diff --git a/inputFiles/chemoMechanics/ChemoMechanicsEigenStrain_CarbonateSystem_1DNoFlow.xml b/inputFiles/chemoMechanics/ChemoMechanicsEigenStrain_CarbonateSystem_1DNoFlow.xml new file mode 100644 index 00000000000..54e3269d2af --- /dev/null +++ b/inputFiles/chemoMechanics/ChemoMechanicsEigenStrain_CarbonateSystem_1DNoFlow.xml @@ -0,0 +1,64 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/chemoMechanics/ChemoMechanicsEigenStrain_UltramaficSystem_base.xml b/inputFiles/chemoMechanics/ChemoMechanicsEigenStrain_UltramaficSystem_base.xml new file mode 100644 index 00000000000..a4e356abc96 --- /dev/null +++ b/inputFiles/chemoMechanics/ChemoMechanicsEigenStrain_UltramaficSystem_base.xml @@ -0,0 +1,119 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/chemoMechanics/ChemoMechanicsEigenStrain_UltramaficSystem_wellbore.xml b/inputFiles/chemoMechanics/ChemoMechanicsEigenStrain_UltramaficSystem_wellbore.xml new file mode 100644 index 00000000000..fccf56f5e5b --- /dev/null +++ b/inputFiles/chemoMechanics/ChemoMechanicsEigenStrain_UltramaficSystem_wellbore.xml @@ -0,0 +1,127 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/chemoMechanics/ChemoMechanicsPoroElastic_UltramaficSystem_1DInjection.xml b/inputFiles/chemoMechanics/ChemoMechanicsPoroElastic_UltramaficSystem_1DInjection.xml new file mode 100644 index 00000000000..4ecb603c3cd --- /dev/null +++ b/inputFiles/chemoMechanics/ChemoMechanicsPoroElastic_UltramaficSystem_1DInjection.xml @@ -0,0 +1,132 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/chemoMechanics/ChemoMechanicsPoroElastic_UltramaficSystem_wellbore.xml b/inputFiles/chemoMechanics/ChemoMechanicsPoroElastic_UltramaficSystem_wellbore.xml new file mode 100644 index 00000000000..9fc7d6c7392 --- /dev/null +++ b/inputFiles/chemoMechanics/ChemoMechanicsPoroElastic_UltramaficSystem_wellbore.xml @@ -0,0 +1,126 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/chemoMechanics/ChemoMechanicsPoroelastic_CarbonateSystem_1DNoFlow.xml b/inputFiles/chemoMechanics/ChemoMechanicsPoroelastic_CarbonateSystem_1DNoFlow.xml new file mode 100644 index 00000000000..878b4a4c9a8 --- /dev/null +++ b/inputFiles/chemoMechanics/ChemoMechanicsPoroelastic_CarbonateSystem_1DNoFlow.xml @@ -0,0 +1,75 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/chemoMechanics/ChemoMechanicsPoroelastic_UltramaficSystem_base.xml b/inputFiles/chemoMechanics/ChemoMechanicsPoroelastic_UltramaficSystem_base.xml new file mode 100644 index 00000000000..57adfd91e71 --- /dev/null +++ b/inputFiles/chemoMechanics/ChemoMechanicsPoroelastic_UltramaficSystem_base.xml @@ -0,0 +1,120 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/chemoMechanics/ChemoMechanics_CarbonateSystem_base.xml b/inputFiles/chemoMechanics/ChemoMechanics_CarbonateSystem_base.xml new file mode 100644 index 00000000000..816582d4fab --- /dev/null +++ b/inputFiles/chemoMechanics/ChemoMechanics_CarbonateSystem_base.xml @@ -0,0 +1,167 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/chemoMechanics/ChemoMechanics_CarbonateSystem_initialAggregate.xml b/inputFiles/chemoMechanics/ChemoMechanics_CarbonateSystem_initialAggregate.xml new file mode 100644 index 00000000000..1d5feb21230 --- /dev/null +++ b/inputFiles/chemoMechanics/ChemoMechanics_CarbonateSystem_initialAggregate.xml @@ -0,0 +1,133 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/chemoMechanics/ChemoMechanics_CarbonateSystem_initialPrimaryConc.xml b/inputFiles/chemoMechanics/ChemoMechanics_CarbonateSystem_initialPrimaryConc.xml new file mode 100644 index 00000000000..8ce64e0159d --- /dev/null +++ b/inputFiles/chemoMechanics/ChemoMechanics_CarbonateSystem_initialPrimaryConc.xml @@ -0,0 +1,133 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/chemoMechanics/ChemoMechanics_UltramaficSystem_initialAggregate.xml b/inputFiles/chemoMechanics/ChemoMechanics_UltramaficSystem_initialAggregate.xml new file mode 100644 index 00000000000..06056cba380 --- /dev/null +++ b/inputFiles/chemoMechanics/ChemoMechanics_UltramaficSystem_initialAggregate.xml @@ -0,0 +1,79 @@ + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/chemoMechanics/ChemoMechanics_UltramaficSystem_initialPrimaryConc.xml b/inputFiles/chemoMechanics/ChemoMechanics_UltramaficSystem_initialPrimaryConc.xml new file mode 100644 index 00000000000..906b0b38dbe --- /dev/null +++ b/inputFiles/chemoMechanics/ChemoMechanics_UltramaficSystem_initialPrimaryConc.xml @@ -0,0 +1,79 @@ + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/chemoMechanics/ChemoMechanics_UltramaficSystem_logConc.xml b/inputFiles/chemoMechanics/ChemoMechanics_UltramaficSystem_logConc.xml new file mode 100644 index 00000000000..bacb1ca3ac9 --- /dev/null +++ b/inputFiles/chemoMechanics/ChemoMechanics_UltramaficSystem_logConc.xml @@ -0,0 +1,38 @@ + + + + + + + + + + + + + + diff --git a/inputFiles/chemoMechanics/ChemoMechanics_UltramaficSystem_surfaceArea.xml b/inputFiles/chemoMechanics/ChemoMechanics_UltramaficSystem_surfaceArea.xml new file mode 100644 index 00000000000..f21fafeb092 --- /dev/null +++ b/inputFiles/chemoMechanics/ChemoMechanics_UltramaficSystem_surfaceArea.xml @@ -0,0 +1,98 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/chemoMechanics/ChemoMechanics_UltramaficSystem_varyingSurfaceArea.xml b/inputFiles/chemoMechanics/ChemoMechanics_UltramaficSystem_varyingSurfaceArea.xml new file mode 100644 index 00000000000..9b43163b118 --- /dev/null +++ b/inputFiles/chemoMechanics/ChemoMechanics_UltramaficSystem_varyingSurfaceArea.xml @@ -0,0 +1,110 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/chemoMechanics/surfaceArea_tables/phi.geos b/inputFiles/chemoMechanics/surfaceArea_tables/phi.geos new file mode 100644 index 00000000000..e29b63b8639 --- /dev/null +++ b/inputFiles/chemoMechanics/surfaceArea_tables/phi.geos @@ -0,0 +1,201 @@ +1.0000000001 +0.9999000050998333 +0.9996000800893344 +0.9991004049785274 +0.9984012794176064 +0.9975031224974601 +0.9964064723309933 +0.9951119855158298 +0.993620436479149 +0.9919327167055711 +0.9900498338491681 +0.9879729107308383 +0.985703184222443 +0.9832420040192554 +0.9805908313024284 +0.9777512372933364 +0.9747249017017939 +0.9715136110702958 +0.9681192570165628 +0.9645438343768155 +0.9607894392523232 +0.9568582669619097 +0.952752609903211 +0.9484748553256046 +0.9440274830178357 +0.9394130629134758 +0.9346342526174487 +0.9296937948569544 +0.9245945148602107 +0.9193393176665182 +0.9139311853712282 +0.908373174309268 +0.9026684121809421 +0.8968200951237868 +0.8908314847343088 +0.8847059050434836 +0.8784467394499313 +0.8720574276147192 +0.865541462321766 +0.8589023863078477 +0.8521437890662114 +0.8452693036278185 +0.8382826033242335 +0.8311873985361709 +0.8239874334317032 +0.8166864826981108 +0.8092883482713321 +0.8017968560669413 +0.7942158527165467 +0.7865492023134549 +0.7788007831714049 +0.7709744846001153 +0.763074203701336 +0.7551038421890235 +0.7470673032371956 +0.7389684883589442 +0.7308112943200039 +0.7225996100901936 +0.7143373138359576 +0.7060282699571397 +0.697676326171031 +0.6892853106466263 +0.6808590291919248 +0.6724012624970027 +0.6639157634354735 +0.6554062544268405 +0.6468764248621298 +0.638329928595075 +0.6297703815010031 +0.6212013591054508 +0.612626394284416 +0.6040489750380253 +0.5954725423392698 +0.586900488059338 +0.5783361529709441 +0.569782824830923 +0.5612437365432349 +0.5527220644033939 +0.5442209264252074 +0.5357433807505847 +0.5272924241430486 +0.5188709905654523 +0.5104819498422892 +0.5021281064068468 +0.4938121981333462 +0.4855368952540795 +0.4773047993614458 +0.46911844249466417 +0.4609802863108336 +0.45289272133989467 +0.4448580663229411 +0.4368785676332222 +0.428956398779073 +0.4210936599879118 +0.4132923778703441 +0.40555450516332053 +0.3978819205512047 +0.3902764285635212 +0.3827397595480691 +0.37527356971800735 +0.36787944127144234 +0.3605588825819756 +0.3533133284586014 +0.3461441404732787 +0.3390526073544422 +0.33203994544466064 +0.3251072992205958 +0.31825574187337075 +0.3114862759474071 +0.3047998340357532 +0.29819727952988734 +0.2916794074219465 +0.28524694515730514 +0.2789005535353997 +0.27264082765667685 +0.26646829791352405 +0.260383431023029 +0.25438663109940385 +0.24847824076390423 +0.24265854229007083 +0.23692775878212177 +0.23128605538432867 +0.2257335405192169 +0.22027026715244175 +0.21489623408220537 +0.20961138725109782 +0.20441562107826397 +0.19930877980982267 +0.19429065888548816 +0.1893610063193741 +0.18451952409298925 +0.1797658695584678 +0.17509965685011336 +0.17052045830237153 +0.16602780587238664 +0.16162119256533922 +0.1573000738608037 +0.1530638691384112 +0.1489119631011477 +0.14484370719466705 +0.14085842102104495 +0.13695539374545324 +0.13313388549428185 +0.12939312874329137 +0.1257323296944279 +0.12215066963999 +0.11864730631288813 +0.11522137522179346 +0.11187199097002817 +0.10859824855710293 +0.10539922466186433 +0.1022739789062699 +0.09922155509886324 +0.09624098245707788 +0.09333127680755256 +0.09049144176369588 +0.08772046987979239 +0.08501734378099508 +0.0823810372686037 +0.07981051640007963 +0.07730474054329971 +0.07486266340460347 +0.07248323403023646 +0.07016539778084316 +0.06790809727870939 +0.06571027332750282 +0.06357086580430447 +0.06148881452377007 +0.059463060074303194 +0.05749254462616458 +0.05557621271148309 +0.05371301197617355 +0.051901893903805396 +0.05014181451150294 +0.04843173501799419 +0.04677062248395898 +0.04515745042486091 +0.0435911993964786 +0.0420708575533823 +0.04059542118063025 +0.03916389519898707 +0.03777529364399197 +0.03642864011922902 +0.03512296822417519 +0.03385732195702312 +0.03263075609289602 +0.03144233653789129 +0.03029114065940633 +0.029176257593216945 +0.028096788527793078 +0.027051846966350403 +0.026040558967148728 +0.025062063362558763 +0.024115511957429073 +0.023200069707293474 +0.022314914876966414 +0.021459239180080407 +0.02063224790012471 +0.019833159993548503 +0.019061208175494868 +0.01831563898873418 diff --git a/inputFiles/chemoMechanics/surfaceArea_tables/x.geos b/inputFiles/chemoMechanics/surfaceArea_tables/x.geos new file mode 100644 index 00000000000..383f4e6e6bd --- /dev/null +++ b/inputFiles/chemoMechanics/surfaceArea_tables/x.geos @@ -0,0 +1,201 @@ +0.0 +0.005 +0.01 +0.015 +0.02 +0.025 +0.03 +0.035 +0.04 +0.045 +0.05 +0.055 +0.06 +0.065 +0.07 +0.075 +0.08 +0.085 +0.09 +0.095 +0.1 +0.105 +0.11 +0.115 +0.12 +0.125 +0.13 +0.135 +0.14 +0.145 +0.15 +0.155 +0.16 +0.165 +0.17 +0.17500000000000002 +0.18 +0.185 +0.19 +0.195 +0.2 +0.20500000000000002 +0.21 +0.215 +0.22 +0.225 +0.23 +0.23500000000000001 +0.24 +0.245 +0.25 +0.255 +0.26 +0.265 +0.27 +0.275 +0.28 +0.28500000000000003 +0.29 +0.295 +0.3 +0.305 +0.31 +0.315 +0.32 +0.325 +0.33 +0.335 +0.34 +0.34500000000000003 +0.35000000000000003 +0.355 +0.36 +0.365 +0.37 +0.375 +0.38 +0.385 +0.39 +0.395 +0.4 +0.405 +0.41000000000000003 +0.41500000000000004 +0.42 +0.425 +0.43 +0.435 +0.44 +0.445 +0.45 +0.455 +0.46 +0.465 +0.47000000000000003 +0.47500000000000003 +0.48 +0.485 +0.49 +0.495 +0.5 +0.505 +0.51 +0.515 +0.52 +0.525 +0.53 +0.535 +0.54 +0.545 +0.55 +0.555 +0.56 +0.5650000000000001 +0.5700000000000001 +0.5750000000000001 +0.58 +0.585 +0.59 +0.595 +0.6 +0.605 +0.61 +0.615 +0.62 +0.625 +0.63 +0.635 +0.64 +0.645 +0.65 +0.655 +0.66 +0.665 +0.67 +0.675 +0.68 +0.685 +0.6900000000000001 +0.6950000000000001 +0.7000000000000001 +0.705 +0.71 +0.715 +0.72 +0.725 +0.73 +0.735 +0.74 +0.745 +0.75 +0.755 +0.76 +0.765 +0.77 +0.775 +0.78 +0.785 +0.79 +0.795 +0.8 +0.805 +0.81 +0.8150000000000001 +0.8200000000000001 +0.8250000000000001 +0.8300000000000001 +0.835 +0.84 +0.845 +0.85 +0.855 +0.86 +0.865 +0.87 +0.875 +0.88 +0.885 +0.89 +0.895 +0.9 +0.905 +0.91 +0.915 +0.92 +0.925 +0.93 +0.935 +0.9400000000000001 +0.9450000000000001 +0.9500000000000001 +0.9550000000000001 +0.96 +0.965 +0.97 +0.975 +0.98 +0.985 +0.99 +0.995 +1.0 diff --git a/inputFiles/chemoMechanics/surfaceArea_tables/y.geos b/inputFiles/chemoMechanics/surfaceArea_tables/y.geos new file mode 100644 index 00000000000..171538eb0b0 --- /dev/null +++ b/inputFiles/chemoMechanics/surfaceArea_tables/y.geos @@ -0,0 +1 @@ +0.0 \ No newline at end of file diff --git a/inputFiles/chemoMechanics/surfaceArea_tables/z.geos b/inputFiles/chemoMechanics/surfaceArea_tables/z.geos new file mode 100644 index 00000000000..171538eb0b0 --- /dev/null +++ b/inputFiles/chemoMechanics/surfaceArea_tables/z.geos @@ -0,0 +1 @@ +0.0 \ No newline at end of file diff --git a/inputFiles/singlePhaseFlow/reactiveTransport/2DUltramaficSystem/mixedReactionUltramaficSystem_base.xml b/inputFiles/singlePhaseFlow/reactiveTransport/2DUltramaficSystem/mixedReactionUltramaficSystem_base.xml new file mode 100644 index 00000000000..41f69bbd504 --- /dev/null +++ b/inputFiles/singlePhaseFlow/reactiveTransport/2DUltramaficSystem/mixedReactionUltramaficSystem_base.xml @@ -0,0 +1,89 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/singlePhaseFlow/reactiveTransport/2DUltramaficSystem/mixedReactionUltramaficSystem_initialAggregate.xml b/inputFiles/singlePhaseFlow/reactiveTransport/2DUltramaficSystem/mixedReactionUltramaficSystem_initialAggregate.xml new file mode 100644 index 00000000000..06056cba380 --- /dev/null +++ b/inputFiles/singlePhaseFlow/reactiveTransport/2DUltramaficSystem/mixedReactionUltramaficSystem_initialAggregate.xml @@ -0,0 +1,79 @@ + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/singlePhaseFlow/reactiveTransport/2DUltramaficSystem/mixedReactionUltramaficSystem_initialPrimaryConc.xml b/inputFiles/singlePhaseFlow/reactiveTransport/2DUltramaficSystem/mixedReactionUltramaficSystem_initialPrimaryConc.xml new file mode 100644 index 00000000000..906b0b38dbe --- /dev/null +++ b/inputFiles/singlePhaseFlow/reactiveTransport/2DUltramaficSystem/mixedReactionUltramaficSystem_initialPrimaryConc.xml @@ -0,0 +1,79 @@ + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/singlePhaseFlow/reactiveTransport/2DUltramaficSystem/mixedReactionUltramaficSystem_logConc.xml b/inputFiles/singlePhaseFlow/reactiveTransport/2DUltramaficSystem/mixedReactionUltramaficSystem_logConc.xml new file mode 100644 index 00000000000..8aefe228ee1 --- /dev/null +++ b/inputFiles/singlePhaseFlow/reactiveTransport/2DUltramaficSystem/mixedReactionUltramaficSystem_logConc.xml @@ -0,0 +1,38 @@ + + + + + + + + + + + + + + diff --git a/inputFiles/singlePhaseFlow/reactiveTransport/2DUltramaficSystem/mixedReactionUltramaficSystem_surfaceArea.xml b/inputFiles/singlePhaseFlow/reactiveTransport/2DUltramaficSystem/mixedReactionUltramaficSystem_surfaceArea.xml new file mode 100644 index 00000000000..589eec19a67 --- /dev/null +++ b/inputFiles/singlePhaseFlow/reactiveTransport/2DUltramaficSystem/mixedReactionUltramaficSystem_surfaceArea.xml @@ -0,0 +1,98 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/singlePhaseFlow/reactiveTransport/2DUltramaficSystem/mixedReactionUltramaficSystem_wellbore.xml b/inputFiles/singlePhaseFlow/reactiveTransport/2DUltramaficSystem/mixedReactionUltramaficSystem_wellbore.xml new file mode 100644 index 00000000000..77f4961de58 --- /dev/null +++ b/inputFiles/singlePhaseFlow/reactiveTransport/2DUltramaficSystem/mixedReactionUltramaficSystem_wellbore.xml @@ -0,0 +1,94 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/singlePhaseWell/thermal_singlePhase_well.xml b/inputFiles/singlePhaseWell/thermal_singlePhase_well.xml new file mode 100755 index 00000000000..78d27f5283c --- /dev/null +++ b/inputFiles/singlePhaseWell/thermal_singlePhase_well.xml @@ -0,0 +1,266 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/coreComponents/constitutive/CMakeLists.txt b/src/coreComponents/constitutive/CMakeLists.txt index 7606cfe25f6..ee3078aa266 100644 --- a/src/coreComponents/constitutive/CMakeLists.txt +++ b/src/coreComponents/constitutive/CMakeLists.txt @@ -155,6 +155,7 @@ set( constitutive_headers permeability/PermeabilityFields.hpp permeability/PressurePermeability.hpp permeability/ProppantPermeability.hpp + permeability/ReactivePressurePermeability.hpp permeability/SlipDependentPermeability.hpp permeability/WillisRichardsPermeability.hpp relativePermeability/BrooksCoreyBakerRelativePermeability.hpp @@ -185,6 +186,7 @@ set( constitutive_headers solid/DruckerPragerExtended.hpp solid/ModifiedCamClay.hpp solid/DelftEgg.hpp + solid/EigenstrainReactiveSolid.hpp solid/ElasticIsotropic.hpp solid/ElasticIsotropicPressureDependent.hpp solid/ElasticTransverseIsotropic.hpp @@ -192,6 +194,7 @@ set( constitutive_headers solid/InvariantDecompositions.hpp solid/PorousDamageSolid.hpp solid/PerfectlyPlastic.hpp + solid/PorousReactiveSolid.hpp solid/PorousSolid.hpp solid/PropertyConversions.hpp solid/ReactiveSolid.hpp @@ -207,10 +210,11 @@ set( constitutive_headers solid/SolidFields.hpp solid/porosity/PorosityFields.hpp solid/porosity/BiotPorosity.hpp + solid/porosity/BiotReactivePorosity.hpp solid/porosity/PorosityBase.hpp solid/porosity/PressurePorosity.hpp solid/porosity/ProppantPorosity.hpp - solid/porosity/ReactivePorosity.hpp + solid/porosity/ReactivePorosityBase.hpp thermalConductivity/MultiPhaseConstantThermalConductivity.hpp thermalConductivity/MultiPhaseThermalConductivityBase.hpp thermalConductivity/MultiPhaseThermalConductivityFields.hpp @@ -306,6 +310,7 @@ set( constitutive_sources permeability/PermeabilityBase.cpp permeability/PressurePermeability.cpp permeability/ProppantPermeability.cpp + permeability/ReactivePressurePermeability.cpp permeability/SlipDependentPermeability.cpp permeability/WillisRichardsPermeability.cpp relativePermeability/BrooksCoreyBakerRelativePermeability.cpp @@ -332,22 +337,25 @@ set( constitutive_sources solid/DruckerPragerExtended.cpp solid/ModifiedCamClay.cpp solid/DelftEgg.cpp + solid/EigenstrainReactiveSolid.cpp solid/ElasticIsotropic.cpp solid/ElasticIsotropicPressureDependent.cpp solid/ElasticTransverseIsotropic.cpp solid/ElasticOrthotropic.cpp solid/PorousDamageSolid.cpp solid/PerfectlyPlastic.cpp + solid/PorousReactiveSolid.cpp solid/PorousSolid.cpp solid/ReactiveSolid.cpp solid/SolidBase.cpp solid/SolidInternalEnergy.cpp solid/CeramicDamage.cpp solid/porosity/BiotPorosity.cpp + solid/porosity/BiotReactivePorosity.cpp solid/porosity/PorosityBase.cpp solid/porosity/PressurePorosity.cpp solid/porosity/ProppantPorosity.cpp - solid/porosity/ReactivePorosity.cpp + solid/porosity/ReactivePorosityBase.cpp thermalConductivity/MultiPhaseConstantThermalConductivity.cpp thermalConductivity/MultiPhaseThermalConductivityBase.cpp thermalConductivity/MultiPhaseVolumeWeightedThermalConductivity.cpp diff --git a/src/coreComponents/constitutive/ConstitutivePassThru.hpp b/src/coreComponents/constitutive/ConstitutivePassThru.hpp index d78d7bb1bf8..0313969eafd 100644 --- a/src/coreComponents/constitutive/ConstitutivePassThru.hpp +++ b/src/coreComponents/constitutive/ConstitutivePassThru.hpp @@ -36,21 +36,24 @@ #include "solid/ElasticIsotropicPressureDependent.hpp" #include "solid/ElasticTransverseIsotropic.hpp" #include "solid/ElasticOrthotropic.hpp" +#include "solid/EigenstrainReactiveSolid.hpp" #include "solid/PorousSolid.hpp" #include "solid/PorousDamageSolid.hpp" +#include "solid/PorousReactiveSolid.hpp" #include "solid/CompressibleSolid.hpp" #include "solid/ProppantSolid.hpp" #include "solid/CeramicDamage.hpp" #include "solid/ReactiveSolid.hpp" #include "solid/porosity/PressurePorosity.hpp" #include "solid/porosity/ProppantPorosity.hpp" -#include "solid/porosity/ReactivePorosity.hpp" +#include "solid/porosity/ReactivePorosityBase.hpp" #include "permeability/ConstantPermeability.hpp" #include "permeability/CarmanKozenyPermeability.hpp" #include "permeability/ExponentialDecayPermeability.hpp" #include "permeability/ParallelPlatesPermeability.hpp" #include "permeability/PressurePermeability.hpp" #include "permeability/ProppantPermeability.hpp" +#include "permeability/ReactivePressurePermeability.hpp" #include "permeability/SlipDependentPermeability.hpp" #include "permeability/WillisRichardsPermeability.hpp" #include "contact/CoulombFriction.hpp" @@ -341,8 +344,9 @@ struct ConstitutivePassThru< PorousSolidBase > PorousSolid< ElasticIsotropic, CarmanKozenyPermeability >, PorousSolid< ElasticTransverseIsotropic, CarmanKozenyPermeability >, PorousSolid< ElasticIsotropicPressureDependent, CarmanKozenyPermeability >, - PorousSolid< ElasticOrthotropic, CarmanKozenyPermeability > >::execute( constitutiveRelation, - std::forward< LAMBDA >( lambda ) ); + PorousSolid< ElasticOrthotropic, CarmanKozenyPermeability >, + PorousSolid< ElasticIsotropic, PressurePermeability > >::execute( constitutiveRelation, + std::forward< LAMBDA >( lambda ) ); } }; @@ -362,6 +366,38 @@ struct ConstitutivePassThru< PorousDamageSolidBase > } }; +/** + * Specialization for the EigenstrainReactiveSolid models. + */ +template<> +struct ConstitutivePassThru< EigenstrainReactiveSolidBase > +{ + template< typename LAMBDA > + static void execute( ConstitutiveBase & constitutiveRelation, LAMBDA && lambda ) + { + ConstitutivePassThruHandler< EigenstrainReactiveSolid< ElasticIsotropic, ConstantPermeability >, + EigenstrainReactiveSolid< ElasticIsotropic, CarmanKozenyPermeability >, + EigenstrainReactiveSolid< ElasticIsotropic, PressurePermeability >, + EigenstrainReactiveSolid< ElasticIsotropic, ReactivePressurePermeability > >::execute( constitutiveRelation, + std::forward< LAMBDA >( lambda ) ); + } +}; + +/** + * Specialization for the PorousReactiveSolid models. + */ +template<> +struct ConstitutivePassThru< PorousReactiveSolidBase > +{ + template< typename LAMBDA > + static void execute( ConstitutiveBase & constitutiveRelation, LAMBDA && lambda ) + { + ConstitutivePassThruHandler< PorousReactiveSolid< ElasticIsotropic, ConstantPermeability >, + PorousReactiveSolid< ElasticIsotropic, CarmanKozenyPermeability > >::execute( constitutiveRelation, + std::forward< LAMBDA >( lambda ) ); + } +}; + /** * Specialization for the CompressibleSolid models. */ @@ -406,9 +442,9 @@ struct ConstitutivePassThru< ReactiveSolidBase > template< typename LAMBDA > static void execute( ConstitutiveBase & constitutiveRelation, LAMBDA && lambda ) { - ConstitutivePassThruHandler< ReactiveSolid< ReactivePorosity, ConstantPermeability >, - ReactiveSolid< ReactivePorosity, CarmanKozenyPermeability >, - ReactiveSolid< ReactivePorosity, PressurePermeability > + ConstitutivePassThruHandler< ReactiveSolid< ReactivePorosityBase, ConstantPermeability >, + ReactiveSolid< ReactivePorosityBase, CarmanKozenyPermeability >, + ReactiveSolid< ReactivePorosityBase, PressurePermeability > >::execute( constitutiveRelation, std::forward< LAMBDA >( lambda ) ); } @@ -416,9 +452,9 @@ struct ConstitutivePassThru< ReactiveSolidBase > template< typename LAMBDA > static void execute( ConstitutiveBase const & constitutiveRelation, LAMBDA && lambda ) { - ConstitutivePassThruHandler< ReactiveSolid< ReactivePorosity, ConstantPermeability >, - ReactiveSolid< ReactivePorosity, CarmanKozenyPermeability >, - ReactiveSolid< ReactivePorosity, PressurePermeability > + ConstitutivePassThruHandler< ReactiveSolid< ReactivePorosityBase, ConstantPermeability >, + ReactiveSolid< ReactivePorosityBase, CarmanKozenyPermeability >, + ReactiveSolid< ReactivePorosityBase, PressurePermeability > >::execute( constitutiveRelation, std::forward< LAMBDA >( lambda ) ); } @@ -465,6 +501,10 @@ struct ConstitutivePassThru< CoupledSolidBase > CompressibleSolid< PressurePorosity, PressurePermeability >, CompressibleSolid< PressurePorosity, SlipDependentPermeability >, CompressibleSolid< PressurePorosity, WillisRichardsPermeability >, + EigenstrainReactiveSolid< ElasticIsotropic, ConstantPermeability >, + EigenstrainReactiveSolid< ElasticIsotropic, CarmanKozenyPermeability >, + EigenstrainReactiveSolid< ElasticIsotropic, PressurePermeability >, + EigenstrainReactiveSolid< ElasticIsotropic, ReactivePressurePermeability >, PorousSolid< DruckerPragerExtended, ConstantPermeability >, PorousSolid< ModifiedCamClay, ConstantPermeability >, PorousSolid< DelftEgg, ConstantPermeability >, @@ -487,13 +527,16 @@ struct ConstitutivePassThru< CoupledSolidBase > PorousSolid< ElasticTransverseIsotropic, CarmanKozenyPermeability >, PorousSolid< ElasticIsotropicPressureDependent, CarmanKozenyPermeability >, PorousSolid< ElasticOrthotropic, CarmanKozenyPermeability >, + PorousSolid< ElasticIsotropic, PressurePermeability >, PorousDamageSolid< DamageSpectral< ElasticIsotropic > >, PorousDamageSolid< DamageVolDev< ElasticIsotropic > >, PorousDamageSolid< Damage< ElasticIsotropic > >, - ReactiveSolid< ReactivePorosity, ConstantPermeability >, - ReactiveSolid< ReactivePorosity, CarmanKozenyPermeability >, - ReactiveSolid< ReactivePorosity, PressurePermeability > >::execute( constitutiveRelation, - std::forward< LAMBDA >( lambda ) ); + PorousReactiveSolid< ElasticIsotropic, ConstantPermeability >, + PorousReactiveSolid< ElasticIsotropic, CarmanKozenyPermeability >, + ReactiveSolid< ReactivePorosityBase, ConstantPermeability >, + ReactiveSolid< ReactivePorosityBase, CarmanKozenyPermeability >, + ReactiveSolid< ReactivePorosityBase, PressurePermeability > >::execute( constitutiveRelation, + std::forward< LAMBDA >( lambda ) ); } template< typename LAMBDA > @@ -506,6 +549,10 @@ struct ConstitutivePassThru< CoupledSolidBase > CompressibleSolid< PressurePorosity, PressurePermeability >, CompressibleSolid< PressurePorosity, SlipDependentPermeability >, CompressibleSolid< PressurePorosity, WillisRichardsPermeability >, + EigenstrainReactiveSolid< ElasticIsotropic, ConstantPermeability >, + EigenstrainReactiveSolid< ElasticIsotropic, CarmanKozenyPermeability >, + EigenstrainReactiveSolid< ElasticIsotropic, PressurePermeability >, + EigenstrainReactiveSolid< ElasticIsotropic, ReactivePressurePermeability >, PorousSolid< DruckerPragerExtended, ConstantPermeability >, PorousSolid< ModifiedCamClay, ConstantPermeability >, PorousSolid< DelftEgg, ConstantPermeability >, @@ -528,13 +575,16 @@ struct ConstitutivePassThru< CoupledSolidBase > PorousSolid< ElasticTransverseIsotropic, CarmanKozenyPermeability >, PorousSolid< ElasticIsotropicPressureDependent, CarmanKozenyPermeability >, PorousSolid< ElasticOrthotropic, CarmanKozenyPermeability >, + PorousSolid< ElasticIsotropic, PressurePermeability >, PorousDamageSolid< DamageSpectral< ElasticIsotropic > >, PorousDamageSolid< DamageVolDev< ElasticIsotropic > >, PorousDamageSolid< Damage< ElasticIsotropic > >, - ReactiveSolid< ReactivePorosity, ConstantPermeability >, - ReactiveSolid< ReactivePorosity, CarmanKozenyPermeability >, - ReactiveSolid< ReactivePorosity, PressurePermeability > >::execute( constitutiveRelation, - std::forward< LAMBDA >( lambda ) ); + PorousReactiveSolid< ElasticIsotropic, ConstantPermeability >, + PorousReactiveSolid< ElasticIsotropic, CarmanKozenyPermeability >, + ReactiveSolid< ReactivePorosityBase, ConstantPermeability >, + ReactiveSolid< ReactivePorosityBase, CarmanKozenyPermeability >, + ReactiveSolid< ReactivePorosityBase, PressurePermeability > >::execute( constitutiveRelation, + std::forward< LAMBDA >( lambda ) ); } }; diff --git a/src/coreComponents/constitutive/HPCReact b/src/coreComponents/constitutive/HPCReact index 5268dc5c2a5..9ec1bf4f1d5 160000 --- a/src/coreComponents/constitutive/HPCReact +++ b/src/coreComponents/constitutive/HPCReact @@ -1 +1 @@ -Subproject commit 5268dc5c2a59e9dae23f435fba9aecf5a5c5de33 +Subproject commit 9ec1bf4f1d5ab35e7d1b25529e06720162fe1b1c diff --git a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.cpp b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.cpp index 69b7b750934..ecb4e3b3340 100644 --- a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.cpp +++ b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.cpp @@ -81,7 +81,7 @@ void ReactiveSinglePhaseFluid< BASE >::postInputInitialization() switch( m_chemicalSystemType ) { case ChemicalSystemType::ultramafic: - m_numPrimarySpecies = 9; + m_numPrimarySpecies = 4; m_numSecondarySpecies = 16; m_numKineticReactions = 5; break; @@ -98,6 +98,12 @@ void ReactiveSinglePhaseFluid< BASE >::postInputInitialization() m_numKineticReactions = 0; break; + case ChemicalSystemType::forge: + m_numPrimarySpecies = 10; + m_numSecondarySpecies = 19; + m_numKineticReactions = 5; + break; + case ChemicalSystemType::chainSerialAllKinetic: m_numPrimarySpecies = 3; m_numSecondarySpecies = 0; diff --git a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.hpp b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.hpp index 6346d030c4f..5e581447672 100644 --- a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.hpp +++ b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.hpp @@ -52,6 +52,7 @@ enum class ChemicalSystemType : integer carbonate, carbonateAllEquilibrium, ultramafic, + forge, momasEasy, momasMedium, chainSerialAllKinetic @@ -228,6 +229,7 @@ class ReactiveSinglePhaseFluid : public BASE typename ReactiveSinglePhaseFluid< BASE >::template ReactionKernelWrapper< hpcReact::geochemistry::ultramaficSystemType >, typename ReactiveSinglePhaseFluid< BASE >::template ReactionKernelWrapper< hpcReact::geochemistry::carbonateSystemType >, typename ReactiveSinglePhaseFluid< BASE >::template ReactionKernelWrapper< hpcReact::geochemistry::carbonateSystemAllEquilibriumType >, + typename ReactiveSinglePhaseFluid< BASE >::template ReactionKernelWrapper< hpcReact::geochemistry::forgeSystemType >, typename ReactiveSinglePhaseFluid< BASE >::template ReactionKernelWrapper< hpcReact::ChainGeneric::serialAllKineticType >, typename ReactiveSinglePhaseFluid< BASE >::template ReactionKernelWrapper< hpcReact::MoMasBenchmark::mediumCaseType >, typename ReactiveSinglePhaseFluid< BASE >::template ReactionKernelWrapper< hpcReact::MoMasBenchmark::easyCaseType > > @@ -282,6 +284,20 @@ class ReactiveSinglePhaseFluid : public BASE m_numSecondarySpecies, m_numKineticReactions, carbonateSystemAllEquilibrium ); + case ChemicalSystemType::forge: + return ReactionKernelWrapper< forgeSystemType >( m_primarySpeciesAggregateConcentration, + m_primarySpeciesMobileAggregateConcentration, + m_dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations, + m_dPrimarySpeciesMobileAggregateConcentration_dLogPrimarySpeciesConcentrations, + m_initialPrimarySpeciesConcentration, + m_secondarySpeciesConcentration, + m_kineticReactionRates, + m_aggregateSpeciesRates, + m_dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations, + m_numPrimarySpecies, + m_numSecondarySpecies, + m_numKineticReactions, + forgeSystem ); case ChemicalSystemType::chainSerialAllKinetic: return ReactionKernelWrapper< serialAllKineticType >( m_primarySpeciesAggregateConcentration, m_primarySpeciesMobileAggregateConcentration, @@ -516,6 +532,7 @@ ENUM_STRINGS( ChemicalSystemType, "carbonate", "carbonateAllEquilibrium", "ultramafic", + "forge", "momasEasy", "momasMedium", "chainSerialAllKinetic" ); diff --git a/src/coreComponents/constitutive/fluid/singlefluid/CompressibleSinglePhaseFluid.hpp b/src/coreComponents/constitutive/fluid/singlefluid/CompressibleSinglePhaseFluid.hpp index e177053795c..47ec5a0da0c 100644 --- a/src/coreComponents/constitutive/fluid/singlefluid/CompressibleSinglePhaseFluid.hpp +++ b/src/coreComponents/constitutive/fluid/singlefluid/CompressibleSinglePhaseFluid.hpp @@ -169,7 +169,7 @@ class CompressibleSinglePhaseFluid : public SingleFluidBase localIndex const numPts ) override; /// Type of kernel wrapper for in-kernel update (TODO: support multiple EAT combinations, not just this combination) - using KernelWrapper = CompressibleSinglePhaseUpdate< ExponentApproximationType::Full, ExponentApproximationType::Linear >; + using KernelWrapper = CompressibleSinglePhaseUpdate< ExponentApproximationType::Linear, ExponentApproximationType::Linear >; /** * @brief Create an update kernel wrapper. diff --git a/src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.hpp b/src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.hpp index 83958834e33..c1e1c1bed49 100644 --- a/src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.hpp +++ b/src/coreComponents/constitutive/fluid/singlefluid/ThermalCompressibleSinglePhaseFluid.hpp @@ -227,7 +227,7 @@ class ThermalCompressibleSinglePhaseFluid : public CompressibleSinglePhaseFluid using CompressibleSinglePhaseFluid::m_densityModelType; /// Type of kernel wrapper for in-kernel update (TODO: support multiple EAT combinations, not just this combination) - using KernelWrapper = ThermalCompressibleSinglePhaseUpdate< ExponentApproximationType::Full, ExponentApproximationType::Linear, ExponentApproximationType::Linear >; + using KernelWrapper = ThermalCompressibleSinglePhaseUpdate< ExponentApproximationType::Linear, ExponentApproximationType::Linear, ExponentApproximationType::Linear >; /** * @brief Create an update kernel wrapper. diff --git a/src/coreComponents/constitutive/permeability/CarmanKozenyPermeability.hpp b/src/coreComponents/constitutive/permeability/CarmanKozenyPermeability.hpp index 2a121e0d36c..785c4c3878b 100644 --- a/src/coreComponents/constitutive/permeability/CarmanKozenyPermeability.hpp +++ b/src/coreComponents/constitutive/permeability/CarmanKozenyPermeability.hpp @@ -54,9 +54,10 @@ class CarmanKozenyPermeabilityUpdate : public PermeabilityBaseUpdate virtual void updateFromPressureAndPorosity( localIndex const k, localIndex const q, real64 const & pressure, + real64 const & pressure_n, real64 const & porosity ) const override { - GEOS_UNUSED_VAR( pressure ); + GEOS_UNUSED_VAR( pressure, pressure_n ); compute( porosity, m_permeability[k][q], diff --git a/src/coreComponents/constitutive/permeability/PermeabilityBase.cpp b/src/coreComponents/constitutive/permeability/PermeabilityBase.cpp index 2333490d7f7..15fa07a673e 100644 --- a/src/coreComponents/constitutive/permeability/PermeabilityBase.cpp +++ b/src/coreComponents/constitutive/permeability/PermeabilityBase.cpp @@ -33,6 +33,7 @@ PermeabilityBase::PermeabilityBase( string const & name, Group * const parent ): ConstitutiveBase( name, parent ) { registerField< fields::permeability::permeability >( &m_permeability ); + registerField< fields::permeability::permeability_n >( &m_permeability_n ); registerField< fields::permeability::dPerm_dPressure >( &m_dPerm_dPressure ); } @@ -55,10 +56,33 @@ void PermeabilityBase::allocateConstitutiveData( Group & parent, { // NOTE: enforcing 1 quadrature point m_permeability.resize( 0, 1, 3 ); + m_permeability_n.resize( 0, 1, 3 ); m_dPerm_dPressure.resize( 0, 1, 3 ); ConstitutiveBase::allocateConstitutiveData( parent, numPts ); } +void PermeabilityBase::saveConvergedState() const +{ + localIndex const numE = m_permeability.size( 0 ); + integer constexpr numQuad = 1; // NOTE: enforcing 1 quadrature point + + auto permView_n = m_permeability_n.toView(); + + forAll< parallelDevicePolicy<> >( numE, [=] GEOS_HOST_DEVICE ( localIndex const ei ) + { + for( localIndex q = 0; q < numQuad; ++q ) + { + real64 const permComponents[3] = { m_permeability[ei][q][0], + m_permeability[ei][q][1], + m_permeability[ei][q][2] }; + for( integer dim=0; dim < 3; ++dim ) + { + permView_n[ei][q][dim] = permComponents[dim]; + } + } + } ); +} + } } /* namespace geos */ diff --git a/src/coreComponents/constitutive/permeability/PermeabilityBase.hpp b/src/coreComponents/constitutive/permeability/PermeabilityBase.hpp index 190b8f1b4f6..0e00afc80f5 100644 --- a/src/coreComponents/constitutive/permeability/PermeabilityBase.hpp +++ b/src/coreComponents/constitutive/permeability/PermeabilityBase.hpp @@ -51,9 +51,10 @@ class PermeabilityBaseUpdate virtual void updateFromPressureAndPorosity( localIndex const k, localIndex const q, real64 const & pressure, + real64 const & pressure_n, real64 const & porosity ) const { - GEOS_UNUSED_VAR( k, q, pressure, porosity ); + GEOS_UNUSED_VAR( k, q, pressure, pressure_n, porosity ); } GEOS_HOST_DEVICE @@ -136,11 +137,17 @@ class PermeabilityBase : public ConstitutiveBase virtual void initializeState() const {} + /// Save state data in preparation for next timestep + virtual void saveConvergedState() const override; + protected: /// Vector of absolute permeability array3d< real64 > m_permeability; + /// Vector of absolute permeability at previous time step + array3d< real64 > m_permeability_n; + /// Vector of derivative of permeability wrt pressure array3d< real64 > m_dPerm_dPressure; }; diff --git a/src/coreComponents/constitutive/permeability/PermeabilityFields.hpp b/src/coreComponents/constitutive/permeability/PermeabilityFields.hpp index e1a48f17618..62511640d68 100644 --- a/src/coreComponents/constitutive/permeability/PermeabilityFields.hpp +++ b/src/coreComponents/constitutive/permeability/PermeabilityFields.hpp @@ -39,6 +39,14 @@ DECLARE_FIELD( permeability, WRITE_AND_READ, "Rock permeability" ); +DECLARE_FIELD( permeability_n, + "permeability_n", + array3d< real64 >, + -1, + LEVEL_0, + WRITE_AND_READ, + "Rock permeability at the previous converged time step" ); + DECLARE_FIELD( dPerm_dPressure, "dPerm_dPressure", array3d< real64 >, diff --git a/src/coreComponents/constitutive/permeability/PressurePermeability.cpp b/src/coreComponents/constitutive/permeability/PressurePermeability.cpp index 838439efd88..7df34504029 100644 --- a/src/coreComponents/constitutive/permeability/PressurePermeability.cpp +++ b/src/coreComponents/constitutive/permeability/PressurePermeability.cpp @@ -40,8 +40,12 @@ PressurePermeability::PressurePermeability( string const & name, Group * const p setInputFlag( InputFlags::REQUIRED ). setDescription( "Pressure dependence coefficients for each permeability component." ); - registerWrapper( viewKeyStruct::referencePressureString(), &m_referencePressure ). + registerWrapper( viewKeyStruct::defaultReferencePressureString(), &m_defaultReferencePressure ). setInputFlag( InputFlags::REQUIRED ). + setDescription( "Default reference pressure for the pressure permeability model" ); + + registerWrapper( viewKeyStruct::referencePressureString(), &m_referencePressure ). + setApplyDefaultValue( -1.0 ). setDescription( "Reference pressure for the pressure permeability model [Pa]" ); registerWrapper( viewKeyStruct::referencePermeabilityString(), &m_referencePermeability ). @@ -56,8 +60,13 @@ PressurePermeability::PressurePermeability( string const & name, Group * const p registerWrapper( viewKeyStruct::pressureModelTypeString(), &m_presModelType ). setInputFlag( InputFlags::OPTIONAL ). - setApplyDefaultValue( PressureModelType::Hyperbolic ). + setApplyDefaultValue( PressureModelType::Exponential ). setDescription( "Type of the pressure dependence model. " ); + + registerWrapper( viewKeyStruct::explicitUpdateFlagString(), &m_explicitFlag ). + setInputFlag( InputFlags::OPTIONAL ). + setApplyDefaultValue( 0 ). + setDescription( "The flag for explicit update. " ); } void PressurePermeability::postInputInitialization() @@ -70,6 +79,9 @@ void PressurePermeability::postInputInitialization() i ), getDataContext() ); } + + this->getWrapper< array1d< real64 > >( viewKeyStruct::referencePressureString() ). + setApplyDefaultValue( m_defaultReferencePressure ); } void PressurePermeability::allocateConstitutiveData( Group & parent, @@ -98,6 +110,7 @@ void PressurePermeability::initializeState() const integer constexpr numQuad = 1; // NOTE: enforcing 1 quadrature point auto permView = m_permeability.toView(); + auto permView_n = m_permeability_n.toView(); real64 const permComponents[3] = { m_referencePermeabilityComponents[0], m_referencePermeabilityComponents[1], m_referencePermeabilityComponents[2] }; @@ -112,6 +125,7 @@ void PressurePermeability::initializeState() const if( permView[ei][q][dim] < 0 ) { permView[ei][q][dim] = permComponents[dim]; + permView_n[ei][q][dim] = permComponents[dim]; } } } diff --git a/src/coreComponents/constitutive/permeability/PressurePermeability.hpp b/src/coreComponents/constitutive/permeability/PressurePermeability.hpp index 9d654048a0e..de39809d832 100644 --- a/src/coreComponents/constitutive/permeability/PressurePermeability.hpp +++ b/src/coreComponents/constitutive/permeability/PressurePermeability.hpp @@ -44,17 +44,21 @@ class PressurePermeabilityUpdate : public PermeabilityBaseUpdate PressurePermeabilityUpdate( PressureModelType const & presModelType, R1Tensor const pressureDependenceConstants, - real64 const & referencePressure, + arrayView1d< real64 const > const & referencePressure, real64 const & maxPermeability, - arrayView3d< real64 > const & referencePermeability, + arrayView3d< real64 const > const & referencePermeability, arrayView3d< real64 > const & permeability, - arrayView3d< real64 > const & dPerm_dPressure ) + arrayView3d< real64 const > const & permeability_n, + arrayView3d< real64 > const & dPerm_dPressure, + integer const & explicitFlag ) : PermeabilityBaseUpdate( permeability, dPerm_dPressure ), m_presModelType( presModelType ), m_pressureDependenceConstants( pressureDependenceConstants ), m_referencePressure( referencePressure ), m_maxPermeability( maxPermeability ), - m_referencePermeability( referencePermeability ) + m_referencePermeability( referencePermeability ), + m_permeability_n( permeability_n ), + m_explicitFlag( explicitFlag ) {} GEOS_HOST_DEVICE @@ -76,11 +80,12 @@ class PressurePermeabilityUpdate : public PermeabilityBaseUpdate virtual void updateFromPressureAndPorosity( localIndex const k, localIndex const q, real64 const & pressure, + real64 const & pressure_n, real64 const & porosity ) const override { GEOS_UNUSED_VAR( q, porosity ); - real64 const deltaPressure = pressure - m_referencePressure; + real64 deltaPressure = (m_explicitFlag)? (pressure_n - m_referencePressure[k]):(pressure - m_referencePressure[k]); real64 referencePermeability[3]; @@ -88,6 +93,12 @@ class PressurePermeabilityUpdate : public PermeabilityBaseUpdate referencePermeability[1] = m_referencePermeability[k][0][1]; referencePermeability[2] = m_referencePermeability[k][0][2]; + real64 permeability_n[3]; + + permeability_n[0] = m_permeability_n[k][0][0]; + permeability_n[1] = m_permeability_n[k][0][1]; + permeability_n[2] = m_permeability_n[k][0][2]; + switch( m_presModelType ) { case PressureModelType::Exponential: @@ -98,6 +109,32 @@ class PressurePermeabilityUpdate : public PermeabilityBaseUpdate m_permeability[k][0], m_dPerm_dPressure[k][0] ); + if( m_explicitFlag ) + { + for( localIndex i=0; i < m_permeability[k][0].size(); i++ ) + { + m_dPerm_dPressure[k][0][i] = 0.0; + + if( m_maxPermeability < 1.0 ) + { + real64 perm = m_permeability[k][0][i]; + + if( perm < referencePermeability[i] ) + { + perm = referencePermeability[i]; + } + + // Ensure permeability non-decreasing + if( perm < permeability_n[i] ) + { + perm = permeability_n[i]; + } + + m_permeability[k][0][i] = (perm > m_maxPermeability)? m_maxPermeability:perm; + } + } + } + break; } case PressureModelType::Hyperbolic: @@ -109,6 +146,20 @@ class PressurePermeabilityUpdate : public PermeabilityBaseUpdate m_permeability[k][0], m_dPerm_dPressure[k][0] ); + if( m_explicitFlag ) + { + for( localIndex i=0; i < m_permeability[k][0].size(); i++ ) + { + m_dPerm_dPressure[k][0][i] = 0.0; + + real64 const perm = m_permeability[k][0][i]; + + if( perm < referencePermeability[i] ) + { + m_permeability[k][0][i] = referencePermeability[i]; + } + } + } break; } default: @@ -127,12 +178,16 @@ class PressurePermeabilityUpdate : public PermeabilityBaseUpdate R1Tensor m_pressureDependenceConstants; /// Reference pressure in the model - real64 const m_referencePressure; + arrayView1d< real64 const > const m_referencePressure; /// Maximum permeability real64 const m_maxPermeability; - arrayView3d< real64 > m_referencePermeability; + arrayView3d< real64 const > const m_referencePermeability; + + arrayView3d< real64 const > const m_permeability_n; + + integer m_explicitFlag; }; @@ -165,17 +220,21 @@ class PressurePermeability : public PermeabilityBase m_maxPermeability, m_referencePermeability, m_permeability, - m_dPerm_dPressure ); + m_permeability_n, + m_dPerm_dPressure, + m_explicitFlag ); } struct viewKeyStruct : public PermeabilityBase::viewKeyStruct { static constexpr char const * referencePermeabilityComponentsString() { return "referencePermeabilityComponents"; } static constexpr char const * pressureDependenceConstantsString() { return "pressureDependenceConstants"; } + static constexpr char const * defaultReferencePressureString() { return "defaultReferencePressure"; } static constexpr char const * referencePressureString() { return "referencePressure"; } static constexpr char const * referencePermeabilityString() { return "referencePermeability"; } static constexpr char const * maxPermeabilityString() { return "maxPermeability"; } static constexpr char const * pressureModelTypeString() { return "pressureModelType"; } + static constexpr char const * explicitUpdateFlagString() { return "explicitUpdateFlag"; } }; virtual void initializeState() const override final; @@ -193,7 +252,8 @@ class PressurePermeability : public PermeabilityBase R1Tensor m_pressureDependenceConstants; /// Reference pressure in the model - real64 m_referencePressure; + real64 m_defaultReferencePressure; + array1d< real64 > m_referencePressure; /// Maximum permeability real64 m_maxPermeability; @@ -203,6 +263,9 @@ class PressurePermeability : public PermeabilityBase /// Pressure dependence model type PressureModelType m_presModelType; + /// Explicit update flag + integer m_explicitFlag; + }; GEOS_HOST_DEVICE diff --git a/src/coreComponents/constitutive/permeability/ReactivePressurePermeability.cpp b/src/coreComponents/constitutive/permeability/ReactivePressurePermeability.cpp new file mode 100644 index 00000000000..30d9a76ec55 --- /dev/null +++ b/src/coreComponents/constitutive/permeability/ReactivePressurePermeability.cpp @@ -0,0 +1,154 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file ReactivePressurePermeability.cpp + */ + +#include "ReactivePressurePermeability.hpp" + +namespace geos +{ + +using namespace dataRepository; + +namespace constitutive +{ + + +ReactivePressurePermeability::ReactivePressurePermeability( string const & name, Group * const parent ): + PermeabilityBase( name, parent ) +{ + registerWrapper( viewKeyStruct::referencePermeabilityComponentsString(), &m_referencePermeabilityComponents ). + setInputFlag( InputFlags::REQUIRED ). + setRestartFlags( RestartFlags::NO_WRITE ). + setDescription( "Reference xx, yy and zz components of a diagonal permeability tensor." ); + + registerWrapper( viewKeyStruct::pressureDependenceConstantsString(), &m_pressureDependenceConstants ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "Pressure dependence coefficients for each permeability component." ); + + registerWrapper( viewKeyStruct::defaultReferencePressureString(), &m_defaultReferencePressure ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "Default reference pressure for the pressure permeability model" ); + + registerWrapper( viewKeyStruct::referencePressureString(), &m_referencePressure ). + setApplyDefaultValue( -1.0 ). + setDescription( "Reference pressure for the pressure permeability model [Pa]" ); + + registerWrapper( viewKeyStruct::referencePermeabilityString(), &m_referencePermeability ). + setApplyDefaultValue( 0.0 ). + setPlotLevel( PlotLevel::LEVEL_0 ). + setDescription( "Reference permeability field" ); + + registerWrapper( viewKeyStruct::defaultReferencePorosityString(), &m_defaultReferencePorosity ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "Default reference porosity" ); + + registerWrapper( viewKeyStruct::referencePorosityString(), &m_referencePorosity ). + setApplyDefaultValue( -1.0 ). + setDescription( "Reference porosity [Pa]" ); + + registerWrapper( viewKeyStruct::maxPermeabilityString(), &m_maxPermeability ). + setApplyDefaultValue( 1.0 ). + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "Max. permeability can be reached." ); + + registerWrapper( viewKeyStruct::minPermeabilityString(), &m_minPermeability ). + setApplyDefaultValue( 1e-20 ). + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "Min. permeability can be reached." ); + + registerWrapper( viewKeyStruct::pressureModelTypeString(), &m_presModelType ). + setInputFlag( InputFlags::OPTIONAL ). + setApplyDefaultValue( PressureModelType::Exponential ). + setDescription( "Type of the pressure dependence model. " ); + + registerWrapper( viewKeyStruct::explicitUpdateFlagString(), &m_explicitFlag ). + setInputFlag( InputFlags::OPTIONAL ). + setApplyDefaultValue( 0 ). + setDescription( "The flag for explicit update. " ); +} + +void ReactivePressurePermeability::postInputInitialization() +{ + for( localIndex i=0; i < 3; i++ ) + { + GEOS_ERROR_IF( std::abs( m_pressureDependenceConstants[i] ) < 1e-15 && m_presModelType == PressureModelType::Hyperbolic, + GEOS_FMT( "The pressure dependent constant at component {} is too close to zero, " + "which is not allowed for the hyperbolic model.", + i ), + getDataContext() ); + } + + this->getWrapper< array1d< real64 > >( viewKeyStruct::referencePressureString() ). + setApplyDefaultValue( m_defaultReferencePressure ); + + this->getWrapper< array1d< real64 > >( viewKeyStruct::referencePorosityString() ). + setApplyDefaultValue( m_defaultReferencePorosity ); +} + +void ReactivePressurePermeability::allocateConstitutiveData( Group & parent, + localIndex const numPts ) +{ + m_referencePermeability.resize( 0, 1, 3 ); // 0 to resize and assign default value later + + PermeabilityBase::allocateConstitutiveData( parent, numPts ); + + integer constexpr numQuad = 1; // NOTE: enforcing 1 quadrature point + + for( localIndex ei = 0; ei < parent.size(); ++ei ) + { + for( localIndex q = 0; q < numQuad; ++q ) + { + m_referencePermeability[ei][q][0] = m_referencePermeabilityComponents[0]; + m_referencePermeability[ei][q][1] = m_referencePermeabilityComponents[1]; + m_referencePermeability[ei][q][2] = m_referencePermeabilityComponents[2]; + } + } +} + +void ReactivePressurePermeability::initializeState() const +{ + localIndex const numE = m_permeability.size( 0 ); + integer constexpr numQuad = 1; // NOTE: enforcing 1 quadrature point + + auto permView = m_permeability.toView(); + auto permView_n = m_permeability_n.toView(); + real64 const permComponents[3] = { m_referencePermeabilityComponents[0], + m_referencePermeabilityComponents[1], + m_referencePermeabilityComponents[2] }; + + forAll< parallelDevicePolicy<> >( numE, [=] GEOS_HOST_DEVICE ( localIndex const ei ) + { + for( localIndex q = 0; q < numQuad; ++q ) + { + for( integer dim=0; dim < 3; ++dim ) + { + // The default value is -1 so if it still -1 it needs to be set to something physical + if( permView[ei][q][dim] < 0 ) + { + permView[ei][q][dim] = permComponents[dim]; + permView_n[ei][q][dim] = permComponents[dim]; + } + } + } + } ); +} + +REGISTER_CATALOG_ENTRY( ConstitutiveBase, ReactivePressurePermeability, string const &, Group * const ) + +} +} /* namespace geos */ diff --git a/src/coreComponents/constitutive/permeability/ReactivePressurePermeability.hpp b/src/coreComponents/constitutive/permeability/ReactivePressurePermeability.hpp new file mode 100644 index 00000000000..14c3f080e19 --- /dev/null +++ b/src/coreComponents/constitutive/permeability/ReactivePressurePermeability.hpp @@ -0,0 +1,355 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file ReactivePressurePermeability.hpp + */ + +#ifndef GEOS_CONSTITUTIVE_PERMEABILITY_REACTIVEPRESSUREPERMEABILITY_HPP_ +#define GEOS_CONSTITUTIVE_PERMEABILITY_REACTIVEPRESSUREPERMEABILITY_HPP_ + +#include "constitutive/permeability/PermeabilityBase.hpp" +#include "constitutive/permeability/PressurePermeability.hpp" + + +namespace geos +{ +namespace constitutive +{ + +class ReactivePressurePermeabilityUpdate : public PermeabilityBaseUpdate +{ +public: + + ReactivePressurePermeabilityUpdate( PressureModelType const & presModelType, + R1Tensor const pressureDependenceConstants, + arrayView1d< real64 const > const & referencePressure, + real64 const & maxPermeability, + real64 const & minPermeability, + arrayView3d< real64 const > const & referencePermeability, + arrayView1d< real64 const > const & referencePorosity, + arrayView3d< real64 > const & permeability, + arrayView3d< real64 const > const & permeability_n, + arrayView3d< real64 > const & dPerm_dPressure, + integer const & explicitFlag ) + : PermeabilityBaseUpdate( permeability, dPerm_dPressure ), + m_presModelType( presModelType ), + m_pressureDependenceConstants( pressureDependenceConstants ), + m_referencePressure( referencePressure ), + m_maxPermeability( maxPermeability ), + m_minPermeability( minPermeability ), + m_referencePermeability( referencePermeability ), + m_referencePorosity( referencePorosity ), + m_permeability_n( permeability_n ), + m_explicitFlag( explicitFlag ) + {} + + GEOS_HOST_DEVICE + void computePressure( real64 const & deltaPressure, + R1Tensor const pressureDependenceConstants, + real64 const (&referencePermeability)[3], + arraySlice1d< real64 > const & permeability, + arraySlice1d< real64 > const & dPerm_dPressure ) const; + + GEOS_HOST_DEVICE + void computePressure( real64 const & deltaPressure, + R1Tensor const pressureDependenceConstants, + real64 const (&referencePermeability)[3], + real64 const maxPermeability, + arraySlice1d< real64 > const & permeability, + arraySlice1d< real64 > const & dPerm_dPressure ) const; + + GEOS_HOST_DEVICE + void computeWithPorosity( real64 const & porosity, + real64 const & minPermeability, + real64 const & referencePorosity, + real64 const (&referencePermeability)[3], + arraySlice1d< real64 > const & permeability, + arraySlice1d< real64 > const & dPerm_dPressure ) const; + + GEOS_HOST_DEVICE + virtual void updateFromPressureAndPorosity( localIndex const k, + localIndex const q, + real64 const & pressure, + real64 const & pressure_n, + real64 const & porosity ) const override + { + GEOS_UNUSED_VAR( q ); + + real64 deltaPressure = (m_explicitFlag)? (pressure_n - m_referencePressure[k]):(pressure - m_referencePressure[k]); + + real64 referencePermeability[3]; + + referencePermeability[0] = m_referencePermeability[k][0][0]; + referencePermeability[1] = m_referencePermeability[k][0][1]; + referencePermeability[2] = m_referencePermeability[k][0][2]; + + real64 permeability_n[3]; + + permeability_n[0] = m_permeability_n[k][0][0]; + permeability_n[1] = m_permeability_n[k][0][1]; + permeability_n[2] = m_permeability_n[k][0][2]; + + switch( m_presModelType ) + { + case PressureModelType::Exponential: + { + // Step 1: Update with pressure + computePressure( deltaPressure, + m_pressureDependenceConstants, + referencePermeability, + m_permeability[k][0], + m_dPerm_dPressure[k][0] ); + + // Step 2: Update with reactive porosity + computeWithPorosity( porosity, m_minPermeability, m_referencePorosity[k], referencePermeability, m_permeability[k][0], m_dPerm_dPressure[k][0] ); + + if( m_explicitFlag ) + { + for( localIndex i=0; i < m_permeability[k][0].size(); i++ ) + { + m_dPerm_dPressure[k][0][i] = 0.0; + + if( m_maxPermeability < 1.0 ) + { + real64 perm = m_permeability[k][0][i]; + + if( perm < referencePermeability[i] ) + { + perm = referencePermeability[i]; + } + + // Ensure permeability non-decreasing + if( perm < permeability_n[i] ) + { + perm = permeability_n[i]; + } + + m_permeability[k][0][i] = (perm > m_maxPermeability)? m_maxPermeability:perm; + } + } + } + + break; + } + case PressureModelType::Hyperbolic: + { + // Step 1: Update with pressure + computePressure( deltaPressure, + m_pressureDependenceConstants, + referencePermeability, + m_maxPermeability, + m_permeability[k][0], + m_dPerm_dPressure[k][0] ); + + // Step 2: Update with reactive porosity + computeWithPorosity( porosity, m_minPermeability, m_referencePorosity[k], referencePermeability, m_permeability[k][0], m_dPerm_dPressure[k][0] ); + + if( m_explicitFlag ) + { + for( localIndex i=0; i < m_permeability[k][0].size(); i++ ) + { + m_dPerm_dPressure[k][0][i] = 0.0; + + real64 const perm = m_permeability[k][0][i]; + + if( perm < referencePermeability[i] ) + { + m_permeability[k][0][i] = referencePermeability[i]; + } + } + } + break; + } + default: + { + GEOS_ERROR( "PressureModelType is invalid! It should be either Exponential or Hyperbolic" ); + } + } + } + +private: + + /// Pressure dependence model type + PressureModelType m_presModelType; + + /// Pressure dependent coefficients for each permeability component + R1Tensor m_pressureDependenceConstants; + + /// Reference pressure in the model + arrayView1d< real64 const > const m_referencePressure; + + /// Maximum and minimum permeability + real64 const m_maxPermeability; + real64 const m_minPermeability; + + arrayView3d< real64 const > const m_referencePermeability; + arrayView1d< real64 const > const m_referencePorosity; + + arrayView3d< real64 const > const m_permeability_n; + + integer m_explicitFlag; + +}; + + +class ReactivePressurePermeability : public PermeabilityBase +{ +public: + + ReactivePressurePermeability( string const & name, dataRepository::Group * const parent ); + + static string catalogName() { return "ReactivePressurePermeability"; } + + virtual string getCatalogName() const override { return catalogName(); } + + virtual void allocateConstitutiveData( dataRepository::Group & parent, + localIndex const numPts ) override; + + /// Type of kernel wrapper for in-kernel update + using KernelWrapper = ReactivePressurePermeabilityUpdate; + + /** + * @brief Create an update kernel wrapper. + * @return the wrapper + */ + KernelWrapper createKernelWrapper() const + { + return KernelWrapper( m_presModelType, + m_pressureDependenceConstants, + m_referencePressure, + m_maxPermeability, + m_minPermeability, + m_referencePermeability, + m_referencePorosity, + m_permeability, + m_permeability_n, + m_dPerm_dPressure, + m_explicitFlag ); + } + + struct viewKeyStruct : public PermeabilityBase::viewKeyStruct + { + static constexpr char const * referencePermeabilityComponentsString() { return "referencePermeabilityComponents"; } + static constexpr char const * pressureDependenceConstantsString() { return "pressureDependenceConstants"; } + static constexpr char const * defaultReferencePressureString() { return "defaultReferencePressure"; } + static constexpr char const * referencePressureString() { return "referencePressure"; } + static constexpr char const * referencePermeabilityString() { return "referencePermeability"; } + static constexpr char const * defaultReferencePorosityString() { return "defaultReferencePorosity"; } + static constexpr char const * referencePorosityString() { return "referencePorosity"; } + static constexpr char const * maxPermeabilityString() { return "maxPermeability"; } + static constexpr char const * minPermeabilityString() { return "minPermeability"; } + static constexpr char const * pressureModelTypeString() { return "pressureModelType"; } + static constexpr char const * explicitUpdateFlagString() { return "explicitUpdateFlag"; } + }; + + virtual void initializeState() const override final; + +protected: + + virtual void postInputInitialization() override; + +private: + + /// Permeability components at the reference pressure + R1Tensor m_referencePermeabilityComponents; + + real64 m_defaultReferencePorosity; + + /// Pressure dependent coefficients for each permeability component + R1Tensor m_pressureDependenceConstants; + + /// Reference pressure in the model + real64 m_defaultReferencePressure; + array1d< real64 > m_referencePressure; + + /// Maximum and minimum permeability + real64 m_maxPermeability; + real64 m_minPermeability; + + array3d< real64 > m_referencePermeability; + array1d< real64 > m_referencePorosity; + + /// Pressure dependence model type + PressureModelType m_presModelType; + + /// Explicit update flag + integer m_explicitFlag; + +}; + +GEOS_HOST_DEVICE +GEOS_FORCE_INLINE +void ReactivePressurePermeabilityUpdate::computePressure( real64 const & deltaPressure, + R1Tensor const pressureDependenceConstants, + real64 const (&referencePermeability)[3], + arraySlice1d< real64 > const & permeability, + arraySlice1d< real64 > const & dPerm_dPressure ) const +{ + for( localIndex i=0; i < permeability.size(); i++ ) + { + real64 const perm = referencePermeability[i] * exp( pressureDependenceConstants[i] * deltaPressure ); + + permeability[i] = perm; + dPerm_dPressure[i] = perm * pressureDependenceConstants[i]; + } +} + +GEOS_HOST_DEVICE +GEOS_FORCE_INLINE +void ReactivePressurePermeabilityUpdate::computePressure( real64 const & deltaPressure, + R1Tensor const pressureDependenceConstants, + real64 const (&referencePermeability)[3], + real64 const maxPermeability, + arraySlice1d< real64 > const & permeability, + arraySlice1d< real64 > const & dPerm_dPressure ) const +{ + for( localIndex i=0; i < permeability.size(); i++ ) + { + real64 const pressureOffSet = log( maxPermeability/referencePermeability[i] - 1 )/pressureDependenceConstants[i]; + + real64 const perm = maxPermeability/( 1 + exp( -pressureDependenceConstants[i]*( deltaPressure - pressureOffSet ) ) ); + permeability[i] = perm; + dPerm_dPressure[i] = perm*perm/maxPermeability*pressureDependenceConstants[i]*exp( -pressureDependenceConstants[i]*deltaPressure ); + } +} + +GEOS_HOST_DEVICE +GEOS_FORCE_INLINE +void ReactivePressurePermeabilityUpdate::computeWithPorosity( real64 const & porosity, + real64 const & minPermeability, + real64 const & referencePorosity, + real64 const (&referencePermeability)[3], + arraySlice1d< real64 > const & permeability, + arraySlice1d< real64 > const & dPerm_dPressure ) const +{ + for( localIndex i=0; i < permeability.size(); i++ ) + { + real64 const permScale = (minPermeability + ( referencePermeability[i] - minPermeability ) * pow( porosity/referencePorosity, 6 )) / referencePermeability[i]; + + real64 const perm = permeability[i]; + permeability[i] = perm * permScale; + + real64 const dPerm_dPres = dPerm_dPressure[i]; + dPerm_dPressure[i] = dPerm_dPres * permScale; + } +} + +}/* namespace constitutive */ + +} /* namespace geos */ + + +#endif //GEOS_CONSTITUTIVE_PERMEABILITY_REACTIVEPRESSUREPERMEABILITY_HPP_ diff --git a/src/coreComponents/constitutive/solid/CompressibleSolid.hpp b/src/coreComponents/constitutive/solid/CompressibleSolid.hpp index fb1832b2804..e835088b4d9 100644 --- a/src/coreComponents/constitutive/solid/CompressibleSolid.hpp +++ b/src/coreComponents/constitutive/solid/CompressibleSolid.hpp @@ -55,11 +55,12 @@ class CompressibleSolidUpdates : public CoupledSolidUpdates< NullModel, PORO_TYP virtual void updateStateFromPressureAndTemperature( localIndex const k, localIndex const q, real64 const & pressure, + real64 const & pressure_n, real64 const & temperature ) const override final { m_porosityUpdate.updateFromPressureAndTemperature( k, q, pressure, temperature ); real64 const porosity = m_porosityUpdate.getPorosity( k, q ); - m_permUpdate.updateFromPressureAndPorosity( k, q, pressure, porosity ); + m_permUpdate.updateFromPressureAndPorosity( k, q, pressure, pressure_n, porosity ); } GEOS_HOST_DEVICE diff --git a/src/coreComponents/constitutive/solid/CoupledSolid.hpp b/src/coreComponents/constitutive/solid/CoupledSolid.hpp index da74fb42ce8..a734eee7442 100644 --- a/src/coreComponents/constitutive/solid/CoupledSolid.hpp +++ b/src/coreComponents/constitutive/solid/CoupledSolid.hpp @@ -22,6 +22,7 @@ #define GEOS_CONSTITUTIVE_SOLID_COUPLEDSOLID_HPP_ #include "constitutive/solid/CoupledSolidBase.hpp" +#include "constitutive/fluid/reactivefluid/ReactiveFluidLayouts.hpp" namespace geos { @@ -85,9 +86,10 @@ class CoupledSolidUpdates virtual void updateStateFromPressureAndTemperature( localIndex const k, localIndex const q, real64 const & pressure, + real64 const & pressure_n, real64 const & temperature ) const { - GEOS_UNUSED_VAR( k, q, pressure, temperature ); + GEOS_UNUSED_VAR( k, q, pressure, pressure_n, temperature ); } GEOS_HOST_DEVICE @@ -105,6 +107,43 @@ class CoupledSolidUpdates temperature, temperature_k, temperature_n ); } + GEOS_HOST_DEVICE + virtual void updateStateReactionsFixedStress( localIndex const k, + localIndex const q, + real64 const & pressure, + real64 const & pressure_k, + real64 const & pressure_n, + real64 const & temperature, + real64 const & temperature_k, + real64 const & temperature_n, + arraySlice1d< real64 const, compflow::USD_COMP - 1 > mineralReactionMolarIncrements ) const + { + GEOS_UNUSED_VAR( k, q, + pressure, pressure_k, pressure_n, + temperature, temperature_k, temperature_n, + mineralReactionMolarIncrements ); + } + + GEOS_HOST_DEVICE + virtual void updateStateFromPressureTemperatureAndReactions( localIndex const k, + localIndex const q, + real64 const & pressure, + real64 const & pressure_n, + real64 const & temperature, + arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & kineticReactionMolarIncrements ) const + { + GEOS_UNUSED_VAR( k, q, pressure, pressure_n, temperature, kineticReactionMolarIncrements ); + } + + GEOS_HOST_DEVICE + virtual void updateSurfaceArea( localIndex const k, + localIndex const q, + arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & initialSurfaceArea, + arraySlice1d< real64, compflow::USD_COMP - 1 > const & surfaceArea ) const + { + GEOS_UNUSED_VAR( k, q, initialSurfaceArea, surfaceArea ); + } + GEOS_HOST_DEVICE virtual real64 getShearModulus( localIndex const k ) const { diff --git a/src/coreComponents/constitutive/solid/CoupledSolidBase.hpp b/src/coreComponents/constitutive/solid/CoupledSolidBase.hpp index a5081d2b7eb..b76e3483026 100644 --- a/src/coreComponents/constitutive/solid/CoupledSolidBase.hpp +++ b/src/coreComponents/constitutive/solid/CoupledSolidBase.hpp @@ -217,6 +217,7 @@ class CoupledSolidBase : public ConstitutiveBase virtual void saveConvergedState() const override final { getBasePorosityModel().saveConvergedState(); + getBasePermModel().saveConvergedState(); if( !m_solidInternalEnergyModelName.empty() ) { /// If the name is provided it has to be saved as well. diff --git a/src/coreComponents/constitutive/solid/EigenstrainReactiveSolid.cpp b/src/coreComponents/constitutive/solid/EigenstrainReactiveSolid.cpp new file mode 100644 index 00000000000..034fc33e011 --- /dev/null +++ b/src/coreComponents/constitutive/solid/EigenstrainReactiveSolid.cpp @@ -0,0 +1,63 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + + +/** + * @file EigenstrainReactiveSolid.cpp + */ + +#include "EigenstrainReactiveSolid.hpp" +#include "ElasticIsotropic.hpp" +#include "constitutive/permeability/ConstantPermeability.hpp" +#include "constitutive/permeability/CarmanKozenyPermeability.hpp" +#include "constitutive/permeability/PressurePermeability.hpp" +#include "constitutive/permeability/ReactivePressurePermeability.hpp" + +namespace geos +{ + +using namespace dataRepository; + +namespace constitutive +{ + +template< typename SOLID_TYPE, + typename PERM_TYPE > +EigenstrainReactiveSolid< SOLID_TYPE, PERM_TYPE >::EigenstrainReactiveSolid( string const & name, Group * const parent ): + CoupledSolid< SOLID_TYPE, ReactivePorosityBase, PERM_TYPE >( name, parent ) +{} + +template< typename SOLID_TYPE, + typename PERM_TYPE > +void EigenstrainReactiveSolid< SOLID_TYPE, PERM_TYPE >::initializeState() const +{ + CoupledSolid< SOLID_TYPE, ReactivePorosityBase, PERM_TYPE >::initializeState(); +} + +// Register all EigenstrainReactiveSolid model types. +typedef EigenstrainReactiveSolid< ElasticIsotropic, ConstantPermeability > EigenStrainReactiveElasticIsotropicConstant; +typedef EigenstrainReactiveSolid< ElasticIsotropic, CarmanKozenyPermeability > EigenStrainReactiveElasticIsotropicCK; +typedef EigenstrainReactiveSolid< ElasticIsotropic, PressurePermeability > EigenStrainReactiveElasticIsotropicPressure; +typedef EigenstrainReactiveSolid< ElasticIsotropic, ReactivePressurePermeability > EigenStrainReactiveElasticIsotropicReactivePressure; + + +REGISTER_CATALOG_ENTRY( ConstitutiveBase, EigenStrainReactiveElasticIsotropicConstant, string const &, Group * const ) +REGISTER_CATALOG_ENTRY( ConstitutiveBase, EigenStrainReactiveElasticIsotropicCK, string const &, Group * const ) +REGISTER_CATALOG_ENTRY( ConstitutiveBase, EigenStrainReactiveElasticIsotropicPressure, string const &, Group * const ) +REGISTER_CATALOG_ENTRY( ConstitutiveBase, EigenStrainReactiveElasticIsotropicReactivePressure, string const &, Group * const ) + + +} +} /* namespace geos */ diff --git a/src/coreComponents/constitutive/solid/EigenstrainReactiveSolid.hpp b/src/coreComponents/constitutive/solid/EigenstrainReactiveSolid.hpp new file mode 100644 index 00000000000..5ee8e3d5fcc --- /dev/null +++ b/src/coreComponents/constitutive/solid/EigenstrainReactiveSolid.hpp @@ -0,0 +1,340 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + + +/** + * @file EigenstrainReactiveSolid.hpp + */ + +#ifndef GEOS_CONSTITUTIVE_SOLID_EIGENSTRAINREACTIVESOLID_HPP_ +#define GEOS_CONSTITUTIVE_SOLID_EIGENSTRAINREACTIVESOLID_HPP_ + +#include "constitutive/solid/CoupledSolid.hpp" +#include "constitutive/solid/porosity/ReactivePorosityBase.hpp" +#include "constitutive/solid/SolidBase.hpp" +#include "constitutive/permeability/ConstantPermeability.hpp" + +#include "constitutive/fluid/reactivefluid/ReactiveFluidLayouts.hpp" + +namespace geos +{ +namespace constitutive +{ + +/** + * @brief Provides kernel-callable constitutive update routines + * + * + * @tparam SOLID_TYPE type of the porosity model + */ +template< typename SOLID_TYPE, + typename PERM_TYPE > +class EigenstrainReactiveSolidUpdates : public CoupledSolidUpdates< SOLID_TYPE, ReactivePorosityBase, PERM_TYPE > +{ +public: + + using DiscretizationOps = typename SOLID_TYPE::KernelWrapper::DiscretizationOps; + + /** + * @brief Constructor + */ + EigenstrainReactiveSolidUpdates( SOLID_TYPE const & solidModel, + ReactivePorosityBase const & porosityModel, + PERM_TYPE const & permModel ): + CoupledSolidUpdates< SOLID_TYPE, ReactivePorosityBase, PERM_TYPE >( solidModel, porosityModel, permModel ) + {} + + GEOS_HOST_DEVICE + virtual void updateStateReactionsFixedStress( localIndex const k, + localIndex const q, + real64 const & pressure, + real64 const & pressure_k, + real64 const & pressure_n, + real64 const & temperature, + real64 const & temperature_k, + real64 const & temperature_n, + arraySlice1d< real64 const, compflow::USD_COMP - 1 > mineralReactionMolarIncrements ) const override final + { + updateSolidBulkModulus( k ); + + m_porosityUpdate.updateFixedStress( k, q, + pressure, pressure_k, pressure_n, + temperature, temperature_k, temperature_n, + mineralReactionMolarIncrements ); + + real64 const porosity = m_porosityUpdate.getPorosity( k, q ); + m_permUpdate.updateFromPressureAndPorosity( k, q, pressure, pressure_n, porosity ); + } + + GEOS_HOST_DEVICE + virtual void updateStateFromPressureTemperatureAndReactions( localIndex const k, + localIndex const q, + real64 const & pressure, + real64 const & pressure_n, + real64 const & temperature, + arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & kineticReactionMolarIncrements ) const override final + { + GEOS_UNUSED_VAR( temperature ); + + m_porosityUpdate.updateFromReactions( k, q, kineticReactionMolarIncrements ); + real64 const porosity = m_porosityUpdate.getPorosity( k, q ); + m_permUpdate.updateFromPressureAndPorosity( k, q, pressure, pressure_n, porosity ); + } + + GEOS_HOST_DEVICE + virtual void updateSurfaceArea( localIndex const k, + localIndex const q, + arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & initialSurfaceArea, + arraySlice1d< real64, compflow::USD_COMP - 1 > const & surfaceArea ) const override final + { + real64 const porosity = m_porosityUpdate.getPorosity( k, q ); + real64 const initialPorosity = m_porosityUpdate.getInitialPorosity( k, q ); + + for( integer r=0; r < initialSurfaceArea.size(); ++r ) + { + real64 const volumeFraction_r = m_porosityUpdate.getVolumeFractionForMineral( k, q, r ); + real64 const initialVolumeFraction_r = m_porosityUpdate.getInitialVolumeFractionForMineral( k, q, r ); + // surfaceArea[r] = initialSurfaceArea[r] * pow( volumeFraction_r / initialVolumeFraction_r, 2.0/3.0 ) + // * pow( porosity / initialPorosity, 2.0/3.0 ); + + if( volumeFraction_r - initialVolumeFraction_r < 0 ) // dissolution + { + surfaceArea[r] = initialSurfaceArea[r] * pow( volumeFraction_r / initialVolumeFraction_r, 2.0/3.0 ) + * pow( porosity / initialPorosity, 2.0/3.0 ); + } + else // precipitation + { + surfaceArea[r] = initialSurfaceArea[r] * pow( porosity / initialPorosity, 6. ); + } + } + } + + GEOS_HOST_DEVICE + void smallStrainUpdateChemoMechanicsFixedStress( localIndex const k, + localIndex const q, + real64 const & timeIncrement, + real64 const & pressure, + real64 const & pressure_n, + real64 const & temperature, + real64 const & temperature_n, + real64 const & referenceTemperature, + arraySlice1d< real64 const, compflow::USD_COMP - 1 > mineralReactionMolarIncrements, + real64 const ( &strainIncrement )[6], + real64 ( & totalStress )[6], + DiscretizationOps & stiffness ) const + { + GEOS_UNUSED_VAR( pressure_n, referenceTemperature ); + + real64 anelasticStrainIncrement = 0.0; + + for( integer r=0; r < mineralReactionMolarIncrements.size(); ++r ) + { + real64 const molarWeight = m_porosityUpdate.getMolarWeights( r ); + real64 const mineralDensity = m_porosityUpdate.getMineralDensities( r ); + + anelasticStrainIncrement -= mineralReactionMolarIncrements[r] * molarWeight/mineralDensity; + } + + // Compute total stress increment and its derivative + real64 const deltaTemperatureFromLastStep = temperature - temperature_n; + computeTotalStress( k, + q, + timeIncrement, + pressure, + deltaTemperatureFromLastStep, + anelasticStrainIncrement, + strainIncrement, + totalStress, + stiffness ); + } + + /** + * @brief Return the stiffness at a given element (small-strain interface) + * + * @note If the material model has a strain-dependent material stiffness (e.g. + * any plasticity, damage, or nonlinear elastic model) then this interface will + * not work. Users should instead use one of the interfaces where a strain + * tensor is provided as input. + * + * @param k the element number + * @param stiffness the stiffness array + */ + GEOS_HOST_DEVICE + inline + void getElasticStiffness( localIndex const k, localIndex const q, real64 ( & stiffness )[6][6] ) const + { + m_solidUpdate.getElasticStiffness( k, q, stiffness ); + } + + /** + * @brief Return the stiffness at a given element (small-strain interface) + * + * @param [in] k the element number + * @param [out] thermalExpansionCoefficient the thermal expansion coefficient + */ + GEOS_HOST_DEVICE + inline + void getThermalExpansionCoefficient( localIndex const k, real64 & thermalExpansionCoefficient ) const + { + thermalExpansionCoefficient = m_solidUpdate.getThermalExpansionCoefficient( k ); + } + +private: + + using CoupledSolidUpdates< SOLID_TYPE, ReactivePorosityBase, PERM_TYPE >::m_solidUpdate; + using CoupledSolidUpdates< SOLID_TYPE, ReactivePorosityBase, PERM_TYPE >::m_porosityUpdate; + using CoupledSolidUpdates< SOLID_TYPE, ReactivePorosityBase, PERM_TYPE >::m_permUpdate; + + GEOS_HOST_DEVICE + inline + void updateSolidBulkModulus( localIndex const k ) const + { + real64 const bulkModulus = m_solidUpdate.getBulkModulus( k ); + + m_porosityUpdate.updateSolidBulkModulus( k, bulkModulus ); + } + + GEOS_HOST_DEVICE + inline + void computeTotalStress( localIndex const k, + localIndex const q, + real64 const & timeIncrement, + real64 const & pressure, + real64 const & deltaTemperatureFromLastStep, + real64 const & anelasticStrainIncrement, + real64 const ( &strainIncrement )[6], + real64 ( & totalStress )[6], + DiscretizationOps & stiffness ) const + { + // For now, the model only used for the sequential fixed stress scheme + // So we ignore the derivatives wrt pressure and temperature + real64 const thermalExpansionCoefficient = m_solidUpdate.getThermalExpansionCoefficient( k ); + + real64 mechanicsStrainIncrement[6]{}; + mechanicsStrainIncrement[0] = strainIncrement[0] - thermalExpansionCoefficient * deltaTemperatureFromLastStep - anelasticStrainIncrement; + mechanicsStrainIncrement[1] = strainIncrement[1] - thermalExpansionCoefficient * deltaTemperatureFromLastStep - anelasticStrainIncrement; + mechanicsStrainIncrement[2] = strainIncrement[2] - thermalExpansionCoefficient * deltaTemperatureFromLastStep - anelasticStrainIncrement; + mechanicsStrainIncrement[3] = strainIncrement[3]; + mechanicsStrainIncrement[4] = strainIncrement[4]; + mechanicsStrainIncrement[5] = strainIncrement[5]; + + // Compute total stress increment and its derivative w.r.t. pressure + m_solidUpdate.smallStrainUpdate( k, + q, + timeIncrement, + mechanicsStrainIncrement, + totalStress, // first effective stress increment accumulated + stiffness ); + + // Add the contributions of pressure to the total stress + LvArray::tensorOps::symAddIdentity< 3 >( totalStress, -pressure ); + + // Compute effective stress increment for the porosity update + real64 const bulkModulus = m_solidUpdate.getBulkModulus( k ); + real64 const meanEffectiveStressIncrement = bulkModulus * ( mechanicsStrainIncrement[0] + mechanicsStrainIncrement[1] + mechanicsStrainIncrement[2] ); + + m_porosityUpdate.updateMeanEffectiveStressIncrement( k, q, meanEffectiveStressIncrement ); + } + +}; + +/** + * @brief EigenstrainReactiveSolidBase class used for dispatch of all Porous solids. + */ +class EigenstrainReactiveSolidBase +{}; + +/** + * @brief Class to represent a porous material for poromechanics simulations. + * It is used as an interface to access all constitutive models relative to the properties of a porous material. + * + * @tparam SOLID_TYPE type of solid model + */ +template< typename SOLID_TYPE, + typename PERM_TYPE > +class EigenstrainReactiveSolid : public CoupledSolid< SOLID_TYPE, ReactivePorosityBase, PERM_TYPE > +{ +public: + + /// Alias for ElasticIsotropicUpdates + using KernelWrapper = EigenstrainReactiveSolidUpdates< SOLID_TYPE, PERM_TYPE >; + + /** + * @brief Constructor + * @param name Object name + * @param parent Object's parent group + */ + EigenstrainReactiveSolid( string const & name, dataRepository::Group * const parent ); + + /** + * @brief Catalog name + * @return Static catalog string + */ + static string catalogName() + { + if constexpr ( std::is_same_v< PERM_TYPE, ConstantPermeability > ) // default case + { + return string( "EigenStrainReactive" ) + SOLID_TYPE::catalogName(); + } + else // special cases + { + return string( "EigenStrainReactive" ) + SOLID_TYPE::catalogName() + PERM_TYPE::catalogName(); + } + } + + /** + * @brief Get catalog name + * @return Catalog name string + */ + virtual string getCatalogName() const override { return catalogName(); } + + /** + * @brief Create a instantiation of the EigenstrainReactiveSolidUpdates class + * that refers to the data in this. + * @return An instantiation of EigenstrainReactiveSolidUpdates. + */ + KernelWrapper createKernelUpdates() const + { + return KernelWrapper( getSolidModel(), + getPorosityModel(), + getPermModel() ); + } + + /** + * @brief initialize the constitutive models fields. + */ + virtual void initializeState() const override final; + + /** + * @brief Const/non-mutable accessor for density + * @return Accessor + */ + arrayView2d< real64 const > const getDensity() const + { + return getSolidModel().getDensity(); + } + +private: + using CoupledSolid< SOLID_TYPE, ReactivePorosityBase, PERM_TYPE >::getSolidModel; + using CoupledSolid< SOLID_TYPE, ReactivePorosityBase, PERM_TYPE >::getPorosityModel; + using CoupledSolid< SOLID_TYPE, ReactivePorosityBase, PERM_TYPE >::getPermModel; +}; + + + +} +} /* namespace geos */ + +#endif /* GEOS_CONSTITUTIVE_SOLID_EIGENSTRAINREACTIVESOLID_HPP_ */ diff --git a/src/coreComponents/constitutive/solid/PorousReactiveSolid.cpp b/src/coreComponents/constitutive/solid/PorousReactiveSolid.cpp new file mode 100644 index 00000000000..95a655f1f77 --- /dev/null +++ b/src/coreComponents/constitutive/solid/PorousReactiveSolid.cpp @@ -0,0 +1,57 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + + +/** + * @file PorousReactiveSolid.cpp + */ + +#include "PorousReactiveSolid.hpp" +#include "ElasticIsotropic.hpp" +#include "constitutive/permeability/ConstantPermeability.hpp" +#include "constitutive/permeability/CarmanKozenyPermeability.hpp" + +namespace geos +{ + +using namespace dataRepository; + +namespace constitutive +{ + +template< typename SOLID_TYPE, + typename PERM_TYPE > +PorousReactiveSolid< SOLID_TYPE, PERM_TYPE >::PorousReactiveSolid( string const & name, Group * const parent ): + CoupledSolid< SOLID_TYPE, BiotReactivePorosity, PERM_TYPE >( name, parent ) +{} + +template< typename SOLID_TYPE, + typename PERM_TYPE > +void PorousReactiveSolid< SOLID_TYPE, PERM_TYPE >::initializeState() const +{ + CoupledSolid< SOLID_TYPE, BiotReactivePorosity, PERM_TYPE >::initializeState(); +} + +// Register all PorousReactiveSolid model types. +typedef PorousReactiveSolid< ElasticIsotropic, ConstantPermeability > PorousReactiveElasticIsotropicConstant; +typedef PorousReactiveSolid< ElasticIsotropic, CarmanKozenyPermeability > PorousReactiveElasticIsotropicCK; + + +REGISTER_CATALOG_ENTRY( ConstitutiveBase, PorousReactiveElasticIsotropicConstant, string const &, Group * const ) +REGISTER_CATALOG_ENTRY( ConstitutiveBase, PorousReactiveElasticIsotropicCK, string const &, Group * const ) + + +} +} /* namespace geos */ diff --git a/src/coreComponents/constitutive/solid/PorousReactiveSolid.hpp b/src/coreComponents/constitutive/solid/PorousReactiveSolid.hpp new file mode 100644 index 00000000000..723b5610d80 --- /dev/null +++ b/src/coreComponents/constitutive/solid/PorousReactiveSolid.hpp @@ -0,0 +1,348 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + + +/** + * @file PorousReactiveSolid.hpp + */ + +#ifndef GEOS_CONSTITUTIVE_SOLID_POROUSREACTIVESOLID_HPP_ +#define GEOS_CONSTITUTIVE_SOLID_POROUSREACTIVESOLID_HPP_ + +#include "constitutive/solid/CoupledSolid.hpp" +#include "constitutive/solid/porosity/BiotReactivePorosity.hpp" +#include "constitutive/solid/SolidBase.hpp" +#include "constitutive/permeability/ConstantPermeability.hpp" + +#include "constitutive/fluid/reactivefluid/ReactiveFluidLayouts.hpp" + +namespace geos +{ +namespace constitutive +{ + +/** + * @brief Provides kernel-callable constitutive update routines + * + * + * @tparam SOLID_TYPE type of the porosity model + */ +template< typename SOLID_TYPE, + typename PERM_TYPE > +class PorousReactiveSolidUpdates : public CoupledSolidUpdates< SOLID_TYPE, BiotReactivePorosity, PERM_TYPE > +{ +public: + + using DiscretizationOps = typename SOLID_TYPE::KernelWrapper::DiscretizationOps; + + /** + * @brief Constructor + */ + PorousReactiveSolidUpdates( SOLID_TYPE const & solidModel, + BiotReactivePorosity const & porosityModel, + PERM_TYPE const & permModel ): + CoupledSolidUpdates< SOLID_TYPE, BiotReactivePorosity, PERM_TYPE >( solidModel, porosityModel, permModel ) + {} + + GEOS_HOST_DEVICE + virtual void updateStateReactionsFixedStress( localIndex const k, + localIndex const q, + real64 const & pressure, + real64 const & pressure_k, + real64 const & pressure_n, + real64 const & temperature, + real64 const & temperature_k, + real64 const & temperature_n, + arraySlice1d< real64 const, compflow::USD_COMP - 1 > mineralReactionMolarIncrements ) const override final + { + updateBiotCoefficientAndAssignModuli( k ); + + m_porosityUpdate.updateFixedStress( k, q, + pressure, pressure_k, pressure_n, + temperature, temperature_k, temperature_n, + mineralReactionMolarIncrements ); + } + + GEOS_HOST_DEVICE + virtual void updateStateFromPressureTemperatureAndReactions( localIndex const k, + localIndex const q, + real64 const & pressure, + real64 const & pressure_n, + real64 const & temperature, + arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & kineticReactionMolarIncrements ) const override final + { + GEOS_UNUSED_VAR( temperature ); + + m_porosityUpdate.updateFromReactions( k, q, kineticReactionMolarIncrements ); + real64 const porosity = m_porosityUpdate.getPorosity( k, q ); + m_permUpdate.updateFromPressureAndPorosity( k, q, pressure, pressure_n, porosity ); + } + + GEOS_HOST_DEVICE + virtual void updateSurfaceArea( localIndex const k, + localIndex const q, + arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & initialSurfaceArea, + arraySlice1d< real64, compflow::USD_COMP - 1 > const & surfaceArea ) const override final + { + real64 const porosity = m_porosityUpdate.getPorosity( k, q ); + real64 const initialPorosity = m_porosityUpdate.getInitialPorosity( k, q ); + + for( integer r=0; r < initialSurfaceArea.size(); ++r ) + { + real64 const volumeFraction_r = m_porosityUpdate.getVolumeFractionForMineral( k, q, r ); + real64 const initialVolumeFraction_r = m_porosityUpdate.getInitialVolumeFractionForMineral( k, q, r ); + surfaceArea[r] = initialSurfaceArea[r] * pow( volumeFraction_r / initialVolumeFraction_r, 2.0/3.0 ) + * pow( porosity / initialPorosity, 2.0/3.0 ); + } + } + + GEOS_HOST_DEVICE + void smallStrainUpdateChemoMechanicsFixedStress( localIndex const k, + localIndex const q, + real64 const & timeIncrement, + real64 const & pressure, + real64 const & pressure_n, + real64 const & temperature, + real64 const & temperature_n, + real64 const & referenceTemperature, + arraySlice1d< real64 const, compflow::USD_COMP - 1 > mineralReactionMolarIncrements, + real64 const ( &strainIncrement )[6], + real64 ( & totalStress )[6], + DiscretizationOps & stiffness ) const + { + GEOS_UNUSED_VAR( referenceTemperature ); + + real64 anelasticStrainIncrement = 0.0; + + for( integer r=0; r < mineralReactionMolarIncrements.size(); ++r ) + { + real64 const molarWeight = m_porosityUpdate.getMolarWeights( r ); + real64 const mineralDensity = m_porosityUpdate.getMineralDensities( r ); + + anelasticStrainIncrement -= mineralReactionMolarIncrements[r] * molarWeight/mineralDensity; + } + + // Compute total stress increment and its derivative + real64 const deltaPressureFromLastStep = pressure - pressure_n; + real64 const deltaTemperatureFromLastStep = temperature - temperature_n; + computeTotalStress( k, + q, + timeIncrement, + pressure, + deltaPressureFromLastStep, + deltaTemperatureFromLastStep, + anelasticStrainIncrement, + strainIncrement, + totalStress, + stiffness ); + } + + /** + * @brief Return the stiffness at a given element (small-strain interface) + * + * @note If the material model has a strain-dependent material stiffness (e.g. + * any plasticity, damage, or nonlinear elastic model) then this interface will + * not work. Users should instead use one of the interfaces where a strain + * tensor is provided as input. + * + * @param k the element number + * @param stiffness the stiffness array + */ + GEOS_HOST_DEVICE + inline + void getElasticStiffness( localIndex const k, localIndex const q, real64 ( & stiffness )[6][6] ) const + { + m_solidUpdate.getElasticStiffness( k, q, stiffness ); + } + + /** + * @brief Return the stiffness at a given element (small-strain interface) + * + * @param [in] k the element number + * @param [out] thermalExpansionCoefficient the thermal expansion coefficient + */ + GEOS_HOST_DEVICE + inline + void getThermalExpansionCoefficient( localIndex const k, real64 & thermalExpansionCoefficient ) const + { + thermalExpansionCoefficient = m_solidUpdate.getThermalExpansionCoefficient( k ); + } + +private: + + using CoupledSolidUpdates< SOLID_TYPE, BiotReactivePorosity, PERM_TYPE >::m_solidUpdate; + using CoupledSolidUpdates< SOLID_TYPE, BiotReactivePorosity, PERM_TYPE >::m_porosityUpdate; + using CoupledSolidUpdates< SOLID_TYPE, BiotReactivePorosity, PERM_TYPE >::m_permUpdate; + + GEOS_HOST_DEVICE + inline + void updateBiotCoefficientAndAssignModuli( localIndex const k ) const + { + // This call is not general like this. + real64 const bulkModulus = m_solidUpdate.getBulkModulus( k ); + + m_porosityUpdate.updateBiotCoefficientAndAssignModuli( k, bulkModulus ); + } + + GEOS_HOST_DEVICE + inline + void computeTotalStress( localIndex const k, + localIndex const q, + real64 const & timeIncrement, + real64 const & pressure, + real64 const & deltaPressureFromLastStep, + real64 const & deltaTemperatureFromLastStep, + real64 const & anelasticStrainIncrement, + real64 const ( &strainIncrement )[6], + real64 ( & totalStress )[6], + DiscretizationOps & stiffness ) const + { + // For now, the model only used for the sequential fixed stress scheme + // So we ignore the derivatives wrt pressure and temperature + real64 const thermalExpansionCoefficient = m_solidUpdate.getThermalExpansionCoefficient( k ); + + real64 mechanicsStrainIncrement[6]{}; + mechanicsStrainIncrement[0] = strainIncrement[0] - thermalExpansionCoefficient * deltaTemperatureFromLastStep; + mechanicsStrainIncrement[1] = strainIncrement[1] - thermalExpansionCoefficient * deltaTemperatureFromLastStep; + mechanicsStrainIncrement[2] = strainIncrement[2] - thermalExpansionCoefficient * deltaTemperatureFromLastStep; + mechanicsStrainIncrement[3] = strainIncrement[3]; + mechanicsStrainIncrement[4] = strainIncrement[4]; + mechanicsStrainIncrement[5] = strainIncrement[5]; + + // Add the contributions of pore material stress/pressure + real64 const biotCoefficient = m_porosityUpdate.getBiotCoefficient( k ); + + // Compute total stress increment and its derivative w.r.t. pressure + m_solidUpdate.smallStrainUpdate( k, + q, + timeIncrement, + mechanicsStrainIncrement, + totalStress, // first effective stress increment accumulated + stiffness ); + + // Compute effective stress increment for the porosity update + real64 const bulkModulus = m_solidUpdate.getBulkModulus( k ); + real64 const meanEffectiveStressIncrement = bulkModulus * ( mechanicsStrainIncrement[0] + mechanicsStrainIncrement[1] + mechanicsStrainIncrement[2] ); + + m_porosityUpdate.updateMeanEffectiveStressIncrement( k, q, meanEffectiveStressIncrement ); + + // Update mineral pressure + real64 dMineralPres_dMeanEffStressIncre = 0.0; + m_porosityUpdate.updatePoreMineralPressure( k, q, + deltaPressureFromLastStep, + meanEffectiveStressIncrement, + anelasticStrainIncrement, + dMineralPres_dMeanEffStressIncre ); + + real64 const mineralPressure = m_porosityUpdate.getPoreMineralPressure( k ); + real64 const totalPorePressure = pressure + mineralPressure; + + // Add the contributions of pressure to the total stress + LvArray::tensorOps::symAddIdentity< 3 >( totalStress, -biotCoefficient * totalPorePressure ); + + // Add the contributions of mineral pressure to the stiffness + stiffness.m_bulkModulus = bulkModulus - biotCoefficient * dMineralPres_dMeanEffStressIncre * bulkModulus; + } + +}; + +/** + * @brief PorousReactiveSolidBase class used for dispatch of all Porous solids. + */ +class PorousReactiveSolidBase +{}; + +/** + * @brief Class to represent a porous material for poromechanics simulations. + * It is used as an interface to access all constitutive models relative to the properties of a porous material. + * + * @tparam SOLID_TYPE type of solid model + */ +template< typename SOLID_TYPE, + typename PERM_TYPE > +class PorousReactiveSolid : public CoupledSolid< SOLID_TYPE, BiotReactivePorosity, PERM_TYPE > +{ +public: + + /// Alias for ElasticIsotropicUpdates + using KernelWrapper = PorousReactiveSolidUpdates< SOLID_TYPE, PERM_TYPE >; + + /** + * @brief Constructor + * @param name Object name + * @param parent Object's parent group + */ + PorousReactiveSolid( string const & name, dataRepository::Group * const parent ); + + /** + * @brief Catalog name + * @return Static catalog string + */ + static string catalogName() + { + if constexpr ( std::is_same_v< PERM_TYPE, ConstantPermeability > ) // default case + { + return string( "PorousReactive" ) + SOLID_TYPE::catalogName(); + } + else // special cases + { + return string( "PorousReactive" ) + SOLID_TYPE::catalogName() + PERM_TYPE::catalogName(); + } + } + + /** + * @brief Get catalog name + * @return Catalog name string + */ + virtual string getCatalogName() const override { return catalogName(); } + + /** + * @brief Create a instantiation of the PorousReactiveSolidUpdates class + * that refers to the data in this. + * @return An instantiation of PorousReactiveSolidUpdates. + */ + KernelWrapper createKernelUpdates() const + { + return KernelWrapper( getSolidModel(), + getPorosityModel(), + getPermModel() ); + } + + /** + * @brief initialize the constitutive models fields. + */ + virtual void initializeState() const override final; + + /** + * @brief Const/non-mutable accessor for density + * @return Accessor + */ + arrayView2d< real64 const > const getDensity() const + { + return getSolidModel().getDensity(); + } + +private: + using CoupledSolid< SOLID_TYPE, BiotReactivePorosity, PERM_TYPE >::getSolidModel; + using CoupledSolid< SOLID_TYPE, BiotReactivePorosity, PERM_TYPE >::getPorosityModel; + using CoupledSolid< SOLID_TYPE, BiotReactivePorosity, PERM_TYPE >::getPermModel; +}; + + + +} +} /* namespace geos */ + +#endif /* GEOS_CONSTITUTIVE_SOLID_POROUSREACTIVESOLID_HPP_ */ diff --git a/src/coreComponents/constitutive/solid/PorousSolid.cpp b/src/coreComponents/constitutive/solid/PorousSolid.cpp index 5b871f186e6..54beec9ecbb 100644 --- a/src/coreComponents/constitutive/solid/PorousSolid.cpp +++ b/src/coreComponents/constitutive/solid/PorousSolid.cpp @@ -29,6 +29,7 @@ #include "DuvautLionsSolid.hpp" #include "constitutive/permeability/ConstantPermeability.hpp" #include "constitutive/permeability/CarmanKozenyPermeability.hpp" +#include "constitutive/permeability/PressurePermeability.hpp" namespace geos { @@ -72,6 +73,16 @@ typedef PorousSolid< DuvautLionsSolid< DruckerPrager >, CarmanKozenyPermeability typedef PorousSolid< DuvautLionsSolid< DruckerPragerExtended >, CarmanKozenyPermeability > PorousViscoDruckerPragerExtendedCK; typedef PorousSolid< DuvautLionsSolid< ModifiedCamClay >, CarmanKozenyPermeability > PorousViscoModifiedCamClayCK; typedef PorousSolid< ModifiedCamClay, CarmanKozenyPermeability > PorousModifiedCamClayCK; +typedef PorousSolid< ElasticIsotropic, PressurePermeability > PorousElasticIsotropicPressure; +typedef PorousSolid< ElasticTransverseIsotropic, PressurePermeability > PorousElasticTransverseIsotropicPressure; +typedef PorousSolid< ElasticOrthotropic, PressurePermeability > PorousElasticOrthotropicPressure; +typedef PorousSolid< DelftEgg, PressurePermeability > PorousDelftEggPressure; +typedef PorousSolid< DruckerPrager, PressurePermeability > PorousDruckerPragerPressure; +typedef PorousSolid< DruckerPragerExtended, PressurePermeability > PorousDruckerPragerExtendedPressure; +typedef PorousSolid< DuvautLionsSolid< DruckerPrager >, PressurePermeability > PorousViscoDruckerPragerPressure; +typedef PorousSolid< DuvautLionsSolid< DruckerPragerExtended >, PressurePermeability > PorousViscoDruckerPragerExtendedPressure; +typedef PorousSolid< DuvautLionsSolid< ModifiedCamClay >, PressurePermeability > PorousViscoModifiedCamClayPressure; +typedef PorousSolid< ModifiedCamClay, PressurePermeability > PorousModifiedCamClayPressure; REGISTER_CATALOG_ENTRY( ConstitutiveBase, PorousElasticIsotropicConstant, string const &, Group * const ) @@ -94,6 +105,16 @@ REGISTER_CATALOG_ENTRY( ConstitutiveBase, PorousModifiedCamClayCK, string const REGISTER_CATALOG_ENTRY( ConstitutiveBase, PorousViscoDruckerPragerCK, string const &, Group * const ) REGISTER_CATALOG_ENTRY( ConstitutiveBase, PorousViscoDruckerPragerExtendedCK, string const &, Group * const ) REGISTER_CATALOG_ENTRY( ConstitutiveBase, PorousViscoModifiedCamClayCK, string const &, Group * const ) +REGISTER_CATALOG_ENTRY( ConstitutiveBase, PorousElasticIsotropicPressure, string const &, Group * const ) +REGISTER_CATALOG_ENTRY( ConstitutiveBase, PorousElasticTransverseIsotropicPressure, string const &, Group * const ) +REGISTER_CATALOG_ENTRY( ConstitutiveBase, PorousElasticOrthotropicPressure, string const &, Group * const ) +REGISTER_CATALOG_ENTRY( ConstitutiveBase, PorousDelftEggPressure, string const &, Group * const ) +REGISTER_CATALOG_ENTRY( ConstitutiveBase, PorousDruckerPragerPressure, string const &, Group * const ) +REGISTER_CATALOG_ENTRY( ConstitutiveBase, PorousDruckerPragerExtendedPressure, string const &, Group * const ) +REGISTER_CATALOG_ENTRY( ConstitutiveBase, PorousModifiedCamClayPressure, string const &, Group * const ) +REGISTER_CATALOG_ENTRY( ConstitutiveBase, PorousViscoDruckerPragerPressure, string const &, Group * const ) +REGISTER_CATALOG_ENTRY( ConstitutiveBase, PorousViscoDruckerPragerExtendedPressure, string const &, Group * const ) +REGISTER_CATALOG_ENTRY( ConstitutiveBase, PorousViscoModifiedCamClayPressure, string const &, Group * const ) } diff --git a/src/coreComponents/constitutive/solid/PorousSolid.hpp b/src/coreComponents/constitutive/solid/PorousSolid.hpp index d7f27433787..d03f1477248 100644 --- a/src/coreComponents/constitutive/solid/PorousSolid.hpp +++ b/src/coreComponents/constitutive/solid/PorousSolid.hpp @@ -69,6 +69,9 @@ class PorousSolidUpdates : public CoupledSolidUpdates< SOLID_TYPE, BiotPorosity, m_porosityUpdate.updateFixedStress( k, q, pressure, pressure_k, pressure_n, temperature, temperature_k, temperature_n ); + + real64 const porosity = m_porosityUpdate.getPorosity( k, q ); // this porosity is actually not being used + m_permUpdate.updateFromPressureAndPorosity( k, q, pressure, pressure_n, porosity ); } GEOS_HOST_DEVICE @@ -119,6 +122,9 @@ class PorousSolidUpdates : public CoupledSolidUpdates< SOLID_TYPE, BiotPorosity, dPorosity_dPressure, dPorosity_dTemperature ); + // Compute permeability + m_permUpdate.updateFromPressureAndPorosity( k, q, pressure, pressure_n, porosity ); + // skip porosity update when doing poromechanics initialization if( performStressInitialization ) { diff --git a/src/coreComponents/constitutive/solid/ReactiveSolid.cpp b/src/coreComponents/constitutive/solid/ReactiveSolid.cpp index eb015de05f4..c8dd448ecab 100644 --- a/src/coreComponents/constitutive/solid/ReactiveSolid.cpp +++ b/src/coreComponents/constitutive/solid/ReactiveSolid.cpp @@ -19,7 +19,7 @@ */ #include "ReactiveSolid.hpp" -#include "porosity/ReactivePorosity.hpp" +#include "porosity/ReactivePorosityBase.hpp" #include "constitutive/permeability/ConstantPermeability.hpp" #include "constitutive/permeability/CarmanKozenyPermeability.hpp" #include "constitutive/permeability/PressurePermeability.hpp" @@ -43,9 +43,9 @@ template< typename PORO_TYPE, ReactiveSolid< PORO_TYPE, PERM_TYPE >::~ReactiveSolid() = default; // Register all ReactiveSolid model types. -typedef ReactiveSolid< ReactivePorosity, ConstantPermeability > ReactiveRockConstant; -typedef ReactiveSolid< ReactivePorosity, CarmanKozenyPermeability > ReactiveRockCK; -typedef ReactiveSolid< ReactivePorosity, PressurePermeability > ReactiveRockPressurePerm; +typedef ReactiveSolid< ReactivePorosityBase, ConstantPermeability > ReactiveRockConstant; +typedef ReactiveSolid< ReactivePorosityBase, CarmanKozenyPermeability > ReactiveRockCK; +typedef ReactiveSolid< ReactivePorosityBase, PressurePermeability > ReactiveRockPressurePerm; REGISTER_CATALOG_ENTRY( ConstitutiveBase, ReactiveRockConstant, string const &, Group * const ) REGISTER_CATALOG_ENTRY( ConstitutiveBase, ReactiveRockCK, string const &, Group * const ) diff --git a/src/coreComponents/constitutive/solid/ReactiveSolid.hpp b/src/coreComponents/constitutive/solid/ReactiveSolid.hpp index 5e06af85ef3..44fcd00858f 100644 --- a/src/coreComponents/constitutive/solid/ReactiveSolid.hpp +++ b/src/coreComponents/constitutive/solid/ReactiveSolid.hpp @@ -22,7 +22,7 @@ #define GEOS_CONSTITUTIVE_SOLID_REACTIVESOLID_HPP_ #include "constitutive/solid/CoupledSolid.hpp" -#include "constitutive/solid/porosity/ReactivePorosity.hpp" +#include "constitutive/solid/porosity/ReactivePorosityBase.hpp" #include "constitutive/NullModel.hpp" #include "constitutive/fluid/reactivefluid/ReactiveFluidLayouts.hpp" @@ -55,21 +55,25 @@ class ReactiveSolidUpdates : public CoupledSolidUpdates< NullModel, PORO_TYPE, P {} GEOS_HOST_DEVICE - void updateStateFromPressureAndReactions( localIndex const k, - localIndex const q, - real64 const & pressure, - arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & kineticReactionMolarIncrements ) const + virtual void updateStateFromPressureTemperatureAndReactions( localIndex const k, + localIndex const q, + real64 const & pressure, + real64 const & pressure_n, + real64 const & temperature, + arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & kineticReactionMolarIncrements ) const override final { + GEOS_UNUSED_VAR( temperature ); + m_porosityUpdate.updateFromReactions( k, q, kineticReactionMolarIncrements ); real64 const porosity = m_porosityUpdate.getPorosity( k, q ); - m_permUpdate.updateFromPressureAndPorosity( k, q, pressure, porosity ); + m_permUpdate.updateFromPressureAndPorosity( k, q, pressure, pressure_n, porosity ); } GEOS_HOST_DEVICE - void updateSurfaceArea( localIndex const k, - localIndex const q, - arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & initialSurfaceArea, - arraySlice1d< real64, compflow::USD_COMP - 1 > const & surfaceArea ) const + virtual void updateSurfaceArea( localIndex const k, + localIndex const q, + arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & initialSurfaceArea, + arraySlice1d< real64, compflow::USD_COMP - 1 > const & surfaceArea ) const override final { real64 const porosity = m_porosityUpdate.getPorosity( k, q ); real64 const initialPorosity = m_porosityUpdate.getInitialPorosity( k, q ); @@ -84,9 +88,9 @@ class ReactiveSolidUpdates : public CoupledSolidUpdates< NullModel, PORO_TYPE, P } private: - using CoupledSolidUpdates< NullModel, ReactivePorosity, PERM_TYPE >::m_solidUpdate; - using CoupledSolidUpdates< NullModel, ReactivePorosity, PERM_TYPE >::m_porosityUpdate; - using CoupledSolidUpdates< NullModel, ReactivePorosity, PERM_TYPE >::m_permUpdate; + using CoupledSolidUpdates< NullModel, ReactivePorosityBase, PERM_TYPE >::m_solidUpdate; + using CoupledSolidUpdates< NullModel, ReactivePorosityBase, PERM_TYPE >::m_porosityUpdate; + using CoupledSolidUpdates< NullModel, ReactivePorosityBase, PERM_TYPE >::m_permUpdate; }; diff --git a/src/coreComponents/constitutive/solid/porosity/BiotReactivePorosity.cpp b/src/coreComponents/constitutive/solid/porosity/BiotReactivePorosity.cpp new file mode 100644 index 00000000000..fa8beac4dde --- /dev/null +++ b/src/coreComponents/constitutive/solid/porosity/BiotReactivePorosity.cpp @@ -0,0 +1,93 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file BiotReactivePorosity.cpp + */ + +#include "BiotReactivePorosity.hpp" +#include "PorosityFields.hpp" +#include "constitutive/solid/SolidBase.hpp" + +namespace geos +{ + +using namespace dataRepository; + +namespace constitutive +{ + +BiotReactivePorosity::BiotReactivePorosity( string const & name, Group * const parent ): + ReactivePorosityBase( name, parent ) +{ + registerWrapper( viewKeyStruct::defaultGrainBulkModulusString(), &m_defaultGrainBulkModulus ). + setInputFlag( InputFlags::REQUIRED ). + setApplyDefaultValue( -1.0 ). + setDescription( "Default grain bulk modulus" ); + + registerWrapper( viewKeyStruct::defaultMineralBulkModulusString(), &m_defaultMineralBulkModulus ). + setInputFlag( InputFlags::REQUIRED ). + setApplyDefaultValue( -1.0 ). + setDescription( "Default mineral bulk modulus" ); + + registerWrapper( viewKeyStruct::mineralBulkModulusString(), &m_mineralBulkModulus ). + setApplyDefaultValue( 0.0 ). + setPlotLevel( PlotLevel::LEVEL_0 ). + setDescription( "Mineral bulk modulus" ); + + registerWrapper( viewKeyStruct::mineralPressureString(), &m_mineralPressure ). + setApplyDefaultValue( 0.0 ). + setPlotLevel( PlotLevel::LEVEL_0 ). + setDescription( "Current mineral pressure" ); + + registerWrapper( viewKeyStruct::mineralPressure_nString(), &m_mineralPressure_n ). + setApplyDefaultValue( 0.0 ). + setPlotLevel( PlotLevel::LEVEL_0 ). + setDescription( "Mineral pressure at last time step" ); + + registerField< fields::porosity::biotCoefficient >( &m_biotCoefficient ); + + registerField< fields::porosity::grainBulkModulus >( &m_grainBulkModulus ); +} + +void BiotReactivePorosity::postInputInitialization() +{ + ReactivePorosityBase::postInputInitialization(); + + // set results as array default values + getWrapper< array1d< real64 > >( fields::porosity::grainBulkModulus::key() ). + setApplyDefaultValue( m_defaultGrainBulkModulus ); + + getWrapper< array1d< real64 > >( viewKeyStruct::mineralBulkModulusString() ). + setApplyDefaultValue( m_defaultMineralBulkModulus ); +} + +void BiotReactivePorosity::initializeState() const +{ + ReactivePorosityBase::initializeState(); + + m_mineralPressure_n.setValues< parallelDevicePolicy<> >( m_mineralPressure.toViewConst() ); +} + +void BiotReactivePorosity::saveConvergedState() const +{ + ReactivePorosityBase::saveConvergedState(); + + m_mineralPressure_n.setValues< parallelDevicePolicy<> >( m_mineralPressure.toViewConst() ); +} + +REGISTER_CATALOG_ENTRY( ConstitutiveBase, BiotReactivePorosity, string const &, Group * const ) +} /* namespace constitutive */ +} /* namespace geos */ diff --git a/src/coreComponents/constitutive/solid/porosity/BiotReactivePorosity.hpp b/src/coreComponents/constitutive/solid/porosity/BiotReactivePorosity.hpp new file mode 100644 index 00000000000..bbb8adfbc07 --- /dev/null +++ b/src/coreComponents/constitutive/solid/porosity/BiotReactivePorosity.hpp @@ -0,0 +1,320 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file BiotReactivePorosity.hpp + */ + +#ifndef GEOS_CONSTITUTIVE_POROSITY_BIOTREACTIVEPOROSITY_HPP_ +#define GEOS_CONSTITUTIVE_POROSITY_BIOTREACTIVEPOROSITY_HPP_ + +#include "ReactivePorosityBase.hpp" +#include "LvArray/src/tensorOps.hpp" + +namespace geos +{ +namespace constitutive +{ + +class BiotReactivePorosityUpdates : public ReactivePorosityBaseUpdates +{ +public: + + BiotReactivePorosityUpdates( arrayView2d< real64 > const & newPorosity, + arrayView2d< real64 > const & porosity_n, + arrayView2d< real64 > const & dPorosity_dPressure, + arrayView2d< real64 > const & dPorosity_dTemperature, + arrayView2d< real64 > const & initialPorosity, + arrayView1d< real64 > const & referencePorosity, + arrayView3d< real64, reactivefluid::USD_SPECIES > const & volumeFractions, + arrayView3d< real64 const, reactivefluid::USD_SPECIES > const & initialVolumeFractions, + arrayView3d< real64 const, reactivefluid::USD_SPECIES > const & volumeFractions_n, + integer const numKineticReactions, + arrayView1d< real64 const > const & molarWeights, + arrayView1d< real64 const > const & mineralDensities, + arrayView1d< real64 > const & biotCoefficient, + arrayView1d< real64 > const & bulkModulus, + arrayView1d< real64 > const & grainBulkModulus, + arrayView1d< real64 > const & mineralBulkModulus, + arrayView1d< real64 > const & mineralPressure, + arrayView1d< real64 > const & mineralPressure_n, + arrayView2d< real64 > const & meanEffectiveStressIncrement_k ): ReactivePorosityBaseUpdates( newPorosity, + porosity_n, + dPorosity_dPressure, + dPorosity_dTemperature, + initialPorosity, + referencePorosity, + volumeFractions, + initialVolumeFractions, + volumeFractions_n, + numKineticReactions, + molarWeights, + mineralDensities, + bulkModulus, + meanEffectiveStressIncrement_k ), + m_grainBulkModulus( grainBulkModulus ), + m_biotCoefficient( biotCoefficient ), + m_mineralBulkModulus( mineralBulkModulus ), + m_mineralPressure( mineralPressure ), + m_mineralPressure_n( mineralPressure_n ) + {} + + GEOS_HOST_DEVICE + real64 getBiotCoefficient( localIndex const k ) const { return m_biotCoefficient[k]; } + + GEOS_HOST_DEVICE + real64 getGrainBulkModulus( localIndex const k ) const { return m_grainBulkModulus[k]; } + + GEOS_HOST_DEVICE + real64 getPoreMineralPressure( localIndex const k ) const { return m_mineralPressure[k]; } + + GEOS_HOST_DEVICE + void computePorosityFixedStress( real64 const & pressure, + real64 const & pressure_k, + real64 const & pressure_n, + real64 const & porosity_n, + real64 const & referencePorosity, + real64 & porosity, + real64 & dPorosity_dPressure, + real64 const & biotCoefficient, + real64 const & meanEffectiveStressIncrement_k, + real64 const & bulkModulus, + real64 const & grainBulkModulus, + real64 const & mineralBulkModulus, + real64 const & reactionPorosityIncrement ) const + { + GEOS_UNUSED_VAR( pressure_k ); + + real64 const biotSkeletonModulusInverse = (biotCoefficient - referencePorosity) / grainBulkModulus; + real64 const porosityMultiplierInverse = 1 / ( 1 + biotSkeletonModulusInverse*mineralBulkModulus/referencePorosity ); + + porosity = porosity_n + // change due to stress increment + + biotCoefficient * meanEffectiveStressIncrement_k / bulkModulus * porosityMultiplierInverse + // change due to pressure increment + + biotSkeletonModulusInverse * ( pressure - pressure_n ) * porosityMultiplierInverse + // change due to mineral pressure increment + + biotSkeletonModulusInverse * reactionPorosityIncrement * mineralBulkModulus * porosityMultiplierInverse; + + dPorosity_dPressure = biotSkeletonModulusInverse * porosityMultiplierInverse; + } + + GEOS_HOST_DEVICE + void computePoreMineralPressure( real64 const & mineralPressure_n, + real64 & mineralPressure, + real64 & dMineralPres_dMeanEffStressIncre, + real64 const & referencePorosity, + real64 const & biotCoefficient, + real64 const & deltaPressureFromLastStep, + real64 const & meanEffectiveStressIncrement, + real64 const & bulkModulus, + real64 const & grainBulkModulus, + real64 const & mineralBulkModulus, + real64 const & anelasticStrainIncrement ) const + { + // GEOS_UNUSED_VAR( meanEffectiveStressIncrement, bulkModulus ); + + real64 const biotSkeletonModulusInverse = (biotCoefficient - referencePorosity) / grainBulkModulus; + real64 const mineralPressureMultiplier = mineralBulkModulus / ( referencePorosity + biotSkeletonModulusInverse*mineralBulkModulus ); + + mineralPressure = mineralPressure_n + // change due to inelastic strain increment + + mineralPressureMultiplier * referencePorosity * anelasticStrainIncrement + // change due to stress increment + - mineralPressureMultiplier * biotCoefficient * meanEffectiveStressIncrement / bulkModulus + // change due to pressure increment + - mineralPressureMultiplier * biotSkeletonModulusInverse * deltaPressureFromLastStep; + + dMineralPres_dMeanEffStressIncre = -mineralPressureMultiplier * biotCoefficient / bulkModulus; + // dMineralPres_dMeanEffStressIncre = 0.0; + } + + GEOS_HOST_DEVICE + void updatePoreMineralPressure( localIndex const k, + localIndex const q, + real64 const & deltaPressureFromLastStep, + real64 const & meanEffectiveStressIncrement, + real64 const & anelasticStrainIncrement, + real64 & dMineralPres_dMeanEffStressIncre ) const + { + GEOS_UNUSED_VAR( q ); + + computePoreMineralPressure( m_mineralPressure_n[k], + m_mineralPressure[k], + dMineralPres_dMeanEffStressIncre, + m_referencePorosity[k], + m_biotCoefficient[k], + deltaPressureFromLastStep, + meanEffectiveStressIncrement, + m_bulkModulus[k], + m_grainBulkModulus[k], + m_mineralBulkModulus[k], + anelasticStrainIncrement ); + } + + // this function is used in flow solver + // it uses average stress increment (element-based) + GEOS_HOST_DEVICE + virtual void updateFixedStress( localIndex const k, + localIndex const q, + real64 const & pressure, + real64 const & pressure_k, + real64 const & pressure_n, + real64 const & temperature, + real64 const & temperature_k, + real64 const & temperature_n, + arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & kineticReactionMolarIncrements ) const override final + { + // 1. Update the porosity due to reactions + real64 reactionPorosityIncrement = 0.0; + + ReactivePorosityBaseUpdates::computePorosityFromReaction( kineticReactionMolarIncrements, + m_volumeFractions[k][q], + m_volumeFractions_n[k][q], + reactionPorosityIncrement, + m_numKineticReactions, + m_molarWeights, + m_mineralDensities ); + + // 2. Update the porosity due to solid, pore mineral, and pore fluid pressure + // Currently ignore thermal effects + GEOS_UNUSED_VAR( temperature, temperature_k, temperature_n ); + computePorosityFixedStress( pressure, pressure_k, pressure_n, + m_porosity_n[k][q], + m_referencePorosity[k], + m_newPorosity[k][q], + m_dPorosity_dPressure[k][q], + m_biotCoefficient[k], + m_meanEffectiveStressIncrement_k[k][q], + m_bulkModulus[k], + m_grainBulkModulus[k], + m_mineralBulkModulus[k], + reactionPorosityIncrement ); + } + + GEOS_HOST_DEVICE + void updateBiotCoefficientAndAssignModuli( localIndex const k, + real64 const bulkModulus ) const + { + m_bulkModulus[k] = bulkModulus; + + m_biotCoefficient[k] = 1.0 - bulkModulus / m_grainBulkModulus[k]; + } + + + +protected: + + /// View on the grain bulk modulus (read from XML) + arrayView1d< real64 > const m_grainBulkModulus; + + /// View on the Biot coefficient (updated by PorousSolid) + arrayView1d< real64 > const m_biotCoefficient; + + /// View on the mineral bulk modulus (read from XML) + arrayView1d< real64 > const m_mineralBulkModulus; + + /// View on the mineral pressure + arrayView1d< real64 > const m_mineralPressure; + + /// View on the mineral pressure at the previous timestep + arrayView1d< real64 > const m_mineralPressure_n; +}; + +class BiotReactivePorosity : public ReactivePorosityBase +{ +public: + BiotReactivePorosity( string const & name, dataRepository::Group * const parent ); + + static string catalogName() { return "BiotReactivePorosity"; } + + virtual string getCatalogName() const override { return catalogName(); } + + struct viewKeyStruct : public ReactivePorosityBase::viewKeyStruct + { + static constexpr char const *defaultMineralBulkModulusString() { return "defaultMineralBulkModulus"; } + + static constexpr char const *mineralBulkModulusString() { return "mineralBulkModulus"; } + + static constexpr char const *mineralPressure_nString() { return "mineralPressure_n"; } + + static constexpr char const *mineralPressureString() { return "mineralPressure"; } + + static constexpr char const *defaultGrainBulkModulusString() { return "defaultGrainBulkModulus"; } + }; + + virtual void initializeState() const override final; + + virtual void saveConvergedState() const override final; + + using KernelWrapper = BiotReactivePorosityUpdates; + + /** + * @brief Create an update kernel wrapper. + * @return the wrapper + */ + KernelWrapper createKernelUpdates() const + { + return KernelWrapper( m_newPorosity, + m_porosity_n, + m_dPorosity_dPressure, + m_dPorosity_dTemperature, + m_initialPorosity, + m_referencePorosity, + m_volumeFractions, + m_initialVolumeFractions, + m_volumeFractions_n, + m_numKineticReactions, + m_molarWeights, + m_mineralDensities, + m_biotCoefficient, + m_bulkModulus, + m_grainBulkModulus, + m_mineralBulkModulus, + m_mineralPressure, + m_mineralPressure_n, + m_meanEffectiveStressIncrement_k ); + } + +protected: + virtual void postInputInitialization() override; + + /// Biot coefficients (update in the update class, not read in input) + array1d< real64 > m_biotCoefficient; + + /// Grain bulk modulus (read from XML) + real64 m_defaultGrainBulkModulus; + + /// Grain bulk modulus (can be specified in XML) + array1d< real64 > m_grainBulkModulus; + + /// Mineral bulk modulus (read from XML) + real64 m_defaultMineralBulkModulus; + + /// Mineral bulk modulus (can be specified in XML) + array1d< real64 > m_mineralBulkModulus; + + /// Mineral pressure + array1d< real64 > m_mineralPressure; + + /// Mineral pressure at the previous timestep + array1d< real64 > m_mineralPressure_n; +}; + +} /* namespace constitutive */ + +} /* namespace geos */ + +#endif //GEOS_CONSTITUTIVE_POROSITY_BIOTREACTIVEPOROSITY_HPP_ diff --git a/src/coreComponents/constitutive/solid/porosity/ReactivePorosity.hpp b/src/coreComponents/constitutive/solid/porosity/ReactivePorosity.hpp deleted file mode 100644 index 54bc69db6a7..00000000000 --- a/src/coreComponents/constitutive/solid/porosity/ReactivePorosity.hpp +++ /dev/null @@ -1,219 +0,0 @@ -/* - * ------------------------------------------------------------------------------------------------------------ - * SPDX-License-Identifier: LGPL-2.1-only - * - * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC - * Copyright (c) 2018-2024 TotalEnergies - * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University - * Copyright (c) 2023-2024 Chevron - * Copyright (c) 2019- GEOS/GEOSX Contributors - * All rights reserved - * - * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. - * ------------------------------------------------------------------------------------------------------------ - */ - -/** - * @file ReactivePorosity.hpp - */ - -#ifndef GEOS_CONSTITUTIVE_POROSITY_REACTIVEPOROSITY_HPP_ -#define GEOS_CONSTITUTIVE_POROSITY_REACTIVEPOROSITY_HPP_ - -#include "PorosityBase.hpp" - -#include "constitutive/fluid/reactivefluid/ReactiveFluidLayouts.hpp" - -namespace geos -{ -namespace constitutive -{ - -class ReactivePorosityUpdates : public PorosityBaseUpdates -{ -public: - - ReactivePorosityUpdates( arrayView2d< real64 > const & newPorosity, - arrayView2d< real64 const > const & porosity_n, - arrayView2d< real64 > const & dPorosity_dPressure, - arrayView2d< real64 > const & dPorosity_dTemperature, - arrayView2d< real64 const > const & initialPorosity, - arrayView1d< real64 const > const & referencePorosity, - arrayView3d< real64, reactivefluid::USD_SPECIES > const & volumeFractions, - arrayView3d< real64 const, reactivefluid::USD_SPECIES > const & initialVolumeFractions, - arrayView3d< real64 const, reactivefluid::USD_SPECIES > const & volumeFractions_n, - integer const numKineticReactions, - arrayView1d< real64 const > const & molarWeights, - arrayView1d< real64 const > const & mineralDensities ): - PorosityBaseUpdates( newPorosity, - porosity_n, - dPorosity_dPressure, - dPorosity_dTemperature, - initialPorosity, - referencePorosity ), - m_volumeFractions( volumeFractions ), - m_initialVolumeFractions( initialVolumeFractions ), - m_volumeFractions_n( volumeFractions_n ), - m_numKineticReactions( numKineticReactions ), - m_molarWeights( molarWeights ), - m_mineralDensities( mineralDensities ) - {} - - GEOS_HOST_DEVICE - void computePorosity( arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & kineticReactionMolarIncrements, - arraySlice1d< real64, reactivefluid::USD_SPECIES - 2 > const & volumeFractions, - arraySlice1d< real64 const, reactivefluid::USD_SPECIES - 2 > const & volumeFractions_n, - real64 & porosity, - real64 const & porosity_n, - integer const & numKineticReactions, - arrayView1d< real64 const > const & molarWeights, - arrayView1d< real64 const > const & mineralDensities ) const - { - real64 porosityIncrement = 0.0; - - for( integer r=0; r < numKineticReactions; ++r ) - { - real64 const volumeFractionIncrement = -kineticReactionMolarIncrements[r] * molarWeights[r]/mineralDensities[r]; - volumeFractions[r] = volumeFractions_n[r] + volumeFractionIncrement; - - porosityIncrement -= volumeFractionIncrement; - } - - porosity = porosity_n + porosityIncrement; - - if( porosity < 0 ) - { - porosity = 0; - } - else if( porosity > 1.0 ) - { - porosity = 1.0; - } - - } - - GEOS_HOST_DEVICE - void updateFromReactions( localIndex const k, - localIndex const q, - arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & kineticReactionMolarIncrements ) const - { - computePorosity( kineticReactionMolarIncrements, - m_volumeFractions[k][q], - m_volumeFractions_n[k][q], - m_newPorosity[k][q], - m_porosity_n[k][q], - m_numKineticReactions, - m_molarWeights, - m_mineralDensities ); - } - - GEOS_HOST_DEVICE - inline - real64 getVolumeFractionForMineral( localIndex const k, - localIndex const q, - localIndex const r ) const - { - return m_volumeFractions[k][q][r]; - } - - GEOS_HOST_DEVICE - inline - real64 getInitialVolumeFractionForMineral( localIndex const k, - localIndex const q, - localIndex const r ) const - { - return m_initialVolumeFractions[k][q][r]; - } - -protected: - - arrayView3d< real64, reactivefluid::USD_SPECIES > m_volumeFractions; - arrayView3d< real64 const, reactivefluid::USD_SPECIES > m_initialVolumeFractions; - arrayView3d< real64 const, reactivefluid::USD_SPECIES > const m_volumeFractions_n; - - integer const m_numKineticReactions; - arrayView1d< real64 const > const m_molarWeights; - arrayView1d< real64 const > const m_mineralDensities; -}; - - -class ReactivePorosity : public PorosityBase -{ -public: - ReactivePorosity( string const & name, Group * const parent ); - - virtual std::unique_ptr< ConstitutiveBase > - deliverClone( string const & name, - dataRepository::Group * const parent ) const override; - - virtual void allocateConstitutiveData( dataRepository::Group & parent, - localIndex const numConstitutivePointsPerParentIndex ) override; - - static string catalogName() { return "ReactivePorosity"; } - - virtual string getCatalogName() const override { return catalogName(); } - - virtual void saveConvergedState() const override; - - integer numKineticReactions() const { return m_numKineticReactions; } - - virtual void initializeState() const override; - - struct viewKeyStruct : public PorosityBase::viewKeyStruct - { - static constexpr char const * defaultInitialVolumeFractionsString() { return "defaultInitialVolumeFractions"; } - static constexpr char const * initialVolumeFractionsString() { return "initialVolumeFractions"; } - static constexpr char const * volumeFractionsString() { return "volumeFractions"; } - static constexpr char const * volumeFractions_nString() { return "volumeFractions_n"; } - static constexpr char const * molarWeightsString() { return "molarWeights"; } - static constexpr char const * mineralDensitiesString() { return "mineralDensities"; } - } viewKeys; - - - using KernelWrapper = ReactivePorosityUpdates; - - /** - * @brief Create an update kernel wrapper. - * @return the wrapper - */ - KernelWrapper createKernelUpdates() const - { - return KernelWrapper( m_newPorosity, - m_porosity_n, - m_dPorosity_dPressure, - m_dPorosity_dTemperature, - m_initialPorosity, - m_referencePorosity, - m_volumeFractions, - m_initialVolumeFractions, - m_volumeFractions_n, - m_numKineticReactions, - m_molarWeights, - m_mineralDensities ); - } - - -private: - virtual void postInputInitialization() override; - - virtual void resizeFields( localIndex const size, localIndex const numPts ); - - array1d< real64 > m_defaultInitialVolumeFractions; - - array3d< real64, constitutive::reactivefluid::LAYOUT_SPECIES > m_volumeFractions; - array3d< real64, constitutive::reactivefluid::LAYOUT_SPECIES > m_initialVolumeFractions; - array3d< real64, constitutive::reactivefluid::LAYOUT_SPECIES > m_volumeFractions_n; - - integer m_numKineticReactions; - array1d< real64 > m_molarWeights; - array1d< real64 > m_mineralDensities; - -}; - - -}/* namespace constitutive */ - -} /* namespace geos */ - - -#endif //GEOS_CONSTITUTIVE_POROSITY_REACTIVEPOROSITY_HPP_ diff --git a/src/coreComponents/constitutive/solid/porosity/ReactivePorosity.cpp b/src/coreComponents/constitutive/solid/porosity/ReactivePorosityBase.cpp similarity index 73% rename from src/coreComponents/constitutive/solid/porosity/ReactivePorosity.cpp rename to src/coreComponents/constitutive/solid/porosity/ReactivePorosityBase.cpp index 7865764fd7c..d75e3a19c48 100644 --- a/src/coreComponents/constitutive/solid/porosity/ReactivePorosity.cpp +++ b/src/coreComponents/constitutive/solid/porosity/ReactivePorosityBase.cpp @@ -14,10 +14,10 @@ */ /** - * @file ReactivePorosity.cpp + * @file ReactivePorosityBase.cpp */ -#include "ReactivePorosity.hpp" +#include "ReactivePorosityBase.hpp" namespace geos { @@ -27,7 +27,7 @@ using namespace dataRepository; namespace constitutive { -ReactivePorosity::ReactivePorosity( string const & name, Group * const parent ): +ReactivePorosityBase::ReactivePorosityBase( string const & name, Group * const parent ): PorosityBase( name, parent ) { registerWrapper( viewKeyStruct::defaultInitialVolumeFractionsString(), &m_defaultInitialVolumeFractions ). @@ -56,20 +56,29 @@ ReactivePorosity::ReactivePorosity( string const & name, Group * const parent ): registerWrapper( viewKeyStruct::mineralDensitiesString(), &m_mineralDensities ). setInputFlag( InputFlags::REQUIRED ). setDescription( "Mineral densities" ); + + registerWrapper( viewKeyStruct::solidBulkModulusString(), &m_bulkModulus ). + setApplyDefaultValue( 1e-6 ). + setDescription( "Solid bulk modulus" ); + + registerWrapper( viewKeyStruct::meanEffectiveStressIncrement_kString(), &m_meanEffectiveStressIncrement_k ). + setApplyDefaultValue( 0.0 ). + setDescription( "Mean effective stress increment at quadrature points at the previous sequential iteration" ); + } -std::unique_ptr< ConstitutiveBase > ReactivePorosity::deliverClone( string const & name, Group * const parent ) const +std::unique_ptr< ConstitutiveBase > ReactivePorosityBase::deliverClone( string const & name, Group * const parent ) const { std::unique_ptr< ConstitutiveBase > clone = ConstitutiveBase::deliverClone( name, parent ); - ReactivePorosity & newConstitutiveRelation = dynamicCast< ReactivePorosity & >( *clone ); + ReactivePorosityBase & newConstitutiveRelation = dynamicCast< ReactivePorosityBase & >( *clone ); newConstitutiveRelation.m_numKineticReactions = m_numKineticReactions; return clone; } -void ReactivePorosity::postInputInitialization() +void ReactivePorosityBase::postInputInitialization() { PorosityBase::postInputInitialization(); @@ -86,8 +95,8 @@ void ReactivePorosity::postInputInitialization() m_numKineticReactions = m_defaultInitialVolumeFractions.size(); } -void ReactivePorosity::allocateConstitutiveData( dataRepository::Group & parent, - localIndex const numConstitutivePointsPerParentIndex ) +void ReactivePorosityBase::allocateConstitutiveData( dataRepository::Group & parent, + localIndex const numConstitutivePointsPerParentIndex ) { PorosityBase::allocateConstitutiveData( parent, numConstitutivePointsPerParentIndex ); @@ -95,23 +104,33 @@ void ReactivePorosity::allocateConstitutiveData( dataRepository::Group & parent, } -void ReactivePorosity::resizeFields( localIndex const size, localIndex const numPts ) +void ReactivePorosityBase::resizeFields( localIndex const size, localIndex const numPts ) { integer const numKineticReactions = this->numKineticReactions(); m_initialVolumeFractions.resize( size, numPts, numKineticReactions ); m_volumeFractions.resize( size, numPts, numKineticReactions ); m_volumeFractions_n.resize( size, numPts, numKineticReactions ); + + m_meanEffectiveStressIncrement_k.resize( 0, numPts ); } -void ReactivePorosity::saveConvergedState() const +void ReactivePorosityBase::saveConvergedState() const { PorosityBase::saveConvergedState(); m_volumeFractions_n.setValues< parallelDevicePolicy<> >( m_volumeFractions.toViewConst() ); + m_meanEffectiveStressIncrement_k.zero(); } -void ReactivePorosity::initializeState() const +void ReactivePorosityBase::ignoreConvergedState() const +{ + PorosityBase::ignoreConvergedState(); + m_meanEffectiveStressIncrement_k.zero(); +} + + +void ReactivePorosityBase::initializeState() const { integer const numKineticReactions = this->numKineticReactions(); @@ -139,6 +158,6 @@ void ReactivePorosity::initializeState() const } } -REGISTER_CATALOG_ENTRY( ConstitutiveBase, ReactivePorosity, string const &, Group * const ) +REGISTER_CATALOG_ENTRY( ConstitutiveBase, ReactivePorosityBase, string const &, Group * const ) } } /* namespace geos */ diff --git a/src/coreComponents/constitutive/solid/porosity/ReactivePorosityBase.hpp b/src/coreComponents/constitutive/solid/porosity/ReactivePorosityBase.hpp new file mode 100644 index 00000000000..cf8be289bfe --- /dev/null +++ b/src/coreComponents/constitutive/solid/porosity/ReactivePorosityBase.hpp @@ -0,0 +1,301 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file ReactivePorosityBase.hpp + */ + +#ifndef GEOS_CONSTITUTIVE_POROSITY_REACTIVEPOROSITYBASE_HPP_ +#define GEOS_CONSTITUTIVE_POROSITY_REACTIVEPOROSITYBASE_HPP_ + +#include "PorosityBase.hpp" + +#include "constitutive/fluid/reactivefluid/ReactiveFluidLayouts.hpp" + +namespace geos +{ +namespace constitutive +{ + +class ReactivePorosityBaseUpdates : public PorosityBaseUpdates +{ +public: + + ReactivePorosityBaseUpdates( arrayView2d< real64 > const & newPorosity, + arrayView2d< real64 const > const & porosity_n, + arrayView2d< real64 > const & dPorosity_dPressure, + arrayView2d< real64 > const & dPorosity_dTemperature, + arrayView2d< real64 const > const & initialPorosity, + arrayView1d< real64 const > const & referencePorosity, + arrayView3d< real64, reactivefluid::USD_SPECIES > const & volumeFractions, + arrayView3d< real64 const, reactivefluid::USD_SPECIES > const & initialVolumeFractions, + arrayView3d< real64 const, reactivefluid::USD_SPECIES > const & volumeFractions_n, + integer const numKineticReactions, + arrayView1d< real64 const > const & molarWeights, + arrayView1d< real64 const > const & mineralDensities, + arrayView1d< real64 > const & bulkModulus, + arrayView2d< real64 > const & meanEffectiveStressIncrement_k ): + PorosityBaseUpdates( newPorosity, + porosity_n, + dPorosity_dPressure, + dPorosity_dTemperature, + initialPorosity, + referencePorosity ), + m_volumeFractions( volumeFractions ), + m_initialVolumeFractions( initialVolumeFractions ), + m_volumeFractions_n( volumeFractions_n ), + m_numKineticReactions( numKineticReactions ), + m_molarWeights( molarWeights ), + m_mineralDensities( mineralDensities ), + m_bulkModulus( bulkModulus ), + m_meanEffectiveStressIncrement_k( meanEffectiveStressIncrement_k ) + {} + + GEOS_HOST_DEVICE + void computePorosityFromReaction( arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & kineticReactionMolarIncrements, + arraySlice1d< real64, reactivefluid::USD_SPECIES - 2 > const & volumeFractions, + arraySlice1d< real64 const, reactivefluid::USD_SPECIES - 2 > const & volumeFractions_n, + real64 & reactionPorosityIncrement, + integer const & numKineticReactions, + arrayView1d< real64 const > const & molarWeights, + arrayView1d< real64 const > const & mineralDensities ) const + { + for( integer r=0; r < numKineticReactions; ++r ) + { + real64 const volumeFractionIncrement = -kineticReactionMolarIncrements[r] * molarWeights[r]/mineralDensities[r]; + volumeFractions[r] = volumeFractions_n[r] + volumeFractionIncrement; + + reactionPorosityIncrement -= volumeFractionIncrement; + } + } + + GEOS_HOST_DEVICE + void updateFromReactions( localIndex const k, + localIndex const q, + arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & kineticReactionMolarIncrements ) const + { + real64 reactionPorosityIncrement = 0.0; + + computePorosityFromReaction( kineticReactionMolarIncrements, + m_volumeFractions[k][q], + m_volumeFractions_n[k][q], + reactionPorosityIncrement, + m_numKineticReactions, + m_molarWeights, + m_mineralDensities ); + + m_newPorosity[k][q] = m_porosity_n[k][q] + reactionPorosityIncrement; + + if( m_newPorosity[k][q] < 0 ) + { + m_newPorosity[k][q] = 0; + } + else if( m_newPorosity[k][q] > 1.0 ) + { + m_newPorosity[k][q] = 1.0; + } + } + + // this function is used in flow solver + // it uses average stress increment (element-based) + GEOS_HOST_DEVICE + virtual void updateFixedStress( localIndex const k, + localIndex const q, + real64 const & pressure, + real64 const & pressure_k, + real64 const & pressure_n, + real64 const & temperature, + real64 const & temperature_k, + real64 const & temperature_n, + arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & kineticReactionMolarIncrements ) const + { + // For now, we ignore the pressure or temperature dependence but just follow Evans et al. for this eigenstrain approach + GEOS_UNUSED_VAR( pressure, pressure_k, pressure_n, temperature, temperature_k, temperature_n ); + + // 1. Update the porosity due to reactions + real64 reactionPorosityIncrement = 0.0; + + computePorosityFromReaction( kineticReactionMolarIncrements, + m_volumeFractions[k][q], + m_volumeFractions_n[k][q], + reactionPorosityIncrement, + m_numKineticReactions, + m_molarWeights, + m_mineralDensities ); + + // 2. Update the porosity due to solid deformation + m_newPorosity[k][q] = m_porosity_n[k][q] + + reactionPorosityIncrement; + // + m_meanEffectiveStressIncrement_k[k][q]/m_bulkModulus[k]; + + if( m_newPorosity[k][q] < 9e-4 ) + { + m_newPorosity[k][q] = 9e-4; + } + else if( m_newPorosity[k][q] > 1.0 ) + { + m_newPorosity[k][q] = 1.0; + } + } + + GEOS_HOST_DEVICE + void updateSolidBulkModulus( localIndex const k, + real64 const bulkModulus ) const + { + m_bulkModulus[k] = bulkModulus; + } + + GEOS_HOST_DEVICE + void updateMeanEffectiveStressIncrement( localIndex const k, + localIndex const q, + real64 const & meanEffectiveStressIncrement ) const + { + m_meanEffectiveStressIncrement_k[k][q] = meanEffectiveStressIncrement; + } + + GEOS_HOST_DEVICE + inline + real64 getVolumeFractionForMineral( localIndex const k, + localIndex const q, + localIndex const r ) const + { + return m_volumeFractions[k][q][r]; + } + + GEOS_HOST_DEVICE + inline + real64 getInitialVolumeFractionForMineral( localIndex const k, + localIndex const q, + localIndex const r ) const + { + return m_initialVolumeFractions[k][q][r]; + } + + GEOS_HOST_DEVICE + inline + real64 getMolarWeights( localIndex const r ) const + { + return m_molarWeights[r]; + } + + GEOS_HOST_DEVICE + inline + real64 getMineralDensities( localIndex const r ) const + { + return m_mineralDensities[r]; + } + +protected: + + arrayView3d< real64, reactivefluid::USD_SPECIES > m_volumeFractions; + arrayView3d< real64 const, reactivefluid::USD_SPECIES > m_initialVolumeFractions; + arrayView3d< real64 const, reactivefluid::USD_SPECIES > const m_volumeFractions_n; + + integer const m_numKineticReactions; + arrayView1d< real64 const > const m_molarWeights; + arrayView1d< real64 const > const m_mineralDensities; + + arrayView1d< real64 > const m_bulkModulus; + arrayView2d< real64 > const m_meanEffectiveStressIncrement_k; +}; + + +class ReactivePorosityBase : public PorosityBase +{ +public: + ReactivePorosityBase( string const & name, Group * const parent ); + + virtual std::unique_ptr< ConstitutiveBase > + deliverClone( string const & name, + dataRepository::Group * const parent ) const override; + + virtual void allocateConstitutiveData( dataRepository::Group & parent, + localIndex const numConstitutivePointsPerParentIndex ) override; + + static string catalogName() { return "ReactivePorosity"; } + + virtual string getCatalogName() const override { return catalogName(); } + + virtual void saveConvergedState() const override; + virtual void ignoreConvergedState() const override; + + integer numKineticReactions() const { return m_numKineticReactions; } + + virtual void initializeState() const override; + + struct viewKeyStruct : public PorosityBase::viewKeyStruct + { + static constexpr char const * defaultInitialVolumeFractionsString() { return "defaultInitialVolumeFractions"; } + static constexpr char const * initialVolumeFractionsString() { return "initialVolumeFractions"; } + static constexpr char const * volumeFractionsString() { return "volumeFractions"; } + static constexpr char const * volumeFractions_nString() { return "volumeFractions_n"; } + static constexpr char const * molarWeightsString() { return "molarWeights"; } + static constexpr char const * mineralDensitiesString() { return "mineralDensities"; } + static constexpr char const * solidBulkModulusString() { return "solidBulkModulus"; } + static constexpr char const * meanEffectiveStressIncrement_kString() { return "meanEffectiveStressIncrement_k"; } + } viewKeys; + + + using KernelWrapper = ReactivePorosityBaseUpdates; + + /** + * @brief Create an update kernel wrapper. + * @return the wrapper + */ + KernelWrapper createKernelUpdates() const + { + return KernelWrapper( m_newPorosity, + m_porosity_n, + m_dPorosity_dPressure, + m_dPorosity_dTemperature, + m_initialPorosity, + m_referencePorosity, + m_volumeFractions, + m_initialVolumeFractions, + m_volumeFractions_n, + m_numKineticReactions, + m_molarWeights, + m_mineralDensities, + m_bulkModulus, + m_meanEffectiveStressIncrement_k ); + } + + +protected: + virtual void postInputInitialization() override; + + virtual void resizeFields( localIndex const size, localIndex const numPts ); + + array1d< real64 > m_defaultInitialVolumeFractions; + + array3d< real64, constitutive::reactivefluid::LAYOUT_SPECIES > m_volumeFractions; + array3d< real64, constitutive::reactivefluid::LAYOUT_SPECIES > m_initialVolumeFractions; + array3d< real64, constitutive::reactivefluid::LAYOUT_SPECIES > m_volumeFractions_n; + + integer m_numKineticReactions; + array1d< real64 > m_molarWeights; + array1d< real64 > m_mineralDensities; + + array1d< real64 > m_bulkModulus; + array2d< real64 > m_meanEffectiveStressIncrement_k; +}; + + +}/* namespace constitutive */ + +} /* namespace geos */ + + +#endif //GEOS_CONSTITUTIVE_POROSITY_REACTIVEPOROSITYBASE_HPP_ diff --git a/src/coreComponents/mesh/PerforationData.cpp b/src/coreComponents/mesh/PerforationData.cpp index 560d728dcd9..e58d1eeb620 100644 --- a/src/coreComponents/mesh/PerforationData.cpp +++ b/src/coreComponents/mesh/PerforationData.cpp @@ -230,6 +230,12 @@ void PerforationData::computeWellTransmissibility( MeshLevel const & mesh, // compute the well Peaceman index m_wellTransmissibility[iperf] = 2 * M_PI * kh / ( std::log( rEq / wellElemRadius[wellElemIndex] ) + m_wellSkinFactor[iperf] ); + { + WellElementRegion const & wellRegion = dynamicCast< WellElementRegion const & >( wellElemSubRegion.getParent().getParent() ); + GEOS_LOG_RANK( "\n \nPerforation " << wellRegion.getWellGeneratorName() << + " has transmissibility of " << m_wellTransmissibility[iperf] << " \n \n" ); + } + if( m_wellTransmissibility[iperf] <= 0 ) { WellElementRegion const & wellRegion = dynamicCast< WellElementRegion const & >( wellElemSubRegion.getParent().getParent() ); diff --git a/src/coreComponents/mesh/generators/WellGeneratorBase.cpp b/src/coreComponents/mesh/generators/WellGeneratorBase.cpp index ee3118133a8..f4c1d2755a9 100644 --- a/src/coreComponents/mesh/generators/WellGeneratorBase.cpp +++ b/src/coreComponents/mesh/generators/WellGeneratorBase.cpp @@ -467,6 +467,24 @@ void WellGeneratorBase::checkPerforationLocationsValidity() mergePerforations( elemToPerfMap ); } + // check that the top well element (wellhead) does not have a perforation + // The top element carries the well control equation (BHP or rate constraint). + // Having a perforation on the same element creates a strong coupling between the control constraint + // and the perforation flux, which degrades Newton convergence (both isothermal and thermal). + // The top segment should be reserved for the well boundary condition constraint only. + for( globalIndex iwelem = 0; iwelem < m_numElems; ++iwelem ) + { + GEOS_THROW_IF( m_nextElemId[iwelem] < 0 && elemToPerfMap[iwelem].size() > 0, + "Well '" << getName() << "': a perforation is placed on the top well element (wellhead). " + << "This is not allowed because the top element is reserved for the well control constraint " + << "(BHP or rate) and must not have any perforation. \n\n" + << "To fix this, extend the well polyline upward by adding a segment above the first perforation. " + << "For example, add a node above the current well head in \"polylineNodeCoords\" and a corresponding " + << "entry in \"polylineSegmentConn\". " + << "This top segment will act as a constraint-only segment with no perforation.", + InputError ); + } + for( globalIndex iwelem = 0; iwelem < m_numElems; ++iwelem ) { // check that there is always a perforation in the last well element (otherwise, the problem is not well posed) @@ -516,7 +534,7 @@ void WellGeneratorBase::mergePerforations( array1d< array1d< localIndex > > cons continue; } - GEOS_LOG_RANK_0( "\n \nThe GEOSX wells currently have the following limitation in parallel: \n" + GEOS_LOG_RANK_0( "\n \nThe GEOS wells currently have the following limitation in parallel: \n" << "We cannot allow an element of the well mesh to have two or more perforations associated with it. \n" << "So, in the present simulation, perforation #" << elemToPerfMap[iwelem][ip] << " of well " << getName() diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CMakeLists.txt b/src/coreComponents/physicsSolvers/fluidFlow/CMakeLists.txt index 5df1a8ce699..7f75939d548 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CMakeLists.txt +++ b/src/coreComponents/physicsSolvers/fluidFlow/CMakeLists.txt @@ -67,6 +67,7 @@ set( fluidFlowSolvers_headers kernels/singlePhase/ThermalAccumulationKernels.hpp kernels/singlePhase/ThermalDirichletFluxComputeKernel.hpp kernels/singlePhase/ThermalFluxComputeKernel.hpp + kernels/singlePhase/ThermalSolutionScalingKernel.hpp kernels/singlePhase/proppant/ProppantBaseKernels.hpp kernels/singlePhase/proppant/ProppantFluxKernels.hpp kernels/singlePhase/reactive/AccumulationKernels.hpp diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp index 3a948b6fee9..5e921ce2165 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp @@ -56,6 +56,7 @@ template< typename POROUSWRAPPER_TYPE > void updatePorosityAndPermeabilityFromPressureAndTemperature( POROUSWRAPPER_TYPE porousWrapper, CellElementSubRegion & subRegion, arrayView1d< real64 const > const & pressure, + arrayView1d< real64 const > const & pressure_n, arrayView1d< real64 const > const & temperature ) { forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_DEVICE ( localIndex const k ) @@ -64,6 +65,7 @@ void updatePorosityAndPermeabilityFromPressureAndTemperature( POROUSWRAPPER_TYPE { porousWrapper.updateStateFromPressureAndTemperature( k, q, pressure[k], + pressure_n[k], temperature[k] ); } } ); @@ -168,6 +170,12 @@ FlowSolverBase::FlowSolverBase( string const & name, setApplyDefaultValue( -1.0 ). // disabled by default setDescription( "Maximum (absolute) pressure change in a Newton iteration" ); + this->registerWrapper( viewKeyStruct::maxAbsoluteTempChangeString(), &m_maxAbsoluteTempChange ). + setSizedFromParent( 0 ). + setInputFlag( InputFlags::OPTIONAL ). + setApplyDefaultValue( -1.0 ). // disabled by default + setDescription( "Maximum (absolute) temperature change in a Newton iteration" ); + this->registerWrapper( viewKeyStruct::maxSequentialPresChangeString(), &m_maxSequentialPresChange ). setSizedFromParent( 0 ). setInputFlag( InputFlags::OPTIONAL ). @@ -667,7 +675,8 @@ void FlowSolverBase::updatePorosityAndPermeability( CellElementSubRegion & subRe } else { - updatePorosityAndPermeabilityFromPressureAndTemperature( porousWrapper, subRegion, pressure, temperature ); + arrayView1d< real64 const > const & pressure_n = subRegion.getField< fields::flow::pressure_n >(); + updatePorosityAndPermeabilityFromPressureAndTemperature( porousWrapper, subRegion, pressure, pressure_n, temperature ); } } ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp index 75afb4fcd65..3c7170ea0fc 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp @@ -77,6 +77,7 @@ class FlowSolverBase : public PhysicsSolverBase static constexpr char const * inputTemperatureString() { return "temperature"; } static constexpr char const * allowNegativePressureString() { return "allowNegativePressure"; } static constexpr char const * maxAbsolutePresChangeString() { return "maxAbsolutePressureChange"; } + static constexpr char const * maxAbsoluteTempChangeString() { return "maxAbsoluteTemperatureChange"; } static constexpr char const * maxSequentialPresChangeString() { return "maxSequentialPressureChange"; } static constexpr char const * maxSequentialTempChangeString() { return "maxSequentialTemperatureChange"; } @@ -301,6 +302,9 @@ class FlowSolverBase : public PhysicsSolverBase /// maximum (absolute) pressure change in a Newton iteration real64 m_maxAbsolutePresChange; + /// maximum (absolute) temperature change in a Newton iteration + real64 m_maxAbsoluteTempChange; + /// maximum (absolute) pressure change in a sequential iteration real64 m_sequentialPresChange; real64 m_maxSequentialPresChange; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp index 562cfae6df4..2f6adac6a37 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp @@ -44,6 +44,7 @@ #include "physicsSolvers/fluidFlow/kernels/singlePhase/MobilityKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/singlePhase/SolutionCheckKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/singlePhase/SolutionScalingKernel.hpp" +#include "physicsSolvers/fluidFlow/kernels/singlePhase/ThermalSolutionScalingKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/singlePhase/StatisticsKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/singlePhase/HydrostaticPressureKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/singlePhase/ThermalHydrostaticPressureKernel.hpp" @@ -103,9 +104,12 @@ void SinglePhaseBase::registerDataOnMesh( Group & meshBodies ) subRegion.registerField< flow::mass_n >( getName() ); subRegion.registerField< flow::dMass >( getName() ).reference().resizeDimension< 1 >( m_numDofPerCell ); + subRegion.registerField< flow::pressureScalingFactor >( getName() ); + if( m_isThermal ) { subRegion.registerField< flow::dEnergy >( getName() ).reference().resizeDimension< 1 >( m_numDofPerCell ); + subRegion.registerField< flow::temperatureScalingFactor >( getName() ); } } ); @@ -1297,34 +1301,89 @@ real64 SinglePhaseBase::scalingForSystemSolution( DomainPartition & domain, string const dofKey = dofManager.getKey( viewKeyStruct::elemDofFieldString() ); real64 scalingFactor = 1.0; real64 maxDeltaPres = 0.0; + real64 maxDeltaTemp = 0.0; - forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, - MeshLevel & mesh, - string_array const & regionNames ) + real64 minPresScalingFactor = 1.0, minTempScalingFactor = 1.0; + + if( m_isThermal ) { - mesh.getElemManager().forElementSubRegions( regionNames, - [&]( localIndex const, - ElementSubRegionBase & subRegion ) + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, + MeshLevel & mesh, + string_array const & regionNames ) { - globalIndex const rankOffset = dofManager.rankOffset(); - arrayView1d< globalIndex const > const dofNumber = subRegion.getReference< array1d< globalIndex > >( dofKey ); - arrayView1d< integer const > const ghostRank = subRegion.ghostRank(); + mesh.getElemManager().forElementSubRegions( regionNames, + [&]( localIndex const, + ElementSubRegionBase & subRegion ) + { + globalIndex const rankOffset = dofManager.rankOffset(); + arrayView1d< globalIndex const > const dofNumber = subRegion.getReference< array1d< globalIndex > >( dofKey ); + arrayView1d< integer const > const ghostRank = subRegion.ghostRank(); + + arrayView1d< real64 > pressureScalingFactor = subRegion.getField< flow::pressureScalingFactor >(); + arrayView1d< real64 > temperatureScalingFactor = subRegion.getField< flow::temperatureScalingFactor >(); - auto const subRegionData = - singlePhaseBaseKernels::SolutionScalingKernel:: - launch< parallelDevicePolicy<> >( localSolution, rankOffset, dofNumber, ghostRank, m_maxAbsolutePresChange ); + auto const subRegionData = thermalSinglePhaseBaseKernels:: + SolutionScalingKernel:: + launch< parallelDevicePolicy<> >( localSolution, rankOffset, 1, dofNumber, ghostRank, + m_maxAbsolutePresChange, m_maxAbsoluteTempChange, + pressureScalingFactor, temperatureScalingFactor ); - scalingFactor = std::min( scalingFactor, subRegionData.first ); - maxDeltaPres = std::max( maxDeltaPres, subRegionData.second ); + scalingFactor = std::min( scalingFactor, std::get< 0 >( subRegionData ) ); + maxDeltaPres = std::max( maxDeltaPres, std::get< 1 >( subRegionData ) ); + maxDeltaTemp = std::max( maxDeltaTemp, std::get< 2 >( subRegionData ) ); + + minPresScalingFactor = std::min( minPresScalingFactor, std::get< 3 >( subRegionData ) ); + minTempScalingFactor = std::min( minTempScalingFactor, std::get< 4 >( subRegionData ) ); + } ); } ); - } ); + } + else + { + GEOS_UNUSED_VAR( maxDeltaTemp ); + + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, + MeshLevel & mesh, + string_array const & regionNames ) + { + mesh.getElemManager().forElementSubRegions( regionNames, + [&]( localIndex const, + ElementSubRegionBase & subRegion ) + { + globalIndex const rankOffset = dofManager.rankOffset(); + arrayView1d< globalIndex const > const dofNumber = subRegion.getReference< array1d< globalIndex > >( dofKey ); + arrayView1d< integer const > const ghostRank = subRegion.ghostRank(); + + auto const subRegionData = singlePhaseBaseKernels:: + SolutionScalingKernel:: + launch< parallelDevicePolicy<> >( localSolution, rankOffset, dofNumber, ghostRank, m_maxAbsolutePresChange ); + + scalingFactor = std::min( scalingFactor, subRegionData.first ); + minPresScalingFactor = std::min( minPresScalingFactor, subRegionData.first ); + maxDeltaPres = std::max( maxDeltaPres, subRegionData.second ); + } ); + } ); + } scalingFactor = MpiWrapper::min( scalingFactor ); + minPresScalingFactor = MpiWrapper::min( minPresScalingFactor ); maxDeltaPres = MpiWrapper::max( maxDeltaPres ); GEOS_LOG_LEVEL_RANK_0( logInfo::Solution, GEOS_FMT( " {}: Max pressure change = {} Pa (before scaling)", getName(), fmt::format( "{:.{}f}", maxDeltaPres, 3 ) ) ); + GEOS_LOG_LEVEL_RANK_0( logInfo::Solution, GEOS_FMT( " {}: Min pressure scaling factor = {}", getName(), minPresScalingFactor ) ); + + if( m_isThermal ) + { + minTempScalingFactor = MpiWrapper::min( minTempScalingFactor ); + maxDeltaTemp = MpiWrapper::max( maxDeltaTemp ); + + GEOS_LOG_LEVEL_RANK_0( logInfo::Solution, GEOS_FMT( " {}: Max temperature change = {} K (before scaling)", + getName(), fmt::format( "{:.{}f}", maxDeltaTemp, 3 ) ) ); + + GEOS_LOG_LEVEL_RANK_0( logInfo::Solution, GEOS_FMT( " {}: Min temperature scaling factor = {}", getName(), minTempScalingFactor ) ); + } + return scalingFactor; } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp index b507ea1b22c..50be0c3a3fc 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp @@ -278,13 +278,13 @@ void SinglePhaseFVM< BASE >::applySystemSolution( DofManager const & dofManager, dofManager.addVectorToField( localSolution, BASE::viewKeyStruct::elemDofFieldString(), flow::pressure::key(), - scalingFactor, + flow::pressureScalingFactor::key(), pressureMask ); dofManager.addVectorToField( localSolution, BASE::viewKeyStruct::elemDofFieldString(), flow::temperature::key(), - scalingFactor, + flow::temperatureScalingFactor::key(), temperatureMask ); } else diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp index b5c850662d3..9c2c6610934 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp @@ -47,18 +47,49 @@ using namespace dataRepository; using namespace constitutive; template< typename POROUSWRAPPER_TYPE > -void updatePorosityAndPermeabilityFromPressureAndReactions( POROUSWRAPPER_TYPE porousWrapper, - ElementSubRegionBase & subRegion, - arrayView1d< real64 const > const & pressure, - arrayView2d< real64 const, compflow::USD_COMP > const & kineticReactionMolarIncrements ) +void updatePorosityAndPermeabilityFromPressureTemperatureAndReactions( POROUSWRAPPER_TYPE porousWrapper, + ElementSubRegionBase & subRegion, + arrayView1d< real64 const > const & pressure, + arrayView1d< real64 const > const & pressure_n, + arrayView1d< real64 const > const & temperature, + arrayView2d< real64 const, compflow::USD_COMP > const & kineticReactionMolarIncrements ) { forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_DEVICE ( localIndex const k ) { for( localIndex q = 0; q < porousWrapper.numGauss(); ++q ) { - porousWrapper.updateStateFromPressureAndReactions( k, q, - pressure[k], - kineticReactionMolarIncrements[k] ); + porousWrapper.updateStateFromPressureTemperatureAndReactions( k, q, + pressure[k], + pressure_n[k], + temperature[k], + kineticReactionMolarIncrements[k] ); + } + } ); +} + +template< typename POROUSWRAPPER_TYPE > +void updatePorosityAndPermeabilityReactionsFixedStress( POROUSWRAPPER_TYPE porousWrapper, + ElementSubRegionBase & subRegion, + arrayView1d< real64 const > const & pressure, + arrayView1d< real64 const > const & pressure_k, + arrayView1d< real64 const > const & pressure_n, + arrayView1d< real64 const > const & temperature, + arrayView1d< real64 const > const & temperature_k, + arrayView1d< real64 const > const & temperature_n, + arrayView2d< real64 const, compflow::USD_COMP > const & kineticReactionMolarIncrements ) +{ + forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_DEVICE ( localIndex const k ) + { + for( localIndex q = 0; q < porousWrapper.numGauss(); ++q ) + { + porousWrapper.updateStateReactionsFixedStress( k, q, + pressure[k], + pressure_k[k], + pressure_n[k], + temperature[k], + temperature_k[k], + temperature_n[k], + kineticReactionMolarIncrements[k] ); } } ); } @@ -234,14 +265,14 @@ void SinglePhaseReactiveTransport::validateConstitutiveModels( DomainPartition & PorosityBase const & porosity = getConstitutiveModel< PorosityBase >( subRegion, porosityModelName ); - GEOS_THROW_IF( m_isUpdateReactivePorosity && (porosity.getCatalogName() != "ReactivePorosity"), + GEOS_THROW_IF( m_isUpdateReactivePorosity && (porosity.getCatalogName() != "ReactivePorosity" && porosity.getCatalogName() != "BiotReactivePorosity"), GEOS_FMT( "SinglePhaseReactiveTransport {}: the reaction porosity update option is enabled in the solver, but the porosity model {} is not for reactive porosity", getDataContext(), porosity.getDataContext() ), InputError ); if( m_isUpdateReactivePorosity ) { - ReactivePorosity const & reactivePorosity = getConstitutiveModel< ReactivePorosity >( subRegion, porosityModelName ); + ReactivePorosityBase const & reactivePorosity = getConstitutiveModel< ReactivePorosityBase >( subRegion, porosityModelName ); GEOS_THROW_IF_NE_MSG( reactivePorosity.numKineticReactions(), m_numKineticReactions, GEOS_FMT( "Mismatch in number of kinetic reactions, check the number of components input in porosity model {}", @@ -662,15 +693,27 @@ void SinglePhaseReactiveTransport::updatePorosityAndPermeability( CellElementSub if( m_isUpdateReactivePorosity ) { arrayView1d< real64 const > const & pressure = subRegion.getField< fields::flow::pressure >(); + arrayView1d< real64 const > const & pressure_n = subRegion.getField< fields::flow::pressure_n >(); + arrayView1d< real64 const > const & temperature = subRegion.getField< fields::flow::temperature >(); arrayView2d< real64 const, compflow::USD_COMP > const kineticReactionMolarIncrements = subRegion.getField< fields::flow::kineticReactionMolarIncrements >(); string const & solidName = subRegion.getReference< string >( viewKeyStruct::solidNamesString() ); CoupledSolidBase & porousSolid = subRegion.template getConstitutiveModel< CoupledSolidBase >( solidName ); - constitutive::ConstitutivePassThru< ReactiveSolidBase >::execute( porousSolid, [=, &subRegion] ( auto & castedPorousSolid ) + constitutive::ConstitutivePassThru< CoupledSolidBase >::execute( porousSolid, [=, &subRegion] ( auto & castedPorousSolid ) { typename TYPEOFREF( castedPorousSolid ) ::KernelWrapper porousWrapper = castedPorousSolid.createKernelUpdates(); - updatePorosityAndPermeabilityFromPressureAndReactions( porousWrapper, subRegion, pressure, kineticReactionMolarIncrements ); + if( m_isFixedStressPoromechanicsUpdate ) + { + arrayView1d< real64 const > const & pressure_k = subRegion.getField< fields::flow::pressure_k >(); + arrayView1d< real64 const > const & temperature_n = subRegion.getField< fields::flow::temperature_n >(); + arrayView1d< real64 const > const & temperature_k = subRegion.getField< fields::flow::temperature_k >(); + updatePorosityAndPermeabilityReactionsFixedStress( porousWrapper, subRegion, pressure, pressure_k, pressure_n, temperature, temperature_k, temperature_n, kineticReactionMolarIncrements ); + } + else + { + updatePorosityAndPermeabilityFromPressureTemperatureAndReactions( porousWrapper, subRegion, pressure, pressure_n, temperature, kineticReactionMolarIncrements ); + } } ); } else @@ -679,6 +722,14 @@ void SinglePhaseReactiveTransport::updatePorosityAndPermeability( CellElementSub } } +// To modify for chemical coupling later +void SinglePhaseReactiveTransport::updatePorosityAndPermeability( SurfaceElementSubRegion & subRegion ) const +{ + GEOS_MARK_FUNCTION; + + FlowSolverBase::updatePorosityAndPermeability( subRegion ); +} + void SinglePhaseReactiveTransport::updateMixedReactionSystem( ElementSubRegionBase & subRegion ) const { GEOS_MARK_FUNCTION; @@ -724,7 +775,7 @@ void SinglePhaseReactiveTransport::updateSurfaceArea( ElementSubRegionBase & sub string const & solidName = subRegion.getReference< string >( viewKeyStruct::solidNamesString() ); CoupledSolidBase & porousSolid = subRegion.template getConstitutiveModel< CoupledSolidBase >( solidName ); - constitutive::ConstitutivePassThru< ReactiveSolidBase >::execute( porousSolid, [=, &subRegion] ( auto & castedPorousSolid ) + constitutive::ConstitutivePassThru< CoupledSolidBase >::execute( porousSolid, [=, &subRegion] ( auto & castedPorousSolid ) { typename TYPEOFREF( castedPorousSolid ) ::KernelWrapper porousWrapper = castedPorousSolid.createKernelUpdates(); updateSurfaceAreaFromReactions( porousWrapper, subRegion, initialSurfaceArea, surfaceArea ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.hpp index 544b0f612f3..352ce1f68ed 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.hpp @@ -159,6 +159,7 @@ class SinglePhaseReactiveTransport : public SinglePhaseBase virtual void updateFluidModel( ObjectManagerBase & dataGroup ) const override; virtual void updatePorosityAndPermeability( CellElementSubRegion & subRegion ) const override; + virtual void updatePorosityAndPermeability( SurfaceElementSubRegion & subRegion ) const override; virtual void initializePostInitialConditionsPreSubGroups() override; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalSolutionScalingKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalSolutionScalingKernel.hpp new file mode 100644 index 00000000000..5ffc3ff7ed7 --- /dev/null +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/ThermalSolutionScalingKernel.hpp @@ -0,0 +1,103 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file ThermalSolutionScalingKernel.hpp + */ + +#ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASE_THERMALSOLUTIONSCALINGKERNEL_HPP +#define GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASE_THERMALSOLUTIONSCALINGKERNEL_HPP + +#include "common/DataTypes.hpp" +#include "common/GEOS_RAJA_Interface.hpp" + +namespace geos +{ + +namespace thermalSinglePhaseBaseKernels +{ + +/******************************** SolutionScalingKernel ********************************/ + +struct SolutionScalingKernel +{ + template< typename POLICY > + static std::tuple< real64, real64, real64, real64, real64 > launch( arrayView1d< real64 const > const & localSolution, + globalIndex const rankOffset, + globalIndex const temperatureOffset, + arrayView1d< globalIndex const > const & dofNumber, + arrayView1d< integer const > const & ghostRank, + real64 const maxAbsolutePresChange, + real64 const maxAbsoluteTempChange, + arrayView1d< real64 > pressureScalingFactor, + arrayView1d< real64 > temperatureScalingFactor ) + { + RAJA::ReduceMin< ReducePolicy< POLICY >, real64 > scalingFactor( 1.0 ); + RAJA::ReduceMax< ReducePolicy< POLICY >, real64 > maxDeltaPres( 0.0 ); + RAJA::ReduceMax< ReducePolicy< POLICY >, real64 > maxDeltaTemp( 0.0 ); + + RAJA::ReduceMin< ReducePolicy< POLICY >, real64 > localMinPresScalingFactor( 1.0 ); + RAJA::ReduceMin< ReducePolicy< POLICY >, real64 > localMinTempScalingFactor( 1.0 ); + + forAll< POLICY >( dofNumber.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei ) mutable + { + if( ghostRank[ei] < 0 && dofNumber[ei] >= 0 ) + { + pressureScalingFactor[ei] = 1.0; + temperatureScalingFactor[ei] = 1.0; + + localIndex const lid = dofNumber[ei] - rankOffset; + + // compute the change in pressure + real64 const absPresChange = LvArray::math::abs( localSolution[lid] ); + maxDeltaPres.max( absPresChange ); + + // compute the change in temperature + real64 const absTempChange = LvArray::math::abs( localSolution[lid + temperatureOffset] ); + maxDeltaTemp.max( absTempChange ); + + // maxAbsolutePresChange <= 0.0 means that scaling is disabled, and we are only collecting maxDeltaPres in this kernel + if( maxAbsolutePresChange > 0.0 && absPresChange > maxAbsolutePresChange ) + { + real64 const presScalingFactor = maxAbsolutePresChange / absPresChange; + pressureScalingFactor[ei] = presScalingFactor; + scalingFactor.min( presScalingFactor ); + + localMinPresScalingFactor.min( presScalingFactor ); + } + + // maxAbsoluteTempChange <= 0.0 means that scaling is disabled, and we are only collecting maxDeltaTemps in this kernel + if( maxAbsoluteTempChange > 0.0 && absTempChange > maxAbsoluteTempChange ) + { + real64 const tempScalingFactor = maxAbsoluteTempChange / absTempChange; + temperatureScalingFactor[ei] = tempScalingFactor; + scalingFactor.min( tempScalingFactor ); + + localMinTempScalingFactor.min( tempScalingFactor ); + } + } + + } ); + + return { scalingFactor.get(), maxDeltaPres.get(), maxDeltaTemp.get(), localMinPresScalingFactor.get(), localMinTempScalingFactor.get() }; + } + +}; + +} // namespace thermalSinglePhaseBaseKernels + +} // namespace geos + +#endif //GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASE_THERMALSOLUTIONSCALINGKERNEL_HPP diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/KernelLaunchSelectors.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/KernelLaunchSelectors.hpp index 7fa4171c9b3..f2298c82df7 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/KernelLaunchSelectors.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/KernelLaunchSelectors.hpp @@ -62,6 +62,8 @@ void kernelLaunchSelectorCompSwitch( T value, LAMBDA && lambda ) { lambda( std::integral_constant< T, 8 >() ); return; } case 9: { lambda( std::integral_constant< T, 9 >() ); return; } + case 10: + { lambda( std::integral_constant< T, 10 >() ); return; } default: { GEOS_ERROR( GEOS_FMT( "Unsupported number of primary species: {}", value ) ); } } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp index 6a72af7a90c..514a1f8869c 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp @@ -87,7 +87,6 @@ CompositionalMultiphaseWell::CompositionalMultiphaseWell( const string & name, m_useTotalMassEquation( 1 ), m_maxCompFracChange( 1.0 ), m_maxRelativePresChange( 0.2 ), - m_maxAbsolutePresChange( -1 ), // disabled by default m_minScalingFactor( 0.01 ), m_allowCompDensChopping( 1 ), m_targetPhaseIndex( -1 ) @@ -120,12 +119,6 @@ CompositionalMultiphaseWell::CompositionalMultiphaseWell( const string & name, setApplyDefaultValue( 1.0 ). setDescription( "Maximum (relative) change in pressure between two Newton iterations (recommended with rate control)" ); - this->registerWrapper( viewKeyStruct::maxAbsolutePresChangeString(), &m_maxAbsolutePresChange ). - setSizedFromParent( 0 ). - setInputFlag( InputFlags::OPTIONAL ). - setApplyDefaultValue( -1.0 ). // disabled by default - setDescription( "Maximum (absolute) pressure change in a Newton iteration" ); - this->registerWrapper( viewKeyStruct::maxRelativeTempChangeString(), &m_maxRelativeTempChange ). setSizedFromParent( 0 ). setInputFlag( InputFlags::OPTIONAL ). diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp index 677ea19d64c..335b0847617 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp @@ -270,8 +270,6 @@ class CompositionalMultiphaseWell : public WellSolverBase static constexpr char const * maxRelativePresChangeString() { return "maxRelativePressureChange"; } - static constexpr char const * maxAbsolutePresChangeString() { return "maxAbsolutePressureChange"; } - static constexpr char const * maxRelativeCompDensChangeString() { return "maxRelativeCompDensChange"; } static constexpr char const * maxRelativeTempChangeString() { return "maxRelativeTemperatureChange"; } @@ -398,9 +396,6 @@ class CompositionalMultiphaseWell : public WellSolverBase /// maximum (relative) change in pressure between two Newton iterations real64 m_maxRelativePresChange; - /// maximum (absolute) change in pressure between two Newton iterations - real64 m_maxAbsolutePresChange; - /// maximum (relative) change in component density between two Newton iterations real64 m_maxRelativeCompDensChange; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp index b2bbfabfe9b..fc7a297ecd0 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp @@ -45,6 +45,8 @@ #include "physicsSolvers/fluidFlow/wells/kernels/SinglePhasePerforationFluxKernels.hpp" #include "physicsSolvers/fluidFlow/kernels/singlePhase/FluidUpdateKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/singlePhase/SolutionCheckKernel.hpp" +#include "physicsSolvers/fluidFlow/kernels/singlePhase/SolutionScalingKernel.hpp" +#include "physicsSolvers/fluidFlow/kernels/singlePhase/ThermalSolutionScalingKernel.hpp" #include "physicsSolvers/fluidFlow/SinglePhaseStatisticsAggregator.hpp" namespace geos @@ -102,12 +104,16 @@ void SinglePhaseWell::registerDataOnMesh( Group & meshBodies ) subRegion.registerField< well::connectionRate_n >( getName() ); subRegion.registerField< well::connectionRate >( getName() ); + subRegion.registerField< well::pressureScalingFactor >( getName() ); + PerforationData & perforationData = *subRegion.getPerforationData(); perforationData.registerField< well::perforationRate >( getName() ); perforationData.registerField< well::dPerforationRate >( getName() ). reference().resizeDimension< 1, 2 >( 2, 2 ); if( isThermal() ) { + subRegion.registerField< well::temperatureScalingFactor >( getName() ); + perforationData.registerField< well::energyPerforationFlux >( getName() ); perforationData.registerField< well::dEnergyPerforationFlux >( getName() ). reference().resizeDimension< 1, 2 >( 2, 2 ); @@ -319,11 +325,28 @@ void SinglePhaseWell::updateVolRateForConstraint( WellElementSubRegion & subRegi // - Surface conditions: using the surface pressure provided by the user // - Reservoir conditions: using the pressure in the top element - fluidWrapper.update( iwelemRef, 0, refConditions.pressure ); + // Update fluid properties with reference conditions + if constexpr ( IS_THERMAL ) + { + fluidWrapper.update( iwelemRef, 0, refConditions.pressure, refConditions.temperature ); + } + else + { + fluidWrapper.update( iwelemRef, 0, refConditions.pressure ); + } + if( useSurfaceConditions && logSurfaceCondition ) { - GEOS_LOG_RANK( GEOS_FMT( "{}: surface density computed with P_surface = {} Pa", - wellControlsName, refConditions.pressure ) ); + if constexpr ( IS_THERMAL ) + { + GEOS_LOG_RANK( GEOS_FMT( "{}: surface density computed with P_surface = {} Pa, T_surface = {} K", + wellControlsName, refConditions.pressure, refConditions.temperature ) ); + } + else + { + GEOS_LOG_RANK( GEOS_FMT( "{}: surface density computed with P_surface = {} Pa", + wellControlsName, refConditions.pressure ) ); + } } #ifdef GEOS_USE_HIP @@ -1059,6 +1082,102 @@ SinglePhaseWell::calculateResidualNorm( real64 const & time_n, return resNorm; } +real64 +SinglePhaseWell::scalingForSystemSolution( DomainPartition & domain, + DofManager const & dofManager, + arrayView1d< real64 const > const & localSolution ) +{ + GEOS_MARK_FUNCTION; + + string const wellDofKey = dofManager.getKey( wellElementDofName() ); + + real64 scalingFactor = 1.0; + real64 maxDeltaPres = 0.0, maxDeltaTemp = 0.0; + + real64 minPresScalingFactor = 1.0, minTempScalingFactor = 1.0; + + if( m_isThermal ) + { + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, + MeshLevel & mesh, + string_array const & regionNames ) + { + mesh.getElemManager().forElementSubRegions( regionNames, + [&]( localIndex const, + ElementSubRegionBase & subRegion ) + { + globalIndex const rankOffset = dofManager.rankOffset(); + arrayView1d< globalIndex const > const dofNumber = subRegion.getReference< array1d< globalIndex > >( wellDofKey ); + arrayView1d< integer const > const ghostRank = subRegion.ghostRank(); + + arrayView1d< real64 > pressureScalingFactor = subRegion.getField< well::pressureScalingFactor >(); + arrayView1d< real64 > temperatureScalingFactor = subRegion.getField< well::temperatureScalingFactor >(); + + auto const subRegionData = thermalSinglePhaseBaseKernels:: + SolutionScalingKernel:: + launch< parallelDevicePolicy<> >( localSolution, rankOffset, 2, dofNumber, ghostRank, + m_maxAbsolutePresChange, m_maxAbsoluteTempChange, + pressureScalingFactor, temperatureScalingFactor ); + + scalingFactor = std::min( scalingFactor, std::get< 0 >( subRegionData ) ); + maxDeltaPres = std::max( maxDeltaPres, std::get< 1 >( subRegionData ) ); + maxDeltaTemp = std::max( maxDeltaTemp, std::get< 2 >( subRegionData ) ); + + minPresScalingFactor = std::min( minPresScalingFactor, std::get< 3 >( subRegionData ) ); + minTempScalingFactor = std::min( minTempScalingFactor, std::get< 4 >( subRegionData ) ); + } ); + } ); + } + else + { + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, + MeshLevel & mesh, + string_array const & regionNames ) + { + mesh.getElemManager().forElementSubRegions( regionNames, + [&]( localIndex const, + ElementSubRegionBase & subRegion ) + { + globalIndex const rankOffset = dofManager.rankOffset(); + arrayView1d< globalIndex const > const dofNumber = subRegion.getReference< array1d< globalIndex > >( wellDofKey ); + arrayView1d< integer const > const ghostRank = subRegion.ghostRank(); + + auto const subRegionData = singlePhaseBaseKernels:: + SolutionScalingKernel:: + launch< parallelDevicePolicy<> >( localSolution, rankOffset, dofNumber, ghostRank, m_maxAbsolutePresChange ); + + scalingFactor = std::min( subRegionData.first, scalingFactor ); + minPresScalingFactor = std::min( minPresScalingFactor, subRegionData.first ); + maxDeltaPres = std::max( maxDeltaPres, subRegionData.second ); + } ); + } ); + } + + scalingFactor = MpiWrapper::min( scalingFactor ); + minPresScalingFactor = MpiWrapper::min( minPresScalingFactor ); + maxDeltaPres = MpiWrapper::max( maxDeltaPres ); + + GEOS_LOG_LEVEL_RANK_0( logInfo::Solution, + GEOS_FMT( " {}: Max well pressure change: {} Pa (before scaling)", + getName(), GEOS_FMT( "{:.{}f}", maxDeltaPres, 3 ) ) ); + + GEOS_LOG_LEVEL_RANK_0( logInfo::Solution, GEOS_FMT( " {}: Min pressure scaling factor = {}", getName(), minPresScalingFactor ) ); + + if( m_isThermal ) + { + minTempScalingFactor = MpiWrapper::min( minTempScalingFactor ); + maxDeltaTemp = MpiWrapper::max( maxDeltaTemp ); + GEOS_LOG_LEVEL_RANK_0( logInfo::Solution, + GEOS_FMT( " {}: Max well temperature change: {} K (before scaling)", + getName(), GEOS_FMT( "{:.{}f}", maxDeltaTemp, 3 ) ) ); + + GEOS_LOG_LEVEL_RANK_0( logInfo::Solution, GEOS_FMT( " {}: Min temperature scaling factor = {}", getName(), minTempScalingFactor ) ); + } + + return scalingFactor; + +} + bool SinglePhaseWell::checkSystemSolution( DomainPartition & domain, DofManager const & dofManager, arrayView1d< real64 const > const & localSolution, @@ -1122,19 +1241,19 @@ SinglePhaseWell::applySystemSolution( DofManager const & dofManager, real64 const dt, DomainPartition & domain ) { - GEOS_UNUSED_VAR( dt ); + GEOS_UNUSED_VAR( dt, scalingFactor ); DofManager::CompMask pressureMask( m_numDofPerWellElement, 0, 1 ); DofManager::CompMask connRateMask( m_numDofPerWellElement, 1, 2 ); dofManager.addVectorToField( localSolution, wellElementDofName(), well::pressure::key(), - scalingFactor, + well::pressureScalingFactor::key(), pressureMask ); dofManager.addVectorToField( localSolution, wellElementDofName(), well::connectionRate::key(), - scalingFactor, + well::pressureScalingFactor::key(), connRateMask ); if( isThermal() ) @@ -1144,7 +1263,7 @@ SinglePhaseWell::applySystemSolution( DofManager const & dofManager, dofManager.addVectorToField( localSolution, wellElementDofName(), fields::well::temperature::key(), - scalingFactor, + well::temperatureScalingFactor::key(), temperatureMask ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.hpp index b82ebe96342..c6aa2b6e721 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.hpp @@ -108,6 +108,11 @@ class SinglePhaseWell : public WellSolverBase DofManager const & dofManager, arrayView1d< real64 const > const & localRhs ) override; + virtual real64 + scalingForSystemSolution( DomainPartition & domain, + DofManager const & dofManager, + arrayView1d< real64 const > const & localSolution ) override; + virtual bool checkSystemSolution( DomainPartition & domain, DofManager const & dofManager, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.cpp index 61b9398d0e6..0d85d268f9e 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.cpp @@ -65,6 +65,18 @@ WellSolverBase::WellSolverBase( string const & name, setInputFlag( dataRepository::InputFlags::OPTIONAL ). setDescription( "Choose time step to honor rates/bhp tables time intervals" ); + this->registerWrapper( viewKeyStruct::maxAbsolutePresChangeString(), &m_maxAbsolutePresChange ). + setSizedFromParent( 0 ). + setInputFlag( InputFlags::OPTIONAL ). + setApplyDefaultValue( -1.0 ). // disabled by default + setDescription( "Maximum (absolute) pressure change in a Newton iteration" ); + + this->registerWrapper( viewKeyStruct::maxAbsoluteTempChangeString(), &m_maxAbsoluteTempChange ). + setSizedFromParent( 0 ). + setInputFlag( InputFlags::OPTIONAL ). + setApplyDefaultValue( -1.0 ). // disabled by default + setDescription( "Maximum (absolute) temperature change in a Newton iteration" ); + addLogLevel< logInfo::WellControl >(); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.hpp index 0f9bd2615b6..90bd499d39e 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.hpp @@ -289,6 +289,9 @@ class WellSolverBase : public PhysicsSolverBase static constexpr char const * timeStepFromTablesFlagString() { return "timeStepFromTables"; } static constexpr char const * fluidNamesString() { return "fluidNames"; } + + static constexpr char const * maxAbsolutePresChangeString() { return "maxAbsolutePressureChange"; } + static constexpr char const * maxAbsoluteTempChangeString() { return "maxAbsoluteTemperatureChange"; } }; private: @@ -357,6 +360,12 @@ class WellSolverBase : public PhysicsSolverBase /// name of the fluid constitutive model used as a reference for component/phase description string m_referenceFluidModelName; + + /// maximum (absolute) change in pressure between two Newton iterations + real64 m_maxAbsolutePresChange; + + /// maximum (absolute) change in temperature between two Newton iterations + real64 m_maxAbsoluteTempChange; }; } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/ThermalSinglePhaseWellKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/ThermalSinglePhaseWellKernels.hpp index f736826561c..2349c6d9042 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/ThermalSinglePhaseWellKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/kernels/ThermalSinglePhaseWellKernels.hpp @@ -122,7 +122,8 @@ class ElementBasedAssemblyKernel : public singlePhaseWellKernels::ElementBasedAs real64 localEnergyAccum; real64 localEnergyAccumDP; real64 localEnergyAccumDT; - if( iwelem == m_iwelemControl && !m_isProducer ) + // if( iwelem == m_iwelemControl && !m_isProducer ) + if( 1 ) { // For top segment energy balance eqn replaced with T(n+1) - T = 0 // No other energy balance derivatives @@ -327,7 +328,8 @@ class FaceBasedAssemblyKernel : public singlePhaseWellKernels::FaceBasedAssembly localEnergyFluxJacobian_dQ[0] = -m_dt *m_enthalpy[iwelem][0]; localEnergyFluxJacobian[FLUID_PROP_COFFSET::dP] = -m_dt * currentConnRate * m_dEnthalpy[iwelem][0][FLUID_PROP_COFFSET::dP]; localEnergyFluxJacobian[FLUID_PROP_COFFSET::dT] = -m_dt * currentConnRate * m_dEnthalpy[iwelem][0][FLUID_PROP_COFFSET::dT]; - if( !m_isProducer && m_globalWellElementIndex[iwelem] == 0 ) + // if( !m_isProducer && m_globalWellElementIndex[iwelem] == 0 ) + if( 1 ) { localEnergyFlux[0]= 0.0; localEnergyFluxJacobian_dQ[0] = 0.0; @@ -373,7 +375,8 @@ class FaceBasedAssemblyKernel : public singlePhaseWellKernels::FaceBasedAssembly real64 dprop_dp = m_dt * currentConnRate *m_dEnthalpy[iwelem][0][FLUID_PROP_COFFSET::dP]; real64 dprop_dt = m_dt * currentConnRate * m_dEnthalpy[iwelem][0][FLUID_PROP_COFFSET::dT]; - if( !m_isProducer && m_globalWellElementIndex[iwelemNext] == 0 ) + // if( !m_isProducer && m_globalWellElementIndex[iwelemNext] == 0 ) + if( 1 ) { localEnergyFlux[TAG::NEXT ] = 0.0; localEnergyFluxJacobian_dQ [TAG::NEXT ][0] = 0.0; @@ -387,10 +390,14 @@ class FaceBasedAssemblyKernel : public singlePhaseWellKernels::FaceBasedAssembly localEnergyFluxJacobian[TAG::NEXT ][FLUID_PROP_COFFSET::dP] = dprop_dp; localEnergyFluxJacobian[TAG::NEXT][FLUID_PROP_COFFSET::dT] = dprop_dt; } - localEnergyFlux[TAG::CURRENT ] = -eflux * currentConnRate; - localEnergyFluxJacobian_dQ [TAG::CURRENT][0] = -eflux_dq; - localEnergyFluxJacobian[TAG::CURRENT][FLUID_PROP_COFFSET::dP] = -dprop_dp; - localEnergyFluxJacobian[TAG::CURRENT][FLUID_PROP_COFFSET::dT] = -dprop_dt; + // localEnergyFlux[TAG::CURRENT ] = -eflux * currentConnRate; + // localEnergyFluxJacobian_dQ [TAG::CURRENT][0] = -eflux_dq; + // localEnergyFluxJacobian[TAG::CURRENT][FLUID_PROP_COFFSET::dP] = -dprop_dp; + // localEnergyFluxJacobian[TAG::CURRENT][FLUID_PROP_COFFSET::dT] = -dprop_dt; + localEnergyFlux[TAG::CURRENT ] = 0.0; + localEnergyFluxJacobian_dQ [TAG::CURRENT][0] = 0.0; + localEnergyFluxJacobian[TAG::CURRENT][FLUID_PROP_COFFSET::dP] = 0.0; + localEnergyFluxJacobian[TAG::CURRENT][FLUID_PROP_COFFSET::dT] = 0.0; // Note this updates diag and offdiag for( integer i = 0; i < 2; ++i ) diff --git a/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndSinglePhaseWellKernels.hpp b/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndSinglePhaseWellKernels.hpp index b16360f690a..c424373a849 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndSinglePhaseWellKernels.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndSinglePhaseWellKernels.hpp @@ -366,14 +366,15 @@ class ThermalSinglePhaseFluxKernel : public IsothermalSinglePhaseFluxKernel< IS_ stackArray1d< globalIndex, 2*resNumDOF > & dofColIndices, localIndex const iwelem ) { - // No energy equation if top element and Injector - // Top element defined by global index == 0 - // Assumption is global index == 0 is top segment with fixed temp BC - if( !m_isProducer ) - { - if( m_globalWellElementIndex[iwelem] == 0 ) - return; - } + // For injector top element (global index == 0), the well energy equation + // is replaced by a Dirichlet BC (T = T_inj), so we must not assemble + // the well-side energy flux. However, we still assemble the + // reservoir-side energy flux so the reservoir cell receives the correct + // enthalpy from the injected mass. + // bool const isTopInjectorElement = !m_isProducer && m_globalWellElementIndex[iwelem] == 0; + GEOS_UNUSED_VAR( iwelem ); + bool const isTopInjectorElement = 1; + // local working variables and arrays stackArray1d< localIndex, 2 > eqnRowIndices( 2 ); @@ -381,22 +382,22 @@ class ThermalSinglePhaseFluxKernel : public IsothermalSinglePhaseFluxKernel< IS_ stackArray2d< real64, 2*2 * resNumDOF > localPerfJacobian( 2, 2 * resNumDOF ); - // equantion offsets - note res and well have different equation lineups + // equation offsets - note res and well have different equation lineups eqnRowIndices[TAG::RES ] = LvArray::integerConversion< localIndex >( resOffset - m_rankOffset ) + 1; eqnRowIndices[TAG::WELL ] = LvArray::integerConversion< localIndex >( wellElemOffset - m_rankOffset ) + WJ_ROFFSET::ENERGYBAL; // populate local flux vector and derivatives localPerf[TAG::RES ] = m_dt * m_energyPerfFlux[iperf]; - localPerf[TAG::WELL ] = -m_dt * m_energyPerfFlux[iperf]; + localPerf[TAG::WELL ] = isTopInjectorElement ? 0.0 : -m_dt * m_energyPerfFlux[iperf]; for( integer ke = 0; ke < 2; ++ke ) { localIndex localDofIndexPres = ke * resNumDOF; localPerfJacobian[TAG::RES ][localDofIndexPres] = m_dt * m_dEnergyPerfFlux[iperf][ke][CP_Deriv::dP]; - localPerfJacobian[TAG::WELL ][localDofIndexPres] = -m_dt * m_dEnergyPerfFlux[iperf][ke][CP_Deriv::dP]; + localPerfJacobian[TAG::WELL ][localDofIndexPres] = isTopInjectorElement ? 0.0 : -m_dt * m_dEnergyPerfFlux[iperf][ke][CP_Deriv::dP]; localPerfJacobian[TAG::RES ][localDofIndexPres+1] = m_dt * m_dEnergyPerfFlux[iperf][ke][CP_Deriv::dT]; - localPerfJacobian[TAG::WELL][localDofIndexPres+1] = -m_dt * m_dEnergyPerfFlux[iperf][ke][CP_Deriv::dT]; + localPerfJacobian[TAG::WELL][localDofIndexPres+1] = isTopInjectorElement ? 0.0 : -m_dt * m_dEnergyPerfFlux[iperf][ke][CP_Deriv::dT]; } diff --git a/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellKernels.hpp b/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellKernels.hpp index 5d106df3156..fd260fd63e8 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellKernels.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellKernels.hpp @@ -432,14 +432,13 @@ class ThermalCompositionalMultiPhaseFluxKernel : public IsothermalCompositionalM stackArray1d< globalIndex, 2*resNumDOF > & dofColIndices, localIndex const iwelem ) { - // No energy equation if top element and Injector - // Top element defined by global index == 0 - // Assumption is global index == 0 is top segment with fixed temp BC - if( !m_isProducer ) - { - if( m_globalWellElementIndex[iwelem] == 0 ) - return; - } + // For injector top element (global index == 0), the well energy equation + // is replaced by a Dirichlet BC (T = T_inj), so we must not assemble + // the well-side energy flux. However, we still assemble the + // reservoir-side energy flux so the reservoir cell receives the correct + // enthalpy from the injected mass. + bool const isTopInjectorElement = !m_isProducer && m_globalWellElementIndex[iwelem] == 0; + // local working variables and arrays stackArray1d< localIndex, 2* numComp > eqnRowIndices( 2 ); @@ -453,23 +452,23 @@ class ThermalCompositionalMultiPhaseFluxKernel : public IsothermalCompositionalM // populate local flux vector and derivatives localPerf[TAG::RES ] = m_dt * m_energyPerfFlux[iperf]; - localPerf[TAG::WELL ] = -m_dt * m_energyPerfFlux[iperf]; + localPerf[TAG::WELL ] = isTopInjectorElement ? 0.0 : -m_dt * m_energyPerfFlux[iperf]; for( integer ke = 0; ke < 2; ++ke ) { localIndex localDofIndexPres = ke * resNumDOF; localPerfJacobian[TAG::RES ][localDofIndexPres] = m_dt * m_dEnergyPerfFlux[iperf][ke][CP_Deriv::dP]; - localPerfJacobian[TAG::WELL ][localDofIndexPres] = -m_dt * m_dEnergyPerfFlux[iperf][ke][CP_Deriv::dP]; + localPerfJacobian[TAG::WELL ][localDofIndexPres] = isTopInjectorElement ? 0.0 : -m_dt * m_dEnergyPerfFlux[iperf][ke][CP_Deriv::dP]; // populate local flux vector and derivatives for( integer ic = 0; ic < numComp; ++ic ) { localIndex const localDofIndexComp = localDofIndexPres + ic + 1; localPerfJacobian[TAG::RES ][localDofIndexComp] = m_dt * m_dEnergyPerfFlux[iperf][ke][CP_Deriv::dC+ic]; - localPerfJacobian[TAG::WELL][localDofIndexComp] = -m_dt * m_dEnergyPerfFlux[iperf][ke][CP_Deriv::dC+ic]; + localPerfJacobian[TAG::WELL][localDofIndexComp] = isTopInjectorElement ? 0.0 : -m_dt * m_dEnergyPerfFlux[iperf][ke][CP_Deriv::dC+ic]; } localPerfJacobian[TAG::RES ][localDofIndexPres+NC+1] = m_dt * m_dEnergyPerfFlux[iperf][ke][CP_Deriv::dT]; - localPerfJacobian[TAG::WELL][localDofIndexPres+NC+1] = -m_dt * m_dEnergyPerfFlux[iperf][ke][CP_Deriv::dT]; + localPerfJacobian[TAG::WELL][localDofIndexPres+NC+1] = isTopInjectorElement ? 0.0 : -m_dt * m_dEnergyPerfFlux[iperf][ke][CP_Deriv::dT]; } diff --git a/src/coreComponents/physicsSolvers/multiphysics/CoupledSolver.hpp b/src/coreComponents/physicsSolvers/multiphysics/CoupledSolver.hpp index c5cddf6acd1..256b817adca 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/CoupledSolver.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/CoupledSolver.hpp @@ -541,7 +541,7 @@ class CoupledSolver : public PhysicsSolverBase solver->saveSequentialIterationState( domain ); } - mapSolutionBetweenSolvers( domain, idx() ); + mapSolutionBetweenSolvers( stepDt, domain, idx() ); if( solverDt < stepDt ) // subsolver had to cut the time step { @@ -618,13 +618,15 @@ class CoupledSolver : public PhysicsSolverBase /** * @brief Maps the solution obtained from one solver to the fields used by the other solver(s) * + * @param dt timestep size * @param domain the domain partition * @param solverType the index of the solver withing this coupled solver. */ - virtual void mapSolutionBetweenSolvers( DomainPartition & domain, + virtual void mapSolutionBetweenSolvers( real64 const & dt, + DomainPartition & domain, integer const solverType ) { - GEOS_UNUSED_VAR( domain, solverType ); + GEOS_UNUSED_VAR( dt, domain, solverType ); } virtual bool checkSequentialConvergence( integer const cycleNumber, diff --git a/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldFractureSolver.cpp b/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldFractureSolver.cpp index 751dd7dea77..905eca73f24 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldFractureSolver.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldFractureSolver.cpp @@ -50,8 +50,9 @@ void PhaseFieldFractureSolver::postInputInitialization() getNonlinearSolverParameters().m_couplingType = NonlinearSolverParameters::CouplingType::Sequential; } -void PhaseFieldFractureSolver::mapSolutionBetweenSolvers( DomainPartition & domain, integer const solverType ) +void PhaseFieldFractureSolver::mapSolutionBetweenSolvers( real64 const & dt, DomainPartition & domain, integer const solverType ) { + GEOS_UNUSED_VAR( dt ); GEOS_MARK_FUNCTION; if( solverType == static_cast< integer >( SolverType::Damage ) ) diff --git a/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldFractureSolver.hpp b/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldFractureSolver.hpp index a1aae598694..af63526d7ca 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldFractureSolver.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldFractureSolver.hpp @@ -86,7 +86,7 @@ class PhaseFieldFractureSolver : public CoupledSolver< SolidMechanicsLagrangianF return std::get< toUnderlying( SolverType::Damage ) >( m_solvers ); } - virtual void mapSolutionBetweenSolvers( DomainPartition & Domain, integer const idx ) override final; + virtual void mapSolutionBetweenSolvers( real64 const & dt, DomainPartition & Domain, integer const idx ) override final; protected: diff --git a/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldPoromechanicsSolver.cpp b/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldPoromechanicsSolver.cpp index 06bb5259f06..a378cc65f83 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldPoromechanicsSolver.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldPoromechanicsSolver.cpp @@ -57,8 +57,10 @@ PhaseFieldPoromechanicsSolver::~PhaseFieldPoromechanicsSolver() // TODO Auto-generated destructor stub } -void PhaseFieldPoromechanicsSolver::mapSolutionBetweenSolvers( DomainPartition & domain, integer const solverType ) +void PhaseFieldPoromechanicsSolver::mapSolutionBetweenSolvers( real64 const & dt, DomainPartition & domain, integer const solverType ) { + GEOS_UNUSED_VAR( dt ); + if( solverType == static_cast< integer >( SolverType::Damage ) ) { GEOS_MARK_FUNCTION; diff --git a/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldPoromechanicsSolver.hpp b/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldPoromechanicsSolver.hpp index 75de5e22b49..33bcc15a738 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldPoromechanicsSolver.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldPoromechanicsSolver.hpp @@ -87,7 +87,7 @@ class PhaseFieldPoromechanicsSolver : public CoupledSolver< SinglePhasePoromecha return std::get< toUnderlying( SolverType::Damage ) >( m_solvers ); } - virtual void mapSolutionBetweenSolvers( DomainPartition & Domain, integer const idx ) override final; + virtual void mapSolutionBetweenSolvers( real64 const & dt, DomainPartition & Domain, integer const idx ) override final; void mapDamageAndGradientToQuadrature( DomainPartition & domain ); diff --git a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsConformingFractures.hpp b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsConformingFractures.hpp index 9e0828e19d5..1d886a18132 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsConformingFractures.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsConformingFractures.hpp @@ -564,7 +564,7 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) = 0; - virtual void mapSolutionBetweenSolvers( DomainPartition & domain, integer const solverType ) override + virtual void mapSolutionBetweenSolvers( real64 const & dt, DomainPartition & domain, integer const solverType ) override { GEOS_MARK_FUNCTION; @@ -581,7 +581,7 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER, this->flowSolver()->updateStencilWeights( domain ); } - Base::mapSolutionBetweenSolvers( domain, solverType ); + Base::mapSolutionBetweenSolvers( dt, domain, solverType ); } void updateHydraulicApertureAndFracturePermeability( DomainPartition & domain ) diff --git a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.cpp b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.cpp index 40cb0eaba82..96db0f63688 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.cpp @@ -162,6 +162,7 @@ typedef PoromechanicsInitialization< SinglePhasePoromechanicsConformingFractures typedef PoromechanicsInitialization< SinglePhasePoromechanicsConformingFracturesALM< SinglePhaseReservoirAndWells<> > > SinglePhaseReservoirPoromechanicsConformingFracturesALMInitialization; typedef PoromechanicsInitialization< SinglePhasePoromechanicsEmbeddedFractures > SinglePhasePoromechanicsEmbeddedFracturesInitialization; typedef PoromechanicsInitialization< SinglePhasePoromechanics< SinglePhaseReservoirAndWells<> > > SinglePhaseReservoirPoromechanicsInitialization; +typedef PoromechanicsInitialization< SinglePhasePoromechanics< SinglePhaseReactiveTransport > > SinglePhaseReactiveTransportPoromechanicsInitialization; typedef PoromechanicsInitialization< HydrofractureSolver< SinglePhasePoromechanics<> > > HydrofractureInitialization; REGISTER_CATALOG_ENTRY( TaskBase, MultiphasePoromechanicsInitialization, string const &, Group * const ) REGISTER_CATALOG_ENTRY( TaskBase, MultiphasePoromechanicsConformingFracturesInitialization, string const &, Group * const ) @@ -176,6 +177,7 @@ REGISTER_CATALOG_ENTRY( TaskBase, SinglePhasePoromechanicsConformingFracturesALM REGISTER_CATALOG_ENTRY( TaskBase, SinglePhaseReservoirPoromechanicsConformingFracturesALMInitialization, string const &, Group * const ) REGISTER_CATALOG_ENTRY( TaskBase, SinglePhasePoromechanicsEmbeddedFracturesInitialization, string const &, Group * const ) REGISTER_CATALOG_ENTRY( TaskBase, SinglePhaseReservoirPoromechanicsInitialization, string const &, Group * const ) +REGISTER_CATALOG_ENTRY( TaskBase, SinglePhaseReactiveTransportPoromechanicsInitialization, string const &, Group * const ) REGISTER_CATALOG_ENTRY( TaskBase, HydrofractureInitialization, string const &, Group * const ) } diff --git a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp index 1b34ab94f65..ff5a8c7bb3e 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp @@ -121,6 +121,12 @@ class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER EnumStrings< SolidMechanicsLagrangianFEM::TimeIntegrationOption >::toString( SolidMechanicsLagrangianFEM::TimeIntegrationOption::QuasiStatic ) ), InputError, this->solidMechanicsSolver()->getDataContext() ); + GEOS_THROW_IF( this->flowSolver()->getCatalogName() == "SinglePhaseReactiveTransport" && + this->getNonlinearSolverParameters().m_couplingType != NonlinearSolverParameters::CouplingType::Sequential, + GEOS_FMT( "{} {}: The coupling type must be Sequential since it is coupled with {}", + this->getCatalogName(), this->getName(), this->flowSolver()->getCatalogName() ), + InputError ); + // Sequential coupling uses the subsolver linear systems directly, so the // coupled solver does not need a top-level MGR strategy. if( this->getNonlinearSolverParameters().couplingType() != NonlinearSolverParameters::CouplingType::Sequential ) @@ -190,10 +196,20 @@ class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER if( this->getNonlinearSolverParameters().m_couplingType == NonlinearSolverParameters::CouplingType::Sequential ) { - // to let the solid mechanics solver that there is a pressure and temperature RHS in the mechanics solve - solidMechanicsSolver()->enableFixedStressPoromechanicsUpdate(); - // to let the flow solver that saving pressure_k and temperature_k is necessary (for the fixed-stress porosity terms) - flowSolver()->enableFixedStressPoromechanicsUpdate(); + if( flowSolver()->getCatalogName() == "SinglePhaseReactiveTransport" ) + { + // to let the solid mechanics solver to account for anelastic strain due to chemistry + solidMechanicsSolver()->enableExplicitChemomechanicsUpdate(); + // to let the flow solver that saving pressure_k and temperature_k is necessary (for the fixed-stress porosity terms) + flowSolver()->enableFixedStressPoromechanicsUpdate(); + } + else + { + // to let the solid mechanics solver that there is a pressure and temperature RHS in the mechanics solve + solidMechanicsSolver()->enableFixedStressPoromechanicsUpdate(); + // to let the flow solver that saving pressure_k and temperature_k is necessary (for the fixed-stress porosity terms) + flowSolver()->enableFixedStressPoromechanicsUpdate(); + } } if( m_stabilizationType == stabilization::StabilizationType::Global || m_stabilizationType == stabilization::StabilizationType::Local ) @@ -652,8 +668,10 @@ class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER } } - virtual void mapSolutionBetweenSolvers( DomainPartition & domain, integer const solverType ) override + virtual void mapSolutionBetweenSolvers( real64 const & dt, DomainPartition & domain, integer const solverType ) override { + GEOS_UNUSED_VAR( dt ); + GEOS_MARK_FUNCTION; /// After the flow solver @@ -666,8 +684,12 @@ class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER if( solverType == static_cast< integer >( SolverType::SolidMechanics ) && !m_performStressInitialization ) // do not update during poromechanics initialization { - // compute the average of the mean total stress increment over quadrature points - averageMeanTotalStressIncrement( domain ); + if( flowSolver()->getCatalogName() != "SinglePhaseReactiveTransport" ) // For now, Biot Poromechanics is not considered for + // ChemoMechanics + { + // compute the average of the mean total stress increment over quadrature points + averageMeanTotalStressIncrement( domain ); + } this->template forDiscretizationOnMeshTargets<>( domain.getMeshBodies(), [&]( string const &, MeshLevel & mesh, @@ -696,7 +718,11 @@ class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER if( solverType == static_cast< integer >( SolverType::SolidMechanics ) && this->getNonlinearSolverParameters().m_nonlinearAccelerationType== NonlinearSolverParameters::NonlinearAccelerationType::Aitken ) { - recordAverageMeanTotalStressIncrement( domain, m_s2_tilde ); + if( flowSolver()->getCatalogName() != "SinglePhaseReactiveTransport" ) // For now, Biot Poromechanics is not considered for + // ChemoMechanics + { + recordAverageMeanTotalStressIncrement( domain, m_s2_tilde ); + } } } diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.cpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.cpp index 468117b09dd..e643617bd35 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.cpp @@ -185,6 +185,48 @@ void SinglePhasePoromechanics< SinglePhaseReservoirAndWells<>, SolidMechanicsLag flowSolver()->assembleCouplingTerms( time_n, dt, domain, dofManager, localMatrix, localRhs ); } +template<> +void SinglePhasePoromechanics< SinglePhaseReactiveTransport, SolidMechanicsLagrangianFEM >::mapSolutionBetweenSolvers( real64 const & dt, + DomainPartition & domain, + integer const solverType ) +{ + GEOS_MARK_FUNCTION; + + Base::mapSolutionBetweenSolvers( dt, domain, solverType ); + + /// After the solid mechanics solver + if( solverType == static_cast< integer >( Base::SolverType::SolidMechanics ) + && !this->m_performStressInitialization ) // do not update during poromechanics initialization + { + this->template forDiscretizationOnMeshTargets<>( domain.getMeshBodies(), [&]( string const &, + MeshLevel & mesh, + string_array const & regionNames ) + { + + mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const, + auto & subRegion ) + { + // update mass after porosity change due to mechanics solve + this->flowSolver()->updateMass( subRegion ); + } ); + } ); + } + + if( solverType == static_cast< integer >( SolverType::Flow ) ) + { + this->template forDiscretizationOnMeshTargets<>( domain.getMeshBodies(), [&]( string const &, + MeshLevel & mesh, + string_array const & regionNames ) + { + mesh.getElemManager().forElementSubRegions( regionNames, [&]( localIndex const, + ElementSubRegionBase & subRegion ) + { + flowSolver()->updateKineticReactionMolarIncrements( dt, subRegion ); + } ); + } ); + } +} + template< typename FLOW_SOLVER, typename MECHANICS_SOLVER > void SinglePhasePoromechanics< FLOW_SOLVER, MECHANICS_SOLVER >::assembleElementBasedTerms( real64 const time_n, real64 const dt, @@ -355,9 +397,12 @@ template class SinglePhasePoromechanics< SinglePhaseReservoirAndWells<> >; template class SinglePhasePoromechanics< SinglePhaseReservoirAndWells<>, SolidMechanicsLagrangeContact >; template class SinglePhasePoromechanics< SinglePhaseReservoirAndWells<>, SolidMechanicsAugmentedLagrangianContact >; //template class SinglePhasePoromechanics< SinglePhaseReservoirAndWells<>, SolidMechanicsEmbeddedFractures >; +template class SinglePhasePoromechanics< SinglePhaseReactiveTransport >; namespace { +typedef SinglePhasePoromechanics< SinglePhaseReactiveTransport > SinglePhaseChemomechanics; +REGISTER_CATALOG_ENTRY( PhysicsSolverBase, SinglePhaseChemomechanics, string const &, Group * const ) typedef SinglePhasePoromechanics< SinglePhaseReservoirAndWells<> > SinglePhaseReservoirPoromechanics; REGISTER_CATALOG_ENTRY( PhysicsSolverBase, SinglePhaseReservoirPoromechanics, string const &, Group * const ) typedef SinglePhasePoromechanics<> SinglePhasePoromechanics; diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp index d7e40171f79..12318c46600 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp @@ -22,6 +22,7 @@ #include "physicsSolvers/multiphysics/PoromechanicsSolver.hpp" #include "physicsSolvers/fluidFlow/SinglePhaseBase.hpp" +#include "physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.hpp" #include "physicsSolvers/multiphysics/SinglePhaseReservoirAndWells.hpp" namespace geos @@ -122,11 +123,11 @@ class SinglePhasePoromechanics : public PoromechanicsSolver< FLOW_SOLVER, MECHAN GEOS_ERROR( GEOS_FMT( "{}: MGR strategy is not implemented for {}", this->getName(), this->getCatalogName())); } - virtual void mapSolutionBetweenSolvers( DomainPartition & domain, integer const solverType ) override + virtual void mapSolutionBetweenSolvers( real64 const & dt, DomainPartition & domain, integer const solverType ) override { GEOS_MARK_FUNCTION; - Base::mapSolutionBetweenSolvers( domain, solverType ); + Base::mapSolutionBetweenSolvers( dt, domain, solverType ); /// After the solid mechanics solver if( solverType == static_cast< integer >( Base::SolverType::SolidMechanics ) diff --git a/src/coreComponents/physicsSolvers/multiphysics/kernelSpecs.json b/src/coreComponents/physicsSolvers/multiphysics/kernelSpecs.json index bf248e898f3..063cf06b507 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/kernelSpecs.json +++ b/src/coreComponents/physicsSolvers/multiphysics/kernelSpecs.json @@ -61,6 +61,7 @@ "PorousSolid, CarmanKozenyPermeability>", "PorousSolid, CarmanKozenyPermeability>", "PorousSolid, CarmanKozenyPermeability>", + "PorousSolid", "PorousDamageSolid>", "PorousDamageSolid>", "PorousDamageSolid>" @@ -136,7 +137,8 @@ "PorousSolid", "PorousSolid", "PorousSolid", - "PorousSolid" + "PorousSolid", + "PorousSolid" ], "FE_TYPE": [ "H1_Hexahedron_Lagrange1_GaussLegendre2", diff --git a/src/coreComponents/physicsSolvers/solidMechanics/CMakeLists.txt b/src/coreComponents/physicsSolvers/solidMechanics/CMakeLists.txt index c6f9a5b0035..52dbfd4ced4 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/CMakeLists.txt +++ b/src/coreComponents/physicsSolvers/solidMechanics/CMakeLists.txt @@ -25,6 +25,8 @@ set( solidMechanicsSolvers_headers kernels/SolidMechanicsLagrangianFEMKernels.hpp SolidMechanicsMPM.hpp MPMSolverFields.hpp + kernels/ExplicitChemoMechanics.hpp + kernels/ExplicitChemoMechanics_impl.hpp kernels/ExplicitFiniteStrain.hpp kernels/ExplicitFiniteStrain_impl.hpp kernels/ExplicitMPM.hpp @@ -85,6 +87,7 @@ set( kernelTemplateFileList "" ) list( APPEND kernelTemplateFileList kernels/SolidMechanicsKernels.cpp.template kernels/SolidMechanicsFixedStressThermoPoromechanicsKernels.cpp.template + kernels/SolidMechanicsExplicitChemoMechanicsKernels.cpp.template contact/kernels/SolidMechanicsALMContactPorousKernels.cpp.template ) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp index 2e2ab2d365d..a8e6b95d600 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp @@ -25,6 +25,7 @@ #include "kernels/ExplicitSmallStrain.hpp" #include "kernels/ExplicitFiniteStrain.hpp" #include "kernels/FixedStressThermoPoromechanics.hpp" +#include "kernels/ExplicitChemoMechanics.hpp" #include "common/GEOS_RAJA_Interface.hpp" #include "constitutive/ConstitutiveManager.hpp" @@ -49,6 +50,7 @@ #include "physicsSolvers/LogLevelsInfo.hpp" #include "physicsSolvers/solidMechanics/kernels/SolidMechanicsKernelsDispatchTypeList.hpp" #include "physicsSolvers/solidMechanics/kernels/SolidMechanicsFixedStressThermoPoromechanicsKernelsDispatchTypeList.hpp" +#include "physicsSolvers/solidMechanics/kernels/SolidMechanicsExplicitChemoMechanicsKernelsDispatchTypeList.hpp" #include "physicsSolvers/fluidFlow/FlowSolverBase.hpp" namespace geos @@ -70,7 +72,8 @@ SolidMechanicsLagrangianFEM::SolidMechanicsLagrangianFEM( const string & name, m_maxNumResolves( 10 ), m_strainTheory( 0 ), m_isFixedStressPoromechanicsUpdate( false ), - m_performStressInitialization( false ) + m_performStressInitialization( false ), + m_isExplicitChemomechanicsUpdate( false ) { registerWrapper( viewKeyStruct::newmarkGammaString(), &m_newmarkGamma ). @@ -1154,7 +1157,7 @@ void SolidMechanicsLagrangianFEM::assembleSystem( real64 const GEOS_UNUSED_PARAM MeshLevel & mesh, string_array const & regionNames ) { - if( m_isFixedStressPoromechanicsUpdate || m_performStressInitialization ) + if( m_isFixedStressPoromechanicsUpdate || (m_performStressInitialization && !m_isExplicitChemomechanicsUpdate) ) { set< string > poromechanicsRegions; set< string > mechanicsRegions; @@ -1207,6 +1210,22 @@ void SolidMechanicsLagrangianFEM::assembleSystem( real64 const GEOS_UNUSED_PARAM m_maxForce = LvArray::math::max( mechanicsMaxForce, poromechanicsMaxForce ); } + else if( m_isExplicitChemomechanicsUpdate ) + { + + // first pass for coupled poromechanics regions + real64 const chemomechanicsMaxForce= assemblyLaunch< SolidMechanicsExplicitChemoMechanicsKernelsDispatchTypeList, + solidMechanicsLagrangianFEMKernels::ExplicitChemoMechanicsFactory >( mesh, + dofManager, + regionNames, + FlowSolverBase::viewKeyStruct::solidNamesString(), + localMatrix, + localRhs, + dt ); + + + m_maxForce = chemomechanicsMaxForce; + } else { if( m_timeIntegrationOption == TimeIntegrationOption::QuasiStatic ) @@ -1587,6 +1606,11 @@ void SolidMechanicsLagrangianFEM::enableFixedStressPoromechanicsUpdate() m_isFixedStressPoromechanicsUpdate = true; } +void SolidMechanicsLagrangianFEM::enableExplicitChemomechanicsUpdate() +{ + m_isExplicitChemomechanicsUpdate = true; +} + void SolidMechanicsLagrangianFEM::saveSequentialIterationState( DomainPartition & GEOS_UNUSED_PARAM( domain ) ) { // nothing to save diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp index 89e149c5ef4..7ce8c1a220e 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp @@ -228,6 +228,8 @@ class SolidMechanicsLagrangianFEM : public PhysicsSolverBase void enableFixedStressPoromechanicsUpdate(); + void enableExplicitChemomechanicsUpdate(); + virtual void saveSequentialIterationState( DomainPartition & domain ) override; struct viewKeyStruct : PhysicsSolverBase::viewKeyStruct @@ -311,6 +313,9 @@ class SolidMechanicsLagrangianFEM : public PhysicsSolverBase /// Flag to indicate that the solver is going to perform stress initialization bool m_performStressInitialization; + /// Flag to indicate that the solver is running with explicit chemomechancis update + bool m_isExplicitChemomechanicsUpdate; + /// Rigid body modes; TODO remove mutable hack mutable array1d< ParallelVector > m_rigidBodyModes; diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernelSpecs.json b/src/coreComponents/physicsSolvers/solidMechanics/kernelSpecs.json index 3caafa8247e..2809914896b 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/kernelSpecs.json +++ b/src/coreComponents/physicsSolvers/solidMechanics/kernelSpecs.json @@ -29,6 +29,7 @@ "constants": [ [ "ExplicitSmallStrainPolicy", "geos::parallelDevicePolicy< GEOS_BLOCK_SIZE >" ], [ "ExplicitFiniteStrainPolicy", "geos::parallelDevicePolicy< GEOS_BLOCK_SIZE >" ], + [ "ExplicitChemoMechanicsPolicy", "geos::parallelDevicePolicy< GEOS_BLOCK_SIZE >" ], [ "FixedStressThermoPoromechanicsPolicy", "geos::parallelDevicePolicy< GEOS_BLOCK_SIZE >" ], [ "ImplicitSmallStrainNewmarkPolicy", "geos::parallelDevicePolicy< GEOS_BLOCK_SIZE >" ], [ "ImplicitSmallStrainQuasiStaticPolicy", "geos::parallelDevicePolicy< GEOS_BLOCK_SIZE >" ] @@ -101,7 +102,8 @@ "PorousSolid, ConstantPermeability>", "PorousSolid, CarmanKozenyPermeability>", "PorousSolid, ConstantPermeability>", - "PorousSolid, CarmanKozenyPermeability>" + "PorousSolid, CarmanKozenyPermeability>", + "PorousSolid" ], "FE_TYPE": [ "H1_Hexahedron_Lagrange1_GaussLegendre2", @@ -139,7 +141,44 @@ ], "CONSTITUTIVE_TYPE": [ "PorousSolid", - "PorousSolid" + "PorousSolid", + "PorousSolid" + ], + "FE_TYPE": [ + "H1_Hexahedron_Lagrange1_GaussLegendre2", + "H1_Wedge_Lagrange1_Gauss6", + "H1_Tetrahedron_Lagrange1_Gauss1", + "H1_Pyramid_Lagrange1_Gauss5" + ] + }, + "explicit": [] + }, + + "SolidMechanicsExplicitChemoMechanicsKernels": { + "vars": [ + "SUBREGION_TYPE", + "CONSTITUTIVE_TYPE", + "FE_TYPE" + ], + "constants": [ + [ "ExplicitSmallStrainPolicy", "geos::parallelDevicePolicy< GEOS_BLOCK_SIZE >" ], + [ "ExplicitFiniteStrainPolicy", "geos::parallelDevicePolicy< GEOS_BLOCK_SIZE >" ], + [ "ExplicitChemoMechanicsPolicy", "geos::parallelDevicePolicy< GEOS_BLOCK_SIZE >" ], + [ "FixedStressThermoPoromechanicsPolicy", "geos::parallelDevicePolicy< GEOS_BLOCK_SIZE >" ], + [ "ImplicitSmallStrainNewmarkPolicy", "geos::parallelDevicePolicy< GEOS_BLOCK_SIZE >" ], + [ "ImplicitSmallStrainQuasiStaticPolicy", "geos::parallelDevicePolicy< GEOS_BLOCK_SIZE >" ] + ], + "combinations": { + "SUBREGION_TYPE": [ + "CellElementSubRegion" + ], + "CONSTITUTIVE_TYPE": [ + "EigenstrainReactiveSolid", + "EigenstrainReactiveSolid", + "EigenstrainReactiveSolid", + "EigenstrainReactiveSolid", + "PorousReactiveSolid", + "PorousReactiveSolid" ], "FE_TYPE": [ "H1_Hexahedron_Lagrange1_GaussLegendre2", diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernels/ExplicitChemoMechanics.hpp b/src/coreComponents/physicsSolvers/solidMechanics/kernels/ExplicitChemoMechanics.hpp new file mode 100644 index 00000000000..e76ce456899 --- /dev/null +++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/ExplicitChemoMechanics.hpp @@ -0,0 +1,234 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file ExplicitChemoMechanics.hpp + */ + +#ifndef GEOS_PHYSICSSOLVERS_SOLIDMECHANICS_KERNELS_EXPLICITCHEMOMECHANICS_HPP_ +#define GEOS_PHYSICSSOLVERS_SOLIDMECHANICS_KERNELS_EXPLICITCHEMOMECHANICS_HPP_ + +#include "finiteElement/kernelInterface/ImplicitKernelBase.hpp" + +namespace geos +{ + +namespace solidMechanicsLagrangianFEMKernels +{ + +/** + * @brief Implements kernels for solving the solid part of the chemomechanics problem. + * @copydoc geos::finiteElement::ImplicitKernelBase + * @tparam NUM_NODES_PER_ELEM The number of nodes per element for the + * @p SUBREGION_TYPE. + * @tparam UNUSED An unused parameter since we are assuming that the test and + * trial space have the same number of support points. + * + * ### ExplicitChemoMechanics Description + * Implements the KernelBase interface functions required for using the + * effective stress for the integration of the stress divergence. This is + * templated on one of the "finite element kernel application" functions + * such as geos::finiteElement::RegionBasedKernelApplication. + */ +template< typename SUBREGION_TYPE, + typename CONSTITUTIVE_TYPE, + typename FE_TYPE > +class ExplicitChemoMechanics : + public finiteElement::ImplicitKernelBase< SUBREGION_TYPE, + CONSTITUTIVE_TYPE, + FE_TYPE, + 3, + 3 > +{ +public: + /// Alias for the base class; + using Base = finiteElement::ImplicitKernelBase< SUBREGION_TYPE, + CONSTITUTIVE_TYPE, + FE_TYPE, + 3, + 3 >; + + /// Maximum number of nodes per element, which is equal to the maxNumTestSupportPointPerElem and + /// maxNumTrialSupportPointPerElem by definition. When the FE_TYPE is not a Virtual Element, this + /// will be the actual number of nodes per element. + static constexpr int numNodesPerElem = Base::maxNumTestSupportPointsPerElem; + using Base::numDofPerTestSupportPoint; + using Base::numDofPerTrialSupportPoint; + using Base::m_dofNumber; + using Base::m_dofRankOffset; + using Base::m_matrix; + using Base::m_rhs; + using Base::m_elemsToNodes; + using Base::m_constitutiveUpdate; + using Base::m_finiteElementSpace; + using Base::m_meshData; + using Base::m_dt; + + /** + * @brief Constructor + * @copydoc geos::finiteElement::ImplicitKernelBase::ImplicitKernelBase + * @param inputGravityVector The gravity vector. + */ + ExplicitChemoMechanics( NodeManager const & nodeManager, + EdgeManager const & edgeManager, + FaceManager const & faceManager, + localIndex const targetRegionIndex, + SUBREGION_TYPE const & elementSubRegion, + FE_TYPE const & finiteElementSpace, + CONSTITUTIVE_TYPE & inputConstitutiveType, + arrayView1d< globalIndex const > const inputDofNumber, + globalIndex const rankOffset, + CRSMatrixView< real64, globalIndex const > const inputMatrix, + arrayView1d< real64 > const inputRhs, + real64 const inputDt, + real64 const (&inputGravityVector)[3] ); + + //***************************************************************************** + /** + * @class StackVariables + * @copydoc geos::finiteElement::ImplicitKernelBase::StackVariables + * + * Adds a stack array for the displacement, incremental displacement, and the + * constitutive stiffness. + */ + struct StackVariables : public Base::StackVariables + { +public: + + /// Constructor. + GEOS_HOST_DEVICE + StackVariables(): + Base::StackVariables(), + xLocal(), + u_local(), + uhat_local(), + constitutiveStiffness() + {} + + /// C-array stack storage for element local the nodal positions. + real64 xLocal[ numNodesPerElem ][ 3 ]; + + /// Stack storage for the element local nodal displacement + real64 u_local[numNodesPerElem][numDofPerTrialSupportPoint]; + + /// Stack storage for the element local nodal incremental displacement + real64 uhat_local[numNodesPerElem][numDofPerTrialSupportPoint]; + + /// Stack storage for the constitutive stiffness at a quadrature point. + real64 constitutiveStiffness[ 6 ][ 6 ]; + }; + //***************************************************************************** + + /** + * @brief Copy global values from primary field to a local stack array. + * @copydoc ::geos::finiteElement::ImplicitKernelBase::setup + * + * For the ExplicitChemoMechanics implementation, global values from the displacement, + * incremental displacement, and degree of freedom numbers are placed into + * element local stack storage. + */ + GEOS_HOST_DEVICE + void setup( localIndex const k, + StackVariables & stack ) const; + + /** + * @copydoc geos::finiteElement::KernelBase::quadraturePointKernel + * For solid mechanics kernels, the strain increment is calculated, and the + * constitutive update is called. In addition, the constitutive stiffness + * stack variable is filled by the constitutive model. + */ + GEOS_HOST_DEVICE + void quadraturePointKernel( localIndex const k, + localIndex const q, + StackVariables & stack ) const; + + /** + * @copydoc geos::finiteElement::ImplicitKernelBase::complete + */ + GEOS_HOST_DEVICE + GEOS_FORCE_INLINE + real64 complete( localIndex const k, + StackVariables & stack ) const; + + /** + * @copydoc geos::finiteElement::KernelBase::kernelLaunch + */ + template< typename POLICY, + typename KERNEL_TYPE > + static real64 + kernelLaunch( localIndex const numElems, + KERNEL_TYPE const & kernelComponent ); + +protected: + /// The array containing the nodal position array. + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const m_X; + + /// The rank-global displacement array. + arrayView2d< real64 const, nodes::TOTAL_DISPLACEMENT_USD > const m_disp; + + /// The rank-global incremental displacement array. + arrayView2d< real64 const, nodes::INCR_DISPLACEMENT_USD > const m_uhat; + + /// The gravity vector. + real64 const m_gravityVector[3]; + + /// The rank global bulk density + arrayView2d< real64 const > const m_bulkDensity; + + /// The rank-global fluid pressure arrays. + arrayView1d< real64 const > const m_pressure; + arrayView1d< real64 const > const m_pressure_n; + + /// The rank-global initial temperature array + arrayView1d< real64 const > const m_initialTemperature; + + /// The rank-global temperature arrays. + arrayView1d< real64 const > const m_temperature; + arrayView1d< real64 const > const m_temperature_n; + + /// The rank-global mineral reaction molar increments arrays. + arrayView2d< real64 const, compflow::USD_COMP > const m_mineralReactionMolarIncrements; + + /** + * @brief Get a parameter representative of the stiffness, used as physical scaling for the + * stabilization matrix. + * @param[in] k Element index. + * @return A parameter representative of the stiffness matrix dstress/dstrain + */ + GEOS_HOST_DEVICE + GEOS_FORCE_INLINE + real64 computeStabilizationScaling( localIndex const k ) const + { + // TODO: generalize this to other constitutive models (currently we assume linear elasticity). + return 2.0 * m_constitutiveUpdate.getShearModulus( k ); + } +}; + +/// The factory used to construct a ExplicitChemoMechanics kernel. +using ExplicitChemoMechanicsFactory = finiteElement::KernelFactory< ExplicitChemoMechanics, + arrayView1d< globalIndex const > const, + globalIndex, + CRSMatrixView< real64, globalIndex const > const, + arrayView1d< real64 > const, + real64 const, + real64 const (&)[3] >; + +} // namespace solidMechanicsLagrangianFEMKernels + +} // namespace geos + +#include "finiteElement/kernelInterface/SparsityKernelBase.hpp" + +#endif // GEOS_PHYSICSSOLVERS_SOLIDMECHANICS_KERNELS_EXPLICITCHEMOMECHANICS_HPP_ diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernels/ExplicitChemoMechanics_impl.hpp b/src/coreComponents/physicsSolvers/solidMechanics/kernels/ExplicitChemoMechanics_impl.hpp new file mode 100644 index 00000000000..78659936f08 --- /dev/null +++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/ExplicitChemoMechanics_impl.hpp @@ -0,0 +1,233 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file ExplicitChemoMechanics_impl.hpp + */ + +#ifndef GEOS_PHYSICSSOLVERS_SOLIDMECHANICS_KERNELS_EXPLICITCHEMOMECHANICS_IMPL_HPP_ +#define GEOS_PHYSICSSOLVERS_SOLIDMECHANICS_KERNELS_EXPLICITCHEMOMECHANICS_IMPL_HPP_ + +#include "ExplicitChemoMechanics.hpp" +#include "finiteElement/elementFormulations/FiniteElementOperators.hpp" +#include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" +#include "physicsSolvers/fluidFlow/SinglePhaseReactiveTransportFields.hpp" +#include "physicsSolvers/multiphysics/PoromechanicsFields.hpp" +#include "physicsSolvers/solidMechanics/SolidMechanicsFields.hpp" + +namespace geos +{ + +namespace solidMechanicsLagrangianFEMKernels +{ + +template< typename SUBREGION_TYPE, + typename CONSTITUTIVE_TYPE, + typename FE_TYPE > + +ExplicitChemoMechanics< SUBREGION_TYPE, CONSTITUTIVE_TYPE, FE_TYPE >:: +ExplicitChemoMechanics( NodeManager const & nodeManager, + EdgeManager const & edgeManager, + FaceManager const & faceManager, + localIndex const targetRegionIndex, + SUBREGION_TYPE const & elementSubRegion, + FE_TYPE const & finiteElementSpace, + CONSTITUTIVE_TYPE & inputConstitutiveType, + arrayView1d< globalIndex const > const inputDofNumber, + globalIndex const rankOffset, + CRSMatrixView< real64, globalIndex const > const inputMatrix, + arrayView1d< real64 > const inputRhs, + real64 const inputDt, + real64 const (&inputGravityVector)[3] ): + Base( nodeManager, + edgeManager, + faceManager, + targetRegionIndex, + elementSubRegion, + finiteElementSpace, + inputConstitutiveType, + inputDofNumber, + rankOffset, + inputMatrix, + inputRhs, + inputDt ), + m_X( nodeManager.referencePosition()), + m_disp( nodeManager.getField< fields::solidMechanics::totalDisplacement >() ), + m_uhat( nodeManager.getField< fields::solidMechanics::incrementalDisplacement >() ), + m_gravityVector{ inputGravityVector[0], inputGravityVector[1], inputGravityVector[2] }, + m_bulkDensity( elementSubRegion.template getField< fields::poromechanics::bulkDensity >() ), + m_pressure( elementSubRegion.template getField< fields::flow::pressure >() ), + m_pressure_n( elementSubRegion.template getField< fields::flow::pressure_n >() ), + m_initialTemperature( elementSubRegion.template getField< fields::flow::initialTemperature >() ), + m_temperature( elementSubRegion.template getField< fields::flow::temperature >() ), + m_temperature_n( elementSubRegion.template getField< fields::flow::temperature_n >() ), + m_mineralReactionMolarIncrements( elementSubRegion.template getField< fields::flow::kineticReactionMolarIncrements >() ) +{} + +template< typename SUBREGION_TYPE, + typename CONSTITUTIVE_TYPE, + typename FE_TYPE > +GEOS_HOST_DEVICE +GEOS_FORCE_INLINE +void ExplicitChemoMechanics< SUBREGION_TYPE, CONSTITUTIVE_TYPE, FE_TYPE >:: +setup( localIndex const k, + StackVariables & stack ) const +{ + m_finiteElementSpace.template setup< FE_TYPE >( k, m_meshData, stack.feStack ); + localIndex const numSupportPoints = + m_finiteElementSpace.getNumSupportPoints( stack.feStack ); + stack.numRows = 3 * numSupportPoints; + stack.numCols = stack.numRows; + + for( localIndex a = 0; a < numSupportPoints; ++a ) + { + localIndex const localNodeIndex = m_elemsToNodes( k, a ); + + for( int i = 0; i < 3; ++i ) + { + stack.xLocal[ a ][ i ] = m_X[ localNodeIndex ][ i ]; + stack.u_local[ a ][i] = m_disp[ localNodeIndex ][i]; + stack.uhat_local[ a ][i] = m_uhat[ localNodeIndex ][i]; + stack.localRowDofIndex[a*3+i] = m_dofNumber[localNodeIndex]+i; + stack.localColDofIndex[a*3+i] = m_dofNumber[localNodeIndex]+i; + } + } + + // Add stabilization to block diagonal parts of the local jacobian + // (this is a no-operation with FEM classes) + real64 const stabilizationScaling = computeStabilizationScaling( k ); + m_finiteElementSpace.template addGradGradStabilizationMatrix + < FE_TYPE, numDofPerTrialSupportPoint, true >( stack.feStack, + stack.localJacobian, + -stabilizationScaling ); +} + +template< typename SUBREGION_TYPE, + typename CONSTITUTIVE_TYPE, + typename FE_TYPE > +GEOS_HOST_DEVICE +GEOS_FORCE_INLINE +void ExplicitChemoMechanics< SUBREGION_TYPE, CONSTITUTIVE_TYPE, FE_TYPE >:: +quadraturePointKernel( localIndex const k, + localIndex const q, + StackVariables & stack ) const +{ + real64 dNdX[ numNodesPerElem ][ 3 ]; + real64 const detJxW = FE_TYPE::calcGradN( q, stack.xLocal, + stack.feStack, dNdX ); + + real64 strainInc[6] = {0}; + real64 totalStress[6] = {0}; + + typename CONSTITUTIVE_TYPE::KernelWrapper::DiscretizationOps stiffness; + + finiteElement::feOps::symmetricGradient( dNdX, stack.uhat_local, strainInc ); + + // Evaluate total stress and its derivatives + // TODO: allow for a customization of the kernel to pass the average pressure to the small strain update (to account for cap pressure + // later) + m_constitutiveUpdate.smallStrainUpdateChemoMechanicsFixedStress( k, q, + m_dt, + m_pressure[k], + m_pressure_n[k], + m_temperature[k], + m_temperature_n[k], + m_initialTemperature[k], + m_mineralReactionMolarIncrements[k], + strainInc, + totalStress, + stiffness ); + + for( localIndex i=0; i<6; ++i ) + { + totalStress[i] *= -detJxW; + } + + // Here we consider the bodyForce is purely from the solid + // Warning: here, we lag (in iteration) the displacement dependence of bulkDensity + real64 const gravityForce[3] = { m_gravityVector[0] * m_bulkDensity( k, q )* detJxW, + m_gravityVector[1] * m_bulkDensity( k, q )* detJxW, + m_gravityVector[2] * m_bulkDensity( k, q )* detJxW }; + + real64 N[numNodesPerElem]; + FE_TYPE::calcN( q, stack.feStack, N ); + finiteElement::feOps::plusGradNajAijPlusNaFi( dNdX, + totalStress, + N, + gravityForce, + reinterpret_cast< real64 (&)[numNodesPerElem][3] >(stack.localResidual) ); + real64 const stabilizationScaling = computeStabilizationScaling( k ); + m_finiteElementSpace.template + addEvaluatedGradGradStabilizationVector< FE_TYPE, + numDofPerTrialSupportPoint >( stack.feStack, + stack.uhat_local, + reinterpret_cast< real64 (&)[numNodesPerElem][3] >(stack.localResidual), + -stabilizationScaling ); + stiffness.template upperBTDB< numNodesPerElem >( dNdX, -detJxW, stack.localJacobian ); +} + +template< typename SUBREGION_TYPE, + typename CONSTITUTIVE_TYPE, + typename FE_TYPE > +GEOS_HOST_DEVICE +GEOS_FORCE_INLINE +real64 ExplicitChemoMechanics< SUBREGION_TYPE, CONSTITUTIVE_TYPE, FE_TYPE >:: +complete( localIndex const k, + StackVariables & stack ) const +{ + GEOS_UNUSED_VAR( k ); + real64 maxForce = 0; + + // TODO: Does this work if BTDB is non-symmetric? + CONSTITUTIVE_TYPE::KernelWrapper::DiscretizationOps::template fillLowerBTDB< numNodesPerElem >( stack.localJacobian ); + localIndex const numSupportPoints = + m_finiteElementSpace.getNumSupportPoints( stack.feStack ); + for( int localNode = 0; localNode < numSupportPoints; ++localNode ) + { + for( int dim = 0; dim < numDofPerTestSupportPoint; ++dim ) + { + localIndex const dof = + LvArray::integerConversion< localIndex >( stack.localRowDofIndex[ numDofPerTestSupportPoint * localNode + dim ] - m_dofRankOffset ); + if( dof < 0 || dof >= m_matrix.numRows() ) + continue; + m_matrix.template addToRowBinarySearchUnsorted< parallelDeviceAtomic >( dof, + stack.localRowDofIndex, + stack.localJacobian[ numDofPerTestSupportPoint * localNode + dim ], + stack.numRows ); + + RAJA::atomicAdd< parallelDeviceAtomic >( &m_rhs[ dof ], stack.localResidual[ numDofPerTestSupportPoint * localNode + dim ] ); + maxForce = fmax( maxForce, fabs( stack.localResidual[ numDofPerTestSupportPoint * localNode + dim ] ) ); + } + } + return maxForce; +} + +template< typename SUBREGION_TYPE, + typename CONSTITUTIVE_TYPE, + typename FE_TYPE > +template< typename POLICY, + typename KERNEL_TYPE > +real64 +ExplicitChemoMechanics< SUBREGION_TYPE, CONSTITUTIVE_TYPE, FE_TYPE >::kernelLaunch( localIndex const numElems, + KERNEL_TYPE const & kernelComponent ) +{ + return Base::template kernelLaunch< POLICY, KERNEL_TYPE >( numElems, kernelComponent ); +} + +} // namespace solidMechanicsLagrangianFEMKernels + +} // namespace geos + +#endif // GEOS_PHYSICSSOLVERS_SOLIDMECHANICS_KERNELS_EXPLICITCHEMOMECHANICS_IMPL_HPP_ diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernels/SolidMechanicsExplicitChemoMechanicsKernels.cpp.template b/src/coreComponents/physicsSolvers/solidMechanics/kernels/SolidMechanicsExplicitChemoMechanicsKernels.cpp.template new file mode 100644 index 00000000000..48196f470bc --- /dev/null +++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/SolidMechanicsExplicitChemoMechanicsKernels.cpp.template @@ -0,0 +1,38 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +#include "physicsSolvers/solidMechanics/kernels/ExplicitChemoMechanics_impl.hpp" + +using ExplicitChemoMechanicsPolicy = @ExplicitChemoMechanicsPolicy@; + +#define INSTANTIATION( NAME )\ +template class NAME < @SUBREGION_TYPE@, @CONSTITUTIVE_TYPE@, @FE_TYPE@_impl >; \ +template real64 NAME < @SUBREGION_TYPE@, @CONSTITUTIVE_TYPE@, @FE_TYPE@_impl >::kernelLaunch< NAME##Policy, \ + NAME < @SUBREGION_TYPE@, @CONSTITUTIVE_TYPE@, @FE_TYPE@_impl > > \ + ( localIndex const, \ + NAME < @SUBREGION_TYPE@, @CONSTITUTIVE_TYPE@, @FE_TYPE@_impl > const & ); \ + + +namespace geos +{ +using namespace constitutive; +using namespace finiteElement; +namespace solidMechanicsLagrangianFEMKernels +{ + INSTANTIATION( ExplicitChemoMechanics ) +} +} + +#undef INSTANTIATION diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernels/policies.hpp.in b/src/coreComponents/physicsSolvers/solidMechanics/kernels/policies.hpp.in index c14654cf0c7..bafc3a064b9 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/kernels/policies.hpp.in +++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/policies.hpp.in @@ -18,6 +18,7 @@ using ExplicitSmallStrainPolicy = @ExplicitSmallStrainPolicy@; using ExplicitFiniteStrainPolicy = @ExplicitFiniteStrainPolicy@; +using ExplicitChemoMechanicsPolicy = @ExplicitChemoMechanicsPolicy@; using FixedStressThermoPoromechanicsPolicy = @FixedStressThermoPoromechanicsPolicy@; using ImplicitSmallStrainNewmarkPolicy = @ImplicitSmallStrainNewmarkPolicy@; using ImplicitSmallStrainQuasiStaticPolicy = @ImplicitSmallStrainQuasiStaticPolicy@;