step-52.h

Go to the documentation of this file.

,
 *                                                      const double          tau,
 *                                                      const Vector<double> &y)
 *     {
 *       SparseDirectUMFPACK inverse_mass_minus_tau_Jacobian;
 *   
 *       mass_minus_tau_Jacobian.copy_from(mass_matrix);
 *       mass_minus_tau_Jacobian.add(-tau, system_matrix);
 *   
 *       inverse_mass_minus_tau_Jacobian.initialize(mass_minus_tau_Jacobian);
 *   
 *       Vector<double> tmp(dof_handler.n_dofs());
 *       mass_matrix.vmult(tmp, y);
 *   
 *       Vector<double> result(y);
 *       inverse_mass_minus_tau_Jacobian.vmult(result, tmp);
 *   
 *       return result;
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="codeDiffusionoutput_resultscode"></a> 
 * <h4><code>Diffusion::output_results</code></h4>
 *   
 
 * 
 * The following function then outputs the solution in vtu files indexed by
 * the number of the time step and the name of the time stepping method. Of
 * course, the (exact) result should really be the same for all time
 * stepping method, but the output here at least allows us to compare them.
 * 
 * @code
 *     void Diffusion::output_results(const double                     time,
 *                                    const unsigned int               time_step,
 *                                    TimeStepping::runge_kutta_method method) const
 *     {
 *       std::string method_name;
 *   
 *       switch (method)
 *         {
 *           case TimeStepping::FORWARD_EULER:
 *             {
 *               method_name = "forward_euler";
 *               break;
 *             }
 *           case TimeStepping::RK_THIRD_ORDER:
 *             {
 *               method_name = "rk3";
 *               break;
 *             }
 *           case TimeStepping::RK_CLASSIC_FOURTH_ORDER:
 *             {
 *               method_name = "rk4";
 *               break;
 *             }
 *           case TimeStepping::BACKWARD_EULER:
 *             {
 *               method_name = "backward_euler";
 *               break;
 *             }
 *           case TimeStepping::IMPLICIT_MIDPOINT:
 *             {
 *               method_name = "implicit_midpoint";
 *               break;
 *             }
 *           case TimeStepping::SDIRK_TWO_STAGES:
 *             {
 *               method_name = "sdirk";
 *               break;
 *             }
 *           case TimeStepping::HEUN_EULER:
 *             {
 *               method_name = "heun_euler";
 *               break;
 *             }
 *           case TimeStepping::BOGACKI_SHAMPINE:
 *             {
 *               method_name = "bogacki_shampine";
 *               break;
 *             }
 *           case TimeStepping::DOPRI:
 *             {
 *               method_name = "dopri";
 *               break;
 *             }
 *           case TimeStepping::FEHLBERG:
 *             {
 *               method_name = "fehlberg";
 *               break;
 *             }
 *           case TimeStepping::CASH_KARP:
 *             {
 *               method_name = "cash_karp";
 *               break;
 *             }
 *           default:
 *             {
 *               break;
 *             }
 *         }
 *   
 *       DataOut<2> data_out;
 *   
 *       data_out.attach_dof_handler(dof_handler);
 *       data_out.add_data_vector(solution, "solution");
 *   
 *       data_out.build_patches();
 *   
 *       data_out.set_flags(DataOutBase::VtkFlags(time, time_step));
 *   
 *       const std::string filename = "solution_" + method_name + "-" +
 *                                    Utilities::int_to_string(time_step, 3) +
 *                                    ".vtu";
 *       std::ofstream output(filename);
 *       data_out.write_vtu(output);
 *   
 *       static std::vector<std::pair<double, std::string>> times_and_names;
 *   
 *       static std::string method_name_prev = "";
 *       static std::string pvd_filename;
 *       if (method_name_prev != method_name)
 *         {
 *           times_and_names.clear();
 *           method_name_prev = method_name;
 *           pvd_filename     = "solution_" + method_name + ".pvd";
 *         }
 *       times_and_names.emplace_back(time, filename);
 *       std::ofstream pvd_output(pvd_filename);
 *       DataOutBase::write_pvd_record(pvd_output, times_and_names);
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="codeDiffusionexplicit_methodcode"></a> 
 * <h4><code>Diffusion::explicit_method</code></h4>
 *   
 
 * 
 * This function is the driver for all the explicit methods. At the
 * top it initializes the time stepping and the solution (by setting
 * it to zero and then ensuring that boundary value and hanging node
 * 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 @f$M^{-1}(f(t,y))@f$, i.e, it needs
 * <code>evaluate_diffusion</code>. Because
 * <code>evaluate_diffusion</code> is a member function, it needs to
 * be bound to <code>this</code>. After each evolution step, we
 * again apply the correct boundary values and hanging node
 * constraints.
 *   
 
 * 
 * Finally, the solution is output
 * every 10 time steps.
 * 
 * @code
 *     void Diffusion::explicit_method(const TimeStepping::runge_kutta_method method,
 *                                     const unsigned int n_time_steps,
 *                                     const double       initial_time,
 *                                     const double       final_time)
 *     {
 *       const double time_step =
 *         (final_time - initial_time) / static_cast<double>(n_time_steps);
 *   
 *       solution = 0.;
 *       constraint_matrix.distribute(solution);
 *   
 *       TimeStepping::ExplicitRungeKutta<Vector<double>> explicit_runge_kutta(
 *         method);
 *       output_results(initial_time, 0, method);
 *       DiscreteTime time(initial_time, final_time, time_step);
 *       while (time.is_at_end() == false)
 *         {
 *           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.advance_time();
 *   
 *           constraint_matrix.distribute(solution);
 *   
 *           if (time.get_step_number() % 10 == 0)
 *             output_results(time.get_current_time(),
 *                            time.get_step_number(),
 *                            method);
 *         }
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="codeDiffusionimplicit_methodcode"></a> 
 * <h4><code>Diffusion::implicit_method</code></h4>
 * This function is equivalent to <code>explicit_method</code> but for
 * implicit methods. When using implicit methods, we need to evaluate
 * @f$M^{-1}(f(t,y))@f$ and @f$\left(I-\tau M^{-1} \frac{\partial f(t,y)}{\partial
 * y}\right)^{-1}@f$ for which we use the two member functions previously
 * introduced.
 * 
 * @code
 *     void Diffusion::implicit_method(const TimeStepping::runge_kutta_method method,
 *                                     const unsigned int n_time_steps,
 *                                     const double       initial_time,
 *                                     const double       final_time)
 *     {
 *       const double time_step =
 *         (final_time - initial_time) / static_cast<double>(n_time_steps);
 *   
 *       solution = 0.;
 *       constraint_matrix.distribute(solution);
 *   
 *       TimeStepping::ImplicitRungeKutta<Vector<double>> implicit_runge_kutta(
 *         method);
 *       output_results(initial_time, 0, method);
 *       DiscreteTime time(initial_time, final_time, time_step);
 *       while (time.is_at_end() == false)
 *         {
 *           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.get_current_time(),
 *             time.get_next_step_size(),
 *             solution);
 *           time.advance_time();
 *   
 *           constraint_matrix.distribute(solution);
 *   
 *           if (time.get_step_number() % 10 == 0)
 *             output_results(time.get_current_time(),
 *                            time.get_step_number(),
 *                            method);
 *         }
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="codeDiffusionembedded_explicit_methodcode"></a> 
 * <h4><code>Diffusion::embedded_explicit_method</code></h4>
 * This function is the driver for the embedded explicit methods. It requires
 * more parameters:
 * - coarsen_param: factor multiplying the current time step when the error
 * is below the threshold.
 * - refine_param: factor multiplying the current time step when the error
 * is above the threshold.
 * - min_delta: smallest time step acceptable.
 * - 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. 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.
 * 
 * @code
 *     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)
 *     {
 *       const double time_step =
 *         (final_time - initial_time) / static_cast<double>(n_time_steps);
 *       const double coarsen_param = 1.2;
 *       const double refine_param  = 0.8;
 *       const double min_delta     = 1e-8;
 *       const double max_delta     = 10 * time_step;
 *       const double refine_tol    = 1e-1;
 *       const double coarsen_tol   = 1e-5;
 *   
 *       solution = 0.;
 *       constraint_matrix.distribute(solution);
 *   
 *       TimeStepping::EmbeddedExplicitRungeKutta<Vector<double>>
 *         embedded_explicit_runge_kutta(method,
 *                                       coarsen_param,
 *                                       refine_param,
 *                                       min_delta,
 *                                       max_delta,
 *                                       refine_tol,
 *                                       coarsen_tol);
 *       output_results(initial_time, 0, method);
 *       DiscreteTime time(initial_time, final_time, time_step);
 *       while (time.is_at_end() == false)
 *         {
 *           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 (time.get_step_number() % 10 == 0)
 *             output_results(time.get_current_time(),
 *                            time.get_step_number(),
 *                            method);
 *   
 *           time.set_desired_next_step_size(
 *             embedded_explicit_runge_kutta.get_status().delta_t_guess);
 *         }
 *   
 *       return time.get_step_number();
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="codeDiffusionruncode"></a> 
 * <h4><code>Diffusion::run</code></h4>
 *   
 
 * 
 * The following is the main function of the program. At the top, we create
 * the grid (a @f$[0,5]\times [0,5]@f$ square) and refine it four times to get a
 * mesh that has 16 by 16 cells, for a total of 256.  We then set the boundary
 * indicator to 1 for those parts of the boundary where @f$x=0@f$ and @f$x=5@f$.
 * 
 * @code
 *     void Diffusion::run()
 *     {
 *       GridGenerator::hyper_cube(triangulation, 0., 5.);
 *       triangulation.refine_global(4);
 *   
 *       for (const auto &cell : triangulation.active_cell_iterators())
 *         for (const auto &face : cell->face_iterators())
 *           if (face->at_boundary())
 *             {
 *               if ((face->center()[0] == 0.) || (face->center()[0] == 5.))
 *                 face->set_boundary_id(1);
 *               else
 *                 face->set_boundary_id(0);
 *             }
 *   
 * @endcode
 * 
 * Next, we set up the linear systems and fill them with content so that
 * they can be used throughout the time stepping process:
 * 
 * @code
 *       setup_system();
 *   
 *       assemble_system();
 *   
 * @endcode
 * 
 * Finally, we solve the diffusion problem using several of the
 * Runge-Kutta methods implemented in namespace TimeStepping, each time
 * outputting the error at the end time. (As explained in the
 * introduction, since the exact solution is zero at the final time, the
 * error equals the numerical solution and can be computed by just taking
 * the @f$l_2@f$ norm of the solution vector.)
 * 
 * @code
 *       unsigned int       n_steps      = 0;
 *       const unsigned int n_time_steps = 200;
 *       const double       initial_time = 0.;
 *       const double       final_time   = 10.;
 *   
 *       std::cout << "Explicit methods:" << std::endl;
 *       explicit_method(TimeStepping::FORWARD_EULER,
 *                       n_time_steps,
 *                       initial_time,
 *                       final_time);
 *       std::cout << "   Forward Euler:            error=" << solution.l2_norm()
 *                 << std::endl;
 *   
 *       explicit_method(TimeStepping::RK_THIRD_ORDER,
 *                       n_time_steps,
 *                       initial_time,
 *                       final_time);
 *       std::cout << "   Third order Runge-Kutta:  error=" << solution.l2_norm()
 *                 << std::endl;
 *   
 *       explicit_method(TimeStepping::RK_CLASSIC_FOURTH_ORDER,
 *                       n_time_steps,
 *                       initial_time,
 *                       final_time);
 *       std::cout << "   Fourth order Runge-Kutta: error=" << solution.l2_norm()
 *                 << std::endl;
 *       std::cout << std::endl;
 *   
 *   
 *       std::cout << "Implicit methods:" << std::endl;
 *       implicit_method(TimeStepping::BACKWARD_EULER,
 *                       n_time_steps,
 *                       initial_time,
 *                       final_time);
 *       std::cout << "   Backward Euler:           error=" << solution.l2_norm()
 *                 << std::endl;
 *   
 *       implicit_method(TimeStepping::IMPLICIT_MIDPOINT,
 *                       n_time_steps,
 *                       initial_time,
 *                       final_time);
 *       std::cout << "   Implicit Midpoint:        error=" << solution.l2_norm()
 *                 << std::endl;
 *   
 *       implicit_method(TimeStepping::CRANK_NICOLSON,
 *                       n_time_steps,
 *                       initial_time,
 *                       final_time);
 *       std::cout << "   Crank-Nicolson:           error=" << solution.l2_norm()
 *                 << std::endl;
 *   
 *       implicit_method(TimeStepping::SDIRK_TWO_STAGES,
 *                       n_time_steps,
 *                       initial_time,
 *                       final_time);
 *       std::cout << "   SDIRK:                    error=" << solution.l2_norm()
 *                 << std::endl;
 *       std::cout << std::endl;
 *   
 *   
 *       std::cout << "Embedded explicit methods:" << std::endl;
 *       n_steps = embedded_explicit_method(TimeStepping::HEUN_EULER,
 *                                          n_time_steps,
 *                                          initial_time,
 *                                          final_time);
 *       std::cout << "   Heun-Euler:               error=" << solution.l2_norm()
 *                 << std::endl;
 *       std::cout << "                   steps performed=" << n_steps << std::endl;
 *   
 *       n_steps = embedded_explicit_method(TimeStepping::BOGACKI_SHAMPINE,
 *                                          n_time_steps,
 *                                          initial_time,
 *                                          final_time);
 *       std::cout << "   Bogacki-Shampine:         error=" << solution.l2_norm()
 *                 << std::endl;
 *       std::cout << "                   steps performed=" << n_steps << std::endl;
 *   
 *       n_steps = embedded_explicit_method(TimeStepping::DOPRI,
 *                                          n_time_steps,
 *                                          initial_time,
 *                                          final_time);
 *       std::cout << "   Dopri:                    error=" << solution.l2_norm()
 *                 << std::endl;
 *       std::cout << "                   steps performed=" << n_steps << std::endl;
 *   
 *       n_steps = embedded_explicit_method(TimeStepping::FEHLBERG,
 *                                          n_time_steps,
 *                                          initial_time,
 *                                          final_time);
 *       std::cout << "   Fehlberg:                 error=" << solution.l2_norm()
 *                 << std::endl;
 *       std::cout << "                   steps performed=" << n_steps << std::endl;
 *   
 *       n_steps = embedded_explicit_method(TimeStepping::CASH_KARP,
 *                                          n_time_steps,
 *                                          initial_time,
 *                                          final_time);
 *       std::cout << "   Cash-Karp:                error=" << solution.l2_norm()
 *                 << std::endl;
 *       std::cout << "                   steps performed=" << n_steps << std::endl;
 *     }
 *   } // namespace Step52
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="Thecodemaincodefunction"></a> 
 * <h3>The <code>main()</code> function</h3>
 * 
 
 * 
 * The following <code>main</code> function is similar to previous examples
 * and need not be commented on.
 * 
 * @code
 *   int main()
 *   {
 *     try
 *       {
 *         Step52::Diffusion diffusion;
 *         diffusion.run();
 *       }
 *     catch (std::exception &exc)
 *       {
 *         std::cerr << std::endl
 *                   << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         std::cerr << "Exception on processing: " << std::endl
 *                   << exc.what() << std::endl
 *                   << "Aborting!" << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         return 1;
 *       }
 *     catch (...)
 *       {
 *         std::cerr << std::endl
 *                   << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         std::cerr << "Unknown exception!" << std::endl
 *                   << "Aborting!" << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         return 1;
 *       };
 *   
 *     return 0;
 *   }
 * @endcode
<a name="Results"></a><h1>Results</h1>
 
 
The point of this program is less to show particular results, but instead to
show how it is done. This we have already demonstrated simply by discussing
the code above. Consequently, the output the program yields is relatively
sparse and consists only of the console output and the solutions given in VTU
format for visualization.
 
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
 
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.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
also see that the (rather inaccurate) Heun-Euler method increased the number of
time steps in order to satisfy the tolerance. On the other hand, the other
embedded methods used a lot less time steps than what was prescribed.
 *
 *
<a name="PlainProg"></a>
<h1> The plain program</h1>
@include "step-52.cc"
*/