Reference documentation for deal.II version 9.5.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
step-77.h
Go to the documentation of this file.
1) const
467 *   {
468 *   return std::sin(2 * numbers::PI * (p[0] + p[1]));
469 *   }
470 *  
471 *  
472 * @endcode
473 *
474 *
475 * <a name="ThecodeMinimalSurfaceProblemcodeclassimplementation"></a>
476 * <h3>The <code>MinimalSurfaceProblem</code> class implementation</h3>
477 *
478
479 *
480 *
481 * <a name="Constructorandsetupfunctions"></a>
482 * <h4>Constructor and set up functions</h4>
483 *
484
485 *
486 * The following few functions are also essentially copies of what
487 * @ref step_15 "step-15" already does, and so there is little to discuss.
488 *
489 * @code
490 *   template <int dim>
491 *   MinimalSurfaceProblem<dim>::MinimalSurfaceProblem()
492 *   : dof_handler(triangulation)
493 *   , fe(1)
494 *   , computing_timer(std::cout, TimerOutput::never, TimerOutput::wall_times)
495 *   {}
496 *  
497 *  
498 *  
499 *   template <int dim>
500 *   void MinimalSurfaceProblem<dim>::setup_system(const bool initial_step)
501 *   {
502 *   TimerOutput::Scope t(computing_timer, "set up");
503 *  
504 *   if (initial_step)
505 *   {
506 *   dof_handler.distribute_dofs(fe);
507 *   current_solution.reinit(dof_handler.n_dofs());
508 *  
509 *   hanging_node_constraints.clear();
511 *   hanging_node_constraints);
512 *   hanging_node_constraints.close();
513 *   }
514 *  
515 *   DynamicSparsityPattern dsp(dof_handler.n_dofs());
516 *   DoFTools::make_sparsity_pattern(dof_handler, dsp);
517 *  
518 *   hanging_node_constraints.condense(dsp);
519 *  
520 *   sparsity_pattern.copy_from(dsp);
521 *   jacobian_matrix.reinit(sparsity_pattern);
522 *   jacobian_matrix_factorization.reset();
523 *   }
524 *  
525 *  
526 *  
527 * @endcode
528 *
529 *
530 * <a name="AssemblingandfactorizingtheJacobianmatrix"></a>
531 * <h4>Assembling and factorizing the Jacobian matrix</h4>
532 *
533
534 *
535 * The following function is then responsible for assembling and factorizing
536 * the Jacobian matrix. The first half of the function is in essence the
537 * `assemble_system()` function of @ref step_15 "step-15", except that it does not deal with
538 * also forming a right hand side vector (i.e., the residual) since we do not
539 * always have to do these operations at the same time.
540 *
541
542 *
543 * We put the whole assembly functionality into a code block enclosed by curly
544 * braces so that we can use a TimerOutput::Scope variable to measure how much
545 * time is spent in this code block, excluding everything that happens in this
546 * function after the matching closing brace `}`.
547 *
548 * @code
549 *   template <int dim>
550 *   void MinimalSurfaceProblem<dim>::compute_and_factorize_jacobian(
551 *   const Vector<double> &evaluation_point)
552 *   {
553 *   {
554 *   TimerOutput::Scope t(computing_timer, "assembling the Jacobian");
555 *  
556 *   std::cout << " Computing Jacobian matrix" << std::endl;
557 *  
558 *   const QGauss<dim> quadrature_formula(fe.degree + 1);
559 *  
560 *   jacobian_matrix = 0;
561 *  
562 *   FEValues<dim> fe_values(fe,
563 *   quadrature_formula,
565 *   update_JxW_values);
566 *  
567 *   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
568 *   const unsigned int n_q_points = quadrature_formula.size();
569 *  
570 *   FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
571 *  
572 *   std::vector<Tensor<1, dim>> evaluation_point_gradients(n_q_points);
573 *  
574 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
575 *  
576 *   for (const auto &cell : dof_handler.active_cell_iterators())
577 *   {
578 *   cell_matrix = 0;
579 *  
580 *   fe_values.reinit(cell);
581 *  
582 *   fe_values.get_function_gradients(evaluation_point,
583 *   evaluation_point_gradients);
584 *  
585 *   for (unsigned int q = 0; q < n_q_points; ++q)
586 *   {
587 *   const double coeff =
588 *   1.0 / std::sqrt(1 + evaluation_point_gradients[q] *
589 *   evaluation_point_gradients[q]);
590 *  
591 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
592 *   {
593 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
594 *   cell_matrix(i, j) +=
595 *   (((fe_values.shape_grad(i, q) // ((\nabla \phi_i
596 *   * coeff // * a_n
597 *   * fe_values.shape_grad(j, q)) // * \nabla \phi_j)
598 *   - // -
599 *   (fe_values.shape_grad(i, q) // (\nabla \phi_i
600 *   * coeff * coeff * coeff // * a_n^3
601 *   *
602 *   (fe_values.shape_grad(j, q) // * (\nabla \phi_j
603 *   * evaluation_point_gradients[q]) // * \nabla u_n)
604 *   * evaluation_point_gradients[q])) // * \nabla u_n)))
605 *   * fe_values.JxW(q)); // * dx
606 *   }
607 *   }
608 *  
609 *   cell->get_dof_indices(local_dof_indices);
610 *   hanging_node_constraints.distribute_local_to_global(cell_matrix,
611 *   local_dof_indices,
612 *   jacobian_matrix);
613 *   }
614 *  
615 *   std::map<types::global_dof_index, double> boundary_values;
617 *   0,
619 *   boundary_values);
620 *   Vector<double> dummy_solution(dof_handler.n_dofs());
621 *   Vector<double> dummy_rhs(dof_handler.n_dofs());
622 *   MatrixTools::apply_boundary_values(boundary_values,
623 *   jacobian_matrix,
624 *   dummy_solution,
625 *   dummy_rhs);
626 *   }
627 *  
628 * @endcode
629 *
630 * The second half of the function then deals with factorizing the
631 * so-computed matrix. To do this, we first create a new SparseDirectUMFPACK
632 * object and by assigning it to the member variable
633 * `jacobian_matrix_factorization`, we also destroy whatever object that
634 * pointer previously pointed to (if any). Then we tell the object to
635 * factorize the Jacobian.
636 *
637
638 *
639 * As above, we enclose this block of code into curly braces and use a timer
640 * to assess how long this part of the program takes.
641 *
642
643 *
644 * (Strictly speaking, we don't actually need the matrix any more after we
645 * are done here, and could throw the matrix object away. A code intended to
646 * be memory efficient would do this, and only create the matrix object in
647 * this function, rather than as a member variable of the surrounding class.
648 * We omit this step here because using the same coding style as in previous
649 * tutorial programs breeds familiarity with the common style and helps make
650 * these tutorial programs easier to read.)
651 *
652 * @code
653 *   {
654 *   TimerOutput::Scope t(computing_timer, "factorizing the Jacobian");
655 *  
656 *   std::cout << " Factorizing Jacobian matrix" << std::endl;
657 *  
658 *   jacobian_matrix_factorization = std::make_unique<SparseDirectUMFPACK>();
659 *   jacobian_matrix_factorization->factorize(jacobian_matrix);
660 *   }
661 *   }
662 *  
663 *  
664 *  
665 * @endcode
666 *
667 *
668 * <a name="Computingtheresidualvector"></a>
669 * <h4>Computing the residual vector</h4>
670 *
671
672 *
673 * The second part of what `assemble_system()` used to do in @ref step_15 "step-15" is
674 * computing the residual vector, i.e., the right hand side vector of the
675 * Newton linear systems. We have broken this out of the previous function,
676 * but the following function will be easy to understand if you understood
677 * what `assemble_system()` in @ref step_15 "step-15" did. Importantly, however, we need to
678 * compute the residual not linearized around the current solution vector, but
679 * whatever we get from KINSOL. This is necessary for operations such as line
680 * search where we want to know what the residual @f$F(U^k + \alpha_k \delta
681 * U^K)@f$ is for different values of @f$\alpha_k@f$; KINSOL in those cases simply
682 * gives us the argument to the function @f$F@f$ and we then compute the residual
683 * @f$F(\cdot)@f$ at this point.
684 *
685
686 *
687 * The function prints the norm of the so-computed residual at the end as a
688 * way for us to follow along the progress of the program.
689 *
690 * @code
691 *   template <int dim>
692 *   void MinimalSurfaceProblem<dim>::compute_residual(
693 *   const Vector<double> &evaluation_point,
694 *   Vector<double> & residual)
695 *   {
696 *   TimerOutput::Scope t(computing_timer, "assembling the residual");
697 *  
698 *   std::cout << " Computing residual vector..." << std::flush;
699 *  
700 *   residual = 0.0;
701 *  
702 *   const QGauss<dim> quadrature_formula(fe.degree + 1);
703 *   FEValues<dim> fe_values(fe,
704 *   quadrature_formula,
705 *   update_gradients | update_quadrature_points |
706 *   update_JxW_values);
707 *  
708 *   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
709 *   const unsigned int n_q_points = quadrature_formula.size();
710 *  
711 *   Vector<double> cell_residual(dofs_per_cell);
712 *   std::vector<Tensor<1, dim>> evaluation_point_gradients(n_q_points);
713 *  
714 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
715 *  
716 *   for (const auto &cell : dof_handler.active_cell_iterators())
717 *   {
718 *   cell_residual = 0;
719 *   fe_values.reinit(cell);
720 *  
721 *   fe_values.get_function_gradients(evaluation_point,
722 *   evaluation_point_gradients);
723 *  
724 *  
725 *   for (unsigned int q = 0; q < n_q_points; ++q)
726 *   {
727 *   const double coeff =
728 *   1.0 / std::sqrt(1 + evaluation_point_gradients[q] *
729 *   evaluation_point_gradients[q]);
730 *  
731 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
732 *   cell_residual(i) +=
733 *   (fe_values.shape_grad(i, q) // \nabla \phi_i
734 *   * coeff // * a_n
735 *   * evaluation_point_gradients[q] // * \nabla u_n
736 *   * fe_values.JxW(q)); // * dx
737 *   }
738 *  
739 *   cell->get_dof_indices(local_dof_indices);
740 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
741 *   residual(local_dof_indices[i]) += cell_residual(i);
742 *   }
743 *  
744 *   hanging_node_constraints.condense(residual);
745 *  
746 *   for (const types::global_dof_index i :
747 *   DoFTools::extract_boundary_dofs(dof_handler))
748 *   residual(i) = 0;
749 *  
750 *   for (const types::global_dof_index i :
751 *   DoFTools::extract_hanging_node_dofs(dof_handler))
752 *   residual(i) = 0;
753 *  
754 *   std::cout << " norm=" << residual.l2_norm() << std::endl;
755 *   }
756 *  
757 *  
758 *  
759 * @endcode
760 *
761 *
762 * <a name="SolvinglinearsystemswiththeJacobianmatrix"></a>
763 * <h4>Solving linear systems with the Jacobian matrix</h4>
764 *
765
766 *
767 * Next up is the function that implements the solution of a linear system
768 * with the Jacobian matrix. Since we have already factored the matrix when we
769 * built the matrix, solving a linear system comes down to applying the
770 * inverse matrix to the given right hand side vector: This is what the
771 * SparseDirectUMFPACK::vmult() function does that we use here. Following
772 * this, we have to make sure that we also address the values of hanging nodes
773 * in the solution vector, and this is done using
774 * AffineConstraints::distribute().
775 *
776
777 *
778 * The function takes an additional, but unused, argument `tolerance` that
779 * indicates how accurately we have to solve the linear system. The meaning of
780 * this argument is discussed in the introduction in the context of the
781 * "Eisenstat Walker trick", but since we are using a direct rather than an
782 * iterative solver, we are not using this opportunity to solve linear systems
783 * only inexactly.
784 *
785 * @code
786 *   template <int dim>
787 *   void MinimalSurfaceProblem<dim>::solve(const Vector<double> &rhs,
788 *   Vector<double> & solution,
789 *   const double /*tolerance*/)
790 *   {
791 *   TimerOutput::Scope t(computing_timer, "linear system solve");
792 *  
793 *   std::cout << " Solving linear system" << std::endl;
794 *  
795 *   jacobian_matrix_factorization->vmult(solution, rhs);
796 *  
797 *   hanging_node_constraints.distribute(solution);
798 *   }
799 *  
800 *  
801 *  
802 * @endcode
803 *
804 *
805 * <a name="Refiningthemeshsettingboundaryvaluesandgeneratinggraphicaloutput"></a>
806 * <h4>Refining the mesh, setting boundary values, and generating graphical output</h4>
807 *
808
809 *
810 * The following three functions are again simply copies of the ones in
811 * @ref step_15 "step-15":
812 *
813 * @code
814 *   template <int dim>
815 *   void MinimalSurfaceProblem<dim>::refine_mesh()
816 *   {
817 *   Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
818 *  
819 *   KellyErrorEstimator<dim>::estimate(
820 *   dof_handler,
821 *   QGauss<dim - 1>(fe.degree + 1),
822 *   std::map<types::boundary_id, const Function<dim> *>(),
823 *   current_solution,
824 *   estimated_error_per_cell);
825 *  
826 *   GridRefinement::refine_and_coarsen_fixed_number(triangulation,
827 *   estimated_error_per_cell,
828 *   0.3,
829 *   0.03);
830 *  
831 *   triangulation.prepare_coarsening_and_refinement();
832 *  
833 *   SolutionTransfer<dim> solution_transfer(dof_handler);
834 *   solution_transfer.prepare_for_coarsening_and_refinement(current_solution);
835 *  
836 *   triangulation.execute_coarsening_and_refinement();
837 *  
838 *   dof_handler.distribute_dofs(fe);
839 *  
840 *   Vector<double> tmp(dof_handler.n_dofs());
841 *   solution_transfer.interpolate(current_solution, tmp);
842 *   current_solution = std::move(tmp);
843 *  
844 *   hanging_node_constraints.clear();
845 *  
846 *   DoFTools::make_hanging_node_constraints(dof_handler,
847 *   hanging_node_constraints);
848 *   hanging_node_constraints.close();
849 *  
850 *   hanging_node_constraints.distribute(current_solution);
851 *  
852 *   set_boundary_values();
853 *  
854 *   setup_system(/*initial_step=*/false);
855 *   }
856 *  
857 *  
858 *  
859 *   template <int dim>
860 *   void MinimalSurfaceProblem<dim>::set_boundary_values()
861 *   {
862 *   std::map<types::global_dof_index, double> boundary_values;
863 *   VectorTools::interpolate_boundary_values(dof_handler,
864 *   0,
865 *   BoundaryValues<dim>(),
866 *   boundary_values);
867 *   for (const auto &boundary_value : boundary_values)
868 *   current_solution(boundary_value.first) = boundary_value.second;
869 *  
870 *   hanging_node_constraints.distribute(current_solution);
871 *   }
872 *  
873 *  
874 *  
875 *   template <int dim>
876 *   void MinimalSurfaceProblem<dim>::output_results(
877 *   const unsigned int refinement_cycle)
878 *   {
879 *   TimerOutput::Scope t(computing_timer, "graphical output");
880 *  
881 *   DataOut<dim> data_out;
882 *  
883 *   data_out.attach_dof_handler(dof_handler);
884 *   data_out.add_data_vector(current_solution, "solution");
885 *   data_out.build_patches();
886 *  
887 *   const std::string filename =
888 *   "solution-" + Utilities::int_to_string(refinement_cycle, 2) + ".vtu";
889 *   std::ofstream output(filename);
890 *   data_out.write_vtu(output);
891 *   }
892 *  
893 *  
894 *  
895 * @endcode
896 *
897 *
898 * <a name="Therunfunctionandtheoveralllogicoftheprogram"></a>
899 * <h4>The run() function and the overall logic of the program</h4>
900 *
901
902 *
903 * The only function that *really* is interesting in this program is the one
904 * that drives the overall algorithm of starting on a coarse mesh, doing some
905 * mesh refinement cycles, and on each mesh using KINSOL to find the solution
906 * of the nonlinear algebraic equation we obtain from discretization on this
907 * mesh. The `refine_mesh()` function above makes sure that the solution on
908 * one mesh is used as the starting guess on the next mesh. We also use a
909 * TimerOutput object to measure how much time every operation on each mesh
910 * costs, and reset the timer at the beginning of each cycle.
911 *
912
913 *
914 * As discussed in the introduction, it is not necessary to solve problems on
915 * coarse meshes particularly accurately since these will only solve as
916 * starting guesses for the next mesh. As a consequence, we will use a target
917 * tolerance of
918 * @f$\tau=10^{-3} \frac{1}{10^k}@f$ for the @f$k@f$th mesh refinement cycle.
919 *
920
921 *
922 * All of this is encoded in the first part of this function:
923 *
924 * @code
925 *   template <int dim>
926 *   void MinimalSurfaceProblem<dim>::run()
927 *   {
928 *   GridGenerator::hyper_ball(triangulation);
929 *   triangulation.refine_global(2);
930 *  
931 *   setup_system(/*initial_step=*/true);
932 *   set_boundary_values();
933 *  
934 *   for (unsigned int refinement_cycle = 0; refinement_cycle < 6;
935 *   ++refinement_cycle)
936 *   {
937 *   computing_timer.reset();
938 *   std::cout << "Mesh refinement step " << refinement_cycle << std::endl;
939 *  
940 *   if (refinement_cycle != 0)
941 *   refine_mesh();
942 *  
943 *   const double target_tolerance = 1e-3 * std::pow(0.1, refinement_cycle);
944 *   std::cout << " Target_tolerance: " << target_tolerance << std::endl
945 *   << std::endl;
946 *  
947 * @endcode
948 *
949 * This is where the fun starts. At the top we create the KINSOL solver
950 * object and feed it with an object that encodes a number of additional
951 * specifics (of which we only change the nonlinear tolerance we want to
952 * reach; but you might want to look into what other members the
953 * SUNDIALS::KINSOL::AdditionalData class has and play with them).
954 *
955 * @code
956 *   {
957 *   typename SUNDIALS::KINSOL<Vector<double>>::AdditionalData
958 *   additional_data;
959 *   additional_data.function_tolerance = target_tolerance;
960 *  
961 *   SUNDIALS::KINSOL<Vector<double>> nonlinear_solver(additional_data);
962 *  
963 * @endcode
964 *
965 * Then we have to describe the operations that were already mentioned
966 * in the introduction. In essence, we have to teach KINSOL how to (i)
967 * resize a vector to the correct size, (ii) compute the residual
968 * vector, (iii) compute the Jacobian matrix (during which we also
969 * compute its factorization), and (iv) solve a linear system with the
970 * Jacobian.
971 *
972
973 *
974 * All four of these operations are represented by member variables of
975 * the SUNDIALS::KINSOL class that are of type `std::function`, i.e.,
976 * they are objects to which we can assign a pointer to a function or,
977 * as we do here, a "lambda function" that takes the appropriate
978 * arguments and returns the appropriate information. It turns out
979 * that we can do all of this in just over 20 lines of code.
980 *
981
982 *
983 * (If you're not familiar what "lambda functions" are, take
984 * a look at @ref step_12 "step-12" or at the
985 * [wikipedia page](https://en.wikipedia.org/wiki/Anonymous_function)
986 * on the subject. The idea of lambda functions is that one
987 * wants to define a function with a certain set of
988 * arguments, but (i) not make it a named functions because,
989 * typically, the function is used in only one place and it
990 * seems unnecessary to give it a global name; and (ii) that
991 * the function has access to some of the variables that
992 * exist at the place where it is defined, including member
993 * variables. The syntax of lambda functions is awkward, but
994 * ultimately quite useful.)
995 *
996
997 *
998 * At the very end of the code block we then tell KINSOL to go to work
999 * and solve our problem. The member functions called from the
1000 * 'residual', 'setup_jacobian', and 'solve_with_jacobian' functions
1001 * will then print output to screen that allows us to follow along
1002 * with the progress of the program.
1003 *
1004 * @code
1005 *   nonlinear_solver.reinit_vector = [&](Vector<double> &x) {
1006 *   x.reinit(dof_handler.n_dofs());
1007 *   };
1008 *  
1009 *   nonlinear_solver.residual =
1010 *   [&](const Vector<double> &evaluation_point,
1011 *   Vector<double> & residual) {
1012 *   compute_residual(evaluation_point, residual);
1013 *   };
1014 *  
1015 *   nonlinear_solver.setup_jacobian =
1016 *   [&](const Vector<double> &current_u,
1017 *   const Vector<double> & /*current_f*/) {
1018 *   compute_and_factorize_jacobian(current_u);
1019 *   };
1020 *  
1021 *   nonlinear_solver.solve_with_jacobian = [&](const Vector<double> &rhs,
1022 *   Vector<double> & dst,
1023 *   const double tolerance) {
1024 *   solve(rhs, dst, tolerance);
1025 *   };
1026 *  
1027 *   nonlinear_solver.solve(current_solution);
1028 *   }
1029 *  
1030 * @endcode
1031 *
1032 * The rest is then just house-keeping: Writing data to a file for
1033 * visualizing, and showing a summary of the timing collected so that we
1034 * can interpret how long each operation has taken, how often it was
1035 * executed, etc:
1036 *
1037 * @code
1038 *   output_results(refinement_cycle);
1039 *  
1040 *   computing_timer.print_summary();
1041 *  
1042 *   std::cout << std::endl;
1043 *   }
1044 *   }
1045 *   } // namespace Step77
1046 *  
1047 *  
1048 *   int main()
1049 *   {
1050 *   try
1051 *   {
1052 *   using namespace Step77;
1053 *  
1054 *   MinimalSurfaceProblem<2> laplace_problem_2d;
1055 *   laplace_problem_2d.run();
1056 *   }
1057 *   catch (std::exception &exc)
1058 *   {
1059 *   std::cerr << std::endl
1060 *   << std::endl
1061 *   << "----------------------------------------------------"
1062 *   << std::endl;
1063 *   std::cerr << "Exception on processing: " << std::endl
1064 *   << exc.what() << std::endl
1065 *   << "Aborting!" << std::endl
1066 *   << "----------------------------------------------------"
1067 *   << std::endl;
1068 *  
1069 *   return 1;
1070 *   }
1071 *   catch (...)
1072 *   {
1073 *   std::cerr << std::endl
1074 *   << std::endl
1075 *   << "----------------------------------------------------"
1076 *   << std::endl;
1077 *   std::cerr << "Unknown exception!" << std::endl
1078 *   << "Aborting!" << std::endl
1079 *   << "----------------------------------------------------"
1080 *   << std::endl;
1081 *   return 1;
1082 *   }
1083 *   return 0;
1084 *   }
1085 * @endcode
1086<a name="Results"></a><h1>Results</h1>
1087
1088
1089When running the program, you get output that looks like this:
1090@code
1091Mesh refinement step 0
1092 Target_tolerance: 0.001
1093
1094 Computing residual vector... norm=0.867975
1095 Computing Jacobian matrix
1096 Factorizing Jacobian matrix
1097 Solving linear system
1098 Computing residual vector... norm=0.867975
1099 Computing residual vector... norm=0.212073
1100 Solving linear system
1101 Computing residual vector... norm=0.212073
1102 Computing residual vector... norm=0.202631
1103 Solving linear system
1104 Computing residual vector... norm=0.202631
1105 Computing residual vector... norm=0.165773
1106 Solving linear system
1107 Computing residual vector... norm=0.165774
1108 Computing residual vector... norm=0.162594
1109 Solving linear system
1110 Computing residual vector... norm=0.162594
1111 Computing residual vector... norm=0.148175
1112 Solving linear system
1113 Computing residual vector... norm=0.148175
1114 Computing residual vector... norm=0.145391
1115 Solving linear system
1116 Computing residual vector... norm=0.145391
1117 Computing residual vector... norm=0.137551
1118 Solving linear system
1119 Computing residual vector... norm=0.137551
1120 Computing residual vector... norm=0.135366
1121 Solving linear system
1122 Computing residual vector... norm=0.135365
1123 Computing residual vector... norm=0.130367
1124 Solving linear system
1125 Computing residual vector... norm=0.130367
1126 Computing residual vector... norm=0.128704
1127 Computing Jacobian matrix
1128 Factorizing Jacobian matrix
1129 Solving linear system
1130 Computing residual vector... norm=0.128704
1131 Computing residual vector... norm=0.0302623
1132 Solving linear system
1133 Computing residual vector... norm=0.0302624
1134 Computing residual vector... norm=0.0126764
1135 Solving linear system
1136 Computing residual vector... norm=0.0126763
1137 Computing residual vector... norm=0.00488315
1138 Solving linear system
1139 Computing residual vector... norm=0.00488322
1140 Computing residual vector... norm=0.00195788
1141 Solving linear system
1142 Computing residual vector... norm=0.00195781
1143 Computing residual vector... norm=0.000773169
1144
1145
1146+---------------------------------------------+------------+------------+
1147| Total wallclock time elapsed since start | 0.121s | |
1148| | | |
1149| Section | no. calls | wall time | % of total |
1150+---------------------------------+-----------+------------+------------+
1151| assembling the Jacobian | 2 | 0.0151s | 12% |
1152| assembling the residual | 31 | 0.0945s | 78% |
1153| factorizing the Jacobian | 2 | 0.00176s | 1.5% |
1154| graphical output | 1 | 0.00504s | 4.2% |
1155| linear system solve | 15 | 0.000893s | 0.74% |
1156+---------------------------------+-----------+------------+------------+
1157
1158
1159Mesh refinement step 1
1160 Target_tolerance: 0.0001
1161
1162 Computing residual vector... norm=0.2467
1163 Computing Jacobian matrix
1164 Factorizing Jacobian matrix
1165 Solving linear system
1166 Computing residual vector... norm=0.246699
1167 Computing residual vector... norm=0.0357783
1168 Solving linear system
1169 Computing residual vector... norm=0.0357784
1170 Computing residual vector... norm=0.0222161
1171 Solving linear system
1172 Computing residual vector... norm=0.022216
1173 Computing residual vector... norm=0.0122148
1174 Solving linear system
1175 Computing residual vector... norm=0.0122149
1176 Computing residual vector... norm=0.00750795
1177 Solving linear system
1178 Computing residual vector... norm=0.00750787
1179 Computing residual vector... norm=0.00439629
1180 Solving linear system
1181 Computing residual vector... norm=0.00439638
1182 Computing residual vector... norm=0.00265093
1183 Solving linear system
1184
1185[...]
1186@endcode
1187
1188The way this should be interpreted is most easily explained by looking at
1189the first few lines of the output on the first mesh:
1190@code
1191Mesh refinement step 0
1192Mesh refinement step 0
1193 Target_tolerance: 0.001
1194
1195 Computing residual vector... norm=0.867975
1196 Computing Jacobian matrix
1197 Factorizing Jacobian matrix
1198 Solving linear system
1199 Computing residual vector... norm=0.867975
1200 Computing residual vector... norm=0.212073
1201 Solving linear system
1202 Computing residual vector... norm=0.212073
1203 Computing residual vector... norm=0.202631
1204 Solving linear system
1205 Computing residual vector... norm=0.202631
1206 Computing residual vector... norm=0.165773
1207 Solving linear system
1208 Computing residual vector... norm=0.165774
1209 Computing residual vector... norm=0.162594
1210 Solving linear system
1211 Computing residual vector... norm=0.162594
1212 Computing residual vector... norm=0.148175
1213 Solving linear system
1214 ...
1215@endcode
1216What is happening is this:
1217- In the first residual computation, KINSOL computes the residual to see whether
1218 the desired tolerance has been reached. The answer is no, so it requests the
1219 user program to compute the Jacobian matrix (and the function then also
1220 factorizes the matrix via SparseDirectUMFPACK).
1221- KINSOL then instructs us to solve a linear system of the form
1222 @f$J_k \, \delta U_k = -F_k@f$ with this matrix and the previously computed
1223 residual vector.
1224- It is then time to determine how far we want to go in this direction,
1225 i.e., do line search. To this end, KINSOL requires us to compute the
1226 residual vector @f$F(U_k + \alpha_k \delta U_k)@f$ for different step lengths
1227 @f$\alpha_k@f$. For the first step above, it finds an acceptable @f$\alpha_k@f$
1228 after two tries, and that's generally what will happen in later line
1229 searches as well.
1230- Having found a suitable updated solution @f$U_{k+1}@f$, the process is
1231 repeated except now KINSOL is happy with the current Jacobian matrix
1232 and does not instruct us to re-build the matrix and its factorization,
1233 instead asking us to solve a linear system with that same matrix. That
1234 will happen several times over, and only after ten solves with the same
1235 matrix are we instructed to build a matrix again, using what is by then an
1236 already substantially improved solution as linearization point.
1237
1238The program also writes the solution to a VTU file at the end
1239of each mesh refinement cycle, and it looks as follows:
1240<table width="60%" align="center">
1241 <tr>
1242 <td align="center">
1243 <img src="https://www.dealii.org/images/steps/developer/step-77.solution.png" alt="">
1244 </td>
1245 </tr>
1246</table>
1247
1248
1249The key takeaway messages of this program are the following:
1250
1251- The solution is the same as the one we computed in @ref step_15 "step-15", i.e., the
1252 interfaces to %SUNDIALS' KINSOL package really did what they were supposed
1253 to do. This should not come as a surprise, but the important point is that
1254 we don't have to spend the time implementing the complex algorithms that
1255 underlie advanced nonlinear solvers ourselves.
1256
1257- KINSOL is able to avoid all sorts of operations such as rebuilding the
1258 Jacobian matrix when that is not actually necessary. Comparing the
1259 number of linear solves in the output above with the number of times
1260 we rebuild the Jacobian and compute its factorization should make it
1261 clear that this leads to very substantial savings in terms of compute
1262 times, without us having to implement the intricacies of algorithms
1263 that determine when we need to rebuild this information.
1264
1265
1266<a name="extensions"></a>
1267<a name="Possibilitiesforextensions"></a><h3> Possibilities for extensions </h3>
1268
1269
1270<a name="Betterlinearsolvers"></a><h4> Better linear solvers </h4>
1271
1272
1273For all but the small problems we consider here, a sparse direct solver
1274requires too much time and memory -- we need an iterative solver like
1275we use in many other programs. The trade-off between constructing an
1276expensive preconditioner (say, a geometric or algebraic multigrid method)
1277is different in the current case, however: Since we can re-use the same
1278matrix for numerous linear solves, we can do the same for the preconditioner
1279and putting more work into building a good preconditioner can more easily
1280be justified than if we used it only for a single linear solve as one
1281does for many other situations.
1282
1283But iterative solvers also afford other opportunities. For example (and as
1284discussed briefly in the introduction), we may not need to solve to
1285very high accuracy (small tolerances) in early nonlinear iterations as long
1286as we are still far away from the actual solution. This was the basis of the
1287Eisenstat-Walker trick mentioned there. (This is also the underlying reason
1288why one can store the matrix in single precision rather than double precision,
1289see the discussion in the "Possibilities for extensions" section of @ref step_15 "step-15".)
1290
1291KINSOL provides the function that does the linear solution with a target
1292tolerance that needs to be reached. We ignore it in the program above
1293because the direct solver we use does not need a tolerance and instead
1294solves the linear system exactly (up to round-off, of course), but iterative
1295solvers could make use of this kind of information -- and, in fact, should.
1296Indeed, the infrastructure is already there: The `solve()` function of this
1297program is declared as
1298@code
1299 template <int dim>
1300 void MinimalSurfaceProblem<dim>::solve(const Vector<double> &rhs,
1301 Vector<double> & solution,
1302 const double /*tolerance*/)
1303@endcode
1304i.e., the `tolerance` parameter already exists, but is unused.
1305
1306
1307<a name="ReplacingSUNDIALSKINSOLbyPETScsSNES"></a><h4> Replacing SUNDIALS' KINSOL by PETSc's SNES </h4>
1308
1309
1310As mentioned in the introduction, SUNDIALS' KINSOL package is not the
1311only player in town. Rather, very similar interfaces exist to the SNES
1312package that is part of PETSc, and the NOX package that is part of
1313Trilinos, via the PETScWrappers::NonlinearSolver and
1314TrilinosWrappers::NOXSolver classes.
1315
1316It is not very difficult to change the program to use either of these
1317two alternatives. Rather than show exactly what needs to be done,
1318let us point out that a version of this program that uses SNES instead
1319of KINSOL is available as part of the test suite, in the file
1320`tests/petsc/step-77-snes.cc`. Setting up the solver for
1321PETScWrappers::NonlinearSolver turns out to be even simpler than
1322for the SUNDIALS::KINSOL class we use here because we don't even
1323need the `reinit` lambda function -- SNES only needs us to set up
1324the remaining three functions `residual`, `setup_jacobian`, and
1325`solve_with_jacobian`. The majority of changes necessary to convert
1326the program to use SNES are related to the fact that SNES can only
1327deal with PETSc vectors and matrices, and these need to be set up
1328slightly differently. On the upside, the test suite program mentioned
1329above already works in parallel.
1330
1331SNES also allows playing with a number of parameters about the
1332solver, and that enables some interesting comparisons between
1333methods. When you run the test program (or a slightly modified
1334version that outputs information to the screen instead of a file),
1335you get output that looks comparable to something like this:
1336@code
1337Mesh refinement step 0
1338 Target_tolerance: 0.001
1339
1340 Computing residual vector
13410 norm=0.867975
1342 Computing Jacobian matrix
1343 Computing residual vector
1344 Computing residual vector
13451 norm=0.212073
1346 Computing Jacobian matrix
1347 Computing residual vector
1348 Computing residual vector
13492 norm=0.0189603
1350 Computing Jacobian matrix
1351 Computing residual vector
1352 Computing residual vector
13533 norm=0.000314854
1354
1355[...]
1356@endcode
1357
1358By default, PETSc uses a Newton solver with cubic backtracking,
1359resampling the Jacobian matrix at each Newton step. That is, we
1360compute and factorize the matrix once per Newton step, and then sample
1361the residual to check for a successul line-search.
1362
1363The attentive reader should have noticed that in this case we are
1364computing one more extra residual per Newton step. This is because
1365the deal.II code is set up to use a Jacobian-free approach, and the
1366extra residual computation pops up when computing a matrix-vector
1367product to test the validity of the Newton solution.
1368
1369PETSc can be configured in many interesting ways via the command line.
1370We can visualize the details of the solver by using the command line
1371argument **-snes_view**, which produces the excerpt below at the end
1372of each solve call:
1373@code
1374Mesh refinement step 0
1375[...]
1376SNES Object: 1 MPI process
1377 type: newtonls
1378 maximum iterations=50, maximum function evaluations=10000
1379 tolerances: relative=1e-08, absolute=0.001, solution=1e-08
1380 total number of linear solver iterations=3
1381 total number of function evaluations=7
1382 norm schedule ALWAYS
1383 Jacobian is applied matrix-free with differencing
1384 Jacobian is applied matrix-free with differencing, no explicit Jacobian
1385 SNESLineSearch Object: 1 MPI process
1386 type: bt
1387 interpolation: cubic
1388 alpha=1.000000e-04
1389 maxstep=1.000000e+08, minlambda=1.000000e-12
1390 tolerances: relative=1.000000e-08, absolute=1.000000e-15, lambda=1.000000e-08
1391 maximum iterations=40
1392 KSP Object: 1 MPI process
1393 type: preonly
1394 maximum iterations=10000, initial guess is zero
1395 tolerances: relative=1e-05, absolute=1e-50, divergence=10000.
1396 left preconditioning
1397 using NONE norm type for convergence test
1398 PC Object: 1 MPI process
1399 type: shell
1400 deal.II user solve
1401 linear system matrix followed by preconditioner matrix:
1402 Mat Object: 1 MPI process
1403 type: mffd
1404 rows=89, cols=89
1405 Matrix-free approximation:
1406 err=1.49012e-08 (relative error in function evaluation)
1407 Using wp compute h routine
1408 Does not compute normU
1409 Mat Object: 1 MPI process
1410 type: seqaij
1411 rows=89, cols=89
1412 total: nonzeros=745, allocated nonzeros=745
1413 total number of mallocs used during MatSetValues calls=0
1414 not using I-node routines
1415[...]
1416@endcode
1417From the above details, we see that we are using the "newtonls" solver
1418type ("Newton line search"), with "bt" ("backtracting") line search.
1419
1420From the output of **-snes_view** we can also get information about
1421the linear solver details; specifically, when using the
1422`solve_with_jacobian` interface, the deal.II interface internally uses
1423a custom solver configuration within a "shell" preconditioner, that
1424wraps the action of `solve_with_jacobian`.
1425
1426We can also see the details of the type of matrices used within the
1427solve: "mffd" (matrix-free finite-differencing) for the action of the
1428linearized operator and "seqaij" for the assembled Jacobian we have
1429used to construct the preconditioner.
1430
1431Diagnostics for the line search procedure can be turned on using the
1432command line **-snes_linesearch_monitor**, producing the excerpt
1433below:
1434@code
1435Mesh refinement step 0
1436 Target_tolerance: 0.001
1437
1438 Computing residual vector
14390 norm=0.867975
1440 Computing Jacobian matrix
1441 Computing residual vector
1442 Computing residual vector
1443 Line search: Using full step: fnorm 8.679748230595e-01 gnorm 2.120728179320e-01
14441 norm=0.212073
1445 Computing Jacobian matrix
1446 Computing residual vector
1447 Computing residual vector
1448 Line search: Using full step: fnorm 2.120728179320e-01 gnorm 1.896033864659e-02
14492 norm=0.0189603
1450 Computing Jacobian matrix
1451 Computing residual vector
1452 Computing residual vector
1453 Line search: Using full step: fnorm 1.896033864659e-02 gnorm 3.148542199408e-04
14543 norm=0.000314854
1455
1456[...]
1457@endcode
1458
1459Within the run, the Jacobian matrix is assembled (and factored) 29 times:
1460@code
1461./step-77-snes | grep "Computing Jacobian" | wc -l
146229
1463@endcode
1464
1465KINSOL internally decided when it was necessary to update the Jacobian
1466matrix (which is when it would call `setup_jacobian`). SNES can do
1467something similar: We can compute the explicit sparse Jacobian matrix
1468only once per refinement step (and reuse the initial factorization) by
1469using the command line **-snes_lag_jacobian -2**, producing:
1470@code
1471./step-77-snes -snes_lag_jacobian -2 | grep "Computing Jacobian" | wc -l
14726
1473@endcode
1474In other words, this dramatically reduces the number of times we have to
1475build the Jacobian matrix, though at a cost to the number of
1476nonlinear steps we have to take.
1477
1478The lagging period can also be decided automatically. For example, if
1479we want to recompute the Jacobian at every other step:
1480@code
1481./step-77-snes -snes_lag_jacobian 2 | grep "Computing Jacobian" | wc -l
148225
1483@endcode
1484Note, however, that we didn't exactly halve the number of Jacobian
1485computations. In this case the solution process will require many more
1486nonlinear iterations since the accuracy of the linear system solve is
1487not enough.
1488
1489If we switch to using the preconditioned conjugate gradient method as
1490a linear solve, still using our initial factorization as
1491preconditioner, we get:
1492@code
1493./step-77-snes -snes_lag_jacobian 2 -ksp_type cg | grep "Computing Jacobian" | wc -l
149417
1495@endcode
1496Note that in this case we use an approximate preconditioner (the LU
1497factorization of the initial approximation) but we use a matrix-free
1498operator for the action of our Jacobian matrix, thus solving for the
1499correct linear system.
1500
1501We can switch to a quasi-Newton method by using the command
1502line **-snes_type qn -snes_qn_scale_type jacobian**, and we can see that
1503our Jacobian is sampled and factored only when needed, at the cost of
1504an increase of the number of steps:
1505@code
1506Mesh refinement step 0
1507 Target_tolerance: 0.001
1508
1509 Computing residual vector
15100 norm=0.867975
1511 Computing Jacobian matrix
1512 Computing residual vector
1513 Computing residual vector
15141 norm=0.166391
1515 Computing residual vector
1516 Computing residual vector
15172 norm=0.0507703
1518 Computing residual vector
1519 Computing residual vector
15203 norm=0.0160007
1521 Computing residual vector
1522 Computing residual vector
1523 Computing residual vector
15244 norm=0.00172425
1525 Computing residual vector
1526 Computing residual vector
1527 Computing residual vector
15285 norm=0.000460486
1529[...]
1530@endcode
1531
1532<a href="https://www.mcs.anl.gov/papers/P2010-0112.pdf">Nonlinear preconditioning</a>
1533can also be used. For example, we can run a right-preconditioned nonlinear
1534GMRES, using one Newton step as a preconditioner, with the command:
1535@code
1536./step-77-snes -snes_type ngmres -npc_snes_type newtonls -snes_monitor -npc_snes_monitor | grep SNES
1537 0 SNES Function norm 8.679748230595e-01
1538 0 SNES Function norm 8.679748230595e-01
1539 1 SNES Function norm 2.120738413585e-01
1540 1 SNES Function norm 1.284613424341e-01
1541 0 SNES Function norm 1.284613424341e-01
1542 1 SNES Function norm 6.539358995036e-03
1543 2 SNES Function norm 5.148828618635e-03
1544 0 SNES Function norm 5.148828618635e-03
1545 1 SNES Function norm 6.048613313899e-06
1546 3 SNES Function norm 3.199913594705e-06
1547 0 SNES Function norm 2.464793634583e-01
1548 0 SNES Function norm 2.464793634583e-01
1549 1 SNES Function norm 3.591625291931e-02
1550 1 SNES Function norm 3.235827289342e-02
1551 0 SNES Function norm 3.235827289342e-02
1552 1 SNES Function norm 1.249214136060e-03
1553 2 SNES Function norm 5.302288687547e-04
1554 0 SNES Function norm 5.302288687547e-04
1555 1 SNES Function norm 1.490247730530e-07
1556 3 SNES Function norm 1.436531309822e-07
1557 0 SNES Function norm 5.044203686086e-01
1558 0 SNES Function norm 5.044203686086e-01
1559 1 SNES Function norm 1.716855756535e-01
1560 1 SNES Function norm 7.770484434662e-02
1561 0 SNES Function norm 7.770484434662e-02
1562 1 SNES Function norm 2.462422395554e-02
1563 2 SNES Function norm 1.438187947066e-02
1564 0 SNES Function norm 1.438187947066e-02
1565 1 SNES Function norm 9.214168343848e-04
1566 3 SNES Function norm 2.268378169625e-04
1567 0 SNES Function norm 2.268378169625e-04
1568 1 SNES Function norm 3.463704776158e-07
1569 4 SNES Function norm 9.964533647277e-08
1570 0 SNES Function norm 1.942213246154e-01
1571 0 SNES Function norm 1.942213246154e-01
1572 1 SNES Function norm 1.125558372384e-01
1573 1 SNES Function norm 1.309880643103e-01
1574 0 SNES Function norm 1.309880643103e-01
1575 1 SNES Function norm 2.595634741967e-02
1576 2 SNES Function norm 1.149616419685e-02
1577 0 SNES Function norm 1.149616419685e-02
1578 1 SNES Function norm 7.204904831783e-04
1579 3 SNES Function norm 6.743539224973e-04
1580 0 SNES Function norm 6.743539224973e-04
1581 1 SNES Function norm 1.521290969181e-05
1582 4 SNES Function norm 8.121151857453e-06
1583 0 SNES Function norm 8.121151857453e-06
1584 1 SNES Function norm 1.460470903719e-09
1585 5 SNES Function norm 9.982794797188e-10
1586 0 SNES Function norm 1.225979459424e-01
1587 0 SNES Function norm 1.225979459424e-01
1588 1 SNES Function norm 4.946412992249e-02
1589 1 SNES Function norm 2.466574163571e-02
1590 0 SNES Function norm 2.466574163571e-02
1591 1 SNES Function norm 8.537739703503e-03
1592 2 SNES Function norm 5.935412895618e-03
1593 0 SNES Function norm 5.935412895618e-03
1594 1 SNES Function norm 3.699307476482e-04
1595 3 SNES Function norm 2.188768476656e-04
1596 0 SNES Function norm 2.188768476656e-04
1597 1 SNES Function norm 9.478344390128e-07
1598 4 SNES Function norm 4.559224590570e-07
1599 0 SNES Function norm 4.559224590570e-07
1600 1 SNES Function norm 1.317127376721e-11
1601 5 SNES Function norm 1.311046524394e-11
1602 0 SNES Function norm 1.011637873732e-01
1603 0 SNES Function norm 1.011637873732e-01
1604 1 SNES Function norm 1.072720108836e-02
1605 1 SNES Function norm 8.985302820531e-03
1606 0 SNES Function norm 8.985302820531e-03
1607 1 SNES Function norm 5.807781788861e-04
1608 2 SNES Function norm 5.594756759727e-04
1609 0 SNES Function norm 5.594756759727e-04
1610 1 SNES Function norm 1.834638371641e-05
1611 3 SNES Function norm 1.408280767367e-05
1612 0 SNES Function norm 1.408280767367e-05
1613 1 SNES Function norm 5.763656314185e-08
1614 4 SNES Function norm 1.702747382189e-08
1615 0 SNES Function norm 1.702747382189e-08
1616 1 SNES Function norm 1.452722802538e-12
1617 5 SNES Function norm 1.444478767837e-12
1618@endcode
1619
1620
1621As also discussed for the KINSOL use above, optimal preconditioners
1622should be used instead of the LU factorization used here by
1623default. This is already possible within this tutorial by playing with
1624the command line options. For example, algebraic multigrid can be
1625used by simply specifying **-pc_type gamg**. When using iterative
1626linear solvers, the "Eisenstat-Walker trick" @cite eiwa96 can be also
1627requested at command line via **-snes_ksp_ew**. Using these options,
1628we can see that the number of nonlinear iterations used by the solver
1629increases as the mesh is refined, and that the number of linear
1630iterations increases as the Newton solver is entering the second-order
1631ball of convergence:
1632@code
1633./step-77-snes -pc_type gamg -ksp_type cg -ksp_converged_reason -snes_converged_reason -snes_ksp_ew | grep CONVERGED
1634 Linear solve converged due to CONVERGED_RTOL iterations 1
1635 Linear solve converged due to CONVERGED_RTOL iterations 2
1636 Linear solve converged due to CONVERGED_RTOL iterations 3
1637Nonlinear solve converged due to CONVERGED_FNORM_ABS iterations 3
1638 Linear solve converged due to CONVERGED_RTOL iterations 1
1639 Linear solve converged due to CONVERGED_RTOL iterations 1
1640 Linear solve converged due to CONVERGED_RTOL iterations 2
1641Nonlinear solve converged due to CONVERGED_FNORM_ABS iterations 3
1642 Linear solve converged due to CONVERGED_RTOL iterations 1
1643 Linear solve converged due to CONVERGED_RTOL iterations 2
1644 Linear solve converged due to CONVERGED_RTOL iterations 2
1645 Linear solve converged due to CONVERGED_RTOL iterations 2
1646 Linear solve converged due to CONVERGED_RTOL iterations 3
1647 Linear solve converged due to CONVERGED_RTOL iterations 4
1648Nonlinear solve converged due to CONVERGED_FNORM_ABS iterations 6
1649 Linear solve converged due to CONVERGED_RTOL iterations 1
1650 Linear solve converged due to CONVERGED_RTOL iterations 1
1651 Linear solve converged due to CONVERGED_RTOL iterations 1
1652 Linear solve converged due to CONVERGED_RTOL iterations 1
1653 Linear solve converged due to CONVERGED_RTOL iterations 1
1654 Linear solve converged due to CONVERGED_RTOL iterations 1
1655 Linear solve converged due to CONVERGED_RTOL iterations 1
1656 Linear solve converged due to CONVERGED_RTOL iterations 1
1657 Linear solve converged due to CONVERGED_RTOL iterations 1
1658 Linear solve converged due to CONVERGED_RTOL iterations 2
1659 Linear solve converged due to CONVERGED_RTOL iterations 4
1660 Linear solve converged due to CONVERGED_RTOL iterations 7
1661Nonlinear solve converged due to CONVERGED_FNORM_ABS iterations 12
1662 Linear solve converged due to CONVERGED_RTOL iterations 1
1663 Linear solve converged due to CONVERGED_RTOL iterations 2
1664 Linear solve converged due to CONVERGED_RTOL iterations 3
1665 Linear solve converged due to CONVERGED_RTOL iterations 4
1666 Linear solve converged due to CONVERGED_RTOL iterations 7
1667Nonlinear solve converged due to CONVERGED_FNORM_ABS iterations 5
1668 Linear solve converged due to CONVERGED_RTOL iterations 2
1669 Linear solve converged due to CONVERGED_RTOL iterations 3
1670 Linear solve converged due to CONVERGED_RTOL iterations 7
1671 Linear solve converged due to CONVERGED_RTOL iterations 6
1672 Linear solve converged due to CONVERGED_RTOL iterations 7
1673 Linear solve converged due to CONVERGED_RTOL iterations 12
1674Nonlinear solve converged due to CONVERGED_FNORM_ABS iterations 6
1675@endcode
1676
1677Finally we describe how to get some diagnostic on the correctness of
1678the computed Jacobian. Deriving the correct linearization is
1679sometimes difficult: It took a page or two in the introduction to
1680derive the exact bilinear form for the Jacobian matrix, and it would
1681be quite nice compute it automatically from the residual of which it
1682is the derivative. (This is what @ref step_72 "step-72" does!) But if one is set on
1683doing things by hand, it would at least be nice if we had a way to
1684check the correctness of the derivation. SNES allows us to do this: we
1685can use the options **-snes_test_jacobian -snes_test_jacobian_view**:
1686@code
1687Mesh refinement step 0
1688 Target_tolerance: 0.001
1689
1690 Computing residual vector
16910 norm=0.867975
1692 Computing Jacobian matrix
1693 ---------- Testing Jacobian -------------
1694 Testing hand-coded Jacobian, if (for double precision runs) ||J - Jfd||_F/||J||_F is
1695 O(1.e-8), the hand-coded Jacobian is probably correct.
1696[...]
1697 ||J - Jfd||_F/||J||_F = 0.0196815, ||J - Jfd||_F = 0.503436
1698[...]
1699 Hand-coded minus finite-difference Jacobian with tolerance 1e-05 ----------
1700Mat Object: 1 MPI process
1701 type: seqaij
1702row 0: (0, 0.125859)
1703row 1: (1, 0.0437112)
1704row 2:
1705row 3:
1706row 4: (4, 0.902232)
1707row 5:
1708row 6:
1709row 7:
1710row 8:
1711row 9: (9, 0.537306)
1712row 10:
1713row 11: (11, 1.38157)
1714row 12:
1715[...]
1716@endcode
1717showing that the only errors we commit in assembling the Jacobian are
1718on the boundary dofs. As discussed in the tutorial, those errors are
1719harmless.
1720
1721The key take-away messages of this modification of the tutorial program are
1722therefore basically the same of what we already found using KINSOL:
1723
1724- The solution is the same as the one we computed in @ref step_15 "step-15", i.e., the
1725 interfaces to PETSc SNES package really did what they were supposed
1726 to do. This should not come as a surprise, but the important point is that
1727 we don't have to spend the time implementing the complex algorithms that
1728 underlie advanced nonlinear solvers ourselves.
1729
1730- SNES offers a wide variety of solvers and line search techniques,
1731 not only Newton. It also allows us to control Jacobian setups;
1732 however, differently from KINSOL, this is not automatically decided
1733 within the library by looking at the residual vector but it needs to
1734 be specified by the user.
1735
1736
1737
1738
1739<a name="ReplacingSUNDIALSKINSOLbyTrilinosNOXpackage"></a><h4> Replacing SUNDIALS' KINSOL by Trilinos' NOX package </h4>
1740
1741
1742Besides KINSOL and SNES, the third option you have is to use the NOX
1743package. As before, rather than showing in detail how that needs to
1744happen, let us simply point out that the test suite program
1745`tests/trilinos/step-77-with-nox.cc` does this. The modifications
1746necessary to use NOX instead of KINSOL are quite minimal; in
1747particular, NOX (unlike SNES) is happy to work with deal.II's own
1748vector and matrix classes.
1749
1750
1751<a name="ReplacingSUNDIALSKINSOLbyagenericnonlinearsolver"></a><h4> Replacing SUNDIALS' KINSOL by a generic nonlinear solver </h4>
1752
1753
1754Having to choose which of these three frameworks (KINSOL, SNES, or NOX)
1755to use at compile time is cumbersome when wanting to compare things. It
1756would be nicer if one could decide the package to use at run time, assuming that one
1757has a copy of deal.II installed that is compiled against all three of these
1758dependencies. It turns out that this is possible, using the class
1759NonlinearSolverSelector that presents a common interface to all three of
1760these solvers, along with the ability to choose which one to use based
1761on run-time parameters.
1762 *
1763 *
1764<a name="PlainProg"></a>
1765<h1> The plain program</h1>
1766@include "step-77.cc"
1767*/
@ wall_times
Definition timer.h:652
Point< 2 > second
Definition grid_out.cc:4616
Point< 2 > first
Definition grid_out.cc:4615
__global__ void set(Number *val, const Number s, const size_type N)
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
@ matrix
Contents is actually a matrix.
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
Definition advection.h:75
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Definition divergence.h:472
void apply_boundary_values(const std::map< types::global_dof_index, number > &boundary_values, SparseMatrix< number > &matrix, Vector< number > &solution, Vector< number > &right_hand_side, const bool eliminate_columns=true)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
constexpr const ReferenceCell Line
VectorType::value_type * end(VectorType &V)
void interpolate_boundary_values(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const std::map< types::boundary_id, const Function< spacedim, number > * > &function_map, std::map< types::global_dof_index, number > &boundary_values, const ComponentMask &component_mask=ComponentMask())
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
int(&) functions(const void *v1, const void *v2)
static constexpr double PI
Definition numbers.h:259
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation