From b0b70b352ac43b6e0e5722e4fc2c560a37533810 Mon Sep 17 00:00:00 2001
From: Christoph Lehmann <christoph.lehmann@ufz.de>
Date: Mon, 11 Sep 2023 15:14:18 +0200
Subject: [PATCH] [PL/TRM] Enabled permeability output

---
 .../ConstitutiveCommon/DarcyLaw.cpp              |  3 ++-
 .../ConstitutiveCommon/DarcyLaw.h                |  1 +
 .../ConstitutiveCommon/EqP.cpp                   |  7 ++++---
 .../ConstitutiveCommon/EqP.h                     |  2 ++
 .../ConstitutiveCommon/PermeabilityData.h        | 12 +++++++++++-
 .../ConstitutiveCommon/PermeabilityModel.cpp     |  3 +--
 .../ConstitutiveCommon/PermeabilityModel.h       |  2 --
 .../ConstitutiveCommon/TRMHeatStorageAndFlux.cpp |  7 ++++---
 .../ConstitutiveCommon/TRMHeatStorageAndFlux.h   |  1 +
 .../ConstitutiveData.h                           |  5 +++--
 .../ConstitutiveSetting.cpp                      | 16 ++++++++--------
 .../ConstitutiveData.h                           |  5 +++--
 .../ConstitutiveSetting.cpp                      | 16 ++++++++--------
 13 files changed, 48 insertions(+), 32 deletions(-)

