OGS
TESProcess.cpp
Go to the documentation of this file.
1 
11 #include "TESProcess.h"
12 
15 
16 namespace ProcessLib
17 {
18 namespace 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 
206 void TESProcess::assembleConcreteProcess(const double t, double const dt,
207  std::vector<GlobalVector*> const& x,
208  std::vector<GlobalVector*> const& xdot,
209  int const process_id, GlobalMatrix& M,
210  GlobalMatrix& K, GlobalVector& b)
211 {
212  DBUG("Assemble TESProcess.");
213 
214  std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
215  dof_table = {std::ref(*_local_to_global_index_map)};
216  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
217 
218  // Call global assembler for each local assembly item.
221  pv.getActiveElementIDs(), dof_table, t, dt, x, xdot, process_id, M, K,
222  b);
223 }
224 
226  const double t, double const dt, std::vector<GlobalVector*> const& x,
227  std::vector<GlobalVector*> const& xdot, int const process_id,
229 {
230  std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
231  dof_table = {std::ref(*_local_to_global_index_map)};
232  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
233 
234  // Call global assembler for each local assembly item.
237  _local_assemblers, pv.getActiveElementIDs(), dof_table, t, dt, x, xdot,
238  process_id, M, K, b, Jac);
239 }
240 
241 void TESProcess::preTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
242  const double t,
243  const double delta_t,
244  const int process_id)
245 {
246  DBUG("new timestep");
247 
248  _assembly_params.delta_t = delta_t;
250  ++_assembly_params.timestep; // TODO remove that
251 
254 }
255 
257  GlobalVector const& /*x*/)
258 {
262 }
263 
265  GlobalVector const& x)
266 {
267  bool check_passed = true;
268 
269  if (!Trafo::constrained)
270  {
271  // bounds checking only has to happen if the vapour mass fraction is
272  // non-logarithmic.
273 
274  std::vector<GlobalIndexType> indices_cache;
275  std::vector<double> local_x_cache;
276  std::vector<double> local_x_prev_ts_cache;
277 
279 
280  auto check_variable_bounds =
281  [&](std::size_t id, TESLocalAssemblerInterface& loc_asm)
282  {
283  auto const r_c_indices = NumLib::getRowColumnIndices(
284  id, *this->_local_to_global_index_map, indices_cache);
285  local_x_cache = x.get(r_c_indices.rows);
286  local_x_prev_ts_cache = _x_previous_timestep->get(r_c_indices.rows);
287 
288  if (!loc_asm.checkBounds(local_x_cache, local_x_prev_ts_cache))
289  {
290  check_passed = false;
291  }
292  };
293 
294  GlobalExecutor::executeDereferenced(check_variable_bounds,
296  }
297 
298  if (!check_passed)
299  {
301  }
302 
303  // TODO remove
304  DBUG("ts {:d} iteration {:d} (in current ts: {:d}) try {:d} accepted",
308 
310 
312 }
313 
315  const double /*t*/,
316  std::vector<GlobalVector*> const& x,
317  std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
318  std::unique_ptr<GlobalVector>& result_cache)
319 {
320  constexpr int process_id = 0; // monolithic scheme.
321  assert(dof_table[process_id] == _local_to_global_index_map.get());
322 
323  auto const& dof_table_single = getSingleComponentDOFTable();
325  {dof_table_single.dofSizeWithoutGhosts(),
326  dof_table_single.dofSizeWithoutGhosts(),
327  &dof_table_single.getGhostIndices(), nullptr});
328 
329  GlobalIndexType const nnodes = _mesh.getNumberOfNodes();
330 
331  for (GlobalIndexType node_id = 0; node_id < nnodes; ++node_id)
332  {
333  auto const p =
334  NumLib::getNodalValue(*x[process_id], _mesh, *dof_table[process_id],
335  node_id, COMPONENT_ID_PRESSURE);
336  auto const x_mV =
337  NumLib::getNodalValue(*x[process_id], _mesh, *dof_table[process_id],
338  node_id, COMPONENT_ID_MASS_FRACTION);
339 
342 
343  result_cache->set(node_id, p * x_nV);
344  }
345 
346  return *result_cache;
347 }
348 
350  double const /*t*/,
351  std::vector<GlobalVector*> const& xs,
352  std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
353  std::unique_ptr<GlobalVector>& result_cache)
354 {
355  constexpr int process_id = 0; // monolithic scheme.
356  assert(dof_table[process_id] == _local_to_global_index_map.get());
357 
358  auto const& dof_table_single = getSingleComponentDOFTable();
360  {dof_table_single.dofSizeWithoutGhosts(),
361  dof_table_single.dofSizeWithoutGhosts(),
362  &dof_table_single.getGhostIndices(), nullptr});
363 
364  GlobalIndexType const nnodes = _mesh.getNumberOfNodes();
365 
366  auto const& x = *xs[0]; // monolithic process
367  for (GlobalIndexType node_id = 0; node_id < nnodes; ++node_id)
368  {
369  auto const p = NumLib::getNodalValue(x, _mesh, *dof_table[process_id],
370  node_id, COMPONENT_ID_PRESSURE);
371  auto const T = NumLib::getNodalValue(x, _mesh, *dof_table[process_id],
372  node_id, COMPONENT_ID_TEMPERATURE);
373  auto const x_mV =
374  NumLib::getNodalValue(x, _mesh, *dof_table[process_id], node_id,
376 
379 
380  auto const p_S =
382 
383  result_cache->set(node_id, p * x_nV / p_S);
384  }
385 
386  return *result_cache;
387 }
388 
390  double const /*t*/,
391  std::vector<GlobalVector*> const& xs,
392  std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
393  std::unique_ptr<GlobalVector>& result_cache)
394 {
395  constexpr int process_id = 0; // monolithic scheme.
396  assert(dof_table[process_id] == _local_to_global_index_map.get());
397 
398  auto const& dof_table_single = getSingleComponentDOFTable();
400  {dof_table_single.dofSizeWithoutGhosts(),
401  dof_table_single.dofSizeWithoutGhosts(),
402  &dof_table_single.getGhostIndices(), nullptr});
403 
404  GlobalIndexType const nnodes = _mesh.getNumberOfNodes();
405 
406  auto const& x = *xs[0]; // monolithic process
407  for (GlobalIndexType node_id = 0; node_id < nnodes; ++node_id)
408  {
409  auto const p = NumLib::getNodalValue(x, _mesh, *dof_table[process_id],
410  node_id, COMPONENT_ID_PRESSURE);
411  auto const T = NumLib::getNodalValue(x, _mesh, *dof_table[process_id],
412  node_id, COMPONENT_ID_TEMPERATURE);
413  auto const x_mV =
414  NumLib::getNodalValue(x, _mesh, *dof_table[process_id], node_id,
416 
419 
420  auto const p_V = p * x_nV;
421  auto const C_eq =
422  (p_V <= 0.0) ? 0.0
423  : _assembly_params.react_sys->getEquilibriumLoading(
424  p_V, T, _assembly_params.M_react);
425 
426  result_cache->set(node_id, C_eq);
427  }
428 
429  return *result_cache;
430 }
431 
432 } // namespace TES
433 
434 } // namespace ProcessLib
GlobalMatrix::IndexType GlobalIndexType
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
Definition: Logging.h:34
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition: Logging.h:29
static double getMolarFraction(double xm, double M_this, double M_other)
Definition: Adsorption.cpp:96
static double getEquilibriumVapourPressure(const double T_Ads)
Definition: Adsorption.cpp:40
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
Definition: ConfigTree.cpp:150
Global vector based on Eigen vector.
Definition: EigenVector.h:28
double get(IndexType rowId) const
get entry
Definition: EigenVector.h:61
bool isAxiallySymmetric() const
Definition: Mesh.h:131
unsigned getDimension() const
Returns the dimension of the mesh (determined by the maximum dimension over all elements).
Definition: Mesh.h:76
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition: Mesh.h:103
std::size_t getNumberOfNodes() const
Get the number of nodes.
Definition: Mesh.h:94
std::vector< std::size_t > const & getActiveElementIDs() const
NumLib::Extrapolator & getExtrapolator() const
Definition: Process.h:184
NumLib::LocalToGlobalIndexMap const & getSingleComponentDOFTable() const
Definition: Process.h:189
MeshLib::Mesh & _mesh
Definition: Process.h:330
SecondaryVariableCollection _secondary_variables
Definition: Process.h:335
std::vector< std::reference_wrapper< ProcessVariable > > const & getProcessVariables(const int process_id) const
Definition: Process.h:148
VectorMatrixAssembler _global_assembler
Definition: Process.h:337
std::unique_ptr< NumLib::LocalToGlobalIndexMap > _local_to_global_index_map
Definition: Process.h:333
Handles configuration of several secondary variables from the project file.
void addSecondaryVariable(std::string const &internal_name, SecondaryVariableFunctions &&fcts)
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 & 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
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 & getIntPtSolidDensity(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:93
void assembleWithJacobianConcreteProcess(const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &xdot, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b, GlobalMatrix &Jac) override
Definition: TESProcess.cpp:225
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)
Definition: TESProcess.cpp:349
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)
Definition: TESProcess.cpp:20
void preTimestepConcreteProcess(std::vector< GlobalVector * > const &x, const double t, const double delta_t, const int process_id) override
Definition: TESProcess.cpp:241
void assembleConcreteProcess(const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &xdot, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b) override
Definition: TESProcess.cpp:206
NumLib::IterationResult postIterationConcreteProcess(GlobalVector const &x) override
Definition: TESProcess.cpp:264
std::unique_ptr< GlobalVector > _x_previous_timestep
Definition: TESProcess.h:98
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)
Definition: TESProcess.cpp:389
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)
Definition: TESProcess.cpp:314
void initializeConcreteProcess(NumLib::LocalToGlobalIndexMap const &dof_table, MeshLib::Mesh const &mesh, unsigned const integration_order) override
Process specific initialization called by initialize().
Definition: TESProcess.cpp:140
void preIterationConcreteProcess(const unsigned iter, GlobalVector const &x) override
Definition: TESProcess.cpp:256
AssemblyParams _assembly_params
Definition: TESProcess.h:95
void assembleWithJacobian(std::size_t const mesh_item_id, LocalAssemblerInterface &local_assembler, std::vector< std::reference_wrapper< NumLib::LocalToGlobalIndexMap >> const &dof_tables, const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &xdot, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b, GlobalMatrix &Jac)
void assemble(std::size_t const mesh_item_id, LocalAssemblerInterface &local_assembler, std::vector< std::reference_wrapper< NumLib::LocalToGlobalIndexMap >> const &dof_tables, double const t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &xdot, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b)
IterationResult
Status flags telling the NonlinearSolver if an iteration succeeded.
void setLocalAccessibleVector(PETScVector const &x)
Definition: LinAlg.cpp:27
static const double t
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