step-77.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="Constructorandsetupfunctions"></a> 
 * <h4>Constructor and set up functions</h4>
 * 
 
 * 
 * The following few functions are also essentially copies of what
 * @ref step_15 "step-15" already does, and so there is little to discuss.
 * 
 * @code
 *     template <int dim>
 *     MinimalSurfaceProblem<dim>::MinimalSurfaceProblem()
 *       : dof_handler(triangulation)
 *       , fe(1)
 *       , computing_timer(std::cout, TimerOutput::never, TimerOutput::wall_times)
 *     {}
 *   
 *   
 *   
 *     template <int dim>
 *     void MinimalSurfaceProblem<dim>::setup_system(const bool initial_step)
 *     {
 *       TimerOutput::Scope t(computing_timer, "set up");
 *   
 *       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();
 *         }
 *   
 *       DynamicSparsityPattern dsp(dof_handler.n_dofs());
 *       DoFTools::make_sparsity_pattern(dof_handler, dsp);
 *   
 *       hanging_node_constraints.condense(dsp);
 *   
 *       sparsity_pattern.copy_from(dsp);
 *       jacobian_matrix.reinit(sparsity_pattern);
 *       jacobian_matrix_factorization.reset();
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="AssemblingandfactorizingtheJacobianmatrix"></a> 
 * <h4>Assembling and factorizing the Jacobian matrix</h4>
 * 
 
 * 
 * The following function is then responsible for assembling and factorizing
 * the Jacobian matrix. The first half of the function is in essence the
 * `assemble_system()` function of @ref step_15 "step-15", except that it does not deal with
 * also forming a right hand side vector (i.e., the residual) since we do not
 * always have to do these operations at the same time.
 *   
 
 * 
 * We put the whole assembly functionality into a code block enclosed by curly
 * braces so that we can use a TimerOutput::Scope variable to measure how much
 * time is spent in this code block, excluding everything that happens in this
 * function after the matching closing brace `}`.
 * 
 * @code
 *     template <int dim>
 *     void MinimalSurfaceProblem<dim>::compute_and_factorize_jacobian(
 *       const Vector<double> &evaluation_point)
 *     {
 *       {
 *         TimerOutput::Scope t(computing_timer, "assembling the Jacobian");
 *   
 *         std::cout << "  Computing Jacobian matrix" << std::endl;
 *   
 *         const QGauss<dim> quadrature_formula(fe.degree + 1);
 *   
 *         jacobian_matrix = 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);
 *   
 *         std::vector<Tensor<1, dim>> evaluation_point_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;
 *   
 *             fe_values.reinit(cell);
 *   
 *             fe_values.get_function_gradients(evaluation_point,
 *                                              evaluation_point_gradients);
 *   
 *             for (unsigned int q = 0; q < n_q_points; ++q)
 *               {
 *                 const double coeff =
 *                   1.0 / std::sqrt(1 + evaluation_point_gradients[q] *
 *                                         evaluation_point_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
 *                             * evaluation_point_gradients[q]) //      * \nabla u_n)
 *                            * evaluation_point_gradients[q])) //   * \nabla u_n)))
 *                          * fe_values.JxW(q));                // * dx
 *                   }
 *               }
 *   
 *             cell->get_dof_indices(local_dof_indices);
 *             hanging_node_constraints.distribute_local_to_global(cell_matrix,
 *                                                                 local_dof_indices,
 *                                                                 jacobian_matrix);
 *           }
 *   
 *         std::map<types::global_dof_index, double> boundary_values;
 *         VectorTools::interpolate_boundary_values(dof_handler,
 *                                                  0,
 *                                                  Functions::ZeroFunction<dim>(),
 *                                                  boundary_values);
 *         Vector<double> dummy_solution(dof_handler.n_dofs());
 *         Vector<double> dummy_rhs(dof_handler.n_dofs());
 *         MatrixTools::apply_boundary_values(boundary_values,
 *                                            jacobian_matrix,
 *                                            dummy_solution,
 *                                            dummy_rhs);
 *       }
 *   
 * @endcode
 * 
 * The second half of the function then deals with factorizing the
 * so-computed matrix. To do this, we first create a new SparseDirectUMFPACK
 * object and by assigning it to the member variable
 * `jacobian_matrix_factorization`, we also destroy whatever object that
 * pointer previously pointed to (if any). Then we tell the object to
 * factorize the Jacobian.
 *     
 
 * 
 * As above, we enclose this block of code into curly braces and use a timer
 * to assess how long this part of the program takes.
 *     
 
 * 
 * (Strictly speaking, we don't actually need the matrix any more after we
 * are done here, and could throw the matrix object away. A code intended to
 * be memory efficient would do this, and only create the matrix object in
 * this function, rather than as a member variable of the surrounding class.
 * We omit this step here because using the same coding style as in previous
 * tutorial programs breeds familiarity with the common style and helps make
 * these tutorial programs easier to read.)
 * 
 * @code
 *       {
 *         TimerOutput::Scope t(computing_timer, "factorizing the Jacobian");
 *   
 *         std::cout << "  Factorizing Jacobian matrix" << std::endl;
 *   
 *         jacobian_matrix_factorization = std::make_unique<SparseDirectUMFPACK>();
 *         jacobian_matrix_factorization->factorize(jacobian_matrix);
 *       }
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="Computingtheresidualvector"></a> 
 * <h4>Computing the residual vector</h4>
 * 
 
 * 
 * The second part of what `assemble_system()` used to do in @ref step_15 "step-15" is
 * computing the residual vector, i.e., the right hand side vector of the
 * Newton linear systems. We have broken this out of the previous function,
 * but the following function will be easy to understand if you understood
 * what `assemble_system()` in @ref step_15 "step-15" did. Importantly, however, we need to
 * compute the residual not linearized around the current solution vector, but
 * whatever we get from KINSOL. This is necessary for operations such as line
 * search where we want to know what the residual @f$F(U^k + \alpha_k \delta
 * U^K)@f$ is for different values of @f$\alpha_k@f$; KINSOL in those cases simply
 * gives us the argument to the function @f$F@f$ and we then compute the residual
 * @f$F(\cdot)@f$ at this point.
 *   
 
 * 
 * The function prints the norm of the so-computed residual at the end as a
 * way for us to follow along the progress of the program.
 * 
 * @code
 *     template <int dim>
 *     void MinimalSurfaceProblem<dim>::compute_residual(
 *       const Vector<double> &evaluation_point,
 *       Vector<double> &      residual)
 *     {
 *       TimerOutput::Scope t(computing_timer, "assembling the residual");
 *   
 *       std::cout << "  Computing residual vector..." << std::flush;
 *   
 *       residual = 0.0;
 *   
 *       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>> evaluation_point_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);
 *   
 *           fe_values.get_function_gradients(evaluation_point,
 *                                            evaluation_point_gradients);
 *   
 *   
 *           for (unsigned int q = 0; q < n_q_points; ++q)
 *             {
 *               const double coeff =
 *                 1.0 / std::sqrt(1 + evaluation_point_gradients[q] *
 *                                       evaluation_point_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
 *                    * evaluation_point_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);
 *         }
 *   
 *       hanging_node_constraints.condense(residual);
 *   
 *       for (const types::global_dof_index i :
 *            DoFTools::extract_boundary_dofs(dof_handler))
 *         residual(i) = 0;
 *   
 *       for (const types::global_dof_index i :
 *            DoFTools::extract_hanging_node_dofs(dof_handler))
 *         residual(i) = 0;
 *   
 *       std::cout << " norm=" << residual.l2_norm() << std::endl;
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="SolvinglinearsystemswiththeJacobianmatrix"></a> 
 * <h4>Solving linear systems with the Jacobian matrix</h4>
 * 
 
 * 
 * Next up is the function that implements the solution of a linear system
 * with the Jacobian matrix. Since we have already factored the matrix when we
 * built the matrix, solving a linear system comes down to applying the
 * inverse matrix to the given right hand side vector: This is what the
 * SparseDirectUMFPACK::vmult() function does that we use here. Following
 * this, we have to make sure that we also address the values of hanging nodes
 * in the solution vector, and this is done using
 * AffineConstraints::distribute().
 *   
 
 * 
 * The function takes an additional, but unused, argument `tolerance` that
 * indicates how accurately we have to solve the linear system. The meaning of
 * this argument is discussed in the introduction in the context of the
 * "Eisenstat Walker trick", but since we are using a direct rather than an
 * iterative solver, we are not using this opportunity to solve linear systems
 * only inexactly.
 * 
 * @code
 *     template <int dim>
 *     void MinimalSurfaceProblem<dim>::solve(const Vector<double> &rhs,
 *                                            Vector<double> &      solution,
 *                                            const double /*tolerance*/)
 *     {
 *       TimerOutput::Scope t(computing_timer, "linear system solve");
 *   
 *       std::cout << "  Solving linear system" << std::endl;
 *   
 *       jacobian_matrix_factorization->vmult(solution, rhs);
 *   
 *       hanging_node_constraints.distribute(solution);
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="Refiningthemeshsettingboundaryvaluesandgeneratinggraphicaloutput"></a> 
 * <h4>Refining the mesh, setting boundary values, and generating graphical output</h4>
 * 
 
 * 
 * The following three functions are again simply copies of the ones in
 * @ref step_15 "step-15":
 * 
 * @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);
 *   
 *       triangulation.prepare_coarsening_and_refinement();
 *   
 *       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);
 *   
 *       Vector<double> tmp(dof_handler.n_dofs());
 *       solution_transfer.interpolate(current_solution, tmp);
 *       current_solution = std::move(tmp);
 *   
 *       hanging_node_constraints.clear();
 *   
 *       DoFTools::make_hanging_node_constraints(dof_handler,
 *                                               hanging_node_constraints);
 *       hanging_node_constraints.close();
 *   
 *       hanging_node_constraints.distribute(current_solution);
 *   
 *       set_boundary_values();
 *   
 *       setup_system(/*initial_step=*/false);
 *     }
 *   
 *   
 *   
 *     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 (const auto &boundary_value : boundary_values)
 *         current_solution(boundary_value.first) = boundary_value.second;
 *   
 *       hanging_node_constraints.distribute(current_solution);
 *     }
 *   
 *   
 *   
 *     template <int dim>
 *     void MinimalSurfaceProblem<dim>::output_results(
 *       const unsigned int refinement_cycle)
 *     {
 *       TimerOutput::Scope t(computing_timer, "graphical output");
 *   
 *       DataOut<dim> data_out;
 *   
 *       data_out.attach_dof_handler(dof_handler);
 *       data_out.add_data_vector(current_solution, "solution");
 *       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="Therunfunctionandtheoveralllogicoftheprogram"></a> 
 * <h4>The run() function and the overall logic of the program</h4>
 * 
 
 * 
 * The only function that *really* is interesting in this program is the one
 * that drives the overall algorithm of starting on a coarse mesh, doing some
 * mesh refinement cycles, and on each mesh using KINSOL to find the solution
 * of the nonlinear algebraic equation we obtain from discretization on this
 * mesh. The `refine_mesh()` function above makes sure that the solution on
 * one mesh is used as the starting guess on the next mesh. We also use a
 * TimerOutput object to measure how much time every operation on each mesh
 * costs, and reset the timer at the beginning of each cycle.
 *   
 
 * 
 * As discussed in the introduction, it is not necessary to solve problems on
 * coarse meshes particularly accurately since these will only solve as
 * starting guesses for the next mesh. As a consequence, we will use a target
 * tolerance of
 * @f$\tau=10^{-3} \frac{1}{10^k}@f$ for the @f$k@f$th mesh refinement cycle.
 *   
 
 * 
 * All of this is encoded in the first part of this function:
 * 
 * @code
 *     template <int dim>
 *     void MinimalSurfaceProblem<dim>::run()
 *     {
 *       GridGenerator::hyper_ball(triangulation);
 *       triangulation.refine_global(2);
 *   
 *       setup_system(/*initial_step=*/true);
 *       set_boundary_values();
 *   
 *       for (unsigned int refinement_cycle = 0; refinement_cycle < 6;
 *            ++refinement_cycle)
 *         {
 *           computing_timer.reset();
 *           std::cout << "Mesh refinement step " << refinement_cycle << std::endl;
 *   
 *           if (refinement_cycle != 0)
 *             refine_mesh();
 *   
 *           const double target_tolerance = 1e-3 * std::pow(0.1, refinement_cycle);
 *           std::cout << "  Target_tolerance: " << target_tolerance << std::endl
 *                     << std::endl;
 *   
 * @endcode
 * 
 * This is where the fun starts. At the top we create the KINSOL solver
 * object and feed it with an object that encodes a number of additional
 * specifics (of which we only change the nonlinear tolerance we want to
 * reach; but you might want to look into what other members the
 * SUNDIALS::KINSOL::AdditionalData class has and play with them).
 * 
 * @code
 *           {
 *             typename SUNDIALS::KINSOL<Vector<double>>::AdditionalData
 *               additional_data;
 *             additional_data.function_tolerance = target_tolerance;
 *   
 *             SUNDIALS::KINSOL<Vector<double>> nonlinear_solver(additional_data);
 *   
 * @endcode
 * 
 * Then we have to describe the operations that were already mentioned
 * in the introduction. In essence, we have to teach KINSOL how to (i)
 * resize a vector to the correct size, (ii) compute the residual
 * vector, (iii) compute the Jacobian matrix (during which we also
 * compute its factorization), and (iv) solve a linear system with the
 * Jacobian.
 *           
 
 * 
 * All four of these operations are represented by member variables of
 * the SUNDIALS::KINSOL class that are of type `std::function`, i.e.,
 * they are objects to which we can assign a pointer to a function or,
 * as we do here, a "lambda function" that takes the appropriate
 * arguments and returns the appropriate information. It turns out
 * that we can do all of this in just over 20 lines of code.
 *           
 
 * 
 * (If you're not familiar what "lambda functions" are, take
 * a look at @ref step_12 "step-12" or at the
 * [wikipedia page](https://en.wikipedia.org/wiki/Anonymous_function)
 * on the subject. The idea of lambda functions is that one
 * wants to define a function with a certain set of
 * arguments, but (i) not make it a named functions because,
 * typically, the function is used in only one place and it
 * seems unnecessary to give it a global name; and (ii) that
 * the function has access to some of the variables that
 * exist at the place where it is defined, including member
 * variables. The syntax of lambda functions is awkward, but
 * ultimately quite useful.)
 *           
 
 * 
 * At the very end of the code block we then tell KINSOL to go to work
 * and solve our problem. The member functions called from the
 * 'residual', 'setup_jacobian', and 'solve_with_jacobian' functions
 * will then print output to screen that allows us to follow along
 * with the progress of the program.
 * 
 * @code
 *             nonlinear_solver.reinit_vector = [&](Vector<double> &x) {
 *               x.reinit(dof_handler.n_dofs());
 *             };
 *   
 *             nonlinear_solver.residual =
 *               [&](const Vector<double> &evaluation_point,
 *                   Vector<double> &      residual) {
 *                 compute_residual(evaluation_point, residual);
 *               };
 *   
 *             nonlinear_solver.setup_jacobian =
 *               [&](const Vector<double> &current_u,
 *                   const Vector<double> & /*current_f*/) {
 *                 compute_and_factorize_jacobian(current_u);
 *               };
 *   
 *             nonlinear_solver.solve_with_jacobian = [&](const Vector<double> &rhs,
 *                                                        Vector<double> &      dst,
 *                                                        const double tolerance) {
 *               solve(rhs, dst, tolerance);
 *             };
 *   
 *             nonlinear_solver.solve(current_solution);
 *           }
 *   
 * @endcode
 * 
 * The rest is then just house-keeping: Writing data to a file for
 * visualizing, and showing a summary of the timing collected so that we
 * can interpret how long each operation has taken, how often it was
 * executed, etc:
 * 
 * @code
 *           output_results(refinement_cycle);
 *   
 *           computing_timer.print_summary();
 *   
 *           std::cout << std::endl;
 *         }
 *     }
 *   } // namespace Step77
 *   
 *   
 *   int main()
 *   {
 *     try
 *       {
 *         using namespace Step77;
 *   
 *         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>
 
 
When running the program, you get output that looks like this:
@code
Mesh refinement step 0
  Target_tolerance: 0.001
 
  Computing residual vector... norm=0.867975
  Computing Jacobian matrix
  Factorizing Jacobian matrix
  Solving linear system
  Computing residual vector... norm=0.867975
  Computing residual vector... norm=0.212073
  Solving linear system
  Computing residual vector... norm=0.212073
  Computing residual vector... norm=0.202631
  Solving linear system
  Computing residual vector... norm=0.202631
  Computing residual vector... norm=0.165773
  Solving linear system
  Computing residual vector... norm=0.165774
  Computing residual vector... norm=0.162594
  Solving linear system
  Computing residual vector... norm=0.162594
  Computing residual vector... norm=0.148175
  Solving linear system
  Computing residual vector... norm=0.148175
  Computing residual vector... norm=0.145391
  Solving linear system
  Computing residual vector... norm=0.145391
  Computing residual vector... norm=0.137551
  Solving linear system
  Computing residual vector... norm=0.137551
  Computing residual vector... norm=0.135366
  Solving linear system
  Computing residual vector... norm=0.135365
  Computing residual vector... norm=0.130367
  Solving linear system
  Computing residual vector... norm=0.130367
  Computing residual vector... norm=0.128704
  Computing Jacobian matrix
  Factorizing Jacobian matrix
  Solving linear system
  Computing residual vector... norm=0.128704
  Computing residual vector... norm=0.0302623
  Solving linear system
  Computing residual vector... norm=0.0302624
  Computing residual vector... norm=0.0126764
  Solving linear system
  Computing residual vector... norm=0.0126763
  Computing residual vector... norm=0.00488315
  Solving linear system
  Computing residual vector... norm=0.00488322
  Computing residual vector... norm=0.00195788
  Solving linear system
  Computing residual vector... norm=0.00195781
  Computing residual vector... norm=0.000773169
 
 
+---------------------------------------------+------------+------------+
| Total wallclock time elapsed since start    |     0.121s |            |
|                                             |            |            |
| Section                         | no. calls |  wall time | % of total |
+---------------------------------+-----------+------------+------------+
| assembling the Jacobian         |         2 |    0.0151s |        12% |
| assembling the residual         |        31 |    0.0945s |        78% |
| factorizing the Jacobian        |         2 |   0.00176s |       1.5% |
| graphical output                |         1 |   0.00504s |       4.2% |
| linear system solve             |        15 |  0.000893s |      0.74% |
+---------------------------------+-----------+------------+------------+
 
 
Mesh refinement step 1
  Target_tolerance: 0.0001
 
  Computing residual vector... norm=0.2467
  Computing Jacobian matrix
  Factorizing Jacobian matrix
  Solving linear system
  Computing residual vector... norm=0.246699
  Computing residual vector... norm=0.0357783
  Solving linear system
  Computing residual vector... norm=0.0357784
  Computing residual vector... norm=0.0222161
  Solving linear system
  Computing residual vector... norm=0.022216
  Computing residual vector... norm=0.0122148
  Solving linear system
  Computing residual vector... norm=0.0122149
  Computing residual vector... norm=0.00750795
  Solving linear system
  Computing residual vector... norm=0.00750787
  Computing residual vector... norm=0.00439629
  Solving linear system
  Computing residual vector... norm=0.00439638
  Computing residual vector... norm=0.00265093
  Solving linear system
 
[...]
@endcode
 
The way this should be interpreted is most easily explained by looking at
the first few lines of the output on the first mesh:
@code
Mesh refinement step 0
Mesh refinement step 0
  Target_tolerance: 0.001
 
  Computing residual vector... norm=0.867975
  Computing Jacobian matrix
  Factorizing Jacobian matrix
  Solving linear system
  Computing residual vector... norm=0.867975
  Computing residual vector... norm=0.212073
  Solving linear system
  Computing residual vector... norm=0.212073
  Computing residual vector... norm=0.202631
  Solving linear system
  Computing residual vector... norm=0.202631
  Computing residual vector... norm=0.165773
  Solving linear system
  Computing residual vector... norm=0.165774
  Computing residual vector... norm=0.162594
  Solving linear system
  Computing residual vector... norm=0.162594
  Computing residual vector... norm=0.148175
  Solving linear system
  ...
@endcode
What is happening is this:
- In the first residual computation, KINSOL computes the residual to see whether
  the desired tolerance has been reached. The answer is no, so it requests the
  user program to compute the Jacobian matrix (and the function then also
  factorizes the matrix via SparseDirectUMFPACK).
- KINSOL then instructs us to solve a linear system of the form
  @f$J_k \, \delta U_k = -F_k@f$ with this matrix and the previously computed
  residual vector.
- It is then time to determine how far we want to go in this direction,
  i.e., do line search. To this end, KINSOL requires us to compute the
  residual vector @f$F(U_k + \alpha_k \delta U_k)@f$ for different step lengths
  @f$\alpha_k@f$. For the first step above, it finds an acceptable @f$\alpha_k@f$
  after two tries, and that's generally what will happen in later line
  searches as well.
- Having found a suitable updated solution @f$U_{k+1}@f$, the process is
  repeated except now KINSOL is happy with the current Jacobian matrix
  and does not instruct us to re-build the matrix and its factorization,
  instead asking us to solve a linear system with that same matrix. That
  will happen several times over, and only after ten solves with the same
  matrix are we instructed to build a matrix again, using what is by then an
  already substantially improved solution as linearization point.
 
The program also writes the solution to a VTU file at the end
of each mesh refinement cycle, and it looks as follows:
<table width="60%" align="center">
  <tr>
    <td align="center">
      <img src="https://www.dealii.org/images/steps/developer/step-77.solution.png" alt="">
    </td>
  </tr>
</table>
 
 
The key takeaway messages of this program are the following:
 
- The solution is the same as the one we computed in @ref step_15 "step-15", i.e., the
  interfaces to %SUNDIALS' KINSOL package really did what they were supposed
  to do. This should not come as a surprise, but the important point is that
  we don't have to spend the time implementing the complex algorithms that
  underlie advanced nonlinear solvers ourselves.
 
- KINSOL is able to avoid all sorts of operations such as rebuilding the
  Jacobian matrix when that is not actually necessary. Comparing the
  number of linear solves in the output above with the number of times
  we rebuild the Jacobian and compute its factorization should make it
  clear that this leads to very substantial savings in terms of compute
  times, without us having to implement the intricacies of algorithms
  that determine when we need to rebuild this information.
 
 
<a name="extensions"></a>
<a name="Possibilitiesforextensions"></a><h3> Possibilities for extensions </h3>
 
 
<a name="Betterlinearsolvers"></a><h4> Better linear solvers </h4>
 
 
For all but the small problems we consider here, a sparse direct solver
requires too much time and memory -- we need an iterative solver like
we use in many other programs. The trade-off between constructing an
expensive preconditioner (say, a geometric or algebraic multigrid method)
is different in the current case, however: Since we can re-use the same
matrix for numerous linear solves, we can do the same for the preconditioner
and putting more work into building a good preconditioner can more easily
be justified than if we used it only for a single linear solve as one
does for many other situations.
 
But iterative solvers also afford other opportunities. For example (and as
discussed briefly in the introduction), we may not need to solve to
very high accuracy (small tolerances) in early nonlinear iterations as long
as we are still far away from the actual solution. This was the basis of the
Eisenstat-Walker trick mentioned there. (This is also the underlying reason
why one can store the matrix in single precision rather than double precision,
see the discussion in the "Possibilities for extensions" section of @ref step_15 "step-15".)
 
KINSOL provides the function that does the linear solution with a target
tolerance that needs to be reached. We ignore it in the program above
because the direct solver we use does not need a tolerance and instead
solves the linear system exactly (up to round-off, of course), but iterative
solvers could make use of this kind of information -- and, in fact, should.
Indeed, the infrastructure is already there: The `solve()` function of this
program is declared as
@code
  template <int dim>
  void MinimalSurfaceProblem<dim>::solve(const Vector<double> &rhs,
                                         Vector<double> &      solution,
                                         const double /*tolerance*/)
@endcode
i.e., the `tolerance` parameter already exists, but is unused.
 
 
<a name="ReplacingSUNDIALSKINSOLbyPETScsSNES"></a><h4> Replacing SUNDIALS' KINSOL by PETSc's SNES </h4>
 
 
As mentioned in the introduction, SUNDIALS' KINSOL package is not the
only player in town. Rather, very similar interfaces exist to the SNES
package that is part of PETSc, and the NOX package that is part of
Trilinos, via the PETScWrappers::NonlinearSolver and
TrilinosWrappers::NOXSolver classes.
 
It is not very difficult to change the program to use either of these
two alternatives. Rather than show exactly what needs to be done,
let us point out that a version of this program that uses SNES instead
of KINSOL is available as part of the test suite, in the file
`tests/petsc/step-77-snes.cc`. Setting up the solver for
PETScWrappers::NonlinearSolver turns out to be even simpler than
for the SUNDIALS::KINSOL class we use here because we don't even
need the `reinit` lambda function -- SNES only needs us to set up
the remaining three functions `residual`, `setup_jacobian`, and
`solve_with_jacobian`. The majority of changes necessary to convert
the program to use SNES are related to the fact that SNES can only
deal with PETSc vectors and matrices, and these need to be set up
slightly differently. On the upside, the test suite program mentioned
above already works in parallel.
 
SNES also allows playing with a number of parameters about the
solver, and that enables some interesting comparisons between
methods. When you run the test program (or a slightly modified
version that outputs information to the screen instead of a file),
you get output that looks comparable to something like this:
@code
Mesh refinement step 0
  Target_tolerance: 0.001
 
  Computing residual vector
0 norm=0.867975
  Computing Jacobian matrix
  Computing residual vector
  Computing residual vector
1 norm=0.212073
  Computing Jacobian matrix
  Computing residual vector
  Computing residual vector
2 norm=0.0189603
  Computing Jacobian matrix
  Computing residual vector
  Computing residual vector
3 norm=0.000314854
 
[...]
@endcode
 
By default, PETSc uses a Newton solver with cubic backtracking,
resampling the Jacobian matrix at each Newton step.  That is, we
compute and factorize the matrix once per Newton step, and then sample
the residual to check for a successul line-search.
 
The attentive reader should have noticed that in this case we are
computing one more extra residual per Newton step.  This is because
the deal.II code is set up to use a Jacobian-free approach, and the
extra residual computation pops up when computing a matrix-vector
product to test the validity of the Newton solution.
 
PETSc can be configured in many interesting ways via the command line.
We can visualize the details of the solver by using the command line
argument **-snes_view**, which produces the excerpt below at the end
of each solve call:
@code
Mesh refinement step 0
[...]
SNES Object: 1 MPI process
  type: newtonls
  maximum iterations=50, maximum function evaluations=10000
  tolerances: relative=1e-08, absolute=0.001, solution=1e-08
  total number of linear solver iterations=3
  total number of function evaluations=7
  norm schedule ALWAYS
  Jacobian is applied matrix-free with differencing
  Jacobian is applied matrix-free with differencing, no explicit Jacobian
  SNESLineSearch Object: 1 MPI process
    type: bt
      interpolation: cubic
      alpha=1.000000e-04
    maxstep=1.000000e+08, minlambda=1.000000e-12
    tolerances: relative=1.000000e-08, absolute=1.000000e-15, lambda=1.000000e-08
    maximum iterations=40
  KSP Object: 1 MPI process
    type: preonly
    maximum iterations=10000, initial guess is zero
    tolerances:  relative=1e-05, absolute=1e-50, divergence=10000.
    left preconditioning
    using NONE norm type for convergence test
  PC Object: 1 MPI process
    type: shell
      deal.II user solve
    linear system matrix followed by preconditioner matrix:
    Mat Object: 1 MPI process
      type: mffd
      rows=89, cols=89
        Matrix-free approximation:
          err=1.49012e-08 (relative error in function evaluation)
          Using wp compute h routine
              Does not compute normU
    Mat Object: 1 MPI process
      type: seqaij
      rows=89, cols=89
      total: nonzeros=745, allocated nonzeros=745
      total number of mallocs used during MatSetValues calls=0
        not using I-node routines
[...]
@endcode
From the above details, we see that we are using the "newtonls" solver
type ("Newton line search"), with "bt" ("backtracting") line search.
 
From the output of **-snes_view** we can also get information about
the linear solver details; specifically, when using the
`solve_with_jacobian` interface, the deal.II interface internally uses
a custom solver configuration within a "shell" preconditioner, that
wraps the action of `solve_with_jacobian`.
 
We can also see the details of the type of matrices used within the
solve: "mffd" (matrix-free finite-differencing) for the action of the
linearized operator and "seqaij" for the assembled Jacobian we have
used to construct the preconditioner.
 
Diagnostics for the line search procedure can be turned on using the
command line **-snes_linesearch_monitor**, producing the excerpt
below:
@code
Mesh refinement step 0
  Target_tolerance: 0.001
 
  Computing residual vector
0 norm=0.867975
  Computing Jacobian matrix
  Computing residual vector
  Computing residual vector
      Line search: Using full step: fnorm 8.679748230595e-01 gnorm 2.120728179320e-01
1 norm=0.212073
  Computing Jacobian matrix
  Computing residual vector
  Computing residual vector
      Line search: Using full step: fnorm 2.120728179320e-01 gnorm 1.896033864659e-02
2 norm=0.0189603
  Computing Jacobian matrix
  Computing residual vector
  Computing residual vector
      Line search: Using full step: fnorm 1.896033864659e-02 gnorm 3.148542199408e-04
3 norm=0.000314854
 
[...]
@endcode
 
Within the run, the Jacobian matrix is assembled (and factored) 29 times:
@code
./step-77-snes | grep "Computing Jacobian" | wc -l
29
@endcode
 
KINSOL internally decided when it was necessary to update the Jacobian
matrix (which is when it would call `setup_jacobian`). SNES can do
something similar: We can compute the explicit sparse Jacobian matrix
only once per refinement step (and reuse the initial factorization) by
using the command line **-snes_lag_jacobian -2**, producing:
@code
./step-77-snes -snes_lag_jacobian -2 | grep "Computing Jacobian" | wc -l
6
@endcode
In other words, this dramatically reduces the number of times we have to
build the Jacobian matrix, though at a cost to the number of
nonlinear steps we have to take.
 
The lagging period can also be decided automatically. For example, if
we want to recompute the Jacobian at every other step:
@code
./step-77-snes -snes_lag_jacobian 2 | grep "Computing Jacobian" | wc -l
25
@endcode
Note, however, that we didn't exactly halve the number of Jacobian
computations. In this case the solution process will require many more
nonlinear iterations since the accuracy of the linear system solve is
not enough.
 
If we switch to using the preconditioned conjugate gradient method as
a linear solve, still using our initial factorization as
preconditioner, we get:
@code
./step-77-snes -snes_lag_jacobian 2 -ksp_type cg | grep "Computing Jacobian" | wc -l
17
@endcode
Note that in this case we use an approximate preconditioner (the LU
factorization of the initial approximation) but we use a matrix-free
operator for the action of our Jacobian matrix, thus solving for the
correct linear system.
 
We can switch to a quasi-Newton method by using the command
line **-snes_type qn -snes_qn_scale_type jacobian**, and we can see that
our Jacobian is sampled and factored only when needed, at the cost of
an increase of the number of steps:
@code
Mesh refinement step 0
  Target_tolerance: 0.001
 
  Computing residual vector
0 norm=0.867975
  Computing Jacobian matrix
  Computing residual vector
  Computing residual vector
1 norm=0.166391
  Computing residual vector
  Computing residual vector
2 norm=0.0507703
  Computing residual vector
  Computing residual vector
3 norm=0.0160007
  Computing residual vector
  Computing residual vector
  Computing residual vector
4 norm=0.00172425
  Computing residual vector
  Computing residual vector
  Computing residual vector
5 norm=0.000460486
[...]
@endcode
 
<a href="https://www.mcs.anl.gov/papers/P2010-0112.pdf">Nonlinear preconditioning</a>
can also be used. For example, we can run a right-preconditioned nonlinear
GMRES, using one Newton step as a preconditioner, with the command:
@code
./step-77-snes -snes_type ngmres -npc_snes_type newtonls -snes_monitor -npc_snes_monitor | grep SNES
  0 SNES Function norm 8.679748230595e-01
    0 SNES Function norm 8.679748230595e-01
    1 SNES Function norm 2.120738413585e-01
  1 SNES Function norm 1.284613424341e-01
    0 SNES Function norm 1.284613424341e-01
    1 SNES Function norm 6.539358995036e-03
  2 SNES Function norm 5.148828618635e-03
    0 SNES Function norm 5.148828618635e-03
    1 SNES Function norm 6.048613313899e-06
  3 SNES Function norm 3.199913594705e-06
  0 SNES Function norm 2.464793634583e-01
    0 SNES Function norm 2.464793634583e-01
    1 SNES Function norm 3.591625291931e-02
  1 SNES Function norm 3.235827289342e-02
    0 SNES Function norm 3.235827289342e-02
    1 SNES Function norm 1.249214136060e-03
  2 SNES Function norm 5.302288687547e-04
    0 SNES Function norm 5.302288687547e-04
    1 SNES Function norm 1.490247730530e-07
  3 SNES Function norm 1.436531309822e-07
  0 SNES Function norm 5.044203686086e-01
    0 SNES Function norm 5.044203686086e-01
    1 SNES Function norm 1.716855756535e-01
  1 SNES Function norm 7.770484434662e-02
    0 SNES Function norm 7.770484434662e-02
    1 SNES Function norm 2.462422395554e-02
  2 SNES Function norm 1.438187947066e-02
    0 SNES Function norm 1.438187947066e-02
    1 SNES Function norm 9.214168343848e-04
  3 SNES Function norm 2.268378169625e-04
    0 SNES Function norm 2.268378169625e-04
    1 SNES Function norm 3.463704776158e-07
  4 SNES Function norm 9.964533647277e-08
  0 SNES Function norm 1.942213246154e-01
    0 SNES Function norm 1.942213246154e-01
    1 SNES Function norm 1.125558372384e-01
  1 SNES Function norm 1.309880643103e-01
    0 SNES Function norm 1.309880643103e-01
    1 SNES Function norm 2.595634741967e-02
  2 SNES Function norm 1.149616419685e-02
    0 SNES Function norm 1.149616419685e-02
    1 SNES Function norm 7.204904831783e-04
  3 SNES Function norm 6.743539224973e-04
    0 SNES Function norm 6.743539224973e-04
    1 SNES Function norm 1.521290969181e-05
  4 SNES Function norm 8.121151857453e-06
    0 SNES Function norm 8.121151857453e-06
    1 SNES Function norm 1.460470903719e-09
  5 SNES Function norm 9.982794797188e-10
  0 SNES Function norm 1.225979459424e-01
    0 SNES Function norm 1.225979459424e-01
    1 SNES Function norm 4.946412992249e-02
  1 SNES Function norm 2.466574163571e-02
    0 SNES Function norm 2.466574163571e-02
    1 SNES Function norm 8.537739703503e-03
  2 SNES Function norm 5.935412895618e-03
    0 SNES Function norm 5.935412895618e-03
    1 SNES Function norm 3.699307476482e-04
  3 SNES Function norm 2.188768476656e-04
    0 SNES Function norm 2.188768476656e-04
    1 SNES Function norm 9.478344390128e-07
  4 SNES Function norm 4.559224590570e-07
    0 SNES Function norm 4.559224590570e-07
    1 SNES Function norm 1.317127376721e-11
  5 SNES Function norm 1.311046524394e-11
  0 SNES Function norm 1.011637873732e-01
    0 SNES Function norm 1.011637873732e-01
    1 SNES Function norm 1.072720108836e-02
  1 SNES Function norm 8.985302820531e-03
    0 SNES Function norm 8.985302820531e-03
    1 SNES Function norm 5.807781788861e-04
  2 SNES Function norm 5.594756759727e-04
    0 SNES Function norm 5.594756759727e-04
    1 SNES Function norm 1.834638371641e-05
  3 SNES Function norm 1.408280767367e-05
    0 SNES Function norm 1.408280767367e-05
    1 SNES Function norm 5.763656314185e-08
  4 SNES Function norm 1.702747382189e-08
    0 SNES Function norm 1.702747382189e-08
    1 SNES Function norm 1.452722802538e-12
  5 SNES Function norm 1.444478767837e-12
@endcode
 
 
As also discussed for the KINSOL use above, optimal preconditioners
should be used instead of the LU factorization used here by
default. This is already possible within this tutorial by playing with
the command line options.  For example, algebraic multigrid can be
used by simply specifying **-pc_type gamg**.  When using iterative
linear solvers, the "Eisenstat-Walker trick" @cite eiwa96 can be also
requested at command line via **-snes_ksp_ew**.  Using these options,
we can see that the number of nonlinear iterations used by the solver
increases as the mesh is refined, and that the number of linear
iterations increases as the Newton solver is entering the second-order
ball of convergence:
@code
./step-77-snes -pc_type gamg -ksp_type cg -ksp_converged_reason -snes_converged_reason -snes_ksp_ew | grep CONVERGED
  Linear solve converged due to CONVERGED_RTOL iterations 1
  Linear solve converged due to CONVERGED_RTOL iterations 2
  Linear solve converged due to CONVERGED_RTOL iterations 3
Nonlinear solve converged due to CONVERGED_FNORM_ABS iterations 3
  Linear solve converged due to CONVERGED_RTOL iterations 1
  Linear solve converged due to CONVERGED_RTOL iterations 1
  Linear solve converged due to CONVERGED_RTOL iterations 2
Nonlinear solve converged due to CONVERGED_FNORM_ABS iterations 3
  Linear solve converged due to CONVERGED_RTOL iterations 1
  Linear solve converged due to CONVERGED_RTOL iterations 2
  Linear solve converged due to CONVERGED_RTOL iterations 2
  Linear solve converged due to CONVERGED_RTOL iterations 2
  Linear solve converged due to CONVERGED_RTOL iterations 3
  Linear solve converged due to CONVERGED_RTOL iterations 4
Nonlinear solve converged due to CONVERGED_FNORM_ABS iterations 6
  Linear solve converged due to CONVERGED_RTOL iterations 1
  Linear solve converged due to CONVERGED_RTOL iterations 1
  Linear solve converged due to CONVERGED_RTOL iterations 1
  Linear solve converged due to CONVERGED_RTOL iterations 1
  Linear solve converged due to CONVERGED_RTOL iterations 1
  Linear solve converged due to CONVERGED_RTOL iterations 1
  Linear solve converged due to CONVERGED_RTOL iterations 1
  Linear solve converged due to CONVERGED_RTOL iterations 1
  Linear solve converged due to CONVERGED_RTOL iterations 1
  Linear solve converged due to CONVERGED_RTOL iterations 2
  Linear solve converged due to CONVERGED_RTOL iterations 4
  Linear solve converged due to CONVERGED_RTOL iterations 7
Nonlinear solve converged due to CONVERGED_FNORM_ABS iterations 12
  Linear solve converged due to CONVERGED_RTOL iterations 1
  Linear solve converged due to CONVERGED_RTOL iterations 2
  Linear solve converged due to CONVERGED_RTOL iterations 3
  Linear solve converged due to CONVERGED_RTOL iterations 4
  Linear solve converged due to CONVERGED_RTOL iterations 7
Nonlinear solve converged due to CONVERGED_FNORM_ABS iterations 5
  Linear solve converged due to CONVERGED_RTOL iterations 2
  Linear solve converged due to CONVERGED_RTOL iterations 3
  Linear solve converged due to CONVERGED_RTOL iterations 7
  Linear solve converged due to CONVERGED_RTOL iterations 6
  Linear solve converged due to CONVERGED_RTOL iterations 7
  Linear solve converged due to CONVERGED_RTOL iterations 12
Nonlinear solve converged due to CONVERGED_FNORM_ABS iterations 6
@endcode
 
Finally we describe how to get some diagnostic on the correctness of
the computed Jacobian.  Deriving the correct linearization is
sometimes difficult: It took a page or two in the introduction to
derive the exact bilinear form for the Jacobian matrix, and it would
be quite nice compute it automatically from the residual of which it
is the derivative.  (This is what @ref step_72 "step-72" does!) But if one is set on
doing things by hand, it would at least be nice if we had a way to
check the correctness of the derivation. SNES allows us to do this: we
can use the options **-snes_test_jacobian -snes_test_jacobian_view**:
@code
Mesh refinement step 0
  Target_tolerance: 0.001
 
  Computing residual vector
0 norm=0.867975
  Computing Jacobian matrix
  ---------- Testing Jacobian -------------
  Testing hand-coded Jacobian, if (for double precision runs) ||J - Jfd||_F/||J||_F is
    O(1.e-8), the hand-coded Jacobian is probably correct.
[...]
  ||J - Jfd||_F/||J||_F = 0.0196815, ||J - Jfd||_F = 0.503436
[...]
  Hand-coded minus finite-difference Jacobian with tolerance 1e-05 ----------
Mat Object: 1 MPI process
  type: seqaij
row 0: (0, 0.125859)
row 1: (1, 0.0437112)
row 2:
row 3:
row 4: (4, 0.902232)
row 5:
row 6:
row 7:
row 8:
row 9: (9, 0.537306)
row 10:
row 11: (11, 1.38157)
row 12:
[...]
@endcode
showing that the only errors we commit in assembling the Jacobian are
on the boundary dofs.  As discussed in the tutorial, those errors are
harmless.
 
The key take-away messages of this modification of the tutorial program are
therefore basically the same of what we already found using KINSOL:
 
- The solution is the same as the one we computed in @ref step_15 "step-15", i.e., the
  interfaces to PETSc SNES package really did what they were supposed
  to do. This should not come as a surprise, but the important point is that
  we don't have to spend the time implementing the complex algorithms that
  underlie advanced nonlinear solvers ourselves.
 
- SNES offers a wide variety of solvers and line search techniques,
  not only Newton.  It also allows us to control Jacobian setups;
  however, differently from KINSOL, this is not automatically decided
  within the library by looking at the residual vector but it needs to
  be specified by the user.
 
 
 
 
<a name="ReplacingSUNDIALSKINSOLbyTrilinosNOXpackage"></a><h4> Replacing SUNDIALS' KINSOL by Trilinos' NOX package </h4>
 
 
Besides KINSOL and SNES, the third option you have is to use the NOX
package. As before, rather than showing in detail how that needs to
happen, let us simply point out that the test suite program
`tests/trilinos/step-77-with-nox.cc` does this. The modifications
necessary to use NOX instead of KINSOL are quite minimal; in
particular, NOX (unlike SNES) is happy to work with deal.II's own
vector and matrix classes.
 
 
<a name="ReplacingSUNDIALSKINSOLbyagenericnonlinearsolver"></a><h4> Replacing SUNDIALS' KINSOL by a generic nonlinear solver </h4>
 
 
Having to choose which of these three frameworks (KINSOL, SNES, or NOX)
to use at compile time is cumbersome when wanting to compare things. It
would be nicer if one could decide the package to use at run time, assuming that one
has a copy of deal.II installed that is compiled against all three of these
dependencies. It turns out that this is possible, using the class
NonlinearSolverSelector that presents a common interface to all three of
these solvers, along with the ability to choose which one to use based
on run-time parameters.
 *
 *
<a name="PlainProg"></a>
<h1> The plain program</h1>
@include "step-77.cc"
*/