Loading [MathJax]/extensions/TeX/newcommand.js
 Reference documentation for deal.II version 9.3.3
\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\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
step-77.h
Go to the documentation of this file.
1) const
442 * {
443 * return std::sin(2 * numbers::PI * (p[0] + p[1]));
444 * }
445 *
446 *
447 * @endcode
448 *
449 *
450 * <a name="ThecodeMinimalSurfaceProblemcodeclassimplementation"></a>
451 * <h3>The <code>MinimalSurfaceProblem</code> class implementation</h3>
452 *
453
454 *
455 *
456 * <a name="Constructorandsetupfunctions"></a>
457 * <h4>Constructor and set up functions</h4>
458 *
459
460 *
461 * The following few functions are also essentially copies of what
462 * @ref step_15 "step-15" already does, and so there is little to discuss.
463 *
464 * @code
465 * template <int dim>
466 * MinimalSurfaceProblem<dim>::MinimalSurfaceProblem()
467 * : dof_handler(triangulation)
468 * , fe(1)
469 * , computing_timer(std::cout, TimerOutput::never, TimerOutput::wall_times)
470 * {}
471 *
472 *
473 *
474 * template <int dim>
475 * void MinimalSurfaceProblem<dim>::setup_system(const bool initial_step)
476 * {
477 * TimerOutput::Scope t(computing_timer, "set up");
478 *
479 * if (initial_step)
480 * {
481 * dof_handler.distribute_dofs(fe);
482 * current_solution.reinit(dof_handler.n_dofs());
483 *
484 * hanging_node_constraints.clear();
486 * hanging_node_constraints);
487 * hanging_node_constraints.close();
488 * }
489 *
490 * DynamicSparsityPattern dsp(dof_handler.n_dofs());
491 * DoFTools::make_sparsity_pattern(dof_handler, dsp);
492 *
493 * hanging_node_constraints.condense(dsp);
494 *
495 * sparsity_pattern.copy_from(dsp);
496 * jacobian_matrix.reinit(sparsity_pattern);
497 * jacobian_matrix_factorization.reset();
498 * }
499 *
500 *
501 *
502 * @endcode
503 *
504 *
505 * <a name="AssemblingandfactorizingtheJacobianmatrix"></a>
506 * <h4>Assembling and factorizing the Jacobian matrix</h4>
507 *
508
509 *
510 * The following function is then responsible for assembling and factorizing
511 * the Jacobian matrix. The first half of the function is in essence the
512 * `assemble_system()` function of @ref step_15 "step-15", except that it does not deal with
513 * also forming a right hand side vector (i.e., the residual) since we do not
514 * always have to do these operations at the same time.
515 *
516
517 *
518 * We put the whole assembly functionality into a code block enclosed by curly
519 * braces so that we can use a TimerOutput::Scope variable to measure how much
520 * time is spent in this code block, excluding everything that happens in this
521 * function after the matching closing brace `}`.
522 *
523 * @code
524 * template <int dim>
525 * void MinimalSurfaceProblem<dim>::compute_and_factorize_jacobian(
526 * const Vector<double> &evaluation_point)
527 * {
528 * {
529 * TimerOutput::Scope t(computing_timer, "assembling the Jacobian");
530 *
531 * std::cout << " Computing Jacobian matrix" << std::endl;
532 *
533 * const QGauss<dim> quadrature_formula(fe.degree + 1);
534 *
535 * jacobian_matrix = 0;
536 *
537 * FEValues<dim> fe_values(fe,
538 * quadrature_formula,
541 *
542 * const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
543 * const unsigned int n_q_points = quadrature_formula.size();
544 *
545 * FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
546 *
547 * std::vector<Tensor<1, dim>> evaluation_point_gradients(n_q_points);
548 *
549 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
550 *
551 * for (const auto &cell : dof_handler.active_cell_iterators())
552 * {
553 * cell_matrix = 0;
554 *
555 * fe_values.reinit(cell);
556 *
557 * fe_values.get_function_gradients(evaluation_point,
558 * evaluation_point_gradients);
559 *
560 * for (unsigned int q = 0; q < n_q_points; ++q)
561 * {
562 * const double coeff =
563 * 1.0 / std::sqrt(1 + evaluation_point_gradients[q] *
564 * evaluation_point_gradients[q]);
565 *
566 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
567 * {
568 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
569 * cell_matrix(i, j) +=
570 * (((fe_values.shape_grad(i, q) // ((\nabla \phi_i
571 * * coeff // * a_n
572 * * fe_values.shape_grad(j, q)) // * \nabla \phi_j)
573 * - // -
574 * (fe_values.shape_grad(i, q) // (\nabla \phi_i
575 * * coeff * coeff * coeff // * a_n^3
576 * *
577 * (fe_values.shape_grad(j, q) // * (\nabla \phi_j
578 * * evaluation_point_gradients[q]) // * \nabla u_n)
579 * * evaluation_point_gradients[q])) // * \nabla u_n)))
580 * * fe_values.JxW(q)); // * dx
581 * }
582 * }
583 *
584 * cell->get_dof_indices(local_dof_indices);
585 * hanging_node_constraints.distribute_local_to_global(cell_matrix,
586 * local_dof_indices,
587 * jacobian_matrix);
588 * }
589 *
590 * std::map<types::global_dof_index, double> boundary_values;
592 * 0,
594 * boundary_values);
595 * Vector<double> dummy_solution(dof_handler.n_dofs());
596 * Vector<double> dummy_rhs(dof_handler.n_dofs());
597 * MatrixTools::apply_boundary_values(boundary_values,
598 * jacobian_matrix,
599 * dummy_solution,
600 * dummy_rhs);
601 * }
602 *
603 * @endcode
604 *
605 * The second half of the function then deals with factorizing the
606 * so-computed matrix. To do this, we first create a new SparseDirectUMFPACK
607 * object and by assigning it to the member variable
608 * `jacobian_matrix_factorization`, we also destroy whatever object that
609 * pointer previously pointed to (if any). Then we tell the object to
610 * factorize the Jacobian.
611 *
612
613 *
614 * As above, we enclose this block of code into curly braces and use a timer
615 * to assess how long this part of the program takes.
616 *
617
618 *
619 * (Strictly speaking, we don't actually need the matrix any more after we
620 * are done here, and could throw the matrix object away. A code intended to
621 * be memory efficient would do this, and only create the matrix object in
622 * this function, rather than as a member variable of the surrounding class.
623 * We omit this step here because using the same coding style as in previous
624 * tutorial programs breeds familiarity with the common style and helps make
625 * these tutorial programs easier to read.)
626 *
627 * @code
628 * {
629 * TimerOutput::Scope t(computing_timer, "factorizing the Jacobian");
630 *
631 * std::cout << " Factorizing Jacobian matrix" << std::endl;
632 *
633 * jacobian_matrix_factorization = std::make_unique<SparseDirectUMFPACK>();
634 * jacobian_matrix_factorization->factorize(jacobian_matrix);
635 * }
636 * }
637 *
638 *
639 *
640 * @endcode
641 *
642 *
643 * <a name="Computingtheresidualvector"></a>
644 * <h4>Computing the residual vector</h4>
645 *
646
647 *
648 * The second part of what `assemble_system()` used to do in @ref step_15 "step-15" is
649 * computing the residual vector, i.e., the right hand side vector of the
650 * Newton linear systems. We have broken this out of the previous function,
651 * but the following function will be easy to understand if you understood
652 * what `assemble_system()` in @ref step_15 "step-15" did. Importantly, however, we need to
653 * compute the residual not linearized around the current solution vector, but
654 * whatever we get from KINSOL. This is necessary for operations such as line
655 * search where we want to know what the residual @f$F(U^k + \alpha_k \delta
656 * U^K)@f$ is for different values of @f$\alpha_k@f$; KINSOL in those cases simply
657 * gives us the argument to the function @f$F@f$ and we then compute the residual
658 * @f$F(\cdot)@f$ at this point.
659 *
660
661 *
662 * The function prints the norm of the so-computed residual at the end as a
663 * way for us to follow along the progress of the program.
664 *
665 * @code
666 * template <int dim>
667 * void MinimalSurfaceProblem<dim>::compute_residual(
668 * const Vector<double> &evaluation_point,
669 * Vector<double> & residual)
670 * {
671 * TimerOutput::Scope t(computing_timer, "assembling the residual");
672 *
673 * std::cout << " Computing residual vector..." << std::flush;
674 *
675 * const QGauss<dim> quadrature_formula(fe.degree + 1);
676 * FEValues<dim> fe_values(fe,
677 * quadrature_formula,
678 * update_gradients | update_quadrature_points |
679 * update_JxW_values);
680 *
681 * const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
682 * const unsigned int n_q_points = quadrature_formula.size();
683 *
684 * Vector<double> cell_residual(dofs_per_cell);
685 * std::vector<Tensor<1, dim>> evaluation_point_gradients(n_q_points);
686 *
687 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
688 *
689 * for (const auto &cell : dof_handler.active_cell_iterators())
690 * {
691 * cell_residual = 0;
692 * fe_values.reinit(cell);
693 *
694 * fe_values.get_function_gradients(evaluation_point,
695 * evaluation_point_gradients);
696 *
697 *
698 * for (unsigned int q = 0; q < n_q_points; ++q)
699 * {
700 * const double coeff =
701 * 1.0 / std::sqrt(1 + evaluation_point_gradients[q] *
702 * evaluation_point_gradients[q]);
703 *
704 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
705 * cell_residual(i) = (fe_values.shape_grad(i, q) // \nabla \phi_i
706 * * coeff // * a_n
707 * * evaluation_point_gradients[q] // * u_n
708 * * fe_values.JxW(q)); // * dx
709 * }
710 *
711 * cell->get_dof_indices(local_dof_indices);
712 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
713 * residual(local_dof_indices[i]) += cell_residual(i);
714 * }
715 *
716 * hanging_node_constraints.condense(residual);
717 *
718 * for (const types::global_dof_index i :
719 * DoFTools::extract_boundary_dofs(dof_handler))
720 * residual(i) = 0;
721 *
722 * for (const types::global_dof_index i :
723 * DoFTools::extract_hanging_node_dofs(dof_handler))
724 * residual(i) = 0;
725 *
726 * std::cout << " norm=" << residual.l2_norm() << std::endl;
727 * }
728 *
729 *
730 *
731 * @endcode
732 *
733 *
734 * <a name="SolvinglinearsystemswiththeJacobianmatrix"></a>
735 * <h4>Solving linear systems with the Jacobian matrix</h4>
736 *
737
738 *
739 * Next up is the function that implements the solution of a linear system
740 * with the Jacobian matrix. Since we have already factored the matrix when we
741 * built the matrix, solving a linear system comes down to applying the
742 * inverse matrix to the given right hand side vector: This is what the
743 * SparseDirectUMFPACK::vmult() function does that we use here. Following
744 * this, we have to make sure that we also address the values of hanging nodes
745 * in the solution vector, and this is done using
746 * AffineConstraints::distribute().
747 *
748
749 *
750 * The function takes an additional, but unused, argument `tolerance` that
751 * indicates how accurately we have to solve the linear system. The meaning of
752 * this argument is discussed in the introduction in the context of the
753 * "Eisenstat Walker trick", but since we are using a direct rather than an
754 * iterative solver, we are not using this opportunity to solve linear systems
755 * only inexactly.
756 *
757 * @code
758 * template <int dim>
759 * void MinimalSurfaceProblem<dim>::solve(const Vector<double> &rhs,
760 * Vector<double> & solution,
761 * const double /*tolerance*/)
762 * {
763 * TimerOutput::Scope t(computing_timer, "linear system solve");
764 *
765 * std::cout << " Solving linear system" << std::endl;
766 *
767 * jacobian_matrix_factorization->vmult(solution, rhs);
768 *
769 * hanging_node_constraints.distribute(solution);
770 * }
771 *
772 *
773 *
774 * @endcode
775 *
776 *
777 * <a name="Refiningthemeshsettingboundaryvaluesandgeneratinggraphicaloutput"></a>
778 * <h4>Refining the mesh, setting boundary values, and generating graphical output</h4>
779 *
780
781 *
782 * The following three functions are again simply copies of the ones in
783 * @ref step_15 "step-15":
784 *
785 * @code
786 * template <int dim>
787 * void MinimalSurfaceProblem<dim>::refine_mesh()
788 * {
789 * Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
790 *
791 * KellyErrorEstimator<dim>::estimate(
792 * dof_handler,
793 * QGauss<dim - 1>(fe.degree + 1),
794 * std::map<types::boundary_id, const Function<dim> *>(),
795 * current_solution,
796 * estimated_error_per_cell);
797 *
798 * GridRefinement::refine_and_coarsen_fixed_number(triangulation,
799 * estimated_error_per_cell,
800 * 0.3,
801 * 0.03);
802 *
803 * triangulation.prepare_coarsening_and_refinement();
804 *
805 * SolutionTransfer<dim> solution_transfer(dof_handler);
806 * solution_transfer.prepare_for_coarsening_and_refinement(current_solution);
807 *
808 * triangulation.execute_coarsening_and_refinement();
809 *
810 * dof_handler.distribute_dofs(fe);
811 *
812 * Vector<double> tmp(dof_handler.n_dofs());
813 * solution_transfer.interpolate(current_solution, tmp);
814 * current_solution = std::move(tmp);
815 *
816 * hanging_node_constraints.clear();
817 *
818 * DoFTools::make_hanging_node_constraints(dof_handler,
819 * hanging_node_constraints);
820 * hanging_node_constraints.close();
821 *
822 * hanging_node_constraints.distribute(current_solution);
823 *
824 * set_boundary_values();
825 *
826 * setup_system(/*initial_step=*/false);
827 * }
828 *
829 *
830 *
831 * template <int dim>
832 * void MinimalSurfaceProblem<dim>::set_boundary_values()
833 * {
834 * std::map<types::global_dof_index, double> boundary_values;
835 * VectorTools::interpolate_boundary_values(dof_handler,
836 * 0,
837 * BoundaryValues<dim>(),
838 * boundary_values);
839 * for (const auto &boundary_value : boundary_values)
840 * current_solution(boundary_value.first) = boundary_value.second;
841 *
842 * hanging_node_constraints.distribute(current_solution);
843 * }
844 *
845 *
846 *
847 * template <int dim>
848 * void MinimalSurfaceProblem<dim>::output_results(
849 * const unsigned int refinement_cycle)
850 * {
851 * TimerOutput::Scope t(computing_timer, "graphical output");
852 *
853 * DataOut<dim> data_out;
854 *
855 * data_out.attach_dof_handler(dof_handler);
856 * data_out.add_data_vector(current_solution, "solution");
857 * data_out.build_patches();
858 *
859 * const std::string filename =
860 * "solution-" + Utilities::int_to_string(refinement_cycle, 2) + ".vtu";
861 * std::ofstream output(filename);
862 * data_out.write_vtu(output);
863 * }
864 *
865 *
866 *
867 * @endcode
868 *
869 *
870 * <a name="Therunfunctionandtheoveralllogicoftheprogram"></a>
871 * <h4>The run() function and the overall logic of the program</h4>
872 *
873
874 *
875 * The only function that *really* is interesting in this program is the one
876 * that drives the overall algorithm of starting on a coarse mesh, doing some
877 * mesh refinement cycles, and on each mesh using KINSOL to find the solution
878 * of the nonlinear algebraic equation we obtain from discretization on this
879 * mesh. The `refine_mesh()` function above makes sure that the solution on
880 * one mesh is used as the starting guess on the next mesh. We also use a
881 * TimerOutput object to measure how much time every operation on each mesh
882 * costs, and reset the timer at the beginning of each cycle.
883 *
884
885 *
886 * As discussed in the introduction, it is not necessary to solve problems on
887 * coarse meshes particularly accurately since these will only solve as
888 * starting guesses for the next mesh. As a consequence, we will use a target
889 * tolerance of
890 * @f$\tau=10^{-3} \frac{1}{10^k}@f$ for the @f$k@f$th mesh refinement cycle.
891 *
892
893 *
894 * All of this is encoded in the first part of this function:
895 *
896 * @code
897 * template <int dim>
898 * void MinimalSurfaceProblem<dim>::run()
899 * {
900 * GridGenerator::hyper_ball(triangulation);
901 * triangulation.refine_global(2);
902 *
903 * setup_system(/*initial_step=*/true);
904 * set_boundary_values();
905 *
906 * for (unsigned int refinement_cycle = 0; refinement_cycle < 6;
907 * ++refinement_cycle)
908 * {
909 * computing_timer.reset();
910 * std::cout << "Mesh refinement step " << refinement_cycle << std::endl;
911 *
912 * if (refinement_cycle != 0)
913 * refine_mesh();
914 *
915 * const double target_tolerance = 1e-3 * std::pow(0.1, refinement_cycle);
916 * std::cout << " Target_tolerance: " << target_tolerance << std::endl
917 * << std::endl;
918 *
919 * @endcode
920 *
921 * This is where the fun starts. At the top we create the KINSOL solver
922 * object and feed it with an object that encodes a number of additional
923 * specifics (of which we only change the nonlinear tolerance we want to
924 * reach; but you might want to look into what other members the
925 * SUNDIALS::KINSOL::AdditionalData class has and play with them).
926 *
927 * @code
928 * {
929 * typename SUNDIALS::KINSOL<Vector<double>>::AdditionalData
930 * additional_data;
931 * additional_data.function_tolerance = target_tolerance;
932 *
933 * SUNDIALS::KINSOL<Vector<double>> nonlinear_solver(additional_data);
934 *
935 * @endcode
936 *
937 * Then we have to describe the operations that were already mentioned
938 * in the introduction. In essence, we have to teach KINSOL how to (i)
939 * resize a vector to the correct size, (ii) compute the residual
940 * vector, (iii) compute the Jacobian matrix (during which we also
941 * compute its factorization), and (iv) solve a linear system with the
942 * Jacobian.
943 *
944
945 *
946 * All four of these operations are represented by member variables of
947 * the SUNDIALS::KINSOL class that are of type `std::function`, i.e.,
948 * they are objects to which we can assign a pointer to a function or,
949 * as we do here, a "lambda function" that takes the appropriate
950 * arguments and returns the appropriate information. By convention,
951 * KINSOL wants that functions doing something nontrivial return an
952 * integer where zero indicates success. It turns out that we can do
953 * all of this in just 25 lines of code.
954 *
955
956 *
957 * (If you're not familiar what "lambda functions" are, take
958 * a look at @ref step_12 "step-12" or at the
959 * [wikipedia page](https://en.wikipedia.org/wiki/Anonymous_function)
960 * on the subject. The idea of lambda functions is that one
961 * wants to define a function with a certain set of
962 * arguments, but (i) not make it a named functions because,
963 * typically, the function is used in only one place and it
964 * seems unnecessary to give it a global name; and (ii) that
965 * the function has access to some of the variables that
966 * exist at the place where it is defined, including member
967 * variables. The syntax of lambda functions is awkward, but
968 * ultimately quite useful.)
969 *
970
971 *
972 * At the very end of the code block we then tell KINSOL to go to work
973 * and solve our problem. The member functions called from the
974 * 'residual', 'setup_jacobian', and 'solve_jacobian_system' functions
975 * will then print output to screen that allows us to follow along
976 * with the progress of the program.
977 *
978 * @code
979 * nonlinear_solver.reinit_vector = [&](Vector<double> &x) {
980 * x.reinit(dof_handler.n_dofs());
981 * };
982 *
983 * nonlinear_solver.residual =
984 * [&](const Vector<double> &evaluation_point,
985 * Vector<double> & residual) {
986 * compute_residual(evaluation_point, residual);
987 *
988 * return 0;
989 * };
990 *
991 * nonlinear_solver.setup_jacobian =
992 * [&](const Vector<double> &current_u,
993 * const Vector<double> & /*current_f*/) {
994 * compute_and_factorize_jacobian(current_u);
995 *
996 * return 0;
997 * };
998 *
999 * nonlinear_solver.solve_with_jacobian = [&](const Vector<double> &rhs,
1000 * Vector<double> & dst,
1001 * const double tolerance) {
1002 * this->solve(rhs, dst, tolerance);
1003 *
1004 * return 0;
1005 * };
1006 *
1007 * nonlinear_solver.solve(current_solution);
1008 * }
1009 *
1010 * @endcode
1011 *
1012 * The rest is then just house-keeping: Writing data to a file for
1013 * visualizing, and showing a summary of the timing collected so that we
1014 * can interpret how long each operation has taken, how often it was
1015 * executed, etc:
1016 *
1017 * @code
1018 * output_results(refinement_cycle);
1019 *
1020 * computing_timer.print_summary();
1021 *
1022 * std::cout << std::endl;
1023 * }
1024 * }
1025 * } // namespace Step77
1026 *
1027 *
1028 * int main()
1029 * {
1030 * try
1031 * {
1032 * using namespace Step77;
1033 *
1034 * MinimalSurfaceProblem<2> laplace_problem_2d;
1035 * laplace_problem_2d.run();
1036 * }
1037 * catch (std::exception &exc)
1038 * {
1039 * std::cerr << std::endl
1040 * << std::endl
1041 * << "----------------------------------------------------"
1042 * << std::endl;
1043 * std::cerr << "Exception on processing: " << std::endl
1044 * << exc.what() << std::endl
1045 * << "Aborting!" << std::endl
1046 * << "----------------------------------------------------"
1047 * << std::endl;
1048 *
1049 * return 1;
1050 * }
1051 * catch (...)
1052 * {
1053 * std::cerr << std::endl
1054 * << std::endl
1055 * << "----------------------------------------------------"
1056 * << std::endl;
1057 * std::cerr << "Unknown exception!" << std::endl
1058 * << "Aborting!" << std::endl
1059 * << "----------------------------------------------------"
1060 * << std::endl;
1061 * return 1;
1062 * }
1063 * return 0;
1064 * }
1065 * @endcode
1066<a name="Results"></a><h1>Results</h1>
1067
1068
1069When running the program, you get output that looks like this:
1070@code
1071Mesh refinement step 0
1072 Target_tolerance: 0.001
1073
1074 Computing residual vector... norm=0.231202
1075 Computing Jacobian matrix
1076 Factorizing Jacobian matrix
1077 Solving linear system
1078 Computing residual vector... norm=0.231202
1079 Computing residual vector... norm=0.171585
1080 Solving linear system
1081 Computing residual vector... norm=0.171585
1082 Computing residual vector... norm=0.127245
1083 Computing residual vector... norm=0.0796471
1084 Solving linear system
1085 Computing residual vector... norm=0.0796471
1086 Computing residual vector... norm=0.0625301
1087 Solving linear system
1088 Computing residual vector... norm=0.0625301
1089 Computing residual vector... norm=0.0498864
1090 Solving linear system
1091 Computing residual vector... norm=0.0498864
1092 Computing residual vector... norm=0.0407765
1093 Solving linear system
1094 Computing residual vector... norm=0.0407765
1095 Computing residual vector... norm=0.0341589
1096 Solving linear system
1097 Computing residual vector... norm=0.0341589
1098 Computing residual vector... norm=0.0292867
1099 Solving linear system
1100 Computing residual vector... norm=0.0292867
1101 Computing residual vector... norm=0.0256309
1102 Computing residual vector... norm=0.0223448
1103 Solving linear system
1104 Computing residual vector... norm=0.0223448
1105 Computing residual vector... norm=0.0202797
1106 Computing residual vector... norm=0.0183817
1107 Solving linear system
1108 Computing residual vector... norm=0.0183817
1109 Computing residual vector... norm=0.0170464
1110 Computing residual vector... norm=0.0157967
1111 Computing Jacobian matrix
1112 Factorizing Jacobian matrix
1113 Solving linear system
1114 Computing residual vector... norm=0.0157967
1115 Computing residual vector... norm=0.0141572
1116 Computing residual vector... norm=0.012657
1117 Solving linear system
1118 Computing residual vector... norm=0.012657
1119 Computing residual vector... norm=0.0116863
1120 Computing residual vector... norm=0.0107696
1121 Solving linear system
1122 Computing residual vector... norm=0.0107696
1123 Computing residual vector... norm=0.0100986
1124 Computing residual vector... norm=0.00944829
1125 Computing residual vector... norm=0.00822576
1126 Solving linear system
1127 Computing residual vector... norm=0.00822576
1128 Computing residual vector... norm=0.00781983
1129 Computing residual vector... norm=0.00741619
1130 Computing residual vector... norm=0.00661792
1131 Solving linear system
1132 Computing residual vector... norm=0.00661792
1133 Computing residual vector... norm=0.00630571
1134 Computing residual vector... norm=0.00599457
1135 Computing residual vector... norm=0.00537663
1136 Solving linear system
1137 Computing residual vector... norm=0.00537663
1138 Computing residual vector... norm=0.00512813
1139 Computing residual vector... norm=0.00488033
1140 Computing residual vector... norm=0.00438751
1141 Computing residual vector... norm=0.00342052
1142 Solving linear system
1143 Computing residual vector... norm=0.00342052
1144 Computing residual vector... norm=0.00326581
1145 Computing residual vector... norm=0.00311176
1146 Computing residual vector... norm=0.00280617
1147 Computing residual vector... norm=0.00220992
1148 Solving linear system
1149 Computing residual vector... norm=0.00220992
1150 Computing residual vector... norm=0.00209976
1151 Computing residual vector... norm=0.00199943
1152 Solving linear system
1153 Computing residual vector... norm=0.00199942
1154 Computing residual vector... norm=0.00190953
1155 Computing residual vector... norm=0.00182005
1156 Computing residual vector... norm=0.00164259
1157 Computing residual vector... norm=0.00129652
1158
1159
1160+---------------------------------------------+------------+------------+
1161| Total wallclock time elapsed since start | 0.192s | |
1162| | | |
1163| Section | no. calls | wall time | % of total |
1164+---------------------------------+-----------+------------+------------+
1165| assembling the Jacobian | 2 | 0.0141s | 7.4% |
1166| assembling the residual | 61 | 0.168s | 88% |
1167| factorizing the Jacobian | 2 | 0.0016s | 0.83% |
1168| graphical output | 1 | 0.00385s | 2% |
1169| linear system solve | 19 | 0.0013s | 0.68% |
1170+---------------------------------+-----------+------------+------------+
1171
1172
1173Mesh refinement step 1
1174 Target_tolerance: 0.0001
1175
1176 Computing residual vector... norm=0.0883422
1177 Computing Jacobian matrix
1178 Factorizing Jacobian matrix
1179 Solving linear system
1180 Computing residual vector... norm=0.0883422
1181 Computing residual vector... norm=0.0607066
1182 Solving linear system
1183 Computing residual vector... norm=0.0607066
1184 Computing residual vector... norm=0.0437266
1185 Solving linear system
1186 Computing residual vector... norm=0.0437266
1187 Computing residual vector... norm=0.0327999
1188 Solving linear system
1189 Computing residual vector... norm=0.0327999
1190 Computing residual vector... norm=0.0255418
1191 Solving linear system
1192 Computing residual vector... norm=0.0255417
1193 Computing residual vector... norm=0.0206042
1194 Solving linear system
1195 Computing residual vector... norm=0.0206042
1196 Computing residual vector... norm=0.0171602
1197 Solving linear system
1198 Computing residual vector... norm=0.0171602
1199 Computing residual vector... norm=0.014689
1200 Solving linear system
1201
1202[...]
1203@endcode
1204
1205The way this should be interpreted is most easily explained by looking at
1206the first few lines of the output on the first mesh:
1207@code
1208Mesh refinement step 0
1209Mesh refinement step 0
1210 Target_tolerance: 0.001
1211
1212 Computing residual vector... norm=0.231202
1213 Computing Jacobian matrix
1214 Factorizing Jacobian matrix
1215 Solving linear system
1216 Computing residual vector... norm=0.231202
1217 Computing residual vector... norm=0.171585
1218 Solving linear system
1219 Computing residual vector... norm=0.171585
1220 Computing residual vector... norm=0.127245
1221 Computing residual vector... norm=0.0796471
1222 Solving linear system
1223 Computing residual vector... norm=0.0796471
1224 ...
1225@endcode
1226What is happening is this:
1227- In the first residual computation, KINSOL computes the residual to see whether
1228 the desired tolerance has been reached. The answer is no, so it requests the
1229 user program to compute the Jacobian matrix (and the function then also
1230 factorizes the matrix via SparseDirectUMFPACK).
1231- KINSOL then instructs us to solve a linear system of the form
1232 @f$J_k \, \delta U_k = -F_k@f$ with this matrix and the previously computed
1233 residual vector.
1234- It is then time to determine how far we want to go in this direction,
1235 i.e., do line search. To this end, KINSOL requires us to compute the
1236 residual vector @f$F(U_k + \alpha_k \delta U_k)@f$ for different step lengths
1237 @f$\alpha_k@f$. For the first step above, it finds an acceptable @f$\alpha_k@f$
1238 after two tries, the second time around it takes three tries.
1239- Having found a suitable updated solution @f$U_{k+1}@f$, the process is
1240 repeated except now KINSOL is happy with the current Jacobian matrix
1241 and does not instruct us to re-build the matrix and its factorization,
1242 and instead asks us to solve a linear system with that same matrix.
1243
1244The program also writes the solution to a VTU file at the end
1245of each mesh refinement cycle, and it looks as follows:
1246<table width="60%" align="center">
1247 <tr>
1248 <td align="center">
1249 <img src="https://www.dealii.org/images/steps/developer/step-77.solution.png" alt="">
1250 </td>
1251 </tr>
1252</table>
1253
1254
1255The key takeaway messages of this program are the following:
1256
1257- The solution is the same as the one we computed in @ref step_15 "step-15", i.e., the
1258 interfaces to %SUNDIALS' KINSOL package really did what they were supposed
1259 to do. This should not come as a surprise, but the important point is that
1260 we don't have to spend the time implementing the complex algorithms that
1261 underlie advanced nonlinear solvers ourselves.
1262
1263- KINSOL is able to avoid all sorts of operations such as rebuilding the
1264 Jacobian matrix when that is not actually necessary. Comparing the
1265 number of linear solves in the output above with the number of times
1266 we rebuild the Jacobian and compute its factorization should make it
1267 clear that this leads to very substantial savings in terms of compute
1268 times, without us having to implement the intricacies of algorithms
1269 that determine when we need to rebuild this information.
1270
1271<a name="extensions"></a>
1272<a name="Possibilitiesforextensions"></a><h3> Possibilities for extensions </h3>
1273
1274
1275For all but the small problems we consider here, a sparse direct solver
1276requires too much time and memory -- we need an iterative solver like
1277we use in many other programs. The trade-off between constructing an
1278expensive preconditioner (say, a geometric or algebraic multigrid method)
1279is different in the current case, however: Since we can re-use the same
1280matrix for numerous linear solves, we can do the same for the preconditioner
1281and putting more work into building a good preconditioner can more easily
1282be justified than if we used it only for a single linear solve as one
1283does for many other situations.
1284
1285But iterative solvers also afford other opportunities. For example (and as
1286discussed briefly in the introduction), we may not need to solve to
1287very high accuracy (small tolerances) in early nonlinear iterations as long
1288as we are still far away from the actual solution. This was the basis of the
1289Eisenstat-Walker trick mentioned there.
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.
1296 *
1297 *
1298<a name="PlainProg"></a>
1299<h1> The plain program</h1>
1300@include "step-77.cc"
1301*/
@ wall_times
Definition: timer.h:653
Definition: vector.h:110
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
Point< 2 > second
Definition: grid_out.cc:4588
Point< 2 > first
Definition: grid_out.cc:4587
__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, SparsityPatternType &sparsity_pattern, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
@ matrix
Contents is actually a matrix.
static const types::blas_int one
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)
Definition: matrix_tools.cc:81
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())
int(&) functions(const void *v1, const void *v2)
static constexpr double PI
Definition: numbers.h:231
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation