Merge branch 'gas_HM' into 'master'

Reduce some redundancy in LiquidFlow and Hydromechanics See merge request ogs/ogs!5132
ufz · Oct 15, 2024 · 17d0669 · 17d0669
2 parents afdb5bd + 9f72c2e
commit 17d0669
Show file tree

Hide file tree

Showing 8 changed files with 143 additions and 97 deletions.
diff --git a/MaterialLib/MPL/Utils/CheckMPLPhasesForSinglePhaseFlow.cpp b/MaterialLib/MPL/Utils/CheckMPLPhasesForSinglePhaseFlow.cpp
@@ -0,0 +1,53 @@
+/**
+ * \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
+ *
+ * Created on October 9, 2024, 6:03 PM
+ */
+
+#include "CheckMPLPhasesForSinglePhaseFlow.h"
+
+#include <range/v3/algorithm/any_of.hpp>
+#include <range/v3/algorithm/for_each.hpp>
+#include <range/v3/view/transform.hpp>
+
+#include "MaterialLib/MPL/MaterialSpatialDistributionMap.h"
+#include "MaterialLib/MPL/Medium.h"
+#include "MeshLib/Elements/Element.h"
+#include "MeshLib/Mesh.h"
+
+namespace MaterialPropertyLib
+{
+void checkMPLPhasesForSinglePhaseFlow(
+    MeshLib::Mesh const& mesh,
+    MaterialPropertyLib::MaterialSpatialDistributionMap const& media_map)
+{
+    // Check all of the elements have a medium defined.
+    ranges::for_each(mesh.getElements() | MeshLib::views::ids,
+                     [&](auto const& element_id)
+                     { media_map.checkElementHasMedium(element_id); });
+
+    // Collect phases of all elements...
+    auto all_phases =
+        media_map.media() |
+        ranges::views::transform([&](auto const& medium)
+                                 { return &fluidPhase(*medium); }) |
+        ranges::to_vector;
+
+    assert(!all_phases.empty());
+
+    // ... and check if any of the phases are different by name.
+    if (ranges::any_of(all_phases,
+                       [p0 = all_phases.front()](auto const* const p)
+                       { return p->name != p0->name; }))
+    {
+        OGS_FATAL(
+            "You are mixing liquid and gas phases in your model domain. OGS "
+            "does not yet know how to handle this.");
+    }
+}
+}  // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Utils/CheckMPLPhasesForSinglePhaseFlow.h b/MaterialLib/MPL/Utils/CheckMPLPhasesForSinglePhaseFlow.h
@@ -0,0 +1,26 @@
+/**
+ * \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
+ *
+ * Created on October 9, 2024, 6:03 PM
+ */
+
+#pragma once
+
+namespace MeshLib
+{
+class Mesh;
+}
+
+namespace MaterialPropertyLib
+{
+class MaterialSpatialDistributionMap;
+
+void checkMPLPhasesForSinglePhaseFlow(
+    MeshLib::Mesh const& mesh,
+    MaterialPropertyLib::MaterialSpatialDistributionMap const& media_map);
+}  // namespace MaterialPropertyLib
diff --git a/ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.cpp b/ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.cpp
@@ -11,13 +11,15 @@
 #include "CreateHydroMechanicsProcess.h"
 
 #include <cassert>
+#include <range/v3/range.hpp>
 #include <string>
 
 #include "HydroMechanicsProcess.h"
 #include "HydroMechanicsProcessData.h"
 #include "MaterialLib/MPL/CreateMaterialSpatialDistributionMap.h"
 #include "MaterialLib/MPL/MaterialSpatialDistributionMap.h"
 #include "MaterialLib/MPL/Medium.h"
+#include "MaterialLib/MPL/Utils/CheckMPLPhasesForSinglePhaseFlow.h"
 #include "MaterialLib/SolidModels/CreateConstitutiveRelation.h"
 #include "MaterialLib/SolidModels/MechanicsBase.h"
 #include "MeshLib/Utils/Is2DMeshOnRotatedVerticalPlane.h"
@@ -215,36 +217,20 @@ std::unique_ptr<Process> createHydroMechanicsProcess(
                                                 MaterialPropertyLib::density};
     std::array const requiredSolidProperties = {MaterialPropertyLib::density};
 
