diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
index 786253238fa852f6652e08314ed4966a2e47e5ed..e4e26e9910d200946f181cb1dfdc6bd4f763077e 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
@@ -25,7 +25,7 @@ namespace LiquidFlow
template <typename ShapeFunction, typename IntegrationMethod,
unsigned GlobalDim>
void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
- assemble(double const /*t*/, std::vector<double> const& local_x,
+ assemble(double const t, std::vector<double> const& local_x,
std::vector<double>& local_M_data,
std::vector<double>& local_K_data,
std::vector<double>& local_b_data)
@@ -43,9 +43,12 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
unsigned const n_integration_points =
_integration_method.getNumberOfPoints();
- const unsigned mat_id = 0; // TODO for heterogeneous medium
+ SpatialPosition pos;
+ pos.setElementID(_element.getID());
+ _material_properties.setMaterialID(pos);
+
const Eigen::MatrixXd& perm =
- _material_properties.intrinsic_permeability[mat_id];
+ _material_properties.getPermeability(t, pos, _element.getDimension());
// 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()
@@ -69,8 +72,8 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
// Assemble mass matrix, M
local_M.noalias() +=
- _material_properties.getMassCoefficient(porosity_variable, storage_variable,
- p, _temperature, mat_id) *
+ _material_properties.getMassCoefficient(
+ t, pos, porosity_variable, storage_variable, p, _temperature) *
sm.N.transpose() * sm.N * integration_factor;
// Compute density:
@@ -81,7 +84,7 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
const double mu = _material_properties.getViscosity(p, _temperature);
// Assemble Laplacian, K, and RHS by the gravitational term
- if (perm.size() == 1) //Save the computing time for isotropic permeability.
+ if (perm.size() == 1)
{
// Use scalar number for isotropic permeability
// to save the computation time.
diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
index 4d92ebb247a6bb5ccb160f1a037d11a3944960b4..acefca848bfc84a8b74894ebae47928fcde9dc28 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
@@ -62,13 +62,12 @@ class LiquidFlowLocalAssembler : public LiquidFlowLocalAssemblerInterface
using GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType;
public:
- LiquidFlowLocalAssembler(
- MeshLib::Element const& element,
- std::size_t const /*local_matrix_size*/,
- bool const is_axially_symmetric,
- unsigned const integration_order,
- bool const compute_gravitational_term,
- LiquidFlowMaterialProperties const& material_propertries)
+ LiquidFlowLocalAssembler(MeshLib::Element const& element,
+ std::size_t const /*local_matrix_size*/,
+ bool const is_axially_symmetric,
+ unsigned const integration_order,
+ bool const compute_gravitational_term,
+ LiquidFlowMaterialProperties& material_propertries)
: _element(element),
_integration_method(integration_order),
_shape_matrices(initShapeMatrices<ShapeFunction, ShapeMatricesType,
@@ -126,7 +125,7 @@ private:
std::vector<double>(_integration_method.getNumberOfPoints()));
const bool _compute_gravitational_term;
- LiquidFlowMaterialProperties const& _material_properties;
+ LiquidFlowMaterialProperties& _material_properties;
double _temperature;
};
diff --git a/ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.cpp b/ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.cpp
index 36e8a3bd42e7094900dae734116a66a83076901b..b6d4d138f5d9fc46bbf2f1746b9e49f7fcc55b27 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.cpp
+++ b/ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.cpp
@@ -29,15 +29,9 @@ namespace ProcessLib
namespace LiquidFlow
{
LiquidFlowMaterialProperties::LiquidFlowMaterialProperties(
- BaseLib::ConfigTree const& config,
- MeshLib::PropertyVector<int> const& material_ids,
- Parameter<double> const& intrinsic_permeability_data,
- Parameter<double> const& porosity_data,
- Parameter<double> const& storage_data)
- : _material_ids(material_ids),
- _intrinsic_permeability_data(intrinsic_permeability_data),
- _porosity_data(porosity_data),
- _storage_data(storage_data)
+ BaseLib::ConfigTree const& config,
+ MeshLib::PropertyVector<int> const& material_ids)
+ : _material_ids(material_ids)
{
DBUG("Reading material properties of liquid flow process.");
@@ -75,6 +69,15 @@ LiquidFlowMaterialProperties::LiquidFlowMaterialProperties(
}
}
+void LiquidFlowMaterialProperties::setMaterialID(const SpatialPosition& pos)
+{
+ if (_material_ids.empty())
+ {
+ assert(pos.getElementID().get() < _material_ids.size());
+ _current_material_id = _material_ids[pos.getElementID().get()];
+ }
+}
+
double LiquidFlowMaterialProperties::getLiquidDensity(const double p,
const double T) const
{
@@ -94,49 +97,28 @@ double LiquidFlowMaterialProperties::getViscosity(const double p,
}
double LiquidFlowMaterialProperties::getMassCoefficient(
- const double t, const SpatialPosition& pos,
- const double p, const double T,
- const double porosity_variable,
- const double storage_variable) const
+ const double /*t*/, const SpatialPosition& /*pos*/, const double p,
+ const double T, const double porosity_variable,
+ const double storage_variable) const
{
ArrayType vars;
vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::T)] = T;
vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::pl)] = p;
- const double drho_dp =
- _liquid_density->getdValue(vars,
- MaterialLib::Fluid::PropertyVariableType::pl);
+ const double drho_dp = _liquid_density->getdValue(
+ vars, MaterialLib::Fluid::PropertyVariableType::pl);
const double rho = _liquid_density->getValue(vars);
- if (_storage_models.size() > 0)
- {
- const int mat_id = _material_ids[pos.getElementID()];
- const double porosity = _porosity_models[mat_id]->getValue(porosity_variable, T);
- const double storage = _storage_models[mat_id]->getValue(storage_variable);
- return porosity * drho_dp / rho + storage;
- }
- else
- {
- const double storage = _storage_data(t, pos)[0];
- const double porosity = _porosity_data(t, pos)[0];
- return porosity * drho_dp / rho + storage;
- }
+ const double porosity =
+ _porosity_models[_current_material_id]->getValue(porosity_variable, T);
+ const double storage =
+ _storage_models[_current_material_id]->getValue(storage_variable);
+ return porosity * drho_dp / rho + storage;
}
-Eigen::MatrixXd const& LiquidFlowMaterialProperties
- ::getPermeability(const double t,
- const SpatialPosition& pos,
- const int dim) const
+Eigen::MatrixXd const& LiquidFlowMaterialProperties::getPermeability(
+ const double /*t*/, const SpatialPosition& /*pos*/, const int /*dim*/) const
{
- if (_intrinsic_permeability_models.size() > 0)
- {
- const int mat_id = _material_ids[pos.getElementID()];
- return _intrinsic_permeability_models[mat_id];
- }
- else
- {
- auto const permeability = _intrinsic_permeability_data(t, pos)[0];
- return MathLib::toMatrix(permeability, dim, dim);
- }
+ return _intrinsic_permeability_models[_current_material_id];
}
} // end of namespace
diff --git a/ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.h b/ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.h
index 0c03a17c4ded8a78097674f1a038366792dcc648..dc59a38c8bf4fc669553cdfb62e723adf382de76 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.h
@@ -29,12 +29,12 @@ class Storage;
namespace MeshLib
{
-template <typename PROP_VAL_TYPE> class PropertyVector;
+template <typename PROP_VAL_TYPE>
+class PropertyVector;
}
namespace ProcessLib
{
-template <typename T> struct Parameter;
class SpatialPosition;
namespace LiquidFlow
@@ -45,11 +45,10 @@ public:
typedef MaterialLib::Fluid::FluidProperty::ArrayType ArrayType;
LiquidFlowMaterialProperties(
- BaseLib::ConfigTree const& config,
- MeshLib::PropertyVector<int> const& material_ids,
- Parameter<double> const& intrinsic_permeability_data,
- Parameter<double> const& porosity_data,
- Parameter<double> const& storage_data);
+ BaseLib::ConfigTree const& config,
+ MeshLib::PropertyVector<int> const& material_ids);
+
+ void setMaterialID(const SpatialPosition& pos);
/**
* \brief Compute the coefficient of the mass term by
@@ -60,7 +59,6 @@ public:
* \f$bata_s\f$ is the storage.
* \param t Time.
* \param pos Position of element.
- * \param dim Dimension of space.
* \param porosity_variable The first variable for porosity model, and it
* passes a double type value that could be
* saturation, and invariant of stress or strain.
@@ -74,8 +72,8 @@ public:
const double storage_variable) const;
Eigen::MatrixXd const& getPermeability(const double t,
- const SpatialPosition& pos,
- const int dim) const;
+ const SpatialPosition& pos,
+ const int dim) const;
double getLiquidDensity(const double p, const double T) const;
@@ -89,14 +87,13 @@ private:
* Material IDs must be given as mesh element properties.
*/
MeshLib::PropertyVector<int> const& _material_ids;
- std::vector<Eigen::MatrixXd> _intrinsic_permeability_models;
- std::vector<std::unique_ptr<MaterialLib::PorousMedium::Porosity>> _porosity_models;
- std::vector<std::unique_ptr<MaterialLib::PorousMedium::Storage>> _storage_models;
- /// Use data for porous medium properties.
- Parameter<double> const& _intrinsic_permeability_data;
- Parameter<double> const& _porosity_data;
- Parameter<double> const& _storage_data;
+ int _current_material_id = 0;
+ std::vector<Eigen::MatrixXd> _intrinsic_permeability_models;
+ std::vector<std::unique_ptr<MaterialLib::PorousMedium::Porosity>>
+ _porosity_models;
+ std::vector<std::unique_ptr<MaterialLib::PorousMedium::Storage>>
+ _storage_models;
};
} // end of namespace
diff --git a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
index c59dfdf12965ab5f6643ac8e583260a13a11a721..aae82265449b4e97f3ec2acfd97c142c6f81c292 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
+++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
@@ -14,6 +14,8 @@
#include <cassert>
+#include "MeshLib/PropertyVector.h"
+
#include "ProcessLib/Utils/CreateLocalAssemblers.h"
#include "LiquidFlowLocalAssembler.h"
@@ -32,14 +34,14 @@ LiquidFlowProcess::LiquidFlowProcess(
std::vector<std::reference_wrapper<ProcessVariable>>&& process_variables,
SecondaryVariableCollection&& secondary_variables,
NumLib::NamedFunctionCaller&& named_function_caller,
+ MeshLib::PropertyVector<int> const& material_ids,
bool const compute_gravitational_term,
BaseLib::ConfigTree const& config)
: Process(mesh, std::move(jacobian_assembler), parameters,
- integration_order,
- std::move(process_variables), std::move(secondary_variables),
- std::move(named_function_caller)),
+ integration_order, std::move(process_variables),
+ std::move(secondary_variables), std::move(named_function_caller)),
_compute_gravitational_term(compute_gravitational_term),
- _material_properties(LiquidFlowMaterialProperties(config))
+ _material_properties(LiquidFlowMaterialProperties(config, material_ids))
{
DBUG("Create Liquid flow process.");
}
diff --git a/ProcessLib/LiquidFlow/LiquidFlowProcess.h b/ProcessLib/LiquidFlow/LiquidFlowProcess.h
index f19556e1294db09962b02974f2c23f2307d1b608..f5b9bbbdf79499144cdcd8a994c16793b24cda96 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowProcess.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.h
@@ -63,6 +63,7 @@ public:
process_variables,
SecondaryVariableCollection&& secondary_variables,
NumLib::NamedFunctionCaller&& named_function_caller,
+ MeshLib::PropertyVector<int> const& material_ids,
bool const compute_gravitational_term,
BaseLib::ConfigTree const& config);
diff --git a/ProcessLib/LiquidFlow/createLiquidFlowProcess.cpp b/ProcessLib/LiquidFlow/createLiquidFlowProcess.cpp
index bf4fd29b84401202cb6d49d4e97850916e0fc7b0..89fe6a2538161092bcb51595e03dffc3e4eb5a02 100644
--- a/ProcessLib/LiquidFlow/createLiquidFlowProcess.cpp
+++ b/ProcessLib/LiquidFlow/createLiquidFlowProcess.cpp
@@ -11,6 +11,8 @@
*/
#include "createLiquidFlowProcess.h"
+#include "MeshLib/MeshGenerators/MeshGenerator.h"
+
#include "ProcessLib/Utils/ParseSecondaryVariables.h"
#include "ProcessLib/Utils/ProcessUtils.h"
@@ -54,11 +56,31 @@ std::unique_ptr<Process> createLiquidFlowProcess(
//! \ogs_file_param{process__LIQUID_FLOW__material_property}
auto const& mat_config = config.getConfigSubtree("material_property");
- return std::unique_ptr<Process>{new LiquidFlowProcess{
- mesh, std::move(jacobian_assembler), parameters,
- integration_order,
- std::move(process_variables), std::move(secondary_variables),
- std::move(named_function_caller), has_gravitational_term, mat_config}};
+ auto const& mat_ids =
+ mesh.getProperties().getPropertyVector<int>("MaterialIDs");
+ if (mat_ids)
+ {
+ INFO("The liquid flow is in heterogeneous porous media.");
+ return std::unique_ptr<Process>{new LiquidFlowProcess{
+ mesh, std::move(jacobian_assembler), parameters, integration_order,
+ std::move(process_variables), std::move(secondary_variables),
+ std::move(named_function_caller), mat_ids.get(),
+ has_gravitational_term, mat_config}};
+ }
+ else
+ {
+ INFO("The liquid flow is in homogeneous porous media.");
+
+ MeshLib::Properties dummy_property;
+ auto const& dummy_property_vector =
+ dummy_property.createNewPropertyVector<int>(
+ "MaterialIDs", MeshLib::MeshItemType::Cell, 1);
+ return std::unique_ptr<Process>{new LiquidFlowProcess{
+ mesh, std::move(jacobian_assembler), parameters, integration_order,
+ std::move(process_variables), std::move(secondary_variables),
+ std::move(named_function_caller), dummy_property_vector.get(),
+ has_gravitational_term, mat_config}};
+ }
}
} // end of namespace
diff --git a/Tests/ProcessLib/LiquidFlow/TestLiquidFlowMaterialProperties.cpp b/Tests/ProcessLib/LiquidFlow/TestLiquidFlowMaterialProperties.cpp
index f948a6bcc5e7226c4147fff6a67f31102f82db91..677d5d80a7c473da5a5579c88a02281fed67c264 100644
--- a/Tests/ProcessLib/LiquidFlow/TestLiquidFlowMaterialProperties.cpp
+++ b/Tests/ProcessLib/LiquidFlow/TestLiquidFlowMaterialProperties.cpp
@@ -11,12 +11,17 @@
*/
#include <gtest/gtest.h>
+#include <memory>
#include "TestTools.h"
+#include "MeshLib/Mesh.h"
+#include "MeshLib/MeshGenerators/MeshGenerator.h"
+
+#include "ProcessLib/Parameter/SpatialPosition.h"
#include "ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.h"
-#include "MaterialLib/Fluid/FluidProperty.h"
+#include "MaterialLib/Fluid/FluidProperty.h"
#include "MaterialLib/PorousMedium/Porosity/Porosity.h"
#include "MaterialLib/PorousMedium/Storage/Storage.h"
@@ -67,41 +72,20 @@ TEST(ProcessLibLiquidFlow, checkLiquidFlowMaterialProperties)
BaseLib::ConfigTree::onwarning);
auto const& sub_config = conf.getConfigSubtree("material_property");
- LiquidFlowMaterialProperties lprop(sub_config);
-
- // Check density
- const ArrayType vars = {273.15 + 60.0, 1.e+6};
- const double T0 = 273.15;
- const double p0 = 1.e+5;
- const double rho0 = 999.8;
- const double K = 2.15e+9;
- const double beta = 2.e-4;
- const double T = vars[0];
- const double p = vars[1];
+ std::unique_ptr<MeshLib::Mesh> mesh(MeshLib::MeshGenerator::generateLineMesh(1u, 1.0));
+ std::vector<int> material_ids({0});
+ MeshLib::addPropertyToMesh(*mesh, "MaterialIDs", MeshLib::MeshItemType::Cell,
+ 1, material_ids);
+
+ auto const& mat_ids = mesh->getProperties().getPropertyVector<int>("MaterialIDs");
- const double fac_T = 1. + beta * (T - T0);
- ASSERT_NEAR(rho0 / fac_T / (1. - (p - p0) / K),
- lprop.liquid_density->getValue(vars), 1.e-10);
+ LiquidFlowMaterialProperties lprop(sub_config, mat_ids.get());
- // Test the derivative with respect to temperature.
- ASSERT_NEAR(-beta * rho0 / (fac_T * fac_T) / (1. - (p - p0) / K),
- lprop.liquid_density->getdValue(vars, PropertyVariableType::T), 1.e-10);
-
- // Test the derivative with respect to pressure.
- const double fac_p = 1. - (p - p0) / K;
- ASSERT_NEAR(rho0 / (1. + beta * (T - T0)) / (fac_p * fac_p * K),
- lprop.liquid_density->getdValue(vars, PropertyVariableType::pl), 1.e-10);
-
- // Check viscosity
- ArrayType vars1;
- vars1[0] = 303.0;
- const auto var_type = MaterialLib::Fluid::PropertyVariableType::T;
- ASSERT_NEAR(0.802657e-3, lprop.viscosity->getValue(vars1), 1.e-5);
- ASSERT_NEAR(-1.87823e-5, lprop.viscosity->getdValue(vars1, var_type),
- 1.e-5);
+ ProcessLib::SpatialPosition pos;
+ pos.setElementID(0);
// Check permeability
- Eigen::MatrixXd& perm = lprop.intrinsic_permeability[0];
+ const Eigen::MatrixXd& perm = lprop.getPermeability(0., pos, 1);
ASSERT_EQ(2.e-10, perm(0, 0));
ASSERT_EQ(0., perm(0, 1));
ASSERT_EQ(0., perm(0, 2));
@@ -111,12 +95,9 @@ TEST(ProcessLibLiquidFlow, checkLiquidFlowMaterialProperties)
ASSERT_EQ(0., perm(2, 0));
ASSERT_EQ(0., perm(2, 1));
ASSERT_EQ(4.e-10, perm(2, 2));
-
- // Check porosity
- const double variable = 0.;
- const double temperature = 0.;
- ASSERT_EQ(0.2, lprop.porosity[0]->getValue(variable, temperature));
-
- // Check storage
- ASSERT_EQ(1.e-4, lprop.storage[0]->getValue(variable));
+
+ const double T = 273.15 + 60.0;
+ const double p = 1.e+6;
+ const double mass_coef = lprop.getMassCoefficient(0., pos, p, T, 0., 0.);
+ ASSERT_NEAR(0.000100000093, mass_coef, 1.e-10);
}