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