13#include <range/v3/algorithm/any_of.hpp>
14#include <range/v3/algorithm/contains.hpp>
28 std::vector<std::unique_ptr<ProcessLib::ProcessData>>
const&
32 for (
auto& process_data : per_process_data)
34 process_data->process.updateDeactivatedSubdomains(
35 t, process_data->process_id);
40 const int timestep,
const double t,
const double end_time)
42 if (std::abs(end_time - t) < std::numeric_limits<double>::epsilon())
48 return ranges::any_of(outputs, [timestep, t](
auto const& output)
49 {
return output.isOutputStep(timestep, t); });
53 int const timestep,
double const t,
double const dt,
const double end_time,
54 std::vector<std::unique_ptr<ProcessLib::ProcessData>>
const&
56 std::vector<GlobalVector*>
const& process_solutions,
57 std::vector<GlobalVector*>
const& process_solutions_prev,
58 std::vector<ProcessLib::Output>
const& outputs)
65 for (
auto& process_data : per_process_data)
67 auto const process_id = process_data->process_id;
68 auto& pcs = process_data->process;
70 pcs.preOutput(t, dt, process_solutions, process_solutions_prev,
79 double const t,
double const dt,
80 std::vector<std::unique_ptr<ProcessData>>
const& per_process_data,
81 std::vector<GlobalVector*>
const& _process_solutions)
83 for (
auto& process_data : per_process_data)
85 auto const process_id = process_data->process_id;
86 auto& pcs = process_data->process;
87 pcs.preTimestep(_process_solutions, t, dt, process_id);
92 double const t,
double const dt,
93 std::vector<std::unique_ptr<ProcessData>>
const& per_process_data,
94 std::vector<GlobalVector*>
const& process_solutions,
95 std::vector<GlobalVector*>
const& process_solutions_prev)
97 for (
auto& process_data : per_process_data)
99 auto const process_id = process_data->process_id;
100 auto& pcs = process_data->process;
102 pcs.computeSecondaryVariable(t, dt, process_solutions,
103 *process_solutions_prev[process_id],
105 pcs.postTimestep(process_solutions, process_solutions_prev, t, dt,
110template <NumLib::ODESystemTag ODETag>
134 else if ((
dynamic_cast<NonlinearSolverNewton*
>(
145 if (
auto* ode_newton =
dynamic_cast<ODENewton*
>(&ode_sys))
154 "You are trying to solve a non-Newton-ready ODE with the"
155 " Newton-Raphson method. Aborting");
160 OGS_FATAL(
"Encountered unknown nonlinear solver type. Aborting");
169std::pair<std::vector<GlobalVector*>, std::vector<GlobalVector*>>
172 std::vector<std::unique_ptr<ProcessData>>
const& per_process_data)
174 std::vector<GlobalVector*> process_solutions;
175 std::vector<GlobalVector*> process_solutions_prev;
177 for (
auto const& process_data : per_process_data)
179 auto const process_id = process_data->process_id;
180 auto& ode_sys = *process_data->tdisc_ode_sys;
183 process_solutions.emplace_back(
185 ode_sys.getMatrixSpecifications(process_id)));
186 process_solutions_prev.emplace_back(
188 ode_sys.getMatrixSpecifications(process_id)));
191 for (
auto const& process_data : per_process_data)
193 auto& pcs = process_data->process;
194 auto const process_id = process_data->process_id;
195 pcs.setInitialConditions(process_solutions, process_solutions_prev, t0,
198 auto& time_disc = *process_data->time_disc;
199 time_disc.setInitialState(t0);
202 return {process_solutions, process_solutions_prev};
206 std::vector<std::unique_ptr<ProcessData>>
const& per_process_data,
207 std::vector<GlobalVector*>
const& process_solutions,
208 std::vector<GlobalVector*>
const& process_solutions_prev)
210 for (
auto const& process_data : per_process_data)
212 auto& nonlinear_solver = process_data->nonlinear_solver;
215 nonlinear_solver.calculateNonEquilibriumInitialResiduum(
216 process_solutions, process_solutions_prev,
217 process_data->process_id);
222 std::vector<GlobalVector*>& x, std::vector<GlobalVector*>
const& x_prev,
223 std::size_t
const timestep,
double const t,
double const delta_t,
224 ProcessData const& process_data, std::vector<Output>
const& outputs)
226 auto& process = process_data.
process;
227 int const process_id = process_data.
process_id;
228 auto& time_disc = *process_data.
time_disc;
238 time_disc.nextTimestep(t, delta_t);
240 auto const post_iteration_callback =
241 [&](
int iteration, std::vector<GlobalVector*>
const& x)
246 for (
auto const& output : outputs)
248 output.doOutputNonlinearIteration(process, process_id, timestep, t,
253 auto const nonlinear_solver_status =
254 nonlinear_solver.solve(x, x_prev, post_iteration_callback, process_id);
256 if (!nonlinear_solver_status.error_norms_met)
258 return nonlinear_solver_status;
261 process.postNonLinearSolver(x, x_prev, t, delta_t, process_id);
263 return nonlinear_solver_status;
267 std::vector<Output>&& outputs,
268 std::vector<std::unique_ptr<ProcessData>>&& per_process_data,
269 std::unique_ptr<NumLib::StaggeredCoupling>&& staggered_coupling,
270 const double start_time,
const double end_time)
271 : _outputs{std::move(outputs)},
272 _per_process_data(std::move(per_process_data)),
273 _start_time(start_time),
275 _staggered_coupling(std::move(staggered_coupling))
284 if (time == timestep_algorithm.
begin())
292 const double prev_dt,
double& t, std::size_t& accepted_steps,
293 std::size_t& rejected_steps,
294 std::vector<std::function<
double(
double,
double)>>
const&
295 time_step_constraints)
297 bool all_process_steps_accepted =
true;
299 double dt = std::numeric_limits<double>::max();
300 constexpr double eps = std::numeric_limits<double>::epsilon();
302 bool const is_initial_step =
304 [](
auto const& ppd) ->
bool
305 { return ppd->timestep_current.timeStepNumber() == 0; });
310 auto& timestep_algorithm = *ppd.timestep_algorithm.get();
315 const double solution_error =
319 ppd.conv_crit.get() ? ppd.conv_crit->getVectorNormType()
323 ppd.timestep_current.setAccepted(
324 ppd.nonlinear_solver_status.error_norms_met);
326 auto [previous_step_accepted, timestepper_dt] = timestep_algorithm.next(
327 solution_error, ppd.nonlinear_solver_status.number_iterations,
328 ppd.timestep_previous, ppd.timestep_current);
330 if (!previous_step_accepted &&
334 t + eps < timestep_algorithm.end())
337 all_process_steps_accepted =
false;
340 if (!ppd.nonlinear_solver_status.error_norms_met)
343 "Time step will be rejected due to nonlinear solver "
345 all_process_steps_accepted =
false;
348 if (timestepper_dt > eps ||
349 std::abs(t - timestep_algorithm.end()) < eps)
351 dt = std::min(timestepper_dt, dt);
355 if (all_process_steps_accepted)
364 bool last_step_rejected =
false;
365 if (!is_initial_step)
367 if (all_process_steps_accepted)
370 last_step_rejected =
false;
378 last_step_rejected =
true;
384 for (
auto const& time_step_constain : time_step_constraints)
386 dt = std::min(dt, time_step_constain(t, dt));
390 if (std::abs(dt - prev_dt) < eps)
392 if (last_step_rejected)
395 "The new step size of {:g} is the same as that of the previous "
396 "rejected time step. \nPlease re-run ogs with a proper "
397 "adjustment in the numerical settings, \ne.g those for time "
398 "stepper, local or global non-linear solver.",
403 DBUG(
"The time stepping is stabilized with the step size of {:g}.",
412 if (all_process_steps_accepted)
416 ppd.timestep_current);
417 auto& timestep_algorithm = ppd.timestep_algorithm;
418 timestep_algorithm->resetCurrentTimeStep(dt, ppd.timestep_previous,
419 ppd.timestep_current);
424 if (all_process_steps_accepted)
433 "Time step {:d} was rejected {:d} times and it will be "
434 "repeated with a reduced step size.",
441 return {dt, last_step_rejected};
444std::vector<std::function<double(
double,
double)>>
446 std::vector<double>&& fixed_times)
const
448 std::vector<std::function<double(
double,
double)>>
const
449 time_step_constraints{
450 [fixed_times = std::move(fixed_times)](
double t,
double dt) {
454 [
this](
double t,
double dt)
462 return time_step_constraints;
470 auto& pcs = process_data->process;
473 output.addProcess(pcs);
478 if (
auto* conv_crit =
480 process_data->conv_crit.get()))
482 int const process_id = process_data->process_id;
483 conv_crit->setDOFTable(pcs.getDOFTable(process_id), pcs.getMesh());
518 time_timestep.
start();
525 "=== Time stepping at step #{:d} and time {:.15g} with step size "
533 INFO(
"[time] Time step #{:d} took {:g} s.", timesteps,
540 const double prev_dt =
_dt;
559 std::numeric_limits<double>::epsilon() ||
565 if (
_dt < std::numeric_limits<double>::epsilon())
568 "Time step size of {:g} is too small.\n"
569 "Time stepping stops at step {:d} and at time of {:g}.",
580 "The whole computation of the time stepping took {:d} steps, in which\n"
581 "\t the accepted steps are {:d}, and the rejected steps are {:d}.\n",
593 std::size_t
const timesteps)
601 nonlinear_solver_status =
606 nonlinear_solver_status =
630 const double t,
const double dt,
const std::size_t timestep_id,
631 ProcessData const& process_data, std::vector<GlobalVector*>& x,
632 std::vector<GlobalVector*>
const& x_prev,
633 std::vector<Output>
const& outputs)
636 time_timestep_process.
start();
639 x, x_prev, timestep_id, t, dt, process_data, outputs);
641 INFO(
"[time] Solving process #{:d} took {:g} s in time step #{:d}",
644 return nonlinear_solver_status;
648 "Time stepper cannot reduce the time step size further.";
651 const double t,
const double dt,
const std::size_t timestep_id)
657 auto const process_id = process_data->process_id;
662 process_data->nonlinear_solver_status = nonlinear_solver_status;
665 ERR(
"The nonlinear solver failed in time step #{:d} at t = {:g} s "
666 "for process #{:d}.",
667 timestep_id, t, process_id);
669 if (!process_data->timestep_algorithm->canReduceTimestepSize(
670 process_data->timestep_current,
671 process_data->timestep_previous))
676 output.doOutputAlways(
677 process_data->process, process_id, timestep_id, t,
678 process_data->nonlinear_solver_status.number_iterations,
684 return nonlinear_solver_status;
688 return nonlinear_solver_status;
693 const double t,
const double dt,
const std::size_t timestep_id)
695 auto const nonlinear_solver_status =
705 auto& pcs = process_data->process;
706 int const process_id = process_data->process_id;
707 auto& ode_sys = *process_data->tdisc_ode_sys;
714 return nonlinear_solver_status;
717template <
typename OutputClassMember>
719 OutputClassMember output_class_member)
const
725 if (!process_data->nonlinear_solver_status.error_norms_met)
730 auto const process_id = process_data->process_id;
731 auto const& pcs = process_data->process;
733 for (
auto const& output_object :
_outputs)
735 (output_object.*output_class_member)(
736 pcs, process_id, timestep, t,
737 process_data->nonlinear_solver_status.number_iterations,
761 if (!process_data->nonlinear_solver_status.error_norms_met)
766 auto const process_id = process_data->process_id;
767 auto& pcs = process_data->process;
772 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 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
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)
TimeLoop(std::vector< Output > &&outputs, std::vector< std::unique_ptr< ProcessData > > &&per_process_data, std::unique_ptr< NumLib::StaggeredCoupling > &&staggered_coupling, const double start_time, const double end_time)
void outputLastTimeStep() const
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< GlobalVector * > _process_solutions
std::unique_ptr< NumLib::StaggeredCoupling > _staggered_coupling
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.
NonlinearSolverTag
Tag used to specify which nonlinear solver will be used.
void copy(PETScVector const &x, PETScVector &y)
double computeRelativeNorm(VectorType const &x, VectorType const &y, MathLib::VecNormType norm_type)
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)
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)
void updateDeactivatedSubdomains(std::vector< std::unique_ptr< ProcessLib::ProcessData > > const &per_process_data, double const t)
void preOutputForAllProcesses(int const timestep, double const t, double const dt, const double end_time, 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)
bool isOutputStep(std::vector< ProcessLib::Output > const &outputs, const int timestep, const double t, const double end_time)
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