36 auto& eq_sys_ =
static_cast<EqSys&
>(eq_sys);
42 nl_solver->setEquationSystem(eq_sys_, conv_crit);
47 "Could not cast nonlinear solver to Picard type solver.");
54 auto& eq_sys_ =
static_cast<EqSys&
>(eq_sys);
61 nl_solver->setEquationSystem(eq_sys_, conv_crit);
64 else if (
auto* nl_solver =
69 nl_solver->setEquationSystem(eq_sys_, conv_crit);
75 "Could not cast nonlinear solver to Newton type solver.");
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)
95 for (
auto& process_data : per_process_data)
97 auto const process_id = process_data->process_id;
98 auto& pcs = process_data->process;
99 pcs.preTimestep(_process_solutions, t, dt, process_id);
104 double const t,
double const dt,
105 std::vector<std::unique_ptr<ProcessData>>
const& per_process_data,
106 std::vector<GlobalVector*>
const& process_solutions,
107 std::vector<GlobalVector*>
const& process_solutions_prev,
108 std::vector<std::size_t>& xdot_vector_ids)
110 std::vector<GlobalVector*> x_dots;
111 x_dots.reserve(per_process_data.size());
112 xdot_vector_ids.resize(per_process_data.size());
115 for (
auto& process_data : per_process_data)
117 auto const process_id = process_data->process_id;
118 auto const& ode_sys = *process_data->tdisc_ode_sys;
119 auto const& time_discretization = *process_data->time_disc;
122 ode_sys.getMatrixSpecifications(process_id), xdot_vector_ids[cnt]));
125 time_discretization.getXdot(*process_solutions[process_id],
126 *process_solutions_prev[process_id],
127 *x_dots[process_id]);
131 bool const is_staggered_coupling =
134 for (
auto& process_data : per_process_data)
136 auto const process_id = process_data->process_id;
137 auto& pcs = process_data->process;
139 if (is_staggered_coupling)
143 pcs.setCoupledSolutionsForStaggeredScheme(&coupled_solutions);
145 auto& x_dot = *x_dots[process_id];
146 pcs.computeSecondaryVariable(t, dt, process_solutions, x_dot,
148 pcs.postTimestep(process_solutions, t, dt, process_id);
150 for (
auto& x_dot : x_dots)
156 template <NumLib::ODESystemTag ODETag>
180 else if ((
dynamic_cast<NonlinearSolverNewton*
>(
191 if (
auto* ode_newton =
dynamic_cast<ODENewton*
>(&ode_sys))
200 "You are trying to solve a non-Newton-ready ODE with the"
201 " Newton-Raphson method. Aborting");
206 OGS_FATAL(
"Encountered unknown nonlinear solver type. Aborting");
215 std::pair<std::vector<GlobalVector*>, std::vector<GlobalVector*>>
218 std::vector<std::unique_ptr<ProcessData>>
const& per_process_data)
220 std::vector<GlobalVector*> process_solutions;
221 std::vector<GlobalVector*> process_solutions_prev;
223 for (
auto& process_data : per_process_data)
225 auto const process_id = process_data->process_id;
226 auto& ode_sys = *process_data->tdisc_ode_sys;
229 process_solutions.emplace_back(
231 ode_sys.getMatrixSpecifications(process_id)));
232 process_solutions_prev.emplace_back(
234 ode_sys.getMatrixSpecifications(process_id)));
237 for (
auto& process_data : per_process_data)
239 auto& pcs = process_data->process;
240 auto const process_id = process_data->process_id;
241 pcs.setInitialConditions(process_solutions, process_solutions_prev, t0,
244 auto& time_disc = *process_data->time_disc;
245 time_disc.setInitialState(t0);
248 return {process_solutions, process_solutions_prev};
252 std::vector<std::unique_ptr<ProcessData>>
const& per_process_data,
253 std::vector<GlobalVector*>
255 std::vector<GlobalVector*>
const& process_solutions_prev)
257 INFO(
"Calculate non-equilibrium initial residuum.");
258 for (
auto& process_data : per_process_data)
260 auto& ode_sys = *process_data->tdisc_ode_sys;
262 auto& time_disc = *process_data->time_disc;
263 auto& nonlinear_solver = process_data->nonlinear_solver;
264 auto& conv_crit = *process_data->conv_crit;
266 auto const nl_tag = process_data->nonlinear_solver_tag;
272 time_disc.nextTimestep(t, dt);
273 nonlinear_solver.calculateNonEquilibriumInitialResiduum(
274 process_solutions, process_solutions_prev,
275 process_data->process_id);
280 std::vector<GlobalVector*>& x, std::vector<GlobalVector*>
const& x_prev,
281 std::size_t
const timestep,
double const t,
double const delta_t,
283 std::size_t& xdot_id)
285 auto& process = process_data.
process;
286 int const process_id = process_data.
process_id;
287 auto& time_disc = *process_data.
time_disc;
288 auto& conv_crit = *process_data.
conv_crit;
300 time_disc.nextTimestep(t, delta_t);
302 auto const post_iteration_callback =
303 [&](
int iteration, std::vector<GlobalVector*>
const& x)
309 auto const nonlinear_solver_status =
310 nonlinear_solver.solve(x, x_prev, post_iteration_callback, process_id);
312 if (!nonlinear_solver_status.error_norms_met)
314 return nonlinear_solver_status;
318 ode_sys.getMatrixSpecifications(process_id), xdot_id);
320 time_disc.getXdot(*x[process_id], *x_prev[process_id], x_dot);
322 process.postNonLinearSolver(*x[process_id], x_dot, t, delta_t, process_id);
325 return nonlinear_solver_status;
329 std::vector<std::unique_ptr<ProcessData>>&& per_process_data,
330 const int global_coupling_max_iterations,
331 std::vector<std::unique_ptr<NumLib::ConvergenceCriterion>>&&
332 global_coupling_conv_crit,
333 const double start_time,
const double end_time)
334 : _output(std::move(output)),
335 _per_process_data(std::move(per_process_data)),
336 _start_time(start_time),
338 _global_coupling_max_iterations(global_coupling_max_iterations),
339 _global_coupling_conv_crit(std::move(global_coupling_conv_crit))
359 std::size_t& accepted_steps,
360 std::size_t& rejected_steps)
362 bool all_process_steps_accepted =
true;
364 double dt = std::numeric_limits<double>::max();
366 bool is_initial_step =
false;
370 const auto& timestepper = ppd.timestepper;
371 if (timestepper->getTimeStep().steps() == 0)
373 is_initial_step =
true;
376 auto& time_disc = ppd.time_disc;
380 auto const& conv_crit = ppd.conv_crit;
382 (conv_crit) ? conv_crit->getVectorNormType()
385 const double solution_error =
386 (timestepper->isSolutionErrorComputationNeeded())
387 ? ((t == timestepper->begin())
389 : time_disc->computeRelativeChangeFromPreviousTimestep(
390 x, x_prev, norm_type))
393 if (!ppd.nonlinear_solver_status.error_norms_met)
395 timestepper->setAcceptedOrNot(
false);
399 timestepper->setAcceptedOrNot(
true);
402 auto [step_accepted, timestepper_dt] = timestepper->next(
403 solution_error, ppd.nonlinear_solver_status.number_iterations);
405 if (!step_accepted &&
408 t + std::numeric_limits<double>::epsilon() < timestepper->end())
411 all_process_steps_accepted =
false;
414 if (!ppd.nonlinear_solver_status.error_norms_met)
417 "Time step will be rejected due to nonlinear solver "
419 all_process_steps_accepted =
false;
422 if (timestepper_dt > std::numeric_limits<double>::epsilon() ||
423 std::abs(t - timestepper->end()) <
424 std::numeric_limits<double>::epsilon())
426 dt = std::min(timestepper_dt, dt);
430 if (all_process_steps_accepted)
439 if (!is_initial_step)
441 if (all_process_steps_accepted)
449 std::numeric_limits<double>::epsilon())
465 _output->getFixedOutputTimes());
467 if (std::abs(dt - prev_dt) < std::numeric_limits<double>::epsilon())
472 "The new step size of {:g} is the same as that of the previous "
473 "rejected time step. \nPlease re-run ogs with a proper "
474 "adjustment in the numerical settings, \ne.g those for time "
475 "stepper, local or global non-linear solver.",
480 DBUG(
"The time stepping is stabilized with the step size of {:g}.",
490 auto& timestepper = ppd.timestepper;
491 if (all_process_steps_accepted)
493 timestepper->resetCurrentTimeStep(dt);
496 if (t == timestepper->begin())
503 if (all_process_steps_accepted)
510 std::numeric_limits<double>::epsilon())
513 "Time step {:d} was rejected {:d} times and it will be "
514 "repeated with a reduced step size.",
529 auto& pcs = process_data->process;
530 int const process_id = process_data->process_id;
535 if (
auto* conv_crit =
537 process_data->conv_crit.get()))
539 conv_crit->setDOFTable(pcs.getDOFTable(process_id), pcs.getMesh());
548 bool const is_staggered_coupling =
551 if (is_staggered_coupling)
558 const bool output_initial_condition =
true;
575 bool const is_staggered_coupling =
578 bool non_equilibrium_initial_residuum_computed =
false;
580 std::size_t accepted_steps = 0;
581 std::size_t rejected_steps = 0;
589 time_timestep.
start();
592 const double prev_dt = dt;
594 const std::size_t timesteps = accepted_steps + 1;
597 "=== Time stepping at step #{:d} and time {:g} with step size {:g}",
603 process_data->process.updateDeactivatedSubdomains(
604 t, process_data->process_id);
607 if (!non_equilibrium_initial_residuum_computed)
611 non_equilibrium_initial_residuum_computed =
true;
617 if (is_staggered_coupling)
619 nonlinear_solver_status =
624 nonlinear_solver_status =
638 INFO(
"[time] Time step #{:d} took {:g} s.", timesteps,
641 double const current_time = t;
647 const bool output_initial_condition =
false;
652 if (std::abs(t -
_end_time) < std::numeric_limits<double>::epsilon() ||
658 if (dt < std::numeric_limits<double>::epsilon())
661 "Time step size of {:g} is too small.\n"
662 "Time stepping stops at step {:d} and at time of {:g}.",
669 "The whole computation of the time stepping took {:d} steps, in which\n"
670 "\t the accepted steps are {:d}, and the rejected steps are {:d}.\n",
671 accepted_steps + rejected_steps, accepted_steps, rejected_steps);
676 const bool output_initial_condition =
false;
678 accepted_steps + rejected_steps, t, *
_output,
686 const double t,
const double dt,
const std::size_t timestep_id,
687 ProcessData const& process_data, std::vector<GlobalVector*>& x,
688 std::vector<GlobalVector*>
const& x_prev,
Output& output,
689 std::size_t& xdot_id)
692 time_timestep_process.
start();
695 x, x_prev, timestep_id, t, dt, process_data, output, xdot_id);
697 INFO(
"[time] Solving process #{:d} took {:g} s in time step #{:d} ",
700 return nonlinear_solver_status;
704 "Time stepper cannot reduce the time step size further.";
707 const double t,
const double dt,
const std::size_t timestep_id)
716 auto const process_id = process_data->process_id;
722 process_data->nonlinear_solver_status = nonlinear_solver_status;
725 ERR(
"The nonlinear solver failed in time step #{:d} at t = {:g} s "
726 "for process #{:d}.",
727 timestep_id, t, process_id);
729 if (!process_data->timestepper->canReduceTimestepSize())
733 process_data->process, process_id, timestep_id, t,
734 process_data->nonlinear_solver_status.number_iterations,
739 return nonlinear_solver_status;
743 return nonlinear_solver_status;
748 const double t,
const double dt,
const std::size_t timestep_id)
756 conv_crit->preFirstIteration();
759 auto resetCouplingConvergenceCriteria = [&]()
768 bool coupling_iteration_converged =
true;
769 for (
int global_coupling_iteration = 0;
771 global_coupling_iteration++, resetCouplingConvergenceCriteria())
774 coupling_iteration_converged =
true;
779 auto const process_id = process_data->process_id;
781 time_timestep_process.
start();
790 process_data->process.setCoupledSolutionsForStaggeredScheme(
797 process_data->nonlinear_solver_status = nonlinear_solver_status;
800 "[time] Solving process #{:d} took {:g} s in time step #{:d} "
801 "coupling iteration #{:d}",
802 process_id, time_timestep_process.
elapsed(), timestep_id,
803 global_coupling_iteration);
805 if (!nonlinear_solver_status.error_norms_met)
808 "The nonlinear solver failed in time step #{:d} at t = "
809 "{:g} s for process #{:d}.",
810 timestep_id, t, process_id);
812 return nonlinear_solver_status;
818 if (global_coupling_iteration > 0)
822 "------- Checking convergence criterion for coupled "
823 "solution of process #{:d} -------",
826 coupling_iteration_converged =
827 coupling_iteration_converged &&
833 if (coupling_iteration_converged && global_coupling_iteration > 0)
838 if (!nonlinear_solver_status.error_norms_met)
840 return nonlinear_solver_status;
844 if (!coupling_iteration_converged)
847 "The coupling iterations reaches its maximum number in time step "
848 "#{:d} at t = {:g} s",
855 auto& pcs = process_data->process;
856 int const process_id = process_data->process_id;
857 auto& ode_sys = *process_data->tdisc_ode_sys;
864 return nonlinear_solver_status;
867 template <
typename OutputClass,
typename OutputClassMember>
869 unsigned timestep,
const double t,
870 OutputClass& output_object,
871 OutputClassMember output_class_member)
const
874 bool const is_staggered_coupling =
881 if (!process_data->nonlinear_solver_status.error_norms_met)
886 auto const process_id = process_data->process_id;
887 auto& pcs = process_data->process;
889 if (!is_staggered_coupling && output_initial_condition)
891 auto const& ode_sys = *process_data->tdisc_ode_sys;
895 process_data->time_disc->nextTimestep(t, dt);
898 ode_sys.getMatrixSpecifications(process_id));
909 if (is_staggered_coupling && output_initial_condition)
914 process_data->process.setCoupledSolutionsForStaggeredScheme(
916 process_data->process
917 .setCoupledTermForTheStaggeredSchemeToLocalAssemblers(
920 auto const& ode_sys = *process_data->tdisc_ode_sys;
924 process_data->time_disc->nextTimestep(t, dt);
927 ode_sys.getMatrixSpecifications(process_id));
938 (output_object.*output_class_member)(
939 pcs, process_id, timestep, t,
940 process_data->nonlinear_solver_status.number_iterations,
void INFO(char const *fmt, Args const &... args)
void ERR(char const *fmt, Args const &... args)
void DBUG(char const *fmt, Args const &... args)
void WARN(char const *fmt, Args const &... args)
Definition of the RunTime class.
double elapsed() const
Get the elapsed time in seconds.
void start()
Start the timer.
Global vector based on Eigen vector.
virtual GlobalVector & getVector(std::size_t &id)=0
Get an uninitialized vector with the given id.
virtual void releaseVector(GlobalVector const &x)=0
void doOutputNonlinearIteration(Process const &process, const int process_id, int const timestep, const double t, const int iteration, std::vector< GlobalVector * > const &x)
void doOutput(Process const &process, const int process_id, int const timestep, const double t, int const iteration, std::vector< GlobalVector * > const &x)
void doOutputLastTimestep(Process const &process, const int process_id, int const timestep, const double t, int const iteration, std::vector< GlobalVector * > const &x)
bool isMonolithicSchemeUsed() const
void outputSolutions(bool const output_initial_condition, unsigned timestep, const double t, OutputClass &output_object, OutputClassMember output_class_member) const
std::vector< GlobalVector * > _solutions_of_last_cpl_iteration
std::vector< std::size_t > _xdot_vector_ids
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,...
std::vector< std::unique_ptr< ProcessData > > _per_process_data
double computeTimeStepping(const double prev_dt, double &t, std::size_t &accepted_steps, std::size_t &rejected_steps)
std::vector< std::unique_ptr< NumLib::ConvergenceCriterion > > _global_coupling_conv_crit
Convergence criteria of processes for the global coupling iterations.
void setCoupledSolutions()
const int _global_coupling_max_iterations
Maximum iterations of the global coupling.
std::vector< GlobalVector * > _process_solutions
std::unique_ptr< Output > _output
void initialize()
initialize output, convergence criterion, etc.
int _repeating_times_of_rejected_step
TimeLoop(std::unique_ptr< Output > &&output, std::vector< std::unique_ptr< ProcessData >> &&per_process_data, const int global_coupling_max_iterations, std::vector< std::unique_ptr< NumLib::ConvergenceCriterion >> &&global_coupling_conv_crit, const double start_time, const double end_time)
std::vector< GlobalVector * > _process_solutions_prev
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.
NonlinearSolverTag
Tag used to specify which nonlinear solver will be used.
void axpy(PETScVector &y, double const a, PETScVector const &x)
void copy(PETScVector const &x, PETScVector &y)
VecNormType
Norm type. Not declared as class type in order to use the members as integers.
double possiblyClampDtToNextFixedTime(double const t, double const dt, std::vector< double > const &fixed_output_times)
void postTimestepForAllProcesses(double const t, double const dt, std::vector< std::unique_ptr< ProcessData >> const &per_process_data, std::vector< GlobalVector * > const &process_solutions, std::vector< GlobalVector * > const &process_solutions_prev, std::vector< std::size_t > &xdot_vector_ids)
static NumLib::NonlinearSolverStatus solveMonolithicProcess(const double t, const double dt, const std::size_t timestep_id, ProcessData const &process_data, std::vector< GlobalVector * > &x, std::vector< GlobalVector * > const &x_prev, Output &output, std::size_t &xdot_id)
static constexpr std::string_view timestepper_cannot_reduce_dt
void setTimeDiscretizedODESystem(ProcessData &process_data, NumLib::ODESystem< ODETag, NumLib::NonlinearSolverTag::Picard > &ode_sys)
void calculateNonEquilibriumInitialResiduum(std::vector< std::unique_ptr< ProcessData >> const &per_process_data, std::vector< GlobalVector * > process_solutions, std::vector< GlobalVector * > const &process_solutions_prev)
void preTimestepForAllProcesses(double const t, double const dt, std::vector< std::unique_ptr< ProcessData >> const &per_process_data, std::vector< GlobalVector * > const &_process_solutions)
std::pair< std::vector< GlobalVector * >, std::vector< GlobalVector * > > setInitialConditions(double const t0, std::vector< std::unique_ptr< ProcessData >> const &per_process_data)
NumLib::NonlinearSolverStatus solveOneTimeStepOneProcess(std::vector< GlobalVector * > &x, std::vector< GlobalVector * > const &x_prev, std::size_t const timestep, double const t, double const delta_t, ProcessData const &process_data, Output &output_control, std::size_t &xdot_id)
bool isMonolithicProcess(ProcessLib::ProcessData const &process_data)
void setEquationSystem(NumLib::NonlinearSolverBase &nonlinear_solver, NumLib::EquationSystem &eq_sys, NumLib::ConvergenceCriterion &conv_crit, NumLib::NonlinearSolverTag nl_tag)
static NUMLIB_EXPORT VectorProvider & provider
Status of the non-linear solver.
std::unique_ptr< NumLib::TimeDiscretization > time_disc
std::unique_ptr< NumLib::ConvergenceCriterion > conv_crit
NumLib::NonlinearSolverBase & nonlinear_solver
std::unique_ptr< NumLib::EquationSystem > tdisc_ode_sys
type-erased time-discretized ODE system
NumLib::NonlinearSolverTag const nonlinear_solver_tag