Merge branch 'H_body_force' into 'master'

[LiquidFlow]  Use body force vector in the local assember

See merge request ogs/ogs!3623

ProcessLib/LiquidFlow/CreateLiquidFlowProcess.cpp

@@ -90,15 +90,12 @@ std::unique_ptr<Process> createLiquidFlowProcess(
             "dimension is {:d}",
             b.size(), mesh.getDimension());
     }
-    int gravity_axis_id = -1; // default: no gravity
-    double const g = -b[mesh.getDimension()-1];
+
+    Eigen::VectorXd specific_body_force(b.size());
     bool const has_gravity = MathLib::toVector(b).norm() > 0;
     if (has_gravity)
     {
-        if (g != 0.0)
-        {
-            gravity_axis_id = mesh.getDimension()-1;
-        }
+        std::copy_n(b.data(), b.size(), specific_body_force.data());
     }
 
     std::unique_ptr<ProcessLib::SurfaceFluxData> surfaceflux;
@@ -118,7 +115,8 @@ std::unique_ptr<Process> createLiquidFlowProcess(
     checkMPLProperties(mesh, *media_map);
     DBUG("Media properties verified.");
 
-    LiquidFlowData process_data{std::move(media_map), gravity_axis_id, g};
+    LiquidFlowData process_data{std::move(media_map),
+                                std::move(specific_body_force), has_gravity};
 
     return std::make_unique<LiquidFlowProcess>(
         std::move(name), mesh, std::move(jacobian_assembler), parameters,

ProcessLib/LiquidFlow/LiquidFlowData.h

@@ -10,6 +10,7 @@
 
 #pragma once
 
+#include <Eigen/Dense>
 #include <memory>
 
 namespace MaterialPropertyLib
@@ -25,8 +26,8 @@ struct LiquidFlowData final
 {
     std::unique_ptr<MaterialPropertyLib::MaterialSpatialDistributionMap>
         media_map;
-    const int gravitational_axis_id;
-    const double gravitational_acceleration;
+    Eigen::VectorXd const specific_body_force;
+    bool const has_gravity;
 };
 
 }  // namespace LiquidFlow

ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h

@@ -190,9 +190,8 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
 
         // Assemble Laplacian, K, and RHS by the gravitational term
         LaplacianGravityVelocityCalculator::calculateLaplacianAndGravityTerm(
-            local_K, local_b, ip_data, permeability, viscosity,
-            fluid_density * _process_data.gravitational_acceleration,
-            _process_data.gravitational_axis_id);
+            local_K, local_b, ip_data, permeability, viscosity, fluid_density,
+            _process_data);
     }
 }
 
@@ -299,9 +298,8 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
                 medium[MaterialPropertyLib::PropertyType::permeability].value(
                     vars, pos, t, dt));
         LaplacianGravityVelocityCalculator::calculateVelocity(
-            ip, local_p_vec, ip_data, permeability, viscosity,
-            fluid_density * _process_data.gravitational_acceleration,
-            _process_data.gravitational_axis_id, darcy_velocity_at_ips);
+            ip, local_p_vec, ip_data, permeability, viscosity, fluid_density,
+            darcy_velocity_at_ips, _process_data);
     }
 }
 
@@ -314,20 +312,16 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
         IntegrationPointData<NodalRowVectorType,
                              GlobalDimNodalMatrixType> const& ip_data,
         Eigen::MatrixXd const& permeability, double const mu,
-        double const rho_g, int const gravitational_axis_id)
+        double const rho_L, const LiquidFlowData& process_data)
 {
     const double K = permeability(0, 0) / mu;
     const double fac = K * ip_data.integration_weight;
     local_K.noalias() += fac * ip_data.dNdx.transpose() * ip_data.dNdx;
 
-    if (gravitational_axis_id >= 0)
+    if (process_data.has_gravity)
     {
-        // Note: Since permeability, K, is a scalar number in this function,
-        // (gradN)*K*V becomes K*(grad N)*V. With the gravity vector of V,
-        // the simplification of (grad N)*V is the gravitational_axis_id th
-        // column of the transpose of (grad N) multiplied with rho_g.
-        local_b.noalias() -=
-            fac * ip_data.dNdx.transpose().col(gravitational_axis_id) * rho_g;
+        local_b.noalias() += (fac * rho_L) * ip_data.dNdx.transpose() *
+                             process_data.specific_body_force;
     }
 }
 
@@ -339,16 +333,18 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
         IntegrationPointData<NodalRowVectorType,
                              GlobalDimNodalMatrixType> const& ip_data,
         Eigen::MatrixXd const& permeability, double const mu,
-        double const rho_g, int const gravitational_axis_id,
-        MatrixOfVelocityAtIntegrationPoints& darcy_velocity_at_ips)
+        double const rho_L,
+        MatrixOfVelocityAtIntegrationPoints& darcy_velocity_at_ips,
+        const LiquidFlowData& process_data)
 {
     const double K = permeability(0, 0) / mu;
     // Compute the velocity
     darcy_velocity_at_ips.col(ip).noalias() = -K * ip_data.dNdx * local_p;
     // gravity term
-    if (gravitational_axis_id >= 0)
+    if (process_data.has_gravity)
     {
-        darcy_velocity_at_ips.col(ip)[gravitational_axis_id] -= K * rho_g;
+        darcy_velocity_at_ips.col(ip).noalias() +=
+            (K * rho_L) * process_data.specific_body_force;
     }
 }
 
