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