step-25.h

Go to the documentation of this file.

 = 0) const override
 *       {
 *         const double t = this->get_time();
 *   
 *         switch (dim)
 *           {
 *             case 1:
 *               {
 *                 const double m  = 0.5;
 *                 const double c1 = 0.;
 *                 const double c2 = 0.;
 *                 return -4. * std::atan(m / std::sqrt(1. - m * m) *
 *                                        std::sin(std::sqrt(1. - m * m) * t + c2) /
 *                                        std::cosh(m * p[0] + c1));
 *               }
 *   
 *             case 2:
 *               {
 *                 const double theta  = numbers::PI / 4.;
 *                 const double lambda = 1.;
 *                 const double a0     = 1.;
 *                 const double s      = 1.;
 *                 const double arg    = p[0] * std::cos(theta) +
 *                                    std::sin(theta) * (p[1] * std::cosh(lambda) +
 *                                                       t * std::sinh(lambda));
 *                 return 4. * std::atan(a0 * std::exp(s * arg));
 *               }
 *   
 *             case 3:
 *               {
 *                 const double theta = numbers::PI / 4;
 *                 const double phi   = numbers::PI / 4;
 *                 const double tau   = 1.;
 *                 const double c0    = 1.;
 *                 const double s     = 1.;
 *                 const double arg   = p[0] * std::cos(theta) +
 *                                    p[1] * std::sin(theta) * std::cos(phi) +
 *                                    std::sin(theta) * std::sin(phi) *
 *                                      (p[2] * std::cosh(tau) + t * std::sinh(tau));
 *                 return 4. * std::atan(c0 * std::exp(s * arg));
 *               }
 *   
 *             default:
 *               DEAL_II_NOT_IMPLEMENTED();
 *               return -1e8;
 *           }
 *       }
 *     };
 *   
 * @endcode
 * 
 * In the second part of this section, we provide the initial conditions. We
 * are lazy (and cautious) and don't want to implement the same functions as
 * above a second time. Rather, if we are queried for initial conditions, we
 * create an object <code>ExactSolution</code>, set it to the correct time,
 * and let it compute whatever values the exact solution has at that time:
 * 
 * @code
 *     template <int dim>
 *     class InitialValues : public Function<dim>
 *     {
 *     public:
 *       InitialValues(const unsigned int n_components = 1, const double time = 0.)
 *         : Function<dim>(n_components, time)
 *       {}
 *   
 *       virtual double value(const Point<dim>  &p,
 *                            const unsigned int component = 0) const override
 *       {
 *         return ExactSolution<dim>(1, this->get_time()).value(p, component);
 *       }
 *     };
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_25-ImplementationofthecodeSineGordonProblemcodeclass"></a> 
 * <h3>Implementation of the <code>SineGordonProblem</code> class</h3>
 * 
 
 * 
 * Let's move on to the implementation of the main class, as it implements
 * the algorithm outlined in the introduction.
 * 
 
 * 
 * 
 * <a name="step_25-SineGordonProblemSineGordonProblem"></a> 
 * <h4>SineGordonProblem::SineGordonProblem</h4>
 * 
 
 * 
 * This is the constructor of the <code>SineGordonProblem</code> class. It
 * specifies the desired polynomial degree of the finite elements,
 * associates a <code>DoFHandler</code> to the <code>triangulation</code>
 * object (just as in the example programs @ref step_3 "step-3" and @ref step_4 "step-4"), initializes
 * the current or initial time, the final time, the time step size, and the
 * value of @f$\theta@f$ for the time stepping scheme. Since the solutions we
 * compute here are time-periodic, the actual value of the start-time
 * doesn't matter, and we choose it so that we start at an interesting time.
 *   
 
 * 
 * Note that if we were to chose the explicit Euler time stepping scheme
 * (@f$\theta = 0@f$), then we must pick a time step @f$k \le h@f$, otherwise the
 * scheme is not stable and oscillations might arise in the solution. The
 * Crank-Nicolson scheme (@f$\theta = \frac{1}{2}@f$) and the implicit Euler
 * scheme (@f$\theta=1@f$) do not suffer from this deficiency, since they are
 * unconditionally stable. However, even then the time step should be chosen
 * to be on the order of @f$h@f$ in order to obtain a good solution. Since we
 * know that our mesh results from the uniform subdivision of a rectangle,
 * we can compute that time step easily; if we had a different domain, the
 * technique in @ref step_24 "step-24" using GridTools::minimal_cell_diameter would work as
 * well.
 * 
 * @code
 *     template <int dim>
 *     SineGordonProblem<dim>::SineGordonProblem()
 *       : fe(1)
 *       , dof_handler(triangulation)
 *       , n_global_refinements(6)
 *       , time(-5.4414)
 *       , final_time(2.7207)
 *       , time_step(10 * 1. / std::pow(2., 1. * n_global_refinements))
 *       , theta(0.5)
 *       , output_timestep_skip(1)
 *     {}
 *   
 * @endcode
 * 
 * 
 * <a name="step_25-SineGordonProblemmake_grid_and_dofs"></a> 
 * <h4>SineGordonProblem::make_grid_and_dofs</h4>
 * 
 
 * 
 * This function creates a rectangular grid in <code>dim</code> dimensions
 * and refines it several times. Also, all matrix and vector members of the
 * <code>SineGordonProblem</code> class are initialized to their appropriate
 * sizes once the degrees of freedom have been assembled. Like @ref step_24 "step-24", we
 * use <code>MatrixCreator</code> functions to generate a mass matrix @f$M@f$
 * and a Laplace matrix @f$A@f$ and store them in the appropriate variables for
 * the remainder of the program's life.
 * 
 * @code
 *     template <int dim>
 *     void SineGordonProblem<dim>::make_grid_and_dofs()
 *     {
 *       GridGenerator::hyper_cube(triangulation, -10, 10);
 *       triangulation.refine_global(n_global_refinements);
 *   
 *       std::cout << "   Number of active cells: " << triangulation.n_active_cells()
 *                 << std::endl
 *                 << "   Total number of cells: " << triangulation.n_cells()
 *                 << std::endl;
 *   
 *       dof_handler.distribute_dofs(fe);
 *   
 *       std::cout << "   Number of degrees of freedom: " << dof_handler.n_dofs()
 *                 << std::endl;
 *   
 *       DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
 *       DoFTools::make_sparsity_pattern(dof_handler, dsp);
 *       sparsity_pattern.copy_from(dsp);
 *   
 *       system_matrix.reinit(sparsity_pattern);
 *       mass_matrix.reinit(sparsity_pattern);
 *       laplace_matrix.reinit(sparsity_pattern);
 *   
 *       MatrixCreator::create_mass_matrix(dof_handler,
 *                                         QGauss<dim>(fe.degree + 1),
 *                                         mass_matrix);
 *       MatrixCreator::create_laplace_matrix(dof_handler,
 *                                            QGauss<dim>(fe.degree + 1),
 *                                            laplace_matrix);
 *   
 *       solution.reinit(dof_handler.n_dofs());
 *       solution_update.reinit(dof_handler.n_dofs());
 *       old_solution.reinit(dof_handler.n_dofs());
 *       M_x_velocity.reinit(dof_handler.n_dofs());
 *       system_rhs.reinit(dof_handler.n_dofs());
 *     }
 *   
 * @endcode
 * 
 * 
 * <a name="step_25-SineGordonProblemassemble_system"></a> 
 * <h4>SineGordonProblem::assemble_system</h4>
 * 
 
 * 
 * This function assembles the system matrix and right-hand side vector for
 * each iteration of Newton's method. The reader should refer to the
 * Introduction for the explicit formulas for the system matrix and
 * right-hand side.
 *   
 
 * 
 * Note that during each time step, we have to add up the various
 * contributions to the matrix and right hand sides. In contrast to @ref step_23 "step-23"
 * and @ref step_24 "step-24", this requires assembling a few more terms, since they depend
 * on the solution of the previous time step or previous nonlinear step. We
 * use the functions <code>compute_nl_matrix</code> and
 * <code>compute_nl_term</code> to do this, while the present function
 * provides the top-level logic.
 * 
 * @code
 *     template <int dim>
 *     void SineGordonProblem<dim>::assemble_system()
 *     {
 * @endcode
 * 
 * First we assemble the Jacobian matrix @f$F'_h(U^{n,l})@f$, where @f$U^{n,l}@f$
 * is stored in the vector <code>solution</code> for convenience.
 * 
 * @code
 *       system_matrix.copy_from(mass_matrix);
 *       system_matrix.add(Utilities::fixed_power<2>(time_step * theta),
 *                         laplace_matrix);
 *   
 *       SparseMatrix<double> tmp_matrix(sparsity_pattern);
 *       compute_nl_matrix(old_solution, solution, tmp_matrix);
 *       system_matrix.add(Utilities::fixed_power<2>(time_step * theta), tmp_matrix);
 *   
 * @endcode
 * 
 * Next we compute the right-hand side vector. This is just the
 * combination of matrix-vector products implied by the description of
 * @f$-F_h(U^{n,l})@f$ in the introduction.
 * 
 * @code
 *       system_rhs = 0.;
 *   
 *       Vector<double> tmp_vector(solution.size());
 *   
 *       mass_matrix.vmult(system_rhs, solution);
 *       laplace_matrix.vmult(tmp_vector, solution);
 *       system_rhs.add(Utilities::fixed_power<2>(time_step * theta), tmp_vector);
 *   
 *       mass_matrix.vmult(tmp_vector, old_solution);
 *       system_rhs.add(-1.0, tmp_vector);
 *       laplace_matrix.vmult(tmp_vector, old_solution);
 *       system_rhs.add(Utilities::fixed_power<2>(time_step) * theta * (1 - theta),
 *                      tmp_vector);
 *   
 *       system_rhs.add(-time_step, M_x_velocity);
 *   
 *       compute_nl_term(old_solution, solution, tmp_vector);
 *       system_rhs.add(Utilities::fixed_power<2>(time_step) * theta, tmp_vector);
 *   
 *       system_rhs *= -1.;
 *     }
 *   
 * @endcode
 * 
 * 
 * <a name="step_25-SineGordonProblemcompute_nl_term"></a> 
 * <h4>SineGordonProblem::compute_nl_term</h4>
 * 
 
 * 
 * This function computes the vector @f$S(\cdot,\cdot)@f$, which appears in the
 * nonlinear term in both equations of the split formulation. This
 * function not only simplifies the repeated computation of this term, but
 * it is also a fundamental part of the nonlinear iterative solver that we
 * use when the time stepping is implicit (i.e. @f$\theta\ne 0@f$). Moreover, we
 * must allow the function to receive as input an "old" and a "new"
 * solution. These may not be the actual solutions of the problem stored in
 * <code>old_solution</code> and <code>solution</code>, but are simply the
 * two functions we linearize about. For the purposes of this function, let
 * us call the first two arguments @f$w_{\mathrm{old}}@f$ and @f$w_{\mathrm{new}}@f$
 * in the documentation of this class below, respectively.
 *   
 
 * 
 * As a side-note, it is perhaps worth investigating what order quadrature
 * formula is best suited for this type of integration. Since @f$\sin(\cdot)@f$
 * is not a polynomial, there are probably no quadrature formulas that can
 * integrate these terms exactly. It is usually sufficient to just make sure
 * that the right hand side is integrated up to the same order of accuracy
 * as the discretization scheme is, but it may be possible to improve on the
 * constant in the asymptotic statement of convergence by choosing a more
 * accurate quadrature formula.
 * 
 * @code
 *     template <int dim>
 *     void SineGordonProblem<dim>::compute_nl_term(const Vector<double> &old_data,
 *                                                  const Vector<double> &new_data,
 *                                                  Vector<double> &nl_term) const
 *     {
 *       nl_term = 0;
 *       const QGauss<dim> quadrature_formula(fe.degree + 1);
 *       FEValues<dim>     fe_values(fe,
 *                               quadrature_formula,
 *                               update_values | update_JxW_values |
 *                                 update_quadrature_points);
 *   
 *       const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
 *       const unsigned int n_q_points    = quadrature_formula.size();
 *   
 *       Vector<double>                       local_nl_term(dofs_per_cell);
 *       std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
 *       std::vector<double>                  old_data_values(n_q_points);
 *       std::vector<double>                  new_data_values(n_q_points);
 *   
 *       for (const auto &cell : dof_handler.active_cell_iterators())
 *         {
 *           local_nl_term = 0;
 * @endcode
 * 
 * Once we re-initialize our <code>FEValues</code> instantiation to
 * the current cell, we make use of the
 * <code>get_function_values</code> routine to get the values of the
 * "old" data (presumably at @f$t=t_{n-1}@f$) and the "new" data
 * (presumably at @f$t=t_n@f$) at the nodes of the chosen quadrature
 * formula.
 * 
 * @code
 *           fe_values.reinit(cell);
 *           fe_values.get_function_values(old_data, old_data_values);
 *           fe_values.get_function_values(new_data, new_data_values);
 *   
 * @endcode
 * 
 * Now, we can evaluate @f$\int_K \sin\left[\theta w_{\mathrm{new}} +
 * (1-\theta) w_{\mathrm{old}}\right] \,\varphi_j\,\mathrm{d}x@f$ using
 * the desired quadrature formula.
 * 
 * @code
 *           for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
 *             for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *               local_nl_term(i) +=
 *                 (std::sin(theta * new_data_values[q_point] +
 *                           (1 - theta) * old_data_values[q_point]) *
 *                  fe_values.shape_value(i, q_point) * fe_values.JxW(q_point));
 *   
 * @endcode
 * 
 * We conclude by adding up the contributions of the integrals over
 * the cells to the global integral.
 * 
 * @code
 *           cell->get_dof_indices(local_dof_indices);
 *   
 *           for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *             nl_term(local_dof_indices[i]) += local_nl_term(i);
 *         }
 *     }
 *   
 * @endcode
 * 
 * 
 * <a name="step_25-SineGordonProblemcompute_nl_matrix"></a> 
 * <h4>SineGordonProblem::compute_nl_matrix</h4>
 * 
 
 * 
 * This is the second function dealing with the nonlinear scheme. It
 * computes the matrix @f$N(\cdot,\cdot)@f$, which appears in the nonlinear
 * term in the Jacobian of @f$F(\cdot)@f$. Just as <code>compute_nl_term</code>,
 * we must allow this function to receive as input an "old" and a "new"
 * solution, which we again call @f$w_{\mathrm{old}}@f$ and @f$w_{\mathrm{new}}@f$
 * below, respectively.
 * 
 * @code
 *     template <int dim>
 *     void SineGordonProblem<dim>::compute_nl_matrix(
 *       const Vector<double> &old_data,
 *       const Vector<double> &new_data,
 *       SparseMatrix<double> &nl_matrix) const
 *     {
 *       const QGauss<dim> quadrature_formula(fe.degree + 1);
 *       FEValues<dim>     fe_values(fe,
 *                               quadrature_formula,
 *                               update_values | update_JxW_values |
 *                                 update_quadrature_points);
 *   
 *       const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
 *       const unsigned int n_q_points    = quadrature_formula.size();
 *   
 *       FullMatrix<double> local_nl_matrix(dofs_per_cell, dofs_per_cell);
 *       std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
 *       std::vector<double>                  old_data_values(n_q_points);
 *       std::vector<double>                  new_data_values(n_q_points);
 *   
 *       for (const auto &cell : dof_handler.active_cell_iterators())
 *         {
 *           local_nl_matrix = 0;
 * @endcode
 * 
 * Again, first we re-initialize our <code>FEValues</code>
 * instantiation to the current cell.
 * 
 * @code
 *           fe_values.reinit(cell);
 *           fe_values.get_function_values(old_data, old_data_values);
 *           fe_values.get_function_values(new_data, new_data_values);
 *   
 * @endcode
 * 
 * Then, we evaluate @f$\int_K \cos\left[\theta w_{\mathrm{new}} +
 * (1-\theta) w_{\mathrm{old}}\right]\, \varphi_i\,
 * \varphi_j\,\mathrm{d}x@f$ using the desired quadrature formula.
 * 
 * @code
 *           for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
 *             for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *               for (unsigned int j = 0; j < dofs_per_cell; ++j)
 *                 local_nl_matrix(i, j) +=
 *                   (std::cos(theta * new_data_values[q_point] +
 *                             (1 - theta) * old_data_values[q_point]) *
 *                    fe_values.shape_value(i, q_point) *
 *                    fe_values.shape_value(j, q_point) * fe_values.JxW(q_point));
 *   
 * @endcode
 * 
 * Finally, we add up the contributions of the integrals over the
 * cells to the global integral.
 * 
 * @code
 *           cell->get_dof_indices(local_dof_indices);
 *   
 *           for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *             for (unsigned int j = 0; j < dofs_per_cell; ++j)
 *               nl_matrix.add(local_dof_indices[i],
 *                             local_dof_indices[j],
 *                             local_nl_matrix(i, j));
 *         }
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_25-SineGordonProblemsolve"></a> 
 * <h4>SineGordonProblem::solve</h4>
 * 
 
 * 
 * As discussed in the Introduction, this function uses the CG iterative
 * solver on the linear system of equations resulting from the finite
 * element spatial discretization of each iteration of Newton's method for
 * the (nonlinear) first equation of the split formulation. The solution to
 * the system is, in fact, @f$\delta U^{n,l}@f$ so it is stored in
 * <code>solution_update</code> and used to update <code>solution</code> in
 * the <code>run</code> function.
 *   
 
 * 
 * Note that we re-set the solution update to zero before solving for
 * it. This is not necessary: iterative solvers can start from any point and
 * converge to the correct solution. If one has a good estimate about the
 * solution of a linear system, it may be worthwhile to start from that
 * vector, but as a general observation it is a fact that the starting point
 * doesn't matter very much: it has to be a very, very good guess to reduce
 * the number of iterations by more than a few. It turns out that for this
 * problem, using the previous nonlinear update as a starting point actually
 * hurts convergence and increases the number of iterations needed, so we
 * simply set it to zero.
 *   
 
 * 
 * The function returns the number of iterations it took to converge to a
 * solution. This number will later be used to generate output on the screen
 * showing how many iterations were needed in each nonlinear iteration.
 * 
 * @code
 *     template <int dim>
 *     unsigned int SineGordonProblem<dim>::solve()
 *     {
 *       SolverControl            solver_control(1000, 1e-12 * system_rhs.l2_norm());
 *       SolverCG<Vector<double>> cg(solver_control);
 *   
 *       PreconditionSSOR<SparseMatrix<double>> preconditioner;
 *       preconditioner.initialize(system_matrix, 1.2);
 *   
 *       cg.solve(system_matrix, solution_update, system_rhs, preconditioner);
 *   
 *       return solver_control.last_step();
 *     }
 *   
 * @endcode
 * 
 * 
 * <a name="step_25-SineGordonProblemoutput_results"></a> 
 * <h4>SineGordonProblem::output_results</h4>
 * 
 
 * 
 * This function outputs the results to a file. It is pretty much identical
 * to the respective functions in @ref step_23 "step-23" and @ref step_24 "step-24":
 * 
 * @code
 *     template <int dim>
 *     void SineGordonProblem<dim>::output_results(
 *       const unsigned int timestep_number) const
 *     {
 *       DataOut<dim> data_out;
 *   
 *       data_out.attach_dof_handler(dof_handler);
 *       data_out.add_data_vector(solution, "u");
 *       data_out.build_patches();
 *   
 *       const std::string filename =
 *         "solution-" + Utilities::int_to_string(timestep_number, 3) + ".vtu";
 *       DataOutBase::VtkFlags vtk_flags;
 *       vtk_flags.compression_level = DataOutBase::CompressionLevel::best_speed;
 *       data_out.set_flags(vtk_flags);
 *       std::ofstream output(filename);
 *       data_out.write_vtu(output);
 *     }
 *   
 * @endcode
 * 
 * 
 * <a name="step_25-SineGordonProblemrun"></a> 
 * <h4>SineGordonProblem::run</h4>
 * 
 
 * 
 * This function has the top-level control over everything: it runs the
 * (outer) time-stepping loop, the (inner) nonlinear-solver loop, and
 * outputs the solution after each time step.
 * 
 * @code
 *     template <int dim>
 *     void SineGordonProblem<dim>::run()
 *     {
 *       make_grid_and_dofs();
 *   
 * @endcode
 * 
 * To acknowledge the initial condition, we must use the function @f$u_0(x)@f$
 * to compute @f$U^0@f$. To this end, below we will create an object of type
 * <code>InitialValues</code>; note that when we create this object (which
 * is derived from the <code>Function</code> class), we set its internal
 * time variable to @f$t_0@f$, to indicate that the initial condition is a
 * function of space and time evaluated at @f$t=t_0@f$.
 *     
 
 * 
 * Then we produce @f$U^0@f$ by projecting @f$u_0(x)@f$ onto the grid using
 * <code>VectorTools::project</code>. We have to use the same construct
 * using hanging node constraints as in @ref step_21 "step-21": the VectorTools::project
 * function requires a hanging node constraints object, but to be used we
 * first need to close it:
 * 
 * @code
 *       {
 *         AffineConstraints<double> constraints;
 *         constraints.close();
 *         VectorTools::project(dof_handler,
 *                              constraints,
 *                              QGauss<dim>(fe.degree + 1),
 *                              InitialValues<dim>(1, time),
 *                              solution);
 *       }
 *   
 * @endcode
 * 
 * For completeness, we output the zeroth time step to a file just like
 * any other time step.
 * 
 * @code
 *       output_results(0);
 *   
 * @endcode
 * 
 * Now we perform the time stepping: at every time step we solve the
 * matrix equation(s) corresponding to the finite element discretization
 * of the problem, and then advance our solution according to the time
 * stepping formulas we discussed in the Introduction.
 * 
 * @code
 *       unsigned int timestep_number = 1;
 *       for (time += time_step; time <= final_time;
 *            time += time_step, ++timestep_number)
 *         {
 *           old_solution = solution;
 *   
 *           std::cout << std::endl
 *                     << "Time step #" << timestep_number << "; "
 *                     << "advancing to t = " << time << '.' << std::endl;
 *   
 * @endcode
 * 
 * At the beginning of each time step we must solve the nonlinear
 * equation in the split formulation via Newton's method ---
 * i.e. solve for @f$\delta U^{n,l}@f$ then compute @f$U^{n,l+1}@f$ and so
 * on. The stopping criterion for this nonlinear iteration is that
 * @f$\|F_h(U^{n,l})\|_2 \le 10^{-6} \|F_h(U^{n,0})\|_2@f$. Consequently,
 * we need to record the norm of the residual in the first iteration.
 *         
 
 * 
 * At the end of each iteration, we output to the console how many
 * linear solver iterations it took us. When the loop below is done,
 * we have (an approximation of) @f$U^n@f$.
 * 
 * @code
 *           double initial_rhs_norm = 0.;
 *           bool   first_iteration  = true;
 *           do
 *             {
 *               assemble_system();
 *   
 *               if (first_iteration == true)
 *                 initial_rhs_norm = system_rhs.l2_norm();
 *   
 *               const unsigned int n_iterations = solve();
 *   
 *               solution += solution_update;
 *   
 *               if (first_iteration == true)
 *                 std::cout << "    " << n_iterations;
 *               else
 *                 std::cout << '+' << n_iterations;
 *               first_iteration = false;
 *             }
 *           while (system_rhs.l2_norm() > 1e-6 * initial_rhs_norm);
 *   
 *           std::cout << " CG iterations per nonlinear step." << std::endl;
 *   
 * @endcode
 * 
 * Upon obtaining the solution to the first equation of the problem at
 * @f$t=t_n@f$, we must update the auxiliary velocity variable
 * @f$V^n@f$. However, we do not compute and store @f$V^n@f$ since it is not a
 * quantity we use directly in the problem. Hence, for simplicity, we
 * update @f$MV^n@f$ directly:
 * 
 * @code
 *           Vector<double> tmp_vector(solution.size());
 *           laplace_matrix.vmult(tmp_vector, solution);
 *           M_x_velocity.add(-time_step * theta, tmp_vector);
 *   
 *           laplace_matrix.vmult(tmp_vector, old_solution);
 *           M_x_velocity.add(-time_step * (1 - theta), tmp_vector);
 *   
 *           compute_nl_term(old_solution, solution, tmp_vector);
 *           M_x_velocity.add(-time_step, tmp_vector);
 *   
 * @endcode
 * 
 * Oftentimes, in particular for fine meshes, we must pick the time
 * step to be quite small in order for the scheme to be
 * stable. Therefore, there are a lot of time steps during which
 * "nothing interesting happens" in the solution. To improve overall
 * efficiency -- in particular, speed up the program and save disk
 * space -- we only output the solution every
 * <code>output_timestep_skip</code> time steps:
 * 
 * @code
 *           if (timestep_number % output_timestep_skip == 0)
 *             output_results(timestep_number);
 *         }
 *     }
 *   } // namespace Step25
 *   
 * @endcode
 * 
 * 
 * <a name="step_25-Thecodemaincodefunction"></a> 
 * <h3>The <code>main</code> function</h3>
 * 
 
 * 
 * This is the main function of the program. It creates an object of top-level
 * class and calls its principal function. If exceptions are thrown during the
 * execution of the run method of the <code>SineGordonProblem</code> class, we
 * catch and report them here. For more information about exceptions the
 * reader should consult @ref step_6 "step-6".
 * 
 * @code
 *   int main()
 *   {
 *     try
 *       {
 *         using namespace Step25;
 *   
 *         SineGordonProblem<1> sg_problem;
 *         sg_problem.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="step_25-Results"></a><h1>Results</h1>
 
The explicit Euler time stepping scheme  (@f$\theta=0@f$) performs adequately for the problems we wish to solve. Unfortunately, a rather small time step has to be chosen due to stability issues --- @f$k\sim h/10@f$ appears to work for most the simulations we performed. On the other hand, the Crank-Nicolson scheme (@f$\theta=\frac{1}{2}@f$) is unconditionally stable, and (at least for the case of the 1D breather) we can pick the time step to be as large as @f$25h@f$ without any ill effects on the solution. The implicit Euler scheme (@f$\theta=1@f$) is "exponentially damped," so it is not a good choice for solving the sine-Gordon equation, which is conservative. However, some of the damped schemes in the continuum that is offered by the @f$\theta@f$-method were useful for eliminating spurious oscillations due to boundary effects.
 
In the simulations below, we solve the sine-Gordon equation on the interval @f$\Omega =
[-10,10]@f$ in 1D and on the square @f$\Omega = [-10,10]\times [-10,10]@f$ in 2D. In
each case, the respective grid is refined uniformly 6 times, i.e. @f$h\sim
2^{-6}@f$.
 
<a name="step_25-An11dSolution"></a><h3>An (1+1)-d Solution</h3>
 
The first example we discuss is the so-called 1D (stationary) breather
solution of the sine-Gordon equation. The breather has the following
closed-form expression, as mentioned in the Introduction:
\f[
u_{\mathrm{breather}}(x,t) = -4\arctan \left(\frac{m}{\sqrt{1-m^2}} \frac{\sin\left(\sqrt{1-m^2}t +c_2\right)}{\cosh(mx+c_1)} \right),
\f]
where @f$c_1@f$, @f$c_2@f$ and @f$m<1@f$ are constants. In the simulation below, we have chosen @f$c_1=0@f$, @f$c_2=0@f$, @f$m=0.5@f$. Moreover, it is know that the period of oscillation of the breather is @f$2\pi\sqrt{1-m^2}@f$, hence we have chosen @f$t_0=-5.4414@f$ and @f$t_f=2.7207@f$ so that we can observe three oscillations of the solution. Then, taking @f$u_0(x) = u_{\mathrm{breather}}(x,t_0)@f$, @f$\theta=0@f$ and @f$k=h/10@f$, the program computed the following solution.
 
<img src="https://www.dealii.org/images/steps/developer/step-25.1d-breather.gif" alt="Animation of the 1D stationary breather.">
 
Though not shown how to do this in the program, another way to visualize the
(1+1)-d solution is to use output generated by the DataOutStack class; it
allows to "stack" the solutions of individual time steps, so that we get
2D space-time graphs from 1D time-dependent
solutions. This produces the space-time plot below instead of the animation
above.
 
<img src="https://www.dealii.org/images/steps/developer/step-25.1d-breather_stp.png" alt="A space-time plot of the 1D stationary breather.">
 
Furthermore, since the breather is an analytical solution of the sine-Gordon
equation, we can use it to validate our code, although we have to assume that
the error introduced by our choice of Neumann boundary conditions is small
compared to the numerical error. Under this assumption, one could use the
VectorTools::integrate_difference function to compute the difference between
the numerical solution and the function described by the
<code>ExactSolution</code> class of this program. For the
simulation shown in the two images above, the @f$L^2@f$ norm of the error in the
finite element solution at each time step remained on the order of
@f$10^{-2}@f$. Hence, we can conclude that the numerical method has been
implemented correctly in the program.
 
 
<a name="step_25-Afew21DSolutions"></a><h3>A few (2+1)D Solutions</h3>
 
 
The only analytical solution to the sine-Gordon equation in (2+1)D that can be found in the literature is the so-called kink solitary wave. It has the following closed-form expression:
  @f[
    u(x,y,t) = 4 \arctan \left[a_0 e^{s\xi}\right]
  @f]
with
  @f[
    \xi = x \cos\vartheta + \sin(\vartheta) (y\cosh\lambda + t\sinh \lambda)
  @f]
where @f$a_0@f$, @f$\vartheta@f$ and @f$\lambda@f$ are constants. In the simulation below
we have chosen @f$a_0=\lambda=1@f$. Notice that if @f$\vartheta=\pi@f$ the kink is
stationary, hence it would make a good solution against which we can
validate the program in 2D because no reflections off the boundary of the
domain occur.
 
The simulation shown below was performed with @f$u_0(x) = u_{\mathrm{kink}}(x,t_0)@f$, @f$\theta=\frac{1}{2}@f$, @f$k=20h@f$, @f$t_0=1@f$ and @f$t_f=500@f$. The @f$L^2@f$ norm of the error of the finite element solution at each time step remained on the order of @f$10^{-2}@f$, showing that the program is working correctly in 2D, as well as 1D. Unfortunately, the solution is not very interesting, nonetheless we have included a snapshot of it below for completeness.
 
<img src="https://www.dealii.org/images/steps/developer/step-25.2d-kink.png" alt="Stationary 2D kink.">
 
Now that we have validated the code in 1D and 2D, we move to a problem where the analytical solution is unknown.
 
To this end, we rotate the kink solution discussed above about the @f$z@f$
axis: we let  @f$\vartheta=\frac{\pi}{4}@f$. The latter results in a
solitary wave that is not aligned with the grid, so reflections occur
at the boundaries of the domain immediately. For the simulation shown
below, we have taken @f$u_0(x)=u_{\mathrm{kink}}(x,t_0)@f$,
@f$\theta=\frac{2}{3}@f$, @f$k=20h@f$, @f$t_0=0@f$ and @f$t_f=20@f$. Moreover, we had
to pick @f$\theta=\frac{2}{3}@f$ because for any @f$\theta\le\frac{1}{2}@f$
oscillations arose at the boundary, which are likely due to the scheme
and not the equation, thus picking a value of @f$\theta@f$ a good bit into
the "exponentially damped" spectrum of the time stepping schemes
assures these oscillations are not created.
 
<img src="https://www.dealii.org/images/steps/developer/step-25.2d-angled_kink.gif" alt="Animation of a moving 2D kink, at 45 degrees to the axes of the grid, showing boundary effects.">
 
Another interesting solution to the sine-Gordon equation (which cannot be
obtained analytically) can be produced by using two 1D breathers to construct
the following separable 2D initial condition:
\f[
  u_0(x) =
  u_{\mathrm{pseudobreather}}(x,t_0) =
  16\arctan \left(
    \frac{m}{\sqrt{1-m^2}}
    \frac{\sin\left(\sqrt{1-m^2}t_0\right)}{\cosh(mx_1)} \right)
  \arctan \left(
    \frac{m}{\sqrt{1-m^2}}
    \frac{\sin\left(\sqrt{1-m^2}t_0\right)}{\cosh(mx_2)} \right),
\f]
where @f$x=(x_1,x_2)\in{R}^2@f$, @f$m=0.5<1@f$ as in the 1D case we discussed
above. For the simulation shown below, we have chosen @f$\theta=\frac{1}{2}@f$,
@f$k=10h@f$, @f$t_0=-5.4414@f$ and @f$t_f=2.7207@f$. The solution is pretty interesting
--- it acts like a breather (as far as the pictures are concerned); however,
it appears to break up and reassemble, rather than just oscillate.
 
<img src="https://www.dealii.org/images/steps/developer/step-25.2d-pseudobreather.gif" alt="Animation of a 2D pseudobreather.">
 
 
<a name="step-25-extensions"></a>
<a name="step_25-Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
 
 
It is instructive to change the initial conditions. Most choices will not lead
to solutions that stay localized (in the soliton community, such
solutions are called "stationary", though the solution does change
with time), but lead to solutions where the wave-like
character of the equation dominates and a wave travels away from the location
of a localized initial condition. For example, it is worth playing around with
the <code>InitialValues</code> class, by replacing the call to the
<code>ExactSolution</code> class by something like this function:
@f[
  u_0(x,y) = \cos\left(\frac x2\right)\cos\left(\frac y2\right)
@f]
if @f$|x|,|y|\le \frac\pi 2@f$, and @f$u_0(x,y)=0@f$ outside this region.
 
A second area would be to investigate whether the scheme is
energy-preserving. For the pure wave equation, discussed in @ref
step_23 "step-23", this is the case if we choose the time stepping
parameter such that we get the Crank-Nicolson scheme. One could do a
similar thing here, noting that the energy in the sine-Gordon solution
is defined as
@f[
  E(t) = \frac 12 \int_\Omega \left(\frac{\partial u}{\partial
  t}\right)^2
  + \left(\nabla u\right)^2 + 2 (1-\cos u) \; dx.
@f]
(We use @f$1-\cos u@f$ instead of @f$-\cos u@f$ in the formula to ensure that all
contributions to the energy are positive, and so that decaying solutions have
finite energy on unbounded domains.)
 
Beyond this, there are two obvious areas:
 
- Clearly, adaptivity (i.e. time-adaptive grids) would be of interest
  to problems like these. Their complexity leads us to leave this out
  of this program again, though the general comments in the
  introduction of @ref step_23 "step-23" remain true.
 
- Faster schemes to solve this problem. While computers today are
  plenty fast enough to solve 2d and, frequently, even 3d stationary
  problems within not too much time, time dependent problems present
  an entirely different class of problems. We address this topic in
  @ref step_48 "step-48" where we show how to solve this problem in parallel and
  without assembling or inverting any matrix at all.
 *
 *
<a name="step_25-PlainProg"></a>
<h1> The plain program</h1>
@include "step-25.cc"
*/