OGS
Process.cpp
Go to the documentation of this file.
1
11#include "Process.h"
12
13#include <range/v3/algorithm/any_of.hpp>
14#include <range/v3/algorithm/set_algorithm.hpp>
15#include <range/v3/view/drop.hpp>
16#include <range/v3/view/transform.hpp>
17
22#include "ProcessVariable.h"
23
24namespace ProcessLib
25{
26const std::string Process::constant_one_parameter_name = "constant_one";
27
29 std::string name_,
30 MeshLib::Mesh& mesh,
31 std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
32 std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
33 unsigned const integration_order,
34 std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
35 process_variables,
36 SecondaryVariableCollection&& secondary_variables,
37 const bool use_monolithic_scheme)
38 : name(std::move(name_)),
39 _mesh(mesh),
40 _secondary_variables(std::move(secondary_variables)),
41 _jacobian_assembler(std::move(jacobian_assembler)),
42 _global_assembler(*_jacobian_assembler),
43 _use_monolithic_scheme(use_monolithic_scheme),
44 _integration_order(integration_order),
45 _process_variables(std::move(process_variables)),
46 _boundary_conditions(
47 [&](const std::size_t number_of_process_variables)
48 -> std::vector<BoundaryConditionCollection>
49 {
50 std::vector<BoundaryConditionCollection> pcs_BCs;
51 pcs_BCs.reserve(number_of_process_variables);
52 for (std::size_t i = 0; i < number_of_process_variables; i++)
53 {
54 pcs_BCs.emplace_back(BoundaryConditionCollection(parameters));
55 }
56 return pcs_BCs;
57 }(_process_variables.size())),
58 _source_term_collections(
59 [&](const std::size_t number_of_processes)
60 -> std::vector<SourceTermCollection>
61 {
62 std::vector<SourceTermCollection> pcs_sts;
63 pcs_sts.reserve(number_of_processes);
64 for (std::size_t i = 0; i < number_of_processes; i++)
65 {
66 pcs_sts.emplace_back(SourceTermCollection(parameters));
67 }
68 return pcs_sts;
69 }(_process_variables.size()))
70{
71}
72
74 const NumLib::LocalToGlobalIndexMap& dof_table, const int process_id,
75 std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media)
76{
77 auto const& per_process_variables = _process_variables[process_id];
78 auto& per_process_BCs = _boundary_conditions[process_id];
79
80 per_process_BCs.addBCsForProcessVariables(per_process_variables, dof_table,
81 _integration_order, *this, media);
82
83 auto& per_process_sts = _source_term_collections[process_id];
84 per_process_sts.addSourceTermsForProcessVariables(
85 per_process_variables, dof_table, _integration_order);
86}
87
89 std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media)
90{
91 // The number of processes is identical to the size of _process_variables,
92 // the vector contains variables for different processes. See the
93 // documentation of _process_variables.
94 const std::size_t number_of_processes = _process_variables.size();
95 for (std::size_t pcs_id = 0; pcs_id < number_of_processes; pcs_id++)
96 {
98 *_local_to_global_index_map, pcs_id, media);
99 }
100}
101
103 std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media)
104{
105 DBUG("Initialize process.");
106
107 DBUG("Construct dof mappings.");
109
110 DBUG("Compute sparsity pattern");
112
113 DBUG("Initialize the extrapolator");
115
118
119 DBUG("Initialize boundary conditions.");
121}
122
124 std::vector<GlobalVector*>& process_solutions,
125 std::vector<GlobalVector*> const& process_solutions_prev,
126 double const t,
127 int const process_id)
128{
129 auto& x = *process_solutions[process_id];
130 auto& x_prev = *process_solutions_prev[process_id];
131
132 // getDOFTableOfProcess can be overloaded by the specific process.
133 auto const& dof_table_of_process = getDOFTable(process_id);
134
135 auto const& per_process_variables = _process_variables[process_id];
136 for (std::size_t variable_id = 0;
137 variable_id < per_process_variables.size();
138 variable_id++)
139 {
142
143 auto const& pv = per_process_variables[variable_id];
144 DBUG("Set the initial condition of variable {:s} of process {:d}.",
145 pv.get().getName().data(), process_id);
146
147 auto const& ic = pv.get().getInitialCondition();
148
149 auto const num_comp = pv.get().getNumberOfGlobalComponents();
150
151 for (int component_id = 0; component_id < num_comp; ++component_id)
152 {
153 auto const& mesh_subset =
154 dof_table_of_process.getMeshSubset(variable_id, component_id);
155 auto const mesh_id = mesh_subset.getMeshID();
156 for (auto const* node : mesh_subset.getNodes())
157 {
159 node->getID());
160
161 pos.setNodeID(node->getID());
162 pos.setCoordinates(*node);
163 auto const& ic_value = ic(t, pos);
164
165 auto global_index =
166 std::abs(dof_table_of_process.getGlobalIndex(l, variable_id,
167 component_id));
168#ifdef USE_PETSC
169 // The global indices of the ghost entries of the global
170 // matrix or the global vectors need to be set as negative
171 // values for equation assembly, however the global indices
172 // start from zero. Therefore, any ghost entry with zero
173 // index is assigned an negative value of the vector size
174 // or the matrix dimension. To assign the initial value for
175 // the ghost entries, the negative indices of the ghost
176 // entries are restored to zero.
177 if (global_index == x.size())
178 global_index = 0;
179#endif
180 x.set(global_index, ic_value[component_id]);
181 }
182 }
183 }
184
186 MathLib::LinAlg::copy(x, x_prev); // pushState
187
190
191 setInitialConditionsConcreteProcess(process_solutions, t, process_id);
192}
193
195 const int /*process_id*/) const
196{
197 auto const& l = *_local_to_global_index_map;
198 return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
199 &l.getGhostIndices(), &_sparsity_pattern};
200}
201
203 const int process_id)
204{
205 auto const& variables_per_process = getProcessVariables(process_id);
206 for (auto const& variable : variables_per_process)
207 {
208 variable.get().updateDeactivatedSubdomains(time);
209 }
210
212
213 auto active_elements_ids = ranges::views::transform(
214 [](auto const& variable)
215 { return variable.get().getActiveElementIDs(); });
216
217 // Early return if there's any process variable with all elements active.
218 if (ranges::any_of(variables_per_process | active_elements_ids,
219 [](auto const& vector) { return vector.empty(); }))
220 {
221 return;
222 }
223
224 // Some process variable has deactivated elements. Create union of active
225 // ids.
226
227 _ids_of_active_elements = // there is at least one process variable.
228 variables_per_process[0].get().getActiveElementIDs();
229
230 for (auto const& pv_active_element_ids :
231 variables_per_process | ranges::views::drop(1) | active_elements_ids)
232 {
233 std::vector<std::size_t> new_active_elements;
234 new_active_elements.reserve(_ids_of_active_elements.size() +
235 pv_active_element_ids.size());
236 ranges::set_union(_ids_of_active_elements, pv_active_element_ids,
237 std::back_inserter(new_active_elements));
238
239 _ids_of_active_elements = std::move(new_active_elements);
240 }
241}
242
243void Process::preAssemble(const double t, double const dt,
244 GlobalVector const& x)
245{
247}
248
249void Process::assemble(const double t, double const dt,
250 std::vector<GlobalVector*> const& x,
251 std::vector<GlobalVector*> const& x_prev,
252 int const process_id, GlobalMatrix& M, GlobalMatrix& K,
253 GlobalVector& b)
254{
255 assert(x.size() == x_prev.size());
256 for (std::size_t i = 0; i < x.size(); i++)
257 {
260 }
261
262 assembleConcreteProcess(t, dt, x, x_prev, process_id, M, K, b);
263
264 // the last argument is for the jacobian, nullptr is for a unused jacobian
265 _boundary_conditions[process_id].applyNaturalBC(t, x, process_id, K, b,
266 nullptr);
267
268 // the last argument is for the jacobian, nullptr is for a unused jacobian
269 _source_term_collections[process_id].integrate(t, *x[process_id], b,
270 nullptr);
271}
272
273void Process::assembleWithJacobian(const double t, double const dt,
274 std::vector<GlobalVector*> const& x,
275 std::vector<GlobalVector*> const& x_prev,
276 int const process_id, GlobalMatrix& M,
278 GlobalMatrix& Jac)
279{
280 assert(x.size() == x_prev.size());
281 for (std::size_t i = 0; i < x.size(); i++)
282 {
285 }
286
287 assembleWithJacobianConcreteProcess(t, dt, x, x_prev, process_id, M, K, b,
288 Jac);
289
290 // TODO: apply BCs to Jacobian.
291 _boundary_conditions[process_id].applyNaturalBC(t, x, process_id, K, b,
292 &Jac);
293
294 // the last argument is for the jacobian, nullptr is for a unused jacobian
295 _source_term_collections[process_id].integrate(t, *x[process_id], b, &Jac);
296}
297
299{
301 {
303
304 return;
305 }
306
307 // For staggered scheme:
308 const int specified_process_id = 0;
310}
311
313{
314 // Create single component dof in every of the mesh nodes.
316 std::make_unique<MeshLib::MeshSubset>(_mesh, _mesh.getNodes());
317
318 // Vector of mesh subsets.
319 std::vector<MeshLib::MeshSubset> all_mesh_subsets;
320
321 // Collect the mesh subsets in a vector for the components of each
322 // variables.
323 for (ProcessVariable const& pv : _process_variables[0])
324 {
325 std::generate_n(std::back_inserter(all_mesh_subsets),
326 pv.getNumberOfGlobalComponents(),
327 [&]() { return *_mesh_subset_all_nodes; });
328 }
329
330 // Create a vector of the number of variable components
331 std::vector<int> vec_var_n_components;
332 transform(cbegin(_process_variables[0]), cend(_process_variables[0]),
333 back_inserter(vec_var_n_components),
334 [](ProcessVariable const& pv)
335 { return pv.getNumberOfGlobalComponents(); });
336
338 std::make_unique<NumLib::LocalToGlobalIndexMap>(
339 std::move(all_mesh_subsets), vec_var_n_components,
341
343}
344
346 const int specified_process_id)
347{
348 // Create single component dof in every of the mesh nodes.
350 std::make_unique<MeshLib::MeshSubset>(_mesh, _mesh.getNodes());
351
352 // Vector of mesh subsets.
353 std::vector<MeshLib::MeshSubset> all_mesh_subsets;
354
355 // Vector of the number of variable components
356 std::vector<int> vec_var_n_components;
357 // Collect the mesh subsets in a vector for each variables' components.
358 std::generate_n(std::back_inserter(all_mesh_subsets),
359 _process_variables[specified_process_id][0]
360 .get()
361 .getNumberOfGlobalComponents(),
362 [&]() { return *_mesh_subset_all_nodes; });
363
364 // Create a vector of the number of variable components.
365 vec_var_n_components.push_back(_process_variables[specified_process_id][0]
366 .get()
367 .getNumberOfGlobalComponents());
369 std::make_unique<NumLib::LocalToGlobalIndexMap>(
370 std::move(all_mesh_subsets), vec_var_n_components,
372
374}
375
376std::tuple<NumLib::LocalToGlobalIndexMap*, bool>
378{
379 if (_local_to_global_index_map->getNumberOfGlobalComponents() == 1)
380 {
381 // For single-variable-single-component processes reuse the existing DOF
382 // table.
383 const bool manage_storage = false;
384 return std::make_tuple(_local_to_global_index_map.get(),
385 manage_storage);
386 }
387
388 // Otherwise construct a new DOF table.
389 std::vector<MeshLib::MeshSubset> all_mesh_subsets_single_component;
390 all_mesh_subsets_single_component.emplace_back(*_mesh_subset_all_nodes);
391
392 const bool manage_storage = true;
393
394 return std::make_tuple(new NumLib::LocalToGlobalIndexMap(
395 std::move(all_mesh_subsets_single_component),
396 // by location order is needed for output
398 manage_storage);
399}
400
402{
403 NumLib::LocalToGlobalIndexMap* dof_table_single_component;
404 bool manage_storage;
405
406 std::tie(dof_table_single_component, manage_storage) =
408
409 std::unique_ptr<NumLib::Extrapolator> extrapolator(
411 *dof_table_single_component));
412
413 // TODO: Later on the DOF table can change during the simulation!
415 std::move(extrapolator), dof_table_single_component, manage_storage);
416}
417
423
424void Process::preTimestep(std::vector<GlobalVector*> const& x, const double t,
425 const double delta_t, const int process_id)
426{
427 for (auto* const solution : x)
429 preTimestepConcreteProcess(x, t, delta_t, process_id);
430
431 _boundary_conditions[process_id].preTimestep(t, x, process_id);
432}
433
434void Process::postTimestep(std::vector<GlobalVector*> const& x,
435 std::vector<GlobalVector*> const& x_prev,
436 const double t, const double delta_t,
437 int const process_id)
438{
439 for (auto* const solution : x)
440 {
442 }
443 for (auto* const solution : x_prev)
444 {
446 }
447
448 postTimestepConcreteProcess(x, x_prev, t, delta_t, process_id);
449
450 _boundary_conditions[process_id].postTimestep(t, x, process_id);
451}
452
453void Process::postNonLinearSolver(std::vector<GlobalVector*> const& x,
454 std::vector<GlobalVector*> const& x_prev,
455 const double t, double const dt,
456 int const process_id)
457{
458 for (auto* const solution : x)
459 {
461 }
462 for (auto* const solution : x_prev)
463 {
465 }
466
467 postNonLinearSolverConcreteProcess(x, x_prev, t, dt, process_id);
468}
469
471 double const dt,
472 std::vector<GlobalVector*> const& x,
473 GlobalVector const& x_prev,
474 int const process_id)
475{
476 for (auto const* solution : x)
479
480 computeSecondaryVariableConcrete(t, dt, x, x_prev, process_id);
481}
482
483void Process::preIteration(const unsigned iter, const GlobalVector& x)
484{
487}
488
494
495void Process::preOutput(const double t, double const dt,
496 std::vector<GlobalVector*> const& x,
497 std::vector<GlobalVector*> const& x_prev,
498 int const process_id)
499{
500 for (auto const* solution : x)
502
503 preOutputConcreteProcess(t, dt, x, x_prev, process_id);
504}
505
506std::vector<GlobalIndexType>
508{
509 std::vector<GlobalIndexType> indices;
510
511 for (std::size_t process_id = 0; process_id < _process_variables.size();
512 process_id++)
513 {
514 auto const& dof_table_of_process = getDOFTable(process_id);
515
516 auto const& per_process_variables = _process_variables[process_id];
517 for (std::size_t variable_id = 0;
518 variable_id < per_process_variables.size();
519 variable_id++)
520 {
521 auto const& pv = per_process_variables[variable_id];
522 DBUG("Set the initial condition of variable {:s} of process {:d}.",
523 pv.get().getName().data(), process_id);
524
525 if ((pv.get().compensateNonEquilibriumInitialResiduum()))
526 {
527 continue;
528 }
529
530 auto const num_comp = pv.get().getNumberOfGlobalComponents();
531
532 for (int component_id = 0; component_id < num_comp; ++component_id)
533 {
534 auto const& mesh_subset = dof_table_of_process.getMeshSubset(
535 variable_id, component_id);
536 auto const mesh_id = mesh_subset.getMeshID();
537 for (auto const* node : mesh_subset.getNodes())
538 {
539 MeshLib::Location const l(
540 mesh_id, MeshLib::MeshItemType::Node, node->getID());
541
542 auto global_index =
543 std::abs(dof_table_of_process.getGlobalIndex(
544 l, variable_id, component_id));
545
546 indices.push_back(global_index);
547 }
548 }
549 }
550 }
551
552 return indices;
553}
554
555} // namespace ProcessLib
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:30
Global vector based on Eigen vector.
Definition EigenVector.h:25
std::vector< Node * > const & getNodes() const
Get the nodes-vector for the mesh.
Definition Mesh.h:106
void setNodeID(std::size_t node_id)
void setCoordinates(MathLib::Point3d const &coordinates)
int getNumberOfGlobalComponents() const
Returns the number of components of the process variable.
virtual NumLib::IterationResult postIterationConcreteProcess(GlobalVector const &)
Definition Process.h:291
virtual void assembleWithJacobianConcreteProcess(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, GlobalMatrix &Jac)=0
virtual void preTimestepConcreteProcess(std::vector< GlobalVector * > const &, const double, const double, const int)
Definition Process.h:253
std::vector< BoundaryConditionCollection > _boundary_conditions
Definition Process.h:392
virtual void preIterationConcreteProcess(const unsigned, GlobalVector const &)
Definition Process.h:277
void constructDofTableOfSpecifiedProcessStaggeredScheme(const int specified_process_id)
Definition Process.cpp:345
static PROCESSLIB_EXPORT const std::string constant_one_parameter_name
Definition Process.h:46
void preAssemble(const double t, double const dt, GlobalVector const &x) final
Definition Process.cpp:243
virtual void initializeBoundaryConditions(std::map< int, std::shared_ptr< MaterialPropertyLib::Medium > > const &media)
Definition Process.cpp:88
void initialize(std::map< int, std::shared_ptr< MaterialPropertyLib::Medium > > const &media)
Definition Process.cpp:102
void initializeProcessBoundaryConditionsAndSourceTerms(const NumLib::LocalToGlobalIndexMap &dof_table, const int process_id, std::map< int, std::shared_ptr< MaterialPropertyLib::Medium > > const &media)
Definition Process.cpp:73
virtual void initializeConcreteProcess(NumLib::LocalToGlobalIndexMap const &dof_table, MeshLib::Mesh const &mesh, unsigned const integration_order)=0
Process specific initialization called by initialize().
std::vector< std::size_t > _ids_of_active_elements
Union of active element ids per process variable.
Definition Process.h:403
void assembleWithJacobian(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, GlobalMatrix &Jac) final
Definition Process.cpp:273
void assemble(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) final
Definition Process.cpp:249
virtual void postTimestepConcreteProcess(std::vector< GlobalVector * > const &, std::vector< GlobalVector * > const &, const double, const double, int const)
Definition Process.h:261
Process(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, const bool use_monolithic_scheme=true)
Definition Process.cpp:28
virtual void computeSecondaryVariableConcrete(double const, double const, std::vector< GlobalVector * > const &, GlobalVector const &, int const)
Definition Process.h:282
std::unique_ptr< MeshLib::MeshSubset const > _mesh_subset_all_nodes
Definition Process.h:355
void preOutput(const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id)
Definition Process.cpp:495
MeshLib::Mesh & _mesh
Definition Process.h:354
virtual 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)=0
std::vector< std::vector< std::reference_wrapper< ProcessVariable > > > _process_variables
Definition Process.h:387
std::vector< std::size_t > const & getActiveElementIDs() const
Definition Process.h:160
virtual NumLib::LocalToGlobalIndexMap const & getDOFTable(const int) const
Definition Process.h:147
virtual void setInitialConditionsConcreteProcess(std::vector< GlobalVector * > &, double const, int const)
Definition Process.h:229
NumLib::IterationResult postIteration(GlobalVector const &x) final
Definition Process.cpp:489
void postTimestep(std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, const double t, const double delta_t, int const process_id)
Postprocessing after a complete timestep.
Definition Process.cpp:434
void preTimestep(std::vector< GlobalVector * > const &x, const double t, const double delta_t, const int process_id)
Preprocessing before starting assembly for new timestep.
Definition Process.cpp:424
void computeSecondaryVariable(double const t, double const dt, std::vector< GlobalVector * > const &x, GlobalVector const &x_prev, int const process_id)
compute secondary variables for the coupled equations or for output.
Definition Process.cpp:470
void initializeExtrapolator()
Definition Process.cpp:401
virtual void postNonLinearSolverConcreteProcess(std::vector< GlobalVector * > const &, std::vector< GlobalVector * > const &, const double, double const, int const)
Definition Process.h:270
std::vector< std::reference_wrapper< ProcessVariable > > const & getProcessVariables(const int process_id) const
Definition Process.h:155
virtual void constructDofTable()
Definition Process.cpp:298
void constructMonolithicProcessDofTable()
Definition Process.cpp:312
virtual void preAssembleConcreteProcess(const double, double const, GlobalVector const &)
Definition Process.h:236
void postNonLinearSolver(std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, const double t, double const dt, int const process_id)
Definition Process.cpp:453
std::vector< GlobalIndexType > getIndicesOfResiduumWithoutInitialCompensation() const override
Definition Process.cpp:507
unsigned const _integration_order
Definition Process.h:371
std::unique_ptr< NumLib::LocalToGlobalIndexMap > _local_to_global_index_map
Definition Process.h:357
std::vector< SourceTermCollection > _source_term_collections
Definition Process.h:398
void updateDeactivatedSubdomains(double const time, const int process_id)
Definition Process.cpp:202
void computeSparsityPattern()
Definition Process.cpp:418
MathLib::MatrixSpecifications getMatrixSpecifications(const int process_id) const override
Definition Process.cpp:194
void setInitialConditions(std::vector< GlobalVector * > &process_solutions, std::vector< GlobalVector * > const &process_solutions_prev, double const t, int const process_id)
Definition Process.cpp:123
virtual void preOutputConcreteProcess(const double, double const, std::vector< GlobalVector * > const &, std::vector< GlobalVector * > const &, int const)
Definition Process.h:297
ExtrapolatorData _extrapolator_data
Definition Process.h:400
GlobalSparsityPattern _sparsity_pattern
Definition Process.h:379
void preIteration(const unsigned iter, GlobalVector const &x) final
Definition Process.cpp:483
const bool _use_monolithic_scheme
Definition Process.h:366
virtual std::tuple< NumLib::LocalToGlobalIndexMap *, bool > getDOFTableForExtrapolatorData() const
Definition Process.cpp:377
Handles configuration of several secondary variables from the project file.
IterationResult
Status flags telling the NonlinearSolver if an iteration succeeded.
void finalizeAssembly(PETScMatrix &A)
Definition LinAlg.cpp:170
void copy(PETScVector const &x, PETScVector &y)
Definition LinAlg.cpp:37
void setLocalAccessibleVector(PETScVector const &x)
Definition LinAlg.cpp:27
@ BY_LOCATION
Ordering data by spatial location.
GlobalSparsityPattern computeSparsityPattern(LocalToGlobalIndexMap const &dof_table, MeshLib::Mesh const &mesh)
Computes a sparsity pattern for the given inputs.