Use DiscreteTime in step-52

diff --git a/examples/step-52/doc/results.dox b/examples/step-52/doc/results.dox
index 8311f7fda07a43d28cd0f299f53a6a5b231c2dde..03ea68c15ed00a4ae8168fe3c9b5f622f0f0a821 100644 (file)
--- a/examples/step-52/doc/results.dox
+++ b/examples/step-52/doc/results.dox
@@ -10,28 +10,27 @@ The console output contains both errors and, for some of the methods, the
 number of steps they performed:
 @code
 Explicit methods:
-Forward Euler:            error=1.00883
-Third order Runge-Kutta:  error=0.000227982
-Fourth order Runge-Kutta: error=1.90541e-06
+   Forward Euler:            error=1.00883
+   Third order Runge-Kutta:  error=0.000227982
+   Fourth order Runge-Kutta: error=1.90541e-06
 
 Implicit methods:
-Backward Euler:           error=1.03428
-Implicit Midpoint:        error=0.00862702
-Crank-Nicolson:           error=0.00862675
-SDIRK:                    error=0.0042349
-
-Embedded %explicit methods:
-Heun-Euler:               error=0.0073012
-                steps performed=284
-Bogacki-Shampine:         error=0.000403281
-                steps performed=181
-Dopri:                    error=0.0165485
-                steps performed=119
-Fehlberg:                 error=0.00104926
-                steps performed=106
-Cash-Karp:                error=8.59366e-07
-                steps performed=107
+   Backward Euler:           error=1.03428
+   Implicit Midpoint:        error=0.00862702
+   Crank-Nicolson:           error=0.00862675
+   SDIRK:                    error=0.0042349
 
+Embedded explicit methods:
+   Heun-Euler:               error=0.0073012
+                   steps performed=284
+   Bogacki-Shampine:         error=0.000408407
+                   steps performed=181
+   Dopri:                    error=0.000836695
+                   steps performed=120
+   Fehlberg:                 error=0.00248922
+                   steps performed=106
+   Cash-Karp:                error=0.0787735
+                   steps performed=106
 @endcode
 
 As expected the higher order methods give (much) more accurate solutions. We

diff --git a/examples/step-52/step-52.cc b/examples/step-52/step-52.cc
index 7e7ef9b63681256b70bf87da15ef71111081b5e0..b1f0aa757710977602e4ea373cc73125d3f2c3b8 100644 (file)
--- a/examples/step-52/step-52.cc
+++ b/examples/step-52/step-52.cc
@@ -21,6 +21,7 @@
 
 // The first task as usual is to include the functionality of these well-known
 // deal.II library files and some C++ header files.
+#include <deal.II/base/discrete_time.h>
 #include <deal.II/base/function.h>
 #include <deal.II/base/quadrature_lib.h>
 
@@ -485,6 +486,7 @@ namespace Step52
   // constraints are respected; of course, with the mesh we use here,
   // hanging node constraints are not in fact an issue). It then calls
   // <code>evolve_one_time_step</code> which performs one time step.
+  // Time is stored and incremented through a DiscreteTime object.
   //
   // For explicit methods, <code>evolve_one_time_step</code> needs to
   // evaluate $M^{-1}(f(t,y))$, i.e, it needs
