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