-    // setting default to suppress -Wmaybe-uninitialized warning
-    MaterialPropertyLib::Variable phase_pressure =
-        MaterialPropertyLib::Variable::liquid_phase_pressure;
-    for (auto const& element_id : mesh.getElements() | MeshLib::views::ids)
+    MaterialPropertyLib::checkMPLPhasesForSinglePhaseFlow(mesh, media_map);
+
+    // The uniqueness of phase has already been checked in
+    // `checkMPLPhasesForSinglePhaseFlow`.
+    MaterialPropertyLib::Variable const phase_variable =
+        (*ranges::begin(media_map.media()))->hasPhase("Gas")
+            ? MaterialPropertyLib::Variable::gas_phase_pressure
+            : MaterialPropertyLib::Variable::liquid_phase_pressure;
+    for (auto const& medium : media_map.media())
     {
-        media_map.checkElementHasMedium(element_id);
-        auto const& medium = *media_map.getMedium(element_id);
-        if (element_id == 0)
-        {
-            phase_pressure =
-                medium.hasPhase("Gas")
-                    ? MaterialPropertyLib::Variable::gas_phase_pressure
-                    : MaterialPropertyLib::Variable::liquid_phase_pressure;
-        }
-        else
-        {
-            auto const phase_pressure_tmp =
-                medium.hasPhase("Gas")
-                    ? MaterialPropertyLib::Variable::gas_phase_pressure
-                    : MaterialPropertyLib::Variable::liquid_phase_pressure;
-            if (phase_pressure != phase_pressure_tmp)
-            {
-                OGS_FATAL(
-                    "You are mixing liquid and gas phases in your model domain."
-                    "OGS does not yet know how to handle this.");
-            }
-        }
-        checkRequiredProperties(medium, requiredMediumProperties);
-        checkRequiredProperties(fluidPhase(medium), requiredFluidProperties);
-        checkRequiredProperties(medium.phase("Solid"), requiredSolidProperties);
+        checkRequiredProperties(*medium, requiredMediumProperties);
+        checkRequiredProperties(fluidPhase(*medium), requiredFluidProperties);
+        checkRequiredProperties(medium->phase("Solid"),
+                                requiredSolidProperties);
     }
     DBUG("Media properties verified.");
 
@@ -266,7 +252,7 @@ std::unique_ptr<Process> createHydroMechanicsProcess(
         hydraulic_process_id,
         mechanics_related_process_id,
         use_taylor_hood_elements,
-        phase_pressure};
+        phase_variable};
 
     SecondaryVariableCollection secondary_variables;
 

diff --git a/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h b/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
@@ -182,8 +182,9 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
     auto const& medium = _process_data.media_map.getMedium(_element.getID());
     auto const& solid = medium->phase("Solid");
     auto const& fluid = fluidPhase(*medium);
-    auto const& phase_pressure = _process_data.phase_pressure;
     MPL::VariableArray vars;
+    auto& phase_pressure = medium->hasPhase("Gas") ? vars.gas_phase_pressure
+                                                   : vars.liquid_phase_pressure;
 
     auto const T_ref =
         medium->property(MPL::PropertyType::reference_temperature)
@@ -219,9 +220,7 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
 
         double const p_int_pt = N_p.dot(p);
 
-        // setting both vars equal to enable MPL access for gas and liquid
-        // properties
-        vars.liquid_phase_pressure = vars.gas_phase_pressure = p_int_pt;
+        phase_pressure = p_int_pt;
 
         auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
             t, x_position, dt, T_ref);
@@ -256,10 +255,11 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
         auto const mu = fluid.property(MPL::PropertyType::viscosity)
                             .template value<double>(vars, x_position, t, dt);
 
-        auto const beta_p = fluid.property(MPL::PropertyType::density)
-                                .template dValue<double>(vars, phase_pressure,
-                                                         x_position, t, dt) /
-                            rho_fr;
+        auto const beta_p =
+            fluid.property(MPL::PropertyType::density)
+                .template dValue<double>(vars, _process_data.phase_variable,
+                                         x_position, t, dt) /
+            rho_fr;
 
         // Set mechanical variables for the intrinsic permeability model
         // For stress dependent permeability.
@@ -411,6 +411,8 @@ std::vector<double> const& HydroMechanicsLocalAssembler<
     auto const& medium = _process_data.media_map.getMedium(_element.getID());
     auto const& fluid = fluidPhase(*medium);
     MPL::VariableArray vars;
+    auto& phase_pressure = medium->hasPhase("Gas") ? vars.gas_phase_pressure
+                                                   : vars.liquid_phase_pressure;
 
     // TODO (naumov) Temporary value not used by current material models. Need
     // extension of secondary variables interface.