@@ -503,28 +505,31 @@ namespace Step52
   {
     const double time_step =
       (final_time - initial_time) / static_cast<double>(n_time_steps);
-    double time = initial_time;
 
     solution = 0.;
     constraint_matrix.distribute(solution);
 
     TimeStepping::ExplicitRungeKutta<Vector<double>> explicit_runge_kutta(
       method);
-    output_results(time, 0, method);
-    for (unsigned int i = 0; i < n_time_steps; ++i)
+    output_results(initial_time, 0, method);
+    DiscreteTime time(initial_time, final_time, time_step);
+    while (time.is_at_end() == false)
       {
-        time = explicit_runge_kutta.evolve_one_time_step(
+        explicit_runge_kutta.evolve_one_time_step(
           [this](const double time, const Vector<double> &y) {
             return this->evaluate_diffusion(time, y);
           },
-          time,
-          time_step,
+          time.get_current_time(),
+          time.get_next_step_size(),
           solution);
+        time.advance_time();
 
         constraint_matrix.distribute(solution);
 
-        if ((i + 1) % 10 == 0)
-          output_results(time, i + 1, method);
+        if (time.get_step_number() % 10 == 0)
+          output_results(time.get_current_time(),
+                         time.get_step_number(),
+                         method);
       }
   }
 
@@ -543,31 +548,34 @@ namespace Step52
   {
     const double time_step =
       (final_time - initial_time) / static_cast<double>(n_time_steps);
-    double time = initial_time;
 
     solution = 0.;
     constraint_matrix.distribute(solution);
 
     TimeStepping::ImplicitRungeKutta<Vector<double>> implicit_runge_kutta(
       method);
-    output_results(time, 0, method);
-    for (unsigned int i = 0; i < n_time_steps; ++i)
+    output_results(initial_time, 0, method);
+    DiscreteTime time(initial_time, final_time, time_step);
+    while (time.is_at_end() == false)
       {
-        time = implicit_runge_kutta.evolve_one_time_step(
+        implicit_runge_kutta.evolve_one_time_step(
           [this](const double time, const Vector<double> &y) {
             return this->evaluate_diffusion(time, y);
           },
           [this](const double time, const double tau, const Vector<double> &y) {
             return this->id_minus_tau_J_inverse(time, tau, y);
           },
-          time,
-          time_step,
+          time.get_current_time(),
+          time.get_next_step_size(),
           solution);
+        time.advance_time();
 
         constraint_matrix.distribute(solution);
 
-        if ((i + 1) % 10 == 0)
-          output_results(time, i + 1, method);
+        if (time.get_step_number() % 10 == 0)
+          output_results(time.get_current_time(),
+                         time.get_step_number(),
+                         method);
       }
   }
 
@@ -584,20 +592,27 @@ namespace Step52
   //   - max_delta: largest time step acceptable.
   //   - refine_tol: threshold above which the time step is refined.
   //   - coarsen_tol: threshold below which the time step is coarsen.
+  //
   // Embedded methods use a guessed time step. If the error using this time step
   // is too large, the time step will be reduced. If the error is below the
   // threshold, a larger time step will be tried for the next time step.
   // <code>delta_t_guess</code> is the guessed time step produced by the
-  // embedded method.
+  // embedded method. In summary, time step size is potentially modified in
+  // three ways:
+  //   - Reducing or increasing time step size within
+  //     TimeStepping::EmbeddedExplicitRungeKutta::evolve_one_time_step().
+  //   - Using the calculated <code>delta_t_guess</code>.
+  //   - Automatically adjusting the step size of the last time step to ensure
+  //     simulation ends precisely at <code>final_time</code>. This adjustment
+  //     is handled inside the DiscreteTime instance.
   unsigned int Diffusion::embedded_explicit_method(
     const TimeStepping::runge_kutta_method method,
     const unsigned int                     n_time_steps,
     const double                           initial_time,
     const double                           final_time)
   {
-    double time_step =
+    const double time_step =
       (final_time - initial_time) / static_cast<double>(n_time_steps);
-    double       time          = initial_time;
     const double coarsen_param = 1.2;
     const double refine_param  = 0.8;
     const double min_delta     = 1e-8;
@@ -616,34 +631,33 @@ namespace Step52
                                     max_delta,
                                     refine_tol,
                                     coarsen_tol);
-    output_results(time, 0, method);
-
-    // Now for the time loop. The last time step is chosen such that the final
-    // time is exactly reached.
-    unsigned int n_steps = 0;
-    while (time < final_time)
+    output_results(initial_time, 0, method);
+    DiscreteTime time(initial_time, final_time, time_step);
+    while (time.is_at_end() == false)
       {
-        if (time + time_step > final_time)
-          time_step = final_time - time;
-
-        time = embedded_explicit_runge_kutta.evolve_one_time_step(
-          [this](const double time, const Vector<double> &y) {
-            return this->evaluate_diffusion(time, y);
-          },
-          time,
-          time_step,
-          solution);
+        const double new_time =
+          embedded_explicit_runge_kutta.evolve_one_time_step(
+            [this](const double time, const Vector<double> &y) {
+              return this->evaluate_diffusion(time, y);
+            },
+            time.get_current_time(),
+            time.get_next_step_size(),
+            solution);
+        time.set_next_step_size(new_time - time.get_current_time());
+        time.advance_time();
 
         constraint_matrix.distribute(solution);
 
-        if ((n_steps + 1) % 10 == 0)
-          output_results(time, n_steps + 1, method);
+        if (time.get_step_number() % 10 == 0)
+          output_results(time.get_current_time(),
+                         time.get_step_number(),
+                         method);
 
-        time_step = embedded_explicit_runge_kutta.get_status().delta_t_guess;
-        ++n_steps;
+        time.set_desired_next_step_size(
+          embedded_explicit_runge_kutta.get_status().delta_t_guess);
       }
 
-    return n_steps;
+    return time.get_step_number();
   }
 
 
@@ -748,7 +762,7 @@ namespace Step52
                                        final_time);
     std::cout << "   Heun-Euler:               error=" << solution.l2_norm()
               << std::endl;
-    std::cout << "                steps performed=" << n_steps << std::endl;
+    std::cout << "                   steps performed=" << n_steps << std::endl;
 
     n_steps = embedded_explicit_method(TimeStepping::BOGACKI_SHAMPINE,
                                        n_time_steps,
@@ -756,7 +770,7 @@ namespace Step52
                                        final_time);
     std::cout << "   Bogacki-Shampine:         error=" << solution.l2_norm()
               << std::endl;
-    std::cout << "                steps performed=" << n_steps << std::endl;
+    std::cout << "                   steps performed=" << n_steps << std::endl;
 
     n_steps = embedded_explicit_method(TimeStepping::DOPRI,
                                        n_time_steps,
@@ -764,7 +778,7 @@ namespace Step52
                                        final_time);
     std::cout << "   Dopri:                    error=" << solution.l2_norm()
               << std::endl;
-    std::cout << "                steps performed=" << n_steps << std::endl;
+    std::cout << "                   steps performed=" << n_steps << std::endl;
 
     n_steps = embedded_explicit_method(TimeStepping::FEHLBERG,
                                        n_time_steps,
@@ -772,7 +786,7 @@ namespace Step52
                                        final_time);
     std::cout << "   Fehlberg:                 error=" << solution.l2_norm()
               << std::endl;
-    std::cout << "                steps performed=" << n_steps << std::endl;
+    std::cout << "                   steps performed=" << n_steps << std::endl;
 
     n_steps = embedded_explicit_method(TimeStepping::CASH_KARP,
                                        n_time_steps,
@@ -780,7 +794,7 @@ namespace Step52
                                        final_time);
     std::cout << "   Cash-Karp:                error=" << solution.l2_norm()
               << std::endl;
-    std::cout << "                steps performed=" << n_steps << std::endl;
+    std::cout << "                   steps performed=" << n_steps << std::endl;
   }
 } // namespace Step52