diff --git a/inputFiles/chemoMechanics/ChemoMechanicsEigenStrain_CarbonateSystem_1DBar.xml b/inputFiles/chemoMechanics/ChemoMechanicsEigenStrain_CarbonateSystem_1DBar.xml new file mode 100644 index 00000000000..9aa0a31e3a7 --- /dev/null +++ b/inputFiles/chemoMechanics/ChemoMechanicsEigenStrain_CarbonateSystem_1DBar.xml @@ -0,0 +1,58 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/chemoMechanics/ChemoMechanicsEigenStrain_UltramaficSystem_wellbore.xml b/inputFiles/chemoMechanics/ChemoMechanicsEigenStrain_UltramaficSystem_wellbore.xml new file mode 100644 index 00000000000..86e43ad62f1 --- /dev/null +++ b/inputFiles/chemoMechanics/ChemoMechanicsEigenStrain_UltramaficSystem_wellbore.xml @@ -0,0 +1,70 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/chemoMechanics/ChemoMechanicsPoroElastic_UltramaficSystem_wellbore.xml b/inputFiles/chemoMechanics/ChemoMechanicsPoroElastic_UltramaficSystem_wellbore.xml new file mode 100644 index 00000000000..607812a9a45 --- /dev/null +++ b/inputFiles/chemoMechanics/ChemoMechanicsPoroElastic_UltramaficSystem_wellbore.xml @@ -0,0 +1,72 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/chemoMechanics/ChemoMechanicsPoroelastic_CarbonateSystem_1DBar.xml b/inputFiles/chemoMechanics/ChemoMechanicsPoroelastic_CarbonateSystem_1DBar.xml new file mode 100644 index 00000000000..1e5dc3b4146 --- /dev/null +++ b/inputFiles/chemoMechanics/ChemoMechanicsPoroelastic_CarbonateSystem_1DBar.xml @@ -0,0 +1,60 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/chemoMechanics/ChemoMechanics_CarbonateSystem_base.xml b/inputFiles/chemoMechanics/ChemoMechanics_CarbonateSystem_base.xml new file mode 100644 index 00000000000..e2345463eb6 --- /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_base.xml b/inputFiles/chemoMechanics/ChemoMechanics_UltramaficSystem_base.xml new file mode 100644 index 00000000000..1f04420e9b0 --- /dev/null +++ b/inputFiles/chemoMechanics/ChemoMechanics_UltramaficSystem_base.xml @@ -0,0 +1,173 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 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..8aefe228ee1 --- /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/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 6af439fe23f..c9e5e94c354 100644 --- a/src/coreComponents/constitutive/CMakeLists.txt +++ b/src/coreComponents/constitutive/CMakeLists.txt @@ -183,6 +183,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 @@ -190,6 +191,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 @@ -205,10 +207,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 @@ -328,22 +331,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 1c63ae20f7f..f5e3667fa28 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" @@ -361,6 +363,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. */ @@ -405,9 +437,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 ) ); } @@ -415,9 +447,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 ) ); } @@ -464,6 +496,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 >, @@ -489,10 +523,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 > @@ -505,6 +541,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 >, @@ -530,10 +568,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..4c0a07c5374 160000 --- a/src/coreComponents/constitutive/HPCReact +++ b/src/coreComponents/constitutive/HPCReact @@ -1 +1 @@ -Subproject commit 5268dc5c2a59e9dae23f435fba9aecf5a5c5de33 +Subproject commit 4c0a07c53742d0ac2de6cc4870264ccb930548e7 diff --git a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.cpp b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.cpp index 69b7b750934..551ad59bd89 100644 --- a/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.cpp +++ b/src/coreComponents/constitutive/fluid/reactivefluid/ReactiveSinglePhaseFluid.cpp @@ -81,7 +81,7 @@ void ReactiveSinglePhaseFluid< BASE >::postInputInitialization() switch( m_chemicalSystemType ) { case ChemicalSystemType::ultramafic: - m_numPrimarySpecies = 9; + m_numPrimarySpecies = 4; m_numSecondarySpecies = 16; m_numKineticReactions = 5; break; diff --git a/src/coreComponents/constitutive/solid/CoupledSolid.hpp b/src/coreComponents/constitutive/solid/CoupledSolid.hpp index 93e26caa84a..52f173acf58 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..95a655f1f77 --- /dev/null +++ b/src/coreComponents/constitutive/solid/PorousReactiveSolid.cpp @@ -0,0 +1,57 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + + +/** + * @file PorousReactiveSolid.cpp + */ + +#include "PorousReactiveSolid.hpp" +#include "ElasticIsotropic.hpp" +#include "constitutive/permeability/ConstantPermeability.hpp" +#include "constitutive/permeability/CarmanKozenyPermeability.hpp" + +namespace geos +{ + +using namespace dataRepository; + +namespace constitutive +{ + +template< typename SOLID_TYPE, + typename PERM_TYPE > +PorousReactiveSolid< SOLID_TYPE, PERM_TYPE >::PorousReactiveSolid( string const & name, Group * const parent ): + CoupledSolid< SOLID_TYPE, BiotReactivePorosity, PERM_TYPE >( name, parent ) +{} + +template< typename SOLID_TYPE, + typename PERM_TYPE > +void PorousReactiveSolid< SOLID_TYPE, PERM_TYPE >::initializeState() const +{ + CoupledSolid< SOLID_TYPE, BiotReactivePorosity, PERM_TYPE >::initializeState(); +} + +// Register all PorousReactiveSolid model types. +typedef PorousReactiveSolid< ElasticIsotropic, ConstantPermeability > PorousReactiveElasticIsotropicConstant; +typedef PorousReactiveSolid< ElasticIsotropic, CarmanKozenyPermeability > PorousReactiveElasticIsotropicCK; + + +REGISTER_CATALOG_ENTRY( ConstitutiveBase, PorousReactiveElasticIsotropicConstant, string const &, Group * const ) +REGISTER_CATALOG_ENTRY( ConstitutiveBase, PorousReactiveElasticIsotropicCK, string const &, Group * const ) + + +} +} /* namespace geos */ diff --git a/src/coreComponents/constitutive/solid/PorousReactiveSolid.hpp b/src/coreComponents/constitutive/solid/PorousReactiveSolid.hpp new file mode 100644 index 00000000000..0df46325518 --- /dev/null +++ b/src/coreComponents/constitutive/solid/PorousReactiveSolid.hpp @@ -0,0 +1,347 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + + +/** + * @file PorousReactiveSolid.hpp + */ + +#ifndef GEOS_CONSTITUTIVE_SOLID_POROUSREACTIVESOLID_HPP_ +#define GEOS_CONSTITUTIVE_SOLID_POROUSREACTIVESOLID_HPP_ + +#include "constitutive/solid/CoupledSolid.hpp" +#include "constitutive/solid/porosity/BiotReactivePorosity.hpp" +#include "constitutive/solid/SolidBase.hpp" +#include "constitutive/permeability/ConstantPermeability.hpp" + +#include "constitutive/fluid/reactivefluid/ReactiveFluidLayouts.hpp" + +namespace geos +{ +namespace constitutive +{ + +/** + * @brief Provides kernel-callable constitutive update routines + * + * + * @tparam SOLID_TYPE type of the porosity model + */ +template< typename SOLID_TYPE, + typename PERM_TYPE > +class PorousReactiveSolidUpdates : public CoupledSolidUpdates< SOLID_TYPE, BiotReactivePorosity, PERM_TYPE > +{ +public: + + using DiscretizationOps = typename SOLID_TYPE::KernelWrapper::DiscretizationOps; + + /** + * @brief Constructor + */ + PorousReactiveSolidUpdates( SOLID_TYPE const & solidModel, + BiotReactivePorosity const & porosityModel, + PERM_TYPE const & permModel ): + CoupledSolidUpdates< SOLID_TYPE, BiotReactivePorosity, PERM_TYPE >( solidModel, porosityModel, permModel ) + {} + + GEOS_HOST_DEVICE + virtual void updateStateReactionsFixedStress( localIndex const k, + localIndex const q, + real64 const & pressure, + real64 const & pressure_k, + real64 const & pressure_n, + real64 const & temperature, + real64 const & temperature_k, + real64 const & temperature_n, + arraySlice1d< real64 const, compflow::USD_COMP - 1 > mineralReactionMolarIncrements ) const override final + { + updateBiotCoefficientAndAssignModuli( k ); + + m_porosityUpdate.updateFixedStress( k, q, + pressure, pressure_k, pressure_n, + temperature, temperature_k, temperature_n, + mineralReactionMolarIncrements ); + } + + GEOS_HOST_DEVICE + virtual void updateStateFromPressureTemperatureAndReactions( localIndex const k, + localIndex const q, + real64 const & pressure, + real64 const & 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; + + /** + * @brief Const/non-mutable accessor for density + * @return Accessor + */ + arrayView2d< real64 const > const getDensity() const + { + return getSolidModel().getDensity(); + } + +private: + using CoupledSolid< SOLID_TYPE, BiotReactivePorosity, PERM_TYPE >::getSolidModel; + using CoupledSolid< SOLID_TYPE, BiotReactivePorosity, PERM_TYPE >::getPorosityModel; + using CoupledSolid< SOLID_TYPE, BiotReactivePorosity, PERM_TYPE >::getPermModel; +}; + + + +} +} /* namespace geos */ + +#endif /* GEOS_CONSTITUTIVE_SOLID_POROUSREACTIVESOLID_HPP_ */ diff --git a/src/coreComponents/constitutive/solid/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..bbb8adfbc07 --- /dev/null +++ b/src/coreComponents/constitutive/solid/porosity/BiotReactivePorosity.hpp @@ -0,0 +1,320 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file BiotReactivePorosity.hpp + */ + +#ifndef GEOS_CONSTITUTIVE_POROSITY_BIOTREACTIVEPOROSITY_HPP_ +#define GEOS_CONSTITUTIVE_POROSITY_BIOTREACTIVEPOROSITY_HPP_ + +#include "ReactivePorosityBase.hpp" +#include "LvArray/src/tensorOps.hpp" + +namespace geos +{ +namespace constitutive +{ + +class BiotReactivePorosityUpdates : public ReactivePorosityBaseUpdates +{ +public: + + BiotReactivePorosityUpdates( arrayView2d< real64 > const & newPorosity, + arrayView2d< real64 > const & porosity_n, + arrayView2d< real64 > const & dPorosity_dPressure, + arrayView2d< real64 > const & dPorosity_dTemperature, + arrayView2d< real64 > const & initialPorosity, + arrayView1d< real64 > const & referencePorosity, + arrayView3d< real64, reactivefluid::USD_SPECIES > const & volumeFractions, + arrayView3d< real64 const, reactivefluid::USD_SPECIES > const & initialVolumeFractions, + arrayView3d< real64 const, reactivefluid::USD_SPECIES > const & volumeFractions_n, + integer const numKineticReactions, + arrayView1d< real64 const > const & molarWeights, + arrayView1d< real64 const > const & mineralDensities, + arrayView1d< real64 > const & biotCoefficient, + arrayView1d< real64 > const & bulkModulus, + arrayView1d< real64 > const & grainBulkModulus, + arrayView1d< real64 > const & mineralBulkModulus, + arrayView1d< real64 > const & mineralPressure, + arrayView1d< real64 > const & mineralPressure_n, + arrayView2d< real64 > const & meanEffectiveStressIncrement_k ): ReactivePorosityBaseUpdates( newPorosity, + porosity_n, + dPorosity_dPressure, + dPorosity_dTemperature, + initialPorosity, + referencePorosity, + volumeFractions, + initialVolumeFractions, + volumeFractions_n, + numKineticReactions, + molarWeights, + mineralDensities, + bulkModulus, + meanEffectiveStressIncrement_k ), + m_grainBulkModulus( grainBulkModulus ), + m_biotCoefficient( biotCoefficient ), + m_mineralBulkModulus( mineralBulkModulus ), + m_mineralPressure( mineralPressure ), + m_mineralPressure_n( mineralPressure_n ) + {} + + GEOS_HOST_DEVICE + real64 getBiotCoefficient( localIndex const k ) const { return m_biotCoefficient[k]; } + + GEOS_HOST_DEVICE + real64 getGrainBulkModulus( localIndex const k ) const { return m_grainBulkModulus[k]; } + + GEOS_HOST_DEVICE + real64 getPoreMineralPressure( localIndex const k ) const { return m_mineralPressure[k]; } + + GEOS_HOST_DEVICE + void computePorosityFixedStress( real64 const & pressure, + real64 const & pressure_k, + real64 const & pressure_n, + real64 const & porosity_n, + real64 const & referencePorosity, + real64 & porosity, + real64 & dPorosity_dPressure, + real64 const & biotCoefficient, + real64 const & meanEffectiveStressIncrement_k, + real64 const & bulkModulus, + real64 const & grainBulkModulus, + real64 const & mineralBulkModulus, + real64 const & reactionPorosityIncrement ) const + { + GEOS_UNUSED_VAR( pressure_k ); + + real64 const biotSkeletonModulusInverse = (biotCoefficient - referencePorosity) / grainBulkModulus; + real64 const porosityMultiplierInverse = 1 / ( 1 + biotSkeletonModulusInverse*mineralBulkModulus/referencePorosity ); + + porosity = porosity_n + // change due to stress increment + + biotCoefficient * meanEffectiveStressIncrement_k / bulkModulus * porosityMultiplierInverse + // change due to pressure increment + + biotSkeletonModulusInverse * ( pressure - pressure_n ) * porosityMultiplierInverse + // change due to mineral pressure increment + + biotSkeletonModulusInverse * reactionPorosityIncrement * mineralBulkModulus * porosityMultiplierInverse; + + dPorosity_dPressure = biotSkeletonModulusInverse * porosityMultiplierInverse; + } + + GEOS_HOST_DEVICE + void computePoreMineralPressure( real64 const & mineralPressure_n, + real64 & mineralPressure, + real64 & dMineralPres_dMeanEffStressIncre, + real64 const & referencePorosity, + real64 const & biotCoefficient, + real64 const & deltaPressureFromLastStep, + real64 const & meanEffectiveStressIncrement, + real64 const & bulkModulus, + real64 const & grainBulkModulus, + real64 const & mineralBulkModulus, + real64 const & anelasticStrainIncrement ) const + { + // GEOS_UNUSED_VAR( meanEffectiveStressIncrement, bulkModulus ); + + real64 const biotSkeletonModulusInverse = (biotCoefficient - referencePorosity) / grainBulkModulus; + real64 const mineralPressureMultiplier = mineralBulkModulus / ( referencePorosity + biotSkeletonModulusInverse*mineralBulkModulus ); + + mineralPressure = mineralPressure_n + // change due to inelastic strain increment + + mineralPressureMultiplier * referencePorosity * anelasticStrainIncrement + // change due to stress increment + - mineralPressureMultiplier * biotCoefficient * meanEffectiveStressIncrement / bulkModulus + // change due to pressure increment + - mineralPressureMultiplier * biotSkeletonModulusInverse * deltaPressureFromLastStep; + + dMineralPres_dMeanEffStressIncre = -mineralPressureMultiplier * biotCoefficient / bulkModulus; + // dMineralPres_dMeanEffStressIncre = 0.0; + } + + GEOS_HOST_DEVICE + void updatePoreMineralPressure( localIndex const k, + localIndex const q, + real64 const & deltaPressureFromLastStep, + real64 const & meanEffectiveStressIncrement, + real64 const & anelasticStrainIncrement, + real64 & dMineralPres_dMeanEffStressIncre ) const + { + GEOS_UNUSED_VAR( q ); + + computePoreMineralPressure( m_mineralPressure_n[k], + m_mineralPressure[k], + dMineralPres_dMeanEffStressIncre, + m_referencePorosity[k], + m_biotCoefficient[k], + deltaPressureFromLastStep, + meanEffectiveStressIncrement, + m_bulkModulus[k], + m_grainBulkModulus[k], + m_mineralBulkModulus[k], + anelasticStrainIncrement ); + } + + // this function is used in flow solver + // it uses average stress increment (element-based) + GEOS_HOST_DEVICE + virtual void updateFixedStress( localIndex const k, + localIndex const q, + real64 const & pressure, + real64 const & pressure_k, + real64 const & pressure_n, + real64 const & temperature, + real64 const & temperature_k, + real64 const & temperature_n, + arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & kineticReactionMolarIncrements ) const override final + { + // 1. Update the porosity due to reactions + real64 reactionPorosityIncrement = 0.0; + + ReactivePorosityBaseUpdates::computePorosityFromReaction( kineticReactionMolarIncrements, + m_volumeFractions[k][q], + m_volumeFractions_n[k][q], + reactionPorosityIncrement, + m_numKineticReactions, + m_molarWeights, + m_mineralDensities ); + + // 2. Update the porosity due to solid, pore mineral, and pore fluid pressure + // Currently ignore thermal effects + GEOS_UNUSED_VAR( temperature, temperature_k, temperature_n ); + computePorosityFixedStress( pressure, pressure_k, pressure_n, + m_porosity_n[k][q], + m_referencePorosity[k], + m_newPorosity[k][q], + m_dPorosity_dPressure[k][q], + m_biotCoefficient[k], + m_meanEffectiveStressIncrement_k[k][q], + m_bulkModulus[k], + m_grainBulkModulus[k], + m_mineralBulkModulus[k], + reactionPorosityIncrement ); + } + + GEOS_HOST_DEVICE + void updateBiotCoefficientAndAssignModuli( localIndex const k, + real64 const bulkModulus ) const + { + m_bulkModulus[k] = bulkModulus; + + m_biotCoefficient[k] = 1.0 - bulkModulus / m_grainBulkModulus[k]; + } + + + +protected: + + /// View on the grain bulk modulus (read from XML) + arrayView1d< real64 > const m_grainBulkModulus; + + /// View on the Biot coefficient (updated by PorousSolid) + arrayView1d< real64 > const m_biotCoefficient; + + /// View on the mineral bulk modulus (read from XML) + arrayView1d< real64 > const m_mineralBulkModulus; + + /// View on the mineral pressure + arrayView1d< real64 > const m_mineralPressure; + + /// View on the mineral pressure at the previous timestep + arrayView1d< real64 > const m_mineralPressure_n; +}; + +class BiotReactivePorosity : public ReactivePorosityBase +{ +public: + BiotReactivePorosity( string const & name, dataRepository::Group * const parent ); + + static string catalogName() { return "BiotReactivePorosity"; } + + virtual string getCatalogName() const override { return catalogName(); } + + struct viewKeyStruct : public ReactivePorosityBase::viewKeyStruct + { + static constexpr char const *defaultMineralBulkModulusString() { return "defaultMineralBulkModulus"; } + + static constexpr char const *mineralBulkModulusString() { return "mineralBulkModulus"; } + + static constexpr char const *mineralPressure_nString() { return "mineralPressure_n"; } + + static constexpr char const *mineralPressureString() { return "mineralPressure"; } + + static constexpr char const *defaultGrainBulkModulusString() { return "defaultGrainBulkModulus"; } + }; + + virtual void initializeState() const override final; + + virtual void saveConvergedState() const override final; + + using KernelWrapper = BiotReactivePorosityUpdates; + + /** + * @brief Create an update kernel wrapper. + * @return the wrapper + */ + KernelWrapper createKernelUpdates() const + { + return KernelWrapper( m_newPorosity, + m_porosity_n, + m_dPorosity_dPressure, + m_dPorosity_dTemperature, + m_initialPorosity, + m_referencePorosity, + m_volumeFractions, + m_initialVolumeFractions, + m_volumeFractions_n, + m_numKineticReactions, + m_molarWeights, + m_mineralDensities, + m_biotCoefficient, + m_bulkModulus, + m_grainBulkModulus, + m_mineralBulkModulus, + m_mineralPressure, + m_mineralPressure_n, + m_meanEffectiveStressIncrement_k ); + } + +protected: + virtual void postInputInitialization() override; + + /// Biot coefficients (update in the update class, not read in input) + array1d< real64 > m_biotCoefficient; + + /// Grain bulk modulus (read from XML) + real64 m_defaultGrainBulkModulus; + + /// Grain bulk modulus (can be specified in XML) + array1d< real64 > m_grainBulkModulus; + + /// Mineral bulk modulus (read from XML) + real64 m_defaultMineralBulkModulus; + + /// Mineral bulk modulus (can be specified in XML) + array1d< real64 > m_mineralBulkModulus; + + /// Mineral pressure + array1d< real64 > m_mineralPressure; + + /// Mineral pressure at the previous timestep + array1d< real64 > m_mineralPressure_n; +}; + +} /* namespace constitutive */ + +} /* namespace geos */ + +#endif //GEOS_CONSTITUTIVE_POROSITY_BIOTREACTIVEPOROSITY_HPP_ diff --git a/src/coreComponents/constitutive/solid/porosity/ReactivePorosityBase.cpp b/src/coreComponents/constitutive/solid/porosity/ReactivePorosityBase.cpp new file mode 100644 index 00000000000..d75e3a19c48 --- /dev/null +++ b/src/coreComponents/constitutive/solid/porosity/ReactivePorosityBase.cpp @@ -0,0 +1,163 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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.cpp + */ + +#include "ReactivePorosityBase.hpp" + +namespace geos +{ + +using namespace dataRepository; + +namespace constitutive +{ + +ReactivePorosityBase::ReactivePorosityBase( string const & name, Group * const parent ): + PorosityBase( name, parent ) +{ + registerWrapper( viewKeyStruct::defaultInitialVolumeFractionsString(), &m_defaultInitialVolumeFractions ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "Default initial volume fractions" ); + + registerWrapper( viewKeyStruct::initialVolumeFractionsString(), &m_initialVolumeFractions ). + setApplyDefaultValue( 0.0 ). + setPlotLevel( PlotLevel::LEVEL_0 ). + setDescription( "Initial volume fractions" ); + + registerWrapper( viewKeyStruct::volumeFractionsString(), &m_volumeFractions ). + setApplyDefaultValue( 0.0 ). + setPlotLevel( PlotLevel::LEVEL_0 ). + setDescription( "Current volume fractions" ); + + registerWrapper( viewKeyStruct::volumeFractions_nString(), &m_volumeFractions_n ). + setApplyDefaultValue( 0.0 ). + setPlotLevel( PlotLevel::LEVEL_0 ). + setDescription( "Volume fractions at last time step" ); + + registerWrapper( viewKeyStruct::molarWeightsString(), &m_molarWeights ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "Mineral molar weights" ); + + registerWrapper( viewKeyStruct::mineralDensitiesString(), &m_mineralDensities ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "Mineral densities" ); + + registerWrapper( viewKeyStruct::solidBulkModulusString(), &m_bulkModulus ). + setApplyDefaultValue( 1e-6 ). + setDescription( "Solid bulk modulus" ); + + registerWrapper( viewKeyStruct::meanEffectiveStressIncrement_kString(), &m_meanEffectiveStressIncrement_k ). + setApplyDefaultValue( 0.0 ). + setDescription( "Mean effective stress increment at quadrature points at the previous sequential iteration" ); + +} + +std::unique_ptr< ConstitutiveBase > ReactivePorosityBase::deliverClone( string const & name, Group * const parent ) const +{ + std::unique_ptr< ConstitutiveBase > clone = ConstitutiveBase::deliverClone( name, parent ); + + ReactivePorosityBase & newConstitutiveRelation = dynamicCast< ReactivePorosityBase & >( *clone ); + + newConstitutiveRelation.m_numKineticReactions = m_numKineticReactions; + + return clone; +} + +void ReactivePorosityBase::postInputInitialization() +{ + PorosityBase::postInputInitialization(); + + GEOS_THROW_IF_NE_MSG( m_molarWeights.size(), m_defaultInitialVolumeFractions.size(), + GEOS_FMT( "{}: mismatch in number of components in porosity model attribute '{}' and '{}", + getFullName(), viewKeyStruct::molarWeightsString(), viewKeyStruct::initialVolumeFractionsString() ), + InputError ); + + GEOS_THROW_IF_NE_MSG( m_mineralDensities.size(), m_defaultInitialVolumeFractions.size(), + GEOS_FMT( "{}: mismatch in number of components in porosity model attribute '{}' and '{}", + getFullName(), viewKeyStruct::mineralDensitiesString(), viewKeyStruct::initialVolumeFractionsString() ), + InputError ); + + m_numKineticReactions = m_defaultInitialVolumeFractions.size(); +} + +void ReactivePorosityBase::allocateConstitutiveData( dataRepository::Group & parent, + localIndex const numConstitutivePointsPerParentIndex ) +{ + PorosityBase::allocateConstitutiveData( parent, numConstitutivePointsPerParentIndex ); + + resizeFields( parent.size(), numConstitutivePointsPerParentIndex ); + +} + +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 ReactivePorosityBase::saveConvergedState() const +{ + PorosityBase::saveConvergedState(); + + m_volumeFractions_n.setValues< parallelDevicePolicy<> >( m_volumeFractions.toViewConst() ); + m_meanEffectiveStressIncrement_k.zero(); +} + +void ReactivePorosityBase::ignoreConvergedState() const +{ + PorosityBase::ignoreConvergedState(); + m_meanEffectiveStressIncrement_k.zero(); +} + + +void ReactivePorosityBase::initializeState() const +{ + integer const numKineticReactions = this->numKineticReactions(); + + for( localIndex ei = 0; ei < m_newPorosity.size( 0 ); ++ei ) + { + for( localIndex q = 0; q < m_newPorosity.size( 1 ); ++q ) + { + m_newPorosity[ei][q] = m_referencePorosity[ei]; + } + } + + PorosityBase::initializeState(); + + for( localIndex ei = 0; ei < m_volumeFractions.size( 0 ); ++ei ) + { + for( localIndex q = 0; q < m_volumeFractions.size( 1 ); ++q ) + { + for( integer r = 0; r < numKineticReactions; ++r ) + { + m_volumeFractions[ei][q][r] = m_defaultInitialVolumeFractions[r]; + m_initialVolumeFractions[ei][q][r] = m_defaultInitialVolumeFractions[r]; + m_volumeFractions_n[ei][q][r] = m_defaultInitialVolumeFractions[r]; + } + } + } +} + +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..eb7a18d975d --- /dev/null +++ b/src/coreComponents/constitutive/solid/porosity/ReactivePorosityBase.hpp @@ -0,0 +1,301 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file ReactivePorosityBase.hpp + */ + +#ifndef GEOS_CONSTITUTIVE_POROSITY_REACTIVEPOROSITYBASE_HPP_ +#define GEOS_CONSTITUTIVE_POROSITY_REACTIVEPOROSITYBASE_HPP_ + +#include "PorosityBase.hpp" + +#include "constitutive/fluid/reactivefluid/ReactiveFluidLayouts.hpp" + +namespace geos +{ +namespace constitutive +{ + +class ReactivePorosityBaseUpdates : public PorosityBaseUpdates +{ +public: + + ReactivePorosityBaseUpdates( arrayView2d< real64 > const & newPorosity, + arrayView2d< real64 const > const & porosity_n, + arrayView2d< real64 > const & dPorosity_dPressure, + arrayView2d< real64 > const & dPorosity_dTemperature, + arrayView2d< real64 const > const & initialPorosity, + arrayView1d< real64 const > const & referencePorosity, + arrayView3d< real64, reactivefluid::USD_SPECIES > const & volumeFractions, + arrayView3d< real64 const, reactivefluid::USD_SPECIES > const & initialVolumeFractions, + arrayView3d< real64 const, reactivefluid::USD_SPECIES > const & volumeFractions_n, + integer const numKineticReactions, + arrayView1d< real64 const > const & molarWeights, + arrayView1d< real64 const > const & mineralDensities, + arrayView1d< real64 > const & bulkModulus, + arrayView2d< real64 > const & meanEffectiveStressIncrement_k ): + PorosityBaseUpdates( newPorosity, + porosity_n, + dPorosity_dPressure, + dPorosity_dTemperature, + initialPorosity, + referencePorosity ), + m_volumeFractions( volumeFractions ), + m_initialVolumeFractions( initialVolumeFractions ), + m_volumeFractions_n( volumeFractions_n ), + m_numKineticReactions( numKineticReactions ), + m_molarWeights( molarWeights ), + m_mineralDensities( mineralDensities ), + m_bulkModulus( bulkModulus ), + m_meanEffectiveStressIncrement_k( meanEffectiveStressIncrement_k ) + {} + + GEOS_HOST_DEVICE + void computePorosityFromReaction( arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & kineticReactionMolarIncrements, + arraySlice1d< real64, reactivefluid::USD_SPECIES - 2 > const & volumeFractions, + arraySlice1d< real64 const, reactivefluid::USD_SPECIES - 2 > const & volumeFractions_n, + real64 & reactionPorosityIncrement, + integer const & numKineticReactions, + arrayView1d< real64 const > const & molarWeights, + arrayView1d< real64 const > const & mineralDensities ) const + { + for( integer r=0; r < numKineticReactions; ++r ) + { + real64 const volumeFractionIncrement = -kineticReactionMolarIncrements[r] * molarWeights[r]/mineralDensities[r]; + volumeFractions[r] = volumeFractions_n[r] + volumeFractionIncrement; + + reactionPorosityIncrement -= volumeFractionIncrement; + } + } + + GEOS_HOST_DEVICE + void updateFromReactions( localIndex const k, + localIndex const q, + arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & kineticReactionMolarIncrements ) const + { + real64 reactionPorosityIncrement = 0.0; + + computePorosityFromReaction( kineticReactionMolarIncrements, + m_volumeFractions[k][q], + m_volumeFractions_n[k][q], + reactionPorosityIncrement, + m_numKineticReactions, + m_molarWeights, + m_mineralDensities ); + + m_newPorosity[k][q] = m_porosity_n[k][q] + reactionPorosityIncrement; + + if( m_newPorosity[k][q] < 0 ) + { + m_newPorosity[k][q] = 0; + } + else if( m_newPorosity[k][q] > 1.0 ) + { + m_newPorosity[k][q] = 1.0; + } + } + + // this function is used in flow solver + // it uses average stress increment (element-based) + GEOS_HOST_DEVICE + virtual void updateFixedStress( localIndex const k, + localIndex const q, + real64 const & pressure, + real64 const & pressure_k, + real64 const & pressure_n, + real64 const & temperature, + real64 const & temperature_k, + real64 const & temperature_n, + arraySlice1d< real64 const, compflow::USD_COMP - 1 > const & kineticReactionMolarIncrements ) const + { + // For now, we ignore the pressure or temperature dependence but just follow Evans et al. for this eigenstrain approach + GEOS_UNUSED_VAR( pressure, pressure_k, pressure_n, temperature, temperature_k, temperature_n ); + + // 1. Update the porosity due to reactions + real64 reactionPorosityIncrement = 0.0; + + computePorosityFromReaction( kineticReactionMolarIncrements, + m_volumeFractions[k][q], + m_volumeFractions_n[k][q], + reactionPorosityIncrement, + m_numKineticReactions, + m_molarWeights, + m_mineralDensities ); + + // 2. Update the porosity due to solid deformation + m_newPorosity[k][q] = m_porosity_n[k][q] + - reactionPorosityIncrement; + // + m_meanEffectiveStressIncrement_k[k][q]/m_bulkModulus[k]; + + if( m_newPorosity[k][q] < 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; +}; + + +class ReactivePorosityBase : public PorosityBase +{ +public: + ReactivePorosityBase( string const & name, Group * const parent ); + + virtual std::unique_ptr< ConstitutiveBase > + deliverClone( string const & name, + dataRepository::Group * const parent ) const override; + + virtual void allocateConstitutiveData( dataRepository::Group & parent, + localIndex const numConstitutivePointsPerParentIndex ) override; + + static string catalogName() { return "ReactivePorosity"; } + + virtual string getCatalogName() const override { return catalogName(); } + + virtual void saveConvergedState() const override; + virtual void ignoreConvergedState() const override; + + integer numKineticReactions() const { return m_numKineticReactions; } + + virtual void initializeState() const override; + + struct viewKeyStruct : public PorosityBase::viewKeyStruct + { + static constexpr char const * defaultInitialVolumeFractionsString() { return "defaultInitialVolumeFractions"; } + static constexpr char const * initialVolumeFractionsString() { return "initialVolumeFractions"; } + static constexpr char const * volumeFractionsString() { return "volumeFractions"; } + static constexpr char const * volumeFractions_nString() { return "volumeFractions_n"; } + static constexpr char const * molarWeightsString() { return "molarWeights"; } + static constexpr char const * mineralDensitiesString() { return "mineralDensities"; } + static constexpr char const * solidBulkModulusString() { return "solidBulkModulus"; } + static constexpr char const * meanEffectiveStressIncrement_kString() { return "meanEffectiveStressIncrement_k"; } + } viewKeys; + + + using KernelWrapper = ReactivePorosityBaseUpdates; + + /** + * @brief Create an update kernel wrapper. + * @return the wrapper + */ + KernelWrapper createKernelUpdates() const + { + return KernelWrapper( m_newPorosity, + m_porosity_n, + m_dPorosity_dPressure, + m_dPorosity_dTemperature, + m_initialPorosity, + m_referencePorosity, + m_volumeFractions, + m_initialVolumeFractions, + m_volumeFractions_n, + m_numKineticReactions, + m_molarWeights, + m_mineralDensities, + m_bulkModulus, + m_meanEffectiveStressIncrement_k ); + } + + +protected: + virtual void postInputInitialization() override; + + virtual void resizeFields( localIndex const size, localIndex const numPts ); + + array1d< real64 > m_defaultInitialVolumeFractions; + + array3d< real64, constitutive::reactivefluid::LAYOUT_SPECIES > m_volumeFractions; + array3d< real64, constitutive::reactivefluid::LAYOUT_SPECIES > m_initialVolumeFractions; + array3d< real64, constitutive::reactivefluid::LAYOUT_SPECIES > m_volumeFractions_n; + + integer m_numKineticReactions; + array1d< real64 > m_molarWeights; + array1d< real64 > m_mineralDensities; + + array1d< real64 > m_bulkModulus; + array2d< real64 > m_meanEffectiveStressIncrement_k; +}; + + +}/* namespace constitutive */ + +} /* namespace geos */ + + +#endif //GEOS_CONSTITUTIVE_POROSITY_REACTIVEPOROSITYBASE_HPP_ diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.cpp index 5d543342aaa..b10fdba1231 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 {}", @@ -662,15 +691,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 +720,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 +773,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 e8c45a7260a..e4f88ad98fb 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseReactiveTransport.hpp @@ -158,6 +158,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/multiphysics/CoupledSolver.hpp b/src/coreComponents/physicsSolvers/multiphysics/CoupledSolver.hpp index c6da8277ae7..5a6fafb857e 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/CoupledSolver.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/CoupledSolver.hpp @@ -538,7 +538,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 { @@ -615,13 +615,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 359f428a3c3..76eac71b4c3 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 2366f915fd4..5a25e5241ec 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 9471da04798..7435f2e86d6 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp @@ -127,6 +127,12 @@ class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER EnumStrings< SolidMechanicsLagrangianFEM::TimeIntegrationOption >::toString( SolidMechanicsLagrangianFEM::TimeIntegrationOption::QuasiStatic ) ), InputError ); + 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 ); + setMGRStrategy(); } @@ -186,10 +192,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 ) @@ -649,8 +665,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 @@ -663,8 +681,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, @@ -693,7 +715,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 15e1250edbe..bb4908ee13d 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.cpp @@ -180,6 +180,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, @@ -350,9 +392,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 ced15c47ba5..c399776c304 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 @@ -82,7 +84,8 @@ set( kernelTemplateFileList "" ) list( APPEND kernelTemplateFileList kernels/SolidMechanicsKernels.cpp.template - kernels/SolidMechanicsFixedStressThermoPoromechanicsKernels.cpp.template ) + kernels/SolidMechanicsFixedStressThermoPoromechanicsKernels.cpp.template + kernels/SolidMechanicsExplicitChemoMechanicsKernels.cpp.template ) foreach( kernelTemplateFile ${kernelTemplateFileList} ) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp index 18e7714e97d..f2680d3cea8 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" @@ -48,6 +49,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 @@ -69,7 +71,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 ). @@ -1136,6 +1139,22 @@ void SolidMechanicsLagrangianFEM::assembleSystem( real64 const GEOS_UNUSED_PARAM m_maxForce = LvArray::math::max( mechanicsMaxForce, poromechanicsMaxForce ); } + else if( m_isExplicitChemomechanicsUpdate ) + { + + // first pass for coupled poromechanics regions + real64 const chemomechanicsMaxForce= assemblyLaunch< SolidMechanicsExplicitChemoMechanicsKernelsDispatchTypeList, + solidMechanicsLagrangianFEMKernels::ExplicitChemoMechanicsFactory >( mesh, + dofManager, + regionNames, + FlowSolverBase::viewKeyStruct::solidNamesString(), + localMatrix, + localRhs, + dt ); + + + m_maxForce = chemomechanicsMaxForce; + } else { if( m_timeIntegrationOption == TimeIntegrationOption::QuasiStatic ) @@ -1512,6 +1531,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 8681e826d44..d4cb0970f9c 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp @@ -226,6 +226,8 @@ class SolidMechanicsLagrangianFEM : public PhysicsSolverBase void enableFixedStressPoromechanicsUpdate(); + void enableExplicitChemomechanicsUpdate(); + virtual void saveSequentialIterationState( DomainPartition & domain ) override; struct viewKeyStruct : PhysicsSolverBase::viewKeyStruct @@ -307,6 +309,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 9eb41224111..72fe97eb4c9 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 >" ] @@ -149,5 +150,39 @@ ] }, "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": [] } } diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernels/ExplicitChemoMechanics.hpp b/src/coreComponents/physicsSolvers/solidMechanics/kernels/ExplicitChemoMechanics.hpp new file mode 100644 index 00000000000..c0321447f26 --- /dev/null +++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/ExplicitChemoMechanics.hpp @@ -0,0 +1,239 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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() + {} + +#if !defined(CALC_FEM_SHAPE_IN_KERNEL) + /// Dummy + int xLocal; +#else + /// C-array stack storage for element local the nodal positions. + real64 xLocal[ numNodesPerElem ][ 3 ]; +#endif + + /// 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..2d58fd3c6f5 --- /dev/null +++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/ExplicitChemoMechanics_impl.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_impl.hpp + */ + +#ifndef GEOS_PHYSICSSOLVERS_SOLIDMECHANICS_KERNELS_EXPLICITCHEMOMECHANICS_IMPL_HPP_ +#define GEOS_PHYSICSSOLVERS_SOLIDMECHANICS_KERNELS_EXPLICITCHEMOMECHANICS_IMPL_HPP_ + +#include "ExplicitChemoMechanics.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.template numSupportPoints< FE_TYPE >( 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 ) + { +#if defined(CALC_FEM_SHAPE_IN_KERNEL) + stack.xLocal[ a ][ i ] = m_X[ localNodeIndex ][ i ]; +#endif + 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 = m_finiteElementSpace.template getGradN< FE_TYPE >( k, q, stack.xLocal, + stack.feStack, dNdX ); + + real64 strainInc[6] = {0}; + real64 totalStress[6] = {0}; + + typename CONSTITUTIVE_TYPE::KernelWrapper::DiscretizationOps stiffness; + + FE_TYPE::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 ); + FE_TYPE::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.template numSupportPoints< FE_TYPE >( 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..a7fc5514a81 --- /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@ >; \ +template real64 NAME < @SUBREGION_TYPE@, @CONSTITUTIVE_TYPE@, @FE_TYPE@ >::kernelLaunch< NAME##Policy, \ + NAME < @SUBREGION_TYPE@, @CONSTITUTIVE_TYPE@, @FE_TYPE@ > > \ + ( localIndex const, \ + NAME < @SUBREGION_TYPE@, @CONSTITUTIVE_TYPE@, @FE_TYPE@ > 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@;