@@ -427,9 +429,7 @@ std::vector<double> const& HydroMechanicsLocalAssembler<
 
         double const p_int_pt = _ip_data[ip].N_p.dot(p);
 
-        // setting both vars equal to enable MPL access for gas and liquid
-        // properties
-        vars.liquid_phase_pressure = vars.gas_phase_pressure = p_int_pt;
+        phase_pressure = p_int_pt;
 
         auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
                                .template value<double>(vars, x_position, t, dt);
@@ -535,8 +535,9 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
 
     auto const& medium = _process_data.media_map.getMedium(_element.getID());
     auto const& fluid = fluidPhase(*medium);
-    auto const& phase_pressure = _process_data.phase_pressure;
     MPL::VariableArray vars;
+    auto& phase_pressure = medium->hasPhase("Gas") ? vars.gas_phase_pressure
+                                                   : vars.liquid_phase_pressure;
 
     auto const T_ref =
         medium->property(MPL::PropertyType::reference_temperature)
@@ -563,9 +564,7 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
 
         double const p_int_pt = N_p.dot(p);
 
-        // setting both vars equal to enable MPL access for gas and liquid
-        // properties
-        vars.liquid_phase_pressure = vars.gas_phase_pressure = p_int_pt;
+        phase_pressure = p_int_pt;
 
         auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
             t, x_position, dt, T_ref);
@@ -615,10 +614,11 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
 
         auto const mu = fluid.property(MPL::PropertyType::viscosity)
                             .template value<double>(vars, x_position, t, dt);
-        auto const beta_p = fluid.property(MPL::PropertyType::density)
-                                .template dValue<double>(vars, phase_pressure,
-                                                         x_position, t, dt) /
-                            rho_fr;
+        auto const beta_p =
+            fluid.property(MPL::PropertyType::density)
+                .template dValue<double>(vars, _process_data.phase_variable,
+                                         x_position, t, dt) /
+            rho_fr;
 
         auto const K_over_mu = K / mu;
 
@@ -701,6 +701,8 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
     auto const& solid = medium->phase("Solid");
     auto const& fluid = fluidPhase(*medium);
     MPL::VariableArray vars;
+    auto& phase_pressure = medium->hasPhase("Gas") ? vars.gas_phase_pressure
+                                                   : vars.liquid_phase_pressure;
 
     auto const T_ref =
         medium->property(MPL::PropertyType::reference_temperature)
@@ -731,9 +733,7 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
         auto& eps = _ip_data[ip].eps;
         auto const& sigma_eff = _ip_data[ip].sigma_eff;
 
-        // setting both vars equal to enable MPL access for gas and liquid
-        // properties
-        vars.liquid_phase_pressure = vars.gas_phase_pressure = N_p.dot(p);
+        phase_pressure = N_p.dot(p);
 
         auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
                                .template value<double>(vars, x_position, t, dt);
@@ -1159,6 +1159,8 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
 
     auto const& medium = _process_data.media_map.getMedium(elem_id);
     MPL::VariableArray vars;
+    auto& phase_pressure = medium->hasPhase("Gas") ? vars.gas_phase_pressure
+                                                   : vars.liquid_phase_pressure;
 
     SymmetricTensor k_sum = SymmetricTensor::Zero(KelvinVectorSize);
     auto sigma_eff_sum = MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
@@ -1177,9 +1179,7 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
                                .template value<double>(vars, x_position, t, dt);
         double const p_int_pt = _ip_data[ip].N_p.dot(p);
 
-        // setting both vars equal to enable MPL access for gas and liquid
-        // properties
-        vars.liquid_phase_pressure = vars.gas_phase_pressure = p_int_pt;
+        phase_pressure = p_int_pt;
 
         // Set mechanical variables for the intrinsic permeability model
         // For stress dependent permeability.

diff --git a/ProcessLib/HydroMechanics/HydroMechanicsProcessData.h b/ProcessLib/HydroMechanics/HydroMechanicsProcessData.h
@@ -90,7 +90,7 @@ struct HydroMechanicsProcessData
 
     const bool use_taylor_hood_elements;
 
