Merge branch 'Th2m-CR-FinishMPL_NaN_CleanupIPData' into 'master'

TH2M; Finalize extraction of MPL evaluations into constitutive relation models

See merge request ogs/ogs!4976

MaterialLib/MPL/Properties/EffectiveThermalConductivityPorosityMixing.h

@@ -39,9 +39,9 @@ class EffectiveThermalConductivityPorosityMixing final : public Property
                            double const dt) const override;
     PropertyDataType dValue(VariableArray const& variable_array,
                             Variable const variable,
-                            ParameterLib::SpatialPosition const& /*pos*/,
-                            double const /*t*/,
-                            double const /*dt*/) const override;
+                            ParameterLib::SpatialPosition const& pos,
+                            double const t,
+                            double const dt) const override;
 
 private:
     ParameterLib::CoordinateSystem const* const local_coordinate_system_;

ProcessLib/TH2M/ConstitutiveRelations/Base.h

@@ -32,6 +32,9 @@ template <int DisplacementDim>
 using GlobalDimMatrix =
     Eigen::Matrix<double, DisplacementDim, DisplacementDim, Eigen::RowMajor>;
 
+template <int DisplacementDim>
+using GlobalDimVector = Eigen::Vector<double, DisplacementDim>;
+
 struct MediaData
 {
     explicit MediaData(MaterialPropertyLib::Medium const& medium)
@@ -47,8 +50,8 @@ struct MediaData
 
 struct TemperatureData
 {
-    double T;
-    double T_prev;
+    double T = nan;
+    double T_prev = nan;
 };
 
 using ReferenceTemperatureData =

ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h

@@ -12,10 +12,13 @@
 #include "Biot.h"
 #include "Bishops.h"
 #include "ConstitutiveDensity.h"
+#include "DarcyVelocity.h"
+#include "DiffusionVelocity.h"
 #include "ElasticTangentStiffnessData.h"
 #include "Enthalpy.h"
 #include "EquivalentPlasticStrainData.h"
 #include "FluidDensity.h"
+#include "InternalEnergy.h"
 #include "MassMoleFractions.h"
 #include "MechanicalStrain.h"
 #include "PermeabilityData.h"
@@ -28,9 +31,11 @@
 #include "Saturation.h"
 #include "SolidCompressibility.h"
 #include "SolidDensity.h"
+#include "SolidHeatCapacity.h"
 #include "SolidMechanics.h"
 #include "SolidThermalExpansion.h"
 #include "Swelling.h"
+#include "ThermalConductivity.h"
 #include "TotalStress.h"
 #include "VapourPartialPressure.h"
 #include "Viscosity.h"
@@ -50,6 +55,7 @@ struct StatefulData
     MechanicalStrainData<DisplacementDim> mechanical_strain_data;
     PureLiquidDensityData rho_W_LR;
     ConstituentDensityData constituent_density_data;
+    InternalEnergyData internal_energy_data;
 
     static auto reflect()
     {
@@ -70,6 +76,7 @@ struct StatefulDataPrev
     PrevState<MechanicalStrainData<DisplacementDim>> mechanical_strain_data;
     PrevState<PureLiquidDensityData> rho_W_LR;
     PrevState<ConstituentDensityData> constituent_density_data;
+    PrevState<InternalEnergyData> internal_energy_data;
 
     StatefulDataPrev<DisplacementDim>& operator=(
         StatefulData<DisplacementDim> const& state)
@@ -80,6 +87,7 @@ struct StatefulDataPrev
         mechanical_strain_data = state.mechanical_strain_data;
         rho_W_LR = state.rho_W_LR;
         constituent_density_data = state.constituent_density_data;
+        internal_energy_data = state.internal_energy_data;
 
         return *this;
     }
@@ -97,6 +105,8 @@ struct OutputData
     VapourPartialPressureData vapour_pressure_data;
     PorosityData porosity_data;
     SolidDensityData solid_density_data;
+    DiffusionVelocityData<DisplacementDim> diffusion_velocity_data;
+    DarcyVelocityData<DisplacementDim> darcy_velocity_data;
 
     static auto reflect()
     {
@@ -109,7 +119,9 @@ struct OutputData
                                               &Self::fluid_density_data,
                                               &Self::vapour_pressure_data,
                                               &Self::porosity_data,
-                                              &Self::solid_density_data);
+                                              &Self::solid_density_data,
+                                              &Self::diffusion_velocity_data,
+                                              &Self::darcy_velocity_data);
     }
 };
 
