[PL/TRM] Extracted transport porosity evaluation to a separate model

ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/PermeabilityModel.cpp

@@ -17,17 +17,12 @@ namespace ProcessLib::ThermoRichardsMechanics
 template <int DisplacementDim>
 void PermeabilityModel<DisplacementDim>::eval(
     SpaceTimeData const& x_t, MediaData const& media_data,
-    SolidCompressibilityData const& solid_compressibility_data,
-    SaturationData const& S_L_data, BishopsData const& bishops_data,
-    BishopsData const& bishops_data_prev,
+    SaturationData const& S_L_data,
     CapillaryPressureData<DisplacementDim> const& p_cap_data,
     TemperatureData<DisplacementDim> const& T_data,
-    PorosityData const& poro_data, LiquidViscosityData const& mu_L_data,
-    TransportPorosityData& transport_poro_data,
-    TransportPorosityData const& transport_poro_data_prev,
+    LiquidViscosityData const& mu_L_data,
+    TransportPorosityData const& transport_poro_data,
     TotalStressData<DisplacementDim> const& total_stress_data,
-    StrainData<DisplacementDim> const& eps_data,
-    StrainData<DisplacementDim> const& eps_prev_data,
     EquivalentPlasticStrainData const& equiv_plast_strain_data,
     PermeabilityData<DisplacementDim>& out) const
 {
@@ -39,44 +34,7 @@ void PermeabilityModel<DisplacementDim>::eval(
     variables.liquid_saturation = S_L_data.S_L;
     variables.temperature = T_data.T;
     variables.capillary_pressure = p_cap_data.p_cap;
-    MPL::VariableArray variables_prev;
-
-    if (medium.hasProperty(MPL::PropertyType::transport_porosity))
-    {
-        static constexpr int kelvin_vector_size =
-            MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
-        using Invariants =
-            MathLib::KelvinVector::Invariants<kelvin_vector_size>;
-        // Used in
-        // MaterialLib/MPL/Properties/PermeabilityOrthotropicPowerLaw.cpp
-        variables_prev.transport_porosity = transport_poro_data_prev.phi;
-
-        // Used in
-        // MaterialLib/MPL/Properties/TransportPorosityFromMassBalance.cpp
-        variables.grain_compressibility = solid_compressibility_data.beta_SR;
-        // Set volumetric strain rate for the general case without swelling.
-        variables.volumetric_strain = Invariants::trace(eps_data.eps);
-        variables_prev.volumetric_strain = Invariants::trace(eps_prev_data.eps);
-        variables.effective_pore_pressure =
-            -bishops_data.chi_S_L * p_cap_data.p_cap;
-        variables.porosity = poro_data.phi;
-
-        // Used in MaterialLib/MPL/Properties/PorosityFromMassBalance.cpp
-        // and MaterialLib/MPL/Properties/TransportPorosityFromMassBalance.cpp
-        variables_prev.effective_pore_pressure =
-            -bishops_data_prev.chi_S_L *
-            (p_cap_data.p_cap - p_cap_data.p_cap_dot * x_t.dt);
-
-        transport_poro_data.phi =
-            medium.property(MPL::PropertyType::transport_porosity)
-                .template value<double>(variables, variables_prev, x_t.x, x_t.t,
-                                        x_t.dt);
-        variables.transport_porosity = transport_poro_data.phi;
-    }
-    else
-    {
-        variables.transport_porosity = poro_data.phi;
-    }
+    variables.transport_porosity = transport_poro_data.phi;
 
     out.k_rel = medium.property(MPL::PropertyType::relative_permeability)
                     .template value<double>(variables, x_t.x, x_t.t, x_t.dt);

ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/PermeabilityModel.h

@@ -13,8 +13,8 @@
 #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/EquivalentPlasticStrainData.h"
 #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/LiquidViscosity.h"
 #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/PermeabilityData.h"
-#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/Porosity.h"
 #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/TotalStressData.h"
+#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/TransportPorosity.h"
 
 namespace ProcessLib::ThermoRichardsMechanics
 {
@@ -22,19 +22,12 @@ template <int DisplacementDim>
 struct PermeabilityModel
 {
     void eval(SpaceTimeData const& x_t, MediaData const& media_data,
-              SolidCompressibilityData const& solid_compressibility_data,
-              SaturationData const& S_L_data, BishopsData const& bishops_data,
-              BishopsData const& bishops_data_prev,
+              SaturationData const& S_L_data,
               CapillaryPressureData<DisplacementDim> const& p_cap_data,
               TemperatureData<DisplacementDim> const& T_data,
-              PorosityData const& poro_data,
               LiquidViscosityData const& mu_L_data,
-              // TODO evaluate transport porosity evolution separately
-              TransportPorosityData& transport_poro_data,
-              TransportPorosityData const& transport_poro_data_prev,
+              TransportPorosityData const& transport_poro_data,
               TotalStressData<DisplacementDim> const& total_stress_data,
-              StrainData<DisplacementDim> const& eps_data,
-              StrainData<DisplacementDim> const& eps_prev_data,
               EquivalentPlasticStrainData const& equiv_plast_strain_data,
               PermeabilityData<DisplacementDim>& out) const;
 };

ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/Porosity.h

@@ -27,17 +27,6 @@ struct PorosityData
     }
 };
 
-struct TransportPorosityData
-{
-    double phi;
-
-    static auto reflect()
-    {
-        return ProcessLib::Reflection::reflectWithName(
-            "transport_porosity", &TransportPorosityData::phi);
-    }
-};
-
 template <int DisplacementDim>
 struct PorosityModel
 {

ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/TransportPorosity.cpp

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2022, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "TransportPorosity.h"
+
+namespace ProcessLib::ThermoRichardsMechanics
+{
+
+template <int DisplacementDim>
+
+void TransportPorosityModel<DisplacementDim>::eval(
+    SpaceTimeData const& x_t, MediaData const& media_data,
+    SolidCompressibilityData const& solid_compressibility_data,
+    BishopsData const& bishops_data, BishopsData const& bishops_data_prev,
+    CapillaryPressureData<DisplacementDim> const& p_cap_data,
+    PorosityData const& poro_data, StrainData<DisplacementDim> const& eps_data,
+    StrainData<DisplacementDim> const& eps_prev_data,
+    TransportPorosityData const& transport_poro_data_prev,
+    TransportPorosityData& transport_poro_data) const
+{
+    namespace MPL = MaterialPropertyLib;
+
+    auto const& medium = media_data.medium;
+
+    if (!medium.hasProperty(MPL::PropertyType::transport_porosity))
+    {
+        transport_poro_data.phi = poro_data.phi;
+        return;
+    }
+
+    MPL::VariableArray variables;
+    MPL::VariableArray variables_prev;
+
+    static constexpr int kelvin_vector_size =
+        MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
+    using Invariants = MathLib::KelvinVector::Invariants<kelvin_vector_size>;
+    // Used in
+    // MaterialLib/MPL/Properties/PermeabilityOrthotropicPowerLaw.cpp
+    variables_prev.transport_porosity = transport_poro_data_prev.phi;
+
+    // Used in
+    // MaterialLib/MPL/Properties/TransportPorosityFromMassBalance.cpp
+    variables.grain_compressibility = solid_compressibility_data.beta_SR;
+    // Set volumetric strain rate for the general case without swelling.
+    variables.volumetric_strain = Invariants::trace(eps_data.eps);
+    variables_prev.volumetric_strain = Invariants::trace(eps_prev_data.eps);
+    variables.effective_pore_pressure =
+        -bishops_data.chi_S_L * p_cap_data.p_cap;
+    variables.porosity = poro_data.phi;
+
+    // Used in MaterialLib/MPL/Properties/PorosityFromMassBalance.cpp
+    // and MaterialLib/MPL/Properties/TransportPorosityFromMassBalance.cpp
+    variables_prev.effective_pore_pressure =
+        -bishops_data_prev.chi_S_L *
+        (p_cap_data.p_cap - p_cap_data.p_cap_dot * x_t.dt);
+
+    transport_poro_data.phi =
+        medium.property(MPL::PropertyType::transport_porosity)
+            .template value<double>(variables, variables_prev, x_t.x, x_t.t,
+                                    x_t.dt);
+}
+
+template struct TransportPorosityModel<2>;
+template struct TransportPorosityModel<3>;
+
+}  // namespace ProcessLib::ThermoRichardsMechanics

ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/TransportPorosity.h

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2022, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
 #pragma once
+
+#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/Bishops.h"
+#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/Porosity.h"
+#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/SolidCompressibility.h"
+
+namespace ProcessLib::ThermoRichardsMechanics
+{
+
+struct TransportPorosityData
+{
+    double phi;
+
+    static auto reflect()
+    {
+        return ProcessLib::Reflection::reflectWithName(
+            "transport_porosity", &TransportPorosityData::phi);
+    }
+};
+
+template <int DisplacementDim>
+struct TransportPorosityModel
+{
+    void eval(SpaceTimeData const& x_t, MediaData const& media_data,
+              SolidCompressibilityData const& solid_compressibility_data,
+              BishopsData const& bishops_data,
+              BishopsData const& bishops_data_prev,
+              CapillaryPressureData<DisplacementDim> const& p_cap_data,
+              PorosityData const& poro_data,
+              StrainData<DisplacementDim> const& eps_data,
+              StrainData<DisplacementDim> const& eps_prev_data,
+              TransportPorosityData const& transport_poro_data_prev,
+              TransportPorosityData& transport_poro_data) const;
+};
+
+extern template struct TransportPorosityModel<2>;
+extern template struct TransportPorosityModel<3>;
+}  // namespace ProcessLib::ThermoRichardsMechanics

ProcessLib/ThermoRichardsMechanics/ConstitutiveOriginal/ConstitutiveData.h

@@ -24,6 +24,7 @@
 #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/TRMStorage.h"
 #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/TRMVaporDiffusion.h"
 #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/TotalStressData.h"
+#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/TransportPorosity.h"
 #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveOriginal/SolidMechanics.h"
 
 namespace ProcessLib::ThermoRichardsMechanics::ConstitutiveOriginal

ProcessLib/ThermoRichardsMechanics/ConstitutiveOriginal/ConstitutiveModels.h

@@ -18,9 +18,11 @@
 #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/Gravity.h"
 #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/PermeabilityModel.h"
 #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/SolidDensity.h"
+#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/Swelling.h"
 #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/TRMHeatStorageAndFlux.h"
 #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/TRMStorage.h"
 #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/TRMVaporDiffusion.h"
+#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveOriginal/SolidMechanics.h"
 
 namespace ProcessLib::ThermoRichardsMechanics::ConstitutiveOriginal
 {
@@ -50,6 +52,7 @@ struct ConstitutiveModels
     SaturationModel<DisplacementDim> S_L_model;
     BishopsModel bishops_model;
     PorosityModel<DisplacementDim> poro_model;
+    TransportPorosityModel<DisplacementDim> transport_poro_model;
     SwellingModel<DisplacementDim> swelling_model;
     SolidThermalExpansionModel<DisplacementDim> s_therm_exp_model;
     SolidMechanicsModel<DisplacementDim> s_mech_model;

ProcessLib/ThermoRichardsMechanics/ConstitutiveOriginal/ConstitutiveSetting.cpp

@@ -123,14 +123,17 @@ void ConstitutiveSetting<DisplacementDim>::eval(
 
     models.mu_L_model.eval(x_t, media_data, T_data, mu_L_data);
 
-    models.perm_model.eval(
-        x_t, media_data, solid_compressibility_data, S_L_data, bishops_data,
-        bishops_data_prev, p_cap_data, T_data, poro_data, mu_L_data,
-        state.transport_poro_data, prev_state.transport_poro_data,
-        cd.total_stress_data, state.eps_data,
+    models.transport_poro_model.eval(
+        x_t, media_data, solid_compressibility_data, bishops_data,
+        bishops_data_prev, p_cap_data, poro_data, state.eps_data,
         StrainData<DisplacementDim>{
             eps_prev_arg} /* TODO why not eqU.eps_prev? */,
-        tmp.equiv_plast_strain_data, perm_data);
+        prev_state.transport_poro_data, state.transport_poro_data);
+
+    models.perm_model.eval(x_t, media_data, S_L_data, p_cap_data, T_data,
+                           mu_L_data, state.transport_poro_data,
+                           cd.total_stress_data, tmp.equiv_plast_strain_data,
+                           perm_data);
 
     models.th_osmosis_model.eval(x_t, media_data, T_data, rho_L_data,
                                  cd.th_osmosis_data);