From 9d37d364e7bba1fd70aacdf0887292a1904a2366 Mon Sep 17 00:00:00 2001 From: Jian HUANG Date: Thu, 18 Sep 2025 15:52:23 -0500 Subject: [PATCH 01/11] add StrainDependentPermeability model --- .../StrainDependentPermeability.cpp | 60 ++++++ .../StrainDependentPermeability.hpp | 191 ++++++++++++++++++ 2 files changed, 251 insertions(+) create mode 100644 src/coreComponents/constitutive/permeability/StrainDependentPermeability.cpp create mode 100644 src/coreComponents/constitutive/permeability/StrainDependentPermeability.hpp diff --git a/src/coreComponents/constitutive/permeability/StrainDependentPermeability.cpp b/src/coreComponents/constitutive/permeability/StrainDependentPermeability.cpp new file mode 100644 index 00000000000..f65745e1930 --- /dev/null +++ b/src/coreComponents/constitutive/permeability/StrainDependentPermeability.cpp @@ -0,0 +1,60 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 StrainDependentPermeability.cpp + */ + +#include "StrainDependentPermeability.hpp" + +namespace geos +{ + +using namespace dataRepository; + +namespace constitutive +{ + + +StrainDependentPermeability::StrainDependentPermeability( string const & name, Group * const parent ): + PermeabilityBase( name, parent ) +{ + registerWrapper( viewKeyStruct::empiricalConstantString(), &m_empiricalConstant ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "an empirical constant." ); + + registerField< fields::permeability::dPerm_dVolStrain >( &m_dPerm_dVolStrain ); +} + +std::unique_ptr< ConstitutiveBase > +StrainDependentPermeability::deliverClone( string const & name, + Group * const parent ) const +{ + return PermeabilityBase::deliverClone( name, parent ); +} + +void StrainDependentPermeability::allocateConstitutiveData( Group & parent, + localIndex const numPts ) +{ + // NOTE: enforcing 1 quadrature point + m_dPerm_dVolStrain.resize( 0, 1, 3 ); + + PermeabilityBase::allocateConstitutiveData( parent, numPts ); +} + +REGISTER_CATALOG_ENTRY( ConstitutiveBase, StrainDependentPermeability, string const &, Group * const ) + +} +} /* namespace geos */ diff --git a/src/coreComponents/constitutive/permeability/StrainDependentPermeability.hpp b/src/coreComponents/constitutive/permeability/StrainDependentPermeability.hpp new file mode 100644 index 00000000000..8509c412b23 --- /dev/null +++ b/src/coreComponents/constitutive/permeability/StrainDependentPermeability.hpp @@ -0,0 +1,191 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 StrainDependentPermeability.hpp + */ + +#ifndef GEOS_CONSTITUTIVE_PERMEABILITY_STRAINDEPENDENTPERMEABILITY_HPP_ +#define GEOS_CONSTITUTIVE_PERMEABILITY_STRAINDEPENDENTPERMEABILITY_HPP_ + +#include "constitutive/permeability/PermeabilityBase.hpp" + + +namespace geos +{ +namespace constitutive +{ + +class StrainDependentPermeabilityUpdate : public PermeabilityBaseUpdate +{ +public: + + StrainDependentPermeabilityUpdate( arrayView3d< real64 > const & permeability, + arrayView3d< real64 > const & dPerm_dPressure, + arrayView3d< real64 > const & dPerm_dPorosity, + arrayView3d< real64 > const & dPerm_dVolStrain, + arrayView1d< real64 > const & referencePorosity, + arrayView3d< real64 > const & referencePermeability, + real64 const empiricalConstant ) + : PermeabilityBaseUpdate( permeability, dPerm_dPressure ), + m_dPerm_dPorosity( dPerm_dPorosity ), + m_dPerm_dVolStrain( dPerm_dVolStrain ), + m_referencePorosity( referencePorosity ), + m_referencePermeability( referencePermeability ), + m_empiricalConstant( empiricalConstant ) + {} + + GEOS_HOST_DEVICE + void compute( real64 const & referencePorosity, + real64 const (&referencePermeability)[3], + real64 const & volStrain, + arraySlice1d< real64 > const & permeability, + arraySlice1d< real64 > const & dPerm_dPorosity, + arraySlice1d< real64 > const & dPerm_dVolStrain ) const; + + GEOS_HOST_DEVICE + virtual void updateFromPressurePorosityAndStrain( localIndex const k, + localIndex const q, + real64 const & volStrain, + real64 const & pressure, + real64 const & porosity ) const override + { + GEOS_UNUSED_VAR( pressure, porosity ); + + real64 referencePermeability[3]; + + referencePermeability[0] = m_referencePermeability[k][0][0]; + referencePermeability[1] = m_referencePermeability[k][0][1]; + referencePermeability[2] = m_referencePermeability[k][0][2]; + + compute( m_referencePorosity[k], + referencePermeability, + volStrain, + m_permeability[k][q], + m_dPerm_dPorosity[k][q], + m_dPerm_dVolStrain[k][q] ); + } + +private: + + /// dPermeability_dPorosity + arrayView3d< real64 > m_dPerm_dPorosity; + + /// dPermeability_dVolumetricStrain + arrayView3d< real64 > m_dPerm_dVolStrain; + + /// An empirical constant + real64 m_empiricalConstant; + + /// Reference permeability + arrayView3d< real64 > m_referencePermeability; + + /// Reference porosity + arrayView1d< real64 > m_referencePorosity; +}; + + +class StrainDependentPermeability : public PermeabilityBase +{ +public: + + StrainDependentPermeability( string const & name, dataRepository::Group * const parent ); + + std::unique_ptr< ConstitutiveBase > deliverClone( string const & name, + dataRepository::Group * const parent ) const override; + + static string catalogName() { return "StrainDependentPermeability"; } + + virtual string getCatalogName() const override { return catalogName(); } + + virtual void allocateConstitutiveData( dataRepository::Group & parent, + localIndex const numPts ) override; + + /// Type of kernel wrapper for in-kernel update + using KernelWrapper = StrainDependentPermeabilityUpdate; + + /** + * @brief Create an update kernel wrapper. + * @return the wrapper + */ + KernelWrapper createKernelWrapper() const + { + return KernelWrapper( m_permeability, + m_dPerm_dPressure, + m_dPerm_dPorosity, + m_dPerm_dVolStrain, + m_empiricalConstant, + m_referencePorosity, + m_referencePermeability ); + } + + + struct viewKeyStruct : public PermeabilityBase::viewKeyStruct + { + static constexpr char const * dPerm_dPorosityString() { return "dPerm_dPorosity"; } + static constexpr char const * dPerm_dVolStrainString() { return "dPerm_dVolStrain"; } + static constexpr char const * empiricalConstantString() { return "empiricalConstant"; } + static constexpr char const * referencePorosityString() { return "referencePorosity"; } + static constexpr char const * referencePermeabilityString() { return "referencePermeability"; } + }; + +private: + + /// dPermeability_dPorosity + array3d< real64 > m_dPerm_dPorosity; + + /// dPermeability_dVolumetricStrain + array3d< real64 > m_dPerm_dVolStrain; + + /// An empirical constant + real64 m_empiricalConstant; + + /// Reference permeability + array3d< real64 > m_referencePermeability; + + /// Reference porosity + array1d< real64 > m_referencePorosity; +}; + + +GEOS_HOST_DEVICE +inline +void StrainDependentPermeabilityUpdate::compute( real64 const & referencePorosity, + real64 const (&referencePermeability)[3], + real64 const & volStrain, + arraySlice1d< real64 > const & permeability, + arraySlice1d< real64 > const & dPerm_dPorosity, + arraySlice1d< real64 > const & dPerm_dVolStrain ) const +{ + real64 const por = 1 - (1 - referencePorosity) * LvArray::math::exp(-volStrain); + + real64 const permMultiplier = LvArray::math::pow( por/referencePorosity, m_empiricalConstant ); + + real64 const dpermMultiplier_dVolStrain = m_empiricalConstant * LvArray::math::pow( por/referencePorosity, m_empiricalConstant-1 ) /referencePorosity * ((1 - referencePorosity) * LvArray::math::exp(-volStrain)); + + for( localIndex i = 0; i < permeability.size(); ++i ) + { + permeability[i] = permMultiplier * referencePermeability[i]; + dPerm_dVolStrain[i] = dpermMultiplier_dVolStrain * referencePermeability[i]; + } +} + + +}/* namespace constitutive */ + +} /* namespace geos */ + + +#endif //GEOS_CONSTITUTIVE_PERMEABILITY_STRAINDEPENDENTPERMEABILITY_HPP_ From ddfc5858d9a75cbac7f633b7ef32edc70dec95b0 Mon Sep 17 00:00:00 2001 From: Jian HUANG Date: Thu, 2 Oct 2025 15:44:27 -0500 Subject: [PATCH 02/11] register in kernels --- .../constitutive/ConstitutivePassThru.hpp | 37 ++++++++++++++++++- .../constitutive/solid/PorousSolid.cpp | 20 ++++++++++ .../multiphysics/kernelSpecs.json | 20 +++++++++- .../solidMechanics/kernelSpecs.json | 8 +++- 4 files changed, 79 insertions(+), 6 deletions(-) diff --git a/src/coreComponents/constitutive/ConstitutivePassThru.hpp b/src/coreComponents/constitutive/ConstitutivePassThru.hpp index 48467e5c1d8..6bc24b9645d 100644 --- a/src/coreComponents/constitutive/ConstitutivePassThru.hpp +++ b/src/coreComponents/constitutive/ConstitutivePassThru.hpp @@ -316,8 +316,19 @@ struct ConstitutivePassThru< PorousSolidBase > PorousSolid< ElasticIsotropic, CarmanKozenyPermeability >, PorousSolid< ElasticTransverseIsotropic, CarmanKozenyPermeability >, PorousSolid< ElasticIsotropicPressureDependent, CarmanKozenyPermeability >, - PorousSolid< ElasticOrthotropic, CarmanKozenyPermeability > >::execute( constitutiveRelation, - std::forward< LAMBDA >( lambda ) ); + PorousSolid< ElasticOrthotropic, CarmanKozenyPermeability >, + PorousSolid< DruckerPragerExtended, StrainDependentPermeability >, + PorousSolid< ModifiedCamClay, StrainDependentPermeability >, + PorousSolid< DelftEgg, StrainDependentPermeability >, + PorousSolid< DruckerPrager, StrainDependentPermeability >, + PorousSolid< DuvautLionsSolid< DruckerPrager >, StrainDependentPermeability >, + PorousSolid< DuvautLionsSolid< DruckerPragerExtended >, StrainDependentPermeability >, + PorousSolid< DuvautLionsSolid< ModifiedCamClay >, StrainDependentPermeability >, + PorousSolid< ElasticIsotropic, StrainDependentPermeability >, + PorousSolid< ElasticTransverseIsotropic, StrainDependentPermeability >, + PorousSolid< ElasticIsotropicPressureDependent, StrainDependentPermeability >, + PorousSolid< ElasticOrthotropic, StrainDependentPermeability > >::execute( constitutiveRelation, + std::forward< LAMBDA >( lambda ) ); } }; @@ -434,6 +445,17 @@ struct ConstitutivePassThru< CoupledSolidBase > PorousSolid< ElasticTransverseIsotropic, CarmanKozenyPermeability >, PorousSolid< ElasticIsotropicPressureDependent, CarmanKozenyPermeability >, PorousSolid< ElasticOrthotropic, CarmanKozenyPermeability >, + PorousSolid< DruckerPragerExtended, StrainDependentPermeability >, + PorousSolid< ModifiedCamClay, StrainDependentPermeability >, + PorousSolid< DelftEgg, StrainDependentPermeability >, + PorousSolid< DruckerPrager, StrainDependentPermeability >, + PorousSolid< DuvautLionsSolid< DruckerPrager >, StrainDependentPermeability >, + PorousSolid< DuvautLionsSolid< DruckerPragerExtended >, StrainDependentPermeability >, + PorousSolid< DuvautLionsSolid< ModifiedCamClay >, StrainDependentPermeability >, + PorousSolid< ElasticIsotropic, StrainDependentPermeability >, + PorousSolid< ElasticTransverseIsotropic, StrainDependentPermeability >, + PorousSolid< ElasticIsotropicPressureDependent, StrainDependentPermeability >, + PorousSolid< ElasticOrthotropic, StrainDependentPermeability >, PorousDamageSolid< DamageSpectral< ElasticIsotropic > >, PorousDamageSolid< DamageVolDev< ElasticIsotropic > >, PorousDamageSolid< Damage< ElasticIsotropic > > >::execute( constitutiveRelation, @@ -472,6 +494,17 @@ struct ConstitutivePassThru< CoupledSolidBase > PorousSolid< ElasticTransverseIsotropic, CarmanKozenyPermeability >, PorousSolid< ElasticIsotropicPressureDependent, CarmanKozenyPermeability >, PorousSolid< ElasticOrthotropic, CarmanKozenyPermeability >, + PorousSolid< DruckerPragerExtended, StrainDependentPermeability >, + PorousSolid< ModifiedCamClay, StrainDependentPermeability >, + PorousSolid< DelftEgg, StrainDependentPermeability >, + PorousSolid< DruckerPrager, StrainDependentPermeability >, + PorousSolid< DuvautLionsSolid< DruckerPrager >, StrainDependentPermeability >, + PorousSolid< DuvautLionsSolid< DruckerPragerExtended >, StrainDependentPermeability >, + PorousSolid< DuvautLionsSolid< ModifiedCamClay >, StrainDependentPermeability >, + PorousSolid< ElasticIsotropic, StrainDependentPermeability >, + PorousSolid< ElasticTransverseIsotropic, StrainDependentPermeability >, + PorousSolid< ElasticIsotropicPressureDependent, StrainDependentPermeability >, + PorousSolid< ElasticOrthotropic, StrainDependentPermeability >, PorousDamageSolid< DamageSpectral< ElasticIsotropic > >, PorousDamageSolid< DamageVolDev< ElasticIsotropic > >, PorousDamageSolid< Damage< ElasticIsotropic > > >::execute( constitutiveRelation, diff --git a/src/coreComponents/constitutive/solid/PorousSolid.cpp b/src/coreComponents/constitutive/solid/PorousSolid.cpp index 5b871f186e6..10e23bc1a9b 100644 --- a/src/coreComponents/constitutive/solid/PorousSolid.cpp +++ b/src/coreComponents/constitutive/solid/PorousSolid.cpp @@ -72,6 +72,16 @@ typedef PorousSolid< DuvautLionsSolid< DruckerPrager >, CarmanKozenyPermeability typedef PorousSolid< DuvautLionsSolid< DruckerPragerExtended >, CarmanKozenyPermeability > PorousViscoDruckerPragerExtendedCK; typedef PorousSolid< DuvautLionsSolid< ModifiedCamClay >, CarmanKozenyPermeability > PorousViscoModifiedCamClayCK; typedef PorousSolid< ModifiedCamClay, CarmanKozenyPermeability > PorousModifiedCamClayCK; +typedef PorousSolid< ElasticIsotropic, StrainDependentPermeability > PorousElasticIsotropicSD; +typedef PorousSolid< ElasticTransverseIsotropic, StrainDependentPermeability > PorousElasticTransverseIsotropicSD; +typedef PorousSolid< ElasticOrthotropic, StrainDependentPermeability > PorousElasticOrthotropicSD; +typedef PorousSolid< DelftEgg, StrainDependentPermeability > PorousDelftEggSD; +typedef PorousSolid< DruckerPrager, StrainDependentPermeability > PorousDruckerPragerSD; +typedef PorousSolid< DruckerPragerExtended, StrainDependentPermeability > PorousDruckerPragerExtendedSD; +typedef PorousSolid< DuvautLionsSolid< DruckerPrager >, StrainDependentPermeability > PorousViscoDruckerPragerSD; +typedef PorousSolid< DuvautLionsSolid< DruckerPragerExtended >, StrainDependentPermeability > PorousViscoDruckerPragerExtendedSD; +typedef PorousSolid< DuvautLionsSolid< ModifiedCamClay >, StrainDependentPermeability > PorousViscoModifiedCamClaySD; +typedef PorousSolid< ModifiedCamClay, StrainDependentPermeability > PorousModifiedCamClaySD; REGISTER_CATALOG_ENTRY( ConstitutiveBase, PorousElasticIsotropicConstant, string const &, Group * const ) @@ -94,6 +104,16 @@ REGISTER_CATALOG_ENTRY( ConstitutiveBase, PorousModifiedCamClayCK, string const REGISTER_CATALOG_ENTRY( ConstitutiveBase, PorousViscoDruckerPragerCK, string const &, Group * const ) REGISTER_CATALOG_ENTRY( ConstitutiveBase, PorousViscoDruckerPragerExtendedCK, string const &, Group * const ) REGISTER_CATALOG_ENTRY( ConstitutiveBase, PorousViscoModifiedCamClayCK, string const &, Group * const ) +REGISTER_CATALOG_ENTRY( ConstitutiveBase, PorousElasticIsotropicSD, string const &, Group * const ) +REGISTER_CATALOG_ENTRY( ConstitutiveBase, PorousElasticTransverseIsotropicSD, string const &, Group * const ) +REGISTER_CATALOG_ENTRY( ConstitutiveBase, PorousElasticOrthotropicSD, string const &, Group * const ) +REGISTER_CATALOG_ENTRY( ConstitutiveBase, PorousDelftEggSD, string const &, Group * const ) +REGISTER_CATALOG_ENTRY( ConstitutiveBase, PorousDruckerPragerSD, string const &, Group * const ) +REGISTER_CATALOG_ENTRY( ConstitutiveBase, PorousDruckerPragerExtendedSD, string const &, Group * const ) +REGISTER_CATALOG_ENTRY( ConstitutiveBase, PorousModifiedCamClaySD, string const &, Group * const ) +REGISTER_CATALOG_ENTRY( ConstitutiveBase, PorousViscoDruckerPragerSD, string const &, Group * const ) +REGISTER_CATALOG_ENTRY( ConstitutiveBase, PorousViscoDruckerPragerExtendedSD, string const &, Group * const ) +REGISTER_CATALOG_ENTRY( ConstitutiveBase, PorousViscoModifiedCamClaySD, string const &, Group * const ) } diff --git a/src/coreComponents/physicsSolvers/multiphysics/kernelSpecs.json b/src/coreComponents/physicsSolvers/multiphysics/kernelSpecs.json index bf248e898f3..b68e4240cfa 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/kernelSpecs.json +++ b/src/coreComponents/physicsSolvers/multiphysics/kernelSpecs.json @@ -61,6 +61,17 @@ "PorousSolid, CarmanKozenyPermeability>", "PorousSolid, CarmanKozenyPermeability>", "PorousSolid, CarmanKozenyPermeability>", + "PorousSolid", + "PorousSolid", + "PorousSolid", + "PorousSolid", + "PorousSolid", + "PorousSolid", + "PorousSolid", + "PorousSolid", + "PorousSolid, StrainDependentPermeability>", + "PorousSolid, StrainDependentPermeability>", + "PorousSolid, StrainDependentPermeability>", "PorousDamageSolid>", "PorousDamageSolid>", "PorousDamageSolid>" @@ -100,7 +111,8 @@ ], "CONSTITUTIVE_TYPE": [ "PorousSolid", - "PorousSolid" + "PorousSolid", + "PorousSolid" ], "FE_TYPE": [ "H1_Hexahedron_Lagrange1_GaussLegendre2", @@ -136,7 +148,11 @@ "PorousSolid", "PorousSolid", "PorousSolid", - "PorousSolid" + "PorousSolid", + "PorousSolid", + "PorousSolid", + "PorousSolid", + "PorousSolid" ], "FE_TYPE": [ "H1_Hexahedron_Lagrange1_GaussLegendre2", diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernelSpecs.json b/src/coreComponents/physicsSolvers/solidMechanics/kernelSpecs.json index ca527443f5a..a964c4c10dc 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/kernelSpecs.json +++ b/src/coreComponents/physicsSolvers/solidMechanics/kernelSpecs.json @@ -94,10 +94,13 @@ "CONSTITUTIVE_TYPE": [ "PorousSolid", "PorousSolid", + "PorousSolid", "PorousSolid", "PorousSolid", + "PorousSolid", "PorousSolid, ConstantPermeability>", - "PorousSolid, CarmanKozenyPermeability>" + "PorousSolid, CarmanKozenyPermeability>", + "PorousSolid, StrainDependentPermeability>" ], "FE_TYPE": [ "H1_Hexahedron_Lagrange1_GaussLegendre2", @@ -134,7 +137,8 @@ ], "CONSTITUTIVE_TYPE": [ "PorousSolid", - "PorousSolid" + "PorousSolid", + "PorousSolid", ], "FE_TYPE": [ "H1_Hexahedron_Lagrange1_GaussLegendre2", From ffca9cd3320245f7d21e9c768f4a2695c6ac4f80 Mon Sep 17 00:00:00 2001 From: Jian HUANG Date: Thu, 2 Oct 2025 15:49:08 -0500 Subject: [PATCH 03/11] update PorousSolid.cpp --- src/coreComponents/constitutive/solid/PorousSolid.cpp | 1 + .../physicsSolvers/solidMechanics/kernelSpecs.json | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/src/coreComponents/constitutive/solid/PorousSolid.cpp b/src/coreComponents/constitutive/solid/PorousSolid.cpp index 10e23bc1a9b..2973c423f41 100644 --- a/src/coreComponents/constitutive/solid/PorousSolid.cpp +++ b/src/coreComponents/constitutive/solid/PorousSolid.cpp @@ -29,6 +29,7 @@ #include "DuvautLionsSolid.hpp" #include "constitutive/permeability/ConstantPermeability.hpp" #include "constitutive/permeability/CarmanKozenyPermeability.hpp" +#include "constitutive/permeability/StrainDependentPermeability.hpp" namespace geos { diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernelSpecs.json b/src/coreComponents/physicsSolvers/solidMechanics/kernelSpecs.json index a964c4c10dc..b9450368482 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/kernelSpecs.json +++ b/src/coreComponents/physicsSolvers/solidMechanics/kernelSpecs.json @@ -138,7 +138,7 @@ "CONSTITUTIVE_TYPE": [ "PorousSolid", "PorousSolid", - "PorousSolid", + "PorousSolid" ], "FE_TYPE": [ "H1_Hexahedron_Lagrange1_GaussLegendre2", From 07d4c6aa104d5430815b384f8116153f2b7cda4e Mon Sep 17 00:00:00 2001 From: Jian HUANG Date: Thu, 2 Oct 2025 16:59:32 -0500 Subject: [PATCH 04/11] add referencePermeability and strainDependenceConstants --- .../StrainDependentPermeability.cpp | 57 ++++++++++++++++++- .../StrainDependentPermeability.hpp | 33 +++++++---- 2 files changed, 76 insertions(+), 14 deletions(-) diff --git a/src/coreComponents/constitutive/permeability/StrainDependentPermeability.cpp b/src/coreComponents/constitutive/permeability/StrainDependentPermeability.cpp index f65745e1930..e1518ffc8d4 100644 --- a/src/coreComponents/constitutive/permeability/StrainDependentPermeability.cpp +++ b/src/coreComponents/constitutive/permeability/StrainDependentPermeability.cpp @@ -31,9 +31,19 @@ namespace constitutive StrainDependentPermeability::StrainDependentPermeability( string const & name, Group * const parent ): PermeabilityBase( name, parent ) { - registerWrapper( viewKeyStruct::empiricalConstantString(), &m_empiricalConstant ). + registerWrapper( viewKeyStruct::referencePermeabilityComponentsString(), &m_referencePermeabilityComponents ). setInputFlag( InputFlags::REQUIRED ). - setDescription( "an empirical constant." ); + setRestartFlags( RestartFlags::NO_WRITE ). + setDescription( "Reference xx, yy and zz components of a diagonal permeability tensor." ); + + registerWrapper( viewKeyStruct::strainDependenceConstantsString(), &m_strainDependenceConstants ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "Volumetric strain dependence coefficients for each permeability component." ); + + registerWrapper( viewKeyStruct::referencePermeabilityString(), &m_referencePermeability ). + setApplyDefaultValue( 0.0 ). + setPlotLevel( PlotLevel::LEVEL_0 ). + setDescription( "Reference permeability field" ); registerField< fields::permeability::dPerm_dVolStrain >( &m_dPerm_dVolStrain ); } @@ -51,9 +61,52 @@ void StrainDependentPermeability::allocateConstitutiveData( Group & parent, // NOTE: enforcing 1 quadrature point m_dPerm_dVolStrain.resize( 0, 1, 3 ); + m_referencePermeability.resize( 0, 1, 3 ); + PermeabilityBase::allocateConstitutiveData( parent, numPts ); + + integer constexpr numQuad = 1; // NOTE: enforcing 1 quadrature point + + for( localIndex ei = 0; ei < parent.size(); ++ei ) + { + for( localIndex q = 0; q < numQuad; ++q ) + { + m_referencePermeability[ei][q][0] = m_referencePermeabilityComponents[0]; + m_referencePermeability[ei][q][1] = m_referencePermeabilityComponents[1]; + m_referencePermeability[ei][q][2] = m_referencePermeabilityComponents[2]; + } + } } +void StrainDependentPermeability::initializeState() const +{ + localIndex const numE = m_permeability.size( 0 ); + integer constexpr numQuad = 1; // NOTE: enforcing 1 quadrature point + + auto permView = m_permeability.toView(); + real64 const permComponents[3] = { m_referencePermeabilityComponents[0], + m_referencePermeabilityComponents[1], + m_referencePermeabilityComponents[2] }; + + forAll< parallelDevicePolicy<> >( numE, [=] GEOS_HOST_DEVICE ( localIndex const ei ) + { + for( localIndex q = 0; q < numQuad; ++q ) + { + for( integer dim=0; dim < 3; ++dim ) + { + // The default value is -1 so if it still -1 it needs to be set to something physical + if( permView[ei][q][dim] < 0 ) + { + permView[ei][q][dim] = permComponents[dim]; + } + } + } + } ); +} + +void StrainDependentPermeability::postInputInitialization() +{} + REGISTER_CATALOG_ENTRY( ConstitutiveBase, StrainDependentPermeability, string const &, Group * const ) } diff --git a/src/coreComponents/constitutive/permeability/StrainDependentPermeability.hpp b/src/coreComponents/constitutive/permeability/StrainDependentPermeability.hpp index 8509c412b23..3190edb7421 100644 --- a/src/coreComponents/constitutive/permeability/StrainDependentPermeability.hpp +++ b/src/coreComponents/constitutive/permeability/StrainDependentPermeability.hpp @@ -38,18 +38,19 @@ class StrainDependentPermeabilityUpdate : public PermeabilityBaseUpdate arrayView3d< real64 > const & dPerm_dVolStrain, arrayView1d< real64 > const & referencePorosity, arrayView3d< real64 > const & referencePermeability, - real64 const empiricalConstant ) + real64 const strainDependenceConstants ) : PermeabilityBaseUpdate( permeability, dPerm_dPressure ), m_dPerm_dPorosity( dPerm_dPorosity ), m_dPerm_dVolStrain( dPerm_dVolStrain ), m_referencePorosity( referencePorosity ), m_referencePermeability( referencePermeability ), - m_empiricalConstant( empiricalConstant ) + m_strainDependenceConstants( strainDependenceConstants ) {} GEOS_HOST_DEVICE void compute( real64 const & referencePorosity, real64 const (&referencePermeability)[3], + R1Tensor const strainDependenceConstants, real64 const & volStrain, arraySlice1d< real64 > const & permeability, arraySlice1d< real64 > const & dPerm_dPorosity, @@ -72,6 +73,7 @@ class StrainDependentPermeabilityUpdate : public PermeabilityBaseUpdate compute( m_referencePorosity[k], referencePermeability, + m_strainDependenceConstants, volStrain, m_permeability[k][q], m_dPerm_dPorosity[k][q], @@ -86,8 +88,8 @@ class StrainDependentPermeabilityUpdate : public PermeabilityBaseUpdate /// dPermeability_dVolumetricStrain arrayView3d< real64 > m_dPerm_dVolStrain; - /// An empirical constant - real64 m_empiricalConstant; + /// Volumetric strain dependence coefficients for each permeability component + R1Tensor m_strainDependenceConstants; /// Reference permeability arrayView3d< real64 > m_referencePermeability; @@ -126,7 +128,7 @@ class StrainDependentPermeability : public PermeabilityBase m_dPerm_dPressure, m_dPerm_dPorosity, m_dPerm_dVolStrain, - m_empiricalConstant, + m_strainDependenceConstants, m_referencePorosity, m_referencePermeability ); } @@ -136,11 +138,17 @@ class StrainDependentPermeability : public PermeabilityBase { static constexpr char const * dPerm_dPorosityString() { return "dPerm_dPorosity"; } static constexpr char const * dPerm_dVolStrainString() { return "dPerm_dVolStrain"; } - static constexpr char const * empiricalConstantString() { return "empiricalConstant"; } + static constexpr char const * strainDependenceConstantsString() { return "strainDependenceConstants"; } static constexpr char const * referencePorosityString() { return "referencePorosity"; } static constexpr char const * referencePermeabilityString() { return "referencePermeability"; } }; + virtual void initializeState() const override final; + +protected: + + virtual void postInputInitialization() override; + private: /// dPermeability_dPorosity @@ -149,8 +157,8 @@ class StrainDependentPermeability : public PermeabilityBase /// dPermeability_dVolumetricStrain array3d< real64 > m_dPerm_dVolStrain; - /// An empirical constant - real64 m_empiricalConstant; + /// Volumetric strain dependence coefficients for each permeability component + R1Tensor m_strainDependenceConstants; /// Reference permeability array3d< real64 > m_referencePermeability; @@ -164,6 +172,7 @@ GEOS_HOST_DEVICE inline void StrainDependentPermeabilityUpdate::compute( real64 const & referencePorosity, real64 const (&referencePermeability)[3], + R1Tensor const strainDependenceConstants, real64 const & volStrain, arraySlice1d< real64 > const & permeability, arraySlice1d< real64 > const & dPerm_dPorosity, @@ -171,12 +180,12 @@ void StrainDependentPermeabilityUpdate::compute( real64 const & referencePorosit { real64 const por = 1 - (1 - referencePorosity) * LvArray::math::exp(-volStrain); - real64 const permMultiplier = LvArray::math::pow( por/referencePorosity, m_empiricalConstant ); - - real64 const dpermMultiplier_dVolStrain = m_empiricalConstant * LvArray::math::pow( por/referencePorosity, m_empiricalConstant-1 ) /referencePorosity * ((1 - referencePorosity) * LvArray::math::exp(-volStrain)); - for( localIndex i = 0; i < permeability.size(); ++i ) { + real64 const permMultiplier = LvArray::math::pow( por/referencePorosity, strainDependenceConstants[i] ); + + real64 const dpermMultiplier_dVolStrain = strainDependenceConstants[i] * LvArray::math::pow( por/referencePorosity, strainDependenceConstants[i]-1 ) /referencePorosity * ((1 - referencePorosity) * LvArray::math::exp(-volStrain)); + permeability[i] = permMultiplier * referencePermeability[i]; dPerm_dVolStrain[i] = dpermMultiplier_dVolStrain * referencePermeability[i]; } From c2b355cdf91c7b93663e972619a473faa6790a4d Mon Sep 17 00:00:00 2001 From: Jian HUANG Date: Thu, 2 Oct 2025 17:17:17 -0500 Subject: [PATCH 05/11] updateFromPorosityAndStrain --- .../constitutive/permeability/PermeabilityBase.hpp | 9 +++++++++ .../permeability/StrainDependentPermeability.hpp | 11 +++++------ 2 files changed, 14 insertions(+), 6 deletions(-) diff --git a/src/coreComponents/constitutive/permeability/PermeabilityBase.hpp b/src/coreComponents/constitutive/permeability/PermeabilityBase.hpp index 190b8f1b4f6..1a68fc3c634 100644 --- a/src/coreComponents/constitutive/permeability/PermeabilityBase.hpp +++ b/src/coreComponents/constitutive/permeability/PermeabilityBase.hpp @@ -56,6 +56,15 @@ class PermeabilityBaseUpdate GEOS_UNUSED_VAR( k, q, pressure, porosity ); } + GEOS_HOST_DEVICE + virtual void updateFromPorosityAndStrain( localIndex const k, + localIndex const q, + real64 const & volStrain, + real64 const & porosity ) const + { + GEOS_UNUSED_VAR( k, q, volStrain, porosity ); + } + GEOS_HOST_DEVICE virtual void updateFromAperture( localIndex const k, localIndex const q, diff --git a/src/coreComponents/constitutive/permeability/StrainDependentPermeability.hpp b/src/coreComponents/constitutive/permeability/StrainDependentPermeability.hpp index 3190edb7421..717e8de2e80 100644 --- a/src/coreComponents/constitutive/permeability/StrainDependentPermeability.hpp +++ b/src/coreComponents/constitutive/permeability/StrainDependentPermeability.hpp @@ -57,13 +57,12 @@ class StrainDependentPermeabilityUpdate : public PermeabilityBaseUpdate arraySlice1d< real64 > const & dPerm_dVolStrain ) const; GEOS_HOST_DEVICE - virtual void updateFromPressurePorosityAndStrain( localIndex const k, - localIndex const q, - real64 const & volStrain, - real64 const & pressure, - real64 const & porosity ) const override + virtual void updateFromPorosityAndStrain( localIndex const k, + localIndex const q, + real64 const & volStrain, + real64 const & porosity ) const override { - GEOS_UNUSED_VAR( pressure, porosity ); + GEOS_UNUSED_VAR( porosity ); real64 referencePermeability[3]; From 8e43665a1faa4d1abaa0165e20eb2496203f8e4c Mon Sep 17 00:00:00 2001 From: Jian HUANG Date: Thu, 2 Oct 2025 17:26:24 -0500 Subject: [PATCH 06/11] get VolStrain --- .../constitutive/permeability/StrainDependentPermeability.hpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/coreComponents/constitutive/permeability/StrainDependentPermeability.hpp b/src/coreComponents/constitutive/permeability/StrainDependentPermeability.hpp index 717e8de2e80..1feb2c97c3a 100644 --- a/src/coreComponents/constitutive/permeability/StrainDependentPermeability.hpp +++ b/src/coreComponents/constitutive/permeability/StrainDependentPermeability.hpp @@ -70,6 +70,8 @@ class StrainDependentPermeabilityUpdate : public PermeabilityBaseUpdate referencePermeability[1] = m_referencePermeability[k][0][1]; referencePermeability[2] = m_referencePermeability[k][0][2]; + volStrain = m_constitutiveUpdate.getVolStrain( k, q ); + compute( m_referencePorosity[k], referencePermeability, m_strainDependenceConstants, From 5c426a29a4a19059256a703c0374822d5fd139bf Mon Sep 17 00:00:00 2001 From: Jian HUANG Date: Thu, 2 Oct 2025 18:16:27 -0500 Subject: [PATCH 07/11] fix compilation errors --- .../constitutive/ConstitutivePassThru.hpp | 1 + .../StrainDependentPermeability.hpp | 34 +++++++++---------- 2 files changed, 18 insertions(+), 17 deletions(-) diff --git a/src/coreComponents/constitutive/ConstitutivePassThru.hpp b/src/coreComponents/constitutive/ConstitutivePassThru.hpp index 6bc24b9645d..2376464c3b6 100644 --- a/src/coreComponents/constitutive/ConstitutivePassThru.hpp +++ b/src/coreComponents/constitutive/ConstitutivePassThru.hpp @@ -45,6 +45,7 @@ #include "solid/porosity/ProppantPorosity.hpp" #include "permeability/ConstantPermeability.hpp" #include "permeability/CarmanKozenyPermeability.hpp" +#include "permeability/StrainDependentPermeability.hpp" #include "permeability/ExponentialDecayPermeability.hpp" #include "permeability/ParallelPlatesPermeability.hpp" #include "permeability/PressurePermeability.hpp" diff --git a/src/coreComponents/constitutive/permeability/StrainDependentPermeability.hpp b/src/coreComponents/constitutive/permeability/StrainDependentPermeability.hpp index 1feb2c97c3a..2f9157e986d 100644 --- a/src/coreComponents/constitutive/permeability/StrainDependentPermeability.hpp +++ b/src/coreComponents/constitutive/permeability/StrainDependentPermeability.hpp @@ -38,7 +38,7 @@ class StrainDependentPermeabilityUpdate : public PermeabilityBaseUpdate arrayView3d< real64 > const & dPerm_dVolStrain, arrayView1d< real64 > const & referencePorosity, arrayView3d< real64 > const & referencePermeability, - real64 const strainDependenceConstants ) + R1Tensor const strainDependenceConstants ) : PermeabilityBaseUpdate( permeability, dPerm_dPressure ), m_dPerm_dPorosity( dPerm_dPorosity ), m_dPerm_dVolStrain( dPerm_dVolStrain ), @@ -70,8 +70,6 @@ class StrainDependentPermeabilityUpdate : public PermeabilityBaseUpdate referencePermeability[1] = m_referencePermeability[k][0][1]; referencePermeability[2] = m_referencePermeability[k][0][2]; - volStrain = m_constitutiveUpdate.getVolStrain( k, q ); - compute( m_referencePorosity[k], referencePermeability, m_strainDependenceConstants, @@ -89,14 +87,14 @@ class StrainDependentPermeabilityUpdate : public PermeabilityBaseUpdate /// dPermeability_dVolumetricStrain arrayView3d< real64 > m_dPerm_dVolStrain; - /// Volumetric strain dependence coefficients for each permeability component - R1Tensor m_strainDependenceConstants; + /// Reference porosity + arrayView1d< real64 > m_referencePorosity; /// Reference permeability arrayView3d< real64 > m_referencePermeability; - /// Reference porosity - arrayView1d< real64 > m_referencePorosity; + /// Volumetric strain dependence coefficients for each permeability component + R1Tensor m_strainDependenceConstants; }; @@ -128,10 +126,10 @@ class StrainDependentPermeability : public PermeabilityBase return KernelWrapper( m_permeability, m_dPerm_dPressure, m_dPerm_dPorosity, - m_dPerm_dVolStrain, - m_strainDependenceConstants, + m_dPerm_dVolStrain, m_referencePorosity, - m_referencePermeability ); + m_referencePermeability, + m_strainDependenceConstants); } @@ -158,14 +156,14 @@ class StrainDependentPermeability : public PermeabilityBase /// dPermeability_dVolumetricStrain array3d< real64 > m_dPerm_dVolStrain; - /// Volumetric strain dependence coefficients for each permeability component - R1Tensor m_strainDependenceConstants; + /// Reference porosity + array1d< real64 > m_referencePorosity; /// Reference permeability array3d< real64 > m_referencePermeability; - /// Reference porosity - array1d< real64 > m_referencePorosity; + /// Volumetric strain dependence coefficients for each permeability component + R1Tensor m_strainDependenceConstants; }; @@ -178,14 +176,16 @@ void StrainDependentPermeabilityUpdate::compute( real64 const & referencePorosit arraySlice1d< real64 > const & permeability, arraySlice1d< real64 > const & dPerm_dPorosity, arraySlice1d< real64 > const & dPerm_dVolStrain ) const -{ +{ + (void)dPerm_dPorosity; + real64 const por = 1 - (1 - referencePorosity) * LvArray::math::exp(-volStrain); for( localIndex i = 0; i < permeability.size(); ++i ) { - real64 const permMultiplier = LvArray::math::pow( por/referencePorosity, strainDependenceConstants[i] ); + real64 const permMultiplier = std::pow( por/referencePorosity, strainDependenceConstants[i] ); - real64 const dpermMultiplier_dVolStrain = strainDependenceConstants[i] * LvArray::math::pow( por/referencePorosity, strainDependenceConstants[i]-1 ) /referencePorosity * ((1 - referencePorosity) * LvArray::math::exp(-volStrain)); + real64 const dpermMultiplier_dVolStrain = strainDependenceConstants[i] * std::pow( por/referencePorosity, strainDependenceConstants[i]-1 ) /referencePorosity * ((1 - referencePorosity) * LvArray::math::exp(-volStrain)); permeability[i] = permMultiplier * referencePermeability[i]; dPerm_dVolStrain[i] = dpermMultiplier_dVolStrain * referencePermeability[i]; From c6ffa703409ba29eddbb7d2d42589402627114e7 Mon Sep 17 00:00:00 2001 From: Jian HUANG Date: Thu, 2 Oct 2025 18:53:39 -0500 Subject: [PATCH 08/11] add to CMakeLists --- src/coreComponents/constitutive/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/src/coreComponents/constitutive/CMakeLists.txt b/src/coreComponents/constitutive/CMakeLists.txt index c835c57966f..868927524df 100644 --- a/src/coreComponents/constitutive/CMakeLists.txt +++ b/src/coreComponents/constitutive/CMakeLists.txt @@ -147,6 +147,7 @@ set( constitutive_headers fluid/twophaseimmisciblefluid/TwoPhaseImmiscibleFluidFields.hpp permeability/CarmanKozenyPermeability.hpp permeability/ConstantPermeability.hpp + permeability/StrainDependentPermeability.hpp permeability/DamagePermeability.hpp permeability/ExponentialDecayPermeability.hpp permeability/ParallelPlatesPermeability.hpp From b172335bf6ad426cf439c3b40edf9015904b15b0 Mon Sep 17 00:00:00 2001 From: Jian HUANG Date: Thu, 2 Oct 2025 19:11:22 -0500 Subject: [PATCH 09/11] update CMakeLists --- src/coreComponents/constitutive/CMakeLists.txt | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/coreComponents/constitutive/CMakeLists.txt b/src/coreComponents/constitutive/CMakeLists.txt index 868927524df..4c00dd610d6 100644 --- a/src/coreComponents/constitutive/CMakeLists.txt +++ b/src/coreComponents/constitutive/CMakeLists.txt @@ -147,7 +147,6 @@ set( constitutive_headers fluid/twophaseimmisciblefluid/TwoPhaseImmiscibleFluidFields.hpp permeability/CarmanKozenyPermeability.hpp permeability/ConstantPermeability.hpp - permeability/StrainDependentPermeability.hpp permeability/DamagePermeability.hpp permeability/ExponentialDecayPermeability.hpp permeability/ParallelPlatesPermeability.hpp @@ -156,6 +155,7 @@ set( constitutive_headers permeability/PressurePermeability.hpp permeability/ProppantPermeability.hpp permeability/SlipDependentPermeability.hpp + permeability/StrainDependentPermeability.hpp permeability/WillisRichardsPermeability.hpp relativePermeability/BrooksCoreyBakerRelativePermeability.hpp relativePermeability/BrooksCoreyStone2RelativePermeability.hpp @@ -305,6 +305,7 @@ set( constitutive_sources permeability/PressurePermeability.cpp permeability/ProppantPermeability.cpp permeability/SlipDependentPermeability.cpp + permeability/StrainDependentPermeability.cpp permeability/WillisRichardsPermeability.cpp relativePermeability/BrooksCoreyBakerRelativePermeability.cpp relativePermeability/BrooksCoreyStone2RelativePermeability.cpp From 43515179c049a7fa367d143ad4aef0bc089a7674 Mon Sep 17 00:00:00 2001 From: Jian HUANG Date: Thu, 2 Oct 2025 20:04:39 -0500 Subject: [PATCH 10/11] handle PermeabilityFields --- .../constitutive/permeability/PermeabilityFields.hpp | 8 ++++++++ .../permeability/StrainDependentPermeability.cpp | 2 ++ .../permeability/StrainDependentPermeability.hpp | 4 ++++ 3 files changed, 14 insertions(+) diff --git a/src/coreComponents/constitutive/permeability/PermeabilityFields.hpp b/src/coreComponents/constitutive/permeability/PermeabilityFields.hpp index e1a48f17618..9ad4377123c 100644 --- a/src/coreComponents/constitutive/permeability/PermeabilityFields.hpp +++ b/src/coreComponents/constitutive/permeability/PermeabilityFields.hpp @@ -71,6 +71,14 @@ DECLARE_FIELD( dPerm_dTraction, WRITE_AND_READ, "Derivative of rock permeability with respect to the traction vector" ); +DECLARE_FIELD( dPerm_dVolStrain, + "dPerm_dVolStrain", + array3d< real64 >, + 0, + NOPLOT, + WRITE_AND_READ, + "Derivative of rock permeability with respect to volumetric strain" ); + DECLARE_FIELD( permeabilityMultiplier, "permeabilityMultiplier", array3d< real64 >, diff --git a/src/coreComponents/constitutive/permeability/StrainDependentPermeability.cpp b/src/coreComponents/constitutive/permeability/StrainDependentPermeability.cpp index e1518ffc8d4..c1d45f2f91a 100644 --- a/src/coreComponents/constitutive/permeability/StrainDependentPermeability.cpp +++ b/src/coreComponents/constitutive/permeability/StrainDependentPermeability.cpp @@ -19,6 +19,8 @@ #include "StrainDependentPermeability.hpp" +#include "constitutive/permeability/PermeabilityFields.hpp" + namespace geos { diff --git a/src/coreComponents/constitutive/permeability/StrainDependentPermeability.hpp b/src/coreComponents/constitutive/permeability/StrainDependentPermeability.hpp index 2f9157e986d..d0197de6d04 100644 --- a/src/coreComponents/constitutive/permeability/StrainDependentPermeability.hpp +++ b/src/coreComponents/constitutive/permeability/StrainDependentPermeability.hpp @@ -135,6 +135,7 @@ class StrainDependentPermeability : public PermeabilityBase struct viewKeyStruct : public PermeabilityBase::viewKeyStruct { + static constexpr char const * referencePermeabilityComponentsString() { return "referencePermeabilityComponents"; } static constexpr char const * dPerm_dPorosityString() { return "dPerm_dPorosity"; } static constexpr char const * dPerm_dVolStrainString() { return "dPerm_dVolStrain"; } static constexpr char const * strainDependenceConstantsString() { return "strainDependenceConstants"; } @@ -150,6 +151,9 @@ class StrainDependentPermeability : public PermeabilityBase private: + /// Permeability components at the reference pressure + R1Tensor m_referencePermeabilityComponents; + /// dPermeability_dPorosity array3d< real64 > m_dPerm_dPorosity; From ec5013363b31a38311a476308893bc85d3e093e3 Mon Sep 17 00:00:00 2001 From: Jian HUANG Date: Wed, 8 Oct 2025 08:37:59 -0500 Subject: [PATCH 11/11] enable m_permUpdate --- .../permeability/PermeabilityBase.hpp | 3 +-- .../StrainDependentPermeability.hpp | 10 ++++--- .../constitutive/solid/PorousSolid.hpp | 27 +++++++++++++++++++ 3 files changed, 34 insertions(+), 6 deletions(-) diff --git a/src/coreComponents/constitutive/permeability/PermeabilityBase.hpp b/src/coreComponents/constitutive/permeability/PermeabilityBase.hpp index 1a68fc3c634..0d81fb65240 100644 --- a/src/coreComponents/constitutive/permeability/PermeabilityBase.hpp +++ b/src/coreComponents/constitutive/permeability/PermeabilityBase.hpp @@ -58,11 +58,10 @@ class PermeabilityBaseUpdate GEOS_HOST_DEVICE virtual void updateFromPorosityAndStrain( localIndex const k, - localIndex const q, real64 const & volStrain, real64 const & porosity ) const { - GEOS_UNUSED_VAR( k, q, volStrain, porosity ); + GEOS_UNUSED_VAR( k, volStrain, porosity ); } GEOS_HOST_DEVICE diff --git a/src/coreComponents/constitutive/permeability/StrainDependentPermeability.hpp b/src/coreComponents/constitutive/permeability/StrainDependentPermeability.hpp index d0197de6d04..eb2d7177ee9 100644 --- a/src/coreComponents/constitutive/permeability/StrainDependentPermeability.hpp +++ b/src/coreComponents/constitutive/permeability/StrainDependentPermeability.hpp @@ -21,6 +21,7 @@ #define GEOS_CONSTITUTIVE_PERMEABILITY_STRAINDEPENDENTPERMEABILITY_HPP_ #include "constitutive/permeability/PermeabilityBase.hpp" +#include "LvArray/src/tensorOps.hpp" namespace geos @@ -58,7 +59,6 @@ class StrainDependentPermeabilityUpdate : public PermeabilityBaseUpdate GEOS_HOST_DEVICE virtual void updateFromPorosityAndStrain( localIndex const k, - localIndex const q, real64 const & volStrain, real64 const & porosity ) const override { @@ -74,9 +74,9 @@ class StrainDependentPermeabilityUpdate : public PermeabilityBaseUpdate referencePermeability, m_strainDependenceConstants, volStrain, - m_permeability[k][q], - m_dPerm_dPorosity[k][q], - m_dPerm_dVolStrain[k][q] ); + m_permeability[k][0], + m_dPerm_dPorosity[k][0], + m_dPerm_dVolStrain[k][0] ); } private: @@ -183,6 +183,8 @@ void StrainDependentPermeabilityUpdate::compute( real64 const & referencePorosit { (void)dPerm_dPorosity; + std::cout << "volStrain = " << volStrain << std::endl; + real64 const por = 1 - (1 - referencePorosity) * LvArray::math::exp(-volStrain); for( localIndex i = 0; i < permeability.size(); ++i ) diff --git a/src/coreComponents/constitutive/solid/PorousSolid.hpp b/src/coreComponents/constitutive/solid/PorousSolid.hpp index 1edd995d78b..c450c3ac02c 100644 --- a/src/coreComponents/constitutive/solid/PorousSolid.hpp +++ b/src/coreComponents/constitutive/solid/PorousSolid.hpp @@ -128,6 +128,8 @@ class PorousSolidUpdates : public CoupledSolidUpdates< SOLID_TYPE, BiotPorosity, dPorosity_dTemperature = 0.0; } + updateMatrixPermeability( k ); + // Save the derivative of solid density wrt pressure for the computation of the body force dSolidDensity_dPressure = m_porosityUpdate.dGrainDensity_dPressure( k ); } @@ -241,6 +243,31 @@ class PorousSolidUpdates : public CoupledSolidUpdates< SOLID_TYPE, BiotPorosity, m_porosityUpdate.updateThermalExpansionCoefficient( k, thermalExpansionCoefficient ); } + GEOS_HOST_DEVICE + void updateMatrixPermeability( localIndex const k ) const + { + // Use the averaged volumetric strain from all quadrature points to get the cell-centered permeability + real64 const avgStrain[6] = m_solidUpdate.m_avgStrain[k]; + + real64 const volStrain = avgStrain[0] + avgStrain[1] + avgStrain[2]; + + // Use the averaged porosity from all quadrature points to get the cell-centered permeability + integer const numQuad = m_solidUpdate.m_strain[k].size(); + + real64 avgPorosity = 0.0; + + for( localIndex q = 0; q < numQuad; ++q ) + { + avgPorosity += m_porosityUpdate.getPorosity( k, q ); + } + + avgPorosity = avgPorosity/numQuad; + + m_permUpdate.updateFromPorosityAndStrain( k, + volStrain, + avgPorosity ); + } + GEOS_HOST_DEVICE void computePorosity( localIndex const k, localIndex const q,