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/src/coreComponents/constitutive/CMakeLists.txt b/src/coreComponents/constitutive/CMakeLists.txt index 7606cfe25f6..43c82a055b3 100644 --- a/src/coreComponents/constitutive/CMakeLists.txt +++ b/src/coreComponents/constitutive/CMakeLists.txt @@ -185,6 +185,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 +193,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 +209,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 @@ -332,22 +335,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..c83bbb8b0f1 100644 --- a/src/coreComponents/constitutive/ConstitutivePassThru.hpp +++ b/src/coreComponents/constitutive/ConstitutivePassThru.hpp @@ -36,15 +36,17 @@ #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" @@ -362,6 +364,36 @@ 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 > >::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 +438,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 +448,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 +497,8 @@ struct ConstitutivePassThru< CoupledSolidBase > CompressibleSolid< PressurePorosity, PressurePermeability >, CompressibleSolid< PressurePorosity, SlipDependentPermeability >, CompressibleSolid< PressurePorosity, WillisRichardsPermeability >, + EigenstrainReactiveSolid< ElasticIsotropic, ConstantPermeability >, + EigenstrainReactiveSolid< ElasticIsotropic, CarmanKozenyPermeability >, PorousSolid< DruckerPragerExtended, ConstantPermeability >, PorousSolid< ModifiedCamClay, ConstantPermeability >, PorousSolid< DelftEgg, ConstantPermeability >, @@ -490,10 +524,12 @@ struct ConstitutivePassThru< CoupledSolidBase > 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 +542,8 @@ struct ConstitutivePassThru< CoupledSolidBase > CompressibleSolid< PressurePorosity, PressurePermeability >, CompressibleSolid< PressurePorosity, SlipDependentPermeability >, CompressibleSolid< PressurePorosity, WillisRichardsPermeability >, + EigenstrainReactiveSolid< ElasticIsotropic, ConstantPermeability >, + EigenstrainReactiveSolid< ElasticIsotropic, CarmanKozenyPermeability >, PorousSolid< DruckerPragerExtended, ConstantPermeability >, PorousSolid< ModifiedCamClay, ConstantPermeability >, PorousSolid< DelftEgg, ConstantPermeability >, @@ -531,10 +569,12 @@ struct ConstitutivePassThru< CoupledSolidBase > 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..2504e76a821 160000 --- a/src/coreComponents/constitutive/HPCReact +++ b/src/coreComponents/constitutive/HPCReact @@ -1 +1 @@ -Subproject commit 5268dc5c2a59e9dae23f435fba9aecf5a5c5de33 +Subproject commit 2504e76a8211f01355044ebddb7005a9bd5d98be diff --git a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveFluidFields.hpp b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveFluidFields.hpp index 54e1c9a1fc2..62a4f9a2f07 100644 --- a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveFluidFields.hpp +++ b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveFluidFields.hpp @@ -39,7 +39,7 @@ DECLARE_FIELD( initialPrimarySpeciesConcentration, "initialPrimarySpeciesConcentration", array3dLayoutSpecies, 1e-16, - LEVEL_0, + LEVEL_1, WRITE_AND_READ, "initialPrimarySpeciesConcentration" ); @@ -55,7 +55,7 @@ DECLARE_FIELD( primarySpeciesAggregateConcentration_n, "primarySpeciesAggregateConcentration_n", array3dLayoutSpecies, 1e-16, - LEVEL_0, + NOPLOT, WRITE_AND_READ, "primarySpeciesAggregateConcentration at the previous timestep" ); @@ -63,7 +63,7 @@ DECLARE_FIELD( primarySpeciesMobileAggregateConcentration, "primarySpeciesMobileAggregateConcentration", array3dLayoutSpecies, 1e-16, - LEVEL_0, + LEVEL_1, WRITE_AND_READ, "primarySpeciesMobileAggregateConcentration" ); @@ -71,7 +71,7 @@ DECLARE_FIELD( dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentra "dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations", array4dLayoutSpecies_dC, 0, - LEVEL_0, + NOPLOT, WRITE_AND_READ, "Deivatives of primarySpeciesAggregateConcentration w.r.t log primary species concentration" ); @@ -79,7 +79,7 @@ DECLARE_FIELD( dPrimarySpeciesMobileAggregateConcentration_dLogPrimarySpeciesCon "dPrimarySpeciesMobileAggregateConcentration_dLogPrimarySpeciesConcentrations", array4dLayoutSpecies_dC, 0, - LEVEL_0, + NOPLOT, WRITE_AND_READ, "Deivatives of primarySpeciesMobileAggregateConcentration w.r.t log primary species concentration" ); @@ -95,7 +95,7 @@ DECLARE_FIELD( kineticReactionRates, "kineticReactionRates", array3dLayoutSpecies, 0, - NOPLOT, + LEVEL_0, WRITE_AND_READ, "kineticReactionRates" ); @@ -103,7 +103,7 @@ DECLARE_FIELD( aggregateSpeciesRates, "aggregateSpeciesRates", array3dLayoutSpecies, 0, - LEVEL_0, + LEVEL_1, WRITE_AND_READ, "aggregateSpeciesRates" ); @@ -111,7 +111,7 @@ DECLARE_FIELD( dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations, "dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations", array4dLayoutSpecies_dC, 0, - LEVEL_0, + NOPLOT, WRITE_AND_READ, "Deivatives of aggregate concentration rates w.r.t log primary species concentration" ); } diff --git a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.cpp b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.cpp index 69b7b750934..108308e2d3e 100644 --- a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.cpp +++ b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.cpp @@ -31,6 +31,9 @@ namespace reactivefluid { using namespace hpcReact::bulkGeneric; +using namespace hpcReact::geochemistry; +using namespace hpcReact::ChainGeneric; +using namespace hpcReact::MoMasBenchmark; template< typename BASE > ReactiveSinglePhaseFluid< BASE >:: @@ -69,6 +72,7 @@ deliverClone( string const & name, Group * const parent ) const newConstitutiveRelation.m_numPrimarySpecies = m_numPrimarySpecies; newConstitutiveRelation.m_numSecondarySpecies = m_numSecondarySpecies; newConstitutiveRelation.m_numKineticReactions = m_numKineticReactions; + newConstitutiveRelation.m_solventDensity = m_solventDensity; return clone; } @@ -81,39 +85,59 @@ void ReactiveSinglePhaseFluid< BASE >::postInputInitialization() switch( m_chemicalSystemType ) { case ChemicalSystemType::ultramafic: - m_numPrimarySpecies = 9; + m_numPrimarySpecies = 4; m_numSecondarySpecies = 16; m_numKineticReactions = 5; + m_solventDensity = ultramaficSystem.getSolventDensity(); break; case ChemicalSystemType::carbonate: m_numPrimarySpecies = 7; m_numSecondarySpecies = 10; m_numKineticReactions = 1; + m_solventDensity = carbonateSystem.getSolventDensity(); break; case ChemicalSystemType::carbonateAllEquilibrium: m_numPrimarySpecies = 7; m_numSecondarySpecies = 11; m_numKineticReactions = 0; + m_solventDensity = carbonateSystemAllEquilibrium.getSolventDensity(); + break; + + case ChemicalSystemType::serpentinization: + m_numPrimarySpecies = 3; + m_numSecondarySpecies = 0; + m_numKineticReactions = 3; + m_solventDensity = serpentinizationSystem.getSolventDensity(); + break; + + case ChemicalSystemType::kineticCarbonate: + m_numPrimarySpecies = 3; + m_numSecondarySpecies = 0; + m_numKineticReactions = 1; + m_solventDensity = kineticCarbonateSystem.getSolventDensity(); break; case ChemicalSystemType::chainSerialAllKinetic: m_numPrimarySpecies = 3; m_numSecondarySpecies = 0; m_numKineticReactions = 3; + m_solventDensity = serialAllKineticParams.getSolventDensity(); break; case ChemicalSystemType::momasMedium: m_numPrimarySpecies = 5; m_numSecondarySpecies = 9; m_numKineticReactions = 1; + m_solventDensity = mediumCaseParams.getSolventDensity(); break; default: m_numPrimarySpecies = 5; m_numSecondarySpecies = 7; m_numKineticReactions = 0; + m_solventDensity = easyCaseParams.getSolventDensity(); break; } } diff --git a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.hpp b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.hpp index 6346d030c4f..ed9ef32a5d7 100644 --- a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.hpp +++ b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.hpp @@ -52,6 +52,8 @@ enum class ChemicalSystemType : integer carbonate, carbonateAllEquilibrium, ultramafic, + serpentinization, + kineticCarbonate, momasEasy, momasMedium, chainSerialAllKinetic @@ -113,6 +115,8 @@ class ReactiveSinglePhaseFluid : public BASE integer numKineticReactions() const { return m_numKineticReactions; } + real64 solventDensity() const { return m_solventDensity; } + /** * @brief Kernel wrapper class for ReactiveSinglePhaseFluid. */ @@ -230,7 +234,8 @@ class ReactiveSinglePhaseFluid : public BASE typename ReactiveSinglePhaseFluid< BASE >::template ReactionKernelWrapper< hpcReact::geochemistry::carbonateSystemAllEquilibriumType >, typename ReactiveSinglePhaseFluid< BASE >::template ReactionKernelWrapper< hpcReact::ChainGeneric::serialAllKineticType >, typename ReactiveSinglePhaseFluid< BASE >::template ReactionKernelWrapper< hpcReact::MoMasBenchmark::mediumCaseType >, - typename ReactiveSinglePhaseFluid< BASE >::template ReactionKernelWrapper< hpcReact::MoMasBenchmark::easyCaseType > > + typename ReactiveSinglePhaseFluid< BASE >::template ReactionKernelWrapper< hpcReact::MoMasBenchmark::easyCaseType >, + typename ReactiveSinglePhaseFluid< BASE >::template ReactionKernelWrapper< hpcReact::geochemistry::kineticCarbonateSystemType > > createReactionKernelWrapper() const { using namespace hpcReact::geochemistry; @@ -282,6 +287,34 @@ class ReactiveSinglePhaseFluid : public BASE m_numSecondarySpecies, m_numKineticReactions, carbonateSystemAllEquilibrium ); + case ChemicalSystemType::serpentinization: + return ReactionKernelWrapper< serpentinizationSystemType >( 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, + serpentinizationSystem ); + case ChemicalSystemType::kineticCarbonate: + return ReactionKernelWrapper< kineticCarbonateSystemType >( 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, + kineticCarbonateSystem ); case ChemicalSystemType::chainSerialAllKinetic: return ReactionKernelWrapper< serialAllKineticType >( m_primarySpeciesAggregateConcentration, m_primarySpeciesMobileAggregateConcentration, @@ -365,6 +398,8 @@ class ReactiveSinglePhaseFluid : public BASE array4d< real64, constitutive::reactivefluid::LAYOUT_SPECIES_DC > m_dAggregateSpeciesRates_dLogPrimarySpeciesConcentrations; ChemicalSystemType m_chemicalSystemType; + + real64 m_solventDensity; }; // these aliases are useful in constitutive dispatch @@ -516,6 +551,8 @@ ENUM_STRINGS( ChemicalSystemType, "carbonate", "carbonateAllEquilibrium", "ultramafic", + "serpentinization", + "kineticCarbonate", "momasEasy", "momasMedium", "chainSerialAllKinetic" ); diff --git a/src/coreComponents/constitutive/fluid/singlefluid/CompressibleSinglePhaseFluid.hpp b/src/coreComponents/constitutive/fluid/singlefluid/CompressibleSinglePhaseFluid.hpp index e177053795c..8b4f0d3c02f 100644 --- a/src/coreComponents/constitutive/fluid/singlefluid/CompressibleSinglePhaseFluid.hpp +++ b/src/coreComponents/constitutive/fluid/singlefluid/CompressibleSinglePhaseFluid.hpp @@ -193,6 +193,9 @@ class CompressibleSinglePhaseFluid : public SingleFluidBase real64 defaultDensity() const override final { return m_defaultDensity; } real64 defaultViscosity() const override final { return m_defaultViscosity; } + /// @return the (constant) fluid compressibility [1/Pa] + real64 compressibility() const { return m_compressibility; } + protected: virtual void postInputInitialization() override; diff --git a/src/coreComponents/constitutive/permeability/CarmanKozenyPermeability.hpp b/src/coreComponents/constitutive/permeability/CarmanKozenyPermeability.hpp index 2a121e0d36c..48c1b18f4ec 100644 --- a/src/coreComponents/constitutive/permeability/CarmanKozenyPermeability.hpp +++ b/src/coreComponents/constitutive/permeability/CarmanKozenyPermeability.hpp @@ -151,9 +151,11 @@ void CarmanKozenyPermeabilityUpdate::compute( real64 const & porosity, { real64 const constant = pow( m_sphericity*m_particleDiameter, 2 ) / 150; - real64 const permValue = constant * pow( porosity, 3 )/ pow( (1 - porosity), 2 ); + real64 const eps = 1e-20; + + real64 const permValue = constant * pow( porosity, 3 )/ pow( (1 - porosity + eps), 2 ); - real64 const dPerm_dPorValue = -constant * ( (porosity - 3) * pow( porosity, 2 ) / pow( (1-porosity), 3 ) ); + real64 const dPerm_dPorValue = -constant * ( (porosity - 3) * pow( porosity, 2 ) / pow( (1-porosity+eps), 3 ) ); for( localIndex i = 0; i < permeability.size(); ++i ) { diff --git a/src/coreComponents/constitutive/solid/CoupledSolid.hpp b/src/coreComponents/constitutive/solid/CoupledSolid.hpp index da74fb42ce8..b01a9a1f0c8 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 { @@ -105,6 +106,42 @@ 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 & temperature, + arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & kineticReactionMolarIncrements ) const + { + GEOS_UNUSED_VAR( k, q, pressure, 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/EigenstrainReactiveSolid.cpp b/src/coreComponents/constitutive/solid/EigenstrainReactiveSolid.cpp new file mode 100644 index 00000000000..c658f6c2942 --- /dev/null +++ b/src/coreComponents/constitutive/solid/EigenstrainReactiveSolid.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 EigenstrainReactiveSolid.cpp + */ + +#include "EigenstrainReactiveSolid.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 > +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; + + +REGISTER_CATALOG_ENTRY( ConstitutiveBase, EigenStrainReactiveElasticIsotropicConstant, string const &, Group * const ) +REGISTER_CATALOG_ENTRY( ConstitutiveBase, EigenStrainReactiveElasticIsotropicCK, 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..b2e2c3f75e9 --- /dev/null +++ b/src/coreComponents/constitutive/solid/EigenstrainReactiveSolid.hpp @@ -0,0 +1,326 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 ); + } + + GEOS_HOST_DEVICE + virtual void updateStateFromPressureTemperatureAndReactions( localIndex const k, + localIndex const q, + real64 const & pressure, + 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 ); + } + + 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( 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..64454643642 --- /dev/null +++ b/src/coreComponents/constitutive/solid/PorousReactiveSolid.cpp @@ -0,0 +1,82 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 ) +{ + this->registerWrapper( "fluidModelName", &m_fluidModelName ). + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "Name of the fluid constitutive model. When set, its (constant) compressibility " + "is handed to the porosity model for the pore-mineral-pressure / porosity coupling." ); +} + +template< typename SOLID_TYPE, + typename PERM_TYPE > +void PorousReactiveSolid< SOLID_TYPE, PERM_TYPE >::initializeState() const +{ + CoupledSolid< SOLID_TYPE, BiotReactivePorosity, PERM_TYPE >::initializeState(); +} + +template< typename SOLID_TYPE, + typename PERM_TYPE > +void PorousReactiveSolid< SOLID_TYPE, PERM_TYPE >::initializePreSubGroups() +{ + CoupledSolid< SOLID_TYPE, BiotReactivePorosity, PERM_TYPE >::initializePreSubGroups(); + + // If a fluid model is specified, read its (constant) compressibility and hand it to the porosity + // model. The porosity model needs the fluid bulk modulus for the pore-mineral-pressure / porosity + // coupling but has no direct handle on the fluid. The fluid is looked up through its + // CompressibleSinglePhaseFluid base, which also covers the reactive and thermal variants. + if( !m_fluidModelName.empty() ) + { + CompressibleSinglePhaseFluid const & fluid = + this->getParent().template getGroup< CompressibleSinglePhaseFluid >( m_fluidModelName ); + BiotReactivePorosity & porosity = + dynamicCast< BiotReactivePorosity & >( this->getBasePorosityModel() ); + porosity.setFluidCompressibility( fluid.compressibility() ); + } +} + +// 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..1f6d390b953 --- /dev/null +++ b/src/coreComponents/constitutive/solid/PorousReactiveSolid.hpp @@ -0,0 +1,354 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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" +#include "constitutive/fluid/singlefluid/CompressibleSinglePhaseFluid.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 & 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 ); + } + + 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; + + virtual void initializePreSubGroups() override; + + /** + * @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; + + /// Name of the fluid model used to read the (constant) fluid compressibility for the + /// pore-mineral-pressure / porosity coupling (optional; empty => fluid treated as incompressible) + string m_fluidModelName; +}; + + + +} +} /* namespace geos */ + +#endif /* GEOS_CONSTITUTIVE_SOLID_POROUSREACTIVESOLID_HPP_ */ 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..f22cb568b28 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,24 @@ 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 & 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 ); } 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 +87,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..53d7c1938a4 --- /dev/null +++ b/src/coreComponents/constitutive/solid/porosity/BiotReactivePorosity.hpp @@ -0,0 +1,359 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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, + real64 const fluidCompressibility, + arrayView2d< real64 > const & meanEffectiveStressIncrement_k, + integer const fixedPorosity ): ReactivePorosityBaseUpdates( newPorosity, + porosity_n, + dPorosity_dPressure, + dPorosity_dTemperature, + initialPorosity, + referencePorosity, + volumeFractions, + initialVolumeFractions, + volumeFractions_n, + numKineticReactions, + molarWeights, + mineralDensities, + bulkModulus, + meanEffectiveStressIncrement_k, + fixedPorosity ), + m_grainBulkModulus( grainBulkModulus ), + m_biotCoefficient( biotCoefficient ), + m_mineralBulkModulus( mineralBulkModulus ), + m_mineralPressure( mineralPressure ), + m_mineralPressure_n( mineralPressure_n ), + m_fluidCompressibility( fluidCompressibility ) + {} + + 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]; } + + /// @return the fluid bulk modulus K_f = 1/compressibility [Pa] (0 if incompressible) + GEOS_HOST_DEVICE + real64 getFluidBulkModulus( localIndex const k ) const + { + GEOS_UNUSED_VAR( k ); + return m_fluidCompressibility > 0.0 ? 1.0 / m_fluidCompressibility : 0.0; + } + + 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 & fluidBulkModulus, + real64 const & reactionPorosityIncrement ) const + { + GEOS_UNUSED_VAR( pressure_k ); + + real64 const biotSkeletonModulusInverse = (biotCoefficient - referencePorosity) / grainBulkModulus; + real64 const porosityMultiplierInverse = 1 / ( 1 + biotSkeletonModulusInverse*mineralBulkModulus/referencePorosity ); + real64 const poreFluidMineralBulkModRatio = ( fluidBulkModulus + mineralBulkModulus ) / fluidBulkModulus; + + real64 const anelasticStrainIncrement = -reactionPorosityIncrement; + + porosity = porosity_n + // change due to stress increment + + biotCoefficient * meanEffectiveStressIncrement_k / bulkModulus * porosityMultiplierInverse + // change due to pressure increment + + biotSkeletonModulusInverse * poreFluidMineralBulkModRatio * ( pressure - pressure_n ) * porosityMultiplierInverse + // change due to mineral pressure increment + + biotSkeletonModulusInverse * 3 * anelasticStrainIncrement * mineralBulkModulus * porosityMultiplierInverse + // change due to mineral volume fraction increment + + reactionPorosityIncrement; + + dPorosity_dPressure = biotSkeletonModulusInverse * poreFluidMineralBulkModRatio * 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 & fluidBulkModulus, + 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 * 3 + // change due to stress increment + - mineralPressureMultiplier * biotCoefficient * meanEffectiveStressIncrement / bulkModulus + // change due to pressure increment + - mineralPressureMultiplier * biotSkeletonModulusInverse * deltaPressureFromLastStep + + mineralPressureMultiplier * referencePorosity / fluidBulkModulus * 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], + getFluidBulkModulus( 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 ); + if( !m_fixedPorosity ) + { + 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], + getFluidBulkModulus( 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; + + /// Fluid compressibility [1/Pa] (constant, pushed in from the flow solver's fluid model) + real64 const m_fluidCompressibility; +}; + +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_fluidCompressibility, + m_meanEffectiveStressIncrement_k, + m_fixedPorosity ); + } + + /** + * @brief Store the (constant) fluid compressibility read from the flow solver's fluid model. + * @param fluidCompressibility the fluid compressibility [1/Pa] + */ + void setFluidCompressibility( real64 const fluidCompressibility ) + { m_fluidCompressibility = fluidCompressibility; } + +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; + + /// Fluid compressibility [1/Pa], pushed in from the flow solver's fluid model (0 until set) + real64 m_fluidCompressibility = 0.0; +}; + +} /* 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 70% rename from src/coreComponents/constitutive/solid/porosity/ReactivePorosity.cpp rename to src/coreComponents/constitutive/solid/porosity/ReactivePorosityBase.cpp index 7865764fd7c..eada144ca67 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,34 @@ 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" ); + + registerWrapper( viewKeyStruct::fixedPorosityString(), &m_fixedPorosity ). + setApplyDefaultValue( 0 ). + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "Flag to update volume fractions without changing porosity (0 = update porosity, 1 = fixed porosity)" ); + } -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 +100,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 +109,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 +163,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..391470de66e --- /dev/null +++ b/src/coreComponents/constitutive/solid/porosity/ReactivePorosityBase.hpp @@ -0,0 +1,315 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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, + integer const fixedPorosity ): + 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 ), + m_fixedPorosity( fixedPorosity ) + {} + + 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 ); + + if( !m_fixedPorosity ) + { + 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 + if( !m_fixedPorosity ) + { + m_newPorosity[k][q] = m_porosity_n[k][q] + + reactionPorosityIncrement; + // + m_meanEffectiveStressIncrement_k[k][q]/m_bulkModulus[k]; + + 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; + } + } + } + + 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; + + integer const m_fixedPorosity; +}; + + +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"; } + static constexpr char const * fixedPorosityString() { return "fixedPorosity"; } + } 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, + m_fixedPorosity ); + } + + +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; + + integer m_fixedPorosity; +}; + + +}/* namespace constitutive */ + +} /* namespace geos */ + + +#endif //GEOS_CONSTITUTIVE_POROSITY_REACTIVEPOROSITYBASE_HPP_ diff --git a/src/coreComponents/dataRepository/unitTests/testErrorHandling.cpp b/src/coreComponents/dataRepository/unitTests/testErrorHandling.cpp index d2a2d16257a..c53161c01f0 100644 --- a/src/coreComponents/dataRepository/unitTests/testErrorHandling.cpp +++ b/src/coreComponents/dataRepository/unitTests/testErrorHandling.cpp @@ -28,6 +28,8 @@ #include #include +#include + using namespace geos; using namespace dataRepository; diff --git a/src/coreComponents/integrationTests/fluidFlowTests/testSinglePhaseReactiveTransport.cpp b/src/coreComponents/integrationTests/fluidFlowTests/testSinglePhaseReactiveTransport.cpp index c4b50dfab30..20e3980fe89 100644 --- a/src/coreComponents/integrationTests/fluidFlowTests/testSinglePhaseReactiveTransport.cpp +++ b/src/coreComponents/integrationTests/fluidFlowTests/testSinglePhaseReactiveTransport.cpp @@ -530,7 +530,7 @@ class SinglePhaseReactiveTransportTest : public ::testing::Test } static real64 constexpr time = 0.0; - static real64 constexpr dt = 1.0; + static real64 constexpr dt = 0.01; static real64 constexpr eps = std::numeric_limits< real64 >::epsilon(); GeosxState state; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp index b5c850662d3..1f63bb1b163 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp @@ -47,18 +47,47 @@ 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 & 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], + 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 +263,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 {}", @@ -463,7 +492,7 @@ void SinglePhaseReactiveTransport::assembleFluxTerms( real64 const dt, forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, MeshLevel const & mesh, - string_array const & ) + string_array const & regionNames ) { NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager(); FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager(); @@ -471,6 +500,25 @@ void SinglePhaseReactiveTransport::assembleFluxTerms( real64 const dt, string const & dofKey = dofManager.getKey( viewKeyStruct::elemDofFieldString() ); + real64 solventDensity = 1.0; + mesh.getElemManager().forElementSubRegions( regionNames, + [&]( localIndex const, + ElementSubRegionBase const & subRegion ) + { + if( m_isThermal ) + { + reactivefluid::ReactiveThermalCompressibleSinglePhaseFluid const & fluid = + getConstitutiveModel< reactivefluid::ReactiveThermalCompressibleSinglePhaseFluid >( subRegion, subRegion.template getReference< string >( viewKeyStruct::fluidNamesString() ) ); + solventDensity = fluid.solventDensity(); + } + else + { + reactivefluid::ReactiveCompressibleSinglePhaseFluid const & fluid = + getConstitutiveModel< reactivefluid::ReactiveCompressibleSinglePhaseFluid >( subRegion, subRegion.template getReference< string >( viewKeyStruct::fluidNamesString() ) ); + solventDensity = fluid.solventDensity(); + } + } ); + fluxApprox.forAllStencils( mesh, [&] ( auto & stencil ) { typename TYPEOFREF( stencil ) ::KernelWrapper stencilWrapper = stencil.createKernelWrapper(); @@ -481,6 +529,7 @@ void SinglePhaseReactiveTransport::assembleFluxTerms( real64 const dt, FluxComputeKernelFactory::createAndLaunch< parallelDevicePolicy<> >( m_numPrimarySpecies, m_hasDiffusion, mobilePrimarySpeciesFlags.toViewConst(), + solventDensity, dofManager.rankOffset(), dofKey, getName(), @@ -496,6 +545,7 @@ void SinglePhaseReactiveTransport::assembleFluxTerms( real64 const dt, FluxComputeKernelFactory::createAndLaunch< parallelDevicePolicy<> >( m_numPrimarySpecies, m_hasDiffusion, mobilePrimarySpeciesFlags.toViewConst(), + solventDensity, dofManager.rankOffset(), dofKey, getName(), @@ -552,15 +602,16 @@ void SinglePhaseReactiveTransport::updateSpeciesAmount( ElementSubRegionBase & s getConstitutiveModel< reactivefluid::ReactiveThermalCompressibleSinglePhaseFluid >( subRegion, subRegion.getReference< string >( viewKeyStruct::fluidNamesString() ) ); arrayView3d< real64 const, reactivefluid::USD_SPECIES > const primarySpeciesAggregateConcentration = fluid.primarySpeciesAggregateConcentration(); arrayView3d< real64 const, reactivefluid::USD_SPECIES > const primarySpeciesAggregateConcentration_n = fluid.primarySpeciesAggregateConcentration_n(); + real64 const solventDensity = fluid.solventDensity(); forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei ) { for( integer is = 0; is < numPrimarySpecies; ++is ) { - primarySpeciesAggregateMole[ei][is] = porosity[ei][0] * ( volume[ei] + deltaVolume[ei] ) * primarySpeciesAggregateConcentration[ei][0][is]; + primarySpeciesAggregateMole[ei][is] = porosity[ei][0] * ( volume[ei] + deltaVolume[ei] ) * primarySpeciesAggregateConcentration[ei][0][is] * solventDensity; if( isZero( primarySpeciesAggregateMole_n[ei][is] ) ) - primarySpeciesAggregateMole_n[ei][is] = porosity_n[ei][0] * volume[ei] * primarySpeciesAggregateConcentration_n[ei][0][is]; + primarySpeciesAggregateMole_n[ei][is] = porosity_n[ei][0] * volume[ei] * primarySpeciesAggregateConcentration_n[ei][0][is] * solventDensity; } } ); } @@ -570,15 +621,16 @@ void SinglePhaseReactiveTransport::updateSpeciesAmount( ElementSubRegionBase & s getConstitutiveModel< reactivefluid::ReactiveCompressibleSinglePhaseFluid >( subRegion, subRegion.getReference< string >( viewKeyStruct::fluidNamesString() ) ); arrayView3d< real64 const, reactivefluid::USD_SPECIES > const primarySpeciesAggregateConcentration = fluid.primarySpeciesAggregateConcentration(); arrayView3d< real64 const, reactivefluid::USD_SPECIES > const primarySpeciesAggregateConcentration_n = fluid.primarySpeciesAggregateConcentration_n(); + real64 const solventDensity = fluid.solventDensity(); forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei ) { for( integer is = 0; is < numPrimarySpecies; ++is ) { - primarySpeciesAggregateMole[ei][is] = porosity[ei][0] * ( volume[ei] + deltaVolume[ei] ) * primarySpeciesAggregateConcentration[ei][0][is]; + primarySpeciesAggregateMole[ei][is] = porosity[ei][0] * ( volume[ei] + deltaVolume[ei] ) * primarySpeciesAggregateConcentration[ei][0][is] * solventDensity; if( isZero( primarySpeciesAggregateMole_n[ei][is] ) ) - primarySpeciesAggregateMole_n[ei][is] = porosity_n[ei][0] * volume[ei] * primarySpeciesAggregateConcentration_n[ei][0][is]; + primarySpeciesAggregateMole_n[ei][is] = porosity_n[ei][0] * volume[ei] * primarySpeciesAggregateConcentration_n[ei][0][is] * solventDensity; } } ); } @@ -598,12 +650,13 @@ void SinglePhaseReactiveTransport::updateKineticReactionMolarIncrements( real64 reactivefluid::ReactiveThermalCompressibleSinglePhaseFluid & fluid = getConstitutiveModel< reactivefluid::ReactiveThermalCompressibleSinglePhaseFluid >( subRegion, subRegion.getReference< string >( viewKeyStruct::fluidNamesString() ) ); arrayView3d< real64 const, reactivefluid::USD_SPECIES > const kineticReactionRates = fluid.kineticReactionRates(); + real64 const solventDensity = fluid.solventDensity(); forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei ) { for( integer r = 0; r < numKineticReactions; ++r ) { - kineticReactionMolarIncrements[ei][r] = dt* kineticReactionRates[ei][0][r]; + kineticReactionMolarIncrements[ei][r] = dt * kineticReactionRates[ei][0][r] * solventDensity; } } ); } @@ -612,12 +665,13 @@ void SinglePhaseReactiveTransport::updateKineticReactionMolarIncrements( real64 reactivefluid::ReactiveCompressibleSinglePhaseFluid & fluid = getConstitutiveModel< reactivefluid::ReactiveCompressibleSinglePhaseFluid >( subRegion, subRegion.getReference< string >( viewKeyStruct::fluidNamesString() ) ); arrayView3d< real64 const, reactivefluid::USD_SPECIES > const kineticReactionRates = fluid.kineticReactionRates(); + real64 const solventDensity = fluid.solventDensity(); forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei ) { for( integer r = 0; r < numKineticReactions; ++r ) { - kineticReactionMolarIncrements[ei][r] = dt* kineticReactionRates[ei][0][r]; + kineticReactionMolarIncrements[ei][r] = dt * kineticReactionRates[ei][0][r] * solventDensity; } } ); } @@ -662,15 +716,27 @@ void SinglePhaseReactiveTransport::updatePorosityAndPermeability( CellElementSub if( m_isUpdateReactivePorosity ) { arrayView1d< real64 const > const & pressure = subRegion.getField< fields::flow::pressure >(); + 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_n = subRegion.getField< fields::flow::pressure_n >(); + 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, temperature, kineticReactionMolarIncrements ); + } } ); } else @@ -679,6 +745,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 +798,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/SinglePhaseReactiveTransportFields.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransportFields.hpp index 28c72f5c170..ea1684ffd4b 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransportFields.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransportFields.hpp @@ -41,29 +41,29 @@ DECLARE_FIELD( logPrimarySpeciesConcentration, 0, LEVEL_0, WRITE_AND_READ, - "Natural log of primary species concentration (molarity)" ); + "Natural log of primary species concentration (molality)" ); DECLARE_FIELD( logPrimarySpeciesConcentration_n, "logPrimarySpeciesConcentration_n", array2dLayoutComp, 0, - LEVEL_0, + NOPLOT, WRITE_AND_READ, - "Natural log of primary species concentration (molarity) at the previous converged time step" ); + "Natural log of primary species concentration (molality) at the previous converged time step" ); DECLARE_FIELD( bcLogPrimarySpeciesConcentration, "bcLogPrimarySpeciesConcentration", array2dLayoutComp, 0, - LEVEL_0, + NOPLOT, WRITE_AND_READ, - "Boundary condition for natural log of primary species concentration (molarity)" ); + "Boundary condition for natural log of primary species concentration (molality)" ); DECLARE_FIELD( primarySpeciesAggregateMole, "primarySpeciesAggregateMole", array2dLayoutComp, 0, - LEVEL_0, + NOPLOT, WRITE_AND_READ, "Aggregate amount of primary species in mole" ); @@ -71,7 +71,7 @@ DECLARE_FIELD( primarySpeciesAggregateMole_n, "primarySpeciesAggregateMole_n", array2dLayoutComp, 0, - LEVEL_0, + NOPLOT, WRITE_AND_READ, "Aggregate amount of primary species in mole at the previous converged time step" ); @@ -95,7 +95,7 @@ DECLARE_FIELD( surfaceArea, "surfaceArea", array2dLayoutComp, 0, - NOPLOT, + LEVEL_0, WRITE_AND_READ, "Surface area for surface reactions." ); @@ -103,7 +103,7 @@ DECLARE_FIELD( initialSurfaceArea, "initialSurfaceArea", array2dLayoutComp, 0, - NOPLOT, + LEVEL_0, WRITE_AND_READ, "Initial surface area for surface reactions." ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/AccumulationKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/AccumulationKernels.hpp index 28c3b1cc181..e3e3dee58c8 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/AccumulationKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/AccumulationKernels.hpp @@ -82,6 +82,7 @@ class AccumulationKernel : public singlePhaseBaseKernels::AccumulationKernel< SU arrayView1d< real64 > const & localRhs ) : Base( rankOffset, dofKey, subRegion, localMatrix, localRhs ), m_dt( dt ), + m_solventDensity( fluid.solventDensity() ), m_volume( subRegion.getElementVolume() ), m_deltaVolume( subRegion.template getField< fields::flow::deltaVolume >() ), m_porosity( solid.getPorosity() ), @@ -160,8 +161,8 @@ class AccumulationKernel : public singlePhaseBaseKernels::AccumulationKernel< SU void computeAccumulation( localIndex const ei, StackVariables & stack ) const { - // Residual[is] += (primarySpeciesAggregateConcentration[is] * stack.poreVolume - primarySpeciesAggregateMole_n[is]) - // - dt * m_volume * primarySpeciesKineticRate[is] // To Check: what's the unit of the kinetic rate + // Residual[is] += (primarySpeciesAggregateConcentration[is] * solventDensity * stack.poreVolume - primarySpeciesAggregateMole_n[is]) + // - dt * m_volume * primarySpeciesKineticRate[is] * solventDensity Base::computeAccumulation( ei, stack ); @@ -183,14 +184,14 @@ class AccumulationKernel : public singlePhaseBaseKernels::AccumulationKernel< SU // Step 2.1: residual // Primary species mole amount in pore volume stack.localResidual[is+numEqn-numSpecies] -= m_primarySpeciesAggregateMole_n[ei][is]; - stack.localResidual[is+numEqn-numSpecies] += m_primarySpeciesAggregateConcentration[ei][0][is] * stack.poreVolume; + stack.localResidual[is+numEqn-numSpecies] += m_primarySpeciesAggregateConcentration[ei][0][is] * m_solventDensity * stack.poreVolume; // Reaction term - stack.localResidual[is+numEqn-numSpecies] -= m_dt * ( m_volume[ei] + m_deltaVolume[ei] ) * m_primarySpeciesAggregateKineticRate[ei][0][is]; + stack.localResidual[is+numEqn-numSpecies] -= m_dt * ( m_volume[ei] + m_deltaVolume[ei] ) * m_primarySpeciesAggregateKineticRate[ei][0][is] * m_solventDensity; // Step 2.1: jacobian // Drivative of primary species amount in pore volume wrt pressure - stack.localJacobian[is+numEqn-numSpecies][0] += stack.dPoreVolume_dPres * m_primarySpeciesAggregateConcentration[ei][0][is] + stack.localJacobian[is+numEqn-numSpecies][0] += stack.dPoreVolume_dPres * m_primarySpeciesAggregateConcentration[ei][0][is] * m_solventDensity /* + stack.poreVolume * m_dTotalPrimarySpeciesConcentration_dPres[ei][is] */; // // Derivative of reaction term wrt pressure // stack.localJacobian[is+numEqn-numSpecies][0] -= m_dt * ( m_volume[ei] + m_deltaVolume[ei] ) * @@ -201,15 +202,9 @@ class AccumulationKernel : public singlePhaseBaseKernels::AccumulationKernel< SU { stack.localJacobian[is+numEqn-numSpecies][js+numDof-numSpecies] = /* stack.dPoreVolume_dLogPrimaryConc[js] * m_primarySpeciesAggregateConcentration[ei][0][is] - + */stack.poreVolume * dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations[is][js]; // To - // check - // if - // the - // permutation - // is - // consistent - - stack.localJacobian[is+numEqn-numSpecies][js+numDof-numSpecies] -= m_dt * ( m_volume[ei] + m_deltaVolume[ei] ) * dPrimarySpeciesAggregateKineticRate_dLogPrimaryConc[is][js]; + + */stack.poreVolume * dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations[is][js] * m_solventDensity; + + stack.localJacobian[is+numEqn-numSpecies][js+numDof-numSpecies] -= m_dt * ( m_volume[ei] + m_deltaVolume[ei] ) * dPrimarySpeciesAggregateKineticRate_dLogPrimaryConc[is][js] * m_solventDensity; } } } @@ -245,6 +240,9 @@ class AccumulationKernel : public singlePhaseBaseKernels::AccumulationKernel< SU /// Time step size real64 const m_dt; + /// Solvent density [kg/m³] used to convert molality [mol/kg] to molarity [mol/m³] + real64 const m_solventDensity; + /// View on the element volumes arrayView1d< real64 const > const m_volume; arrayView1d< real64 const > const m_deltaVolume; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/DirichletFluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/DirichletFluxComputeKernel.hpp index 80ff9bce9ba..fac64cc70da 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/DirichletFluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/DirichletFluxComputeKernel.hpp @@ -117,6 +117,7 @@ class DirichletFluxComputeKernel : public singlePhaseFVMKernels::DirichletFluxCo ReactiveSinglePhaseFluidAccessors const & reactiveSinglePhaseFluidAccessors, PermeabilityAccessors const & permeabilityAccessors, arrayView1d< integer const > const & mobilePrimarySpeciesFlags, + real64 const & solventDensity, real64 const & dt, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) @@ -135,7 +136,8 @@ class DirichletFluxComputeKernel : public singlePhaseFVMKernels::DirichletFluxCo m_primarySpeciesMobileAggregateConc( reactiveSinglePhaseFluidAccessors.get( fields::reactivefluid::primarySpeciesMobileAggregateConcentration {} ) ), m_dPrimarySpeciesMobileAggregateConc_dLogPrimaryConc( reactiveSinglePhaseFluidAccessors.get( fields::reactivefluid::dPrimarySpeciesMobileAggregateConcentration_dLogPrimarySpeciesConcentrations {} ) ), - m_mobilePrimarySpeciesFlags( mobilePrimarySpeciesFlags ) + m_mobilePrimarySpeciesFlags( mobilePrimarySpeciesFlags ), + m_solventDensity( solventDensity ) {} /** @@ -194,17 +196,17 @@ class DirichletFluxComputeKernel : public singlePhaseFVMKernels::DirichletFluxCo for( integer is = 0; is < numSpecies; ++is ) { - real64 const aggregateConc_i = m_primarySpeciesMobileAggregateConc[seri][sesri][sei][0][is]; - speciesFlux[is] = aggregateConc_i / dens_up * fluxVal * mobility_up; + real64 const aggregateConcMolarity_i = m_primarySpeciesMobileAggregateConc[seri][sesri][sei][0][is] * m_solventDensity; + speciesFlux[is] = aggregateConcMolarity_i / dens_up * fluxVal * mobility_up; - dSpeciesFlux_dP[is] = aggregateConc_i / dens_up * dFlux_dP * mobility_up - + aggregateConc_i / dens_up * fluxVal * dMobility_dP_up - - aggregateConc_i * fluxVal * mobility_up * dDens_dP_up / (dens_up * dens_up); + dSpeciesFlux_dP[is] = aggregateConcMolarity_i / dens_up * dFlux_dP * mobility_up + + aggregateConcMolarity_i / dens_up * fluxVal * dMobility_dP_up + - aggregateConcMolarity_i * fluxVal * mobility_up * dDens_dP_up / (dens_up * dens_up); for( integer js = 0; js < numSpecies; ++js ) { - real64 const dAggregateConc_i_dLogConc_j = m_dPrimarySpeciesMobileAggregateConc_dLogPrimaryConc[seri][sesri][sei][0][is][js]; - dSpeciesFlux_dLogConc[is][js] += dAggregateConc_i_dLogConc_j / dens_up * fluxVal * mobility_up; + real64 const dAggregateConcMolarity_i_dLogConc_j = m_dPrimarySpeciesMobileAggregateConc_dLogPrimaryConc[seri][sesri][sei][0][is][js] * m_solventDensity; + dSpeciesFlux_dLogConc[is][js] += dAggregateConcMolarity_i_dLogConc_j / dens_up * fluxVal * mobility_up; } } @@ -267,6 +269,9 @@ class DirichletFluxComputeKernel : public singlePhaseFVMKernels::DirichletFluxCo /// Array of flags to indicate mobile primary species arrayView1d< integer const > const m_mobilePrimarySpeciesFlags; + /// Solvent density [kg/m³] used to convert molality [mol/kg] to molarity [mol/m³] + real64 const m_solventDensity; + }; @@ -344,6 +349,7 @@ class DirichletFluxComputeKernelFactory reactiveFluidAccessors, permeabilityAccessors, mobilePrimarySpeciesFlags, + reactiveFluid.solventDensity(), dt, localMatrix, localRhs ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/FluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/FluxComputeKernel.hpp index a408c29aa0d..2fdcc47d656 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/FluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/FluxComputeKernel.hpp @@ -122,6 +122,7 @@ class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_E * @param[in] porosityAccessors * @param[in] hasDiffusion the flag to turn on diffusion calculation * @param[in] mobilePrimarySpeciesFlags the array of flags to indicate mobile primary species + * @param[in] solventDensity the density of the solvent (e.g., water) [kg/m3] * @param[in] dt time step size * @param[inout] localMatrix the local CRS matrix * @param[inout] localRhs the local right-hand side vector @@ -138,6 +139,7 @@ class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_E PorosityAccessors const & porosityAccessors, integer const & hasDiffusion, arrayView1d< integer const > const & mobilePrimarySpeciesFlags, + real64 const & solventDensity, real64 const & dt, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) @@ -159,7 +161,8 @@ class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_E m_dDiffusivity_dTemp( diffusionAccessors.get( fields::diffusion::dDiffusivity_dTemperature {} ) ), m_referencePorosity( porosityAccessors.get( fields::porosity::referencePorosity {} ) ), m_hasDiffusion( hasDiffusion ), - m_mobilePrimarySpeciesFlags( mobilePrimarySpeciesFlags ) + m_mobilePrimarySpeciesFlags( mobilePrimarySpeciesFlags ), + m_solventDensity( solventDensity ) {} /** @@ -249,20 +252,22 @@ class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_E // compute species fluxes and derivatives using upstream cell concentration for( integer is = 0; is < numSpecies; ++is ) { - real64 const aggregateConc_i = m_primarySpeciesMobileAggregateConc[er_up][esr_up][ei_up][0][is]; - speciesFlux[is] = aggregateConc_i / fluidDens_up * fluxVal; + real64 const aggregateConcMolarity_i = m_primarySpeciesMobileAggregateConc[er_up][esr_up][ei_up][0][is] + * m_solventDensity; // convert from mol/kg to mol/m3 using solvent density + speciesFlux[is] = aggregateConcMolarity_i / fluidDens_up * fluxVal; for( integer ke = 0; ke < numFluxSupportPoints; ++ke ) { - dSpeciesFlux_dP[ke][is] += aggregateConc_i / fluidDens_up * dFlux_dP[ke]; + dSpeciesFlux_dP[ke][is] += aggregateConcMolarity_i / fluidDens_up * dFlux_dP[ke]; } - dSpeciesFlux_dP[k_up][is] += -aggregateConc_i * fluxVal * dDens_dPres / (fluidDens_up * fluidDens_up); + dSpeciesFlux_dP[k_up][is] += -aggregateConcMolarity_i * fluxVal * dDens_dPres / (fluidDens_up * fluidDens_up); for( integer js = 0; js < numSpecies; ++js ) { - real64 const dAggregateConc_i_dLogConc_j = m_dPrimarySpeciesMobileAggregateConc_dLogPrimaryConc[er_up][esr_up][ei_up][0][is][js]; - dSpeciesFlux_dLogConc[k_up][is][js] += dAggregateConc_i_dLogConc_j / fluidDens_up * fluxVal; + real64 const dAggregateConcMolarity_i_dLogConc_j = m_dPrimarySpeciesMobileAggregateConc_dLogPrimaryConc[er_up][esr_up][ei_up][0][is][js] + * m_solventDensity; // convert from mol/kg to mol/m3 using solvent density + dSpeciesFlux_dLogConc[k_up][is][js] += dAggregateConcMolarity_i_dLogConc_j / fluidDens_up * fluxVal; } } @@ -354,15 +359,16 @@ class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_E localIndex const esr = sesri[ke]; localIndex const ei = sei[ke]; - real64 const aggregateConc_i = m_primarySpeciesMobileAggregateConc[er][esr][ei][0][is]; + real64 const aggregateConcMolarity_i = m_primarySpeciesMobileAggregateConc[er][esr][ei][0][is] + * m_solventDensity; // convert from mol/kg to mol/m3 using solvent density - speciesGrad[is] += diffusionTrans[ke] * aggregateConc_i; + speciesGrad[is] += diffusionTrans[ke] * aggregateConcMolarity_i; for( integer js = 0; js < numSpecies; ++js ) { - real64 const dAggregateConc_i_dLogConc_j = m_dPrimarySpeciesMobileAggregateConc_dLogPrimaryConc[er][esr][ei][0][is][js]; + real64 const dAggregateConcMolarity_i_dLogConc_j = m_dPrimarySpeciesMobileAggregateConc_dLogPrimaryConc[er][esr][ei][0][is][js] * m_solventDensity; - dSpeciesGrad_i_dLogConc[ke][js] += diffusionTrans[ke] * dAggregateConc_i_dLogConc_j; + dSpeciesGrad_i_dLogConc[ke][js] += diffusionTrans[ke] * dAggregateConcMolarity_i_dLogConc_j; } } @@ -500,6 +506,9 @@ class FluxComputeKernel : public singlePhaseFVMKernels::FluxComputeKernel< NUM_E /// Array of flags to indicate mobile primary species arrayView1d< integer const > const m_mobilePrimarySpeciesFlags; + + /// Density of the solvent (e.g., water) [kg/m3] + real64 const m_solventDensity; }; /** @@ -516,6 +525,7 @@ class FluxComputeKernelFactory * @param[in] numSpecies the number of primary species * @param[in] hasDiffusion the flag of adding diffusion term * @param[in] mobilePrimarySpeciesFlags the array of flags to indicate mobile primary species + * @param[in] solventDensity the density of the solvent (e.g., water) [kg/m3] * @param[in] rankOffset the offset of my MPI rank * @param[in] dofKey string to get the element degrees of freedom numbers * @param[in] solverName name of the solver (to name accessors) @@ -530,6 +540,7 @@ class FluxComputeKernelFactory createAndLaunch( integer const numSpecies, integer const hasDiffusion, arrayView1d< integer const > const mobilePrimarySpeciesFlags, + real64 const solventDensity, globalIndex const rankOffset, string const & dofKey, string const & solverName, @@ -561,7 +572,7 @@ class FluxComputeKernelFactory KernelType kernel( rankOffset, stencilWrapper, dofNumberAccessor, flowAccessors, reactiveFlowAccessors, fluidAccessors, reactiveFluidAccessors, permAccessors, diffusionAccessors, porosityAccessors, hasDiffusion, mobilePrimarySpeciesFlags, - dt, localMatrix, localRhs ); + solventDensity, dt, localMatrix, localRhs ); KernelType::template launch< POLICY >( stencilWrapper.size(), kernel ); } ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/SourceFluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/SourceFluxComputeKernel.hpp index 1ef58f942aa..8f2897a451f 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/SourceFluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/SourceFluxComputeKernel.hpp @@ -74,6 +74,7 @@ class SourceFluxComputeKernel m_elemGhostRank( elemGhostRank ), m_rhsContributionArrayView( rhsContributionArrayView ), m_sizeScalingFactor( sizeScalingFactor ), + m_solventDensity( fluid.solventDensity() ), m_primarySpeciesAggregateConcentration( fluid.primarySpeciesAggregateConcentration() ), m_dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations( fluid.dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations() ), m_density( fluid.density() ), @@ -151,12 +152,12 @@ class SourceFluxComputeKernel for( integer i = 0; i < numSpecies; ++i ) { - stack.localSpeciesRhs[i] += m_primarySpeciesAggregateConcentration[ei][0][i] / m_density[ei][0] * scaledInflowMass; - stack.localSpeciesJacobian[i][0] += -m_primarySpeciesAggregateConcentration[ei][0][i] * m_dDensity[ei][0][DerivOffset::dP] / (m_density[ei][0] * m_density[ei][0]) * scaledInflowMass; + stack.localSpeciesRhs[i] += m_primarySpeciesAggregateConcentration[ei][0][i] * m_solventDensity / m_density[ei][0] * scaledInflowMass; + stack.localSpeciesJacobian[i][0] += -m_primarySpeciesAggregateConcentration[ei][0][i] * m_solventDensity * m_dDensity[ei][0][DerivOffset::dP] / (m_density[ei][0] * m_density[ei][0]) * scaledInflowMass; for( integer j = 0; j < numSpecies; ++j ) { - stack.localSpeciesJacobian[i][j+numDof-numSpecies] += m_dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations[ei][0][i][j] / m_density[ei][0] * scaledInflowMass; + stack.localSpeciesJacobian[i][j+numDof-numSpecies] += m_dPrimarySpeciesAggregateConcentration_dLogPrimarySpeciesConcentrations[ei][0][i][j] * m_solventDensity / m_density[ei][0] * scaledInflowMass; } } } @@ -237,6 +238,9 @@ class SourceFluxComputeKernel /// size scaling factor real64 const m_sizeScalingFactor; + /// Solvent density [kg/m³] used to convert molality [mol/kg] to molarity [mol/m³] + real64 const m_solventDensity; + // View on the total concentration of ions that contain the primary species arrayView3d< real64 const, constitutive::reactivefluid::USD_SPECIES > const m_primarySpeciesAggregateConcentration; // View on the derivatives of total ion concentration for the primary species wrt log of primary species concentration diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalAccumulationKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalAccumulationKernels.hpp index 33fa2ac0f3b..276ea4bea5d 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalAccumulationKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalAccumulationKernels.hpp @@ -53,6 +53,7 @@ class AccumulationKernel : public singlePhaseReactiveBaseKernels::AccumulationKe using Base::m_volume; using Base::m_deltaVolume; using Base::m_primarySpeciesAggregateConcentration; + using Base::m_solventDensity; /// Note: Derivative lineup only supports dP & dT, not component terms using DerivOffset = constitutive::singlefluid::DerivativeOffsetC< 1 >; @@ -161,9 +162,9 @@ class AccumulationKernel : public singlePhaseReactiveBaseKernels::AccumulationKe for( integer is = 0; is < numSpecies; ++is ) { // Drivative of primary species amount in pore volume wrt temperature - stack.localJacobian[is+numEqn-numSpecies][numDof-numSpecies-1] += stack.dPoreVolume_dTemp * m_primarySpeciesAggregateConcentration[ei][0][is] + stack.localJacobian[is+numEqn-numSpecies][numDof-numSpecies-1] += stack.dPoreVolume_dTemp * m_primarySpeciesAggregateConcentration[ei][0][is] * m_solventDensity /* + stack.poreVolume * - m_dPrimarySpeciesAggregateConcentration_dTemp[ei][is] */; + m_dPrimarySpeciesAggregateConcentration_dTemp[ei][is] * m_solventDensity */; // // Derivative of reaction term wrt temperature // stack.localJacobian[is+numEqn-numSpecies][numDof-numSpecies-1] -= m_dt * ( m_volume[ei] + m_deltaVolume[ei] ) * // m_dPrimarySpeciesTotalKineticRate_dTemp[is]; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalDirichletFluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalDirichletFluxComputeKernel.hpp index 7a87bdfbc50..73d7135862c 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalDirichletFluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalDirichletFluxComputeKernel.hpp @@ -86,6 +86,7 @@ class DirichletFluxComputeKernel : public singlePhaseReactiveFVMKernels::Dirichl using Base::m_sei; using Base::m_facePres; using Base::m_faceGravCoef; + using Base::m_solventDensity; using ReactiveSinglePhaseFlowAccessors = typename Base::ReactiveSinglePhaseFlowAccessors; using ReactiveSinglePhaseFluidAccessors = typename Base::ReactiveSinglePhaseFluidAccessors; @@ -137,6 +138,7 @@ class DirichletFluxComputeKernel : public singlePhaseReactiveFVMKernels::Dirichl PermeabilityAccessors const & permeabilityAccessors, ThermalConductivityAccessors const & thermalConductivityAccessors, arrayView1d< integer const > const & mobilePrimarySpeciesFlags, + real64 const & solventDensity, real64 const & dt, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) @@ -152,6 +154,7 @@ class DirichletFluxComputeKernel : public singlePhaseReactiveFVMKernels::Dirichl reactiveSinglePhaseFluidAccessors, permeabilityAccessors, mobilePrimarySpeciesFlags, + solventDensity, dt, localMatrix, localRhs ), @@ -392,6 +395,7 @@ class DirichletFluxComputeKernelFactory permeabilityAccessors, thermalConductivityAccessors, mobilePrimarySpeciesFlags, + reactiveFluid.solventDensity(), dt, localMatrix, localRhs ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalFluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalFluxComputeKernel.hpp index d436b261628..80c262efcd4 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalFluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalFluxComputeKernel.hpp @@ -89,6 +89,7 @@ class FluxComputeKernel : public singlePhaseReactiveFVMKernels::FluxComputeKerne using Base::m_primarySpeciesMobileAggregateConc; using Base::m_referencePorosity; using Base::m_mobilePrimarySpeciesFlags; + using Base::m_solventDensity; using ThermalSinglePhaseFlowAccessors = StencilAccessors< fields::flow::temperature >; @@ -121,6 +122,7 @@ class FluxComputeKernel : public singlePhaseReactiveFVMKernels::FluxComputeKerne * @param[in] thermalConductivityAccessors accessor for wrappers registered by the thermal conductivity model * @param[in] hasDiffusion the flag to turn on diffusion calculation * @param[in] mobilePrimarySpeciesFlags the array of flags to indicate mobile primary species + * @param[in] solventDensity the density of the solvent (e.g., water) [kg/m3] * @param[in] dt time step size * @param[inout] localMatrix the local CRS matrix * @param[inout] localRhs the local right-hand side vector @@ -140,6 +142,7 @@ class FluxComputeKernel : public singlePhaseReactiveFVMKernels::FluxComputeKerne ThermalConductivityAccessors const & thermalConductivityAccessors, integer const & hasDiffusion, arrayView1d< integer const > const & mobilePrimarySpeciesFlags, + real64 const & solventDensity, real64 const & dt, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) @@ -155,6 +158,7 @@ class FluxComputeKernel : public singlePhaseReactiveFVMKernels::FluxComputeKerne porosityAccessors, hasDiffusion, mobilePrimarySpeciesFlags, + solventDensity, dt, localMatrix, localRhs ), @@ -328,15 +332,15 @@ class FluxComputeKernel : public singlePhaseReactiveFVMKernels::FluxComputeKerne // Step 2.2: compute speciesFlux derivative wrt temperature for( integer is = 0; is < numSpecies; ++is ) { - real64 const aggregateConc_i = m_primarySpeciesMobileAggregateConc[er_up][esr_up][ei_up][0][is]; + real64 const aggregateConcMolarity_i = m_primarySpeciesMobileAggregateConc[er_up][esr_up][ei_up][0][is] * m_solventDensity; // real64 const dAggregateConc_i_dTemp = m_dPrimarySpeciesMobileAggregateConcentration_dTemp[er_up][esr_up][ei_up][is]; - // dSpeciesFlux_dT[k_up][is] += dAggregateConc_i_dTemp * fluxVal / fluidDens_up; - dSpeciesFlux_dT[k_up][is] += -aggregateConc_i * fluxVal * dDens_dTemp / (fluidDens_up * fluidDens_up); + // dSpeciesFlux_dT[k_up][is] += dAggregateConc_i_dTemp * m_solventDensity * fluxVal / fluidDens_up; + dSpeciesFlux_dT[k_up][is] += -aggregateConcMolarity_i * fluxVal * dDens_dTemp / (fluidDens_up * fluidDens_up); for( integer ke = 0; ke < numFluxSupportPoints; ++ke ) { - dSpeciesFlux_dT[ke][is] += aggregateConc_i / fluidDens_up * dFlux_dT[ke]; + dSpeciesFlux_dT[ke][is] += aggregateConcMolarity_i / fluidDens_up * dFlux_dT[ke]; } } } @@ -489,7 +493,7 @@ class FluxComputeKernel : public singlePhaseReactiveFVMKernels::FluxComputeKerne // dSpeciesGrad_dT[ke] += stack.diffusionTransmissibility[connectionIndex][ke] // * m_dPrimarySpeciesMobileAggregateConcentration_dTemp[er][esr][ei][is]; - dSpeciesGrad_dT[ke] += stack.dDiffusionTrans_dT[connectionIndex][ke] * m_primarySpeciesMobileAggregateConc[er][esr][ei][0][is]; + dSpeciesGrad_dT[ke] += stack.dDiffusionTrans_dT[connectionIndex][ke] * m_primarySpeciesMobileAggregateConc[er][esr][ei][0][is] * m_solventDensity; } for( integer ke = 0; ke < numFluxSupportPoints; ke++ ) @@ -588,6 +592,7 @@ class FluxComputeKernelFactory createAndLaunch( integer const numSpecies, integer const hasDiffusion, arrayView1d< integer const > const mobilePrimarySpeciesFlags, + real64 const solventDensity, globalIndex const rankOffset, string const & dofKey, string const & solverName, @@ -622,7 +627,7 @@ class FluxComputeKernelFactory KernelType kernel( rankOffset, stencilWrapper, dofNumberAccessor, flowAccessors, reactiveFlowAccessors, thermalFlowAccessors, fluidAccessors, reactiveFluidAccessors, thermalFluidAccessors, permAccessors, diffusionAccessors, porosityAccessors, thermalConductivityAccessors, - hasDiffusion, mobilePrimarySpeciesFlags, dt, localMatrix, localRhs ); + hasDiffusion, mobilePrimarySpeciesFlags, solventDensity, dt, localMatrix, localRhs ); KernelType::template launch< POLICY >( stencilWrapper.size(), kernel ); } ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalSourceFluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalSourceFluxComputeKernel.hpp index 5fa7f2c201c..21f542c262a 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalSourceFluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/reactive/ThermalSourceFluxComputeKernel.hpp @@ -45,6 +45,7 @@ class SourceFluxComputeKernel : public singlePhaseReactiveBaseKernels::SourceFlu using Base::numDof; using Base::numEqn; using Base::m_sizeScalingFactor; + using Base::m_solventDensity; using Base::m_primarySpeciesAggregateConcentration; using Base::m_density; using Base::m_dDensity; @@ -122,7 +123,7 @@ class SourceFluxComputeKernel : public singlePhaseReactiveBaseKernels::SourceFlu for( integer i = 0; i < numSpecies; ++i ) { - stack.localSpeciesJacobian[i][numDof-numSpecies-1] += -m_primarySpeciesAggregateConcentration[ei][0][i] * m_dDensity[ei][0][DerivOffset::dT] / (m_density[ei][0] * m_density[ei][0]) * + stack.localSpeciesJacobian[i][numDof-numSpecies-1] += -m_primarySpeciesAggregateConcentration[ei][0][i] * m_solventDensity * m_dDensity[ei][0][DerivOffset::dT] / (m_density[ei][0] * m_density[ei][0]) * scaledInflowMass; } } 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/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/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..688962c91e2 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 ). @@ -1207,6 +1210,59 @@ void SolidMechanicsLagrangianFEM::assembleSystem( real64 const GEOS_UNUSED_PARAM m_maxForce = LvArray::math::max( mechanicsMaxForce, poromechanicsMaxForce ); } + else if( m_isExplicitChemomechanicsUpdate ) + { + set< string > chemomechanicsRegions; + set< string > mechanicsRegions; + ElementRegionManager const & elementRegionManager = mesh.getElemManager(); + elementRegionManager.forElementSubRegions< CellElementSubRegion >( regionNames, + [&] + ( localIndex const regionIndex, auto & elementSubRegion ) + { + if( elementSubRegion.template hasWrapper< string >( FlowSolverBase::viewKeyStruct::solidNamesString() ) ) + { + chemomechanicsRegions.insert( regionNames[regionIndex] ); + } + else + { + mechanicsRegions.insert( regionNames[regionIndex] ); + } + } ); + + string_array chemomechanicsRegionNames; + chemomechanicsRegionNames.reserve( chemomechanicsRegions.size() ); + for( auto const & region : chemomechanicsRegions ) + { + chemomechanicsRegionNames.emplace_back( region ); + } + string_array mechanicsRegionNames; + mechanicsRegionNames.reserve( mechanicsRegions.size() ); + for( auto const & region : mechanicsRegions ) + { + mechanicsRegionNames.emplace_back( region ); + } + + // first pass for coupled chemomechanics regions + real64 const chemomechanicsMaxForce = assemblyLaunch< SolidMechanicsExplicitChemoMechanicsKernelsDispatchTypeList, + solidMechanicsLagrangianFEMKernels::ExplicitChemoMechanicsFactory >( mesh, + dofManager, + chemomechanicsRegionNames, + FlowSolverBase::viewKeyStruct::solidNamesString(), + localMatrix, + localRhs, + dt ); + // second pass for pure mechanics regions + real64 const mechanicsMaxForce = assemblyLaunch< SolidMechanicsKernelsDispatchTypeList, + solidMechanicsLagrangianFEMKernels::QuasiStaticFactory >( mesh, + dofManager, + mechanicsRegionNames, + viewKeyStruct::solidMaterialNamesString(), + localMatrix, + localRhs, + dt ); + + m_maxForce = LvArray::math::max( chemomechanicsMaxForce, mechanicsMaxForce ); + } else { if( m_timeIntegrationOption == TimeIntegrationOption::QuasiStatic ) @@ -1587,6 +1643,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..3bd4f0e8041 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 >" ] @@ -151,6 +152,40 @@ "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", + "PorousReactiveSolid", + "PorousReactiveSolid" + ], + "FE_TYPE": [ + "H1_Hexahedron_Lagrange1_GaussLegendre2", + "H1_Wedge_Lagrange1_Gauss6", + "H1_Tetrahedron_Lagrange1_Gauss1", + "H1_Pyramid_Lagrange1_Gauss5" + ] + }, + "explicit": [] + }, + "SolidMechanicsALMContactPorousKernels": { "vars": [ "SUBREGION_TYPE", 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@;