deal.II version GIT relicensing-3325-gccc124ab5a 2025-05-17 05:00:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
time_stepping.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2014 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_time_stepping_h
16#define dealii_time_stepping_h
17
18
19#include <deal.II/base/config.h>
20
22
23#include <functional>
24#include <vector>
25
27
32namespace TimeStepping
33{
152
153
154
174
175
176
181 template <typename VectorType>
183 {
184 public:
188 virtual ~TimeStepping() = default;
189
200 virtual double
202 std::vector<std::function<VectorType(const double, const VectorType &)>>
203 &F,
204 std::vector<std::function<
205 VectorType(const double, const double, const VectorType &)>> &J_inverse,
206 double t,
207 double delta_t,
208 VectorType &y) = 0;
209
213 struct Status
214 {};
215
219 virtual const Status &
220 get_status() const = 0;
221 };
222
223
224
228 template <typename VectorType>
229 class RungeKutta : public TimeStepping<VectorType>
230 {
231 public:
235 virtual ~RungeKutta() override = default;
236
240 virtual void
242
254 double
256 std::vector<std::function<VectorType(const double, const VectorType &)>>
257 &F,
258 std::vector<std::function<
259 VectorType(const double, const double, const VectorType &)>> &J_inverse,
260 double t,
261 double delta_t,
262 VectorType &y) override;
263
275 virtual double
277 const std::function<VectorType(const double, const VectorType &)> &f,
278 const std::function<
279 VectorType(const double, const double, const VectorType &)>
280 &id_minus_tau_J_inverse,
281 double t,
282 double delta_t,
283 VectorType &y) = 0;
284
285 protected:
289 unsigned int n_stages;
290
294 std::vector<double> b;
295
299 std::vector<double> c;
300
304 std::vector<std::vector<double>> a;
305 };
306
307
308
313 template <typename VectorType>
314 class ExplicitRungeKutta : public RungeKutta<VectorType>
315 {
316 public:
317 using RungeKutta<VectorType>::evolve_one_time_step;
318
325
330
334 void
335 initialize(const runge_kutta_method method) override;
336
350 double
352 const std::function<VectorType(const double, const VectorType &)> &f,
353 const std::function<
354 VectorType(const double, const double, const VectorType &)>
355 &id_minus_tau_J_inverse,
356 double t,
357 double delta_t,
358 VectorType &y) override;
359
367 double
369 const std::function<VectorType(const double, const VectorType &)> &f,
370 double t,
371 double delta_t,
372 VectorType &y);
373
377 struct Status : public TimeStepping<VectorType>::Status
378 {
380 : method(invalid)
381 {}
382
384 };
385
389 const Status &
390 get_status() const override;
391
392 private:
396 void
398 const std::function<VectorType(const double, const VectorType &)> &f,
399 const double t,
400 const double delta_t,
401 const VectorType &y,
402 std::vector<VectorType> &f_stages) const;
403
408 };
409
410
411
417 template <typename VectorType>
418 class LowStorageRungeKutta : public RungeKutta<VectorType>
419 {
420 public:
421 using RungeKutta<VectorType>::evolve_one_time_step;
422
429
434
438 void
439 initialize(const runge_kutta_method method) override;
440
452 double
454 const std::function<VectorType(const double, const VectorType &)> &f,
455 const std::function<
456 VectorType(const double, const double, const VectorType &)>
457 &id_minus_tau_J_inverse,
458 double t,
459 double delta_t,
460 VectorType &y) override;
461
471 double
473 const std::function<VectorType(const double, const VectorType &)> &f,
474 double t,
475 double delta_t,
476 VectorType &solution,
477 VectorType &vec_ri,
478 VectorType &vec_ki);
479
486 void
487 get_coefficients(std::vector<double> &a,
488 std::vector<double> &b,
489 std::vector<double> &c) const;
490
494 struct Status : public TimeStepping<VectorType>::Status
495 {
497 : method(invalid)
498 {}
499
501 };
502
506 const Status &
507 get_status() const override;
508
509 private:
513 void
515 const std::function<VectorType(const double, const VectorType &)> &f,
516 const double t,
517 const double factor_solution,
518 const double factor_ai,
519 const VectorType &current_ri,
520 VectorType &vec_ki,
521 VectorType &solution,
522 VectorType &next_ri) const;
523
528 };
529
530
531
536 template <typename VectorType>
537 class ImplicitRungeKutta : public RungeKutta<VectorType>
538 {
539 public:
540 using RungeKutta<VectorType>::evolve_one_time_step;
541
548
555 const unsigned int max_it = 100,
556 const double tolerance = 1e-6);
557
561 void
562 initialize(const runge_kutta_method method) override;
563
575 double
577 const std::function<VectorType(const double, const VectorType &)> &f,
578 const std::function<
579 VectorType(const double, const double, const VectorType &)>
580 &id_minus_tau_J_inverse,
581 double t,
582 double delta_t,
583 VectorType &y) override;
584
589 void
591 const double tolerance);
592
597 struct Status : public TimeStepping<VectorType>::Status
598 {
600 : method(invalid)
601 , n_iterations(numbers::invalid_unsigned_int)
602 , norm_residual(numbers::signaling_nan<double>())
603 {}
604
606 unsigned int n_iterations;
608 };
609
613 const Status &
614 get_status() const override;
615
616 private:
620 void
622 const std::function<VectorType(const double, const VectorType &)> &f,
623 const std::function<
624 VectorType(const double, const double, const VectorType &)>
625 &id_minus_tau_J_inverse,
626 double t,
627 double delta_t,
628 VectorType &y,
629 std::vector<VectorType> &f_stages);
630
634 void
636 const std::function<void(const VectorType &, VectorType &)> &get_residual,
637 const std::function<VectorType(const VectorType &)>
638 &id_minus_tau_J_inverse,
639 VectorType &y);
640
644 void
646 const std::function<VectorType(const double, const VectorType &)> &f,
647 double t,
648 double delta_t,
649 const VectorType &new_y,
650 const VectorType &y,
651 VectorType &tendency,
652 VectorType &residual) const;
653
657 unsigned int max_it;
658
662 double tolerance;
663
668 };
669
670
671
676 template <typename VectorType>
677 class EmbeddedExplicitRungeKutta : public RungeKutta<VectorType>
678 {
679 public:
680 using RungeKutta<VectorType>::evolve_one_time_step;
681
688
694 const double coarsen_param = 1.2,
695 const double refine_param = 0.8,
696 const double min_delta = 1e-14,
697 const double max_delta = 1e100,
698 const double refine_tol = 1e-8,
699 const double coarsen_tol = 1e-12);
700
705 {
706 free_memory();
707 }
708
712 void
714
718 void
719 initialize(const runge_kutta_method method) override;
720
734 double
736 const std::function<VectorType(const double, const VectorType &)> &f,
737 const std::function<
738 VectorType(const double, const double, const VectorType &)>
739 &id_minus_tau_J_inverse,
740 double t,
741 double delta_t,
742 VectorType &y) override;
743
751 double
753 const std::function<VectorType(const double, const VectorType &)> &f,
754 double t,
755 double delta_t,
756 VectorType &y);
757
761 void
763 const double refine_param,
764 const double min_delta,
765 const double max_delta,
766 const double refine_tol,
767 const double coarsen_tol);
768
775 struct Status : public TimeStepping<VectorType>::Status
776 {
778 : method(invalid)
779 {}
780
783 unsigned int n_iterations;
786 };
787
791 const Status &
792 get_status() const override;
793
794 private:
798 void
800 const std::function<VectorType(const double, const VectorType &)> &f,
801 const double t,
802 const double delta_t,
803 const VectorType &y,
804 std::vector<VectorType> &f_stages);
805
811
817
822
827
833
839
844 bool last_same_as_first = false;
845
849 std::vector<double> b1;
850
854 std::vector<double> b2;
855
860 VectorType *last_stage = nullptr;
861
866 };
867} // namespace TimeStepping
868
870
871#endif
double evolve_one_time_step(const std::function< VectorType(const double, const VectorType &)> &f, double t, double delta_t, VectorType &y)
void compute_stages(const std::function< VectorType(const double, const VectorType &)> &f, const double t, const double delta_t, const VectorType &y, std::vector< VectorType > &f_stages)
const Status & get_status() const override
void set_time_adaptation_parameters(const double coarsen_param, const double refine_param, const double min_delta, const double max_delta, const double refine_tol, const double coarsen_tol)
void initialize(const runge_kutta_method method) override
double evolve_one_time_step(const std::function< VectorType(const double, const VectorType &)> &f, const std::function< VectorType(const double, const double, const VectorType &)> &id_minus_tau_J_inverse, double t, double delta_t, VectorType &y) override
EmbeddedExplicitRungeKutta(const runge_kutta_method method, const double coarsen_param=1.2, const double refine_param=0.8, const double min_delta=1e-14, const double max_delta=1e100, const double refine_tol=1e-8, const double coarsen_tol=1e-12)
const Status & get_status() const override
void initialize(const runge_kutta_method method) override
void compute_stages(const std::function< VectorType(const double, const VectorType &)> &f, const double t, const double delta_t, const VectorType &y, std::vector< VectorType > &f_stages) const
double evolve_one_time_step(const std::function< VectorType(const double, const VectorType &)> &f, double t, double delta_t, VectorType &y)
double evolve_one_time_step(const std::function< VectorType(const double, const VectorType &)> &f, const std::function< VectorType(const double, const double, const VectorType &)> &id_minus_tau_J_inverse, double t, double delta_t, VectorType &y) override
ExplicitRungeKutta(const runge_kutta_method method)
void compute_stages(const std::function< VectorType(const double, const VectorType &)> &f, const std::function< VectorType(const double, const double, const VectorType &)> &id_minus_tau_J_inverse, double t, double delta_t, VectorType &y, std::vector< VectorType > &f_stages)
void newton_solve(const std::function< void(const VectorType &, VectorType &)> &get_residual, const std::function< VectorType(const VectorType &)> &id_minus_tau_J_inverse, VectorType &y)
const Status & get_status() const override
void initialize(const runge_kutta_method method) override
double evolve_one_time_step(const std::function< VectorType(const double, const VectorType &)> &f, const std::function< VectorType(const double, const double, const VectorType &)> &id_minus_tau_J_inverse, double t, double delta_t, VectorType &y) override
void compute_residual(const std::function< VectorType(const double, const VectorType &)> &f, double t, double delta_t, const VectorType &new_y, const VectorType &y, VectorType &tendency, VectorType &residual) const
void set_newton_solver_parameters(const unsigned int max_it, const double tolerance)
ImplicitRungeKutta(const runge_kutta_method method, const unsigned int max_it=100, const double tolerance=1e-6)
void compute_one_stage(const std::function< VectorType(const double, const VectorType &)> &f, const double t, const double factor_solution, const double factor_ai, const VectorType &current_ri, VectorType &vec_ki, VectorType &solution, VectorType &next_ri) const
void initialize(const runge_kutta_method method) override
LowStorageRungeKutta(const runge_kutta_method method)
void get_coefficients(std::vector< double > &a, std::vector< double > &b, std::vector< double > &c) const
double evolve_one_time_step(const std::function< VectorType(const double, const VectorType &)> &f, double t, double delta_t, VectorType &solution, VectorType &vec_ri, VectorType &vec_ki)
const Status & get_status() const override
double evolve_one_time_step(const std::function< VectorType(const double, const VectorType &)> &f, const std::function< VectorType(const double, const double, const VectorType &)> &id_minus_tau_J_inverse, double t, double delta_t, VectorType &y) override
std::vector< double > b
double evolve_one_time_step(std::vector< std::function< VectorType(const double, const VectorType &)> > &F, std::vector< std::function< VectorType(const double, const double, const VectorType &)> > &J_inverse, double t, double delta_t, VectorType &y) override
virtual ~RungeKutta() override=default
virtual void initialize(const runge_kutta_method method)=0
std::vector< std::vector< double > > a
std::vector< double > c
virtual double evolve_one_time_step(const std::function< VectorType(const double, const VectorType &)> &f, const std::function< VectorType(const double, const double, const VectorType &)> &id_minus_tau_J_inverse, double t, double delta_t, VectorType &y)=0
virtual double evolve_one_time_step(std::vector< std::function< VectorType(const double, const VectorType &)> > &F, std::vector< std::function< VectorType(const double, const double, const VectorType &)> > &J_inverse, double t, double delta_t, VectorType &y)=0
virtual ~TimeStepping()=default
virtual const Status & get_status() const =0
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
@ LOW_STORAGE_RK_STAGE9_ORDER5
@ LOW_STORAGE_RK_STAGE3_ORDER3
@ LOW_STORAGE_RK_STAGE7_ORDER4
@ LOW_STORAGE_RK_STAGE5_ORDER4
embedded_runge_kutta_time_step exit_delta_t