From 0ca170b8f3ea2372f03e4ff2a1bfad3e946ccc65 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=B6rg=20Buchwald?= <joerg.buchwald@ufz.de>
Date: Tue, 15 Jun 2021 18:24:41 +0200
Subject: [PATCH] add ThermoRichardsFlow process
---
Applications/ApplicationsLib/ProjectData.cpp | 15 +-
ProcessLib/ThermoRichardsFlow/CMakeLists.txt | 8 +
.../CreateThermoRichardsFlowProcess.cpp | 170 +++
.../CreateThermoRichardsFlowProcess.h | 54 +
.../ThermoRichardsFlow/IntegrationPointData.h | 52 +
.../LocalAssemblerInterface.h | 58 +
.../ThermoRichardsFlowFEM-impl.h | 1352 +++++++++++++++++
.../ThermoRichardsFlowFEM.h | 172 +++
.../ThermoRichardsFlowProcess.cpp | 305 ++++
.../ThermoRichardsFlowProcess.h | 156 ++
.../ThermoRichardsFlowProcessData.h | 53 +
scripts/cmake/ProcessesSetup.cmake | 1 +
12 files changed, 2395 insertions(+), 1 deletion(-)
create mode 100644 ProcessLib/ThermoRichardsFlow/CMakeLists.txt
create mode 100644 ProcessLib/ThermoRichardsFlow/CreateThermoRichardsFlowProcess.cpp
create mode 100644 ProcessLib/ThermoRichardsFlow/CreateThermoRichardsFlowProcess.h
create mode 100644 ProcessLib/ThermoRichardsFlow/IntegrationPointData.h
create mode 100644 ProcessLib/ThermoRichardsFlow/LocalAssemblerInterface.h
create mode 100644 ProcessLib/ThermoRichardsFlow/ThermoRichardsFlowFEM-impl.h
create mode 100644 ProcessLib/ThermoRichardsFlow/ThermoRichardsFlowFEM.h
create mode 100644 ProcessLib/ThermoRichardsFlow/ThermoRichardsFlowProcess.cpp
create mode 100644 ProcessLib/ThermoRichardsFlow/ThermoRichardsFlowProcess.h
create mode 100644 ProcessLib/ThermoRichardsFlow/ThermoRichardsFlowProcessData.h
diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp
index f5c400293ad..c0747fc788f 100644
--- a/Applications/ApplicationsLib/ProjectData.cpp
+++ b/Applications/ApplicationsLib/ProjectData.cpp
@@ -118,6 +118,9 @@
#ifdef OGS_BUILD_PROCESS_THERMOMECHANICS
#include "ProcessLib/ThermoMechanics/CreateThermoMechanicsProcess.h"
#endif
+#ifdef OGS_BUILD_PROCESS_THERMORICHARDSFLOW
+#include "ProcessLib/ThermoRichardsFlow/CreateThermoRichardsFlowProcess.h"
+#endif
#ifdef OGS_BUILD_PROCESS_TWOPHASEFLOWWITHPP
#include "ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPProcess.h"
#endif
@@ -1047,7 +1050,17 @@ void ProjectData::parseProcesses(
}
else
#endif
-
+#ifdef OGS_BUILD_PROCESS_THERMORICHARDSFLOW
+ if (type == "THERMO_RICHARDS_FLOW")
+ {
+ process =
+ ProcessLib::ThermoRichardsFlow::createThermoRichardsFlowProcess(
+ name, *_mesh_vec[0], std::move(jacobian_assembler),
+ _process_variables, _parameters, integration_order,
+ process_config, _media);
+ }
+ else
+#endif
#ifdef OGS_BUILD_PROCESS_THERMORICHARDSMECHANICS
if (type == "THERMO_RICHARDS_MECHANICS")
{
diff --git a/ProcessLib/ThermoRichardsFlow/CMakeLists.txt b/ProcessLib/ThermoRichardsFlow/CMakeLists.txt
new file mode 100644
index 00000000000..aeb8302d1dc
--- /dev/null
+++ b/ProcessLib/ThermoRichardsFlow/CMakeLists.txt
@@ -0,0 +1,8 @@
+get_source_files(SOURCES)
+
+ogs_add_library(ThermoRichardsFlow ${SOURCES})
+target_link_libraries(ThermoRichardsFlow PUBLIC ProcessLib PRIVATE ParameterLib)
+
+if(OGS_BUILD_TESTING)
+ include(Tests.cmake)
+endif()
diff --git a/ProcessLib/ThermoRichardsFlow/CreateThermoRichardsFlowProcess.cpp b/ProcessLib/ThermoRichardsFlow/CreateThermoRichardsFlowProcess.cpp
new file mode 100644
index 00000000000..2ac2c523e5a
--- /dev/null
+++ b/ProcessLib/ThermoRichardsFlow/CreateThermoRichardsFlowProcess.cpp
@@ -0,0 +1,170 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2021, 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 "CreateThermoRichardsFlowProcess.h"
+
+#include <cassert>
+
+#include "CreateSimplifiedElasticityModel.h"
+#include "LocalAssemblerInterface.h"
+#include "MaterialLib/MPL/CreateMaterialSpatialDistributionMap.h"
+#include "MaterialLib/MPL/MaterialSpatialDistributionMap.h"
+#include "MaterialLib/MPL/Medium.h"
+#include "MaterialLib/SolidModels/CreateConstitutiveRelation.h"
+#include "MaterialLib/SolidModels/MechanicsBase.h"
+#include "ParameterLib/Utils.h"
+#include "ProcessLib/Output/CreateSecondaryVariables.h"
+#include "ProcessLib/Utils/ProcessUtils.h"
+#include "SimplifiedElasticityModel.h"
+#include "ThermoRichardsFlowProcess.h"
+#include "ThermoRichardsFlowProcessData.h"
+
+namespace ProcessLib
+{
+namespace ThermoRichardsFlow
+{
+void checkMPLProperties(
+ std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media)
+{
+ std::array const required_medium_properties = {
+ MaterialPropertyLib::permeability,
+ MaterialPropertyLib::porosity, MaterialPropertyLib::biot_coefficient,
+ MaterialPropertyLib::relative_permeability,
+ MaterialPropertyLib::saturation};
+ std::array const required_liquid_properties = {
+ MaterialPropertyLib::viscosity,
+ MaterialPropertyLib::density,
+ };
+ std::array const required_solid_properties = {
+ MaterialPropertyLib::density};
+
+ // Thermal properties are not checked because they can be phase property or
+ // meduim property (will be enabled later).
+ for (auto const& m : media)
+ {
+ checkRequiredProperties(*m.second, required_medium_properties);
+ checkRequiredProperties(m.second->phase("AqueousLiquid"),
+ required_liquid_properties);
+ checkRequiredProperties(m.second->phase("Solid"),
+ required_solid_properties);
+ }
+}
+
+void checkProcessVariableComponents(ProcessVariable const& variable)
+{
+ if (variable.getNumberOfGlobalComponents() != 1)
+ {
+ OGS_FATAL(
+ "Number of components of the process variable '{:s}' is different "
+ "from one: got {:d}.",
+ variable.getName(),
+ variable.getNumberOfGlobalComponents());
+ }
+}
+
+std::unique_ptr<Process> createThermoRichardsFlowProcess(
+ 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,
+ unsigned const integration_order,
+ BaseLib::ConfigTree const& config,
+ std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media)
+{
+ //! \ogs_file_param{prj__processes__process__type}
+ config.checkConfigParameter("type", "THERMO_RICHARDS_FLOW");
+ DBUG("Create ThermoRichardsFlowProcess.");
+
+ auto const staggered_scheme =
+ //! \ogs_file_param{prj__processes__process__THERMO_RICHARDS_FLOW__coupling_scheme}
+ config.getConfigParameterOptional<std::string>("coupling_scheme");
+ const bool use_monolithic_scheme =
+ !(staggered_scheme && (*staggered_scheme == "staggered"));
+
+ // Process variable.
+
+ //! \ogs_file_param{prj__processes__process__THERMO_RICHARDS_FLOW__process_variables}
+ auto const pv_config = config.getConfigSubtree("process_variables");
+
+ ProcessVariable* variable_T;
+ ProcessVariable* variable_p;
+ std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>
+ process_variables;
+ if (use_monolithic_scheme) // monolithic scheme.
+ {
+ auto per_process_variables = findProcessVariables(
+ variables, pv_config,
+ {//! \ogs_file_param_special{prj__processes__process__THERMO_RICHARDS_FLOW__process_variables__temperature}
+ "temperature",
+ //! \ogs_file_param_special{prj__processes__process__THERMO_RICHARDS_FLOW__process_variables__pressure}
+ "pressure"});
+ variable_T = &per_process_variables[0].get();
+ variable_p = &per_process_variables[1].get();
+ process_variables.push_back(std::move(per_process_variables));
+ }
+ else // staggered scheme.
+ {
+ OGS_FATAL(
+ "So far, only the monolithic scheme is implemented for "
+ "THERMO_RICHARDS_FLOW");
+ }
+
+ checkProcessVariableComponents(*variable_T);
+ checkProcessVariableComponents(*variable_p);
+
+ // Specific body force parameter.
+ Eigen::VectorXd specific_body_force;
+ {
+ std::vector<double> const b =
+ //! \ogs_file_param{prj__processes__process__THERMO_RICHARDS_FLOW__specific_body_force}
+ config.getConfigParameter<std::vector<double>>(
+ "specific_body_force");
+ if (b.size() != mesh.getDimension())
+ {
+ OGS_FATAL(
+ "specific body force (gravity vector) has {:d} components, "
+ "but mesh dimension is {:d}",
+ b.size(), mesh.getDimension());
+ }
+ specific_body_force.resize(b.size());
+ std::copy_n(b.data(), b.size(), specific_body_force.data());
+ }
+
+ auto media_map =
+ MaterialPropertyLib::createMaterialSpatialDistributionMap(media, mesh);
+ DBUG(
+ "Check the media properties of ThermoRichardsFlow process "
+ "...");
+ checkMPLProperties(media);
+ DBUG("Media properties verified.");
+
+ //! \ogs_file_param{prj__processes__process__THERMO_RICHARDS_FLOW__mass_lumping}
+ bool const mass_lumping = config.getConfigParameter<bool>("mass_lumping", false);
+
+ std::unique_ptr<SimplifiedElasticityModel> simplified_elasticity =
+ createElasticityModel(config);
+
+ ThermoRichardsFlowProcessData process_data{
+ std::move(media_map), std::move(specific_body_force),
+ mass_lumping, std::move(simplified_elasticity)};
+
+ SecondaryVariableCollection secondary_variables;
+
+ ProcessLib::createSecondaryVariables(config, secondary_variables);
+
+ return std::make_unique<ThermoRichardsFlowProcess>(
+ std::move(name), mesh, std::move(jacobian_assembler), parameters,
+ integration_order, std::move(process_variables),
+ std::move(process_data), std::move(secondary_variables),
+ use_monolithic_scheme);
+}
+} // namespace ThermoRichardsFlow
+} // namespace ProcessLib
diff --git a/ProcessLib/ThermoRichardsFlow/CreateThermoRichardsFlowProcess.h b/ProcessLib/ThermoRichardsFlow/CreateThermoRichardsFlowProcess.h
new file mode 100644
index 00000000000..2abf3aabdc0
--- /dev/null
+++ b/ProcessLib/ThermoRichardsFlow/CreateThermoRichardsFlowProcess.h
@@ -0,0 +1,54 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2021, 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 <map>
+#include <memory>
+#include <vector>
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+namespace MeshLib
+{
+class Mesh;
+}
+namespace MaterialPropertyLib
+{
+class Medium;
+}
+namespace ParameterLib
+{
+struct ParameterBase;
+} // namespace ParameterLib
+namespace ProcessLib
+{
+class AbstractJacobianAssembler;
+class Process;
+class ProcessVariable;
+} // namespace ProcessLib
+
+namespace ProcessLib
+{
+namespace ThermoRichardsFlow
+{
+std::unique_ptr<Process> createThermoRichardsFlowProcess(
+ 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,
+ unsigned const integration_order,
+ BaseLib::ConfigTree const& config,
+ std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media);
+} // namespace ThermoRichardsFlow
+} // namespace ProcessLib
diff --git a/ProcessLib/ThermoRichardsFlow/IntegrationPointData.h b/ProcessLib/ThermoRichardsFlow/IntegrationPointData.h
new file mode 100644
index 00000000000..55eeaf3e8e4
--- /dev/null
+++ b/ProcessLib/ThermoRichardsFlow/IntegrationPointData.h
@@ -0,0 +1,52 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2021, 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>
+
+namespace ProcessLib
+{
+namespace ThermoRichardsFlow
+{
+template <typename ShapeMatricesType>
+struct IntegrationPointData final
+{
+ typename ShapeMatricesType::NodalRowVectorType N;
+ typename ShapeMatricesType::GlobalDimNodalMatrixType dNdx;
+
+ typename ShapeMatricesType::GlobalDimVectorType v_darcy;
+
+ double saturation = std::numeric_limits<double>::quiet_NaN();
+ double saturation_prev = std::numeric_limits<double>::quiet_NaN();
+ double porosity = std::numeric_limits<double>::quiet_NaN();
+ double porosity_prev = std::numeric_limits<double>::quiet_NaN();
+ double transport_porosity = std::numeric_limits<double>::quiet_NaN();
+ double transport_porosity_prev = std::numeric_limits<double>::quiet_NaN();
+ double dry_density_solid = std::numeric_limits<double>::quiet_NaN();
+ double dry_density_pellet_saturated =
+ std::numeric_limits<double>::quiet_NaN();
+ double dry_density_pellet_unsaturated =
+ std::numeric_limits<double>::quiet_NaN();
+
+ double integration_weight;
+
+ void pushBackState()
+ {
+ saturation_prev = saturation;
+ porosity_prev = porosity;
+ transport_porosity_prev = transport_porosity;
+ }
+
+ EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
+};
+
+} // namespace ThermoRichardsFlow
+} // namespace ProcessLib
diff --git a/ProcessLib/ThermoRichardsFlow/LocalAssemblerInterface.h b/ProcessLib/ThermoRichardsFlow/LocalAssemblerInterface.h
new file mode 100644
index 00000000000..d7579227601
--- /dev/null
+++ b/ProcessLib/ThermoRichardsFlow/LocalAssemblerInterface.h
@@ -0,0 +1,58 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2021, 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 <vector>
+
+#include "NumLib/Extrapolation/ExtrapolatableElement.h"
+#include "ProcessLib/LocalAssemblerInterface.h"
+
+namespace ProcessLib
+{
+namespace ThermoRichardsFlow
+{
+struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface,
+ public NumLib::ExtrapolatableElement
+{
+ virtual std::size_t setIPDataInitialConditions(
+ std::string const& name, double const* values,
+ int const integration_order) = 0;
+
+ virtual std::vector<double> const& getIntPtDarcyVelocity(
+ const double t,
+ std::vector<GlobalVector*> const& x,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
+ std::vector<double>& cache) const = 0;
+
+ virtual std::vector<double> getSaturation() const = 0;
+ virtual std::vector<double> const& getIntPtSaturation(
+ const double t,
+ std::vector<GlobalVector*> const& x,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
+ std::vector<double>& cache) const = 0;
+
+ virtual std::vector<double> getPorosity() const = 0;
+ virtual std::vector<double> const& getIntPtPorosity(
+ const double t,
+ std::vector<GlobalVector*> const& x,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
+ std::vector<double>& cache) const = 0;
+
+ virtual std::vector<double> const& getIntPtDryDensitySolid(
+ const double t,
+ std::vector<GlobalVector*> const& x,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
+ std::vector<double>& cache) const = 0;
+
+ // TODO move to NumLib::ExtrapolatableElement
+ virtual unsigned getNumberOfIntegrationPoints() const = 0;
+};
+} // namespace ThermoRichardsFlow
+} // namespace ProcessLib
diff --git a/ProcessLib/ThermoRichardsFlow/ThermoRichardsFlowFEM-impl.h b/ProcessLib/ThermoRichardsFlow/ThermoRichardsFlowFEM-impl.h
new file mode 100644
index 00000000000..255c0b12afe
--- /dev/null
+++ b/ProcessLib/ThermoRichardsFlow/ThermoRichardsFlowFEM-impl.h
@@ -0,0 +1,1352 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2021, 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 <cassert>
+
+#include "HydrostaticElasticityModel.h"
+#include "MaterialLib/MPL/Medium.h"
+#include "MaterialLib/MPL/Utils/FormEffectiveThermalConductivity.h"
+#include "MaterialLib/MPL/Utils/FormEigenTensor.h"
+#include "MaterialLib/MPL/Utils/FormEigenVector.h"
+#include "MaterialLib/MPL/Utils/GetLiquidThermalExpansivity.h"
+#include "MaterialLib/PhysicalConstant.h"
+#include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
+#include "MaterialLib/MPL/MaterialSpatialDistributionMap.h"
+#include "NumLib/Function/Interpolation.h"
+#include "ProcessLib/Utils/SetOrGetIntegrationPointData.h"
+#include "RigidElasticityModel.h"
+#include "UniaxialElasticityModel.h"
+#include "UserDefinedElasticityModel.h"
+
+namespace ProcessLib
+{
+namespace ThermoRichardsFlow
+{
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+ThermoRichardsFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
+ ThermoRichardsFlowLocalAssembler(
+ MeshLib::Element const& e,
+ std::size_t const /*local_matrix_size*/,
+ bool const is_axially_symmetric,
+ unsigned const integration_order,
+ ThermoRichardsFlowProcessData& process_data)
+ : _process_data(process_data),
+ _integration_method(integration_order),
+ _element(e),
+ _is_axially_symmetric(is_axially_symmetric)
+{
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ _ip_data.reserve(n_integration_points);
+
+ auto const shape_matrices =
+ NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType, GlobalDim>(
+ e, is_axially_symmetric, _integration_method);
+
+ auto const& medium = *_process_data.media_map->getMedium(_element.getID());
+
+ ParameterLib::SpatialPosition x_position;
+ x_position.setElementID(_element.getID());
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ auto const& sm = shape_matrices[ip];
+ x_position.setIntegrationPoint(ip);
+ _ip_data.emplace_back();
+ auto& ip_data = _ip_data[ip];
+ _ip_data[ip].integration_weight =
+ _integration_method.getWeightedPoint(ip).getWeight() *
+ sm.integralMeasure * sm.detJ;
+
+ ip_data.N = sm.N;
+ ip_data.dNdx = sm.dNdx;
+
+ // Initial porosity. Could be read from integration point data or mesh.
+ ip_data.porosity =
+ medium[MPL::porosity]
+ .template initialValue<double>(
+ x_position,
+ std::numeric_limits<
+ double>::quiet_NaN() /* t independent */);
+ }
+}
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+std::size_t ThermoRichardsFlowLocalAssembler<
+ ShapeFunction, IntegrationMethod,
+ GlobalDim>::setIPDataInitialConditions(std::string const& name,
+ double const* values,
+ int const integration_order)
+{
+ if (integration_order !=
+ static_cast<int>(_integration_method.getIntegrationOrder()))
+ {
+ OGS_FATAL(
+ "Setting integration point initial conditions; The integration "
+ "order of the local assembler for element {:d} is different "
+ "from the integration order in the initial condition.",
+ _element.getID());
+ }
+
+ if (name == "saturation_ip")
+ {
+ return ProcessLib::setIntegrationPointScalarData(values, _ip_data,
+ &IpData::saturation);
+ }
+ if (name == "porosity_ip")
+ {
+ return ProcessLib::setIntegrationPointScalarData(values, _ip_data,
+ &IpData::porosity);
+ }
+ return 0;
+}
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+void ThermoRichardsFlowLocalAssembler<ShapeFunction, IntegrationMethod,
+ GlobalDim>::
+ setInitialConditionsConcrete(std::vector<double> const& local_x,
+ double const t,
+ bool const /*use_monolithic_scheme*/,
+ int const /*process_id*/)
+{
+ assert(local_x.size() == temperature_size + pressure_size);
+
+ auto p_L = Eigen::Map<
+ typename ShapeMatricesType::template VectorType<pressure_size> const>(
+ local_x.data() + pressure_index, pressure_size);
+
+ auto const& medium = *_process_data.media_map->getMedium(_element.getID());
+ MPL::VariableArray variables;
+
+ ParameterLib::SpatialPosition x_position;
+ x_position.setElementID(_element.getID());
+
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ x_position.setIntegrationPoint(ip);
+
+ auto const& N = _ip_data[ip].N;
+
+ double p_cap_ip;
+ NumLib::shapeFunctionInterpolate(-p_L, N, p_cap_ip);
+
+ variables[static_cast<int>(MPL::Variable::capillary_pressure)] =
+ p_cap_ip;
+ variables[static_cast<int>(MPL::Variable::phase_pressure)] = -p_cap_ip;
+
+ // Note: temperature dependent saturation model is not considered so
+ // far.
+ _ip_data[ip].saturation_prev =
+ medium[MPL::PropertyType::saturation]
+ .template value<double>(
+ variables, x_position, t,
+ std::numeric_limits<double>::quiet_NaN());
+ }
+}
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+void ThermoRichardsFlowLocalAssembler<
+ ShapeFunction, IntegrationMethod,
+ GlobalDim>::assembleWithJacobian(double const t, double const dt,
+ std::vector<double> const& local_x,
+ std::vector<double> const& local_xdot,
+ const double /*dxdot_dx*/,
+ const double /*dx_dx*/,
+ std::vector<double>& /*local_M_data*/,
+ std::vector<double>& /*local_K_data*/,
+ std::vector<double>& local_rhs_data,
+ std::vector<double>& local_Jac_data)
+{
+ auto const local_matrix_dim = pressure_size + temperature_size;
+ assert(local_x.size() == local_matrix_dim);
+
+ auto const T = Eigen::Map<typename ShapeMatricesType::template VectorType<
+ temperature_size> const>(local_x.data() + temperature_index,
+ temperature_size);
+ auto const p_L = Eigen::Map<
+ typename ShapeMatricesType::template VectorType<pressure_size> const>(
+ local_x.data() + pressure_index, pressure_size);
+
+ auto const T_dot =
+ Eigen::Map<typename ShapeMatricesType::template VectorType<
+ temperature_size> const>(local_xdot.data() + temperature_index,
+ temperature_size);
+ auto const p_L_dot = Eigen::Map<
+ typename ShapeMatricesType::template VectorType<pressure_size> const>(
+ local_xdot.data() + pressure_index, pressure_size);
+
+ auto local_Jac = MathLib::createZeroedMatrix<
+ typename ShapeMatricesType::template MatrixType<local_matrix_dim,
+ local_matrix_dim>>(
+ local_Jac_data, local_matrix_dim, local_matrix_dim);
+
+ auto local_rhs = MathLib::createZeroedVector<
+ typename ShapeMatricesType::template VectorType<local_matrix_dim>>(
+ local_rhs_data, local_matrix_dim);
+
+ typename ShapeMatricesType::NodalMatrixType M_TT =
+ ShapeMatricesType::NodalMatrixType::Zero(temperature_size,
+ temperature_size);
+ typename ShapeMatricesType::NodalMatrixType M_Tp =
+ ShapeMatricesType::NodalMatrixType::Zero(temperature_size,
+ pressure_size);
+ typename ShapeMatricesType::NodalMatrixType K_TT =
+ ShapeMatricesType::NodalMatrixType::Zero(temperature_size,
+ temperature_size);
+ typename ShapeMatricesType::NodalMatrixType K_Tp =
+ ShapeMatricesType::NodalMatrixType::Zero(temperature_size,
+ pressure_size);
+ typename ShapeMatricesType::NodalMatrixType M_pT =
+ ShapeMatricesType::NodalMatrixType::Zero(pressure_size,
+ temperature_size);
+ typename ShapeMatricesType::NodalMatrixType laplace_p =
+ ShapeMatricesType::NodalMatrixType::Zero(pressure_size, pressure_size);
+
+ typename ShapeMatricesType::NodalMatrixType storage_p_a_p =
+ ShapeMatricesType::NodalMatrixType::Zero(pressure_size, pressure_size);
+
+ typename ShapeMatricesType::NodalMatrixType storage_p_a_S_Jpp =
+ ShapeMatricesType::NodalMatrixType::Zero(pressure_size, pressure_size);
+
+ typename ShapeMatricesType::NodalMatrixType storage_p_a_S =
+ ShapeMatricesType::NodalMatrixType::Zero(pressure_size, pressure_size);
+
+ auto const& medium = *_process_data.media_map->getMedium(_element.getID());
+ auto const& liquid_phase = medium.phase("AqueousLiquid");
+ auto const& solid_phase = medium.phase("Solid");
+ MPL::Phase const* gas_phase = medium.hasPhase("Gas") ?
+ &medium.phase("Gas") : nullptr;
+ MPL::VariableArray variables;
+ MPL::VariableArray variables_prev;
+
+ ParameterLib::SpatialPosition x_position;
+ x_position.setElementID(_element.getID());
+
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ x_position.setIntegrationPoint(ip);
+ auto const& w = _ip_data[ip].integration_weight;
+
+ auto const& N = _ip_data[ip].N;
+ auto const& dNdx = _ip_data[ip].dNdx;
+
+ double T_ip;
+ NumLib::shapeFunctionInterpolate(T, N, T_ip);
+ double T_dot_ip;
+ NumLib::shapeFunctionInterpolate(T_dot, N, T_dot_ip);
+
+ double p_cap_ip;
+ NumLib::shapeFunctionInterpolate(-p_L, N, p_cap_ip);
+
+ double p_cap_dot_ip;
+ NumLib::shapeFunctionInterpolate(-p_L_dot, N, p_cap_dot_ip);
+
+ variables[static_cast<int>(MPL::Variable::capillary_pressure)] =
+ p_cap_ip;
+ variables[static_cast<int>(MPL::Variable::phase_pressure)] = -p_cap_ip;
+ variables[static_cast<int>(MPL::Variable::temperature)] = T_ip;
+
+ auto& S_L = _ip_data[ip].saturation;
+ auto const S_L_prev = _ip_data[ip].saturation_prev;
+ auto const alpha =
+ medium[MPL::PropertyType::biot_coefficient]
+ .template value<double>(variables, x_position, t, dt);
+
+ auto& solid_elasticity = *_process_data.simplified_elasticity;
+ // TODO (buchwaldj)
+ // is bulk_modulus good name for bulk modulus of solid skeleton?
+ auto const beta_S =
+ solid_elasticity.bulkCompressibilityFromYoungsModulus(
+ solid_phase, variables, x_position, t, dt);
+ auto const beta_SR = (1 - alpha) * beta_S;
+ variables[static_cast<int>(MPL::Variable::grain_compressibility)] =
+ beta_SR;
+
+ auto const rho_LR =
+ liquid_phase[MPL::PropertyType::density]
+ .template value<double>(variables, x_position, t, dt);
+ auto const& b = _process_data.specific_body_force;
+
+ double const drho_LR_dp =
+ liquid_phase[MPL::PropertyType::density]
+ .template dValue<double>(variables,
+ MPL::Variable::phase_pressure,
+ x_position, t, dt);
+ auto const beta_LR = drho_LR_dp / rho_LR;
+
+ S_L = medium[MPL::PropertyType::saturation]
+ .template value<double>(variables, x_position, t, dt);
+ variables[static_cast<int>(MPL::Variable::liquid_saturation)] = S_L;
+ variables_prev[static_cast<int>(MPL::Variable::liquid_saturation)] =
+ S_L_prev;
+
+ // tangent derivative for Jacobian
+ double const dS_L_dp_cap = medium[MPL::PropertyType::saturation]
+ .template dValue<double>(variables,
+ MPL::Variable::capillary_pressure,
+ x_position, t, dt);
+ // secant derivative from time discretization for storage
+ // use tangent, if secant is not available
+ double const DeltaS_L_Deltap_cap =
+ (p_cap_dot_ip == 0) ? dS_L_dp_cap
+ : (S_L - S_L_prev) / (dt * p_cap_dot_ip);
+
+ auto chi_S_L = S_L;
+ auto chi_S_L_prev = S_L_prev;
+ auto dchi_dS_L = 1.0;
+ if (medium.hasProperty(MPL::PropertyType::bishops_effective_stress))
+ {
+ auto const chi = [&medium, x_position, t, dt](double const S_L) {
+ MPL::VariableArray variables;
+ variables[static_cast<int>(MPL::Variable::liquid_saturation)] =
+ S_L;
+ return medium[MPL::PropertyType::bishops_effective_stress]
+ .template value<double>(variables, x_position, t, dt);
+ };
+ chi_S_L = chi(S_L);
+ chi_S_L_prev = chi(S_L_prev);
+
+ dchi_dS_L = medium[MPL::PropertyType::bishops_effective_stress]
+ .template dValue<double>(variables,
+ MPL::Variable::liquid_saturation,
+ x_position, t, dt);
+ }
+ // TODO (buchwaldj)
+ // should solid_grain_pressure or effective_pore_pressure remain?
+ // double const p_FR = -chi_S_L * p_cap_ip;
+ // variables[static_cast<int>(MPL::Variable::solid_grain_pressure)] =
+ // p_FR;
+
+ variables[static_cast<int>(MPL::Variable::effective_pore_pressure)] =
+ -chi_S_L * p_cap_ip;
+ variables_prev[static_cast<int>(
+ MPL::Variable::effective_pore_pressure)] =
+ -chi_S_L_prev * (p_cap_ip - p_cap_dot_ip * dt);
+
+ auto& phi = _ip_data[ip].porosity;
+ { // Porosity update
+
+ variables_prev[static_cast<int>(MPL::Variable::porosity)] =
+ _ip_data[ip].porosity_prev;
+ phi = medium[MPL::PropertyType::porosity]
+ .template value<double>(variables, variables_prev,
+ x_position, t, dt);
+ variables[static_cast<int>(MPL::Variable::porosity)] = phi;
+ }
+
+ if (alpha < phi)
+ {
+ OGS_FATAL(
+ "ThermoRichardsFlow: Biot-coefficient {} is smaller than "
+ "porosity {} in element/integration point {}/{}.",
+ alpha, phi, _element.getID(), ip);
+ }
+
+ double const k_rel =
+ medium[MPL::PropertyType::relative_permeability]
+ .template value<double>(variables, x_position, t, dt);
+ auto const mu =
+ liquid_phase[MPL::PropertyType::viscosity]
+ .template value<double>(variables, x_position, t, dt);
+
+ auto const K_intrinsic = MPL::formEigenTensor<GlobalDim>(
+ medium[MPL::PropertyType::permeability]
+ .value(variables, x_position, t, dt));
+
+ GlobalDimMatrixType const Ki_over_mu = K_intrinsic / mu;
+ GlobalDimMatrixType const rho_Ki_over_mu = rho_LR * Ki_over_mu;
+
+ // Consider anisotropic thermal expansion.
+ // Read in 3x3 tensor. 2D case also requires expansion coeff. for z-
+ // component.
+ Eigen::Matrix<double, 3,
+ 3> const solid_linear_thermal_expansion_coefficient =
+ MaterialPropertyLib::formEigenTensor<3>(
+ solid_phase[
+ MaterialPropertyLib::PropertyType::thermal_expansivity]
+ .value(variables, x_position, t, dt));
+
+ auto const rho_SR =
+ solid_phase[MPL::PropertyType::density]
+ .template value<double>(variables, x_position, t, dt);
+
+ //
+ // pressure equation, pressure part.
+ //
+ laplace_p.noalias() +=
+ dNdx.transpose() * k_rel * rho_Ki_over_mu * dNdx * w;
+
+ const double alphaB_minus_phi = alpha - phi;
+ double const a0 = alphaB_minus_phi * beta_SR;
+ double const specific_storage_a_p =
+ S_L * (phi * beta_LR + S_L * a0 +
+ chi_S_L * alpha * alpha *
+ solid_elasticity.storageContribution(
+ solid_phase, variables, x_position, t, dt));
+ double const specific_storage_a_S = phi - p_cap_ip * S_L * a0;
+
+ double const dspecific_storage_a_p_dp_cap =
+ dS_L_dp_cap * (phi * beta_LR + 2 * S_L * a0 +
+ alpha * alpha *
+ solid_elasticity.storageContribution(
+ solid_phase, variables, x_position, t, dt) *
+ (chi_S_L + dchi_dS_L * S_L));
+ double const dspecific_storage_a_S_dp_cap =
+ -a0 * (S_L + p_cap_ip * dS_L_dp_cap);
+
+ storage_p_a_p.noalias() +=
+ N.transpose() * rho_LR * specific_storage_a_p * N * w;
+
+ storage_p_a_S.noalias() -= N.transpose() * rho_LR *
+ specific_storage_a_S * DeltaS_L_Deltap_cap *
+ N * w;
+
+ local_Jac
+ .template block<pressure_size, pressure_size>(pressure_index,
+ pressure_index)
+ .noalias() += N.transpose() * p_cap_dot_ip * rho_LR *
+ dspecific_storage_a_p_dp_cap * N * w;
+
+ storage_p_a_S_Jpp.noalias() -=
+ N.transpose() * rho_LR *
+ ((S_L - S_L_prev) * dspecific_storage_a_S_dp_cap +
+ specific_storage_a_S * dS_L_dp_cap) /
+ dt * N * w;
+
+ double const dk_rel_dS_L =
+ medium[MPL::PropertyType::relative_permeability]
+ .template dValue<double>(variables,
+ MPL::Variable::liquid_saturation,
+ x_position, t, dt);
+ GlobalDimVectorType const grad_p_cap = -dNdx * p_L;
+ local_Jac
+ .template block<pressure_size, pressure_size>(pressure_index,
+ pressure_index)
+ .noalias() += dNdx.transpose() * rho_Ki_over_mu * grad_p_cap *
+ dk_rel_dS_L * dS_L_dp_cap * N * w;
+
+ local_Jac
+ .template block<pressure_size, pressure_size>(pressure_index,
+ pressure_index)
+ .noalias() += dNdx.transpose() * rho_LR * rho_Ki_over_mu * b *
+ dk_rel_dS_L * dS_L_dp_cap * N * w;
+
+ local_rhs.template segment<pressure_size>(pressure_index).noalias() +=
+ dNdx.transpose() * rho_LR * k_rel * rho_Ki_over_mu * b * w;
+
+ //
+ // pressure equation, temperature part.
+ //
+ double const fluid_volumetric_thermal_expansion_coefficient =
+ MPL::getLiquidThermalExpansivity(liquid_phase, variables, rho_LR,
+ x_position, t, dt);
+ const double eff_thermal_expansion =
+ S_L * (alphaB_minus_phi *
+ solid_linear_thermal_expansion_coefficient.trace() +
+ phi * fluid_volumetric_thermal_expansion_coefficient +
+ alpha * solid_elasticity.thermalExpansivityContribution(
+ solid_linear_thermal_expansion_coefficient,
+ solid_phase, variables, x_position, t, dt));
+ M_pT.noalias() -=
+ N.transpose() * rho_LR * eff_thermal_expansion * N * w;
+
+ //
+ // temperature equation.
+ //
+ {
+ auto const specific_heat_capacity_fluid =
+ liquid_phase[MaterialPropertyLib::specific_heat_capacity]
+ .template value<double>(variables, x_position, t, dt);
+
+ auto const specific_heat_capacity_solid =
+ solid_phase
+ [MaterialPropertyLib::PropertyType::
+ specific_heat_capacity]
+ .template value<double>(variables, x_position, t, dt);
+
+ M_TT.noalias() +=
+ w *
+ (rho_SR * specific_heat_capacity_solid * (1 - phi) +
+ (S_L * rho_LR * specific_heat_capacity_fluid) * phi) *
+ N.transpose() * N;
+
+ auto const thermal_conductivity =
+ MaterialPropertyLib::formEigenTensor<GlobalDim>(
+ medium[MaterialPropertyLib::PropertyType::
+ thermal_conductivity]
+ .value(variables, x_position, t, dt));
+
+ GlobalDimVectorType const velocity_L =
+ GlobalDimVectorType(-Ki_over_mu * (dNdx * p_L - rho_LR * b));
+
+ K_TT.noalias() +=
+ (dNdx.transpose() * thermal_conductivity * dNdx +
+ N.transpose() * velocity_L.transpose() * dNdx * rho_LR *
+ specific_heat_capacity_fluid) *
+ w;
+
+ //
+ // temperature equation, pressure part
+ //
+ K_Tp.noalias() -= rho_LR * specific_heat_capacity_fluid *
+ N.transpose() * (dNdx * T).transpose() *
+ Ki_over_mu * dNdx * w;
+ }
+ if (liquid_phase.hasProperty(MPL::PropertyType::vapour_diffusion) &&
+ S_L < 1.0)
+ {
+ variables[static_cast<int>(MPL::Variable::density)] = rho_LR;
+
+ double const rho_wv =
+ liquid_phase[MaterialPropertyLib::vapour_density]
+ .template value<double>(variables, x_position, t, dt);
+
+ double const drho_wv_dT =
+ liquid_phase[MaterialPropertyLib::vapour_density]
+ .template dValue<double>(variables,
+ MPL::Variable::temperature,
+ x_position, t, dt);
+ double const drho_wv_dp =
+ liquid_phase[MaterialPropertyLib::vapour_density]
+ .template dValue<double>(variables,
+ MPL::Variable::phase_pressure,
+ x_position, t, dt);
+ auto const f_Tv =
+ liquid_phase[
+ MPL::PropertyType::thermal_diffusion_enhancement_factor]
+ .template value<double>(variables, x_position, t, dt);
+
+ variables[static_cast<int>(MPL::Variable::porosity)] = phi;
+ double const D_v =
+ liquid_phase[MPL::PropertyType::vapour_diffusion]
+ .template value<double>(variables, x_position, t, dt);
+
+ double const f_Tv_D_Tv = f_Tv * D_v * drho_wv_dT;
+ double const D_pv = D_v * drho_wv_dp;
+
+ if (gas_phase &&
+ gas_phase->hasProperty(MPL::PropertyType::heat_capacity))
+ {
+ GlobalDimVectorType const grad_T = dNdx * T;
+ // Vapour velocity
+ GlobalDimVectorType const q_v =
+ -(f_Tv_D_Tv * grad_T - D_pv * grad_p_cap) / rho_LR;
+ double const specific_heat_capacity_vapour =
+ gas_phase->property(MaterialPropertyLib::PropertyType::
+ specific_heat_capacity)
+ .template value<double>(variables, x_position, t, dt);
+
+ M_TT.noalias() +=
+ w *
+ (rho_wv * specific_heat_capacity_vapour * (1 - S_L) * phi) *
+ N.transpose() * N;
+
+ K_TT.noalias() += N.transpose() * q_v.transpose() * dNdx *
+ rho_wv * specific_heat_capacity_vapour * w;
+ }
+
+ double const storage_coefficient_by_water_vapor =
+ phi * (rho_wv * dS_L_dp_cap + (1 - S_L) * drho_wv_dp);
+
+ storage_p_a_p.noalias() +=
+ N.transpose() * storage_coefficient_by_water_vapor * N * w;
+
+ double const vapor_expansion_factor = phi * (1 - S_L) * drho_wv_dT;
+ M_pT.noalias() += N.transpose() * vapor_expansion_factor * N * w;
+
+ local_Jac
+ .template block<pressure_size, temperature_size>(
+ pressure_index, temperature_index)
+ .noalias() += dNdx.transpose() * f_Tv_D_Tv * dNdx * w;
+
+ local_rhs.template segment<pressure_size>(pressure_index)
+ .noalias() -= f_Tv_D_Tv * dNdx.transpose() * (dNdx * T) * w;
+
+ laplace_p.noalias() += dNdx.transpose() * D_pv * dNdx * w;
+
+ //
+ // Latent heat term
+ //
+ if (liquid_phase.hasProperty(MPL::PropertyType::latent_heat))
+ {
+ double const factor = phi * (1 - S_L) / rho_LR;
+ // The volumetric latent heat of vaporization of liquid water
+ double const L0 =
+ liquid_phase[MPL::PropertyType::latent_heat]
+ .template value<double>(variables, x_position, t, dt) *
+ rho_LR;
+
+ double const drho_LR_dT =
+ liquid_phase[MPL::PropertyType::density]
+ .template dValue<double>(variables,
+ MPL::Variable::temperature,
+ x_position, t, dt);
+
+ double const rho_wv_over_rho_L = rho_wv / rho_LR;
+ M_TT.noalias() +=
+ factor * L0 *
+ (drho_wv_dT - rho_wv_over_rho_L * drho_LR_dT) *
+ N.transpose() * N * w;
+
+ M_Tp.noalias() +=
+ (factor * L0 *
+ (drho_wv_dp - rho_wv_over_rho_L * drho_LR_dp) +
+ L0 * phi * rho_wv_over_rho_L * dS_L_dp_cap) *
+ N.transpose() * N * w;
+
+ // temperature equation, temperature part
+ K_TT.noalias() +=
+ L0 * f_Tv_D_Tv * dNdx.transpose() * dNdx * w / rho_LR;
+ // temperature equation, pressure part
+ K_Tp.noalias() +=
+ L0 * D_pv * dNdx.transpose() * dNdx * w / rho_LR;
+ }
+ }
+ }
+
+ if (_process_data.apply_mass_lumping)
+ {
+ storage_p_a_p = storage_p_a_p.colwise().sum().eval().asDiagonal();
+ storage_p_a_S = storage_p_a_S.colwise().sum().eval().asDiagonal();
+ storage_p_a_S_Jpp =
+ storage_p_a_S_Jpp.colwise().sum().eval().asDiagonal();
+ }
+
+ //
+ // -- Jacobian
+ //
+ // temperature equation.
+ local_Jac
+ .template block<temperature_size, temperature_size>(temperature_index,
+ temperature_index)
+ .noalias() += M_TT / dt + K_TT;
+ // temperature equation, pressure part
+ local_Jac
+ .template block<temperature_size, pressure_size>(temperature_index,
+ pressure_index)
+ .noalias() += K_Tp;
+
+ // pressure equation, pressure part.
+ local_Jac
+ .template block<pressure_size, pressure_size>(pressure_index,
+ pressure_index)
+ .noalias() += laplace_p + storage_p_a_p / dt + storage_p_a_S_Jpp;
+
+ // pressure equation, temperature part (contributed by thermal expansion).
+ local_Jac
+ .template block<pressure_size, temperature_size>(pressure_index,
+ temperature_index)
+ .noalias() += M_pT / dt;
+
+ //
+ // -- Residual
+ //
+ // temperature equation
+ local_rhs.template segment<temperature_size>(temperature_index).noalias() -=
+ M_TT * T_dot + K_TT * T;
+
+ // pressure equation
+ local_rhs.template segment<pressure_size>(pressure_index).noalias() -=
+ laplace_p * p_L + (storage_p_a_p + storage_p_a_S) * p_L_dot +
+ M_pT * T_dot;
+ if (liquid_phase.hasProperty(MPL::PropertyType::vapour_diffusion) &&
+ liquid_phase.hasProperty(MPL::PropertyType::latent_heat))
+ {
+ // Jacobian: temperature equation, pressure part
+ local_Jac
+ .template block<temperature_size, pressure_size>(temperature_index,
+ pressure_index)
+ .noalias() += M_Tp / dt;
+ // RHS: temperature part
+ local_rhs.template segment<temperature_size>(temperature_index)
+ .noalias() -= M_Tp * p_L_dot;
+ }
+}
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+void ThermoRichardsFlowLocalAssembler<
+ ShapeFunction, IntegrationMethod,
+ GlobalDim>::assemble(double const t, double const dt,
+ std::vector<double> const& local_x,
+ std::vector<double> const& local_xdot,
+ std::vector<double>& local_M_data,
+ std::vector<double>& local_K_data,
+ std::vector<double>& local_rhs_data)
+{
+ auto const local_matrix_dim = pressure_size + temperature_size;
+ assert(local_x.size() == local_matrix_dim);
+
+ auto const T = Eigen::Map<typename ShapeMatricesType::template VectorType<
+ temperature_size> const>(local_x.data() + temperature_index,
+ temperature_size);
+ auto const p_L = Eigen::Map<
+ typename ShapeMatricesType::template VectorType<pressure_size> const>(
+ local_x.data() + pressure_index, pressure_size);
+
+ auto const T_dot =
+ Eigen::Map<typename ShapeMatricesType::template VectorType<
+ temperature_size> const>(local_xdot.data() + temperature_index,
+ temperature_size);
+ auto const p_L_dot = Eigen::Map<
+ typename ShapeMatricesType::template VectorType<pressure_size> const>(
+ local_xdot.data() + pressure_index, pressure_size);
+
+ auto local_K = MathLib::createZeroedMatrix<
+ typename ShapeMatricesType::template MatrixType<
+ local_matrix_dim, local_matrix_dim>>(
+ local_K_data, local_matrix_dim, local_matrix_dim);
+
+ auto local_M = MathLib::createZeroedMatrix<
+ typename ShapeMatricesType::template MatrixType<
+ local_matrix_dim, local_matrix_dim>>(
+ local_M_data, local_matrix_dim, local_matrix_dim);
+
+ auto local_rhs = MathLib::createZeroedVector<
+ typename ShapeMatricesType::template VectorType<local_matrix_dim>>(
+ local_rhs_data, local_matrix_dim);
+
+ auto const& medium = *_process_data.media_map->getMedium(_element.getID());
+ auto const& liquid_phase = medium.phase("AqueousLiquid");
+ auto const& solid_phase = medium.phase("Solid");
+ MPL::Phase const* gas_phase = medium.hasPhase("Gas") ?
+ &medium.phase("Gas") : nullptr;
+ MPL::VariableArray variables;
+ MPL::VariableArray variables_prev;
+
+ ParameterLib::SpatialPosition x_position;
+ x_position.setElementID(_element.getID());
+
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ x_position.setIntegrationPoint(ip);
+ auto const& w = _ip_data[ip].integration_weight;
+
+ auto const& N = _ip_data[ip].N;
+ auto const& dNdx = _ip_data[ip].dNdx;
+
+ double T_ip;
+ NumLib::shapeFunctionInterpolate(T, N, T_ip);
+ double T_dot_ip;
+ NumLib::shapeFunctionInterpolate(T_dot, N, T_dot_ip);
+
+ double p_cap_ip;
+ NumLib::shapeFunctionInterpolate(-p_L, N, p_cap_ip);
+
+ double p_cap_dot_ip;
+ NumLib::shapeFunctionInterpolate(-p_L_dot, N, p_cap_dot_ip);
+
+ variables[static_cast<int>(MPL::Variable::capillary_pressure)] =
+ p_cap_ip;
+ variables[static_cast<int>(MPL::Variable::phase_pressure)] = -p_cap_ip;
+ variables[static_cast<int>(MPL::Variable::temperature)] = T_ip;
+
+ auto& S_L = _ip_data[ip].saturation;
+ auto const S_L_prev = _ip_data[ip].saturation_prev;
+ auto const alpha =
+ medium[MPL::PropertyType::biot_coefficient]
+ .template value<double>(variables, x_position, t, dt);
+
+ auto& solid_elasticity = *_process_data.simplified_elasticity;
+ // TODO (buchwaldj)
+ // is bulk_modulus good name for bulk modulus of solid skeleton?
+ auto const beta_S =
+ solid_elasticity.bulkCompressibilityFromYoungsModulus(
+ solid_phase, variables, x_position, t, dt);
+ auto const beta_SR = (1 - alpha) * beta_S;
+ variables[static_cast<int>(MPL::Variable::grain_compressibility)] =
+ beta_SR;
+
+ auto const rho_LR =
+ liquid_phase[MPL::PropertyType::density]
+ .template value<double>(variables, x_position, t, dt);
+ auto const& b = _process_data.specific_body_force;
+
+ double const drho_LR_dp =
+ liquid_phase[MPL::PropertyType::density]
+ .template dValue<double>(variables,
+ MPL::Variable::phase_pressure,
+ x_position, t, dt);
+ auto const beta_LR = drho_LR_dp / rho_LR;
+
+ S_L = medium[MPL::PropertyType::saturation]
+ .template value<double>(variables, x_position, t, dt);
+ variables[static_cast<int>(MPL::Variable::liquid_saturation)] = S_L;
+ variables_prev[static_cast<int>(MPL::Variable::liquid_saturation)] =
+ S_L_prev;
+
+ // tangent derivative for Jacobian
+ double const dS_L_dp_cap = medium[MPL::PropertyType::saturation]
+ .template dValue<double>(variables,
+ MPL::Variable::capillary_pressure,
+ x_position, t, dt);
+ // secant derivative from time discretization for storage
+ // use tangent, if secant is not available
+ double const DeltaS_L_Deltap_cap =
+ (p_cap_dot_ip == 0) ? dS_L_dp_cap
+ : (S_L - S_L_prev) / (dt * p_cap_dot_ip);
+
+ auto chi_S_L = S_L;
+ auto chi_S_L_prev = S_L_prev;
+ if (medium.hasProperty(MPL::PropertyType::bishops_effective_stress))
+ {
+ auto const chi = [&medium, x_position, t, dt](double const S_L) {
+ MPL::VariableArray variables;
+ variables[static_cast<int>(MPL::Variable::liquid_saturation)] =
+ S_L;
+ return medium[MPL::PropertyType::bishops_effective_stress]
+ .template value<double>(variables, x_position, t, dt);
+ };
+ chi_S_L = chi(S_L);
+ chi_S_L_prev = chi(S_L_prev);
+
+ }
+ // TODO (buchwaldj)
+ // should solid_grain_pressure or effective_pore_pressure remain?
+ // double const p_FR = -chi_S_L * p_cap_ip;
+ // variables[static_cast<int>(MPL::Variable::solid_grain_pressure)] =
+ // p_FR;
+
+ variables[static_cast<int>(MPL::Variable::effective_pore_pressure)] =
+ -chi_S_L * p_cap_ip;
+ variables_prev[static_cast<int>(
+ MPL::Variable::effective_pore_pressure)] =
+ -chi_S_L_prev * (p_cap_ip - p_cap_dot_ip * dt);
+
+ auto& phi = _ip_data[ip].porosity;
+ { // Porosity update
+
+ variables_prev[static_cast<int>(MPL::Variable::porosity)] =
+ _ip_data[ip].porosity_prev;
+ phi = medium[MPL::PropertyType::porosity]
+ .template value<double>(variables, variables_prev,
+ x_position, t, dt);
+ variables[static_cast<int>(MPL::Variable::porosity)] = phi;
+ }
+
+ if (alpha < phi)
+ {
+ OGS_FATAL(
+ "ThermoRichardsFlow: Biot-coefficient {} is smaller than "
+ "porosity {} in element/integration point {}/{}.",
+ alpha, phi, _element.getID(), ip);
+ }
+
+ double const k_rel =
+ medium[MPL::PropertyType::relative_permeability]
+ .template value<double>(variables, x_position, t, dt);
+ auto const mu =
+ liquid_phase[MPL::PropertyType::viscosity]
+ .template value<double>(variables, x_position, t, dt);
+
+ auto const K_intrinsic = MPL::formEigenTensor<GlobalDim>(
+ medium[MPL::PropertyType::permeability]
+ .value(variables, x_position, t, dt));
+
+ GlobalDimMatrixType const Ki_over_mu = K_intrinsic / mu;
+ GlobalDimMatrixType const rho_Ki_over_mu = rho_LR * Ki_over_mu;
+
+ // Consider anisotropic thermal expansion.
+ // Read in 3x3 tensor. 2D case also requires expansion coeff. for z-
+ // component.
+ Eigen::Matrix<double, 3,
+ 3> const solid_linear_thermal_expansion_coefficient =
+ MaterialPropertyLib::formEigenTensor<3>(
+ solid_phase[
+ MaterialPropertyLib::PropertyType::thermal_expansivity]
+ .value(variables, x_position, t, dt));
+
+ auto const rho_SR =
+ solid_phase[MPL::PropertyType::density]
+ .template value<double>(variables, x_position, t, dt);
+
+ //
+ // pressure equation, pressure part.
+ //
+ local_K.template block<pressure_size, pressure_size>(pressure_index,
+ pressure_index)
+ .noalias() += dNdx.transpose() * k_rel * rho_Ki_over_mu * dNdx * w;
+
+ const double alphaB_minus_phi = alpha - phi;
+ double const a0 = alphaB_minus_phi * beta_SR;
+ double const specific_storage_a_p =
+ S_L * (phi * beta_LR + S_L * a0 +
+ chi_S_L * alpha * alpha *
+ solid_elasticity.storageContribution(
+ solid_phase, variables, x_position, t, dt));
+ double const specific_storage_a_S = phi - p_cap_ip * S_L * a0;
+
+
+ local_M.template block<pressure_size, pressure_size>(pressure_index,
+ pressure_index)
+ .noalias() += N.transpose() * rho_LR * (specific_storage_a_p -
+ specific_storage_a_S * DeltaS_L_Deltap_cap) * N * w;
+
+
+ local_rhs.template segment<pressure_size>(pressure_index).noalias() +=
+ dNdx.transpose() * rho_LR * k_rel * rho_Ki_over_mu * b * w;
+
+ //
+ // pressure equation, temperature part.
+ //
+ double const fluid_volumetric_thermal_expansion_coefficient =
+ MPL::getLiquidThermalExpansivity(liquid_phase, variables, rho_LR,
+ x_position, t, dt);
+ const double eff_thermal_expansion =
+ S_L * (alphaB_minus_phi *
+ solid_linear_thermal_expansion_coefficient.trace() +
+ phi * fluid_volumetric_thermal_expansion_coefficient +
+ alpha * solid_elasticity.thermalExpansivityContribution(
+ solid_linear_thermal_expansion_coefficient,
+ solid_phase, variables, x_position, t, dt));
+
+ local_M.template block<pressure_size, temperature_size>(pressure_index,
+ temperature_index)
+ .noalias() -=
+ N.transpose() * rho_LR * eff_thermal_expansion * N * w;
+
+ //
+ // temperature equation.
+ //
+ {
+ auto const specific_heat_capacity_fluid =
+ liquid_phase[MaterialPropertyLib::specific_heat_capacity]
+ .template value<double>(variables, x_position, t, dt);
+
+ auto const specific_heat_capacity_solid =
+ solid_phase[MaterialPropertyLib::PropertyType::
+ specific_heat_capacity]
+ .template value<double>(variables, x_position, t, dt);
+
+ local_M.template block<temperature_size, temperature_size>(temperature_index,
+ temperature_index)
+ .noalias() +=
+ w *
+ (rho_SR * specific_heat_capacity_solid * (1 - phi) +
+ (S_L * rho_LR * specific_heat_capacity_fluid) * phi) *
+ N.transpose() * N;
+
+ auto const thermal_conductivity =
+ MaterialPropertyLib::formEigenTensor<GlobalDim>(
+ medium[MaterialPropertyLib::PropertyType::
+ thermal_conductivity]
+ .value(variables, x_position, t, dt));
+
+ GlobalDimVectorType const velocity_L =
+ GlobalDimVectorType(-Ki_over_mu * (dNdx * p_L - rho_LR * b));
+
+ local_K.template block<temperature_size, temperature_size>(temperature_index,
+ temperature_index)
+ .noalias() +=
+ (dNdx.transpose() * thermal_conductivity * dNdx +
+ N.transpose() * velocity_L.transpose() * dNdx * rho_LR *
+ specific_heat_capacity_fluid) *
+ w;
+
+ //
+ // temperature equation, pressure part
+ //
+ local_K.template block<temperature_size, pressure_size>(temperature_index,
+ pressure_index)
+ .noalias() -=
+ rho_LR * specific_heat_capacity_fluid *
+ N.transpose() * (dNdx * T).transpose() *
+ Ki_over_mu * dNdx * w;
+ }
+ if (liquid_phase.hasProperty(MPL::PropertyType::vapour_diffusion) &&
+ S_L < 1.0)
+ {
+ variables[static_cast<int>(MPL::Variable::density)] = rho_LR;
+
+ double const rho_wv =
+ liquid_phase[MaterialPropertyLib::vapour_density]
+ .template value<double>(variables, x_position, t, dt);
+
+ double const drho_wv_dT =
+ liquid_phase[MaterialPropertyLib::vapour_density]
+ .template dValue<double>(variables,
+ MPL::Variable::temperature,
+ x_position, t, dt);
+ double const drho_wv_dp =
+ liquid_phase[MaterialPropertyLib::vapour_density]
+ .template dValue<double>(variables,
+ MPL::Variable::phase_pressure,
+ x_position, t, dt);
+ auto const f_Tv =
+ liquid_phase[
+ MPL::PropertyType::thermal_diffusion_enhancement_factor]
+ .template value<double>(variables, x_position, t, dt);
+
+ variables[static_cast<int>(MPL::Variable::porosity)] = phi;
+ double const D_v =
+ liquid_phase[MPL::PropertyType::vapour_diffusion]
+ .template value<double>(variables, x_position, t, dt);
+
+ double const f_Tv_D_Tv = f_Tv * D_v * drho_wv_dT;
+ double const D_pv = D_v * drho_wv_dp;
+
+ if (gas_phase &&
+ gas_phase->hasProperty(MPL::PropertyType::heat_capacity))
+ {
+ GlobalDimVectorType const grad_T = dNdx * T;
+ GlobalDimVectorType const grad_p_cap = -dNdx * p_L;
+ // Vapour velocity
+ GlobalDimVectorType const q_v =
+ -(f_Tv_D_Tv * grad_T - D_pv * grad_p_cap) / rho_LR;
+ double const specific_heat_capacity_vapour =
+ gas_phase->property(
+ MaterialPropertyLib::PropertyType::
+ specific_heat_capacity)
+ .template value<double>(variables, x_position, t, dt);
+
+ local_M.template block<temperature_size, temperature_size>(temperature_index,
+ temperature_index)
+ .noalias() +=
+ w *
+ (rho_wv * specific_heat_capacity_vapour * (1 - S_L) * phi) *
+ N.transpose() * N;
+
+ local_K.template block<temperature_size, temperature_size>(temperature_index,
+ temperature_index)
+ .noalias() += N.transpose() * q_v.transpose() * dNdx *
+ rho_wv * specific_heat_capacity_vapour * w;
+ }
+
+ double const storage_coefficient_by_water_vapor =
+ phi * (rho_wv * dS_L_dp_cap + (1 - S_L) * drho_wv_dp);
+ local_M.template block<pressure_size, pressure_size>(pressure_index,
+ pressure_index)
+ .noalias() +=
+ N.transpose() * storage_coefficient_by_water_vapor * N * w;
+
+ double const vapor_expansion_factor = phi * (1 - S_L) * drho_wv_dT;
+ local_M.template block<pressure_size, temperature_size>(pressure_index,
+ temperature_index)
+ .noalias() += N.transpose() * vapor_expansion_factor * N * w;
+
+ local_rhs.template segment<pressure_size>(pressure_index)
+ .noalias() -= f_Tv_D_Tv * dNdx.transpose() * (dNdx * T) * w;
+
+ local_K.template block<pressure_size, pressure_size>(pressure_index,
+ pressure_index)
+ .noalias() += dNdx.transpose() * D_pv * dNdx * w;
+
+ //
+ // Latent heat term
+ //
+ if (liquid_phase.hasProperty(MPL::PropertyType::latent_heat))
+ {
+ double const factor = phi * (1 - S_L) / rho_LR;
+ // The volumetric latent heat of vaporization of liquid water
+ double const L0 =
+ liquid_phase[MPL::PropertyType::latent_heat]
+ .template value<double>(variables, x_position, t, dt) *
+ rho_LR;
+
+ double const drho_LR_dT =
+ liquid_phase[MPL::PropertyType::density]
+ .template dValue<double>(variables,
+ MPL::Variable::temperature,
+ x_position, t, dt);
+
+ double const rho_wv_over_rho_L = rho_wv / rho_LR;
+ local_M.template block<temperature_size, temperature_size>(temperature_index,
+ temperature_index)
+ .noalias() +=
+ factor * L0 *
+ (drho_wv_dT - rho_wv_over_rho_L * drho_LR_dT) *
+ N.transpose() * N * w;
+
+ local_M.template block<temperature_size, pressure_size>(temperature_index,
+ pressure_index)
+ .noalias() +=
+ (factor * L0 *
+ (drho_wv_dp - rho_wv_over_rho_L * drho_LR_dp) +
+ L0 * phi * rho_wv_over_rho_L * dS_L_dp_cap) *
+ N.transpose() * N * w;
+
+ // temperature equation, temperature part
+ local_K.template block<temperature_size, temperature_size>(temperature_index,
+ temperature_index)
+ .noalias() +=
+ L0 * f_Tv_D_Tv * dNdx.transpose() * dNdx * w / rho_LR;
+ // temperature equation, pressure part
+ local_K.template block<temperature_size, pressure_size>(temperature_index,
+ pressure_index)
+ .noalias() +=
+ L0 * D_pv * dNdx.transpose() * dNdx * w / rho_LR;
+ }
+ }
+ }
+
+ if (_process_data.apply_mass_lumping)
+ {
+ auto Mpp = local_M.template block<pressure_size, pressure_size>(
+ pressure_index, pressure_index);
+ Mpp = Mpp.colwise().sum().eval().asDiagonal();
+ }
+
+
+}
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+std::vector<double> const&
+ThermoRichardsFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
+ getIntPtDarcyVelocity(
+ const double /*t*/,
+ std::vector<GlobalVector*> const& /*x*/,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
+ std::vector<double>& cache) const
+{
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ cache.clear();
+ auto cache_matrix = MathLib::createZeroedMatrix<
+ Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
+ cache, GlobalDim, n_integration_points);
+
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ cache_matrix.col(ip).noalias() = _ip_data[ip].v_darcy;
+ }
+
+ return cache;
+}
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+std::vector<double> ThermoRichardsFlowLocalAssembler<
+ ShapeFunction, IntegrationMethod, GlobalDim>::getSaturation() const
+{
+ std::vector<double> result;
+ getIntPtSaturation(0, {}, {}, result);
+ return result;
+}
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+std::vector<double> const&
+ThermoRichardsFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
+ getIntPtSaturation(
+ const double /*t*/,
+ std::vector<GlobalVector*> const& /*x*/,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
+ std::vector<double>& cache) const
+{
+ return ProcessLib::getIntegrationPointScalarData(
+ _ip_data, &IpData::saturation, cache);
+}
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+std::vector<double> ThermoRichardsFlowLocalAssembler<
+ ShapeFunction, IntegrationMethod, GlobalDim>::getPorosity() const
+{
+ std::vector<double> result;
+ getIntPtPorosity(0, {}, {}, result);
+ return result;
+}
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+std::vector<double> const&
+ThermoRichardsFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
+ getIntPtPorosity(
+ const double /*t*/,
+ std::vector<GlobalVector*> const& /*x*/,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
+ std::vector<double>& cache) const
+{
+ return ProcessLib::getIntegrationPointScalarData(_ip_data,
+ &IpData::porosity, cache);
+}
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+std::vector<double> const&
+ThermoRichardsFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
+ getIntPtDryDensitySolid(
+ const double /*t*/,
+ std::vector<GlobalVector*> const& /*x*/,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
+ std::vector<double>& cache) const
+{
+ return ProcessLib::getIntegrationPointScalarData(
+ _ip_data, &IpData::dry_density_solid, cache);
+}
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+void ThermoRichardsFlowLocalAssembler<ShapeFunction, IntegrationMethod,
+ GlobalDim>::
+ computeSecondaryVariableConcrete(double const t, double const dt,
+ Eigen::VectorXd const& local_x,
+ Eigen::VectorXd const& local_x_dot)
+{
+ auto const T =
+ local_x.template segment<temperature_size>(temperature_index);
+
+ auto const T_dot =
+ local_x_dot.template segment<temperature_size>(temperature_index);
+
+ auto const p_L = local_x.template segment<pressure_size>(pressure_index);
+
+ auto p_L_dot = local_x_dot.template segment<pressure_size>(pressure_index);
+
+ auto const& medium = *_process_data.media_map->getMedium(_element.getID());
+ auto const& liquid_phase = medium.phase("AqueousLiquid");
+ auto const& solid_phase = medium.phase("Solid");
+ MPL::VariableArray variables;
+ MPL::VariableArray variables_prev;
+
+ ParameterLib::SpatialPosition x_position;
+ x_position.setElementID(_element.getID());
+
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ double saturation_avg = 0;
+ double porosity_avg = 0;
+
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ x_position.setIntegrationPoint(ip);
+ auto const& N = _ip_data[ip].N;
+
+ double T_ip;
+ NumLib::shapeFunctionInterpolate(T, N, T_ip);
+ double T_dot_ip;
+ NumLib::shapeFunctionInterpolate(T_dot, N, T_dot_ip);
+
+ double p_cap_ip;
+ NumLib::shapeFunctionInterpolate(-p_L, N, p_cap_ip);
+
+ double p_cap_dot_ip;
+ NumLib::shapeFunctionInterpolate(-p_L_dot, N, p_cap_dot_ip);
+
+ variables[static_cast<int>(MPL::Variable::capillary_pressure)] =
+ p_cap_ip;
+ variables[static_cast<int>(MPL::Variable::phase_pressure)] = -p_cap_ip;
+
+ variables[static_cast<int>(MPL::Variable::temperature)] = T_ip;
+
+ auto& S_L = _ip_data[ip].saturation;
+ auto const S_L_prev = _ip_data[ip].saturation_prev;
+ S_L = medium[MPL::PropertyType::saturation]
+ .template value<double>(variables, x_position, t, dt);
+ variables[static_cast<int>(MPL::Variable::liquid_saturation)] = S_L;
+ variables_prev[static_cast<int>(MPL::Variable::liquid_saturation)] =
+ S_L_prev;
+
+ auto chi_S_L = S_L;
+ auto chi_S_L_prev = S_L_prev;
+ if (medium.hasProperty(MPL::PropertyType::bishops_effective_stress))
+ {
+ auto const chi = [&medium, x_position, t, dt](double const S_L) {
+ MPL::VariableArray variables;
+ variables[static_cast<int>(MPL::Variable::liquid_saturation)] =
+ S_L;
+ return medium[MPL::PropertyType::bishops_effective_stress]
+ .template value<double>(variables, x_position, t, dt);
+ };
+ chi_S_L = chi(S_L);
+ chi_S_L_prev = chi(S_L_prev);
+ }
+ variables[static_cast<int>(MPL::Variable::effective_pore_pressure)] =
+ -chi_S_L * p_cap_ip;
+ variables_prev[static_cast<int>(
+ MPL::Variable::effective_pore_pressure)] =
+ -chi_S_L_prev * (p_cap_ip - p_cap_dot_ip * dt);
+
+ auto const alpha =
+ medium[MPL::PropertyType::biot_coefficient]
+ .template value<double>(variables, x_position, t, dt);
+
+ auto& solid_elasticity = *_process_data.simplified_elasticity;
+ auto const beta_S =
+ solid_elasticity.bulkCompressibilityFromYoungsModulus(
+ solid_phase, variables, x_position, t, dt);
+ auto const beta_SR = (1 - alpha) * beta_S;
+ variables[static_cast<int>(MPL::Variable::grain_compressibility)] =
+ beta_SR;
+
+ auto& phi = _ip_data[ip].porosity;
+ { // Porosity update
+ variables_prev[static_cast<int>(MPL::Variable::porosity)] =
+ _ip_data[ip].porosity_prev;
+ phi = medium[MPL::PropertyType::porosity]
+ .template value<double>(variables, variables_prev,
+ x_position, t, dt);
+ variables[static_cast<int>(MPL::Variable::porosity)] = phi;
+ }
+
+ auto const mu =
+ liquid_phase[MPL::PropertyType::viscosity]
+ .template value<double>(variables, x_position, t, dt);
+ auto const rho_LR =
+ liquid_phase[MPL::PropertyType::density]
+ .template value<double>(variables, x_position, t, dt);
+
+ auto const K_intrinsic = MPL::formEigenTensor<GlobalDim>(
+ medium[MPL::PropertyType::permeability]
+ .value(variables, x_position, t, dt));
+
+ double const k_rel =
+ medium[MPL::PropertyType::relative_permeability]
+ .template value<double>(variables, x_position, t, dt);
+
+ GlobalDimMatrixType const K_over_mu = k_rel * K_intrinsic / mu;
+
+ auto const rho_SR =
+ solid_phase[MPL::PropertyType::density]
+ .template value<double>(variables, x_position, t, dt);
+ _ip_data[ip].dry_density_solid = (1 - phi) * rho_SR;
+
+ auto const& b = _process_data.specific_body_force;
+
+ // Compute the velocity
+ auto const& dNdx = _ip_data[ip].dNdx;
+ _ip_data[ip].v_darcy.noalias() =
+ -K_over_mu * dNdx * p_L + rho_LR * K_over_mu * b;
+
+ saturation_avg += S_L;
+ porosity_avg += phi;
+ }
+ saturation_avg /= n_integration_points;
+ porosity_avg /= n_integration_points;
+
+ (*_process_data.element_saturation)[_element.getID()] = saturation_avg;
+ (*_process_data.element_porosity)[_element.getID()] = porosity_avg;
+}
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+unsigned ThermoRichardsFlowLocalAssembler<
+ ShapeFunction, IntegrationMethod, GlobalDim>::getNumberOfIntegrationPoints()
+ const
+{
+ return _integration_method.getNumberOfPoints();
+}
+
+} // namespace ThermoRichardsFlow
+} // namespace ProcessLib
diff --git a/ProcessLib/ThermoRichardsFlow/ThermoRichardsFlowFEM.h b/ProcessLib/ThermoRichardsFlow/ThermoRichardsFlowFEM.h
new file mode 100644
index 00000000000..9330d99849b
--- /dev/null
+++ b/ProcessLib/ThermoRichardsFlow/ThermoRichardsFlowFEM.h
@@ -0,0 +1,172 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2021, 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 "IntegrationPointData.h"
+#include "LocalAssemblerInterface.h"
+#include "MaterialLib/SolidModels/LinearElasticIsotropic.h"
+#include "NumLib/DOF/DOFTableUtil.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
+#include "NumLib/Fem/ShapeMatrixPolicy.h"
+#include "ParameterLib/Parameter.h"
+#include "ProcessLib/Deformation/BMatrixPolicy.h"
+#include "ProcessLib/Deformation/LinearBMatrix.h"
+#include "ProcessLib/LocalAssemblerTraits.h"
+#include "ThermoRichardsFlowProcessData.h"
+
+namespace ProcessLib
+{
+namespace ThermoRichardsFlow
+{
+namespace MPL = MaterialPropertyLib;
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+class ThermoRichardsFlowLocalAssembler : public LocalAssemblerInterface
+{
+public:
+ // Note: temperature variable uses the same shape functions as that are used
+ // by pressure variable.
+ using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
+
+ using GlobalDimMatrixType = typename ShapeMatricesType::GlobalDimMatrixType;
+ using GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType;
+
+ using IpData = IntegrationPointData<ShapeMatricesType>;
+
+ ThermoRichardsFlowLocalAssembler(ThermoRichardsFlowLocalAssembler const&) =
+ delete;
+ ThermoRichardsFlowLocalAssembler(ThermoRichardsFlowLocalAssembler&&) =
+ delete;
+
+ ThermoRichardsFlowLocalAssembler(
+ MeshLib::Element const& e,
+ std::size_t const /*local_matrix_size*/,
+ bool const is_axially_symmetric,
+ unsigned const integration_order,
+ ThermoRichardsFlowProcessData& process_data);
+
+ /// \return the number of read integration points.
+ std::size_t setIPDataInitialConditions(
+ std::string const& name,
+ double const* values,
+ int const integration_order) override;
+
+ void setInitialConditionsConcrete(std::vector<double> const& local_x,
+ double const t,
+ bool const /*use_monolithic_scheme*/,
+ int const /*process_id*/) override;
+
+ void assembleWithJacobian(double const t, double const dt,
+ std::vector<double> const& local_x,
+ std::vector<double> const& local_xdot,
+ const double /*dxdot_dx*/, const double /*dx_dx*/,
+ std::vector<double>& /*local_M_data*/,
+ std::vector<double>& /*local_K_data*/,
+ std::vector<double>& local_rhs_data,
+ std::vector<double>& local_Jac_data) override;
+
+ void assemble(double const t, double const dt,
+ std::vector<double> const& local_x,
+ std::vector<double> const& local_xdot,
+ std::vector<double>& local_M_data,
+ std::vector<double>& local_K_data,
+ std::vector<double>& local_rhs_data) override;
+
+ void initializeConcrete() override
+ {
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ auto& ip_data = _ip_data[ip];
+ ip_data.pushBackState();
+ }
+ }
+
+ void postTimestepConcrete(Eigen::VectorXd const& /*local_x*/,
+ double const /*t*/,
+ double const /*dt*/) override
+ {
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ _ip_data[ip].pushBackState();
+ }
+ }
+
+ void computeSecondaryVariableConcrete(
+ double const t, double const dt, Eigen::VectorXd const& local_x,
+ Eigen::VectorXd const& local_x_dot) override;
+
+ Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
+ const unsigned integration_point) const override
+ {
+ 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());
+ }
+
+ std::vector<double> const& getIntPtDarcyVelocity(
+ const double t,
+ std::vector<GlobalVector*> const& x,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
+ std::vector<double>& cache) const override;
+
+ std::vector<double> getSaturation() const override;
+ std::vector<double> const& getIntPtSaturation(
+ const double t,
+ std::vector<GlobalVector*> const& x,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
+ std::vector<double>& cache) const override;
+
+ std::vector<double> getPorosity() const override;
+ std::vector<double> const& getIntPtPorosity(
+ const double t,
+ std::vector<GlobalVector*> const& x,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
+ std::vector<double>& cache) const override;
+
+ std::vector<double> const& getIntPtDryDensitySolid(
+ const double t,
+ std::vector<GlobalVector*> const& x,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
+ std::vector<double>& cache) const override;
+
+private:
+ unsigned getNumberOfIntegrationPoints() const override;
+
+private:
+ ThermoRichardsFlowProcessData& _process_data;
+
+ std::vector<IpData, Eigen::aligned_allocator<IpData>> _ip_data;
+
+ IntegrationMethod _integration_method;
+ MeshLib::Element const& _element;
+ bool const _is_axially_symmetric;
+
+ static const int temperature_index = 0;
+ static const int temperature_size = ShapeFunction::NPOINTS;
+ static const int pressure_index = temperature_size;
+ static const int pressure_size = ShapeFunction::NPOINTS;
+};
+
+} // namespace ThermoRichardsFlow
+} // namespace ProcessLib
+
+#include "ThermoRichardsFlowFEM-impl.h"
diff --git a/ProcessLib/ThermoRichardsFlow/ThermoRichardsFlowProcess.cpp b/ProcessLib/ThermoRichardsFlow/ThermoRichardsFlowProcess.cpp
new file mode 100644
index 00000000000..7aca7870df0
--- /dev/null
+++ b/ProcessLib/ThermoRichardsFlow/ThermoRichardsFlowProcess.cpp
@@ -0,0 +1,305 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2021, 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 "ThermoRichardsFlowProcess.h"
+
+#include <cassert>
+
+#include "BaseLib/Error.h"
+#include "MeshLib/Elements/Utils.h"
+#include "NumLib/DOF/ComputeSparsityPattern.h"
+#include "ProcessLib/Process.h"
+#include "ProcessLib/Utils/CreateLocalAssemblers.h"
+#include "ThermoRichardsFlowFEM.h"
+
+namespace ProcessLib
+{
+namespace ThermoRichardsFlow
+{
+ThermoRichardsFlowProcess::ThermoRichardsFlowProcess(
+ std::string name,
+ MeshLib::Mesh& mesh,
+ std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+ std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
+ unsigned const integration_order,
+ std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
+ process_variables,
+ ThermoRichardsFlowProcessData&& process_data,
+ SecondaryVariableCollection&& secondary_variables,
+ bool const use_monolithic_scheme)
+ : Process(std::move(name), mesh, std::move(jacobian_assembler), parameters,
+ integration_order, std::move(process_variables),
+ std::move(secondary_variables), use_monolithic_scheme),
+ _process_data(std::move(process_data))
+{
+ _heat_flux = MeshLib::getOrCreateMeshProperty<double>(
+ mesh, "HeatFlux", MeshLib::MeshItemType::Node, 1);
+
+ _hydraulic_flow = MeshLib::getOrCreateMeshProperty<double>(
+ mesh, "HydraulicFlow", MeshLib::MeshItemType::Node, 1);
+
+ // TODO (naumov) remove ip suffix. Probably needs modification of the mesh
+ // properties, s.t. there is no "overlapping" with cell/point data.
+ // See getOrCreateMeshProperty.
+ _integration_point_writer.emplace_back(
+ std::make_unique<IntegrationPointWriter>(
+ "saturation_ip", 1 /*n components*/, integration_order, [this]() {
+ // Result containing integration point data for each local
+ // assembler.
+ std::vector<std::vector<double>> result;
+ result.resize(_local_assemblers.size());
+
+ for (std::size_t i = 0; i < _local_assemblers.size(); ++i)
+ {
+ auto const& local_asm = *_local_assemblers[i];
+ result[i] = local_asm.getSaturation();
+ }
+
+ return result;
+ }));
+
+ _integration_point_writer.emplace_back(
+ std::make_unique<IntegrationPointWriter>(
+ "porosity_ip", 1 /*n components*/, integration_order, [this]() {
+ // Result containing integration point data for each local
+ // assembler.
+ std::vector<std::vector<double>> result;
+ result.resize(_local_assemblers.size());
+
+ for (std::size_t i = 0; i < _local_assemblers.size(); ++i)
+ {
+ auto const& local_asm = *_local_assemblers[i];
+ result[i] = local_asm.getPorosity();
+ }
+
+ return result;
+ }));
+}
+
+void ThermoRichardsFlowProcess::initializeConcreteProcess(
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ MeshLib::Mesh const& mesh,
+ unsigned const integration_order)
+{
+ using nlohmann::json;
+
+ const int process_id = 0;
+ const int variable_id = 0;
+ ProcessLib::createLocalAssemblers<ThermoRichardsFlowLocalAssembler>(
+ mesh.getDimension(), mesh.getElements(), dof_table,
+ getProcessVariables(process_id)[variable_id]
+ .get()
+ .getShapeFunctionOrder(),
+ _local_assemblers, mesh.isAxiallySymmetric(), integration_order,
+ _process_data);
+
+ auto add_secondary_variable = [&](std::string const& name,
+ int const num_components,
+ auto get_ip_values_function) {
+ _secondary_variables.addSecondaryVariable(
+ name,
+ makeExtrapolator(num_components, getExtrapolator(),
+ _local_assemblers,
+ std::move(get_ip_values_function)));
+ };
+
+ add_secondary_variable("velocity", mesh.getDimension(),
+ &LocalAssemblerIF::getIntPtDarcyVelocity);
+
+ add_secondary_variable("saturation", 1,
+ &LocalAssemblerIF::getIntPtSaturation);
+
+ add_secondary_variable("porosity", 1, &LocalAssemblerIF::getIntPtPorosity);
+
+ add_secondary_variable("dry_density_solid", 1,
+ &LocalAssemblerIF::getIntPtDryDensitySolid);
+
+ _process_data.element_saturation = MeshLib::getOrCreateMeshProperty<double>(
+ const_cast<MeshLib::Mesh&>(mesh), "saturation_avg",
+ MeshLib::MeshItemType::Cell, 1);
+
+ _process_data.element_porosity = MeshLib::getOrCreateMeshProperty<double>(
+ const_cast<MeshLib::Mesh&>(mesh), "porosity_avg",
+ MeshLib::MeshItemType::Cell, 1);
+
+ // Set initial conditions for integration point data.
+ for (auto const& ip_writer : _integration_point_writer)
+ {
+ // Find the mesh property with integration point writer's name.
+ auto const& name = ip_writer->name();
+ if (!mesh.getProperties().existsPropertyVector<double>(name))
+ {
+ continue;
+ }
+ auto const& mesh_property =
+ *mesh.getProperties().template getPropertyVector<double>(name);
+
+ // The mesh property must be defined on integration points.
+ if (mesh_property.getMeshItemType() !=
+ MeshLib::MeshItemType::IntegrationPoint)
+ {
+ continue;
+ }
+
+ auto const ip_meta_data = getIntegrationPointMetaData(mesh, name);
+
+ // Check the number of components.
+ if (ip_meta_data.n_components !=
+ mesh_property.getNumberOfGlobalComponents())
+ {
+ OGS_FATAL(
+ "Different number of components in meta data ({:d}) than in "
+ "the integration point field data for '{:s}': {:d}.",
+ ip_meta_data.n_components, name,
+ mesh_property.getNumberOfGlobalComponents());
+ }
+
+ // Now we have a properly named vtk's field data array and the
+ // corresponding meta data.
+ std::size_t position = 0;
+ for (auto& local_asm : _local_assemblers)
+ {
+ std::size_t const integration_points_read =
+ local_asm->setIPDataInitialConditions(
+ name, &mesh_property[position],
+ ip_meta_data.integration_order);
+ if (integration_points_read == 0)
+ {
+ OGS_FATAL(
+ "No integration points read in the integration point "
+ "initial conditions set function.");
+ }
+ position += integration_points_read * ip_meta_data.n_components;
+ }
+ }
+
+ // Initialize local assemblers after all variables have been set.
+ GlobalExecutor::executeMemberOnDereferenced(&LocalAssemblerIF::initialize,
+ _local_assemblers,
+ *_local_to_global_index_map);
+}
+
+void ThermoRichardsFlowProcess::setInitialConditionsConcreteProcess(
+ std::vector<GlobalVector*>& x, double const t, int const process_id)
+{
+ if (process_id != 0)
+ {
+ return;
+ }
+ DBUG("SetInitialConditions ThermoRichardsFlowProcess.");
+
+ GlobalExecutor::executeMemberOnDereferenced(
+ &LocalAssemblerIF::setInitialConditions, _local_assemblers,
+ *_local_to_global_index_map, *x[process_id], t, _use_monolithic_scheme,
+ process_id);
+}
+
+void ThermoRichardsFlowProcess::assembleConcreteProcess(
+ const double t, double const dt,
+ std::vector<GlobalVector*> const& x,
+ std::vector<GlobalVector*> const& xdot, int const process_id,
+ GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b)
+{
+ DBUG("Assemble the equations for ThermoRichardsFlowProcess.");
+
+ std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+ dof_table = {std::ref(*_local_to_global_index_map)};
+ ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
+
+ // Call global assembler for each local assembly item.
+ GlobalExecutor::executeSelectedMemberDereferenced(
+ _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
+ pv.getActiveElementIDs(), dof_table, t, dt, x, xdot, process_id, M, K,
+ b);
+}
+
+void ThermoRichardsFlowProcess::assembleWithJacobianConcreteProcess(
+ const double t, double const dt, std::vector<GlobalVector*> const& x,
+ std::vector<GlobalVector*> const& xdot, const double dxdot_dx,
+ const double dx_dx, int const process_id, GlobalMatrix& M, GlobalMatrix& K,
+ GlobalVector& b, GlobalMatrix& Jac)
+{
+ std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+ dof_tables;
+
+ DBUG(
+ "Assemble the Jacobian of ThermoRichardsFlow for the monolithic "
+ "scheme.");
+ dof_tables.emplace_back(*_local_to_global_index_map);
+
+ ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
+
+ GlobalExecutor::executeSelectedMemberDereferenced(
+ _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
+ _local_assemblers, pv.getActiveElementIDs(), dof_tables, t, dt, x, xdot,
+ dxdot_dx, dx_dx, process_id, M, K, b, Jac);
+
+ auto copyRhs = [&](int const variable_id, auto& output_vector) {
+ transformVariableFromGlobalVector(b, variable_id, dof_tables[0],
+ output_vector, std::negate<double>());
+ };
+
+ copyRhs(0, *_heat_flux);
+ copyRhs(1, *_hydraulic_flow);
+}
+
+void ThermoRichardsFlowProcess::postTimestepConcreteProcess(
+ std::vector<GlobalVector*> const& x, double const t, double const dt,
+ const int process_id)
+{
+ if (process_id != 0)
+ {
+ return;
+ }
+
+ DBUG("PostTimestep ThermoRichardsFlowProcess.");
+
+ auto const dof_tables = getDOFTables(x.size());
+
+ ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
+ GlobalExecutor::executeSelectedMemberOnDereferenced(
+ &LocalAssemblerIF::postTimestep, _local_assemblers,
+ pv.getActiveElementIDs(), dof_tables, x, t, dt);
+}
+
+void ThermoRichardsFlowProcess::computeSecondaryVariableConcrete(
+ const double t, const double dt, std::vector<GlobalVector*> const& x,
+ GlobalVector const& x_dot, int const process_id)
+{
+ if (process_id != 0)
+ {
+ return;
+ }
+ DBUG(
+ "Compute the secondary variables for "
+ "ThermoRichardsFlowProcess.");
+ auto const dof_tables = getDOFTables(x.size());
+ ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
+
+ GlobalExecutor::executeSelectedMemberOnDereferenced(
+ &LocalAssemblerIF::computeSecondaryVariable, _local_assemblers,
+ pv.getActiveElementIDs(), dof_tables, t, dt, x, x_dot, process_id);
+}
+
+
+std::vector<NumLib::LocalToGlobalIndexMap const*>
+ThermoRichardsFlowProcess::getDOFTables(
+ int const number_of_processes) const
+{
+ std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
+ dof_tables.reserve(number_of_processes);
+ std::generate_n(std::back_inserter(dof_tables), number_of_processes,
+ [&]() { return _local_to_global_index_map.get(); });
+ return dof_tables;
+}
+
+
+} // namespace ThermoRichardsFlow
+} // namespace ProcessLib
diff --git a/ProcessLib/ThermoRichardsFlow/ThermoRichardsFlowProcess.h b/ProcessLib/ThermoRichardsFlow/ThermoRichardsFlowProcess.h
new file mode 100644
index 00000000000..ed5b3a48809
--- /dev/null
+++ b/ProcessLib/ThermoRichardsFlow/ThermoRichardsFlowProcess.h
@@ -0,0 +1,156 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2021, 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 "LocalAssemblerInterface.h"
+#include "ProcessLib/Process.h"
+#include "ThermoRichardsFlowProcessData.h"
+
+namespace ProcessLib
+{
+namespace ThermoRichardsFlow
+{
+/**
+ * \brief Global assembler for the monolithic scheme of the non-isothermal
+ * Richards flow.
+ *
+ * <b>Governing equations without vapor diffusion</b>
+ *
+ * The energy balance equation is given by
+ * \f[
+ * (\rho c_p)^{eff}\dot T -
+ * \nabla (\mathbf{k}_T^{eff} \nabla T)+\rho^l c_p^l \nabla T \cdot
+ * \mathbf{v}^l
+ * = Q_T
+ * \f]
+ * with\f$T\f$ the temperature, \f$(\rho c_p)^{eff}\f$ the effective
+ * volumetric heat
+ * capacity, \f$\mathbf{k}_T^{eff} \f$
+ * the effective thermal conductivity, \f$\rho^l\f$ the density of liquid,
+ * \f$c_p^l\f$ the specific heat capacity of liquid, \f$\mathbf{v}^l\f$ the
+ * liquid velocity, and \f$Q_T\f$ the point heat source. The effective
+ * volumetric heat can be considered as a composite of the contributions of
+ * solid phase and the liquid phase as
+ * \f[
+ * (\rho c_p)^{eff} = (1-\phi) \rho^s c_p^s + S^l \phi \rho^l c_p^l
+ * \f]
+ * with \f$\phi\f$ the porosity, \f$S^l\f$ the liquid saturation, \f$\rho^s \f$
+ * the solid density, and \f$c_p^s\f$ the specific heat capacity of solid.
+ * Similarly, the effective thermal conductivity is given by
+ * \f[
+ * \mathbf{k}_T^{eff} = (1-\phi) \mathbf{k}_T^s + S^l \phi k_T^l \mathbf I
+ * \f]
+ * where \f$\mathbf{k}_T^s\f$ is the thermal conductivity tensor of solid, \f$
+ * k_T^l\f$ is the thermal conductivity of liquid, and \f$\mathbf I\f$ is the
+ * identity tensor.
+ *
+ * The mass balance equation is given by
+ * \f{eqnarray*}{
+ * \left(S^l\beta - \phi\frac{\partial S}{\partial p_c}\right) \rho^l\dot p
+ * - S \left( \frac{\partial \rho^l}{\partial T}
+ * +\rho^l(\alpha_B -S)
+ * \alpha_T^s
+ * \right)\dot T\\
+ * +\nabla (\rho^l \mathbf{v}^l) + S \alpha_B \rho^l \nabla \cdot \dot {\mathbf
+ * u}= Q_H
+ * \f}
+ * where \f$p\f$ is the pore pressure, \f$p_c\f$ is the
+ * capillary pressure, which is \f$-p\f$ under the single phase assumption,
+ * \f$\beta\f$ is a composite coefficient by the liquid compressibility and
+ * solid compressibility, \f$\alpha_B\f$ is the Biot's constant,
+ * \f$\alpha_T^s\f$ is the linear thermal expansivity of solid, \f$Q_H\f$
+ * is the point source or sink term, \f$H(S-1)\f$ is the Heaviside function, and
+ * \f$ \mathbf u\f$ is the displacement. While this process does not contain a fully
+ * mechanical coupling, simplfied expressions can be given to approximate the latter
+ * term under certain stress conditions.
+ * The liquid velocity \f$\mathbf{v}^l\f$ is
+ * described by the Darcy's law as
+ * \f[
+ * \mathbf{v}^l=-\frac{{\mathbf k} k_{ref}}{\mu} (\nabla p - \rho^l \mathbf g)
+ * \f]
+ * with \f${\mathbf k}\f$ the intrinsic permeability, \f$k_{ref}\f$ the relative
+ * permeability, \f$\mathbf g\f$ the gravitational force.
+ */
+class ThermoRichardsFlowProcess final : public Process
+{
+public:
+ ThermoRichardsFlowProcess(
+ std::string name,
+ MeshLib::Mesh& mesh,
+ std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&&
+ jacobian_assembler,
+ std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const&
+ parameters,
+ unsigned const integration_order,
+ std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
+ process_variables,
+ ThermoRichardsFlowProcessData&& process_data,
+ SecondaryVariableCollection&& secondary_variables,
+ bool const use_monolithic_scheme);
+
+ //! \name ODESystem interface
+ //! @{
+
+ bool isLinear() const override { return false; }
+ //! @}
+
+private:
+ using LocalAssemblerIF = LocalAssemblerInterface;
+
+ void initializeConcreteProcess(
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ MeshLib::Mesh const& mesh,
+ unsigned const integration_order) override;
+
+ void setInitialConditionsConcreteProcess(std::vector<GlobalVector*>& x,
+ double const t,
+ int const /*process_id*/) override;
+
+ void assembleConcreteProcess(const double t, double const dt,
+ std::vector<GlobalVector*> const& x,
+ std::vector<GlobalVector*> const& xdot,
+ int const process_id, GlobalMatrix& M,
+ GlobalMatrix& K, GlobalVector& b) override;
+
+ void assembleWithJacobianConcreteProcess(
+ const double t, double const dt, std::vector<GlobalVector*> const& x,
+ std::vector<GlobalVector*> const& xdot, const double dxdot_dx,
+ const double dx_dx, int const process_id, GlobalMatrix& M,
+ GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac) override;
+
+ void postTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
+ double const t, double const dt,
+ const int process_id) override;
+
+private:
+ std::vector<MeshLib::Node*> _base_nodes;
+ std::unique_ptr<MeshLib::MeshSubset const> _mesh_subset_base_nodes;
+ ThermoRichardsFlowProcessData _process_data;
+
+ std::vector<std::unique_ptr<LocalAssemblerIF>> _local_assemblers;
+
+ /// Sparsity pattern for the flow equation, and it is initialized only if
+ /// the staggered scheme is used.
+ GlobalSparsityPattern _sparsity_pattern_with_linear_element;
+
+ void computeSecondaryVariableConcrete(double const t, double const dt,
+ std::vector<GlobalVector*> const& x,
+ GlobalVector const& x_dot,
+ int const process_id) override;
+ std::vector<NumLib::LocalToGlobalIndexMap const*> getDOFTables(
+ const int number_of_processes) const;
+
+ MeshLib::PropertyVector<double>* _heat_flux = nullptr;
+ MeshLib::PropertyVector<double>* _hydraulic_flow = nullptr;
+};
+
+} // namespace ThermoRichardsFlow
+} // namespace ProcessLib
diff --git a/ProcessLib/ThermoRichardsFlow/ThermoRichardsFlowProcessData.h b/ProcessLib/ThermoRichardsFlow/ThermoRichardsFlowProcessData.h
new file mode 100644
index 00000000000..a355697363e
--- /dev/null
+++ b/ProcessLib/ThermoRichardsFlow/ThermoRichardsFlowProcessData.h
@@ -0,0 +1,53 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2021, 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 <Eigen/Dense>
+#include <memory>
+
+namespace ProcessLib
+{
+namespace ThermoRichardsFlow
+{
+struct SimplifiedElasticityModel;
+}
+}
+
+namespace MaterialPropertyLib
+{
+class MaterialSpatialDistributionMap;
+}
+
+namespace ProcessLib
+{
+namespace ThermoRichardsFlow
+{
+struct ThermoRichardsFlowProcessData
+{
+ std::unique_ptr<MaterialPropertyLib::MaterialSpatialDistributionMap>
+ media_map = nullptr;
+
+ /// Specific body forces applied to solid and fluid.
+ /// It is usually used to apply gravitational forces.
+ /// A vector of global mesh dimension's length.
+ Eigen::VectorXd const specific_body_force;
+
+ bool const apply_mass_lumping;
+ std::unique_ptr<SimplifiedElasticityModel> simplified_elasticity = nullptr;
+
+ MeshLib::PropertyVector<double>* element_saturation = nullptr;
+ MeshLib::PropertyVector<double>* element_porosity = nullptr;
+
+ EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
+};
+
+} // namespace ThermoRichardsFlow
+} // namespace ProcessLib
diff --git a/scripts/cmake/ProcessesSetup.cmake b/scripts/cmake/ProcessesSetup.cmake
index 7ce7c258314..21f9407d562 100644
--- a/scripts/cmake/ProcessesSetup.cmake
+++ b/scripts/cmake/ProcessesSetup.cmake
@@ -23,6 +23,7 @@ set(_processes_list
ThermoHydroMechanics
ThermoMechanicalPhaseField
ThermoMechanics
+ ThermoRichardsFlow
TwoPhaseFlowWithPP
TwoPhaseFlowWithPrho
)
--
GitLab