Reference documentation for deal.II version GIT relicensing-489-g2d48aca8cc 2024-04-28 17:30:02+00:00
|
This tutorial depends on step-15.
Table of contents | |
---|---|
This program was contributed by Wolfgang Bangerth, Colorado State University.
This material is based upon work partially supported by National Science Foundation grants OAC-1835673, DMS-1821210, and EAR-1925595; and by the Computational Infrastructure in Geodynamics initiative (CIG), through the National Science Foundation under Award No. EAR-1550901 and The University of California-Davis.
Stefano Zampini (King Abdullah University of Science and Technology) contributed the results obtained with the PETSc variant of this program discussed in the results section below.
The step-15 program solved the following, nonlinear equation describing the minimal surface problem:
\begin{align*} -\nabla \cdot \left( \frac{1}{\sqrt{1+|\nabla u|^{2}}}\nabla u \right) &= 0 \qquad \qquad &&\textrm{in} ~ \Omega \\ u&=g \qquad\qquad &&\textrm{on} ~ \partial \Omega. \end{align*}
step-15 uses a Newton method, and Newton's method works by repeatedly solving a linearized problem for an update \(\delta u_k\) – called the "search direction" –, computing a "step length" \(\alpha_k\), and then combining them to compute the new guess for the solution via
\begin{align*} u_{k+1} = u_k + \alpha_k \, \delta u_k. \end{align*}
In the course of the discussions in step-15, we found that it is awkward to compute the step length, and so just settled for simple choice: Always choose \(\alpha_k=0.1\). This is of course not efficient: We know that we can only realize Newton's quadratic convergence rate if we eventually are able to choose \(\alpha_k=1\), though we may have to choose it smaller for the first few iterations where we are still too far away to use this long a step length.
Among the goals of this program is therefore to address this shortcoming. Since line search algorithms are not entirely trivial to implement, one does as one should do anyway: Import complicated functionality from an external library. To this end, we will make use of the interfaces deal.II has to one of the big nonlinear solver packages, namely the KINSOL sub-package of the SUNDIALS suite. SUNDIALS is, at its heart, a package meant to solve complex ordinary differential equations (ODEs) and differential-algebraic equations (DAEs), and the deal.II interfaces allow for this via the classes in the SUNDIALS namespace: Notably the SUNDIALS::ARKode and SUNDIALS::IDA classes. But, because that is an important step in the solution of ODEs and DAEs with implicit methods, SUNDIALS also has a solver for nonlinear problems called KINSOL, and deal.II has an interface to it in the form of the SUNDIALS::KINSOL class. This is what we will use for the solution of our problem.
But SUNDIALS isn't just a convenient way for us to avoid writing a line search algorithm. In general, the solution of nonlinear problems is quite expensive, and one typically wants to save as much compute time as possible. One way one can achieve this is as follows: The algorithm in step-15 discretizes the problem and then in every iteration solves a linear system of the form
\begin{align*} J_k \, \delta U_k = -F_k \end{align*}
where \(F_k\) is the residual vector computed using the current vector of nodal values \(U_k\), \(J_k\) is its derivative (called the "Jacobian"), and \(\delta U_k\) is the update vector that corresponds to the function \(\delta u_k\) mentioned above. The construction of \(J_k,F_k\) has been thoroughly discussed in step-15, as has the way to solve the linear system in each Newton iteration. So let us focus on another aspect of the nonlinear solution procedure: Computing \(F_k\) is expensive, and assembling the matrix \(J_k\) even more so. Do we actually need to do that in every iteration? It turns out that in many applications, this is not actually necessary: These methods often converge even if we replace \(J_k\) by an approximation \(\tilde J_k\) and solve
\begin{align*} \tilde J_k \, \widetilde{\delta U}_k = -F_k \end{align*}
instead, then update
\begin{align*} U_{k+1} = U_k + \alpha_k \, \widetilde{\delta U}_k. \end{align*}
This may require an iteration or two more because our update \(\widetilde{\delta U}_k\) is not quite as good as \(\delta U_k\), but it may still be a win because we don't have to assemble \(J_k\) quite as often.
What kind of approximation \(\tilde J_k\) would we like for \(J_k\)? Theory says that as \(U_k\) converges to the exact solution \(U^\ast\), we need to ensure that \(\tilde J_k\) needs to converge to \(J^\ast = \nabla F(U^\ast)\). In particular, since \(J_k\rightarrow J^\ast\), a valid choice is \(\tilde J_k = J_k\). But so is choosing \(\tilde J_k = J_k\) every, say, fifth iteration \(k=0,5,10,\ldots\) and for the other iterations, we choose \(\tilde J_k\) equal to the last computed \(J_{k'}\). This is what we will do here: we will just re-use \(\tilde J_{k-1}\) from the previous iteration, which may again be what we had used in the iteration before that, \(\tilde J_{k-2}\).
This scheme becomes even more interesting if, for the solution of the linear system with \(J_k\), we don't just have to assemble a matrix, but also compute a good preconditioner. For example, if we were to use a sparse LU decomposition via the SparseDirectUMFPACK class, or used a geometric or algebraic multigrid. In those cases, we would also not have to update the preconditioner, whose computation may have taken about as long or longer than the assembly of the matrix in the first place. Indeed, with this mindset, we should probably think about using the best preconditioner we can think of, even though their construction is typically quite expensive: We will hope to amortize the cost of computing this preconditioner by applying it to more than one just one linear solve.
The big question is, of course: By what criterion do we decide whether we can get away with the approximation \(\tilde J_k\) based on a previously computed Jacobian matrix \(J_{k-s}\) that goes back \(s\) steps, or whether we need to – at least in this iteration – actually re-compute the Jacobian \(J_k\) and the corresponding preconditioner? This is, like the issue with line search, one that requires a non-trivial amount of code that monitors the convergence of the overall algorithm. We could implement these sorts of things ourselves, but we probably shouldn't: KINSOL already does that for us. It will tell our code when to "update" the Jacobian matrix.
One last consideration if we were to use an iterative solver instead of the sparse direct one mentioned above: Not only is it possible to get away with replacing \(J_k\) by some approximation \(\tilde J_k\) when solving for the update \(\delta U_k\), but one can also ask whether it is necessary to solve the linear system
\begin{align*} \tilde J_k \widetilde{\delta U}_k = -F_k \end{align*}
to high accuracy. The thinking goes like this: While our current solution \(U_k\) is still far away from \(U^\ast\), why would we solve this linear system particularly accurately? The update \(U_{k+1}=U_k + \widetilde{\delta U}_k\) is likely still going to be far away from the exact solution, so why spend much time on solving the linear system to great accuracy? This is the kind of thinking that underlies algorithms such as the "Eisenstat-Walker trick" [184] in which one is given a tolerance to which the linear system above in iteration \(k\) has to be solved, with this tolerance dependent on the progress in the overall nonlinear solver. As before, one could try to implement this oneself, but KINSOL already provides this kind of information for us – though we will not use it in this program since we use a direct solver that requires no solver tolerance and just solves the linear system exactly up to round-off.
As a summary of all of these considerations, we could say the following: There is no need to reinvent the wheel. Just like deal.II provides a vast amount of finite-element functionality, SUNDIALS' KINSOL package provides a vast amount of nonlinear solver functionality, and we better use it.
KINSOL, like many similar packages, works in a pretty abstract way. At its core, it sees a nonlinear problem of the form
\begin{align*} F(U) = 0 \end{align*}
and constructs a sequence of iterates \(U_k\) which, in general, are vectors of the same length as the vector returned by the function \(F\). To do this, there are a few things it needs from the user:
All of these operations need to be provided to KINSOL by std::function objects that take the appropriate set of arguments and that generally return an integer that indicates success (a zero return value) or failure (a nonzero return value). Specifically, the objects we will access are the SUNDIALS::KINSOL::reinit_vector, SUNDIALS::KINSOL::residual, SUNDIALS::KINSOL::setup_jacobian, and SUNDIALS::KINSOL::solve_with_jacobian member variables. (See the documentation of these variables for their details.) In our implementation, we will use lambda functions to implement these "callbacks" that in turn can call member functions; KINSOL will then call these callbacks whenever its internal algorithms think it is useful.
The majority of the code of this tutorial program is as in step-15, and we will not comment on it in much detail. There is really just one aspect one has to pay some attention to, namely how to compute \(F(U)\) given a vector \(U\) on the one hand, and \(J(U)\) given a vector \(U\) separately. At first, this seems trivial: We just take the assemble_system()
function and in the one case throw out all code that deals with the matrix and in the other case with the right hand side vector. There: Problem solved.
But it isn't quite as simple. That's because the two are not independent if we have nonzero Dirichlet boundary values, as we do here. The linear system we want to solve contains both interior and boundary degrees of freedom, and when eliminating those degrees of freedom from those that are truly "free", using for example AffineConstraints::distribute_local_to_global(), we need to know the matrix when assembling the right hand side vector.
Of course, this completely contravenes the original intent: To not assemble the matrix if we can get away without it. We solve this problem as follows:
There is an assumption here that whenever KINSOL asks for a linear solver with the (approximation of the) Jacobian, that this will be for an update \(\delta U\) (which has zero boundary values), a multiple of which will be added to the solution (which already has the right boundary values). This may not be true and if so, we might have to rethink our approach. That said, it turns out that in practice this is exactly what KINSOL does when using a Newton method, and so our approach is successful.
This program starts out like most others with well known include files. Compared to the step-15 program from which most of what we do here is copied, the only difference is the include of the header files from which we import the SparseDirectUMFPACK class and the actual interface to KINSOL:
MinimalSurfaceProblem
class templateSimilarly, the main class of this program is essentially a copy of the one in step-15. The class does, however, split the computation of the Jacobian (system) matrix (and its factorization using a direct solver) and residual into separate functions for the reasons outlined in the introduction. For the same reason, the class also has a pointer to a factorization of the Jacobian matrix that is reset every time we update the Jacobian matrix.
(If you are wondering why the program uses a direct object for the Jacobian matrix but a pointer for the factorization: Every time KINSOL requests that the Jacobian be updated, we can simply write jacobian_matrix=0;
to reset it to an empty matrix that we can then fill again. On the other hand, the SparseDirectUMFPACK class does not have any way to throw away its content or to replace it with a new factorization, and so we use a pointer: We just throw away the whole object and create a new one whenever we have a new Jacobian matrix to factor.)
Finally, the class has a timer variable that we will use to assess how long the different parts of the program take so that we can assess whether KINSOL's tendency to not rebuild the matrix and its factorization makes sense. We will discuss this in the "Results" section below.
The classes implementing boundary values are a copy from step-15:
MinimalSurfaceProblem
class implementationThe following few functions are also essentially copies of what step-15 already does, and so there is little to discuss.
The following function is then responsible for assembling and factorizing the Jacobian matrix. The first half of the function is in essence the assemble_system()
function of step-15, except that it does not deal with also forming a right hand side vector (i.e., the residual) since we do not always have to do these operations at the same time.
We put the whole assembly functionality into a code block enclosed by curly braces so that we can use a TimerOutput::Scope variable to measure how much time is spent in this code block, excluding everything that happens in this function after the matching closing brace }
.
The second half of the function then deals with factorizing the so-computed matrix. To do this, we first create a new SparseDirectUMFPACK object and by assigning it to the member variable jacobian_matrix_factorization
, we also destroy whatever object that pointer previously pointed to (if any). Then we tell the object to factorize the Jacobian.
As above, we enclose this block of code into curly braces and use a timer to assess how long this part of the program takes.
(Strictly speaking, we don't actually need the matrix any more after we are done here, and could throw the matrix object away. A code intended to be memory efficient would do this, and only create the matrix object in this function, rather than as a member variable of the surrounding class. We omit this step here because using the same coding style as in previous tutorial programs breeds familiarity with the common style and helps make these tutorial programs easier to read.)
The second part of what assemble_system()
used to do in step-15 is computing the residual vector, i.e., the right hand side vector of the Newton linear systems. We have broken this out of the previous function, but the following function will be easy to understand if you understood what assemble_system()
in step-15 did. Importantly, however, we need to compute the residual not linearized around the current solution vector, but whatever we get from KINSOL. This is necessary for operations such as line search where we want to know what the residual \(F(U^k + \alpha_k \delta
U^K)\) is for different values of \(\alpha_k\); KINSOL in those cases simply gives us the argument to the function \(F\) and we then compute the residual \(F(\cdot)\) at this point.
The function prints the norm of the so-computed residual at the end as a way for us to follow along the progress of the program.
Next up is the function that implements the solution of a linear system with the Jacobian matrix. Since we have already factored the matrix when we built the matrix, solving a linear system comes down to applying the inverse matrix to the given right hand side vector: This is what the SparseDirectUMFPACK::vmult() function does that we use here. Following this, we have to make sure that we also address the values of hanging nodes in the solution vector, and this is done using AffineConstraints::distribute().
The function takes an additional, but unused, argument tolerance
that indicates how accurately we have to solve the linear system. The meaning of this argument is discussed in the introduction in the context of the "Eisenstat Walker trick", but since we are using a direct rather than an iterative solver, we are not using this opportunity to solve linear systems only inexactly.
The following three functions are again simply copies of the ones in step-15:
The only function that really is interesting in this program is the one that drives the overall algorithm of starting on a coarse mesh, doing some mesh refinement cycles, and on each mesh using KINSOL to find the solution of the nonlinear algebraic equation we obtain from discretization on this mesh. The refine_mesh()
function above makes sure that the solution on one mesh is used as the starting guess on the next mesh. We also use a TimerOutput object to measure how much time every operation on each mesh costs, and reset the timer at the beginning of each cycle.
As discussed in the introduction, it is not necessary to solve problems on coarse meshes particularly accurately since these will only solve as starting guesses for the next mesh. As a consequence, we will use a target tolerance of \(\tau=10^{-3} \frac{1}{10^k}\) for the \(k\)th mesh refinement cycle.
All of this is encoded in the first part of this function:
This is where the fun starts. At the top we create the KINSOL solver object and feed it with an object that encodes a number of additional specifics (of which we only change the nonlinear tolerance we want to reach; but you might want to look into what other members the SUNDIALS::KINSOL::AdditionalData class has and play with them).
Then we have to describe the operations that were already mentioned in the introduction. In essence, we have to teach KINSOL how to (i) resize a vector to the correct size, (ii) compute the residual vector, (iii) compute the Jacobian matrix (during which we also compute its factorization), and (iv) solve a linear system with the Jacobian.
All four of these operations are represented by member variables of the SUNDIALS::KINSOL class that are of type std::function
, i.e., they are objects to which we can assign a pointer to a function or, as we do here, a "lambda function" that takes the appropriate arguments and returns the appropriate information. It turns out that we can do all of this in just over 20 lines of code.
(If you're not familiar what "lambda functions" are, take a look at step-12 or at the wikipedia page on the subject. The idea of lambda functions is that one wants to define a function with a certain set of arguments, but (i) not make it a named functions because, typically, the function is used in only one place and it seems unnecessary to give it a global name; and (ii) that the function has access to some of the variables that exist at the place where it is defined, including member variables. The syntax of lambda functions is awkward, but ultimately quite useful.)
At the very end of the code block we then tell KINSOL to go to work and solve our problem. The member functions called from the 'residual', 'setup_jacobian', and 'solve_with_jacobian' functions will then print output to screen that allows us to follow along with the progress of the program.
The rest is then just house-keeping: Writing data to a file for visualizing, and showing a summary of the timing collected so that we can interpret how long each operation has taken, how often it was executed, etc:
When running the program, you get output that looks like this:
The way this should be interpreted is most easily explained by looking at the first few lines of the output on the first mesh:
What is happening is this:
The program also writes the solution to a VTU file at the end of each mesh refinement cycle, and it looks as follows:
The key takeaway messages of this program are the following:
For all but the small problems we consider here, a sparse direct solver requires too much time and memory – we need an iterative solver like we use in many other programs. The trade-off between constructing an expensive preconditioner (say, a geometric or algebraic multigrid method) is different in the current case, however: Since we can re-use the same matrix for numerous linear solves, we can do the same for the preconditioner and putting more work into building a good preconditioner can more easily be justified than if we used it only for a single linear solve as one does for many other situations.
But iterative solvers also afford other opportunities. For example (and as discussed briefly in the introduction), we may not need to solve to very high accuracy (small tolerances) in early nonlinear iterations as long as we are still far away from the actual solution. This was the basis of the Eisenstat-Walker trick mentioned there. (This is also the underlying reason why one can store the matrix in single precision rather than double precision, see the discussion in the "Possibilities for extensions" section of step-15.)
KINSOL provides the function that does the linear solution with a target tolerance that needs to be reached. We ignore it in the program above because the direct solver we use does not need a tolerance and instead solves the linear system exactly (up to round-off, of course), but iterative solvers could make use of this kind of information – and, in fact, should. Indeed, the infrastructure is already there: The solve()
function of this program is declared as
i.e., the tolerance
parameter already exists, but is unused.
As mentioned in the introduction, SUNDIALS' KINSOL package is not the only player in town. Rather, very similar interfaces exist to the SNES package that is part of PETSc, and the NOX package that is part of Trilinos, via the PETScWrappers::NonlinearSolver and TrilinosWrappers::NOXSolver classes.
It is not very difficult to change the program to use either of these two alternatives. Rather than show exactly what needs to be done, let us point out that a version of this program that uses SNES instead of KINSOL is available as part of the test suite, in the file tests/petsc/step-77-snes.cc
. Setting up the solver for PETScWrappers::NonlinearSolver turns out to be even simpler than for the SUNDIALS::KINSOL class we use here because we don't even need the reinit
lambda function – SNES only needs us to set up the remaining three functions residual
, setup_jacobian
, and solve_with_jacobian
. The majority of changes necessary to convert the program to use SNES are related to the fact that SNES can only deal with PETSc vectors and matrices, and these need to be set up slightly differently. On the upside, the test suite program mentioned above already works in parallel.
SNES also allows playing with a number of parameters about the solver, and that enables some interesting comparisons between methods. When you run the test program (or a slightly modified version that outputs information to the screen instead of a file), you get output that looks comparable to something like this:
By default, PETSc uses a Newton solver with cubic backtracking, resampling the Jacobian matrix at each Newton step. That is, we compute and factorize the matrix once per Newton step, and then sample the residual to check for a successful line-search.
The attentive reader should have noticed that in this case we are computing one more extra residual per Newton step. This is because the deal.II code is set up to use a Jacobian-free approach, and the extra residual computation pops up when computing a matrix-vector product to test the validity of the Newton solution.
PETSc can be configured in many interesting ways via the command line. We can visualize the details of the solver by using the command line argument -snes_view, which produces the excerpt below at the end of each solve call:
From the above details, we see that we are using the "newtonls" solver type ("Newton line search"), with "bt" ("backtracting") line search.
From the output of -snes_view we can also get information about the linear solver details; specifically, when using the solve_with_jacobian
interface, the deal.II interface internally uses a custom solver configuration within a "shell" preconditioner, that wraps the action of solve_with_jacobian
.
We can also see the details of the type of matrices used within the solve: "mffd" (matrix-free finite-differencing) for the action of the linearized operator and "seqaij" for the assembled Jacobian we have used to construct the preconditioner.
Diagnostics for the line search procedure can be turned on using the command line -snes_linesearch_monitor, producing the excerpt below:
Within the run, the Jacobian matrix is assembled (and factored) 29 times:
KINSOL internally decided when it was necessary to update the Jacobian matrix (which is when it would call setup_jacobian
). SNES can do something similar: We can compute the explicit sparse Jacobian matrix only once per refinement step (and reuse the initial factorization) by using the command line -snes_lag_jacobian -2, producing:
In other words, this dramatically reduces the number of times we have to build the Jacobian matrix, though at a cost to the number of nonlinear steps we have to take.
The lagging period can also be decided automatically. For example, if we want to recompute the Jacobian at every other step:
Note, however, that we didn't exactly halve the number of Jacobian computations. In this case the solution process will require many more nonlinear iterations since the accuracy of the linear system solve is not enough.
If we switch to using the preconditioned conjugate gradient method as a linear solve, still using our initial factorization as preconditioner, we get:
Note that in this case we use an approximate preconditioner (the LU factorization of the initial approximation) but we use a matrix-free operator for the action of our Jacobian matrix, thus solving for the correct linear system.
We can switch to a quasi-Newton method by using the command line -snes_type qn -snes_qn_scale_type jacobian, and we can see that our Jacobian is sampled and factored only when needed, at the cost of an increase of the number of steps:
Nonlinear preconditioning can also be used. For example, we can run a right-preconditioned nonlinear GMRES, using one Newton step as a preconditioner, with the command:
As also discussed for the KINSOL use above, optimal preconditioners should be used instead of the LU factorization used here by default. This is already possible within this tutorial by playing with the command line options. For example, algebraic multigrid can be used by simply specifying -pc_type gamg. When using iterative linear solvers, the "Eisenstat-Walker trick" [184] can be also requested at command line via -snes_ksp_ew. Using these options, we can see that the number of nonlinear iterations used by the solver increases as the mesh is refined, and that the number of linear iterations increases as the Newton solver is entering the second-order ball of convergence:
Finally we describe how to get some diagnostic on the correctness of the computed Jacobian. Deriving the correct linearization is sometimes difficult: It took a page or two in the introduction to derive the exact bilinear form for the Jacobian matrix, and it would be quite nice compute it automatically from the residual of which it is the derivative. (This is what step-72 does!) But if one is set on doing things by hand, it would at least be nice if we had a way to check the correctness of the derivation. SNES allows us to do this: we can use the options -snes_test_jacobian -snes_test_jacobian_view:
showing that the only errors we commit in assembling the Jacobian are on the boundary dofs. As discussed in the tutorial, those errors are harmless.
The key take-away messages of this modification of the tutorial program are therefore basically the same of what we already found using KINSOL:
Besides KINSOL and SNES, the third option you have is to use the NOX package. As before, rather than showing in detail how that needs to happen, let us simply point out that the test suite program tests/trilinos/step-77-with-nox.cc
does this. The modifications necessary to use NOX instead of KINSOL are quite minimal; in particular, NOX (unlike SNES) is happy to work with deal.II's own vector and matrix classes.
Having to choose which of these three frameworks (KINSOL, SNES, or NOX) to use at compile time is cumbersome when wanting to compare things. It would be nicer if one could decide the package to use at run time, assuming that one has a copy of deal.II installed that is compiled against all three of these dependencies. It turns out that this is possible, using the class NonlinearSolverSelector that presents a common interface to all three of these solvers, along with the ability to choose which one to use based on run-time parameters.