@@ -361,19 +357,16 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
         IntegrationPointData<NodalRowVectorType,
                              GlobalDimNodalMatrixType> const& ip_data,
         Eigen::MatrixXd const& permeability, double const mu,
-        double const rho_g, int const gravitational_axis_id)
+        double const rho_L, const LiquidFlowData& process_data)
 {
     const double fac = ip_data.integration_weight / mu;
     local_K.noalias() +=
         fac * ip_data.dNdx.transpose() * permeability * ip_data.dNdx;
-    if (gravitational_axis_id >= 0)
+
+    if (process_data.has_gravity)
     {
-        // Note: permeability * gravity_vector = rho_g *
-        // permeability.col(gravitational_axis_id)
-        //       i.e the result is the gravitational_axis_id th column of
-        //       permeability multiplied with rho_g.
-        local_b.noalias() -= fac * rho_g * ip_data.dNdx.transpose() *
-                             permeability.col(gravitational_axis_id);
+        local_b.noalias() += (fac * rho_L) * ip_data.dNdx.transpose() *
+                             permeability * process_data.specific_body_force;
     }
 }
 
@@ -385,16 +378,18 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
         IntegrationPointData<NodalRowVectorType,
                              GlobalDimNodalMatrixType> const& ip_data,
         Eigen::MatrixXd const& permeability, double const mu,
-        double const rho_g, int const gravitational_axis_id,
-        MatrixOfVelocityAtIntegrationPoints& darcy_velocity_at_ips)
+        double const rho_L,
+        MatrixOfVelocityAtIntegrationPoints& darcy_velocity_at_ips,
+        const LiquidFlowData& process_data)
 {
     // Compute the velocity
     darcy_velocity_at_ips.col(ip).noalias() =
         -permeability * ip_data.dNdx * local_p / mu;
-    if (gravitational_axis_id >= 0)
+    // gravity term
+    if (process_data.has_gravity)
     {
-        darcy_velocity_at_ips.col(ip).noalias() -=
-            rho_g * permeability.col(gravitational_axis_id) / mu;
+        darcy_velocity_at_ips.col(ip).noalias() +=
+            (rho_L / mu) * permeability * process_data.specific_body_force;
     }
 }
 

ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h

@@ -37,10 +37,9 @@ struct IntegrationPointData final
     explicit IntegrationPointData(NodalRowVectorType const& N_,
                                   GlobalDimNodalMatrixType const& dNdx_,
                                   double const& integration_weight_)
-        : N(N_),
-          dNdx(dNdx_),
-          integration_weight(integration_weight_)
-    {}
+        : N(N_), dNdx(dNdx_), integration_weight(integration_weight_)
+    {
+    }
 
     NodalRowVectorType const N;
     GlobalDimNodalMatrixType const dNdx;
@@ -83,12 +82,11 @@ class LiquidFlowLocalAssembler : public LiquidFlowLocalAssemblerInterface
         Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>;
 
 public:
-    LiquidFlowLocalAssembler(
-        MeshLib::Element const& element,
-        std::size_t const /*local_matrix_size*/,
-        bool const is_axially_symmetric,
-        unsigned const integration_order,
-        LiquidFlowData const& process_data)
+    LiquidFlowLocalAssembler(MeshLib::Element const& element,
+                             std::size_t const /*local_matrix_size*/,
+                             bool const is_axially_symmetric,
+                             unsigned const integration_order,
+                             LiquidFlowData const& process_data)
         : _element(element),
           _integration_method(integration_order),
           _process_data(process_data)
@@ -162,15 +160,16 @@ private:
             IntegrationPointData<NodalRowVectorType,
                                  GlobalDimNodalMatrixType> const& ip_data,
             Eigen::MatrixXd const& permeability, double const mu,
-            double const rho_g, int const gravitational_axis_id);
+            double const rho_L, const LiquidFlowData& process_data);
 
         static void calculateVelocity(
             unsigned const ip, Eigen::Map<const NodalVectorType> const& local_p,
             IntegrationPointData<NodalRowVectorType,
                                  GlobalDimNodalMatrixType> const& ip_data,
             Eigen::MatrixXd const& permeability, double const mu,
-            double const rho_g, int const gravitational_axis_id,
-            MatrixOfVelocityAtIntegrationPoints& darcy_velocity_at_ips);
+            double const rho_L,
+            MatrixOfVelocityAtIntegrationPoints& darcy_velocity_at_ips,
+            const LiquidFlowData& process_data);
     };
 
     /**
@@ -185,15 +184,16 @@ private:
             IntegrationPointData<NodalRowVectorType,
                                  GlobalDimNodalMatrixType> const& ip_data,
             Eigen::MatrixXd const& permeability, double const mu,
-            double const rho_g, int const gravitational_axis_id);
+            double const rho_L, const LiquidFlowData& process_data);
 
         static void calculateVelocity(
             unsigned const ip, Eigen::Map<const NodalVectorType> const& local_p,
             IntegrationPointData<NodalRowVectorType,
                                  GlobalDimNodalMatrixType> const& ip_data,
             Eigen::MatrixXd const& permeability, double const mu,
-            double const rho_g, int const gravitational_axis_id,
-            MatrixOfVelocityAtIntegrationPoints& darcy_velocity_at_ips);
+            double const rho_L,
+            MatrixOfVelocityAtIntegrationPoints& darcy_velocity_at_ips,
+            const LiquidFlowData& process_data);
     };
 
     template <typename LaplacianGravityVelocityCalculator>