475 * <a name=
"step_77-ThecodeMinimalSurfaceProblemcodeclassimplementation"></a>
481 * <a name=
"step_77-Constructorandsetupfunctions"></a>
506 *
dof_handler.distribute_dofs(fe);
520 *
sparsity_pattern.copy_from(
dsp);
530 * <a name=
"step_77-AssemblingandfactorizingtheJacobianmatrix"></a>
556 *
std::cout <<
" Computing Jacobian matrix" << std::endl;
567 *
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
574 *
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
576 *
for (
const auto &cell : dof_handler.active_cell_iterators())
580 *
fe_values.reinit(cell);
585 *
for (
unsigned int q = 0;
q < n_q_points; ++
q)
587 *
const double coeff =
591 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
593 *
for (
unsigned int j = 0;
j < dofs_per_cell; ++
j)
595 *
(((fe_values.shape_grad(i,
q)
597 *
* fe_values.shape_grad(
j,
q))
599 *
(fe_values.shape_grad(i,
q)
602 *
(fe_values.shape_grad(
j,
q)
605 *
* fe_values.JxW(
q));
609 *
cell->get_dof_indices(local_dof_indices);
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.)
654 * TimerOutput::Scope t(computing_timer, "factorizing the Jacobian");
656 * std::cout << " Factorizing Jacobian matrix" << std::endl;
658 * jacobian_matrix_factorization = std::make_unique<SparseDirectUMFPACK>();
659 * jacobian_matrix_factorization->factorize(jacobian_matrix);
668 * <a name="step_77-Computingtheresidualvector"></a>
669 * <h4>Computing the residual vector</h4>
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.
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.
692 * void MinimalSurfaceProblem<dim>::compute_residual(
693 * const Vector<double> &evaluation_point,
694 * Vector<double> &residual)
696 * TimerOutput::Scope t(computing_timer, "assembling the residual");
698 * std::cout << " Computing residual vector..." << std::flush;
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);
708 * const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
709 * const unsigned int n_q_points = quadrature_formula.size();
711 * Vector<double> cell_residual(dofs_per_cell);
712 * std::vector<Tensor<1, dim>> evaluation_point_gradients(n_q_points);
714 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
716 * for (const auto &cell : dof_handler.active_cell_iterators())
719 * fe_values.reinit(cell);
721 * fe_values.get_function_gradients(evaluation_point,
722 * evaluation_point_gradients);
725 * for (unsigned int q = 0; q < n_q_points; ++q)
727 * const double coeff =
728 * 1.0 / std::sqrt(1 + evaluation_point_gradients[q] *
729 * evaluation_point_gradients[q]);
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
735 * * evaluation_point_gradients[q] // * \nabla u_n
736 * * fe_values.JxW(q)); // * dx
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);
744 * hanging_node_constraints.condense(residual);
746 * for (const types::global_dof_index i :
747 * DoFTools::extract_boundary_dofs(dof_handler))
750 * for (const types::global_dof_index i :
751 * DoFTools::extract_hanging_node_dofs(dof_handler))
754 * std::cout << " norm=" << residual.l2_norm() << std::endl;
762 * <a name="step_77-SolvinglinearsystemswiththeJacobianmatrix"></a>
763 * <h4>Solving linear systems with the Jacobian matrix</h4>
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().
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
787 * void MinimalSurfaceProblem<dim>::solve(const Vector<double> &rhs,
788 * Vector<double> &solution,
789 * const double /*tolerance*/)
791 * TimerOutput::Scope t(computing_timer, "linear system solve");
793 * std::cout << " Solving linear system" << std::endl;
795 * jacobian_matrix_factorization->vmult(solution, rhs);
797 * hanging_node_constraints.distribute(solution);
805 * <a name="step_77-Refiningthemeshsettingboundaryvaluesandgeneratinggraphicaloutput"></a>
806 * <h4>Refining the mesh, setting boundary values, and generating graphical output</h4>
810 * The following three functions are again simply copies of the ones in
811 * @ref step_15 "step-15":
815 * void MinimalSurfaceProblem<dim>::refine_mesh()
817 * Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
819 * KellyErrorEstimator<dim>::estimate(
821 * QGauss<dim - 1>(fe.degree + 1),
822 * std::map<types::boundary_id, const Function<dim> *>(),
824 * estimated_error_per_cell);
826 * GridRefinement::refine_and_coarsen_fixed_number(triangulation,
827 * estimated_error_per_cell,
831 * triangulation.prepare_coarsening_and_refinement();
833 * SolutionTransfer<dim> solution_transfer(dof_handler);
834 * solution_transfer.prepare_for_coarsening_and_refinement(current_solution);
836 * triangulation.execute_coarsening_and_refinement();
838 * dof_handler.distribute_dofs(fe);
840 * Vector<double> tmp(dof_handler.n_dofs());
841 * solution_transfer.interpolate(current_solution, tmp);
842 * current_solution = std::move(tmp);
844 * hanging_node_constraints.clear();
846 * DoFTools::make_hanging_node_constraints(dof_handler,
847 * hanging_node_constraints);
848 * hanging_node_constraints.close();
850 * hanging_node_constraints.distribute(current_solution);
852 * set_boundary_values();
854 * setup_system(/*initial_step=*/false);
860 * void MinimalSurfaceProblem<dim>::set_boundary_values()
862 * std::map<types::global_dof_index, double> boundary_values;
863 * VectorTools::interpolate_boundary_values(dof_handler,
865 * BoundaryValues<dim>(),
867 * for (const auto &boundary_value : boundary_values)
868 * current_solution(boundary_value.first) = boundary_value.second;
870 * hanging_node_constraints.distribute(current_solution);
876 * void MinimalSurfaceProblem<dim>::output_results(
877 * const unsigned int refinement_cycle)
879 * TimerOutput::Scope t(computing_timer, "graphical output");
881 * DataOut<dim> data_out;
883 * data_out.attach_dof_handler(dof_handler);
884 * data_out.add_data_vector(current_solution, "solution");
885 * data_out.build_patches();
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);
898 * <a name="step_77-Therunfunctionandtheoveralllogicoftheprogram"></a>
899 * <h4>The run() function and the overall logic of the program</h4>
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.
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
918 * @f$\tau=10^{-3} \frac{1}{10^k}@f$ for the @f$k@f$th mesh refinement cycle.
922 * All of this is encoded in the first part of this function:
926 * void MinimalSurfaceProblem<dim>::run()
928 * GridGenerator::hyper_ball(triangulation);
929 * triangulation.refine_global(2);
931 * setup_system(/*initial_step=*/true);
932 * set_boundary_values();
934 * for (unsigned int refinement_cycle = 0; refinement_cycle < 6;
935 * ++refinement_cycle)
937 * computing_timer.reset();
938 * std::cout << "Mesh refinement step " << refinement_cycle << std::endl;
940 * if (refinement_cycle != 0)
943 * const double target_tolerance = 1e-3 * std::pow(0.1, refinement_cycle);
944 * std::cout << " Target_tolerance: " << target_tolerance << std::endl
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).
957 * typename SUNDIALS::KINSOL<Vector<double>>::AdditionalData
959 * additional_data.function_tolerance = target_tolerance;
961 * SUNDIALS::KINSOL<Vector<double>> nonlinear_solver(additional_data);
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
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.
1000 *
'residual',
'setup_jacobian',
and 'solve_with_jacobian' functions
1006 *
x.reinit(dof_handler.n_dofs());
1023 *
const double tolerance) {
1024 *
solve(
rhs, dst, tolerance);
1042 *
std::cout << std::endl;
1052 *
using namespace Step77;
1057 *
catch (std::exception &exc)
1059 *
std::cerr << std::endl
1061 *
<<
"----------------------------------------------------"
1063 *
std::cerr <<
"Exception on processing: " << std::endl
1064 *
<< exc.what() << std::endl
1065 *
<<
"Aborting!" << std::endl
1066 *
<<
"----------------------------------------------------"
1073 *
std::cerr << std::endl
1075 *
<<
"----------------------------------------------------"
1077 *
std::cerr <<
"Unknown exception!" << std::endl
1078 *
<<
"Aborting!" << std::endl
1079 *
<<
"----------------------------------------------------"
1146+---------------------------------------------+------------+------------+
1150+---------------------------------+-----------+------------+------------+
1154|
graphical output | 1 | 0.00504s | 4.2% |
1155| linear
system solve | 15 | 0.000893s | 0.74% |
1156+---------------------------------+-----------+------------+------------+
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.
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">
1243 <img src="https://www.dealii.org/images/steps/developer/step-77.solution.png" alt="">
1249The key takeaway messages of this program are the following:
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.
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.
1266<a name="step-77-extensions"></a>
1267<a name="step_77-Possibilitiesforextensions"></a><h3> Possibilities for extensions </h3>
1270<a name="step_77-Betterlinearsolvers"></a><h4> Better linear solvers </h4>
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.
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".)
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
1300 void MinimalSurfaceProblem<dim>::solve(const Vector<double> &rhs,
1301 Vector<double> & solution,
1304i.e., the `tolerance` parameter already exists, but is unused.
1307<a name="step_77-ReplacingSUNDIALSKINSOLbyPETScsSNES"></a><h4> Replacing SUNDIALS' KINSOL by PETSc's SNES </h4>
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.
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.
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:
1337Mesh refinement step 0
1338 Target_tolerance: 0.001
1340 Computing residual vector
1342 Computing Jacobian matrix
1343 Computing residual vector
1344 Computing residual vector
1346 Computing Jacobian matrix
1347 Computing residual vector
1348 Computing residual vector
1350 Computing Jacobian matrix
1351 Computing residual vector
1352 Computing residual vector
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 successful line-search.
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.
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
1374Mesh refinement step 0
1376SNES Object: 1 MPI process
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
1387 interpolation: cubic
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
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
1401 linear system matrix followed by preconditioner matrix:
1402 Mat Object: 1 MPI process
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
1412 total: nonzeros=745, allocated nonzeros=745
1413 total number of mallocs used during MatSetValues calls=0
1414 not using I-node routines
1417From the above details, we see that we are using the "newtonls" solver
1418type ("Newton line search"), with "bt" ("backtracting") line search.
1420From the output of **-snes_view** we can also get information about
1485computations. In this case the solution process will require many more
1486nonlinear iterations since the accuracy of the linear system solve is
1489If we switch to using the preconditioned conjugate gradient method as
1490a linear solve, still using our initial factorization as
1491preconditioner, we get:
1493./step-77-snes -snes_lag_jacobian 2 -ksp_type cg | grep "Computing Jacobian" | wc -l
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.
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:
1506Mesh refinement step 0
1507 Target_tolerance: 0.001
1509 Computing residual vector
1511 Computing Jacobian matrix
1512 Computing residual vector
1513 Computing residual vector
1515 Computing residual vector
1516 Computing residual vector
1518 Computing residual vector
1519 Computing residual vector
1521 Computing residual vector
1522 Computing residual vector
1523 Computing residual vector
1525 Computing residual vector
1526 Computing residual vector
1527 Computing residual vector
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:
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
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
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
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**:
1687Mesh refinement step 0
1688 Target_tolerance: 0.001
1690 Computing residual vector
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.
1697 ||J - Jfd||_F/||J||_F = 0.0196815, ||J - Jfd||_F = 0.503436
1699 Hand-coded minus finite-difference Jacobian with tolerance 1e-05 ----------
1700Mat Object: 1 MPI process
1703row 1: (1, 0.0437112)
1713row 11: (11, 1.38157)
1717showing that the only errors we commit in assembling the Jacobian are
1718on the boundary dofs. As discussed in the tutorial, those errors are
1721The key take-away messages of this modification of the tutorial program are
1722therefore basically the same of what we already found using KINSOL:
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
1739<a name=
"step_77-ReplacingSUNDIALSKINSOLbyTrilinosNOXpackage"></a><
h4>
Replacing SUNDIALS' KINSOL by Trilinos' NOX package </h4>
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
1748vector and matrix classes.
1751<a name="step_77-ReplacingSUNDIALSKINSOLbyagenericnonlinearsolver"></a><h4> Replacing SUNDIALS' KINSOL
by a
generic nonlinear solver </
h4>
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.
1764<a name="step_77-PlainProg"></a>
1765<h1> The plain program</h1>
1766@include "step-77.cc"
numbers::NumberTraits< Number >::real_type norm() const
__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={}, 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.)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
constexpr const ReferenceCell Line
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
const Iterator const std_cxx20::type_identity_t< Iterator > & end
const InputIterator OutputIterator const Function & function
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation