OGS
TimeLoop.cpp
Go to the documentation of this file.
1
11#include "TimeLoop.h"
12
13#include "BaseLib/Error.h"
14#include "BaseLib/RunTime.h"
20#include "ProcessData.h"
22namespace
23{
25{
26 return process_data.process.isMonolithicSchemeUsed();
27}
28
30 std::vector<std::unique_ptr<ProcessLib::ProcessData>> const&
31 per_process_data,
32 double const t)
33{
34 for (auto& process_data : per_process_data)
35 {
36 process_data->process.updateDeactivatedSubdomains(
37 t, process_data->process_id);
38 }
39}
40
41} // namespace
42
43namespace ProcessLib
44{
46 double const t, double const dt,
47 std::vector<std::unique_ptr<ProcessData>> const& per_process_data,
48 std::vector<GlobalVector*> const& _process_solutions)
49{
50 for (auto& process_data : per_process_data)
51 {
52 auto const process_id = process_data->process_id;
53 auto& pcs = process_data->process;
54 pcs.preTimestep(_process_solutions, t, dt, process_id);
55 }
56}
57
59 double const t, double const dt,
60 std::vector<std::unique_ptr<ProcessData>> const& per_process_data,
61 std::vector<GlobalVector*> const& process_solutions,
62 std::vector<GlobalVector*> const& process_solutions_prev,
63 std::vector<std::size_t>& xdot_vector_ids)
64{
65 std::vector<GlobalVector*> x_dots;
66 x_dots.reserve(per_process_data.size());
67 xdot_vector_ids.resize(per_process_data.size());
68
69 std::size_t cnt = 0;
70 for (auto& process_data : per_process_data)
71 {
72 auto const process_id = process_data->process_id;
73 auto const& ode_sys = *process_data->tdisc_ode_sys;
74 auto const& time_discretization = *process_data->time_disc;
75
76 x_dots.emplace_back(&NumLib::GlobalVectorProvider::provider.getVector(
77 ode_sys.getMatrixSpecifications(process_id), xdot_vector_ids[cnt]));
78 cnt++;
79
80 time_discretization.getXdot(*process_solutions[process_id],
81 *process_solutions_prev[process_id],
82 *x_dots[process_id]);
83 }
84
85 // All _per_process_data share the first process.
86 bool const is_staggered_coupling =
87 !isMonolithicProcess(*per_process_data[0]);
88
89 for (auto& process_data : per_process_data)
90 {
91 auto const process_id = process_data->process_id;
92 auto& pcs = process_data->process;
93
94 if (is_staggered_coupling)
95 {
96 CoupledSolutionsForStaggeredScheme coupled_solutions(
97 process_solutions);
98 pcs.setCoupledSolutionsForStaggeredScheme(&coupled_solutions);
99 }
100 auto& x_dot = *x_dots[process_id];
101 pcs.computeSecondaryVariable(t, dt, process_solutions, x_dot,
102 process_id);
103 pcs.postTimestep(process_solutions, t, dt, process_id);
104 }
105 for (auto& x_dot : x_dots)
106 {
108 }
109}
110
111template <NumLib::ODESystemTag ODETag>
113 ProcessData& process_data,
115{
116 using Tag = NumLib::NonlinearSolverTag;
117 // A concrete Picard solver
118 using NonlinearSolverPicard = NumLib::NonlinearSolver<Tag::Picard>;
119 // A concrete Newton solver
120 using NonlinearSolverNewton = NumLib::NonlinearSolver<Tag::Newton>;
121
122 if (dynamic_cast<NonlinearSolverPicard*>(&process_data.nonlinear_solver))
123 {
124 // The Picard solver can also work with a Newton-ready ODE,
125 // because the Newton ODESystem derives from the Picard ODESystem.
126 // So no further checks are needed here.
127
128 process_data.tdisc_ode_sys = std::make_unique<
130 process_data.process_id, ode_sys, *process_data.time_disc);
131 }
132 // TODO (naumov) Provide a function to nonlinear_solver to distinguish the
133 // types. Could be handy, because a nonlinear solver could handle both types
134 // like PETScSNES.
135 else if ((dynamic_cast<NonlinearSolverNewton*>(
136 &process_data.nonlinear_solver) != nullptr)
137#ifdef USE_PETSC
138 || (dynamic_cast<NumLib::PETScNonlinearSolver*>(
139 &process_data.nonlinear_solver) != nullptr)
140#endif // USE_PETSC
141 )
142 {
143 // The Newton-Raphson method needs a Newton-ready ODE.
144
146 if (auto* ode_newton = dynamic_cast<ODENewton*>(&ode_sys))
147 {
148 process_data.tdisc_ode_sys = std::make_unique<
150 process_data.process_id, *ode_newton, *process_data.time_disc);
151 }
152 else
153 {
154 OGS_FATAL(
155 "You are trying to solve a non-Newton-ready ODE with the"
156 " Newton-Raphson method. Aborting");
157 }
158 }
159 else
160 {
161 OGS_FATAL("Encountered unknown nonlinear solver type. Aborting");
162 }
163}
164
166{
167 setTimeDiscretizedODESystem(process_data, process_data.process);
168}
169
170std::pair<std::vector<GlobalVector*>, std::vector<GlobalVector*>>
172 double const t0,
173 std::vector<std::unique_ptr<ProcessData>> const& per_process_data)
174{
175 std::vector<GlobalVector*> process_solutions;
176 std::vector<GlobalVector*> process_solutions_prev;
177
178 for (auto& process_data : per_process_data)
179 {
180 auto const process_id = process_data->process_id;
181 auto& ode_sys = *process_data->tdisc_ode_sys;
182
183 // append a solution vector of suitable size
184 process_solutions.emplace_back(
186 ode_sys.getMatrixSpecifications(process_id)));
187 process_solutions_prev.emplace_back(
189 ode_sys.getMatrixSpecifications(process_id)));
190 }
191
192 for (auto& process_data : per_process_data)
193 {
194 auto& pcs = process_data->process;
195 auto const process_id = process_data->process_id;
196 pcs.setInitialConditions(process_solutions, process_solutions_prev, t0,
197 process_id);
198
199 auto& time_disc = *process_data->time_disc;
200 time_disc.setInitialState(t0); // push IC
201 }
202
203 return {process_solutions, process_solutions_prev};
204}
205
207 std::vector<std::unique_ptr<ProcessData>> const& per_process_data,
208 std::vector<GlobalVector*>
209 process_solutions,
210 std::vector<GlobalVector*> const& process_solutions_prev)
211{
212 for (auto& process_data : per_process_data)
213 {
214 auto& time_disc = *process_data->time_disc;
215 auto& nonlinear_solver = process_data->nonlinear_solver;
216
217 setEquationSystem(*process_data);
218 // dummy values to handle the time derivative terms more or less
219 // correctly, i.e. to ignore them.
220 double const t = 0;
221 double const dt = 1;
222 time_disc.nextTimestep(t, dt);
223 nonlinear_solver.calculateNonEquilibriumInitialResiduum(
224 process_solutions, process_solutions_prev,
225 process_data->process_id);
226 }
227}
228
230 std::vector<GlobalVector*>& x, std::vector<GlobalVector*> const& x_prev,
231 std::size_t const timestep, double const t, double const delta_t,
232 ProcessData const& process_data, Output& output_control,
233 std::size_t& xdot_id)
234{
235 auto& process = process_data.process;
236 int const process_id = process_data.process_id;
237 auto& time_disc = *process_data.time_disc;
238 auto& ode_sys = *process_data.tdisc_ode_sys;
239 auto& nonlinear_solver = process_data.nonlinear_solver;
240
241 setEquationSystem(process_data);
242
243 // Note: Order matters!
244 // First advance to the next timestep, then set known solutions at that
245 // time, afterwards pass the right solution vector and time to the
246 // preTimestep() hook.
247
248 time_disc.nextTimestep(t, delta_t);
249
250 auto const post_iteration_callback =
251 [&](int iteration, std::vector<GlobalVector*> const& x)
252 {
253 output_control.doOutputNonlinearIteration(process, process_id, timestep,
254 t, iteration, x);
255 };
256
257 auto const nonlinear_solver_status =
258 nonlinear_solver.solve(x, x_prev, post_iteration_callback, process_id);
259
260 if (!nonlinear_solver_status.error_norms_met)
261 {
262 return nonlinear_solver_status;
263 }
264
266 ode_sys.getMatrixSpecifications(process_id), xdot_id);
267
268 time_disc.getXdot(*x[process_id], *x_prev[process_id], x_dot);
269
270 process.postNonLinearSolver(*x[process_id], x_dot, t, delta_t, process_id);
272
273 return nonlinear_solver_status;
274}
275
276TimeLoop::TimeLoop(std::unique_ptr<Output>&& output,
277 std::vector<std::unique_ptr<ProcessData>>&& per_process_data,
278 const int global_coupling_max_iterations,
279 std::vector<std::unique_ptr<NumLib::ConvergenceCriterion>>&&
280 global_coupling_conv_crit,
281 const double start_time, const double end_time)
282 : _output(std::move(output)),
283 _per_process_data(std::move(per_process_data)),
284 _start_time(start_time),
285 _end_time(end_time),
286 _global_coupling_max_iterations(global_coupling_max_iterations),
287 _global_coupling_conv_crit(std::move(global_coupling_conv_crit))
288{
289}
290
292{
293 for (auto& process_data : _per_process_data)
294 {
295 auto const& x = *_process_solutions[process_data->process_id];
296
297 // Create a vector to store the solution of the last coupling iteration
300
301 // append a solution vector of suitable size
302 _solutions_of_last_cpl_iteration.emplace_back(&x0);
303 }
304}
305
306std::pair<double, bool> TimeLoop::computeTimeStepping(
307 const double prev_dt, double& t, std::size_t& accepted_steps,
308 std::size_t& rejected_steps,
309 std::vector<std::function<double(double, double)>> const&
310 time_step_constraints)
311{
312 bool all_process_steps_accepted = true;
313 // Get minimum time step size among step sizes of all processes.
314 double dt = std::numeric_limits<double>::max();
315 constexpr double eps = std::numeric_limits<double>::epsilon();
316
317 bool const is_initial_step =
319 [](auto const& ppd) -> bool
320 { return ppd->timestep_current.timeStepNumber() == 0; });
321
322 auto computeSolutionError = [this, t](auto const i) -> double
323 {
324 auto const& ppd = *_per_process_data[i];
325 const auto& timestep_algorithm = ppd.timestep_algorithm;
326 if (!timestep_algorithm->isSolutionErrorComputationNeeded())
327 {
328 return 0.0;
329 }
330 if (t == timestep_algorithm->begin())
331 {
332 // Always accepts the zeroth step
333 return 0.0;
334 }
335
336 auto const& x = *_process_solutions[i];
337 auto const& x_prev = *_process_solutions_prev[i];
338
339 auto const& conv_crit = ppd.conv_crit;
340 const MathLib::VecNormType norm_type =
341 (conv_crit) ? conv_crit->getVectorNormType()
343
344 const double solution_error =
346 norm_type);
347 return solution_error;
348 };
349
350 for (std::size_t i = 0; i < _per_process_data.size(); i++)
351 {
352 auto& ppd = *_per_process_data[i];
353 const auto& timestep_algorithm = ppd.timestep_algorithm;
354
355 const double solution_error = computeSolutionError(i);
356
357 ppd.timestep_current.setAccepted(
358 ppd.nonlinear_solver_status.error_norms_met);
359
360 auto [previous_step_accepted, timestepper_dt] =
361 timestep_algorithm->next(
362 solution_error, ppd.nonlinear_solver_status.number_iterations,
363 ppd.timestep_previous, ppd.timestep_current);
364
365 if (!previous_step_accepted &&
366 // In case of FixedTimeStepping, which makes
367 // timestep_algorithm->next(...) return false when the ending time
368 // is reached.
369 t + eps < timestep_algorithm->end())
370 {
371 // Not all processes have accepted steps.
372 all_process_steps_accepted = false;
373 }
374
375 if (!ppd.nonlinear_solver_status.error_norms_met)
376 {
377 WARN(
378 "Time step will be rejected due to nonlinear solver "
379 "divergence.");
380 all_process_steps_accepted = false;
381 }
382
383 if (timestepper_dt > eps ||
384 std::abs(t - timestep_algorithm->end()) < eps)
385 {
386 dt = std::min(timestepper_dt, dt);
387 }
388 }
389
390 if (all_process_steps_accepted)
391 {
393 }
394 else
395 {
397 }
398
399 bool last_step_rejected = false;
400 if (!is_initial_step)
401 {
402 if (all_process_steps_accepted)
403 {
404 accepted_steps++;
405 last_step_rejected = false;
406 }
407 else
408 {
409 if (t < _end_time || std::abs(t - _end_time) < eps)
410 {
411 t -= prev_dt;
412 rejected_steps++;
413 last_step_rejected = true;
414 }
415 }
416 }
417
418 // adjust step size considering external communciation_point_calculators
419 for (auto const& time_step_constain : time_step_constraints)
420 {
421 dt = std::min(dt, time_step_constain(t, dt));
422 }
423
424 // Check whether the time stepping is stabilized
425 if (std::abs(dt - prev_dt) < eps)
426 {
427 if (last_step_rejected)
428 {
429 OGS_FATAL(
430 "The new step size of {:g} is the same as that of the previous "
431 "rejected time step. \nPlease re-run ogs with a proper "
432 "adjustment in the numerical settings, \ne.g those for time "
433 "stepper, local or global non-linear solver.",
434 dt);
435 }
436 else
437 {
438 DBUG("The time stepping is stabilized with the step size of {:g}.",
439 dt);
440 }
441 }
442
443 // Reset the time step with the minimum step size, dt
444 // Update the solution of the previous time step.
445 for (std::size_t i = 0; i < _per_process_data.size(); i++)
446 {
447 auto& ppd = *_per_process_data[i];
448 auto& timestep_algorithm = ppd.timestep_algorithm;
449 if (all_process_steps_accepted)
450 {
451 NumLib::updateTimeSteps(dt, ppd.timestep_previous,
452 ppd.timestep_current);
453 timestep_algorithm->resetCurrentTimeStep(dt, ppd.timestep_previous,
454 ppd.timestep_current);
455 }
456
457 if (t == timestep_algorithm->begin())
458 {
459 continue;
460 }
461
462 auto& x = *_process_solutions[i];
463 auto& x_prev = *_process_solutions_prev[i];
464 if (all_process_steps_accepted)
465 {
466 MathLib::LinAlg::copy(x, x_prev); // pushState
467 }
468 else
469 {
470 if (t < _end_time || std::abs(t - _end_time) < eps)
471 {
472 WARN(
473 "Time step {:d} was rejected {:d} times and it will be "
474 "repeated with a reduced step size.",
475 accepted_steps + 1, _repeating_times_of_rejected_step);
476 MathLib::LinAlg::copy(x_prev, x); // popState
477 }
478 }
479 }
480
481 return {dt, last_step_rejected};
482}
483
486{
487 for (auto& process_data : _per_process_data)
488 {
489 auto& pcs = process_data->process;
490 int const process_id = process_data->process_id;
491 _output->addProcess(pcs);
492
493 setTimeDiscretizedODESystem(*process_data);
494
495 if (auto* conv_crit =
497 process_data->conv_crit.get()))
498 {
499 conv_crit->setDOFTable(pcs.getDOFTable(process_id), pcs.getMesh());
500 }
501 }
502
503 // initial solution storage
506
507 // All _per_process_data share the first process.
508 bool const is_staggered_coupling =
510
511 if (is_staggered_coupling)
512 {
514 }
515
516 // Output initial conditions
517 {
518 const bool output_initial_condition = true;
519 outputSolutions(output_initial_condition, 0, _start_time, *_output,
521 }
522
523 auto const& fixed_times = _output->getFixedOutputTimes();
524 std::vector<std::function<double(double, double)>> time_step_constraints{
525 [&fixed_times](double t, double dt)
526 { return NumLib::possiblyClampDtToNextFixedTime(t, dt, fixed_times); },
527 [this](double t, double dt)
528 {
529 if (t < _end_time && t + dt > _end_time)
530 {
531 return _end_time - t;
532 }
533 return dt;
534 }};
535 std::tie(_dt, _last_step_rejected) =
537 _rejected_steps, time_step_constraints);
538
540
543}
544
546{
547 BaseLib::RunTime time_timestep;
548 time_timestep.start();
549
551
552 const std::size_t timesteps = _accepted_steps + 1;
553 // TODO(wenqing): , input option for time unit.
554 INFO("=== Time stepping at step #{:d} and time {:g} with step size {:g}",
555 timesteps, _current_time, _dt);
556
558
560 INFO("[time] Time step #{:d} took {:g} s.", timesteps,
561 time_timestep.elapsed());
563}
564
566{
567 const double prev_dt = _dt;
568 double const current_time = _current_time;
569
570 const std::size_t timesteps = _accepted_steps + 1;
571
572 auto const& fixed_times = _output->getFixedOutputTimes();
573 std::vector<std::function<double(double, double)>> time_step_constraints{
574 [&fixed_times](double t, double dt)
575 { return NumLib::possiblyClampDtToNextFixedTime(t, dt, fixed_times); },
576 [this](double t, double dt)
577 {
578 if (t < _end_time && t + dt > _end_time)
579 {
580 return _end_time - t;
581 }
582 return dt;
583 }};
584
585 // _last_step_rejected is also checked in computeTimeStepping.
586 std::tie(_dt, _last_step_rejected) =
588 _rejected_steps, time_step_constraints);
589
591 {
592 const bool output_initial_condition = false;
593 outputSolutions(output_initial_condition, timesteps, current_time,
595 }
596
597 if (std::abs(_current_time - _end_time) <
598 std::numeric_limits<double>::epsilon() ||
600 {
601 return false;
602 }
603
604 if (_dt < std::numeric_limits<double>::epsilon())
605 {
606 WARN(
607 "Time step size of {:g} is too small.\n"
608 "Time stepping stops at step {:d} and at time of {:g}.",
609 _dt, timesteps, _current_time);
610 return false;
611 }
612
613 return true;
614}
615
617{
618 INFO(
619 "The whole computation of the time stepping took {:d} steps, in which\n"
620 "\t the accepted steps are {:d}, and the rejected steps are {:d}.\n",
622
623 // output last time step
625 {
626 const bool output_initial_condition = false;
627 outputSolutions(output_initial_condition,
630 }
631}
632
633bool TimeLoop::doNonlinearIteration(double const t, double const dt,
634 std::size_t const timesteps)
635{
637
638 // All _per_process_data share the first process.
639 bool const is_staggered_coupling =
641 NumLib::NonlinearSolverStatus nonlinear_solver_status;
642
643 if (is_staggered_coupling)
644 {
645 nonlinear_solver_status =
647 }
648 else
649 {
650 nonlinear_solver_status =
651 solveUncoupledEquationSystems(t, dt, timesteps);
652 }
653
654 // Run post time step only if the last iteration was successful.
655 // Otherwise it runs the risks to get the same errors as in the last
656 // iteration, an exception thrown in assembly, for example.
657 if (nonlinear_solver_status.error_norms_met)
658 {
662 }
663 return nonlinear_solver_status.error_norms_met;
664}
665
667 const double t, const double dt, const std::size_t timestep_id,
668 ProcessData const& process_data, std::vector<GlobalVector*>& x,
669 std::vector<GlobalVector*> const& x_prev, Output& output,
670 std::size_t& xdot_id)
671{
672 BaseLib::RunTime time_timestep_process;
673 time_timestep_process.start();
674
675 auto const nonlinear_solver_status = solveOneTimeStepOneProcess(
676 x, x_prev, timestep_id, t, dt, process_data, output, xdot_id);
677
678 INFO("[time] Solving process #{:d} took {:g} s in time step #{:d} ",
679 process_data.process_id, time_timestep_process.elapsed(), timestep_id);
680
681 return nonlinear_solver_status;
682}
683
684static constexpr std::string_view timestepper_cannot_reduce_dt =
685 "Time stepper cannot reduce the time step size further.";
686
688 const double t, const double dt, const std::size_t timestep_id)
689{
690 NumLib::NonlinearSolverStatus nonlinear_solver_status;
691
692 _xdot_vector_ids.resize(_per_process_data.size());
693 std::size_t cnt = 0;
694
695 for (auto& process_data : _per_process_data)
696 {
697 auto const process_id = process_data->process_id;
698 nonlinear_solver_status = solveMonolithicProcess(
699 t, dt, timestep_id, *process_data, _process_solutions,
701 cnt++;
702
703 process_data->nonlinear_solver_status = nonlinear_solver_status;
704 if (!nonlinear_solver_status.error_norms_met)
705 {
706 ERR("The nonlinear solver failed in time step #{:d} at t = {:g} s "
707 "for process #{:d}.",
708 timestep_id, t, process_id);
709
710 if (!process_data->timestep_algorithm->canReduceTimestepSize(
711 process_data->timestep_current,
712 process_data->timestep_previous))
713 {
714 // save unsuccessful solution
715 _output->doOutputAlways(
716 process_data->process, process_id, timestep_id, t,
717 process_data->nonlinear_solver_status.number_iterations,
720 }
721
722 return nonlinear_solver_status;
723 }
724 }
725
726 return nonlinear_solver_status;
727}
728
731 const double t, const double dt, const std::size_t timestep_id)
732{
733 // Coupling iteration
735 {
736 // Set the flag of the first iteration be true.
737 for (auto& conv_crit : _global_coupling_conv_crit)
738 {
739 conv_crit->preFirstIteration();
740 }
741 }
742 auto resetCouplingConvergenceCriteria = [&]()
743 {
744 for (auto& conv_crit : _global_coupling_conv_crit)
745 {
746 conv_crit->reset();
747 }
748 };
749
750 NumLib::NonlinearSolverStatus nonlinear_solver_status{false, -1};
751 bool coupling_iteration_converged = true;
752 for (int global_coupling_iteration = 0;
753 global_coupling_iteration < _global_coupling_max_iterations;
754 global_coupling_iteration++, resetCouplingConvergenceCriteria())
755 {
756 // TODO(wenqing): use process name
757 coupling_iteration_converged = true;
758 _xdot_vector_ids.resize(_per_process_data.size());
759 std::size_t cnt = 0;
760 for (auto& process_data : _per_process_data)
761 {
762 auto const process_id = process_data->process_id;
763 BaseLib::RunTime time_timestep_process;
764 time_timestep_process.start();
765
766 // The following setting of coupled_solutions can be removed only if
767 // the CoupledSolutionsForStaggeredScheme and related functions are
768 // removed totally from the computation of the secondary variable
769 // and from post-time functions.
770 CoupledSolutionsForStaggeredScheme coupled_solutions(
772
773 process_data->process.setCoupledSolutionsForStaggeredScheme(
774 &coupled_solutions);
775
776 nonlinear_solver_status = solveOneTimeStepOneProcess(
778 *process_data, *_output, _xdot_vector_ids[cnt]);
779 cnt++;
780 process_data->nonlinear_solver_status = nonlinear_solver_status;
781
782 INFO(
783 "[time] Solving process #{:d} took {:g} s in time step #{:d} "
784 "coupling iteration #{:d}",
785 process_id, time_timestep_process.elapsed(), timestep_id,
786 global_coupling_iteration);
787
788 if (!nonlinear_solver_status.error_norms_met)
789 {
790 WARN(
791 "The nonlinear solver failed in time step #{:d} at t = "
792 "{:g} s for process #{:d}.",
793 timestep_id, t, process_id);
794 _last_step_rejected = true;
795 return nonlinear_solver_status;
796 }
797
798 // Check the convergence of the coupling iteration
799 auto& x = *_process_solutions[process_id];
800 auto& x_old = *_solutions_of_last_cpl_iteration[process_id];
801 if (global_coupling_iteration > 0)
802 {
803 MathLib::LinAlg::axpy(x_old, -1.0, x); // save dx to x_old
804 INFO(
805 "------- Checking convergence criterion for coupled "
806 "solution of process #{:d} -------",
807 process_id);
808 _global_coupling_conv_crit[process_id]->checkDeltaX(x_old, x);
809 coupling_iteration_converged =
810 coupling_iteration_converged &&
811 _global_coupling_conv_crit[process_id]->isSatisfied();
812 }
813 MathLib::LinAlg::copy(x, x_old);
814 } // end of for (auto& process_data : _per_process_data)
815
816 if (coupling_iteration_converged && global_coupling_iteration > 0)
817 {
818 break;
819 }
820
821 if (!nonlinear_solver_status.error_norms_met)
822 {
823 return nonlinear_solver_status;
824 }
825 }
826
827 if (!coupling_iteration_converged)
828 {
829 WARN(
830 "The coupling iterations reaches its maximum number in time step "
831 "#{:d} at t = {:g} s",
832 timestep_id, t);
833 }
834
835 {
836 for (auto& process_data : _per_process_data)
837 {
838 auto& pcs = process_data->process;
839 int const process_id = process_data->process_id;
840 auto& ode_sys = *process_data->tdisc_ode_sys;
841 pcs.solveReactionEquation(_process_solutions,
842 _process_solutions_prev, t, dt, ode_sys,
843 process_id);
844 }
845 }
846
847 return nonlinear_solver_status;
848}
849
850template <typename OutputClass, typename OutputClassMember>
851void TimeLoop::outputSolutions(bool const output_initial_condition,
852 unsigned timestep, const double t,
853 OutputClass& output_object,
854 OutputClassMember output_class_member) const
855{
856 // All _per_process_data share the first process.
857 bool const is_staggered_coupling =
859
860 for (auto& process_data : _per_process_data)
861 {
862 // If nonlinear solver diverged, the solution has already been
863 // saved.
864 if (!process_data->nonlinear_solver_status.error_norms_met)
865 {
866 continue;
867 }
868
869 auto const process_id = process_data->process_id;
870 auto& pcs = process_data->process;
871
872 if (!is_staggered_coupling && output_initial_condition)
873 {
874 auto const& ode_sys = *process_data->tdisc_ode_sys;
875 // dummy value to handle the time derivative terms more or less
876 // correctly, i.e. to ignore them.
877 double const dt = 1;
878 process_data->time_disc->nextTimestep(t, dt);
879
881 ode_sys.getMatrixSpecifications(process_id));
882 x_dot.setZero();
883
884 pcs.preTimestep(_process_solutions, _start_time, dt, process_id);
885 // Update secondary variables, which might be uninitialized, before
886 // output.
887 pcs.computeSecondaryVariable(_start_time, dt, _process_solutions,
888 x_dot, process_id);
889
891 }
892 if (is_staggered_coupling && output_initial_condition)
893 {
894 CoupledSolutionsForStaggeredScheme coupled_solutions(
896
897 process_data->process.setCoupledSolutionsForStaggeredScheme(
898 &coupled_solutions);
899 process_data->process
900 .setCoupledTermForTheStaggeredSchemeToLocalAssemblers(
901 process_id);
902
903 auto const& ode_sys = *process_data->tdisc_ode_sys;
904 // dummy value to handle the time derivative terms more or less
905 // correctly, i.e. to ignore them.
906 double const dt = 1;
907 process_data->time_disc->nextTimestep(t, dt);
908
910 ode_sys.getMatrixSpecifications(process_id));
911 x_dot.setZero();
912
913 pcs.preTimestep(_process_solutions, _start_time, dt, process_id);
914 // Update secondary variables, which might be uninitialized, before
915 // output.
916 pcs.computeSecondaryVariable(_start_time, dt, _process_solutions,
917 x_dot, process_id);
918
920 }
921 (output_object.*output_class_member)(
922 pcs, process_id, timestep, t,
923 process_data->nonlinear_solver_status.number_iterations,
925 }
926}
927
929{
930 for (auto* x : _process_solutions)
931 {
933 }
934 for (auto* x : _process_solutions_prev)
935 {
937 }
938
940 {
942 }
943}
944
945} // namespace ProcessLib
#define OGS_FATAL(...)
Definition: Error.h:26
void INFO(char const *fmt, Args const &... args)
Definition: Logging.h:34
void ERR(char const *fmt, Args const &... args)
Definition: Logging.h:44
void DBUG(char const *fmt, Args const &... args)
Definition: Logging.h:29
void WARN(char const *fmt, Args const &... args)
Definition: Logging.h:39
Definition of the RunTime class.
Count the running time.
Definition: RunTime.h:29
double elapsed() const
Get the elapsed time in seconds.
Definition: RunTime.h:42
void start()
Start the timer.
Definition: RunTime.h:32
Global vector based on Eigen vector.
Definition: EigenVector.h:28
virtual GlobalVector & getVector(std::size_t &id)=0
Get an uninitialized vector with the given id.
virtual void releaseVector(GlobalVector const &x)=0
void doOutputLastTimestep(Process const &process, const int process_id, int const timestep, const double t, int const iteration, std::vector< GlobalVector * > const &xs)
Definition: Output.cpp:497
void doOutputNonlinearIteration(Process const &process, const int process_id, int const timestep, const double t, const int iteration, std::vector< GlobalVector * > const &xs)
Definition: Output.cpp:514
void doOutput(Process const &process, const int process_id, int const timestep, const double t, int const iteration, std::vector< GlobalVector * > const &xs)
Definition: Output.cpp:478
bool isMonolithicSchemeUsed() const
Definition: Process.h:103
std::pair< double, bool > computeTimeStepping(const double prev_dt, double &t, std::size_t &accepted_steps, std::size_t &rejected_steps, std::vector< std::function< double(double, double)> > const &time_step_constraints)
Definition: TimeLoop.cpp:306
void outputSolutions(bool const output_initial_condition, unsigned timestep, const double t, OutputClass &output_object, OutputClassMember output_class_member) const
Definition: TimeLoop.cpp:851
bool doNonlinearIteration(double const t, double const dt, std::size_t const timesteps)
Definition: TimeLoop.cpp:633
TimeLoop(std::unique_ptr< Output > &&output, std::vector< std::unique_ptr< ProcessData > > &&per_process_data, const int global_coupling_max_iterations, std::vector< std::unique_ptr< NumLib::ConvergenceCriterion > > &&global_coupling_conv_crit, const double start_time, const double end_time)
Definition: TimeLoop.cpp:276
void outputLastTimeStep() const
Definition: TimeLoop.cpp:616
std::vector< GlobalVector * > _solutions_of_last_cpl_iteration
Definition: TimeLoop.h:144
const double _start_time
Definition: TimeLoop.h:127
std::vector< std::size_t > _xdot_vector_ids
Definition: TimeLoop.h:149
NumLib::NonlinearSolverStatus solveUncoupledEquationSystems(const double t, const double dt, const std::size_t timestep_id)
Member to solver non coupled systems of equations, which can be a single system of equations,...
Definition: TimeLoop.cpp:687
std::vector< std::unique_ptr< ProcessData > > _per_process_data
Definition: TimeLoop.h:125
const double _end_time
Definition: TimeLoop.h:128
std::size_t _accepted_steps
Definition: TimeLoop.h:130
bool successful_time_step
Definition: TimeLoop.h:55
std::vector< std::unique_ptr< NumLib::ConvergenceCriterion > > _global_coupling_conv_crit
Convergence criteria of processes for the global coupling iterations.
Definition: TimeLoop.h:140
void setCoupledSolutions()
Definition: TimeLoop.cpp:291
const int _global_coupling_max_iterations
Maximum iterations of the global coupling.
Definition: TimeLoop.h:137
std::vector< GlobalVector * > _process_solutions
Definition: TimeLoop.h:122
std::unique_ptr< Output > _output
Definition: TimeLoop.h:124
void initialize()
initialize output, convergence criterion, etc.
Definition: TimeLoop.cpp:485
int _repeating_times_of_rejected_step
Definition: TimeLoop.h:133
std::size_t _rejected_steps
Definition: TimeLoop.h:131
bool calculateNextTimeStep()
Definition: TimeLoop.cpp:565
std::vector< GlobalVector * > _process_solutions_prev
Definition: TimeLoop.h:123
NumLib::NonlinearSolverStatus solveCoupledEquationSystemsByStaggeredScheme(const double t, const double dt, const std::size_t timestep_id)
Member to solver coupled systems of equations by the staggered scheme.
Definition: TimeLoop.cpp:730
NonlinearSolverTag
Tag used to specify which nonlinear solver will be used.
Definition: Types.h:20
constexpr bool any_of(List const &values)
Checks if any of the elements in the given list is true.
Definition: Algorithm.h:325
void copy(PETScVector const &x, PETScVector &y)
Definition: LinAlg.cpp:37
void axpy(PETScVector &y, PetscScalar const a, PETScVector const &x)
Definition: LinAlg.cpp:57
VecNormType
Norm type. Not declared as class type in order to use the members as integers.
Definition: LinAlgEnums.h:22
static const double t
void updateTimeSteps(double const dt, TimeStep &previous_timestep, TimeStep &current_timestep)
Definition: TimeStep.h:112
double computeRelativeChangeFromPreviousTimestep(GlobalVector const &x, GlobalVector const &x_old, MathLib::VecNormType norm_type)
double possiblyClampDtToNextFixedTime(double const t, double const dt, std::vector< double > const &fixed_output_times)
static NumLib::NonlinearSolverStatus solveMonolithicProcess(const double t, const double dt, const std::size_t timestep_id, ProcessData const &process_data, std::vector< GlobalVector * > &x, std::vector< GlobalVector * > const &x_prev, Output &output, std::size_t &xdot_id)
Definition: TimeLoop.cpp:666
static constexpr std::string_view timestepper_cannot_reduce_dt
Definition: TimeLoop.cpp:684
void setTimeDiscretizedODESystem(ProcessData &process_data, NumLib::ODESystem< ODETag, NumLib::NonlinearSolverTag::Picard > &ode_sys)
Definition: TimeLoop.cpp:112
void preTimestepForAllProcesses(double const t, double const dt, std::vector< std::unique_ptr< ProcessData > > const &per_process_data, std::vector< GlobalVector * > const &_process_solutions)
Definition: TimeLoop.cpp:45
NumLib::NonlinearSolverStatus solveOneTimeStepOneProcess(std::vector< GlobalVector * > &x, std::vector< GlobalVector * > const &x_prev, std::size_t const timestep, double const t, double const delta_t, ProcessData const &process_data, Output &output_control, std::size_t &xdot_id)
Definition: TimeLoop.cpp:229
void postTimestepForAllProcesses(double const t, double const dt, std::vector< std::unique_ptr< ProcessData > > const &per_process_data, std::vector< GlobalVector * > const &process_solutions, std::vector< GlobalVector * > const &process_solutions_prev, std::vector< std::size_t > &xdot_vector_ids)
Definition: TimeLoop.cpp:58
void setEquationSystem(ProcessData const &process_data)
Definition: ProcessData.cpp:17
void calculateNonEquilibriumInitialResiduum(std::vector< std::unique_ptr< ProcessData > > const &per_process_data, std::vector< GlobalVector * > process_solutions, std::vector< GlobalVector * > const &process_solutions_prev)
Definition: TimeLoop.cpp:206
std::pair< std::vector< GlobalVector * >, std::vector< GlobalVector * > > setInitialConditions(double const t0, std::vector< std::unique_ptr< ProcessData > > const &per_process_data)
Definition: TimeLoop.cpp:171
bool isMonolithicProcess(ProcessLib::ProcessData const &process_data)
Definition: TimeLoop.cpp:24
void updateDeactivatedSubdomains(std::vector< std::unique_ptr< ProcessLib::ProcessData > > const &per_process_data, double const t)
Definition: TimeLoop.cpp:29
static NUMLIB_EXPORT VectorProvider & provider
Status of the non-linear solver.
std::unique_ptr< NumLib::TimeDiscretization > time_disc
Definition: ProcessData.h:72
NumLib::NonlinearSolverBase & nonlinear_solver
Definition: ProcessData.h:68
std::unique_ptr< NumLib::EquationSystem > tdisc_ode_sys
type-erased time-discretized ODE system
Definition: ProcessData.h:74