Merge branch 'liquidflow_gas' into 'master'

Gas flow calculation in LiquidFLow

See merge request ogs/ogs!5112

Documentation/ProjectFile/prj/processes/process/LIQUID_FLOW/t_equation_balance_type.md

@@ -0,0 +1 @@
+\copydoc ProcessLib::LiquidFlow::LiquidFlowData::equation_balance_type

MaterialLib/MPL/MaterialSpatialDistributionMap.h

@@ -13,6 +13,7 @@
 
 #include <map>
 #include <memory>
+#include <range/v3/view.hpp>
 #include <vector>
 
 namespace MeshLib
@@ -35,6 +36,8 @@ class MaterialSpatialDistributionMap
  {
  }
 
+ auto media() const { return media_ | ranges::views::values; }
+
  Medium* getMedium(std::size_t element_id);
  Medium const* getMedium(std::size_t element_id) const;
  void checkElementHasMedium(std::size_t const element_id) const;

ProcessLib/LiquidFlow/CreateLiquidFlowProcess.cpp

@@ -12,10 +12,16 @@
 #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 "MeshLib/Utils/GetElementRotationMatrices.h"
 #include "MeshLib/Utils/GetSpaceDimension.h"
@@ -29,8 +35,33 @@ namespace LiquidFlow
 {
 void checkMPLProperties(
  MeshLib::Mesh const& mesh,
- MaterialPropertyLib::MaterialSpatialDistributionMap const& media_map)
+ 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.");
+ }
+
  std::array const required_medium_properties = {
  MaterialPropertyLib::reference_temperature,
  MaterialPropertyLib::PropertyType::permeability,
@@ -41,17 +72,50 @@ void checkMPLProperties(
  MaterialPropertyLib::PropertyType::viscosity,
  MaterialPropertyLib::PropertyType::density};
 
- std::array<MaterialPropertyLib::PropertyType, 0> const
- required_solid_properties{};
-
- std::array<MaterialPropertyLib::PropertyType, 0> const
- required_gas_properties{};
+ // Check Constant-type density.
+ if (is_equation_type_volume)
+ {
+ for (auto const& medium : media_map.media())
+ {
+ // auto const& medium = *media_map.getMedium(element_id);
+ auto const& fluid_phase_density =
+ fluidPhase(*medium)[MaterialPropertyLib::PropertyType::density];
+ if (typeid(fluid_phase_density) !=
+ typeid(MaterialPropertyLib::Constant))
+ {
+ OGS_FATAL(
+ "Since `equation_balance_type` is set to `volume`,the "
+ "phase density type must be `Constant`. Note: by "
+ "default, `equation_balance_type` is set to `volume`.");
+ }
+ }
+ }
 
- MaterialPropertyLib::checkMaterialSpatialDistributionMap(
- mesh, media_map, required_medium_properties, required_solid_properties,
- required_liquid_properties, required_gas_properties);
+ for (auto const& medium : media_map.media())
+ {
+ checkRequiredProperties(*medium, required_medium_properties);
+ checkRequiredProperties(fluidPhase(*medium),
+ required_liquid_properties);
+ }
+ DBUG("Media properties verified.");
 }
 
