step-15.h

Go to the documentation of this file.

) const
 *     {
 *       return std::sin(2 * numbers::PI * (p[0] + p[1]));
 *     }
 *   
 * @endcode
 * 
 * 
 * <a name="ThecodeMinimalSurfaceProblemcodeclassimplementation"></a> 
 * <h3>The <code>MinimalSurfaceProblem</code> class implementation</h3>
 * 
 
 * 
 * 
 * <a name="MinimalSurfaceProblemMinimalSurfaceProblem"></a> 
 * <h4>MinimalSurfaceProblem::MinimalSurfaceProblem</h4>
 * 
 
 * 
 * The constructor and destructor of the class are the same as in the first
 * few tutorials.
 * 
 
 * 
 * 
 * @code
 *     template <int dim>
 *     MinimalSurfaceProblem<dim>::MinimalSurfaceProblem()
 *       : dof_handler(triangulation)
 *       , fe(2)
 *     {}
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="MinimalSurfaceProblemsetup_system"></a> 
 * <h4>MinimalSurfaceProblem::setup_system</h4>
 * 
 
 * 
 * As always in the setup-system function, we set up the variables of the
 * finite element method. There are same differences to @ref step_6 "step-6", because
 * there we start solving the PDE from scratch in every refinement cycle
 * whereas here we need to take the solution from the previous mesh onto the
 * current mesh. Consequently, we can't just reset solution vectors. The
 * argument passed to this function thus indicates whether we can
 * distributed degrees of freedom (plus compute constraints) and set the
 * solution vector to zero or whether this has happened elsewhere already
 * (specifically, in <code>refine_mesh()</code>).
 * 
 
 * 
 * 
 * @code
 *     template <int dim>
 *     void MinimalSurfaceProblem<dim>::setup_system(const bool initial_step)
 *     {
 *       if (initial_step)
 *         {
 *           dof_handler.distribute_dofs(fe);
 *           current_solution.reinit(dof_handler.n_dofs());
 *   
 *           hanging_node_constraints.clear();
 *           DoFTools::make_hanging_node_constraints(dof_handler,
 *                                                   hanging_node_constraints);
 *           hanging_node_constraints.close();
 *         }
 *   
 *   
 * @endcode
 * 
 * The remaining parts of the function are the same as in @ref step_6 "step-6".
 * 
 
 * 
 * 
 * @code
 *       newton_update.reinit(dof_handler.n_dofs());
 *       system_rhs.reinit(dof_handler.n_dofs());
 *   
 *       DynamicSparsityPattern dsp(dof_handler.n_dofs());
 *       DoFTools::make_sparsity_pattern(dof_handler, dsp);
 *   
 *       hanging_node_constraints.condense(dsp);
 *   
 *       sparsity_pattern.copy_from(dsp);
 *       system_matrix.reinit(sparsity_pattern);
 *     }
 *   
 * @endcode
 * 
 * 
 * <a name="MinimalSurfaceProblemassemble_system"></a> 
 * <h4>MinimalSurfaceProblem::assemble_system</h4>
 * 
 
 * 
 * This function does the same as in the previous tutorials except that now,
 * of course, the matrix and right hand side functions depend on the
 * previous iteration's solution. As discussed in the introduction, we need
 * to use zero boundary values for the Newton updates; we compute them at
 * the end of this function.
 *   
 
 * 
 * The top of the function contains the usual boilerplate code, setting up
 * the objects that allow us to evaluate shape functions at quadrature
 * points and temporary storage locations for the local matrices and
 * vectors, as well as for the gradients of the previous solution at the
 * quadrature points. We then start the loop over all cells:
 * 
 * @code
 *     template <int dim>
 *     void MinimalSurfaceProblem<dim>::assemble_system()
 *     {
 *       const QGauss<dim> quadrature_formula(fe.degree + 1);
 *   
 *       system_matrix = 0;
 *       system_rhs    = 0;
 *   
 *       FEValues<dim> fe_values(fe,
 *                               quadrature_formula,
 *                               update_gradients | update_quadrature_points |
 *                                 update_JxW_values);
 *   
 *       const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
 *       const unsigned int n_q_points    = quadrature_formula.size();
 *   
 *       FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
 *       Vector<double>     cell_rhs(dofs_per_cell);
 *   
 *       std::vector<Tensor<1, dim>> old_solution_gradients(n_q_points);
 *   
 *       std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
 *   
 *       for (const auto &cell : dof_handler.active_cell_iterators())
 *         {
 *           cell_matrix = 0;
 *           cell_rhs    = 0;
 *   
 *           fe_values.reinit(cell);
 *   
 * @endcode
 * 
 * For the assembly of the linear system, we have to obtain the values
 * of the previous solution's gradients at the quadrature
 * points. There is a standard way of doing this: the
 * FEValues::get_function_gradients function takes a vector that
 * represents a finite element field defined on a DoFHandler, and
 * evaluates the gradients of this field at the quadrature points of the
 * cell with which the FEValues object has last been reinitialized.
 * The values of the gradients at all quadrature points are then written
 * into the second argument:
 * 
 * @code
 *           fe_values.get_function_gradients(current_solution,
 *                                            old_solution_gradients);
 *   
 * @endcode
 * 
 * With this, we can then do the integration loop over all quadrature
 * points and shape functions.  Having just computed the gradients of
 * the old solution in the quadrature points, we are able to compute
 * the coefficients @f$a_{n}@f$ in these points.  The assembly of the
 * system itself then looks similar to what we always do with the
 * exception of the nonlinear terms, as does copying the results from
 * the local objects into the global ones:
 * 
 * @code
 *           for (unsigned int q = 0; q < n_q_points; ++q)
 *             {
 *               const double coeff =
 *                 1.0 / std::sqrt(1 + old_solution_gradients[q] *
 *                                       old_solution_gradients[q]);
 *   
 *               for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *                 {
 *                   for (unsigned int j = 0; j < dofs_per_cell; ++j)
 *                     cell_matrix(i, j) +=
 *                       (((fe_values.shape_grad(i, q)      // ((\nabla \phi_i
 *                          * coeff                         //   * a_n
 *                          * fe_values.shape_grad(j, q))   //   * \nabla \phi_j)
 *                         -                                //  -
 *                         (fe_values.shape_grad(i, q)      //  (\nabla \phi_i
 *                          * coeff * coeff * coeff         //   * a_n^3
 *                          * (fe_values.shape_grad(j, q)   //   * (\nabla \phi_j
 *                             * old_solution_gradients[q]) //      * \nabla u_n)
 *                          * old_solution_gradients[q]))   //   * \nabla u_n)))
 *                        * fe_values.JxW(q));              // * dx
 *   
 *                   cell_rhs(i) -= (fe_values.shape_grad(i, q)  // \nabla \phi_i
 *                                   * coeff                     // * a_n
 *                                   * old_solution_gradients[q] // * \nabla u_n
 *                                   * fe_values.JxW(q));        // * dx
 *                 }
 *             }
 *   
 *           cell->get_dof_indices(local_dof_indices);
 *           for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *             {
 *               for (unsigned int j = 0; j < dofs_per_cell; ++j)
 *                 system_matrix.add(local_dof_indices[i],
 *                                   local_dof_indices[j],
 *                                   cell_matrix(i, j));
 *   
 *               system_rhs(local_dof_indices[i]) += cell_rhs(i);
 *             }
 *         }
 *   
 * @endcode
 * 
 * Finally, we remove hanging nodes from the system and apply zero
 * boundary values to the linear system that defines the Newton updates
 * @f$\delta u^n@f$:
 * 
 * @code
 *       hanging_node_constraints.condense(system_matrix);
 *       hanging_node_constraints.condense(system_rhs);
 *   
 *       std::map<types::global_dof_index, double> boundary_values;
 *       VectorTools::interpolate_boundary_values(dof_handler,
 *                                                0,
 *                                                Functions::ZeroFunction<dim>(),
 *                                                boundary_values);
 *       MatrixTools::apply_boundary_values(boundary_values,
 *                                          system_matrix,
 *                                          newton_update,
 *                                          system_rhs);
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="MinimalSurfaceProblemsolve"></a> 
 * <h4>MinimalSurfaceProblem::solve</h4>
 * 
 
 * 
 * The solve function is the same as always. At the end of the solution
 * process we update the current solution by setting
 * @f$u^{n+1}=u^n+\alpha^n\;\delta u^n@f$.
 * 
 * @code
 *     template <int dim>
 *     void MinimalSurfaceProblem<dim>::solve()
 *     {
 *       SolverControl            solver_control(system_rhs.size(),
 *                                    system_rhs.l2_norm() * 1e-6);
 *       SolverCG<Vector<double>> solver(solver_control);
 *   
 *       PreconditionSSOR<SparseMatrix<double>> preconditioner;
 *       preconditioner.initialize(system_matrix, 1.2);
 *   
 *       solver.solve(system_matrix, newton_update, system_rhs, preconditioner);
 *   
 *       hanging_node_constraints.distribute(newton_update);
 *   
 *       const double alpha = determine_step_length();
 *       current_solution.add(alpha, newton_update);
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="MinimalSurfaceProblemrefine_mesh"></a> 
 * <h4>MinimalSurfaceProblem::refine_mesh</h4>
 * 
 
 * 
 * The first part of this function is the same as in @ref step_6 "step-6"... However,
 * after refining the mesh we have to transfer the old solution to the new
 * one which we do with the help of the SolutionTransfer class. The process
 * is slightly convoluted, so let us describe it in detail:
 * 
 * @code
 *     template <int dim>
 *     void MinimalSurfaceProblem<dim>::refine_mesh()
 *     {
 *       Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
 *   
 *       KellyErrorEstimator<dim>::estimate(
 *         dof_handler,
 *         QGauss<dim - 1>(fe.degree + 1),
 *         std::map<types::boundary_id, const Function<dim> *>(),
 *         current_solution,
 *         estimated_error_per_cell);
 *   
 *       GridRefinement::refine_and_coarsen_fixed_number(triangulation,
 *                                                       estimated_error_per_cell,
 *                                                       0.3,
 *                                                       0.03);
 *   
 * @endcode
 * 
 * Then we need an additional step: if, for example, you flag a cell that
 * is once more refined than its neighbor, and that neighbor is not
 * flagged for refinement, we would end up with a jump of two refinement
 * levels across a cell interface.  To avoid these situations, the library
 * will silently also have to refine the neighbor cell once. It does so by
 * calling the Triangulation::prepare_coarsening_and_refinement function
 * before actually doing the refinement and coarsening.  This function
 * flags a set of additional cells for refinement or coarsening, to
 * enforce rules like the one-hanging-node rule.  The cells that are
 * flagged for refinement and coarsening after calling this function are
 * exactly the ones that will actually be refined or coarsened. Usually,
 * you don't have to do this by hand
 * (Triangulation::execute_coarsening_and_refinement does this for
 * you). However, we need to initialize the SolutionTransfer class and it
 * needs to know the final set of cells that will be coarsened or refined
 * in order to store the data from the old mesh and transfer to the new
 * one. Thus, we call the function by hand:
 * 
 * @code
 *       triangulation.prepare_coarsening_and_refinement();
 *   
 * @endcode
 * 
 * With this out of the way, we initialize a SolutionTransfer object with
 * the present DoFHandler and attach the solution vector to it, followed
 * by doing the actual refinement and distribution of degrees of freedom
 * on the new mesh
 * 
 * @code
 *       SolutionTransfer<dim> solution_transfer(dof_handler);
 *       solution_transfer.prepare_for_coarsening_and_refinement(current_solution);
 *   
 *       triangulation.execute_coarsening_and_refinement();
 *   
 *       dof_handler.distribute_dofs(fe);
 *   
 * @endcode
 * 
 * Finally, we retrieve the old solution interpolated to the new
 * mesh. Since the SolutionTransfer function does not actually store the
 * values of the old solution, but rather indices, we need to preserve the
 * old solution vector until we have gotten the new interpolated
 * values. Thus, we have the new values written into a temporary vector,
 * and only afterwards write them into the solution vector object:
 * 
 * @code
 *       Vector<double> tmp(dof_handler.n_dofs());
 *       solution_transfer.interpolate(current_solution, tmp);
 *       current_solution = tmp;
 *   
 * @endcode
 * 
 * On the new mesh, there are different hanging nodes, for which we have to
 * compute constraints again, after throwing away previous content of the
 * object. To be on the safe side, we should then also make sure that the
 * current solution's vector entries satisfy the hanging node constraints
 * (see the discussion in the documentation of the SolutionTransfer class
 * for why this is necessary). We could do this by calling
 * `hanging_node_constraints.distribute(current_solution)` explicitly; we
 * omit this step because this will happen at the end of the call to
 * `set_boundary_values()` below, and it is not necessary to do it twice.
 * 
 * @code
 *       hanging_node_constraints.clear();
 *   
 *       DoFTools::make_hanging_node_constraints(dof_handler,
 *                                               hanging_node_constraints);
 *       hanging_node_constraints.close();
 *   
 * @endcode
 * 
 * Once we have the interpolated solution and all information about
 * hanging nodes, we have to make sure that the @f$u^n@f$ we now have
 * actually has the correct boundary values. As explained at the end of
 * the introduction, this is not automatically the case even if the
 * solution before refinement had the correct boundary values, and so we
 * have to explicitly make sure that it now has:
 * 
 * @code
 *       set_boundary_values();
 *   
 * @endcode
 * 
 * We end the function by updating all the remaining data structures,
 * indicating to <code>setup_dofs()</code> that this is not the first
 * go-around and that it needs to preserve the content of the solution
 * vector:
 * 
 * @code
 *       setup_system(false);
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="MinimalSurfaceProblemset_boundary_values"></a> 
 * <h4>MinimalSurfaceProblem::set_boundary_values</h4>
 * 
 
 * 
 * The next function ensures that the solution vector's entries respect the
 * boundary values for our problem.  Having refined the mesh (or just
 * started computations), there might be new nodal points on the
 * boundary. These have values that are simply interpolated from the
 * previous mesh in `refine_mesh()`, instead of the correct boundary
 * values. This is fixed up by setting all boundary nodes of the current
 * solution vector explicit to the right value.
 *   
 
 * 
 * There is one issue we have to pay attention to, though: If we have
 * a hanging node right next to a new boundary node, then its value
 * must also be adjusted to make sure that the finite element field
 * remains continuous. This is what the call in the last line of this
 * function does.
 * 
 * @code
 *     template <int dim>
 *     void MinimalSurfaceProblem<dim>::set_boundary_values()
 *     {
 *       std::map<types::global_dof_index, double> boundary_values;
 *       VectorTools::interpolate_boundary_values(dof_handler,
 *                                                0,
 *                                                BoundaryValues<dim>(),
 *                                                boundary_values);
 *       for (auto &boundary_value : boundary_values)
 *         current_solution(boundary_value.first) = boundary_value.second;
 *   
 *       hanging_node_constraints.distribute(current_solution);
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="MinimalSurfaceProblemcompute_residual"></a> 
 * <h4>MinimalSurfaceProblem::compute_residual</h4>
 * 
 
 * 
 * In order to monitor convergence, we need a way to compute the norm of the
 * (discrete) residual, i.e., the norm of the vector
 * @f$\left<F(u^n),\varphi_i\right>@f$ with @f$F(u)=-\nabla \cdot \left(
 * \frac{1}{\sqrt{1+|\nabla u|^{2}}}\nabla u \right)@f$ as discussed in the
 * introduction. It turns out that (although we don't use this feature in
 * the current version of the program) one needs to compute the residual
 * @f$\left<F(u^n+\alpha^n\;\delta u^n),\varphi_i\right>@f$ when determining
 * optimal step lengths, and so this is what we implement here: the function
 * takes the step length @f$\alpha^n@f$ as an argument. The original
 * functionality is of course obtained by passing a zero as argument.
 *   
 
 * 
 * In the function below, we first set up a vector for the residual, and
 * then a vector for the evaluation point @f$u^n+\alpha^n\;\delta u^n@f$. This
 * is followed by the same boilerplate code we use for all integration
 * operations:
 * 
 * @code
 *     template <int dim>
 *     double MinimalSurfaceProblem<dim>::compute_residual(const double alpha) const
 *     {
 *       Vector<double> residual(dof_handler.n_dofs());
 *   
 *       Vector<double> evaluation_point(dof_handler.n_dofs());
 *       evaluation_point = current_solution;
 *       evaluation_point.add(alpha, newton_update);
 *   
 *       const QGauss<dim> quadrature_formula(fe.degree + 1);
 *       FEValues<dim>     fe_values(fe,
 *                               quadrature_formula,
 *                               update_gradients | update_quadrature_points |
 *                                 update_JxW_values);
 *   
 *       const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
 *       const unsigned int n_q_points    = quadrature_formula.size();
 *   
 *       Vector<double>              cell_residual(dofs_per_cell);
 *       std::vector<Tensor<1, dim>> gradients(n_q_points);
 *   
 *       std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
 *   
 *       for (const auto &cell : dof_handler.active_cell_iterators())
 *         {
 *           cell_residual = 0;
 *           fe_values.reinit(cell);
 *   
 * @endcode
 * 
 * The actual computation is much as in
 * <code>assemble_system()</code>. We first evaluate the gradients of
 * @f$u^n+\alpha^n\,\delta u^n@f$ at the quadrature points, then compute
 * the coefficient @f$a_n@f$, and then plug it all into the formula for
 * the residual:
 * 
 * @code
 *           fe_values.get_function_gradients(evaluation_point, gradients);
 *   
 *   
 *           for (unsigned int q = 0; q < n_q_points; ++q)
 *             {
 *               const double coeff =
 *                 1. / std::sqrt(1 + gradients[q] * gradients[q]);
 *   
 *               for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *                 cell_residual(i) -= (fe_values.shape_grad(i, q) // \nabla \phi_i
 *                                      * coeff                    // * a_n
 *                                      * gradients[q]             // * \nabla u_n
 *                                      * fe_values.JxW(q));       // * dx
 *             }
 *   
 *           cell->get_dof_indices(local_dof_indices);
 *           for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *             residual(local_dof_indices[i]) += cell_residual(i);
 *         }
 *   
 * @endcode
 * 
 * At the end of this function we also have to deal with the hanging node
 * constraints and with the issue of boundary values. With regard to the
 * latter, we have to set to zero the elements of the residual vector for
 * all entries that correspond to degrees of freedom that sit at the
 * boundary. The reason is that because the value of the solution there is
 * fixed, they are of course no "real" degrees of freedom and so, strictly
 * speaking, we shouldn't have assembled entries in the residual vector
 * for them. However, as we always do, we want to do exactly the same
 * thing on every cell and so we didn't want to deal with the question
 * of whether a particular degree of freedom sits at the boundary in the
 * integration above. Rather, we will simply set to zero these entries
 * after the fact. To this end, we need to determine which degrees
 * of freedom do in fact belong to the boundary and then loop over all of
 * those and set the residual entry to zero. This happens in the following
 * lines which we have already seen used in @ref step_11 "step-11", using the appropriate
 * function from namespace DoFTools:
 * 
 * @code
 *       hanging_node_constraints.condense(residual);
 *   
 *       for (const types::global_dof_index i :
 *            DoFTools::extract_boundary_dofs(dof_handler))
 *         residual(i) = 0;
 *   
 * @endcode
 * 
 * At the end of the function, we return the norm of the residual:
 * 
 * @code
 *       return residual.l2_norm();
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="MinimalSurfaceProblemdetermine_step_length"></a> 
 * <h4>MinimalSurfaceProblem::determine_step_length</h4>
 * 
 
 * 
 * As discussed in the introduction, Newton's method frequently does not
 * converge if we always take full steps, i.e., compute @f$u^{n+1}=u^n+\delta
 * u^n@f$. Rather, one needs a damping parameter (step length) @f$\alpha^n@f$ and
 * set @f$u^{n+1}=u^n+\alpha^n\delta u^n@f$. This function is the one called
 * to compute @f$\alpha^n@f$.
 *   
 
 * 
 * Here, we simply always return 0.1. This is of course a sub-optimal
 * choice: ideally, what one wants is that the step size goes to one as we
 * get closer to the solution, so that we get to enjoy the rapid quadratic
 * convergence of Newton's method. We will discuss better strategies below
 * in the results section, and @ref step_77 "step-77" also covers this aspect.
 * 
 * @code
 *     template <int dim>
 *     double MinimalSurfaceProblem<dim>::determine_step_length() const
 *     {
 *       return 0.1;
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="MinimalSurfaceProblemoutput_results"></a> 
 * <h4>MinimalSurfaceProblem::output_results</h4>
 * 
 
 * 
 * This last function to be called from `run()` outputs the current solution
 * (and the Newton update) in graphical form as a VTU file. It is entirely the
 * same as what has been used in previous tutorials.
 * 
 * @code
 *     template <int dim>
 *     void MinimalSurfaceProblem<dim>::output_results(
 *       const unsigned int refinement_cycle) const
 *     {
 *       DataOut<dim> data_out;
 *   
 *       data_out.attach_dof_handler(dof_handler);
 *       data_out.add_data_vector(current_solution, "solution");
 *       data_out.add_data_vector(newton_update, "update");
 *       data_out.build_patches();
 *   
 *       const std::string filename =
 *         "solution-" + Utilities::int_to_string(refinement_cycle, 2) + ".vtu";
 *       std::ofstream output(filename);
 *       data_out.write_vtu(output);
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="MinimalSurfaceProblemrun"></a> 
 * <h4>MinimalSurfaceProblem::run</h4>
 * 
 
 * 
 * In the run function, we build the first grid and then have the top-level
 * logic for the Newton iteration.
 *   
 
 * 
 * As described in the introduction, the domain is the unit disk around
 * the origin, created in the same way as shown in @ref step_6 "step-6". The mesh is
 * globally refined twice followed later on by several adaptive cycles.
 *   
 
 * 
 * Before starting the Newton loop, we also need to do a bit of
 * setup work: We need to create the basic data structures and
 * ensure that the first Newton iterate already has the correct
 * boundary values, as discussed in the introduction.
 * 
 * @code
 *     template <int dim>
 *     void MinimalSurfaceProblem<dim>::run()
 *     {
 *       GridGenerator::hyper_ball(triangulation);
 *       triangulation.refine_global(2);
 *   
 *       setup_system(/*first time=*/true);
 *       set_boundary_values();
 *   
 * @endcode
 * 
 * The Newton iteration starts next. We iterate until the (norm of the)
 * residual computed at the end of the previous iteration is less than
 * @f$10^{-3}@f$, as checked at the end of the `do { ... } while` loop that
 * starts here. Because we don't have a reasonable value to initialize
 * the variable, we just use the largest value that can be represented
 * as a `double`.
 * 
 * @code
 *       double       last_residual_norm = std::numeric_limits<double>::max();
 *       unsigned int refinement_cycle   = 0;
 *       do
 *         {
 *           std::cout << "Mesh refinement step " << refinement_cycle << std::endl;
 *   
 *           if (refinement_cycle != 0)
 *             refine_mesh();
 *   
 * @endcode
 * 
 * On every mesh we do exactly five Newton steps. We print the initial
 * residual here and then start the iterations on this mesh.
 *         
 
 * 
 * In every Newton step the system matrix and the right hand side have
 * to be computed first, after which we store the norm of the right
 * hand side as the residual to check against when deciding whether to
 * stop the iterations. We then solve the linear system (the function
 * also updates @f$u^{n+1}=u^n+\alpha^n\;\delta u^n@f$) and output the
 * norm of the residual at the end of this Newton step.
 *         
 
 * 
 * After the end of this loop, we then also output the solution on the
 * current mesh in graphical form and increment the counter for the
 * mesh refinement cycle.
 * 
 * @code
 *           std::cout << "  Initial residual: " << compute_residual(0) << std::endl;
 *   
 *           for (unsigned int inner_iteration = 0; inner_iteration < 5;
 *                ++inner_iteration)
 *             {
 *               assemble_system();
 *               last_residual_norm = system_rhs.l2_norm();
 *   
 *               solve();
 *   
 *               std::cout << "  Residual: " << compute_residual(0) << std::endl;
 *             }
 *   
 *           output_results(refinement_cycle);
 *   
 *           ++refinement_cycle;
 *           std::cout << std::endl;
 *         }
 *       while (last_residual_norm > 1e-2);
 *     }
 *   } // namespace Step15
 *   
 * @endcode
 * 
 * 
 * <a name="Themainfunction"></a> 
 * <h4>The main function</h4>
 * 
 
 * 
 * Finally the main function. This follows the scheme of all other main
 * functions:
 * 
 * @code
 *   int main()
 *   {
 *     try
 *       {
 *         using namespace Step15;
 *   
 *         MinimalSurfaceProblem<2> laplace_problem_2d;
 *         laplace_problem_2d.run();
 *       }
 *     catch (std::exception &exc)
 *       {
 *         std::cerr << std::endl
 *                   << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         std::cerr << "Exception on processing: " << std::endl
 *                   << exc.what() << std::endl
 *                   << "Aborting!" << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *   
 *         return 1;
 *       }
 *     catch (...)
 *       {
 *         std::cerr << std::endl
 *                   << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         std::cerr << "Unknown exception!" << std::endl
 *                   << "Aborting!" << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         return 1;
 *       }
 *     return 0;
 *   }
 * @endcode
<a name="Results"></a><h1>Results</h1>
 
 
 
The output of the program looks as follows:
@code
Mesh refinement step 0
  Initial residual: 1.53143
  Residual: 1.08746
  Residual: 0.966748
  Residual: 0.859602
  Residual: 0.766462
  Residual: 0.685475
 
Mesh refinement step 1
  Initial residual: 0.868959
  Residual: 0.762125
  Residual: 0.677792
  Residual: 0.605762
  Residual: 0.542748
  Residual: 0.48704
 
Mesh refinement step 2
  Initial residual: 0.426445
  Residual: 0.382731
  Residual: 0.343865
  Residual: 0.30918
  Residual: 0.278147
  Residual: 0.250327
 
Mesh refinement step 3
  Initial residual: 0.282026
  Residual: 0.253146
  Residual: 0.227414
  Residual: 0.20441
  Residual: 0.183803
  Residual: 0.165319
 
Mesh refinement step 4
  Initial residual: 0.154404
  Residual: 0.138723
  Residual: 0.124694
  Residual: 0.112124
  Residual: 0.100847
  Residual: 0.0907222
 
....
@endcode
 
Obviously, the scheme converges, if not very fast. We will come back to
strategies for accelerating the method below.
 
One can visualize the solution after each set of five Newton
iterations, i.e., on each of the meshes on which we approximate the
solution. This yields the following set of images:
 
<div class="twocolumn" style="width: 80%">
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step_15_solution_1.png"
         alt="Solution after zero cycles with contour lines." width="230" height="273">
  </div>
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step_15_solution_2.png"
         alt="Solution after one cycle with contour lines." width="230" height="273">
  </div>
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step_15_solution_3.png"
         alt="Solution after two cycles with contour lines." width="230" height="273">
  </div>
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step_15_solution_4.png"
         alt="Solution after three cycles with contour lines." width="230" height="273">
  </div>
</div>
 
It is clearly visible, that the solution minimizes the surface
after each refinement. The solution converges to a picture one
would imagine a soap bubble to be that is located inside a wire loop
that is bent like
the boundary. Also it is visible, how the boundary
is smoothed out after each refinement. On the coarse mesh,
the boundary doesn't look like a sine, whereas it does the
finer the mesh gets.
 
The mesh is mostly refined near the boundary, where the solution
increases or decreases strongly, whereas it is coarsened on
the inside of the domain, where nothing interesting happens,
because there isn't much change in the solution. The ninth
solution and mesh are shown here:
 
<div class="onecolumn" style="width: 60%">
  <div>
    <img src="https://www.dealii.org/images/steps/developer/step_15_solution_9.png"
         alt="Grid and solution of the ninth cycle with contour lines." width="507" height="507">
  </div>
</div>
 
 
 
<a name="extensions"></a>
<a name="Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
 
 
The program shows the basic structure of a solver for a nonlinear, stationary
problem. However, it does not converge particularly fast, for good reasons:
 
- The program always takes a step size of 0.1. This precludes the rapid,
  quadratic convergence for which Newton's method is typically chosen.
- It does not connect the nonlinear iteration with the mesh refinement
  iteration.
 
Obviously, a better program would have to address these two points.
We will discuss them in the following.
 
 
<a name="Steplengthcontrol"></a><h4> Step length control </h4>
 
 
Newton's method has two well known properties:
- It may not converge from arbitrarily chosen starting points. Rather, a
  starting point has to be close enough to the solution to guarantee
  convergence. However, we can enlarge the area from which Newton's method
  converges by damping the iteration using a <i>step length</i> 0<@f$\alpha^n\le
  1@f$.
- It exhibits rapid convergence of quadratic order if (i) the step length is
  chosen as @f$\alpha^n=1@f$, and (ii) it does in fact converge with this choice
  of step length.
 
A consequence of these two observations is that a successful strategy is to
choose @f$\alpha^n<1@f$ for the initial iterations until the iterate has come
close enough to allow for convergence with full step length, at which point we
want to switch to @f$\alpha^n=1@f$. The question is how to choose @f$\alpha^n@f$ in an
automatic fashion that satisfies these criteria.
 
We do not want to review the literature on this topic here, but only briefly
mention that there are two fundamental approaches to the problem: backtracking
line search and trust region methods. The former is more widely used for
partial differential equations and essentially does the following:
- Compute a search direction
- See if the resulting residual of @f$u^n + \alpha^n\;\delta u^n@f$ with
  @f$\alpha^n=1@f$ is "substantially smaller" than that of @f$u^n@f$ alone.
- If so, then take @f$\alpha^n=1@f$.
- If not, try whether the residual is "substantially smaller" with
  @f$\alpha^n=2/3@f$.
- If so, then take @f$\alpha^n=2/3@f$.
- If not, try whether the residual is "substantially smaller" with
  @f$\alpha^n=(2/3)^2@f$.
- Etc.
One can of course choose other factors @f$r, r^2, \ldots@f$ than the @f$2/3,
(2/3)^2, \ldots@f$ chosen above, for @f$0<r<1@f$. It is obvious where the term
"backtracking" comes from: we try a long step, but if that doesn't work we try
a shorter step, and ever shorter step, etc. The function
<code>determine_step_length()</code> is written the way it is to support
exactly this kind of use case.
 
Whether we accept a particular step length @f$\alpha^n@f$ depends on how we define
"substantially smaller". There are a number of ways to do so, but without
going into detail let us just mention that the most common ones are to use the
Wolfe and Armijo-Goldstein conditions. For these, one can show the following:
- There is always a step length @f$\alpha^n@f$ for which the conditions are
  satisfied, i.e., the iteration never gets stuck as long as the problem is
  convex.
- If we are close enough to the solution, then the conditions allow for
  @f$\alpha^n=1@f$, thereby enabling quadratic convergence.
 
We will not dwell on this here any further but leave the implementation of
such algorithms as an exercise. We note, however, that when implemented
correctly then it is a common observation that most reasonably nonlinear
problems can be solved in anywhere between 5 and 15 Newton iterations to
engineering accuracy &mdash; substantially fewer than we need with the current
version of the program.
 
More details on globalization methods including backtracking can be found,
for example, in @cite GNS08 and @cite NW99.
 
A separate point, very much worthwhile making, however, is that in practice
the implementation of efficient nonlinear solvers is about as complicated as
the implementation of efficient finite element methods. One should not
attempt to reinvent the wheel by implementing all of the necessary steps
oneself. Substantial pieces of the puzzle are already available in
the LineMinimization::line_search() function and could be used to this end.
But, instead, just like building finite element solvers on libraries
such as deal.II, one should be building nonlinear solvers on libraries such
as [SUNDIALS](https://computing.llnl.gov/projects/sundials). In fact,
deal.II has interfaces to SUNDIALS and in particular to its nonlinear solver
sub-package KINSOL through the SUNDIALS::KINSOL class. It would not be
very difficult to base the current problem on that interface --
indeed, that is what @ref step_77 "step-77" does.
 
 
 
<a name="Integratingmeshrefinementandnonlinearandlinearsolvers"></a><h4> Integrating mesh refinement and nonlinear and linear solvers </h4>
 
 
We currently do exactly 5 iterations on each mesh. But is this optimal? One
could ask the following questions:
- Maybe it is worthwhile doing more iterations on the initial meshes since
  there, computations are cheap.
- On the other hand, we do not want to do too many iterations on every mesh:
  yes, we could drive the residual to zero on every mesh, but that would only
  mean that the nonlinear iteration error is far smaller than the
  discretization error.
- Should we use solve the linear systems in each Newton step with higher or
  lower accuracy?
 
Ultimately, what this boils down to is that we somehow need to couple the
discretization error on the current mesh with the nonlinear residual we want
to achieve with the Newton iterations on a given mesh, and to the linear
iteration we want to achieve with the CG method within each Newton
iterations.
 
How to do this is, again, not entirely trivial, and we again leave it as a
future exercise.
 
 
 
<a name="UsingautomaticdifferentiationtocomputetheJacobianmatrix"></a><h4> Using automatic differentiation to compute the Jacobian matrix </h4>
 
 
As outlined in the introduction, when solving a nonlinear problem of
the form
  @f[
    F(u) \dealcoloneq
    -\nabla \cdot \left( \frac{1}{\sqrt{1+|\nabla u|^{2}}}\nabla u \right)
    = 0
  @f]
we use a Newton iteration that requires us to repeatedly solve the
linear partial differential equation
  @f{align*}
    F'(u^{n},\delta u^{n}) &=- F(u^{n})
  @f}
so that we can compute the update
  @f{align*}
    u^{n+1}&=u^{n}+\alpha^n \delta u^{n}
  @f}
with the solution @f$\delta u^{n}@f$ of the Newton step. For the problem
here, we could compute the derivative @f$F'(u,\delta u)@f$ by hand and
obtained
  @f[
  F'(u,\delta u)
  =
  - \nabla \cdot \left( \frac{1}{\left(1+|\nabla u|^{2}\right)^{\frac{1}{2}}}\nabla
  \delta u \right) +
  \nabla \cdot \left( \frac{\nabla u \cdot
  \nabla \delta u}{\left(1+|\nabla u|^{2}\right)^{\frac{3}{2}}} \nabla u
  \right).
  @f]
But this is already a sizable expression that is cumbersome both to
derive and to implement. It is also, in some sense, duplicative: If we
implement what @f$F(u)@f$ is somewhere in the code, then @f$F'(u,\delta u)@f$
is not an independent piece of information but is something that, at
least in principle, a computer should be able to infer itself.
Wouldn't it be nice if that could actually happen? That is, if we
really only had to implement @f$F(u)@f$, and @f$F'(u,\delta u)@f$ was then somehow
done implicitly? That is in fact possible, and runs under the name
"automatic differentiation". @ref step_71 "step-71" discusses this very
concept in general terms, and @ref step_72 "step-72" illustrates how this can be
applied in practice for the very problem we are considering here.
 
 
<a name="StoringtheJacobianmatrixinlowerprecisionfloatingpointvariables"></a><h4> Storing the Jacobian matrix in lower-precision floating point variables </h4>
 
 
On modern computer systems, *accessing* data in main memory takes far
longer than *actually doing* something with it: We can do many floating
point operations for the time it takes to load one floating point
number from memory onto the processor. Unfortunately, when we do things
such as matrix-vector products, we only multiply each matrix entry once
with another number (the corresponding entry of the vector) and then we
add it to something else -- so two floating point operations for one
load. (Strictly speaking, we also have to load the corresponding vector
entry, but at least sometimes we get to re-use that vector entry in
doing the products that correspond to the next row of the matrix.) This
is a fairly low "arithmetic intensity", and consequently we spend most
of our time during matrix-vector products waiting for data to arrive
from memory rather than actually doing floating point operations.
 
This is of course one of the rationales for the "matrix-free" approach to
solving linear systems (see @ref step_37 "step-37", for example). But if you don't quite
want to go all that way to change the structure of the program, then
here is a different approach: Storing the system matrix (the "Jacobian")
in single-precision instead of double precision floating point numbers
(i.e., using `float` instead of `double` as the data type). This reduces
the amount of memory necessary by a factor of 1.5 (each matrix entry
in a SparseMatrix object requires storing the column index -- 4 bytes --
and the actual value -- either 4 or 8 bytes), and consequently
will speed up matrix-vector products by a factor of around 1.5 as well because,
as pointed out above, most of the time is spent loading data from memory
and loading 2/3 the amount of data should be roughly 3/2 times as fast. All
of this could be done using SparseMatrix<float> as the data type
for the system matrix. (In principle, we would then also like it if
the SparseDirectUMFPACK solver we use in this program computes and
stores its sparse decomposition in `float` arithmetic. This is not
currently implemented, though could be done.)
 
Of course, there is a downside to this: Lower precision data storage
also implies that we will not solve the linear system of the Newton
step as accurately as we might with `double` precision. At least
while we are far away from the solution of the nonlinear problem,
this may not be terrible: If we can do a Newton iteration in half
the time, we can afford to do a couple more Newton steps if the
search directions aren't as good.
But it turns out that even that isn't typically necessary: Both
theory and computational experience shows that it is entirely
sufficient to store the Jacobian matrix in single precision
*as long as one stores the right hand side in double precision*.
A great overview of why this is so, along with numerical
experiments that also consider "half precision" floating point
numbers, can be found in @cite Kelley2022 .
 *
 *
<a name="PlainProg"></a>
<h1> The plain program</h1>
@include "step-15.cc"
*/