From 2bd86efa9ea27ce461a3b333b0dc8045c522c134 Mon Sep 17 00:00:00 2001
From: Christoph Lehmann <christoph.lehmann@ufz.de>
Date: Mon, 30 May 2016 18:12:16 +0200
Subject: [PATCH] [PL] added TES process
---
Applications/ApplicationsLib/ProjectData.cpp | 8 +
ProcessLib/CMakeLists.txt | 3 +
ProcessLib/TES/.gitignore | 1 +
ProcessLib/TES/CMakeLists.txt | 11 +
ProcessLib/TES/TESAssemblyParams.h | 84 ++++
ProcessLib/TES/TESLocalAssembler-impl.h | 224 +++++++++
ProcessLib/TES/TESLocalAssembler.h | 111 +++++
ProcessLib/TES/TESLocalAssemblerData.cpp | 31 ++
ProcessLib/TES/TESLocalAssemblerData.h | 59 +++
ProcessLib/TES/TESLocalAssemblerInner-fwd.h | 33 ++
...LocalAssemblerInner-impl-incl-dynamic.h.in | 6 +
...ESLocalAssemblerInner-impl-incl-fixed.h.in | 7 +
ProcessLib/TES/TESLocalAssemblerInner-impl.h | 372 +++++++++++++++
ProcessLib/TES/TESLocalAssemblerInner.cpp | 35 ++
ProcessLib/TES/TESLocalAssemblerInner.h | 98 ++++
ProcessLib/TES/TESOGS5MaterialModels.cpp | 35 ++
ProcessLib/TES/TESOGS5MaterialModels.h | 411 ++++++++++++++++
ProcessLib/TES/TESProcess.cpp | 451 ++++++++++++++++++
ProcessLib/TES/TESProcess.h | 149 ++++++
ProcessLib/TES/TESReactionAdaptor.cpp | 373 +++++++++++++++
ProcessLib/TES/TESReactionAdaptor.h | 133 ++++++
21 files changed, 2635 insertions(+)
create mode 100644 ProcessLib/TES/.gitignore
create mode 100644 ProcessLib/TES/CMakeLists.txt
create mode 100644 ProcessLib/TES/TESAssemblyParams.h
create mode 100644 ProcessLib/TES/TESLocalAssembler-impl.h
create mode 100644 ProcessLib/TES/TESLocalAssembler.h
create mode 100644 ProcessLib/TES/TESLocalAssemblerData.cpp
create mode 100644 ProcessLib/TES/TESLocalAssemblerData.h
create mode 100644 ProcessLib/TES/TESLocalAssemblerInner-fwd.h
create mode 100644 ProcessLib/TES/TESLocalAssemblerInner-impl-incl-dynamic.h.in
create mode 100644 ProcessLib/TES/TESLocalAssemblerInner-impl-incl-fixed.h.in
create mode 100644 ProcessLib/TES/TESLocalAssemblerInner-impl.h
create mode 100644 ProcessLib/TES/TESLocalAssemblerInner.cpp
create mode 100644 ProcessLib/TES/TESLocalAssemblerInner.h
create mode 100644 ProcessLib/TES/TESOGS5MaterialModels.cpp
create mode 100644 ProcessLib/TES/TESOGS5MaterialModels.h
create mode 100644 ProcessLib/TES/TESProcess.cpp
create mode 100644 ProcessLib/TES/TESProcess.h
create mode 100644 ProcessLib/TES/TESReactionAdaptor.cpp
create mode 100644 ProcessLib/TES/TESReactionAdaptor.h
diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp
index 9d8586a0013..83dca05b08f 100644
--- a/Applications/ApplicationsLib/ProjectData.cpp
+++ b/Applications/ApplicationsLib/ProjectData.cpp
@@ -35,6 +35,7 @@
#include "UncoupledProcessesTimeLoop.h"
#include "ProcessLib/GroundwaterFlow/GroundwaterFlowProcess-fwd.h"
+#include "ProcessLib/TES/TESProcess.h"
namespace detail
@@ -178,6 +179,13 @@ void ProjectData::buildProcesses()
*_mesh_vec[0], *nl_slv, std::move(time_disc),
_process_variables, _parameters, pc));
}
+ else if (type == "TES")
+ {
+ _processes.emplace_back(
+ ProcessLib::TES::createTESProcess<GlobalSetupType>(
+ *_mesh_vec[0], *nl_slv, std::move(time_disc),
+ _process_variables, _parameters, pc));
+ }
else
{
ERR("Unknown process type: %s", type.c_str());
diff --git a/ProcessLib/CMakeLists.txt b/ProcessLib/CMakeLists.txt
index e1738ab9840..9c22324d5e3 100644
--- a/ProcessLib/CMakeLists.txt
+++ b/ProcessLib/CMakeLists.txt
@@ -6,6 +6,9 @@ APPEND_SOURCE_FILES(SOURCES)
add_subdirectory(GroundwaterFlow)
APPEND_SOURCE_FILES(SOURCES GroundwaterFlow)
+add_subdirectory(TES)
+APPEND_SOURCE_FILES(SOURCES TES)
+
add_library(ProcessLib ${SOURCES})
target_link_libraries(ProcessLib
diff --git a/ProcessLib/TES/.gitignore b/ProcessLib/TES/.gitignore
new file mode 100644
index 00000000000..238847fa1e3
--- /dev/null
+++ b/ProcessLib/TES/.gitignore
@@ -0,0 +1 @@
+TESLocalAssemblerInner-impl-incl.h
diff --git a/ProcessLib/TES/CMakeLists.txt b/ProcessLib/TES/CMakeLists.txt
new file mode 100644
index 00000000000..fac50da00c8
--- /dev/null
+++ b/ProcessLib/TES/CMakeLists.txt
@@ -0,0 +1,11 @@
+if(OGS_EIGEN_DYNAMIC_SHAPE_MATRICES)
+ configure_file(
+ "${CMAKE_CURRENT_SOURCE_DIR}/TESLocalAssemblerInner-impl-incl-dynamic.h.in"
+ "${CMAKE_CURRENT_SOURCE_DIR}/TESLocalAssemblerInner-impl-incl.h" COPYONLY
+ )
+else()
+ configure_file(
+ "${CMAKE_CURRENT_SOURCE_DIR}/TESLocalAssemblerInner-impl-incl-fixed.h.in"
+ "${CMAKE_CURRENT_SOURCE_DIR}/TESLocalAssemblerInner-impl-incl.h" COPYONLY
+ )
+endif()
diff --git a/ProcessLib/TES/TESAssemblyParams.h b/ProcessLib/TES/TESAssemblyParams.h
new file mode 100644
index 00000000000..352bc93a2dc
--- /dev/null
+++ b/ProcessLib/TES/TESAssemblyParams.h
@@ -0,0 +1,84 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef PROCESS_LIB_TESPROCESS_NOTPL_H_
+#define PROCESS_LIB_TESPROCESS_NOTPL_H_
+
+#include <Eigen/Eigen>
+#include <Eigen/Sparse>
+
+#include "MaterialsLib/Adsorption/Reaction.h"
+
+#include "ProcessLib/VariableTransformation.h"
+
+namespace ProcessLib
+{
+namespace TES
+{
+const unsigned NODAL_DOF = 3;
+const unsigned COMPONENT_ID_PRESSURE = 0;
+const unsigned COMPONENT_ID_TEMPERATURE = 1;
+const unsigned COMPONENT_ID_MASS_FRACTION = 2;
+
+const double M_N2 = 0.028013;
+const double M_H2O = 0.018016;
+
+struct AssemblyParams
+{
+ Trafo trafo_p{1.0};
+ Trafo trafo_T{1.0};
+ Trafo trafo_x{1.0};
+
+ std::unique_ptr<Adsorption::Reaction> react_sys;
+
+ double fluid_specific_heat_source =
+ std::numeric_limits<double>::quiet_NaN();
+ double cpG = std::numeric_limits<
+ double>::quiet_NaN(); // specific isobaric fluid heat capacity
+
+ Eigen::MatrixXd solid_perm_tensor = Eigen::MatrixXd::Constant(
+ 3, 3, std::numeric_limits<double>::quiet_NaN()); // TODO get dimensions
+ double solid_specific_heat_source =
+ std::numeric_limits<double>::quiet_NaN();
+ double solid_heat_cond = std::numeric_limits<double>::quiet_NaN();
+ double cpS = std::numeric_limits<
+ double>::quiet_NaN(); // specific isobaric solid heat capacity
+
+ double tortuosity = std::numeric_limits<double>::quiet_NaN();
+ double diffusion_coefficient_component =
+ std::numeric_limits<double>::quiet_NaN();
+
+ double poro = std::numeric_limits<double>::quiet_NaN();
+
+ double rho_SR_dry = std::numeric_limits<double>::quiet_NaN();
+
+ const double M_inert = M_N2; // N2
+ const double M_react = M_H2O;
+
+ // TODO unify variable names
+ double initial_solid_density = std::numeric_limits<double>::quiet_NaN();
+
+ double delta_t = std::numeric_limits<double>::quiet_NaN();
+ unsigned iteration_in_current_timestep = 0;
+
+ bool output_element_matrices = false;
+
+ unsigned number_of_try_of_iteration = 0;
+ double current_time = std::numeric_limits<double>::quiet_NaN();
+
+ //! Output global matrix/rhs after first iteration.
+ std::size_t timestep = 0;
+ std::size_t total_iteration = 0;
+};
+
+} // namespace TES
+
+} // namespace ProcessLib
+
+#endif // PROCESS_LIB_TESPROCESS_NOTPL_H_
diff --git a/ProcessLib/TES/TESLocalAssembler-impl.h b/ProcessLib/TES/TESLocalAssembler-impl.h
new file mode 100644
index 00000000000..3b9ad99e5c2
--- /dev/null
+++ b/ProcessLib/TES/TESLocalAssembler-impl.h
@@ -0,0 +1,224 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ */
+
+#ifndef PROCESS_LIB_TES_FEM_IMPL_H_
+#define PROCESS_LIB_TES_FEM_IMPL_H_
+
+#include "MaterialsLib/Adsorption/Adsorption.h"
+#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "NumLib/Fem/ShapeMatrixPolicy.h"
+#include "NumLib/Function/Interpolation.h"
+#include "ProcessLib/Utils/InitShapeMatrices.h"
+
+#include "TESLocalAssembler.h"
+#include "TESReactionAdaptor.h"
+
+namespace
+{
+enum class MatOutType
+{
+ OGS5,
+ PYTHON
+};
+
+const MatOutType MATRIX_OUTPUT_FORMAT = MatOutType::PYTHON;
+
+// TODO move to some location in the OGS core.
+template <typename Mat>
+void ogs5OutMat(const Mat& mat)
+{
+ for (unsigned r = 0; r < mat.rows(); ++r)
+ {
+ switch (MATRIX_OUTPUT_FORMAT)
+ {
+ case MatOutType::OGS5:
+ if (r != 0)
+ std::printf("\n");
+ std::printf("|");
+ break;
+ case MatOutType::PYTHON:
+ if (r != 0)
+ std::printf(",\n");
+ std::printf("[");
+ break;
+ }
+
+ for (unsigned c = 0; c < mat.cols(); ++c)
+ {
+ switch (MATRIX_OUTPUT_FORMAT)
+ {
+ case MatOutType::OGS5:
+ std::printf(" %.16e", mat(r, c));
+ break;
+ case MatOutType::PYTHON:
+ if (c != 0)
+ std::printf(",");
+ std::printf(" %23.16g", mat(r, c));
+ break;
+ }
+ }
+
+ switch (MATRIX_OUTPUT_FORMAT)
+ {
+ case MatOutType::OGS5:
+ std::printf(" | ");
+ break;
+ case MatOutType::PYTHON:
+ std::printf(" ]");
+ break;
+ }
+ }
+ std::printf("\n");
+}
+
+template <typename Vec>
+void ogs5OutVec(const Vec& vec)
+{
+ for (unsigned r = 0; r < vec.size(); ++r)
+ {
+ switch (MATRIX_OUTPUT_FORMAT)
+ {
+ case MatOutType::OGS5:
+ if (r != 0)
+ std::printf("\n");
+ std::printf("| %.16e | ", vec[r]);
+ break;
+ case MatOutType::PYTHON:
+ if (r != 0)
+ std::printf(",\n");
+ std::printf("[ %23.16g ]", vec[r]);
+ break;
+ }
+ }
+ std::printf("\n");
+}
+
+} // anonymous namespace
+
+namespace ProcessLib
+{
+namespace TES
+{
+template <typename ShapeFunction_, typename IntegrationMethod_,
+ typename GlobalMatrix, typename GlobalVector, unsigned GlobalDim>
+TESLocalAssembler<
+ ShapeFunction_, IntegrationMethod_, GlobalMatrix, GlobalVector,
+ GlobalDim>::TESLocalAssembler(MeshLib::Element const& e,
+ std::size_t const /*local_matrix_size*/,
+ unsigned const integration_order,
+ AssemblyParams const& asm_params)
+ : _shape_matrices(initShapeMatrices<ShapeFunction, ShapeMatricesType,
+ IntegrationMethod_, GlobalDim>(
+ e, integration_order)),
+ _d(asm_params,
+ // TODO narrowing conversion
+ static_cast<const unsigned>(
+ _shape_matrices.front()
+ .N.rows()) /* number of integration points */,
+ GlobalDim),
+ _local_M(ShapeFunction::NPOINTS * NODAL_DOF,
+ ShapeFunction::NPOINTS * NODAL_DOF),
+ _local_K(ShapeFunction::NPOINTS * NODAL_DOF,
+ ShapeFunction::NPOINTS * NODAL_DOF),
+ _local_b(ShapeFunction::NPOINTS * NODAL_DOF),
+ _integration_order(integration_order)
+{
+}
+
+template <typename ShapeFunction_, typename IntegrationMethod_,
+ typename GlobalMatrix, typename GlobalVector, unsigned GlobalDim>
+void TESLocalAssembler<ShapeFunction_, IntegrationMethod_, GlobalMatrix,
+ GlobalVector,
+ GlobalDim>::assemble(const double /*t*/,
+ std::vector<double> const& local_x)
+{
+ _local_M.setZero();
+ _local_K.setZero();
+ _local_b.setZero();
+
+ IntegrationMethod_ integration_method(_integration_order);
+ unsigned const n_integration_points = integration_method.getNPoints();
+
+ _d.preEachAssemble();
+
+ for (std::size_t ip(0); ip < n_integration_points; ip++)
+ {
+ auto const& sm = _shape_matrices[ip];
+ auto const& wp = integration_method.getWeightedPoint(ip);
+ auto const weight = wp.getWeight();
+
+ _d.assembleIntegrationPoint(ip, local_x, sm.N, sm.dNdx, sm.J, sm.detJ,
+ weight, _local_M, _local_K, _local_b);
+ }
+
+ if (_d.getAssemblyParameters().output_element_matrices)
+ {
+ std::puts("### Element: ?");
+
+ std::puts("---Velocity of water");
+ for (auto const& vs : _d.getData().velocity)
+ {
+ std::printf("| ");
+ for (auto v : vs)
+ {
+ std::printf("%23.16e ", v);
+ }
+ std::printf("|\n");
+ }
+
+ std::printf("\n---Mass matrix: \n");
+ ogs5OutMat(_local_M);
+ std::printf("\n");
+
+ std::printf("---Laplacian + Advective + Content matrix: \n");
+ ogs5OutMat(_local_K);
+ std::printf("\n");
+
+ std::printf("---RHS: \n");
+ ogs5OutVec(_local_b);
+ std::printf("\n");
+ }
+}
+
+template <typename ShapeFunction_, typename IntegrationMethod_,
+ typename GlobalMatrix, typename GlobalVector, unsigned GlobalDim>
+void TESLocalAssembler<ShapeFunction_, IntegrationMethod_, GlobalMatrix,
+ GlobalVector, GlobalDim>::
+ addToGlobal(
+ AssemblerLib::LocalToGlobalIndexMap::RowColumnIndices const& indices,
+ GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b) const
+{
+ M.add(indices, _local_M);
+ K.add(indices, _local_K);
+ b.add(indices.rows, _local_b);
+}
+
+template <typename ShapeFunction_, typename IntegrationMethod_,
+ typename GlobalMatrix, typename GlobalVector, unsigned GlobalDim>
+std::vector<double> const& TESLocalAssembler<
+ ShapeFunction_, IntegrationMethod_, GlobalMatrix, GlobalVector,
+ GlobalDim>::getIntegrationPointValues(TESIntPtVariables const var,
+ std::vector<double>& cache) const
+{
+ return _d.getIntegrationPointValues(var, cache);
+}
+
+template <typename ShapeFunction_, typename IntegrationMethod_,
+ typename GlobalMatrix, typename GlobalVector, unsigned GlobalDim>
+bool TESLocalAssembler<
+ ShapeFunction_, IntegrationMethod_, GlobalMatrix, GlobalVector,
+ GlobalDim>::checkBounds(std::vector<double> const& local_x,
+ std::vector<double> const& local_x_prev_ts)
+{
+ return _d.getReactionAdaptor().checkBounds(local_x, local_x_prev_ts);
+}
+
+} // namespace TES
+} // namespace ProcessLib
+
+#endif // PROCESS_LIB_TES_FEM_IMPL_H_
diff --git a/ProcessLib/TES/TESLocalAssembler.h b/ProcessLib/TES/TESLocalAssembler.h
new file mode 100644
index 00000000000..68da0989069
--- /dev/null
+++ b/ProcessLib/TES/TESLocalAssembler.h
@@ -0,0 +1,111 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef PROCESS_LIB_TES_FEM_H_
+#define PROCESS_LIB_TES_FEM_H_
+
+#include <memory>
+#include <vector>
+
+#include "TESAssemblyParams.h"
+#include "TESLocalAssemblerInner-fwd.h"
+
+#include "NumLib/Extrapolation/Extrapolator.h"
+
+namespace ProcessLib
+{
+namespace TES
+{
+template <typename GlobalMatrix, typename GlobalVector>
+class TESLocalAssemblerInterface
+ : public NumLib::Extrapolatable<GlobalVector, TESIntPtVariables>
+{
+public:
+ virtual ~TESLocalAssemblerInterface() = default;
+
+ virtual void assemble(double const t,
+ std::vector<double> const& local_x) = 0;
+
+ virtual void addToGlobal(
+ AssemblerLib::LocalToGlobalIndexMap::RowColumnIndices const&,
+ GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b) const = 0;
+
+ virtual bool checkBounds(std::vector<double> const& local_x,
+ std::vector<double> const& local_x_prev_ts) = 0;
+};
+
+template <typename ShapeFunction_, typename IntegrationMethod_,
+ typename GlobalMatrix, typename GlobalVector, unsigned GlobalDim>
+class TESLocalAssembler final
+ : public TESLocalAssemblerInterface<GlobalMatrix, GlobalVector>
+{
+public:
+ using ShapeFunction = ShapeFunction_;
+ using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
+ using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices;
+
+ TESLocalAssembler(MeshLib::Element const& e,
+ std::size_t const local_matrix_size,
+ unsigned const integration_order,
+ AssemblyParams const& asm_params);
+
+ void assemble(double const t, std::vector<double> const& local_x) override;
+
+ void addToGlobal(
+ AssemblerLib::LocalToGlobalIndexMap::RowColumnIndices const& indices,
+ GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b) const override;
+
+ Eigen::Map<const Eigen::VectorXd> getShapeMatrix(
+ const unsigned integration_point) const override
+ {
+ auto const& N = _shape_matrices[integration_point].N;
+
+ // assumes N is stored contiguously in memory
+ return Eigen::Map<const Eigen::VectorXd>(N.data(), N.size());
+ }
+
+ bool checkBounds(std::vector<double> const& local_x,
+ std::vector<double> const& local_x_prev_ts) override;
+
+ std::vector<double> const& getIntegrationPointValues(
+ TESIntPtVariables const var, std::vector<double>& cache) const override;
+
+private:
+ std::vector<ShapeMatrices> _shape_matrices;
+
+ using LAT = LocalAssemblerTraits<ShapeMatricesType, ShapeFunction::NPOINTS,
+ NODAL_DOF, GlobalDim>;
+
+ TESLocalAssemblerInner<LAT> _d;
+
+ using NodalMatrixType = typename LAT::LocalMatrix;
+ using NodalVectorType = typename LAT::LocalVector;
+
+ static_assert(
+ std::is_same<NodalMatrixType, typename LAT::LocalMatrix>::value,
+ "local matrix and data traits matrix do not coincide");
+ static_assert(
+ std::is_same<NodalVectorType, typename LAT::LocalVector>::value,
+ "local vector and data traits vector do not coincide");
+
+ // TODO Change VectorMatrixAssembler s.t. these can be omitted.
+ NodalMatrixType _local_M;
+ NodalMatrixType _local_K;
+ NodalVectorType _local_b;
+
+ // TODO Use the value from Process
+ unsigned const _integration_order;
+};
+
+} // namespace TES
+} // namespace ProcessLib
+
+#include "TESLocalAssembler-impl.h"
+
+#endif // PROCESS_LIB_TES_FEM_H_
diff --git a/ProcessLib/TES/TESLocalAssemblerData.cpp b/ProcessLib/TES/TESLocalAssemblerData.cpp
new file mode 100644
index 00000000000..9a1ae0d91fb
--- /dev/null
+++ b/ProcessLib/TES/TESLocalAssemblerData.cpp
@@ -0,0 +1,31 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, 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 "TESLocalAssemblerData.h"
+#include "TESReactionAdaptor.h"
+
+namespace ProcessLib
+{
+namespace TES
+{
+TESLocalAssemblerData::TESLocalAssemblerData(AssemblyParams const& ap_,
+ const unsigned num_int_pts,
+ const unsigned dimension)
+ : ap(ap_),
+ solid_density(num_int_pts, ap_.initial_solid_density),
+ reaction_rate(num_int_pts),
+ velocity(dimension, std::vector<double>(num_int_pts)),
+ reaction_adaptor(TESFEMReactionAdaptor::newInstance(*this)),
+ solid_density_prev_ts(num_int_pts, ap_.initial_solid_density),
+ reaction_rate_prev_ts(num_int_pts)
+{
+}
+
+TESLocalAssemblerData::~TESLocalAssemblerData() = default;
+}
+} // namespaces
diff --git a/ProcessLib/TES/TESLocalAssemblerData.h b/ProcessLib/TES/TESLocalAssemblerData.h
new file mode 100644
index 00000000000..cd86950c22b
--- /dev/null
+++ b/ProcessLib/TES/TESLocalAssemblerData.h
@@ -0,0 +1,59 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef PROCESSLIB_TES_TESLOCALASSEMBLERDATA_H
+#define PROCESSLIB_TES_TESLOCALASSEMBLERDATA_H
+
+#include "TESAssemblyParams.h"
+
+namespace ProcessLib
+{
+namespace TES
+{
+class TESFEMReactionAdaptor;
+
+struct TESLocalAssemblerData
+{
+ TESLocalAssemblerData(AssemblyParams const& ap_, const unsigned num_int_pts,
+ const unsigned dimension);
+
+ ~TESLocalAssemblerData();
+
+ AssemblyParams const& ap;
+
+ // integration point quantities
+ std::vector<double> solid_density;
+ std::vector<double> reaction_rate; // dC/dt * rho_SR_dry
+ std::vector<std::vector<double>>
+ velocity; // vector of velocities for each integration point
+
+ // integration point values of unknowns -- temporary storage
+ double p = std::numeric_limits<double>::quiet_NaN(); // gas pressure
+ double T = std::numeric_limits<double>::quiet_NaN(); // temperature
+ double vapour_mass_fraction = std::numeric_limits<
+ double>::quiet_NaN(); // fluid mass fraction of the second component
+
+ // temporary storage for some properties
+ // values do not change during the assembly of one integration point
+ double rho_GR = std::numeric_limits<double>::quiet_NaN();
+ double p_V =
+ std::numeric_limits<double>::quiet_NaN(); // vapour partial pressure
+ double qR = std::numeric_limits<
+ double>::quiet_NaN(); // reaction rate, use this in assembly!!!
+
+ std::unique_ptr<TESFEMReactionAdaptor> const reaction_adaptor;
+
+ // variables at previous timestep
+ std::vector<double> solid_density_prev_ts;
+ std::vector<double> reaction_rate_prev_ts; // could also be calculated from
+ // solid_density_prev_ts
+};
+}
+} // namespaces
+#endif // PROCESSLIB_TES_TESLOCALASSEMBLERDATA_H
diff --git a/ProcessLib/TES/TESLocalAssemblerInner-fwd.h b/ProcessLib/TES/TESLocalAssemblerInner-fwd.h
new file mode 100644
index 00000000000..b15365a1520
--- /dev/null
+++ b/ProcessLib/TES/TESLocalAssemblerInner-fwd.h
@@ -0,0 +1,33 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef PROCESS_LIB_TESFEM_DATA_FWD_H_
+#define PROCESS_LIB_TESFEM_DATA_FWD_H_
+
+#include "TESLocalAssemblerInner.h"
+
+namespace ProcessLib
+{
+namespace TES
+{
+#ifdef OGS_EIGEN_DYNAMIC_SHAPE_MATRICES
+
+extern template class TESLocalAssemblerInner<
+ LocalAssemblerTraits<ShapeMatrixPolicyType<void, 0>, 0, 0, 0>>;
+
+static_assert(OGS_EIGEN_DYNAMIC_SHAPE_MATRICES_FLAG == 1,
+ "inconsistent use of macros");
+#else
+static_assert(OGS_EIGEN_DYNAMIC_SHAPE_MATRICES_FLAG == 0,
+ "inconsistent use of macros");
+#endif
+}
+}
+
+#endif // PROCESS_LIB_TESFEM_DATA_FWD_H_
diff --git a/ProcessLib/TES/TESLocalAssemblerInner-impl-incl-dynamic.h.in b/ProcessLib/TES/TESLocalAssemblerInner-impl-incl-dynamic.h.in
new file mode 100644
index 00000000000..a905c565945
--- /dev/null
+++ b/ProcessLib/TES/TESLocalAssemblerInner-impl-incl-dynamic.h.in
@@ -0,0 +1,6 @@
+#ifndef PROCESSLIB_TES_FEM_DATA_IMPL_INCL
+#define PROCESSLIB_TES_FEM_DATA_IMPL_INCL
+
+static_assert(OGS_EIGEN_DYNAMIC_SHAPE_MATRICES_FLAG == 1, "inconsistent use of macros");
+
+#endif // PROCESSLIB_TES_FEM_DATA_IMPL_INCL
diff --git a/ProcessLib/TES/TESLocalAssemblerInner-impl-incl-fixed.h.in b/ProcessLib/TES/TESLocalAssemblerInner-impl-incl-fixed.h.in
new file mode 100644
index 00000000000..4da6ccdcae2
--- /dev/null
+++ b/ProcessLib/TES/TESLocalAssemblerInner-impl-incl-fixed.h.in
@@ -0,0 +1,7 @@
+#ifndef PROCESSLIB_TES_FEM_DATA_IMPL_INCL
+#define PROCESSLIB_TES_FEM_DATA_IMPL_INCL
+
+static_assert(OGS_EIGEN_DYNAMIC_SHAPE_MATRICES_FLAG == 0, "inconsistent use of macros");
+#include "TESLocalAssemblerInner-impl.h"
+
+#endif // PROCESSLIB_TES_FEM_DATA_IMPL_INCL
diff --git a/ProcessLib/TES/TESLocalAssemblerInner-impl.h b/ProcessLib/TES/TESLocalAssemblerInner-impl.h
new file mode 100644
index 00000000000..f58e1cb4c5d
--- /dev/null
+++ b/ProcessLib/TES/TESLocalAssemblerInner-impl.h
@@ -0,0 +1,372 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * The code of this file is used to decouple the evaluation of matrix elements
+ * from the rest of OGS6,
+ * not all of OGS6 has to be recompiled every time a small change is done.
+ */
+
+#ifndef PROCESS_LIB_TESFEM_DATA_IMPL_H_
+#define PROCESS_LIB_TESFEM_DATA_IMPL_H_
+
+#include <cstdio>
+#include <iostream>
+
+#include <logog/include/logog.hpp>
+
+#include "NumLib/Function/Interpolation.h"
+
+#include "TESLocalAssemblerInner-fwd.h"
+#include "TESOGS5MaterialModels.h"
+#include "TESReactionAdaptor.h"
+
+namespace ProcessLib
+{
+namespace TES
+{
+template <typename Traits>
+TESLocalAssemblerInner<Traits>::TESLocalAssemblerInner(
+ const AssemblyParams& ap, const unsigned num_int_pts,
+ const unsigned dimension)
+ : _d{ap, num_int_pts, dimension}
+{
+}
+
+template <typename Traits>
+Eigen::Matrix3d TESLocalAssemblerInner<Traits>::getMassCoeffMatrix(
+ const unsigned int_pt)
+{
+ // TODO: Dalton's law property
+ const double dxn_dxm = Adsorption::AdsorptionReaction::dMolarFraction(
+ _d.vapour_mass_fraction, _d.ap.M_react, _d.ap.M_inert);
+
+ const double M_pp = _d.ap.poro / _d.p * _d.rho_GR;
+ const double M_pT = -_d.ap.poro / _d.T * _d.rho_GR;
+ const double M_px = (_d.ap.M_react - _d.ap.M_inert) * _d.p /
+ (GAS_CONST * _d.T) * dxn_dxm * _d.ap.poro;
+
+ const double M_Tp = -_d.ap.poro;
+ const double M_TT =
+ _d.ap.poro * _d.rho_GR * _d.ap.cpG // TODO: vapour heat capacity
+ +
+ (1.0 - _d.ap.poro) * _d.solid_density[int_pt] *
+ _d.ap.cpS; // TODO: adsorbate heat capacity
+ const double M_Tx = 0.0;
+
+ const double M_xp = 0.0;
+ const double M_xT = 0.0;
+ const double M_xx = _d.ap.poro * _d.rho_GR;
+
+ Eigen::Matrix3d M;
+ M << M_pp, M_pT, M_px, M_Tp, M_TT, M_Tx, M_xp, M_xT, M_xx;
+
+ return M;
+}
+
+template <typename Traits>
+typename Traits::LaplaceMatrix
+TESLocalAssemblerInner<Traits>::getLaplaceCoeffMatrix(const unsigned /*int_pt*/,
+ const unsigned dim)
+{
+ const double eta_GR = fluid_viscosity(_d.p, _d.T, _d.vapour_mass_fraction);
+
+ const double lambda_F =
+ fluid_heat_conductivity(_d.p, _d.T, _d.vapour_mass_fraction);
+ const double lambda_S = _d.ap.solid_heat_cond;
+
+ using Mat = typename Traits::MatrixDimDim;
+
+ typename Traits::LaplaceMatrix L =
+ Traits::LaplaceMatrix::Zero(dim * NODAL_DOF, dim * NODAL_DOF);
+
+ // TODO: k_rel
+ // L_pp
+ Traits::blockDimDim(L, 0, 0, dim, dim) =
+ Traits::blockDimDim(_d.ap.solid_perm_tensor, 0, 0, dim, dim) *
+ _d.rho_GR / eta_GR;
+
+ // TODO: add zeolite part
+ // L_TT
+ Traits::blockDimDim(L, dim, dim, dim, dim) =
+ Mat::Identity(dim, dim) *
+ (_d.ap.poro * lambda_F + (1.0 - _d.ap.poro) * lambda_S);
+
+ // L_xx
+ Traits::blockDimDim(L, 2 * dim, 2 * dim, dim, dim) =
+ Mat::Identity(dim, dim) * (_d.ap.tortuosity * _d.ap.poro * _d.rho_GR *
+ _d.ap.diffusion_coefficient_component);
+
+ return L;
+}
+
+template <typename Traits>
+Eigen::Matrix3d TESLocalAssemblerInner<Traits>::getAdvectionCoeffMatrix(
+ const unsigned /*int_pt*/)
+{
+ const double A_pp = 0.0;
+ const double A_pT = 0.0;
+
+ const double A_px = 0.0;
+
+ const double A_Tp = 0.0;
+
+ const double A_TT = _d.rho_GR * _d.ap.cpG; // porosity?
+ const double A_Tx = 0.0;
+
+ const double A_xp = 0.0;
+ const double A_xT = 0.0;
+ const double A_xx = _d.rho_GR; // porosity?
+
+ Eigen::Matrix3d A;
+ A << A_pp, A_pT, A_px, A_Tp, A_TT, A_Tx, A_xp, A_xT, A_xx;
+
+ return A;
+}
+
+template <typename Traits>
+Eigen::Matrix3d TESLocalAssemblerInner<Traits>::getContentCoeffMatrix(
+ const unsigned /*int_pt*/)
+{
+ const double C_pp = 0.0;
+ const double C_pT = 0.0;
+
+ const double C_px = 0.0;
+
+ const double C_Tp = 0.0;
+
+ const double C_TT = 0.0;
+ const double C_Tx = 0.0;
+
+ const double C_xp = 0.0;
+ const double C_xT = 0.0;
+ const double C_xx = (_d.ap.poro - 1.0) * _d.qR;
+
+ Eigen::Matrix3d C;
+ C << C_pp, C_pT, C_px, C_Tp, C_TT, C_Tx, C_xp, C_xT, C_xx;
+
+ return C;
+}
+
+template <typename Traits>
+Eigen::Vector3d TESLocalAssemblerInner<Traits>::getRHSCoeffVector(
+ const unsigned int_pt)
+{
+ const double reaction_enthalpy =
+ _d.ap.react_sys->getEnthalpy(_d.p_V, _d.T, _d.ap.M_react);
+
+ const double rhs_p =
+ (_d.ap.poro - 1.0) * _d.qR; // TODO [CL] body force term
+
+ const double rhs_T =
+ _d.rho_GR * _d.ap.poro * _d.ap.fluid_specific_heat_source +
+ (1.0 - _d.ap.poro) * _d.qR * reaction_enthalpy +
+ _d.solid_density[int_pt] * (1.0 - _d.ap.poro) *
+ _d.ap.solid_specific_heat_source;
+ // TODO [CL] momentum production term
+
+ const double rhs_x =
+ (_d.ap.poro - 1.0) * _d.qR; // TODO [CL] what if x < 0.0
+
+ Eigen::Vector3d rhs;
+ rhs << rhs_p, rhs_T, rhs_x;
+
+ return rhs;
+}
+
+template <typename Traits>
+void TESLocalAssemblerInner<Traits>::initReaction(const unsigned int_pt)
+{
+ auto const& rate = _d.reaction_adaptor->initReaction(int_pt);
+ _d.qR = rate.reaction_rate;
+ _d.reaction_rate[int_pt] = rate.reaction_rate;
+ _d.solid_density[int_pt] = rate.solid_density;
+}
+
+template <typename Traits>
+void TESLocalAssemblerInner<Traits>::preEachAssembleIntegrationPoint(
+ const unsigned int_pt,
+ const std::vector<double>& localX,
+ typename Traits::ShapeMatrices::ShapeType const& smN,
+ typename Traits::ShapeMatrices::DxShapeType const& /*smDNdx*/,
+ typename Traits::ShapeMatrices::JacobianType const& /*smJ*/,
+ const double /*smDetJ*/)
+{
+#ifndef NDEBUG
+ // fill local data with garbage to aid in debugging
+ _d.p = _d.T = _d.vapour_mass_fraction = _d.p_V = _d.rho_GR = _d.qR =
+ std::numeric_limits<double>::quiet_NaN();
+#endif
+
+ NumLib::shapeFunctionInterpolate(localX, smN, _d.p, _d.T,
+ _d.vapour_mass_fraction);
+
+ // pre-compute certain properties
+ _d.p_V = _d.p * Adsorption::AdsorptionReaction::getMolarFraction(
+ _d.vapour_mass_fraction, _d.ap.M_react, _d.ap.M_inert);
+
+ initReaction(int_pt);
+
+ assert(_d.p > 0.0);
+ assert(_d.T > 0.0);
+ assert(0.0 <= _d.vapour_mass_fraction && _d.vapour_mass_fraction <= 1.0);
+
+ _d.rho_GR = fluid_density(_d.p, _d.T, _d.vapour_mass_fraction);
+}
+
+template <typename Traits>
+std::vector<double> const&
+TESLocalAssemblerInner<Traits>::getIntegrationPointValues(
+ TESIntPtVariables const var, std::vector<double>& cache) const
+{
+ switch (var)
+ {
+ case TESIntPtVariables::REACTION_RATE:
+ return _d.reaction_rate;
+ case TESIntPtVariables::SOLID_DENSITY:
+ return _d.solid_density;
+ case TESIntPtVariables::VELOCITY_X:
+ return _d.velocity[0];
+ case TESIntPtVariables::VELOCITY_Y:
+ assert(_d.velocity.size() >= 2);
+ return _d.velocity[1];
+ case TESIntPtVariables::VELOCITY_Z:
+ assert(_d.velocity.size() >= 3);
+ return _d.velocity[2];
+
+ case TESIntPtVariables::LOADING:
+ {
+ auto& Cs = cache;
+ Cs.clear();
+ Cs.reserve(_d.solid_density.size());
+
+ for (auto rho_SR : _d.solid_density)
+ {
+ Cs.push_back(rho_SR / _d.ap.rho_SR_dry - 1.0);
+ }
+
+ return Cs;
+ }
+
+ // TODO that's an element value, ain't it?
+ case TESIntPtVariables::REACTION_DAMPING_FACTOR:
+ {
+ auto const num_integration_points = _d.solid_density.size();
+ auto& alphas = cache;
+ alphas.clear();
+ alphas.resize(num_integration_points,
+ _d.reaction_adaptor->getReactionDampingFactor());
+
+ return alphas;
+ }
+ }
+
+ cache.clear();
+ return cache;
+}
+
+template <typename Traits>
+void TESLocalAssemblerInner<Traits>::assembleIntegrationPoint(
+ unsigned integration_point,
+ std::vector<double> const& localX,
+ const typename Traits::ShapeMatrices::ShapeType& smN,
+ const typename Traits::ShapeMatrices::DxShapeType& smDNdx,
+ const typename Traits::ShapeMatrices::JacobianType& smJ,
+ const double smDetJ,
+ const double weight,
+ typename Traits::LocalMatrix& local_M,
+ typename Traits::LocalMatrix& local_K,
+ typename Traits::LocalVector& local_b)
+{
+ preEachAssembleIntegrationPoint(integration_point, localX, smN, smDNdx, smJ,
+ smDetJ);
+
+ auto const N = smDNdx.cols(); // number of integration points
+ auto const D = smDNdx.rows(); // global dimension: 1, 2 or 3
+
+ assert(N * NODAL_DOF == local_M.cols());
+
+ auto const laplaceCoeffMat = getLaplaceCoeffMatrix(integration_point, D);
+ assert(laplaceCoeffMat.cols() == D * NODAL_DOF);
+ auto const massCoeffMat = getMassCoeffMatrix(integration_point);
+ auto const advCoeffMat = getAdvectionCoeffMatrix(integration_point);
+ auto const contentCoeffMat = getContentCoeffMatrix(integration_point);
+
+ // calculate velocity
+ assert((unsigned)smDNdx.rows() == _d.velocity.size() &&
+ (unsigned)smDNdx.cols() == _d.velocity[0].size());
+
+ auto const velocity =
+ (Traits::blockDimDim(laplaceCoeffMat, 0, 0, D, D) *
+ (smDNdx * Eigen::Map<const typename Traits::Vector1Comp>(localX.data(),
+ N) // grad_p
+ /
+ -_d.rho_GR))
+ .eval();
+ assert(velocity.size() == D);
+
+ for (unsigned d = 0; d < D; ++d)
+ {
+ _d.velocity[d][integration_point] = velocity[d];
+ }
+
+ auto const detJ_w_N = (smDetJ * weight * smN).eval();
+ auto const detJ_w_N_NT = (detJ_w_N * smN.transpose()).eval();
+ assert(detJ_w_N_NT.rows() == N && detJ_w_N_NT.cols() == N);
+
+ auto const detJ_w_N_vT_dNdx =
+ (detJ_w_N * velocity.transpose() * smDNdx).eval();
+ assert(detJ_w_N_vT_dNdx.rows() == N && detJ_w_N_vT_dNdx.cols() == N);
+
+ for (unsigned r = 0; r < NODAL_DOF; ++r)
+ {
+ for (unsigned c = 0; c < NODAL_DOF; ++c)
+ {
+ Traits::blockShpShp(local_K, N * r, N * c, N, N).noalias() +=
+ smDetJ * weight * smDNdx.transpose() *
+ Traits::blockDimDim(laplaceCoeffMat, D * r, D * c, D, D) *
+ smDNdx // end Laplacian part
+ + detJ_w_N_NT * contentCoeffMat(r, c) +
+ detJ_w_N_vT_dNdx * advCoeffMat(r, c);
+ Traits::blockShpShp(local_M, N * r, N * c, N, N).noalias() +=
+ detJ_w_N_NT * massCoeffMat(r, c);
+ }
+ }
+
+ auto const rhsCoeffVector = getRHSCoeffVector(integration_point);
+
+ for (unsigned r = 0; r < NODAL_DOF; ++r)
+ {
+ Traits::blockShp(local_b, N * r, N).noalias() +=
+ rhsCoeffVector(r) * smN * smDetJ * weight;
+ }
+}
+
+template <typename Traits>
+void TESLocalAssemblerInner<Traits>::preEachAssemble()
+{
+ if (_d.ap.iteration_in_current_timestep == 1)
+ {
+ if (_d.ap.number_of_try_of_iteration ==
+ 1) // TODO has to hold if the above holds.
+ {
+ _d.solid_density_prev_ts = _d.solid_density;
+ _d.reaction_rate_prev_ts = _d.reaction_rate;
+
+ _d.reaction_adaptor->preZerothTryAssemble();
+ }
+ else
+ {
+ _d.solid_density = _d.solid_density_prev_ts;
+ }
+ }
+}
+
+} // namespace TES
+
+} // namespace ProcessLib
+
+#endif // PROCESS_LIB_TESFEM_DATA_IMPL_H_
diff --git a/ProcessLib/TES/TESLocalAssemblerInner.cpp b/ProcessLib/TES/TESLocalAssemblerInner.cpp
new file mode 100644
index 00000000000..12893b6d2fa
--- /dev/null
+++ b/ProcessLib/TES/TESLocalAssemblerInner.cpp
@@ -0,0 +1,35 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * The code of this file is used to decouple the evaluation of matrix elements
+ * from the rest of OGS6,
+ * not all of OGS6 has to be recompiled every time a small change is done.
+ */
+
+#include "TESLocalAssemblerInner-fwd.h"
+
+#ifdef OGS_EIGEN_DYNAMIC_SHAPE_MATRICES
+#include "TESLocalAssemblerInner-impl.h"
+#endif
+
+namespace ProcessLib
+{
+namespace TES
+{
+#ifdef OGS_EIGEN_DYNAMIC_SHAPE_MATRICES
+
+template class TESLocalAssemblerInner<
+ LocalAssemblerTraits<ShapeMatrixPolicyType<void, 0>, 0, 0, 0>>;
+
+static_assert(OGS_EIGEN_DYNAMIC_SHAPE_MATRICES_FLAG == 1,
+ "inconsistent use of macros");
+#else
+static_assert(OGS_EIGEN_DYNAMIC_SHAPE_MATRICES_FLAG == 0,
+ "inconsistent use of macros");
+#endif
+}
+}
diff --git a/ProcessLib/TES/TESLocalAssemblerInner.h b/ProcessLib/TES/TESLocalAssemblerInner.h
new file mode 100644
index 00000000000..3e36cb47180
--- /dev/null
+++ b/ProcessLib/TES/TESLocalAssemblerInner.h
@@ -0,0 +1,98 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * The code of this file is used to decouple the evaluation of matrix elements
+ * from the rest of OGS6,
+ * not all of OGS6 has to be recompiled every time a small change is done.
+ */
+
+#ifndef PROCESS_LIB_TES_FEM_NOTPL_H_
+#define PROCESS_LIB_TES_FEM_NOTPL_H_
+
+#include "NumLib/Fem/ShapeMatrixPolicy.h"
+#include "ProcessLib/LocalAssemblerTraits.h"
+#include "ProcessLib/VariableTransformation.h"
+
+#include "TESLocalAssemblerData.h"
+
+namespace ProcessLib
+{
+namespace TES
+{
+enum class TESIntPtVariables : unsigned
+{
+ SOLID_DENSITY,
+ REACTION_RATE,
+ VELOCITY_X,
+ VELOCITY_Y,
+ VELOCITY_Z,
+ LOADING,
+ REACTION_DAMPING_FACTOR
+};
+
+template <typename Traits>
+class TESLocalAssemblerInner
+{
+public:
+ explicit TESLocalAssemblerInner(AssemblyParams const& ap,
+ const unsigned num_int_pts,
+ const unsigned dimension);
+
+ void assembleIntegrationPoint(
+ unsigned integration_point,
+ std::vector<double> const& localX,
+ typename Traits::ShapeMatrices::ShapeType const& smN,
+ typename Traits::ShapeMatrices::DxShapeType const& smDNdx,
+ typename Traits::ShapeMatrices::JacobianType const& smJ,
+ const double smDetJ,
+ const double weight,
+ typename Traits::LocalMatrix& local_M,
+ typename Traits::LocalMatrix& local_K,
+ typename Traits::LocalVector& local_b);
+
+ void preEachAssemble();
+
+ std::vector<double> const& getIntegrationPointValues(
+ TESIntPtVariables const var, std::vector<double>& cache) const;
+
+ // TODO better encapsulation
+ AssemblyParams const& getAssemblyParameters() const { return _d.ap; }
+ TESFEMReactionAdaptor const& getReactionAdaptor() const
+ {
+ return *_d.reaction_adaptor;
+ }
+ TESFEMReactionAdaptor& getReactionAdaptor() { return *_d.reaction_adaptor; }
+ TESLocalAssemblerData const& getData() const { return _d; }
+private:
+ Eigen::Matrix3d getMassCoeffMatrix(const unsigned int_pt);
+ typename Traits::LaplaceMatrix getLaplaceCoeffMatrix(const unsigned int_pt,
+ const unsigned dim);
+ Eigen::Matrix3d getAdvectionCoeffMatrix(const unsigned int_pt);
+ Eigen::Matrix3d getContentCoeffMatrix(const unsigned int_pt);
+ Eigen::Vector3d getRHSCoeffVector(const unsigned int_pt);
+
+ void preEachAssembleIntegrationPoint(
+ const unsigned int_pt,
+ std::vector<double> const& localX,
+ typename Traits::ShapeMatrices::ShapeType const& smN,
+ typename Traits::ShapeMatrices::DxShapeType const& smDNdx,
+ typename Traits::ShapeMatrices::JacobianType const& smJ,
+ const double smDetJ);
+
+ void initReaction(const unsigned int_pt);
+
+ TESLocalAssemblerData _d;
+};
+
+} // namespace TES
+
+} // namespace ProcessLib
+
+// tricking cmake dependency checker
+#include "TESLocalAssemblerInner-impl-incl.h"
+
+#endif // PROCESS_LIB_TES_FEM_NOTPL_H_
diff --git a/ProcessLib/TES/TESOGS5MaterialModels.cpp b/ProcessLib/TES/TESOGS5MaterialModels.cpp
new file mode 100644
index 00000000000..1a234d686b3
--- /dev/null
+++ b/ProcessLib/TES/TESOGS5MaterialModels.cpp
@@ -0,0 +1,35 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, 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 "TESOGS5MaterialModels.h"
+
+namespace ProcessLib
+{
+namespace TES
+{
+const double FluidHeatConductivityN2::A[5] = {0.46649, -0.57015, 0.19164,
+ -0.03708, 0.00241};
+
+const double FluidHeatConductivityN2::f[9] = {
+ -0.837079888737e3, 0.37914711487e2, -0.601737844275,
+ 0.350418363823e1, -0.874955653028e-5, 0.148968607239e-7,
+ -0.256370354277e-11, 0.100773735767e1, 0.335340610e4};
+
+const double FluidHeatConductivityN2::C[4] = {3.3373542, 0.37098251, 0.89913456,
+ 0.16972505};
+
+const double FluidViscosityN2::A[5] = {0.46649, -0.57015, 0.19164, -0.03708,
+ 0.00241};
+const double FluidViscosityN2::C[5] = {-20.09997, 3.4376416, -1.4470051,
+ -0.027766561, -0.21662362};
+
+const double FluidHeatConductivityH2O::a[4] = {0.0102811, 0.0299621, 0.0156146,
+ -0.00422464};
+
+} // TES
+} // ProcessLib
diff --git a/ProcessLib/TES/TESOGS5MaterialModels.h b/ProcessLib/TES/TESOGS5MaterialModels.h
new file mode 100644
index 00000000000..c4e366d9826
--- /dev/null
+++ b/ProcessLib/TES/TESOGS5MaterialModels.h
@@ -0,0 +1,411 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ */
+
+#ifndef PROCESSLIB_TES_OGS5MATERIALMODELS
+#define PROCESSLIB_TES_OGS5MATERIALMODELS
+
+#include "TESAssemblyParams.h"
+
+namespace
+{
+const double GAS_CONST = 8.3144621;
+}
+
+namespace ProcessLib
+{
+namespace TES
+{
+inline double fluid_density(const double p, const double T, const double x)
+{
+ // OGS-5 density model 26
+
+ const double M0 = ProcessLib::TES::M_N2;
+ const double M1 = ProcessLib::TES::M_H2O;
+
+ const double xn = M0 * x / (M0 * x + M1 * (1.0 - x));
+
+ return p / (GAS_CONST * T) * (M1 * xn + M0 * (1.0 - xn));
+ ;
+}
+
+template <int i>
+double mypow(const double x)
+{
+ if (i < 0)
+ {
+ return 1.0 / mypow<-i>(x);
+ }
+ else
+ {
+ const double p = mypow<(i >> 1)>(x);
+ return (i & 1) ? p * p * x : p * p;
+ }
+}
+
+template <>
+inline double mypow<0>(const double /*x*/)
+{
+ return 1.0;
+}
+
+struct FluidViscosityN2
+{
+ static double get(double rho, double T)
+ {
+ const double rho_c = 314; // [kg/m3]
+ const double CVF = 14.058; // [1e-3 Pa-s]
+
+ const double sigma = 0.36502496e-09;
+ const double k = 1.38062e-23;
+ const double eps = 138.08483e-23;
+ const double c1 = 0.3125;
+ const double c2 = 2.0442e-49;
+
+ const double T_star = T * k / eps;
+ rho = rho / rho_c;
+
+ double Omega = loop1_term<0>(T_star);
+ Omega += loop1_term<1>(T_star);
+ Omega += loop1_term<2>(T_star);
+ Omega += loop1_term<3>(T_star);
+ Omega += loop1_term<4>(T_star);
+
+ Omega = std::exp(Omega);
+
+ // eta in [Pa*s]
+ const double eta_0 = c1 * std::sqrt(c2 * T) / (sigma * sigma * Omega);
+
+ double sum = loop2_term<2>(rho);
+ sum += loop2_term<3>(rho);
+ sum += loop2_term<4>(rho);
+
+ //
+ const double eta_r =
+ CVF * 1e-6 * (C[0] / (rho - C[1]) + C[0] / C[1] + sum);
+
+ return eta_0 + eta_r; // [Pa*s]
+ }
+
+private:
+ template <unsigned i>
+ static double loop1_term(double T_star)
+ {
+ return A[i] * mypow<i>(log(T_star));
+ }
+
+ template <unsigned i>
+ static double loop2_term(double rho)
+ {
+ return C[i] * mypow<i - 1>(rho);
+ }
+
+ static const double A[5];
+ static const double C[5];
+};
+
+struct FluidViscosityH2O
+{
+ static double get(double rho, double T)
+ {
+ double my, my_0, my_1;
+ double H[4];
+
+ T = T / 647.096;
+ rho = rho / 322.0;
+
+ H[0] = 1.67752;
+ H[1] = 2.20462;
+ H[2] = 0.6366564;
+ H[3] = -0.241605;
+
+ double h[6][7] = {0.0};
+ h[0][0] = 0.520094000;
+ h[1][0] = 0.085089500;
+ h[2][0] = -1.083740000;
+ h[3][0] = -0.289555000;
+ h[0][1] = 0.222531000;
+ h[1][1] = 0.999115000;
+ h[2][1] = 1.887970000;
+ h[3][1] = 1.266130000;
+ h[5][1] = 0.120573000;
+ h[0][2] = -0.281378000;
+ h[1][2] = -0.906851000;
+ h[2][2] = -0.772479000;
+ h[3][2] = -0.489837000;
+ h[4][2] = -0.257040000;
+ h[0][3] = 0.161913000;
+ h[1][3] = 0.257399000;
+ h[0][4] = -0.032537200;
+ h[3][4] = 0.069845200;
+ h[4][5] = 0.008721020;
+ h[3][6] = -0.004356730;
+ h[5][6] = -0.000593264;
+
+ double sum1 = H[0] / mypow<0>(T);
+ sum1 += H[1] / mypow<1>(T);
+ sum1 += H[2] / mypow<2>(T);
+ sum1 += H[3] / mypow<3>(T);
+
+ my_0 = 100 * std::sqrt(T) / sum1;
+
+ double sum2 = inner_loop<0>(rho, T, h);
+ sum2 += inner_loop<1>(rho, T, h);
+ sum2 += inner_loop<2>(rho, T, h);
+ sum2 += inner_loop<3>(rho, T, h);
+ sum2 += inner_loop<4>(rho, T, h);
+ sum2 += inner_loop<5>(rho, T, h);
+
+ my_1 = std::exp(rho * sum2);
+
+ my = (my_0 * my_1) / 1e6;
+ return my;
+ }
+
+private:
+ template <int i>
+ static double inner_loop(const double rho,
+ const double T,
+ const double (&h)[6][7])
+ {
+ const double base = rho - 1.0;
+
+ double sum3 = h[i][0] * mypow<0>(base);
+ sum3 += h[i][1] * mypow<1>(base);
+ sum3 += h[i][2] * mypow<2>(base);
+ sum3 += h[i][3] * mypow<3>(base);
+ sum3 += h[i][4] * mypow<4>(base);
+ sum3 += h[i][5] * mypow<5>(base);
+ sum3 += h[i][6] * mypow<6>(base);
+
+ return mypow<i>(1 / T - 1) * sum3;
+ }
+};
+
+inline double fluid_viscosity(const double p, const double T, const double x)
+{
+ // OGS 5 viscosity model 26
+
+ const double M0 = ProcessLib::TES::M_N2;
+ const double M1 = ProcessLib::TES::M_H2O;
+
+ // reactive component
+ const double x0 =
+ M0 * x / (M0 * x + M1 * (1.0 - x)); // mass in mole fraction
+ const double V0 = FluidViscosityH2O::get(M1 * p / (GAS_CONST * T), T);
+ // inert component
+ const double x1 = 1.0 - x0;
+ const double V1 = FluidViscosityN2::get(M0 * p / (GAS_CONST * T), T);
+
+ const double M0_over_M1(M1 / M0); // reactive over inert
+ const double V0_over_V1(V0 / V1);
+
+ const double phi_12 =
+ mypow<2>(1.0 +
+ std::sqrt(V0_over_V1) * std::pow(1.0 / M0_over_M1, 0.25)) /
+ std::sqrt(8.0 * (1.0 + M0_over_M1));
+
+ const double phi_21 = phi_12 * M0_over_M1 / V0_over_V1;
+
+ return V0 * x0 / (x0 + x1 * phi_12) + V1 * x1 / (x1 + x0 * phi_21);
+}
+
+struct FluidHeatConductivityN2
+{
+ static double get(double rho, double T)
+ {
+ const double X1 = 0.95185202;
+ const double X2 = 1.0205422;
+
+ const double rho_c = 314; // [kg/m3]
+ const double M = 28.013;
+ const double k = 1.38062e-23;
+ const double eps = 138.08483e-23;
+ const double N_A = 6.02213E26;
+ const double R = 8.31434;
+ // const double R = GAS_CONST;
+ const double CCF = 4.173; // mW/m/K
+
+ const double c1 = 0.3125;
+ const double c2 = 2.0442e-49;
+ const double sigma = 0.36502496e-09;
+
+ rho /= rho_c;
+
+ // dilute heat conductivity
+ const double sum1 = loop1_term<0>(T) + loop1_term<1>(T) +
+ loop1_term<2>(T) + loop1_term<3>(T) +
+ loop1_term<4>(T) + loop1_term<5>(T) +
+ loop1_term<6>(T);
+ const double temp(std::exp((f[8] / T)) - 1);
+ const double c_v0 =
+ R *
+ (sum1 + ((f[7] * (f[8] / T) * (f[8] / T) * (std::exp((f[8] / T)))) /
+ (temp * temp) -
+ 1));
+
+ double cvint;
+ cvint = c_v0 * 1000 / N_A;
+
+ // dilute gas viscosity
+ const double log_T_star = std::log(T * k / eps);
+
+ const double Omega =
+ std::exp(loop2_term<0>(log_T_star) + loop2_term<1>(log_T_star) +
+ loop2_term<2>(log_T_star) + loop2_term<3>(log_T_star) +
+ loop2_term<4>(log_T_star));
+
+ // eta in [Pa*s]
+ const double eta_0 =
+ 1e6 * (c1 * std::sqrt(c2 * T) / (sigma * sigma * Omega));
+
+ const double F = eta_0 * k * N_A / (M * 1000);
+
+ const double lambda_tr = 2.5 * (1.5 - X1);
+ const double lambda_in = X2 * (cvint / k + X1);
+
+ const double lambda_0 = F * (lambda_tr + lambda_in);
+
+ const double sum2 = loop3_term<0>(rho) + loop3_term<1>(rho) +
+ loop3_term<2>(rho) + loop3_term<3>(rho);
+ const double lambda_r = sum2 * CCF;
+
+ return (lambda_0 + lambda_r) / 1000; // lambda in [W/m/K]
+ }
+
+private:
+ template <int i>
+ static double loop1_term(const double T)
+ {
+ return f[i] * mypow<i - 3>(T);
+ }
+
+ template <int i>
+ static double loop2_term(const double log_T_star)
+ {
+ return A[i] * mypow<i>(log_T_star);
+ }
+
+ template <int i>
+ static double loop3_term(const double rho)
+ {
+ return C[i] * mypow<i + 1>(rho);
+ }
+
+ const static double A[5];
+ const static double f[9];
+ const static double C[4];
+};
+
+struct FluidHeatConductivityH2O
+{
+ static double get(double rho, double T)
+ {
+ double S, Q;
+ double b[3], B[2], d[4], C[6];
+
+ T /= 647.096;
+ rho /= 317.11;
+
+ b[0] = -0.397070;
+ b[1] = 0.400302;
+ b[2] = 1.060000;
+
+ B[0] = -0.171587;
+ B[1] = 2.392190;
+
+ d[0] = 0.0701309;
+ d[1] = 0.0118520;
+ d[2] = 0.00169937;
+ d[3] = -1.0200;
+
+ C[0] = 0.642857;
+ C[1] = -4.11717;
+ C[2] = -6.17937;
+ C[3] = 0.00308976;
+ C[4] = 0.0822994;
+ C[5] = 10.0932;
+
+ const double sum1 = loop_term<0>(T) + loop_term<1>(T) +
+ loop_term<2>(T) + loop_term<3>(T);
+
+ const double lambda_0 = std::sqrt(T) * sum1;
+ const double lambda_1 =
+ b[0] + b[1] * rho +
+ b[2] * std::exp(B[0] * (rho + B[1]) * (rho + B[1]));
+
+ const double dT = fabs(T - 1) + C[3];
+ const double dT_pow_3_5 = std::pow(dT, 3. / 5.);
+ Q = 2 + (C[4] / dT_pow_3_5);
+
+ if (T >= 1)
+ S = 1 / dT;
+ else
+ S = C[5] / dT_pow_3_5;
+
+ const double rho_pow_9_5 = std::pow(rho, 9. / 5.);
+ const double rho_pow_Q = std::pow(rho, Q);
+ const double T_pow_3_2 = T * std::sqrt(T);
+ const double lambda_2 =
+ (d[0] / mypow<10>(T) + d[1]) * rho_pow_9_5 *
+ std::exp(C[0] * (1 - rho * rho_pow_9_5)) +
+ d[2] * S * rho_pow_Q *
+ std::exp((Q / (1. + Q)) * (1 - rho * rho_pow_Q)) +
+ d[3] * std::exp(C[1] * T_pow_3_2 + C[2] / mypow<5>(rho));
+
+ return lambda_0 + lambda_1 + lambda_2; // lambda in [W/m/K]
+ }
+
+private:
+ template <unsigned i>
+ static double loop_term(const double T)
+ {
+ return a[i] * mypow<i>(T);
+ }
+
+ static const double a[4];
+};
+
+inline double fluid_heat_conductivity(const double p,
+ const double T,
+ const double x)
+{
+ // OGS 5 fluid heat conductivity model 11
+
+ const double M0 = ProcessLib::TES::M_N2;
+ const double M1 = ProcessLib::TES::M_H2O;
+
+ // TODO [CL] max() is redundant if the fraction is guaranteed to be between
+ // 0 and 1.
+ // reactive component
+ const double x0 = std::max(M0 * x / (M0 * x + M1 * (1.0 - x)),
+ 0.); // convert mass to mole fraction
+ const double k0 =
+ FluidHeatConductivityH2O::get(M1 * p / (GAS_CONST * T), T);
+ // inert component
+ const double x1 = 1.0 - x0;
+ const double k1 = FluidHeatConductivityN2::get(M0 * p / (GAS_CONST * T), T);
+
+ const double M1_over_M2 = M1 / M0; // reactive over inert
+ const double V1_over_V2 =
+ FluidViscosityH2O::get(M1 * p / (GAS_CONST * T), T) /
+ FluidViscosityN2::get(M0 * p / (GAS_CONST * T), T);
+ const double L1_over_L2 = V1_over_V2 / M1_over_M2;
+
+ const double M12_pow_mquarter = std::pow(M1_over_M2, -0.25);
+ const double phi_12 = (1.0 + std::sqrt(L1_over_L2) * M12_pow_mquarter) *
+ (1.0 + std::sqrt(V1_over_V2) * M12_pow_mquarter) /
+ std::sqrt(8.0 * (1.0 + M1_over_M2));
+ const double phi_21 = phi_12 * M1_over_M2 / V1_over_V2;
+
+ return k0 * x0 / (x0 + x1 * phi_12) + k1 * x1 / (x1 + x0 * phi_21);
+}
+
+} // TES
+} // ProcessLib
+
+#endif // PROCESSLIB_TES_OGS5MATERIALMODELS
diff --git a/ProcessLib/TES/TESProcess.cpp b/ProcessLib/TES/TESProcess.cpp
new file mode 100644
index 00000000000..930bc4d4c92
--- /dev/null
+++ b/ProcessLib/TES/TESProcess.cpp
@@ -0,0 +1,451 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, 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 "TESProcess.h"
+
+// TODO Copied from VectorMatrixAssembler. Could be provided by the DOF table.
+inline AssemblerLib::LocalToGlobalIndexMap::RowColumnIndices
+getRowColumnIndices_(std::size_t const id,
+ AssemblerLib::LocalToGlobalIndexMap const& dof_table,
+ std::vector<GlobalIndexType>& indices)
+{
+ assert(dof_table.size() > id);
+ indices.clear();
+
+ // Local matrices and vectors will always be ordered by component,
+ // no matter what the order of the global matrix is.
+ for (unsigned c = 0; c < dof_table.getNumComponents(); ++c)
+ {
+ auto const& idcs = dof_table(id, c).rows;
+ indices.reserve(indices.size() + idcs.size());
+ indices.insert(indices.end(), idcs.begin(), idcs.end());
+ }
+
+ return AssemblerLib::LocalToGlobalIndexMap::RowColumnIndices(indices,
+ indices);
+}
+
+template <typename GlobalVector>
+void getVectorValues(
+ GlobalVector const& x,
+ AssemblerLib::LocalToGlobalIndexMap::RowColumnIndices const& r_c_indices,
+ std::vector<double>& local_x)
+{
+ auto const& indices = r_c_indices.rows;
+ local_x.clear();
+ local_x.reserve(indices.size());
+
+ for (auto i : indices)
+ {
+ local_x.emplace_back(x.get(i));
+ }
+}
+
+// TODO that essentially duplicates code which is also present in ProcessOutput.
+template <typename GlobalVector>
+double getNodalValue(GlobalVector const& x, MeshLib::Mesh const& mesh,
+ AssemblerLib::LocalToGlobalIndexMap const& dof_table,
+ std::size_t const node_id,
+ std::size_t const global_component_id)
+{
+ MeshLib::Location const l{mesh.getID(), MeshLib::MeshItemType::Node,
+ node_id};
+
+ auto const index = dof_table.getLocalIndex(
+ l, global_component_id, x.getRangeBegin(), x.getRangeEnd());
+
+ return x.get(index);
+}
+
+namespace ProcessLib
+{
+namespace TES
+{
+template <typename GlobalSetup>
+TESProcess<GlobalSetup>::TESProcess(
+ MeshLib::Mesh& mesh,
+ typename Process<GlobalSetup>::NonlinearSolver& nonlinear_solver,
+ std::unique_ptr<typename Process<GlobalSetup>::TimeDiscretization>&&
+ time_discretization,
+ std::vector<std::reference_wrapper<ProcessVariable>>&& process_variables,
+ SecondaryVariableCollection<GlobalVector>&& secondary_variables,
+ ProcessOutput<GlobalVector>&& process_output,
+ const BaseLib::ConfigTree& config)
+ : Process<GlobalSetup>(
+ mesh, nonlinear_solver, std::move(time_discretization),
+ std::move(process_variables), std::move(secondary_variables),
+ std::move(process_output))
+{
+ DBUG("Create TESProcess.");
+
+ // physical parameters for local assembly
+ {
+ std::vector<std::pair<std::string, double*>> params{
+ {"fluid_specific_heat_source",
+ &_assembly_params.fluid_specific_heat_source},
+ {"fluid_specific_isobaric_heat_capacity", &_assembly_params.cpG},
+ {"solid_specific_heat_source",
+ &_assembly_params.solid_specific_heat_source},
+ {"solid_heat_conductivity", &_assembly_params.solid_heat_cond},
+ {"solid_specific_isobaric_heat_capacity", &_assembly_params.cpS},
+ {"tortuosity", &_assembly_params.tortuosity},
+ {"diffusion_coefficient",
+ &_assembly_params.diffusion_coefficient_component},
+ {"porosity", &_assembly_params.poro},
+ {"solid_density_dry", &_assembly_params.rho_SR_dry},
+ {"solid_density_initial", &_assembly_params.initial_solid_density}};
+
+ for (auto const& p : params)
+ {
+ if (auto const par = config.getConfParamOptional<double>(p.first))
+ {
+ DBUG("setting parameter `%s' to value `%g'", p.first.c_str(),
+ *par);
+ *p.second = *par;
+ }
+ }
+ }
+
+ // characteristic values of primary variables
+ {
+ std::vector<std::pair<std::string, Trafo*>> const params{
+ {"characteristic_pressure", &_assembly_params.trafo_p},
+ {"characteristic_temperature", &_assembly_params.trafo_T},
+ {"characteristic_vapour_mass_fraction", &_assembly_params.trafo_x}};
+
+ for (auto const& p : params)
+ {
+ if (auto const par = config.getConfParamOptional<double>(p.first))
+ {
+ INFO("setting parameter `%s' to value `%g'", p.first.c_str(),
+ *par);
+ *p.second = Trafo{*par};
+ }
+ }
+ }
+
+ // permeability
+ if (auto par =
+ config.getConfParamOptional<double>("solid_hydraulic_permeability"))
+ {
+ DBUG(
+ "setting parameter `solid_hydraulic_permeability' to isotropic "
+ "value `%g'",
+ *par);
+ const auto dim = mesh.getDimension();
+ _assembly_params.solid_perm_tensor =
+ Eigen::MatrixXd::Identity(dim, dim) * (*par);
+ }
+
+ // reactive system
+ _assembly_params.react_sys = Adsorption::AdsorptionReaction::newInstance(
+ config.getConfSubtree("reactive_system"));
+
+ // debug output
+ if (auto const param =
+ config.getConfParamOptional<bool>("output_element_matrices"))
+ {
+ DBUG("output_element_matrices: %s", (*param) ? "true" : "false");
+
+ _assembly_params.output_element_matrices = *param;
+ }
+
+ // TODO somewhere else
+ /*
+ if (auto const param =
+ config.getConfParamOptional<bool>("output_global_matrix"))
+ {
+ DBUG("output_global_matrix: %s", (*param) ? "true" : "false");
+
+ this->_process_output.output_global_matrix = *param;
+ }
+ */
+}
+
+template <typename GlobalSetup>
+void TESProcess<GlobalSetup>::initializeConcreteProcess(
+ AssemblerLib::LocalToGlobalIndexMap const& dof_table,
+ MeshLib::Mesh const& mesh, unsigned const integration_order)
+{
+ DBUG("Create global assembler.");
+ _global_assembler.reset(new GlobalAssembler(dof_table));
+
+ ProcessLib::createLocalAssemblers<GlobalSetup, TESLocalAssembler>(
+ mesh.getDimension(), mesh.getElements(), dof_table, integration_order,
+ _local_assemblers, _assembly_params);
+
+ // TODO move the two data members somewhere else.
+ // for extrapolation of secondary variables
+ std::vector<std::unique_ptr<MeshLib::MeshSubsets>>
+ all_mesh_subsets_single_component;
+ all_mesh_subsets_single_component.emplace_back(
+ new MeshLib::MeshSubsets(this->_mesh_subset_all_nodes.get()));
+ _local_to_global_index_map_single_component.reset(
+ new AssemblerLib::LocalToGlobalIndexMap(
+ std::move(all_mesh_subsets_single_component),
+ // by location order is needed for output
+ AssemblerLib::ComponentOrder::BY_LOCATION));
+
+ _extrapolator.reset(
+ new ExtrapolatorImplementation(MathLib::MatrixSpecifications(
+ 0u, 0u, nullptr, _local_to_global_index_map_single_component.get(),
+ &mesh)));
+
+ // secondary variables
+ auto add2nd = [&](std::string const& var_name, unsigned const n_components,
+ SecondaryVariableFunctions<GlobalVector>&& fcts) {
+ this->_secondary_variables.addSecondaryVariable(var_name, n_components,
+ std::move(fcts));
+ };
+ auto makeEx =
+ [&](TESIntPtVariables var) -> SecondaryVariableFunctions<GlobalVector> {
+ return ProcessLib::makeExtrapolator(var, *_extrapolator,
+ _local_assemblers);
+ };
+
+ add2nd("solid_density", 1, makeEx(TESIntPtVariables::SOLID_DENSITY));
+ add2nd("reaction_rate", 1, makeEx(TESIntPtVariables::REACTION_RATE));
+ add2nd("velocity_x", 1, makeEx(TESIntPtVariables::VELOCITY_X));
+ if (mesh.getDimension() >= 2)
+ add2nd("velocity_y", 1, makeEx(TESIntPtVariables::VELOCITY_Y));
+ if (mesh.getDimension() >= 3)
+ add2nd("velocity_z", 1, makeEx(TESIntPtVariables::VELOCITY_Z));
+
+ add2nd("loading", 1, makeEx(TESIntPtVariables::LOADING));
+ add2nd("reaction_damping_factor", 1,
+ makeEx(TESIntPtVariables::REACTION_DAMPING_FACTOR));
+
+ namespace PH = std::placeholders;
+ using Self = TESProcess<GlobalSetup>;
+
+ add2nd("vapour_partial_pressure", 1,
+ {std::bind(&Self::computeVapourPartialPressure, this, PH::_1, PH::_2,
+ PH::_3),
+ nullptr});
+ add2nd("relative_humidity", 1, {std::bind(&Self::computeRelativeHumidity,
+ this, PH::_1, PH::_2, PH::_3),
+ nullptr});
+ add2nd("equilibrium_loading", 1,
+ {std::bind(&Self::computeEquilibriumLoading, this, PH::_1, PH::_2,
+ PH::_3),
+ nullptr});
+}
+
+template <typename GlobalSetup>
+void TESProcess<GlobalSetup>::assembleConcreteProcess(const double t,
+ GlobalVector const& x,
+ GlobalMatrix& M,
+ GlobalMatrix& K,
+ GlobalVector& b)
+{
+ DBUG("Assemble TESProcess.");
+
+ // Call global assembler for each local assembly item.
+ GlobalSetup::executeMemberDereferenced(*_global_assembler,
+ &GlobalAssembler::assemble,
+ _local_assemblers, t, x, M, K, b);
+}
+
+template <typename GlobalSetup>
+void TESProcess<GlobalSetup>::preTimestep(GlobalVector const& x, const double t,
+ const double delta_t)
+{
+ DBUG("new timestep");
+
+ _assembly_params.delta_t = delta_t;
+ _assembly_params.current_time = t;
+ ++_assembly_params.timestep; // TODO remove that
+
+ _x_previous_timestep =
+ MathLib::MatrixVectorTraits<GlobalVector>::newInstance(x);
+}
+
+template <typename GlobalSetup>
+void TESProcess<GlobalSetup>::preIteration(const unsigned iter,
+ GlobalVector const& /*x*/)
+{
+ _assembly_params.iteration_in_current_timestep = iter;
+ ++_assembly_params.total_iteration;
+ ++_assembly_params.number_of_try_of_iteration;
+}
+
+template <typename GlobalSetup>
+NumLib::IterationResult TESProcess<GlobalSetup>::postIteration(
+ GlobalVector const& x)
+{
+ if (this->_process_output.output_iteration_results)
+ {
+ DBUG("output results of iteration %li",
+ _assembly_params.total_iteration);
+ std::string fn =
+ "tes_iter_" + std::to_string(_assembly_params.total_iteration) +
+ +"_ts_" + std::to_string(_assembly_params.timestep) + "_" +
+ std::to_string(_assembly_params.iteration_in_current_timestep) +
+ "_" + std::to_string(_assembly_params.number_of_try_of_iteration) +
+ ".vtu";
+
+ this->output(fn, 0, x);
+ }
+
+ bool check_passed = true;
+
+ if (!Trafo::constrained)
+ {
+ // bounds checking only has to happen if the vapour mass fraction is
+ // non-logarithmic.
+
+ std::vector<GlobalIndexType> indices_cache;
+ std::vector<double> local_x_cache;
+ std::vector<double> local_x_prev_ts_cache;
+
+ auto check_variable_bounds = [&](std::size_t id,
+ LocalAssembler& loc_asm) {
+ auto const r_c_indices = getRowColumnIndices_(
+ id, *this->_local_to_global_index_map, indices_cache);
+ getVectorValues(x, r_c_indices, local_x_cache);
+ getVectorValues(*_x_previous_timestep, r_c_indices,
+ local_x_prev_ts_cache);
+
+ if (!loc_asm.checkBounds(local_x_cache, local_x_prev_ts_cache))
+ check_passed = false;
+ };
+
+ GlobalSetup::executeDereferenced(check_variable_bounds,
+ _local_assemblers);
+ }
+
+ if (!check_passed)
+ return NumLib::IterationResult::REPEAT_ITERATION;
+
+ // TODO remove
+ DBUG("ts %lu iteration %lu (in current ts: %lu) try %u accepted",
+ _assembly_params.timestep, _assembly_params.total_iteration,
+ _assembly_params.iteration_in_current_timestep,
+ _assembly_params.number_of_try_of_iteration);
+
+ _assembly_params.number_of_try_of_iteration = 0;
+
+ return NumLib::IterationResult::SUCCESS;
+}
+
+template <typename GlobalSetup>
+typename TESProcess<GlobalSetup>::GlobalVector const&
+TESProcess<GlobalSetup>::computeVapourPartialPressure(
+ typename TESProcess::GlobalVector const& x,
+ AssemblerLib::LocalToGlobalIndexMap const& dof_table,
+ std::unique_ptr<GlobalVector>& result_cache)
+{
+ assert(&dof_table == this->_local_to_global_index_map.get());
+
+ result_cache = MathLib::MatrixVectorTraits<GlobalVector>::newInstance(
+ {0, 0, nullptr, _local_to_global_index_map_single_component.get(),
+ nullptr});
+
+ GlobalIndexType nnodes = this->_mesh.getNNodes();
+
+ for (GlobalIndexType node_id = 0; node_id < nnodes; ++node_id)
+ {
+ auto const p = getNodalValue(x, this->_mesh, dof_table, node_id,
+ COMPONENT_ID_PRESSURE);
+ auto const x_mV = getNodalValue(x, this->_mesh, dof_table, node_id,
+ COMPONENT_ID_MASS_FRACTION);
+
+ auto const x_nV = Adsorption::AdsorptionReaction::getMolarFraction(
+ x_mV, _assembly_params.M_react, _assembly_params.M_inert);
+
+ // TODO Problems with PETSc? (local vs. global index)
+ result_cache->set(node_id, p * x_nV);
+ }
+
+ return *result_cache;
+}
+
+template <typename GlobalSetup>
+typename TESProcess<GlobalSetup>::GlobalVector const&
+TESProcess<GlobalSetup>::computeRelativeHumidity(
+ typename TESProcess::GlobalVector const& x,
+ AssemblerLib::LocalToGlobalIndexMap const& dof_table,
+ std::unique_ptr<GlobalVector>& result_cache)
+{
+ assert(&dof_table == this->_local_to_global_index_map.get());
+
+ result_cache = MathLib::MatrixVectorTraits<GlobalVector>::newInstance(
+ {0, 0, nullptr, _local_to_global_index_map_single_component.get(),
+ nullptr});
+
+ GlobalIndexType nnodes = this->_mesh.getNNodes();
+
+ for (GlobalIndexType node_id = 0; node_id < nnodes; ++node_id)
+ {
+ auto const p = getNodalValue(x, this->_mesh, dof_table, node_id,
+ COMPONENT_ID_PRESSURE);
+ auto const T = getNodalValue(x, this->_mesh, dof_table, node_id,
+ COMPONENT_ID_TEMPERATURE);
+ auto const x_mV = getNodalValue(x, this->_mesh, dof_table, node_id,
+ COMPONENT_ID_MASS_FRACTION);
+
+ auto const x_nV = Adsorption::AdsorptionReaction::getMolarFraction(
+ x_mV, _assembly_params.M_react, _assembly_params.M_inert);
+
+ auto const p_S =
+ Adsorption::AdsorptionReaction::getEquilibriumVapourPressure(T);
+
+ // TODO Problems with PETSc? (local vs. global index)
+ result_cache->set(node_id, p * x_nV / p_S);
+ }
+
+ return *result_cache;
+}
+
+template <typename GlobalSetup>
+typename TESProcess<GlobalSetup>::GlobalVector const&
+TESProcess<GlobalSetup>::computeEquilibriumLoading(
+ typename TESProcess::GlobalVector const& x,
+ AssemblerLib::LocalToGlobalIndexMap const& dof_table,
+ std::unique_ptr<GlobalVector>& result_cache)
+{
+ assert(&dof_table == this->_local_to_global_index_map.get());
+
+ result_cache = MathLib::MatrixVectorTraits<GlobalVector>::newInstance(
+ {0, 0, nullptr, _local_to_global_index_map_single_component.get(),
+ nullptr});
+
+ GlobalIndexType nnodes = this->_mesh.getNNodes();
+
+ for (GlobalIndexType node_id = 0; node_id < nnodes; ++node_id)
+ {
+ auto const p = getNodalValue(x, this->_mesh, dof_table, node_id,
+ COMPONENT_ID_PRESSURE);
+ auto const T = getNodalValue(x, this->_mesh, dof_table, node_id,
+ COMPONENT_ID_TEMPERATURE);
+ auto const x_mV = getNodalValue(x, this->_mesh, dof_table, node_id,
+ COMPONENT_ID_MASS_FRACTION);
+
+ auto const x_nV = Adsorption::AdsorptionReaction::getMolarFraction(
+ x_mV, _assembly_params.M_react, _assembly_params.M_inert);
+
+ auto const p_V = p * x_nV;
+ auto const C_eq =
+ (p_V <= 0.0) ? 0.0
+ : _assembly_params.react_sys->getEquilibriumLoading(
+ p_V, T, _assembly_params.M_react);
+
+ // TODO Problems with PETSc? (local vs. global index)
+ result_cache->set(node_id, C_eq);
+ }
+
+ return *result_cache;
+}
+
+// Explicitly instantiate TESProcess for GlobalSetupType.
+template class TESProcess<GlobalSetupType>;
+
+} // namespace TES
+
+} // namespace ProcessLib
diff --git a/ProcessLib/TES/TESProcess.h b/ProcessLib/TES/TESProcess.h
new file mode 100644
index 00000000000..5149b329af5
--- /dev/null
+++ b/ProcessLib/TES/TESProcess.h
@@ -0,0 +1,149 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef PROCESS_LIB_TESPROCESS_H_
+#define PROCESS_LIB_TESPROCESS_H_
+
+#include "AssemblerLib/LocalToGlobalIndexMap.h"
+#include "NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h"
+#include "ProcessLib/Process.h"
+
+#include "TESAssemblyParams.h"
+#include "TESLocalAssembler.h"
+
+namespace MeshLib
+{
+class Element;
+class Mesh;
+template <typename PROP_VAL_TYPE>
+class PropertyVector;
+}
+
+namespace ProcessLib
+{
+namespace TES
+{
+template <typename GlobalSetup>
+class TESProcess final : public Process<GlobalSetup>
+{
+ using BP = Process<GlobalSetup>; //!< "Base Process"
+
+public:
+ using GlobalVector = typename GlobalSetup::VectorType;
+ using GlobalMatrix = typename GlobalSetup::MatrixType;
+
+ TESProcess(
+ MeshLib::Mesh& mesh,
+ typename Process<GlobalSetup>::NonlinearSolver& nonlinear_solver,
+ std::unique_ptr<typename Process<GlobalSetup>::TimeDiscretization>&&
+ time_discretization,
+ std::vector<std::reference_wrapper<ProcessVariable>>&&
+ process_variables,
+ SecondaryVariableCollection<GlobalVector>&& secondary_variables,
+ ProcessOutput<GlobalVector>&& process_output,
+ BaseLib::ConfigTree const& config);
+
+ void preTimestep(GlobalVector const& x, const double t,
+ const double delta_t) override;
+ void preIteration(const unsigned iter, GlobalVector const& x) override;
+ NumLib::IterationResult postIteration(GlobalVector const& x) override;
+
+ bool isLinear() const override { return false; }
+private:
+ using LocalAssembler =
+ TESLocalAssemblerInterface<GlobalMatrix, GlobalVector>;
+
+ using GlobalAssembler = AssemblerLib::VectorMatrixAssembler<
+ GlobalMatrix, GlobalVector, LocalAssembler,
+ NumLib::ODESystemTag::FirstOrderImplicitQuasilinear>;
+
+ using ExtrapolatorInterface =
+ NumLib::Extrapolator<GlobalVector, TESIntPtVariables, LocalAssembler>;
+ using ExtrapolatorImplementation =
+ NumLib::LocalLinearLeastSquaresExtrapolator<
+ GlobalVector, TESIntPtVariables, LocalAssembler>;
+
+ void initializeConcreteProcess(
+ AssemblerLib::LocalToGlobalIndexMap const& dof_table,
+ MeshLib::Mesh const& mesh, unsigned const integration_order) override;
+
+ void assembleConcreteProcess(const double t, GlobalVector const& x,
+ GlobalMatrix& M, GlobalMatrix& K,
+ GlobalVector& b) override;
+
+ GlobalVector const& computeVapourPartialPressure(
+ GlobalVector const& x,
+ AssemblerLib::LocalToGlobalIndexMap const& dof_table,
+ std::unique_ptr<GlobalVector>& result_cache);
+
+ GlobalVector const& computeRelativeHumidity(
+ GlobalVector const& x,
+ AssemblerLib::LocalToGlobalIndexMap const& dof_table,
+ std::unique_ptr<GlobalVector>& result_cache);
+
+ GlobalVector const& computeEquilibriumLoading(
+ GlobalVector const& x,
+ AssemblerLib::LocalToGlobalIndexMap const& dof_table,
+ std::unique_ptr<GlobalVector>& result_cache);
+
+ std::unique_ptr<GlobalAssembler> _global_assembler;
+ std::vector<std::unique_ptr<LocalAssembler>> _local_assemblers;
+
+ AssemblyParams _assembly_params;
+
+ std::unique_ptr<AssemblerLib::LocalToGlobalIndexMap>
+ _local_to_global_index_map_single_component;
+
+ std::unique_ptr<ExtrapolatorInterface> _extrapolator;
+
+ // used for checkBounds()
+ std::unique_ptr<GlobalVector> _x_previous_timestep;
+};
+
+template <typename GlobalSetup>
+std::unique_ptr<TESProcess<GlobalSetup>> createTESProcess(
+ MeshLib::Mesh& mesh,
+ typename Process<GlobalSetup>::NonlinearSolver& nonlinear_solver,
+ std::unique_ptr<typename Process<GlobalSetup>::TimeDiscretization>&&
+ time_discretization,
+ std::vector<ProcessVariable> const& variables,
+ std::vector<std::unique_ptr<ParameterBase>> const& /*parameters*/,
+ BaseLib::ConfigTree const& config)
+{
+ config.checkConfParam("type", "TES");
+
+ DBUG("Create TESProcess.");
+
+ auto process_variables = findProcessVariables(
+ variables, config,
+ {"fluid_pressure", "temperature", "vapour_mass_fraction"});
+
+ SecondaryVariableCollection<typename GlobalSetup::VectorType>
+ secondary_variables{
+ config.getConfSubtreeOptional("secondary_variables"),
+ {"solid_density", "reaction_rate", "velocity_x", "velocity_y",
+ "velocity_z", "loading", "reaction_damping_factor",
+ "vapour_partial_pressure", "relative_humidity",
+ "equilibrium_loading"}};
+
+ ProcessOutput<typename GlobalSetup::VectorType> process_output{
+ config.getConfSubtree("output"), process_variables,
+ secondary_variables};
+
+ return std::unique_ptr<TESProcess<GlobalSetup>>{new TESProcess<GlobalSetup>{
+ mesh, nonlinear_solver, std::move(time_discretization),
+ std::move(process_variables), std::move(secondary_variables),
+ std::move(process_output), config}};
+}
+
+} // namespace TES
+
+} // namespace ProcessLib
+
+#endif // PROCESS_LIB_TESPROCESS_H_
diff --git a/ProcessLib/TES/TESReactionAdaptor.cpp b/ProcessLib/TES/TESReactionAdaptor.cpp
new file mode 100644
index 00000000000..8259ba24bb3
--- /dev/null
+++ b/ProcessLib/TES/TESReactionAdaptor.cpp
@@ -0,0 +1,373 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, 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 <cassert>
+
+#include <logog/include/logog.hpp>
+
+#include "MathLib/Nonlinear/Root1D.h"
+#include "MathLib/ODE/ODESolverBuilder.h"
+
+#include "MaterialsLib/Adsorption/Adsorption.h"
+#include "MaterialsLib/Adsorption/ReactionInert.h"
+#include "MaterialsLib/Adsorption/ReactionSinusoidal.h"
+
+#include "TESLocalAssemblerInner.h"
+#include "TESReactionAdaptor.h"
+
+namespace ProcessLib
+{
+namespace TES
+{
+std::unique_ptr<TESFEMReactionAdaptor> TESFEMReactionAdaptor::newInstance(
+ TESLocalAssemblerData const& data)
+{
+ auto const* ads = data.ap.react_sys.get();
+ if (dynamic_cast<Adsorption::AdsorptionReaction const*>(ads) != nullptr)
+ {
+ return std::unique_ptr<TESFEMReactionAdaptor>(
+ new TESFEMReactionAdaptorAdsorption(data));
+ }
+ else if (dynamic_cast<Adsorption::ReactionInert const*>(ads) != nullptr)
+ {
+ return std::unique_ptr<TESFEMReactionAdaptor>(
+ new TESFEMReactionAdaptorInert(data));
+ }
+ else if (dynamic_cast<Adsorption::ReactionSinusoidal const*>(ads) !=
+ nullptr)
+ {
+ return std::unique_ptr<TESFEMReactionAdaptor>(
+ new TESFEMReactionAdaptorSinusoidal(data));
+ }
+ else if (dynamic_cast<Adsorption::ReactionCaOH2 const*>(ads) != nullptr)
+ {
+ return std::unique_ptr<TESFEMReactionAdaptor>(
+ new TESFEMReactionAdaptorCaOH2(data));
+ }
+
+ ERR("No suitable TESFEMReactionAdaptor found. Aborting.");
+ std::abort();
+ return std::unique_ptr<TESFEMReactionAdaptor>(nullptr);
+}
+
+TESFEMReactionAdaptorAdsorption::TESFEMReactionAdaptorAdsorption(
+ TESLocalAssemblerData const& data)
+ // caution fragile: this relies in this constructor b eing called __after__
+ // data.solid_density has been properly set up!
+ : _bounds_violation(data.solid_density.size(), false),
+ _d(data)
+{
+ assert(dynamic_cast<Adsorption::AdsorptionReaction const*>(
+ data.ap.react_sys.get()) != nullptr &&
+ "Reactive system has wrong type.");
+ assert(_bounds_violation.size() != 0);
+}
+
+ReactionRate
+TESFEMReactionAdaptorAdsorption::initReaction_slowDownUndershootStrategy(
+ const unsigned int_pt)
+{
+ assert(_d.ap.number_of_try_of_iteration <= 20);
+
+ const double loading = Adsorption::AdsorptionReaction::getLoading(
+ _d.solid_density_prev_ts[int_pt], _d.ap.rho_SR_dry);
+
+ double react_rate_R =
+ _d.ap.react_sys->getReactionRate(_d.p_V, _d.T, _d.ap.M_react, loading) *
+ _d.ap.rho_SR_dry;
+
+ // set reaction rate based on current damping factor
+ react_rate_R = (_reaction_damping_factor > 1e-3)
+ ? _reaction_damping_factor * react_rate_R
+ : 0.0;
+
+ if (_d.p_V <
+ 0.01 * Adsorption::AdsorptionReaction::getEquilibriumVapourPressure(
+ _d.T) &&
+ react_rate_R > 0.0)
+ {
+ react_rate_R = 0.0;
+ }
+ else if (_d.p_V < 100.0 ||
+ _d.p_V < 0.05 * Adsorption::AdsorptionReaction::
+ getEquilibriumVapourPressure(_d.T))
+ {
+ // use equilibrium reaction for dry regime
+
+ // in the case of zeroth try in zeroth iteration: _p_V and loading are
+ // the values
+ // at the end of the previous timestep
+
+ const double pV_eq = estimateAdsorptionEquilibrium(_d.p_V, loading);
+ // TODO [CL]: it would be more correct to subtract pV from the previous
+ // timestep here
+ const double delta_pV = pV_eq - _d.p_V;
+ const double delta_rhoV = delta_pV * _d.ap.M_react /
+ Adsorption::GAS_CONST / _d.T * _d.ap.poro;
+ const double delta_rhoSR = delta_rhoV / (_d.ap.poro - 1.0);
+ double react_rate_R2 = delta_rhoSR / _d.ap.delta_t;
+
+ if (_bounds_violation[int_pt])
+ {
+ react_rate_R2 *= 0.5;
+ }
+
+ // 0th try: make sure reaction is not slower than allowed by local
+ // estimation
+ // nth try: make sure reaction is not made faster by local estimation
+ if ((_d.ap.number_of_try_of_iteration == 1 &&
+ std::abs(react_rate_R2) > std::abs(react_rate_R)) ||
+ (_d.ap.number_of_try_of_iteration > 1 &&
+ std::abs(react_rate_R2) < std::abs(react_rate_R)))
+ {
+ react_rate_R = react_rate_R2;
+ }
+ }
+
+ // smooth out readjustment of reaction rate
+ if (_d.ap.iteration_in_current_timestep > 4)
+ {
+ if (_d.ap.iteration_in_current_timestep <= 9)
+ {
+ // update reaction rate for for five iterations
+ const auto N = _d.ap.iteration_in_current_timestep - 4;
+
+ // take average s.t. does not oscillate so much
+ react_rate_R = 1.0 / (1.0 + N) *
+ (N * _d.reaction_rate[int_pt] + 1.0 * react_rate_R);
+ }
+ else
+ {
+ // afterwards no update anymore
+ react_rate_R = _d.reaction_rate[int_pt];
+ }
+ }
+
+ if (_d.ap.number_of_try_of_iteration > 1)
+ {
+ // assert that within tries reaction does not get faster
+ // (e.g. due to switch equilibrium reaction <--> kinetic reaction)
+
+ // factor of 0.9*N: in fact, even slow down reaction over tries
+ const double r = std::pow(0.9, _d.ap.number_of_try_of_iteration - 1) *
+ _d.reaction_rate[int_pt];
+ if (std::abs(react_rate_R) > std::abs(r))
+ {
+ react_rate_R = r;
+ }
+ }
+
+ return {react_rate_R,
+ _d.solid_density_prev_ts[int_pt] + react_rate_R * _d.ap.delta_t};
+}
+
+double TESFEMReactionAdaptorAdsorption::estimateAdsorptionEquilibrium(
+ const double p_V0, const double C0) const
+{
+ auto f = [this, p_V0, C0](double pV) -> double {
+ // pV0 := _p_V
+ const double C_eq =
+ _d.ap.react_sys->getEquilibriumLoading(pV, _d.T, _d.ap.M_react);
+ return (pV - p_V0) * _d.ap.M_react / Adsorption::GAS_CONST / _d.T *
+ _d.ap.poro +
+ (1.0 - _d.ap.poro) * (C_eq - C0) * _d.ap.rho_SR_dry;
+ };
+
+ // range where to search for roots of f
+ const double C_eq0 =
+ _d.ap.react_sys->getEquilibriumLoading(p_V0, _d.T, _d.ap.M_react);
+ const double limit =
+ (C_eq0 > C0)
+ ? 1e-8
+ : Adsorption::AdsorptionReaction::getEquilibriumVapourPressure(
+ _d.T);
+
+ // search for roots
+ auto rf = MathLib::Nonlinear::makeRegulaFalsi<MathLib::Nonlinear::Pegasus>(
+ f, p_V0, limit);
+ rf.step(3);
+
+ // set vapour pressure
+ return rf.getResult();
+}
+
+bool TESFEMReactionAdaptorAdsorption::checkBounds(
+ std::vector<double> const& local_x,
+ std::vector<double> const& local_x_prev_ts)
+{
+ double alpha = 1.0;
+
+ const double min_xmV = 1e-6;
+ const std::size_t nnodes = local_x.size() / NODAL_DOF;
+ const std::size_t xmV_offset = COMPONENT_ID_MASS_FRACTION * nnodes;
+
+ for (std::size_t i = 0; i < nnodes; ++i)
+ {
+ auto const xnew = local_x[xmV_offset + i];
+ auto const xold = local_x_prev_ts[xmV_offset + i];
+
+ if (xnew < min_xmV)
+ {
+ const auto a = xold / (xold - xnew);
+ alpha = std::min(alpha, a);
+ _bounds_violation[i] = true;
+ }
+ else if (xnew > 1.0)
+ {
+ const auto a = xold / (xnew - xold);
+ alpha = std::min(alpha, a);
+ _bounds_violation[i] = true;
+ }
+ else
+ {
+ _bounds_violation[i] = false;
+ }
+ }
+
+ assert(alpha > 0.0);
+
+ if (alpha != 1.0)
+ {
+ if (alpha > 0.5)
+ alpha = 0.5;
+ if (alpha < 0.05)
+ alpha = 0.05;
+
+ if (_d.ap.number_of_try_of_iteration <= 3)
+ {
+ _reaction_damping_factor *= sqrt(alpha);
+ }
+ else
+ {
+ _reaction_damping_factor *= alpha;
+ }
+ }
+
+ return alpha == 1.0;
+}
+
+void TESFEMReactionAdaptorAdsorption::preZerothTryAssemble()
+{
+ if (_reaction_damping_factor < 1e-3)
+ _reaction_damping_factor = 1e-3;
+
+ _reaction_damping_factor = std::min(std::sqrt(_reaction_damping_factor),
+ 10.0 * _reaction_damping_factor);
+}
+
+TESFEMReactionAdaptorInert::TESFEMReactionAdaptorInert(
+ TESLocalAssemblerData const& data)
+ : _d(data)
+{
+}
+
+ReactionRate TESFEMReactionAdaptorInert::initReaction(const unsigned int_pt)
+{
+ return {0.0, _d.solid_density_prev_ts[int_pt]};
+ // _d._qR = 0.0;
+ // _d._reaction_rate[int_pt] = 0.0;
+}
+
+TESFEMReactionAdaptorSinusoidal::TESFEMReactionAdaptorSinusoidal(
+ TESLocalAssemblerData const& data)
+ : _d(data)
+{
+ assert(dynamic_cast<Adsorption::ReactionSinusoidal const*>(
+ data.ap.react_sys.get()) != nullptr &&
+ "Reactive system has wrong type.");
+}
+
+ReactionRate TESFEMReactionAdaptorSinusoidal::initReaction(
+ const unsigned /*int_pt*/)
+{
+ const double t = _d.ap.current_time;
+
+ // Cf. OGS5
+ const double rhoSR0 = 1.0;
+ const double rhoTil = 0.1;
+ const double omega = 2.0 * 3.1416;
+ const double poro = _d.ap.poro;
+
+ return {rhoTil * omega * cos(omega * t) / (1.0 - poro),
+ rhoSR0 + rhoTil * std::sin(omega * t) / (1.0 - poro)};
+}
+
+TESFEMReactionAdaptorCaOH2::TESFEMReactionAdaptorCaOH2(
+ TESLocalAssemblerData const& data)
+ : _d(data),
+ _react(dynamic_cast<Adsorption::ReactionCaOH2&>(*data.ap.react_sys.get()))
+{
+ _ode_solver = MathLib::ODE::createODESolver<1>(_react.getOdeSolverConfig());
+ // TODO invalidate config
+
+ _ode_solver->setTolerance(1e-10, 1e-10);
+
+ auto f = [this](const double /*t*/,
+ MathLib::ODE::MappedConstVector<1> const y,
+ MathLib::ODE::MappedVector<1>
+ ydot) -> bool {
+ ydot[0] = _react.getReactionRate(y[0]);
+ return true;
+ };
+
+ _ode_solver->setFunction(f, nullptr);
+}
+
+ReactionRate TESFEMReactionAdaptorCaOH2::initReaction(const unsigned int int_pt)
+{
+ // TODO if the first holds, the second also has to hold
+ if (_d.ap.iteration_in_current_timestep > 1 ||
+ _d.ap.number_of_try_of_iteration > 1)
+ {
+ return {_d.reaction_rate[int_pt], _d.solid_density[int_pt]};
+ }
+
+ // TODO: double check!
+ // const double xv_NR = SolidProp->non_reactive_solid_volume_fraction;
+ // const double rho_NR = SolidProp->non_reactive_solid_density;
+ const double xv_NR = 0.0;
+ const double rho_NR = 0.0;
+
+ const double t0 = 0.0;
+ const double y0 =
+ (_d.solid_density_prev_ts[int_pt] - xv_NR * rho_NR) / (1.0 - xv_NR);
+
+ const double t_end = _d.ap.delta_t;
+
+ _react.updateParam(_d.T, _d.p, _d.vapour_mass_fraction,
+ _d.solid_density_prev_ts[int_pt]);
+
+ _ode_solver->setIC(t0, {y0});
+ _ode_solver->preSolve();
+ _ode_solver->solve(t_end);
+
+ const double time_reached = _ode_solver->getTime();
+ (void)time_reached;
+ assert(std::abs(t_end - time_reached) <
+ std::numeric_limits<double>::epsilon());
+
+ auto const& y_new = _ode_solver->getSolution();
+ auto const& y_dot_new = _ode_solver->getYDot(t_end, y_new);
+
+ double rho_react;
+
+ // cut off when limits are reached
+ if (y_new[0] < _react.rho_low)
+ rho_react = _react.rho_low;
+ else if (y_new[0] > _react.rho_up)
+ rho_react = _react.rho_up;
+ else
+ rho_react = y_new[0];
+
+ return {y_dot_new[0] * (1.0 - xv_NR),
+ (1.0 - xv_NR) * rho_react + xv_NR * rho_NR};
+}
+
+} // namespace TES
+} // namespace ProcessLib
diff --git a/ProcessLib/TES/TESReactionAdaptor.h b/ProcessLib/TES/TESReactionAdaptor.h
new file mode 100644
index 00000000000..a1742dc78bb
--- /dev/null
+++ b/ProcessLib/TES/TESReactionAdaptor.h
@@ -0,0 +1,133 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef PROCESSLIB_TES_TESREACTIONADAPTOR_H
+#define PROCESSLIB_TES_TESREACTIONADAPTOR_H
+
+#include <memory>
+#include <vector>
+
+#include "MaterialsLib/Adsorption/ReactionCaOH2.h"
+#include "MathLib/ODE/ODESolver.h"
+
+namespace ProcessLib
+{
+namespace TES
+{
+struct TESLocalAssemblerData;
+
+struct ReactionRate
+{
+ const double reaction_rate;
+ const double solid_density;
+};
+
+class TESFEMReactionAdaptor
+{
+public:
+ virtual bool checkBounds(std::vector<double> const& /*local_x*/,
+ std::vector<double> const& /*local_x_prev_ts*/)
+ {
+ return true; // by default accept everything
+ }
+
+ virtual ReactionRate initReaction(const unsigned int_pt) = 0;
+
+ virtual void preZerothTryAssemble() {}
+ // TODO: remove
+ virtual double getReactionDampingFactor() const { return -1.0; }
+ virtual ~TESFEMReactionAdaptor() = default;
+
+ static std::unique_ptr<TESFEMReactionAdaptor> newInstance(
+ TESLocalAssemblerData const& data);
+};
+
+class TESFEMReactionAdaptorAdsorption final : public TESFEMReactionAdaptor
+{
+public:
+ explicit TESFEMReactionAdaptorAdsorption(TESLocalAssemblerData const& data);
+
+ bool checkBounds(std::vector<double> const& local_x,
+ std::vector<double> const& local_x_prev_ts) override;
+
+ ReactionRate initReaction(const unsigned int_pt) override
+ {
+ return initReaction_slowDownUndershootStrategy(int_pt);
+ }
+
+ void preZerothTryAssemble() override;
+
+ // TODO: get rid of
+ double getReactionDampingFactor() const override
+ {
+ return _reaction_damping_factor;
+ }
+
+private:
+ ReactionRate initReaction_slowDownUndershootStrategy(const unsigned int_pt);
+
+ /// returns estimated equilibrium vapour pressure
+ /// based on a local (i.e. no diffusion/advection) balance
+ /// of adsorbate loading and vapour partial pressure
+ double estimateAdsorptionEquilibrium(const double p_V0,
+ const double C0) const;
+
+ double _reaction_damping_factor = 1.0;
+ std::vector<bool> _bounds_violation;
+
+ TESLocalAssemblerData const& _d;
+};
+
+class TESFEMReactionAdaptorInert final : public TESFEMReactionAdaptor
+{
+public:
+ explicit TESFEMReactionAdaptorInert(TESLocalAssemblerData const&);
+
+ ReactionRate initReaction(const unsigned int_pt) override;
+
+private:
+ TESLocalAssemblerData const& _d;
+};
+
+class TESFEMReactionAdaptorSinusoidal final : public TESFEMReactionAdaptor
+{
+public:
+ explicit TESFEMReactionAdaptorSinusoidal(TESLocalAssemblerData const& data);
+
+ ReactionRate initReaction(const unsigned) override;
+
+private:
+ TESLocalAssemblerData const& _d;
+};
+
+class TESFEMReactionAdaptorCaOH2 final : public TESFEMReactionAdaptor
+{
+public:
+ explicit TESFEMReactionAdaptorCaOH2(TESLocalAssemblerData const& data);
+
+ ReactionRate initReaction(const unsigned) override;
+
+private:
+ using Data = TESLocalAssemblerData;
+ using React = Adsorption::ReactionCaOH2;
+ Data const& _d;
+ React& _react;
+
+ std::unique_ptr<MathLib::ODE::ODESolver<1>> _ode_solver;
+
+ static bool odeRhs(const double /*t*/,
+ MathLib::ODE::MappedConstVector<1> const y,
+ MathLib::ODE::MappedVector<1>
+ ydot,
+ React& reaction);
+};
+
+} // namespace TES
+} // namespace ProcessLib
+#endif // PROCESSLIB_TES_TESREACTIONADAPTOR_H
--
GitLab