+EquationBalanceType covertEquationBalanceTypeFromString(
+ std::string_view const type_in_string)
+{
+ if (type_in_string == "volume")
+ {
+ return EquationBalanceType::volume;
+ }
+ if (type_in_string == "mass")
+ {
+ return EquationBalanceType::mass;
+ }
+ OGS_FATAL(
+ "The value of `equation_balance_type` must be either `volume` or "
+ "`mass`");
+};
+
 std::unique_ptr<Process> createLiquidFlowProcess(
  std::string const& name,
  MeshLib::Mesh& mesh,
@@ -128,8 +192,13 @@ std::unique_ptr<Process> createLiquidFlowProcess(
  INFO("LiquidFlow process is set to be linear.");
  }
 
+ auto const equation_balance_type_str =
+ //! \ogs_file_param{prj__processes__process__LIQUID_FLOW__equation_balance_type}
+ config.getConfigParameter<std::string>("equation_balance_type",
+ "volume");
+
  DBUG("Check the media properties of LiquidFlow process ...");
- checkMPLProperties(mesh, media_map);
+ checkMPLProperties(mesh, media_map, equation_balance_type_str == "volume");
  DBUG("Media properties verified.");
 
  auto const* aperture_size_parameter = &ParameterLib::findParameter<double>(
@@ -146,6 +215,7 @@ std::unique_ptr<Process> createLiquidFlowProcess(
  }
 
  LiquidFlowData process_data{
+ covertEquationBalanceTypeFromString(equation_balance_type_str),
  std::move(media_map),
  MeshLib::getElementRotationMatrices(
  mesh_space_dimension, mesh.getDimension(), mesh.getElements()),

ProcessLib/LiquidFlow/LiquidFlowData.h

@@ -20,8 +20,24 @@ namespace ProcessLib
 {
 namespace LiquidFlow
 {
+
+/// Governing equation balance type
+enum class EquationBalanceType
+{
+ volume,
+ mass
+};
+
 struct LiquidFlowData final
 {
+ /// This indicates whether the governing equation is a volume balance or a
+ /// mass balance. Its value can be `volume` or `mass`. If it is set to
+ /// `volume`, note that the phase density must be constant, and the unit of
+ /// the Neumann boundary condition is m/s. Otherwise, the unit of the
+ /// Neumann boundary condition is kg/m³·m/s = kg/m²/s. By default, it is set
+ /// to `volume`.
+ EquationBalanceType const equation_balance_type;
+
  MaterialPropertyLib::MaterialSpatialDistributionMap media_map;
 
  /// A vector of the rotation matrices for all elements.

ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h

@@ -81,7 +81,7 @@ Eigen::Vector3d LiquidFlowLocalAssembler<ShapeFunction, GlobalDim>::getFlux(
  pos.setElementID(_element.getID());
 
  auto const& medium = *_process_data.media_map.getMedium(_element.getID());
- auto const& liquid_phase = medium.phase("AqueousLiquid");
+ auto const& fluid_phase = fluidPhase(medium);
 
  MaterialPropertyLib::VariableArray vars;
 
@@ -94,7 +94,7 @@ Eigen::Vector3d LiquidFlowLocalAssembler<ShapeFunction, GlobalDim>::getFlux(
  medium[MaterialPropertyLib::PropertyType::permeability].value(
  vars, pos, t, dt));
  auto const viscosity =
- liquid_phase[MaterialPropertyLib::PropertyType::viscosity]
+ fluid_phase[MaterialPropertyLib::PropertyType::viscosity]
  .template value<double>(vars, pos, t, dt);
 
  Eigen::Vector3d flux(0.0, 0.0, 0.0);
@@ -123,6 +123,8 @@ void LiquidFlowLocalAssembler<ShapeFunction, GlobalDim>::
  local_K_data, local_matrix_size, local_matrix_size);
  auto local_b = MathLib::createZeroedVector<NodalVectorType>(
  local_b_data, local_matrix_size);
+ const auto local_p_vec =
+ MathLib::toVector<NodalVectorType>(local_x, local_matrix_size);
 
  unsigned const n_integration_points =
  _integration_method.getNumberOfPoints();
@@ -131,9 +133,16 @@ void LiquidFlowLocalAssembler<ShapeFunction, GlobalDim>::
  pos.setElementID(_element.getID());
 
  auto const& medium = *_process_data.media_map.getMedium(_element.getID());
- auto const& liquid_phase = medium.phase("AqueousLiquid");
+ auto const& fluid_phase = fluidPhase(medium);
+ bool is_gas_phase = medium.hasPhase("Gas");
+ auto const fluid_pressure_variable =
+ is_gas_phase ? MaterialPropertyLib::Variable::gas_phase_pressure
+ : MaterialPropertyLib::Variable::liquid_phase_pressure;
 
  MaterialPropertyLib::VariableArray vars;
+ auto& phase_pressue =
+ is_gas_phase ? vars.gas_phase_pressure : vars.liquid_phase_pressure;
+
  vars.temperature =
  medium[MaterialPropertyLib::PropertyType::reference_temperature]
  .template value<double>(vars, pos, t, dt);
@@ -151,38 +160,36 @@ void LiquidFlowLocalAssembler<ShapeFunction, GlobalDim>::
  auto const& ip_data = _ip_data[ip];
  auto const& N = Ns[ip];
 
- double p = 0.;
- NumLib::shapeFunctionInterpolate(local_x, N, p);
- vars.liquid_phase_pressure = p;
-
- // Compute density:
- auto const fluid_density =
- liquid_phase[MaterialPropertyLib::PropertyType::density]
- .template value<double>(vars, pos, t, dt);
- assert(fluid_density > 0.);
- vars.density = fluid_density;
+ phase_pressue = N.dot(local_p_vec);
 
- auto const ddensity_dpressure =
- liquid_phase[MaterialPropertyLib::PropertyType::density]
- .template dValue<double>(
- vars, MaterialPropertyLib::Variable::liquid_phase_pressure,
- pos, t, dt);
+ auto const [fluid_density, viscosity] =
+ getFluidDensityAndViscosity(t, dt, pos, fluid_phase, vars);
 
  auto const porosity =
  medium[MaterialPropertyLib::PropertyType::porosity]
  .template value<double>(vars, pos, t, dt);
- auto const storage = medium[MaterialPropertyLib::PropertyType::storage]
- .template value<double>(vars, pos, t, dt);
+ auto const specific_storage =
+ medium[MaterialPropertyLib::PropertyType::storage]
+ .template value<double>(vars, pos, t, dt);
 
+ auto const get_drho_dp = [&]()
+ {
+ return fluid_phase[MaterialPropertyLib::PropertyType::density]
+ .template dValue<double>(vars, fluid_pressure_variable, pos, t,
+ dt);
+ };
+ auto const storage =
+ _process_data.equation_balance_type == EquationBalanceType::volume
+ ? specific_storage
+ : specific_storage + porosity * get_drho_dp() / fluid_density;
+
+ double const scaling_factor =
+ _process_data.equation_balance_type == EquationBalanceType::volume
+ ? 1.0
+ : fluid_density;
  // Assemble mass matrix, M
- local_M.noalias() +=
- (porosity * ddensity_dpressure / fluid_density + storage) *
- N.transpose() * N * ip_data.integration_weight;
-
- // Compute viscosity:
- auto const viscosity =
- liquid_phase[MaterialPropertyLib::PropertyType::viscosity]
- .template value<double>(vars, pos, t, dt);
+ local_M.noalias() += scaling_factor * storage * N.transpose() * N *
+ ip_data.integration_weight;
 
  pos.setIntegrationPoint(ip);
  GlobalDimMatrixType const permeability =
@@ -192,8 +199,9 @@ void LiquidFlowLocalAssembler<ShapeFunction, GlobalDim>::
 
  // Assemble Laplacian, K, and RHS by the gravitational term
  LaplacianGravityVelocityCalculator::calculateLaplacianAndGravityTerm(
- local_K, local_b, ip_data, permeability, viscosity, fluid_density,
- projected_body_force_vector, _process_data.has_gravity);
+ local_K, local_b, ip_data, scaling_factor * permeability, viscosity,
+ fluid_density, projected_body_force_vector,
+ _process_data.has_gravity);
  }
 }
 
@@ -284,9 +292,12 @@ void LiquidFlowLocalAssembler<ShapeFunction, GlobalDim>::
  _integration_method.getNumberOfPoints();
 
  auto const& medium = *_process_data.media_map.getMedium(_element.getID());
- auto const& liquid_phase = medium.phase("AqueousLiquid");
+ auto const& fluid_phase = fluidPhase(medium);
 
  MaterialPropertyLib::VariableArray vars;
+ auto& phase_pressue = medium.hasPhase("Gas") ? vars.gas_phase_pressure
+ : vars.liquid_phase_pressure;
+
  vars.temperature =
  medium[MaterialPropertyLib::PropertyType::reference_temperature]
  .template value<double>(vars, pos, t, dt);
@@ -303,20 +314,11 @@ void LiquidFlowLocalAssembler<ShapeFunction, GlobalDim>::
  {
  auto const& ip_data = _ip_data[ip];
  auto const& N = Ns[ip];
- double p = 0.;
- NumLib::shapeFunctionInterpolate(local_x, N, p);
- vars.liquid_phase_pressure = p;
 
- // Compute density:
- auto const fluid_density =
- liquid_phase[MaterialPropertyLib::PropertyType::density]
- .template value<double>(vars, pos, t, dt);
- vars.density = fluid_density;
+ phase_pressue = N.dot(local_p_vec);
 
- // Compute viscosity:
- auto const viscosity =
- liquid_phase[MaterialPropertyLib::PropertyType::viscosity]
- .template value<double>(vars, pos, t, dt);
+ auto const [fluid_density, viscosity] =
+ getFluidDensityAndViscosity(t, dt, pos, fluid_phase, vars);
 
  GlobalDimMatrixType const permeability =
  MaterialPropertyLib::formEigenTensor<GlobalDim>(
@@ -336,11 +338,11 @@ void LiquidFlowLocalAssembler<ShapeFunction, GlobalDim>::IsotropicCalculator::
  Eigen::Map<NodalMatrixType>& local_K,
  Eigen::Map<NodalVectorType>& local_b,
  IntegrationPointData<GlobalDimNodalMatrixType> const& ip_data,
- GlobalDimMatrixType const& permeability, double const mu,
- double const rho_L, GlobalDimVectorType const& specific_body_force,
- bool const has_gravity)
+ GlobalDimMatrixType const& permeability_with_density_factor,
+ double const mu, double const rho_L,
+ GlobalDimVectorType const& specific_body_force, bool const has_gravity)
 {
- const double K = permeability(0, 0) / mu;
+ const double K = permeability_with_density_factor(0, 0) / mu;
  const double fac = K * ip_data.integration_weight;
  local_K.noalias() += fac * ip_data.dNdx.transpose() * ip_data.dNdx;
 
@@ -378,18 +380,19 @@ void LiquidFlowLocalAssembler<ShapeFunction, GlobalDim>::AnisotropicCalculator::
  Eigen::Map<NodalMatrixType>& local_K,
  Eigen::Map<NodalVectorType>& local_b,
  IntegrationPointData<GlobalDimNodalMatrixType> const& ip_data,
- GlobalDimMatrixType const& permeability, double const mu,
- double const rho_L, GlobalDimVectorType const& specific_body_force,
- bool const has_gravity)
+ GlobalDimMatrixType const& permeability_with_density_factor,
+ double const mu, double const rho_L,
+ GlobalDimVectorType const& specific_body_force, bool const has_gravity)
 {
  const double fac = ip_data.integration_weight / mu;
- local_K.noalias() +=
- fac * ip_data.dNdx.transpose() * permeability * ip_data.dNdx;
+ local_K.noalias() += fac * ip_data.dNdx.transpose() *
+  permeability_with_density_factor * ip_data.dNdx;
 
  if (has_gravity)
  {
  local_b.noalias() += (fac * rho_L) * ip_data.dNdx.transpose() *
- permeability * specific_body_force;
+ permeability_with_density_factor *
+ specific_body_force;
  }
 }