OGS 6.2.2-87-g988ee9c30.dirty.20200123122242
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"
16 #include "MathLib/LinAlg/LinAlg.h"
21 
23 #include "ProcessData.h"
24 
25 namespace
26 {
29 template <NumLib::NonlinearSolverTag NLTag>
31  NumLib::EquationSystem& eq_sys,
33 {
34  using Solver = NumLib::NonlinearSolver<NLTag>;
35  using EqSys = NumLib::NonlinearSystem<NLTag>;
36 
37  assert(dynamic_cast<Solver*>(&nonlinear_solver) != nullptr);
38  assert(dynamic_cast<EqSys*>(&eq_sys) != nullptr);
39 
40  auto& nl_solver_ = static_cast<Solver&>(nonlinear_solver);
41  auto& eq_sys_ = static_cast<EqSys&>(eq_sys);
42 
43  nl_solver_.setEquationSystem(eq_sys_, conv_crit);
44 }
45 
49  NumLib::EquationSystem& eq_sys,
52 {
53  using Tag = NumLib::NonlinearSolverTag;
54  switch (nl_tag)
55  {
56  case Tag::Picard:
57  setEquationSystem<Tag::Picard>(nonlinear_solver, eq_sys, conv_crit);
58  break;
59  case Tag::Newton:
60  setEquationSystem<Tag::Newton>(nonlinear_solver, eq_sys, conv_crit);
61  break;
62  }
63 }
64 
66 {
67  return process_data.process.isMonolithicSchemeUsed();
68 }
69 } // namespace
70 
71 namespace ProcessLib
72 {
73 template <NumLib::ODESystemTag ODETag>
75  ProcessData& process_data,
77 {
78  using Tag = NumLib::NonlinearSolverTag;
79  // A concrete Picard solver
80  using NonlinearSolverPicard = NumLib::NonlinearSolver<Tag::Picard>;
81  // A concrete Newton solver
82  using NonlinearSolverNewton = NumLib::NonlinearSolver<Tag::Newton>;
83 
84  if (dynamic_cast<NonlinearSolverPicard*>(&process_data.nonlinear_solver))
85  {
86  // The Picard solver can also work with a Newton-ready ODE,
87  // because the Newton ODESystem derives from the Picard ODESystem.
88  // So no further checks are needed here.
89 
90  process_data.tdisc_ode_sys = std::make_unique<
92  process_data.process_id, ode_sys, *process_data.time_disc);
93  }
94  else if (dynamic_cast<NonlinearSolverNewton*>(
95  &process_data.nonlinear_solver))
96  {
97  // The Newton-Raphson method needs a Newton-ready ODE.
98 
100  if (auto* ode_newton = dynamic_cast<ODENewton*>(&ode_sys))
101  {
102  process_data.tdisc_ode_sys = std::make_unique<
104  process_data.process_id, *ode_newton, *process_data.time_disc);
105  }
106  else
107  {
108  OGS_FATAL(
109  "You are trying to solve a non-Newton-ready ODE with the"
110  " Newton-Raphson method. Aborting");
111  }
112  }
113  else
114  {
115  OGS_FATAL("Encountered unknown nonlinear solver type. Aborting");
116  }
117 
118  process_data.mat_strg = dynamic_cast<NumLib::InternalMatrixStorage*>(
119  process_data.tdisc_ode_sys.get());
120 }
121 
123 {
124  setTimeDiscretizedODESystem(process_data, process_data.process);
125 }
126 
127 std::vector<GlobalVector*> setInitialConditions(
128  double const t0,
129  std::vector<std::unique_ptr<ProcessData>> const& per_process_data)
130 {
131  std::vector<GlobalVector*> process_solutions;
132 
133  for (auto& process_data : per_process_data)
134  {
135  auto& pcs = process_data->process;
136  auto const process_id = process_data->process_id;
137  auto& time_disc = *process_data->time_disc;
138 
139  auto& ode_sys = *process_data->tdisc_ode_sys;
140 
141  // append a solution vector of suitable size
142  process_solutions.emplace_back(
144  ode_sys.getMatrixSpecifications(process_id)));
145 
146  auto& x0 = *process_solutions[process_id];
147  pcs.setInitialConditions(process_id, t0, x0);
149 
150  time_disc.setInitialState(t0, x0); // push IC
151 
152  if (time_disc.needsPreload())
153  {
154  auto& nonlinear_solver = process_data->nonlinear_solver;
155  auto& mat_strg = *process_data->mat_strg;
156  auto& conv_crit = *process_data->conv_crit;
157 
158  auto const nl_tag = process_data->nonlinear_solver_tag;
159  setEquationSystem(nonlinear_solver, ode_sys, conv_crit, nl_tag);
160  nonlinear_solver.assemble(process_solutions, process_id);
161  time_disc.pushState(
162  t0, x0,
163  mat_strg); // TODO: that might do duplicate work
164  }
165  }
166 
167  return process_solutions;
168 }
169 
171  std::vector<std::unique_ptr<ProcessData>> const& per_process_data,
172  std::vector<GlobalVector*>
173  process_solutions)
174 {
175  INFO("Calculate non-equilibrium initial residuum.");
176  for (auto& process_data : per_process_data)
177  {
178  auto& ode_sys = *process_data->tdisc_ode_sys;
179 
180  auto& time_disc = *process_data->time_disc;
181  auto& nonlinear_solver = process_data->nonlinear_solver;
182  auto& conv_crit = *process_data->conv_crit;
183 
184  auto const nl_tag = process_data->nonlinear_solver_tag;
185  setEquationSystem(nonlinear_solver, ode_sys, conv_crit, nl_tag);
186  // dummy values to handle the time derivative terms more or less
187  // correctly, i.e. to ignore them.
188  double const t = 0;
189  double const dt = 1;
190  time_disc.nextTimestep(t, dt);
191  nonlinear_solver.calculateNonEquilibriumInitialResiduum(
192  process_solutions, process_data->process_id);
193  }
194 }
195 
197  std::vector<GlobalVector*>& x, std::size_t const timestep, double const t,
198  double const delta_t, ProcessData const& process_data,
199  Output& output_control)
200 {
201  auto& process = process_data.process;
202  int const process_id = process_data.process_id;
203  auto& time_disc = *process_data.time_disc;
204  auto& conv_crit = *process_data.conv_crit;
205  auto& ode_sys = *process_data.tdisc_ode_sys;
206  auto& nonlinear_solver = process_data.nonlinear_solver;
207  auto const nl_tag = process_data.nonlinear_solver_tag;
208 
209  setEquationSystem(nonlinear_solver, ode_sys, conv_crit, nl_tag);
210 
211  // Note: Order matters!
212  // First advance to the next timestep, then set known solutions at that
213  // time, afterwards pass the right solution vector and time to the
214  // preTimestep() hook.
215 
216  time_disc.nextTimestep(t, delta_t);
217 
218  auto const post_iteration_callback =
219  [&](int iteration, std::vector<GlobalVector*> const& x) {
220  output_control.doOutputNonlinearIteration(
221  process, process_id, timestep, t, x, iteration);
222  };
223 
224  auto const nonlinear_solver_status =
225  nonlinear_solver.solve(x, post_iteration_callback, process_id);
226 
227  if (nonlinear_solver_status.error_norms_met)
228  {
229  process.postNonLinearSolver(*x[process_id], t, delta_t, process_id);
230  }
231 
232  return nonlinear_solver_status;
233 }
234 
236  std::unique_ptr<Output>&& output,
237  std::vector<std::unique_ptr<ProcessData>>&& per_process_data,
238  const int global_coupling_max_iterations,
239  std::vector<std::unique_ptr<NumLib::ConvergenceCriterion>>&&
240  global_coupling_conv_crit,
241  std::unique_ptr<ChemistryLib::ChemicalSolverInterface>&& chemical_system,
242  const double start_time, const double end_time)
243  : _output(std::move(output)),
244  _per_process_data(std::move(per_process_data)),
245  _start_time(start_time),
246  _end_time(end_time),
247  _global_coupling_max_iterations(global_coupling_max_iterations),
248  _global_coupling_conv_crit(std::move(global_coupling_conv_crit)),
249  _chemical_system(std::move(chemical_system))
250 {
251 }
252 
254 {
255  for (auto& process_data : _per_process_data)
256  {
257  auto const& x = *_process_solutions[process_data->process_id];
258 
259  // Create a vector to store the solution of the last coupling iteration
261  MathLib::LinAlg::copy(x, x0);
262 
263  // append a solution vector of suitable size
264  _solutions_of_last_cpl_iteration.emplace_back(&x0);
265  }
266 }
267 
268 double TimeLoop::computeTimeStepping(const double prev_dt, double& t,
269  std::size_t& accepted_steps,
270  std::size_t& rejected_steps)
271 {
272  bool all_process_steps_accepted = true;
273  // Get minimum time step size among step sizes of all processes.
274  double dt = std::numeric_limits<double>::max();
275  for (std::size_t i = 0; i < _per_process_data.size(); i++)
276  {
277  auto& ppd = *_per_process_data[i];
278  const auto& timestepper = ppd.timestepper;
279 
280  auto& time_disc = ppd.time_disc;
281  auto const& x = *_process_solutions[i];
282 
283  auto const& conv_crit = ppd.conv_crit;
284  const MathLib::VecNormType norm_type =
285  (conv_crit) ? conv_crit->getVectorNormType()
287 
288  const double solution_error =
289  (timestepper->isSolutionErrorComputationNeeded())
290  ? ((t == timestepper->begin())
291  ? 0. // Always accepts the zeroth step
292  : time_disc->getRelativeChangeFromPreviousTimestep(
293  x, norm_type))
294  : 0.;
295 
296  if (!ppd.nonlinear_solver_status.error_norms_met)
297  {
298  timestepper->setAcceptedOrNot(false);
299  }
300  else
301  {
302  timestepper->setAcceptedOrNot(true);
303  }
304 
305  if (!timestepper->next(solution_error,
306  ppd.nonlinear_solver_status.number_iterations) &&
307  // In case of FixedTimeStepping, which makes timestepper->next(...)
308  // return false when the ending time is reached.
309  t + std::numeric_limits<double>::epsilon() < timestepper->end())
310  {
311  // Not all processes have accepted steps.
312  all_process_steps_accepted = false;
313  }
314 
315  if (!ppd.nonlinear_solver_status.error_norms_met)
316  {
317  WARN("Time step will be rejected due to nonlinear solver diverged");
318  all_process_steps_accepted = false;
319  }
320 
321  if (timestepper->getTimeStep().dt() >
322  std::numeric_limits<double>::min() ||
323  std::abs(t - timestepper->end()) <
324  std::numeric_limits<double>::epsilon())
325  {
326  if (timestepper->getTimeStep().dt() < dt)
327  {
328  dt = timestepper->getTimeStep().dt();
329  }
330  }
331  }
332 
333  if (all_process_steps_accepted)
334  {
336  }
337  else
338  {
340  }
341 
342  bool is_initial_step = false;
343  // Reset the time step with the minimum step size, dt
344  // Update the solution of the previous time step in time_disc.
345  for (std::size_t i = 0; i < _per_process_data.size(); i++)
346  {
347  const auto& ppd = *_per_process_data[i];
348  auto& timestepper = ppd.timestepper;
349  timestepper->resetCurrentTimeStep(dt);
350 
351  if (t == timestepper->begin())
352  {
353  is_initial_step = true;
354  continue;
355  }
356 
357  auto& time_disc = ppd.time_disc;
358  auto& x = *_process_solutions[i];
359  if (all_process_steps_accepted)
360  {
361  time_disc->pushState(t, x, *ppd.mat_strg);
362  }
363  else
364  {
365  if (t < _end_time || std::abs(t - _end_time) <
366  std::numeric_limits<double>::epsilon())
367  {
368  WARN(
369  "Time step %d was rejected %d times "
370  "and it will be repeated with a reduced step size.",
371  accepted_steps + 1, _repeating_times_of_rejected_step);
372  time_disc->popState(x);
373  }
374  }
375  }
376 
377  if (!is_initial_step)
378  {
379  if (all_process_steps_accepted)
380  {
381  accepted_steps++;
382  _last_step_rejected = false;
383  }
384  else
385  {
386  if (t < _end_time || std::abs(t - _end_time) <
387  std::numeric_limits<double>::epsilon())
388  {
389  t -= prev_dt;
390  rejected_steps++;
391  _last_step_rejected = true;
392  }
393  }
394  }
395 
396  // Adjust step size if t < _end_time, while t+dt exceeds the end time
397  if (t < _end_time && t + dt > _end_time)
398  {
399  dt = _end_time - t;
400  }
401 
402  return dt;
403 }
404 
407 {
408  for (auto& process_data : _per_process_data)
409  {
410  auto& pcs = process_data->process;
411  int const process_id = process_data->process_id;
412  _output->addProcess(pcs, process_id);
413 
414  setTimeDiscretizedODESystem(*process_data);
415 
416  if (auto* conv_crit =
417  dynamic_cast<NumLib::ConvergenceCriterionPerComponent*>(
418  process_data->conv_crit.get()))
419  {
420  conv_crit->setDOFTable(pcs.getDOFTable(process_id), pcs.getMesh());
421  }
422 
423  // Add the fixed times of output to time stepper in order that
424  // the time stepping is performed and the results are output at
425  // these times. Note: only the adaptive time steppers can have the
426  // the fixed times.
427  auto& timestepper = process_data->timestepper;
428  timestepper->addFixedOutputTimes(_output->getFixedOutputTimes());
429  }
430 
431  // init solution storage
433 
434  if (_chemical_system != nullptr)
435  {
436  BaseLib::RunTime time_phreeqc;
437  time_phreeqc.start();
438  _chemical_system->executeInitialCalculation(_process_solutions);
439  INFO("[time] Phreeqc took %g s.", time_phreeqc.elapsed());
440  }
441 
442  // All _per_process_data share the first process.
443  bool const is_staggered_coupling =
444  !isMonolithicProcess(*_per_process_data[0]);
445 
446  if (is_staggered_coupling)
447  {
449  }
450 
451  // Output initial conditions
452  {
453  const bool output_initial_condition = true;
454  outputSolutions(output_initial_condition, 0, _start_time, *_output,
456  }
457 
458 }
459 
460 /*
461  * TODO:
462  * Now we have a structure inside the time loop which is very similar to the
463  * nonlinear solver. And admittedly, the control flow inside the nonlinear
464  * solver is rather complicated. Maybe in the future one can introduce an
465  * abstraction that can do both the convergence checks of the coupling loop and
466  * of the nonlinear solver.
467  */
469 {
470  // All _per_process_data share the first process.
471  bool const is_staggered_coupling =
473 
474  bool non_equilibrium_initial_residuum_computed = false;
475  double t = _start_time;
476  std::size_t accepted_steps = 0;
477  std::size_t rejected_steps = 0;
478  NumLib::NonlinearSolverStatus nonlinear_solver_status;
479 
480  double dt = computeTimeStepping(0.0, t, accepted_steps, rejected_steps);
481 
482  while (t < _end_time)
483  {
484  BaseLib::RunTime time_timestep;
485  time_timestep.start();
486 
487  t += dt;
488  const double prev_dt = dt;
489 
490  const std::size_t timesteps = accepted_steps + 1;
491  // TODO(wenqing): , input option for time unit.
492  INFO("=== Time stepping at step #%u and time %g with step size %g",
493  timesteps, t, dt);
494 
495  // Check element deactivation:
496  for (auto& process_data : _per_process_data)
497  {
498  process_data->process.updateDeactivatedSubdomains(
499  t, process_data->process_id);
500  }
501 
502  if (!non_equilibrium_initial_residuum_computed)
503  {
504  calculateNonEquilibriumInitialResiduum(_per_process_data,
506  non_equilibrium_initial_residuum_computed = true;
507  }
508 
509  if (is_staggered_coupling)
510  {
511  nonlinear_solver_status =
513  }
514  else
515  {
516  nonlinear_solver_status =
517  solveUncoupledEquationSystems(t, dt, timesteps);
518  }
519 
520  INFO("[time] Time step #%u took %g s.", timesteps,
521  time_timestep.elapsed());
522 
523  dt = computeTimeStepping(prev_dt, t, accepted_steps, rejected_steps);
524 
525  if (!_last_step_rejected)
526  {
527  const bool output_initial_condition = false;
528  outputSolutions(output_initial_condition, timesteps, t, *_output,
530  }
531 
532  if (t == _end_time || t + dt > _end_time ||
533  t + std::numeric_limits<double>::epsilon() > _end_time)
534  {
535  break;
536  }
537 
538  if (dt < std::numeric_limits<double>::epsilon())
539  {
540  WARN(
541  "Time step size of %g is too small.\n"
542  "Time stepping stops at step %u and at time of %g.",
543  dt, timesteps, t);
544  break;
545  }
546  }
547 
548  INFO(
549  "The whole computation of the time stepping took %u steps, in which\n"
550  "\t the accepted steps are %u, and the rejected steps are %u.\n",
551  accepted_steps + rejected_steps, accepted_steps, rejected_steps);
552 
553  // output last time step
554  if (nonlinear_solver_status.error_norms_met)
555  {
556  const bool output_initial_condition = false;
557  outputSolutions(output_initial_condition,
558  accepted_steps + rejected_steps, t, *_output,
560  }
561 
562  return nonlinear_solver_status.error_norms_met;
563 }
564 
566  double const t, double const dt,
567  std::vector<std::unique_ptr<ProcessData>> const& per_process_data,
568  std::vector<GlobalVector*> const& _process_solutions)
569 {
570  for (auto& process_data : per_process_data)
571  {
572  auto const process_id = process_data->process_id;
573  auto& pcs = process_data->process;
574  pcs.preTimestep(_process_solutions, t, dt, process_id);
575  }
576 }
577 
579  const double t, const double dt, const std::size_t timestep_id,
580  ProcessData const& process_data, std::vector<GlobalVector*>& x,
581  Output& output)
582 {
583  BaseLib::RunTime time_timestep_process;
584  time_timestep_process.start();
585 
586  auto const nonlinear_solver_status =
587  solveOneTimeStepOneProcess(x, timestep_id, t, dt, process_data, output);
588 
589  INFO("[time] Solving process #%u took %g s in time step #%u ",
590  process_data.process_id, time_timestep_process.elapsed(), timestep_id);
591 
592  return nonlinear_solver_status;
593 }
594 
595 static constexpr std::string_view timestepper_cannot_reduce_dt =
596  "Time stepper cannot reduce the time step size further.";
597 
599  double const t, double const dt,
600  std::vector<std::unique_ptr<ProcessData>> const& per_process_data,
601  std::vector<GlobalVector*> const& process_solutions)
602 {
603  // All _per_process_data share the first process.
604  bool const is_staggered_coupling =
605  !isMonolithicProcess(*per_process_data[0]);
606 
607  for (auto& process_data : per_process_data)
608  {
609  auto const process_id = process_data->process_id;
610  auto& pcs = process_data->process;
611 
612  if (is_staggered_coupling)
613  {
614  CoupledSolutionsForStaggeredScheme coupled_solutions(
615  process_solutions);
616  pcs.setCoupledSolutionsForStaggeredScheme(&coupled_solutions);
617  }
618  auto& x = *process_solutions[process_id];
619  pcs.postTimestep(process_solutions, t, dt, process_id);
620  pcs.computeSecondaryVariable(t, x, process_id);
621  }
622 }
623 
625  const double t, const double dt, const std::size_t timestep_id)
626 {
628 
629  NumLib::NonlinearSolverStatus nonlinear_solver_status;
630  for (auto& process_data : _per_process_data)
631  {
632  auto const process_id = process_data->process_id;
633  nonlinear_solver_status = solveMonolithicProcess(
634  t, dt, timestep_id, *process_data, _process_solutions, *_output);
635 
636  process_data->nonlinear_solver_status = nonlinear_solver_status;
637  if (!nonlinear_solver_status.error_norms_met)
638  {
639  ERR("The nonlinear solver failed in time step #%u at t = %g s for "
640  "process #%u.",
641  timestep_id, t, process_id);
642 
643  if (!process_data->timestepper->canReduceTimestepSize())
644  {
645  // save unsuccessful solution
646  _output->doOutputAlways(process_data->process, process_id,
647  timestep_id, t, _process_solutions);
648  OGS_FATAL(timestepper_cannot_reduce_dt.data());
649  }
650 
651  return nonlinear_solver_status;
652  }
653  }
654 
655  postTimestepForAllProcesses(t, dt, _per_process_data, _process_solutions);
656 
657  return nonlinear_solver_status;
658 }
659 
662  const double t, const double dt, const std::size_t timestep_id)
663 {
664  // Coupling iteration
666  {
667  // Set the flag of the first iteration be true.
668  for (auto& conv_crit : _global_coupling_conv_crit)
669  {
670  conv_crit->preFirstIteration();
671  }
672  }
673  auto resetCouplingConvergenceCriteria = [&]() {
674  for (auto& conv_crit : _global_coupling_conv_crit)
675  {
676  conv_crit->reset();
677  }
678  };
679 
681 
682  NumLib::NonlinearSolverStatus nonlinear_solver_status{false, -1};
683  bool coupling_iteration_converged = true;
684  for (int global_coupling_iteration = 0;
685  global_coupling_iteration < _global_coupling_max_iterations;
686  global_coupling_iteration++, resetCouplingConvergenceCriteria())
687  {
688  // TODO(wenqing): use process name
689  coupling_iteration_converged = true;
690  int const last_process_id = _per_process_data.size() - 1;
691  for (auto& process_data : _per_process_data)
692  {
693  auto const process_id = process_data->process_id;
694  BaseLib::RunTime time_timestep_process;
695  time_timestep_process.start();
696 
697  CoupledSolutionsForStaggeredScheme coupled_solutions(
699 
700  process_data->process.setCoupledSolutionsForStaggeredScheme(
701  &coupled_solutions);
702 
703  nonlinear_solver_status =
705  dt, *process_data, *_output);
706  process_data->nonlinear_solver_status = nonlinear_solver_status;
707 
708  INFO(
709  "[time] Solving process #%u took %g s in time step #%u "
710  " coupling iteration #%u",
711  process_id, time_timestep_process.elapsed(), timestep_id,
712  global_coupling_iteration);
713 
714  if (!nonlinear_solver_status.error_norms_met)
715  {
716  ERR("The nonlinear solver failed in time step #%u at t = %g s "
717  "for process #%u.",
718  timestep_id, t, process_id);
719 
720  if (!process_data->timestepper->canReduceTimestepSize())
721  {
722  // save unsuccessful solution
723  _output->doOutputAlways(process_data->process, process_id,
724  timestep_id, t, _process_solutions);
725  OGS_FATAL(timestepper_cannot_reduce_dt.data());
726  }
727  break;
728  }
729 
730  // Check the convergence of the coupling iteration
731  auto& x = *_process_solutions[process_id];
732  auto& x_old = *_solutions_of_last_cpl_iteration[process_id];
733  if (global_coupling_iteration > 0)
734  {
735  MathLib::LinAlg::axpy(x_old, -1.0, x); // save dx to x_old
736  if (process_id == last_process_id)
737  {
738  INFO(
739  "------- Checking convergence criterion for coupled "
740  "solution -------");
741  _global_coupling_conv_crit[process_id]->checkDeltaX(x_old,
742  x);
743  coupling_iteration_converged =
744  coupling_iteration_converged &&
745  _global_coupling_conv_crit[process_id]->isSatisfied();
746  }
747  }
748  MathLib::LinAlg::copy(x, x_old);
749  } // end of for (auto& process_data : _per_process_data)
750 
751  if (coupling_iteration_converged && global_coupling_iteration > 0)
752  {
753  break;
754  }
755 
756  if (!nonlinear_solver_status.error_norms_met)
757  {
758  return nonlinear_solver_status;
759  }
760  }
761 
762  if (!coupling_iteration_converged)
763  {
764  WARN(
765  "The coupling iterations reaches its maximum number in time step"
766  "#%u at t = %g s",
767  timestep_id, t);
768  }
769 
770  if (_chemical_system != nullptr)
771  {
772  // Sequential non-iterative approach applied here to perform water
773  // chemistry calculation followed by resolving component transport
774  // process.
775  // TODO: move into a global loop to consider both mass balance over
776  // space and localized chemical equilibrium between solutes.
777  BaseLib::RunTime time_phreeqc;
778  time_phreeqc.start();
779  _chemical_system->doWaterChemistryCalculation(_process_solutions, dt);
780  INFO("[time] Phreeqc took %g s.", time_phreeqc.elapsed());
781  }
782 
784 
785  return nonlinear_solver_status;
786 }
787 
788 template <typename OutputClass, typename OutputClassMember>
789 void TimeLoop::outputSolutions(bool const output_initial_condition,
790  unsigned timestep, const double t,
791  OutputClass& output_object,
792  OutputClassMember output_class_member) const
793 {
794  // All _per_process_data share the first process.
795  bool const is_staggered_coupling =
797 
798  for (auto& process_data : _per_process_data)
799  {
800  auto const process_id = process_data->process_id;
801  auto& pcs = process_data->process;
802  // If nonlinear solver diverged, the solution has already been
803  // saved.
804  if (!process_data->nonlinear_solver_status.error_norms_met)
805  {
806  continue;
807  }
808 
809  auto const& x = *_process_solutions[process_id];
810 
811  if (output_initial_condition)
812  {
813  pcs.preTimestep(_process_solutions, _start_time,
814  process_data->timestepper->getTimeStep().dt(),
815  process_id);
816  // Update secondary variables, which might be uninitialized, before
817  // output.
818  pcs.computeSecondaryVariable(_start_time, x, process_id);
819  }
820  if (is_staggered_coupling)
821  {
822  CoupledSolutionsForStaggeredScheme coupled_solutions(
824 
825  process_data->process.setCoupledSolutionsForStaggeredScheme(
826  &coupled_solutions);
827  process_data->process
828  .setCoupledTermForTheStaggeredSchemeToLocalAssemblers(
829  process_id);
830  }
831  (output_object.*output_class_member)(pcs, process_id, timestep, t,
833  }
834 }
835 
837 {
838  for (auto* x : _process_solutions)
839  {
841  }
842 
843  for (auto* x : _solutions_of_last_cpl_iteration)
844  {
846  }
847 }
848 
849 } // namespace ProcessLib
const double _start_time
Definition: TimeLoop.h:119
void setEquationSystem(NumLib::NonlinearSolverBase &nonlinear_solver, NumLib::EquationSystem &eq_sys, NumLib::ConvergenceCriterion &conv_crit, NumLib::NonlinearSolverTag nl_tag)
Definition: TimeLoop.cpp:48
std::unique_ptr< NumLib::TimeDiscretization > time_disc
Definition: ProcessData.h:67
void doOutput(Process const &process, const int process_id, const int timestep, const double t, std::vector< GlobalVector *> const &x)
Definition: Output.cpp:304
void doOutputNonlinearIteration(Process const &process, const int process_id, const int timestep, const double t, std::vector< GlobalVector *> const &x, const int iteration)
Definition: Output.cpp:336
void finalizeAssembly(Matrix &)
NumLib::NonlinearSolverStatus solveOneTimeStepOneProcess(std::vector< GlobalVector *> &x, std::size_t const timestep, double const t, double const delta_t, ProcessData const &process_data, Output &output_control)
Definition: TimeLoop.cpp:196
std::unique_ptr< NumLib::EquationSystem > tdisc_ode_sys
type-erased time-discretized ODE system
Definition: ProcessData.h:69
const int _global_coupling_max_iterations
Maximum iterations of the global coupling.
Definition: TimeLoop.h:123
VecNormType
Norm type. Not declared as class type in order to use the members as integers.
Definition: LinAlgEnums.h:21
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:624
std::vector< GlobalVector * > setInitialConditions(double const t0, std::vector< std::unique_ptr< ProcessData >> const &per_process_data)
Definition: TimeLoop.cpp:127
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:565
virtual GlobalVector & getVector()=0
Get an uninitialized vector.
bool isMonolithicSchemeUsed() const
Definition: Process.h:95
int _repeating_times_of_rejected_step
Definition: TimeLoop.h:118
void setCoupledSolutions()
Definition: TimeLoop.cpp:253
std::vector< GlobalVector * > _process_solutions
Definition: TimeLoop.h:113
NumLib::NonlinearSolverBase & nonlinear_solver
Definition: ProcessData.h:63
double computeTimeStepping(const double prev_dt, double &t, std::size_t &accepted_steps, std::size_t &rejected_steps)
Definition: TimeLoop.cpp:268
void doOutputLastTimestep(Process const &process, const int process_id, const int timestep, const double t, std::vector< GlobalVector *> const &x)
Definition: Output.cpp:321
static NUMLIB_EXPORT VectorProvider & provider
Definition of the RunTime class.
std::vector< GlobalVector * > _solutions_of_last_cpl_iteration
Definition: TimeLoop.h:132
NonlinearSolverTag
Tag used to specify which nonlinear solver will be used.
Definition: Types.h:20
const double _end_time
Definition: TimeLoop.h:120
std::unique_ptr< Output > _output
Definition: TimeLoop.h:114
NumLib::NonlinearSolverTag const nonlinear_solver_tag
Definition: ProcessData.h:62
Status of the non-linear solver.
std::unique_ptr< NumLib::ConvergenceCriterion > conv_crit
Definition: ProcessData.h:65
bool isMonolithicProcess(ProcessLib::ProcessData const &process_data)
Definition: TimeLoop.cpp:65
std::vector< std::unique_ptr< NumLib::ConvergenceCriterion > > _global_coupling_conv_crit
Convergence criteria of processes for the global coupling iterations.
Definition: TimeLoop.h:126
NumLib::InternalMatrixStorage * mat_strg
cast of tdisc_ode_sys to NumLib::InternalMatrixStorage
Definition: ProcessData.h:71
std::vector< std::unique_ptr< ProcessData > > _per_process_data
Definition: TimeLoop.h:115
std::unique_ptr< ChemistryLib::ChemicalSolverInterface > _chemical_system
Definition: TimeLoop.h:128
void setTimeDiscretizedODESystem(ProcessData &process_data, NumLib::ODESystem< ODETag, NumLib::NonlinearSolverTag::Picard > &ode_sys)
Definition: TimeLoop.cpp:74
Count the running time.
Definition: RunTime.h:31
static constexpr std::string_view timestepper_cannot_reduce_dt
Definition: TimeLoop.cpp:595
void start()
Start the timer.
Definition: RunTime.h:35
void calculateNonEquilibriumInitialResiduum(std::vector< std::unique_ptr< ProcessData >> const &per_process_data, std::vector< GlobalVector *> process_solutions)
Definition: TimeLoop.cpp:170
void postTimestepForAllProcesses(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:598
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:661
#define OGS_FATAL(fmt,...)
Definition: Error.h:64
virtual void releaseVector(GlobalVector const &x)=0
void initialize()
initialize output, convergence criterion, etc.
Definition: TimeLoop.cpp:406
void outputSolutions(bool const output_initial_condition, unsigned timestep, const double t, OutputClass &output_object, OutputClassMember output_class_member) const
Definition: TimeLoop.cpp:789
void axpy(MatrixOrVector &y, double const a, MatrixOrVector const &x)
Computes .
Definition: LinAlg.h:58
void copy(MatrixOrVector const &x, MatrixOrVector &y)
Copies x to y.
Definition: LinAlg.h:37
double elapsed() const
Get the elapsed time after started.
Definition: RunTime.h:51
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, std::unique_ptr< ChemistryLib::ChemicalSolverInterface > &&chemical_system, const double start_time, const double end_time)
Definition: TimeLoop.cpp:235
static NumLib::NonlinearSolverStatus solveMonolithicProcess(const double t, const double dt, const std::size_t timestep_id, ProcessData const &process_data, std::vector< GlobalVector *> &x, Output &output)
Definition: TimeLoop.cpp:578