@@ -138,14 +150,17 @@ struct ConstitutiveTempData
     PhaseTransitionData phase_transition_data;
     PorosityDerivativeData porosity_d_data;
     SolidDensityDerivativeData solid_density_d_data;
+    SolidHeatCapacityData solid_heat_capacity_data;
+    ThermalConductivityData<DisplacementDim> thermal_conductivity_data;
+    EffectiveVolumetricEnthalpy effective_volumetric_enthalpy_data;
+    EffectiveVolumetricEnthalpyDerivatives effective_volumetric_enthalpy_d_data;
+    EffectiveVolumetricInternalEnergyDerivatives
+        effective_volumetric_internal_energy_d_data;
 
     using DisplacementDimVector = Eigen::Matrix<double, DisplacementDim, 1>;
     using DisplacementDimMatrix =
         Eigen::Matrix<double, DisplacementDim, DisplacementDim>;
 
-    DisplacementDimMatrix dlambda_dp_GR;
-    DisplacementDimMatrix dlambda_dp_cap;
-    DisplacementDimMatrix dlambda_dT;
     DisplacementDimVector drho_GR_h_w_eff_dp_GR_Npart;
     DisplacementDimMatrix drho_GR_h_w_eff_dp_GR_gradNpart;
     DisplacementDimVector drho_LR_h_w_eff_dp_cap_Npart;
@@ -174,14 +189,6 @@ struct ConstitutiveTempData
     DisplacementDimMatrix dadvection_C_dp_cap;
     DisplacementDimMatrix dk_over_mu_G_dp_cap;
     DisplacementDimMatrix dk_over_mu_L_dp_cap;
-    double drho_u_eff_dT = std::numeric_limits<double>::quiet_NaN();
-    double drho_u_eff_dp_GR = std::numeric_limits<double>::quiet_NaN();
-    double drho_u_eff_dp_cap = std::numeric_limits<double>::quiet_NaN();
-    double drho_h_eff_dT = std::numeric_limits<double>::quiet_NaN();
-    double drho_h_eff_dp_GR = std::numeric_limits<double>::quiet_NaN();
-    double drho_h_eff_dp_cap = std::numeric_limits<double>::quiet_NaN();
-    double dh_G_dT = std::numeric_limits<double>::quiet_NaN();
-    double dh_L_dT = std::numeric_limits<double>::quiet_NaN();
     double dfC_4_MCpG_dp_GR = std::numeric_limits<double>::quiet_NaN();
     double dfC_4_MCpG_dT = std::numeric_limits<double>::quiet_NaN();
     double dfC_4_MCT_dT = std::numeric_limits<double>::quiet_NaN();

ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h

@@ -20,9 +20,11 @@
 #include "Saturation.h"
 #include "SolidCompressibility.h"
 #include "SolidDensity.h"
+#include "SolidHeatCapacity.h"
 #include "SolidMechanics.h"
 #include "SolidThermalExpansion.h"
 #include "Swelling.h"
+#include "ThermalConductivity.h"
 #include "TotalStress.h"
 #include "Viscosity.h"
 
@@ -60,6 +62,7 @@ struct ConstitutiveModels
     PermeabilityModel<DisplacementDim> permeability_model;
     PureLiquidDensityModel pure_liquid_density_model;
     PhaseTransitionModel const& phase_transition_model;
+    ViscosityModel viscosity_model;
 #ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
     PorosityModelNonConstantSolidPhaseVolumeFraction<DisplacementDim>
         porosity_model;
@@ -69,6 +72,8 @@ struct ConstitutiveModels
     PorosityModel porosity_model;
     SolidDensityModel solid_density_model;
 #endif  // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
+    SolidHeatCapacityModel solid_heat_capacity_model;
+    ThermalConductivityModel<DisplacementDim> thermal_conductivity_model;
 };
 }  // namespace ConstitutiveRelations
 }  // namespace ProcessLib::TH2M

ProcessLib/TH2M/ConstitutiveRelations/DarcyVelocity.h

@@ -0,0 +1,36 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2024, 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 "Base.h"
+#include "ProcessLib/Reflection/ReflectionData.h"
+
+namespace ProcessLib::TH2M
+{
+namespace ConstitutiveRelations
+{
+template <int DisplacementDim>
+struct DarcyVelocityData
+{
+    GlobalDimVector<DisplacementDim> w_GS;
+    GlobalDimVector<DisplacementDim> w_LS;
+
+    static auto reflect()
+    {
+        using Self = DarcyVelocityData<DisplacementDim>;
+        namespace R = ProcessLib::Reflection;
+
+        return std::tuple{
+            R::makeReflectionData("velocity_gas", &Self::w_GS),
+            R::makeReflectionData("velocity_liquid", &Self::w_LS)};
+    }
+};
+}  // namespace ConstitutiveRelations
+}  // namespace ProcessLib::TH2M

