Reference documentation for deal.II version 9.6.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\}}\)
Loading...
Searching...
No Matches
step-15.h
Go to the documentation of this file.
1
548) const
549 *   {
550 *   return std::sin(2 * numbers::PI * (p[0] + p[1]));
551 *   }
552 *  
553 * @endcode
554 *
555 *
556 * <a name="step_15-ThecodeMinimalSurfaceProblemcodeclassimplementation"></a>
557 * <h3>The <code>MinimalSurfaceProblem</code> class implementation</h3>
558 *
559
560 *
561 *
562 * <a name="step_15-MinimalSurfaceProblemMinimalSurfaceProblem"></a>
563 * <h4>MinimalSurfaceProblem::MinimalSurfaceProblem</h4>
564 *
565
566 *
567 * The constructor and destructor of the class are the same as in the first
568 * few tutorials.
569 *
570
571 *
572 *
573 * @code
574 *   template <int dim>
575 *   MinimalSurfaceProblem<dim>::MinimalSurfaceProblem()
576 *   : dof_handler(triangulation)
577 *   , fe(2)
578 *   {}
579 *  
580 *  
581 * @endcode
582 *
583 *
584 * <a name="step_15-MinimalSurfaceProblemsetup_system"></a>
585 * <h4>MinimalSurfaceProblem::setup_system</h4>
586 *
587
588 *
589 * As always in the setup-system function, we set up the variables of the
590 * finite element method. There are some differences to @ref step_6 "step-6", because
591 * we need to construct two AffineConstraint<> objects.
592 *
593 * @code
594 *   template <int dim>
595 *   void MinimalSurfaceProblem<dim>::setup_system()
596 *   {
597 *   dof_handler.distribute_dofs(fe);
598 *   current_solution.reinit(dof_handler.n_dofs());
599 *  
600 *   zero_constraints.clear();
602 *   0,
604 *   zero_constraints);
605 *   DoFTools::make_hanging_node_constraints(dof_handler, zero_constraints);
606 *   zero_constraints.close();
607 *  
608 *   nonzero_constraints.clear();
610 *   0,
611 *   BoundaryValues<dim>(),
612 *   nonzero_constraints);
613 *  
614 *   DoFTools::make_hanging_node_constraints(dof_handler, nonzero_constraints);
615 *   nonzero_constraints.close();
616 *  
617 *   newton_update.reinit(dof_handler.n_dofs());
618 *   system_rhs.reinit(dof_handler.n_dofs());
619 *  
620 *   DynamicSparsityPattern dsp(dof_handler.n_dofs());
621 *   DoFTools::make_sparsity_pattern(dof_handler, dsp, zero_constraints);
622 *  
623 *   sparsity_pattern.copy_from(dsp);
624 *   system_matrix.reinit(sparsity_pattern);
625 *   }
626 *  
627 * @endcode
628 *
629 *
630 * <a name="step_15-MinimalSurfaceProblemassemble_system"></a>
631 * <h4>MinimalSurfaceProblem::assemble_system</h4>
632 *
633
634 *
635 * This function does the same as in the previous tutorials except that now,
636 * of course, the matrix and right hand side functions depend on the
637 * previous iteration's solution. As discussed in the introduction, we need
638 * to use zero boundary values for the Newton updates; this is done by using
639 * the `zero_constraints` object when assembling into the global matrix and
640 * vector.
641 *
642
643 *
644 * The top of the function contains the usual boilerplate code, setting up
645 * the objects that allow us to evaluate shape functions at quadrature
646 * points and temporary storage locations for the local matrices and
647 * vectors, as well as for the gradients of the previous solution at the
648 * quadrature points. We then start the loop over all cells:
649 *
650 * @code
651 *   template <int dim>
652 *   void MinimalSurfaceProblem<dim>::assemble_system()
653 *   {
654 *   const QGauss<dim> quadrature_formula(fe.degree + 1);
655 *  
656 *   system_matrix = 0;
657 *   system_rhs = 0;
658 *  
659 *   FEValues<dim> fe_values(fe,
660 *   quadrature_formula,
661 *   update_gradients | update_quadrature_points |
662 *   update_JxW_values);
663 *  
664 *   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
665 *   const unsigned int n_q_points = quadrature_formula.size();
666 *  
667 *   FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
668 *   Vector<double> cell_rhs(dofs_per_cell);
669 *  
670 *   std::vector<Tensor<1, dim>> old_solution_gradients(n_q_points);
671 *  
672 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
673 *  
674 *   for (const auto &cell : dof_handler.active_cell_iterators())
675 *   {
676 *   cell_matrix = 0;
677 *   cell_rhs = 0;
678 *  
679 *   fe_values.reinit(cell);
680 *  
681 * @endcode
682 *
683 * For the assembly of the linear system, we have to obtain the values
684 * of the previous solution's gradients at the quadrature
685 * points. There is a standard way of doing this: the
686 * FEValues::get_function_gradients function takes a vector that
687 * represents a finite element field defined on a DoFHandler, and
688 * evaluates the gradients of this field at the quadrature points of the
689 * cell with which the FEValues object has last been reinitialized.
690 * The values of the gradients at all quadrature points are then written
691 * into the second argument:
692 *
693 * @code
694 *   fe_values.get_function_gradients(current_solution,
695 *   old_solution_gradients);
696 *  
697 * @endcode
698 *
699 * With this, we can then do the integration loop over all quadrature
700 * points and shape functions. Having just computed the gradients of
701 * the old solution in the quadrature points, we are able to compute
702 * the coefficients @f$a_{n}@f$ in these points. The assembly of the
703 * system itself then looks similar to what we always do with the
704 * exception of the nonlinear terms, as does copying the results from
705 * the local objects into the global ones:
706 *
707 * @code
708 *   for (unsigned int q = 0; q < n_q_points; ++q)
709 *   {
710 *   const double coeff =
711 *   1.0 / std::sqrt(1 + old_solution_gradients[q] *
712 *   old_solution_gradients[q]);
713 *  
714 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
715 *   {
716 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
717 *   cell_matrix(i, j) +=
718 *   (((fe_values.shape_grad(i, q) // ((\nabla \phi_i
719 *   * coeff // * a_n
720 *   * fe_values.shape_grad(j, q)) // * \nabla \phi_j)
721 *   - // -
722 *   (fe_values.shape_grad(i, q) // (\nabla \phi_i
723 *   * coeff * coeff * coeff // * a_n^3
724 *   * (fe_values.shape_grad(j, q) // * (\nabla \phi_j
725 *   * old_solution_gradients[q]) // * \nabla u_n)
726 *   * old_solution_gradients[q])) // * \nabla u_n)))
727 *   * fe_values.JxW(q)); // * dx
728 *  
729 *   cell_rhs(i) -= (fe_values.shape_grad(i, q) // \nabla \phi_i
730 *   * coeff // * a_n
731 *   * old_solution_gradients[q] // * \nabla u_n
732 *   * fe_values.JxW(q)); // * dx
733 *   }
734 *   }
735 *  
736 *   cell->get_dof_indices(local_dof_indices);
737 *   zero_constraints.distribute_local_to_global(
738 *   cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
739 *   }
740 *   }
741 *  
742 *  
743 *  
744 * @endcode
745 *
746 *
747 * <a name="step_15-MinimalSurfaceProblemsolve"></a>
748 * <h4>MinimalSurfaceProblem::solve</h4>
749 *
750
751 *
752 * The solve function is the same as always. At the end of the solution
753 * process we update the current solution by setting
754 * @f$u^{n+1}=u^n+\alpha^n\;\delta u^n@f$.
755 *
756 * @code
757 *   template <int dim>
758 *   void MinimalSurfaceProblem<dim>::solve()
759 *   {
760 *   SolverControl solver_control(system_rhs.size(),
761 *   system_rhs.l2_norm() * 1e-6);
762 *   SolverCG<Vector<double>> solver(solver_control);
763 *  
764 *   PreconditionSSOR<SparseMatrix<double>> preconditioner;
765 *   preconditioner.initialize(system_matrix, 1.2);
766 *  
767 *   solver.solve(system_matrix, newton_update, system_rhs, preconditioner);
768 *  
769 *   zero_constraints.distribute(newton_update);
770 *  
771 *   const double alpha = determine_step_length();
772 *   current_solution.add(alpha, newton_update);
773 *   }
774 *  
775 *  
776 * @endcode
777 *
778 *
779 * <a name="step_15-MinimalSurfaceProblemrefine_mesh"></a>
780 * <h4>MinimalSurfaceProblem::refine_mesh</h4>
781 *
782
783 *
784 * The first part of this function is the same as in @ref step_6 "step-6"... However,
785 * after refining the mesh we have to transfer the old solution to the new
786 * one which we do with the help of the SolutionTransfer class. The process
787 * is slightly convoluted, so let us describe it in detail:
788 *
789 * @code
790 *   template <int dim>
791 *   void MinimalSurfaceProblem<dim>::refine_mesh()
792 *   {
793 *   Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
794 *  
796 *   dof_handler,
797 *   QGauss<dim - 1>(fe.degree + 1),
798 *   std::map<types::boundary_id, const Function<dim> *>(),
799 *   current_solution,
800 *   estimated_error_per_cell);
801 *  
803 *   estimated_error_per_cell,
804 *   0.3,
805 *   0.03);
806 *  
807 * @endcode
808 *
809 * Then we need an additional step: if, for example, you flag a cell that
810 * is once more refined than its neighbor, and that neighbor is not
811 * flagged for refinement, we would end up with a jump of two refinement
812 * levels across a cell interface. To avoid these situations, the library
813 * will silently also have to refine the neighbor cell once. It does so by
815 * before actually doing the refinement and coarsening. This function
816 * flags a set of additional cells for refinement or coarsening, to
817 * enforce rules like the one-hanging-node rule. The cells that are
818 * flagged for refinement and coarsening after calling this function are
819 * exactly the ones that will actually be refined or coarsened. Usually,
820 * you don't have to do this by hand
821 * (Triangulation::execute_coarsening_and_refinement does this for
822 * you). However, we need to initialize the SolutionTransfer class and it
823 * needs to know the final set of cells that will be coarsened or refined
824 * in order to store the data from the old mesh and transfer to the new
825 * one. Thus, we call the function by hand:
826 *
827 * @code
828 *   triangulation.prepare_coarsening_and_refinement();
829 *  
830 * @endcode
831 *
832 * With this out of the way, we initialize a SolutionTransfer object with
833 * the present DoFHandler. We make a copy of the solution vector and attach
834 * it to the SolutionTransfer. Now we can actually execute the refinement
835 * and create the new matrices and vectors including the vector
836 * `current_solution`, that will hold the current solution on the new mesh
837 * after calling SolutionTransfer::interpolate():
838 *
839 * @code
840 *   SolutionTransfer<dim> solution_transfer(dof_handler);
841 *   const Vector<double> coarse_solution = current_solution;
842 *   solution_transfer.prepare_for_coarsening_and_refinement(coarse_solution);
843 *  
844 *   triangulation.execute_coarsening_and_refinement();
845 *  
846 *   setup_system();
847 *  
848 *   solution_transfer.interpolate(coarse_solution, current_solution);
849 *  
850 * @endcode
851 *
852 * On the new mesh, there are different hanging nodes, computed in
853 * `setup_system()` above. To be on the safe side, we should make sure that
854 * the current solution's vector entries satisfy the hanging node
855 * constraints (see the discussion in the documentation of the
856 * SolutionTransfer class for why this is necessary) and boundary values. As
857 * explained at the end of the introduction, the interpolated solution does
858 * not automatically satisfy the boundary values even if the solution before
859 * refinement had the correct boundary values.
860 *
861 * @code
862 *   nonzero_constraints.distribute(current_solution);
863 *   }
864 *  
865 *  
866 *  
867 * @endcode
868 *
869 *
870 * <a name="step_15-MinimalSurfaceProblemcompute_residual"></a>
871 * <h4>MinimalSurfaceProblem::compute_residual</h4>
872 *
873
874 *
875 * In order to monitor convergence, we need a way to compute the norm of the
876 * (discrete) residual, i.e., the norm of the vector
877 * @f$\left<F(u^n),\varphi_i\right>@f$ with @f$F(u)=-\nabla \cdot \left(
878 * \frac{1}{\sqrt{1+|\nabla u|^{2}}}\nabla u \right)@f$ as discussed in the
879 * introduction. It turns out that (although we don't use this feature in
880 * the current version of the program) one needs to compute the residual
881 * @f$\left<F(u^n+\alpha^n\;\delta u^n),\varphi_i\right>@f$ when determining
882 * optimal step lengths, and so this is what we implement here: the function
883 * takes the step length @f$\alpha^n@f$ as an argument. The original
884 * functionality is of course obtained by passing a zero as argument.
885 *
886
887 *
888 * In the function below, we first set up a vector for the residual, and
889 * then a vector for the evaluation point @f$u^n+\alpha^n\;\delta u^n@f$. This
890 * is followed by the same boilerplate code we use for all integration
891 * operations:
892 *
893 * @code
894 *   template <int dim>
895 *   double MinimalSurfaceProblem<dim>::compute_residual(const double alpha) const
896 *   {
897 *   Vector<double> residual(dof_handler.n_dofs());
898 *  
899 *   Vector<double> evaluation_point(dof_handler.n_dofs());
900 *   evaluation_point = current_solution;
901 *   evaluation_point.add(alpha, newton_update);
902 *  
903 *   const QGauss<dim> quadrature_formula(fe.degree + 1);
904 *   FEValues<dim> fe_values(fe,
905 *   quadrature_formula,
906 *   update_gradients | update_quadrature_points |
907 *   update_JxW_values);
908 *  
909 *   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
910 *   const unsigned int n_q_points = quadrature_formula.size();
911 *  
912 *   Vector<double> cell_residual(dofs_per_cell);
913 *   std::vector<Tensor<1, dim>> gradients(n_q_points);
914 *  
915 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
916 *  
917 *   for (const auto &cell : dof_handler.active_cell_iterators())
918 *   {
919 *   cell_residual = 0;
920 *   fe_values.reinit(cell);
921 *  
922 * @endcode
923 *
924 * The actual computation is much as in
925 * <code>assemble_system()</code>. We first evaluate the gradients of
926 * @f$u^n+\alpha^n\,\delta u^n@f$ at the quadrature points, then compute
927 * the coefficient @f$a_n@f$, and then plug it all into the formula for
928 * the residual:
929 *
930 * @code
931 *   fe_values.get_function_gradients(evaluation_point, gradients);
932 *  
933 *  
934 *   for (unsigned int q = 0; q < n_q_points; ++q)
935 *   {
936 *   const double coeff =
937 *   1. / std::sqrt(1 + gradients[q] * gradients[q]);
938 *  
939 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
940 *   cell_residual(i) -= (fe_values.shape_grad(i, q) // \nabla \phi_i
941 *   * coeff // * a_n
942 *   * gradients[q] // * \nabla u_n
943 *   * fe_values.JxW(q)); // * dx
944 *   }
945 *  
946 *   cell->get_dof_indices(local_dof_indices);
947 *   zero_constraints.distribute_local_to_global(cell_residual,
948 *   local_dof_indices,
949 *   residual);
950 *   }
951 *  
952 *   return residual.l2_norm();
953 *   }
954 *  
955 *  
956 *  
957 * @endcode
958 *
959 *
960 * <a name="step_15-MinimalSurfaceProblemdetermine_step_length"></a>
961 * <h4>MinimalSurfaceProblem::determine_step_length</h4>
962 *
963
964 *
965 * As discussed in the introduction, Newton's method frequently does not
966 * converge if we always take full steps, i.e., compute @f$u^{n+1}=u^n+\delta
967 * u^n@f$. Rather, one needs a damping parameter (step length) @f$\alpha^n@f$ and
968 * set @f$u^{n+1}=u^n+\alpha^n\delta u^n@f$. This function is the one called
969 * to compute @f$\alpha^n@f$.
970 *
971
972 *
973 * Here, we simply always return 0.1. This is of course a sub-optimal
974 * choice: ideally, what one wants is that the step size goes to one as we
975 * get closer to the solution, so that we get to enjoy the rapid quadratic
976 * convergence of Newton's method. We will discuss better strategies below
977 * in the results section, and @ref step_77 "step-77" also covers this aspect.
978 *
979 * @code
980 *   template <int dim>
981 *   double MinimalSurfaceProblem<dim>::determine_step_length() const
982 *   {
983 *   return 0.1;
984 *   }
985 *  
986 *  
987 *  
988 * @endcode
989 *
990 *
991 * <a name="step_15-MinimalSurfaceProblemoutput_results"></a>
992 * <h4>MinimalSurfaceProblem::output_results</h4>
993 *
994
995 *
996 * This last function to be called from `run()` outputs the current solution
997 * (and the Newton update) in graphical form as a VTU file. It is entirely the
998 * same as what has been used in previous tutorials.
999 *
1000 * @code
1001 *   template <int dim>
1002 *   void MinimalSurfaceProblem<dim>::output_results(
1003 *   const unsigned int refinement_cycle) const
1004 *   {
1005 *   DataOut<dim> data_out;
1006 *  
1007 *   data_out.attach_dof_handler(dof_handler);
1008 *   data_out.add_data_vector(current_solution, "solution");
1009 *   data_out.add_data_vector(newton_update, "update");
1010 *   data_out.build_patches();
1011 *  
1012 *   const std::string filename =
1013 *   "solution-" + Utilities::int_to_string(refinement_cycle, 2) + ".vtu";
1014 *   std::ofstream output(filename);
1015 *   data_out.write_vtu(output);
1016 *   }
1017 *  
1018 *  
1019 * @endcode
1020 *
1021 *
1022 * <a name="step_15-MinimalSurfaceProblemrun"></a>
1023 * <h4>MinimalSurfaceProblem::run</h4>
1024 *
1025
1026 *
1027 * In the run function, we build the first grid and then have the top-level
1028 * logic for the Newton iteration.
1029 *
1030
1031 *
1032 * As described in the introduction, the domain is the unit disk around
1033 * the origin, created in the same way as shown in @ref step_6 "step-6". The mesh is
1034 * globally refined twice followed later on by several adaptive cycles.
1035 *
1036
1037 *
1038 * Before starting the Newton loop, we also need to do
1039 * ensure that the first Newton iterate already has the correct
1040 * boundary values, as discussed in the introduction.
1041 *
1042 * @code
1043 *   template <int dim>
1044 *   void MinimalSurfaceProblem<dim>::run()
1045 *   {
1046 *   GridGenerator::hyper_ball(triangulation);
1047 *   triangulation.refine_global(2);
1048 *  
1049 *   setup_system();
1050 *   nonzero_constraints.distribute(current_solution);
1051 *  
1052 * @endcode
1053 *
1054 * The Newton iteration starts next. We iterate until the (norm of the)
1055 * residual computed at the end of the previous iteration is less than
1056 * @f$10^{-3}@f$, as checked at the end of the `do { ... } while` loop that
1057 * starts here. Because we don't have a reasonable value to initialize
1058 * the variable, we just use the largest value that can be represented
1059 * as a `double`.
1060 *
1061 * @code
1062 *   double last_residual_norm = std::numeric_limits<double>::max();
1063 *   unsigned int refinement_cycle = 0;
1064 *   do
1065 *   {
1066 *   std::cout << "Mesh refinement step " << refinement_cycle << std::endl;
1067 *  
1068 *   if (refinement_cycle != 0)
1069 *   refine_mesh();
1070 *  
1071 * @endcode
1072 *
1073 * On every mesh we do exactly five Newton steps. We print the initial
1074 * residual here and then start the iterations on this mesh.
1075 *
1076
1077 *
1078 * In every Newton step the system matrix and the right hand side have
1079 * to be computed first, after which we store the norm of the right
1080 * hand side as the residual to check against when deciding whether to
1081 * stop the iterations. We then solve the linear system (the function
1082 * also updates @f$u^{n+1}=u^n+\alpha^n\;\delta u^n@f$) and output the
1083 * norm of the residual at the end of this Newton step.
1084 *
1085
1086 *
1087 * After the end of this loop, we then also output the solution on the
1088 * current mesh in graphical form and increment the counter for the
1089 * mesh refinement cycle.
1090 *
1091 * @code
1092 *   std::cout << " Initial residual: " << compute_residual(0) << std::endl;
1093 *  
1094 *   for (unsigned int inner_iteration = 0; inner_iteration < 5;
1095 *   ++inner_iteration)
1096 *   {
1097 *   assemble_system();
1098 *   last_residual_norm = system_rhs.l2_norm();
1099 *  
1100 *   solve();
1101 *  
1102 *   std::cout << " Residual: " << compute_residual(0) << std::endl;
1103 *   }
1104 *  
1105 *   output_results(refinement_cycle);
1106 *  
1107 *   ++refinement_cycle;
1108 *   std::cout << std::endl;
1109 *   }
1110 *   while (last_residual_norm > 1e-2);
1111 *   }
1112 *   } // namespace Step15
1113 *  
1114 * @endcode
1115 *
1116 *
1117 * <a name="step_15-Themainfunction"></a>
1118 * <h4>The main function</h4>
1119 *
1120
1121 *
1122 * Finally the main function. This follows the scheme of all other main
1123 * functions:
1124 *
1125 * @code
1126 *   int main()
1127 *   {
1128 *   try
1129 *   {
1130 *   using namespace Step15;
1131 *  
1132 *   MinimalSurfaceProblem<2> problem;
1133 *   problem.run();
1134 *   }
1135 *   catch (std::exception &exc)
1136 *   {
1137 *   std::cerr << std::endl
1138 *   << std::endl
1139 *   << "----------------------------------------------------"
1140 *   << std::endl;
1141 *   std::cerr << "Exception on processing: " << std::endl
1142 *   << exc.what() << std::endl
1143 *   << "Aborting!" << std::endl
1144 *   << "----------------------------------------------------"
1145 *   << std::endl;
1146 *  
1147 *   return 1;
1148 *   }
1149 *   catch (...)
1150 *   {
1151 *   std::cerr << std::endl
1152 *   << std::endl
1153 *   << "----------------------------------------------------"
1154 *   << std::endl;
1155 *   std::cerr << "Unknown exception!" << std::endl
1156 *   << "Aborting!" << std::endl
1157 *   << "----------------------------------------------------"
1158 *   << std::endl;
1159 *   return 1;
1160 *   }
1161 *   return 0;
1162 *   }
1163 * @endcode
1164<a name="step_15-Results"></a><h1>Results</h1>
1165
1166
1167
1168The output of the program looks as follows:
1169@code
1170Mesh refinement step 0
1171 Initial residual: 1.53143
1172 Residual: 1.08746
1173 Residual: 0.966748
1174 Residual: 0.859602
1175 Residual: 0.766462
1176 Residual: 0.685475
1177
1178Mesh refinement step 1
1179 Initial residual: 0.868959
1180 Residual: 0.762125
1181 Residual: 0.677792
1182 Residual: 0.605762
1183 Residual: 0.542748
1184 Residual: 0.48704
1185
1186Mesh refinement step 2
1187 Initial residual: 0.426445
1188 Residual: 0.382731
1189 Residual: 0.343865
1190 Residual: 0.30918
1191 Residual: 0.278147
1192 Residual: 0.250327
1193
1194Mesh refinement step 3
1195 Initial residual: 0.282026
1196 Residual: 0.253146
1197 Residual: 0.227414
1198 Residual: 0.20441
1199 Residual: 0.183803
1200 Residual: 0.165319
1201
1202Mesh refinement step 4
1203 Initial residual: 0.154404
1204 Residual: 0.138723
1205 Residual: 0.124694
1206 Residual: 0.112124
1207 Residual: 0.100847
1208 Residual: 0.0907222
1209
1210....
1211@endcode
1212
1213Obviously, the scheme converges, if not very fast. We will come back to
1214strategies for accelerating the method below.
1215
1216One can visualize the solution after each set of five Newton
1217iterations, i.e., on each of the meshes on which we approximate the
1218solution. This yields the following set of images:
1219
1220<div class="twocolumn" style="width: 80%">
1221 <div>
1222 <img src="https://www.dealii.org/images/steps/developer/step_15_solution_1.png"
1223 alt="Solution after zero cycles with contour lines." width="230" height="273">
1224 </div>
1225 <div>
1226 <img src="https://www.dealii.org/images/steps/developer/step_15_solution_2.png"
1227 alt="Solution after one cycle with contour lines." width="230" height="273">
1228 </div>
1229 <div>
1230 <img src="https://www.dealii.org/images/steps/developer/step_15_solution_3.png"
1231 alt="Solution after two cycles with contour lines." width="230" height="273">
1232 </div>
1233 <div>
1234 <img src="https://www.dealii.org/images/steps/developer/step_15_solution_4.png"
1235 alt="Solution after three cycles with contour lines." width="230" height="273">
1236 </div>
1237</div>
1238
1239It is clearly visible, that the solution minimizes the surface
1240after each refinement. The solution converges to a picture one
1241would imagine a soap bubble to be that is located inside a wire loop
1242that is bent like
1243the boundary. Also it is visible, how the boundary
1244is smoothed out after each refinement. On the coarse mesh,
1245the boundary doesn't look like a sine, whereas it does the
1246finer the mesh gets.
1247
1248The mesh is mostly refined near the boundary, where the solution
1249increases or decreases strongly, whereas it is coarsened on
1250the inside of the domain, where nothing interesting happens,
1251because there isn't much change in the solution. The ninth
1252solution and mesh are shown here:
1253
1254<div class="onecolumn" style="width: 60%">
1255 <div>
1256 <img src="https://www.dealii.org/images/steps/developer/step_15_solution_9.png"
1257 alt="Grid and solution of the ninth cycle with contour lines." width="507" height="507">
1258 </div>
1259</div>
1260
1261
1262
1263<a name="step-15-extensions"></a>
1264<a name="step_15-Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
1265
1266
1267The program shows the basic structure of a solver for a nonlinear, stationary
1268problem. However, it does not converge particularly fast, for good reasons:
1269
1270- The program always takes a step size of 0.1. This precludes the rapid,
1271 quadratic convergence for which Newton's method is typically chosen.
1272- It does not connect the nonlinear iteration with the mesh refinement
1273 iteration.
1274
1275Obviously, a better program would have to address these two points.
1276We will discuss them in the following.
1277
1278
1279<a name="step_15-Steplengthcontrol"></a><h4> Step length control </h4>
1280
1281
1282Newton's method has two well known properties:
1283- It may not converge from arbitrarily chosen starting points. Rather, a
1284 starting point has to be close enough to the solution to guarantee
1285 convergence. However, we can enlarge the area from which Newton's method
1286 converges by damping the iteration using a <i>step length</i> 0<@f$\alpha^n\le
1287 1@f$.
1288- It exhibits rapid convergence of quadratic order if (i) the step length is
1289 chosen as @f$\alpha^n=1@f$, and (ii) it does in fact converge with this choice
1290 of step length.
1291
1292A consequence of these two observations is that a successful strategy is to
1293choose @f$\alpha^n<1@f$ for the initial iterations until the iterate has come
1294close enough to allow for convergence with full step length, at which point we
1295want to switch to @f$\alpha^n=1@f$. The question is how to choose @f$\alpha^n@f$ in an
1296automatic fashion that satisfies these criteria.
1297
1298We do not want to review the literature on this topic here, but only briefly
1299mention that there are two fundamental approaches to the problem: backtracking
1300line search and trust region methods. The former is more widely used for
1301partial differential equations and essentially does the following:
1302- Compute a search direction
1303- See if the resulting residual of @f$u^n + \alpha^n\;\delta u^n@f$ with
1304 @f$\alpha^n=1@f$ is "substantially smaller" than that of @f$u^n@f$ alone.
1305- If so, then take @f$\alpha^n=1@f$.
1306- If not, try whether the residual is "substantially smaller" with
1307 @f$\alpha^n=2/3@f$.
1308- If so, then take @f$\alpha^n=2/3@f$.
1309- If not, try whether the residual is "substantially smaller" with
1310 @f$\alpha^n=(2/3)^2@f$.
1311- Etc.
1312One can of course choose other factors @f$r, r^2, \ldots@f$ than the @f$2/3,
1313(2/3)^2, \ldots@f$ chosen above, for @f$0<r<1@f$. It is obvious where the term
1314"backtracking" comes from: we try a long step, but if that doesn't work we try
1315a shorter step, and ever shorter step, etc. The function
1316<code>determine_step_length()</code> is written the way it is to support
1317exactly this kind of use case.
1318
1319Whether we accept a particular step length @f$\alpha^n@f$ depends on how we define
1320"substantially smaller". There are a number of ways to do so, but without
1321going into detail let us just mention that the most common ones are to use the
1322Wolfe and Armijo-Goldstein conditions. For these, one can show the following:
1323- There is always a step length @f$\alpha^n@f$ for which the conditions are
1324 satisfied, i.e., the iteration never gets stuck as long as the problem is
1325 convex.
1326- If we are close enough to the solution, then the conditions allow for
1327 @f$\alpha^n=1@f$, thereby enabling quadratic convergence.
1328
1329We will not dwell on this here any further but leave the implementation of
1330such algorithms as an exercise. We note, however, that when implemented
1331correctly then it is a common observation that most reasonably nonlinear
1332problems can be solved in anywhere between 5 and 15 Newton iterations to
1333engineering accuracy &mdash; substantially fewer than we need with the current
1334version of the program.
1335
1336More details on globalization methods including backtracking can be found,
1337for example, in @cite GNS08 and @cite NW99.
1338
1339A separate point, very much worthwhile making, however, is that in practice
1340the implementation of efficient nonlinear solvers is about as complicated as
1341the implementation of efficient finite element methods. One should not
1342attempt to reinvent the wheel by implementing all of the necessary steps
1343oneself. Substantial pieces of the puzzle are already available in
1344the LineMinimization::line_search() function and could be used to this end.
1345But, instead, just like building finite element solvers on libraries
1346such as deal.II, one should be building nonlinear solvers on libraries such
1347as [SUNDIALS](https://computing.llnl.gov/projects/sundials). In fact,
1348deal.II has interfaces to SUNDIALS and in particular to its nonlinear solver
1349sub-package KINSOL through the SUNDIALS::KINSOL class. It would not be
1350very difficult to base the current problem on that interface --
1351indeed, that is what @ref step_77 "step-77" does.
1352
1353
1354
1355<a name="step_15-Integratingmeshrefinementandnonlinearandlinearsolvers"></a><h4> Integrating mesh refinement and nonlinear and linear solvers </h4>
1356
1357
1358We currently do exactly 5 iterations on each mesh. But is this optimal? One
1359could ask the following questions:
1360- Maybe it is worthwhile doing more iterations on the initial meshes since
1361 there, computations are cheap.
1362- On the other hand, we do not want to do too many iterations on every mesh:
1363 yes, we could drive the residual to zero on every mesh, but that would only
1364 mean that the nonlinear iteration error is far smaller than the
1365 discretization error.
1366- Should we use solve the linear systems in each Newton step with higher or
1367 lower accuracy?
1368
1369Ultimately, what this boils down to is that we somehow need to couple the
1370discretization error on the current mesh with the nonlinear residual we want
1371to achieve with the Newton iterations on a given mesh, and to the linear
1372iteration we want to achieve with the CG method within each Newton
1373iterations.
1374
1375How to do this is, again, not entirely trivial, and we again leave it as a
1376future exercise.
1377
1378
1379
1380<a name="step_15-UsingautomaticdifferentiationtocomputetheJacobianmatrix"></a><h4> Using automatic differentiation to compute the Jacobian matrix </h4>
1381
1382
1383As outlined in the introduction, when solving a nonlinear problem of
1384the form
1385 @f[
1386 F(u) \dealcoloneq
1387 -\nabla \cdot \left( \frac{1}{\sqrt{1+|\nabla u|^{2}}}\nabla u \right)
1388 = 0
1389 @f]
1390we use a Newton iteration that requires us to repeatedly solve the
1391linear partial differential equation
1392 @f{align*}{
1393 F'(u^{n},\delta u^{n}) &=- F(u^{n})
1394 @f}
1395so that we can compute the update
1396 @f{align*}{
1397 u^{n+1}&=u^{n}+\alpha^n \delta u^{n}
1398 @f}
1399with the solution @f$\delta u^{n}@f$ of the Newton step. For the problem
1400here, we could compute the derivative @f$F'(u,\delta u)@f$ by hand and
1401obtained
1402 @f[
1403 F'(u,\delta u)
1404 =
1405 - \nabla \cdot \left( \frac{1}{\left(1+|\nabla u|^{2}\right)^{\frac{1}{2}}}\nabla
1406 \delta u \right) +
1407 \nabla \cdot \left( \frac{\nabla u \cdot
1408 \nabla \delta u}{\left(1+|\nabla u|^{2}\right)^{\frac{3}{2}}} \nabla u
1409 \right).
1410 @f]
1411But this is already a sizable expression that is cumbersome both to
1412derive and to implement. It is also, in some sense, duplicative: If we
1413implement what @f$F(u)@f$ is somewhere in the code, then @f$F'(u,\delta u)@f$
1414is not an independent piece of information but is something that, at
1415least in principle, a computer should be able to infer itself.
1416Wouldn't it be nice if that could actually happen? That is, if we
1417really only had to implement @f$F(u)@f$, and @f$F'(u,\delta u)@f$ was then somehow
1418done implicitly? That is in fact possible, and runs under the name
1419"automatic differentiation". @ref step_71 "step-71" discusses this very
1420concept in general terms, and @ref step_72 "step-72" illustrates how this can be
1421applied in practice for the very problem we are considering here.
1422
1423
1424<a name="step_15-StoringtheJacobianmatrixinlowerprecisionfloatingpointvariables"></a><h4> Storing the Jacobian matrix in lower-precision floating point variables </h4>
1425
1426
1427On modern computer systems, *accessing* data in main memory takes far
1428longer than *actually doing* something with it: We can do many floating
1429point operations for the time it takes to load one floating point
1430number from memory onto the processor. Unfortunately, when we do things
1431such as matrix-vector products, we only multiply each matrix entry once
1432with another number (the corresponding entry of the vector) and then we
1433add it to something else -- so two floating point operations for one
1434load. (Strictly speaking, we also have to load the corresponding vector
1435entry, but at least sometimes we get to re-use that vector entry in
1436doing the products that correspond to the next row of the matrix.) This
1437is a fairly low "arithmetic intensity", and consequently we spend most
1438of our time during matrix-vector products waiting for data to arrive
1439from memory rather than actually doing floating point operations.
1440
1441This is of course one of the rationales for the "matrix-free" approach to
1442solving linear systems (see @ref step_37 "step-37", for example). But if you don't quite
1443want to go all that way to change the structure of the program, then
1444here is a different approach: Storing the system matrix (the "Jacobian")
1445in single-precision instead of double precision floating point numbers
1446(i.e., using `float` instead of `double` as the data type). This reduces
1447the amount of memory necessary by a factor of 1.5 (each matrix entry
1448in a SparseMatrix object requires storing the column index -- 4 bytes --
1449and the actual value -- either 4 or 8 bytes), and consequently
1450will speed up matrix-vector products by a factor of around 1.5 as well because,
1451as pointed out above, most of the time is spent loading data from memory
1452and loading 2/3 the amount of data should be roughly 3/2 times as fast. All
1453of this could be done using SparseMatrix<float> as the data type
1454for the system matrix. (In principle, we would then also like it if
1455the SparseDirectUMFPACK solver we use in this program computes and
1456stores its sparse decomposition in `float` arithmetic. This is not
1457currently implemented, though could be done.)
1458
1459Of course, there is a downside to this: Lower precision data storage
1460also implies that we will not solve the linear system of the Newton
1461step as accurately as we might with `double` precision. At least
1462while we are far away from the solution of the nonlinear problem,
1463this may not be terrible: If we can do a Newton iteration in half
1464the time, we can afford to do a couple more Newton steps if the
1465search directions aren't as good.
1466But it turns out that even that isn't typically necessary: Both
1467theory and computational experience shows that it is entirely
1468sufficient to store the Jacobian matrix in single precision
1469*as long as one stores the right hand side in double precision*.
1470A great overview of why this is so, along with numerical
1471experiments that also consider "half precision" floating point
1472numbers, can be found in @cite Kelley2022 .
1473 *
1474 *
1475<a name="step_15-PlainProg"></a>
1476<h1> The plain program</h1>
1477@include "step-15.cc"
1478*/
void get_function_gradients(const ReadVector< Number > &fe_function, std::vector< Tensor< 1, spacedim, Number > > &gradients) const
static void estimate(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Quadrature< dim - 1 > &quadrature, const std::map< types::boundary_id, const Function< spacedim, Number > * > &neumann_bc, const ReadVector< Number > &solution, Vector< float > &error, const ComponentMask &component_mask={}, const Function< spacedim > *coefficients=nullptr, const unsigned int n_threads=numbers::invalid_unsigned_int, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id, const types::material_id material_id=numbers::invalid_material_id, const Strategy strategy=cell_diameter_over_24)
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
unsigned int n_active_cells() const
virtual bool prepare_coarsening_and_refinement()
Point< 2 > second
Definition grid_out.cc:4624
Point< 2 > first
Definition grid_out.cc:4623
__global__ void set(Number *val, const Number s, const size_type N)
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_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
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)
const Event initial
Definition event.cc:64
void approximate(const SynchronousIterators< std::tuple< typename DoFHandler< dim, spacedim >::active_cell_iterator, Vector< float >::iterator > > &cell, const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof_handler, const InputVector &solution, const unsigned int component)
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
void refine_and_coarsen_fixed_number(Triangulation< dim, spacedim > &triangulation, const Vector< Number > &criteria, const double top_fraction_of_cells, const double bottom_fraction_of_cells, const unsigned int max_n_cells=std::numeric_limits< unsigned int >::max())
@ matrix
Contents is actually a matrix.
std::pair< NumberType, unsigned int > line_search(const std::function< std::pair< NumberType, NumberType >(const NumberType x)> &func, const NumberType f0, const NumberType g0, const std::function< NumberType(const NumberType x_low, const NumberType f_low, const NumberType g_low, const NumberType x_hi, const NumberType f_hi, const NumberType g_hi, const FiniteSizeHistory< NumberType > &x_rec, const FiniteSizeHistory< NumberType > &f_rec, const FiniteSizeHistory< NumberType > &g_rec, const std::pair< NumberType, NumberType > bounds)> &interpolate, const NumberType a1, const NumberType eta=0.9, const NumberType mu=0.01, const NumberType a_max=std::numeric_limits< NumberType >::max(), const unsigned int max_evaluations=20, const bool debug_output=false)
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
Definition advection.h:74
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Definition divergence.h:471
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:191
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
VectorType::value_type * end(VectorType &V)
void interpolate_boundary_values(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const std::map< types::boundary_id, const Function< spacedim, number > * > &function_map, std::map< types::global_dof_index, number > &boundary_values, const ComponentMask &component_mask={})
bool check(const ConstraintKinds kind_in, const unsigned int dim)
int(& functions)(const void *v1, const void *v2)
static constexpr double PI
Definition numbers.h:259
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation