OGS
TESProcess.cpp
Go to the documentation of this file.
1
11#include "TESProcess.h"
12
15
16namespace ProcessLib
17{
18namespace TES
19{
21 std::string name,
22 MeshLib::Mesh& mesh,
23 std::unique_ptr<AbstractJacobianAssembler>&& jacobian_assembler,
24 std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
25 unsigned const integration_order,
26 std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
27 process_variables,
28 SecondaryVariableCollection&& secondary_variables,
29 const BaseLib::ConfigTree& config)
30 : Process(std::move(name), mesh, std::move(jacobian_assembler), parameters,
31 integration_order, std::move(process_variables),
32 std::move(secondary_variables))
33{
34 DBUG("Create TESProcess.");
35
36 // physical parameters for local assembly
37 {
38 std::vector<std::pair<std::string, double*>> params{
40 {"fluid_specific_heat_source",
43 {"fluid_specific_isobaric_heat_capacity", &_assembly_params.cpG},
45 {"solid_specific_heat_source",
48 {"solid_heat_conductivity", &_assembly_params.solid_heat_cond},
50 {"solid_specific_isobaric_heat_capacity", &_assembly_params.cpS},
52 {"tortuosity", &_assembly_params.tortuosity},
54 {"diffusion_coefficient",
57 {"porosity", &_assembly_params.poro},
59 {"solid_density_dry", &_assembly_params.rho_SR_dry},
61 {"solid_density_initial", &_assembly_params.initial_solid_density}};
62
63 for (auto const& p : params)
64 {
65 if (auto const par =
67 config.getConfigParameterOptional<double>(p.first))
68 {
69 DBUG("setting parameter `{:s}' to value `{:g}'", p.first, *par);
70 *p.second = *par;
71 }
72 }
73 }
74
75 // characteristic values of primary variables
76 {
77 std::vector<std::pair<std::string, Trafo*>> const params{
79 {"characteristic_pressure", &_assembly_params.trafo_p},
81 {"characteristic_temperature", &_assembly_params.trafo_T},
83 {"characteristic_vapour_mass_fraction", &_assembly_params.trafo_x}};
84
85 for (auto const& p : params)
86 {
87 if (auto const par =
89 config.getConfigParameterOptional<double>(p.first))
90 {
91 INFO("setting parameter `{:s}' to value `{:g}'", p.first, *par);
92 *p.second = Trafo{*par};
93 }
94 }
95 }
96
97 // permeability
98 if (auto par =
100 config.getConfigParameterOptional<double>(
101 "solid_hydraulic_permeability"))
102 {
103 DBUG(
104 "setting parameter `solid_hydraulic_permeability' to isotropic "
105 "value `{:g}'",
106 *par);
107 const auto dim = mesh.getDimension();
109 Eigen::MatrixXd::Identity(dim, dim) * (*par);
110 }
111
112 // reactive system
115 config.getConfigSubtree("reactive_system"));
116
117 // debug output
118 if (auto const param =
120 config.getConfigParameterOptional<bool>("output_element_matrices"))
121 {
122 DBUG("output_element_matrices: {:s}", (*param) ? "true" : "false");
123
125 }
126
127 // TODO somewhere else
128 /*
129 if (auto const param =
131 config.getConfigParameterOptional<bool>("output_global_matrix"))
132 {
133 DBUG("output_global_matrix: {:s}", (*param) ? "true" : "false");
134
135 this->_process_output.output_global_matrix = *param;
136 }
137 */
138}
139
141 NumLib::LocalToGlobalIndexMap const& dof_table, MeshLib::Mesh const& mesh,
142 unsigned const integration_order)
143{
144 ProcessLib::createLocalAssemblers<TESLocalAssembler>(
145 mesh.getDimension(), mesh.getElements(), dof_table, _local_assemblers,
146 NumLib::IntegrationOrder{integration_order}, mesh.isAxiallySymmetric(),
148
150}
151
153{
154 // adds a secondary variables to the collection of all secondary variables.
155 auto add2nd =
156 [&](std::string const& var_name, SecondaryVariableFunctions&& fcts)
157 { _secondary_variables.addSecondaryVariable(var_name, std::move(fcts)); };
158
159 // creates an extrapolator
160 auto makeEx =
161 [&](unsigned const n_components,
162 std::vector<double> const& (TESLocalAssemblerInterface::*method)(
163 const double /*t*/,
164 std::vector<GlobalVector*> const& /*x*/,
165 std::vector<
166 NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
167 std::vector<double>& /*cache*/)
169 {
170 return ProcessLib::makeExtrapolator(n_components, getExtrapolator(),
171 _local_assemblers, method);
172 };
173
174 add2nd("solid_density",
176
177 add2nd("reaction_rate",
179
180 add2nd("darcy_velocity",
181 makeEx(_mesh.getDimension(),
183
184 add2nd("loading", makeEx(1, &TESLocalAssemblerInterface::getIntPtLoading));
185 add2nd(
186 "reaction_damping_factor",
188
189 add2nd("vapour_partial_pressure",
190 {1,
191 [&](auto&&... args) -> GlobalVector const&
192 { return computeVapourPartialPressure(args...); },
193 nullptr});
194 add2nd("relative_humidity",
195 {1,
196 [&](auto&&... args) -> GlobalVector const&
197 { return computeRelativeHumidity(args...); },
198 nullptr});
199 add2nd("equilibrium_loading",
200 {1,
201 [&](auto&&... args) -> GlobalVector const&
202 { return computeEquilibriumLoading(args...); },
203 nullptr});
204}
205
207 const double t, double const dt, std::vector<GlobalVector*> const& x,
208 std::vector<GlobalVector*> const& x_prev, int const process_id,
210{
211 DBUG("Assemble TESProcess.");
212
213 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_table = {
215
216 // Call global assembler for each local assembly item.
219 getActiveElementIDs(), dof_table, t, dt, x, x_prev, process_id, &M, &K,
220 &b);
221}
222
224 const double t, double const dt, std::vector<GlobalVector*> const& x,
225 std::vector<GlobalVector*> const& x_prev, int const process_id,
226 GlobalVector& b, GlobalMatrix& Jac)
227{
228 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_table = {
230
231 // Call global assembler for each local assembly item.
234 _local_assemblers, getActiveElementIDs(), dof_table, t, dt, x, x_prev,
235 process_id, &b, &Jac);
236}
237
238void TESProcess::preTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
239 const double t,
240 const double delta_t,
241 const int process_id)
242{
243 DBUG("new timestep");
244
245 _assembly_params.delta_t = delta_t;
247 ++_assembly_params.timestep; // TODO remove that
248
251}
252
260
262 GlobalVector const& x)
263{
264 bool check_passed = true;
265
267 {
268 // bounds checking only has to happen if the vapour mass fraction is
269 // non-logarithmic.
270
271 std::vector<GlobalIndexType> indices_cache;
272 std::vector<double> local_x_cache;
273 std::vector<double> local_x_prev_ts_cache;
274
276
277 auto check_variable_bounds =
278 [&](std::size_t id, TESLocalAssemblerInterface& loc_asm)
279 {
280 auto const r_c_indices = NumLib::getRowColumnIndices(
281 id, *this->_local_to_global_index_map, indices_cache);
282 local_x_cache = x.get(r_c_indices.rows);
283 local_x_prev_ts_cache = _x_previous_timestep->get(r_c_indices.rows);
284
285 if (!loc_asm.checkBounds(local_x_cache, local_x_prev_ts_cache))
286 {
287 check_passed = false;
288 }
289 };
290
291 GlobalExecutor::executeDereferenced(check_variable_bounds,
293 }
294
295 if (!check_passed)
296 {
298 }
299
300 // TODO remove
301 DBUG("ts {:d} iteration {:d} (in current ts: {:d}) try {:d} accepted",
305
307
309}
310
312 const double /*t*/,
313 std::vector<GlobalVector*> const& x,
314 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
315 std::unique_ptr<GlobalVector>& result_cache)
316{
317 constexpr int process_id = 0; // monolithic scheme.
318 assert(dof_table[process_id] == _local_to_global_index_map.get());
319
320 auto const& dof_table_single = getSingleComponentDOFTable();
322 {dof_table_single.dofSizeWithoutGhosts(),
323 dof_table_single.dofSizeWithoutGhosts(),
324 &dof_table_single.getGhostIndices(), nullptr});
325
326 GlobalIndexType const nnodes = _mesh.getNumberOfNodes();
327
328 for (GlobalIndexType node_id = 0; node_id < nnodes; ++node_id)
329 {
330 auto const p =
331 NumLib::getNodalValue(*x[process_id], _mesh, *dof_table[process_id],
332 node_id, COMPONENT_ID_PRESSURE);
333 auto const x_mV =
334 NumLib::getNodalValue(*x[process_id], _mesh, *dof_table[process_id],
336
339
340 result_cache->set(node_id, p * x_nV);
341 }
342
343 return *result_cache;
344}
345
347 double const /*t*/,
348 std::vector<GlobalVector*> const& xs,
349 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
350 std::unique_ptr<GlobalVector>& result_cache)
351{
352 constexpr int process_id = 0; // monolithic scheme.
353 assert(dof_table[process_id] == _local_to_global_index_map.get());
354
355 auto const& dof_table_single = getSingleComponentDOFTable();
357 {dof_table_single.dofSizeWithoutGhosts(),
358 dof_table_single.dofSizeWithoutGhosts(),
359 &dof_table_single.getGhostIndices(), nullptr});
360
361 GlobalIndexType const nnodes = _mesh.getNumberOfNodes();
362
363 auto const& x = *xs[0]; // monolithic process
364 for (GlobalIndexType node_id = 0; node_id < nnodes; ++node_id)
365 {
366 auto const p = NumLib::getNodalValue(x, _mesh, *dof_table[process_id],
367 node_id, COMPONENT_ID_PRESSURE);
368 auto const T = NumLib::getNodalValue(x, _mesh, *dof_table[process_id],
369 node_id, COMPONENT_ID_TEMPERATURE);
370 auto const x_mV =
371 NumLib::getNodalValue(x, _mesh, *dof_table[process_id], node_id,
373
376
377 auto const p_S =
379
380 result_cache->set(node_id, p * x_nV / p_S);
381 }
382
383 return *result_cache;
384}
385
387 double const /*t*/,
388 std::vector<GlobalVector*> const& xs,
389 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
390 std::unique_ptr<GlobalVector>& result_cache)
391{
392 constexpr int process_id = 0; // monolithic scheme.
393 assert(dof_table[process_id] == _local_to_global_index_map.get());
394
395 auto const& dof_table_single = getSingleComponentDOFTable();
397 {dof_table_single.dofSizeWithoutGhosts(),
398 dof_table_single.dofSizeWithoutGhosts(),
399 &dof_table_single.getGhostIndices(), nullptr});
400
401 GlobalIndexType const nnodes = _mesh.getNumberOfNodes();
402
403 auto const& x = *xs[0]; // monolithic process
404 for (GlobalIndexType node_id = 0; node_id < nnodes; ++node_id)
405 {
406 auto const p = NumLib::getNodalValue(x, _mesh, *dof_table[process_id],
407 node_id, COMPONENT_ID_PRESSURE);
408 auto const T = NumLib::getNodalValue(x, _mesh, *dof_table[process_id],
409 node_id, COMPONENT_ID_TEMPERATURE);
410 auto const x_mV =
411 NumLib::getNodalValue(x, _mesh, *dof_table[process_id], node_id,
413
416
417 auto const p_V = p * x_nV;
418 auto const C_eq =
419 (p_V <= 0.0) ? 0.0
420 : _assembly_params.react_sys->getEquilibriumLoading(
421 p_V, T, _assembly_params.M_react);
422
423 result_cache->set(node_id, C_eq);
424 }
425
426 return *result_cache;
427}
428
429} // namespace TES
430
431} // namespace ProcessLib
GlobalMatrix::IndexType GlobalIndexType
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:35
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:30
static double getMolarFraction(double xm, double M_this, double M_other)
static double getEquilibriumVapourPressure(const double T_Ads)
static std::unique_ptr< Reaction > newInstance(BaseLib::ConfigTree const &conf)
Definition Reaction.cpp:27
std::optional< T > getConfigParameterOptional(std::string const &param) const
ConfigTree getConfigSubtree(std::string const &root) const
Global vector based on Eigen vector.
Definition EigenVector.h:25
double get(IndexType rowId) const
get entry
Definition EigenVector.h:58
bool isAxiallySymmetric() const
Definition Mesh.h:137
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition Mesh.h:109
unsigned getDimension() const
Returns the dimension of the mesh (determined by the maximum dimension over all elements).
Definition Mesh.h:88
std::size_t getNumberOfNodes() const
Get the number of nodes.
Definition Mesh.h:100
MeshLib::Mesh & _mesh
Definition Process.h:365
std::vector< std::size_t > const & getActiveElementIDs() const
Definition Process.h:167
NumLib::LocalToGlobalIndexMap const & getSingleComponentDOFTable() const
Definition Process.h:213
SecondaryVariableCollection _secondary_variables
Definition Process.h:370
VectorMatrixAssembler _global_assembler
Definition Process.h:377
std::unique_ptr< NumLib::LocalToGlobalIndexMap > _local_to_global_index_map
Definition Process.h:368
NumLib::Extrapolator & getExtrapolator() const
Definition Process.h:208
Handles configuration of several secondary variables from the project file.
void addSecondaryVariable(std::string const &internal_name, SecondaryVariableFunctions &&fcts)
virtual std::vector< double > const & getIntPtSolidDensity(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
virtual std::vector< double > const & getIntPtLoading(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
virtual std::vector< double > const & getIntPtReactionDampingFactor(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
virtual std::vector< double > const & getIntPtDarcyVelocity(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
virtual std::vector< double > const & getIntPtReactionRate(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
std::vector< std::unique_ptr< TESLocalAssemblerInterface > > _local_assemblers
Definition TESProcess.h:91
TESProcess(std::string name, MeshLib::Mesh &mesh, std::unique_ptr< AbstractJacobianAssembler > &&jacobian_assembler, std::vector< std::unique_ptr< ParameterLib::ParameterBase > > const &parameters, unsigned const integration_order, std::vector< std::vector< std::reference_wrapper< ProcessVariable > > > &&process_variables, SecondaryVariableCollection &&secondary_variables, BaseLib::ConfigTree const &config)
GlobalVector const & computeRelativeHumidity(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::unique_ptr< GlobalVector > &result_cache)
void preTimestepConcreteProcess(std::vector< GlobalVector * > const &x, const double t, const double delta_t, const int process_id) override
NumLib::IterationResult postIterationConcreteProcess(GlobalVector const &x) override
std::unique_ptr< GlobalVector > _x_previous_timestep
Definition TESProcess.h:96
GlobalVector const & computeEquilibriumLoading(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::unique_ptr< GlobalVector > &result_cache)
void assembleConcreteProcess(const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b) override
void assembleWithJacobianConcreteProcess(const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id, GlobalVector &b, GlobalMatrix &Jac) override
GlobalVector const & computeVapourPartialPressure(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::unique_ptr< GlobalVector > &result_cache)
void initializeConcreteProcess(NumLib::LocalToGlobalIndexMap const &dof_table, MeshLib::Mesh const &mesh, unsigned const integration_order) override
Process specific initialization called by initialize().
void preIterationConcreteProcess(const unsigned iter, GlobalVector const &x) override
AssemblyParams _assembly_params
Definition TESProcess.h:93
void assemble(std::size_t const mesh_item_id, LocalAssemblerInterface &local_assembler, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, double const t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id, GlobalMatrix *M, GlobalMatrix *K, GlobalVector *b)
void assembleWithJacobian(std::size_t const mesh_item_id, LocalAssemblerInterface &local_assembler, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id, GlobalVector *b, GlobalMatrix *Jac)
IterationResult
Status flags telling the NonlinearSolver if an iteration succeeded.
void setLocalAccessibleVector(PETScVector const &x)
Definition LinAlg.cpp:27
double getNodalValue(GlobalVector const &x, MeshLib::Mesh const &mesh, NumLib::LocalToGlobalIndexMap const &dof_table, std::size_t const node_id, std::size_t const global_component_id)
NumLib::LocalToGlobalIndexMap::RowColumnIndices getRowColumnIndices(std::size_t const id, NumLib::LocalToGlobalIndexMap const &dof_table, std::vector< GlobalIndexType > &indices)
const unsigned COMPONENT_ID_MASS_FRACTION
const unsigned COMPONENT_ID_TEMPERATURE
const unsigned COMPONENT_ID_PRESSURE
SecondaryVariableFunctions makeExtrapolator(const unsigned num_components, NumLib::Extrapolator &extrapolator, LocalAssemblerCollection const &local_assemblers, typename NumLib::ExtrapolatableLocalAssemblerCollection< LocalAssemblerCollection >::IntegrationPointValuesMethod integration_point_values_method)
static void executeSelectedMemberDereferenced(Object &object, Method method, Container const &container, std::vector< std::size_t > const &active_container_ids, Args &&... args)
static void executeDereferenced(F const &f, C const &c, Args_ &&... args)
std::size_t timestep
Output global matrix/rhs after first iteration.
std::unique_ptr< Adsorption::Reaction > react_sys