From 4639a1733bf267574e29b8c38a6523475119cdbe Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <dmitri.naumov@ufz.de>
Date: Mon, 22 Jan 2024 11:39:36 +0100
Subject: [PATCH] [PL/CT] Optional temperature variable.
---
.../ComponentTransportFEM.h | 86 ++++++++++++++-----
.../ComponentTransportProcess.cpp | 42 ++++++---
.../ComponentTransportProcessData.h | 6 +-
.../CreateComponentTransportProcess.cpp | 47 ++++++++--
4 files changed, 137 insertions(+), 44 deletions(-)
diff --git a/ProcessLib/ComponentTransport/ComponentTransportFEM.h b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
index 1982fa9d57a..215cbfc9748 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportFEM.h
+++ b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
@@ -210,7 +210,8 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface
// When monolithic scheme is adopted, nodal pressure and nodal concentration
// are accessed by vector index.
static const int pressure_index = 0;
- static const int first_concentration_index = ShapeFunction::NPOINTS;
+ const int temperature_index = -1;
+ const int first_concentration_index = -1;
static const int pressure_size = ShapeFunction::NPOINTS;
static const int temperature_size = ShapeFunction::NPOINTS;
@@ -247,7 +248,12 @@ public:
ComponentTransportProcessData const& process_data,
std::vector<std::reference_wrapper<ProcessVariable>> const&
transport_process_variables)
- : _element(element),
+ : temperature_index(process_data.isothermal ? -1
+ : ShapeFunction::NPOINTS),
+ first_concentration_index(process_data.isothermal
+ ? ShapeFunction::NPOINTS
+ : 2 * ShapeFunction::NPOINTS),
+ _element(element),
_process_data(process_data),
_integration_method(integration_method),
_transport_process_variables(transport_process_variables)
@@ -811,6 +817,25 @@ public:
auto const local_C_prev =
local_x_prev.segment<concentration_size>(first_concentration_index);
+ NodalVectorType local_T;
+ if (_process_data.isothermal)
+ {
+ if (_process_data.temperature)
+ {
+ local_T = _process_data.temperature->getNodalValuesOnElement(
+ _element, t);
+ }
+ else
+ {
+ local_T = NodalVectorType::Zero(temperature_size);
+ }
+ }
+ else
+ {
+ local_T =
+ local_x.template segment<temperature_size>(temperature_index);
+ }
+
auto local_M = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
local_M_data, pressure_size, pressure_size);
auto local_K = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
@@ -846,14 +871,13 @@ public:
auto& porosity = ip_data.porosity;
auto const& porosity_prev = ip_data.porosity_prev;
- double C_int_pt = 0.0;
- double p_int_pt = 0.0;
-
- NumLib::shapeFunctionInterpolate(local_C, N, C_int_pt);
- NumLib::shapeFunctionInterpolate(local_p, N, p_int_pt);
+ double const C_int_pt = N.dot(local_C);
+ double const p_int_pt = N.dot(local_p);
+ double const T_int_pt = N.dot(local_T);
vars.concentration = C_int_pt;
vars.liquid_phase_pressure = p_int_pt;
+ vars.temperature = T_int_pt;
// porosity
{
@@ -1048,20 +1072,38 @@ public:
std::vector<double>& local_K_data,
std::vector<double>& /*local_b_data*/, int const transport_process_id)
{
+ assert(static_cast<int>(local_x.size()) ==
+ pressure_size +
+ concentration_size *
+ static_cast<int>(_transport_process_variables.size()) +
+ (_process_data.isothermal ? 0 : temperature_size));
+
auto const local_p =
local_x.template segment<pressure_size>(pressure_index);
- auto const local_C = local_x.template segment<concentration_size>(
- first_concentration_index +
- (transport_process_id - 1) * concentration_size);
- auto const local_p_prev =
- local_x_prev.segment<pressure_size>(pressure_index);
-
NodalVectorType local_T;
- if (_process_data.temperature)
+ if (_process_data.isothermal)
+ {
+ if (_process_data.temperature)
+ {
+ local_T = _process_data.temperature->getNodalValuesOnElement(
+ _element, t);
+ }
+ else
+ {
+ local_T = NodalVectorType::Zero(temperature_size);
+ }
+ }
+ else
{
local_T =
- _process_data.temperature->getNodalValuesOnElement(_element, t);
+ local_x.template segment<temperature_size>(temperature_index);
}
+ auto const local_C = local_x.template segment<concentration_size>(
+ first_concentration_index +
+ (transport_process_id - (_process_data.isothermal ? 1 : 2)) *
+ concentration_size);
+ auto const local_p_prev =
+ local_x_prev.segment<pressure_size>(pressure_index);
auto local_M = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
local_M_data, concentration_size, concentration_size);
@@ -1094,9 +1136,8 @@ public:
auto const& medium =
*_process_data.media_map.getMedium(_element.getID());
auto const& phase = medium.phase("AqueousLiquid");
- // Hydraulic process id is 0 and thus transport process id starts
- // from 1.
- auto const component_id = transport_process_id - 1;
+ auto const component_id =
+ transport_process_id - (_process_data.isothermal ? 1 : 2);
auto const& component = phase.component(
_transport_process_variables[component_id].get().getName());
@@ -1111,14 +1152,13 @@ public:
auto& porosity = ip_data.porosity;
auto const& porosity_prev = ip_data.porosity_prev;
- double C_int_pt = 0.0;
- double p_int_pt = 0.0;
-
- NumLib::shapeFunctionInterpolate(local_C, N, C_int_pt);
- NumLib::shapeFunctionInterpolate(local_p, N, p_int_pt);
+ double const C_int_pt = N.dot(local_C);
+ double const p_int_pt = N.dot(local_p);
+ double const T_int_pt = N.dot(local_T);
vars.concentration = C_int_pt;
vars.liquid_phase_pressure = p_int_pt;
+ vars.temperature = T_int_pt;
if (_process_data.temperature)
{
diff --git a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
index 5e816c4a883..6a15ad753ba 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
+++ b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
@@ -11,6 +11,7 @@
#include "ComponentTransportProcess.h"
#include <cassert>
+#include <range/v3/algorithm/copy.hpp>
#include <range/v3/view/drop.hpp>
#include "BaseLib/RunTime.h"
@@ -88,27 +89,36 @@ ComponentTransportProcess::ComponentTransportProcess(
"probably expect.");
}
- _residua.push_back(MeshLib::getOrCreateMeshProperty<double>(
- mesh, "LiquidMassFlowRate", MeshLib::MeshItemType::Node, 1));
+ auto residuum_name = [&](auto const& pv) -> std::string
+ {
+ std::string const& pv_name = pv.getName();
+ if (pv_name == "pressure")
+ {
+ return "LiquidMassFlowRate";
+ }
+ if (pv_name == "temperature")
+ {
+ return "HeatFlowRate";
+ }
+ return pv_name + "FlowRate";
+ };
if (_use_monolithic_scheme)
{
- const int process_id = 0;
- for (auto const& pv :
- _process_variables[process_id] | ranges::views::drop(1))
+ int const process_id = 0;
+ for (auto const& pv : _process_variables[process_id])
{
_residua.push_back(MeshLib::getOrCreateMeshProperty<double>(
- mesh, pv.get().getName() + "FlowRate",
- MeshLib::MeshItemType::Node, 1));
+ mesh, residuum_name(pv.get()), MeshLib::MeshItemType::Node, 1));
}
}
else
{
- for (auto const& pv : _process_variables | ranges::views::drop(1))
+ for (auto const& pv : _process_variables)
{
_residua.push_back(MeshLib::getOrCreateMeshProperty<double>(
- mesh, pv[0].get().getName() + "FlowRate",
- MeshLib::MeshItemType::Node, 1));
+ mesh, residuum_name(pv[0].get()), MeshLib::MeshItemType::Node,
+ 1));
}
}
}
@@ -142,11 +152,13 @@ void ComponentTransportProcess::initializeConcreteProcess(
}
else
{
- for (auto pv_iter = std::next(_process_variables.begin());
- pv_iter != _process_variables.end();
- ++pv_iter)
+ // All process variables but the pressure and optionally the temperature
+ // are transport variables.
+ for (auto const& pv :
+ _process_variables |
+ ranges::views::drop(_process_data.isothermal ? 1 : 2))
{
- transport_process_variables.push_back((*pv_iter)[0]);
+ transport_process_variables.push_back(pv[0]);
}
}
@@ -232,6 +244,8 @@ void ComponentTransportProcess::assembleConcreteProcess(
_global_assembler, _local_assemblers,
pv.getActiveElementIDs());
+ // TODO (naumov) What about temperature? A test with FCT would reveal any
+ // problems.
if (process_id == _process_data.hydraulic_process_id)
{
return;
diff --git a/ProcessLib/ComponentTransport/ComponentTransportProcessData.h b/ProcessLib/ComponentTransport/ComponentTransportProcessData.h
index 3aa65da281d..2ac5512eed9 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportProcessData.h
+++ b/ProcessLib/ComponentTransport/ComponentTransportProcessData.h
@@ -82,13 +82,17 @@ struct ComponentTransportProcessData
ParameterLib::Parameter<double> const& aperture_size =
ParameterLib::ConstantParameter<double>("constant_one", 1.0);
+ bool const isothermal;
+
static const int hydraulic_process_id = 0;
+ // Thermal process is optional, indicated by -1. If present, it is positive.
+ const int thermal_process_id = isothermal ? -1 : 1;
// TODO (renchao-lu): This variable is used in the calculation of the
// fluid's density and flux, indicating the transport process id. For now it
// is assumed that these quantities depend on the first occurring transport
// process only. The density and flux calculations have to be extended to
// all processes.
- static const int first_transport_process_id = 1;
+ const int first_transport_process_id = isothermal ? 1 : 2;
MeshLib::PropertyVector<double>* mesh_prop_velocity = nullptr;
MeshLib::PropertyVector<double>* mesh_prop_porosity = nullptr;
diff --git a/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp b/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp
index 2aea148a7bd..bb69b4fe94e 100644
--- a/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp
+++ b/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp
@@ -107,14 +107,27 @@ std::unique_ptr<Process> createComponentTransportProcess(
std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>
process_variables;
- /// Primary process variables as they appear in the global component vector:
- auto const collected_process_variables = findProcessVariables(
+ /// Primary process variables as they appear in the global component vector.
+ std::vector collected_process_variables = findProcessVariables(
variables, pv_config,
{//! \ogs_file_param_special{prj__processes__process__ComponentTransport__process_variables__pressure}
"pressure",
//! \ogs_file_param_special{prj__processes__process__ComponentTransport__process_variables__concentration}
"concentration"});
+ /// Temperature is optional.
+ std::vector const temperature_variable = findProcessVariables(
+ variables, pv_config,
+ //! \ogs_file_param_special{prj__processes__process__ComponentTransport__process_variables__temperature}
+ "temperature", true /*temperature is optional*/);
+ bool const isothermal = temperature_variable.empty();
+ if (!isothermal)
+ {
+ assert(temperature_variable.size() == 1);
+ collected_process_variables.insert(
+ ++collected_process_variables.begin(), temperature_variable[0]);
+ }
+
// Check number of components for each process variable
auto it = std::find_if(
collected_process_variables.cbegin(),
@@ -136,6 +149,13 @@ std::unique_ptr<Process> createComponentTransportProcess(
// depending on what scheme is adopted
if (use_monolithic_scheme) // monolithic scheme.
{
+ if (!isothermal)
+ {
+ OGS_FATAL(
+ "Currently, non-isothermal component transport process can "
+ "only be simulated in staggered scheme.");
+ }
+
process_variables.push_back(std::move(collected_process_variables));
}
else // staggered scheme.
@@ -183,8 +203,16 @@ std::unique_ptr<Process> createComponentTransportProcess(
// pressure
per_process_variable.emplace_back(collected_process_variables[0]);
process_variables.push_back(std::move(per_process_variable));
+ // temperature
+ if (!isothermal)
+ {
+ per_process_variable.emplace_back(
+ collected_process_variables[1]);
+ process_variables.push_back(std::move(per_process_variable));
+ }
// concentration
- assert(components.size() + 1 == collected_process_variables.size());
+ assert(components.size() + (isothermal ? 1 : 2) ==
+ collected_process_variables.size());
std::transform(components.begin(), components.end(),
std::back_inserter(process_variables),
sort_by_component);
@@ -222,9 +250,14 @@ std::unique_ptr<Process> createComponentTransportProcess(
config.getConfigParameter<bool>("chemically_induced_porosity_change",
false);
- auto const temperature = ParameterLib::findOptionalTagParameter<double>(
+ auto const temperature_field = ParameterLib::findOptionalTagParameter<
+ double>(
//! \ogs_file_param_special{prj__processes__process__ComponentTransport__temperature_field}
config, "temperature_field", parameters, 1, &mesh);
+ if (!isothermal && temperature_field != nullptr)
+ {
+ OGS_FATAL("Temperature field is set for non-isothermal setup.")
+ }
auto media_map =
MaterialPropertyLib::createMaterialSpatialDistributionMap(media, mesh);
@@ -275,14 +308,16 @@ std::unique_ptr<Process> createComponentTransportProcess(
std::move(media_map),
has_gravity,
non_advective_form,
- temperature,
+ temperature_field,
chemically_induced_porosity_change,
chemical_solver_interface.get(),
std::move(lookup_table),
std::move(stabilizer),
projected_specific_body_force_vectors,
mesh_space_dimension,
- *aperture_size_parameter};
+ *aperture_size_parameter,
+ isothermal,
+ };
SecondaryVariableCollection secondary_variables;
--
GitLab