diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml index 440853cc61d..14a5dd73076 100644 --- a/.integrated_tests.yaml +++ b/.integrated_tests.yaml @@ -1,6 +1,6 @@ baselines: bucket: geosx - baseline: integratedTests/baseline_integratedTests-pr3202-6636-2821c93 + baseline: integratedTests/baseline_integratedTests-pr3197-6679-a003f43 allow_fail: all: '' diff --git a/BASELINE_NOTES.md b/BASELINE_NOTES.md index dd686fc7877..be5e9d7d934 100644 --- a/BASELINE_NOTES.md +++ b/BASELINE_NOTES.md @@ -6,6 +6,10 @@ This file is designed to track changes to the integrated test baselines. Any developer who updates the baseline ID in the .integrated_tests.yaml file is expected to create an entry in this file with the pull request number, date, and their justification for rebaselining. These notes should be in reverse-chronological order, and use the following time format: (YYYY-MM-DD). +PR 3197 (2024-08-07) +==================== +Separation of contact and friction laws. + PR #3202 (2024-08-03) ====================== Acoustic VTI tests needed rebaselining after update in source and receiver location algorithm. diff --git a/examples/pygeosxExamples/modifyBoundaryCondition/pkn_example.xml b/examples/pygeosxExamples/modifyBoundaryCondition/pkn_example.xml index 18702b7b27f..57f6d330ce2 100644 --- a/examples/pygeosxExamples/modifyBoundaryCondition/pkn_example.xml +++ b/examples/pygeosxExamples/modifyBoundaryCondition/pkn_example.xml @@ -22,7 +22,6 @@ name="hydrofracture" solidSolverName="lagsolve" logLevel="1" - contactRelationName="fractureContact" flowSolverName="SinglePhaseFlow" surfaceGeneratorName="SurfaceGen" targetRegions="{ Domain, Fracture }"> @@ -39,8 +38,8 @@ timeIntegrationOption="QuasiStatic" discretization="FE1" name="lagsolve" - contactRelationName="fractureContact" - targetRegions="{ Domain, Fracture }"> + targetRegions="{ Domain, Fracture }" + contactPenaltyStiffness="1.0e0"> + materialList="{ water, rock, fractureFilling, fracturePorosity, fracturePerm, nullSolid, hApertureTable }"/> @@ -112,9 +111,11 @@ compressibility="0.0"/> + name="fractureContact"/> + + diff --git a/host-configs/apple/macOS_customized.cmake b/host-configs/apple/macOS_customized.cmake deleted file mode 100644 index c1707a0b945..00000000000 --- a/host-configs/apple/macOS_customized.cmake +++ /dev/null @@ -1,7 +0,0 @@ -set( HOMEBREW_DIR $ENV{HOMEBREW_DIR} ) -set( CONFIG_NAME $ENV{GEOS_CONFIG_NAME} ) - -set(Python3_ROOT_DIR $ENV{GEOS_PYTHON_DIR} CACHE PATH "") -set(Python3_EXECUTABLE ${Python3_ROOT_DIR}/bin/python3 CACHE PATH "") - -include(${CMAKE_CURRENT_LIST_DIR}/macOS_base.cmake) diff --git a/inputFiles/almContactMechanics/ALM_stickFault_base.xml b/inputFiles/almContactMechanics/ALM_stickFault_base.xml index 444c2cee4c7..7a8c79aa7d8 100644 --- a/inputFiles/almContactMechanics/ALM_stickFault_base.xml +++ b/inputFiles/almContactMechanics/ALM_stickFault_base.xml @@ -92,8 +92,11 @@ + frictionCoefficient="0.1"/> + + diff --git a/inputFiles/efemFractureMechanics/EmbFrac_Compression_CoulombFriction_base.xml b/inputFiles/efemFractureMechanics/EmbFrac_Compression_CoulombFriction_base.xml index b6e12ad2dd5..2e2568a0250 100644 --- a/inputFiles/efemFractureMechanics/EmbFrac_Compression_CoulombFriction_base.xml +++ b/inputFiles/efemFractureMechanics/EmbFrac_Compression_CoulombFriction_base.xml @@ -9,7 +9,8 @@ initialDt="10" discretization="FE1" timeIntegrationOption="QuasiStatic" - logLevel="1"> + logLevel="1" + contactPenaltyStiffness="1.0e12"> @@ -60,12 +61,10 @@ defaultShearModulus="1.0e10"/> + frictionCoefficient = "0.8"/> @@ -116,13 +115,6 @@ scale="0.0" setNames="{ zneg, zpos }"/> - - - - + logLevel="1" + contactPenaltyStiffness="1.0e8"> @@ -60,18 +61,9 @@ defaultShearModulus="1.0e10"/> + name="frictionLaw"/> - - - - + useStaticCondensation="1" + contactPenaltyStiffness="0.0e8"> + useStaticCondensation="1" + contactPenaltyStiffness="0.0e8"> + useStaticCondensation="1" + contactPenaltyStiffness="0.0e8"> + logLevel="1" + contactPenaltyStiffness="0.0e8"> @@ -85,19 +86,10 @@ defaultShearModulus="1.0e10"/> + name="frictionLaw"/> - - - - @@ -35,19 +35,10 @@ defaultShearModulus="1.0e10"/> + name="frictionLaw"/> - - - - + logLevel="1" + contactPenaltyStiffness="0.0e8"> + logLevel="1" + contactPenaltyStiffness="0.0e8"> + useStaticCondensation="1" + contactPenaltyStiffness="0.0e8"> + useStaticCondensation="1" + contactPenaltyStiffness="0.0e8"> + logLevel="1" + contactPenaltyStiffness="0.0e8"> + contactRelationName="fractureContact" + contactPenaltyStiffness="1.0e0"/> + materialList="{ water, fractureFilling, fractureContact, hApertureModel }"/> @@ -126,9 +126,11 @@ name="fracturePerm"/> + name="fractureContact"/> + + diff --git a/inputFiles/hydraulicFracturing/heterogeneousInSitu_base.xml b/inputFiles/hydraulicFracturing/heterogeneousInSitu_base.xml index bbf01a6a976..8c55433af13 100644 --- a/inputFiles/hydraulicFracturing/heterogeneousInSitu_base.xml +++ b/inputFiles/hydraulicFracturing/heterogeneousInSitu_base.xml @@ -11,7 +11,6 @@ surfaceGeneratorName="SurfaceGen" logLevel="1" targetRegions="{ Fracture }" - contactRelationName="fractureContact" maxNumResolves="2"> + contactRelationName="fractureContact" + contactPenaltyStiffness="1e12"> + materialList="{ water, fractureFilling, fractureContact, hApertureModel }"/> @@ -157,9 +157,11 @@ compressibility="0.0"/> + name="fractureContact"/> + + diff --git a/inputFiles/hydraulicFracturing/heterogeneousInSitu_smoke.xml b/inputFiles/hydraulicFracturing/heterogeneousInSitu_smoke.xml index b75048cd2c2..10a8b07477f 100644 --- a/inputFiles/hydraulicFracturing/heterogeneousInSitu_smoke.xml +++ b/inputFiles/hydraulicFracturing/heterogeneousInSitu_smoke.xml @@ -44,7 +44,6 @@ surfaceGeneratorName="SurfaceGen" logLevel="1" targetRegions="{ Fracture }" - contactRelationName="fractureContact" maxNumResolves="2"> + contactRelationName="fractureContact" + contactPenaltyStiffness="0.0e12"> + materialList="{ water, fractureFilling, fractureContact, hApertureModel }"/> @@ -167,9 +167,11 @@ compressibility="0.0"/> + name="fractureContact"/> + + diff --git a/inputFiles/hydraulicFracturing/hydrofractureSinglePhase2d.xml b/inputFiles/hydraulicFracturing/hydrofractureSinglePhase2d.xml index 5d1e865b642..40cb8c20c2c 100644 --- a/inputFiles/hydraulicFracturing/hydrofractureSinglePhase2d.xml +++ b/inputFiles/hydraulicFracturing/hydrofractureSinglePhase2d.xml @@ -9,8 +9,7 @@ solidSolverName="lagsolve" flowSolverName="SinglePhaseFlow" surfaceGeneratorName="SurfaceGen" - targetRegions="{ Fracture }" - contactRelationName="fractureContact"> + targetRegions="{ Fracture }"> + contactRelationName="fractureContact" + contactPenaltyStiffness="0.0e12"/> @@ -160,9 +160,11 @@ name="fracturePerm"/> + name="fractureContact"/> + + diff --git a/inputFiles/hydraulicFracturing/kgdBase_C3D6_base.xml b/inputFiles/hydraulicFracturing/kgdBase_C3D6_base.xml index 19f102b76c5..3efb892cf10 100644 --- a/inputFiles/hydraulicFracturing/kgdBase_C3D6_base.xml +++ b/inputFiles/hydraulicFracturing/kgdBase_C3D6_base.xml @@ -24,7 +24,7 @@ + materialList="{ water, fractureFilling, fractureContact, hApertureModel }"/> @@ -62,10 +62,11 @@ name="fracturePerm"/> + name="fractureContact"/> + diff --git a/inputFiles/hydraulicFracturing/kgdEdgeBased_C3D6_base.xml b/inputFiles/hydraulicFracturing/kgdEdgeBased_C3D6_base.xml index 298684e5927..c775276de50 100644 --- a/inputFiles/hydraulicFracturing/kgdEdgeBased_C3D6_base.xml +++ b/inputFiles/hydraulicFracturing/kgdEdgeBased_C3D6_base.xml @@ -13,8 +13,7 @@ solidSolverName="lagsolve" flowSolverName="SinglePhaseFlow" surfaceGeneratorName="SurfaceGen" - targetRegions="{ Fracture }" - contactRelationName="fractureContact"> + targetRegions="{ Fracture }"> + contactRelationName="fractureContact" + contactPenaltyStiffness="0.0e8"> + targetRegions="{ Fracture }"> + contactRelationName="fractureContact" + contactPenaltyStiffness="0.0e8"/> @@ -30,7 +29,8 @@ timeIntegrationOption="QuasiStatic" discretization="FE1" targetRegions="{ Domain, Fracture }" - contactRelationName="fractureContact"/> + contactRelationName="fractureContact" + contactPenaltyStiffness="1.0"/> @@ -72,7 +72,7 @@ + materialList="{ water, fractureFilling, fractureContact, hApertureModel }"/> @@ -114,9 +114,11 @@ name="fracturePerm"/> + name="fractureContact"/> + + diff --git a/inputFiles/hydraulicFracturing/kgdToughnessDominated_poroelastic_base.xml b/inputFiles/hydraulicFracturing/kgdToughnessDominated_poroelastic_base.xml index 2aba63788b4..9624bd12259 100644 --- a/inputFiles/hydraulicFracturing/kgdToughnessDominated_poroelastic_base.xml +++ b/inputFiles/hydraulicFracturing/kgdToughnessDominated_poroelastic_base.xml @@ -11,7 +11,6 @@ logLevel="1" targetRegions="{ Domain, Fracture }" isMatrixPoroelastic="1" - contactRelationName="fractureContact" maxNumResolves="2"> + contactRelationName="fractureContact" + contactPenaltyStiffness="1.0"/> + materialList="{ water, fractureFilling, fractureContact, hApertureModel }"/> @@ -117,9 +117,11 @@ name="fracturePerm"/> + name="fractureContact"/> + + diff --git a/inputFiles/hydraulicFracturing/kgdValidation_base.xml b/inputFiles/hydraulicFracturing/kgdValidation_base.xml index 59c48a60734..1cb74613cac 100644 --- a/inputFiles/hydraulicFracturing/kgdValidation_base.xml +++ b/inputFiles/hydraulicFracturing/kgdValidation_base.xml @@ -11,7 +11,6 @@ surfaceGeneratorName="SurfaceGen" logLevel="1" targetRegions="{ Fracture }" - contactRelationName="fractureContact" maxNumResolves="2"> + contactRelationName="fractureContact" + contactPenaltyStiffness="1.0"/> @@ -75,7 +75,7 @@ + materialList="{ water, fractureFilling, fractureContact, hApertureModel }"/> @@ -115,9 +115,11 @@ name="fracturePerm"/> + name="fractureContact"/> + + diff --git a/inputFiles/hydraulicFracturing/kgdViscosityDominated_base.xml b/inputFiles/hydraulicFracturing/kgdViscosityDominated_base.xml index 28cb49e72b3..8bc09cdc194 100644 --- a/inputFiles/hydraulicFracturing/kgdViscosityDominated_base.xml +++ b/inputFiles/hydraulicFracturing/kgdViscosityDominated_base.xml @@ -12,7 +12,6 @@ logLevel="1" initialDt="0.1" targetRegions="{ Fracture }" - contactRelationName="fractureContact" maxNumResolves="2"> @@ -30,7 +29,8 @@ timeIntegrationOption="QuasiStatic" discretization="FE1" targetRegions="{ Domain, Fracture }" - contactRelationName="fractureContact"/> + contactRelationName="fractureContact" + contactPenaltyStiffness="1.0e0"/> @@ -73,7 +73,7 @@ + materialList="{ water, fractureFilling, fractureContact, hApertureModel }"/> @@ -115,9 +115,11 @@ name="fracturePerm"/> + name="fractureContact"/> + + diff --git a/inputFiles/hydraulicFracturing/kgdViscosityDominated_poroelastic_base.xml b/inputFiles/hydraulicFracturing/kgdViscosityDominated_poroelastic_base.xml index 8161813edac..c369c1c522f 100644 --- a/inputFiles/hydraulicFracturing/kgdViscosityDominated_poroelastic_base.xml +++ b/inputFiles/hydraulicFracturing/kgdViscosityDominated_poroelastic_base.xml @@ -11,7 +11,6 @@ logLevel="1" targetRegions="{ Domain, Fracture }" isMatrixPoroelastic="1" - contactRelationName="fractureContact" maxNumResolves="2"> + contactRelationName="fractureContact" + contactPenaltyStiffness="1.0e0"/> + materialList="{ water, fractureFilling, fractureContact, hApertureModel }"/> @@ -118,9 +118,11 @@ name="fracturePerm"/> + name="fractureContact"/> + + diff --git a/inputFiles/hydraulicFracturing/pennyShapedToughnessDominated_base.xml b/inputFiles/hydraulicFracturing/pennyShapedToughnessDominated_base.xml index bf1560af82a..726c70762b7 100644 --- a/inputFiles/hydraulicFracturing/pennyShapedToughnessDominated_base.xml +++ b/inputFiles/hydraulicFracturing/pennyShapedToughnessDominated_base.xml @@ -11,7 +11,7 @@ + materialList="{ water, fractureFilling, fractureContact, hApertureModel }"/> @@ -51,9 +51,12 @@ name="fracturePerm"/> + name="fractureContact"/> + + + diff --git a/inputFiles/hydraulicFracturing/pennyShapedToughnessDominated_benchmark.xml b/inputFiles/hydraulicFracturing/pennyShapedToughnessDominated_benchmark.xml index be67d03b8e2..59abb35871a 100644 --- a/inputFiles/hydraulicFracturing/pennyShapedToughnessDominated_benchmark.xml +++ b/inputFiles/hydraulicFracturing/pennyShapedToughnessDominated_benchmark.xml @@ -16,7 +16,6 @@ surfaceGeneratorName="SurfaceGen" logLevel="1" targetRegions="{ Fracture }" - contactRelationName="fractureContact" maxNumResolves="5" initialDt="0.1"> + contactRelationName="fractureContact" + contactPenaltyStiffness="1.0e0"> + materialList="{ water, fractureFilling, fractureContact, hApertureModel }"/> @@ -66,9 +66,11 @@ name="fracturePerm"/> + name="fractureContact"/> + + @@ -82,8 +84,7 @@ + name="singlePhaseTPFA"/> diff --git a/inputFiles/hydraulicFracturing/pennyShapedToughnessDominated_poroelastic_benchmark.xml b/inputFiles/hydraulicFracturing/pennyShapedToughnessDominated_poroelastic_benchmark.xml index dff5dd1605c..e830a52c5fd 100644 --- a/inputFiles/hydraulicFracturing/pennyShapedToughnessDominated_poroelastic_benchmark.xml +++ b/inputFiles/hydraulicFracturing/pennyShapedToughnessDominated_poroelastic_benchmark.xml @@ -17,7 +17,6 @@ logLevel="1" targetRegions="{ Domain, Fracture }" isMatrixPoroelastic="1" - contactRelationName="fractureContact" maxNumResolves="5" initialDt="0.1"> + contactRelationName="fractureContact" + contactPenaltyStiffness="1.0e0"> + contactRelationName="fractureContact" + contactPenaltyStiffness="1.0e0"> + materialList="{ water, fractureFilling, fractureContact, hApertureModel }"/> @@ -51,9 +51,11 @@ name="fracturePerm"/> + name="fractureContact"/> + + diff --git a/inputFiles/hydraulicFracturing/pennyShapedViscosityDominated_benchmark.xml b/inputFiles/hydraulicFracturing/pennyShapedViscosityDominated_benchmark.xml index 06b52ff40af..a8f78adf606 100644 --- a/inputFiles/hydraulicFracturing/pennyShapedViscosityDominated_benchmark.xml +++ b/inputFiles/hydraulicFracturing/pennyShapedViscosityDominated_benchmark.xml @@ -16,7 +16,6 @@ surfaceGeneratorName="SurfaceGen" logLevel="1" targetRegions="{ Fracture }" - contactRelationName="fractureContact" maxNumResolves="1" initialDt="0.1"> + contactRelationName="fractureContact" + contactPenaltyStiffness="1.0e0"> + materialList="{ water, fractureFilling, fractureContact, hApertureModel }"/> @@ -66,9 +66,11 @@ name="fracturePerm"/> + name="fractureContact"/> + + diff --git a/inputFiles/hydraulicFracturing/pennyShapedViscosityDominated_poroelastic_benchmark.xml b/inputFiles/hydraulicFracturing/pennyShapedViscosityDominated_poroelastic_benchmark.xml index cf731e4fe59..0d1b9485772 100644 --- a/inputFiles/hydraulicFracturing/pennyShapedViscosityDominated_poroelastic_benchmark.xml +++ b/inputFiles/hydraulicFracturing/pennyShapedViscosityDominated_poroelastic_benchmark.xml @@ -17,7 +17,6 @@ logLevel="1" targetRegions="{ Domain, Fracture }" isMatrixPoroelastic="1" - contactRelationName="fractureContact" maxNumResolves="1" initialDt="0.1"> + contactRelationName="fractureContact" + contactPenaltyStiffness="1.0e0"> + contactRelationName="fractureContact" + contactPenaltyStiffness="1.0e0"> + contactRelationName="fractureContact" + contactPenaltyStiffness="1.0e0"> + materialList="{ water, fractureFilling, fractureContact, hApertureModel }"/> @@ -51,9 +51,12 @@ name="fracturePerm"/> + name="fractureContact"/> + + + diff --git a/inputFiles/hydraulicFracturing/pknViscosityDominated_benchmark.xml b/inputFiles/hydraulicFracturing/pknViscosityDominated_benchmark.xml index 9b4d96a8d12..2b8dce7a719 100644 --- a/inputFiles/hydraulicFracturing/pknViscosityDominated_benchmark.xml +++ b/inputFiles/hydraulicFracturing/pknViscosityDominated_benchmark.xml @@ -16,7 +16,6 @@ surfaceGeneratorName="SurfaceGen" logLevel="1" targetRegions="{ Fracture }" - contactRelationName="fractureContact" maxNumResolves="5" initialDt="0.1"> + contactRelationName="fractureContact" + contactPenaltyStiffness="1.0e0"> + materialList="{ water, fractureFilling, fractureContact, hApertureModel }"/> @@ -66,9 +66,11 @@ name="fracturePerm"/> + name="fractureContact"/> + + diff --git a/inputFiles/hydraulicFracturing/pknViscosityDominated_poroelastic_benchmark.xml b/inputFiles/hydraulicFracturing/pknViscosityDominated_poroelastic_benchmark.xml index 4dde5d644d2..72b1be2f4d8 100644 --- a/inputFiles/hydraulicFracturing/pknViscosityDominated_poroelastic_benchmark.xml +++ b/inputFiles/hydraulicFracturing/pknViscosityDominated_poroelastic_benchmark.xml @@ -17,7 +17,6 @@ logLevel="1" targetRegions="{ Domain, Fracture }" isMatrixPoroelastic="1" - contactRelationName="fractureContact" maxNumResolves="5" initialDt="0.1"> + contactRelationName="fractureContact" + contactPenaltyStiffness="1.0e0"> + contactRelationName="fractureContact" + contactPenaltyStiffness="1.0e0"> + contactRelationName="fractureContact" + contactPenaltyStiffness="1.0e0"> + materialList="{water, rock, fractureFilling, fractureContact, hApertureModel}"/> @@ -74,9 +74,11 @@ name="fracturePerm"/> + name="fractureContact"/> + + diff --git a/inputFiles/hydraulicFracturing/walshQuarterNoChombo_benchmark.xml b/inputFiles/hydraulicFracturing/walshQuarterNoChombo_benchmark.xml index 0a715d4089a..957b63bd278 100644 --- a/inputFiles/hydraulicFracturing/walshQuarterNoChombo_benchmark.xml +++ b/inputFiles/hydraulicFracturing/walshQuarterNoChombo_benchmark.xml @@ -28,7 +28,6 @@ surfaceGeneratorName="SurfaceGen" logLevel="1" targetRegions="{Fracture}" - contactRelationName="fractureContact" maxNumResolves="5" initialDt="0.1"> diff --git a/inputFiles/hydraulicFracturing/walshQuarterNoChombo_smoke.xml b/inputFiles/hydraulicFracturing/walshQuarterNoChombo_smoke.xml index 3d0d8a5c3e8..9990dd5c291 100644 --- a/inputFiles/hydraulicFracturing/walshQuarterNoChombo_smoke.xml +++ b/inputFiles/hydraulicFracturing/walshQuarterNoChombo_smoke.xml @@ -29,7 +29,6 @@ surfaceGeneratorName="SurfaceGen" logLevel="1" targetRegions="{Fracture}" - contactRelationName="fractureContact" maxNumResolves="5" initialDt="0.1"> diff --git a/inputFiles/lagrangianContactMechanics/ContactMechanics_PEBICrack_base.xml b/inputFiles/lagrangianContactMechanics/ContactMechanics_PEBICrack_base.xml index c08000ed8c1..09667b5e88f 100644 --- a/inputFiles/lagrangianContactMechanics/ContactMechanics_PEBICrack_base.xml +++ b/inputFiles/lagrangianContactMechanics/ContactMechanics_PEBICrack_base.xml @@ -26,7 +26,7 @@ - + diff --git a/inputFiles/lagrangianContactMechanics/ContactMechanics_PassingCrack_base.xml b/inputFiles/lagrangianContactMechanics/ContactMechanics_PassingCrack_base.xml index 58afadaed04..c309080025d 100644 --- a/inputFiles/lagrangianContactMechanics/ContactMechanics_PassingCrack_base.xml +++ b/inputFiles/lagrangianContactMechanics/ContactMechanics_PassingCrack_base.xml @@ -48,8 +48,8 @@ + defaultAperture="1.0e-6" + materialList="{ frictionLaw, rock }"/> @@ -60,9 +60,8 @@ defaultShearModulus="1.0e4"/> @@ -123,13 +122,6 @@ scale="0.0" setNames="{ rear }"/> - - - - + defaultAperture="1.0e-6" + materialList="{ frictionLaw }"/> @@ -62,10 +62,9 @@ defaultShearModulus="2.0e3"/> + frictionCoefficient="0.577350269"/> @@ -156,10 +155,6 @@ inputVarNames="{ time }" coordinates="{ 0.0, 5.0, 10.0 }" values="{ 0.0, 5.e0, -5.e0 }"/> - diff --git a/inputFiles/lagrangianContactMechanics/ContactMechanics_SingleFracCompression_base.xml b/inputFiles/lagrangianContactMechanics/ContactMechanics_SingleFracCompression_base.xml index fa6999ef17a..28b80df0368 100755 --- a/inputFiles/lagrangianContactMechanics/ContactMechanics_SingleFracCompression_base.xml +++ b/inputFiles/lagrangianContactMechanics/ContactMechanics_SingleFracCompression_base.xml @@ -57,8 +57,8 @@ + defaultAperture="1.0e-6" + materialList="{ frictionLaw }"/> @@ -70,21 +70,12 @@ defaultShearModulus="1.0e10"/> + frictionCoefficient="0.577350269"/> - - - - - + defaultAperture="1.0e-6" + materialList="{ frictionLaw }"/> @@ -72,10 +72,9 @@ defaultShearModulus="29.17e9"/> + frictionCoefficient="0.577350269"/> @@ -160,11 +159,6 @@ inputVarNames="{ time }" coordinates="{ 0.0, 1.0 }" values="{ 0.0, 1.e0 }"/> - - diff --git a/inputFiles/lagrangianContactMechanics/ContactMechanics_UnstructuredCrack_base.xml b/inputFiles/lagrangianContactMechanics/ContactMechanics_UnstructuredCrack_base.xml index c47369c9a73..4c4c2bb1379 100644 --- a/inputFiles/lagrangianContactMechanics/ContactMechanics_UnstructuredCrack_base.xml +++ b/inputFiles/lagrangianContactMechanics/ContactMechanics_UnstructuredCrack_base.xml @@ -48,8 +48,8 @@ + defaultAperture="1.0e-6" + materialList="{ frictionLaw }"/> @@ -60,10 +60,9 @@ defaultShearModulus="1.0e4"/> + frictionCoefficient="0.577350269"/> @@ -137,12 +136,7 @@ name="ForceTimeFunction" inputVarNames="{ time }" coordinates="{ 0.0, 2.0 }" - values="{ 0.0, 2.e0 }"/> - - + values="{ 0.0, 2.e0 }"/> diff --git a/inputFiles/lagrangianContactMechanics/ContactMechanics_slippingFault_base.xml b/inputFiles/lagrangianContactMechanics/ContactMechanics_slippingFault_base.xml index 7abd1e6c74b..0c5350192c1 100644 --- a/inputFiles/lagrangianContactMechanics/ContactMechanics_slippingFault_base.xml +++ b/inputFiles/lagrangianContactMechanics/ContactMechanics_slippingFault_base.xml @@ -82,8 +82,8 @@ + defaultAperture="1.0e-6" + materialList="{ frictionLaw, rock }"/> @@ -95,18 +95,12 @@ defaultShearModulus="9e9"/> + frictionCoefficient="0.1"/> - - + defaultAperture="1.0e-6" + materialList="{ frictionLaw }"/> @@ -58,10 +58,9 @@ defaultShearModulus="1.0e10"/> + frictionCoefficient="0.577350269"/> @@ -115,13 +114,6 @@ setNames="{ core }"/> - - - - diff --git a/inputFiles/poromechanicsFractures/ExponentialDecayPermeability_conformingFracture_base.xml b/inputFiles/poromechanicsFractures/ExponentialDecayPermeability_conformingFracture_base.xml index 0a1704fd9c5..7c53d0d4981 100644 --- a/inputFiles/poromechanicsFractures/ExponentialDecayPermeability_conformingFracture_base.xml +++ b/inputFiles/poromechanicsFractures/ExponentialDecayPermeability_conformingFracture_base.xml @@ -88,7 +88,7 @@ @@ -146,8 +146,11 @@ + frictionCoefficient="0.01"/> + + diff --git a/inputFiles/poromechanicsFractures/ExponentialDecayPermeability_edfm_base.xml b/inputFiles/poromechanicsFractures/ExponentialDecayPermeability_edfm_base.xml index f43116ee181..480796bc952 100755 --- a/inputFiles/poromechanicsFractures/ExponentialDecayPermeability_edfm_base.xml +++ b/inputFiles/poromechanicsFractures/ExponentialDecayPermeability_edfm_base.xml @@ -22,13 +22,13 @@ name="fractureMechSolver" targetRegions="{ Domain, Fracture }" timeIntegrationOption="QuasiStatic" - discretization="FE1"/> + discretization="FE1" + contactPenaltyStiffness="2.0e11"/> + targetRegions="{ Domain, Fracture }"/> @@ -120,9 +120,11 @@ initialPermeability="{1.0e-12, 1.0e-12, 1.0e-12}"/> + name="frictionLaw"/> + + diff --git a/inputFiles/poromechanicsFractures/PoroElastic_conformingFracture_2d_faultSlip_base.xml b/inputFiles/poromechanicsFractures/PoroElastic_conformingFracture_2d_faultSlip_base.xml index 4c674d7bacc..390f1cec5a2 100644 --- a/inputFiles/poromechanicsFractures/PoroElastic_conformingFracture_2d_faultSlip_base.xml +++ b/inputFiles/poromechanicsFractures/PoroElastic_conformingFracture_2d_faultSlip_base.xml @@ -31,7 +31,7 @@ @@ -87,8 +87,11 @@ + frictionCoefficient="0.1"/> + + diff --git a/inputFiles/poromechanicsFractures/PoroElastic_conformingFracture_2d_openingFrac_base.xml b/inputFiles/poromechanicsFractures/PoroElastic_conformingFracture_2d_openingFrac_base.xml index 9aa53c702f6..29509e3a858 100644 --- a/inputFiles/poromechanicsFractures/PoroElastic_conformingFracture_2d_openingFrac_base.xml +++ b/inputFiles/poromechanicsFractures/PoroElastic_conformingFracture_2d_openingFrac_base.xml @@ -24,7 +24,7 @@ @@ -78,10 +78,13 @@ name="fracturePerm"/> + frictionCoefficient="0.577350269"/> + + diff --git a/inputFiles/poromechanicsFractures/PoroElastic_dfm_PEBICrack_base.xml b/inputFiles/poromechanicsFractures/PoroElastic_dfm_PEBICrack_base.xml index 5f5aa06d519..1aa468ae144 100644 --- a/inputFiles/poromechanicsFractures/PoroElastic_dfm_PEBICrack_base.xml +++ b/inputFiles/poromechanicsFractures/PoroElastic_dfm_PEBICrack_base.xml @@ -11,7 +11,7 @@ @@ -118,10 +118,13 @@ name="fracturePerm"/> + frictionCoefficient="0.01"/> + + diff --git a/inputFiles/poromechanicsFractures/PoroElastic_efem-edfm_base.xml b/inputFiles/poromechanicsFractures/PoroElastic_efem-edfm_base.xml index 8242a94caed..3d8a8b10309 100644 --- a/inputFiles/poromechanicsFractures/PoroElastic_efem-edfm_base.xml +++ b/inputFiles/poromechanicsFractures/PoroElastic_efem-edfm_base.xml @@ -23,7 +23,8 @@ name="fractureMechSolver" targetRegions="{ Domain, Fracture }" timeIntegrationOption="QuasiStatic" - discretization="FE1"/> + discretization="FE1" + contactPenaltyStiffness="4.0e9"/> @@ -120,11 +121,15 @@ name="fracturePerm"/> + name="fractureContact"/> + + + + + discretization="FE1" + contactPenaltyStiffness="1e12"> @@ -247,9 +248,11 @@ name="fracturePerm"/> + name="fractureContact"/> + + diff --git a/inputFiles/poromechanicsFractures/PoroElastic_efem-edfm_eggModel_small.xml b/inputFiles/poromechanicsFractures/PoroElastic_efem-edfm_eggModel_small.xml index a758775ec56..074d11adb4e 100644 --- a/inputFiles/poromechanicsFractures/PoroElastic_efem-edfm_eggModel_small.xml +++ b/inputFiles/poromechanicsFractures/PoroElastic_efem-edfm_eggModel_small.xml @@ -25,7 +25,8 @@ name="fractureMechSolver" targetRegions="{ Reservoir, Overburden, Underburden, Sideburden, Fracture }" timeIntegrationOption="QuasiStatic" - discretization="FE1"> + discretization="FE1" + contactPenaltyStiffness="1e12"> @@ -247,9 +248,11 @@ name="fracturePerm"/> + name="fractureContact"/> + + diff --git a/inputFiles/poromechanicsFractures/PoroElastic_efem-edfm_pennyCrack_base.xml b/inputFiles/poromechanicsFractures/PoroElastic_efem-edfm_pennyCrack_base.xml index c17733c43ed..7e4c4ad9393 100644 --- a/inputFiles/poromechanicsFractures/PoroElastic_efem-edfm_pennyCrack_base.xml +++ b/inputFiles/poromechanicsFractures/PoroElastic_efem-edfm_pennyCrack_base.xml @@ -38,7 +38,7 @@ name="Fracture" faceBlock="embeddedSurfaceSubRegion" subRegionType="embeddedElement" - materialList="{ water, fractureFilling, fractureContact }" + materialList="{ water, fractureFilling, fractureContact, hApertureModel }" defaultAperture="1e-3"/> @@ -92,9 +92,11 @@ name="fracturePerm"/> + name="fractureContact"/> + + diff --git a/inputFiles/poromechanicsFractures/PoroElastic_efem-edfm_pennyCrack_benchmark.xml b/inputFiles/poromechanicsFractures/PoroElastic_efem-edfm_pennyCrack_benchmark.xml index a2979a5b39a..3ca8d754e85 100644 --- a/inputFiles/poromechanicsFractures/PoroElastic_efem-edfm_pennyCrack_benchmark.xml +++ b/inputFiles/poromechanicsFractures/PoroElastic_efem-edfm_pennyCrack_benchmark.xml @@ -28,7 +28,8 @@ name="fractureMechSolver" targetRegions="{ Domain, Fracture }" timeIntegrationOption="QuasiStatic" - discretization="FE1"/> + discretization="FE1" + contactPenaltyStiffness="1e12"/> + discretization="FE1" + contactPenaltyStiffness="1e12"/> + discretization="FE1" + contactPenaltyStiffness="0.0e8"/> @@ -163,9 +164,11 @@ name="fracturePerm"/> + name="fractureContact"/> + + diff --git a/inputFiles/poromechanicsFractures/SlipPermeability_embeddedFrac.xml b/inputFiles/poromechanicsFractures/SlipPermeability_embeddedFrac.xml index 81aa4adae38..ce28512cb48 100755 --- a/inputFiles/poromechanicsFractures/SlipPermeability_embeddedFrac.xml +++ b/inputFiles/poromechanicsFractures/SlipPermeability_embeddedFrac.xml @@ -22,7 +22,8 @@ name="fractureMechSolver" targetRegions="{ Domain, Fracture }" timeIntegrationOption="QuasiStatic" - discretization="FE1"/> + discretization="FE1" + contactPenaltyStiffness="2.0e11"/> @@ -237,9 +238,11 @@ name="fracturePerm2"/> + name="fractureContact"/> + + diff --git a/inputFiles/poromechanicsFractures/SlipPermeability_pEDFM_base.xml b/inputFiles/poromechanicsFractures/SlipPermeability_pEDFM_base.xml index 5d915f898c5..8be34031377 100755 --- a/inputFiles/poromechanicsFractures/SlipPermeability_pEDFM_base.xml +++ b/inputFiles/poromechanicsFractures/SlipPermeability_pEDFM_base.xml @@ -22,7 +22,8 @@ name="fractureMechSolver" targetRegions="{ Domain, Fracture }" timeIntegrationOption="QuasiStatic" - discretization="FE1"/> + discretization="FE1" + contactPenaltyStiffness="2.0e11"/> + usePEDFM="1"/> @@ -64,7 +64,7 @@ name="Fracture" faceBlock="embeddedSurfaceSubRegion" subRegionType="embeddedElement" - materialList="{ water, fractureFilling, fractureContact }" + materialList="{ water, fractureFilling, fractureContact, hApertureModel}" defaultAperture="1e-3"/> @@ -121,9 +121,11 @@ initialPermeability="{1e-15, 1e-15, 1e-15}"/> + name="fractureContact"/> + + diff --git a/inputFiles/poromechanicsFractures/WillisRichardsPermeability_efem-edfm_base.xml b/inputFiles/poromechanicsFractures/WillisRichardsPermeability_efem-edfm_base.xml index 826edbe9e27..f5fee5d296b 100755 --- a/inputFiles/poromechanicsFractures/WillisRichardsPermeability_efem-edfm_base.xml +++ b/inputFiles/poromechanicsFractures/WillisRichardsPermeability_efem-edfm_base.xml @@ -22,7 +22,8 @@ name="fractureMechSolver" targetRegions="{ Domain, Fracture }" timeIntegrationOption="QuasiStatic" - discretization="FE1"/> + discretization="FE1" + contactPenaltyStiffness="2.0e11"/> @@ -123,9 +124,11 @@ refClosureStress="1.0e7"/> + name="fractureContact"/> + + diff --git a/inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_base.xml b/inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_base.xml index bd3bba82795..0623d767f54 100644 --- a/inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_base.xml +++ b/inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_base.xml @@ -49,8 +49,11 @@ + frictionCoefficient="0.01"/> + + @@ -147,7 +150,7 @@ diff --git a/inputFiles/thermoPoromechanicsFractures/ThermoPoroElastic_base.xml b/inputFiles/thermoPoromechanicsFractures/ThermoPoroElastic_base.xml index cb00829f403..1ba0e49f27f 100644 --- a/inputFiles/thermoPoromechanicsFractures/ThermoPoroElastic_base.xml +++ b/inputFiles/thermoPoromechanicsFractures/ThermoPoroElastic_base.xml @@ -66,9 +66,11 @@ name="fracturePerm"/> + name="fractureContact"/> + + diff --git a/inputFiles/thermoPoromechanicsFractures/ThermoPoroElastic_conforming_base.xml b/inputFiles/thermoPoromechanicsFractures/ThermoPoroElastic_conforming_base.xml index e8ca36eaade..6f855810a38 100644 --- a/inputFiles/thermoPoromechanicsFractures/ThermoPoroElastic_conforming_base.xml +++ b/inputFiles/thermoPoromechanicsFractures/ThermoPoroElastic_conforming_base.xml @@ -74,7 +74,7 @@ diff --git a/inputFiles/thermoPoromechanicsFractures/ThermoPoroElastic_efem-edfm_base.xml b/inputFiles/thermoPoromechanicsFractures/ThermoPoroElastic_efem-edfm_base.xml index 3d8cf3e4a03..47df41a1533 100644 --- a/inputFiles/thermoPoromechanicsFractures/ThermoPoroElastic_efem-edfm_base.xml +++ b/inputFiles/thermoPoromechanicsFractures/ThermoPoroElastic_efem-edfm_base.xml @@ -29,7 +29,8 @@ name="fractureMechSolver" targetRegions="{ RockMatrix, Fracture }" timeIntegrationOption="QuasiStatic" - discretization="FE1"/> + discretization="FE1" + contactPenaltyStiffness="4.0e9"/> diff --git a/inputFiles/thermoPoromechanicsFractures/ThermoPoroElastic_efem-edfm_eggModel_small.xml b/inputFiles/thermoPoromechanicsFractures/ThermoPoroElastic_efem-edfm_eggModel_small.xml index f268a9c3ee1..d2b84ef7de2 100644 --- a/inputFiles/thermoPoromechanicsFractures/ThermoPoroElastic_efem-edfm_eggModel_small.xml +++ b/inputFiles/thermoPoromechanicsFractures/ThermoPoroElastic_efem-edfm_eggModel_small.xml @@ -27,7 +27,8 @@ name="fractureMechSolver" targetRegions="{ Reservoir, Overburden, Underburden, Sideburden, Fracture }" timeIntegrationOption="QuasiStatic" - discretization="FE1"> + discretization="FE1" + contactPenaltyStiffness="0.0e8"> + name="fractureContact"/> + + diff --git a/inputFiles/wellbore/CasedElasticWellbore_ImperfectInterfaces_base.xml b/inputFiles/wellbore/CasedElasticWellbore_ImperfectInterfaces_base.xml index afb8cd86e99..0f782918ec9 100644 --- a/inputFiles/wellbore/CasedElasticWellbore_ImperfectInterfaces_base.xml +++ b/inputFiles/wellbore/CasedElasticWellbore_ImperfectInterfaces_base.xml @@ -74,9 +74,7 @@ + frictionCoefficient="0.5"/> - diff --git a/inputFiles/wellbore/CasedThermoElasticWellbore_ImperfectInterfaces_base.xml b/inputFiles/wellbore/CasedThermoElasticWellbore_ImperfectInterfaces_base.xml index 3a6d0653f9e..4cc92cc6c6c 100644 --- a/inputFiles/wellbore/CasedThermoElasticWellbore_ImperfectInterfaces_base.xml +++ b/inputFiles/wellbore/CasedThermoElasticWellbore_ImperfectInterfaces_base.xml @@ -78,7 +78,7 @@ @@ -100,11 +100,13 @@ + frictionCoefficient="0.5"/> + + 0 + */ + GEOS_HOST_DEVICE + real64 computeHydraulicAperture( real64 const aperture, + real64 const normalTraction, + real64 & dHydraulicAperture_aperture, + real64 & dHydraulicAperture_dNormalStress ) const; + +private: + real64 m_aperture0; + + real64 m_referenceNormalStress; +}; + + +/** + * @class BartonBandis + * + * This class serves as the interface for implementing contact enforcement constitutive relations. + * This does not include the actual enforcement algorithm, but only the constitutive relations that + * govern the behavior of the contact. So things like penalty, or friction, or kinematic constraint. + */ +class BartonBandis : public HydraulicApertureBase +{ +public: + + /** + * @brief The standard data repository constructor + * @param name The name of the relation in the data repository + * @param parent The name of the parent Group that holds this relation object. + */ + BartonBandis( string const & name, + Group * const parent ); + + /** + * @brief default destructor + */ + virtual ~BartonBandis() override; + + + /// Type of kernel wrapper for in-kernel update + using KernelWrapper = BartonBandisUpdates; + + /** + * @brief Create an update kernel wrapper. + * @return the wrapper + */ + KernelWrapper createKernelWrapper() const; + +private: + + struct viewKeyStruct : public HydraulicApertureBase::viewKeyStruct + { + /// string/key for reference normal stress + static constexpr char const * referenceNormalStressString() { return "referenceNormalStress"; } + }; + /// Reference normal stress + real64 m_referenceNormalStress; +}; + +GEOS_HOST_DEVICE +GEOS_FORCE_INLINE +real64 BartonBandisUpdates::computeHydraulicAperture( real64 const aperture, + real64 const normalTraction, + real64 & dHydraulicAperture_aperture, + real64 & dHydraulicAperture_dNormalStress ) const +{ + real64 const hydraulicAperture = ( aperture >= 0.0 ) ? (aperture + m_aperture0) : m_aperture0 / ( 1 + 9*normalTraction/m_referenceNormalStress ); + dHydraulicAperture_dNormalStress = ( aperture >= 0.0 ) ? 0.0 : -hydraulicAperture / ( 1 + 9*normalTraction/m_referenceNormalStress ) * 9/m_referenceNormalStress; + dHydraulicAperture_aperture = ( aperture >= 0.0 ) ? 1.0 : 0.0; + + return hydraulicAperture; ///It would be nice to change this to return a tuple. +} + +} /* namespace constitutive */ + +} /* namespace geos */ + +#endif /* GEOS_CONSTITUTIVE_CONTACT_HYDRAULICAPERTURETABLE_HPP_ */ diff --git a/src/coreComponents/constitutive/contact/ContactBase.hpp b/src/coreComponents/constitutive/contact/ContactBase.hpp deleted file mode 100644 index 7b088ed5e17..00000000000 --- a/src/coreComponents/constitutive/contact/ContactBase.hpp +++ /dev/null @@ -1,269 +0,0 @@ -/* - * ------------------------------------------------------------------------------------------------------------ - * SPDX-License-Identifier: LGPL-2.1-only - * - * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC - * Copyright (c) 2018-2024 Total, S.A - * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University - * Copyright (c) 2018-2024 Chevron - * Copyright (c) 2019- GEOS/GEOSX Contributors - * All rights reserved - * - * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. - * ------------------------------------------------------------------------------------------------------------ - */ - -/** - * @file ContactBase.hpp - */ - -#ifndef GEOS_CONSTITUTIVE_CONTACT_CONTACTBASE_HPP_ -#define GEOS_CONSTITUTIVE_CONTACT_CONTACTBASE_HPP_ - -#include "constitutive/ConstitutiveBase.hpp" -#include "functions/TableFunction.hpp" -#include "physicsSolvers/contact/ContactFields.hpp" - - -namespace geos -{ - -namespace constitutive -{ - -/** - * @class ContactBaseUpdates - * - * This class is used for in-kernel contact relation updates - */ -class ContactBaseUpdates -{ -public: - - ContactBaseUpdates( real64 const & penaltyStiffness, - real64 const & shearStiffness, - real64 const & displacementJumpThreshold, - TableFunction const & apertureTable ) - : m_penaltyStiffness( penaltyStiffness ), - m_shearStiffness( shearStiffness ), - m_displacementJumpThreshold( displacementJumpThreshold ), - m_apertureTable( apertureTable.createKernelWrapper() ) - {} - - /// Default copy constructor - ContactBaseUpdates( ContactBaseUpdates const & ) = default; - - /// Default move constructor - ContactBaseUpdates( ContactBaseUpdates && ) = default; - - /// Deleted default constructor - ContactBaseUpdates() = default; - - /// Deleted copy assignment operator - ContactBaseUpdates & operator=( ContactBaseUpdates const & ) = delete; - - /// Deleted move assignment operator - ContactBaseUpdates & operator=( ContactBaseUpdates && ) = delete; - - /** - * @brief Evaluate the effective aperture, and its derivative wrt aperture - * @param[in] aperture the model aperture/gap - * @param[out] dHydraulicAperture_dAperture the derivative of the effective aperture wrt aperture - * @return The hydraulic aperture that is always > 0 - */ - GEOS_HOST_DEVICE - virtual real64 computeHydraulicAperture( real64 const aperture, - real64 & dHydraulicAperture_dAperture ) const; - - - /** - * @brief Evaluate the traction vector and its derivatives wrt to pressure and jump - * @param[in] dispJump the displacement jump - * @param[in] fractureState the fracture state - * @param[out] tractionVector the traction vector - * @param[out] dTractionVector_dJump the derivative of the traction vector wrt displacement jump - */ - GEOS_HOST_DEVICE - inline - virtual void computeTraction( localIndex const k, - arraySlice1d< real64 const > const & oldDispJump, - arraySlice1d< real64 const > const & dispJump, - integer const & fractureState, - arraySlice1d< real64 > const & tractionVector, - arraySlice2d< real64 > const & dTractionVector_dJump ) const - {GEOS_UNUSED_VAR( k, oldDispJump, dispJump, tractionVector, dTractionVector_dJump, fractureState );} - - /** - * @brief Evaluate the traction vector and its derivatives wrt to pressure and jump - * @param[in] dispJump the displacement jump - * @param[in] tractionVector the traction vector - * @param[out] fractureState the fracture state - */ - GEOS_HOST_DEVICE - inline - virtual void updateFractureState( localIndex const k, - arraySlice1d< real64 const > const & dispJump, - arraySlice1d< real64 const > const & tractionVector, - integer & fractureState ) const - { GEOS_UNUSED_VAR( k, dispJump, tractionVector, fractureState ); } - - /** - * @brief Update the traction with the pressure term - * @param[in] pressure the pressure term - * @param[in] isOpen a flag specifying whether the fracture is open or closed - * @param[inout] traction the current tractionVector - * @param[out] dTraction_dPressure the derivative of the fist component of traction wrt pressure - * @return the updated traction - */ - GEOS_HOST_DEVICE - virtual void addPressureToTraction( real64 const & pressure, - arraySlice1d< real64 >const & tractionVector, - real64 & dTraction_dPressure ) const; - - /** - * @brief Evaluate the limit tangential traction norm and return the derivative wrt normal traction - * @param[in] normalTraction the normal traction - * @param[out] dLimitTangentialTractionNorm_dTraction the derivative of the limit tangential traction norm wrt normal traction - * @return the limit tangential traction norm - */ - GEOS_HOST_DEVICE - inline - virtual real64 computeLimitTangentialTractionNorm( real64 const & normalTraction, - real64 & dLimitTangentialTractionNorm_dTraction ) const - { GEOS_UNUSED_VAR( normalTraction, dLimitTangentialTractionNorm_dTraction ); return 0; }; - -protected: - - /// The penalty stiffness - real64 m_penaltyStiffness; - - /// The shear stiffness - real64 m_shearStiffness; - - /// A threshold valued to determine whether a fracture is open or not. - real64 m_displacementJumpThreshold; - - /// The aperture table function wrapper - TableFunction::KernelWrapper m_apertureTable; -}; - - -/** - * @class ContactBase - * - * This class serves as the interface for implementing contact enforcement constitutive relations. - * This does not include the actual enforcement algorithm, but only the constitutive relations that - * govern the behavior of the contact. So things like penalty, or friction, or kinematic constraint. - */ -class ContactBase : public ConstitutiveBase -{ -public: - - /** - * @brief The standard data repository constructor - * @param name The name of the relation in the data repository - * @param parent The name of the parent Group that holds this relation object. - */ - ContactBase( string const & name, - Group * const parent ); - - /** - * @brief default destructor - */ - virtual ~ContactBase() override; - - virtual void allocateConstitutiveData( dataRepository::Group & parent, - localIndex const numConstitutivePointsPerParentIndex ) override; - - /** - * @brief accessor for penalty stiffness - * @return the stiffness - */ - real64 stiffness() const { return m_penaltyStiffness; } - - - /// Type of kernel wrapper for in-kernel update - using KernelWrapper = ContactBaseUpdates; - - /** - * @brief Create an update kernel wrapper. - * @return the wrapper - */ - KernelWrapper createKernelWrapper() const; - - /** - * @struct Structure to hold scoped key names - */ - struct viewKeyStruct : public ConstitutiveBase::viewKeyStruct - { - /// string/key for penalty stiffness - static constexpr char const * penaltyStiffnessString() { return "penaltyStiffness"; } - - /// string/key for shear stiffness - static constexpr char const * shearStiffnessString() { return "shearStiffness"; } - - /// string/key for the displacement jump threshold value - static constexpr char const * displacementJumpThresholdString() { return "displacementJumpThreshold"; } - - /// string/key for aperture tolerance - static constexpr char const * apertureToleranceString() { return "apertureTolerance"; } - - /// string/key for aperture table name - static constexpr char const * apertureTableNameString() { return "apertureTableName"; } - }; - -protected: - virtual void postInputInitialization() override; - - virtual void initializePreSubGroups() override; - - /** - * @brief Validate the values provided in the aperture table - * @param[in] apertureTable the effective aperture vs aperture table - */ - void validateApertureTable( TableFunction const & apertureTable ) const; - - - /// The value of penalty to penetration - real64 m_penaltyStiffness; - - /// The value of the shear stiffness - real64 m_shearStiffness; - - /// The aperture tolerance to avoid floating point errors in expressions involving aperture - real64 m_apertureTolerance; - - /// The name of the aperture table, if any - string m_apertureTableName; - - /// A threshold valued to determine whether a fracture is open or not. - real64 m_displacementJumpThreshold; - - /// Pointer to the function that limits the model aperture to a physically admissible value. - TableFunction const * m_apertureTable; -}; - -GEOS_HOST_DEVICE -GEOS_FORCE_INLINE -real64 ContactBaseUpdates::computeHydraulicAperture( real64 const aperture, - real64 & dHydraulicAperture_dAperture ) const -{ - return m_apertureTable.compute( &aperture, &dHydraulicAperture_dAperture ); -} - -GEOS_HOST_DEVICE -GEOS_FORCE_INLINE -void ContactBaseUpdates::addPressureToTraction( real64 const & pressure, - arraySlice1d< real64 > const & tractionVector, - real64 & dTraction_dPressure ) const -{ - tractionVector[0] -= pressure; - dTraction_dPressure = -1.0; -} - - -} /* namespace constitutive */ - -} /* namespace geos */ - -#endif /* GEOS_CONSTITUTIVE_CONTACT_CONTACTBASE_HPP_ */ diff --git a/src/coreComponents/constitutive/contact/CoulombContact.hpp b/src/coreComponents/constitutive/contact/CoulombContact.hpp deleted file mode 100644 index 89b39211b0d..00000000000 --- a/src/coreComponents/constitutive/contact/CoulombContact.hpp +++ /dev/null @@ -1,328 +0,0 @@ -/* - * ------------------------------------------------------------------------------------------------------------ - * SPDX-License-Identifier: LGPL-2.1-only - * - * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC - * Copyright (c) 2018-2024 Total, S.A - * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University - * Copyright (c) 2018-2024 Chevron - * Copyright (c) 2019- GEOS/GEOSX Contributors - * All rights reserved - * - * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. - * ------------------------------------------------------------------------------------------------------------ - */ - -/** - * @file CoulombContact.hpp - */ - -#ifndef GEOS_CONSTITUTIVE_CONTACT_COULOMBCONTACT_HPP_ -#define GEOS_CONSTITUTIVE_CONTACT_COULOMBCONTACT_HPP_ - -#include "ContactBase.hpp" - -namespace geos -{ - -namespace constitutive -{ - -/** - * @class CoulombContactUpdates - * - * This class is used for in-kernel contact relation updates - */ -class CoulombContactUpdates : public ContactBaseUpdates -{ -public: - - CoulombContactUpdates( real64 const & penaltyStiffness, - real64 const & shearStiffness, - real64 const & displacementJumpThreshold, - TableFunction const & apertureTable, - real64 const & cohesion, - real64 const & frictionCoefficient, - arrayView2d< real64 > const & elasticSlip ) - : ContactBaseUpdates( penaltyStiffness, shearStiffness, displacementJumpThreshold, apertureTable ), - m_cohesion( cohesion ), - m_frictionCoefficient( frictionCoefficient ), - m_elasticSlip( elasticSlip ) - {} - - /// Default copy constructor - CoulombContactUpdates( CoulombContactUpdates const & ) = default; - - /// Default move constructor - CoulombContactUpdates( CoulombContactUpdates && ) = default; - - /// Deleted default constructor - CoulombContactUpdates() = delete; - - /// Deleted copy assignment operator - CoulombContactUpdates & operator=( CoulombContactUpdates const & ) = delete; - - /// Deleted move assignment operator - CoulombContactUpdates & operator=( CoulombContactUpdates && ) = delete; - - /** - * @brief Evaluate the limit tangential traction norm and return the derivative wrt normal traction - * @param[in] normalTraction the normal traction - * @param[out] dLimitTangentialTractionNorm_dTraction the derivative of the limit tangential traction norm wrt normal traction - * @return the limit tangential traction norm - */ - GEOS_HOST_DEVICE - inline - virtual real64 computeLimitTangentialTractionNorm( real64 const & normalTraction, - real64 & dLimitTangentialTractionNorm_dTraction ) const override final; - - GEOS_HOST_DEVICE - inline - virtual void computeTraction( localIndex const k, - arraySlice1d< real64 const > const & oldDispJump, - arraySlice1d< real64 const > const & dispJump, - integer const & fractureState, - arraySlice1d< real64 > const & tractionVector, - arraySlice2d< real64 > const & dTractionVector_dJump ) const override final; - - GEOS_HOST_DEVICE - inline - virtual void updateFractureState( localIndex const k, - arraySlice1d< real64 const > const & dispJump, - arraySlice1d< real64 const > const & tractionVector, - integer & fractureState ) const override final; - -private: - - /// The cohesion for each upper level dimension (i.e. cell) of *this - real64 m_cohesion; - - /// The friction coefficient for each upper level dimension (i.e. cell) of *this - real64 m_frictionCoefficient; - - arrayView2d< real64 > m_elasticSlip; -}; - - -/** - * @class CoulombContact - * - * Class to provide a CoulombContact friction model. - */ -class CoulombContact : public ContactBase -{ -public: - - /** - * constructor - * @param[in] name name of the instance in the catalog - * @param[in] parent the group which contains this instance - */ - CoulombContact( string const & name, Group * const parent ); - - /** - * Default Destructor - */ - virtual ~CoulombContact() override; - - /** - * @name Static Factory Catalog members and functions - */ - ///@{ - - /** - * @return A string that is used to register/lookup this class in the registry - */ - static string catalogName() { return "Coulomb"; } - - virtual string getCatalogName() const override { return catalogName(); } - - ///@} - - virtual void allocateConstitutiveData( dataRepository::Group & parent, - localIndex const numConstitutivePointsPerParentIndex ) override final; - - /** - * @brief Const accessor for cohesion - * @return A const reference to arrayView1d containing the - * cohesions (at every element). - */ - real64 const & cohesion() const { return m_cohesion; } - - /** - * @brief Const accessor for friction angle - * @return A const reference to arrayView1d containing the - * friction coefficient (at every element). - */ - real64 const & frictionCoefficient() const { return m_frictionCoefficient; } - - /// Type of kernel wrapper for in-kernel update - using KernelWrapper = CoulombContactUpdates; - - /** - * @brief Create an update kernel wrapper. - * @return the wrapper - */ - KernelWrapper createKernelWrapper() const; - -protected: - - virtual void postInputInitialization() override; - -private: - - /// The cohesion for each upper level dimension (i.e. cell) of *this - real64 m_cohesion; - - /// The friction coefficient for each upper level dimension (i.e. cell) of *this - real64 m_frictionCoefficient; - - /// Elastic slip - array2d< real64 > m_elasticSlip; - -/** - * @struct Set of "char const *" and keys for data specified in this class. - */ - struct viewKeyStruct : public ContactBase::viewKeyStruct - { - /// string/key for cohesion - static constexpr char const * cohesionString() { return "cohesion"; } - - /// string/key for friction coefficient - static constexpr char const * frictionCoefficientString() { return "frictionCoefficient"; } - - /// string/key for the elastic slip - static constexpr char const * elasticSlipString() { return "elasticSlip"; } - }; - -}; - - -GEOS_HOST_DEVICE -real64 CoulombContactUpdates::computeLimitTangentialTractionNorm( real64 const & normalTraction, - real64 & dLimitTangentialTractionNorm_dTraction ) const -{ - dLimitTangentialTractionNorm_dTraction = m_frictionCoefficient; - return ( m_cohesion - normalTraction * m_frictionCoefficient ); -} - - -GEOS_HOST_DEVICE -inline void CoulombContactUpdates::computeTraction( localIndex const k, - arraySlice1d< real64 const > const & oldDispJump, - arraySlice1d< real64 const > const & dispJump, - integer const & fractureState, - arraySlice1d< real64 > const & tractionVector, - arraySlice2d< real64 > const & dTractionVector_dJump ) const -{ - - bool const isOpen = fractureState == fields::contact::FractureState::Open; - - // Initialize everyting to 0 - tractionVector[0] = 0.0; - tractionVector[1] = 0.0; - tractionVector[2] = 0.0; - LvArray::forValuesInSlice( dTractionVector_dJump, []( real64 & val ){ val = 0.0; } ); - // If the fracture is open the traction is 0 and so are its derivatives so there is nothing to do - if( !isOpen ) - { - // normal component of the traction - tractionVector[0] = m_penaltyStiffness * dispJump[0]; - // derivative of the normal component w.r.t. to the nomral dispJump - dTractionVector_dJump[0][0] = m_penaltyStiffness; - - // Compute the slip - real64 const slip[2] = { dispJump[1] - oldDispJump[1], - dispJump[2] - oldDispJump[2] }; - - - real64 const tau[2] = { m_shearStiffness * ( slip[0] + m_elasticSlip[k][0] ), - m_shearStiffness * ( slip[1] + m_elasticSlip[k][1] ) }; - - switch( fractureState ) - { - case fields::contact::FractureState::Stick: - { - // Elastic slip case - // Tangential components of the traction are equal to tau - tractionVector[1] = tau[0]; - tractionVector[2] = tau[1]; - - dTractionVector_dJump[1][1] = m_shearStiffness; - dTractionVector_dJump[2][2] = m_shearStiffness; - - // The slip is only elastic: we add the full slip to the elastic one - LvArray::tensorOps::add< 2 >( m_elasticSlip[k], slip ); - - break; - } - case fields::contact::FractureState::Slip: - { - // Plastic slip case - real64 dLimitTau_dNormalTraction; - real64 const limitTau = computeLimitTangentialTractionNorm( tractionVector[0], - dLimitTau_dNormalTraction ); - - real64 const slipNorm = LvArray::tensorOps::l2Norm< 2 >( slip ); - - // Tangential components of the traction computed based on the limitTau - tractionVector[1] = limitTau * slip[0] / slipNorm; - tractionVector[2] = limitTau * slip[1] / slipNorm; - - dTractionVector_dJump[1][0] = m_penaltyStiffness * dLimitTau_dNormalTraction * slip[0] / slipNorm; - dTractionVector_dJump[1][1] = limitTau * pow( slip[1], 2 ) / pow( LvArray::tensorOps::l2NormSquared< 2 >( slip ), 1.5 ); - dTractionVector_dJump[1][2] = limitTau * slip[0] * slip[1] / pow( LvArray::tensorOps::l2NormSquared< 2 >( slip ), 1.5 ); - - dTractionVector_dJump[2][0] = m_penaltyStiffness * dLimitTau_dNormalTraction * slip[1] / slipNorm; - dTractionVector_dJump[2][1] = limitTau * slip[0] * slip[1] / pow( LvArray::tensorOps::l2NormSquared< 2 >( slip ), 1.5 ); - dTractionVector_dJump[2][2] = limitTau * pow( slip[0], 2 ) / pow( LvArray::tensorOps::l2NormSquared< 2 >( slip ), 1.5 ); - - // Compute elastic component of the slip for this case - real64 const plasticSlip[2] = { tractionVector[1] / m_shearStiffness, - tractionVector[2] / m_shearStiffness }; - - LvArray::tensorOps::copy< 2 >( m_elasticSlip[k], slip ); - LvArray::tensorOps::subtract< 2 >( m_elasticSlip[k], plasticSlip ); - - break; - } - } - } -} - -GEOS_HOST_DEVICE -inline void CoulombContactUpdates::updateFractureState( localIndex const k, - arraySlice1d< real64 const > const & dispJump, - arraySlice1d< real64 const > const & tractionVector, - integer & fractureState ) const -{ - using namespace fields::contact; - - if( dispJump[0] > -m_displacementJumpThreshold ) - { - fractureState = FractureState::Open; - m_elasticSlip[k][0] = 0.0; - m_elasticSlip[k][1] = 0.0; - } - else - { - real64 const tau[2] = { tractionVector[1], - tractionVector[2] }; - real64 const tauNorm = LvArray::tensorOps::l2Norm< 2 >( tau ); - - real64 dLimitTau_dNormalTraction; - real64 const limitTau = computeLimitTangentialTractionNorm( tractionVector[0], - dLimitTau_dNormalTraction ); - - // Yield function (not necessary but makes it clearer) - real64 const yield = tauNorm - limitTau; - - fractureState = yield < 0 ? FractureState::Stick : FractureState::Slip; - } -} - -} /* namespace constitutive */ - -} /* namespace geos */ - -#endif /* GEOS_CONSTITUTIVE_CONTACT_COULOMBCONTACT_HPP_ */ diff --git a/src/coreComponents/constitutive/contact/CoulombContact.cpp b/src/coreComponents/constitutive/contact/CoulombFriction.cpp similarity index 59% rename from src/coreComponents/constitutive/contact/CoulombContact.cpp rename to src/coreComponents/constitutive/contact/CoulombFriction.cpp index 95ef7d5f8ea..bf9d89cb073 100644 --- a/src/coreComponents/constitutive/contact/CoulombContact.cpp +++ b/src/coreComponents/constitutive/contact/CoulombFriction.cpp @@ -14,10 +14,10 @@ */ /** - * @file CoulombContact.cpp + * @file CoulombFriction.cpp */ -#include "CoulombContact.hpp" +#include "CoulombFriction.hpp" namespace geos { @@ -27,12 +27,16 @@ using namespace dataRepository; namespace constitutive { -CoulombContact::CoulombContact( string const & name, Group * const parent ): - ContactBase( name, parent ), +CoulombFriction::CoulombFriction( string const & name, Group * const parent ): + FrictionBase( name, parent ), m_cohesion(), m_frictionCoefficient(), m_elasticSlip() { + registerWrapper( viewKeyStruct::shearStiffnessString(), &m_shearStiffness ). + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "Value of the shear elastic stiffness. Units of Pressure/length" ); + registerWrapper( viewKeyStruct::cohesionString(), &m_cohesion ). setApplyDefaultValue( -1 ). setInputFlag( InputFlags::REQUIRED ). @@ -49,10 +53,10 @@ CoulombContact::CoulombContact( string const & name, Group * const parent ): } -CoulombContact::~CoulombContact() +CoulombFriction::~CoulombFriction() {} -void CoulombContact::postInputInitialization() +void CoulombFriction::postInputInitialization() { GEOS_THROW_IF( m_frictionCoefficient < 0.0, getFullName() << ": The provided friction coefficient is less than zero. Value: " << m_frictionCoefficient, @@ -60,27 +64,25 @@ void CoulombContact::postInputInitialization() } -void CoulombContact::allocateConstitutiveData( Group & parent, - localIndex const numConstitutivePointsPerParentIndex ) +void CoulombFriction::allocateConstitutiveData( Group & parent, + localIndex const numConstitutivePointsPerParentIndex ) { m_elasticSlip.resize( 0, 2 ); - ContactBase::allocateConstitutiveData( parent, numConstitutivePointsPerParentIndex ); + FrictionBase::allocateConstitutiveData( parent, numConstitutivePointsPerParentIndex ); } -CoulombContactUpdates CoulombContact::createKernelWrapper() const +CoulombFrictionUpdates CoulombFriction::createKernelWrapper() const { - return CoulombContactUpdates( m_penaltyStiffness, - m_shearStiffness, - m_displacementJumpThreshold, - *m_apertureTable, - m_cohesion, - m_frictionCoefficient, - m_elasticSlip ); + return CoulombFrictionUpdates( m_displacementJumpThreshold, + m_shearStiffness, + m_cohesion, + m_frictionCoefficient, + m_elasticSlip ); } -REGISTER_CATALOG_ENTRY( ConstitutiveBase, CoulombContact, string const &, Group * const ) +REGISTER_CATALOG_ENTRY( ConstitutiveBase, CoulombFriction, string const &, Group * const ) } /* namespace constitutive */ diff --git a/src/coreComponents/constitutive/contact/CoulombFriction.hpp b/src/coreComponents/constitutive/contact/CoulombFriction.hpp new file mode 100644 index 00000000000..e2cd87f024f --- /dev/null +++ b/src/coreComponents/constitutive/contact/CoulombFriction.hpp @@ -0,0 +1,309 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 Total, S.A + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file CoulombFriction.hpp + */ + +#ifndef GEOS_CONSTITUTIVE_CONTACT_COULOMBFRICTION_HPP_ +#define GEOS_CONSTITUTIVE_CONTACT_COULOMBFRICTION_HPP_ + +#include "FrictionBase.hpp" + +namespace geos +{ + +namespace constitutive +{ + +/** + * @class CoulombFrictionUpdates + * + * This class is used for in-kernel contact relation updates + */ +class CoulombFrictionUpdates : public FrictionBaseUpdates +{ +public: + CoulombFrictionUpdates( real64 const & displacementJumpThreshold, + real64 const & shearStiffness, + real64 const & cohesion, + real64 const & frictionCoefficient, + arrayView2d< real64 > const & elasticSlip ) + : FrictionBaseUpdates( displacementJumpThreshold ), + m_shearStiffness( shearStiffness ), + m_cohesion( cohesion ), + m_frictionCoefficient( frictionCoefficient ), + m_elasticSlip( elasticSlip ) + {} + + /// Default copy constructor + CoulombFrictionUpdates( CoulombFrictionUpdates const & ) = default; + + /// Default move constructor + CoulombFrictionUpdates( CoulombFrictionUpdates && ) = default; + + /// Deleted default constructor + CoulombFrictionUpdates() = delete; + + /// Deleted copy assignment operator + CoulombFrictionUpdates & operator=( CoulombFrictionUpdates const & ) = delete; + + /// Deleted move assignment operator + CoulombFrictionUpdates & operator=( CoulombFrictionUpdates && ) = delete; + + /** + * @brief Evaluate the limit tangential traction norm and return the derivative wrt normal traction + * @param[in] normalTraction the normal traction + * @param[out] dLimitTangentialTractionNorm_dTraction the derivative of the limit tangential traction norm wrt normal traction + * @return the limit tangential traction norm + */ + GEOS_HOST_DEVICE + inline + virtual real64 computeLimitTangentialTractionNorm( real64 const & normalTraction, + real64 & dLimitTangentialTractionNorm_dTraction ) const override final; + + GEOS_HOST_DEVICE + inline + virtual void computeShearTraction( localIndex const k, + arraySlice1d< real64 const > const & oldDispJump, + arraySlice1d< real64 const > const & dispJump, + integer const & fractureState, + arraySlice1d< real64 > const & tractionVector, + arraySlice2d< real64 > const & dTractionVector_dJump ) const override final; + + GEOS_HOST_DEVICE + inline + virtual void updateFractureState( localIndex const k, + arraySlice1d< real64 const > const & dispJump, + arraySlice1d< real64 const > const & tractionVector, + integer & fractureState ) const override final; + +private: + /// The shear stiffness + real64 m_shearStiffness; + + /// The cohesion for each upper level dimension (i.e. cell) of *this + real64 m_cohesion; + + /// The friction coefficient for each upper level dimension (i.e. cell) of *this + real64 m_frictionCoefficient; + + arrayView2d< real64 > m_elasticSlip; +}; + + +/** + * @class CoulombFriction + * + * Class to provide a CoulombFriction friction model. + */ +class CoulombFriction : public FrictionBase +{ +public: + + /** + * constructor + * @param[in] name name of the instance in the catalog + * @param[in] parent the group which contains this instance + */ + CoulombFriction( string const & name, Group * const parent ); + + /** + * Default Destructor + */ + virtual ~CoulombFriction() override; + + static string catalogName() { return "Coulomb"; } + + virtual string getCatalogName() const override { return catalogName(); } + + ///@} + + virtual void allocateConstitutiveData( dataRepository::Group & parent, + localIndex const numConstitutivePointsPerParentIndex ) override final; + + /** + * @brief Const accessor for cohesion + * @return A const reference to arrayView1d containing the + * cohesions (at every element). + */ + real64 const & cohesion() const { return m_cohesion; } + + /** + * @brief Const accessor for friction angle + * @return A const reference to arrayView1d containing the + * friction coefficient (at every element). + */ + real64 const & frictionCoefficient() const { return m_frictionCoefficient; } + + /// Type of kernel wrapper for in-kernel update + using KernelWrapper = CoulombFrictionUpdates; + + /** + * @brief Create an update kernel wrapper. + * @return the wrapper + */ + KernelWrapper createKernelWrapper() const; + +protected: + + virtual void postInputInitialization() override; + +private: + + /// The shear stiffness + real64 m_shearStiffness; + + /// The cohesion for each upper level dimension (i.e. cell) of *this + real64 m_cohesion; + + /// The friction coefficient for each upper level dimension (i.e. cell) of *this + real64 m_frictionCoefficient; + + /// Elastic slip + array2d< real64 > m_elasticSlip; + +/** + * @struct Set of "char const *" and keys for data specified in this class. + */ + struct viewKeyStruct : public FrictionBase::viewKeyStruct + { + /// string/key for shear stiffness + static constexpr char const * shearStiffnessString() { return "shearStiffness"; } + + /// string/key for cohesion + static constexpr char const * cohesionString() { return "cohesion"; } + + /// string/key for friction coefficient + static constexpr char const * frictionCoefficientString() { return "frictionCoefficient"; } + + /// string/key for the elastic slip + static constexpr char const * elasticSlipString() { return "elasticSlip"; } + }; + +}; + + +GEOS_HOST_DEVICE +real64 CoulombFrictionUpdates::computeLimitTangentialTractionNorm( real64 const & normalTraction, + real64 & dLimitTangentialTractionNorm_dTraction ) const +{ + dLimitTangentialTractionNorm_dTraction = m_frictionCoefficient; + return ( m_cohesion - normalTraction * m_frictionCoefficient ); +} + + +GEOS_HOST_DEVICE +inline void CoulombFrictionUpdates::computeShearTraction( localIndex const k, + arraySlice1d< real64 const > const & oldDispJump, + arraySlice1d< real64 const > const & dispJump, + integer const & fractureState, + arraySlice1d< real64 > const & tractionVector, + arraySlice2d< real64 > const & dTractionVector_dJump ) const +{ + // Compute the slip + real64 const slip[2] = { dispJump[1] - oldDispJump[1], + dispJump[2] - oldDispJump[2] }; + + + real64 const tau[2] = { m_shearStiffness * ( slip[0] + m_elasticSlip[k][0] ), + m_shearStiffness * ( slip[1] + m_elasticSlip[k][1] ) }; + + switch( fractureState ) + { + case fields::contact::FractureState::Stick: + { + // Elastic slip case + // Tangential components of the traction are equal to tau + tractionVector[1] = tau[0]; + tractionVector[2] = tau[1]; + + dTractionVector_dJump[1][1] = m_shearStiffness; + dTractionVector_dJump[2][2] = m_shearStiffness; + + // The slip is only elastic: we add the full slip to the elastic one + LvArray::tensorOps::add< 2 >( m_elasticSlip[k], slip ); + + break; + } + case fields::contact::FractureState::Slip: + { + // Plastic slip case + real64 dLimitTau_dNormalTraction; + real64 const limitTau = computeLimitTangentialTractionNorm( tractionVector[0], + dLimitTau_dNormalTraction ); + + real64 const slipNorm = LvArray::tensorOps::l2Norm< 2 >( slip ); + + // Tangential components of the traction computed based on the limitTau + tractionVector[1] = limitTau * slip[0] / slipNorm; + tractionVector[2] = limitTau * slip[1] / slipNorm; + + dTractionVector_dJump[1][0] = dTractionVector_dJump[0][0] * dLimitTau_dNormalTraction * slip[0] / slipNorm; + dTractionVector_dJump[1][1] = limitTau * pow( slip[1], 2 ) / pow( LvArray::tensorOps::l2NormSquared< 2 >( slip ), 1.5 ); + dTractionVector_dJump[1][2] = limitTau * slip[0] * slip[1] / pow( LvArray::tensorOps::l2NormSquared< 2 >( slip ), 1.5 ); + + dTractionVector_dJump[2][0] = dTractionVector_dJump[0][0] * dLimitTau_dNormalTraction * slip[1] / slipNorm; + dTractionVector_dJump[2][1] = limitTau * slip[0] * slip[1] / pow( LvArray::tensorOps::l2NormSquared< 2 >( slip ), 1.5 ); + dTractionVector_dJump[2][2] = limitTau * pow( slip[0], 2 ) / pow( LvArray::tensorOps::l2NormSquared< 2 >( slip ), 1.5 ); + + // Compute elastic component of the slip for this case + real64 const plasticSlip[2] = { tractionVector[1] / m_shearStiffness, + tractionVector[2] / m_shearStiffness }; + + LvArray::tensorOps::copy< 2 >( m_elasticSlip[k], slip ); + LvArray::tensorOps::subtract< 2 >( m_elasticSlip[k], plasticSlip ); + + break; + } + } +} + +GEOS_HOST_DEVICE +inline void CoulombFrictionUpdates::updateFractureState( localIndex const k, + arraySlice1d< real64 const > const & dispJump, + arraySlice1d< real64 const > const & tractionVector, + integer & fractureState ) const +{ + using namespace fields::contact; + + if( dispJump[0] > -m_displacementJumpThreshold ) + { + fractureState = FractureState::Open; + m_elasticSlip[k][0] = 0.0; + m_elasticSlip[k][1] = 0.0; + } + else + { + real64 const tau[2] = { tractionVector[1], + tractionVector[2] }; + real64 const tauNorm = LvArray::tensorOps::l2Norm< 2 >( tau ); + + real64 dLimitTau_dNormalTraction; + real64 const limitTau = computeLimitTangentialTractionNorm( tractionVector[0], + dLimitTau_dNormalTraction ); + + // Yield function (not necessary but makes it clearer) + real64 const yield = tauNorm - limitTau; + + fractureState = yield < 0 ? FractureState::Stick : FractureState::Slip; + } +} + +} /* namespace constitutive */ + +} /* namespace geos */ + +#endif /* GEOS_CONSTITUTIVE_CONTACT_COULOMBFRICTION_HPP_ */ diff --git a/src/coreComponents/constitutive/contact/FrictionBase.cpp b/src/coreComponents/constitutive/contact/FrictionBase.cpp new file mode 100644 index 00000000000..a4443bcf335 --- /dev/null +++ b/src/coreComponents/constitutive/contact/FrictionBase.cpp @@ -0,0 +1,51 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2020 TotalEnergies + * Copyright (c) 2019- GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file FrictionBase.cpp + */ + +#include "FrictionBase.hpp" +#include "functions/FunctionManager.hpp" +#include "functions/TableFunction.hpp" + +namespace geos +{ + +using namespace dataRepository; + +namespace constitutive +{ + +FrictionBase::FrictionBase( string const & name, + Group * const parent ): + ConstitutiveBase( name, parent ) +{ + registerWrapper( viewKeyStruct::displacementJumpThresholdString(), &m_displacementJumpThreshold ). + setApplyDefaultValue( std::numeric_limits< real64 >::epsilon() ). + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "A threshold valued to determine whether a fracture is open or not." ); +} + +FrictionBase::~FrictionBase() +{} + +FrictionBaseUpdates FrictionBase::createKernelWrapper() const +{ + return FrictionBaseUpdates( m_displacementJumpThreshold ); +} + +} /* end namespace constitutive */ + +} /* end namespace geos */ diff --git a/src/coreComponents/constitutive/contact/FrictionBase.hpp b/src/coreComponents/constitutive/contact/FrictionBase.hpp new file mode 100644 index 00000000000..13186378e42 --- /dev/null +++ b/src/coreComponents/constitutive/contact/FrictionBase.hpp @@ -0,0 +1,167 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2020 TotalEnergies + * Copyright (c) 2019- GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file FrictionBase.hpp + */ + +#ifndef GEOS_CONSTITUTIVE_CONTACT_FRICTIONBASE_HPP_ +#define GEOS_CONSTITUTIVE_CONTACT_FRICTIONBASE_HPP_ + +#include "constitutive/ConstitutiveBase.hpp" +#include "functions/TableFunction.hpp" +#include "physicsSolvers/contact/ContactFields.hpp" + + +namespace geos +{ + +namespace constitutive +{ + +/** + * @class FrictionBaseUpdates + * + * This class is used for in-kernel contact relation updates + */ +class FrictionBaseUpdates +{ +public: + + FrictionBaseUpdates( real64 const & displacementJumpThreshold ) + : m_displacementJumpThreshold( displacementJumpThreshold ) + {} + + /// Default copy constructor + FrictionBaseUpdates( FrictionBaseUpdates const & ) = default; + + /// Default move constructor + FrictionBaseUpdates( FrictionBaseUpdates && ) = default; + + /// Deleted default constructor + FrictionBaseUpdates() = default; + + /// Deleted copy assignment operator + FrictionBaseUpdates & operator=( FrictionBaseUpdates const & ) = delete; + + /// Deleted move assignment operator + FrictionBaseUpdates & operator=( FrictionBaseUpdates && ) = delete; + + /** + * @brief Evaluate the traction vector and its derivatives wrt to pressure and jump + * @param[in] dispJump the displacement jump + * @param[in] fractureState the fracture state + * @param[out] tractionVector the traction vector + * @param[out] dTractionVector_dJump the derivative of the traction vector wrt displacement jump + */ + GEOS_HOST_DEVICE + inline + virtual void computeShearTraction( localIndex const k, + arraySlice1d< real64 const > const & oldDispJump, + arraySlice1d< real64 const > const & dispJump, + integer const & fractureState, + arraySlice1d< real64 > const & tractionVector, + arraySlice2d< real64 > const & dTractionVector_dJump ) const + {GEOS_UNUSED_VAR( k, oldDispJump, dispJump, tractionVector, dTractionVector_dJump, fractureState );} + + /** + * @brief Evaluate the traction vector and its derivatives wrt to pressure and jump + * @param[in] dispJump the displacement jump + * @param[in] tractionVector the traction vector + * @param[out] fractureState the fracture state + */ + GEOS_HOST_DEVICE + inline + virtual void updateFractureState( localIndex const k, + arraySlice1d< real64 const > const & dispJump, + arraySlice1d< real64 const > const & tractionVector, + integer & fractureState ) const + { GEOS_UNUSED_VAR( k, dispJump, tractionVector, fractureState ); } + + + + /** + * @brief Evaluate the limit tangential traction norm and return the derivative wrt normal traction + * @param[in] normalTraction the normal traction + * @param[out] dLimitTangentialTractionNorm_dTraction the derivative of the limit tangential traction norm wrt normal traction + * @return the limit tangential traction norm + */ + GEOS_HOST_DEVICE + inline + virtual real64 computeLimitTangentialTractionNorm( real64 const & normalTraction, + real64 & dLimitTangentialTractionNorm_dTraction ) const + { GEOS_UNUSED_VAR( normalTraction, dLimitTangentialTractionNorm_dTraction ); return 0; }; + +protected: + + /// A threshold valued to determine whether a fracture is open or not. + real64 m_displacementJumpThreshold; + +}; + + +/** + * @class FrictionBase + * + * This class serves as the interface for implementing contact enforcement constitutive relations. + * This does not include the actual enforcement algorithm, but only the constitutive relations that + * govern the behavior of the contact. So things like penalty, or friction, or kinematic constraint. + */ +class FrictionBase : public ConstitutiveBase +{ +public: + + /** + * @brief The standard data repository constructor + * @param name The name of the relation in the data repository + * @param parent The name of the parent Group that holds this relation object. + */ + FrictionBase( string const & name, + Group * const parent ); + + /** + * @brief default destructor + */ + virtual ~FrictionBase() override; + + /// Type of kernel wrapper for in-kernel update + using KernelWrapper = FrictionBaseUpdates; + + /** + * @brief Create an update kernel wrapper. + * @return the wrapper + */ + KernelWrapper createKernelWrapper() const; + + /** + * @struct Structure to hold scoped key names + */ + struct viewKeyStruct : public ConstitutiveBase::viewKeyStruct + { + /// string/key for the displacement jump threshold value + static constexpr char const * displacementJumpThresholdString() { return "displacementJumpThreshold"; } + }; + +protected: + + /// A threshold valued to determine whether a fracture is open or not. + real64 m_displacementJumpThreshold; + +}; + +} /* namespace constitutive */ + +} /* namespace geos */ + +#endif /* GEOS_CONSTITUTIVE_CONTACT_FRICTIONBASE_HPP_ */ diff --git a/src/coreComponents/constitutive/contact/ContactSelector.hpp b/src/coreComponents/constitutive/contact/FrictionSelector.hpp similarity index 68% rename from src/coreComponents/constitutive/contact/ContactSelector.hpp rename to src/coreComponents/constitutive/contact/FrictionSelector.hpp index e067ddfc90b..70f510eb0a7 100644 --- a/src/coreComponents/constitutive/contact/ContactSelector.hpp +++ b/src/coreComponents/constitutive/contact/FrictionSelector.hpp @@ -6,7 +6,7 @@ * Copyright (c) 2018-2024 Total, S.A * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University * Copyright (c) 2018-2024 Chevron - * Copyright (c) 2019- GEOS/GEOSX Contributors + * Copyright (c) 2019- GEOS/GEOSX Contributors * All rights reserved * * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. @@ -14,14 +14,14 @@ */ /** - * @file ContactSelector.hpp + * @file FrictionSelector.hpp */ #ifndef GEOS_CONSTITUTIVE_CONTACT_CONTACTSELECTOR_HPP_ #define GEOS_CONSTITUTIVE_CONTACT_CONTACTSELECTOR_HPP_ #include "constitutive/ConstitutivePassThruHandler.hpp" -#include "constitutive/contact/CoulombContact.hpp" +#include "constitutive/contact/CoulombFriction.hpp" #include "constitutive/contact/FrictionlessContact.hpp" namespace geos @@ -31,19 +31,19 @@ namespace constitutive { template< typename LAMBDA > -void constitutiveUpdatePassThru( ContactBase const & contact, +void constitutiveUpdatePassThru( FrictionBase const & contact, LAMBDA && lambda ) { - ConstitutivePassThruHandler< CoulombContact, - FrictionlessContact >::execute( contact, std::forward< LAMBDA >( lambda ) ); + ConstitutivePassThruHandler< FrictionlessContact, + CoulombFriction >::execute( contact, std::forward< LAMBDA >( lambda ) ); } template< typename LAMBDA > -void constitutiveUpdatePassThru( ContactBase & contact, +void constitutiveUpdatePassThru( FrictionBase & contact, LAMBDA && lambda ) { - ConstitutivePassThruHandler< CoulombContact, - FrictionlessContact >::execute( contact, std::forward< LAMBDA >( lambda ) ); + ConstitutivePassThruHandler< FrictionlessContact, + CoulombFriction >::execute( contact, std::forward< LAMBDA >( lambda ) ); } } /* namespace constitutive */ diff --git a/src/coreComponents/constitutive/contact/FrictionlessContact.cpp b/src/coreComponents/constitutive/contact/FrictionlessContact.cpp index 5a30413f633..a74b36a76fb 100644 --- a/src/coreComponents/constitutive/contact/FrictionlessContact.cpp +++ b/src/coreComponents/constitutive/contact/FrictionlessContact.cpp @@ -29,24 +29,15 @@ namespace constitutive FrictionlessContact::FrictionlessContact( string const & name, Group * const parent ): - ContactBase( name, parent ) + FrictionBase( name, parent ) {} FrictionlessContact::~FrictionlessContact() {} -void FrictionlessContact::allocateConstitutiveData( Group & parent, - localIndex const numConstitutivePointsPerParentIndex ) -{ - ContactBase::allocateConstitutiveData( parent, numConstitutivePointsPerParentIndex ); -} - FrictionlessContactUpdates FrictionlessContact::createKernelWrapper() const { - return FrictionlessContactUpdates( m_penaltyStiffness, - m_shearStiffness, - m_displacementJumpThreshold, - *m_apertureTable ); + return FrictionlessContactUpdates( m_displacementJumpThreshold ); } REGISTER_CATALOG_ENTRY( ConstitutiveBase, FrictionlessContact, string const &, Group * const ) diff --git a/src/coreComponents/constitutive/contact/FrictionlessContact.hpp b/src/coreComponents/constitutive/contact/FrictionlessContact.hpp index 2b99908609e..81aeea036d3 100644 --- a/src/coreComponents/constitutive/contact/FrictionlessContact.hpp +++ b/src/coreComponents/constitutive/contact/FrictionlessContact.hpp @@ -20,7 +20,7 @@ #ifndef GEOS_CONSTITUTIVE_CONTACT_FRICTIONLESSCONTACT_HPP_ #define GEOS_CONSTITUTIVE_CONTACT_FRICTIONLESSCONTACT_HPP_ -#include "constitutive/contact/ContactBase.hpp" +#include "constitutive/contact/FrictionBase.hpp" namespace geos { @@ -33,15 +33,12 @@ namespace constitutive * * This class is used for in-kernel contact relation updates */ -class FrictionlessContactUpdates : public ContactBaseUpdates +class FrictionlessContactUpdates : public FrictionBaseUpdates { public: - FrictionlessContactUpdates( real64 const & penaltyStiffness, - real64 const & shearStiffness, - real64 const & displacementJumpThreshold, - TableFunction const & apertureTable ) - : ContactBaseUpdates( penaltyStiffness, shearStiffness, displacementJumpThreshold, apertureTable ) + FrictionlessContactUpdates( real64 const & displacementJumpThreshold ) + : FrictionBaseUpdates( displacementJumpThreshold ) {} /// Default copy constructor @@ -59,16 +56,6 @@ class FrictionlessContactUpdates : public ContactBaseUpdates /// Deleted move assignment operator FrictionlessContactUpdates & operator=( FrictionlessContactUpdates && ) = delete; - GEOS_HOST_DEVICE - inline - virtual void computeTraction( localIndex const k, - arraySlice1d< real64 const > const & oldDispJump, - arraySlice1d< real64 const > const & dispJump, - integer const & fractureState, - arraySlice1d< real64 > const & tractionVector, - arraySlice2d< real64 > const & dTractionVector_dJump ) const override final; - - GEOS_HOST_DEVICE inline virtual void updateFractureState( localIndex const k, @@ -100,7 +87,7 @@ class FrictionlessContactUpdates : public ContactBaseUpdates * This does not include the actual enforcement algorithm, but only the constitutive relations that * govern the behavior of the contact. So things like penalty, or friction, or kinematic constraint. */ -class FrictionlessContact : public ContactBase +class FrictionlessContact : public FrictionBase { public: @@ -124,9 +111,6 @@ class FrictionlessContact : public ContactBase virtual string getCatalogName() const override { return catalogName(); } - virtual void allocateConstitutiveData( dataRepository::Group & parent, - localIndex const numConstitutivePointsPerParentIndex ) override; - /// Type of kernel wrapper for in-kernel update using KernelWrapper = FrictionlessContactUpdates; @@ -146,25 +130,6 @@ class FrictionlessContact : public ContactBase }; -GEOS_HOST_DEVICE -inline void FrictionlessContactUpdates::computeTraction( localIndex const k, - arraySlice1d< real64 const > const & oldDispJump, - arraySlice1d< real64 const > const & dispJump, - integer const & fractureState, - arraySlice1d< real64 > const & tractionVector, - arraySlice2d< real64 > const & dTractionVector_dJump ) const -{ - GEOS_UNUSED_VAR( k, oldDispJump ); - - bool const isOpen = fractureState == fields::contact::FractureState::Open; - - tractionVector[0] = isOpen ? 0.0 : m_penaltyStiffness * dispJump[0]; - tractionVector[1] = 0.0; - tractionVector[2] = 0.0; - - LvArray::forValuesInSlice( dTractionVector_dJump, []( real64 & val ){ val = 0.0; } ); - dTractionVector_dJump( 0, 0 ) = isOpen ? 0.0 : m_penaltyStiffness; -} GEOS_HOST_DEVICE inline void FrictionlessContactUpdates::updateFractureState( localIndex const k, diff --git a/src/coreComponents/constitutive/contact/HydraulicApertureBase.cpp b/src/coreComponents/constitutive/contact/HydraulicApertureBase.cpp new file mode 100644 index 00000000000..81ca7a55097 --- /dev/null +++ b/src/coreComponents/constitutive/contact/HydraulicApertureBase.cpp @@ -0,0 +1,48 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2020 TotalEnergies + * Copyright (c) 2019- GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file HydraulicApertureBase.cpp + */ + + +#include "HydraulicApertureBase.hpp" + + +namespace geos +{ + +namespace constitutive +{ + +HydraulicApertureBase::HydraulicApertureBase( string const & name, + Group * const parent ): + ConstitutiveBase( name, parent ), + m_aperture0( 0.0 ) +{ + /// TODO: must become a required parameter. + registerWrapper( viewKeyStruct::apertureZeroString(), &m_aperture0 ). + setInputFlag( dataRepository::InputFlags::OPTIONAL ). + setApplyDefaultValue( 1e-6 ). + setDescription( "Reference hydraulic aperture. It is the aperture at zero normal stress." ); +} + +HydraulicApertureBase::~HydraulicApertureBase() +{} + + + +} /* namespace constitutive */ + +} /* namespace geos */ diff --git a/src/coreComponents/constitutive/contact/HydraulicApertureBase.hpp b/src/coreComponents/constitutive/contact/HydraulicApertureBase.hpp new file mode 100644 index 00000000000..dab7b35957d --- /dev/null +++ b/src/coreComponents/constitutive/contact/HydraulicApertureBase.hpp @@ -0,0 +1,73 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2020 TotalEnergies + * Copyright (c) 2019- GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file HydraulicApertureBase.hpp + */ + +#ifndef GEOS_CONSTITUTIVE_CONTACT_HYDRAULICAPERTUREBASE_HPP_ +#define GEOS_CONSTITUTIVE_CONTACT_HYDRAULICAPERTUREBASE_HPP_ + +#include "constitutive/ConstitutiveBase.hpp" +#include "functions/TableFunction.hpp" +#include "physicsSolvers/contact/ContactFields.hpp" + + +namespace geos +{ + +namespace constitutive +{ + +/** + * @class HydraulicApertureBase + * + * This class serves as the interface for implementing contact enforcement constitutive relations. + * This does not include the actual enforcement algorithm, but only the constitutive relations that + * govern the behavior of the contact. So things like penalty, or friction, or kinematic constraint. + */ +class HydraulicApertureBase : public ConstitutiveBase +{ +public: + + /** + * @brief The standard data repository constructor + * @param name The name of the relation in the data repository + * @param parent The name of the parent Group that holds this relation object. + */ + HydraulicApertureBase( string const & name, + Group * const parent ); + + /** + * @brief default destructor + */ + virtual ~HydraulicApertureBase() override; + +protected: + + struct viewKeyStruct : public ConstitutiveBase::viewKeyStruct + { + /// string/key for aperture under zero normal stress + static constexpr char const * apertureZeroString() { return "referenceAperture"; } + }; + + /// Reference hydraulic aperture. Aperture at zero normal stress + real64 m_aperture0; /// TODO: this will replace what is currently called defaultAperture. +}; + +} /* namespace constitutive */ + +} /* namespace geos */ + +#endif /* GEOS_CONSTITUTIVE_CONTACT_HYDRAULICAPERTURETABLE_HPP_ */ diff --git a/src/coreComponents/constitutive/contact/HydraulicApertureRelationSelector.hpp b/src/coreComponents/constitutive/contact/HydraulicApertureRelationSelector.hpp new file mode 100644 index 00000000000..0d32794c1eb --- /dev/null +++ b/src/coreComponents/constitutive/contact/HydraulicApertureRelationSelector.hpp @@ -0,0 +1,52 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2020 TotalEnergies + * Copyright (c) 2019- GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file FrictionSelector.hpp + */ + +#ifndef GEOS_CONSTITUTIVE_CONTACT_CONTACTSELECTOR_HPP_ +#define GEOS_CONSTITUTIVE_CONTACT_CONTACTSELECTOR_HPP_ + +#include "constitutive/ConstitutivePassThruHandler.hpp" +#include "constitutive/contact/HydraulicApertureTable.hpp" +#include "constitutive/contact/BartonBandis.hpp" + +namespace geos +{ + +namespace constitutive +{ + +template< typename LAMBDA > +void constitutiveUpdatePassThru( HydraulicApertureBase const & contact, + LAMBDA && lambda ) +{ + ConstitutivePassThruHandler< HydraulicApertureTable, + BartonBandis >::execute( contact, std::forward< LAMBDA >( lambda ) ); +} + +template< typename LAMBDA > +void constitutiveUpdatePassThru( HydraulicApertureBase & contact, + LAMBDA && lambda ) +{ + ConstitutivePassThruHandler< HydraulicApertureTable, + BartonBandis >::execute( contact, std::forward< LAMBDA >( lambda ) ); +} + +} /* namespace constitutive */ + +} /* namespace geos */ + +#endif // GEOS_CONSTITUTIVE_CONTACT_CONTACTSELECTOR_HPP_ diff --git a/src/coreComponents/constitutive/contact/ContactBase.cpp b/src/coreComponents/constitutive/contact/HydraulicApertureTable.cpp similarity index 78% rename from src/coreComponents/constitutive/contact/ContactBase.cpp rename to src/coreComponents/constitutive/contact/HydraulicApertureTable.cpp index fa4de0fb776..60c3f0ab807 100644 --- a/src/coreComponents/constitutive/contact/ContactBase.cpp +++ b/src/coreComponents/constitutive/contact/HydraulicApertureTable.cpp @@ -14,10 +14,10 @@ */ /** - * @file ContactBase.cpp + * @file HydraulicApertureTable.cpp */ -#include "ContactBase.hpp" +#include "HydraulicApertureTable.hpp" #include "functions/FunctionManager.hpp" #include "functions/TableFunction.hpp" @@ -29,19 +29,11 @@ using namespace dataRepository; namespace constitutive { -ContactBase::ContactBase( string const & name, - Group * const parent ): - ConstitutiveBase( name, parent ), +HydraulicApertureTable::HydraulicApertureTable( string const & name, + Group * const parent ): + HydraulicApertureBase( name, parent ), m_apertureTable( nullptr ) { - registerWrapper( viewKeyStruct::penaltyStiffnessString(), &m_penaltyStiffness ). - setInputFlag( InputFlags::OPTIONAL ). - setDescription( "Value of the penetration penalty stiffness. Units of Pressure/length" ); - - registerWrapper( viewKeyStruct::shearStiffnessString(), &m_shearStiffness ). - setInputFlag( InputFlags::OPTIONAL ). - setDescription( "Value of the shear elastic stiffness. Units of Pressure/length" ); - registerWrapper( viewKeyStruct::apertureToleranceString(), &m_apertureTolerance ). setApplyDefaultValue( 1.0e-9 ). setInputFlag( InputFlags::OPTIONAL ). @@ -52,23 +44,18 @@ ContactBase::ContactBase( string const & name, "to smooth out highly nonlinear behavior associated with 1/0 in addition to avoiding the " "1/0 error." ); - registerWrapper( viewKeyStruct::displacementJumpThresholdString(), &m_displacementJumpThreshold ). - setApplyDefaultValue( std::numeric_limits< real64 >::epsilon() ). - setInputFlag( InputFlags::OPTIONAL ). - setDescription( "A threshold valued to determine whether a fracture is open or not." ); - - registerWrapper( viewKeyStruct::apertureTableNameString(), &m_apertureTableName ). setRTTypeName( rtTypes::CustomTypes::groupNameRef ). setInputFlag( InputFlags::REQUIRED ). setDescription( "Name of the aperture table" ); } -ContactBase::~ContactBase() +HydraulicApertureTable::~HydraulicApertureTable() {} -void ContactBase::postInputInitialization() + +void HydraulicApertureTable::postInputInitialization() { GEOS_THROW_IF( m_apertureTableName.empty(), @@ -76,12 +63,9 @@ void ContactBase::postInputInitialization() } -void ContactBase::initializePreSubGroups() -{} - -void ContactBase::allocateConstitutiveData( Group & parent, - localIndex const numConstitutivePointsPerParentIndex ) +void HydraulicApertureTable::allocateConstitutiveData( Group & parent, + localIndex const numConstitutivePointsPerParentIndex ) { ConstitutiveBase::allocateConstitutiveData( parent, numConstitutivePointsPerParentIndex ); @@ -133,7 +117,7 @@ void ContactBase::allocateConstitutiveData( Group & parent, } -void ContactBase::validateApertureTable( TableFunction const & apertureTable ) const +void HydraulicApertureTable::validateApertureTable( TableFunction const & apertureTable ) const { ArrayOfArraysView< real64 const > const coords = apertureTable.getCoordinates(); arrayView1d< real64 const > const & hydraulicApertureValues = apertureTable.getValues(); @@ -163,14 +147,13 @@ void ContactBase::validateApertureTable( TableFunction const & apertureTable ) c -ContactBaseUpdates ContactBase::createKernelWrapper() const +HydraulicApertureTableUpdates HydraulicApertureTable::createKernelWrapper() const { - return ContactBaseUpdates( m_penaltyStiffness, - m_shearStiffness, - m_displacementJumpThreshold, - *m_apertureTable ); + return HydraulicApertureTableUpdates( *m_apertureTable ); } +REGISTER_CATALOG_ENTRY( ConstitutiveBase, HydraulicApertureTable, string const &, Group * const ) + } /* end namespace constitutive */ } /* end namespace geos */ diff --git a/src/coreComponents/constitutive/contact/HydraulicApertureTable.hpp b/src/coreComponents/constitutive/contact/HydraulicApertureTable.hpp new file mode 100644 index 00000000000..c1e9c47c1ce --- /dev/null +++ b/src/coreComponents/constitutive/contact/HydraulicApertureTable.hpp @@ -0,0 +1,170 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 Total, S.A + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file HydraulicApertureTable.hpp + */ + +#ifndef GEOS_CONSTITUTIVE_CONTACT_HYDRAULICAPERTURETABLE_HPP_ +#define GEOS_CONSTITUTIVE_CONTACT_HYDRAULICAPERTURETABLE_HPP_ + +#include "constitutive/contact/HydraulicApertureBase.hpp" +#include "functions/TableFunction.hpp" + + +namespace geos +{ + +namespace constitutive +{ + +/** + * @class HydraulicApertureTableUpdates + * + * This class is used for in-kernel contact relation updates + */ +class HydraulicApertureTableUpdates +{ +public: + + HydraulicApertureTableUpdates( TableFunction const & apertureTable ) + : m_apertureTable( apertureTable.createKernelWrapper() ) + {} + + /// Default copy constructor + HydraulicApertureTableUpdates( HydraulicApertureTableUpdates const & ) = default; + + /// Default move constructor + HydraulicApertureTableUpdates( HydraulicApertureTableUpdates && ) = default; + + /// Deleted default constructor + HydraulicApertureTableUpdates() = default; + + /// Deleted copy assignment operator + HydraulicApertureTableUpdates & operator=( HydraulicApertureTableUpdates const & ) = delete; + + /// Deleted move assignment operator + HydraulicApertureTableUpdates & operator=( HydraulicApertureTableUpdates && ) = delete; + + /** + * @brief Evaluate the effective aperture, and its derivative wrt aperture + * @param[in] aperture the model aperture/gap + * @param[out] dHydraulicAperture_dAperture the derivative of the effective aperture wrt aperture + * @return The hydraulic aperture that is always > 0 + */ + GEOS_HOST_DEVICE + real64 computeHydraulicAperture( real64 const aperture, + real64 const normalTraction, + real64 & dHydraulicAperture_daperture, + real64 & dHydraulicAperture_dNormalStress ) const; + +protected: + + /// The aperture table function wrapper + TableFunction::KernelWrapper m_apertureTable; +}; + + +/** + * @class HydraulicApertureTable + * + * This class serves as the interface for implementing contact enforcement constitutive relations. + * This does not include the actual enforcement algorithm, but only the constitutive relations that + * govern the behavior of the contact. So things like penalty, or friction, or kinematic constraint. + */ +class HydraulicApertureTable : public HydraulicApertureBase +{ +public: + + /** + * @brief The standard data repository constructor + * @param name The name of the relation in the data repository + * @param parent The name of the parent Group that holds this relation object. + */ + HydraulicApertureTable( string const & name, + Group * const parent ); + + /** + * @brief default destructor + */ + virtual ~HydraulicApertureTable() override; + + static string catalogName() { return "HydraulicApertureTable"; } + + virtual string getCatalogName() const override { return catalogName(); } + + virtual void allocateConstitutiveData( dataRepository::Group & parent, + localIndex const numConstitutivePointsPerParentIndex ) override final; + + + + /// Type of kernel wrapper for in-kernel update + using KernelWrapper = HydraulicApertureTableUpdates; + + /** + * @brief Create an update kernel wrapper. + * @return the wrapper + */ + KernelWrapper createKernelWrapper() const; + + /** + * @struct Structure to hold scoped key names + */ + struct viewKeyStruct : public ConstitutiveBase::viewKeyStruct + { + /// string/key for aperture tolerance + static constexpr char const * apertureToleranceString() { return "apertureTolerance"; } + + /// string/key for aperture table name + static constexpr char const * apertureTableNameString() { return "apertureTableName"; } + }; + + +protected: + + virtual void postInputInitialization() override; + + /** + * @brief Validate the values provided in the aperture table + * @param[in] apertureTable the effective aperture vs aperture table + */ + void validateApertureTable( TableFunction const & apertureTable ) const; + + /// The aperture tolerance to avoid floating point errors in expressions involving aperture + real64 m_apertureTolerance; + + /// The name of the aperture table, if any + string m_apertureTableName; + + /// Pointer to the function that limits the model aperture to a physically admissible value. + TableFunction const * m_apertureTable; +}; + +GEOS_HOST_DEVICE +GEOS_FORCE_INLINE +real64 HydraulicApertureTableUpdates::computeHydraulicAperture( real64 const aperture, + real64 const normalTraction, + real64 & dHydraulicAperture_dAperture, + real64 & dHydraulicAperture_dNormalStress ) const +{ + GEOS_UNUSED_VAR( normalTraction, dHydraulicAperture_dNormalStress ); + return m_apertureTable.compute( &aperture, &dHydraulicAperture_dAperture ); +} + +} /* namespace constitutive */ + +} /* namespace geos */ + +#endif /* GEOS_CONSTITUTIVE_CONTACT_HYDRAULICAPERTURETABLE_HPP_ */ diff --git a/src/coreComponents/constitutive/contact/deprecated/ApertureTableContact.cpp b/src/coreComponents/constitutive/contact/deprecated/ApertureTableContact.cpp deleted file mode 100644 index 3f6c6ba15a5..00000000000 --- a/src/coreComponents/constitutive/contact/deprecated/ApertureTableContact.cpp +++ /dev/null @@ -1,132 +0,0 @@ -/* - * ------------------------------------------------------------------------------------------------------------ - * SPDX-License-Identifier: LGPL-2.1-only - * - * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC - * Copyright (c) 2018-2024 Total, S.A - * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University - * Copyright (c) 2018-2024 Chevron - * Copyright (c) 2019- GEOS/GEOSX Contributors - * All rights reserved - * - * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. - * ------------------------------------------------------------------------------------------------------------ - */ - -/** - * @file ApertureTableContact.cpp - */ - -#include "../deprecated/ApertureTableContact.hpp" - -#include "functions/FunctionManager.hpp" -#include "functions/TableFunction.hpp" - -namespace geos -{ - -using namespace dataRepository; - -namespace constitutive -{ - -ApertureTableContact::ApertureTableContact( string const & name, - Group * const parent ) - : ContactBase( name, parent ), - m_apertureTolerance( 1.0e-99 ), - m_apertureTable( nullptr ) -{ - registerWrapper( viewKeyStruct::apertureToleranceString(), &m_apertureTolerance ). - setApplyDefaultValue( 1.0e-9 ). - setInputFlag( InputFlags::OPTIONAL ). - setDescription( "Value to be used to avoid floating point errors in expressions involving aperture. " - "For example in the case of dividing by the actual aperture (not the effective aperture " - "that results from the aperture function) this value may be used to avoid the 1/0 error. " - "Note that this value may have some physical significance in its usage, as it may be used " - "to smooth out highly nonlinear behavior associated with 1/0 in addition to avoiding the " - "1/0 error." ); - - registerWrapper( viewKeyStruct::apertureTableNameString(), &m_apertureTableName ). - setRTTypeName( rtTypes::CustomTypes::groupNameRef ). - setInputFlag( InputFlags::REQUIRED ). - setDescription( "Name of the aperture table" ); -} - -ApertureTableContact::~ApertureTableContact() -{} - -void ApertureTableContact::postInputInitialization() -{ - FunctionManager const & functionManager = FunctionManager::getInstance(); - - GEOS_THROW_IF( m_apertureTableName.empty(), - getFullName() << ": the aperture table name " << m_apertureTableName << " is empty", - InputError ); - - GEOS_THROW_IF( !functionManager.hasGroup( m_apertureTableName ), - getFullName() << ": the aperture table named " << m_apertureTableName << " could not be found", - InputError ); -} - -void ApertureTableContact::initializePreSubGroups() -{ - FunctionManager & functionManager = FunctionManager::getInstance(); - TableFunction & apertureTable = functionManager.getGroup< TableFunction >( m_apertureTableName ); - validateApertureTable( apertureTable ); - - ArrayOfArraysView< real64 > coords = apertureTable.getCoordinates(); - arraySlice1d< real64 const > apertureValues = coords[0]; - array1d< real64 > & effectiveApertureValues = apertureTable.getValues(); - - localIndex const n = apertureValues.size()-1; - real64 const slope = ( effectiveApertureValues[n] - effectiveApertureValues[n-1] ) / ( apertureValues[n] - apertureValues[n-1] ); - real64 const apertureTransition = ( effectiveApertureValues[n] - slope * apertureValues[n] ) / ( 1.0 - slope ); - - coords.emplaceBack( 0, apertureTransition ); - effectiveApertureValues.emplace_back( apertureTransition ); - coords.emplaceBack( 0, apertureTransition*10e9 ); - effectiveApertureValues.emplace_back( apertureTransition*10e9 ); - apertureTable.reInitializeFunction(); - - m_apertureTable = &apertureTable; -} - -void ApertureTableContact::validateApertureTable( TableFunction const & apertureTable ) const -{ - ArrayOfArraysView< real64 const > const coords = apertureTable.getCoordinates(); - arrayView1d< real64 const > const & effectiveApertureValues = apertureTable.getValues(); - - GEOS_THROW_IF( coords.size() > 1, - getFullName() << ": Aperture limiter table cannot be greater than a 1D table.", - InputError ); - - arraySlice1d< real64 const > apertureValues = coords[0]; - localIndex const size = apertureValues.size(); - - GEOS_THROW_IF( coords( 0, size-1 ) > 0.0 || coords( 0, size-1 ) < 0.0, - getFullName() << ": Invalid aperture limiter table. Last coordinate must be zero!", - InputError ); - - GEOS_THROW_IF( apertureValues.size() < 2, - getFullName() << ": Invalid aperture limiter table. Must have more than two points specified", - InputError ); - - localIndex const n = apertureValues.size()-1; - real64 const slope = ( effectiveApertureValues[n] - effectiveApertureValues[n-1] ) / ( apertureValues[n] - apertureValues[n-1] ); - - GEOS_THROW_IF( slope >= 1.0, - getFullName() << ": Invalid aperture table. The slope of the last two points >= 1 is invalid.", - InputError ); -} - -ApertureTableContactUpdates ApertureTableContact::createKernelWrapper() const -{ - return ApertureTableContactUpdates( m_penaltyStiffness, - *m_apertureTable ); -} - -REGISTER_CATALOG_ENTRY( ConstitutiveBase, ApertureTableContact, string const &, Group * const ) - -} /* namespace constitutive */ - -} /* namespace geos */ diff --git a/src/coreComponents/constitutive/contact/deprecated/ApertureTableContact.hpp b/src/coreComponents/constitutive/contact/deprecated/ApertureTableContact.hpp deleted file mode 100644 index 5f82cdf00d2..00000000000 --- a/src/coreComponents/constitutive/contact/deprecated/ApertureTableContact.hpp +++ /dev/null @@ -1,178 +0,0 @@ -/* - * ------------------------------------------------------------------------------------------------------------ - * SPDX-License-Identifier: LGPL-2.1-only - * - * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC - * Copyright (c) 2018-2024 Total, S.A - * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University - * Copyright (c) 2018-2024 Chevron - * Copyright (c) 2019- GEOS/GEOSX Contributors - * All rights reserved - * - * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. - * ------------------------------------------------------------------------------------------------------------ - */ - -/** - * @file ApertureTableContact.hpp - */ - -#ifndef GEOS_CONSTITUTIVE_CONTACT_APERTURETABLECONTACT_HPP_ -#define GEOS_CONSTITUTIVE_CONTACT_APERTURETABLECONTACT_HPP_ - -#include "constitutive/contact/ContactBase.hpp" -#include "functions/TableFunction.hpp" - -namespace geos -{ - -namespace constitutive -{ - -/** - * @class ApertureTableContactUpdates - * - * This class is used for in-kernel contact relation updates - */ -class ApertureTableContactUpdates : public ContactBaseUpdates -{ -public: - - ApertureTableContactUpdates( real64 const & penaltyStiffness, - TableFunction const & apertureTable ) - : ContactBaseUpdates( penaltyStiffness ), - m_apertureTable( apertureTable.createKernelWrapper() ) - {} - - /// Default copy constructor - ApertureTableContactUpdates( ApertureTableContactUpdates const & ) = default; - - /// Default move constructor - ApertureTableContactUpdates( ApertureTableContactUpdates && ) = default; - - /// Deleted default constructor - ApertureTableContactUpdates() = delete; - - /// Deleted copy assignment operator - ApertureTableContactUpdates & operator=( ApertureTableContactUpdates const & ) = delete; - - /// Deleted move assignment operator - ApertureTableContactUpdates & operator=( ApertureTableContactUpdates && ) = delete; - - GEOS_HOST_DEVICE - inline - virtual real64 computeEffectiveAperture( real64 const aperture, - real64 & dEffectiveAperture_dAperture ) const; - - GEOS_HOST_DEVICE - inline - virtual void computeTraction( arraySlice1d< real64 const > const & dispJump, - arraySlice1d< real64 > const & tractionVector, - arraySlice2d< real64 > const & dTractionVector_dJump ) const; - - GEOS_HOST_DEVICE - inline - void addPressureToTraction( real64 const & pressure, - bool const isOpen, - arraySlice1d< real64 >const & tractionVector, - real64 & dTraction_dPressure ) const; - -private: - - /// The aperture table function wrapper - TableFunction::KernelWrapper m_apertureTable; - -}; - - -/** - * @class ApertureTableContact - */ -class ApertureTableContact : public ContactBase -{ -public: - - /** - * @brief The standard data repository constructor - * @param name The name of the relation in the data repository - * @param parent The name of the parent Group that holds this relation object. - */ - ApertureTableContact( string const & name, - Group * const parent ); - - /** - * @brief default destructor - */ - virtual ~ApertureTableContact() override; - - /** - * @brief Name that is used to register this a type of "ApertureTableContact" in the object catalog - * @return See description - */ - static string catalogName() { return "Contact"; } - - virtual string getCatalogName() const override { return catalogName(); } - - /** - * @brief accessor for aperture tolerance - * @return the aperture tolerance - */ - real64 apertureTolerance() const { return m_apertureTolerance; } - - /// Type of kernel wrapper for in-kernel update - using KernelWrapper = ApertureTableContactUpdates; - - /** - * @brief Create an update kernel wrapper. - * @return the wrapper - */ - KernelWrapper createKernelWrapper() const; - - /** - * @struct Structure to hold scoped key names - */ - struct viewKeyStruct : public ConstitutiveBase::viewKeyStruct - { - /// string/key for aperture tolerance - static constexpr char const * apertureToleranceString() { return "apertureTolerance"; } - - /// string/key for aperture table name - static constexpr char const * apertureTableNameString() { return "apertureTableName"; } - }; - -protected: - - virtual void postInputInitialization() override; - - virtual void initializePreSubGroups() override; - - /** - * @brief Validate the values provided in the aperture table - * @param[in] apertureTable the effective aperture vs aperture table - */ - void validateApertureTable( TableFunction const & apertureTable ) const; - - /// The aperture tolerance to avoid floating point errors in expressions involving aperture - real64 m_apertureTolerance; - - /// The name of the aperture table, if any - string m_apertureTableName; - - /// Pointer to the function that limits the model aperture to a physically admissible value. - TableFunction const * m_apertureTable; - -}; - -GEOS_HOST_DEVICE -real64 ApertureTableContactUpdates::computeEffectiveAperture( real64 const aperture, - real64 & dEffectiveAperture_dAperture ) const -{ - return m_apertureTable.compute( &aperture, &dEffectiveAperture_dAperture ); -} - - -} /* namespace constitutive */ - -} /* namespace geos */ - -#endif /* GEOS_CONSTITUTIVE_CONTACT_APERTURETABLECONTACT_HPP_ */ diff --git a/src/coreComponents/fileIO/silo/SiloFile.cpp b/src/coreComponents/fileIO/silo/SiloFile.cpp index 4abcb5894d4..68f3d7289eb 100644 --- a/src/coreComponents/fileIO/silo/SiloFile.cpp +++ b/src/coreComponents/fileIO/silo/SiloFile.cpp @@ -31,7 +31,7 @@ #include "constitutive/ConstitutiveManager.hpp" #include "constitutive/fluid/singlefluid/SingleFluidBase.hpp" #include "constitutive/fluid/multifluid/MultiFluidBase.hpp" -#include "constitutive/contact/ContactBase.hpp" +#include "constitutive/contact/FrictionBase.hpp" #include "constitutive/NullModel.hpp" #include "fileIO/Outputs/OutputUtilities.hpp" #include "mesh/DomainPartition.hpp" @@ -1468,7 +1468,7 @@ void SiloFile::writeElementMesh( ElementRegionBase const & elementRegion, localIndex const numFluids = regionFluidMaterialList.size(); string_array - fractureContactMaterialList = elementRegion.getConstitutiveNames< constitutive::ContactBase >(); + fractureContactMaterialList = elementRegion.getConstitutiveNames< constitutive::FrictionBase >(); localIndex const numContacts = fractureContactMaterialList.size(); diff --git a/src/coreComponents/physicsSolvers/contact/ContactSolverBase.cpp b/src/coreComponents/physicsSolvers/contact/ContactSolverBase.cpp index 2cef13dd221..7851dbdaa4d 100644 --- a/src/coreComponents/physicsSolvers/contact/ContactSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/contact/ContactSolverBase.cpp @@ -20,14 +20,8 @@ #include "ContactSolverBase.hpp" #include "common/TimingMacros.hpp" -#include "constitutive/ConstitutiveManager.hpp" -#include "constitutive/contact/ContactSelector.hpp" -#include "constitutive/solid/ElasticIsotropic.hpp" -#include "finiteElement/elementFormulations/FiniteElementBase.hpp" -#include "linearAlgebra/utilities/LAIHelperFunctions.hpp" +#include "constitutive/contact/FrictionBase.hpp" #include "mesh/DomainPartition.hpp" -#include "fieldSpecification/FieldSpecificationManager.hpp" -#include "mesh/NodeManager.hpp" #include "mesh/SurfaceElementRegion.hpp" #include "physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp" #include "common/GEOS_RAJA_Interface.hpp" @@ -244,15 +238,15 @@ void ContactSolverBase::setConstitutiveNamesCallSuper( ElementSubRegionBase & su } else if( dynamic_cast< SurfaceElementSubRegion * >( &subRegion ) ) { - subRegion.registerWrapper< string >( viewKeyStruct::contactRelationNameString() ). + subRegion.registerWrapper< string >( viewKeyStruct::frictionLawNameString() ). setPlotLevel( PlotLevel::NOPLOT ). setRestartFlags( RestartFlags::NO_WRITE ). setSizedFromParent( 0 ); - string & contactRelationName = subRegion.getReference< string >( viewKeyStruct::contactRelationNameString() ); - contactRelationName = SolverBase::getConstitutiveName< ContactBase >( subRegion ); - GEOS_ERROR_IF( contactRelationName.empty(), GEOS_FMT( "{}: ContactBase model not found on subregion {}", - getDataContext(), subRegion.getDataContext() ) ); + string & frictionLawName = subRegion.getReference< string >( viewKeyStruct::frictionLawNameString() ); + frictionLawName = SolverBase::getConstitutiveName< FrictionBase >( subRegion ); + GEOS_ERROR_IF( frictionLawName.empty(), GEOS_FMT( "{}: FrictionBase model not found on subregion {}", + getDataContext(), subRegion.getDataContext() ) ); } } diff --git a/src/coreComponents/physicsSolvers/contact/ContactSolverBase.hpp b/src/coreComponents/physicsSolvers/contact/ContactSolverBase.hpp index 1c4cf695ae6..d0d1b471ef3 100644 --- a/src/coreComponents/physicsSolvers/contact/ContactSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/contact/ContactSolverBase.hpp @@ -54,6 +54,9 @@ class ContactSolverBase : public SolidMechanicsLagrangianFEM constexpr static char const * fractureStateString() { return "fractureState"; } constexpr static char const * oldFractureStateString() { return "oldFractureState"; } + + constexpr static char const * frictionLawNameString() { return "frictionLawName"; } + }; protected: diff --git a/src/coreComponents/physicsSolvers/contact/SolidMechanicsEFEMKernels.hpp b/src/coreComponents/physicsSolvers/contact/SolidMechanicsEFEMKernels.hpp index 9be74bf112f..b5e530db12f 100644 --- a/src/coreComponents/physicsSolvers/contact/SolidMechanicsEFEMKernels.hpp +++ b/src/coreComponents/physicsSolvers/contact/SolidMechanicsEFEMKernels.hpp @@ -291,17 +291,18 @@ struct StateUpdateKernel /** * @brief Launch the kernel function doing fracture traction updates * @tparam POLICY the type of policy used in the kernel launch - * @tparam CONTACT_WRAPPER the type of contact wrapper doing the fracture traction updates + * @tparam FRICTION_WRAPPER the type of contact wrapper doing the fracture traction updates * @param[in] size the size of the subregion - * @param[in] contactWrapper the wrapper implementing the contact relationship + * @param[in] frictionWrapper the wrapper implementing the contact relationship * @param[in] jump the displacement jump * @param[out] fractureTraction the fracture traction * @param[out] dFractureTraction_dJump the derivative of the fracture traction wrt displacement jump */ - template< typename POLICY, typename CONTACT_WRAPPER > + template< typename POLICY, typename FRICTION_WRAPPER > static void launch( localIndex const size, - CONTACT_WRAPPER const & contactWrapper, + FRICTION_WRAPPER const & frictionWrapper, + real64 const contactPenaltyStiffness, arrayView2d< real64 const > const & oldJump, arrayView2d< real64 const > const & jump, arrayView2d< real64 > const & fractureTraction, @@ -311,10 +312,26 @@ struct StateUpdateKernel { forAll< POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const k ) { - contactWrapper.computeTraction( k, oldJump[k], jump[k], - fractureState[k], - fractureTraction[k], - dFractureTraction_dJump[k] ); + // Initialize traction and derivatives to 0 + LvArray::forValuesInSlice( fractureTraction[k], []( real64 & val ){ val = 0.0; } ); + LvArray::forValuesInSlice( dFractureTraction_dJump[k], []( real64 & val ){ val = 0.0; } ); + + // If the fracture is open the traction is 0 and so are its derivatives so there is nothing to do + bool const isOpen = fractureState[k] == fields::contact::FractureState::Open; + + if( !isOpen ) + { + // normal component of the traction + fractureTraction[k][0] = contactPenaltyStiffness * jump[k][0]; + + // derivative of the normal component w.r.t. to the normal dispJump + dFractureTraction_dJump[k][0][0] = contactPenaltyStiffness; + + frictionWrapper.computeShearTraction( k, oldJump[k], jump[k], + fractureState[k], + fractureTraction[k], + dFractureTraction_dJump[k] ); + } slip[ k ] = LvArray::math::sqrt( LvArray::math::square( jump( k, 1 ) ) + LvArray::math::square( jump( k, 2 ) ) ); diff --git a/src/coreComponents/physicsSolvers/contact/SolidMechanicsEmbeddedFractures.cpp b/src/coreComponents/physicsSolvers/contact/SolidMechanicsEmbeddedFractures.cpp index 6e1653861c2..0d6af0ee41e 100644 --- a/src/coreComponents/physicsSolvers/contact/SolidMechanicsEmbeddedFractures.cpp +++ b/src/coreComponents/physicsSolvers/contact/SolidMechanicsEmbeddedFractures.cpp @@ -22,7 +22,7 @@ #include "common/TimingMacros.hpp" #include "common/GEOS_RAJA_Interface.hpp" #include "constitutive/ConstitutiveManager.hpp" -#include "constitutive/contact/ContactSelector.hpp" +#include "constitutive/contact/FrictionSelector.hpp" #include "constitutive/solid/ElasticIsotropic.hpp" #include "fieldSpecification/FieldSpecificationManager.hpp" #include "finiteElement/elementFormulations/FiniteElementBase.hpp" @@ -51,6 +51,10 @@ SolidMechanicsEmbeddedFractures::SolidMechanicsEmbeddedFractures( const string & setInputFlag( InputFlags::OPTIONAL ). setApplyDefaultValue( 0 ). setDescription( "Defines whether to use static condensation or not." ); + + getWrapperBase( viewKeyStruct::contactPenaltyStiffnessString() ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "Value of the penetration penalty stiffness. Units of Pressure/length" ); } SolidMechanicsEmbeddedFractures::~SolidMechanicsEmbeddedFractures() @@ -718,8 +722,8 @@ void SolidMechanicsEmbeddedFractures::updateState( DomainPartition & domain ) { fractureRegion.forElementSubRegions< SurfaceElementSubRegion >( [&]( SurfaceElementSubRegion & subRegion ) { - string const & contactRelationName = subRegion.template getReference< string >( viewKeyStruct::contactRelationNameString() ); - ContactBase const & contact = getConstitutiveModel< ContactBase >( subRegion, contactRelationName ); + string const & frictionLawName = subRegion.template getReference< string >( viewKeyStruct::frictionLawNameString() ); + FrictionBase const & frictionLaw = getConstitutiveModel< FrictionBase >( subRegion, frictionLawName ); arrayView2d< real64 const > const & jump = subRegion.getField< contact::dispJump >(); @@ -733,14 +737,15 @@ void SolidMechanicsEmbeddedFractures::updateState( DomainPartition & domain ) arrayView1d< real64 > const & slip = subRegion.getField< fields::contact::slip >(); - constitutiveUpdatePassThru( contact, [&] ( auto & castedContact ) + constitutiveUpdatePassThru( frictionLaw, [&] ( auto & castedFrictionLaw ) { - using ContactType = TYPEOFREF( castedContact ); - typename ContactType::KernelWrapper contactWrapper = castedContact.createKernelWrapper(); + using FrictionType = TYPEOFREF( castedFrictionLaw ); + typename FrictionType::KernelWrapper frictionWrapper = castedFrictionLaw.createKernelWrapper(); solidMechanicsEFEMKernels::StateUpdateKernel:: launch< parallelDevicePolicy<> >( subRegion.size(), - contactWrapper, + frictionWrapper, + m_contactPenaltyStiffness, oldJump, jump, fractureTraction, @@ -765,13 +770,13 @@ bool SolidMechanicsEmbeddedFractures::updateConfiguration( DomainPartition & dom arrayView2d< real64 const > const & traction = subRegion.getField< fields::contact::traction >(); arrayView1d< integer > const & fractureState = subRegion.getField< fields::contact::fractureState >(); - string const & contactRelationName = subRegion.template getReference< string >( viewKeyStruct::contactRelationNameString() ); - ContactBase const & contact = getConstitutiveModel< ContactBase >( subRegion, contactRelationName ); + string const & frictionLawName = subRegion.template getReference< string >( viewKeyStruct::frictionLawNameString() ); + FrictionBase const & frictionLaw = getConstitutiveModel< FrictionBase >( subRegion, frictionLawName ); - constitutiveUpdatePassThru( contact, [&] ( auto & castedContact ) + constitutiveUpdatePassThru( frictionLaw, [&] ( auto & castedFrictionLaw ) { - using ContactType = TYPEOFREF( castedContact ); - typename ContactType::KernelWrapper contactWrapper = castedContact.createKernelWrapper(); + using FrictionType = TYPEOFREF( castedFrictionLaw ); + typename FrictionType::KernelWrapper frictionWrapper = castedFrictionLaw.createKernelWrapper(); RAJA::ReduceMin< parallelHostReduce, integer > checkActiveSetSub( 1 ); @@ -780,7 +785,7 @@ bool SolidMechanicsEmbeddedFractures::updateConfiguration( DomainPartition & dom if( ghostRank[kfe] < 0 ) { integer const originalFractureState = fractureState[kfe]; - contactWrapper.updateFractureState( kfe, dispJump[kfe], traction[kfe], fractureState[kfe] ); + frictionWrapper.updateFractureState( kfe, dispJump[kfe], traction[kfe], fractureState[kfe] ); checkActiveSetSub.min( compareFractureStates( originalFractureState, fractureState[kfe] ) ); } } ); diff --git a/src/coreComponents/physicsSolvers/contact/SolidMechanicsEmbeddedFractures.hpp b/src/coreComponents/physicsSolvers/contact/SolidMechanicsEmbeddedFractures.hpp index e4cf9d874b5..d601ccdf50d 100644 --- a/src/coreComponents/physicsSolvers/contact/SolidMechanicsEmbeddedFractures.hpp +++ b/src/coreComponents/physicsSolvers/contact/SolidMechanicsEmbeddedFractures.hpp @@ -118,6 +118,8 @@ class SolidMechanicsEmbeddedFractures : public ContactSolverBase struct viewKeyStruct : ContactSolverBase::viewKeyStruct { constexpr static char const * useStaticCondensationString() { return "useStaticCondensation"; } + + constexpr static char const * contactPenaltyStiffnessString() { return "contactPenaltyStiffness"; } }; protected: @@ -135,6 +137,9 @@ class SolidMechanicsEmbeddedFractures : public ContactSolverBase /// decide whether to use static condensation or not integer m_useStaticCondensation; + // TODO: activate when solidMechanicsPenalty contact is used and this is removed from base solver. + // real64 m_contactPenaltyStiffness; + }; diff --git a/src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContact.cpp b/src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContact.cpp index d4d537ab405..75ac7e16f0e 100644 --- a/src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContact.cpp +++ b/src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContact.cpp @@ -22,7 +22,7 @@ #include "common/TimingMacros.hpp" #include "constitutive/ConstitutiveManager.hpp" -#include "constitutive/contact/ContactSelector.hpp" +#include "constitutive/contact/FrictionSelector.hpp" #include "constitutive/fluid/singlefluid/SingleFluidBase.hpp" #include "finiteVolume/FiniteVolumeManager.hpp" #include "finiteVolume/FluxApproximationBase.hpp" @@ -1458,8 +1458,8 @@ void SolidMechanicsLagrangeContact:: [&]( localIndex const, FaceElementSubRegion const & subRegion ) { - string const & contactRelationName = subRegion.template getReference< string >( viewKeyStruct::contactRelationNameString() ); - ContactBase const & contact = getConstitutiveModel< ContactBase >( subRegion, contactRelationName ); + string const & frictionLawName = subRegion.template getReference< string >( viewKeyStruct::frictionLawNameString() ); + FrictionBase const & frictionLaw = getConstitutiveModel< FrictionBase >( subRegion, frictionLawName ); arrayView1d< globalIndex const > const & tracDofNumber = subRegion.getReference< globalIndex_array >( tracDofKey ); arrayView1d< integer const > const & ghostRank = subRegion.ghostRank(); @@ -1473,10 +1473,10 @@ void SolidMechanicsLagrangeContact:: arrayView2d< real64 const > const & previousDispJump = subRegion.getField< contact::oldDispJump >(); arrayView1d< real64 const > const & slidingTolerance = subRegion.getReference< array1d< real64 > >( viewKeyStruct::slidingToleranceString() ); - constitutiveUpdatePassThru( contact, [&] ( auto & castedContact ) + constitutiveUpdatePassThru( frictionLaw, [&] ( auto & castedFrictionLaw ) { - using ContactType = TYPEOFREF( castedContact ); - typename ContactType::KernelWrapper contactWrapper = castedContact.createKernelWrapper(); + using FrictionType = TYPEOFREF( castedFrictionLaw ); + typename FrictionType::KernelWrapper frictionWrapper = castedFrictionLaw.createKernelWrapper(); forAll< parallelHostPolicy >( subRegion.size(), [=] ( localIndex const kfe ) { @@ -1577,8 +1577,8 @@ void SolidMechanicsLagrangeContact:: } real64 dLimitTau_dNormalTraction = 0; - real64 const limitTau = contactWrapper.computeLimitTangentialTractionNorm( traction[kfe][0], - dLimitTau_dNormalTraction ); + real64 const limitTau = frictionWrapper.computeLimitTangentialTractionNorm( traction[kfe][0], + dLimitTau_dNormalTraction ); real64 sliding[ 2 ] = { dispJump[kfe][1] - previousDispJump[kfe][1], dispJump[kfe][2] - previousDispJump[kfe][2] }; real64 slidingNorm = sqrt( sliding[ 0 ]*sliding[ 0 ] + sliding[ 1 ]*sliding[ 1 ] ); @@ -2230,8 +2230,8 @@ bool SolidMechanicsLagrangeContact::updateConfiguration( DomainPartition & domai elemManager.forElementSubRegions< FaceElementSubRegion >( regionNames, [&]( localIndex const, FaceElementSubRegion & subRegion ) { - string const & contactRelationName = subRegion.template getReference< string >( viewKeyStruct::contactRelationNameString() ); - ContactBase const & contact = getConstitutiveModel< ContactBase >( subRegion, contactRelationName ); + string const & fricitonLawName = subRegion.template getReference< string >( viewKeyStruct::frictionLawNameString() ); + FrictionBase const & frictionLaw = getConstitutiveModel< FrictionBase >( subRegion, fricitonLawName ); arrayView1d< integer const > const & ghostRank = subRegion.ghostRank(); arrayView2d< real64 const > const & traction = subRegion.getField< contact::traction >(); @@ -2247,10 +2247,10 @@ bool SolidMechanicsLagrangeContact::updateConfiguration( DomainPartition & domai RAJA::ReduceSum< parallelHostReduce, real64 > changed( 0 ); RAJA::ReduceSum< parallelHostReduce, real64 > total( 0 ); - constitutiveUpdatePassThru( contact, [&] ( auto & castedContact ) + constitutiveUpdatePassThru( frictionLaw, [&] ( auto & castedFrictionLaw ) { - using ContactType = TYPEOFREF( castedContact ); - typename ContactType::KernelWrapper contactWrapper = castedContact.createKernelWrapper(); + using FrictionType = TYPEOFREF( castedFrictionLaw ); + typename FrictionType::KernelWrapper frictionWrapper = castedFrictionLaw.createKernelWrapper(); forAll< parallelHostPolicy >( subRegion.size(), [=] ( localIndex const kfe ) { @@ -2282,8 +2282,8 @@ bool SolidMechanicsLagrangeContact::updateConfiguration( DomainPartition & domai real64 dLimitTangentialTractionNorm_dTraction = 0.0; real64 const limitTau = - contactWrapper.computeLimitTangentialTractionNorm( traction[kfe][0], - dLimitTangentialTractionNorm_dTraction ); + frictionWrapper.computeLimitTangentialTractionNorm( traction[kfe][0], + dLimitTangentialTractionNorm_dTraction ); if( originalFractureState == FractureState::Stick && currentTau >= limitTau ) { diff --git a/src/coreComponents/physicsSolvers/contact/SolidMechanicsPenaltyContact.cpp b/src/coreComponents/physicsSolvers/contact/SolidMechanicsPenaltyContact.cpp new file mode 100644 index 00000000000..b635abb1b52 --- /dev/null +++ b/src/coreComponents/physicsSolvers/contact/SolidMechanicsPenaltyContact.cpp @@ -0,0 +1,235 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2020 TotalEnergies + * Copyright (c) 2019- GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file SolidMechanicsPenaltyContact.cpp + * + */ + +#include "SolidMechanicsPenaltyContact.hpp" + +#include "common/TimingMacros.hpp" +#include "constitutive/ConstitutiveManager.hpp" +#include "constitutive/contact/FrictionSelector.hpp" + + +#if defined( __INTEL_COMPILER ) +#pragma GCC optimize "O0" +#endif + +namespace geos +{ + +using namespace constitutive; +using namespace dataRepository; +using namespace fields; +using namespace finiteElement; + +SolidMechanicsPenaltyContact::SolidMechanicsPenaltyContact( const string & name, + Group * const parent ): + ContactSolverBase( name, parent ) +{ + registerWrapper( viewKeyStruct::contactPenaltyStiffnessString(), &m_contactPenaltyStiffness ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "Value of the penetration penalty stiffness. Units of Pressure/length" ); +} + +SolidMechanicsPenaltyContact::~SolidMechanicsPenaltyContact() +{ + // TODO Auto-generated destructor stub +} + +void SolidMechanicsPenaltyContact::setupSystem( DomainPartition & domain, + DofManager & dofManager, + CRSMatrix< real64, globalIndex > & localMatrix, + ParallelVector & rhs, + ParallelVector & solution, + bool const setSparsity ) +{ + GEOS_MARK_FUNCTION; + SolverBase::setupSystem( domain, dofManager, localMatrix, rhs, solution, false ); + + SparsityPattern< globalIndex > sparsityPattern( dofManager.numLocalDofs(), + dofManager.numGlobalDofs(), + 8*8*3*1.2 ); + + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & regionNames ) + { + NodeManager const & nodeManager = mesh.getNodeManager(); + arrayView1d< globalIndex const > const + dofNumber = nodeManager.getReference< globalIndex_array >( dofManager.getKey( solidMechanics::totalDisplacement::key() ) ); + + + ElementRegionManager const & elemManager = mesh.getElemManager(); + array1d< string > allFaceElementRegions; + elemManager.forElementRegions< SurfaceElementRegion >( [&]( SurfaceElementRegion const & elemRegion ) + { + allFaceElementRegions.emplace_back( elemRegion.getName() ); + } ); + + finiteElement:: + fillSparsity< FaceElementSubRegion, + solidMechanicsLagrangianFEMKernels::ImplicitSmallStrainQuasiStatic >( mesh, + allFaceElementRegions, + this->getDiscretizationName(), + dofNumber, + dofManager.rankOffset(), + sparsityPattern ); + + finiteElement::fillSparsity< CellElementSubRegion, + solidMechanicsLagrangianFEMKernels::ImplicitSmallStrainQuasiStatic >( mesh, + regionNames, + this->getDiscretizationName(), + dofNumber, + dofManager.rankOffset(), + sparsityPattern ); + + + } ); + + sparsityPattern.compress(); + localMatrix.assimilate< parallelDevicePolicy<> >( std::move( sparsityPattern ) ); +} + +void SolidMechanicsPenaltyContact::assembleSystem( real64 const time, + real64 const dt, + DomainPartition & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) +{ + GEOS_MARK_FUNCTION; + + synchronizeFractureState( domain ); + + SolidMechanicsLagrangianFEM::assembleSystem( time, + dt, + domain, + dofManager, + localMatrix, + localRhs ); + + assembleContact( domain, dofManager, localMatrix, localRhs ); +} + +void SolidMechanicsPenaltyContact::assembleContact( DomainPartition & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) +{ + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & ) + { + FaceManager const & faceManager = mesh.getFaceManager(); + NodeManager & nodeManager = mesh.getNodeManager(); + ElementRegionManager & elemManager = mesh.getElemManager(); + + solidMechanics::arrayViewConst2dLayoutTotalDisplacement const u = + nodeManager.getField< solidMechanics::totalDisplacement >(); + arrayView2d< real64 > const fc = nodeManager.getField< solidMechanics::contactForce >(); + fc.zero(); + + arrayView2d< real64 const > const faceNormal = faceManager.faceNormal(); + ArrayOfArraysView< localIndex const > const facesToNodes = faceManager.nodeList().toViewConst(); + + string const dofKey = dofManager.getKey( solidMechanics::totalDisplacement::key() ); + arrayView1d< globalIndex > const nodeDofNumber = nodeManager.getReference< globalIndex_array >( dofKey ); + globalIndex const rankOffset = dofManager.rankOffset(); + + // TODO: this bound may need to change + constexpr localIndex maxNodexPerFace = 4; + constexpr localIndex maxDofPerElem = maxNodexPerFace * 3 * 2; + + elemManager.forElementSubRegions< FaceElementSubRegion >( [&]( FaceElementSubRegion & subRegion ) + { + real64 const contactStiffness = m_contactPenaltyStiffness; + + arrayView1d< real64 > const area = subRegion.getElementArea(); + ArrayOfArraysView< localIndex const > const elemsToFaces = subRegion.faceList().toViewConst(); + + // TODO: use parallel policy? + forAll< serialPolicy >( subRegion.size(), [=] ( localIndex const kfe ) + { + localIndex const kf0 = elemsToFaces[kfe][0], kf1 = elemsToFaces[kfe][1]; + real64 Nbar[ 3 ] = { faceNormal[kf0][0] - faceNormal[kf1][0], + faceNormal[kf0][1] - faceNormal[kf1][1], + faceNormal[kf0][2] - faceNormal[kf1][2] }; + + LvArray::tensorOps::normalize< 3 >( Nbar ); + + localIndex const numNodesPerFace=facesToNodes.sizeOfArray( kf0 ); + real64 const Ja = area[kfe] / numNodesPerFace; + + stackArray1d< globalIndex, maxDofPerElem > rowDOF( numNodesPerFace*3*2 ); + stackArray1d< real64, maxDofPerElem > nodeRHS( numNodesPerFace*3*2 ); + stackArray2d< real64, maxDofPerElem *maxDofPerElem > dRdP( numNodesPerFace*3*2, numNodesPerFace*3*2 ); + + for( localIndex a=0; a( gap, u[node0] ); + real64 const gapNormal = LvArray::tensorOps::AiBi< 3 >( gap, Nbar ); + + for( int i=0; i<3; ++i ) + { + rowDOF[3*a+i] = nodeDofNumber[node0]+i; + rowDOF[3*(numNodesPerFace + a)+i] = nodeDofNumber[node1]+i; + } + + if( gapNormal < 0 ) + { + LvArray::tensorOps::scale< 3 >( penaltyForce, -contactStiffness * gapNormal * Ja ); + for( int i=0; i<3; ++i ) + { + LvArray::tensorOps::subtract< 3 >( fc[node0], penaltyForce ); + LvArray::tensorOps::add< 3 >( fc[node1], penaltyForce ); + nodeRHS[3*a+i] -= penaltyForce[i]; + nodeRHS[3*(numNodesPerFace + a)+i] += penaltyForce[i]; + + dRdP( 3*a+i, 3*a+i ) -= contactStiffness * Ja * Nbar[i] * Nbar[i]; + dRdP( 3*a+i, 3*(numNodesPerFace + a)+i ) += contactStiffness * Ja * Nbar[i] * Nbar[i]; + dRdP( 3*(numNodesPerFace + a)+i, 3*a+i ) += contactStiffness * Ja * Nbar[i] * Nbar[i]; + dRdP( 3*(numNodesPerFace + a)+i, 3*(numNodesPerFace + a)+i ) -= contactStiffness * Ja * Nbar[i] * Nbar[i]; + } + } + } + + for( localIndex idof = 0; idof < numNodesPerFace*3*2; ++idof ) + { + localIndex const localRow = LvArray::integerConversion< localIndex >( rowDOF[idof] - rankOffset ); + + if( localRow >= 0 && localRow < localMatrix.numRows() ) + { + localMatrix.addToRowBinarySearchUnsorted< serialAtomic >( localRow, + rowDOF.data(), + dRdP[idof].dataIfContiguous(), + numNodesPerFace*3*2 ); + RAJA::atomicAdd( serialAtomic{}, &localRhs[localRow], nodeRHS[idof] ); + } + } + } ); + } ); + } ); +} + + +REGISTER_CATALOG_ENTRY( SolverBase, SolidMechanicsPenaltyContact, string const &, Group * const ) + +} /* namespace geos */ diff --git a/src/coreComponents/physicsSolvers/contact/SolidMechanicsPenaltyContact.hpp b/src/coreComponents/physicsSolvers/contact/SolidMechanicsPenaltyContact.hpp new file mode 100644 index 00000000000..d368fc15ebc --- /dev/null +++ b/src/coreComponents/physicsSolvers/contact/SolidMechanicsPenaltyContact.hpp @@ -0,0 +1,92 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2020 TotalEnergies + * Copyright (c) 2019- GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file SolidMechanicsPenaltyContact.hpp + * + */ + +#ifndef GEOS_PHYSICSSOLVERS_CONTACT_SOLIDMECHANICSPENALTYCONTACT_HPP_ +#define GEOS_PHYSICSSOLVERS_CONTACT_SOLIDMECHANICSPENALTYCONTACT_HPP_ + +#include "physicsSolvers/contact/ContactSolverBase.hpp" +#include "../../linearAlgebra/DofManager.hpp" +#include "../../common/DataTypes.hpp" + +namespace geos +{ + +class SolidMechanicsLagrangianFEM; + +class SolidMechanicsPenaltyContact : public ContactSolverBase +{ +public: + + SolidMechanicsPenaltyContact( const string & name, + Group * const parent ); + + ~SolidMechanicsPenaltyContact() override; + + /** + * @brief name of the node manager in the object catalog + * @return string that contains the catalog name to generate a new NodeManager object through the object catalog. + */ + static string catalogName() + { + return "SolidMechanicsPenaltyContact"; + } + /** + * @copydoc SolverBase::getCatalogName() + */ + string getCatalogName() const override { return catalogName(); } + + /// String used to form the solverName used to register single-physics solvers in CoupledSolver + static string coupledSolverAttributePrefix() { return "PenaltyContact"; } + + virtual void + setupSystem( DomainPartition & domain, + DofManager & dofManager, + CRSMatrix< real64, globalIndex > & localMatrix, + ParallelVector & rhs, + ParallelVector & solution, + bool const setSparsity = true ) override final; + + virtual void + assembleSystem( real64 const time, + real64 const dt, + DomainPartition & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) override; + + void assembleContact( DomainPartition & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ); + +protected: + +private: + + struct viewKeyStruct : ContactSolverBase::viewKeyStruct + { + constexpr static char const * contactPenaltyStiffnessString() { return "contactPenaltyStiffness"; } + }; + + real64 m_contactPenaltyStiffness; +}; + +} /* namespace geos */ + +#endif /* GEOS_PHYSICSSOLVERS_CONTACT_SOLIDMECHANICSPENALTYCONTACT_HPP_ */ diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp index 8c8acc2ebcb..3e63e5e00c2 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp @@ -22,6 +22,7 @@ #include "constitutive/ConstitutivePassThru.hpp" #include "constitutive/permeability/PermeabilityFields.hpp" #include "constitutive/solid/SolidInternalEnergy.hpp" +#include "constitutive/contact/HydraulicApertureBase.hpp" #include "discretizationMethods/NumericalMethodsManager.hpp" #include "fieldSpecification/AquiferBoundaryCondition.hpp" #include "fieldSpecification/EquilibriumInitialCondition.hpp" diff --git a/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.cpp b/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.cpp index 5743772a144..3c108c21daa 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.cpp @@ -6,7 +6,7 @@ * Copyright (c) 2018-2024 Total, S.A * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University * Copyright (c) 2018-2024 Chevron - * Copyright (c) 2019- GEOS/GEOSX Contributors + * Copyright (c) 2019- GEOS/GEOS Contributors * All rights reserved * * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. @@ -17,10 +17,9 @@ * @file HydrofractureSolver.cpp */ - #include "HydrofractureSolver.hpp" -#include "constitutive/contact/ContactSelector.hpp" +#include "constitutive/contact/HydraulicApertureRelationSelector.hpp" #include "constitutive/fluid/singlefluid/SingleFluidBase.hpp" #include "physicsSolvers/multiphysics/HydrofractureSolverKernels.hpp" #include "physicsSolvers/solidMechanics/SolidMechanicsFields.hpp" @@ -41,7 +40,6 @@ template< typename POROMECHANICS_SOLVER > HydrofractureSolver< POROMECHANICS_SOLVER >::HydrofractureSolver( const string & name, Group * const parent ) : Base( name, parent ), - m_contactRelationName(), m_surfaceGeneratorName(), m_surfaceGenerator( nullptr ), m_maxNumResolves( 10 ), @@ -53,11 +51,6 @@ HydrofractureSolver< POROMECHANICS_SOLVER >::HydrofractureSolver( const string & setInputFlag( InputFlags::REQUIRED ). setDescription( "Name of the surface generator to use in the hydrofracture solver" ); - registerWrapper( viewKeyStruct::contactRelationNameString(), &m_contactRelationName ). - setRTTypeName( rtTypes::CustomTypes::groupNameRef ). - setInputFlag( InputFlags::REQUIRED ). - setDescription( "Name of contact relation to enforce constraints on fracture boundary." ); - registerWrapper( viewKeyStruct::maxNumResolvesString(), &m_maxNumResolves ). setApplyDefaultValue( 10 ). setInputFlag( InputFlags::OPTIONAL ). @@ -288,7 +281,8 @@ void HydrofractureSolver< POROMECHANICS_SOLVER >::updateHydraulicApertureAndFrac elemManager.forElementSubRegions< FaceElementSubRegion >( [&]( FaceElementSubRegion & subRegion ) { - ContactBase const & contact = this->template getConstitutiveModel< ContactBase >( subRegion, m_contactRelationName ); + string const & hydraulicApertureRelationName = subRegion.template getReference< string >( viewKeyStruct::hydraulicApertureRelationNameString() ); + HydraulicApertureBase const & hydraulicApertureModel = this->template getConstitutiveModel< HydraulicApertureBase >( subRegion, hydraulicApertureRelationName ); arrayView1d< real64 > const aperture = subRegion.getElementAperture(); arrayView1d< real64 > const hydraulicAperture = subRegion.getField< flow::hydraulicAperture >(); @@ -316,14 +310,14 @@ void HydrofractureSolver< POROMECHANICS_SOLVER >::updateHydraulicApertureAndFrac { typename TYPEOFREF( castedPorousSolid ) ::KernelWrapper porousMaterialWrapper = castedPorousSolid.createKernelUpdates(); - constitutiveUpdatePassThru( contact, [&] ( auto & castedContact ) + constitutiveUpdatePassThru( hydraulicApertureModel, [&] ( auto & castedHydraulicApertureModel ) { - using ContactType = TYPEOFREF( castedContact ); - typename ContactType::KernelWrapper contactWrapper = castedContact.createKernelWrapper(); + using HydraulicApertureModelType = TYPEOFREF( castedHydraulicApertureModel ); + typename HydraulicApertureModelType::KernelWrapper hydraulicApertureWrapper = castedHydraulicApertureModel.createKernelWrapper(); auto const statistics = hydrofractureSolverKernels::DeformationUpdateKernel ::launch< parallelDevicePolicy<> >( subRegion.size(), - contactWrapper, + hydraulicApertureWrapper, porousMaterialWrapper, u, faceNormal, @@ -809,7 +803,8 @@ assembleFluidMassResidualDerivativeWrtDisplacement( DomainPartition const & doma [&]( localIndex const, FaceElementSubRegion const & subRegion ) { - ContactBase const & contact = this->template getConstitutiveModel< ContactBase >( subRegion, m_contactRelationName ); + string const & hydraulicApertureRelationName = subRegion.template getReference< string >( viewKeyStruct::hydraulicApertureRelationNameString() ); + HydraulicApertureBase const & hydraulicApertureModel = this->template getConstitutiveModel< HydraulicApertureBase >( subRegion, hydraulicApertureRelationName ); string const & fluidName = subRegion.getReference< string >( FlowSolverBase::viewKeyStruct::fluidNamesString() ); SingleFluidBase const & fluid = this->template getConstitutiveModel< SingleFluidBase >( subRegion, fluidName ); @@ -827,15 +822,15 @@ assembleFluidMassResidualDerivativeWrtDisplacement( DomainPartition const & doma arrayView2d< real64 const > const faceNormal = faceManager.faceNormal(); - constitutiveUpdatePassThru( contact, [&] ( auto & castedContact ) + constitutiveUpdatePassThru( hydraulicApertureModel, [&] ( auto & castedHydraulicApertureModel ) { - using ContactType = TYPEOFREF( castedContact ); - typename ContactType::KernelWrapper contactWrapper = castedContact.createKernelWrapper(); + using HydraulicApertureModelType = TYPEOFREF( castedHydraulicApertureModel ); + typename HydraulicApertureModelType::KernelWrapper hydraulicApertureWrapper = castedHydraulicApertureModel.createKernelWrapper(); hydrofractureSolverKernels::FluidMassResidualDerivativeAssemblyKernel:: launch< parallelDevicePolicy<> >( subRegion.size(), rankOffset, - contactWrapper, + hydraulicApertureWrapper, m_useQuasiNewton, elemsToFaces, faceToNodeMap, @@ -1049,7 +1044,7 @@ void HydrofractureSolver< POROMECHANICS_SOLVER >::initializeNewFractureFields( D localIndex const newElemIndex = newFractureElements[k]; real64 initialPressure = 1.0e99; real64 initialAperture = 1.0e99; - #ifdef GEOSX_USE_SEPARATION_COEFFICIENT + #ifdef GEOS_USE_SEPARATION_COEFFICIENT apertureF[newElemIndex] = aperture[newElemIndex]; #endif if( m_newFractureInitializationType == InitializationType::Displacement ) diff --git a/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.hpp b/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.hpp index 9df2e20b199..c9598839401 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.hpp @@ -163,8 +163,6 @@ class HydrofractureSolver : public POROMECHANICS_SOLVER struct viewKeyStruct : Base::viewKeyStruct { - constexpr static char const * contactRelationNameString() { return "contactRelationName"; } - constexpr static char const * surfaceGeneratorNameString() { return "surfaceGeneratorName"; } constexpr static char const * maxNumResolvesString() { return "maxNumResolves"; } @@ -175,7 +173,7 @@ class HydrofractureSolver : public POROMECHANICS_SOLVER constexpr static char const * useQuasiNewtonString() { return "useQuasiNewton"; } - static constexpr char const * isLaggingFractureStencilWeightsUpdateString() { return "isLaggingFractureStencilWeightsUpdate"; } + constexpr static char const * isLaggingFractureStencilWeightsUpdateString() { return "isLaggingFractureStencilWeightsUpdate"; } #ifdef GEOS_USE_SEPARATION_COEFFICIENT constexpr static char const * separationCoeff0String() { return "separationCoeff0"; } diff --git a/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolverKernels.hpp b/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolverKernels.hpp index cd193042239..e555a51d551 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolverKernels.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolverKernels.hpp @@ -32,10 +32,10 @@ namespace hydrofractureSolverKernels struct DeformationUpdateKernel { - template< typename POLICY, typename CONTACT_WRAPPER, typename POROUS_WRAPPER > + template< typename POLICY, typename HYDRAULICAPERTURE_WRAPPER, typename POROUS_WRAPPER > static std::tuple< double, double, double, double, double, double > launch( localIndex const size, - CONTACT_WRAPPER const & contactWrapper, + HYDRAULICAPERTURE_WRAPPER const & hydraulicApertureWrapper, POROUS_WRAPPER const & porousMaterialWrapper, arrayView2d< real64 const, nodes::TOTAL_DISPLACEMENT_USD > const & u, arrayView2d< real64 const > const & faceNormal, @@ -86,8 +86,14 @@ struct DeformationUpdateKernel minAperture.min( aperture[kfe] ); maxAperture.max( aperture[kfe] ); - real64 dHydraulicAperture_dNormalJump = 0; - real64 const newHydraulicAperture = contactWrapper.computeHydraulicAperture( aperture[kfe], dHydraulicAperture_dNormalJump ); + real64 normalTraction = 0.0; /// TODO: must be changed to use actual traction + real64 dHydraulicAperture_dNormalTraction = 0.0; + real64 dHydraulicAperture_dNormalJump = 0.0; + real64 const newHydraulicAperture = hydraulicApertureWrapper.computeHydraulicAperture( aperture[kfe], + normalTraction, + dHydraulicAperture_dNormalJump, + dHydraulicAperture_dNormalTraction ); + maxHydraulicApertureChange.max( std::fabs( newHydraulicAperture - hydraulicAperture[kfe] )); real64 const oldHydraulicAperture = hydraulicAperture[kfe]; hydraulicAperture[kfe] = newHydraulicAperture; @@ -127,11 +133,11 @@ struct DeformationUpdateKernel struct FluidMassResidualDerivativeAssemblyKernel { - template< typename CONTACT_WRAPPER > + template< typename HYDRAULICAPERTURE_WRAPPER > GEOS_HOST_DEVICE inline static void - computeAccumulationDerivative( CONTACT_WRAPPER const & contactWrapper, + computeAccumulationDerivative( HYDRAULICAPERTURE_WRAPPER const & hydraulicApertureWrapper, localIndex const numNodesPerFace, arraySlice1d< localIndex const > const elemsToFaces, ArrayOfArraysView< localIndex const > const faceToNodeMap, @@ -143,8 +149,13 @@ struct FluidMassResidualDerivativeAssemblyKernel globalIndex (& nodeDOF)[8 * 3], arraySlice1d< real64 > const dRdU ) { - real64 dHydraulicAperture_dNormalJump = 0; - real64 const hydraulicAperture = contactWrapper.computeHydraulicAperture( aperture, dHydraulicAperture_dNormalJump ); + real64 dHydraulicAperture_dNormalJump = 0.0; + real64 dHydraulicAperture_dTraction = 0.0; + real64 fractureTraction = 0.0; + real64 const hydraulicAperture = hydraulicApertureWrapper.computeHydraulicAperture( aperture, + fractureTraction, + dHydraulicAperture_dNormalJump, + dHydraulicAperture_dTraction ); GEOS_UNUSED_VAR( hydraulicAperture ); constexpr integer kfSign[2] = { -1, 1 }; @@ -201,11 +212,11 @@ struct FluidMassResidualDerivativeAssemblyKernel } } - template< typename POLICY, typename CONTACT_WRAPPER > + template< typename POLICY, typename HYDRAULICAPERTURE_WRAPPER > static void launch( localIndex const size, globalIndex const rankOffset, - CONTACT_WRAPPER const & contactWrapper, + HYDRAULICAPERTURE_WRAPPER const & hydraulicApertureWrapper, integer const useQuasiNewton, ArrayOfArraysView< localIndex const > const elemsToFaces, ArrayOfArraysView< localIndex const > const faceToNodeMap, @@ -230,7 +241,7 @@ struct FluidMassResidualDerivativeAssemblyKernel globalIndex nodeDOF[8 * 3]; stackArray1d< real64, 24 > dRdU( 2 * numNodesPerFace * 3 ); // - computeAccumulationDerivative( contactWrapper, + computeAccumulationDerivative( hydraulicApertureWrapper, numNodesPerFace, elemsToFaces[ei], faceToNodeMap, diff --git a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp index 9030f995508..bf6446aa0a6 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp @@ -27,6 +27,7 @@ #include "physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp" #include "constitutive/solid/CoupledSolidBase.hpp" #include "constitutive/solid/PorousSolid.hpp" +#include "constitutive/contact/HydraulicApertureBase.hpp" #include "mesh/DomainPartition.hpp" #include "mesh/utilities/AverageOverQuadraturePointsKernel.hpp" #include "codingUtilities/Utilities.hpp" @@ -123,6 +124,23 @@ class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER InputError ); } + virtual void setConstitutiveNamesCallSuper( ElementSubRegionBase & subRegion ) const override final + { + if( dynamic_cast< SurfaceElementSubRegion * >( &subRegion ) ) + { + subRegion.registerWrapper< string >( viewKeyStruct::hydraulicApertureRelationNameString() ). + setPlotLevel( PlotLevel::NOPLOT ). + setRestartFlags( RestartFlags::NO_WRITE ). + setSizedFromParent( 0 ); + + string & hydraulicApertureModelName = subRegion.getReference< string >( viewKeyStruct::hydraulicApertureRelationNameString() ); + hydraulicApertureModelName = SolverBase::getConstitutiveName< HydraulicApertureBase >( subRegion ); + GEOS_ERROR_IF( hydraulicApertureModelName.empty(), GEOS_FMT( "{}: HydraulicApertureBase model not found on subregion {}", + this->getDataContext(), subRegion.getDataContext() ) ); + } + + } + virtual void initializePreSubGroups() override { Base::initializePreSubGroups(); @@ -311,6 +329,10 @@ class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER /// Multiplier on stabilization strength constexpr static const char * stabilizationMultiplierString() {return "stabilizationMultiplier"; } + + /// Name of the hydraulicApertureRelationName + static constexpr char const * hydraulicApertureRelationNameString() {return "hydraulicApertureRelationName"; } + }; void updateStabilizationParameters( DomainPartition & domain ) const diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.cpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.cpp index abd76e855ad..2f0a16dc9ef 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.cpp @@ -20,12 +20,11 @@ #include "SinglePhasePoromechanicsConformingFractures.hpp" #include "constitutive/solid/PorousSolid.hpp" -#include "constitutive/contact/ContactSelector.hpp" #include "constitutive/fluid/singlefluid/SingleFluidBase.hpp" #include "linearAlgebra/solvers/BlockPreconditioner.hpp" #include "linearAlgebra/solvers/SeparateComponentPreconditioner.hpp" +#include "constitutive/contact/HydraulicApertureRelationSelector.hpp" #include "physicsSolvers/fluidFlow/SinglePhaseBase.hpp" -//#include "physicsSolvers/multiphysics/SinglePhaseReservoirAndWells.hpp" #include "physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanics.hpp" #include "physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanics.hpp" #include "physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsFractures.hpp" @@ -769,13 +768,13 @@ void SinglePhasePoromechanicsConformingFractures< FLOW_SOLVER >::updateHydraulic string const porousSolidName = subRegion.getReference< string >( FlowSolverBase::viewKeyStruct::solidNamesString() ); CoupledSolidBase & porousSolid = subRegion.getConstitutiveModel< CoupledSolidBase >( porousSolidName ); - string const & contactRelationName = subRegion.template getReference< string >( SolidMechanicsLagrangianFEM::viewKeyStruct::contactRelationNameString() ); - ContactBase const & contact = subRegion.getConstitutiveModel< ContactBase >( contactRelationName ); + string const & hydraulicApertureRelationName = subRegion.template getReference< string >( viewKeyStruct::hydraulicApertureRelationNameString() ); + HydraulicApertureBase const & hydraulicApertureModel = this->template getConstitutiveModel< HydraulicApertureBase >( subRegion, hydraulicApertureRelationName ); - constitutiveUpdatePassThru( contact, [&] ( auto & castedContact ) + constitutiveUpdatePassThru( hydraulicApertureModel, [&] ( auto & castedHydraulicAperture ) { - using ContactType = TYPEOFREF( castedContact ); - typename ContactType::KernelWrapper contactWrapper = castedContact.createKernelWrapper(); + using HydraulicApertureType = TYPEOFREF( castedHydraulicAperture ); + typename HydraulicApertureType::KernelWrapper hydraulicApertureWrapper = castedHydraulicAperture.createKernelWrapper(); constitutive::ConstitutivePassThru< CompressibleSolidBase >::execute( porousSolid, [=, &subRegion] ( auto & castedPorousSolid ) { @@ -784,7 +783,7 @@ void SinglePhasePoromechanicsConformingFractures< FLOW_SOLVER >::updateHydraulic poromechanicsFracturesKernels::StateUpdateKernel:: launch< parallelDevicePolicy<> >( subRegion.size(), porousMaterialWrapper, - contactWrapper, + hydraulicApertureWrapper, dispJump, pressure, area, diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.hpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.hpp index 124c2c261dd..72cfd3d68ca 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.hpp @@ -107,6 +107,9 @@ class SinglePhasePoromechanicsConformingFractures : public SinglePhasePoromechan private: + struct viewKeyStruct : public Base::viewKeyStruct + {}; + static const localIndex m_maxFaceNodes=11; // Maximum number of nodes on a contact face void assembleElementBasedContributions( real64 const time_n, diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.cpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.cpp index deb60634c18..79dba986aa4 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.cpp @@ -18,7 +18,7 @@ */ #include "SinglePhasePoromechanicsEmbeddedFractures.hpp" -#include "constitutive/contact/ContactSelector.hpp" +#include "constitutive/contact/HydraulicApertureRelationSelector.hpp" #include "constitutive/fluid/singlefluid/SingleFluidBase.hpp" #include "physicsSolvers/contact/SolidMechanicsEFEMKernelsHelper.hpp" #include "physicsSolvers/fluidFlow/SinglePhaseBase.hpp" @@ -496,33 +496,38 @@ void SinglePhasePoromechanicsEmbeddedFractures::updateState( DomainPartition & d arrayView1d< real64 const > const & pressure = subRegion.template getField< fields::flow::pressure >(); - string const & contactRelationName = subRegion.template getReference< string >( ContactSolverBase::viewKeyStruct::contactRelationNameString() ); - ContactBase const & contact = getConstitutiveModel< ContactBase >( subRegion, contactRelationName ); - - ContactBase::KernelWrapper contactWrapper = contact.createKernelWrapper(); + string const & hydraulicApertureRelationName = subRegion.template getReference< string >( viewKeyStruct::hydraulicApertureRelationNameString() ); + HydraulicApertureBase const & hydraulicApertureModel = this->template getConstitutiveModel< HydraulicApertureBase >( subRegion, hydraulicApertureRelationName ); string const porousSolidName = subRegion.template getReference< string >( FlowSolverBase::viewKeyStruct::solidNamesString() ); CoupledSolidBase & porousSolid = subRegion.template getConstitutiveModel< CoupledSolidBase >( porousSolidName ); - constitutive::ConstitutivePassThru< CompressibleSolidBase >::execute( porousSolid, [=, &subRegion] ( auto & castedPorousSolid ) + constitutive::ConstitutivePassThru< CompressibleSolidBase >::execute( porousSolid, [=, &subRegion, &hydraulicApertureModel] ( auto & castedPorousSolid ) { typename TYPEOFREF( castedPorousSolid ) ::KernelWrapper porousMaterialWrapper = castedPorousSolid.createKernelUpdates(); - poromechanicsEFEMKernels::StateUpdateKernel:: - launch< parallelDevicePolicy<> >( subRegion.size(), - contactWrapper, - porousMaterialWrapper, - dispJump, - pressure, - area, - volume, - deltaVolume, - aperture, - oldHydraulicAperture, - hydraulicAperture, - fractureTraction, - dTdpf ); + constitutiveUpdatePassThru( hydraulicApertureModel, [=, &subRegion] ( auto & castedHydraulicApertureModel ) + { + using HydraulicApertureModelType = TYPEOFREF( castedHydraulicApertureModel ); + typename HydraulicApertureModelType::KernelWrapper hydraulicApertureModelWrapper = castedHydraulicApertureModel.createKernelWrapper(); + + poromechanicsEFEMKernels::StateUpdateKernel:: + launch< parallelDevicePolicy<> >( subRegion.size(), + hydraulicApertureModelWrapper, + porousMaterialWrapper, + dispJump, + pressure, + area, + volume, + deltaVolume, + aperture, + oldHydraulicAperture, + hydraulicAperture, + fractureTraction, + dTdpf ); + + } ); } ); // update the stencil weights using the updated hydraulic aperture diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM.hpp index 3afa2f8e64d..ed4f3c07813 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsEFEM.hpp @@ -20,7 +20,6 @@ #ifndef GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSKERNELS_SINGLEPHASEPOROMECHANICSEFEM_HPP_ #define GEOS_PHYSICSSOLVERS_MULTIPHYSICS_POROMECHANICSKERNELS_SINGLEPHASEPOROMECHANICSEFEM_HPP_ -#include "constitutive/contact/ContactBase.hpp" #include "finiteElement/kernelInterface/ImplicitKernelBase.hpp" namespace geos @@ -331,10 +330,10 @@ struct StateUpdateKernel * @param[out] fractureTraction the fracture traction * @param[out] dFractureTraction_dPressure the derivative of the fracture traction wrt pressure */ - template< typename POLICY, typename POROUS_WRAPPER > + template< typename POLICY, typename POROUS_WRAPPER, typename CONTACT_WRAPPER > static void launch( localIndex const size, - constitutive::ContactBase::KernelWrapper const & contactWrapper, + CONTACT_WRAPPER const & contactWrapper, POROUS_WRAPPER const & porousMaterialWrapper, arrayView2d< real64 const > const & dispJump, arrayView1d< real64 const > const & pressure, @@ -352,16 +351,18 @@ struct StateUpdateKernel // update aperture to be equal to the normal displacement jump aperture[k] = dispJump[k][0]; // the first component of the jump is the normal one. - real64 dHydraulicAperture_dNormalJump = 0; + real64 dHydraulicAperture_dNormalJump = 0.0; + real64 dHydraulicAperture_dNormalTraction = 0.0; hydraulicAperture[k] = contactWrapper.computeHydraulicAperture( aperture[k], - dHydraulicAperture_dNormalJump ); + fractureTraction[k][0], + dHydraulicAperture_dNormalJump, + dHydraulicAperture_dNormalTraction ); deltaVolume[k] = hydraulicAperture[k] * area[k] - volume[k]; // traction on the fracture to include the pressure contribution - contactWrapper.addPressureToTraction( pressure[k], - fractureTraction[k], - dFractureTraction_dPressure[k] ); + fractureTraction[k][0] -= pressure[k]; + dFractureTraction_dPressure[k] = -1.0; real64 const jump[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3 ( dispJump[k] ); real64 const traction[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3 ( fractureTraction[k] ); diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsFractures.hpp b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsFractures.hpp index 58b6ab6d141..82bfd7a76b6 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsFractures.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanicsFractures.hpp @@ -35,7 +35,7 @@ struct StateUpdateKernel /** * @brief Launch the kernel function doing volume, aperture and fracture traction updates * @tparam POLICY the type of policy used in the kernel launch - * @tparam CONTACT_WRAPPER the type of contact wrapper doing the fracture traction updates + * @tparam HYDRAULIC_APERTURE_WRAPPER the type of hydraulic aperture model wrapper doing the fracture traction updates * @param[in] size the size of the subregion * @param[in] contactWrapper the wrapper implementing the contact relationship * @param[in] dispJump the displacement jump @@ -48,11 +48,11 @@ struct StateUpdateKernel * @param[out] hydraulicAperture the effecture aperture * @param[in] fractureTraction the fracture traction */ - template< typename POLICY, typename POROUS_WRAPPER, typename CONTACT_WRAPPER > + template< typename POLICY, typename POROUS_WRAPPER, typename HYDRAULIC_APERTURE_WRAPPER > static void launch( localIndex const size, POROUS_WRAPPER const & porousMaterialWrapper, - CONTACT_WRAPPER const & contactWrapper, + HYDRAULIC_APERTURE_WRAPPER const & contactWrapper, arrayView2d< real64 const > const & dispJump, arrayView1d< real64 const > const & pressure, arrayView1d< real64 const > const & area, @@ -70,7 +70,11 @@ struct StateUpdateKernel aperture[k] = dispJump[k][0]; // the first component of the jump is the normal one. real64 dHydraulicAperture_dNormalJump = 0.0; - hydraulicAperture[k] = contactWrapper.computeHydraulicAperture( aperture[k], dHydraulicAperture_dNormalJump ); + real64 dHydraulicAperture_dNormalTraction = 0.0; + hydraulicAperture[k] = contactWrapper.computeHydraulicAperture( aperture[k], + fractureTraction[k][0], + dHydraulicAperture_dNormalJump, + dHydraulicAperture_dNormalTraction ); deltaVolume[k] = hydraulicAperture[k] * area[k] - volume[k]; @@ -78,9 +82,11 @@ struct StateUpdateKernel real64 const traction[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3 ( fractureTraction[k] ); porousMaterialWrapper.updateStateFromPressureApertureJumpAndTraction( k, 0, pressure[k], - oldHydraulicAperture[k], hydraulicAperture[k], + oldHydraulicAperture[k], + hydraulicAperture[k], dHydraulicAperture_dNormalJump, - jump, traction ); + jump, + traction ); } ); } diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp index 23477cbf51c..70756a953d7 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp @@ -28,7 +28,6 @@ #include "codingUtilities/Utilities.hpp" #include "constitutive/ConstitutiveManager.hpp" -#include "constitutive/contact/ContactBase.hpp" #include "common/GEOS_RAJA_Interface.hpp" #include "discretizationMethods/NumericalMethodsManager.hpp" #include "fieldSpecification/FieldSpecificationManager.hpp" @@ -117,6 +116,11 @@ SolidMechanicsLagrangianFEM::SolidMechanicsLagrangianFEM( const string & name, setInputFlag( InputFlags::OPTIONAL ). setDescription( "Name of contact relation to enforce constraints on fracture boundary." ); + registerWrapper( viewKeyStruct::contactPenaltyStiffnessString(), &m_contactPenaltyStiffness ). + setInputFlag( InputFlags::OPTIONAL ). + setApplyDefaultValue( 0.0 ). + setDescription( "Value of the penetration penalty stiffness. Units of Pressure/length" ); + registerWrapper( viewKeyStruct::maxForceString(), &m_maxForce ). setInputFlag( InputFlags::FALSE ). setDescription( "The maximum force contribution in the problem domain." ); @@ -1378,9 +1382,7 @@ void SolidMechanicsLagrangianFEM::applyContactConstraint( DofManager const & dof elemManager.forElementSubRegions< FaceElementSubRegion >( [&]( FaceElementSubRegion & subRegion ) { - ContactBase const & contact = getConstitutiveModel< ContactBase >( subRegion, m_contactRelationName ); - - real64 const contactStiffness = contact.stiffness(); + real64 const contactStiffness = m_contactPenaltyStiffness; arrayView1d< real64 > const area = subRegion.getElementArea(); ArrayOfArraysView< localIndex const > const elemsToFaces = subRegion.faceList().toViewConst(); diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp index 9deae4f3272..9417db2282f 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp @@ -251,6 +251,9 @@ class SolidMechanicsLagrangianFEM : public SolverBase static constexpr char const * nonSendOrReceiveNodesString() { return "nonSendOrReceiveNodes";} static constexpr char const * targetNodesString() { return "targetNodes";} static constexpr char const * forceString() { return "Force";} + + static constexpr char const * contactPenaltyStiffnessString() { return "contactPenaltyStiffness"; } + }; SortedArray< localIndex > & getElemsAttachedToSendOrReceiveNodes( ElementSubRegionBase & subRegion ) @@ -297,6 +300,8 @@ class SolidMechanicsLagrangianFEM : public SolverBase /// Rigid body modes array1d< ParallelVector > m_rigidBodyModes; + real64 m_contactPenaltyStiffness; + private: string m_contactRelationName; diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsMPM.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsMPM.cpp index 08640b56332..16e8377c0f1 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsMPM.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsMPM.cpp @@ -29,7 +29,7 @@ #include "codingUtilities/Utilities.hpp" #include "common/TimingMacros.hpp" #include "constitutive/ConstitutiveManager.hpp" -#include "constitutive/contact/ContactBase.hpp" +#include "constitutive/contact/FrictionBase.hpp" #include "finiteElement/FiniteElementDiscretizationManager.hpp" #include "finiteElement/Kinematics.h" #include "LvArray/src/output.hpp" diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernels/ExplicitMPM.hpp b/src/coreComponents/physicsSolvers/solidMechanics/kernels/ExplicitMPM.hpp index 495872c8038..05f39528d30 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/kernels/ExplicitMPM.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/ExplicitMPM.hpp @@ -85,7 +85,7 @@ struct StateUpdateKernel /** * @brief Launch the kernel function doing constitutive updates * @tparam POLICY the type of policy used in the kernel launch - * @tparam CONTACT_WRAPPER the type of contact wrapper doing the constitutive updates + * @tparam CONSTITUTIVE_WRAPPER the type of consitutive wrapper doing the constitutive updates * @param[in] size the size of the subregion * @param[in] constitutiveWrapper the wrapper implementing the constitutive model * @param[in] deformationGradient the deformation gradient diff --git a/src/coreComponents/schema/docs/Constitutive.rst b/src/coreComponents/schema/docs/Constitutive.rst index 09b1aadeaff..c62cb08d281 100644 --- a/src/coreComponents/schema/docs/Constitutive.rst +++ b/src/coreComponents/schema/docs/Constitutive.rst @@ -42,6 +42,7 @@ ElasticTransverseIsotropic node :ref:`XML_ElasticTran ExponentialDecayPermeability node :ref:`XML_ExponentialDecayPermeability` ExtendedDruckerPrager node :ref:`XML_ExtendedDruckerPrager` FrictionlessContact node :ref:`XML_FrictionlessContact` +HydraulicApertureTable node :ref:`XML_HydraulicApertureTable` JFunctionCapillaryPressure node :ref:`XML_JFunctionCapillaryPressure` LinearIsotropicDispersion node :ref:`XML_LinearIsotropicDispersion` ModifiedCamClay node :ref:`XML_ModifiedCamClay` diff --git a/src/coreComponents/schema/docs/Constitutive_other.rst b/src/coreComponents/schema/docs/Constitutive_other.rst index d95527f1483..1195cddf026 100644 --- a/src/coreComponents/schema/docs/Constitutive_other.rst +++ b/src/coreComponents/schema/docs/Constitutive_other.rst @@ -42,6 +42,7 @@ ElasticTransverseIsotropic node :ref:`DATASTRUCTURE_ElasticTr ExponentialDecayPermeability node :ref:`DATASTRUCTURE_ExponentialDecayPermeability` ExtendedDruckerPrager node :ref:`DATASTRUCTURE_ExtendedDruckerPrager` FrictionlessContact node :ref:`DATASTRUCTURE_FrictionlessContact` +HydraulicApertureTable node :ref:`DATASTRUCTURE_HydraulicApertureTable` JFunctionCapillaryPressure node :ref:`DATASTRUCTURE_JFunctionCapillaryPressure` LinearIsotropicDispersion node :ref:`DATASTRUCTURE_LinearIsotropicDispersion` ModifiedCamClay node :ref:`DATASTRUCTURE_ModifiedCamClay` diff --git a/src/coreComponents/schema/docs/Coulomb.rst b/src/coreComponents/schema/docs/Coulomb.rst index 7b1a581e860..2cbf7fdeffc 100644 --- a/src/coreComponents/schema/docs/Coulomb.rst +++ b/src/coreComponents/schema/docs/Coulomb.rst @@ -1,16 +1,13 @@ -========================= ============ =========== ============================================================================================================================================================================================================================================================================================================================================================================================================================================================= -Name Type Default Description -========================= ============ =========== ============================================================================================================================================================================================================================================================================================================================================================================================================================================================= -apertureTableName groupNameRef required Name of the aperture table -apertureTolerance real64 1e-09 Value to be used to avoid floating point errors in expressions involving aperture. For example in the case of dividing by the actual aperture (not the effective aperture that results from the aperture function) this value may be used to avoid the 1/0 error. Note that this value may have some physical significance in its usage, as it may be used to smooth out highly nonlinear behavior associated with 1/0 in addition to avoiding the 1/0 error. -cohesion real64 required Cohesion -displacementJumpThreshold real64 2.22045e-16 A threshold valued to determine whether a fracture is open or not. -frictionCoefficient real64 required Friction coefficient -name groupName required A name is required for any non-unique nodes -penaltyStiffness real64 0 Value of the penetration penalty stiffness. Units of Pressure/length -shearStiffness real64 0 Value of the shear elastic stiffness. Units of Pressure/length -========================= ============ =========== ============================================================================================================================================================================================================================================================================================================================================================================================================================================================= +========================= ========= =========== ================================================================== +Name Type Default Description +========================= ========= =========== ================================================================== +cohesion real64 required Cohesion +displacementJumpThreshold real64 2.22045e-16 A threshold valued to determine whether a fracture is open or not. +frictionCoefficient real64 required Friction coefficient +name groupName required A name is required for any non-unique nodes +shearStiffness real64 0 Value of the shear elastic stiffness. Units of Pressure/length +========================= ========= =========== ================================================================== diff --git a/src/coreComponents/schema/docs/FrictionlessContact.rst b/src/coreComponents/schema/docs/FrictionlessContact.rst index 385d3738b2a..e150b24c06a 100644 --- a/src/coreComponents/schema/docs/FrictionlessContact.rst +++ b/src/coreComponents/schema/docs/FrictionlessContact.rst @@ -1,14 +1,10 @@ -========================= ============ =========== ============================================================================================================================================================================================================================================================================================================================================================================================================================================================= -Name Type Default Description -========================= ============ =========== ============================================================================================================================================================================================================================================================================================================================================================================================================================================================= -apertureTableName groupNameRef required Name of the aperture table -apertureTolerance real64 1e-09 Value to be used to avoid floating point errors in expressions involving aperture. For example in the case of dividing by the actual aperture (not the effective aperture that results from the aperture function) this value may be used to avoid the 1/0 error. Note that this value may have some physical significance in its usage, as it may be used to smooth out highly nonlinear behavior associated with 1/0 in addition to avoiding the 1/0 error. -displacementJumpThreshold real64 2.22045e-16 A threshold valued to determine whether a fracture is open or not. -name groupName required A name is required for any non-unique nodes -penaltyStiffness real64 0 Value of the penetration penalty stiffness. Units of Pressure/length -shearStiffness real64 0 Value of the shear elastic stiffness. Units of Pressure/length -========================= ============ =========== ============================================================================================================================================================================================================================================================================================================================================================================================================================================================= +========================= ========= =========== ================================================================== +Name Type Default Description +========================= ========= =========== ================================================================== +displacementJumpThreshold real64 2.22045e-16 A threshold valued to determine whether a fracture is open or not. +name groupName required A name is required for any non-unique nodes +========================= ========= =========== ================================================================== diff --git a/src/coreComponents/schema/docs/Hydrofracture.rst b/src/coreComponents/schema/docs/Hydrofracture.rst index 896a5cb77b0..d36978900e6 100644 --- a/src/coreComponents/schema/docs/Hydrofracture.rst +++ b/src/coreComponents/schema/docs/Hydrofracture.rst @@ -4,7 +4,6 @@ Name Type Default Description ===================================== =================================================================================================================================== ======== ====================================================================================================================================================================================================================================================================================================================== cflFactor real64 0.5 Factor to apply to the `CFL condition `_ when calculating the maximum allowable time step. Values should be in the interval (0,1] -contactRelationName groupNameRef required Name of contact relation to enforce constraints on fracture boundary. flowSolverName groupNameRef required Name of the flow solver used by the coupled solver initialDt real64 1e+99 Initial time-step value required by the solver to the event manager. isLaggingFractureStencilWeightsUpdate integer 0 Flag to determine whether or not to apply lagging update for the fracture stencil weights. diff --git a/src/coreComponents/schema/docs/SolidMechanicsAugmentedLagrangianContact.rst b/src/coreComponents/schema/docs/SolidMechanicsAugmentedLagrangianContact.rst index 9b736b1b4d6..4ea2de12d3a 100644 --- a/src/coreComponents/schema/docs/SolidMechanicsAugmentedLagrangianContact.rst +++ b/src/coreComponents/schema/docs/SolidMechanicsAugmentedLagrangianContact.rst @@ -4,6 +4,7 @@ Name Type Default Description ========================= ====================================================== =============== ======================================================================================================================================================================================================================================================================================================================== cflFactor real64 0.5 Factor to apply to the `CFL condition `_ when calculating the maximum allowable time step. Values should be in the interval (0,1] +contactPenaltyStiffness real64 0 Value of the penetration penalty stiffness. Units of Pressure/length discretization groupNameRef required Name of discretization object (defined in the :ref:`NumericalMethodsManager`) to use for this solver. For instance, if this is a Finite Element Solver, the name of a :ref:`FiniteElement` should be specified. If this is a Finite Volume Method, the name of a :ref:`FiniteVolume` discretization should be specified. initialDt real64 1e+99 Initial time-step value required by the solver to the event manager. logLevel integer 0 Log level diff --git a/src/coreComponents/schema/docs/SolidMechanicsEmbeddedFractures.rst b/src/coreComponents/schema/docs/SolidMechanicsEmbeddedFractures.rst index eca5d02f2ed..f1eeeefdbf1 100644 --- a/src/coreComponents/schema/docs/SolidMechanicsEmbeddedFractures.rst +++ b/src/coreComponents/schema/docs/SolidMechanicsEmbeddedFractures.rst @@ -4,6 +4,7 @@ Name Type Default Description ========================= ====================================================== =============== ======================================================================================================================================================================================================================================================================================================================== cflFactor real64 0.5 Factor to apply to the `CFL condition `_ when calculating the maximum allowable time step. Values should be in the interval (0,1] +contactPenaltyStiffness real64 required Value of the penetration penalty stiffness. Units of Pressure/length discretization groupNameRef required Name of discretization object (defined in the :ref:`NumericalMethodsManager`) to use for this solver. For instance, if this is a Finite Element Solver, the name of a :ref:`FiniteElement` should be specified. If this is a Finite Volume Method, the name of a :ref:`FiniteVolume` discretization should be specified. initialDt real64 1e+99 Initial time-step value required by the solver to the event manager. logLevel integer 0 Log level diff --git a/src/coreComponents/schema/docs/SolidMechanicsLagrangeContact.rst b/src/coreComponents/schema/docs/SolidMechanicsLagrangeContact.rst index ea949482327..e801acd20e0 100644 --- a/src/coreComponents/schema/docs/SolidMechanicsLagrangeContact.rst +++ b/src/coreComponents/schema/docs/SolidMechanicsLagrangeContact.rst @@ -4,6 +4,7 @@ Name Type Default Description =============================== ====================================================== =============== ======================================================================================================================================================================================================================================================================================================================== cflFactor real64 0.5 Factor to apply to the `CFL condition `_ when calculating the maximum allowable time step. Values should be in the interval (0,1] +contactPenaltyStiffness real64 0 Value of the penetration penalty stiffness. Units of Pressure/length discretization groupNameRef required Name of discretization object (defined in the :ref:`NumericalMethodsManager`) to use for this solver. For instance, if this is a Finite Element Solver, the name of a :ref:`FiniteElement` should be specified. If this is a Finite Volume Method, the name of a :ref:`FiniteVolume` discretization should be specified. initialDt real64 1e+99 Initial time-step value required by the solver to the event manager. logLevel integer 0 Log level diff --git a/src/coreComponents/schema/docs/SolidMechanicsLagrangianSSLE.rst b/src/coreComponents/schema/docs/SolidMechanicsLagrangianSSLE.rst index bed0562d496..1bccbf3c052 100644 --- a/src/coreComponents/schema/docs/SolidMechanicsLagrangianSSLE.rst +++ b/src/coreComponents/schema/docs/SolidMechanicsLagrangianSSLE.rst @@ -4,6 +4,7 @@ Name Type Default Description ========================= ====================================================== =============== ======================================================================================================================================================================================================================================================================================================================== cflFactor real64 0.5 Factor to apply to the `CFL condition `_ when calculating the maximum allowable time step. Values should be in the interval (0,1] +contactPenaltyStiffness real64 0 Value of the penetration penalty stiffness. Units of Pressure/length contactRelationName groupNameRef NOCONTACT Name of contact relation to enforce constraints on fracture boundary. discretization groupNameRef required Name of discretization object (defined in the :ref:`NumericalMethodsManager`) to use for this solver. For instance, if this is a Finite Element Solver, the name of a :ref:`FiniteElement` should be specified. If this is a Finite Volume Method, the name of a :ref:`FiniteVolume` discretization should be specified. initialDt real64 1e+99 Initial time-step value required by the solver to the event manager. diff --git a/src/coreComponents/schema/docs/SolidMechanics_LagrangianFEM.rst b/src/coreComponents/schema/docs/SolidMechanics_LagrangianFEM.rst index bed0562d496..1bccbf3c052 100644 --- a/src/coreComponents/schema/docs/SolidMechanics_LagrangianFEM.rst +++ b/src/coreComponents/schema/docs/SolidMechanics_LagrangianFEM.rst @@ -4,6 +4,7 @@ Name Type Default Description ========================= ====================================================== =============== ======================================================================================================================================================================================================================================================================================================================== cflFactor real64 0.5 Factor to apply to the `CFL condition `_ when calculating the maximum allowable time step. Values should be in the interval (0,1] +contactPenaltyStiffness real64 0 Value of the penetration penalty stiffness. Units of Pressure/length contactRelationName groupNameRef NOCONTACT Name of contact relation to enforce constraints on fracture boundary. discretization groupNameRef required Name of discretization object (defined in the :ref:`NumericalMethodsManager`) to use for this solver. For instance, if this is a Finite Element Solver, the name of a :ref:`FiniteElement` should be specified. If this is a Finite Volume Method, the name of a :ref:`FiniteVolume` discretization should be specified. initialDt real64 1e+99 Initial time-step value required by the solver to the event manager. diff --git a/src/coreComponents/schema/schema.xsd b/src/coreComponents/schema/schema.xsd index d02adab186c..2fe20af996e 100644 --- a/src/coreComponents/schema/schema.xsd +++ b/src/coreComponents/schema/schema.xsd @@ -701,6 +701,10 @@ + + + + @@ -2919,8 +2923,6 @@ Equal to 1 for surface conditions, and to 0 for reservoir conditions--> - - @@ -3568,6 +3570,8 @@ Local- Add jump stabilization on interior of macro elements--> + + @@ -3612,6 +3616,8 @@ Local- Add jump stabilization on interior of macro elements--> + + @@ -3653,6 +3659,8 @@ Local- Add jump stabilization on interior of macro elements--> + + @@ -3696,6 +3704,8 @@ Local- Add jump stabilization on interior of macro elements--> + + @@ -3739,6 +3749,8 @@ Local- Add jump stabilization on interior of macro elements--> + + @@ -4180,6 +4192,7 @@ Local- Add jump stabilization on interior of macro elements--> + @@ -4671,18 +4684,12 @@ The expected format is "{ waterMax, oilMax }", in that order--> - - - - - - @@ -5005,16 +5012,18 @@ For instance, if "oil" is before "gas" in "phaseNames", the table order should b + + + + + + - - - - - - + + diff --git a/src/coreComponents/schema/schema.xsd.other b/src/coreComponents/schema/schema.xsd.other index 73552a29993..fac7e939187 100644 --- a/src/coreComponents/schema/schema.xsd.other +++ b/src/coreComponents/schema/schema.xsd.other @@ -1449,6 +1449,7 @@ + @@ -2335,6 +2336,7 @@ + diff --git a/src/docs/sphinx/CompleteXMLSchema.rst b/src/docs/sphinx/CompleteXMLSchema.rst index 93fd92c0af4..92595791056 100644 --- a/src/docs/sphinx/CompleteXMLSchema.rst +++ b/src/docs/sphinx/CompleteXMLSchema.rst @@ -563,6 +563,13 @@ Element: HybridMimeticDiscretization .. include:: ../../coreComponents/schema/docs/HybridMimeticDiscretization.rst +.. _XML_HydraulicApertureTable: + +Element: HydraulicApertureTable +=============================== +.. include:: ../../coreComponents/schema/docs/HydraulicApertureTable.rst + + .. _XML_Hydrofracture: Element: Hydrofracture @@ -2037,6 +2044,13 @@ Datastructure: HybridMimeticDiscretization .. include:: ../../coreComponents/schema/docs/HybridMimeticDiscretization_other.rst +.. _DATASTRUCTURE_HydraulicApertureTable: + +Datastructure: HydraulicApertureTable +===================================== +.. include:: ../../coreComponents/schema/docs/HydraulicApertureTable_other.rst + + .. _DATASTRUCTURE_Hydrofracture: Datastructure: Hydrofracture diff --git a/src/docs/sphinx/developerGuide/KeyComponents/AddingNewSolver.rst b/src/docs/sphinx/developerGuide/KeyComponents/AddingNewSolver.rst index 00bf5cd0a90..38c122242b0 100644 --- a/src/docs/sphinx/developerGuide/KeyComponents/AddingNewSolver.rst +++ b/src/docs/sphinx/developerGuide/KeyComponents/AddingNewSolver.rst @@ -117,7 +117,7 @@ when ``applyBoundaryConditions()`` is called in this particular class override. Browsing the base class ``SolverBase``, it can be noted that most of the solver interface functions are called during either ``SolverBase::linearImplicitStep()`` or ``SolverBase::nonlinearImplicitStep()`` depending on the solver strategy chosen. -Switching to protected members, ``postProcessInput()`` is a central member function and +Switching to protected members, ``postInputInitialization()`` is a central member function and will be called by ``Group`` object after input is read from XML entry file. It will set and dispatch solver variables from the base class ``SolverBase`` to the most derived class. For ``LaplaceFEM``, it will allow us to set the right time integration scheme based on the XML value @@ -202,7 +202,7 @@ to writing our new *LaplaceDiffFEM* solver. .. note:: - We might want to remove final keyword from ``postProcessInput()`` as it will prevent you from overriding it. + We might want to remove final keyword from ``postInputInitialization()`` as it will prevent you from overriding it. Start doing your own Physic solver ================================== @@ -249,7 +249,7 @@ commented afterwards. } laplaceDiffFEMViewKeys; protected: - virtual void postProcessInput() override final; + virtual void postInputInitialization() override final; private: real64 m_diffusion; @@ -268,7 +268,7 @@ Then as mentioned in :ref:`Implementation`, the diffusion coefficient is used wh we will have to override the ``assembleSystem()`` function as detailed below. Moreover, if we want to introduce a new binding between the input XML and the code we will have to work on the three -``struct viewKeyStruct`` , ``postProcessInput()`` and the constructor. +``struct viewKeyStruct`` , ``postInputInitialization()`` and the constructor. Our new solver ``viewKeyStruct`` will have its own structure inheriting from the *LaplaceFEM* one to have the ``timeIntegrationOption`` and ``fieldName`` field. It will also create a ``diffusionCoeff`` field to be bound to the user defined homogeneous coefficient on one hand @@ -293,13 +293,13 @@ an "input uniform diffusion coefficient for the Laplace equation". setDescription("input uniform diffusion coeff for the laplace equation"); } -Another important spot for binding the value of the XML read parameter to our ``m_diffusion`` is in ``postProcessInput()``. +Another important spot for binding the value of the XML read parameter to our ``m_diffusion`` is in ``postInputInitialization()``. .. code-block:: c++ - void LaplaceDiffFEM::postProcessInput() + void LaplaceDiffFEM::postInputInitialization() { - LaplaceFEM::postProcessInput(); + LaplaceFEM::postInputInitialization(); string sDiffCoeff = this->getReference(laplaceDiffFEMViewKeys.diffusionCoeff); this->m_diffusion = std::stof(sDiffCoeff);