ProcessLib/TH2M/ConstitutiveRelations/DiffusionVelocity.h

@@ -0,0 +1,42 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2024, 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 "Base.h"
+#include "ProcessLib/Reflection/ReflectionData.h"
+
+namespace ProcessLib::TH2M
+{
+namespace ConstitutiveRelations
+{
+template <int DisplacementDim>
+struct DiffusionVelocityData
+{
+    GlobalDimVector<DisplacementDim> d_CG;
+    GlobalDimVector<DisplacementDim> d_WG;
+    GlobalDimVector<DisplacementDim> d_CL;
+    GlobalDimVector<DisplacementDim> d_WL;
+
+    static auto reflect()
+    {
+        using Self = DiffusionVelocityData<DisplacementDim>;
+        namespace R = ProcessLib::Reflection;
+
+        return std::tuple{
+            R::makeReflectionData("diffusion_velocity_gas_gas", &Self::d_CG),
+            R::makeReflectionData("diffusion_velocity_vapour_gas", &Self::d_WG),
+            R::makeReflectionData("diffusion_velocity_solute_liquid",
+                                  &Self::d_CL),
+            R::makeReflectionData("diffusion_velocity_liquid_liquid",
+                                  &Self::d_WL)};
+    }
+};
+}  // namespace ConstitutiveRelations
+}  // namespace ProcessLib::TH2M

ProcessLib/TH2M/ConstitutiveRelations/Enthalpy.h

@@ -20,16 +20,30 @@ struct EnthalpyData
 {
     double h_G = nan;
     double h_L = nan;
+    double h_S = nan;
 
     static auto reflect()
     {
         using Self = EnthalpyData;
         namespace R = ProcessLib::Reflection;
 
         return std::tuple{R::makeReflectionData("enthalpy_gas", &Self::h_G),
-                          R::makeReflectionData("enthalpy_liquid", &Self::h_L)};
+                          R::makeReflectionData("enthalpy_liquid", &Self::h_L),
+                          R::makeReflectionData("enthalpy_solid", &Self::h_S)};
     }
 };
 
+struct EffectiveVolumetricEnthalpy
+{
+    double rho_h_eff = nan;
+};
+
+struct EffectiveVolumetricEnthalpyDerivatives
+{
+    double drho_h_eff_dT = nan;
+    double drho_h_eff_dp_GR = nan;
+    double drho_h_eff_dp_cap = nan;
+};
+
 }  // namespace ConstitutiveRelations
 }  // namespace ProcessLib::TH2M

ProcessLib/TH2M/ConstitutiveRelations/FluidDensity.h

@@ -19,10 +19,10 @@ namespace ConstitutiveRelations
 struct FluidDensityData
 {
     // gas phase density
-    double rho_GR = 0.;
+    double rho_GR = nan;
 
     // liquid phase density
-    double rho_LR = 0.;
+    double rho_LR = nan;
 
     static auto reflect()
     {

ProcessLib/TH2M/ConstitutiveRelations/InternalEnergy.h

@@ -0,0 +1,29 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2024, 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 "Base.h"
+#include "BaseLib/StrongType.h"
+
+namespace ProcessLib::TH2M
+{
+namespace ConstitutiveRelations
+{
+using InternalEnergyData =
+    BaseLib::StrongType<double, struct InternalEnergyTag>;
+
+struct EffectiveVolumetricInternalEnergyDerivatives
+{
+    double drho_u_eff_dT = nan;
+    double drho_u_eff_dp_GR = nan;
+    double drho_u_eff_dp_cap = nan;
+};
+}  // namespace ConstitutiveRelations
+}  // namespace ProcessLib::TH2M

ProcessLib/TH2M/ConstitutiveRelations/MassMoleFractions.h

@@ -18,10 +18,10 @@ namespace ConstitutiveRelations
 {
 struct MassMoleFractionsData
 {
-    double xnCG = 0.;
-    double xmCG = 0.;
-    double xnWL = 0.;
-    double xmWL = 0.;
+    double xnCG = nan;
+    double xmCG = nan;
+    double xnWL = nan;
+    double xmWL = nan;
 
     static auto reflect()
     {