From 1e6446b4a71a8526b2deaea641fba74d02a7bf8c Mon Sep 17 00:00:00 2001 From: Matteo Frigo <43966286+matteofrigo5@users.noreply.github.com> Date: Tue, 17 Sep 2024 10:53:18 -0700 Subject: [PATCH] feat: Augmented Lagrangian (slip and open modes) (#3217) * Implementation of slip and open modes for ALM * Adding a new inputFile for ALM and bug (bubble functions) fixed * Bug bubble functions gradient fixed * Adding initial stress contribution for bubble functions - updating lambda functions * Adding nested ALM * Adding simultaneous and nested ALM * Moving the traction update into the friction model * Moved contact constitutive behavior from ALM to the friction law --- .integrated_tests.yaml | 2 +- BASELINE_NOTES.md | 4 + ...rlock-gcc10-ompi4.1.2-openblas0.3.10.cmake | 1 - .../ALM_stickFault_base.xml | 115 ---- inputFiles/almContactMechanics/alm.ats | 19 - .../ALM_PassingCrack_smoke.xml | 55 ++ .../ALM_SimpleCubes_smoke.xml | 58 ++ .../ALM_SingleFracCompression_benchmark.xml | 77 +++ .../ALM_SingleFracCompression_smoke.xml | 77 +++ .../ALM_Sneddon_benchmark.xml | 67 ++ .../ALM_Sneddon_smoke.xml | 67 ++ .../ALM_TFrac_benchmark.xml | 77 +++ .../ALM_TFrac_smoke.xml | 77 +++ .../ALM_UnstructuredCrack_benchmark.xml | 57 ++ .../ALM_UnstructuredCrack_smoke.xml | 58 ++ .../ALM_slippingFault_horizontal_smoke.xml | 55 ++ .../ALM_slippingFault_vertical_smoke.xml | 55 ++ .../ContactMechanics_PassingCrack_smoke.xml | 95 +-- .../ContactMechanics_SimpleCubes_smoke.xml | 76 +-- ...hanics_SingleFracCompression_benchmark.xml | 106 ++- ...tMechanics_SingleFracCompression_smoke.xml | 100 +-- ...=> ContactMechanics_Sneddon_benchmark.xml} | 30 +- ...xml => ContactMechanics_Sneddon_smoke.xml} | 57 +- .../ContactMechanics_TFrac_benchmark.xml | 76 +-- .../ContactMechanics_TFrac_smoke.xml | 33 +- ...tMechanics_UnstructuredCrack_benchmark.xml | 95 +-- ...ntactMechanics_UnstructuredCrack_smoke.xml | 95 +-- .../ContactMechanics_slippingFault_base.xml | 111 --- ...chanics_slippingFault_horizontal_smoke.xml | 174 ++--- ...Mechanics_slippingFault_vertical_smoke.xml | 147 ++-- ...ngCrack_base.xml => PassingCrack_base.xml} | 27 +- .../PassingCrack_smoke.xml | 58 ++ ...pleCubes_base.xml => SimpleCubes_base.xml} | 24 - .../SimpleCubes_smoke.xml | 54 ++ ...ase.xml => SingleFracCompression_base.xml} | 26 - .../SingleFracCompression_benchmark.xml | 75 ++ .../SingleFracCompression_smoke.xml | 75 ++ .../SlippingFault_base.xml | 59 ++ .../SlippingFault_horizontal_smoke.xml | 134 ++++ .../SlippingFault_vertical_smoke.xml} | 55 +- ...actMechanics_base.xml => Sneddon_base.xml} | 7 +- .../Sneddon_benchmark.xml | 34 + .../Sneddon_smoke.xml | 35 + ...echanics_TFrac_base.xml => TFrac_base.xml} | 6 +- .../TFrac_benchmark.xml | 36 + .../TFrac_smoke.xml | 42 ++ ...ck_base.xml => UnstructuredCrack_base.xml} | 24 +- .../UnstructuredCrack_benchmark.xml | 74 ++ .../UnstructuredCrack_smoke.xml | 74 ++ .../contactMechanics.ats | 41 +- .../constitutive/ConstitutivePassThru.hpp | 15 + .../constitutive/contact/CoulombFriction.cpp | 2 +- .../constitutive/contact/CoulombFriction.hpp | 283 +++++++- .../constitutive/contact/FrictionBase.hpp | 74 ++ .../contact/FrictionlessContact.cpp | 2 +- .../contact/FrictionlessContact.hpp | 2 +- .../elementFormulations/LagrangeBasis1.hpp | 127 ++-- .../physicsSolvers/CMakeLists.txt | 3 +- .../physicsSolvers/contact/ContactFields.hpp | 8 +- .../SolidMechanicsALMBubbleKernels.hpp | 41 +- .../contact/SolidMechanicsALMKernels.hpp | 90 +-- .../contact/SolidMechanicsALMKernelsBase.hpp | 178 ++++- .../SolidMechanicsALMSimultaneousKernels.hpp | 435 ++++++++++++ ...lidMechanicsAugmentedLagrangianContact.cpp | 643 ++++++++++++++---- ...lidMechanicsAugmentedLagrangianContact.hpp | 86 ++- .../SolidMechanicsEmbeddedFractures.cpp | 4 +- .../contact/SolidMechanicsLagrangeContact.cpp | 4 +- .../faultMechanics/intersectFrac/Example.rst | 18 +- .../intersectFrac/intersectFracFigure.py | 2 +- .../singleFracCompression/Example.rst | 18 +- .../singleFracCompressionFigure.py | 4 +- .../faultMechanics/sneddon/Example.rst | 15 +- 72 files changed, 3740 insertions(+), 1390 deletions(-) delete mode 100644 inputFiles/almContactMechanics/ALM_stickFault_base.xml delete mode 100644 inputFiles/almContactMechanics/alm.ats create mode 100644 inputFiles/lagrangianContactMechanics/ALM_PassingCrack_smoke.xml create mode 100644 inputFiles/lagrangianContactMechanics/ALM_SimpleCubes_smoke.xml create mode 100644 inputFiles/lagrangianContactMechanics/ALM_SingleFracCompression_benchmark.xml create mode 100644 inputFiles/lagrangianContactMechanics/ALM_SingleFracCompression_smoke.xml create mode 100644 inputFiles/lagrangianContactMechanics/ALM_Sneddon_benchmark.xml create mode 100644 inputFiles/lagrangianContactMechanics/ALM_Sneddon_smoke.xml create mode 100644 inputFiles/lagrangianContactMechanics/ALM_TFrac_benchmark.xml create mode 100644 inputFiles/lagrangianContactMechanics/ALM_TFrac_smoke.xml create mode 100644 inputFiles/lagrangianContactMechanics/ALM_UnstructuredCrack_benchmark.xml create mode 100644 inputFiles/lagrangianContactMechanics/ALM_UnstructuredCrack_smoke.xml create mode 100644 inputFiles/lagrangianContactMechanics/ALM_slippingFault_horizontal_smoke.xml create mode 100644 inputFiles/lagrangianContactMechanics/ALM_slippingFault_vertical_smoke.xml rename inputFiles/lagrangianContactMechanics/{Sneddon_contactMechanics_smoke.xml => ContactMechanics_Sneddon_benchmark.xml} (75%) mode change 100755 => 100644 rename inputFiles/lagrangianContactMechanics/{Sneddon_contactMechanics_benchmark.xml => ContactMechanics_Sneddon_smoke.xml} (58%) mode change 100755 => 100644 delete mode 100644 inputFiles/lagrangianContactMechanics/ContactMechanics_slippingFault_base.xml rename inputFiles/lagrangianContactMechanics/{ContactMechanics_PassingCrack_base.xml => PassingCrack_base.xml} (78%) create mode 100644 inputFiles/lagrangianContactMechanics/PassingCrack_smoke.xml rename inputFiles/lagrangianContactMechanics/{ContactMechanics_SimpleCubes_base.xml => SimpleCubes_base.xml} (81%) create mode 100644 inputFiles/lagrangianContactMechanics/SimpleCubes_smoke.xml rename inputFiles/lagrangianContactMechanics/{ContactMechanics_SingleFracCompression_base.xml => SingleFracCompression_base.xml} (81%) mode change 100755 => 100644 create mode 100644 inputFiles/lagrangianContactMechanics/SingleFracCompression_benchmark.xml create mode 100644 inputFiles/lagrangianContactMechanics/SingleFracCompression_smoke.xml create mode 100644 inputFiles/lagrangianContactMechanics/SlippingFault_base.xml create mode 100644 inputFiles/lagrangianContactMechanics/SlippingFault_horizontal_smoke.xml rename inputFiles/{almContactMechanics/ALM_stickFault_vertical_smoke.xml => lagrangianContactMechanics/SlippingFault_vertical_smoke.xml} (73%) rename inputFiles/lagrangianContactMechanics/{Sneddon_contactMechanics_base.xml => Sneddon_base.xml} (96%) mode change 100755 => 100644 create mode 100644 inputFiles/lagrangianContactMechanics/Sneddon_benchmark.xml create mode 100644 inputFiles/lagrangianContactMechanics/Sneddon_smoke.xml rename inputFiles/lagrangianContactMechanics/{ContactMechanics_TFrac_base.xml => TFrac_base.xml} (97%) create mode 100644 inputFiles/lagrangianContactMechanics/TFrac_benchmark.xml create mode 100644 inputFiles/lagrangianContactMechanics/TFrac_smoke.xml rename inputFiles/lagrangianContactMechanics/{ContactMechanics_UnstructuredCrack_base.xml => UnstructuredCrack_base.xml} (80%) create mode 100644 inputFiles/lagrangianContactMechanics/UnstructuredCrack_benchmark.xml create mode 100644 inputFiles/lagrangianContactMechanics/UnstructuredCrack_smoke.xml create mode 100644 src/coreComponents/physicsSolvers/contact/SolidMechanicsALMSimultaneousKernels.hpp diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml index 85a23286bf6..53c078f12d0 100644 --- a/.integrated_tests.yaml +++ b/.integrated_tests.yaml @@ -1,6 +1,6 @@ baselines: bucket: geosx - baseline: integratedTests/baseline_integratedTests-pr3318-7578-14f13e2 + baseline: integratedTests/baseline_integratedTests-pr3217-7661-9b4a604 allow_fail: all: '' diff --git a/BASELINE_NOTES.md b/BASELINE_NOTES.md index 7c26b2ef528..a01e4749868 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 #3217 (2024-09-16) +====================== +ALM slip and open modes with relative tests. + PR #3318 (2024-09-12) ====================== Modified SeismicityRate poroelastic case. diff --git a/host-configs/Stanford/sherlock-gcc10-ompi4.1.2-openblas0.3.10.cmake b/host-configs/Stanford/sherlock-gcc10-ompi4.1.2-openblas0.3.10.cmake index 86013187457..d29b6180e11 100644 --- a/host-configs/Stanford/sherlock-gcc10-ompi4.1.2-openblas0.3.10.cmake +++ b/host-configs/Stanford/sherlock-gcc10-ompi4.1.2-openblas0.3.10.cmake @@ -40,5 +40,4 @@ set(LAPACK_LIBRARIES "${OPENBLAS_ROOT}/lib/libopenblas.so" CACHE STRING "") set(ENABLE_VALGRIND OFF CACHE BOOL "") set(ENABLE_CALIPER ON CACHE BOOL "") -set(GEOS_TPL_DIR "$ENV{GEOSX_TPL_DIR}" CACHE PATH "" FORCE) include(${CMAKE_CURRENT_LIST_DIR}/../tpls.cmake) diff --git a/inputFiles/almContactMechanics/ALM_stickFault_base.xml b/inputFiles/almContactMechanics/ALM_stickFault_base.xml deleted file mode 100644 index 7a8c79aa7d8..00000000000 --- a/inputFiles/almContactMechanics/ALM_stickFault_base.xml +++ /dev/null @@ -1,115 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/inputFiles/almContactMechanics/alm.ats b/inputFiles/almContactMechanics/alm.ats deleted file mode 100644 index dc6417df533..00000000000 --- a/inputFiles/almContactMechanics/alm.ats +++ /dev/null @@ -1,19 +0,0 @@ -import os -from geos.ats.test_builder import TestDeck, RestartcheckParameters, generate_geos_tests - -restartcheck_params = {} -restartcheck_params["atol"] = 2.0E-4 -restartcheck_params["rtol"] = 1.0E-7 - -decks = [ - TestDeck( - name="ALM_stickFault_vertical_smoke", - description= - "Cube with a vertical fracture (structured grid)", - partitions=((1, 1, 1), (2, 2, 2)), - restart_step=1, - check_step=2, - restartcheck_params=RestartcheckParameters(**restartcheck_params)), -] - -generate_geos_tests(decks) diff --git a/inputFiles/lagrangianContactMechanics/ALM_PassingCrack_smoke.xml b/inputFiles/lagrangianContactMechanics/ALM_PassingCrack_smoke.xml new file mode 100644 index 00000000000..2025d1e4e19 --- /dev/null +++ b/inputFiles/lagrangianContactMechanics/ALM_PassingCrack_smoke.xml @@ -0,0 +1,55 @@ + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/lagrangianContactMechanics/ALM_SimpleCubes_smoke.xml b/inputFiles/lagrangianContactMechanics/ALM_SimpleCubes_smoke.xml new file mode 100644 index 00000000000..d0c48b78bdf --- /dev/null +++ b/inputFiles/lagrangianContactMechanics/ALM_SimpleCubes_smoke.xml @@ -0,0 +1,58 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/lagrangianContactMechanics/ALM_SingleFracCompression_benchmark.xml b/inputFiles/lagrangianContactMechanics/ALM_SingleFracCompression_benchmark.xml new file mode 100644 index 00000000000..2202aa2f739 --- /dev/null +++ b/inputFiles/lagrangianContactMechanics/ALM_SingleFracCompression_benchmark.xml @@ -0,0 +1,77 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/lagrangianContactMechanics/ALM_SingleFracCompression_smoke.xml b/inputFiles/lagrangianContactMechanics/ALM_SingleFracCompression_smoke.xml new file mode 100644 index 00000000000..8b09acbea67 --- /dev/null +++ b/inputFiles/lagrangianContactMechanics/ALM_SingleFracCompression_smoke.xml @@ -0,0 +1,77 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/lagrangianContactMechanics/ALM_Sneddon_benchmark.xml b/inputFiles/lagrangianContactMechanics/ALM_Sneddon_benchmark.xml new file mode 100644 index 00000000000..ba0462cbe05 --- /dev/null +++ b/inputFiles/lagrangianContactMechanics/ALM_Sneddon_benchmark.xml @@ -0,0 +1,67 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/lagrangianContactMechanics/ALM_Sneddon_smoke.xml b/inputFiles/lagrangianContactMechanics/ALM_Sneddon_smoke.xml new file mode 100644 index 00000000000..41c125ef0e5 --- /dev/null +++ b/inputFiles/lagrangianContactMechanics/ALM_Sneddon_smoke.xml @@ -0,0 +1,67 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/lagrangianContactMechanics/ALM_TFrac_benchmark.xml b/inputFiles/lagrangianContactMechanics/ALM_TFrac_benchmark.xml new file mode 100644 index 00000000000..8acc35c09b1 --- /dev/null +++ b/inputFiles/lagrangianContactMechanics/ALM_TFrac_benchmark.xml @@ -0,0 +1,77 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/lagrangianContactMechanics/ALM_TFrac_smoke.xml b/inputFiles/lagrangianContactMechanics/ALM_TFrac_smoke.xml new file mode 100644 index 00000000000..c7a5bf66543 --- /dev/null +++ b/inputFiles/lagrangianContactMechanics/ALM_TFrac_smoke.xml @@ -0,0 +1,77 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/lagrangianContactMechanics/ALM_UnstructuredCrack_benchmark.xml b/inputFiles/lagrangianContactMechanics/ALM_UnstructuredCrack_benchmark.xml new file mode 100644 index 00000000000..4a893d5dc0a --- /dev/null +++ b/inputFiles/lagrangianContactMechanics/ALM_UnstructuredCrack_benchmark.xml @@ -0,0 +1,57 @@ + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/lagrangianContactMechanics/ALM_UnstructuredCrack_smoke.xml b/inputFiles/lagrangianContactMechanics/ALM_UnstructuredCrack_smoke.xml new file mode 100644 index 00000000000..4b230c006ce --- /dev/null +++ b/inputFiles/lagrangianContactMechanics/ALM_UnstructuredCrack_smoke.xml @@ -0,0 +1,58 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/lagrangianContactMechanics/ALM_slippingFault_horizontal_smoke.xml b/inputFiles/lagrangianContactMechanics/ALM_slippingFault_horizontal_smoke.xml new file mode 100644 index 00000000000..4060a64af2e --- /dev/null +++ b/inputFiles/lagrangianContactMechanics/ALM_slippingFault_horizontal_smoke.xml @@ -0,0 +1,55 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/lagrangianContactMechanics/ALM_slippingFault_vertical_smoke.xml b/inputFiles/lagrangianContactMechanics/ALM_slippingFault_vertical_smoke.xml new file mode 100644 index 00000000000..7beefc7386d --- /dev/null +++ b/inputFiles/lagrangianContactMechanics/ALM_slippingFault_vertical_smoke.xml @@ -0,0 +1,55 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/lagrangianContactMechanics/ContactMechanics_PassingCrack_smoke.xml b/inputFiles/lagrangianContactMechanics/ContactMechanics_PassingCrack_smoke.xml index 4ba2c33f6d8..cccc821fc54 100644 --- a/inputFiles/lagrangianContactMechanics/ContactMechanics_PassingCrack_smoke.xml +++ b/inputFiles/lagrangianContactMechanics/ContactMechanics_PassingCrack_smoke.xml @@ -3,57 +3,35 @@ + name="./PassingCrack_smoke.xml"/> - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + @@ -61,26 +39,23 @@ name="preFracture" target="/Solvers/SurfaceGen"/> - - - + target="/Outputs/vtkOutput"/> + + diff --git a/inputFiles/lagrangianContactMechanics/ContactMechanics_SimpleCubes_smoke.xml b/inputFiles/lagrangianContactMechanics/ContactMechanics_SimpleCubes_smoke.xml index 74ef1533500..733a439db2e 100644 --- a/inputFiles/lagrangianContactMechanics/ContactMechanics_SimpleCubes_smoke.xml +++ b/inputFiles/lagrangianContactMechanics/ContactMechanics_SimpleCubes_smoke.xml @@ -3,54 +3,37 @@ + name="./SimpleCubes_smoke.xml"/> + + + + + + - - - - - - - - - - - - - - - - - + + + + + @@ -74,9 +57,6 @@ targetExactTimestep="0" target="/Outputs/restartOutput"/> - - + - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + @@ -81,7 +48,6 @@ @@ -90,6 +56,12 @@ timeFrequency="1" targetExactTimestep="0" target="/Outputs/vtkOutput"/> + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + diff --git a/inputFiles/lagrangianContactMechanics/Sneddon_contactMechanics_smoke.xml b/inputFiles/lagrangianContactMechanics/ContactMechanics_Sneddon_benchmark.xml old mode 100755 new mode 100644 similarity index 75% rename from inputFiles/lagrangianContactMechanics/Sneddon_contactMechanics_smoke.xml rename to inputFiles/lagrangianContactMechanics/ContactMechanics_Sneddon_benchmark.xml index e7962f070c1..cefdf59664b --- a/inputFiles/lagrangianContactMechanics/Sneddon_contactMechanics_smoke.xml +++ b/inputFiles/lagrangianContactMechanics/ContactMechanics_Sneddon_benchmark.xml @@ -3,7 +3,7 @@ + name="./Sneddon_benchmark.xml"/> - - - - - - - + + + + + + name="./Sneddon_smoke.xml"/> - + - - - + - - - - - + + + + + + + + diff --git a/inputFiles/lagrangianContactMechanics/ContactMechanics_TFrac_benchmark.xml b/inputFiles/lagrangianContactMechanics/ContactMechanics_TFrac_benchmark.xml index 83e9561f2c4..89e2f4854df 100644 --- a/inputFiles/lagrangianContactMechanics/ContactMechanics_TFrac_benchmark.xml +++ b/inputFiles/lagrangianContactMechanics/ContactMechanics_TFrac_benchmark.xml @@ -2,60 +2,44 @@ - + - - + + targetRegions="{ Region, Fracture }"> + newtonTol="1.0e-8" + logLevel="2" + maxNumConfigurationAttempts="10" + newtonMaxIter="10" + lineSearchAction="Require" + lineSearchMaxCuts="2" + maxTimeStepCuts="2"/> - + solverType="direct" + directParallel="0" + logLevel="0"/> + + + - - - + + + + + - - - - - - + maxTime="0.2"> @@ -63,7 +47,7 @@ @@ -93,6 +77,12 @@ name="displacementHistoryOutput" timeFrequency="0.2" targetExactTimestep="0" - target="/Outputs/displacementOutput"/> + target="/Outputs/displacementOutput"/> + + diff --git a/inputFiles/lagrangianContactMechanics/ContactMechanics_TFrac_smoke.xml b/inputFiles/lagrangianContactMechanics/ContactMechanics_TFrac_smoke.xml index b788e0ac183..d3a1a91319c 100644 --- a/inputFiles/lagrangianContactMechanics/ContactMechanics_TFrac_smoke.xml +++ b/inputFiles/lagrangianContactMechanics/ContactMechanics_TFrac_smoke.xml @@ -2,7 +2,7 @@ - + - - - - - - - - - + + + + + + - + name="./UnstructuredCrack_benchmark.xml"/> - - - - - - - - - - - - + + + + + + + + + + + + - - - - - - - - - - - - - + diff --git a/inputFiles/lagrangianContactMechanics/ContactMechanics_UnstructuredCrack_smoke.xml b/inputFiles/lagrangianContactMechanics/ContactMechanics_UnstructuredCrack_smoke.xml index 71b7da13653..ea9c2adfe64 100644 --- a/inputFiles/lagrangianContactMechanics/ContactMechanics_UnstructuredCrack_smoke.xml +++ b/inputFiles/lagrangianContactMechanics/ContactMechanics_UnstructuredCrack_smoke.xml @@ -1,75 +1,38 @@ - + name="./UnstructuredCrack_smoke.xml"/> - - - - - - - - - - - - - - + + + + + + + + + + + + - - - - - - - - - - @@ -90,14 +53,12 @@ targetExactTimestep="0" target="/Outputs/restartOutput"/> - + diff --git a/inputFiles/lagrangianContactMechanics/ContactMechanics_slippingFault_base.xml b/inputFiles/lagrangianContactMechanics/ContactMechanics_slippingFault_base.xml deleted file mode 100644 index 0c5350192c1..00000000000 --- a/inputFiles/lagrangianContactMechanics/ContactMechanics_slippingFault_base.xml +++ /dev/null @@ -1,111 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/inputFiles/lagrangianContactMechanics/ContactMechanics_slippingFault_horizontal_smoke.xml b/inputFiles/lagrangianContactMechanics/ContactMechanics_slippingFault_horizontal_smoke.xml index 6b59abebbbf..7584144aec3 100644 --- a/inputFiles/lagrangianContactMechanics/ContactMechanics_slippingFault_horizontal_smoke.xml +++ b/inputFiles/lagrangianContactMechanics/ContactMechanics_slippingFault_horizontal_smoke.xml @@ -3,132 +3,64 @@ - + - - - - - - - - - - - - - - - - - + + + + + + - + + + + + - + + - + - + - - - - - - - - - - - - - + + + diff --git a/inputFiles/lagrangianContactMechanics/ContactMechanics_slippingFault_vertical_smoke.xml b/inputFiles/lagrangianContactMechanics/ContactMechanics_slippingFault_vertical_smoke.xml index b520f3498c6..9a6ca75ddba 100644 --- a/inputFiles/lagrangianContactMechanics/ContactMechanics_slippingFault_vertical_smoke.xml +++ b/inputFiles/lagrangianContactMechanics/ContactMechanics_slippingFault_vertical_smoke.xml @@ -3,107 +3,64 @@ - + - - - + + + + + + - - + + + + + - - - + + - - + - + - + + - - - - - - - - - - - - - - diff --git a/inputFiles/lagrangianContactMechanics/ContactMechanics_PassingCrack_base.xml b/inputFiles/lagrangianContactMechanics/PassingCrack_base.xml similarity index 78% rename from inputFiles/lagrangianContactMechanics/ContactMechanics_PassingCrack_base.xml rename to inputFiles/lagrangianContactMechanics/PassingCrack_base.xml index c309080025d..14ccc3cd8f8 100644 --- a/inputFiles/lagrangianContactMechanics/ContactMechanics_PassingCrack_base.xml +++ b/inputFiles/lagrangianContactMechanics/PassingCrack_base.xml @@ -3,22 +3,6 @@ - - - - - - - - @@ -124,9 +103,9 @@ - + diff --git a/inputFiles/lagrangianContactMechanics/PassingCrack_smoke.xml b/inputFiles/lagrangianContactMechanics/PassingCrack_smoke.xml new file mode 100644 index 00000000000..5af8f7534a5 --- /dev/null +++ b/inputFiles/lagrangianContactMechanics/PassingCrack_smoke.xml @@ -0,0 +1,58 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/lagrangianContactMechanics/ContactMechanics_SimpleCubes_base.xml b/inputFiles/lagrangianContactMechanics/SimpleCubes_base.xml similarity index 81% rename from inputFiles/lagrangianContactMechanics/ContactMechanics_SimpleCubes_base.xml rename to inputFiles/lagrangianContactMechanics/SimpleCubes_base.xml index 690ecf823a3..2098ad04708 100644 --- a/inputFiles/lagrangianContactMechanics/ContactMechanics_SimpleCubes_base.xml +++ b/inputFiles/lagrangianContactMechanics/SimpleCubes_base.xml @@ -3,25 +3,6 @@ - - - - - - - - - diff --git a/inputFiles/lagrangianContactMechanics/SimpleCubes_smoke.xml b/inputFiles/lagrangianContactMechanics/SimpleCubes_smoke.xml new file mode 100644 index 00000000000..157f01ceda1 --- /dev/null +++ b/inputFiles/lagrangianContactMechanics/SimpleCubes_smoke.xml @@ -0,0 +1,54 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/lagrangianContactMechanics/ContactMechanics_SingleFracCompression_base.xml b/inputFiles/lagrangianContactMechanics/SingleFracCompression_base.xml old mode 100755 new mode 100644 similarity index 81% rename from inputFiles/lagrangianContactMechanics/ContactMechanics_SingleFracCompression_base.xml rename to inputFiles/lagrangianContactMechanics/SingleFracCompression_base.xml index 28b80df0368..2da12a025df --- a/inputFiles/lagrangianContactMechanics/ContactMechanics_SingleFracCompression_base.xml +++ b/inputFiles/lagrangianContactMechanics/SingleFracCompression_base.xml @@ -4,27 +4,6 @@ - - - - - - - - - diff --git a/inputFiles/lagrangianContactMechanics/SingleFracCompression_benchmark.xml b/inputFiles/lagrangianContactMechanics/SingleFracCompression_benchmark.xml new file mode 100644 index 00000000000..346528e4bea --- /dev/null +++ b/inputFiles/lagrangianContactMechanics/SingleFracCompression_benchmark.xml @@ -0,0 +1,75 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/lagrangianContactMechanics/SingleFracCompression_smoke.xml b/inputFiles/lagrangianContactMechanics/SingleFracCompression_smoke.xml new file mode 100644 index 00000000000..a724d175616 --- /dev/null +++ b/inputFiles/lagrangianContactMechanics/SingleFracCompression_smoke.xml @@ -0,0 +1,75 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/lagrangianContactMechanics/SlippingFault_base.xml b/inputFiles/lagrangianContactMechanics/SlippingFault_base.xml new file mode 100644 index 00000000000..5299794fe0f --- /dev/null +++ b/inputFiles/lagrangianContactMechanics/SlippingFault_base.xml @@ -0,0 +1,59 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/lagrangianContactMechanics/SlippingFault_horizontal_smoke.xml b/inputFiles/lagrangianContactMechanics/SlippingFault_horizontal_smoke.xml new file mode 100644 index 00000000000..2d81f9a0593 --- /dev/null +++ b/inputFiles/lagrangianContactMechanics/SlippingFault_horizontal_smoke.xml @@ -0,0 +1,134 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/almContactMechanics/ALM_stickFault_vertical_smoke.xml b/inputFiles/lagrangianContactMechanics/SlippingFault_vertical_smoke.xml similarity index 73% rename from inputFiles/almContactMechanics/ALM_stickFault_vertical_smoke.xml rename to inputFiles/lagrangianContactMechanics/SlippingFault_vertical_smoke.xml index 4eb31ff467b..8c70450485d 100644 --- a/inputFiles/almContactMechanics/ALM_stickFault_vertical_smoke.xml +++ b/inputFiles/lagrangianContactMechanics/SlippingFault_vertical_smoke.xml @@ -3,19 +3,19 @@ - + @@ -26,7 +26,7 @@ origin="{0.0, 0.0, 0.0}" lengthVector="{0.0, 1.0, 0.0}" widthVector="{0.0, 0.0, 1.0}" - dimensions="{ 500, 200 }"/> + dimensions="{ 180, 10 }"/> + dimensions="{ 180, 10 }"/> @@ -56,20 +56,20 @@ scale="1"/> + setNames="{ xneg }"/> + setNames="{ yneg, ypos, xneg, xpos }"/> + setNames="{ zneg, zpos }"/> - + scale="-4.25e6"/> - + plotFileRoot="faultSlip_vertical"/> - + name="restartOutput"/> + diff --git a/inputFiles/lagrangianContactMechanics/Sneddon_contactMechanics_base.xml b/inputFiles/lagrangianContactMechanics/Sneddon_base.xml old mode 100755 new mode 100644 similarity index 96% rename from inputFiles/lagrangianContactMechanics/Sneddon_contactMechanics_base.xml rename to inputFiles/lagrangianContactMechanics/Sneddon_base.xml index 64b949a396e..b4d2d0a3a22 --- a/inputFiles/lagrangianContactMechanics/Sneddon_contactMechanics_base.xml +++ b/inputFiles/lagrangianContactMechanics/Sneddon_base.xml @@ -27,11 +27,6 @@ name="FE1" order="1"/> - - - - @@ -114,7 +109,7 @@ setNames="{ core }"/> - + + + + + + + + + + + + + + + + + diff --git a/inputFiles/lagrangianContactMechanics/Sneddon_smoke.xml b/inputFiles/lagrangianContactMechanics/Sneddon_smoke.xml new file mode 100644 index 00000000000..5e64b204ad1 --- /dev/null +++ b/inputFiles/lagrangianContactMechanics/Sneddon_smoke.xml @@ -0,0 +1,35 @@ + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/lagrangianContactMechanics/ContactMechanics_TFrac_base.xml b/inputFiles/lagrangianContactMechanics/TFrac_base.xml similarity index 97% rename from inputFiles/lagrangianContactMechanics/ContactMechanics_TFrac_base.xml rename to inputFiles/lagrangianContactMechanics/TFrac_base.xml index d8d82afdd12..9c3a4f54cdd 100644 --- a/inputFiles/lagrangianContactMechanics/ContactMechanics_TFrac_base.xml +++ b/inputFiles/lagrangianContactMechanics/TFrac_base.xml @@ -43,11 +43,6 @@ name="FE1" order="1"/> - - - - @@ -159,6 +154,7 @@ inputVarNames="{ time }" coordinates="{ 0.0, 1.0 }" values="{ 0.0, 1.e0 }"/> + diff --git a/inputFiles/lagrangianContactMechanics/TFrac_benchmark.xml b/inputFiles/lagrangianContactMechanics/TFrac_benchmark.xml new file mode 100644 index 00000000000..60f1a1f537a --- /dev/null +++ b/inputFiles/lagrangianContactMechanics/TFrac_benchmark.xml @@ -0,0 +1,36 @@ + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/lagrangianContactMechanics/TFrac_smoke.xml b/inputFiles/lagrangianContactMechanics/TFrac_smoke.xml new file mode 100644 index 00000000000..5ab56b69af9 --- /dev/null +++ b/inputFiles/lagrangianContactMechanics/TFrac_smoke.xml @@ -0,0 +1,42 @@ + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/lagrangianContactMechanics/ContactMechanics_UnstructuredCrack_base.xml b/inputFiles/lagrangianContactMechanics/UnstructuredCrack_base.xml similarity index 80% rename from inputFiles/lagrangianContactMechanics/ContactMechanics_UnstructuredCrack_base.xml rename to inputFiles/lagrangianContactMechanics/UnstructuredCrack_base.xml index 4c4c2bb1379..b8de0e0eaea 100644 --- a/inputFiles/lagrangianContactMechanics/ContactMechanics_UnstructuredCrack_base.xml +++ b/inputFiles/lagrangianContactMechanics/UnstructuredCrack_base.xml @@ -3,23 +3,6 @@ - - - - - - - - - @@ -136,7 +114,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/UnstructuredCrack_benchmark.xml b/inputFiles/lagrangianContactMechanics/UnstructuredCrack_benchmark.xml new file mode 100644 index 00000000000..bb32766c947 --- /dev/null +++ b/inputFiles/lagrangianContactMechanics/UnstructuredCrack_benchmark.xml @@ -0,0 +1,74 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/lagrangianContactMechanics/UnstructuredCrack_smoke.xml b/inputFiles/lagrangianContactMechanics/UnstructuredCrack_smoke.xml new file mode 100644 index 00000000000..faf2455cbf1 --- /dev/null +++ b/inputFiles/lagrangianContactMechanics/UnstructuredCrack_smoke.xml @@ -0,0 +1,74 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/lagrangianContactMechanics/contactMechanics.ats b/inputFiles/lagrangianContactMechanics/contactMechanics.ats index b522d942a63..bc4a2e43d0e 100644 --- a/inputFiles/lagrangianContactMechanics/contactMechanics.ats +++ b/inputFiles/lagrangianContactMechanics/contactMechanics.ats @@ -21,7 +21,7 @@ decks = [ check_step=2, restartcheck_params=RestartcheckParameters(**restartcheck_params)), TestDeck( - name="Sneddon_contactMechanics_smoke", + name="ContactMechanics_Sneddon_smoke", description= "Testing Sneddon problem using contact mechanics (structured grid)", partitions=((1, 1, 1), ), @@ -51,6 +51,45 @@ decks = [ partitions=((1, 1, 1), (2, 1, 1)), restart_step=5, check_step=10, + restartcheck_params=RestartcheckParameters(**restartcheck_params)), + TestDeck( + name="ALM_SimpleCubes_smoke", + description= + "Two cubes with a fracture separating them (structured grid)", + partitions=((1, 1, 1), (2, 2, 2), (1, 3, 3)), + restart_step=10, + check_step=20, + restartcheck_params=RestartcheckParameters(**restartcheck_params)), + TestDeck( + name="ALM_UnstructuredCrack_smoke", + description="A thick plane with a crack in it (unstructured grid)", + partitions=((1, 1, 1), ), + restart_step=1, + check_step=2, + restartcheck_params=RestartcheckParameters(**restartcheck_params)), + TestDeck( + name="ALM_Sneddon_smoke", + description= + "Testing Sneddon problem using contact mechanics (structured grid)", + partitions=((1, 1, 1), ), + restart_step=1, + check_step=2, + restartcheck_params=RestartcheckParameters(**restartcheck_params)), + TestDeck( + name="ALM_TFrac_smoke", + description= + "Two fractures intersecting at a right angle (structured grid)", + partitions=((1, 1, 1), (2, 2, 1)), + restart_step=1, + check_step=2, + restartcheck_params=RestartcheckParameters(**restartcheck_params)), + TestDeck( + name="ALM_SingleFracCompression_smoke", + description= + "Single tilted fracture subjected to remote compression (unstructured grid)", + partitions=((1, 1, 1), ), + restart_step=1, + check_step=2, restartcheck_params=RestartcheckParameters(**restartcheck_params)) ] diff --git a/src/coreComponents/constitutive/ConstitutivePassThru.hpp b/src/coreComponents/constitutive/ConstitutivePassThru.hpp index bef43e6e7ce..decabffbb6a 100644 --- a/src/coreComponents/constitutive/ConstitutivePassThru.hpp +++ b/src/coreComponents/constitutive/ConstitutivePassThru.hpp @@ -49,6 +49,7 @@ #include "permeability/ProppantPermeability.hpp" #include "permeability/SlipDependentPermeability.hpp" #include "permeability/WillisRichardsPermeability.hpp" +#include "contact/CoulombFriction.hpp" namespace geos @@ -83,6 +84,20 @@ struct ConstitutivePassThru< ElasticIsotropic > } }; +/** + * Specialization for models that derive from CoulombFriction. + */ +template<> +struct ConstitutivePassThru< CoulombFriction > +{ + template< typename LAMBDA > + static + void execute( ConstitutiveBase & constitutiveRelation, LAMBDA && lambda ) + { + ConstitutivePassThruHandler< CoulombFriction >::execute( constitutiveRelation, + std::forward< LAMBDA >( lambda ) ); + } +}; /** * Specialization for models that derive from SolidBase. diff --git a/src/coreComponents/constitutive/contact/CoulombFriction.cpp b/src/coreComponents/constitutive/contact/CoulombFriction.cpp index 32a5ee4f765..1199f652d7e 100644 --- a/src/coreComponents/constitutive/contact/CoulombFriction.cpp +++ b/src/coreComponents/constitutive/contact/CoulombFriction.cpp @@ -73,7 +73,7 @@ void CoulombFriction::allocateConstitutiveData( Group & parent, } -CoulombFrictionUpdates CoulombFriction::createKernelWrapper() const +CoulombFrictionUpdates CoulombFriction::createKernelUpdates() const { return CoulombFrictionUpdates( m_displacementJumpThreshold, m_shearStiffness, diff --git a/src/coreComponents/constitutive/contact/CoulombFriction.hpp b/src/coreComponents/constitutive/contact/CoulombFriction.hpp index 02c56ed0926..2ed7d182536 100644 --- a/src/coreComponents/constitutive/contact/CoulombFriction.hpp +++ b/src/coreComponents/constitutive/contact/CoulombFriction.hpp @@ -90,6 +90,41 @@ class CoulombFrictionUpdates : public FrictionBaseUpdates arraySlice1d< real64 const > const & tractionVector, integer & fractureState ) const override final; + GEOS_HOST_DEVICE + inline + virtual void updateTraction( arraySlice1d< real64 const > const & oldDispJump, + arraySlice1d< real64 const > const & dispJump, + arraySlice1d< real64 const > const & penalty, + arraySlice1d< real64 const > const & traction, + bool const symmetric, + bool const fixedLimitTau, + real64 const normalTractionTolerance, + real64 const tangentialTractionTolerance, + real64 ( &dTraction_dDispJump )[3][3], + real64 ( &tractionNew )[3], + integer & fractureState ) const override final; + + + GEOS_HOST_DEVICE + inline + virtual void updateTractionOnly( arraySlice1d< real64 const > const & dispJump, + arraySlice1d< real64 const > const & deltaDispJump, + arraySlice1d< real64 const > const & penalty, + arraySlice1d< real64 const > const & traction, + arraySlice1d< real64 > const & tractionNew ) const override final; + + GEOS_HOST_DEVICE + inline + virtual void constraintCheck( arraySlice1d< real64 const > const & dispJump, + arraySlice1d< real64 const > const & deltaDispJump, + arraySlice1d< real64 > const & tractionVector, + integer const fractureState, + real64 const normalTractionTolerance, + real64 const normalDisplacementTolerance, + real64 const slidingTolerance, + real64 const slidingCheckTolerance, + integer & condConv ) const override final; + private: /// The shear stiffness real64 m_shearStiffness; @@ -155,7 +190,7 @@ class CoulombFriction : public FrictionBase * @brief Create an update kernel wrapper. * @return the wrapper */ - KernelWrapper createKernelWrapper() const; + KernelWrapper createKernelUpdates() const; protected: @@ -302,6 +337,252 @@ inline void CoulombFrictionUpdates::updateFractureState( localIndex const k, } } +GEOS_HOST_DEVICE +inline void CoulombFrictionUpdates::updateTraction( arraySlice1d< real64 const > const & oldDispJump, + arraySlice1d< real64 const > const & dispJump, + arraySlice1d< real64 const > const & penalty, + arraySlice1d< real64 const > const & traction, + bool const symmetric, + bool const fixedLimitTau, + real64 const normalTractionTolerance, + real64 const tangentialTractionTolerance, + real64 ( & dTraction_dDispJump )[3][3], + real64 ( & tractionNew ) [3], + integer & fractureState ) const +{ + + using namespace fields::contact; + + real64 dLimitTangentialTractionNorm_dTraction = 0.0; + real64 limitTau = 0.0; + + // Compute the trial traction + real64 tractionTrial[ 3 ]; + tractionTrial[ 0 ] = traction[0] + penalty[0] * dispJump[0]; + tractionTrial[ 1 ] = traction[1] + penalty[1] * (dispJump[1] - oldDispJump[1]); + tractionTrial[ 2 ] = traction[2] + penalty[1] * (dispJump[2] - oldDispJump[2]); + + // Compute tangential trial traction norm + real64 const tau[2] = { tractionTrial[1], + tractionTrial[2] }; + real64 const tractionTrialNorm = LvArray::tensorOps::l2Norm< 2 >( tau ); + + // If normal tangential trial is positive (opening) + fractureState = FractureState::Stick; + if( tractionTrial[ 0 ] > normalTractionTolerance ) + { + tractionNew[0] = 0.0; + dTraction_dDispJump[0][0] = 0.0; + fractureState = FractureState::Open; + } + else + { + tractionNew[0] = tractionTrial[0]; + dTraction_dDispJump[0][0] = -penalty[ 0 ]; + } + + // Compute limit Tau + if( fixedLimitTau ) + { + limitTau = computeLimitTangentialTractionNorm( traction[0], + dLimitTangentialTractionNorm_dTraction ); + } + else + { + limitTau = computeLimitTangentialTractionNorm( tractionNew[0], + dLimitTangentialTractionNorm_dTraction ); + } + + if( tractionTrialNorm <= tangentialTractionTolerance ) + { + // It is needed for the first iteration (both t and u are equal to zero) + dTraction_dDispJump[1][1] = -penalty[1]; + dTraction_dDispJump[2][2] = -penalty[1]; + + tractionNew[1] = tractionTrial[1]; + tractionNew[2] = tractionTrial[2]; + + if( fractureState != FractureState::Open ) + fractureState = FractureState::Stick; + } + else if( limitTau <= tangentialTractionTolerance ) + { + dTraction_dDispJump[1][1] = 0.0; + dTraction_dDispJump[2][2] = 0.0; + + tractionNew[1] = (fixedLimitTau) ? tractionTrial[1] : 0.0; + tractionNew[2] = (fixedLimitTau) ? tractionTrial[2] : 0.0; + + if( fractureState != FractureState::Open ) + fractureState = FractureState::Slip; + } + else + { + // Compute psi and dpsi + //real64 const psi = std::tanh( tractionTrialNorm/limitTau ); + //real64 const dpsi = 1.0-std::pow(psi,2); + real64 const psi = ( tractionTrialNorm > limitTau) ? 1.0 : tractionTrialNorm/limitTau; + real64 const dpsi = ( tractionTrialNorm > limitTau) ? 0.0 : 1.0; + + if( fractureState != FractureState::Open ) + { + fractureState = ( tractionTrialNorm > limitTau) ? FractureState::Slip : FractureState::Stick; + } + + // Two symmetric 2x2 matrices + real64 dNormTTdgT[ 3 ]; + dNormTTdgT[ 0 ] = tractionTrial[ 1 ] * tractionTrial[ 1 ]; + dNormTTdgT[ 1 ] = tractionTrial[ 2 ] * tractionTrial[ 2 ]; + dNormTTdgT[ 2 ] = tractionTrial[ 1 ] * tractionTrial[ 2 ]; + + real64 dTdgT[ 3 ]; + dTdgT[ 0 ] = (tractionTrialNorm * tractionTrialNorm - dNormTTdgT[0]); + dTdgT[ 1 ] = (tractionTrialNorm * tractionTrialNorm - dNormTTdgT[1]); + dTdgT[ 2 ] = -dNormTTdgT[2]; + + LvArray::tensorOps::scale< 3 >( dNormTTdgT, 1. / std::pow( tractionTrialNorm, 2 ) ); + LvArray::tensorOps::scale< 3 >( dTdgT, 1. / std::pow( tractionTrialNorm, 3 ) ); + + // Compute dTdDispJump + dTraction_dDispJump[1][1] = -penalty[1] * ( + dpsi * dNormTTdgT[0] + + psi * dTdgT[0] * limitTau ); + dTraction_dDispJump[2][2] = -penalty[1] * ( + dpsi * dNormTTdgT[1] + + psi * dTdgT[1] * limitTau ); + dTraction_dDispJump[1][2] = -penalty[1] * ( + dpsi * dNormTTdgT[2] + + psi * dTdgT[2] * limitTau ); + dTraction_dDispJump[2][1] = dTraction_dDispJump[1][2]; + + if( !symmetric ) + { + // Nonsymetric term + dTraction_dDispJump[1][0] = -dTraction_dDispJump[0][0] * m_frictionCoefficient * + tractionTrial[1] * (psi/tractionTrialNorm - dpsi/limitTau); + dTraction_dDispJump[2][0] = -dTraction_dDispJump[0][0] * m_frictionCoefficient * + tractionTrial[2] * (psi/tractionTrialNorm - dpsi/limitTau); + } + + LvArray::tensorOps::scale< 3 >( tractionTrial, (psi * limitTau)/tractionTrialNorm ); + tractionNew[1] = tractionTrial[1]; + tractionNew[2] = tractionTrial[2]; + } + +} + +GEOS_HOST_DEVICE +inline void CoulombFrictionUpdates::updateTractionOnly( arraySlice1d< real64 const > const & dispJump, + arraySlice1d< real64 const > const & deltaDispJump, + arraySlice1d< real64 const > const & penalty, + arraySlice1d< real64 const > const & traction, + arraySlice1d< real64 > const & tractionNew ) const +{ + + // TODO: Pass this tol as an argument or define a new class member + real64 const zero = LvArray::NumericLimits< real64 >::epsilon; + + tractionNew[0] = traction[0] + penalty[0] * dispJump[0]; + tractionNew[1] = traction[1] + penalty[1] * deltaDispJump[1]; + tractionNew[2] = traction[2] + penalty[1] * deltaDispJump[2]; + + real64 const tau[2] = { tractionNew[1], + tractionNew[2] }; + real64 const currentTau = LvArray::tensorOps::l2Norm< 2 >( tau ); + + real64 dLimitTangentialTractionNorm_dTraction = 0.0; + real64 const limitTau = computeLimitTangentialTractionNorm( tractionNew[0], + dLimitTangentialTractionNorm_dTraction ); + + // Compute psi + real64 psi; + if( limitTau < zero ) + { + psi = 1.0; + } + else + { + //psi = std::tanh(currentTau / limitTau); + psi = (currentTau > limitTau ) ? 1.0 : currentTau/limitTau; + } + + // Compute the new tangential traction + if( limitTau > zero && currentTau > zero ) + { + tractionNew[1] *= limitTau * psi / currentTau; + tractionNew[2] *= limitTau * psi / currentTau; + } + else + { + tractionNew[1] *= 0.0; + tractionNew[2] *= 0.0; + } +} + +GEOS_HOST_DEVICE +inline void CoulombFrictionUpdates::constraintCheck( arraySlice1d< real64 const > const & dispJump, + arraySlice1d< real64 const > const & deltaDispJump, + arraySlice1d< real64 > const & tractionVector, + integer const fractureState, + real64 const normalTractionTolerance, + real64 const normalDisplacementTolerance, + real64 const slidingTolerance, + real64 const slidingCheckTolerance, + integer & condConv ) const +{ + + using namespace fields::contact; + + // Compute the slip + real64 const deltaDisp[2] = { deltaDispJump[1], + deltaDispJump[2] }; + + real64 const deltaDispNorm = LvArray::tensorOps::l2Norm< 2 >( deltaDisp ); + + // Compute current Tau and limit Tau + real64 const tau[2] = { tractionVector[1], + tractionVector[2] }; + real64 const currentTau = LvArray::tensorOps::l2Norm< 2 >( tau ); + + real64 dLimitTangentialTractionNorm_dTraction = 0.0; + real64 const limitTau = computeLimitTangentialTractionNorm( tractionVector[0], + dLimitTangentialTractionNorm_dTraction ); + + condConv = 0; + // Case 1: if it is open + if( tractionVector[0] >= normalTractionTolerance ) + { + if( fractureState != FractureState::Open ) + { + condConv = 1; + } + tractionVector[0] = 0.0; + tractionVector[1] = 0.0; + tractionVector[2] = 0.0; + } + else + { + // Case 2: compenetration + if(( LvArray::math::abs( dispJump[0] ) > normalDisplacementTolerance ) && + (fractureState != FractureState::Open)) + { + condConv = 2; + } + // Case 3: it is stick and dg is greater than 0 + if( fractureState == FractureState::Stick && + deltaDispNorm > slidingTolerance ) + { + condConv = 3; + } + + // Case 4: the elastic tangential traction is greater than the limit + if( currentTau > (LvArray::math::abs( limitTau ) * (1.0 + slidingCheckTolerance)) ) + { + condConv = 4; + } + } +} + } /* namespace constitutive */ } /* namespace geos */ diff --git a/src/coreComponents/constitutive/contact/FrictionBase.hpp b/src/coreComponents/constitutive/contact/FrictionBase.hpp index 13186378e42..e05923f2998 100644 --- a/src/coreComponents/constitutive/contact/FrictionBase.hpp +++ b/src/coreComponents/constitutive/contact/FrictionBase.hpp @@ -89,7 +89,81 @@ class FrictionBaseUpdates integer & fractureState ) const { GEOS_UNUSED_VAR( k, dispJump, tractionVector, fractureState ); } + /** + * @brief Update the trial traction vector ( return mapping ) + * @param[in] oldDispJump the displacement jump of the previous time step + * @param[in] dispJump the displacement jump of the current time step + * @param[in] penalty the penalty coefficients + * @param[in] traction the traction vector + * @param[in] symmetric flag to compute only the symmetric part of dTraction_dDispJump + * @param[in] fixedLimitTau flag to keep fixed the tangential stress + * @param[in] normalTractionTolerance normal traction tolerance (if tn > tol => fracture is open) + * @param[in] tangentialTractionTolerance tangential traction tolerance (if tau < tol => tau = 0) + * @param[out] dTraction_dDispJump matrix containing the derivatives of traction over the displacement jump + * @param[out] tractionNew the new traction vector + * @param[out] fractureState the fracture state + */ + GEOS_HOST_DEVICE + inline + virtual void updateTraction( arraySlice1d< real64 const > const & oldDispJump, + arraySlice1d< real64 const > const & dispJump, + arraySlice1d< real64 const > const & penalty, + arraySlice1d< real64 const > const & traction, + bool const symmetric, + bool const fixedLimitTau, + real64 const normalTractionTolerance, + real64 const tangentialTractionTolerance, + real64 ( & dTraction_dDispJump )[3][3], + real64 ( & tractionNew )[3], + integer & fractureState ) const + { + GEOS_UNUSED_VAR( oldDispJump, dispJump, penalty, traction, symmetric, fixedLimitTau, + normalTractionTolerance, tangentialTractionTolerance, + dTraction_dDispJump, tractionNew, fractureState ); + } + /** + * @brief Update the traction vector only ( return mapping ) + * @param[in] dispJump the displacement jump of the current time step + * @param[in] deltaDispJump the delta displacement jump of the current time step + * @param[in] penalty the penalty coefficients + * @param[in] traction the traction vector + * @param[out] tractionNew the new traction vector + */ + GEOS_HOST_DEVICE + inline + virtual void updateTractionOnly( arraySlice1d< real64 const > const & dispJump, + arraySlice1d< real64 const > const & deltaDispJump, + arraySlice1d< real64 const > const & penalty, + arraySlice1d< real64 const > const & traction, + arraySlice1d< real64 > const & tractionNew ) const + { GEOS_UNUSED_VAR( dispJump, deltaDispJump, penalty, traction, tractionNew ); } + + /** + * @brief Check for the constraint satisfaction + * @param[in] dispJump the displacement jump of the current time step + * @param[in] deltaDispJump the delta displacement jump of the current time step + * @param[in] tractionVector the traction vector + * @param[in] fractureState the fracture check + * @param[in] tolerance the tolerance + * @param[out] condConv flag indicating the result of the check + */ + GEOS_HOST_DEVICE + inline + virtual void constraintCheck( arraySlice1d< real64 const > const & dispJump, + arraySlice1d< real64 const > const & deltaDispJump, + arraySlice1d< real64 > const & tractionVector, + integer const fractureState, + real64 const normalTractionTolerance, + real64 const normalDisplacementTolerance, + real64 const slidingTolerance, + real64 const slidingCheckTolerance, + integer & condConv ) const + { + GEOS_UNUSED_VAR( dispJump, deltaDispJump, tractionVector, fractureState, + normalTractionTolerance, normalDisplacementTolerance, slidingTolerance, + slidingCheckTolerance, condConv ); + } /** * @brief Evaluate the limit tangential traction norm and return the derivative wrt normal traction diff --git a/src/coreComponents/constitutive/contact/FrictionlessContact.cpp b/src/coreComponents/constitutive/contact/FrictionlessContact.cpp index 557bdd81a0d..e7fb6bb8ba6 100644 --- a/src/coreComponents/constitutive/contact/FrictionlessContact.cpp +++ b/src/coreComponents/constitutive/contact/FrictionlessContact.cpp @@ -35,7 +35,7 @@ FrictionlessContact::FrictionlessContact( string const & name, FrictionlessContact::~FrictionlessContact() {} -FrictionlessContactUpdates FrictionlessContact::createKernelWrapper() const +FrictionlessContactUpdates FrictionlessContact::createKernelUpdates() const { return FrictionlessContactUpdates( m_displacementJumpThreshold ); } diff --git a/src/coreComponents/constitutive/contact/FrictionlessContact.hpp b/src/coreComponents/constitutive/contact/FrictionlessContact.hpp index b4ce8dba853..83e53fd8e14 100644 --- a/src/coreComponents/constitutive/contact/FrictionlessContact.hpp +++ b/src/coreComponents/constitutive/contact/FrictionlessContact.hpp @@ -118,7 +118,7 @@ class FrictionlessContact : public FrictionBase * @brief Create an update kernel wrapper. * @return the wrapper */ - KernelWrapper createKernelWrapper() const; + KernelWrapper createKernelUpdates() const; /** * @struct Structure to hold scoped key names diff --git a/src/coreComponents/finiteElement/elementFormulations/LagrangeBasis1.hpp b/src/coreComponents/finiteElement/elementFormulations/LagrangeBasis1.hpp index 81d1762bcdb..2ea2d584a79 100644 --- a/src/coreComponents/finiteElement/elementFormulations/LagrangeBasis1.hpp +++ b/src/coreComponents/finiteElement/elementFormulations/LagrangeBasis1.hpp @@ -176,8 +176,7 @@ class LagrangeBasis1 GEOS_FORCE_INLINE constexpr static real64 gradientBubble( const real64 xi ) { - GEOS_UNUSED_VAR( xi ); - return -0.5*xi; + return -2.0*xi; } /** @@ -420,20 +419,29 @@ class LagrangeBasis1 static void valueFaceBubble( real64 const (&coords)[3], real64 (& N)[numSupportFaces] ) { - for( int a=0; a<2; ++a ) - { - N[ a*5 ] = LagrangeBasis1::valueBubble( coords[0] ) * - LagrangeBasis1::valueBubble( coords[1] ) * - LagrangeBasis1::value( a, coords[2] ); + N[ 0 ] = LagrangeBasis1::valueBubble( coords[0] ) * + LagrangeBasis1::value( 0, coords[1] ) * + LagrangeBasis1::valueBubble( coords[2] ); - N[ a*3+1 ] = LagrangeBasis1::valueBubble( coords[0] ) * - LagrangeBasis1::value( a, coords[1] ) * - LagrangeBasis1::valueBubble( coords[2] ); + N[ 1 ] = LagrangeBasis1::valueBubble( coords[0] ) * + LagrangeBasis1::valueBubble( coords[1] ) * + LagrangeBasis1::value( 0, coords[2] ); - N[ a+2 ] = LagrangeBasis1::value( a, coords[0] ) * - LagrangeBasis1::valueBubble( coords[1] ) * - LagrangeBasis1::valueBubble( coords[2] ); - } + N[ 2 ] = LagrangeBasis1::value( 0, coords[0] ) * + LagrangeBasis1::valueBubble( coords[1] ) * + LagrangeBasis1::valueBubble( coords[2] ); + + N[ 3 ] = LagrangeBasis1::value( 1, coords[0] ) * + LagrangeBasis1::valueBubble( coords[1] ) * + LagrangeBasis1::valueBubble( coords[2] ); + + N[ 4 ] = LagrangeBasis1::valueBubble( coords[0] ) * + LagrangeBasis1::value( 1, coords[1] ) * + LagrangeBasis1::valueBubble( coords[2] ); + + N[ 5 ] = LagrangeBasis1::valueBubble( coords[0] ) * + LagrangeBasis1::valueBubble( coords[1] ) * + LagrangeBasis1::value( 1, coords[2] ); } /** @@ -448,38 +456,65 @@ class LagrangeBasis1 static void gradientFaceBubble( real64 const (&coords)[3], real64 (& dNdXi)[numSupportFaces][3] ) { - for( int a=0; a<2; ++a ) - { - dNdXi[ a*5 ][0] = LagrangeBasis1::gradientBubble( coords[0] ) * - LagrangeBasis1::valueBubble( coords[1] ) * - LagrangeBasis1::value( a, coords[2] ); - dNdXi[ a*5 ][1] = LagrangeBasis1::valueBubble( coords[0] ) * - LagrangeBasis1::gradientBubble( coords[1] ) * - LagrangeBasis1::value( a, coords[2] ); - dNdXi[ a*5 ][2] = LagrangeBasis1::valueBubble( coords[0] ) * - LagrangeBasis1::valueBubble( coords[1] ) * - LagrangeBasis1::gradient( a, coords[2] ); - - dNdXi[ a*3+1 ][0] = LagrangeBasis1::gradientBubble( coords[0] ) * - LagrangeBasis1::value( a, coords[1] ) * - LagrangeBasis1::valueBubble( coords[2] ); - dNdXi[ a*3+1 ][1] = LagrangeBasis1::valueBubble( coords[0] ) * - LagrangeBasis1::gradient( a, coords[1] ) * - LagrangeBasis1::valueBubble( coords[2] ); - dNdXi[ a*3+1 ][2] = LagrangeBasis1::valueBubble( coords[0] ) * - LagrangeBasis1::value( a, coords[1] ) * - LagrangeBasis1::gradientBubble( coords[2] ); - - dNdXi[ a+2 ][0] = LagrangeBasis1::gradient( a, coords[0] ) * - LagrangeBasis1::valueBubble( coords[1] ) * - LagrangeBasis1::valueBubble( coords[2] ); - dNdXi[ a+2 ][1] = LagrangeBasis1::value( a, coords[0] ) * - LagrangeBasis1::gradientBubble( coords[1] ) * - LagrangeBasis1::valueBubble( coords[2] ); - dNdXi[ a+2 ][2] = LagrangeBasis1::value( a, coords[0] ) * - LagrangeBasis1::valueBubble( coords[1] ) * - LagrangeBasis1::gradientBubble( coords[2] ); - } + dNdXi[0][0] = LagrangeBasis1::gradientBubble( coords[0] ) * + LagrangeBasis1::value( 0, coords[1] ) * + LagrangeBasis1::valueBubble( coords[2] ); + dNdXi[0][1] = LagrangeBasis1::valueBubble( coords[0] ) * + LagrangeBasis1::gradient( 0, coords[1] ) * + LagrangeBasis1::valueBubble( coords[2] ); + dNdXi[0][2] = LagrangeBasis1::valueBubble( coords[0] ) * + LagrangeBasis1::value( 0, coords[1] ) * + LagrangeBasis1::gradientBubble( coords[2] ); + + dNdXi[1][0] = LagrangeBasis1::gradientBubble( coords[0] ) * + LagrangeBasis1::valueBubble( coords[1] ) * + LagrangeBasis1::value( 0, coords[2] ); + dNdXi[1][1] = LagrangeBasis1::valueBubble( coords[0] ) * + LagrangeBasis1::gradientBubble( coords[1] ) * + LagrangeBasis1::value( 0, coords[2] ); + dNdXi[1][2] = LagrangeBasis1::valueBubble( coords[0] ) * + LagrangeBasis1::valueBubble( coords[1] ) * + LagrangeBasis1::gradient( 0, coords[2] ); + + dNdXi[2][0] = LagrangeBasis1::gradient( 0, coords[0] ) * + LagrangeBasis1::valueBubble( coords[1] ) * + LagrangeBasis1::valueBubble( coords[2] ); + dNdXi[2][1] = LagrangeBasis1::value( 0, coords[0] ) * + LagrangeBasis1::gradientBubble( coords[1] ) * + LagrangeBasis1::valueBubble( coords[2] ); + dNdXi[2][2] = LagrangeBasis1::value( 0, coords[0] ) * + LagrangeBasis1::valueBubble( coords[1] ) * + LagrangeBasis1::gradientBubble( coords[2] ); + + dNdXi[3][0] = LagrangeBasis1::gradient( 1, coords[0] ) * + LagrangeBasis1::valueBubble( coords[1] ) * + LagrangeBasis1::valueBubble( coords[2] ); + dNdXi[3][1] = LagrangeBasis1::value( 1, coords[0] ) * + LagrangeBasis1::gradientBubble( coords[1] ) * + LagrangeBasis1::valueBubble( coords[2] ); + dNdXi[3][2] = LagrangeBasis1::value( 1, coords[0] ) * + LagrangeBasis1::valueBubble( coords[1] ) * + LagrangeBasis1::gradientBubble( coords[2] ); + + dNdXi[4][0] = LagrangeBasis1::gradientBubble( coords[0] ) * + LagrangeBasis1::value( 1, coords[1] ) * + LagrangeBasis1::valueBubble( coords[2] ); + dNdXi[4][1] = LagrangeBasis1::valueBubble( coords[0] ) * + LagrangeBasis1::gradient( 1, coords[1] ) * + LagrangeBasis1::valueBubble( coords[2] ); + dNdXi[4][2] = LagrangeBasis1::valueBubble( coords[0] ) * + LagrangeBasis1::value( 1, coords[1] ) * + LagrangeBasis1::gradientBubble( coords[2] ); + + dNdXi[5][0] = LagrangeBasis1::gradientBubble( coords[0] ) * + LagrangeBasis1::valueBubble( coords[1] ) * + LagrangeBasis1::value( 1, coords[2] ); + dNdXi[5][1] = LagrangeBasis1::valueBubble( coords[0] ) * + LagrangeBasis1::gradientBubble( coords[1] ) * + LagrangeBasis1::value( 1, coords[2] ); + dNdXi[5][2] = LagrangeBasis1::valueBubble( coords[0] ) * + LagrangeBasis1::valueBubble( coords[1] ) * + LagrangeBasis1::gradient( 1, coords[2] ); } /** diff --git a/src/coreComponents/physicsSolvers/CMakeLists.txt b/src/coreComponents/physicsSolvers/CMakeLists.txt index 387979f7a41..ed3ea6de356 100644 --- a/src/coreComponents/physicsSolvers/CMakeLists.txt +++ b/src/coreComponents/physicsSolvers/CMakeLists.txt @@ -6,8 +6,7 @@ set( physicsSolvers_headers SolverBase.hpp SolverBaseKernels.hpp SolverStatistics.hpp - FieldStatisticsBase.hpp - ) + FieldStatisticsBase.hpp ) # Specify solver sources set( physicsSolvers_sources diff --git a/src/coreComponents/physicsSolvers/contact/ContactFields.hpp b/src/coreComponents/physicsSolvers/contact/ContactFields.hpp index 17291a2dbf7..4792ca2d019 100644 --- a/src/coreComponents/physicsSolvers/contact/ContactFields.hpp +++ b/src/coreComponents/physicsSolvers/contact/ContactFields.hpp @@ -50,13 +50,13 @@ struct FractureState }; }; -DECLARE_FIELD( penalty, - "penalty", +DECLARE_FIELD( iterativePenalty, + "iterativePenalty", array2d< real64 >, - 0, + 1.e5, LEVEL_0, WRITE_AND_READ, - "Penalty coefficients" ); + "Penalty coefficients used in the iterative procedure of the Augmented Lagrangian Method" ); DECLARE_FIELD( rotationMatrix, "rotationMatrix", diff --git a/src/coreComponents/physicsSolvers/contact/SolidMechanicsALMBubbleKernels.hpp b/src/coreComponents/physicsSolvers/contact/SolidMechanicsALMBubbleKernels.hpp index 7764d9e0a50..4b76b7a4a3d 100644 --- a/src/coreComponents/physicsSolvers/contact/SolidMechanicsALMBubbleKernels.hpp +++ b/src/coreComponents/physicsSolvers/contact/SolidMechanicsALMBubbleKernels.hpp @@ -285,30 +285,7 @@ class ALMBubbleKernels : real64 strainBubbleMatrix[6][nBubbleUdof]; solidMechanicsALMKernelsHelper::assembleStrainOperator< 6, nBubbleUdof, numFacesPerElem >( strainBubbleMatrix, dBubbleNdX ); - // TODO: Use the following functions - //using namespace PDEUtilities; - //constexpr FunctionSpace displacementTrialSpace = FE_TYPE::template getFunctionSpace< numDofPerTrialSupportPoint >(); - //constexpr FunctionSpace displacementTestSpace = displacementTrialSpace; - //real64 Abb_bilinear[nBubbleUdof][nBubbleUdof]; - //BilinearFormUtilities::compute< displacementTrialSpace, - // displacementTestSpace, - // DifferentialOperator::SymmetricGradient, - // DifferentialOperator::SymmetricGradient > - //( - // Abb_bilinear, - // dBubbleNdX, - // stack.constitutiveStiffness, // fourth-order tensor handled via DiscretizationOps - // dBubbleNdX, - // -detJ ); - - - //LinearFormUtilities::compute< displacementTestSpace, - // DifferentialOperator::Identity > - //( - //stack.localResidualMomentum, - //N, - //stack.bodyForce, - //detJxW ); + // TODO: It would be nice use BilinearFormUtilities::compute real64 matBD[nBubbleUdof][6]; real64 Abb_gauss[nBubbleUdof][nBubbleUdof], Abu_gauss[nBubbleUdof][nUdof], Aub_gauss[nUdof][nBubbleUdof]; @@ -330,6 +307,22 @@ class ALMBubbleKernels : LvArray::tensorOps::scaledAdd< nBubbleUdof, nUdof >( stack.localAbu, Abu_gauss, -detJ ); LvArray::tensorOps::scaledAdd< nUdof, nBubbleUdof >( stack.localAub, Aub_gauss, -detJ ); + // Compute the initial stress + // The following block assumes a linear elastic constitutive model + real64 rb_gauss[nBubbleUdof]{}; + real64 strain[6] = {0}; + LvArray::tensorOps::Ri_eq_AijBj< 6, nUdof >( strain, strainMatrix, stack.uLocal ); + + real64 initStressLocal[ 6 ] = {0}; + LvArray::tensorOps::Ri_eq_AijBj< 6, 6 >( initStressLocal, stack.constitutiveStiffness, strain ); + for( localIndex c = 0; c < 6; ++c ) + { + initStressLocal[ c ] -= m_constitutiveUpdate.m_newStress( k, q, c ); + } + + LvArray::tensorOps::Ri_eq_AjiBj< nBubbleUdof, 6 >( rb_gauss, strainBubbleMatrix, initStressLocal ); + LvArray::tensorOps::scaledAdd< nBubbleUdof >( stack.localRb, rb_gauss, detJ ); + } GEOS_HOST_DEVICE diff --git a/src/coreComponents/physicsSolvers/contact/SolidMechanicsALMKernels.hpp b/src/coreComponents/physicsSolvers/contact/SolidMechanicsALMKernels.hpp index 7c63181a6c9..9f5c53f5132 100644 --- a/src/coreComponents/physicsSolvers/contact/SolidMechanicsALMKernels.hpp +++ b/src/coreComponents/physicsSolvers/contact/SolidMechanicsALMKernels.hpp @@ -52,6 +52,7 @@ class ALM : using Base::m_elemsToFaces; using Base::m_faceToNodes; using Base::m_finiteElementSpace; + using Base::m_constitutiveUpdate; using Base::m_dofNumber; using Base::m_bDofNumber; using Base::m_dofRankOffset; @@ -80,7 +81,8 @@ class ALM : CRSMatrixView< real64, globalIndex const > const inputMatrix, arrayView1d< real64 > const inputRhs, real64 const inputDt, - arrayView1d< localIndex const > const & faceElementList ): + arrayView1d< localIndex const > const & faceElementList, + bool const isSymmetric ): Base( nodeManager, edgeManager, faceManager, @@ -95,7 +97,8 @@ class ALM : inputRhs, inputDt, faceElementList ), - m_traction( elementSubRegion.getField< fields::contact::traction >().toViewConst()) + m_traction( elementSubRegion.getField< fields::contact::traction >().toViewConst()), + m_symmetric( isSymmetric ) {} //*************************************************************************** @@ -223,11 +226,6 @@ class ALM : } } - // The minus sign is consistent with the sign of the Jacobian - stack.localPenalty[0][0] = -m_penalty( k, 0 ); - stack.localPenalty[1][1] = -m_penalty( k, 1 ); - stack.localPenalty[2][2] = -m_penalty( k, 1 ); - for( int i=0; i::epsilon; constexpr int numUdofs = numNodesPerElem * 3 * 2; @@ -259,13 +258,24 @@ class ALM : real64 matRRtAtu[3][numUdofs], matDRtAtu[3][numUdofs]; real64 matRRtAtb[3][numBdofs], matDRtAtb[3][numBdofs]; - real64 dispJumpR[numUdofs]; - real64 oldDispJumpR[numUdofs]; real64 tractionR[numUdofs]; - real64 dispJumpRb[numBdofs]; - real64 oldDispJumpRb[numBdofs]; real64 tractionRb[numBdofs]; + real64 tractionNew[3]; + + integer fractureState; + m_constitutiveUpdate.updateTraction( m_oldDispJump[k], + m_dispJump[k], + m_penalty[k], + m_traction[k], + m_symmetric, + m_symmetric, + zero, + zero, + stack.localPenalty, + tractionNew, + fractureState ); + // transp(R) * Atu LvArray::tensorOps::Rij_eq_AkiBkj< 3, numUdofs, 3 >( matRRtAtu, stack.localRotationMatrix, stack.localAtu ); @@ -274,8 +284,8 @@ class ALM : stack.localAtb ); // Compute the traction contribute of the local residuals - LvArray::tensorOps::Ri_eq_AjiBj< numUdofs, 3 >( tractionR, matRRtAtu, stack.tLocal ); - LvArray::tensorOps::Ri_eq_AjiBj< numBdofs, 3 >( tractionRb, matRRtAtb, stack.tLocal ); + LvArray::tensorOps::Ri_eq_AjiBj< numUdofs, 3 >( tractionR, matRRtAtu, tractionNew ); + LvArray::tensorOps::Ri_eq_AjiBj< numBdofs, 3 >( tractionRb, matRRtAtb, tractionNew ); // D*RtAtu LvArray::tensorOps::Rij_eq_AikBkj< 3, numUdofs, 3 >( matDRtAtu, stack.localPenalty, @@ -307,18 +317,9 @@ class ALM : matRRtAtb ); // Compute the local residuals - LvArray::tensorOps::Ri_eq_AjiBj< numUdofs, 3 >( dispJumpR, matDRtAtu, stack.dispJumpLocal ); - LvArray::tensorOps::Ri_eq_AjiBj< numUdofs, 3 >( oldDispJumpR, matDRtAtu, stack.oldDispJumpLocal ); - LvArray::tensorOps::Ri_eq_AjiBj< numBdofs, 3 >( dispJumpRb, matDRtAtb, stack.dispJumpLocal ); - LvArray::tensorOps::Ri_eq_AjiBj< numBdofs, 3 >( oldDispJumpRb, matDRtAtb, stack.oldDispJumpLocal ); - LvArray::tensorOps::scaledAdd< numUdofs >( stack.localRu, tractionR, -1 ); - LvArray::tensorOps::scaledAdd< numUdofs >( stack.localRu, dispJumpR, 1 ); - LvArray::tensorOps::scaledAdd< numUdofs >( stack.localRu, oldDispJumpR, -1 ); LvArray::tensorOps::scaledAdd< numBdofs >( stack.localRb, tractionRb, -1 ); - LvArray::tensorOps::scaledAdd< numBdofs >( stack.localRb, dispJumpRb, 1 ); - LvArray::tensorOps::scaledAdd< numBdofs >( stack.localRb, oldDispJumpRb, -1 ); for( localIndex i=0; i < numUdofs; ++i ) { @@ -369,6 +370,8 @@ class ALM : arrayView2d< real64 const > const m_traction; + bool const m_symmetric; + }; /// The factory used to construct the kernel. @@ -379,48 +382,45 @@ using ALMFactory = finiteElement::InterfaceKernelFactory< ALM, CRSMatrixView< real64, globalIndex const > const, arrayView1d< real64 > const, real64 const, - arrayView1d< localIndex const > const >; - + arrayView1d< localIndex const > const, + bool const >; /** - * @brief A struct to compute rotation matrices + * @brief A struct to compute the traction after nonlinear solve */ -struct ComputeRotationMatricesKernel +struct ComputeTractionKernel { /** - * @brief Launch the kernel function to comute rotation matrices + * @brief Launch the kernel function to compute the traction * @tparam POLICY the type of policy used in the kernel launch + * @tparam CONTACT_WRAPPER the type of contact wrapper doing the fracture traction updates * @param[in] size the size of the subregion - * @param[in] faceNormal the array of array containing the face to nodes map - * @param[in] elemsToFaces the array of array containing the element to faces map - * @param[out] rotationMatrix the array containing the rotation matrices + * @param[in] penalty the array containing the tangential penalty matrix + * @param[in] traction the array containing the current traction + * @param[in] dispJump the array containing the displacement jump + * @param[in] deltaDispJump the array containing the delta displacement jump + * @param[out] tractionNew the array containing the new traction */ - template< typename POLICY > + template< typename POLICY, typename CONTACT_WRAPPER > static void launch( localIndex const size, - arrayView2d< real64 const > const & faceNormal, - ArrayOfArraysView< localIndex const > const & elemsToFaces, - arrayView3d< real64 > const & rotationMatrix ) + CONTACT_WRAPPER const & contactWrapper, + arrayView2d< real64 const > const & penalty, + arrayView2d< real64 const > const & traction, + arrayView2d< real64 const > const & dispJump, + arrayView2d< real64 const > const & deltaDispJump, + arrayView2d< real64 > const & tractionNew ) { forAll< POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const k ) { - localIndex const f0 = elemsToFaces[k][0]; - localIndex const f1 = elemsToFaces[k][1]; - - real64 Nbar[3]; - Nbar[0] = faceNormal[f0][0] - faceNormal[f1][0]; - Nbar[1] = faceNormal[f0][1] - faceNormal[f1][1]; - Nbar[2] = faceNormal[f0][2] - faceNormal[f1][2]; - - LvArray::tensorOps::normalize< 3 >( Nbar ); - computationalGeometry::RotationMatrix_3D( Nbar, rotationMatrix[k] ); + contactWrapper.updateTractionOnly( dispJump[k], deltaDispJump[k], + penalty[k], traction[k], tractionNew[k] ); } ); } - }; } // namespace SolidMechanicsALMKernels diff --git a/src/coreComponents/physicsSolvers/contact/SolidMechanicsALMKernelsBase.hpp b/src/coreComponents/physicsSolvers/contact/SolidMechanicsALMKernelsBase.hpp index d1debe35a5d..4f3504e3d5f 100644 --- a/src/coreComponents/physicsSolvers/contact/SolidMechanicsALMKernelsBase.hpp +++ b/src/coreComponents/physicsSolvers/contact/SolidMechanicsALMKernelsBase.hpp @@ -98,7 +98,7 @@ class ALMKernelsBase : m_rotationMatrix( elementSubRegion.getField< fields::contact::rotationMatrix >().toViewConst()), m_dispJump( elementSubRegion.getField< fields::contact::dispJump >().toView() ), m_oldDispJump( elementSubRegion.getField< fields::contact::oldDispJump >().toViewConst() ), - m_penalty( elementSubRegion.getField< fields::contact::penalty >().toViewConst() ) + m_penalty( elementSubRegion.getField< fields::contact::iterativePenalty >().toViewConst() ) {} //*************************************************************************** @@ -266,6 +266,182 @@ class ALMKernelsBase : }; +/** + * @brief A struct to compute rotation matrices + */ +struct ComputeRotationMatricesKernel +{ + + /** + * @brief Launch the kernel function to comute rotation matrices + * @tparam POLICY the type of policy used in the kernel launch + * @param[in] size the size of the subregion + * @param[in] faceNormal the array of array containing the face to nodes map + * @param[in] elemsToFaces the array of array containing the element to faces map + * @param[out] rotationMatrix the array containing the rotation matrices + */ + template< typename POLICY > + static void + launch( localIndex const size, + arrayView2d< real64 const > const & faceNormal, + ArrayOfArraysView< localIndex const > const & elemsToFaces, + arrayView3d< real64 > const & rotationMatrix ) + { + + forAll< POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const k ) + { + + localIndex const & f0 = elemsToFaces[k][0]; + localIndex const & f1 = elemsToFaces[k][1]; + + real64 Nbar[3]; + Nbar[0] = faceNormal[f0][0] - faceNormal[f1][0]; + Nbar[1] = faceNormal[f0][1] - faceNormal[f1][1]; + Nbar[2] = faceNormal[f0][2] - faceNormal[f1][2]; + + LvArray::tensorOps::normalize< 3 >( Nbar ); + computationalGeometry::RotationMatrix_3D( Nbar, rotationMatrix[k] ); + + } ); + } + +}; + +/** + * @brief A struct to check for constraint satisfaction + */ +struct ConstraintCheckKernel +{ + + /** + * @brief Launch the kernel function to check the constraint satisfaction + * @tparam POLICY the type of policy used in the kernel launch + * @tparam CONTACT_WRAPPER the type of contact wrapper doing the fracture traction updates + * @param[in] size the size of the subregion + * @param[in] traction the array containing the current traction + * @param[in] dispJump the array containing the displacement jump + * @param[in] deltaDispJump the array containing the delta displacement jump + * @param[in] normalTractionTolerance Check tolerance (normal traction) + * @param[in] normalDisplacementTolerance Check tolerance (compenetration) + * @param[in] slidingTolerance Check tolerance (sliding) + * @param[in] slidingCheckTolerance Check tolerance (if shear strass exceeds tauLim) + * @param[in] area interface element area + * @param[in] fractureState the array containing the fracture state + * @param[out] condConv the array containing the convergence flag: + * 0: Constraint conditions satisfied + * 1: Open + * 2: Compenetration + * 3: Slip exceeds sliding tolerance + * 4: Shear stress exceeds tauLim + */ + template< typename POLICY, typename CONTACT_WRAPPER > + static void + launch( localIndex const size, + CONTACT_WRAPPER const & contactWrapper, + arrayView1d< integer const > const & ghostRank, + arrayView2d< real64 > const & traction, + arrayView2d< real64 const > const & dispJump, + arrayView2d< real64 const > const & deltaDispJump, + arrayView1d< real64 const > const & normalTractionTolerance, + arrayView1d< real64 const > const & normalDisplacementTolerance, + arrayView1d< real64 const > const & slidingTolerance, + real64 const slidingCheckTolerance, + arrayView1d< real64 const > const & area, + arrayView1d< integer const > const & fractureState, + arrayView1d< integer > const & condConv ) + { + + forAll< POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const k ) + { + if( ghostRank[k] < 0 ) + { + contactWrapper.constraintCheck( dispJump[k], + deltaDispJump[k], + traction[k], + fractureState[k], + normalTractionTolerance[k], + normalDisplacementTolerance[k]*area[k], + slidingTolerance[k]*area[k], + slidingCheckTolerance, + condConv[k] ); + } + + } ); + } +}; + +/** + * @brief A struct to check for constraint satisfaction + */ +struct UpdateStateKernel +{ + + /** + * @brief Launch the kernel function to check the constraint satisfaction + * @tparam POLICY the type of policy used in the kernel launch + * @tparam CONTACT_WRAPPER the type of contact wrapper doing the fracture traction updates + * @param[in] size the size of the subregion + * @param[in] oldDispJump the array containing the old displacement jump (previous time step) + * @param[in] dispJump the array containing the displacement jump + * @param[in] penalty the array containing the penalty coefficients + * @param[in] symmetric flag to compute symmetric penalty matrix + * @param[in] normalTractionTolerance Check tolerance (normal traction) + * @param[in] traction the array containing the current traction + * @param[in] fractureState the array containing the fracture state + */ + template< typename POLICY, typename CONTACT_WRAPPER > + static void + launch( localIndex const size, + CONTACT_WRAPPER const & contactWrapper, + arrayView2d< real64 const > const & oldDispJump, + arrayView2d< real64 const > const & dispJump, + arrayView2d< real64 > const & penalty, + bool const symmetric, + arrayView1d< real64 const > const & normalTractionTolerance, + arrayView2d< real64 > const & traction, + arrayView1d< integer > const & fractureState ) + + { + forAll< POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const k ) + { + + real64 const zero = LvArray::NumericLimits< real64 >::epsilon; + + real64 localPenalty[3][3]{}; + real64 localTractionNew[3]{}; + contactWrapper.updateTraction( oldDispJump[k], + dispJump[k], + penalty[k], + traction[k], + symmetric, + false, + normalTractionTolerance[k], + zero, + localPenalty, + localTractionNew, + fractureState[k] ); + + if( fractureState[k] == fields::contact::FractureState::Open ) + { + + LvArray::tensorOps::fill< 3 >( localTractionNew, 0.0 ); + } + else if( LvArray::math::abs( localTractionNew[ 0 ] ) < normalTractionTolerance[k] ) + { + LvArray::tensorOps::fill< 3 >( localTractionNew, 0.0 ); + fractureState[k] = fields::contact::FractureState::Slip; + } + + LvArray::tensorOps::copy< 3 >( traction[k], localTractionNew ); + penalty[k][2] = -localPenalty[1][1]; + penalty[k][3] = -localPenalty[2][2]; + penalty[k][4] = -localPenalty[1][2]; + + } ); + } + +}; + } // namespace SolidMechanicsALMKernels } // namespace geos diff --git a/src/coreComponents/physicsSolvers/contact/SolidMechanicsALMSimultaneousKernels.hpp b/src/coreComponents/physicsSolvers/contact/SolidMechanicsALMSimultaneousKernels.hpp new file mode 100644 index 00000000000..f9cd307b6fe --- /dev/null +++ b/src/coreComponents/physicsSolvers/contact/SolidMechanicsALMSimultaneousKernels.hpp @@ -0,0 +1,435 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 Total, S.A + * Copyright (c) 2019- GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file SolidMechanicsALMSimultaneousKernels.hpp + */ + +#ifndef GEOS_PHYSICSSOLVERS_CONTACT_SOLIDMECHANICSALMSIMULTANEOUSKERNELS_HPP_ +#define GEOS_PHYSICSSOLVERS_CONTACT_SOLIDMECHANICSALMSIMULTANEOUSKERNELS_HPP_ + +#include "SolidMechanicsALMKernelsBase.hpp" + +namespace geos +{ + +namespace solidMechanicsALMKernels +{ + +/** + * @copydoc geos::finiteElement::ImplicitKernelBase + */ +template< typename CONSTITUTIVE_TYPE, + typename FE_TYPE > +class ALMSimultaneous : + public ALMKernelsBase< CONSTITUTIVE_TYPE, + FE_TYPE > +{ +public: + /// Alias for the base class. + using Base = ALMKernelsBase< CONSTITUTIVE_TYPE, + FE_TYPE >; + + /// Maximum number of nodes per element, which is equal to the maxNumTestSupportPointPerElem and + /// maxNumTrialSupportPointPerElem by definition. + static constexpr int numNodesPerElem = Base::maxNumTestSupportPointsPerElem; + + /// Compile time value for the number of quadrature points per element. + static constexpr int numQuadraturePointsPerElem = FE_TYPE::numQuadraturePoints; + + using Base::m_elemsToFaces; + using Base::m_faceToNodes; + using Base::m_finiteElementSpace; + using Base::m_dofNumber; + using Base::m_bDofNumber; + using Base::m_dofRankOffset; + using Base::m_X; + using Base::m_rotationMatrix; + using Base::m_penalty; + using Base::m_dispJump; + using Base::m_oldDispJump; + using Base::m_matrix; + using Base::m_rhs; + + /** + * @brief Constructor + * @copydoc geos::finiteElement::InterfaceKernelBase::InterfaceKernelBase + */ + ALMSimultaneous( NodeManager const & nodeManager, + EdgeManager const & edgeManager, + FaceManager const & faceManager, + localIndex const targetRegionIndex, + FaceElementSubRegion & elementSubRegion, + FE_TYPE const & finiteElementSpace, + CONSTITUTIVE_TYPE & inputConstitutiveType, + arrayView1d< globalIndex const > const uDofNumber, + arrayView1d< globalIndex const > const bDofNumber, + globalIndex const rankOffset, + CRSMatrixView< real64, globalIndex const > const inputMatrix, + arrayView1d< real64 > const inputRhs, + real64 const inputDt, + arrayView1d< localIndex const > const & faceElementList ): + Base( nodeManager, + edgeManager, + faceManager, + targetRegionIndex, + elementSubRegion, + finiteElementSpace, + inputConstitutiveType, + uDofNumber, + bDofNumber, + rankOffset, + inputMatrix, + inputRhs, + inputDt, + faceElementList ), + m_traction( elementSubRegion.getField< fields::contact::traction >().toViewConst()) + {} + + //*************************************************************************** + + /** + * @copydoc finiteElement::KernelBase::StackVariables + */ + struct StackVariables : public Base::StackVariables + { + + /// The number of displacement dofs per element. + static constexpr int numUdofs = numNodesPerElem * 3 * 2; + + /// The number of lagrange multiplier dofs per element. + static constexpr int numTdofs = 3; + + /// The number of bubble dofs per element. + static constexpr int numBdofs = 3*2; + +public: + + GEOS_HOST_DEVICE + StackVariables(): + Base::StackVariables(), + dispEqnRowIndices{}, + dispColIndices{}, + bEqnRowIndices{}, + bColIndices{}, + localRu{}, + localRb{}, + localAutAtu{ {} }, + localAbtAtb{ {} }, + localAbtAtu{ {} }, + localAutAtb{ {} }, + tLocal{} + {} + + /// C-array storage for the element local row degrees of freedom. + globalIndex dispEqnRowIndices[numUdofs]; + + /// C-array storage for the element local column degrees of freedom. + globalIndex dispColIndices[numUdofs]; + + /// C-array storage for the element local row degrees of freedom. + globalIndex bEqnRowIndices[numBdofs]; + + /// C-array storage for the element local column degrees of freedom. + globalIndex bColIndices[numBdofs]; + + /// C-array storage for the element local Ru residual vector. + real64 localRu[numUdofs]; + + /// C-array storage for the element local Rb residual vector. + real64 localRb[numBdofs]; + + /// C-array storage for the element local AutAtu matrix. + real64 localAutAtu[numUdofs][numUdofs]; + + /// C-array storage for the element local AbtAtb matrix. + real64 localAbtAtb[numBdofs][numBdofs]; + + /// C-array storage for the element local AbtAtu matrix. + real64 localAbtAtu[numBdofs][numUdofs]; + + /// C-array storage for the element local AbtAtu matrix. + real64 localAutAtb[numUdofs][numBdofs]; + + /// Stack storage for the element local lagrange multiplier vector + real64 tLocal[numTdofs]; + + }; + + //*************************************************************************** + + //START_kernelLauncher + template< typename POLICY, + typename KERNEL_TYPE > + static + real64 + kernelLaunch( localIndex const numElems, + KERNEL_TYPE const & kernelComponent ) + { + return Base::template kernelLaunch< POLICY, KERNEL_TYPE >( numElems, kernelComponent ); + } + //END_kernelLauncher + + /** + * @brief Copy global values from primary field to a local stack array. + * @copydoc ::geos::finiteElement::InterfaceKernelBase::setup + */ + GEOS_HOST_DEVICE + inline + void setup( localIndex const k, + StackVariables & stack ) const + { + constexpr int shift = numNodesPerElem * 3; + + constexpr int numTdofs = 3; + + int permutation[numNodesPerElem]; + m_finiteElementSpace.getPermutation( permutation ); + + localIndex const kf0 = m_elemsToFaces[k][0]; + localIndex const kf1 = m_elemsToFaces[k][1]; + for( localIndex a=0; a( tractionNew, stack.tLocal, -1.0 ); + LvArray::tensorOps::Ri_add_AijBj< 3, 3 >( tractionNew, stack.localPenalty, dispJump ); + + // If normal tangential trial is positive (opening) + //if (tractionNew[ 0 ] < -zero) + //{ + // tractionNew[0] = 0.0; + // stack.localPenalty[0][0] = 0.0; + //} + + // transp(R) * Atu + LvArray::tensorOps::Rij_eq_AkiBkj< 3, numUdofs, 3 >( matRRtAtu, stack.localRotationMatrix, + stack.localAtu ); + // transp(R) * Atb + LvArray::tensorOps::Rij_eq_AkiBkj< 3, numBdofs, 3 >( matRRtAtb, stack.localRotationMatrix, + stack.localAtb ); + + // Compute the traction contribute of the local residuals + LvArray::tensorOps::Ri_eq_AjiBj< numUdofs, 3 >( tractionR, matRRtAtu, tractionNew ); + LvArray::tensorOps::Ri_eq_AjiBj< numBdofs, 3 >( tractionRb, matRRtAtb, tractionNew ); + + // D*RtAtu + LvArray::tensorOps::Rij_eq_AikBkj< 3, numUdofs, 3 >( matDRtAtu, stack.localPenalty, + matRRtAtu ); + // D*RtAtb + LvArray::tensorOps::Rij_eq_AikBkj< 3, numBdofs, 3 >( matDRtAtb, stack.localPenalty, + matRRtAtb ); + + // R*RtAtu + LvArray::tensorOps::Rij_eq_AikBkj< 3, numUdofs, 3 >( matRRtAtu, stack.localRotationMatrix, + matDRtAtu ); + // R*RtAtb + LvArray::tensorOps::Rij_eq_AikBkj< 3, numBdofs, 3 >( matRRtAtb, stack.localRotationMatrix, + matDRtAtb ); + + // transp(Atu)*RRtAtu + LvArray::tensorOps::Rij_eq_AkiBkj< numUdofs, numUdofs, 3 >( stack.localAutAtu, stack.localAtu, + matRRtAtu ); + // transp(Atb)*RRtAtb + LvArray::tensorOps::Rij_eq_AkiBkj< numBdofs, numBdofs, 3 >( stack.localAbtAtb, stack.localAtb, + matRRtAtb ); + + // transp(Atb)*RRtAtu + LvArray::tensorOps::Rij_eq_AkiBkj< numBdofs, numUdofs, 3 >( stack.localAbtAtu, stack.localAtb, + matRRtAtu ); + + // transp(Atu)*RRtAtb + LvArray::tensorOps::Rij_eq_AkiBkj< numUdofs, numBdofs, 3 >( stack.localAutAtb, stack.localAtu, + matRRtAtb ); + + // Compute the local residuals + LvArray::tensorOps::scaledAdd< numUdofs >( stack.localRu, tractionR, 1 ); + + LvArray::tensorOps::scaledAdd< numBdofs >( stack.localRb, tractionRb, 1 ); + + for( localIndex i=0; i < numUdofs; ++i ) + { + localIndex const dof = LvArray::integerConversion< localIndex >( stack.dispEqnRowIndices[ i ] ); + + if( dof < 0 || dof >= m_matrix.numRows() ) continue; + + // Is it necessary? Each row should be indepenedent + RAJA::atomicAdd< parallelDeviceAtomic >( &m_rhs[dof], stack.localRu[i] ); + + // Fill in matrix + m_matrix.template addToRowBinarySearchUnsorted< parallelDeviceAtomic >( dof, + stack.dispColIndices, + stack.localAutAtu[i], + numUdofs ); + + m_matrix.template addToRowBinarySearchUnsorted< parallelDeviceAtomic >( dof, + stack.bColIndices, + stack.localAutAtb[i], + numBdofs ); + } + + for( localIndex i=0; i < numBdofs; ++i ) + { + localIndex const dof = LvArray::integerConversion< localIndex >( stack.bEqnRowIndices[ i ] ); + + if( dof < 0 || dof >= m_matrix.numRows() ) continue; + + // Is it necessary? Each row should be indepenedent + RAJA::atomicAdd< parallelDeviceAtomic >( &m_rhs[dof], stack.localRb[i] ); + + // Fill in matrix + m_matrix.template addToRowBinarySearchUnsorted< parallelDeviceAtomic >( dof, + stack.bColIndices, + stack.localAbtAtb[i], + numBdofs ); + + m_matrix.template addToRowBinarySearchUnsorted< parallelDeviceAtomic >( dof, + stack.dispColIndices, + stack.localAbtAtu[i], + numUdofs ); + } + + return 0.0; + } + +protected: + + arrayView2d< real64 const > const m_traction; +}; + +/// The factory used to construct the kernel. +using ALMSimultaneousFactory = finiteElement::InterfaceKernelFactory< ALMSimultaneous, + arrayView1d< globalIndex const > const, + arrayView1d< globalIndex const > const, + globalIndex const, + CRSMatrixView< real64, globalIndex const > const, + arrayView1d< real64 > const, + real64 const, + arrayView1d< localIndex const > const >; + +/** + * @brief A struct to compute the traction after nonlinear solve + */ +struct ComputeTractionSimultaneousKernel +{ + + /** + * @brief Launch the kernel function to comute traction + * @tparam POLICY the type of policy used in the kernel launch + * @param[in] size the size of the subregion + * @param[in] penalty the array containing the tangential penalty matrix + * @param[in] traction the array containing the current traction + * @param[in] dispJump the array containing the displacement jump + * @param[in] deltaDispJump the array containing the delta displacement jump + * @param[out] tractionNew the array containing the new traction + */ + template< typename POLICY > + static void + launch( localIndex const size, + arrayView2d< real64 const > const & penalty, + arrayView2d< real64 const > const & traction, + arrayView2d< real64 const > const & dispJump, + arrayView2d< real64 const > const & deltaDispJump, + arrayView2d< real64 > const & tractionNew ) + { + + forAll< POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const kfe ) + { + tractionNew[kfe][0] = traction[kfe][0] + penalty[kfe][0] * dispJump[kfe][0]; + tractionNew[kfe][1] = traction[kfe][1] + penalty[kfe][2] * deltaDispJump[kfe][1] + + penalty[kfe][4] * deltaDispJump[kfe][2]; + tractionNew[kfe][2] = traction[kfe][2] + penalty[kfe][3] * deltaDispJump[kfe][2] + + penalty[kfe][4] * deltaDispJump[kfe][1]; + } ); + } + +}; + +} // namespace SolidMechanicsALMKernels + +} // namespace geos + + +#endif /* GEOS_PHYSICSSOLVERS_CONTACT_SOLIDMECHANICSALMKERNELS_HPP_ */ diff --git a/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.cpp b/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.cpp index f37d5e8d05d..eb9fc81049e 100644 --- a/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.cpp +++ b/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.cpp @@ -21,9 +21,13 @@ #include "SolidMechanicsAugmentedLagrangianContact.hpp" #include "physicsSolvers/contact/SolidMechanicsALMKernels.hpp" +#include "physicsSolvers/contact/SolidMechanicsALMSimultaneousKernels.hpp" #include "physicsSolvers/contact/SolidMechanicsALMJumpUpdateKernels.hpp" #include "physicsSolvers/contact/SolidMechanicsALMBubbleKernels.hpp" +#include "constitutive/ConstitutiveManager.hpp" +#include "constitutive/contact/FrictionSelector.hpp" + namespace geos { @@ -95,9 +99,9 @@ void SolidMechanicsAugmentedLagrangianContact::registerDataOnMesh( dataRepositor subRegion.registerField< contact::oldDispJump >( this->getName() ). reference().resizeDimension< 1 >( 3 ); - // Register the displacement jump - subRegion.registerField< contact::penalty >( this->getName() ). - reference().resizeDimension< 1 >( 2 ); + // Register the penalty coefficients for the iterative procedure + subRegion.registerField< contact::iterativePenalty >( this->getName() ). + reference().resizeDimension< 1 >( 5 ); subRegion.registerWrapper< array1d< real64 > >( viewKeyStruct::normalTractionToleranceString() ). setPlotLevel( PlotLevel::NOPLOT ). @@ -114,6 +118,12 @@ void SolidMechanicsAugmentedLagrangianContact::registerDataOnMesh( dataRepositor setRegisteringObjects( this->getName()). setDescription( "An array that holds the sliding tolerance." ); + subRegion.registerWrapper< array2d< real64 > >( viewKeyStruct::dispJumpUpdPenaltyString() ). + setPlotLevel( PlotLevel::NOPLOT ). + setRegisteringObjects( this->getName()). + setDescription( "An array that stores the displacement jumps used to update the penalty coefficients." ). + reference().resizeDimension< 1 >( 3 ); + } ); } ); @@ -147,14 +157,6 @@ void SolidMechanicsAugmentedLagrangianContact::setupDofs( DomainPartition const solidMechanics::totalBubbleDisplacement::key(), DofManager::Connector::Elem ); - // The dofManager can not create the connection due to the coupling - // between totalDisplacement and totalBubbleDisplacement - // These connection are created using the two functions - // addCouplingNumNonzeros and addCouplingSparsityPattern - // dofManager.addCoupling( solidMechanics::totalDisplacement::key(), - // solidMechanics::totalBubbleDisplacement::key(), - // DofManager::Connector::Elem); - } void SolidMechanicsAugmentedLagrangianContact::setupSystem( DomainPartition & domain, @@ -168,9 +170,12 @@ void SolidMechanicsAugmentedLagrangianContact::setupSystem( DomainPartition & do GEOS_MARK_FUNCTION; GEOS_UNUSED_VAR( setSparsity ); - // Create the list of interface elements that have same type. + // Create the lists of interface elements that have same type. createFaceTypeList( domain ); + // Create the lists of interface elements that have same type and same fracture state. + updateStickSlipList( domain ); + // Create the list of cell elements that they are enriched with bubble functions. createBubbleCellList( domain ); @@ -236,9 +241,12 @@ void SolidMechanicsAugmentedLagrangianContact::implicitStepSetup( real64 const & SurfaceElementRegion & region = elemManager.getRegion< SurfaceElementRegion >( getUniqueFractureRegionName() ); FaceElementSubRegion & subRegion = region.getUniqueSubRegion< FaceElementSubRegion >(); - arrayView2d< real64 const > const & faceNormal = faceManager.faceNormal(); + arrayView2d< real64 const > const faceNormal = faceManager.faceNormal(); ArrayOfArraysView< localIndex const > const elemsToFaces = subRegion.faceList().toViewConst(); + arrayView2d< real64 > const incrBubbleDisp = + faceManager.getField< fields::solidMechanics::incrementalBubbleDisplacement >(); + arrayView3d< real64 > const rotationMatrix = subRegion.getField< fields::contact::rotationMatrix >().toView(); @@ -249,21 +257,67 @@ void SolidMechanicsAugmentedLagrangianContact::implicitStepSetup( real64 const & elemsToFaces, rotationMatrix ); + // Set the tollerances + computeTolerances( domain ); + + // Set array to update penalty coefficients + arrayView2d< real64 > const dispJumpUpdPenalty = + subRegion.getReference< array2d< real64 > >( viewKeyStruct::dispJumpUpdPenaltyString() ); + arrayView2d< real64 > const - penalty = subRegion.getField< fields::contact::penalty >().toView(); + iterativePenalty = subRegion.getField< fields::contact::iterativePenalty >().toView(); + arrayView1d< integer const > const fractureState = subRegion.getField< contact::fractureState >(); + + if( m_simultaneous ) + { + // Set the iterative penalty coefficients + forAll< parallelDevicePolicy<> >( subRegion.size(), + [=] + GEOS_HOST_DEVICE ( localIndex const k ) + { + if( fractureState[k] == contact::FractureState::Stick ) + { + iterativePenalty[k][2] = iterativePenalty[k][1]; + iterativePenalty[k][3] = iterativePenalty[k][1]; + iterativePenalty[k][4] = 0.0; + } + else + { + iterativePenalty[k][2] = 0.0; + iterativePenalty[k][3] = 0.0; + iterativePenalty[k][4] = 0.0; + } + } ); + } - // Set the penalty coefficients - // TODO: This is only temporary. The setting of penalty will be adaptive in the final version. - forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) + forAll< parallelDevicePolicy<> >( subRegion.size(), + [ = ] + GEOS_HOST_DEVICE ( localIndex const k ) { - penalty[k] [0] = 1.e+7; - penalty[k] [1] = 1.e+7; + LvArray::tensorOps::fill< 3 >( dispJumpUpdPenalty[k], 0.0 ); + localIndex const kf0 = elemsToFaces[k][0]; + localIndex const kf1 = elemsToFaces[k][1]; + LvArray::tensorOps::fill< 3 >( incrBubbleDisp[kf0], 0.0 ); + LvArray::tensorOps::fill< 3 >( incrBubbleDisp[kf1], 0.0 ); } ); + } ); + // Sync iterativePenalty + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & ) + { + FieldIdentifiers fieldsToBeSync; + + fieldsToBeSync.addElementFields( { contact::iterativePenalty::key() }, + { getUniqueFractureRegionName() } ); + + CommunicationTools::getInstance().synchronizeFields( fieldsToBeSync, + mesh, + domain.getNeighbors(), + true ); } ); - // Set the tollerances - computeTolerances( domain ); } void SolidMechanicsAugmentedLagrangianContact::assembleSystem( real64 const time, @@ -288,7 +342,6 @@ void SolidMechanicsAugmentedLagrangianContact::assembleSystem( real64 const time //ParallelMatrix parallel_matrix; //parallel_matrix.create( localMatrix.toViewConst(), dofManager.numLocalDofs(), MPI_COMM_GEOS ); //parallel_matrix.write("mech.mtx"); - //abort(); // Loop for assembling contributes from interface elements (Aut*eps^-1*Atu and Aub*eps^-1*Abu) forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const & meshName, @@ -308,31 +361,113 @@ void SolidMechanicsAugmentedLagrangianContact::assembleSystem( real64 const time string const & fractureRegionName = this->getUniqueFractureRegionName(); - forFiniteElementOnFractureSubRegions( meshName, [&] ( string const & finiteElementName, - arrayView1d< localIndex const > const & faceElementList ) + forFiniteElementOnStickFractureSubRegions( meshName, [&] ( string const &, + finiteElement::FiniteElementBase const & subRegionFE, + arrayView1d< localIndex const > const & faceElementList, + bool const ) { - finiteElement::FiniteElementBase & subRegionFE = *(m_faceTypeToFiniteElements[finiteElementName]); + if( m_simultaneous ) + { + solidMechanicsALMKernels::ALMSimultaneousFactory kernelFactory( dispDofNumber, + bubbleDofNumber, + dofManager.rankOffset(), + localMatrix, + localRhs, + dt, + faceElementList ); + + real64 maxTraction = finiteElement:: + interfaceBasedKernelApplication + < parallelDevicePolicy< >, + constitutive::CoulombFriction >( mesh, + fractureRegionName, + faceElementList, + subRegionFE, + viewKeyStruct::frictionLawNameString(), + kernelFactory ); + + GEOS_UNUSED_VAR( maxTraction ); - solidMechanicsALMKernels::ALMFactory kernelFactory( dispDofNumber, - bubbleDofNumber, - dofManager.rankOffset(), - localMatrix, - localRhs, - dt, - faceElementList ); + } + else + { + solidMechanicsALMKernels::ALMFactory kernelFactory( dispDofNumber, + bubbleDofNumber, + dofManager.rankOffset(), + localMatrix, + localRhs, + dt, + faceElementList, + m_symmetric ); + + real64 maxTraction = finiteElement:: + interfaceBasedKernelApplication + < parallelDevicePolicy< >, + constitutive::CoulombFriction >( mesh, + fractureRegionName, + faceElementList, + subRegionFE, + viewKeyStruct::frictionLawNameString(), + kernelFactory ); + + GEOS_UNUSED_VAR( maxTraction ); + } - real64 maxTraction = finiteElement:: - interfaceBasedKernelApplication - < parallelDevicePolicy< >, - constitutive::NullModel >( mesh, - fractureRegionName, - faceElementList, - subRegionFE, - "", - kernelFactory ); + } ); - GEOS_UNUSED_VAR( maxTraction ); + forFiniteElementOnSlipFractureSubRegions( meshName, [&] ( string const &, + finiteElement::FiniteElementBase const & subRegionFE, + arrayView1d< localIndex const > const & faceElementList, + bool const ) + { + + if( m_simultaneous ) + { + solidMechanicsALMKernels::ALMSimultaneousFactory kernelFactory( dispDofNumber, + bubbleDofNumber, + dofManager.rankOffset(), + localMatrix, + localRhs, + dt, + faceElementList ); + + real64 maxTraction = finiteElement:: + interfaceBasedKernelApplication + < parallelDevicePolicy< >, + constitutive::CoulombFriction >( mesh, + fractureRegionName, + faceElementList, + subRegionFE, + viewKeyStruct::frictionLawNameString(), + kernelFactory ); + + GEOS_UNUSED_VAR( maxTraction ); + + } + else + { + solidMechanicsALMKernels::ALMFactory kernelFactory( dispDofNumber, + bubbleDofNumber, + dofManager.rankOffset(), + localMatrix, + localRhs, + dt, + faceElementList, + m_symmetric ); + + real64 maxTraction = finiteElement:: + interfaceBasedKernelApplication + < parallelDevicePolicy< >, + constitutive::CoulombFriction >( mesh, + fractureRegionName, + faceElementList, + subRegionFE, + viewKeyStruct::frictionLawNameString(), + kernelFactory ); + + GEOS_UNUSED_VAR( maxTraction ); + } } ); @@ -384,7 +519,6 @@ void SolidMechanicsAugmentedLagrangianContact::implicitStepComplete( real64 cons DomainPartition & domain ) { - SolidMechanicsLagrangianFEM::implicitStepComplete( time_n, dt, domain ); forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, @@ -399,10 +533,16 @@ void SolidMechanicsAugmentedLagrangianContact::implicitStepComplete( real64 cons arrayView2d< real64 > const oldDispJump = subRegion.getField< contact::oldDispJump >(); arrayView2d< real64 > const deltaDispJump = subRegion.getField< contact::deltaDispJump >(); - forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const kfe ) + arrayView1d< integer const > const fractureState = subRegion.getField< contact::fractureState >(); + arrayView1d< integer > const oldFractureState = subRegion.getField< contact::oldFractureState >(); + + forAll< parallelDevicePolicy<> >( subRegion.size(), + [ = ] + GEOS_HOST_DEVICE ( localIndex const kfe ) { LvArray::tensorOps::fill< 3 >( deltaDispJump[kfe], 0.0 ); LvArray::tensorOps::copy< 3 >( oldDispJump[kfe], dispJump[kfe] ); + oldFractureState[kfe] = fractureState[kfe]; } ); } ); @@ -441,13 +581,15 @@ real64 SolidMechanicsAugmentedLagrangianContact::calculateResidualNorm( real64 c SurfaceElementRegion const & region = elemManager.getRegion< SurfaceElementRegion >( getUniqueFractureRegionName() ); FaceElementSubRegion const & subRegion = region.getUniqueSubRegion< FaceElementSubRegion >(); - arrayView1d< integer const > const & ghostRank = subRegion.ghostRank(); + arrayView1d< integer const > const ghostRank = subRegion.ghostRank(); ArrayOfArraysView< localIndex const > const elemsToFaces = subRegion.faceList().toViewConst(); arrayView1d< globalIndex const > const bubbleDofNumber = faceManager.getReference< globalIndex_array >( bubbleDofKey ); - forAll< parallelDevicePolicy<> >( subRegion.size(), [ elemsToFaces, localRhs, localSum, bubbleDofNumber, rankOffset, ghostRank] GEOS_HOST_DEVICE ( localIndex const kfe ) + forAll< parallelDevicePolicy<> >( subRegion.size(), + [ = ] + GEOS_HOST_DEVICE ( localIndex const kfe ) { if( ghostRank[kfe] < 0 ) @@ -550,12 +692,11 @@ void SolidMechanicsAugmentedLagrangianContact::applySystemSolution( DofManager c CRSMatrix< real64, globalIndex > const voidMatrix; array1d< real64 > const voidRhs; - forFiniteElementOnFractureSubRegions( meshName, [&] ( string const & finiteElementName, + forFiniteElementOnFractureSubRegions( meshName, [&] ( string const &, + finiteElement::FiniteElementBase const & subRegionFE, arrayView1d< localIndex const > const & faceElementList ) { - finiteElement::FiniteElementBase & subRegionFE = *(m_faceTypeToFiniteElements[finiteElementName]); - solidMechanicsALMKernels::ALMJumpUpdateFactory kernelFactory( dispDofNumber, bubbleDofNumber, dofManager.rankOffset(), @@ -611,10 +752,9 @@ bool SolidMechanicsAugmentedLagrangianContact::updateConfiguration( DomainPartit { GEOS_MARK_FUNCTION; - int hasConfigurationConverged = true; - int condConv = true; - // TODO: This function's design is temporary and intended solely for testing the stick mode. - // In the final version the parallelHostPolicy will be substitute with the parallelDevicePolicy<>. + array1d< int > condConv; + localIndex globalCondConv[5] = {0, 0, 0, 0, 0}; + array2d< real64 > traction_new; forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, @@ -627,12 +767,16 @@ bool SolidMechanicsAugmentedLagrangianContact::updateConfiguration( DomainPartit FaceElementSubRegion & subRegion ) { - arrayView1d< integer const > const & ghostRank = subRegion.ghostRank(); + string const & frictionLawName = subRegion.template getReference< string >( viewKeyStruct::frictionLawNameString() ); + FrictionBase const & frictionLaw = getConstitutiveModel< FrictionBase >( subRegion, frictionLawName ); + + arrayView1d< integer const > const ghostRank = subRegion.ghostRank(); arrayView2d< real64 const > const traction = subRegion.getField< contact::traction >(); arrayView2d< real64 const > const dispJump = subRegion.getField< contact::dispJump >(); arrayView2d< real64 const > const deltaDispJump = subRegion.getField< contact::deltaDispJump >(); - arrayView2d< real64 const > const penalty = subRegion.getField< contact::penalty >(); + arrayView2d< real64 const > const iterativePenalty = subRegion.getField< contact::iterativePenalty >(); + arrayView1d< integer const > const fractureState = subRegion.getField< contact::fractureState >(); arrayView1d< real64 const > const normalDisplacementTolerance = subRegion.getReference< array1d< real64 > >( viewKeyStruct::normalDisplacementToleranceString() ); @@ -641,95 +785,317 @@ bool SolidMechanicsAugmentedLagrangianContact::updateConfiguration( DomainPartit arrayView1d< real64 const > const & normalTractionTolerance = subRegion.getReference< array1d< real64 > >( viewKeyStruct::normalTractionToleranceString() ); + arrayView1d< real64 const > const area = subRegion.getElementArea().toViewConst(); + std::ptrdiff_t const sizes[ 2 ] = {subRegion.size(), 3}; traction_new.resize( 2, sizes ); arrayView2d< real64 > const traction_new_v = traction_new.toView(); - // TODO: Create a Kernel to update Traction - forAll< parallelHostPolicy >( subRegion.size(), - [ traction_new_v, traction, penalty, dispJump, deltaDispJump ] ( localIndex const kfe ) + condConv.resize( subRegion.size()); + arrayView1d< int > const condConv_v = condConv.toView(); + + // Update the traction field based on the displacement results from the nonlinear solve + constitutiveUpdatePassThru( frictionLaw, [&] ( auto & castedFrictionLaw ) { - real64 eps_N = penalty[kfe][0]; - real64 eps_T = penalty[kfe][1]; - traction_new_v[kfe][0] = traction[kfe][0] + eps_N * dispJump[kfe][0]; - traction_new_v[kfe][1] = traction[kfe][1] + eps_T * deltaDispJump[kfe][1]; - traction_new_v[kfe][2] = traction[kfe][2] + eps_T * deltaDispJump[kfe][2]; + using FrictionType = TYPEOFREF( castedFrictionLaw ); + typename FrictionType::KernelWrapper frictionWrapper = castedFrictionLaw.createKernelUpdates(); + + if( m_simultaneous ) + { + solidMechanicsALMKernels::ComputeTractionSimultaneousKernel:: + launch< parallelDevicePolicy<> >( subRegion.size(), + iterativePenalty, + traction, + dispJump, + deltaDispJump, + traction_new_v ); + } + else + { + solidMechanicsALMKernels::ComputeTractionKernel:: + launch< parallelDevicePolicy<> >( subRegion.size(), + frictionWrapper, + iterativePenalty, + traction, + dispJump, + deltaDispJump, + traction_new_v ); + } + } ); + + real64 const slidingCheckTolerance = m_slidingCheckTolerance; + + constitutiveUpdatePassThru( frictionLaw, [&] ( auto & castedFrictionLaw ) + { + using FrictionType = TYPEOFREF( castedFrictionLaw ); + typename FrictionType::KernelWrapper frictionWrapper = castedFrictionLaw.createKernelUpdates(); + + solidMechanicsALMKernels::ConstraintCheckKernel:: + launch< parallelDevicePolicy<> >( subRegion.size(), + frictionWrapper, + ghostRank, + traction_new_v, + dispJump, + deltaDispJump, + normalTractionTolerance, + normalDisplacementTolerance, + slidingTolerance, + slidingCheckTolerance, + area, + fractureState, + condConv_v ); } ); - forAll< parallelHostPolicy >( subRegion.size(), - [ &condConv, ghostRank, traction_new_v, normalTractionTolerance, normalDisplacementTolerance, slidingTolerance, deltaDispJump, dispJump ] ( localIndex const kfe ) + RAJA::ReduceSum< parallelDeviceReduce, localIndex > localSum[5] = + { RAJA::ReduceSum< parallelDeviceReduce, localIndex >( 0 ), + RAJA::ReduceSum< parallelDeviceReduce, localIndex >( 0 ), + RAJA::ReduceSum< parallelDeviceReduce, localIndex >( 0 ), + RAJA::ReduceSum< parallelDeviceReduce, localIndex >( 0 ), + RAJA::ReduceSum< parallelDeviceReduce, localIndex >( 0 ) }; + forAll< parallelDevicePolicy<> >( subRegion.size(), [ = ] GEOS_HOST_DEVICE ( localIndex const kfe ) { if( ghostRank[kfe] < 0 ) { - if( traction_new_v[kfe][0] >= normalTractionTolerance[kfe] ) - { - GEOS_ERROR( "Unsuported open mode! Only stick mode is supported with ALM" ); - } - else + localSum[condConv_v[kfe]] += 1; + } + } ); + + localIndex const localConvCond[5] = { static_cast< localIndex >( localSum[0].get()), + static_cast< localIndex >( localSum[1].get()), + static_cast< localIndex >( localSum[2].get()), + static_cast< localIndex >( localSum[3].get()), + static_cast< localIndex >( localSum[4].get()) }; + + int const rank = MpiWrapper::commRank( MPI_COMM_GEOS ); + int const numRanks = MpiWrapper::commSize( MPI_COMM_GEOS ); + array1d< localIndex > globalValues( numRanks * 5 ); + + // Everything is done on rank 0 + MpiWrapper::gather( localConvCond, + 5, + globalValues.data(), + 5, + 0, + MPI_COMM_GEOS ); + + if( rank==0 ) + { + for( int r=0; r normalDisplacementTolerance[kfe] ) - { - //GEOS_LOG_LEVEL(2, - //GEOS_FMT( " Element: {}, Stick mode and g_n > tol1 => compenetration, g_n: {:4.2e} tol1: {:4.2e}", - //kfe,std::abs(dispJump[kfe][0]), normalDisplacementTolerance[kfe] )) - condConv = false; - } - if( deltaDisp > slidingTolerance[kfe] ) - { - //GEOS_LOG_LEVEL(2, - //GEOS_FMT( " Element: {}, Stick and dg_t > tol2, dg_t: {:4.2e} tol2: {:4.2e}", - //kfe, deltaDisp, slidingTolerance[kfe] )) - condConv = false; - } + globalCondConv[i] += globalValues[r*5+i]; } } - } ); + } + + MpiWrapper::bcast( globalCondConv, 5, 0, MPI_COMM_GEOS ); } ); } ); - forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, - MeshLevel & mesh, - arrayView1d< string const > const & regionNames ) + localIndex totCondNotConv = 0; + for( int i=0; i<4; ++i ) { - ElementRegionManager & elemManager = mesh.getElemManager(); + totCondNotConv+=globalCondConv[i+1]; + } - elemManager.forElementSubRegions< FaceElementSubRegion >( regionNames, [&]( localIndex const, - FaceElementSubRegion & subRegion ) + int hasConfigurationConvergedGlobally = (totCondNotConv == 0) ? true : false; + + GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( " ALM convergence summary:" + " converged: {:6} | stick & gn>0: {:6} | compenetration: {:6} | stick & gt>lim: {:6} | tau>tauLim: {:6}\n", + globalCondConv[0], globalCondConv[1], globalCondConv[2], + globalCondConv[3], globalCondConv[4] )); + + if( hasConfigurationConvergedGlobally ) + { + + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & regionNames ) { + ElementRegionManager & elemManager = mesh.getElemManager(); - arrayView2d< real64 > const traction_new_v = traction_new.toView(); - arrayView2d< real64 > const traction = subRegion.getField< contact::traction >(); - forAll< parallelHostPolicy >( subRegion.size(), [ traction, traction_new_v ] ( localIndex const kfe ) + elemManager.forElementSubRegions< FaceElementSubRegion >( regionNames, [&]( localIndex const, + FaceElementSubRegion & subRegion ) { - traction[kfe][0] = traction_new_v( kfe, 0 ); - traction[kfe][1] = traction_new_v( kfe, 1 ); - traction[kfe][2] = traction_new_v( kfe, 2 ); + + arrayView2d< real64 > const traction_new_v = traction_new.toView(); + arrayView2d< real64 > const traction = subRegion.getField< contact::traction >(); + + forAll< parallelDevicePolicy<> >( subRegion.size(), [ = ] GEOS_HOST_DEVICE ( localIndex const kfe ) + { + LvArray::tensorOps::copy< 3 >( traction[kfe], traction_new_v[kfe] ); + } ); } ); } ); - } ); - - if( !condConv ) + } + else { - hasConfigurationConverged = false; + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & regionNames ) + { + ElementRegionManager & elemManager = mesh.getElemManager(); + + elemManager.forElementSubRegions< FaceElementSubRegion >( regionNames, [&]( localIndex const, + FaceElementSubRegion & subRegion ) + { + + string const & frictionLawName = subRegion.template getReference< string >( viewKeyStruct::frictionLawNameString() ); + FrictionBase const & frictionLaw = getConstitutiveModel< FrictionBase >( subRegion, frictionLawName ); + + arrayView2d< real64 > const traction = subRegion.getField< contact::traction >(); + + arrayView1d< real64 const > const normalTractionTolerance = + subRegion.getReference< array1d< real64 > >( viewKeyStruct::normalTractionToleranceString() ); + + arrayView2d< real64 > const iterativePenalty = subRegion.getField< fields::contact::iterativePenalty >().toView(); + + arrayView2d< real64 > const dispJumpUpdPenalty = + subRegion.getReference< array2d< real64 > >( viewKeyStruct::dispJumpUpdPenaltyString() ); + + arrayView1d< integer > const fractureState = subRegion.getField< contact::fractureState >(); + + arrayView2d< real64 const > const dispJump = subRegion.getField< contact::dispJump >(); + + arrayView2d< real64 const > const oldDispJump = subRegion.getField< contact::oldDispJump >(); + + arrayView2d< real64 const > const deltaDispJump = subRegion.getField< contact::deltaDispJump >(); + + constitutiveUpdatePassThru( frictionLaw, [&] ( auto & castedFrictionLaw ) + { + using FrictionType = TYPEOFREF( castedFrictionLaw ); + typename FrictionType::KernelWrapper frictionWrapper = castedFrictionLaw.createKernelUpdates(); + + solidMechanicsALMKernels::UpdateStateKernel:: + launch< parallelDevicePolicy<> >( subRegion.size(), + frictionWrapper, + oldDispJump, + dispJump, + iterativePenalty, + m_symmetric, + normalTractionTolerance, + traction, + fractureState ); + } ); + + forAll< parallelDevicePolicy<> >( subRegion.size(), [ = ] + GEOS_HOST_DEVICE ( localIndex const kfe ) + { + dispJumpUpdPenalty[kfe][0] = dispJump[kfe][0]; + dispJumpUpdPenalty[kfe][1] = deltaDispJump[kfe][1]; + dispJumpUpdPenalty[kfe][2] = deltaDispJump[kfe][2]; + } ); + } ); + } ); } // Need to synchronize the fracture state due to the use will be made of in AssemblyStabilization synchronizeFractureState( domain ); - // Compute if globally the fracture state has changed - int hasConfigurationConvergedGlobally; - MpiWrapper::allReduce( &hasConfigurationConverged, - &hasConfigurationConvergedGlobally, - 1, - MPI_LAND, - MPI_COMM_GEOS ); + // Update lists of stick and slip elements + if( !hasConfigurationConvergedGlobally ) + { + updateStickSlipList( domain ); + } return hasConfigurationConvergedGlobally; } +void SolidMechanicsAugmentedLagrangianContact::updateStickSlipList( DomainPartition const & domain ) +{ + + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const & meshName, + MeshLevel const & mesh, + arrayView1d< string const > const & ) + + { + + ElementRegionManager const & elemManager = mesh.getElemManager(); + SurfaceElementRegion const & region = elemManager.getRegion< SurfaceElementRegion >( getUniqueFractureRegionName() ); + FaceElementSubRegion const & subRegion = region.getUniqueSubRegion< FaceElementSubRegion >(); + + arrayView1d< integer const > const fractureState = subRegion.getField< contact::fractureState >(); + + forFiniteElementOnFractureSubRegions( meshName, [&] ( string const & finiteElementName, + finiteElement::FiniteElementBase const &, + arrayView1d< localIndex const > const & faceElementList ) + { + + array1d< localIndex > keys( subRegion.size()); + array1d< localIndex > vals( subRegion.size()); + array1d< localIndex > stickList; + array1d< localIndex > slipList; + RAJA::ReduceSum< ReducePolicy< parallelDevicePolicy<> >, localIndex > nStick_r( 0 ); + RAJA::ReduceSum< ReducePolicy< parallelDevicePolicy<> >, localIndex > nSlip_r( 0 ); + + arrayView1d< localIndex > const keys_v = keys.toView(); + arrayView1d< localIndex > const vals_v = vals.toView(); + forAll< parallelDevicePolicy<> >( faceElementList.size(), + [ = ] + GEOS_HOST_DEVICE ( localIndex const kfe ) + { + + localIndex const faceIndex = faceElementList[kfe]; + if( fractureState[faceIndex] == contact::FractureState::Stick ) + { + keys_v[kfe]=0; + vals_v[kfe]=faceIndex; + nStick_r += 1; + } + else if(( fractureState[faceIndex] == contact::FractureState::Slip ) || + (fractureState[faceIndex] == contact::FractureState::NewSlip)) + { + keys_v[kfe]=1; + vals_v[kfe]=faceIndex; + nSlip_r += 1; + } + else + { + keys_v[kfe] = 2; + } + + } ); + + localIndex nStick = static_cast< localIndex >(nStick_r.get()); + localIndex nSlip = static_cast< localIndex >(nSlip_r.get()); + + // Sort vals according to keys to ensure that + // elements of the same type are adjacent in the vals list. + // This arrangement allows for efficient copying into the container + // by leveraging parallelism. + RAJA::sort_pairs< parallelDevicePolicy<> >( keys_v, vals_v ); + + stickList.resize( nStick ); + slipList.resize( nSlip ); + arrayView1d< localIndex > const stickList_v = stickList.toView(); + arrayView1d< localIndex > const slipList_v = slipList.toView(); + + forAll< parallelDevicePolicy<> >( nStick, [ = ] + GEOS_HOST_DEVICE ( localIndex const kfe ) + { + stickList_v[kfe] = vals_v[kfe]; + } ); + + forAll< parallelDevicePolicy<> >( nSlip, [ = ] + GEOS_HOST_DEVICE ( localIndex const kfe ) + { + slipList_v[kfe] = vals_v[nStick+kfe]; + } ); + + this->m_faceTypesToFaceElementsStick[meshName][finiteElementName] = stickList; + this->m_faceTypesToFaceElementsSlip[meshName][finiteElementName] = slipList; + + GEOS_LOG_LEVEL( 2, + GEOS_FMT( "# stick elements: {}, # slip elements: {}", nStick, nSlip )) + } ); + } ); + +} + void SolidMechanicsAugmentedLagrangianContact::createFaceTypeList( DomainPartition const & domain ) { @@ -756,7 +1122,8 @@ void SolidMechanicsAugmentedLagrangianContact::createFaceTypeList( DomainPartiti arrayView1d< localIndex > const vals_v = vals.toView(); // Determine the size of the lists and generate the vector keys and vals for parallel indexing into lists. // (With RAJA, parallelizing this operation seems the most viable approach.) - forAll< parallelDevicePolicy<> >( subRegion.size(), [ nTri_r, nQuad_r, faceToNodeMap, keys_v, vals_v ] GEOS_HOST_DEVICE ( localIndex const kfe ) + forAll< parallelDevicePolicy<> >( subRegion.size(), + [ = ] GEOS_HOST_DEVICE ( localIndex const kfe ) { localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( kfe ); @@ -792,12 +1159,12 @@ void SolidMechanicsAugmentedLagrangianContact::createFaceTypeList( DomainPartiti arrayView1d< localIndex > const quadList_v = quadList.toView(); arrayView1d< localIndex > const triList_v = triList.toView(); - forAll< parallelDevicePolicy<> >( nTri, [triList_v, vals_v] GEOS_HOST_DEVICE ( localIndex const kfe ) + forAll< parallelDevicePolicy<> >( nTri, [ = ] GEOS_HOST_DEVICE ( localIndex const kfe ) { triList_v[kfe] = vals_v[kfe]; } ); - forAll< parallelDevicePolicy<> >( nQuad, [quadList_v, vals_v, nTri] GEOS_HOST_DEVICE ( localIndex const kfe ) + forAll< parallelDevicePolicy<> >( nQuad, [ = ] GEOS_HOST_DEVICE ( localIndex const kfe ) { quadList_v[kfe] = vals_v[nTri+kfe]; } ); @@ -819,6 +1186,7 @@ void SolidMechanicsAugmentedLagrangianContact::createBubbleCellList( DomainParti SurfaceElementRegion const & region = elemManager.getRegion< SurfaceElementRegion >( getUniqueFractureRegionName() ); FaceElementSubRegion const & subRegion = region.getUniqueSubRegion< FaceElementSubRegion >(); + // Array to store face indexes array1d< localIndex > tmpSpace( 2*subRegion.size()); SortedArray< localIndex > faceIdList; @@ -827,7 +1195,7 @@ void SolidMechanicsAugmentedLagrangianContact::createBubbleCellList( DomainParti { ArrayOfArraysView< localIndex const > const elemsToFaces = subRegion.faceList().toViewConst(); - forAll< parallelDevicePolicy<> >( subRegion.size(), [ tmpSpace_v, elemsToFaces ] GEOS_HOST_DEVICE ( localIndex const kfe ) + forAll< parallelDevicePolicy<> >( subRegion.size(), [ = ] GEOS_HOST_DEVICE ( localIndex const kfe ) { localIndex const kf0 = elemsToFaces[kfe][0], kf1 = elemsToFaces[kfe][1]; @@ -862,7 +1230,7 @@ void SolidMechanicsAugmentedLagrangianContact::createBubbleCellList( DomainParti SortedArrayView< localIndex const > const faceIdList_v = faceIdList.toViewConst(); forAll< parallelDevicePolicy<> >( cellElementSubRegion.size(), - [ keys_v, vals_v, perms_v, localFaceIds_v, nBubElems_r, elemsToFaces, faceIdList_v ] + [ = ] GEOS_HOST_DEVICE ( localIndex const kfe ) { for( int i=0; i < elemsToFaces.size( 1 ); ++i ) @@ -896,18 +1264,18 @@ void SolidMechanicsAugmentedLagrangianContact::createBubbleCellList( DomainParti arrayView1d< localIndex > const bubbleElemsList_v = bubbleElemsList.toView(); - forAll< parallelDevicePolicy<> >( n_max, [ keys_v, vals_v, perms_v ] GEOS_HOST_DEVICE ( localIndex const k ) + forAll< parallelDevicePolicy<> >( n_max, [ = ] GEOS_HOST_DEVICE ( localIndex const k ) { keys_v[k] = vals_v[perms_v[k]]; } ); - forAll< parallelDevicePolicy<> >( nBubElems, [ bubbleElemsList_v, keys_v ] GEOS_HOST_DEVICE ( localIndex const k ) + forAll< parallelDevicePolicy<> >( nBubElems, [ = ] GEOS_HOST_DEVICE ( localIndex const k ) { bubbleElemsList_v[k] = keys_v[k]; } ); cellElementSubRegion.setBubbleElementsList( bubbleElemsList.toViewConst()); - forAll< parallelDevicePolicy<> >( n_max, [ keys_v, localFaceIds_v, perms_v ] GEOS_HOST_DEVICE ( localIndex const k ) + forAll< parallelDevicePolicy<> >( n_max, [ = ] GEOS_HOST_DEVICE ( localIndex const k ) { keys_v[k] = localFaceIds_v[perms_v[k]]; } ); @@ -918,7 +1286,7 @@ void SolidMechanicsAugmentedLagrangianContact::createBubbleCellList( DomainParti arrayView2d< localIndex > const faceElemsList_v = faceElemsList.toView(); forAll< parallelDevicePolicy<> >( nBubElems, - [ faceElemsList_v, keys_v, elemsToFaces, bubbleElemsList_v ] + [ = ] GEOS_HOST_DEVICE ( localIndex const k ) { localIndex const kfe = bubbleElemsList_v[k]; @@ -1193,7 +1561,6 @@ void SolidMechanicsAugmentedLagrangianContact::addCouplingSparsityPattern( Domai } -// TODO: Is it possible to define this method once? Similar to SolidMechanicsLagrangeContact void SolidMechanicsAugmentedLagrangianContact::computeTolerances( DomainPartition & domain ) const { GEOS_MARK_FUNCTION; @@ -1208,16 +1575,16 @@ void SolidMechanicsAugmentedLagrangianContact::computeTolerances( DomainPartitio // Get the "face to element" map (valid for the entire mesh) FaceManager::ElemMapType const & faceToElem = faceManager.toElementRelation(); - arrayView2d< localIndex const > const & faceToElemRegion = faceToElem.m_toElementRegion; - arrayView2d< localIndex const > const & faceToElemSubRegion = faceToElem.m_toElementSubRegion; - arrayView2d< localIndex const > const & faceToElemIndex = faceToElem.m_toElementIndex; + arrayView2d< localIndex const > const faceToElemRegion = faceToElem.m_toElementRegion; + arrayView2d< localIndex const > const faceToElemSubRegion = faceToElem.m_toElementSubRegion; + arrayView2d< localIndex const > const faceToElemIndex = faceToElem.m_toElementIndex; // Get the volume for all elements ElementRegionManager::ElementViewAccessor< arrayView1d< real64 const > > const elemVolume = elemManager.constructViewAccessor< array1d< real64 >, arrayView1d< real64 const > >( ElementSubRegionBase::viewKeyStruct::elementVolumeString() ); // Get the coordinates for all nodes - arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition = nodeManager.referencePosition(); + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const nodePosition = nodeManager.referencePosition(); // Bulk modulus accessor ElementRegionManager::ElementViewAccessor< arrayView1d< real64 const > > const bulkModulus = @@ -1235,20 +1602,25 @@ void SolidMechanicsAugmentedLagrangianContact::computeTolerances( DomainPartitio { if( subRegion.hasField< contact::traction >() ) { - arrayView1d< integer const > const & ghostRank = subRegion.ghostRank(); - arrayView1d< real64 const > const & faceArea = subRegion.getElementArea().toViewConst(); - arrayView3d< real64 const > const & faceRotationMatrix = subRegion.getField< fields::contact::rotationMatrix >().toView(); - ArrayOfArraysView< localIndex const > const & elemsToFaces = subRegion.faceList().toViewConst(); + arrayView1d< real64 const > const faceArea = subRegion.getElementArea().toViewConst(); + arrayView3d< real64 const > const faceRotationMatrix = subRegion.getField< fields::contact::rotationMatrix >().toView(); + ArrayOfArraysView< localIndex const > const elemsToFaces = subRegion.faceList().toViewConst(); - arrayView1d< real64 > const & normalTractionTolerance = + arrayView1d< real64 > const normalTractionTolerance = subRegion.getReference< array1d< real64 > >( viewKeyStruct::normalTractionToleranceString() ); - arrayView1d< real64 > const & normalDisplacementTolerance = + arrayView1d< real64 > const normalDisplacementTolerance = subRegion.getReference< array1d< real64 > >( viewKeyStruct::normalDisplacementToleranceString() ); - arrayView1d< real64 > const & slidingTolerance = + arrayView1d< real64 > const slidingTolerance = subRegion.getReference< array1d< real64 > >( viewKeyStruct::slidingToleranceString() ); + arrayView2d< real64 > const + iterativePenalty = subRegion.getField< fields::contact::iterativePenalty >().toView(); + + arrayView1d< integer const > const ghostRank = subRegion.ghostRank(); + forAll< parallelHostPolicy >( subRegion.size(), [=] ( localIndex const kfe ) { + if( ghostRank[kfe] < 0 ) { real64 const area = faceArea[kfe]; @@ -1330,10 +1702,13 @@ void SolidMechanicsAugmentedLagrangianContact::computeTolerances( DomainPartitio LvArray::tensorOps::scale< 3, 3 >( rotatedInvStiffApprox, area ); // Finally, compute tolerances for the given fracture element - normalDisplacementTolerance[kfe] = rotatedInvStiffApprox[ 0 ][ 0 ] * averageYoungModulus / 2.e+8; - slidingTolerance[kfe] = sqrt( rotatedInvStiffApprox[ 1 ][ 1 ] * rotatedInvStiffApprox[ 1 ][ 1 ] + - rotatedInvStiffApprox[ 2 ][ 2 ] * rotatedInvStiffApprox[ 2 ][ 2 ] ) * averageYoungModulus / 2.e+8; - normalTractionTolerance[kfe] = 1.0 / 2.0 * averageConstrainedModulus / averageBoxSize0 * normalDisplacementTolerance[kfe]; + normalDisplacementTolerance[kfe] = rotatedInvStiffApprox[ 0 ][ 0 ] * averageYoungModulus / 2.e+7; + slidingTolerance[kfe] = sqrt( pow( rotatedInvStiffApprox[ 1 ][ 1 ], 2 ) + + pow( rotatedInvStiffApprox[ 2 ][ 2 ], 2 )) * averageYoungModulus / 2.e+5; + normalTractionTolerance[kfe] = (1.0 / 2.0) * (averageConstrainedModulus / averageBoxSize0) * + (normalDisplacementTolerance[kfe]); + iterativePenalty[kfe][0] = 1e+01*averageConstrainedModulus/(averageBoxSize0*area); + iterativePenalty[kfe][1] = 1e-01*averageConstrainedModulus/(averageBoxSize0*area); } } ); } diff --git a/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.hpp b/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.hpp index c11bd493159..54ffa40bab7 100644 --- a/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.hpp +++ b/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.hpp @@ -108,11 +108,76 @@ class SolidMechanicsAugmentedLagrangianContact : public ContactSolverBase for( const auto & [finiteElementName, faceElementList] : faceTypesToFaceElements ) { arrayView1d< localIndex const > const faceElemList = faceElementList.toViewConst(); - lambda( finiteElementName, faceElemList ); + + finiteElement::FiniteElementBase const & subRegionFE = *(m_faceTypeToFiniteElements.at( finiteElementName )); + + lambda( finiteElementName, subRegionFE, faceElemList ); } } + /** + * @brief Loop over the finite element type on the stick fracture subregions of meshName and apply callback. + * @tparam LAMBDA The callback function type + * @param meshName The mesh name. + * @param lambda The callback function. Take the finite element type name and + * the list of face element of the same type. + */ + template< typename LAMBDA > + void forFiniteElementOnStickFractureSubRegions( string const & meshName, LAMBDA && lambda ) const + { + + bool const isStickState = true; + + std::map< string, array1d< localIndex > > const & + faceTypesToFaceElements = m_faceTypesToFaceElementsStick.at( meshName ); + + for( const auto & [finiteElementName, faceElementList] : faceTypesToFaceElements ) + { + arrayView1d< localIndex const > const faceElemList = faceElementList.toViewConst(); + + finiteElement::FiniteElementBase const & subRegionFE = *(m_faceTypeToFiniteElements.at( finiteElementName )); + + lambda( finiteElementName, subRegionFE, faceElemList, isStickState ); + } + + } + + /** + * @brief Loop over the finite element type on the slip fracture subregions of meshName and apply callback. + * @tparam LAMBDA The callback function type + * @param meshName The mesh name. + * @param lambda The callback function. Take the finite element type name and + * the list of face element of the same type. + */ + template< typename LAMBDA > + void forFiniteElementOnSlipFractureSubRegions( string const & meshName, LAMBDA && lambda ) const + { + + bool const isStickState = false; + + std::map< string, array1d< localIndex > > const & + faceTypesToFaceElements = m_faceTypesToFaceElementsSlip.at( meshName ); + + for( const auto & [finiteElementName, faceElementList] : faceTypesToFaceElements ) + { + arrayView1d< localIndex const > const faceElemList = faceElementList.toViewConst(); + + finiteElement::FiniteElementBase const & subRegionFE = *(m_faceTypeToFiniteElements.at( finiteElementName )); + + lambda( finiteElementName, subRegionFE, faceElemList, isStickState ); + } + + } + + /** + * @brief Create the list of finite elements of the same type + * for each FaceElementSubRegion (Triangle or Quadrilateral) + * and of the same fracture state (Stick or Slip). + * @param domain The physical domain object + */ + void updateStickSlipList( DomainPartition const & domain ); + /** * @brief Create the list of finite elements of the same type * for each FaceElementSubRegion (Triangle or Quadrilateral). @@ -155,6 +220,12 @@ class SolidMechanicsAugmentedLagrangianContact : public ContactSolverBase /// Finite element type to face element index map std::map< string, std::map< string, array1d< localIndex > > > m_faceTypesToFaceElements; + /// Finite element type to face element index map (stick mode) + std::map< string, std::map< string, array1d< localIndex > > > m_faceTypesToFaceElementsStick; + + /// Finite element type to face element index map (slip mode) + std::map< string, std::map< string, array1d< localIndex > > > m_faceTypesToFaceElementsSlip; + /// Finite element type to finite element object map std::map< string, std::unique_ptr< geos::finiteElement::FiniteElementBase > > m_faceTypeToFiniteElements; @@ -166,8 +237,21 @@ class SolidMechanicsAugmentedLagrangianContact : public ContactSolverBase constexpr static char const * normalTractionToleranceString() { return "normalTractionTolerance"; } constexpr static char const * slidingToleranceString() { return "slidingTolerance"; } + + constexpr static char const * dispJumpUpdPenaltyString() { return "dispJumpUpdPenalty"; } }; + /// Tolerance for the sliding check: the tangential traction must exceed (1 + m_slidingCheckTolerance) * t_lim to activate the sliding + /// condition + real64 const m_slidingCheckTolerance = 0.05; + + /// Flag to update the Lagrange multiplier at each Newton iteration (true), or only after the Newton loop has converged (false) + bool m_simultaneous = true; + + /// Flag to neglect the non-symmetric contribution in the tangential matrix, i.e., the derivative of tangential traction with respect to + /// the normal displacement is neglected + bool m_symmetric = true; + }; } /* namespace geos */ diff --git a/src/coreComponents/physicsSolvers/contact/SolidMechanicsEmbeddedFractures.cpp b/src/coreComponents/physicsSolvers/contact/SolidMechanicsEmbeddedFractures.cpp index df8976ed422..90f33a650e0 100644 --- a/src/coreComponents/physicsSolvers/contact/SolidMechanicsEmbeddedFractures.cpp +++ b/src/coreComponents/physicsSolvers/contact/SolidMechanicsEmbeddedFractures.cpp @@ -740,7 +740,7 @@ void SolidMechanicsEmbeddedFractures::updateState( DomainPartition & domain ) constitutiveUpdatePassThru( frictionLaw, [&] ( auto & castedFrictionLaw ) { using FrictionType = TYPEOFREF( castedFrictionLaw ); - typename FrictionType::KernelWrapper frictionWrapper = castedFrictionLaw.createKernelWrapper(); + typename FrictionType::KernelWrapper frictionWrapper = castedFrictionLaw.createKernelUpdates(); solidMechanicsEFEMKernels::StateUpdateKernel:: launch< parallelDevicePolicy<> >( subRegion.size(), @@ -776,7 +776,7 @@ bool SolidMechanicsEmbeddedFractures::updateConfiguration( DomainPartition & dom constitutiveUpdatePassThru( frictionLaw, [&] ( auto & castedFrictionLaw ) { using FrictionType = TYPEOFREF( castedFrictionLaw ); - typename FrictionType::KernelWrapper frictionWrapper = castedFrictionLaw.createKernelWrapper(); + typename FrictionType::KernelWrapper frictionWrapper = castedFrictionLaw.createKernelUpdates(); RAJA::ReduceMin< parallelHostReduce, integer > checkActiveSetSub( 1 ); diff --git a/src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContact.cpp b/src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContact.cpp index e943728e7f0..2b66ae59da0 100644 --- a/src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContact.cpp +++ b/src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContact.cpp @@ -1476,7 +1476,7 @@ void SolidMechanicsLagrangeContact:: constitutiveUpdatePassThru( frictionLaw, [&] ( auto & castedFrictionLaw ) { using FrictionType = TYPEOFREF( castedFrictionLaw ); - typename FrictionType::KernelWrapper frictionWrapper = castedFrictionLaw.createKernelWrapper(); + typename FrictionType::KernelWrapper frictionWrapper = castedFrictionLaw.createKernelUpdates(); forAll< parallelHostPolicy >( subRegion.size(), [=] ( localIndex const kfe ) { @@ -2250,7 +2250,7 @@ bool SolidMechanicsLagrangeContact::updateConfiguration( DomainPartition & domai constitutiveUpdatePassThru( frictionLaw, [&] ( auto & castedFrictionLaw ) { using FrictionType = TYPEOFREF( castedFrictionLaw ); - typename FrictionType::KernelWrapper frictionWrapper = castedFrictionLaw.createKernelWrapper(); + typename FrictionType::KernelWrapper frictionWrapper = castedFrictionLaw.createKernelUpdates(); forAll< parallelHostPolicy >( subRegion.size(), [=] ( localIndex const kfe ) { diff --git a/src/docs/sphinx/advancedExamples/validationStudies/faultMechanics/intersectFrac/Example.rst b/src/docs/sphinx/advancedExamples/validationStudies/faultMechanics/intersectFrac/Example.rst index a1019bac817..6714806e81f 100644 --- a/src/docs/sphinx/advancedExamples/validationStudies/faultMechanics/intersectFrac/Example.rst +++ b/src/docs/sphinx/advancedExamples/validationStudies/faultMechanics/intersectFrac/Example.rst @@ -14,11 +14,15 @@ method `(Phan et al., 2003) :end-before: @@ -78,7 +82,7 @@ with one group of cell blocks named here ``cb1``. Refinement is necessary to conform with the fracture geometry specified in the ``Geometry`` section. -.. literalinclude:: ../../../../../../../inputFiles/lagrangianContactMechanics/ContactMechanics_TFrac_base.xml +.. literalinclude:: ../../../../../../../inputFiles/lagrangianContactMechanics/TFrac_base.xml :language: xml :start-after: :end-before: @@ -124,7 +128,7 @@ A homogeneous and isotropic domain with one solid material is assumed, and its m Fracture surface slippage is assumed to be governed by the Coulomb failure criterion. The contact constitutive behavior is named ``fractureMaterial`` in the ``Coulomb`` block, where cohesion ``cohesion="0.0"`` and friction coefficient ``frictionCoefficient="0.577350269"`` are specified. -.. literalinclude:: ../../../../../../../inputFiles/lagrangianContactMechanics/ContactMechanics_TFrac_base.xml +.. literalinclude:: ../../../../../../../inputFiles/lagrangianContactMechanics/TFrac_base.xml :language: xml :start-after: :end-before: @@ -145,7 +149,7 @@ In the ``Tasks`` section, ``PackCollection`` tasks are defined to collect time h Either the entire field or specified named sets of indices in the field can be collected. In this example, ``tractionCollection`` and ``displacementJumpCollection`` tasks are specified to output the local traction ``fieldName="traction"`` and relative displacement ``fieldName="displacementJump"`` on the fracture surface. -.. literalinclude:: ../../../../../../../inputFiles/lagrangianContactMechanics/ContactMechanics_TFrac_base.xml +.. literalinclude:: ../../../../../../../inputFiles/lagrangianContactMechanics/TFrac_base.xml :language: xml :start-after: :end-before: @@ -171,7 +175,7 @@ The remaining parts of the outer boundaries are subjected to roller constraints. These boundary conditions are set up through the ``FieldSpecifications`` section. -.. literalinclude:: ../../../../../../../inputFiles/lagrangianContactMechanics/ContactMechanics_TFrac_base.xml +.. literalinclude:: ../../../../../../../inputFiles/lagrangianContactMechanics/TFrac_base.xml :language: xml :start-after: :end-before: diff --git a/src/docs/sphinx/advancedExamples/validationStudies/faultMechanics/intersectFrac/intersectFracFigure.py b/src/docs/sphinx/advancedExamples/validationStudies/faultMechanics/intersectFrac/intersectFracFigure.py index 2aa2ffa8fdf..ac6f796ed04 100644 --- a/src/docs/sphinx/advancedExamples/validationStudies/faultMechanics/intersectFrac/intersectFracFigure.py +++ b/src/docs/sphinx/advancedExamples/validationStudies/faultMechanics/intersectFrac/intersectFracFigure.py @@ -105,7 +105,7 @@ def main(): geosDir = args.geosDir hdf5File1Path = outputDir + "/traction_history.hdf5" hdf5File2Path = outputDir + "/displacementJump_history.hdf5" - xmlFilePath = geosDir + "/inputFiles/lagrangianContactMechanics/ContactMechanics_TFrac_base.xml" + xmlFilePath = geosDir + "/inputFiles/lagrangianContactMechanics/TFrac_base.xml" # Read HDF5 # Global Coordinate of Fracture Element Center diff --git a/src/docs/sphinx/advancedExamples/validationStudies/faultMechanics/singleFracCompression/Example.rst b/src/docs/sphinx/advancedExamples/validationStudies/faultMechanics/singleFracCompression/Example.rst index c34a6cf30cb..4b29541ce9d 100644 --- a/src/docs/sphinx/advancedExamples/validationStudies/faultMechanics/singleFracCompression/Example.rst +++ b/src/docs/sphinx/advancedExamples/validationStudies/faultMechanics/singleFracCompression/Example.rst @@ -13,11 +13,15 @@ In this example, a single fracture is simulated using a Lagrange contact model i **Input file** -Everything required is contained within two GEOS input files and one mesh file located at: +Everything required is contained within these GEOS input files and one mesh file located at: .. code-block:: console - inputFiles/lagrangianContactMechanics/ContactMechanics_SingleFracCompression_base.xml + inputFiles/lagrangianContactMechanics/SingleFracCompression_base.xml + +.. code-block:: console + + inputFiles/lagrangianContactMechanics/SingleFracCompression_benchmark.xml .. code-block:: console @@ -79,7 +83,7 @@ The syntax to import external meshes is simple: in the XML file, the mesh file ``crackInPlane_benchmark.vtu`` is included with its relative or absolute path to the location of the GEOS XML file and a user-specified label (here ``CubeHex``) is given to the mesh object. This unstructured mesh contains quadrilaterals elements and interface elements. Refinement is performed to conform with the fracture geometry specified in the ``Geometry`` section. -.. literalinclude:: ../../../../../../../inputFiles/lagrangianContactMechanics/ContactMechanics_SingleFracCompression_benchmark.xml +.. literalinclude:: ../../../../../../../inputFiles/lagrangianContactMechanics/SingleFracCompression_benchmark.xml :language: xml :start-after: :end-before: @@ -110,7 +114,7 @@ To setup a coupling between rock and fracture deformations, we define three diff - The solver ``SurfaceGenerator`` defines the fracture region and rock toughness. -.. literalinclude:: ../../../../../../../inputFiles/lagrangianContactMechanics/ContactMechanics_SingleFracCompression_base.xml +.. literalinclude:: ../../../../../../../inputFiles/lagrangianContactMechanics/ContactMechanics_SingleFracCompression_benchmark.xml :language: xml :start-after: :end-before: @@ -124,7 +128,7 @@ For this specific problem, we simulate the elastic deformation and fracture slip Fracture surface slippage is assumed to be governed by the Coulomb failure criterion. The contact constitutive behavior is named ``fractureMaterial`` in the ``Coulomb`` block, where cohesion ``cohesion="0.0"`` and friction coefficient ``frictionCoefficient="0.577350269"`` are specified. -.. literalinclude:: ../../../../../../../inputFiles/lagrangianContactMechanics/ContactMechanics_SingleFracCompression_base.xml +.. literalinclude:: ../../../../../../../inputFiles/lagrangianContactMechanics/SingleFracCompression_base.xml :language: xml :start-after: :end-before: @@ -145,7 +149,7 @@ In the ``Tasks`` section, ``PackCollection`` tasks are defined to collect time h Either the entire field or specified named sets of indices in the field can be collected. In this example, ``tractionCollection`` and ``displacementJumpCollection`` tasks are specified to output the local traction ``fieldName="traction"`` and relative displacement ``fieldName="displacementJump"`` on the fracture surface. -.. literalinclude:: ../../../../../../../inputFiles/lagrangianContactMechanics/ContactMechanics_SingleFracCompression_base.xml +.. literalinclude:: ../../../../../../../inputFiles/lagrangianContactMechanics/SingleFracCompression_base.xml :language: xml :start-after: :end-before: @@ -170,7 +174,7 @@ The remaining parts of the outer boundaries are subjected to roller constraints. These boundary conditions are set up through the ``FieldSpecifications`` section. -.. literalinclude:: ../../../../../../../inputFiles/lagrangianContactMechanics/ContactMechanics_SingleFracCompression_base.xml +.. literalinclude:: ../../../../../../../inputFiles/lagrangianContactMechanics/SingleFracCompression_base.xml :language: xml :start-after: :end-before: diff --git a/src/docs/sphinx/advancedExamples/validationStudies/faultMechanics/singleFracCompression/singleFracCompressionFigure.py b/src/docs/sphinx/advancedExamples/validationStudies/faultMechanics/singleFracCompression/singleFracCompressionFigure.py index 01513ae946d..32afe5d42b2 100644 --- a/src/docs/sphinx/advancedExamples/validationStudies/faultMechanics/singleFracCompression/singleFracCompressionFigure.py +++ b/src/docs/sphinx/advancedExamples/validationStudies/faultMechanics/singleFracCompression/singleFracCompressionFigure.py @@ -96,8 +96,8 @@ def main(): geosDir = args.geosDir hdf5File1Path = outputDir + "/traction_history.hdf5" hdf5File2Path = outputDir + "/displacementJump_history.hdf5" - xmlFile1Path = geosDir + "/inputFiles/lagrangianContactMechanics/ContactMechanics_SingleFracCompression_base.xml" - xmlFile2Path = geosDir + "/inputFiles/lagrangianContactMechanics/ContactMechanics_SingleFracCompression_benchmark.xml" + xmlFile1Path = geosDir + "/inputFiles/lagrangianContactMechanics/SingleFracCompression_base.xml" + xmlFile2Path = geosDir + "/inputFiles/lagrangianContactMechanics/SingleFracCompression_benchmark.xml" # Read HDF5 # Global Coordinate of Fracture Element Center diff --git a/src/docs/sphinx/advancedExamples/validationStudies/faultMechanics/sneddon/Example.rst b/src/docs/sphinx/advancedExamples/validationStudies/faultMechanics/sneddon/Example.rst index 65a1ad48d5e..c8a81833630 100644 --- a/src/docs/sphinx/advancedExamples/validationStudies/faultMechanics/sneddon/Example.rst +++ b/src/docs/sphinx/advancedExamples/validationStudies/faultMechanics/sneddon/Example.rst @@ -29,8 +29,9 @@ The xml input files for the case with LagrangianContact solver are located at: .. code-block:: console - inputFiles/lagrangianContactMechanics/Sneddon_contactMechanics_base.xml - inputFiles/lagrangianContactMechanics/Sneddon_contactMechanics_benchmark.xml + inputFiles/lagrangianContactMechanics/Sneddon_base.xml + inputFiles/lagrangianContactMechanics/Sneddon_benchmark.xml + inputFiles/lagrangianContactMechanics/ContactMechanics_Sneddon_benchmark.xml The xml input files for the case with HydroFracture solver are located at: @@ -99,7 +100,7 @@ To setup a coupling between rock and fracture deformations in LagrangianContact - The solver ``SurfaceGenerator`` defines the fracture region and rock toughness. -.. literalinclude:: ../../../../../../../inputFiles/lagrangianContactMechanics/Sneddon_contactMechanics_benchmark.xml +.. literalinclude:: ../../../../../../../inputFiles/lagrangianContactMechanics/ContactMechanics_Sneddon_benchmark.xml :language: xml :start-after: :end-before: @@ -155,7 +156,7 @@ along the Z axes, 121 elements along the X axis and 921 elements along the Y axi The mesh for the case with LagrangianContact solver was also created using the internal mesh generator, as parametrized in the ``InternalMesh`` XML tag. The mesh discretizes the same compational domain (:math:`40\, m \, \times 40 \, m \, \times 1 \, m`) with 300 x 300 x 2 eight-node brick elements in the x, y, and z directions respectively. -.. literalinclude:: ../../../../../../../inputFiles/lagrangianContactMechanics/Sneddon_contactMechanics_benchmark.xml +.. literalinclude:: ../../../../../../../inputFiles/lagrangianContactMechanics/Sneddon_benchmark.xml :language: xml :start-after: :end-before: @@ -209,7 +210,7 @@ The static fracture is defined by a nodeset occupying a small region within the - The test case with LagrangianContact solver: -.. literalinclude:: ../../../../../../../inputFiles/lagrangianContactMechanics/Sneddon_contactMechanics_base.xml +.. literalinclude:: ../../../../../../../inputFiles/lagrangianContactMechanics/Sneddon_base.xml :language: xml :start-after: :end-before: @@ -242,7 +243,7 @@ In this example, a task is specified to output fracture aperture (normal opening - The test case with LagrangianContact solver: -.. literalinclude:: ../../../../../../../inputFiles/lagrangianContactMechanics/Sneddon_contactMechanics_base.xml +.. literalinclude:: ../../../../../../../inputFiles/lagrangianContactMechanics/Sneddon_base.xml :language: xml :start-after: :end-before: @@ -271,7 +272,7 @@ To run these three cases, use the following commands: ``path/to/geos -i inputFiles/efemFractureMechanics/Sneddon_embeddedFrac_verification.xml`` -``path/to/geos -i inputFiles/lagrangianContactMechanics/Sneddon_contactMechanics_benchmark.xml`` +``path/to/geos -i inputFiles/lagrangianContactMechanics/ContactMechanics_Sneddon_benchmark.xml`` ``path/to/geos -i inputFiles/hydraulicFracturing/Sneddon_hydroFrac_benchmark.xml``