[PL] Enable arbitrary perm. laws in LiquidFlow.

ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h

@@ -32,8 +32,10 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
     pos.setElementID(_element.getID());
     const int material_id = _material_properties.getMaterialID(pos);
 
+    double const pressure = std::nan("");
     const Eigen::MatrixXd& permeability = _material_properties.getPermeability(
-        material_id, t, pos, _element.getDimension());
+        material_id, t, pos, _element.getDimension(), pressure,
+        _reference_temperature);
     // Note: For Inclined 1D in 2D/3D or 2D element in 3D, the first item in
     //  the assert must be changed to permeability.rows() ==
     //  _element->getDimension()
@@ -42,14 +44,12 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
     if (permeability.size() == 1)
     {  // isotropic or 1D problem.
         assembleMatrixAndVector<IsotropicCalculator>(
-            material_id, t, local_x, local_M_data, local_K_data, local_b_data,
-            pos, permeability);
+            material_id, t, local_x, local_M_data, local_K_data, local_b_data);
     }
     else
     {
         assembleMatrixAndVector<AnisotropicCalculator>(
-            material_id, t, local_x, local_M_data, local_K_data, local_b_data,
-            pos, permeability);
+            material_id, t, local_x, local_M_data, local_K_data, local_b_data);
     }
 }
 
@@ -61,9 +61,7 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
                             std::vector<double> const& local_x,
                             std::vector<double>& local_M_data,
                             std::vector<double>& local_K_data,
-                            std::vector<double>& local_b_data,
-                            ParameterLib::SpatialPosition const& pos,
-                            Eigen::MatrixXd const& permeability)
+                            std::vector<double>& local_b_data)
 {
     auto const local_matrix_size = local_x.size();
     assert(local_matrix_size == ShapeFunction::NPOINTS);
@@ -78,6 +76,9 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
     unsigned const n_integration_points =
         _integration_method.getNumberOfPoints();
 
+    ParameterLib::SpatialPosition pos;
+    pos.setElementID(_element.getID());
+
     // TODO: The following two variables should be calculated inside the
     //       the integration loop for non-constant porosity and storage models.
     double porosity_variable = 0.;
@@ -104,6 +105,11 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
         const double mu =
             _material_properties.getViscosity(p, _reference_temperature);
 
+        pos.setIntegrationPoint(ip);
+        auto const& permeability = _material_properties.getPermeability(
+            material_id, t, pos, _element.getDimension(), p,
+            _reference_temperature);
+
         // Assemble Laplacian, K, and RHS by the gravitational term
         LaplacianGravityVelocityCalculator::calculateLaplacianAndGravityTerm(
             local_K, local_b, ip_data, permeability, mu, rho_g,
@@ -131,8 +137,12 @@ LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
     ParameterLib::SpatialPosition pos;
     pos.setElementID(_element.getID());
     const int material_id = _material_properties.getMaterialID(pos);
+    // evaluate the permeability to distinguish which computeDarcyVelocity
+    // method should be used
+    double const pressure = std::nan("");
     const Eigen::MatrixXd& permeability = _material_properties.getPermeability(
-        material_id, t, pos, _element.getDimension());
+        material_id, t, pos, _element.getDimension(), pressure,
+        _reference_temperature);
 
     // Note: For Inclined 1D in 2D/3D or 2D element in 3D, the first item in
     //  the assert must be changed to perm.rows() == _element->getDimension()
@@ -140,13 +150,13 @@ LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
 
     if (permeability.size() == 1)
     {  // isotropic or 1D problem.
-        computeDarcyVelocityLocal<IsotropicCalculator>(local_x, permeability,
-                                                       velocity_cache_vectors);
+        computeDarcyVelocityLocal<IsotropicCalculator>(
+            material_id, t, local_x, pos, velocity_cache_vectors);
     }
     else
     {
         computeDarcyVelocityLocal<AnisotropicCalculator>(
-            local_x, permeability, velocity_cache_vectors);
+            material_id, t, local_x, pos, velocity_cache_vectors);
     }
     return velocity_cache;
 }
@@ -156,8 +166,10 @@ template <typename ShapeFunction, typename IntegrationMethod,
 template <typename LaplacianGravityVelocityCalculator>
 void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
     computeDarcyVelocityLocal(
+        const int material_id,
+        const double t,
         std::vector<double> const& local_x,
-        Eigen::MatrixXd const& permeability,
+        ParameterLib::SpatialPosition const& pos,
         MatrixOfVelocityAtIntegrationPoints& darcy_velocity_at_ips) const
 {
     auto const local_matrix_size = local_x.size();
@@ -182,6 +194,9 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
         const double mu =
             _material_properties.getViscosity(p, _reference_temperature);
 
+        auto const& permeability = _material_properties.getPermeability(
+            material_id, t, pos, _element.getDimension(), p,
+            _reference_temperature);
         LaplacianGravityVelocityCalculator::calculateVelocity(
             ip, local_p_vec, ip_data, permeability, mu, rho_g,
             _gravitational_axis_id, darcy_velocity_at_ips);

ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h

@@ -204,14 +204,14 @@ private:
                                  std::vector<double> const& local_x,
                                  std::vector<double>& local_M_data,
                                  std::vector<double>& local_K_data,
-                                 std::vector<double>& local_b_data,
-                                 ParameterLib::SpatialPosition const& pos,
-                                 Eigen::MatrixXd const& permeability);
+                                 std::vector<double>& local_b_data);
 
     template <typename LaplacianGravityVelocityCalculator>
     void computeDarcyVelocityLocal(
+        const int material_id,
+        const double t,
         std::vector<double> const& local_x,
-        Eigen::MatrixXd const& permeability,
+        ParameterLib::SpatialPosition const& pos,
         MatrixOfVelocityAtIntegrationPoints& darcy_velocity_at_ips) const;
 
     const int _gravitational_axis_id;

ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.cpp

@@ -97,10 +97,10 @@ double LiquidFlowMaterialProperties::getMassCoefficient(
 
 Eigen::MatrixXd LiquidFlowMaterialProperties::getPermeability(
     const int material_id, const double t,
-    const ParameterLib::SpatialPosition& pos, const int /*dim*/) const
+    const ParameterLib::SpatialPosition& pos, const int /*dim*/, double const p,
+    double const T) const
 {
-    return _intrinsic_permeability_models[material_id]->getValue(t, pos, 0.0,
-                                                                 0.0);
+    return _intrinsic_permeability_models[material_id]->getValue(t, pos, p, T);
 }
 
 }  // namespace LiquidFlow

ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.h

@@ -91,10 +91,10 @@ public:
                               const double porosity_variable,
                               const double storage_variable) const;
 
-    Eigen::MatrixXd getPermeability(const int material_id,
-                                    const double t,
+    Eigen::MatrixXd getPermeability(const int material_id, const double t,
                                     const ParameterLib::SpatialPosition& pos,
-                                    const int dim) const;
+                                    const int dim, const double p,
+                                    const double T) const;
 
     double getLiquidDensity(const double p, const double T) const;