Reference documentation for deal.II version 9.2.0
\(\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\}}\)
step-25.h
Go to the documentation of this file.
1  = 0) const override
512  * {
513  * const double t = this->get_time();
514  *
515  * switch (dim)
516  * {
517  * case 1:
518  * {
519  * const double m = 0.5;
520  * const double c1 = 0.;
521  * const double c2 = 0.;
522  * return -4. * std::atan(m / std::sqrt(1. - m * m) *
523  * std::sin(std::sqrt(1. - m * m) * t + c2) /
524  * std::cosh(m * p[0] + c1));
525  * }
526  *
527  * case 2:
528  * {
529  * const double theta = numbers::PI / 4.;
530  * const double lambda = 1.;
531  * const double a0 = 1.;
532  * const double s = 1.;
533  * const double arg = p[0] * std::cos(theta) +
534  * std::sin(theta) * (p[1] * std::cosh(lambda) +
535  * t * std::sinh(lambda));
536  * return 4. * std::atan(a0 * std::exp(s * arg));
537  * }
538  *
539  * case 3:
540  * {
541  * const double theta = numbers::PI / 4;
542  * const double phi = numbers::PI / 4;
543  * const double tau = 1.;
544  * const double c0 = 1.;
545  * const double s = 1.;
546  * const double arg = p[0] * std::cos(theta) +
547  * p[1] * std::sin(theta) * std::cos(phi) +
548  * std::sin(theta) * std::sin(phi) *
549  * (p[2] * std::cosh(tau) + t * std::sinh(tau));
550  * return 4. * std::atan(c0 * std::exp(s * arg));
551  * }
552  *
553  * default:
554  * Assert(false, ExcNotImplemented());
555  * return -1e8;
556  * }
557  * }
558  * };
559  *
560  * @endcode
561  *
562  * In the second part of this section, we provide the initial conditions. We
563  * are lazy (and cautious) and don't want to implement the same functions as
564  * above a second time. Rather, if we are queried for initial conditions, we
565  * create an object <code>ExactSolution</code>, set it to the correct time,
566  * and let it compute whatever values the exact solution has at that time:
567  *
568  * @code
569  * template <int dim>
570  * class InitialValues : public Function<dim>
571  * {
572  * public:
573  * InitialValues(const unsigned int n_components = 1, const double time = 0.)
574  * : Function<dim>(n_components, time)
575  * {}
576  *
577  * virtual double value(const Point<dim> & p,
578  * const unsigned int component = 0) const override
579  * {
580  * return ExactSolution<dim>(1, this->get_time()).value(p, component);
581  * }
582  * };
583  *
584  *
585  * @endcode
586  *
587  *
588  * <a name="ImplementationofthecodeSineGordonProblemcodeclass"></a>
589  * <h3>Implementation of the <code>SineGordonProblem</code> class</h3>
590  *
591 
592  *
593  * Let's move on to the implementation of the main class, as it implements
594  * the algorithm outlined in the introduction.
595  *
596 
597  *
598  *
599  * <a name="SineGordonProblemSineGordonProblem"></a>
600  * <h4>SineGordonProblem::SineGordonProblem</h4>
601  *
602 
603  *
604  * This is the constructor of the <code>SineGordonProblem</code> class. It
605  * specifies the desired polynomial degree of the finite elements,
606  * associates a <code>DoFHandler</code> to the <code>triangulation</code>
607  * object (just as in the example programs @ref step_3 "step-3" and @ref step_4 "step-4"), initializes
608  * the current or initial time, the final time, the time step size, and the
609  * value of @f$\theta@f$ for the time stepping scheme. Since the solutions we
610  * compute here are time-periodic, the actual value of the start-time
611  * doesn't matter, and we choose it so that we start at an interesting time.
612  *
613 
614  *
615  * Note that if we were to chose the explicit Euler time stepping scheme
616  * (@f$\theta = 0@f$), then we must pick a time step @f$k \le h@f$, otherwise the
617  * scheme is not stable and oscillations might arise in the solution. The
618  * Crank-Nicolson scheme (@f$\theta = \frac{1}{2}@f$) and the implicit Euler
619  * scheme (@f$\theta=1@f$) do not suffer from this deficiency, since they are
620  * unconditionally stable. However, even then the time step should be chosen
621  * to be on the order of @f$h@f$ in order to obtain a good solution. Since we
622  * know that our mesh results from the uniform subdivision of a rectangle,
623  * we can compute that time step easily; if we had a different domain, the
624  * technique in @ref step_24 "step-24" using GridTools::minimal_cell_diameter would work as
625  * well.
626  *
627  * @code
628  * template <int dim>
629  * SineGordonProblem<dim>::SineGordonProblem()
630  * : fe(1)
631  * , dof_handler(triangulation)
632  * , n_global_refinements(6)
633  * , time(-5.4414)
634  * , final_time(2.7207)
635  * , time_step(10 * 1. / std::pow(2., 1. * n_global_refinements))
636  * , theta(0.5)
637  * , output_timestep_skip(1)
638  * {}
639  *
640  * @endcode
641  *
642  *
643  * <a name="SineGordonProblemmake_grid_and_dofs"></a>
644  * <h4>SineGordonProblem::make_grid_and_dofs</h4>
645  *
646 
647  *
648  * This function creates a rectangular grid in <code>dim</code> dimensions
649  * and refines it several times. Also, all matrix and vector members of the
650  * <code>SineGordonProblem</code> class are initialized to their appropriate
651  * sizes once the degrees of freedom have been assembled. Like @ref step_24 "step-24", we
652  * use <code>MatrixCreator</code> functions to generate a mass matrix @f$M@f$
653  * and a Laplace matrix @f$A@f$ and store them in the appropriate variables for
654  * the remainder of the program's life.
655  *
656  * @code
657  * template <int dim>
658  * void SineGordonProblem<dim>::make_grid_and_dofs()
659  * {
661  * triangulation.refine_global(n_global_refinements);
662  *
663  * std::cout << " Number of active cells: " << triangulation.n_active_cells()
664  * << std::endl
665  * << " Total number of cells: " << triangulation.n_cells()
666  * << std::endl;
667  *
668  * dof_handler.distribute_dofs(fe);
669  *
670  * std::cout << " Number of degrees of freedom: " << dof_handler.n_dofs()
671  * << std::endl;
672  *
673  * DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
674  * DoFTools::make_sparsity_pattern(dof_handler, dsp);
675  * sparsity_pattern.copy_from(dsp);
676  *
677  * system_matrix.reinit(sparsity_pattern);
678  * mass_matrix.reinit(sparsity_pattern);
679  * laplace_matrix.reinit(sparsity_pattern);
680  *
681  * MatrixCreator::create_mass_matrix(dof_handler,
682  * QGauss<dim>(fe.degree + 1),
683  * mass_matrix);
685  * QGauss<dim>(fe.degree + 1),
686  * laplace_matrix);
687  *
688  * solution.reinit(dof_handler.n_dofs());
689  * solution_update.reinit(dof_handler.n_dofs());
690  * old_solution.reinit(dof_handler.n_dofs());
691  * M_x_velocity.reinit(dof_handler.n_dofs());
692  * system_rhs.reinit(dof_handler.n_dofs());
693  * }
694  *
695  * @endcode
696  *
697  *
698  * <a name="SineGordonProblemassemble_system"></a>
699  * <h4>SineGordonProblem::assemble_system</h4>
700  *
701 
702  *
703  * This function assembles the system matrix and right-hand side vector for
704  * each iteration of Newton's method. The reader should refer to the
705  * Introduction for the explicit formulas for the system matrix and
706  * right-hand side.
707  *
708 
709  *
710  * Note that during each time step, we have to add up the various
711  * contributions to the matrix and right hand sides. In contrast to @ref step_23 "step-23"
712  * and @ref step_24 "step-24", this requires assembling a few more terms, since they depend
713  * on the solution of the previous time step or previous nonlinear step. We
714  * use the functions <code>compute_nl_matrix</code> and
715  * <code>compute_nl_term</code> to do this, while the present function
716  * provides the top-level logic.
717  *
718  * @code
719  * template <int dim>
720  * void SineGordonProblem<dim>::assemble_system()
721  * {
722  * @endcode
723  *
724  * First we assemble the Jacobian matrix @f$F'_h(U^{n,l})@f$, where @f$U^{n,l}@f$
725  * is stored in the vector <code>solution</code> for convenience.
726  *
727  * @code
728  * system_matrix.copy_from(mass_matrix);
729  * system_matrix.add(std::pow(time_step * theta, 2), laplace_matrix);
730  *
731  * SparseMatrix<double> tmp_matrix(sparsity_pattern);
732  * compute_nl_matrix(old_solution, solution, tmp_matrix);
733  * system_matrix.add(std::pow(time_step * theta, 2), tmp_matrix);
734  *
735  * @endcode
736  *
737  * Next we compute the right-hand side vector. This is just the
738  * combination of matrix-vector products implied by the description of
739  * @f$-F_h(U^{n,l})@f$ in the introduction.
740  *
741  * @code
742  * system_rhs = 0.;
743  *
744  * Vector<double> tmp_vector(solution.size());
745  *
746  * mass_matrix.vmult(system_rhs, solution);
747  * laplace_matrix.vmult(tmp_vector, solution);
748  * system_rhs.add(std::pow(time_step * theta, 2), tmp_vector);
749  *
750  * mass_matrix.vmult(tmp_vector, old_solution);
751  * system_rhs.add(-1.0, tmp_vector);
752  * laplace_matrix.vmult(tmp_vector, old_solution);
753  * system_rhs.add(std::pow(time_step, 2) * theta * (1 - theta), tmp_vector);
754  *
755  * system_rhs.add(-time_step, M_x_velocity);
756  *
757  * compute_nl_term(old_solution, solution, tmp_vector);
758  * system_rhs.add(std::pow(time_step, 2) * theta, tmp_vector);
759  *
760  * system_rhs *= -1.;
761  * }
762  *
763  * @endcode
764  *
765  *
766  * <a name="SineGordonProblemcompute_nl_term"></a>
767  * <h4>SineGordonProblem::compute_nl_term</h4>
768  *
769 
770  *
771  * This function computes the vector @f$S(\cdot,\cdot)@f$, which appears in the
772  * nonlinear term in both equations of the split formulation. This
773  * function not only simplifies the repeated computation of this term, but
774  * it is also a fundamental part of the nonlinear iterative solver that we
775  * use when the time stepping is implicit (i.e. @f$\theta\ne 0@f$). Moreover, we
776  * must allow the function to receive as input an "old" and a "new"
777  * solution. These may not be the actual solutions of the problem stored in
778  * <code>old_solution</code> and <code>solution</code>, but are simply the
779  * two functions we linearize about. For the purposes of this function, let
780  * us call the first two arguments @f$w_{\mathrm{old}}@f$ and @f$w_{\mathrm{new}}@f$
781  * in the documentation of this class below, respectively.
782  *
783 
784  *
785  * As a side-note, it is perhaps worth investigating what order quadrature
786  * formula is best suited for this type of integration. Since @f$\sin(\cdot)@f$
787  * is not a polynomial, there are probably no quadrature formulas that can
788  * integrate these terms exactly. It is usually sufficient to just make sure
789  * that the right hand side is integrated up to the same order of accuracy
790  * as the discretization scheme is, but it may be possible to improve on the
791  * constant in the asymptotic statement of convergence by choosing a more
792  * accurate quadrature formula.
793  *
794  * @code
795  * template <int dim>
796  * void SineGordonProblem<dim>::compute_nl_term(const Vector<double> &old_data,
797  * const Vector<double> &new_data,
798  * Vector<double> &nl_term) const
799  * {
800  * nl_term = 0;
801  * const QGauss<dim> quadrature_formula(fe.degree + 1);
802  * FEValues<dim> fe_values(fe,
803  * quadrature_formula,
806  *
807  * const unsigned int dofs_per_cell = fe.dofs_per_cell;
808  * const unsigned int n_q_points = quadrature_formula.size();
809  *
810  * Vector<double> local_nl_term(dofs_per_cell);
811  * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
812  * std::vector<double> old_data_values(n_q_points);
813  * std::vector<double> new_data_values(n_q_points);
814  *
815  * for (const auto &cell : dof_handler.active_cell_iterators())
816  * {
817  * local_nl_term = 0;
818  * @endcode
819  *
820  * Once we re-initialize our <code>FEValues</code> instantiation to
821  * the current cell, we make use of the
822  * <code>get_function_values</code> routine to get the values of the
823  * "old" data (presumably at @f$t=t_{n-1}@f$) and the "new" data
824  * (presumably at @f$t=t_n@f$) at the nodes of the chosen quadrature
825  * formula.
826  *
827  * @code
828  * fe_values.reinit(cell);
829  * fe_values.get_function_values(old_data, old_data_values);
830  * fe_values.get_function_values(new_data, new_data_values);
831  *
832  * @endcode
833  *
834  * Now, we can evaluate @f$\int_K \sin\left[\theta w_{\mathrm{new}} +
835  * (1-\theta) w_{\mathrm{old}}\right] \,\varphi_j\,\mathrm{d}x@f$ using
836  * the desired quadrature formula.
837  *
838  * @code
839  * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
840  * for (unsigned int i = 0; i < dofs_per_cell; ++i)
841  * local_nl_term(i) +=
842  * (std::sin(theta * new_data_values[q_point] +
843  * (1 - theta) * old_data_values[q_point]) *
844  * fe_values.shape_value(i, q_point) * fe_values.JxW(q_point));
845  *
846  * @endcode
847  *
848  * We conclude by adding up the contributions of the integrals over
849  * the cells to the global integral.
850  *
851  * @code
852  * cell->get_dof_indices(local_dof_indices);
853  *
854  * for (unsigned int i = 0; i < dofs_per_cell; ++i)
855  * nl_term(local_dof_indices[i]) += local_nl_term(i);
856  * }
857  * }
858  *
859  * @endcode
860  *
861  *
862  * <a name="SineGordonProblemcompute_nl_matrix"></a>
863  * <h4>SineGordonProblem::compute_nl_matrix</h4>
864  *
865 
866  *
867  * This is the second function dealing with the nonlinear scheme. It
868  * computes the matrix @f$N(\cdot,\cdot)@f$, which appears in the nonlinear
869  * term in the Jacobian of @f$F(\cdot)@f$. Just as <code>compute_nl_term</code>,
870  * we must allow this function to receive as input an "old" and a "new"
871  * solution, which we again call @f$w_{\mathrm{old}}@f$ and @f$w_{\mathrm{new}}@f$
872  * below, respectively.
873  *
874  * @code
875  * template <int dim>
876  * void SineGordonProblem<dim>::compute_nl_matrix(
877  * const Vector<double> &old_data,
878  * const Vector<double> &new_data,
879  * SparseMatrix<double> &nl_matrix) const
880  * {
881  * QGauss<dim> quadrature_formula(fe.degree + 1);
882  * FEValues<dim> fe_values(fe,
883  * quadrature_formula,
886  *
887  * const unsigned int dofs_per_cell = fe.dofs_per_cell;
888  * const unsigned int n_q_points = quadrature_formula.size();
889  *
890  * FullMatrix<double> local_nl_matrix(dofs_per_cell, dofs_per_cell);
891  * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
892  * std::vector<double> old_data_values(n_q_points);
893  * std::vector<double> new_data_values(n_q_points);
894  *
895  * for (const auto &cell : dof_handler.active_cell_iterators())
896  * {
897  * local_nl_matrix = 0;
898  * @endcode
899  *
900  * Again, first we re-initialize our <code>FEValues</code>
901  * instantiation to the current cell.
902  *
903  * @code
904  * fe_values.reinit(cell);
905  * fe_values.get_function_values(old_data, old_data_values);
906  * fe_values.get_function_values(new_data, new_data_values);
907  *
908  * @endcode
909  *
910  * Then, we evaluate @f$\int_K \cos\left[\theta w_{\mathrm{new}} +
911  * (1-\theta) w_{\mathrm{old}}\right]\, \varphi_i\,
912  * \varphi_j\,\mathrm{d}x@f$ using the desired quadrature formula.
913  *
914  * @code
915  * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
916  * for (unsigned int i = 0; i < dofs_per_cell; ++i)
917  * for (unsigned int j = 0; j < dofs_per_cell; ++j)
918  * local_nl_matrix(i, j) +=
919  * (std::cos(theta * new_data_values[q_point] +
920  * (1 - theta) * old_data_values[q_point]) *
921  * fe_values.shape_value(i, q_point) *
922  * fe_values.shape_value(j, q_point) * fe_values.JxW(q_point));
923  *
924  * @endcode
925  *
926  * Finally, we add up the contributions of the integrals over the
927  * cells to the global integral.
928  *
929  * @code
930  * cell->get_dof_indices(local_dof_indices);
931  *
932  * for (unsigned int i = 0; i < dofs_per_cell; ++i)
933  * for (unsigned int j = 0; j < dofs_per_cell; ++j)
934  * nl_matrix.add(local_dof_indices[i],
935  * local_dof_indices[j],
936  * local_nl_matrix(i, j));
937  * }
938  * }
939  *
940  *
941  *
942  * @endcode
943  *
944  *
945  * <a name="SineGordonProblemsolve"></a>
946  * <h4>SineGordonProblem::solve</h4>
947  *
948 
949  *
950  * As discussed in the Introduction, this function uses the CG iterative
951  * solver on the linear system of equations resulting from the finite
952  * element spatial discretization of each iteration of Newton's method for
953  * the (nonlinear) first equation of the split formulation. The solution to
954  * the system is, in fact, @f$\delta U^{n,l}@f$ so it is stored in
955  * <code>solution_update</code> and used to update <code>solution</code> in
956  * the <code>run</code> function.
957  *
958 
959  *
960  * Note that we re-set the solution update to zero before solving for
961  * it. This is not necessary: iterative solvers can start from any point and
962  * converge to the correct solution. If one has a good estimate about the
963  * solution of a linear system, it may be worthwhile to start from that
964  * vector, but as a general observation it is a fact that the starting point
965  * doesn't matter very much: it has to be a very, very good guess to reduce
966  * the number of iterations by more than a few. It turns out that for this
967  * problem, using the previous nonlinear update as a starting point actually
968  * hurts convergence and increases the number of iterations needed, so we
969  * simply set it to zero.
970  *
971 
972  *
973  * The function returns the number of iterations it took to converge to a
974  * solution. This number will later be used to generate output on the screen
975  * showing how many iterations were needed in each nonlinear iteration.
976  *
977  * @code
978  * template <int dim>
979  * unsigned int SineGordonProblem<dim>::solve()
980  * {
981  * SolverControl solver_control(1000, 1e-12 * system_rhs.l2_norm());
982  * SolverCG<Vector<double>> cg(solver_control);
983  *
984  * PreconditionSSOR<SparseMatrix<double>> preconditioner;
985  * preconditioner.initialize(system_matrix, 1.2);
986  *
987  * cg.solve(system_matrix, solution_update, system_rhs, preconditioner);
988  *
989  * return solver_control.last_step();
990  * }
991  *
992  * @endcode
993  *
994  *
995  * <a name="SineGordonProblemoutput_results"></a>
996  * <h4>SineGordonProblem::output_results</h4>
997  *
998 
999  *
1000  * This function outputs the results to a file. It is pretty much identical
1001  * to the respective functions in @ref step_23 "step-23" and @ref step_24 "step-24":
1002  *
1003  * @code
1004  * template <int dim>
1005  * void SineGordonProblem<dim>::output_results(
1006  * const unsigned int timestep_number) const
1007  * {
1008  * DataOut<dim> data_out;
1009  *
1010  * data_out.attach_dof_handler(dof_handler);
1011  * data_out.add_data_vector(solution, "u");
1012  * data_out.build_patches();
1013  *
1014  * const std::string filename =
1015  * "solution-" + Utilities::int_to_string(timestep_number, 3) + ".vtu";
1016  * DataOutBase::VtkFlags vtk_flags;
1017  * vtk_flags.compression_level =
1018  * DataOutBase::VtkFlags::ZlibCompressionLevel::best_speed;
1019  * data_out.set_flags(vtk_flags);
1020  * std::ofstream output(filename);
1021  * data_out.write_vtu(output);
1022  * }
1023  *
1024  * @endcode
1025  *
1026  *
1027  * <a name="SineGordonProblemrun"></a>
1028  * <h4>SineGordonProblem::run</h4>
1029  *
1030 
1031  *
1032  * This function has the top-level control over everything: it runs the
1033  * (outer) time-stepping loop, the (inner) nonlinear-solver loop, and
1034  * outputs the solution after each time step.
1035  *
1036  * @code
1037  * template <int dim>
1039  * {
1040  * make_grid_and_dofs();
1041  *
1042  * @endcode
1043  *
1044  * To acknowledge the initial condition, we must use the function @f$u_0(x)@f$
1045  * to compute @f$U^0@f$. To this end, below we will create an object of type
1046  * <code>InitialValues</code>; note that when we create this object (which
1047  * is derived from the <code>Function</code> class), we set its internal
1048  * time variable to @f$t_0@f$, to indicate that the initial condition is a
1049  * function of space and time evaluated at @f$t=t_0@f$.
1050  *
1051 
1052  *
1053  * Then we produce @f$U^0@f$ by projecting @f$u_0(x)@f$ onto the grid using
1054  * <code>VectorTools::project</code>. We have to use the same construct
1055  * using hanging node constraints as in @ref step_21 "step-21": the VectorTools::project
1056  * function requires a hanging node constraints object, but to be used we
1057  * first need to close it:
1058  *
1059  * @code
1060  * {
1061  * AffineConstraints<double> constraints;
1062  * constraints.close();
1063  * VectorTools::project(dof_handler,
1064  * constraints,
1065  * QGauss<dim>(fe.degree + 1),
1066  * InitialValues<dim>(1, time),
1067  * solution);
1068  * }
1069  *
1070  * @endcode
1071  *
1072  * For completeness, we output the zeroth time step to a file just like
1073  * any other time step.
1074  *
1075  * @code
1076  * output_results(0);
1077  *
1078  * @endcode
1079  *
1080  * Now we perform the time stepping: at every time step we solve the
1081  * matrix equation(s) corresponding to the finite element discretization
1082  * of the problem, and then advance our solution according to the time
1083  * stepping formulas we discussed in the Introduction.
1084  *
1085  * @code
1086  * unsigned int timestep_number = 1;
1087  * for (time += time_step; time <= final_time;
1088  * time += time_step, ++timestep_number)
1089  * {
1090  * old_solution = solution;
1091  *
1092  * std::cout << std::endl
1093  * << "Time step #" << timestep_number << "; "
1094  * << "advancing to t = " << time << "." << std::endl;
1095  *
1096  * @endcode
1097  *
1098  * At the beginning of each time step we must solve the nonlinear
1099  * equation in the split formulation via Newton's method ---
1100  * i.e. solve for @f$\delta U^{n,l}@f$ then compute @f$U^{n,l+1}@f$ and so
1101  * on. The stopping criterion for this nonlinear iteration is that
1102  * @f$\|F_h(U^{n,l})\|_2 \le 10^{-6} \|F_h(U^{n,0})\|_2@f$. Consequently,
1103  * we need to record the norm of the residual in the first iteration.
1104  *
1105 
1106  *
1107  * At the end of each iteration, we output to the console how many
1108  * linear solver iterations it took us. When the loop below is done,
1109  * we have (an approximation of) @f$U^n@f$.
1110  *
1111  * @code
1112  * double initial_rhs_norm = 0.;
1113  * bool first_iteration = true;
1114  * do
1115  * {
1116  * assemble_system();
1117  *
1118  * if (first_iteration == true)
1119  * initial_rhs_norm = system_rhs.l2_norm();
1120  *
1121  * const unsigned int n_iterations = solve();
1122  *
1123  * solution += solution_update;
1124  *
1125  * if (first_iteration == true)
1126  * std::cout << " " << n_iterations;
1127  * else
1128  * std::cout << '+' << n_iterations;
1129  * first_iteration = false;
1130  * }
1131  * while (system_rhs.l2_norm() > 1e-6 * initial_rhs_norm);
1132  *
1133  * std::cout << " CG iterations per nonlinear step." << std::endl;
1134  *
1135  * @endcode
1136  *
1137  * Upon obtaining the solution to the first equation of the problem at
1138  * @f$t=t_n@f$, we must update the auxiliary velocity variable
1139  * @f$V^n@f$. However, we do not compute and store @f$V^n@f$ since it is not a
1140  * quantity we use directly in the problem. Hence, for simplicity, we
1141  * update @f$MV^n@f$ directly:
1142  *
1143  * @code
1144  * Vector<double> tmp_vector(solution.size());
1145  * laplace_matrix.vmult(tmp_vector, solution);
1146  * M_x_velocity.add(-time_step * theta, tmp_vector);
1147  *
1148  * laplace_matrix.vmult(tmp_vector, old_solution);
1149  * M_x_velocity.add(-time_step * (1 - theta), tmp_vector);
1150  *
1151  * compute_nl_term(old_solution, solution, tmp_vector);
1152  * M_x_velocity.add(-time_step, tmp_vector);
1153  *
1154  * @endcode
1155  *
1156  * Oftentimes, in particular for fine meshes, we must pick the time
1157  * step to be quite small in order for the scheme to be
1158  * stable. Therefore, there are a lot of time steps during which
1159  * "nothing interesting happens" in the solution. To improve overall
1160  * efficiency -- in particular, speed up the program and save disk
1161  * space -- we only output the solution every
1162  * <code>output_timestep_skip</code> time steps:
1163  *
1164  * @code
1165  * if (timestep_number % output_timestep_skip == 0)
1166  * output_results(timestep_number);
1167  * }
1168  * }
1169  * } // namespace Step25
1170  *
1171  * @endcode
1172  *
1173  *
1174  * <a name="Thecodemaincodefunction"></a>
1175  * <h3>The <code>main</code> function</h3>
1176  *
1177 
1178  *
1179  * This is the main function of the program. It creates an object of top-level
1180  * class and calls its principal function. If exceptions are thrown during the
1181  * execution of the run method of the <code>SineGordonProblem</code> class, we
1182  * catch and report them here. For more information about exceptions the
1183  * reader should consult @ref step_6 "step-6".
1184  *
1185  * @code
1186  * int main()
1187  * {
1188  * try
1189  * {
1190  * using namespace Step25;
1191  *
1192  * SineGordonProblem<1> sg_problem;
1193  * sg_problem.run();
1194  * }
1195  * catch (std::exception &exc)
1196  * {
1197  * std::cerr << std::endl
1198  * << std::endl
1199  * << "----------------------------------------------------"
1200  * << std::endl;
1201  * std::cerr << "Exception on processing: " << std::endl
1202  * << exc.what() << std::endl
1203  * << "Aborting!" << std::endl
1204  * << "----------------------------------------------------"
1205  * << std::endl;
1206  *
1207  * return 1;
1208  * }
1209  * catch (...)
1210  * {
1211  * std::cerr << std::endl
1212  * << std::endl
1213  * << "----------------------------------------------------"
1214  * << std::endl;
1215  * std::cerr << "Unknown exception!" << std::endl
1216  * << "Aborting!" << std::endl
1217  * << "----------------------------------------------------"
1218  * << std::endl;
1219  * return 1;
1220  * }
1221  *
1222  * return 0;
1223  * }
1224  * @endcode
1225 <a name="Results"></a><h1>Results</h1>
1226 
1227 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.
1228 
1229 In the simulations below, we solve the sine-Gordon equation on the interval @f$\Omega =
1230 [-10,10]@f$ in 1D and on the square @f$\Omega = [-10,10]\times [-10,10]@f$ in 2D. In
1231 each case, the respective grid is refined uniformly 6 times, i.e. @f$h\sim
1232 2^{-6}@f$.
1233 
1234 <a name="An11dSolution"></a><h3>An (1+1)-d Solution</h3>
1235 
1236 The first example we discuss is the so-called 1D (stationary) breather
1237 solution of the sine-Gordon equation. The breather has the following
1238 closed-form expression, as mentioned in the Introduction:
1239 \f[
1240 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),
1241 \f]
1242 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.
1243 
1244 <img src="https://www.dealii.org/images/steps/developer/step-25.1d-breather.gif" alt="Animation of the 1D stationary breather.">
1245 
1246 Though not shown how to do this in the program, another way to visualize the
1247 (1+1)-d solution is to use output generated by the DataOutStack class; it
1248 allows to "stack" the solutions of individual time steps, so that we get
1249 2D space-time graphs from 1D time-dependent
1250 solutions. This produces the space-time plot below instead of the animation
1251 above.
1252 
1253 <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.">
1254 
1255 Furthermore, since the breather is an analytical solution of the sine-Gordon
1256 equation, we can use it to validate our code, although we have to assume that
1257 the error introduced by our choice of Neumann boundary conditions is small
1258 compared to the numerical error. Under this assumption, one could use the
1259 VectorTools::integrate_difference function to compute the difference between
1260 the numerical solution and the function described by the
1261 <code>ExactSolution</code> class of this program. For the
1262 simulation shown in the two images above, the @f$L^2@f$ norm of the error in the
1263 finite element solution at each time step remained on the order of
1264 @f$10^{-2}@f$. Hence, we can conclude that the numerical method has been
1265 implemented correctly in the program.
1266 
1267 
1268 <a name="Afew21DSolutions"></a><h3>A few (2+1)D Solutions</h3>
1269 
1270 
1271 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:
1272  @f[
1273  u(x,y,t) = 4 \arctan \left[a_0 e^{s\xi}\right]
1274  @f]
1275 with
1276  @f[
1277  \xi = x \cos\vartheta + \sin(\vartheta) (y\cosh\lambda + t\sinh \lambda)
1278  @f]
1279 where @f$a_0@f$, @f$\vartheta@f$ and @f$\lambda@f$ are constants. In the simulation below
1280 we have chosen @f$a_0=\lambda=1@f$. Notice that if @f$\vartheta=\pi@f$ the kink is
1281 stationary, hence it would make a good solution against which we can
1282 validate the program in 2D because no reflections off the boundary of the
1283 domain occur.
1284 
1285 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.
1286 
1287 <img src="https://www.dealii.org/images/steps/developer/step-25.2d-kink.png" alt="Stationary 2D kink.">
1288 
1289 Now that we have validated the code in 1D and 2D, we move to a problem where the analytical solution is unknown.
1290 
1291 To this end, we rotate the kink solution discussed above about the @f$z@f$
1292 axis: we let @f$\vartheta=\frac{\pi}{4}@f$. The latter results in a
1293 solitary wave that is not aligned with the grid, so reflections occur
1294 at the boundaries of the domain immediately. For the simulation shown
1295 below, we have taken @f$u_0(x)=u_{\mathrm{kink}}(x,t_0)@f$,
1296 @f$\theta=\frac{2}{3}@f$, @f$k=20h@f$, @f$t_0=0@f$ and @f$t_f=20@f$. Moreover, we had
1297 to pick @f$\theta=\frac{2}{3}@f$ because for any @f$\theta\le\frac{1}{2}@f$
1298 oscillations arose at the boundary, which are likely due to the scheme
1299 and not the equation, thus picking a value of @f$\theta@f$ a good bit into
1300 the "exponentially damped" spectrum of the time stepping schemes
1301 assures these oscillations are not created.
1302 
1303 <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.">
1304 
1305 Another interesting solution to the sine-Gordon equation (which cannot be
1306 obtained analytically) can be produced by using two 1D breathers to construct
1307 the following separable 2D initial condition:
1308 \f[
1309  u_0(x) =
1310  u_{\mathrm{pseudobreather}}(x,t_0) =
1311  16\arctan \left(
1312  \frac{m}{\sqrt{1-m^2}}
1313  \frac{\sin\left(\sqrt{1-m^2}t_0\right)}{\cosh(mx_1)} \right)
1314  \arctan \left(
1315  \frac{m}{\sqrt{1-m^2}}
1316  \frac{\sin\left(\sqrt{1-m^2}t_0\right)}{\cosh(mx_2)} \right),
1317 \f]
1318 where @f$x=(x_1,x_2)\in{R}^2@f$, @f$m=0.5<1@f$ as in the 1D case we discussed
1319 above. For the simulation shown below, we have chosen @f$\theta=\frac{1}{2}@f$,
1320 @f$k=10h@f$, @f$t_0=-5.4414@f$ and @f$t_f=2.7207@f$. The solution is pretty interesting
1321 --- it acts like a breather (as far as the pictures are concerned); however,
1322 it appears to break up and reassemble, rather than just oscillate.
1323 
1324 <img src="https://www.dealii.org/images/steps/developer/step-25.2d-pseudobreather.gif" alt="Animation of a 2D pseudobreather.">
1325 
1326 
1327 <a name="extensions"></a>
1328 <a name="Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
1329 
1330 
1331 It is instructive to change the initial conditions. Most choices will not lead
1332 to solutions that stay localized (in the soliton community, such
1333 solutions are called "stationary", though the solution does change
1334 with time), but lead to solutions where the wave-like
1335 character of the equation dominates and a wave travels away from the location
1336 of a localized initial condition. For example, it is worth playing around with
1337 the <code>InitialValues</code> class, by replacing the call to the
1338 <code>ExactSolution</code> class by something like this function:
1339 @f[
1340  u_0(x,y) = \cos\left(\frac x2\right)\cos\left(\frac y2\right)
1341 @f]
1342 if @f$|x|,|y|\le \frac\pi 2@f$, and @f$u_0(x,y)=0@f$ outside this region.
1343 
1344 A second area would be to investigate whether the scheme is
1345 energy-preserving. For the pure wave equation, discussed in @ref
1346 step_23 "step-23", this is the case if we choose the time stepping
1347 parameter such that we get the Crank-Nicolson scheme. One could do a
1348 similar thing here, noting that the energy in the sine-Gordon solution
1349 is defined as
1350 @f[
1351  E(t) = \frac 12 \int_\Omega \left(\frac{\partial u}{\partial
1352  t}\right)^2
1353  + \left(\nabla u\right)^2 + 2 (1-\cos u) \; dx.
1354 @f]
1355 (We use @f$1-\cos u@f$ instead of @f$-\cos u@f$ in the formula to ensure that all
1356 contributions to the energy are positive, and so that decaying solutions have
1357 finite energy on unbounded domains.)
1358 
1359 Beyond this, there are two obvious areas:
1360 
1361 - Clearly, adaptivity (i.e. time-adaptive grids) would be of interest
1362  to problems like these. Their complexity leads us to leave this out
1363  of this program again, though the general comments in the
1364  introduction of @ref step_23 "step-23" remain true.
1365 
1366 - Faster schemes to solve this problem. While computers today are
1367  plenty fast enough to solve 2d and, frequently, even 3d stationary
1368  problems within not too much time, time dependent problems present
1369  an entirely different class of problems. We address this topic in
1370  @ref step_48 "step-48" where we show how to solve this problem in parallel and
1371  without assembling or inverting any matrix at all.
1372  *
1373  *
1374 <a name="PlainProg"></a>
1375 <h1> The plain program</h1>
1376 @include "step-25.cc"
1377 */
Utilities::System::get_time
std::string get_time()
Definition: utilities.cc:1020
update_quadrature_points
@ update_quadrature_points
Transformed quadrature points.
Definition: fe_update_flags.h:122
LocalIntegrators::L2::mass_matrix
void mass_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
Definition: l2.h:63
DataOutInterface::write_vtu
void write_vtu(std::ostream &out) const
Definition: data_out_base.cc:6864
SolverCG
Definition: solver_cg.h:98
Utilities::int_to_string
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:474
Differentiation::SD::atan
Expression atan(const Expression &x)
Definition: symengine_math.cc:147
Differentiation::SD::OptimizerType::lambda
@ lambda
MatrixCreator::create_laplace_matrix
void create_laplace_matrix(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Quadrature< dim > &q, SparseMatrix< double > &matrix, const Function< spacedim > *const a=nullptr, const AffineConstraints< double > &constraints=AffineConstraints< double >())
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
SparseMatrix< double >
LAPACKSupport::U
static const char U
Definition: lapack_support.h:167
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
VectorTools::project
void project(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const AffineConstraints< typename VectorType::value_type > &constraints, const Quadrature< dim > &quadrature, const Function< spacedim, typename VectorType::value_type > &function, VectorType &vec, const bool enforce_zero_boundary=false, const Quadrature< dim - 1 > &q_boundary=(dim > 1 ? QGauss< dim - 1 >(2) :Quadrature< dim - 1 >(0)), const bool project_to_boundary_first=false)
DataOutBase::VtkFlags::compression_level
ZlibCompressionLevel compression_level
Definition: data_out_base.h:1159
Differentiation::SD::sinh
Expression sinh(const Expression &x)
Definition: symengine_math.cc:182
Differentiation::SD::cosh
Expression cosh(const Expression &x)
Definition: symengine_math.cc:189
update_values
@ update_values
Shape function values.
Definition: fe_update_flags.h:78
second
Point< 2 > second
Definition: grid_out.cc:4353
VectorizedArray::sin
VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &x)
Definition: vectorization.h:5297
DataOut::build_patches
virtual void build_patches(const unsigned int n_subdivisions=0)
Definition: data_out.cc:1071
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
internal::p4est::functions
int(&) functions(const void *v1, const void *v2)
Definition: p4est_wrappers.cc:339
DoFHandler
Definition: dof_handler.h:205
LinearAlgebra::CUDAWrappers::kernel::set
__global__ void set(Number *val, const Number s, const size_type N)
std::cos
inline ::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &x)
Definition: vectorization.h:5324
OpenCASCADE::point
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:188
FEValues< dim >
WorkStream::run
void run(const std::vector< std::vector< Iterator >> &colored_iterators, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
Definition: work_stream.h:1185
level
unsigned int level
Definition: grid_out.cc:4355
PreconditionSSOR::initialize
void initialize(const MatrixType &A, const typename BaseClass::AdditionalData &parameters=typename BaseClass::AdditionalData())
DoFTools::make_sparsity_pattern
void make_sparsity_pattern(const DoFHandlerType &dof_handler, SparsityPatternType &sparsity_pattern, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
Definition: dof_tools_sparsity.cc:63
SynchronousIterators::advance
void advance(std::tuple< I1, I2 > &t, const unsigned int n)
Definition: synchronous_iterator.h:146
Algorithms::Events::initial
const Event initial
Definition: event.cc:65
Physics::Elasticity::Kinematics::l
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
FEValues::reinit
void reinit(const TriaIterator< DoFCellAccessor< DoFHandlerType< dim, spacedim >, level_dof_access >> &cell)
MatrixCreator::create_mass_matrix
void create_mass_matrix(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Quadrature< dim > &q, SparseMatrix< number > &matrix, const Function< spacedim, number > *const a=nullptr, const AffineConstraints< number > &constraints=AffineConstraints< number >())
LAPACKSupport::matrix
@ matrix
Contents is actually a matrix.
Definition: lapack_support.h:60
DynamicSparsityPattern
Definition: dynamic_sparsity_pattern.h:323
Threads::internal::call
void call(const std::function< RT()> &function, internal::return_value< RT > &ret_val)
Definition: thread_management.h:607
TrilinosWrappers::internal::end
VectorType::value_type * end(VectorType &V)
Definition: trilinos_sparse_matrix.cc:65
DataOutBase::VtkFlags
Definition: data_out_base.h:1095
QGauss
Definition: quadrature_lib.h:40
DataOut_DoFData::attach_dof_handler
void attach_dof_handler(const DoFHandlerType &)
value
static const bool value
Definition: dof_tools_constraints.cc:433
AffineConstraints< double >
update_JxW_values
@ update_JxW_values
Transformed quadrature weights.
Definition: fe_update_flags.h:129
AdaptationStrategies::Refinement::split
std::vector< value_type > split(const typename ::Triangulation< dim, spacedim >::cell_iterator &parent, const value_type parent_value)
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1419
DataOutInterface::set_flags
void set_flags(const FlagType &flags)
Definition: data_out_base.cc:7837
GridGenerator::hyper_cube
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
LAPACKSupport::zero
static const types::blas_int zero
Definition: lapack_support.h:179
PreconditionSSOR
Definition: precondition.h:665
FullMatrix< double >
Function
Definition: function.h:151
triangulation
const typename ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
Definition: p4est_wrappers.cc:69
SolverControl
Definition: solver_control.h:67
MeshWorker::loop
void loop(ITERATOR begin, typename identity< ITERATOR >::type end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(DOFINFO &, DOFINFO &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, ASSEMBLER &assembler, const LoopControl &lctrl=LoopControl())
Definition: loop.h:443
first
Point< 2 > first
Definition: grid_out.cc:4352
DataOut< dim >
numbers::PI
static constexpr double PI
Definition: numbers.h:237
Vector< double >
AffineConstraints::close
void close()
std::sin
inline ::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &x)
Definition: vectorization.h:5297
DataOut_DoFData::add_data_vector
void add_data_vector(const VectorType &data, const std::vector< std::string > &names, const DataVectorType type=type_automatic, const std::vector< DataComponentInterpretation::DataComponentInterpretation > &data_component_interpretation=std::vector< DataComponentInterpretation::DataComponentInterpretation >())
Definition: data_out_dof_data.h:1090