From c9afadc9ecaddf708488d80e8648a5dd7598f986 Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Wed, 12 Oct 2016 15:51:26 +0200
Subject: [PATCH] [Mat] Alternative implementation of porous medium properties.
---
.../LiquidFlowMaterialProperties.cpp | 75 +++++++++++++++----
.../LiquidFlow/LiquidFlowMaterialProperties.h | 62 +++++++++------
2 files changed, 99 insertions(+), 38 deletions(-)
diff --git a/ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.cpp b/ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.cpp
index 265f8d8b23a..36e8a3bd42e 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.cpp
+++ b/ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.cpp
@@ -15,6 +15,10 @@
#include <logog/include/logog.hpp>
#include "MeshLib/Mesh.h"
+#include "MeshLib/PropertyVector.h"
+
+#include "ProcessLib/Parameter/Parameter.h"
+#include "ProcessLib/Parameter/SpatialPosition.h"
#include "MaterialLib/Fluid/FluidProperty.h"
#include "MaterialLib/PorousMedium/Porosity/Porosity.h"
@@ -25,7 +29,15 @@ namespace ProcessLib
namespace LiquidFlow
{
LiquidFlowMaterialProperties::LiquidFlowMaterialProperties(
- BaseLib::ConfigTree const& config)
+ 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)
{
DBUG("Reading material properties of liquid flow process.");
@@ -35,10 +47,10 @@ LiquidFlowMaterialProperties::LiquidFlowMaterialProperties(
// Get fluid properties
//! \ogs_file_param{prj__material_property__fluid__density}
auto const& rho_conf = fluid_config.getConfigSubtree("density");
- liquid_density = MaterialLib::Fluid::createFluidDensityModel(rho_conf);
+ _liquid_density = MaterialLib::Fluid::createFluidDensityModel(rho_conf);
//! \ogs_file_param{prj__material_property__fluid__viscosity}
auto const& mu_conf = fluid_config.getConfigSubtree("viscosity");
- viscosity = MaterialLib::Fluid::createViscosityModel(mu_conf);
+ _viscosity = MaterialLib::Fluid::createViscosityModel(mu_conf);
// Get porous properties
//! \ogs_file_param{prj__material_property__porous_medium}
@@ -48,18 +60,18 @@ LiquidFlowMaterialProperties::LiquidFlowMaterialProperties(
{
//! \ogs_file_param{prj__material_property__porous_medium__porous_medium__permeability}
auto const& perm_conf = conf.getConfigSubtree("permeability");
- intrinsic_permeability.emplace_back(
+ _intrinsic_permeability_models.emplace_back(
MaterialLib::PorousMedium::createPermeabilityModel(perm_conf));
//! \ogs_file_param{prj__material_property__porous_medium__porous_medium__porosity}
auto const& poro_conf = conf.getConfigSubtree("porosity");
auto n = MaterialLib::PorousMedium::createPorosityModel(poro_conf);
- porosity.emplace_back(std::move(n));
+ _porosity_models.emplace_back(std::move(n));
//! \ogs_file_param{prj__material_property__porous_medium__porous_medium__storage}
auto const& stora_conf = conf.getConfigSubtree("storage");
auto beta = MaterialLib::PorousMedium::createStorageModel(stora_conf);
- storage.emplace_back(std::move(beta));
+ _storage_models.emplace_back(std::move(beta));
}
}
@@ -69,7 +81,7 @@ double LiquidFlowMaterialProperties::getLiquidDensity(const double p,
ArrayType vars;
vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::T)] = T;
vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::pl)] = p;
- return liquid_density->getValue(vars);
+ return _liquid_density->getValue(vars);
}
double LiquidFlowMaterialProperties::getViscosity(const double p,
@@ -78,22 +90,53 @@ double LiquidFlowMaterialProperties::getViscosity(const double p,
ArrayType vars;
vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::T)] = T;
vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::pl)] = p;
- return viscosity->getValue(vars);
+ return _viscosity->getValue(vars);
}
double LiquidFlowMaterialProperties::getMassCoefficient(
- const double porosity_variable, const double storage_variable,
- const double p, const double T, const unsigned material_group_id) 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;
- return porosity[material_group_id]->getValue(porosity_variable, T) *
- liquid_density->getdValue(
- vars, MaterialLib::Fluid::PropertyVariableType::pl)
- // Divided by rho_l because the PDE is scaled with rho_l
- / liquid_density->getValue(vars) +
- storage[material_group_id]->getValue(storage_variable);
+ 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;
+ }
+}
+
+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);
+ }
}
} // end of namespace
diff --git a/ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.h b/ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.h
index 3a3f6d9696d..0c03a17c4de 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.h
@@ -18,13 +18,6 @@
#include "MaterialLib/Fluid/FluidPropertyHeaders.h"
#include "MaterialLib/PorousMedium/PorousPropertyHeaders.h"
-#include "MeshLib/PropertyVector.h"
-
-namespace MeshLib
-{
-class Mesh;
-}
-
namespace MaterialLib
{
namespace PorousMedium
@@ -34,15 +27,29 @@ class Storage;
}
}
+namespace MeshLib
+{
+template <typename PROP_VAL_TYPE> class PropertyVector;
+}
+
namespace ProcessLib
{
+template <typename T> struct Parameter;
+class SpatialPosition;
+
namespace LiquidFlow
{
-struct LiquidFlowMaterialProperties
+class LiquidFlowMaterialProperties
{
+public:
typedef MaterialLib::Fluid::FluidProperty::ArrayType ArrayType;
- explicit LiquidFlowMaterialProperties(BaseLib::ConfigTree const& config);
+ 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);
/**
* \brief Compute the coefficient of the mass term by
@@ -51,34 +58,45 @@ struct LiquidFlowMaterialProperties
* \f]
* where \f$n\f$ is the porosity, \f$rho_l\f$ is the liquid density,
* \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.
* \param storage_variable Variable for storage model.
* \param p Pressure value
* \param T Temperature value
- * \param material_group_id Material ID of the element
*/
- double getMassCoefficient(const double p, const double T,
+ double getMassCoefficient(const double t, const SpatialPosition& pos,
+ const double p, const double T,
const double porosity_variable,
- const double storage_variable,
- const unsigned material_group_id = 0) const;
+ const double storage_variable) const;
+
+ Eigen::MatrixXd const& getPermeability(const double t,
+ const SpatialPosition& pos,
+ const int dim) const;
double getLiquidDensity(const double p, const double T) const;
double getViscosity(const double p, const double T) const;
- std::unique_ptr<MaterialLib::Fluid::FluidProperty> liquid_density;
- std::unique_ptr<MaterialLib::Fluid::FluidProperty> viscosity;
+private:
+ std::unique_ptr<MaterialLib::Fluid::FluidProperty> _liquid_density;
+ std::unique_ptr<MaterialLib::Fluid::FluidProperty> _viscosity;
- /// Porous medium properties of different material zones.
- /// The vector is left empty if the property data are given in vtu file,
- /// e.g for heterogeneous medium.
- std::vector<Eigen::MatrixXd> intrinsic_permeability;
- std::vector<std::unique_ptr<MaterialLib::PorousMedium::Porosity>> porosity;
- std::vector<std::unique_ptr<MaterialLib::PorousMedium::Storage>> storage;
+ /** Use porous medium models for different material zones.
+ * 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;
- // TODO: heterogeneous medium.
+ /// Use data for porous medium properties.
+ Parameter<double> const& _intrinsic_permeability_data;
+ Parameter<double> const& _porosity_data;
+ Parameter<double> const& _storage_data;
};
} // end of namespace
--
GitLab