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