From 399eacc3ee7cad75b6e390414cdba4906e98d1e7 Mon Sep 17 00:00:00 2001
From: Florian Zill <florian.zill@ufz.de>
Date: Thu, 23 Jan 2020 17:23:21 +0100
Subject: [PATCH] [HM] Linkage to MPL
for now porosity and viscosity only,
the rest follows later
---
Applications/ApplicationsLib/ProjectData.cpp | 4 +-
.../CreateHydroMechanicsProcess.cpp | 40 ++++++-------
.../CreateHydroMechanicsProcess.h | 14 ++++-
.../HydroMechanics/HydroMechanicsFEM-impl.h | 56 ++++++++++++++++---
.../HydroMechanicsProcessData.h | 10 ++--
5 files changed, 81 insertions(+), 43 deletions(-)
diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp
index e9d3bc5bc60..1f52451511b 100644
--- a/Applications/ApplicationsLib/ProjectData.cpp
+++ b/Applications/ApplicationsLib/ProjectData.cpp
@@ -574,7 +574,7 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config,
std::move(jacobian_assembler),
_process_variables, _parameters,
_local_coordinate_system, integration_order,
- process_config);
+ process_config, _media);
break;
case 3:
process =
@@ -583,7 +583,7 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config,
std::move(jacobian_assembler),
_process_variables, _parameters,
_local_coordinate_system, integration_order,
- process_config);
+ process_config, _media);
break;
default:
OGS_FATAL(
diff --git a/ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.cpp b/ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.cpp
index 563c2393112..ea1b29bd97d 100644
--- a/ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.cpp
+++ b/ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.cpp
@@ -12,6 +12,9 @@
#include <cassert>
+#include "MaterialLib/MPL/CreateMaterialSpatialDistributionMap.h"
+#include "MaterialLib/MPL/MaterialSpatialDistributionMap.h"
+
#include "MaterialLib/SolidModels/CreateConstitutiveRelation.h"
#include "MaterialLib/SolidModels/MechanicsBase.h"
#include "ParameterLib/Utils.h"
@@ -27,15 +30,14 @@ namespace HydroMechanics
{
template <int DisplacementDim>
std::unique_ptr<Process> createHydroMechanicsProcess(
- std::string name,
- MeshLib::Mesh& mesh,
+ std::string name, MeshLib::Mesh& mesh,
std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
std::vector<ProcessVariable> const& variables,
std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
boost::optional<ParameterLib::CoordinateSystem> const&
local_coordinate_system,
- unsigned const integration_order,
- BaseLib::ConfigTree const& config)
+ unsigned const integration_order, BaseLib::ConfigTree const& config,
+ std::map<int, std::unique_ptr<MaterialPropertyLib::Medium>> const& media)
{
//! \ogs_file_param{prj__processes__process__type}
config.checkConfigParameter("type", "HYDRO_MECHANICS");
@@ -118,14 +120,6 @@ std::unique_ptr<Process> createHydroMechanicsProcess(
DBUG("Use '%s' as intrinsic conductivity parameter.",
intrinsic_permeability.name.c_str());
- // Fluid viscosity
- auto& fluid_viscosity = ParameterLib::findParameter<double>(
- config,
- //! \ogs_file_param_special{prj__processes__process__HYDRO_MECHANICS__fluid_viscosity}
- "fluid_viscosity", parameters, 1, &mesh);
- DBUG("Use '%s' as fluid viscosity parameter.",
- fluid_viscosity.name.c_str());
-
// Fluid density
auto& fluid_density = ParameterLib::findParameter<double>(
config,
@@ -141,13 +135,6 @@ std::unique_ptr<Process> createHydroMechanicsProcess(
DBUG("Use '%s' as Biot coefficient parameter.",
biot_coefficient.name.c_str());
- // Porosity
- auto& porosity = ParameterLib::findParameter<double>(
- config,
- //! \ogs_file_param_special{prj__processes__process__HYDRO_MECHANICS__porosity}
- "porosity", parameters, 1, &mesh);
- DBUG("Use '%s' as porosity parameter.", porosity.name.c_str());
-
// Solid density
auto& solid_density = ParameterLib::findParameter<double>(
config,
@@ -174,6 +161,9 @@ std::unique_ptr<Process> createHydroMechanicsProcess(
std::copy_n(b.data(), b.size(), specific_body_force.data());
}
+ auto media_map =
+ MaterialPropertyLib::createMaterialSpatialDistributionMap(media, mesh);
+
// Initial stress conditions
auto const initial_stress = ParameterLib::findOptionalTagParameter<double>(
//! \ogs_file_param_special{prj__processes__process__HYDRO_MECHANICS__initial_stress}
@@ -217,10 +207,10 @@ std::unique_ptr<Process> createHydroMechanicsProcess(
}
HydroMechanicsProcessData<DisplacementDim> process_data{
- materialIDs(mesh), std::move(solid_constitutive_relations),
+ materialIDs(mesh), std::move(media_map),
+ std::move(solid_constitutive_relations),
initial_stress, intrinsic_permeability,
- fluid_viscosity, fluid_density,
- biot_coefficient, porosity,
+ fluid_density, biot_coefficient,
solid_density, specific_body_force,
fluid_compressibility, reference_temperature,
specific_gas_constant, fluid_type};
@@ -245,7 +235,8 @@ template std::unique_ptr<Process> createHydroMechanicsProcess<2>(
boost::optional<ParameterLib::CoordinateSystem> const&
local_coordinate_system,
unsigned const integration_order,
- BaseLib::ConfigTree const& config);
+ BaseLib::ConfigTree const& config,
+ std::map<int, std::unique_ptr<MaterialPropertyLib::Medium>> const& media);
template std::unique_ptr<Process> createHydroMechanicsProcess<3>(
std::string name,
@@ -256,7 +247,8 @@ template std::unique_ptr<Process> createHydroMechanicsProcess<3>(
boost::optional<ParameterLib::CoordinateSystem> const&
local_coordinate_system,
unsigned const integration_order,
- BaseLib::ConfigTree const& config);
+ BaseLib::ConfigTree const& config,
+ std::map<int, std::unique_ptr<MaterialPropertyLib::Medium>> const& media);
} // namespace HydroMechanics
} // namespace ProcessLib
diff --git a/ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.h b/ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.h
index 3119ebd8d28..7fb0b496022 100644
--- a/ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.h
+++ b/ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.h
@@ -11,6 +11,7 @@
#pragma once
#include <boost/optional.hpp>
+#include <map>
#include <memory>
#include <vector>
@@ -18,6 +19,10 @@ namespace BaseLib
{
class ConfigTree;
}
+namespace MaterialPropertyLib
+{
+class Medium;
+}
namespace MeshLib
{
class Mesh;
@@ -48,7 +53,8 @@ std::unique_ptr<Process> createHydroMechanicsProcess(
boost::optional<ParameterLib::CoordinateSystem> const&
local_coordinate_system,
unsigned const integration_order,
- BaseLib::ConfigTree const& config);
+ BaseLib::ConfigTree const& config,
+ std::map<int, std::unique_ptr<MaterialPropertyLib::Medium>> const& media);
extern template std::unique_ptr<Process> createHydroMechanicsProcess<2>(
std::string name,
@@ -59,7 +65,8 @@ extern template std::unique_ptr<Process> createHydroMechanicsProcess<2>(
boost::optional<ParameterLib::CoordinateSystem> const&
local_coordinate_system,
unsigned const integration_order,
- BaseLib::ConfigTree const& config);
+ BaseLib::ConfigTree const& config,
+ std::map<int, std::unique_ptr<MaterialPropertyLib::Medium>> const& media);
extern template std::unique_ptr<Process> createHydroMechanicsProcess<3>(
std::string name,
@@ -70,6 +77,7 @@ extern template std::unique_ptr<Process> createHydroMechanicsProcess<3>(
boost::optional<ParameterLib::CoordinateSystem> const&
local_coordinate_system,
unsigned const integration_order,
- BaseLib::ConfigTree const& config);
+ BaseLib::ConfigTree const& config,
+ std::map<int, std::unique_ptr<MaterialPropertyLib::Medium>> const& media);
} // namespace HydroMechanics
} // namespace ProcessLib
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h b/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
index 5d542a5f981..fe040e8db2b 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
@@ -14,6 +14,8 @@
#include "HydroMechanicsFEM.h"
#include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
+#include "MaterialLib/MPL/Medium.h"
+#include "MaterialLib/MPL/Property.h"
#include "MathLib/KelvinVector.h"
#include "NumLib/Function/Interpolation.h"
#include "ProcessLib/CoupledSolutionsForStaggeredScheme.h"
@@ -180,6 +182,11 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
unsigned const n_integration_points =
_integration_method.getNumberOfPoints();
+
+ auto const& medium = _process_data.media_map->getMedium(_element.getID());
+ auto const& gas_phase = medium->phase("Gas");
+ auto const& solid_phase = medium->phase("Solid");
+ MaterialPropertyLib::VariableArray vars;
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
x_position.setIntegrationPoint(ip);
@@ -206,9 +213,13 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
auto& eps = _ip_data[ip].eps;
auto const& sigma_eff = _ip_data[ip].sigma_eff;
+ auto const mu =
+ gas_phase.property(MaterialPropertyLib::PropertyType::viscosity)
+ .template value<double>(vars, x_position, t);
+
double const K_over_mu =
- _process_data.intrinsic_permeability(t, x_position)[0] /
- _process_data.fluid_viscosity(t, x_position)[0];
+ _process_data.intrinsic_permeability(t, x_position)[0] / mu;
+
auto const alpha = _process_data.biot_coefficient(t, x_position)[0];
auto const K_S = solid_material.getBulkModulus(t, x_position);
auto const rho_sr = _process_data.solid_density(t, x_position)[0];
@@ -219,7 +230,9 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
: std::numeric_limits<double>::quiet_NaN();
double const rho_fr = _process_data.getFluidDensity(t, x_position, p_fr);
double const beta_p = _process_data.getFluidCompressibility(p_fr);
- auto const porosity = _process_data.porosity(t, x_position)[0];
+ auto const porosity =
+ solid_phase.property(MaterialPropertyLib::PropertyType::porosity)
+ .template value<double>(vars, x_position, t);
auto const& identity2 = MathLib::KelvinVector::Invariants<
MathLib::KelvinVector::KelvinVectorDimensions<
DisplacementDim>::value>::identity2;
@@ -326,12 +339,21 @@ HydroMechanicsLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
ParameterLib::SpatialPosition x_position;
x_position.setElementID(_element.getID());
+
+ auto const& medium = _process_data.media_map->getMedium(_element.getID());
+ auto const& gas_phase = medium->phase("Gas");
+ MaterialPropertyLib::VariableArray vars;
+
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
x_position.setIntegrationPoint(ip);
+
+ auto const mu =
+ gas_phase.property(MaterialPropertyLib::PropertyType::viscosity)
+ .template value<double>(vars, x_position, t);
+
double const K_over_mu =
- _process_data.intrinsic_permeability(t, x_position)[0] /
- _process_data.fluid_viscosity(t, x_position)[0];
+ _process_data.intrinsic_permeability(t, x_position)[0] / mu;
auto const rho_fr = _process_data.fluid_density(t, x_position)[0];
auto const& b = _process_data.specific_body_force;
@@ -391,6 +413,11 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
ParameterLib::SpatialPosition x_position;
x_position.setElementID(_element.getID());
+ auto const& medium = _process_data.media_map->getMedium(_element.getID());
+ auto const& gas_phase = medium->phase("Gas");
+ auto const& solid_phase = medium->phase("Solid");
+ MaterialPropertyLib::VariableArray vars;
+
int const n_integration_points = _integration_method.getNumberOfPoints();
for (int ip = 0; ip < n_integration_points; ip++)
{
@@ -403,9 +430,12 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
auto const& N_p = _ip_data[ip].N_p;
auto const& dNdx_p = _ip_data[ip].dNdx_p;
+ auto const mu =
+ gas_phase.property(MaterialPropertyLib::PropertyType::viscosity)
+ .template value<double>(vars, x_position, t);
+
double const K_over_mu =
- _process_data.intrinsic_permeability(t, x_position)[0] /
- _process_data.fluid_viscosity(t, x_position)[0];
+ _process_data.intrinsic_permeability(t, x_position)[0] / mu;
auto const alpha_b = _process_data.biot_coefficient(t, x_position)[0];
// TODO (FZill) get fluid properties from GPML
double const p_fr =
@@ -414,7 +444,9 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
: std::numeric_limits<double>::quiet_NaN();
double const rho_fr = _process_data.getFluidDensity(t, x_position, p_fr);
double const beta_p = _process_data.getFluidCompressibility(p_fr);
- auto const porosity = _process_data.porosity(t, x_position)[0];
+ auto const porosity =
+ solid_phase.property(MaterialPropertyLib::PropertyType::porosity)
+ .template value<double>(vars, x_position, t);
auto const K_S = solid_material.getBulkModulus(t, x_position);
laplace.noalias() += dNdx_p.transpose() * K_over_mu * dNdx_p * w;
@@ -477,6 +509,10 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
ParameterLib::SpatialPosition x_position;
x_position.setElementID(_element.getID());
+ auto const& medium = _process_data.media_map->getMedium(_element.getID());
+ auto const& solid_phase = medium->phase("Solid");
+ MaterialPropertyLib::VariableArray vars;
+
int const n_integration_points = _integration_method.getNumberOfPoints();
for (int ip = 0; ip < n_integration_points; ip++)
{
@@ -506,7 +542,9 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
auto const alpha = _process_data.biot_coefficient(t, x_position)[0];
auto const rho_sr = _process_data.solid_density(t, x_position)[0];
auto const rho_fr = _process_data.fluid_density(t, x_position)[0];
- auto const porosity = _process_data.porosity(t, x_position)[0];
+ auto const porosity =
+ solid_phase.property(MaterialPropertyLib::PropertyType::porosity)
+ .template value<double>(vars, x_position, t);
auto const& b = _process_data.specific_body_force;
auto const& identity2 = MathLib::KelvinVector::Invariants<
MathLib::KelvinVector::KelvinVectorDimensions<
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsProcessData.h b/ProcessLib/HydroMechanics/HydroMechanicsProcessData.h
index 66db57df0d8..7f1d59f54f9 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsProcessData.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsProcessData.h
@@ -18,6 +18,8 @@
#include <Eigen/Dense>
+#include "MaterialLib/MPL/MaterialSpatialDistributionMap.h"
+
namespace MaterialLib
{
namespace Solids
@@ -35,6 +37,9 @@ struct HydroMechanicsProcessData
{
MeshLib::PropertyVector<int> const* const material_ids = nullptr;
+ std::unique_ptr<MaterialPropertyLib::MaterialSpatialDistributionMap>
+ media_map = nullptr;
+
/// The constitutive relation for the mechanical part.
/// \note Linear elasticity is the only supported one in the moment.
std::map<
@@ -49,15 +54,10 @@ struct HydroMechanicsProcessData
/// Permeability of the solid. A scalar quantity,
/// ParameterLib::Parameter<double>.
ParameterLib::Parameter<double> const& intrinsic_permeability;
- /// Fluid's viscosity. A scalar quantity, ParameterLib::Parameter<double>.
- ParameterLib::Parameter<double> const& fluid_viscosity;
/// Fluid's density. A scalar quantity, ParameterLib::Parameter<double>.
ParameterLib::Parameter<double> const& fluid_density;
/// Biot coefficient. A scalar quantity, ParameterLib::Parameter<double>.
ParameterLib::Parameter<double> const& biot_coefficient;
- /// Porosity of the solid. A scalar quantity,
- /// ParameterLib::Parameter<double>.
- ParameterLib::Parameter<double> const& porosity;
/// Solid's density. A scalar quantity, ParameterLib::Parameter<double>.
ParameterLib::Parameter<double> const& solid_density;
/// Specific body forces applied to solid and fluid.
--
GitLab