[PL/HT] Use porosity model from MaterialLib.

ProcessLib/HT/CreateHTProcess.cpp

@@ -17,6 +17,7 @@
 
 #include "MaterialLib/Fluid/Density/createFluidDensityModel.h"
 #include "MaterialLib/Fluid/Viscosity/createViscosityModel.h"
+#include "MaterialLib/PorousMedium/Porosity/createPorosityModel.h"
 
 namespace ProcessLib
 {
@@ -48,12 +49,11 @@ std::unique_ptr<Process> createHTProcess(
         //! \ogs_file_param_special{prj__processes__process__HT__process_variables__pressure}
         "pressure"});
 
-    // Porosity parameter.
-    auto& porosity = findParameter<double>(
-        config,
-        //! \ogs_file_param_special{prj__processes__process__HT__porosity}
-        "porosity", parameters, 1);
-    DBUG("Use \'%s\' as porosity parameter.", porosity.name.c_str());
+    auto const& porous_medium_config = config.getConfigSubtree("porous_medium");
+    auto const& porosity_conf =
+        porous_medium_config.getConfigSubtree("porosity");
+    auto porosity_model =
+        MaterialLib::PorousMedium::createPorosityModel(porosity_conf);
 
     // Parameter for the intrinsic permeability (only one scalar per element,
     // i.e., the isotropic case is handled at the moment)
@@ -157,7 +157,7 @@ std::unique_ptr<Process> createHTProcess(
         std::copy_n(b.data(), b.size(), specific_body_force.data());
 
     HTProcessData process_data{
-        porosity,
+        std::move(porosity_model),
         intrinsic_permeability,
         specific_storage,
         std::move(viscosity_model),

ProcessLib/HT/HTFEM.h

@@ -127,31 +127,36 @@ public:
             auto const intrinsic_permeability =
                 _process_data.intrinsic_permeability(t, pos)[0];
 
-            auto const porosity = _process_data.porosity(t, pos)[0];
-
-            auto const specific_heat_capacity_solid =
-                _process_data.specific_heat_capacity_solid(t, pos)[0];
-            auto const specific_heat_capacity_fluid =
-                _process_data.specific_heat_capacity_fluid(t, pos)[0];
-            double const heat_capacity =
-                density_solid * specific_heat_capacity_solid * (1 - porosity) +
-                fluid_reference_density * specific_heat_capacity_fluid * porosity;
-
             auto const thermal_conductivity_solid =
                 _process_data.thermal_conductivity_solid(t, pos)[0];
             auto const thermal_conductivity_fluid =
                 _process_data.thermal_conductivity_fluid(t, pos)[0];
 
+            auto const& sm = _shape_matrices[ip];
+            double T_int_pt = 0.0;
+            double p_int_pt = 0.0;
+            // Order matters: First T, then P!
+            NumLib::shapeFunctionInterpolate(local_x, sm.N, T_int_pt, p_int_pt);
+
+            // TODO the first argument has to be changed for non constant
+            // porosity model
+            auto const porosity =
+                _process_data.porosity_model->getValue(0.0, T_int_pt);
+
             Eigen::MatrixXd const thermal_conductivity =
                 (thermal_conductivity_solid * (1 - porosity) +
                 thermal_conductivity_fluid * porosity) * unit_mat;
 
+            auto const specific_heat_capacity_solid =
+                _process_data.specific_heat_capacity_solid(t, pos)[0];
+            auto const specific_heat_capacity_fluid =
+                _process_data.specific_heat_capacity_fluid(t, pos)[0];
+
             auto const thermal_dispersivity_longitudinal =
                 _process_data.thermal_dispersivity_longitudinal(t, pos)[0];
             auto const thermal_dispersivity_transversal =
                 _process_data.thermal_dispersivity_transversal(t, pos)[0];
 
-            auto const& sm = _shape_matrices[ip];
             auto const& wp = _integration_method.getWeightedPoint(ip);
             auto Ktt = local_K.template block<num_nodes, num_nodes>(0, 0);
             auto Mtt = local_M.template block<num_nodes, num_nodes>(0, 0);
@@ -161,11 +166,6 @@ public:
                                                                     num_nodes);
             auto Bp = local_b.template block<num_nodes, 1>(num_nodes, 0);
 
-            double T_int_pt = 0.0;
-            double p_int_pt = 0.0;
-            // Order matters: First T, then P!
-            NumLib::shapeFunctionInterpolate(local_x, sm.N, T_int_pt, p_int_pt);
-
             // Use the fluid density model to compute the density
             vars[static_cast<int>(
                 MaterialLib::Fluid::PropertyVariableType::T)] = T_int_pt;
@@ -194,6 +194,10 @@ public:
             auto const hydrodynamic_thermodispersion =
                 thermal_conductivity + thermal_dispersivity;
 
+            double const heat_capacity =
+                density_solid * specific_heat_capacity_solid * (1 - porosity) +
+                fluid_reference_density * specific_heat_capacity_fluid * porosity;
+
             auto const integral_term =
                 sm.integralMeasure * sm.detJ * wp.getWeight();
             // matrix assembly

ProcessLib/HT/HTProcessData.h

@@ -12,6 +12,7 @@
 #include <memory>
 
 #include "MaterialLib/Fluid/FluidProperty.h"
+#include "MaterialLib/PorousMedium/Porosity/Porosity.h"
 
 namespace MeshLib
 {
@@ -28,7 +29,7 @@ namespace HT
 struct HTProcessData
 {
     HTProcessData(
-        ProcessLib::Parameter<double> const& porosity_,
+        std::unique_ptr<MaterialLib::PorousMedium::Porosity>&& porosity_model_,
         ProcessLib::Parameter<double> const& intrinsic_permeability_,
         ProcessLib::Parameter<double> const& specific_storage_,
         std::unique_ptr<MaterialLib::Fluid::FluidProperty>&& viscosity_model_,
@@ -43,7 +44,7 @@ struct HTProcessData
         ProcessLib::Parameter<double> const& thermal_conductivity_fluid_,
         Eigen::Vector3d const& specific_body_force_,
         bool const has_gravity_)
-        : porosity(porosity_),
+        : porosity_model(std::move(porosity_model_)),
           intrinsic_permeability(intrinsic_permeability_),
           specific_storage(specific_storage_),
           viscosity_model(std::move(viscosity_model_)),
@@ -62,7 +63,7 @@ struct HTProcessData
     }
 
     HTProcessData(HTProcessData&& other)
-        : porosity(other.porosity),
+        : porosity_model(other.porosity_model.release()),
           intrinsic_permeability(other.intrinsic_permeability),
           specific_storage(other.specific_storage),
           viscosity_model(other.viscosity_model.release()),
@@ -91,7 +92,7 @@ struct HTProcessData
     //! Assignments are not needed.
     void operator=(HTProcessData&&) = delete;
 
-    Parameter<double> const& porosity;
+    std::unique_ptr<MaterialLib::PorousMedium::Porosity> porosity_model;
     Parameter<double> const& intrinsic_permeability;
     Parameter<double> const& specific_storage;
     std::unique_ptr<MaterialLib::Fluid::FluidProperty> viscosity_model;