13#include <range/v3/algorithm/any_of.hpp>
14#include <range/v3/algorithm/contains.hpp>
15#include <range/v3/view/map.hpp>
33 std::vector<std::unique_ptr<ProcessLib::ProcessData>>
const&
37 for (
auto& process_data : per_process_data)
39 process_data->process.updateDeactivatedSubdomains(
40 t, process_data->process_id);
45 const int timestep,
const double t)
47 return ranges::any_of(outputs, [timestep, t](
auto const& output)
48 {
return output.isOutputStep(timestep, t); });
52 int const timestep,
double const t,
double const dt,
53 std::vector<std::unique_ptr<ProcessLib::ProcessData>>
const&
55 std::vector<GlobalVector*>
const& process_solutions,
56 std::vector<GlobalVector*>
const& process_solutions_prev,
57 std::vector<ProcessLib::Output>
const& outputs)
64 for (
auto& process_data : per_process_data)
66 auto const process_id = process_data->process_id;
67 auto& pcs = process_data->process;
69 pcs.preOutput(t, dt, process_solutions, process_solutions_prev,
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)
82 for (
auto& process_data : per_process_data)
84 auto const process_id = process_data->process_id;
85 auto& pcs = process_data->process;
86 pcs.preTimestep(_process_solutions, t, dt, process_id);
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)
96 for (
auto& process_data : per_process_data)
98 auto const process_id = process_data->process_id;
99 auto& pcs = process_data->process;
101 pcs.computeSecondaryVariable(t, dt, process_solutions,
102 *process_solutions_prev[process_id],
104 pcs.postTimestep(process_solutions, process_solutions_prev, t, dt,
109template <NumLib::ODESystemTag ODETag>
133 else if ((
dynamic_cast<NonlinearSolverNewton*
>(
144 if (
auto* ode_newton =
dynamic_cast<ODENewton*
>(&ode_sys))
153 "You are trying to solve a non-Newton-ready ODE with the"
154 " Newton-Raphson method. Aborting");
159 OGS_FATAL(
"Encountered unknown nonlinear solver type. Aborting");
168std::pair<std::vector<GlobalVector*>, std::vector<GlobalVector*>>
171 std::vector<std::unique_ptr<ProcessData>>
const& per_process_data)
173 std::vector<GlobalVector*> process_solutions;
174 std::vector<GlobalVector*> process_solutions_prev;
176 for (
auto const& process_data : per_process_data)
178 auto const process_id = process_data->process_id;
179 auto& ode_sys = *process_data->tdisc_ode_sys;
182 process_solutions.emplace_back(
184 ode_sys.getMatrixSpecifications(process_id)));
185 process_solutions_prev.emplace_back(
187 ode_sys.getMatrixSpecifications(process_id)));
190 for (
auto const& process_data : per_process_data)
192 auto& pcs = process_data->process;
193 auto const process_id = process_data->process_id;
194 pcs.setInitialConditions(process_solutions, process_solutions_prev, t0,
197 auto& time_disc = *process_data->time_disc;
198 time_disc.setInitialState(t0);
201 return {process_solutions, process_solutions_prev};
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)
209 for (
auto const& process_data : per_process_data)
211 auto& nonlinear_solver = process_data->nonlinear_solver;
214 nonlinear_solver.calculateNonEquilibriumInitialResiduum(
215 process_solutions, process_solutions_prev,
216 process_data->process_id);
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)
225 auto& process = process_data.
process;
226 int const process_id = process_data.
process_id;
227 auto& time_disc = *process_data.
time_disc;
237 time_disc.nextTimestep(t, delta_t);
239 auto const post_iteration_callback =
240 [&](
int iteration, std::vector<GlobalVector*>
const& x)
245 for (
auto const& output : outputs)
247 output.doOutputNonlinearIteration(process, process_id, timestep, t,
252 auto const nonlinear_solver_status =
253 nonlinear_solver.solve(x, x_prev, post_iteration_callback, process_id);
255 if (!nonlinear_solver_status.error_norms_met)
257 return nonlinear_solver_status;
260 process.postNonLinearSolver(x, x_prev, t, delta_t, process_id);
262 return nonlinear_solver_status;
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),
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))
302 if (time == timestep_algorithm.
begin())
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)
315 bool all_process_steps_accepted =
true;
317 double dt = std::numeric_limits<double>::max();
318 constexpr double eps = std::numeric_limits<double>::epsilon();
320 bool const is_initial_step =
322 [](
auto const& ppd) ->
bool
323 { return ppd->timestep_current.timeStepNumber() == 0; });
328 auto& timestep_algorithm = *ppd.timestep_algorithm.get();
333 const double solution_error =
337 ppd.conv_crit.get() ? ppd.conv_crit->getVectorNormType()
341 ppd.timestep_current.setAccepted(
342 ppd.nonlinear_solver_status.error_norms_met);
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);
348 if (!previous_step_accepted &&
352 t + eps < timestep_algorithm.end())
355 all_process_steps_accepted =
false;
358 if (!ppd.nonlinear_solver_status.error_norms_met)
361 "Time step will be rejected due to nonlinear solver "
363 all_process_steps_accepted =
false;
366 if (timestepper_dt > eps ||
367 std::abs(t - timestep_algorithm.end()) < eps)
369 dt = std::min(timestepper_dt, dt);
373 if (all_process_steps_accepted)
382 bool last_step_rejected =
false;
383 if (!is_initial_step)
385 if (all_process_steps_accepted)
388 last_step_rejected =
false;
396 last_step_rejected =
true;
402 for (
auto const& time_step_constain : time_step_constraints)
404 dt = std::min(dt, time_step_constain(t, dt));
408 if (std::abs(dt - prev_dt) < eps)
410 if (last_step_rejected)
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.",
421 DBUG(
"The time stepping is stabilized with the step size of {:g}.",
430 if (all_process_steps_accepted)
434 ppd.timestep_current);
435 auto& timestep_algorithm = ppd.timestep_algorithm;
436 timestep_algorithm->resetCurrentTimeStep(dt, ppd.timestep_previous,
437 ppd.timestep_current);
442 if (all_process_steps_accepted)
451 "Time step {:d} was rejected {:d} times and it will be "
452 "repeated with a reduced step size.",
459 return {dt, last_step_rejected};
462std::vector<std::function<double(
double,
double)>>
464 std::vector<double>&& fixed_times)
const
466 std::vector<std::function<double(
double,
double)>>
const
467 time_step_constraints{
468 [fixed_times = std::move(fixed_times)](
double t,
double dt) {
472 [
this](
double t,
double dt)
480 return time_step_constraints;
488 auto& pcs = process_data->process;
491 output.addProcess(pcs);
496 if (
auto* conv_crit =
498 process_data->conv_crit.get()))
500 int const process_id = process_data->process_id;
501 conv_crit->setDOFTable(pcs.getDOFTable(process_id), pcs.getMesh());
510 bool const is_staggered_coupling =
513 if (is_staggered_coupling)
540 time_timestep.
start();
547 "=== Time stepping at step #{:d} and time {:.15g} with step size "
555 INFO(
"[time] Time step #{:d} took {:g} s.", timesteps,
562 const double prev_dt =
_dt;
581 std::numeric_limits<double>::epsilon() ||
587 if (
_dt < std::numeric_limits<double>::epsilon())
590 "Time step size of {:g} is too small.\n"
591 "Time stepping stops at step {:d} and at time of {:g}.",
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",
615 std::size_t
const timesteps)
620 bool const is_staggered_coupling =
624 if (is_staggered_coupling)
626 nonlinear_solver_status =
631 nonlinear_solver_status =
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)
661 time_timestep_process.
start();
664 x, x_prev, timestep_id, t, dt, process_data, outputs);
666 INFO(
"[time] Solving process #{:d} took {:g} s in time step #{:d}",
669 return nonlinear_solver_status;
673 "Time stepper cannot reduce the time step size further.";
676 const double t,
const double dt,
const std::size_t timestep_id)
682 auto const process_id = process_data->process_id;
687 process_data->nonlinear_solver_status = nonlinear_solver_status;
690 ERR(
"The nonlinear solver failed in time step #{:d} at t = {:g} s "
691 "for process #{:d}.",
692 timestep_id, t, process_id);
694 if (!process_data->timestep_algorithm->canReduceTimestepSize(
695 process_data->timestep_current,
696 process_data->timestep_previous))
701 output.doOutputAlways(
702 process_data->process, process_id, timestep_id, t,
703 process_data->nonlinear_solver_status.number_iterations,
709 return nonlinear_solver_status;
713 return nonlinear_solver_status;
718 const double t,
const double dt,
const std::size_t timestep_id)
726 conv_crit->preFirstIteration();
729 auto resetCouplingConvergenceCriteria = [&]()
738 bool coupling_iteration_converged =
true;
739 bool local_coupling_iteration_converged =
true;
740 for (
int global_coupling_iteration = 0;
742 global_coupling_iteration++, resetCouplingConvergenceCriteria())
745 coupling_iteration_converged =
true;
746 bool local_iteration_converged =
true;
749 auto const process_id = process_data->process_id;
751 bool const isLocalCouplingProcess = ranges::contains(
754 if (!local_coupling_iteration_converged && !isLocalCouplingProcess)
756 coupling_iteration_converged =
false;
761 time_timestep_process.
start();
766 process_data->nonlinear_solver_status = nonlinear_solver_status;
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);
774 if (!nonlinear_solver_status.error_norms_met)
777 "The nonlinear solver failed in time step #{:d} at t = "
778 "{:g} s for process #{:d}.",
779 timestep_id, t, process_id);
781 return nonlinear_solver_status;
787 if (global_coupling_iteration > 0)
791 "------- Checking convergence criterion for coupled "
792 "solution of process #{:d} -------",
795 coupling_iteration_converged =
796 coupling_iteration_converged &&
798 if (isLocalCouplingProcess)
800 local_iteration_converged =
801 local_iteration_converged &&
808 local_coupling_iteration_converged = local_iteration_converged;
809 if (local_coupling_iteration_converged &&
810 coupling_iteration_converged && global_coupling_iteration > 0)
815 if (!nonlinear_solver_status.error_norms_met)
817 return nonlinear_solver_status;
821 if (!coupling_iteration_converged || !local_coupling_iteration_converged)
824 "The coupling iterations reaches its maximum number in time step "
825 "#{:d} at t = {:g} s",
832 auto& pcs = process_data->process;
833 int const process_id = process_data->process_id;
834 auto& ode_sys = *process_data->tdisc_ode_sys;
841 return nonlinear_solver_status;
844template <
typename OutputClassMember>
846 OutputClassMember output_class_member)
const
852 if (!process_data->nonlinear_solver_status.error_norms_met)
857 auto const process_id = process_data->process_id;
858 auto const& pcs = process_data->process;
860 for (
auto const& output_object :
_outputs)
862 (output_object.*output_class_member)(
863 pcs, process_id, timestep, t,
864 process_data->nonlinear_solver_status.number_iterations,
893 if (!process_data->nonlinear_solver_status.error_norms_met)
898 auto const process_id = process_data->process_id;
899 auto& pcs = process_data->process;
904 process_data->time_disc->nextTimestep(t, dt);
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
void WARN(fmt::format_string< Args... > fmt, Args &&... args)
Definition of the RunTime class.
double elapsed() const
Get the elapsed time in seconds.
void start()
Start the timer.
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
void doOutput(Process const &process, const int process_id, int const timestep, const double t, int const iteration, std::vector< GlobalVector * > const &xs) const
virtual bool isMonolithicSchemeUsed() const
std::vector< std::function< double(double, double)> > generateOutputTimeStepConstraints(std::vector< double > &&fixed_times) const
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)
void outputLastTimeStep() const
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)
std::vector< GlobalVector * > _solutions_of_last_cpl_iteration
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
void outputSolutions(unsigned timestep, const double t, OutputClassMember output_class_member) const
std::vector< Output > _outputs
std::size_t _accepted_steps
bool successful_time_step
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
void initialize()
initialize output, convergence criterion, etc.
bool preTsNonlinearSolvePostTs(double const t, double const dt, std::size_t const timesteps)
int _repeating_times_of_rejected_step
std::size_t _rejected_steps
bool calculateNextTimeStep()
void preOutputInitialConditions(const double t) const
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.
std::map< std::string, int > _local_coupling_processes
Processes that will be solved in a local iteration.
NonlinearSolverTag
Tag used to specify which nonlinear solver will be used.
void copy(PETScVector const &x, PETScVector &y)
void axpy(PETScVector &y, PetscScalar const a, PETScVector const &x)
void updateTimeSteps(double const dt, TimeStep &previous_timestep, TimeStep ¤t_timestep)
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
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 * > const &process_solutions, std::vector< GlobalVector * > const &process_solutions_prev)
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)
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)
void preTimestepForAllProcesses(double const t, double const dt, std::vector< std::unique_ptr< ProcessData > > const &per_process_data, std::vector< GlobalVector * > const &_process_solutions)
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)
void setEquationSystem(ProcessData const &process_data)
bool computationOfChangeNeeded(NumLib::TimeStepAlgorithm const ×tep_algorithm, double const time)
std::vector< double > calculateUniqueFixedTimesForAllOutputs(std::vector< Output > const &outputs)
std::pair< std::vector< GlobalVector * >, std::vector< GlobalVector * > > setInitialConditions(double const t0, std::vector< std::unique_ptr< ProcessData > > const &per_process_data)
bool isMonolithicProcess(ProcessLib::ProcessData const &process_data)
bool isOutputStep(std::vector< ProcessLib::Output > const &outputs, const int timestep, const double t)
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)
void updateDeactivatedSubdomains(std::vector< std::unique_ptr< ProcessLib::ProcessData > > const &per_process_data, double const t)
static NUMLIB_EXPORT VectorProvider & provider
Status of the non-linear solver.
std::unique_ptr< NumLib::TimeDiscretization > time_disc
NumLib::NonlinearSolverBase & nonlinear_solver
std::unique_ptr< NumLib::EquationSystem > tdisc_ode_sys
type-erased time-discretized ODE system