From: Reza Rastak Date: Fri, 21 Aug 2020 17:09:58 +0000 (-0700) Subject: Use DiscreteTime in step-52 X-Git-Tag: v9.3.0-rc1~1084^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=6d9b5f1ff9ec0d5969ee175b62c21e820368d681;p=dealii.git Use DiscreteTime in step-52 --- diff --git a/examples/step-52/doc/results.dox b/examples/step-52/doc/results.dox index 8311f7fda0..03ea68c15e 100644 --- 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 7e7ef9b636..b1f0aa7577 100644 --- 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 #include #include @@ -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 // evolve_one_time_step which performs one time step. + // Time is stored and incremented through a DiscreteTime object. // // For explicit methods, evolve_one_time_step 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(n_time_steps); - double time = initial_time; solution = 0.; constraint_matrix.distribute(solution); TimeStepping::ExplicitRungeKutta> 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 &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(n_time_steps); - double time = initial_time; solution = 0.; constraint_matrix.distribute(solution); TimeStepping::ImplicitRungeKutta> 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 &y) { return this->evaluate_diffusion(time, y); }, [this](const double time, const double tau, const Vector &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. // delta_t_guess 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 delta_t_guess. + // - Automatically adjusting the step size of the last time step to ensure + // simulation ends precisely at final_time. 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(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 &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 &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