diff --git a/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/DarcyLaw.cpp b/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/DarcyLaw.cpp
index 5c89ff49e5e..a56725c47a4 100644
--- a/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/DarcyLaw.cpp
+++ b/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/DarcyLaw.cpp
@@ -18,11 +18,12 @@ template <int DisplacementDim>
 void DarcyLawModel<DisplacementDim>::eval(
     CapillaryPressureData<DisplacementDim> const& p_cap_data,
     LiquidDensityData const& rho_L_data,
+    LiquidViscosityData const& mu_L_data,
     PermeabilityData<DisplacementDim> const& perm_data,
     ThermoOsmosisData<DisplacementDim> const& th_osmosis_data,
     DarcyLawData<DisplacementDim>& out) const
 {
-    out.v_darcy = perm_data.Ki_over_mu *
+    out.v_darcy = perm_data.Ki / mu_L_data.viscosity *
                       (perm_data.k_rel *
                        (p_cap_data.grad_p_cap + rho_L_data.rho_LR * b_)) +
                   th_osmosis_data.seepage_velocity_contribution;
diff --git a/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/DarcyLaw.h b/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/DarcyLaw.h
index d2384690f08..fc64acaa541 100644
--- a/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/DarcyLaw.h
+++ b/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/DarcyLaw.h
@@ -42,6 +42,7 @@ struct DarcyLawModel
 
     void eval(CapillaryPressureData<DisplacementDim> const& p_cap_data,
               LiquidDensityData const& rho_L_data,
+              LiquidViscosityData const& mu_L_data,
               PermeabilityData<DisplacementDim> const& perm_data,
               ThermoOsmosisData<DisplacementDim> const& th_osmosis_data,
               DarcyLawData<DisplacementDim>& out) const;
diff --git a/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/EqP.cpp b/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/EqP.cpp
index 802f7813cf2..b8d7a09315f 100644
--- a/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/EqP.cpp
+++ b/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/EqP.cpp
@@ -20,6 +20,7 @@ void EqPModel<DisplacementDim>::eval(
     SaturationDataDeriv const& dS_L_data,
     BiotData const& biot_data,
     LiquidDensityData const& rho_L_data,
+    LiquidViscosityData const& mu_L_data,
     PermeabilityData<DisplacementDim> const& perm_data,
     FluidThermalExpansionData const& f_therm_exp_data,
     TRMVaporDiffusionData<DisplacementDim> const& vap_data,
@@ -28,14 +29,14 @@ void EqPModel<DisplacementDim>::eval(
 {
     out.M_pu_X_BTI2N = S_L_data.S_L * rho_L_data.rho_LR * biot_data();
 
-    out.K_pp_Laplace =
-        perm_data.k_rel * rho_L_data.rho_LR * perm_data.Ki_over_mu;
+    out.K_pp_Laplace = perm_data.k_rel * rho_L_data.rho_LR * perm_data.Ki /
+                       mu_L_data.viscosity;
 
     out.J_pp_X_BTI2NT_u_dot_N =
         -rho_L_data.rho_LR * dS_L_data.dS_L_dp_cap * biot_data();
 
     out.J_pp_dNT_V_N =
-        perm_data.Ki_over_mu *
+        perm_data.Ki / mu_L_data.viscosity *
         (rho_L_data.rho_LR * perm_data.dk_rel_dS_L * dS_L_data.dS_L_dp_cap *
          (p_cap_data.grad_p_cap + rho_L_data.rho_LR * b_));
 
diff --git a/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/EqP.h b/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/EqP.h
index 6ed67d921a9..6b4d4192069 100644
--- a/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/EqP.h
+++ b/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/EqP.h
@@ -13,6 +13,7 @@
 #include "Biot.h"
 #include "FluidThermalExpansion.h"
 #include "LiquidDensity.h"
+#include "LiquidViscosity.h"
 #include "PermeabilityData.h"
 #include "Saturation.h"
 #include "TRMStorage.h"
@@ -51,6 +52,7 @@ struct EqPModel
               SaturationDataDeriv const& dS_L_data,
               BiotData const& biot_data,
               LiquidDensityData const& rho_L_data,
+              LiquidViscosityData const& mu_L_data,
               PermeabilityData<DisplacementDim> const& perm_data,
               FluidThermalExpansionData const& f_therm_exp_data,
               TRMVaporDiffusionData<DisplacementDim> const& vap_data,
diff --git a/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/PermeabilityData.h b/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/PermeabilityData.h
index 53a09055e30..c5edff77036 100644
--- a/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/PermeabilityData.h
+++ b/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/PermeabilityData.h
@@ -19,6 +19,16 @@ struct PermeabilityData
 {
     double k_rel;
     double dk_rel_dS_L;
-    GlobalDimMatrix<DisplacementDim> Ki_over_mu;
+    GlobalDimMatrix<DisplacementDim> Ki;
+
+    static auto reflect()
+    {
+        using Self = PermeabilityData<DisplacementDim>;
+        namespace R = ProcessLib::Reflection;
+
+        return std::tuple{
+            R::makeReflectionData("intrinsic_permeability", &Self::Ki),
+            R::makeReflectionData("relative_permeability", &Self::k_rel)};
+    }
 };
 }  // namespace ProcessLib::ThermoRichardsMechanics
diff --git a/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/PermeabilityModel.cpp b/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/PermeabilityModel.cpp
index 00772d4b0de..c5a1300a1d3 100644
--- a/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/PermeabilityModel.cpp
+++ b/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/PermeabilityModel.cpp
@@ -20,7 +20,6 @@ void PermeabilityModel<DisplacementDim>::eval(
     SaturationData const& S_L_data,
     CapillaryPressureData<DisplacementDim> const& p_cap_data,
     TemperatureData<DisplacementDim> const& T_data,
-    LiquidViscosityData const& mu_L_data,
     TransportPorosityData const& transport_poro_data,
     TotalStressData<DisplacementDim> const& total_stress_data,
     StrainData<DisplacementDim> const& eps_data,
@@ -68,7 +67,7 @@ void PermeabilityModel<DisplacementDim>::eval(
         medium.property(MPL::PropertyType::permeability)
             .value(variables, x_t.x, x_t.t, x_t.dt));
 
-    out.Ki_over_mu = K_intrinsic / mu_L_data.viscosity;
+    out.Ki = K_intrinsic;
 }
 
 template struct PermeabilityModel<2>;
diff --git a/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/PermeabilityModel.h b/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/PermeabilityModel.h
index acb1ee0800b..be6272b7c23 100644
--- a/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/PermeabilityModel.h
+++ b/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/PermeabilityModel.h
@@ -12,7 +12,6 @@
 
 #include "ProcessLib/ConstitutiveRelations/StrainData.h"
 #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/EquivalentPlasticStrainData.h"
-#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/LiquidViscosity.h"
 #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/PermeabilityData.h"
 #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/TotalStressData.h"
 #include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/TransportPorosity.h"
@@ -26,7 +25,6 @@ struct PermeabilityModel
               SaturationData const& S_L_data,
               CapillaryPressureData<DisplacementDim> const& p_cap_data,
               TemperatureData<DisplacementDim> const& T_data,
-              LiquidViscosityData const& mu_L_data,
               TransportPorosityData const& transport_poro_data,
               TotalStressData<DisplacementDim> const& total_stress_data,
               StrainData<DisplacementDim> const& eps_data,
diff --git a/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/TRMHeatStorageAndFlux.cpp b/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/TRMHeatStorageAndFlux.cpp
index d155c24b3af..0f81bea6959 100644
--- a/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/TRMHeatStorageAndFlux.cpp
+++ b/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/TRMHeatStorageAndFlux.cpp
@@ -19,7 +19,7 @@ void TRMHeatStorageAndFluxModel<DisplacementDim>::eval(
     SpaceTimeData const& x_t, MediaData const& media_data,
     LiquidDensityData const& rho_L_data, SolidDensityData const& rho_S_data,
     SaturationData const& S_L_data, SaturationDataDeriv const& dS_L_data,
-    PorosityData const& poro_data,
+    PorosityData const& poro_data, LiquidViscosityData const& mu_L_data,
     PermeabilityData<DisplacementDim> const& perm,
     TemperatureData<DisplacementDim> const& T_data,
     DarcyLawData<DisplacementDim> const& darcy_data,
@@ -68,8 +68,9 @@ void TRMHeatStorageAndFluxModel<DisplacementDim>::eval(
     //
     // temperature equation, pressure part
     //
-    out.K_Tp_NT_V_dN = -volumetric_heat_capacity_liquid * perm.k_rel *
-                       (perm.Ki_over_mu.transpose() * T_data.grad_T);
+    out.K_Tp_NT_V_dN = -volumetric_heat_capacity_liquid * perm.k_rel /
+                       mu_L_data.viscosity *
+                       (perm.Ki.transpose() * T_data.grad_T);
     out.K_Tp_X_NTN = -volumetric_heat_capacity_liquid *
                      darcy_data.v_darcy.dot(T_data.grad_T) / perm.k_rel *
                      perm.dk_rel_dS_L * dS_L_data.dS_L_dp_cap;
diff --git a/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/TRMHeatStorageAndFlux.h b/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/TRMHeatStorageAndFlux.h
index 09649b1239f..4602e0f5b16 100644
--- a/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/TRMHeatStorageAndFlux.h
+++ b/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/TRMHeatStorageAndFlux.h
@@ -38,6 +38,7 @@ struct TRMHeatStorageAndFluxModel
               SaturationData const& S_L_data,
               SaturationDataDeriv const& dS_L_data,
               PorosityData const& poro_data,
+              LiquidViscosityData const& mu_L_data,
               PermeabilityData<DisplacementDim> const& perm,
               TemperatureData<DisplacementDim> const& T_data,
               DarcyLawData<DisplacementDim> const& darcy_data,
diff --git a/ProcessLib/ThermoRichardsMechanics/ConstitutiveStressSaturation_StrainPressureTemperature/ConstitutiveData.h b/ProcessLib/ThermoRichardsMechanics/ConstitutiveStressSaturation_StrainPressureTemperature/ConstitutiveData.h
index 2c6613177a5..d8260cf5be2 100644
--- a/ProcessLib/ThermoRichardsMechanics/ConstitutiveStressSaturation_StrainPressureTemperature/ConstitutiveData.h
+++ b/ProcessLib/ThermoRichardsMechanics/ConstitutiveStressSaturation_StrainPressureTemperature/ConstitutiveData.h
@@ -95,6 +95,7 @@ struct OutputData
     LiquidDensityData rho_L_data;
     LiquidViscosityData mu_L_data;
     SolidDensityData rho_S_data;
+    PermeabilityData<DisplacementDim> perm_data;
 
     static auto reflect()
     {
@@ -103,7 +104,8 @@ struct OutputData
         return Reflection::reflectWithoutName(&Self::darcy_data,
                                               &Self::rho_L_data,
                                               &Self::mu_L_data,
-                                              &Self::rho_S_data);
+                                              &Self::rho_S_data,
+                                              &Self::perm_data);
     }
 };
 
@@ -134,7 +136,6 @@ struct ConstitutiveTempData
     // TODO why not usual state tracking for that?
     PrevState<BishopsData> bishops_data_prev;
     SolidThermalExpansionData<DisplacementDim> s_therm_exp_data;
-    PermeabilityData<DisplacementDim> perm_data;
     FluidThermalExpansionData f_therm_exp_data;
     EquivalentPlasticStrainData equiv_plast_strain_data;
 };
diff --git a/ProcessLib/ThermoRichardsMechanics/ConstitutiveStressSaturation_StrainPressureTemperature/ConstitutiveSetting.cpp b/ProcessLib/ThermoRichardsMechanics/ConstitutiveStressSaturation_StrainPressureTemperature/ConstitutiveSetting.cpp
index b312ae8cd3b..dfa134a890c 100644
--- a/ProcessLib/ThermoRichardsMechanics/ConstitutiveStressSaturation_StrainPressureTemperature/ConstitutiveSetting.cpp
+++ b/ProcessLib/ThermoRichardsMechanics/ConstitutiveStressSaturation_StrainPressureTemperature/ConstitutiveSetting.cpp
@@ -43,7 +43,7 @@ void ConstitutiveSetting<DisplacementDim>::eval(
     auto& rho_L_data = out.rho_L_data;
     auto& rho_S_data = out.rho_S_data;
     auto& mu_L_data = out.mu_L_data;
-    auto& perm_data = tmp.perm_data;
+    auto& perm_data = out.perm_data;
     auto& darcy_data = out.darcy_data;
     auto& f_therm_exp_data = tmp.f_therm_exp_data;
 
@@ -126,22 +126,22 @@ void ConstitutiveSetting<DisplacementDim>::eval(
 
     assertEvalArgsUnique(models.perm_model);
     models.perm_model.eval(x_t, media_data, S_L_data, p_cap_data, T_data,
-                           mu_L_data, state.transport_poro_data,
-                           state.total_stress_data, state.eps_data,
-                           tmp.equiv_plast_strain_data, perm_data);
+                           state.transport_poro_data, state.total_stress_data,
+                           state.eps_data, tmp.equiv_plast_strain_data,
+                           perm_data);
 
     assertEvalArgsUnique(models.th_osmosis_model);
     models.th_osmosis_model.eval(x_t, media_data, T_data, rho_L_data,
                                  cd.th_osmosis_data);
 
     assertEvalArgsUnique(models.darcy_model);
-    models.darcy_model.eval(p_cap_data, rho_L_data, perm_data,
+    models.darcy_model.eval(p_cap_data, rho_L_data, mu_L_data, perm_data,
                             cd.th_osmosis_data, darcy_data);
 
     assertEvalArgsUnique(models.heat_storage_and_flux_model);
     models.heat_storage_and_flux_model.eval(
         x_t, media_data, rho_L_data, rho_S_data, S_L_data, dS_L_data, poro_data,
-        perm_data, T_data, darcy_data, heat_data);
+        mu_L_data, perm_data, T_data, darcy_data, heat_data);
 
     assertEvalArgsUnique(models.vapor_diffusion_model);
     models.vapor_diffusion_model.eval(x_t, media_data, rho_L_data, S_L_data,
@@ -166,8 +166,8 @@ void ConstitutiveSetting<DisplacementDim>::eval(
 
     assertEvalArgsUnique(models.eq_p_model);
     models.eq_p_model.eval(p_cap_data, T_data, S_L_data, dS_L_data, biot_data,
-                           rho_L_data, perm_data, f_therm_exp_data, vap_data,
-                           storage_data, cd.eq_p_data);
+                           rho_L_data, mu_L_data, perm_data, f_therm_exp_data,
+                           vap_data, storage_data, cd.eq_p_data);
 
     assertEvalArgsUnique(models.eq_T_model);
     models.eq_T_model.eval(heat_data, vap_data, cd.eq_T_data);
diff --git a/ProcessLib/ThermoRichardsMechanics/ConstitutiveStress_StrainTemperature/ConstitutiveData.h b/ProcessLib/ThermoRichardsMechanics/ConstitutiveStress_StrainTemperature/ConstitutiveData.h
index 36f0f35f11d..bf5485816f2 100644
--- a/ProcessLib/ThermoRichardsMechanics/ConstitutiveStress_StrainTemperature/ConstitutiveData.h
+++ b/ProcessLib/ThermoRichardsMechanics/ConstitutiveStress_StrainTemperature/ConstitutiveData.h
@@ -87,6 +87,7 @@ struct OutputData
     LiquidDensityData rho_L_data;
     LiquidViscosityData mu_L_data;
     SolidDensityData rho_S_data;
+    PermeabilityData<DisplacementDim> perm_data;
 
     static auto reflect()
     {
@@ -95,7 +96,8 @@ struct OutputData
         return Reflection::reflectWithoutName(&Self::darcy_data,
                                               &Self::rho_L_data,
                                               &Self::mu_L_data,
-                                              &Self::rho_S_data);
+                                              &Self::rho_S_data,
+                                              &Self::perm_data);
     }
 };
 
@@ -128,7 +130,6 @@ struct ConstitutiveTempData
     // TODO why not usual state tracking for that?
     PrevState<BishopsData> bishops_data_prev;
     SolidThermalExpansionData<DisplacementDim> s_therm_exp_data;
-    PermeabilityData<DisplacementDim> perm_data;
     FluidThermalExpansionData f_therm_exp_data;
     EquivalentPlasticStrainData equiv_plast_strain_data;
 };
diff --git a/ProcessLib/ThermoRichardsMechanics/ConstitutiveStress_StrainTemperature/ConstitutiveSetting.cpp b/ProcessLib/ThermoRichardsMechanics/ConstitutiveStress_StrainTemperature/ConstitutiveSetting.cpp
index 5b7f044187d..ca99347b666 100644
--- a/ProcessLib/ThermoRichardsMechanics/ConstitutiveStress_StrainTemperature/ConstitutiveSetting.cpp
+++ b/ProcessLib/ThermoRichardsMechanics/ConstitutiveStress_StrainTemperature/ConstitutiveSetting.cpp
@@ -43,7 +43,7 @@ void ConstitutiveSetting<DisplacementDim>::eval(
     auto& rho_L_data = out.rho_L_data;
     auto& rho_S_data = out.rho_S_data;
     auto& mu_L_data = out.mu_L_data;
-    auto& perm_data = tmp.perm_data;
+    auto& perm_data = out.perm_data;
     auto& darcy_data = out.darcy_data;
     auto& f_therm_exp_data = tmp.f_therm_exp_data;
 
@@ -148,22 +148,22 @@ void ConstitutiveSetting<DisplacementDim>::eval(
 
     assertEvalArgsUnique(models.perm_model);
     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, state.eps_data,
-                           tmp.equiv_plast_strain_data, perm_data);
+                           state.transport_poro_data, cd.total_stress_data,
+                           state.eps_data, tmp.equiv_plast_strain_data,
+                           perm_data);
 
     assertEvalArgsUnique(models.th_osmosis_model);
     models.th_osmosis_model.eval(x_t, media_data, T_data, rho_L_data,
                                  cd.th_osmosis_data);
 
     assertEvalArgsUnique(models.darcy_model);
-    models.darcy_model.eval(p_cap_data, rho_L_data, perm_data,
+    models.darcy_model.eval(p_cap_data, rho_L_data, mu_L_data, perm_data,
                             cd.th_osmosis_data, darcy_data);
 
     assertEvalArgsUnique(models.heat_storage_and_flux_model);
     models.heat_storage_and_flux_model.eval(
         x_t, media_data, rho_L_data, rho_S_data, S_L_data, dS_L_data, poro_data,
-        perm_data, T_data, darcy_data, heat_data);
+        mu_L_data, perm_data, T_data, darcy_data, heat_data);
 
     assertEvalArgsUnique(models.vapor_diffusion_model);
     models.vapor_diffusion_model.eval(x_t, media_data, rho_L_data, S_L_data,
@@ -182,8 +182,8 @@ void ConstitutiveSetting<DisplacementDim>::eval(
 
     assertEvalArgsUnique(models.eq_p_model);
     models.eq_p_model.eval(p_cap_data, T_data, S_L_data, dS_L_data, biot_data,
-                           rho_L_data, perm_data, f_therm_exp_data, vap_data,
-                           storage_data, cd.eq_p_data);
+                           rho_L_data, mu_L_data, perm_data, f_therm_exp_data,
+                           vap_data, storage_data, cd.eq_p_data);
 
     assertEvalArgsUnique(models.eq_T_model);
     models.eq_T_model.eval(heat_data, vap_data, cd.eq_T_data);
-- 
GitLab