diff --git a/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/c_RICHARDS_FLOW.md b/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/c_RICHARDS_FLOW.md
index 576add64edccf5d980a0ff86143dc43a7109a788..f5925e17d6fb4fd5bcaa304e814d1bf3f53731c8 100644
--- a/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/c_RICHARDS_FLOW.md
+++ b/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/c_RICHARDS_FLOW.md
@@ -1 +1 @@
-\ogs_missing_documentation
+\copydoc ProcessLib::RichardsFlow::RichardsFlowProcess
diff --git a/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/material_property/fluid/i_fluid.md b/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/material_property/fluid/i_fluid.md
new file mode 100644
index 0000000000000000000000000000000000000000..00aef4b95e01a85c9fb21ee4e2f36860cefd7dc1
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/material_property/fluid/i_fluid.md
@@ -0,0 +1 @@
+Defines fluid properties.
diff --git a/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/material_property/i_material_property.md b/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/material_property/i_material_property.md
new file mode 100644
index 0000000000000000000000000000000000000000..c0c96b339623410ba707bd6ddf4a2037c07162d5
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/material_property/i_material_property.md
@@ -0,0 +1 @@
+A tag for material properties of the RichardsFlow process.
diff --git a/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/material_property/porous_medium/i_porous_medium.md b/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/material_property/porous_medium/i_porous_medium.md
new file mode 100644
index 0000000000000000000000000000000000000000..0ff6eac587d8f6af9cd9f8087fcdb7ab184d6bc7
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/material_property/porous_medium/i_porous_medium.md
@@ -0,0 +1,13 @@
+A tag for the properties of porous media, which contains nested tags for
+different porous media, which are distinguished by material IDs. For example
+
+```
+ <porous_medium>
+ <porous_medium id="0">
+ ...
+ </porous_medium>
+ <porous_medium id="1">
+ ...
+ </porous_medium>
+ </porous_medium>
+```
\ No newline at end of file
diff --git a/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/material_property/porous_medium/porous_medium/a_id.md b/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/material_property/porous_medium/porous_medium/a_id.md
new file mode 100644
index 0000000000000000000000000000000000000000..c1dae9f80b40a7c5ba5eb2fc7fae385bf250a298
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/material_property/porous_medium/porous_medium/a_id.md
@@ -0,0 +1,2 @@
+It specifies the material ID of a porous medium. The material IDs of elements
+must be given in the mesh's cell data array "MaterialIDs".
diff --git a/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/material_property/porous_medium/porous_medium/i_porous_medium.md b/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/material_property/porous_medium/porous_medium/i_porous_medium.md
new file mode 100644
index 0000000000000000000000000000000000000000..bb330dc056ceff8400e948dfb52a53fb0f44252a
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/material_property/porous_medium/porous_medium/i_porous_medium.md
@@ -0,0 +1 @@
+A tag for the properties of a porous medium with material ID.
\ No newline at end of file
diff --git a/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/material_property/porous_medium/porous_medium/t_capillary_pressure.md b/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/material_property/porous_medium/porous_medium/t_capillary_pressure.md
new file mode 100644
index 0000000000000000000000000000000000000000..19d854870d0c3537197cc8edb349b03503fcbede
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/material_property/porous_medium/porous_medium/t_capillary_pressure.md
@@ -0,0 +1 @@
+Capillary pressure model.
diff --git a/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/material_property/porous_medium/porous_medium/t_permeability.md b/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/material_property/porous_medium/porous_medium/t_permeability.md
new file mode 100644
index 0000000000000000000000000000000000000000..075c8a48fde6b9a1d7682d9cc08e6bcaa9feeb93
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/material_property/porous_medium/porous_medium/t_permeability.md
@@ -0,0 +1 @@
+A tag for the relative permeability model of a porous medium.
\ No newline at end of file
diff --git a/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/material_property/porous_medium/porous_medium/t_porosity.md b/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/material_property/porous_medium/porous_medium/t_porosity.md
new file mode 100644
index 0000000000000000000000000000000000000000..d90fcba994c498eb7f7c82b6052cf4619ef7735b
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/material_property/porous_medium/porous_medium/t_porosity.md
@@ -0,0 +1 @@
+A tag for the porosity model of a porous medium.
\ No newline at end of file
diff --git a/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/material_property/porous_medium/porous_medium/t_relative_permeability.md b/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/material_property/porous_medium/porous_medium/t_relative_permeability.md
new file mode 100644
index 0000000000000000000000000000000000000000..fd2a51343fb9e4a7a8ce38bdc88c31f400c1b7e7
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/material_property/porous_medium/porous_medium/t_relative_permeability.md
@@ -0,0 +1 @@
+A tag for the relative permeability model.
diff --git a/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/material_property/porous_medium/porous_medium/t_storage.md b/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/material_property/porous_medium/porous_medium/t_storage.md
new file mode 100644
index 0000000000000000000000000000000000000000..756a75cbc74b625cca84e3b6b15dcc29b36026b5
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/material_property/porous_medium/porous_medium/t_storage.md
@@ -0,0 +1 @@
+A tag for the storage model of a porous medium.
\ No newline at end of file
diff --git a/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/process_variables/i_process_variables.md b/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/process_variables/i_process_variables.md
index 576add64edccf5d980a0ff86143dc43a7109a788..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 100644
--- a/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/process_variables/i_process_variables.md
+++ b/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/process_variables/i_process_variables.md
@@ -1 +0,0 @@
-\ogs_missing_documentation
diff --git a/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/process_variables/t_process_variable.md b/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/process_variables/t_process_variable.md
index 576add64edccf5d980a0ff86143dc43a7109a788..b41ba4a237742a7e452d9f8a769a412d3086c487 100644
--- a/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/process_variables/t_process_variable.md
+++ b/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/process_variables/t_process_variable.md
@@ -1 +1 @@
-\ogs_missing_documentation
+A tag for the primary variable in RichardsFlow Process. In this implementation, fluid pressure is used as primary variable.
diff --git a/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/t_intrinsic_permeability.md b/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/t_intrinsic_permeability.md
deleted file mode 100644
index 576add64edccf5d980a0ff86143dc43a7109a788..0000000000000000000000000000000000000000
--- a/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/t_intrinsic_permeability.md
+++ /dev/null
@@ -1 +0,0 @@
-\ogs_missing_documentation
diff --git a/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/t_mass_lumping.md b/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/t_mass_lumping.md
index 576add64edccf5d980a0ff86143dc43a7109a788..b7844df80113ac416e571e335b34b301474f9835 100644
--- a/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/t_mass_lumping.md
+++ b/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/t_mass_lumping.md
@@ -1 +1 @@
-\ogs_missing_documentation
+A tag for enabling mass-lumping scheme
diff --git a/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/t_porosity.md b/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/t_porosity.md
deleted file mode 100644
index 576add64edccf5d980a0ff86143dc43a7109a788..0000000000000000000000000000000000000000
--- a/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/t_porosity.md
+++ /dev/null
@@ -1 +0,0 @@
-\ogs_missing_documentation
diff --git a/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/t_specific_body_force.md b/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/t_specific_body_force.md
index 576add64edccf5d980a0ff86143dc43a7109a788..ed337a7b8e0a63e758207f44be45ffead89e5177 100644
--- a/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/t_specific_body_force.md
+++ b/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/t_specific_body_force.md
@@ -1 +1 @@
-\ogs_missing_documentation
+Specific body forces applied to solid and fluid. It is usually used to apply gravitational forces.
diff --git a/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/t_storage.md b/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/t_storage.md
deleted file mode 100644
index 576add64edccf5d980a0ff86143dc43a7109a788..0000000000000000000000000000000000000000
--- a/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/t_storage.md
+++ /dev/null
@@ -1 +0,0 @@
-\ogs_missing_documentation
diff --git a/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/t_temperature.md b/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/t_temperature.md
new file mode 100644
index 0000000000000000000000000000000000000000..a41a95c52ff9cd3047303ca8f4f5f21e9ea57da7
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/t_temperature.md
@@ -0,0 +1,2 @@
+An input of a reference temperature (K) for RichardsFlow process.
+It is fixed to be constant for the simulation due to an isothermal assumption.
diff --git a/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/t_viscosity.md b/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/t_viscosity.md
deleted file mode 100644
index 576add64edccf5d980a0ff86143dc43a7109a788..0000000000000000000000000000000000000000
--- a/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/t_viscosity.md
+++ /dev/null
@@ -1 +0,0 @@
-\ogs_missing_documentation
diff --git a/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/t_water_density.md b/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/t_water_density.md
deleted file mode 100644
index 576add64edccf5d980a0ff86143dc43a7109a788..0000000000000000000000000000000000000000
--- a/Documentation/ProjectFile/prj/processes/process/RICHARDS_FLOW/t_water_density.md
+++ /dev/null
@@ -1 +0,0 @@
-\ogs_missing_documentation
diff --git a/ProcessLib/RichardsFlow/CreateRichardsFlowMaterialProperties.cpp b/ProcessLib/RichardsFlow/CreateRichardsFlowMaterialProperties.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..da3485782c507565f231676c2075874dc2770ab8
--- /dev/null
+++ b/ProcessLib/RichardsFlow/CreateRichardsFlowMaterialProperties.cpp
@@ -0,0 +1,117 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "CreateRichardsFlowMaterialProperties.h"
+
+#include <logog/include/logog.hpp>
+
+#include "BaseLib/reorderVector.h"
+#include "MaterialLib/Fluid/FluidProperty.h"
+#include "MaterialLib/Fluid/FluidProperties/CreateFluidProperties.h"
+#include "MaterialLib/PorousMedium/Porosity/Porosity.h"
+#include "MaterialLib/PorousMedium/Storage/Storage.h"
+#include "MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CapillaryPressureSaturation.h"
+#include "MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CreateCapillaryPressureModel.h"
+#include "MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/CreateRelativePermeabilityModel.h"
+#include "MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/RelativePermeability.h"
+#include "MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h"
+#include "MeshLib/Mesh.h"
+#include "MeshLib/PropertyVector.h"
+#include "ProcessLib/Parameter/Parameter.h"
+#include "ProcessLib/Parameter/SpatialPosition.h"
+
+#include "RichardsFlowMaterialProperties.h"
+
+namespace ProcessLib
+{
+namespace RichardsFlow
+{
+std::unique_ptr<RichardsFlowMaterialProperties>
+createRichardsFlowMaterialProperties(
+ BaseLib::ConfigTree const& config,
+ boost::optional<MeshLib::PropertyVector<int> const&>
+ material_ids)
+{
+ DBUG("Reading material properties of Richards flow process.");
+
+ //! \ogs_file_param{prj__processes__process__RICHARDS_FLOW__material_property__fluid}
+ auto const& fluid_config = config.getConfigSubtree("fluid");
+ auto fluid_properties =
+ MaterialLib::Fluid::createFluidProperties(fluid_config);
+
+ // Get porous properties
+ std::vector<int> mat_ids;
+ std::vector<int> mat_krel_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;
+ std::vector<
+ std::unique_ptr<MaterialLib::PorousMedium::CapillaryPressureSaturation>>
+ capillary_pressure_models;
+ std::vector<
+ std::unique_ptr<MaterialLib::PorousMedium::RelativePermeability>>
+ relative_permeability_models;
+
+ //! \ogs_file_param{prj__processes__process__RICHARDS_FLOW__material_property__porous_medium}
+ auto const& poro_config = config.getConfigSubtree("porous_medium");
+ //! \ogs_file_param{prj__processes__process__RICHARDS_FLOW__material_property__porous_medium__porous_medium}
+ for (auto const& conf : poro_config.getConfigSubtreeList("porous_medium"))
+ {
+ //! \ogs_file_attr{prj__processes__process__RICHARDS_FLOW__material_property__porous_medium__porous_medium__id}
+ auto const id = conf.getConfigAttributeOptional<int>("id");
+ mat_ids.push_back(*id);
+
+ //! \ogs_file_param{prj__processes__process__RICHARDS_FLOW__material_property__porous_medium__porous_medium__permeability}
+ auto const& permeability_conf = conf.getConfigSubtree("permeability");
+ intrinsic_permeability_models.emplace_back(
+ MaterialLib::PorousMedium::createPermeabilityModel(
+ permeability_conf));
+
+ //! \ogs_file_param{prj__processes__process__RICHARDS_FLOW__material_property__porous_medium__porous_medium__porosity}
+ auto const& porosity_conf = conf.getConfigSubtree("porosity");
+ auto n = MaterialLib::PorousMedium::createPorosityModel(porosity_conf);
+ porosity_models.emplace_back(std::move(n));
+
+ //! \ogs_file_param{prj__processes__process__RICHARDS_FLOW__material_property__porous_medium__porous_medium__storage}
+ auto const& storage_conf = conf.getConfigSubtree("storage");
+ auto beta = MaterialLib::PorousMedium::createStorageModel(storage_conf);
+ storage_models.emplace_back(std::move(beta));
+
+ auto const& capillary_pressure_conf =
+ //! \ogs_file_param{prj__processes__process__RICHARDS_FLOW__material_property__porous_medium__porous_medium__capillary_pressure}
+ conf.getConfigSubtree("capillary_pressure");
+ auto pc = MaterialLib::PorousMedium::createCapillaryPressureModel(
+ capillary_pressure_conf);
+ capillary_pressure_models.emplace_back(std::move(pc));
+
+ auto const& krel_config =
+ //! \ogs_file_param{prj__processes__process__RICHARDS_FLOW__material_property__porous_medium__porous_medium__relative_permeability}
+ conf.getConfigSubtree("relative_permeability");
+ auto krel = MaterialLib::PorousMedium::createRelativePermeabilityModel(
+ krel_config);
+ relative_permeability_models.emplace_back(std::move(krel));
+ }
+
+ BaseLib::reorderVector(intrinsic_permeability_models, mat_ids);
+ BaseLib::reorderVector(porosity_models, mat_ids);
+ BaseLib::reorderVector(storage_models, mat_ids);
+
+ return std::unique_ptr<RichardsFlowMaterialProperties>{
+ new RichardsFlowMaterialProperties{
+ material_ids, std::move(fluid_properties),
+ std::move(intrinsic_permeability_models),
+ std::move(porosity_models), std::move(storage_models),
+ std::move(capillary_pressure_models),
+ std::move(relative_permeability_models)}};
+}
+
+} // end namespace
+} // end namespace
diff --git a/ProcessLib/RichardsFlow/CreateRichardsFlowMaterialProperties.h b/ProcessLib/RichardsFlow/CreateRichardsFlowMaterialProperties.h
new file mode 100644
index 0000000000000000000000000000000000000000..f069f286f291ef3991ff6ce24f78bffffdcfec16
--- /dev/null
+++ b/ProcessLib/RichardsFlow/CreateRichardsFlowMaterialProperties.h
@@ -0,0 +1,29 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include "RichardsFlowMaterialProperties.h"
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace ProcessLib
+{
+namespace RichardsFlow
+{
+std::unique_ptr<RichardsFlowMaterialProperties>
+createRichardsFlowMaterialProperties(
+ BaseLib::ConfigTree const& config,
+ boost::optional<MeshLib::PropertyVector<int> const&>
+ material_ids);
+
+} // end namespace
+} // end namespace
diff --git a/ProcessLib/RichardsFlow/CreateRichardsFlowProcess.cpp b/ProcessLib/RichardsFlow/CreateRichardsFlowProcess.cpp
index 5355698ef95b0ef27d211196af52f686dc000a20..bdcde323b286a701d5cb58b6b8bea000dddebe38 100644
--- a/ProcessLib/RichardsFlow/CreateRichardsFlowProcess.cpp
+++ b/ProcessLib/RichardsFlow/CreateRichardsFlowProcess.cpp
@@ -9,9 +9,11 @@
#include "CreateRichardsFlowProcess.h"
+#include "ProcessLib/Parameter/ConstantParameter.h"
#include "ProcessLib/Utils/ParseSecondaryVariables.h"
#include "ProcessLib/Utils/ProcessUtils.h"
-#include "ProcessLib/Parameter/ConstantParameter.h"
+
+#include "CreateRichardsFlowMaterialProperties.h"
#include "RichardsFlowProcess.h"
#include "RichardsFlowProcessData.h"
@@ -45,75 +47,6 @@ std::unique_ptr<Process> createRichardsFlowProcess(
{//! \ogs_file_param_special{prj__processes__process__RICHARDS_FLOW__process_variables__process_variable}
"process_variable"});
- // Hydraulic conductivity parameter.
- auto& intrinsic_permeability = findParameter<double>(
- config,
- //! \ogs_file_param_special{prj__processes__process__RICHARDS_FLOW__intrinsic_permeability}
- "intrinsic_permeability", parameters, 1);
-
- DBUG("Use \'%s\' as intrinsic permeability parameter.",
- intrinsic_permeability.name.c_str());
-
- // Porosity parameter.
- auto& porosity = findParameter<double>(
- config,
- //! \ogs_file_param_special{prj__processes__process__RICHARDS_FLOW__porosity}
- "porosity", parameters, 1);
-
- DBUG("Use \'%s\' as porosity parameter.", porosity.name.c_str());
-
- // Viscosity parameter.
- auto& viscosity = findParameter<double>(
- config,
- //! \ogs_file_param_special{prj__processes__process__RICHARDS_FLOW__viscosity}
- "viscosity", parameters, 1);
-
- DBUG("Use \'%s\' as porosity parameter.", viscosity.name.c_str());
-
- // storage parameter.
- auto& storage = findParameter<double>(
- config,
- //! \ogs_file_param_special{prj__processes__process__RICHARDS_FLOW__storage}
- "storage", parameters, 1);
-
- DBUG("Use \'%s\' as storage parameter.", storage.name.c_str());
-
- // water_density parameter.
- auto& water_density = findParameter<double>(
- config,
- //! \ogs_file_param_special{prj__processes__process__RICHARDS_FLOW__water_density}
- "water_density", parameters, 1);
-
- DBUG("Use \'%s\' as storage parameter.", water_density.name.c_str());
-
- // Specific body force parameter.
- auto& specific_body_force = findParameter<double>(
- config,
- //! \ogs_file_param_special{prj__processes__process__RICHARDS_FLOW__specific_body_force}
- "specific_body_force", parameters, mesh.getDimension());
- DBUG("Use \'%s\' as specific body force parameter.",
- specific_body_force.name.c_str());
-
- // Assume constant parameter, then check the norm at arbitrary
- // SpatialPosition and time.
- assert(dynamic_cast<ConstantParameter<double>*>(&specific_body_force));
- bool const has_gravity =
- MathLib::toVector(specific_body_force(0, SpatialPosition{})).norm() > 0;
-
- // has mass lumping
- //! \ogs_file_param{prj__processes__process__RICHARDS_FLOW__mass_lumping}
- auto mass_lump = config.getConfigParameter<bool>("mass_lumping");
-
- RichardsFlowProcessData process_data{intrinsic_permeability,
- porosity,
- viscosity,
- storage,
- water_density,
- specific_body_force,
- has_gravity,
- mass_lump,
- *curves.at("curve_PC_S"),
- *curves.at("curve_S_Krel")};
SecondaryVariableCollection secondary_variables;
NumLib::NamedFunctionCaller named_function_caller(
@@ -122,10 +55,50 @@ std::unique_ptr<Process> createRichardsFlowProcess(
ProcessLib::parseSecondaryVariables(config, secondary_variables,
named_function_caller);
+ // Specific body force
+ std::vector<double> const b =
+ //! \ogs_file_param{prj__processes__process__RICHARDS_FLOW__specific_body_force}
+ config.getConfigParameter<std::vector<double>>("specific_body_force");
+ assert(b.size() > 0 && b.size() < 4);
+ Eigen::VectorXd specific_body_force(b.size());
+ bool const has_gravity = MathLib::toVector(b).norm() > 0;
+ if (has_gravity)
+ std::copy_n(b.data(), b.size(), specific_body_force.data());
+
+ // has mass lumping
+ //! \ogs_file_param{prj__processes__process__RICHARDS_FLOW__mass_lumping}
+ auto mass_lumping = config.getConfigParameter<bool>("mass_lumping");
+
+ auto& temperature = findParameter<double>(
+ config,
+ //! \ogs_file_param_special{prj__processes__process__RICHARDS_FLOW__temperature}
+ "temperature", parameters, 1);
+
+ //! \ogs_file_param{prj__processes__process__RICHARDS_FLOW__material_property}
+ auto const& mat_config = config.getConfigSubtree("material_property");
+
+ boost::optional<MeshLib::PropertyVector<int> const&> material_ids;
+ if (mesh.getProperties().existsPropertyVector<int>("MaterialIDs"))
+ {
+ INFO("The Richards flow is in heterogeneous porous media.");
+ material_ids =
+ *mesh.getProperties().getPropertyVector<int>("MaterialIDs");
+ }
+ else
+ {
+ INFO("The Richards flow is in homogeneous porous media.");
+ }
+ std::unique_ptr<RichardsFlowMaterialProperties> material =
+ createRichardsFlowMaterialProperties(mat_config, material_ids);
+ RichardsFlowProcessData process_data{std::move(material),
+ specific_body_force, has_gravity,
+ mass_lumping, temperature};
+
return std::make_unique<RichardsFlowProcess>(
mesh, std::move(jacobian_assembler), parameters, integration_order,
std::move(process_variables), std::move(process_data),
- std::move(secondary_variables), std::move(named_function_caller));
+ std::move(secondary_variables), std::move(named_function_caller),
+ mat_config, curves);
}
} // namespace RichardsFlow
diff --git a/ProcessLib/RichardsFlow/RichardsFlowFEM.h b/ProcessLib/RichardsFlow/RichardsFlowFEM.h
index caf5a78f60533223ae539f9068584dd42d8a4c7a..319bbacd977d968385ed72d5bb0d855934450262 100644
--- a/ProcessLib/RichardsFlow/RichardsFlowFEM.h
+++ b/ProcessLib/RichardsFlow/RichardsFlowFEM.h
@@ -28,6 +28,26 @@ namespace ProcessLib
{
namespace RichardsFlow
{
+template <typename NodalRowVectorType, typename GlobalDimNodalMatrixType,
+ typename NodalMatrixType>
+struct IntegrationPointData final
+{
+ IntegrationPointData(NodalRowVectorType const& N_,
+ GlobalDimNodalMatrixType const& dNdx_,
+ double const& integration_weight_,
+ NodalMatrixType const mass_operator_)
+ : N(N_),
+ dNdx(dNdx_),
+ integration_weight(integration_weight_),
+ mass_operator(mass_operator_)
+ {
+ }
+
+ NodalRowVectorType const N;
+ GlobalDimNodalMatrixType const dNdx;
+ double const integration_weight;
+ NodalMatrixType const mass_operator;
+};
const unsigned NUM_NODAL_DOF = 1;
class RichardsFlowLocalAssemblerInterface
@@ -58,28 +78,54 @@ class LocalAssemblerData : public RichardsFlowLocalAssemblerInterface
using LocalAssemblerTraits = ProcessLib::LocalAssemblerTraits<
ShapeMatricesType, ShapeFunction::NPOINTS, NUM_NODAL_DOF, GlobalDim>;
- using NodalMatrixType = typename LocalAssemblerTraits::LocalMatrix;
- using NodalVectorType = typename LocalAssemblerTraits::LocalVector;
+ using NodalRowVectorType = typename ShapeMatricesType::NodalRowVectorType;
+ using GlobalDimNodalMatrixType =
+ typename ShapeMatricesType::GlobalDimNodalMatrixType;
+ using NodalMatrixType = typename ShapeMatricesType::NodalMatrixType;
+ using NodalVectorType = typename ShapeMatricesType::NodalVectorType;
+ using GlobalDimMatrixType = typename ShapeMatricesType::GlobalDimMatrixType;
using GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType;
public:
LocalAssemblerData(MeshLib::Element const& element,
- std::size_t const /*local_matrix_size*/,
+ std::size_t const local_matrix_size,
bool is_axially_symmetric,
unsigned const integration_order,
RichardsFlowProcessData const& process_data)
: _element(element),
_process_data(process_data),
_integration_method(integration_order),
- _shape_matrices(initShapeMatrices<ShapeFunction, ShapeMatricesType,
- IntegrationMethod, GlobalDim>(
- element, is_axially_symmetric, _integration_method)),
_saturation(
std::vector<double>(_integration_method.getNumberOfPoints())),
_darcy_velocities(
GlobalDim,
std::vector<double>(_integration_method.getNumberOfPoints()))
{
+ // This assertion is valid only if all nodal d.o.f. use the same shape
+ // matrices.
+ assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF);
+ (void)local_matrix_size;
+
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+ _ip_data.reserve(n_integration_points);
+
+ auto const shape_matrices =
+ initShapeMatrices<ShapeFunction, ShapeMatricesType,
+ IntegrationMethod, GlobalDim>(
+ element, is_axially_symmetric, _integration_method);
+
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ auto const& sm = shape_matrices[ip];
+ const double integration_factor = sm.integralMeasure * sm.detJ;
+ _ip_data.emplace_back(
+ sm.N, sm.dNdx,
+ _integration_method.getWeightedPoint(ip).getWeight() *
+ integration_factor,
+ sm.N.transpose() * sm.N * integration_factor *
+ _integration_method.getWeightedPoint(ip).getWeight());
+ }
}
void assemble(double const t, std::vector<double> const& local_x,
@@ -92,11 +138,6 @@ public:
// matrices.
assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF);
- MathLib::PiecewiseLinearInterpolation const& interpolated_Pc =
- _process_data.interpolated_Pc;
- MathLib::PiecewiseLinearInterpolation const& interpolated_Kr =
- _process_data.interpolated_Kr;
-
auto local_M = MathLib::createZeroedMatrix<NodalMatrixType>(
local_M_data, local_matrix_size, local_matrix_size);
auto local_K = MathLib::createZeroedMatrix<NodalMatrixType>(
@@ -108,59 +149,72 @@ public:
_integration_method.getNumberOfPoints();
SpatialPosition pos;
pos.setElementID(_element.getID());
+ const int material_id =
+ _process_data.material->getMaterialID(_element.getID());
+ const Eigen::MatrixXd& perm = _process_data.material->getPermeability(
+ material_id, t, pos, _element.getDimension());
+ assert(perm.rows() == _element.getDimension() || perm.rows() == 1);
+ GlobalDimMatrixType permeability = GlobalDimMatrixType::Zero(
+ _element.getDimension(), _element.getDimension());
+ if (perm.rows() == _element.getDimension())
+ permeability = perm;
+ else if (perm.rows() == 1)
+ permeability.diagonal().setConstant(perm(0, 0));
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
pos.setIntegrationPoint(ip);
- auto const& sm = _shape_matrices[ip];
- auto const& wp = _integration_method.getWeightedPoint(ip);
-
- double P_int_pt = 0.0;
- NumLib::shapeFunctionInterpolate(local_x, sm.N, P_int_pt);
-
- auto const K = _process_data.intrinsic_permeability(t, pos)[0];
- auto const poro = _process_data.porosity(t, pos)[0];
- auto const mu = _process_data.viscosity(t, pos)[0];
- auto const storage = _process_data.storage(t, pos)[0];
- double const Pc = -P_int_pt;
-
- double dSwdPc = interpolated_Pc.getDerivative(Pc);
- if (Pc > interpolated_Pc.getSupportMax())
- dSwdPc = interpolated_Pc.getDerivative(
- interpolated_Pc.getSupportMax());
- else if (Pc < interpolated_Pc.getSupportMin())
- dSwdPc = interpolated_Pc.getDerivative(
- interpolated_Pc.getSupportMin());
-
- double const Sw = interpolated_Pc.getValue(Pc);
+ double p_int_pt = 0.0;
+ NumLib::shapeFunctionInterpolate(local_x, _ip_data[ip].N, p_int_pt);
+ const double& temperature = _process_data.temperature(t, pos)[0];
+ auto const porosity = _process_data.material->getPorosity(
+ material_id, t, pos, p_int_pt, temperature, 0);
+
+ double const pc_int_pt = -p_int_pt;
+
+ double const Sw =
+ (pc_int_pt > 0)
+ ? _process_data.material->getSaturation(
+ material_id, t, pos, p_int_pt, temperature, pc_int_pt)
+ : 1.0;
_saturation[ip] = Sw;
- double const k_rel = interpolated_Kr.getValue(Sw);
+ double const dSw_dpc =
+ (pc_int_pt > 0)
+ ? _process_data.material->getSaturationDerivative(
+ material_id, t, pos, p_int_pt, temperature, Sw)
+ : 0.;
+
// \TODO Extend to pressure dependent density.
double const drhow_dp(0.0);
-
+ auto const storage = _process_data.material->getStorage(
+ material_id, t, pos, p_int_pt, temperature, 0);
double const mass_mat_coeff =
- storage * Sw + poro * Sw * drhow_dp - poro * dSwdPc;
- double const K_mat_coeff = K * k_rel / mu;
+ storage * Sw + porosity * Sw * drhow_dp - porosity * dSw_dpc;
- local_K.noalias() += sm.dNdx.transpose() * K_mat_coeff * sm.dNdx *
- sm.detJ * sm.integralMeasure * wp.getWeight();
+ local_M.noalias() += mass_mat_coeff * _ip_data[ip].mass_operator;
+
+ double const k_rel =
+ _process_data.material->getRelativePermeability(
+ t, pos, p_int_pt, temperature, Sw);
+ auto const mu = _process_data.material->getFluidViscosity(
+ p_int_pt, temperature);
+ local_K.noalias() += _ip_data[ip].dNdx.transpose() * permeability *
+ _ip_data[ip].dNdx *
+ _ip_data[ip].integration_weight * (k_rel / mu);
- local_M.noalias() += sm.N.transpose() * mass_mat_coeff * sm.N *
- sm.detJ * sm.integralMeasure * wp.getWeight();
if (_process_data.has_gravity)
{
- auto const rho_w = _process_data.water_density(t, pos)[0];
- auto const body_force =
- _process_data.specific_body_force(t, pos);
+ auto const rho_w = _process_data.material->getFluidDensity(
+ p_int_pt, temperature);
+ auto const& body_force = _process_data.specific_body_force;
assert(body_force.size() == GlobalDim);
- auto const b = MathLib::toVector<GlobalDimVectorType>(
- body_force, GlobalDim);
- local_b.noalias() += sm.dNdx.transpose() * K_mat_coeff * rho_w *
- b * sm.detJ * sm.integralMeasure *
- wp.getWeight();
+ NodalVectorType gravity_operator =
+ _ip_data[ip].dNdx.transpose() * permeability * body_force *
+ _ip_data[ip].integration_weight;
+ local_b.noalias() += (k_rel / mu) * rho_w * gravity_operator;
}
- } // end of GP
+ }
if (_process_data.has_mass_lumping)
{
for (int idx_ml = 0; idx_ml < local_M.cols(); idx_ml++)
@@ -177,13 +231,18 @@ public:
{
SpatialPosition pos;
pos.setElementID(_element.getID());
-
- auto const K = _process_data.intrinsic_permeability(t, pos)[0];
- MathLib::PiecewiseLinearInterpolation const& interpolated_Pc =
- _process_data.interpolated_Pc;
- MathLib::PiecewiseLinearInterpolation const& interpolated_Kr =
- _process_data.interpolated_Kr;
- auto const mu = _process_data.viscosity(t, pos)[0];
+ const int material_id =
+ _process_data.material->getMaterialID(_element.getID());
+
+ const Eigen::MatrixXd& perm = _process_data.material->getPermeability(
+ material_id, t, pos, _element.getDimension());
+ assert(perm.rows() == _element.getDimension() || perm.rows() == 1);
+ GlobalDimMatrixType permeability = GlobalDimMatrixType::Zero(
+ _element.getDimension(), _element.getDimension());
+ if (perm.rows() == _element.getDimension())
+ permeability = perm;
+ else if (perm.rows() == 1)
+ permeability.diagonal().setConstant(perm(0, 0));
unsigned const n_integration_points =
_integration_method.getNumberOfPoints();
@@ -193,27 +252,29 @@ public:
for (unsigned ip = 0; ip < n_integration_points; ++ip)
{
- auto const& sm = _shape_matrices[ip];
-
- double P_int_pt = 0.0;
- NumLib::shapeFunctionInterpolate(local_x, sm.N, P_int_pt);
- double const Pc = -P_int_pt;
- double const Sw = interpolated_Pc.getValue(Pc);
- double const k_rel = interpolated_Kr.getValue(Sw);
- double const K_mat_coeff = K * k_rel / mu;
-
+ double p_int_pt = 0.0;
+ NumLib::shapeFunctionInterpolate(local_x, _ip_data[ip].N, p_int_pt);
+ double const pc_int_pt = -p_int_pt;
+ const double& temperature = _process_data.temperature(t, pos)[0];
+ double const Sw = _process_data.material->getSaturation(
+ material_id, t, pos, p_int_pt, temperature, pc_int_pt);
+ double const k_rel =
+ _process_data.material->getRelativePermeability(
+ t, pos, p_int_pt, temperature, Sw);
+ auto const mu = _process_data.material->getFluidViscosity(
+ p_int_pt, temperature);
+ auto const K_mat_coeff = permeability * (k_rel / mu);
GlobalDimVectorType velocity =
- -K_mat_coeff * sm.dNdx * p_nodal_values;
+ -K_mat_coeff * _ip_data[ip].dNdx * p_nodal_values;
if (_process_data.has_gravity)
{
- auto const rho_w = _process_data.water_density(t, pos)[0];
- auto const body_force =
- _process_data.specific_body_force(t, pos);
+ auto const rho_w = _process_data.material->getFluidDensity(
+ p_int_pt, temperature);
+ auto const& body_force = _process_data.specific_body_force;
assert(body_force.size() == GlobalDim);
- auto const b = MathLib::toVector<GlobalDimVectorType>(
- body_force, GlobalDim);
- // here it is assumed that the vector b is directed 'downwards'
- velocity += K_mat_coeff * rho_w * b;
+ // here it is assumed that the vector body_force is directed
+ // 'downwards'
+ velocity += K_mat_coeff * rho_w * body_force;
}
for (unsigned d = 0; d < GlobalDim; ++d)
@@ -226,7 +287,7 @@ public:
Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
const unsigned integration_point) const override
{
- auto const& N = _shape_matrices[integration_point].N;
+ auto const& N = _ip_data[integration_point].N;
// assumes N is stored contiguously in memory
return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
@@ -267,7 +328,12 @@ private:
IntegrationMethod const _integration_method;
std::vector<ShapeMatrices, Eigen::aligned_allocator<ShapeMatrices>>
_shape_matrices;
-
+ std::vector<
+ IntegrationPointData<NodalRowVectorType, GlobalDimNodalMatrixType,
+ NodalMatrixType>,
+ Eigen::aligned_allocator<IntegrationPointData<
+ NodalRowVectorType, GlobalDimNodalMatrixType, NodalMatrixType>>>
+ _ip_data;
std::vector<double> _saturation;
std::vector<std::vector<double>> _darcy_velocities;
};
diff --git a/ProcessLib/RichardsFlow/RichardsFlowMaterialProperties.cpp b/ProcessLib/RichardsFlow/RichardsFlowMaterialProperties.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..3c5f64280b7077a777ff0e79ce5d5f4f8fb7d809
--- /dev/null
+++ b/ProcessLib/RichardsFlow/RichardsFlowMaterialProperties.cpp
@@ -0,0 +1,136 @@
+/**
+* \copyright
+* Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+* Distributed under a Modified BSD License.
+* See accompanying file LICENSE.txt or
+* http://www.opengeosys.org/project/license
+*
+*/
+
+#include "RichardsFlowMaterialProperties.h"
+
+#include <logog/include/logog.hpp>
+
+#include "BaseLib/reorderVector.h"
+#include "MaterialLib/Fluid/FluidProperty.h"
+#include "MaterialLib/PorousMedium/Porosity/Porosity.h"
+#include "MaterialLib/PorousMedium/Storage/Storage.h"
+#include "MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CapillaryPressureSaturation.h"
+#include "MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CreateCapillaryPressureModel.h"
+#include "MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/CreateRelativePermeabilityModel.h"
+#include "MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/RelativePermeability.h"
+#include "MeshLib/Mesh.h"
+#include "MeshLib/PropertyVector.h"
+#include "ProcessLib/Parameter/Parameter.h"
+#include "ProcessLib/Parameter/SpatialPosition.h"
+
+namespace ProcessLib
+{
+namespace RichardsFlow
+{
+RichardsFlowMaterialProperties::RichardsFlowMaterialProperties(
+ boost::optional<MeshLib::PropertyVector<int> const&> const material_ids,
+ std::unique_ptr<MaterialLib::Fluid::FluidProperties>&& fluid_properties,
+ 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,
+ std::vector<std::unique_ptr<
+ MaterialLib::PorousMedium::CapillaryPressureSaturation>>&&
+ capillary_pressure_models,
+ std::vector<
+ std::unique_ptr<MaterialLib::PorousMedium::RelativePermeability>>&&
+ relative_permeability_models)
+ : _material_ids(material_ids),
+ _fluid_properties(std::move(fluid_properties)),
+ _intrinsic_permeability_models(intrinsic_permeability_models),
+ _porosity_models(std::move(porosity_models)),
+ _storage_models(std::move(storage_models)),
+ _capillary_pressure_models(std::move(capillary_pressure_models)),
+ _relative_permeability_models(std::move(relative_permeability_models))
+{
+ DBUG("Create material properties for Richards flow.");
+}
+
+int RichardsFlowMaterialProperties::getMaterialID(const std::size_t element_id)
+{
+ if (!_material_ids)
+ {
+ return 0;
+ }
+
+ assert(element_id < _material_ids->size());
+ return (*_material_ids)[element_id];
+}
+
+double RichardsFlowMaterialProperties::getFluidDensity(const double p,
+ const double T) const
+{
+ ArrayType vars;
+ vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::T)] = T;
+ vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::p)] = p;
+ return _fluid_properties->getValue(
+ MaterialLib::Fluid::FluidPropertyType::Density, vars);
+}
+
+double RichardsFlowMaterialProperties::getFluidViscosity(const double p,
+ const double T) const
+{
+ ArrayType vars;
+ vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::T)] = T;
+ vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::p)] = p;
+ return _fluid_properties->getValue(
+ MaterialLib::Fluid::FluidPropertyType::Viscosity, vars);
+}
+
+Eigen::MatrixXd const& RichardsFlowMaterialProperties::getPermeability(
+ const int material_id, const double /*t*/,
+ const ProcessLib::SpatialPosition& /*pos*/, const int /*dim*/) const
+{
+ return _intrinsic_permeability_models[material_id];
+}
+
+double RichardsFlowMaterialProperties::getPorosity(
+ const int material_id, const double /*t*/,
+ const ProcessLib::SpatialPosition& /*pos*/, const double /*p*/,
+ const double T, const double porosity_variable) const
+{
+ return _porosity_models[material_id]->getValue(porosity_variable, T);
+}
+
+double RichardsFlowMaterialProperties::getStorage(
+ const int material_id, const double /*t*/,
+ const ProcessLib::SpatialPosition& /*pos*/, const double /*p*/,
+ const double /*T*/, const double storage_variable) const
+{
+ // \todo getValue() can be extended for non
+ // constant storage model
+ return _storage_models[material_id]->getValue(storage_variable);
+}
+
+double RichardsFlowMaterialProperties::getRelativePermeability(
+ const double /*t*/, const ProcessLib::SpatialPosition& /*pos*/,
+ const double /*p*/, const double /*T*/, const double saturation) const
+{
+ return _relative_permeability_models[0]->getValue(saturation);
+}
+
+double RichardsFlowMaterialProperties::getSaturation(
+ const int material_id, const double /*t*/,
+ const ProcessLib::SpatialPosition& /*pos*/, const double /*p*/,
+ const double /*T*/, const double pc) const
+{
+ return _capillary_pressure_models[material_id]->getSaturation(pc);
+}
+double RichardsFlowMaterialProperties::getSaturationDerivative(
+ const int material_id, const double /*t*/,
+ const ProcessLib::SpatialPosition& /*pos*/, const double /*p*/,
+ const double /*T*/, const double saturation) const
+{
+ const double dpcdsw =
+ _capillary_pressure_models[material_id]->getdPcdS(saturation);
+ return 1 / dpcdsw;
+}
+} // end of namespace
+} // end of namespace
diff --git a/ProcessLib/RichardsFlow/RichardsFlowMaterialProperties.h b/ProcessLib/RichardsFlow/RichardsFlowMaterialProperties.h
new file mode 100644
index 0000000000000000000000000000000000000000..cf6db124c526c25b9dc8a6d93c1a9092afaaa1ca
--- /dev/null
+++ b/ProcessLib/RichardsFlow/RichardsFlowMaterialProperties.h
@@ -0,0 +1,121 @@
+/**
+* \copyright
+* Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+* Distributed under a Modified BSD License.
+* See accompanying file LICENSE.txt or
+* http://www.opengeosys.org/project/license
+*
+*/
+
+#pragma once
+
+#include <memory>
+#include <vector>
+
+#include "MaterialLib/Fluid/FluidPropertyHeaders.h"
+#include "MaterialLib/Fluid/FluidProperties/FluidProperties.h"
+#include "MaterialLib/PhysicalConstant.h"
+#include "MaterialLib/PorousMedium/Porosity/Porosity.h"
+#include "MaterialLib/PorousMedium/PorousPropertyHeaders.h"
+#include "MaterialLib/PorousMedium/Storage/Storage.h"
+#include "MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CapillaryPressureSaturation.h"
+#include "MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CreateCapillaryPressureModel.h"
+#include "MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/CreateRelativePermeabilityModel.h"
+#include "MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/RelativePermeability.h"
+
+namespace MeshLib
+{
+template <typename PROP_VAL_TYPE>
+class PropertyVector;
+}
+
+namespace ProcessLib
+{
+class SpatialPosition;
+namespace RichardsFlow
+{
+/** This class has a collection of material properties for Richards flow model
+* and it makes description of the properties of unsaturated porous media
+* i.e. the fluid density and viscosity models
+* the relative permeability models,
+* the capillary pressure-saturation relationships.
+* It generally provides the computation of the PDE coefficients for Richards
+* flow.
+*/
+
+class RichardsFlowMaterialProperties final
+{
+public:
+ using ArrayType = MaterialLib::Fluid::FluidProperty::ArrayType;
+
+ RichardsFlowMaterialProperties(
+ boost::optional<MeshLib::PropertyVector<int> const&> const material_ids,
+ std::unique_ptr<MaterialLib::Fluid::FluidProperties>&& fluid_properties,
+ 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,
+ std::vector<std::unique_ptr<
+ MaterialLib::PorousMedium::CapillaryPressureSaturation>>&&
+ capillary_pressure_models,
+ std::vector<
+ std::unique_ptr<MaterialLib::PorousMedium::RelativePermeability>>&&
+ relative_permeability_models);
+
+ int getMaterialID(const std::size_t element_id);
+
+ Eigen::MatrixXd const& getPermeability(
+ const int material_id,
+ const double t,
+ const ProcessLib::SpatialPosition& pos,
+ const int dim) const;
+
+ double getPorosity(const int material_id, const double t,
+ const ProcessLib::SpatialPosition& pos, const double p,
+ const double T, const double porosity_variable) const;
+
+ double getStorage(
+ const int material_id, const double t,
+ const ProcessLib::SpatialPosition& pos, const double p,
+ const double T, const double storage_variable) const;
+
+ double getRelativePermeability(const double t,
+ const ProcessLib::SpatialPosition& pos,
+ const double p, const double T,
+ const double saturation) const;
+
+ double getSaturation(const int material_id, const double t,
+ const ProcessLib::SpatialPosition& pos, const double p,
+ const double T, const double pc) const;
+ double getSaturationDerivative(const int material_id, const double t,
+ const ProcessLib::SpatialPosition& pos,
+ const double p, const double T,
+ const double saturation) const;
+ double getFluidDensity(const double p, const double T) const;
+ double getFluidViscosity(const double p, const double T) const;
+
+private:
+ /**
+ * Material IDs must be given as mesh element properties.
+ */
+ boost::optional<MeshLib::PropertyVector<int> const&> const _material_ids;
+
+ const std::unique_ptr<MaterialLib::Fluid::FluidProperties>
+ _fluid_properties;
+
+ std::vector<Eigen::MatrixXd> const _intrinsic_permeability_models;
+ std::vector<std::unique_ptr<MaterialLib::PorousMedium::Porosity>> const
+ _porosity_models;
+ std::vector<std::unique_ptr<MaterialLib::PorousMedium::Storage>> const
+ _storage_models;
+ std::vector<std::unique_ptr<
+ MaterialLib::PorousMedium::CapillaryPressureSaturation>> const
+ _capillary_pressure_models;
+ std::vector<
+ std::unique_ptr<MaterialLib::PorousMedium::RelativePermeability>> const
+ _relative_permeability_models;
+};
+
+} // end of namespace
+} // end of namespace
diff --git a/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp b/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp
index 2b3d2a741792628ca8a926fc3fbeac4ad816efcc..c87d3dd6cdf843931c0a104d42f4ce22b5e3b5e8 100644
--- a/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp
+++ b/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp
@@ -25,7 +25,11 @@ RichardsFlowProcess::RichardsFlowProcess(
std::vector<std::reference_wrapper<ProcessVariable>>&& process_variables,
RichardsFlowProcessData&& process_data,
SecondaryVariableCollection&& secondary_variables,
- NumLib::NamedFunctionCaller&& named_function_caller)
+ NumLib::NamedFunctionCaller&& named_function_caller,
+ BaseLib::ConfigTree const& /*config*/,
+ std::map<std::string,
+ std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&
+ /*curves*/)
: Process(mesh, std::move(jacobian_assembler), parameters,
integration_order, std::move(process_variables),
std::move(secondary_variables), std::move(named_function_caller)),
@@ -74,13 +78,13 @@ void RichardsFlowProcess::initializeConcreteProcess(
}
}
-void RichardsFlowProcess::assembleConcreteProcess(const double t,
- GlobalVector const& x,
- GlobalMatrix& M,
- GlobalMatrix& K,
- GlobalVector& b,
- StaggeredCouplingTerm const&
- coupling_term)
+void RichardsFlowProcess::assembleConcreteProcess(
+ const double t,
+ GlobalVector const& x,
+ GlobalMatrix& M,
+ GlobalMatrix& K,
+ GlobalVector& b,
+ StaggeredCouplingTerm const& coupling_term)
{
DBUG("Assemble RichardsFlowProcess.");
diff --git a/ProcessLib/RichardsFlow/RichardsFlowProcess.h b/ProcessLib/RichardsFlow/RichardsFlowProcess.h
index 6832699082d830df9fde213a7be616a3689654a1..f1b44fa2272d61ddcdb5e666c538dc29e4349cba 100644
--- a/ProcessLib/RichardsFlow/RichardsFlowProcess.h
+++ b/ProcessLib/RichardsFlow/RichardsFlowProcess.h
@@ -11,6 +11,7 @@
#include "NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h"
#include "ProcessLib/Process.h"
+
#include "RichardsFlowFEM.h"
#include "RichardsFlowProcessData.h"
@@ -33,7 +34,11 @@ public:
process_variables,
RichardsFlowProcessData&& process_data,
SecondaryVariableCollection&& secondary_variables,
- NumLib::NamedFunctionCaller&& named_function_caller);
+ NumLib::NamedFunctionCaller&& named_function_caller,
+ BaseLib::ConfigTree const& config,
+ std::map<std::string,
+ std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&
+ curves);
//! \name ODESystem interface
//! @{
@@ -51,11 +56,9 @@ private:
MeshLib::Mesh const& mesh,
unsigned const integration_order) override;
- void assembleConcreteProcess(const double t, GlobalVector const& x,
- GlobalMatrix& M, GlobalMatrix& K,
- GlobalVector& b,
- StaggeredCouplingTerm const& coupling_term
- ) override;
+ void assembleConcreteProcess(
+ const double t, GlobalVector const& x, GlobalMatrix& M, GlobalMatrix& K,
+ GlobalVector& b, StaggeredCouplingTerm const& coupling_term) override;
void assembleWithJacobianConcreteProcess(
const double t, GlobalVector const& x, GlobalVector const& xdot,
diff --git a/ProcessLib/RichardsFlow/RichardsFlowProcessData.h b/ProcessLib/RichardsFlow/RichardsFlowProcessData.h
index 4d7ea9ebb51750fc428e8dd9bacd576b557c3f1d..a2c0a9419c54e0392ed9bad6553c3bdc47b8db5e 100644
--- a/ProcessLib/RichardsFlow/RichardsFlowProcessData.h
+++ b/ProcessLib/RichardsFlow/RichardsFlowProcessData.h
@@ -8,11 +8,7 @@
*/
#pragma once
-
-namespace MeshLib
-{
-class Element;
-}
+#include "RichardsFlowMaterialProperties.h"
namespace ProcessLib
{
@@ -24,40 +20,25 @@ namespace RichardsFlow
struct RichardsFlowProcessData
{
RichardsFlowProcessData(
- Parameter<double> const& intrinsic_permeability_,
- Parameter<double> const& porosity_,
- Parameter<double> const& viscosity_,
- Parameter<double> const& storage_,
- Parameter<double> const& water_density_,
- Parameter<double> const& specific_body_force_,
+ std::unique_ptr<RichardsFlowMaterialProperties>&& material_,
+ Eigen::VectorXd const specific_body_force_,
bool const has_gravity_,
bool const has_mass_lumping_,
- MathLib::PiecewiseLinearInterpolation const& interpolated_Pc_,
- MathLib::PiecewiseLinearInterpolation const& interpolated_Kr_)
- : intrinsic_permeability(intrinsic_permeability_),
- porosity(porosity_),
- viscosity(viscosity_),
- storage(storage_),
- water_density(water_density_),
+ Parameter<double> const& temperature_)
+ : material(std::move(material_)),
specific_body_force(specific_body_force_),
has_gravity(has_gravity_),
has_mass_lumping(has_mass_lumping_),
- interpolated_Pc(interpolated_Pc_),
- interpolated_Kr(interpolated_Kr_)
+ temperature(temperature_)
{
}
RichardsFlowProcessData(RichardsFlowProcessData&& other)
- : intrinsic_permeability(other.intrinsic_permeability),
- porosity(other.porosity),
- viscosity(other.viscosity),
- storage(other.storage),
- water_density(other.water_density),
+ : material(std::move(other.material)),
specific_body_force(other.specific_body_force),
has_gravity(other.has_gravity),
has_mass_lumping(other.has_mass_lumping),
- interpolated_Pc(other.interpolated_Pc),
- interpolated_Kr(other.interpolated_Kr)
+ temperature(other.temperature)
{
}
@@ -70,16 +51,11 @@ struct RichardsFlowProcessData
//! Assignments are not needed.
void operator=(RichardsFlowProcessData&&) = delete;
- Parameter<double> const& intrinsic_permeability;
- Parameter<double> const& porosity;
- Parameter<double> const& viscosity;
- Parameter<double> const& storage;
- Parameter<double> const& water_density;
- Parameter<double> const& specific_body_force;
+ std::unique_ptr<RichardsFlowMaterialProperties> material;
+ Eigen::VectorXd const specific_body_force;
bool const has_gravity;
bool const has_mass_lumping;
- MathLib::PiecewiseLinearInterpolation const& interpolated_Pc;
- MathLib::PiecewiseLinearInterpolation const& interpolated_Kr;
+ Parameter<double> const& temperature;
};
} // namespace RichardsFlow
diff --git a/ProcessLib/RichardsFlow/Tests.cmake b/ProcessLib/RichardsFlow/Tests.cmake
index 46511bb561d41e3afbb08c151ea7fdc2eecc1008..35b569a8d578bc4b292ee4d98cadba46c1c6acb0 100644
--- a/ProcessLib/RichardsFlow/Tests.cmake
+++ b/ProcessLib/RichardsFlow/Tests.cmake
@@ -1,11 +1,33 @@
AddTest(
- NAME LARGE_2D_RichardsFlow_h_us_quad
+ NAME 2D_RichardsFlow_h_us_quad_ogs5
PATH Parabolic/Richards
EXECUTABLE ogs
- EXECUTABLE_ARGS RichardsFlow_2d.prj
+ EXECUTABLE_ARGS RichardsFlow_2d_compare_ogs5.prj
TESTER vtkdiff
ABSTOL 1e-1 RELTOL 1e-1
DIFF_DATA
h_us_quad_1000.vtu richards_pcs_0_ts_100_t_100.000000.vtu PRESSURE1 pressure
REQUIREMENTS NOT OGS_USE_MPI
)
+AddTest(
+ NAME 2D_RichardsFlow_h_us_quad_small
+ PATH Parabolic/Richards
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS RichardsFlow_2d_small.prj
+ TESTER vtkdiff
+ ABSTOL 1e-8 RELTOL 1e-8
+ DIFF_DATA
+ ref_t_1600.000000.vtu richards_pcs_0_ts_1100_t_1600.000000.vtu pressure pressure
+ REQUIREMENTS NOT OGS_USE_MPI
+)
+AddTest(
+ NAME LARGE_2D_RichardsFlow_h_us_quad
+ PATH Parabolic/Richards
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS RichardsFlow_2d_large.prj
+ TESTER vtkdiff
+ ABSTOL 1e-8 RELTOL 1e-8
+ DIFF_DATA
+ ref_t_20000.000000.vtu richards_pcs_0_ts_18200_t_20000.000000.vtu pressure pressure
+ REQUIREMENTS NOT OGS_USE_MPI
+)
diff --git a/Tests/Data b/Tests/Data
index a456c5b48e931239647852196feef1ca6fb0dece..1b30c432bbe45271710b54e83f84ad4f5d75a033 160000
--- a/Tests/Data
+++ b/Tests/Data
@@ -1 +1 @@
-Subproject commit a456c5b48e931239647852196feef1ca6fb0dece
+Subproject commit 1b30c432bbe45271710b54e83f84ad4f5d75a033