OGS
Process.cpp
Go to the documentation of this file.
1
11#include "Process.h"
12
18#include "ProcessVariable.h"
19
20namespace ProcessLib
21{
23 std::string name_,
24 MeshLib::Mesh& mesh,
25 std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
26 std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
27 unsigned const integration_order,
28 std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
29 process_variables,
30 SecondaryVariableCollection&& secondary_variables,
31 const bool use_monolithic_scheme)
32 : name(std::move(name_)),
33 _mesh(mesh),
34 _secondary_variables(std::move(secondary_variables)),
35 _global_assembler(std::move(jacobian_assembler)),
36 _use_monolithic_scheme(use_monolithic_scheme),
37 _coupled_solutions(nullptr),
38 _integration_order(integration_order),
39 _process_variables(std::move(process_variables)),
40 _boundary_conditions(
41 [&](const std::size_t number_of_process_variables)
42 -> std::vector<BoundaryConditionCollection>
43 {
44 std::vector<BoundaryConditionCollection> pcs_BCs;
45 pcs_BCs.reserve(number_of_process_variables);
46 for (std::size_t i = 0; i < number_of_process_variables; i++)
47 {
48 pcs_BCs.emplace_back(BoundaryConditionCollection(parameters));
49 }
50 return pcs_BCs;
51 }(_process_variables.size())),
52 _source_term_collections(
53 [&](const std::size_t number_of_processes)
54 -> std::vector<SourceTermCollection>
55 {
56 std::vector<SourceTermCollection> pcs_sts;
57 pcs_sts.reserve(number_of_processes);
58 for (std::size_t i = 0; i < number_of_processes; i++)
59 {
60 pcs_sts.emplace_back(SourceTermCollection(parameters));
61 }
62 return pcs_sts;
63 }(_process_variables.size()))
64{
65}
66
68 const NumLib::LocalToGlobalIndexMap& dof_table, const int process_id)
69{
70 auto const& per_process_variables = _process_variables[process_id];
71 auto& per_process_BCs = _boundary_conditions[process_id];
72
73 per_process_BCs.addBCsForProcessVariables(per_process_variables, dof_table,
74 _integration_order, *this);
75
76 auto& per_process_sts = _source_term_collections[process_id];
77 per_process_sts.addSourceTermsForProcessVariables(
78 per_process_variables, dof_table, _integration_order);
79}
80
82{
83 // The number of processes is identical to the size of _process_variables,
84 // the vector contains variables for different processes. See the
85 // documentation of _process_variables.
86 const std::size_t number_of_processes = _process_variables.size();
87 for (std::size_t pcs_id = 0; pcs_id < number_of_processes; pcs_id++)
88 {
91 }
92}
93
95{
96 DBUG("Initialize process.");
97
98 DBUG("Construct dof mappings.");
100
101 DBUG("Compute sparsity pattern");
103
104 DBUG("Initialize the extrapolator");
106
109
110 DBUG("Initialize boundary conditions.");
112}
113
115 std::vector<GlobalVector*>& process_solutions,
116 std::vector<GlobalVector*> const& process_solutions_prev,
117 double const t,
118 int const process_id)
119{
120 auto& x = *process_solutions[process_id];
121 auto& x_prev = *process_solutions_prev[process_id];
122
123 // getDOFTableOfProcess can be overloaded by the specific process.
124 auto const& dof_table_of_process = getDOFTable(process_id);
125
126 auto const& per_process_variables = _process_variables[process_id];
127 for (std::size_t variable_id = 0;
128 variable_id < per_process_variables.size();
129 variable_id++)
130 {
133
134 auto const& pv = per_process_variables[variable_id];
135 DBUG("Set the initial condition of variable {:s} of process {:d}.",
136 pv.get().getName().data(), process_id);
137
138 auto const& ic = pv.get().getInitialCondition();
139
140 auto const num_comp = pv.get().getNumberOfGlobalComponents();
141
142 for (int component_id = 0; component_id < num_comp; ++component_id)
143 {
144 auto const& mesh_subset =
145 dof_table_of_process.getMeshSubset(variable_id, component_id);
146 auto const mesh_id = mesh_subset.getMeshID();
147 for (auto const* node : mesh_subset.getNodes())
148 {
150 node->getID());
151
152 pos.setNodeID(node->getID());
153 pos.setCoordinates(*node);
154 auto const& ic_value = ic(t, pos);
155
156 auto global_index =
157 std::abs(dof_table_of_process.getGlobalIndex(l, variable_id,
158 component_id));
159#ifdef USE_PETSC
160 // The global indices of the ghost entries of the global
161 // matrix or the global vectors need to be set as negative
162 // values for equation assembly, however the global indices
163 // start from zero. Therefore, any ghost entry with zero
164 // index is assigned an negative value of the vector size
165 // or the matrix dimension. To assign the initial value for
166 // the ghost entries, the negative indices of the ghost
167 // entries are restored to zero.
168 if (global_index == x.size())
169 global_index = 0;
170#endif
171 x.set(global_index, ic_value[component_id]);
172 }
173 }
174 }
175
177 MathLib::LinAlg::copy(x, x_prev); // pushState
178
181
182 setInitialConditionsConcreteProcess(process_solutions, t, process_id);
183}
184
186 const int /*process_id*/) const
187{
188 auto const& l = *_local_to_global_index_map;
189 return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
190 &l.getGhostIndices(), &_sparsity_pattern};
191}
192
194 const int process_id)
195{
196 auto const& variables_per_process = getProcessVariables(process_id);
197 for (auto const& variable : variables_per_process)
198 {
199 variable.get().updateDeactivatedSubdomains(time);
200 }
201}
202
203void Process::preAssemble(const double t, double const dt,
204 GlobalVector const& x)
205{
207}
208
209void Process::assemble(const double t, double const dt,
210 std::vector<GlobalVector*> const& x,
211 std::vector<GlobalVector*> const& xdot,
212 int const process_id, GlobalMatrix& M, GlobalMatrix& K,
213 GlobalVector& b)
214{
215 assert(x.size() == xdot.size());
216 for (std::size_t i = 0; i < x.size(); i++)
217 {
220 }
221
222 assembleConcreteProcess(t, dt, x, xdot, process_id, M, K, b);
223
224 // the last argument is for the jacobian, nullptr is for a unused jacobian
225 _boundary_conditions[process_id].applyNaturalBC(t, x, process_id, K, b,
226 nullptr);
227
228 // the last argument is for the jacobian, nullptr is for a unused jacobian
229 _source_term_collections[process_id].integrate(t, *x[process_id], b,
230 nullptr);
231}
232
233void Process::assembleWithJacobian(const double t, double const dt,
234 std::vector<GlobalVector*> const& x,
235 std::vector<GlobalVector*> const& xdot,
236 int const process_id, GlobalMatrix& M,
238 GlobalMatrix& Jac)
239{
240 assert(x.size() == xdot.size());
241 for (std::size_t i = 0; i < x.size(); i++)
242 {
245 }
246
247 assembleWithJacobianConcreteProcess(t, dt, x, xdot, process_id, M, K, b,
248 Jac);
249
250 // TODO: apply BCs to Jacobian.
251 _boundary_conditions[process_id].applyNaturalBC(t, x, process_id, K, b,
252 &Jac);
253
254 // the last argument is for the jacobian, nullptr is for a unused jacobian
255 _source_term_collections[process_id].integrate(t, *x[process_id], b, &Jac);
256}
257
259{
261 {
263
264 return;
265 }
266
267 // For staggered scheme:
268 const int specified_process_id = 0;
270}
271
273{
274 // Create single component dof in every of the mesh nodes.
276 std::make_unique<MeshLib::MeshSubset>(_mesh, _mesh.getNodes());
277
278 // Vector of mesh subsets.
279 std::vector<MeshLib::MeshSubset> all_mesh_subsets;
280
281 // Collect the mesh subsets in a vector for the components of each
282 // variables.
283 for (ProcessVariable const& pv : _process_variables[0])
284 {
285 std::generate_n(std::back_inserter(all_mesh_subsets),
286 pv.getNumberOfGlobalComponents(),
287 [&]() { return *_mesh_subset_all_nodes; });
288 }
289
290 // Create a vector of the number of variable components
291 std::vector<int> vec_var_n_components;
292 transform(cbegin(_process_variables[0]), cend(_process_variables[0]),
293 back_inserter(vec_var_n_components),
294 [](ProcessVariable const& pv)
295 { return pv.getNumberOfGlobalComponents(); });
296
298 std::make_unique<NumLib::LocalToGlobalIndexMap>(
299 std::move(all_mesh_subsets), vec_var_n_components,
301
303}
304
306 const int specified_process_id)
307{
308 // Create single component dof in every of the mesh nodes.
310 std::make_unique<MeshLib::MeshSubset>(_mesh, _mesh.getNodes());
311
312 // Vector of mesh subsets.
313 std::vector<MeshLib::MeshSubset> all_mesh_subsets;
314
315 // Vector of the number of variable components
316 std::vector<int> vec_var_n_components;
317 // Collect the mesh subsets in a vector for each variables' components.
318 std::generate_n(std::back_inserter(all_mesh_subsets),
319 _process_variables[specified_process_id][0]
320 .get()
321 .getNumberOfGlobalComponents(),
322 [&]() { return *_mesh_subset_all_nodes; });
323
324 // Create a vector of the number of variable components.
325 vec_var_n_components.push_back(_process_variables[specified_process_id][0]
326 .get()
327 .getNumberOfGlobalComponents());
329 std::make_unique<NumLib::LocalToGlobalIndexMap>(
330 std::move(all_mesh_subsets), vec_var_n_components,
332
334}
335
336std::tuple<NumLib::LocalToGlobalIndexMap*, bool>
338{
339 if (_local_to_global_index_map->getNumberOfGlobalComponents() == 1)
340 {
341 // For single-variable-single-component processes reuse the existing DOF
342 // table.
343 const bool manage_storage = false;
344 return std::make_tuple(_local_to_global_index_map.get(),
345 manage_storage);
346 }
347
348 // Otherwise construct a new DOF table.
349 std::vector<MeshLib::MeshSubset> all_mesh_subsets_single_component;
350 all_mesh_subsets_single_component.emplace_back(*_mesh_subset_all_nodes);
351
352 const bool manage_storage = true;
353
354 return std::make_tuple(new NumLib::LocalToGlobalIndexMap(
355 std::move(all_mesh_subsets_single_component),
356 // by location order is needed for output
358 manage_storage);
359}
360
362{
363 NumLib::LocalToGlobalIndexMap* dof_table_single_component;
364 bool manage_storage;
365
366 std::tie(dof_table_single_component, manage_storage) =
368
369 std::unique_ptr<NumLib::Extrapolator> extrapolator(
371 *dof_table_single_component));
372
373 // TODO: Later on the DOF table can change during the simulation!
375 std::move(extrapolator), dof_table_single_component, manage_storage);
376}
377
379{
382}
383
384void Process::preTimestep(std::vector<GlobalVector*> const& x, const double t,
385 const double delta_t, const int process_id)
386{
387 for (auto* const solution : x)
389 preTimestepConcreteProcess(x, t, delta_t, process_id);
390
391 _boundary_conditions[process_id].preTimestep(t, x, process_id);
392}
393
394void Process::postTimestep(std::vector<GlobalVector*> const& x, const double t,
395 const double delta_t, int const process_id)
396{
397 for (auto* const solution : x)
399 postTimestepConcreteProcess(x, t, delta_t, process_id);
400
401 _boundary_conditions[process_id].postTimestep(t, x, process_id);
402}
403
405 GlobalVector const& xdot, const double t,
406 double const dt, int const process_id)
407{
410 postNonLinearSolverConcreteProcess(x, xdot, t, dt, process_id);
411}
412
414 double const dt,
415 std::vector<GlobalVector*> const& x,
416 GlobalVector const& x_dot,
417 int const process_id)
418{
419 for (auto const* solution : x)
422
423 computeSecondaryVariableConcrete(t, dt, x, x_dot, process_id);
424}
425
426void Process::preIteration(const unsigned iter, const GlobalVector& x)
427{
430}
431
433{
436}
437
438std::vector<GlobalIndexType>
440{
441 std::vector<GlobalIndexType> indices;
442
443 for (std::size_t process_id = 0; process_id < _process_variables.size();
444 process_id++)
445 {
446 auto const& dof_table_of_process = getDOFTable(process_id);
447
448 auto const& per_process_variables = _process_variables[process_id];
449 for (std::size_t variable_id = 0;
450 variable_id < per_process_variables.size();
451 variable_id++)
452 {
453 auto const& pv = per_process_variables[variable_id];
454 DBUG("Set the initial condition of variable {:s} of process {:d}.",
455 pv.get().getName().data(), process_id);
456
457 if ((pv.get().compensateNonEquilibriumInitialResiduum()))
458 {
459 continue;
460 }
461
462 auto const num_comp = pv.get().getNumberOfGlobalComponents();
463
464 for (int component_id = 0; component_id < num_comp; ++component_id)
465 {
466 auto const& mesh_subset = dof_table_of_process.getMeshSubset(
467 variable_id, component_id);
468 auto const mesh_id = mesh_subset.getMeshID();
469 for (auto const* node : mesh_subset.getNodes())
470 {
471 MeshLib::Location const l(
472 mesh_id, MeshLib::MeshItemType::Node, node->getID());
473
474 auto global_index =
475 std::abs(dof_table_of_process.getGlobalIndex(
476 l, variable_id, component_id));
477
478 indices.push_back(global_index);
479 }
480 }
481 }
482 }
483
484 return indices;
485}
486
487} // namespace ProcessLib
void DBUG(char const *fmt, Args const &... args)
Definition: Logging.h:29
Global vector based on Eigen vector.
Definition: EigenVector.h:28
std::vector< Node * > const & getNodes() const
Get the nodes-vector for the mesh.
Definition: Mesh.h:100
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 void initializeBoundaryConditions()
Definition: Process.cpp:81
virtual NumLib::IterationResult postIterationConcreteProcess(GlobalVector const &)
Definition: Process.h:272
virtual void preTimestepConcreteProcess(std::vector< GlobalVector * > const &, const double, const double, const int)
Definition: Process.h:236
std::vector< BoundaryConditionCollection > _boundary_conditions
Definition: Process.h:367
void constructDofTableOfSpecifiedProcessStaggeredScheme(const int specified_prosess_id)
Definition: Process.cpp:305
virtual void postTimestepConcreteProcess(std::vector< GlobalVector * > const &, const double, const double, int const)
Definition: Process.h:244
void postNonLinearSolver(GlobalVector const &x, GlobalVector const &xdot, const double t, double const dt, int const process_id)
Definition: Process.cpp:404
virtual void preIterationConcreteProcess(const unsigned, GlobalVector const &)
Definition: Process.h:258
virtual 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)=0
void assembleWithJacobian(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) final
Definition: Process.cpp:233
void preAssemble(const double t, double const dt, GlobalVector const &x) final
Definition: Process.cpp:203
virtual void initializeConcreteProcess(NumLib::LocalToGlobalIndexMap const &dof_table, MeshLib::Mesh const &mesh, unsigned const integration_order)=0
Process specific initialization called by initialize().
void initializeProcessBoundaryConditionsAndSourceTerms(const NumLib::LocalToGlobalIndexMap &dof_table, const int process_id)
Definition: Process.cpp:67
void computeSecondaryVariable(double const t, double const dt, std::vector< GlobalVector * > const &x, GlobalVector const &x_dot, int const process_id)
compute secondary variables for the coupled equations or for output.
Definition: Process.cpp:413
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:22
virtual void computeSecondaryVariableConcrete(double const, double const, std::vector< GlobalVector * > const &, GlobalVector const &, int const)
Definition: Process.h:263
std::unique_ptr< MeshLib::MeshSubset const > _mesh_subset_all_nodes
Definition: Process.h:329
MeshLib::Mesh & _mesh
Definition: Process.h:328
std::vector< std::vector< std::reference_wrapper< ProcessVariable > > > _process_variables
Definition: Process.h:362
virtual void postNonLinearSolverConcreteProcess(GlobalVector const &, GlobalVector const &, const double, double const, int const)
Definition: Process.h:252
virtual NumLib::LocalToGlobalIndexMap const & getDOFTable(const int) const
Definition: Process.h:138
virtual 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)=0
virtual void setInitialConditionsConcreteProcess(std::vector< GlobalVector * > &, double const, int const)
Definition: Process.h:211
NumLib::IterationResult postIteration(GlobalVector const &x) final
Definition: Process.cpp:432
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:384
void initializeExtrapolator()
Definition: Process.cpp:361
std::vector< std::reference_wrapper< ProcessVariable > > const & getProcessVariables(const int process_id) const
Definition: Process.h:146
virtual void constructDofTable()
Definition: Process.cpp:258
void constructMonolithicProcessDofTable()
Definition: Process.cpp:272
virtual void preAssembleConcreteProcess(const double, double const, GlobalVector const &)
Definition: Process.h:218
std::vector< GlobalIndexType > getIndicesOfResiduumWithoutInitialCompensation() const override
Definition: Process.cpp:439
unsigned const _integration_order
Definition: Process.h:346
std::unique_ptr< NumLib::LocalToGlobalIndexMap > _local_to_global_index_map
Definition: Process.h:331
std::vector< SourceTermCollection > _source_term_collections
Definition: Process.h:373
void updateDeactivatedSubdomains(double const time, const int process_id)
Definition: Process.cpp:193
void computeSparsityPattern()
Definition: Process.cpp:378
MathLib::MatrixSpecifications getMatrixSpecifications(const int process_id) const override
Definition: Process.cpp:185
void setInitialConditions(std::vector< GlobalVector * > &process_solutions, std::vector< GlobalVector * > const &process_solutions_prev, double const t, int const process_id)
Definition: Process.cpp:114
void assemble(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) final
Definition: Process.cpp:209
void postTimestep(std::vector< GlobalVector * > const &x, const double t, const double delta_t, int const process_id)
Postprocessing after a complete timestep.
Definition: Process.cpp:394
ExtrapolatorData _extrapolator_data
Definition: Process.h:375
GlobalSparsityPattern _sparsity_pattern
Definition: Process.h:354
void preIteration(const unsigned iter, GlobalVector const &x) final
Definition: Process.cpp:426
const bool _use_monolithic_scheme
Definition: Process.h:337
virtual std::tuple< NumLib::LocalToGlobalIndexMap *, bool > getDOFTableForExtrapolatorData() const
Definition: Process.cpp:337
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:163
void copy(PETScVector const &x, PETScVector &y)
Definition: LinAlg.cpp:37
void setLocalAccessibleVector(PETScVector const &x)
Definition: LinAlg.cpp:27
static const double t
@ 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.