-    MaterialPropertyLib::Variable phase_pressure;
+    MaterialPropertyLib::Variable const phase_variable;
 
     MeshLib::PropertyVector<double>* pressure_interpolated = nullptr;
     std::array<MeshLib::PropertyVector<double>*, 3> principal_stress_vector = {

diff --git a/ProcessLib/LiquidFlow/CreateLiquidFlowProcess.cpp b/ProcessLib/LiquidFlow/CreateLiquidFlowProcess.cpp
@@ -12,17 +12,13 @@
 #include "CreateLiquidFlowProcess.h"
 
 #include <algorithm>
-#include <range/v3/algorithm/any_of.hpp>
-#include <range/v3/algorithm/for_each.hpp>
-#include <range/v3/range/conversion.hpp>
-#include <range/v3/view/transform.hpp>
 #include <typeinfo>
 
 #include "LiquidFlowProcess.h"
 #include "MaterialLib/MPL/CheckMaterialSpatialDistributionMap.h"
 #include "MaterialLib/MPL/CreateMaterialSpatialDistributionMap.h"
 #include "MaterialLib/MPL/Properties/Constant.h"
-#include "MaterialLib/PhysicalConstant.h"
+#include "MaterialLib/MPL/Utils/CheckMPLPhasesForSinglePhaseFlow.h"
 #include "MeshLib/Utils/GetElementRotationMatrices.h"
 #include "MeshLib/Utils/GetSpaceDimension.h"
 #include "ParameterLib/Utils.h"
@@ -38,29 +34,7 @@ void checkMPLProperties(
     MaterialPropertyLib::MaterialSpatialDistributionMap const& media_map,
     bool const is_equation_type_volume)
 {
-    // Check all of the elements have a medium defined.
-    ranges::for_each(mesh.getElements() | MeshLib::views::ids,
-                     [&](auto const& element_id)
-                     { media_map.checkElementHasMedium(element_id); });
-
-    // Collect phases of all elements...
-    auto all_phases =
-        media_map.media() |
-        ranges::views::transform([&](auto const& medium)
-                                 { return &fluidPhase(*medium); }) |
-        ranges::to_vector;
-
-    assert(!all_phases.empty());
-
-    // ... and check if any of the phases are different by name.
-    if (ranges::any_of(all_phases,
-                       [p0 = all_phases.front()](auto const* const p)
-                       { return p->name != p0->name; }))
-    {
-        OGS_FATAL(
-            "You are mixing liquid and gas phases in your model domain. OGS "
-            "does not yet know how to handle this.");
-    }
+    MaterialPropertyLib::checkMPLPhasesForSinglePhaseFlow(mesh, media_map);
 
     std::array const required_medium_properties = {
         MaterialPropertyLib::reference_temperature,
@@ -214,6 +188,13 @@ std::unique_ptr<Process> createLiquidFlowProcess(
             *aperture_config, "parameter", parameters, 1);
     }
 
+    // The uniqueness of phase has already been checked in
+    // `checkMPLProperties`.
+    MaterialPropertyLib::Variable const phase_variable =
+        (*ranges::begin(media_map.media()))->hasPhase("Gas")
+            ? MaterialPropertyLib::Variable::gas_phase_pressure
+            : MaterialPropertyLib::Variable::liquid_phase_pressure;
+
     LiquidFlowData process_data{
         covertEquationBalanceTypeFromString(equation_balance_type_str),
         std::move(media_map),
@@ -223,7 +204,8 @@ std::unique_ptr<Process> createLiquidFlowProcess(
         std::move(specific_body_force),
         has_gravity,
         *aperture_size_parameter,
-        NumLib::ShapeMatrixCache{integration_order}};
+        NumLib::ShapeMatrixCache{integration_order},
+        phase_variable};
 
     return std::make_unique<LiquidFlowProcess>(
         std::move(name), mesh, std::move(jacobian_assembler), parameters,

diff --git a/ProcessLib/LiquidFlow/LiquidFlowData.h b/ProcessLib/LiquidFlow/LiquidFlowData.h
@@ -13,6 +13,7 @@
 #include <Eigen/Core>
 
 #include "MaterialLib/MPL/MaterialSpatialDistributionMap.h"
+#include "MaterialLib/MPL/VariableType.h"
 #include "NumLib/Fem/ShapeMatrixCache.h"
 #include "ParameterLib/Parameter.h"
 
@@ -54,6 +55,8 @@ struct LiquidFlowData final
 
     /// caches for each mesh element type the shape matrix
     NumLib::ShapeMatrixCache shape_matrix_cache;
+
+    MaterialPropertyLib::Variable const phase_variable;
 };
 
 }  // namespace LiquidFlow