 |
deal.II version GIT relicensing-2287-g6548a49e0a 2024-12-20 18:30:00+00:00
|
Loading...
Searching...
No Matches
|
deal.II version GIT relicensing-2287-g6548a49e0a 2024-12-20 18:30:00+00:00
|
This tutorial depends on step-26, step-40.
This program was contributed by Wolfgang Bangerth (Colorado State University), Stefano Zampini (King Abdullah University of Science and Technology), and Luca Heltai (University of Pisa).
This material is based upon work partially supported by National Science Foundation grants OAC-1835673, DMS-1821210, and EAR-1925595.
step-26 solved the simple heat equation, one of the prototypical examples of time dependent problems:
\begin{align*} \frac{\partial u(\mathbf x, t)}{\partial t} - \Delta u(\mathbf x, t) &= f(\mathbf x, t), \qquad\qquad && \forall \mathbf x \in \Omega, t\in (0,T), \\ u(\mathbf x, 0) &= u_0(\mathbf x) && \forall \mathbf x \in \Omega, \\ u(\mathbf x, t) &= g(\mathbf x,t) && \forall \mathbf x \in \partial\Omega, t \in (0,T). \end{align*}
While that program showed a number of advanced techniques such as using adaptive mesh refinement for time-dependent problems, it did not address one big issue: It hand-rolls its own time stepping scheme, which in that program is the simple Crank-Nicolson method with a fixed time step. This is neither accurate nor efficient: We should be using a higher-order time stepping algorithm, and we should use one of the many ways to efficiently and automatically choose the length of the time step in response to the accuracy obtained.
This would of course require quite a lot of development effort – unless, of course, you do what we always advise: You build on what others have already done and have likely done in a way far superior to what one can do by oneself. In the current case, deal.II has interfaces to two such libraries: SUNDIALS, the SUite of Nonlinear and DIfferential/ALgebraic equation Solvers (and here specifically the Runge-Kutta-type solvers wrapped in the SUNDIALS::ARKode class), and PETSc's TS sub-package (wrapped in the PETScWrappers::TimeStepper class). In this program, we will use the PETSc TS interfaces.
While we're at it with updating step-26, we will also make the program run in parallel – a minor change if you've read step-40, for example.
Both the PETSc TS and the SUNDIALS interfaces require that we first write the partial differential equation in the form of a system of ordinary differential equations. To this end, let us turn around the approach we used in step-26. There, we first discretized in time, obtaining a PDE to be solved at each time step that we could then discretize using the finite element method. This approach is called the "Rothe method". Instead, here, we use what's called the "Method of Lines" where we first discretize in space, obtaining a system of ordinary differential equations to which we can apply traditional time steppers. (There are some trade-offs between these two strategies, principally around using dynamically changing meshes; we will get back to this issue later on.)
To get this started, we take the equation above and multiply it by a test function \(\varphi(\mathbf x)\) and integrate by parts to get a weak form: We seek a function \(u(\mathbf x, t)\) that for all test functions \(\varphi \in H^1_0(\Omega)\) satisfies
\begin{align*} \left(\varphi(\mathbf x), \frac{\partial u(\mathbf x, t)}{\partial t} \right)_\Omega + \left(\nabla \varphi(\mathbf x), \nabla u(\mathbf x, t) \right)_\Omega &= \left(\varphi(\mathbf x), f(\mathbf x, t) \right)_\Omega, \\ \left(\varphi(\mathbf x), u(\mathbf x, 0)\right)_\Omega &= \left(\varphi(\mathbf x), u_0(\mathbf x)\right)_\Omega, && \\ u(\mathbf x, t) &= g(\mathbf x,t) && \forall \mathbf x \in \partial\Omega, t \in (0,T). \end{align*}
(Integration by parts would ordinarily result in boundary terms unless one has Dirichlet boundary conditions – possibly nonzero – all around the boundary. We will assume that this is indeed so herein.)
We then discretize by restricting ourself to finite element functions of the form
\begin{align*} u_h(\mathbf x,t) = \sum_j U_j(t) \varphi_j(\mathbf x), \end{align*}
which leads to the problem of finding a function \(u_h(\mathbf x, t)\) that for all discrete test functions \(\varphi \in V_h(\Omega)\) satisfies
\begin{align*} \left(\varphi_i(\mathbf x), \frac{\partial u_h(\mathbf x, t)}{\partial t} \right)_\Omega + \left(\nabla \varphi_j(\mathbf x), \nabla u_h(\mathbf x, t) \right)_\Omega &= \left(\varphi_i(\mathbf x), f(\mathbf x, t) \right)_\Omega, \\ \left(\varphi_i(\mathbf x), u_h(\mathbf x, 0)\right)_\Omega &= \left(\varphi_i(\mathbf x), u_0(\mathbf x)\right)_\Omega, && \\ u_h(\mathbf x, t) &= g_h(\mathbf x,t) && \forall \mathbf x \in \partial\Omega, t \in (0,T), \end{align*}
where \(g_h\) is an interpolant of the function \(g\) on the boundary.
This equation can be rewritten in matrix form in the usual way, by expanding \(u_h\) into its coefficients times shape function form, pulling the sum over \(j\) out of the integrals, and then considering that choosing test function \(\varphi_i\) leads to the \(i\)th row of the linear system. This then gives us
\begin{align*} M \frac{\partial U(t)}{\partial t} + AU(t) &= F(t), \\ U(0) = U_0, \end{align*}
plus appropriate boundary conditions.
There are now two perspectives on how one should proceed. If we were to use the SUNDIALS::ARKode wrappers to solve this linear system, we would bring the \(AU\) term to the right hand side and consider the ODE
\begin{align*} M \frac{\partial U(t)}{\partial t} &= - AU(t) + F(t), \end{align*}
which matches the form stated in the documentation of SUNDIALS::ARKode. In particular, ARKode is able to deal with the fact that the time derivative is multiplied by the mass matrix \(M\), which is always there when using finite elements.
On the other hand, when using the PETScWrappers::TimeStepper class, we can solve ODEs that are stated in a general "implicit" form, and in that case we simply bring everything to the left hand side and obtain
\begin{align*} \underbrace{ M \frac{\partial U(t)}{\partial t} + AU(t) - F(t) }_{=:R(t,U,\dot U)} = 0. \end{align*}
This matches the form \(R(t,U,\dot U) = 0\) you can find in the documentation of PETScWrappers::TimeStepper if you identify the time dependent function \(y=y(t)\) used there with our solution vector \(U(t)\), and our notation \(R(t,U,\dot U)\) instead of the \(F(t,y,\dot y)\) used there and which we rename because we want to use \(F\) as the right hand side vector of the ODE indicating forcing terms.
This program uses the PETScWrappers::TimeStepper class, and so we will take the latter viewpoint. (It is worth noting that SUNDIALS also has a package that can solve ODEs in implicit form, wrapped by the SUNDIALS::IDA class.) In what follows, we will continue to use \(U(t)\) as the function we seek, even though the documentation of the class uses \(y(t)\).
Having identified how we want to see the problem (namely, as an "implicit" ODE), the question is how we describe the problem to the time stepper. Conceptually, all of the wrappers for time stepping packages we support in deal.II only requires us to provide them with a very limited set of operations. Specifically, for the implicit formulation used by PETScWrappers::TimeStepper, all we need to implement are functions that provide the following:
That's really it. If we can provide these three functions, PETSc will do the rest (as would, for example, SUNDIALS::ARKode or, if you prefer the implicit form, SUNDIALS::IDA). It will not be very difficult to set these things up. In practice, the way this will work is that inside the run() function, we will set up lambda functions that can access the information of the surrounding scopes and that return the requested information.
In practice, we often want to provide a fourth function:
While we like to say that all nodes in a finite element mesh are "degrees of freedom", this is not actually true if we have Dirichlet boundary conditions: Degrees of "freedom" located along Dirichlet boundaries have specific values that are prescribed by the boundary conditions and are, consequently, not "free". Moreover, while the form
\begin{align*} M \frac{\partial U(t)}{\partial t} &= - AU(t) + F(t) \end{align*}
suggests that all elements of the vector \(U(t)\) satisfy a differential equation, this is not actually true for those components of \(U\) that correspond to boundary nodes; rather, their values are prescribed and will in general not satisfy the equation. On second thought, you will also find that the same sort of issue happens as well with handing nodes: These, too, are not really "free" but constrained to values that are derived from the values of neighboring nodes.
For the same reason as we will also discuss in the next section, all of this is easier using the Rothe method we have used in all previous tutorials. There, we end up with a PDE at every time step for which we can independently prescribe boundary conditions and hanging node constraints, and we then deal with those by modifying the matrix and right hand side appropriately. Here, with the method of lines, things are slightly more complicated.
Not too complicated, however: Like with the mesh refinement issues of the next section, Dirichlet boundary conditions (and constrained degrees of freedom in general) are something every PDE solver has to deal with, and because the people who write ODE solvers also have PDEs in mind, they needed to address these cases too and so the interfaces we use are prepared for it. Specifically, what we need to do is mark which entries of the solution vector (i.e., which degrees of freedom) are "algebraic" – that is, satisfy an algebraic, rather than a differential, equation. The way we will do this is that the ODE integrator interface requires us to provide a "callback" function that it can call and that needs to return an IndexSet object in which all the algebraic degrees of freedom are listed. At the end of each solution stage, a second callback then needs to provide the ability to take a solution vector and set its constrained (algebraic) entries to their correct values; for us, that will mean setting boundary values and hanging node constraints correctly.
(We note that when dealing with Differential-Algebraic Equations (DAEs), the algebraic components identified by the first callback mentioned above contain hanging nodes, Dirichlet boundary nodes, and degrees of freedom for which the equation does not provide time derivatives. This is not of concern to us here, and so we will not dwell on the issue – or the fact that the second callback in the case of DAEs should set only a subset of the algebraic components – and instead refer to the examples shown in the SUNDIALS::IDA class documentation.)
When stating an ODE in the form
\begin{align*} M \frac{\partial U(t)}{\partial t} &= - AU(t) + F(t), \end{align*}
or one of the reformulations discussed above, there is an implicit assumption that the number of entries in the vector \(U\) stays constant and that each entry continues to correspond to the same quantity. But if you use mesh refinement, this is not the case: The number of unknowns will go up or down whenever you refine or coarsen the mesh, and the 42nd (or any other) degree of freedom may be located at an entirely different physical location after mesh refinement than where it was located below. In other words, the size of vectors and what individual vector entries mean changes when we do mesh refinement. The ODE form we derived above after spatial discretization simply ceases to be meaningful at these times.
This was precisely why, in all previous time-dependent tutorial programs, we have adopted the Rothe approach. There, one first discretizes in time, and obtains a PDE at each time step. This PDE can then be discretized anew – if one really wanted to, with an entirely different mesh in each time step, though we typically don't go that far. On the other hand, being able to use external ODE integrators is* undoubtedly very useful, and so let us see if we can shoehorn the mesh refinement complication into what external ODE integrators do. This is, in practice, not as difficult as it may at first sound because, perhaps not surprisingly, ODE integrators are written by people who want to solve problems like ours, and so they have had to deal with the same kinds of complications we are discussing here.
The way we approach the situation from a conceptual perspective is that we break things into "time slabs". Let's say we want to solve on the time interval \([0,T]\), then we break things into slabs \([0=\tau_0,\tau_1], [\tau_1,\tau_2], \ldots [\tau_{n-1},\tau_n=T]\) where the break points satisfy \(\tau_{k-1}<\tau_k\). On each time slab, we keep the mesh the same, and so we can call into our time integrator. At the end of a time slab, we then save the solution, refine the mesh, set up other data structures, and restore the solution on the new mesh; then we start the time integrator again at the start of the new time slab. This approach guarantees that for the purposes of ODE solvers, we really only ever give them something that can rightfully be considered an ODE system. A disadvantage is that we typically want to refine or coarsen the mesh relatively frequently (in large-scale codes one often chooses to refine and coarsen the mesh every 10-20 time steps), and that limits the efficiency of time integrators: They gain much of their advantage from being able to choose the time step length automatically, but there is often a cost associated with starting up; if the slabs are too short, then neither the start-up cost nor the benefit of potentially long time steps are realized.
In practice, good integrators such as those in PETSc TS can deal with this transparently. We just have to give them a way to call back into our code at the end of each time step to ask whether we want to refine the mesh and do some prep work; and a second function that the integrator can then call to do the actual refinement and interpolate solution vectors from old to new mesh. You will see in the run() function that none of this is difficult to explain to the ODE integrator.
Compared to the code structure of step-26, the program the current tutorial is principally based on, there are essentially two sets of changes:
run() function: You can see the loop over all time steps, and which functions are called where within the loop. Here, however, this is no longer the case. In essence, in run(), we create an object of type PETScWrappers::TimeStepper, and after some set-up, we turn over control to that object's PETScWrappers::TimeStepper::solve() function that contains the loop over time steps and the logic that decides how large the time step needs to be, what needs to happen when, etc. In other words, the details of the program's logic are no longer visible. Instead, what you have to provide to the PETScWrappers::TimeStepper object is a series of "callbacks": Functions that the time stepper can call whenever appropriate. These callbacks are typically small lambda functions that, if the functionality required only takes a few lines of code do exactly that or, otherwise, call larger member functions of the main class.The program solves the heat equation, which with all right hand sides, initial, and boundary values reads as
\begin{align*} \frac{\partial u(\mathbf x, t)}{\partial t} - \Delta u(\mathbf x, t) &= f(\mathbf x, t), \qquad\qquad && \forall \mathbf x \in \Omega, t\in (0,T), \\ u(\mathbf x, 0) &= u_0(\mathbf x) && \forall \mathbf x \in \Omega, \\ u(\mathbf x, t) &= g(\mathbf x,t) && \forall \mathbf x \in \partial\Omega, t \in (0,T). \end{align*}
The right hand side \(f\), initial conditions \(u_0\), and Dirichlet boundary values \(g\) are all specified in an input file heat_equation.prm in which these functions are provided as expressions that are parsed and evaluated at run time using the Functions::ParsedFunction<dim> class. The version of this file that is distributed with the library uses
\begin{align*} f(\mathbf x,t) &= 0, \\ u_0(\mathbf x) &= 0, \\ g(\mathbf x,t) &= \begin{cases} \cos(4\pi t) & \text{if @f$x=-1@f$}, \\ -\cos(4\pi t) & \text{if @f$x=1@f$}, \\ 0 & \text{otherwise} \end{cases} \end{align*}
but this is easily changed.
The program's input file also contains two sections that control the time stepper:
The first of these two sections describes things such as the end time up to which we want to run the program, the initial time step size, and the type of the time stepper (where beuler indicates "backward Euler"; other choices are listed here. We will play with some of these parameters in the results section. As usual when using PETSc solvers, these runtime configuration options can always be complemented (or overridden) via command line options.
The program starts with the usual prerequisites, all of which you will know by now from either step-26 (the heat equation solver) or step-40 (the parallel Laplace solver):
The only new include file relevant here is the one that provides us with the PETScWrappers::TimerStepper class:
At its core, this program's principal structure can be understood quite easily if you know step-26 (for the heat equation) and step-40 (for how a parallel solver looks like). It has many of the usual member functions and member variables that for convenience of documentation we will list towards the top of the class so that we can document the remainder separately below.
We derive the main class from ParameterAcceptor to make dealing with run time parameters and reading them a parameter file easier. step-60 and step-70 have already explained how this works.
At this point, we start to deviate from the "manual" way of solving time dependent problems. High level packages for the solution of ODEs usually expect the user to provide a function that computes the residual of the equation and the time derivative of the solution.
This allows those packages to abstract away the details of how the time derivative is defined (e.g., backward Euler, Crank-Nicolson, etc.) and allow the user to provide a uniform interface, irrespective of the time stepper used. PETSc TS and SUNDIALS are two examples of such packages, and they require the user to provide C style call-backs to compute residuals, Jacobians, etc. In deal.II, we wrap these C style libraries with our own wrappers that use a style closer to c++, and that allows us to simply define the callbacks via lambda functions. Several of these lambda functions will simply call member functions to do the actual work.
To make it clear what we are doing here, we start by defining functions with the same name of the interface that we will use later on (and documented both in the introduction of this program as well in the documentation of the PETScWrappers::TimeStepper class). These are the functions that compute the right hand side, the "Jacobian matrix" (see the introduction for how that is defined), and that can solve a linear system with the Jacobian matrix. At the bottom of the following block, we also declare the matrix object that will store this Jacobian.
Note that all of these functions receive the current solution (and, where necessary, its time derivative) as inputs. This is because the solution vector – i.e., the variable that stores the current state of the solution – is kept inside the time integrator object. This is useful: If we kept a copy of the solution vector as a member variable of the current class, one would continuously have to wonder whether it is still in sync with the version of the solution vector the time integrator object internally believes is the currently correct version of this vector. As a consequence, we do not store such a copy here: Whenever a function requires access to the current value of the solution, it receives it from the time integrator as a const argument. The same observation can be made about the variable that stores the current time, or the current length of the time step, or the number of time steps performed so far: They are all kept inside the time stepping object.
In this tutorial program, similar to what we did in step-26, we want to adapt the mesh at regular time intervals. However, if we are using an external time stepper, we need to make sure that the time stepper is aware of the mesh changes – see the discussion on this topic in the introduction. In particular, the time stepper must support the fact that the number of degrees of freedom may change, and it must support the fact that all internal vectors (e.g., the solution vector and all intermediate stages used to compute the current time derivative) will need to be transferred to the new mesh. For reasons that will be discussed in the implementation below, we split the mesh refinement operation into two functions: The first will mark which cells to refine, based on a given solution vector from which we can compute error indicators; the second one will do the actual mesh refinement and transfer a set of vectors from the old to the new mesh.
As also discussed in the introduction, we also have to deal with "algebraic" solution components, i.e., degrees of freedom that are not really free but instead have values determined either by boundary conditions or by hanging node constraints. While the values of the boundary conditions can change from time step to time step (because the function \(g(\mathbf x,t)\) may indeed depend on time), hanging node constraints remain the same as long as the mesh remains the same. As a consequence, we will keep an AffineConstraints object that stores the hanging node constraints and that is only updated when the mesh changes, and then an AffineConstraints object current_constraints that we will initialize with the hanging node constraints and then add the constraints due to Dirichlet boundary values at a specified time. This evaluation of boundary values and combining of constraints happens in the update_current_constraints() function.
At one place in the program, we will also need an object that constrains the same degrees of freedom, but with zero values even if the boundary values for the solution are non-zero; we will keep this modified set of constraints in homogeneous_constraints.
The remainder of the class is simply consumed with objects that describe either the functioning of the time stepper, when and how to do mesh refinement, and the objects that describe right hand side, initial conditions, and boundary conditions for the PDE. Since we already derive our class from ParameterAcceptor, we exploit its facilities to parse also the parameters of the functions that define the initial value, the right hand side, and the boundary values. We therefore create three ParameterAcceptorProxy objects, which wrap the actual ParsedFunction class into objects that ParameterAcceptor can handle.
The constructor is responsible for initializing all member variables of the class. This is relatively straightforward, and includes setting up parameters to be set upon reading from an input file via the ParameterAcceptor mechanism previously detailed in step-60 and step-70.
This function is not very different from what we do in many other programs (including step-26). We enumerate degrees of freedom, output some information about then, build constraint objects (recalling that we put hanging node constraints into their separate object), and then also build an AffineConstraint object that contains both the hanging node constraints as well as constraints corresponding to zero Dirichlet boundary conditions. This last object is needed since we impose the constraints through algebraic equations. While technically it would be possible to use the time derivative of the boundary function as a boundary conditions for the time derivative of the solution, this is not done here. Instead, we impose the boundary conditions through algebraic equations, and therefore the time derivative of the boundary conditions is not part of the algebraic system, and we need zero boundary conditions on the time derivative of the solution when computing the residual. We use the homogeneous_constraints object for this purpose.
Finally, we create the actual non-homogeneous current_constraints by calling `update_current_constraints). These are also used during the assembly and during the residual evaluation.
The final block of code resets and initializes the matrix object with the appropriate sparsity pattern. Recall that we do not store solution vectors in this class (the time integrator object does that internally) and so do not have to resize and initialize them either.
This function is called from "monitor" function that is called in turns by the time stepper in each time step. We use it to write the solution to a file, and provide graphical output through paraview or visit. We also write a pvd file, which groups all metadata about the .vtu files into a single file that can be used to load the full time dependent solution in paraview.
As discussed in the introduction, we describe the ODE system to the time stepper via its residual,
\[ R(t,U,\dot U) = M \frac{\partial U(t)}{\partial t} + AU(t) - F(t). \]
The following function computes it, given vectors for \(U,\dot U\) that we will denote by y and y_dot because that's how they are called in the documentation of the PETScWrappers::TimeStepper class.
At the top of the function, we do the usual set up when computing integrals. We face two minor difficulties here: the first is that the y and y_dot vectors we get as input are read only, but we need to make sure they satisfy the correct boundary conditions and so have to set elements in these vectors. The second is that we need to compute the residual, and therefore in general we need to evaluate solution values and gradients inside locally owned cells, and for this need access to degrees of freedom which may be owned by neighboring processors. To address these issues, we create (non-ghosted) writable copies of the input vectors, apply boundary conditions and hanging node current_constraints; and then copy these vectors to ghosted vectors before we can do anything sensible with them.
Now for computing the actual residual. Recall that we wan to compute the vector
\[ R(t,U,\dot U) = M \frac{\partial U(t)}{\partial t} + AU(t) - F(t). \]
We could do that by actually forming the matrices \(M\) and \(A\), but this is not efficient. Instead, recall (by writing out how the elements of \(M\) and \(A\) are defined, and exchanging integrals and sums) that the \(i\)th element of the residual vector is given by
\begin{align*} R(t,U,\dot U)_i &= \sum_j \int_\Omega \varphi_i(\mathbf x, t) \varphi_j(\mathbf x, t) {\partial U_j(t)}{\partial t} + \sum_j \int_\Omega \nabla \varphi_i(\mathbf x, t) \cdot \nabla \varphi_j(\mathbf x, t) U_j(t) - \int_\Omega \varphi_i f(\mathbf x, t) \\ &= \int_\Omega \varphi_i(\mathbf x, t) u_h(\mathbf x, t) + \int_\Omega \nabla \varphi_i(\mathbf x, t) \cdot \nabla u_h(\mathbf x, t) - \int_\Omega \varphi_i f(\mathbf x, t). \end{align*}
We can compute these integrals efficiently by breaking them up into a sum over all cells and then applying quadrature. For the integrand, we need to evaluate the solution and its gradient at the quadrature points within each locally owned cell, and for this, we need also access to degrees of freedom that may be owned by neighboring processors. We therefore use the locally_relevant_solution and and locally_relevant_solution_dot vectors.
The end result of the operations above is a vector that contains the residual vector, having taken into account the constraints due to hanging nodes and Dirichlet boundary conditions (by virtue of having used current_constraints.distribute_local_to_global() to add the local contributions to the global vector. At the end of the day, the residual vector \(r\) will be used in the solution of linear systems of the form \(J z = r\) with the "Jacobian" matrix that we define below. We want to achieve that for algebraic components, the algebraic components of \(z\) have predictable values that achieve the purposes discussed in the following. We do this by ensuring that the entries corresponding to algebraic components in the residual \(r\) have specific values, and then we will do the same in the next function for the matrix; for this, you will have to know that the rows and columns of the matrix corresponding to constrained entries are zero with the exception of the diagonal entries. We will manually set that diagonal entry to one, and so \(z_i=r_i\) for algebraic components.
From the point of view of the residual vector, if the input y vector does not contain the correct values on constrained degrees of freedom (hanging nodes or boundary conditions), we need to communicate this to the time stepper, and we do so by setting the residual to the actual difference between the input y vector and the our local copy of it, in which we have applied the constraints (see the top of the function where we called current_constraints.distribute(tmp_solution) and a similar operation on the time derivative). Since we have made a copy of the input vector for this purpose, we use it to compute the residual value. However, there is a difference between hanging nodes constraints and boundary conditions: we do not want to make hanging node constraints actually depend on their dependent degrees of freedom, since this would imply that we are actually solving for the dependent degrees of freedom. This is not what we are actually doing, however, since hanging nodes are not actually solved for. They are eliminated from the system by the call to AffineConstraints::distribute_local_to_global() above. From the point of view of the Jacobian matrix, we are effectively setting hanging nodes to an artificial value (usually zero), and we simply want to make sure that we solve for those degrees of freedom a posteriori, by calling the function AffineConstraints::distribute().
Here we therefore check that the residual is equal to the input value on the constrained dofs corresponding to hanging nodes (i.e., those for which the lines of the current_constraints contain at least one other entry), and to the difference between the input vector and the actual solution on those constraints that correspond to boundary conditions.
The next operation is to compute the "Jacobian", which PETSc TS defines as the matrix
\[ J_\alpha = \dfrac{\partial R}{\partial y} + \alpha \dfrac{\partial R}{\partial \dot y} \]
which, for the current linear problem, is simply
\[ J_\alpha = A + \alpha M \]
and which is in particular independent of time and the current solution vectors \(y\) and \(\dot y\).
Having seen the assembly of matrices before, there is little that should surprise you in the actual assembly here:
The only interesting part is the following. Recall that we modified the residual vector's entries corresponding to the algebraic components of the solution in the previous function. The outcome of calling current_constraints.distribute_local_to_global() a few lines above is that the global matrix has zero rows and columns for the algebraic (constrained) components of the solution; the function puts a value on the diagonal that is nonzero and has about the same size as the remaining diagonal entries of the matrix. What this diagonal value is is unknown to us – in other cases where we call current_constraints.distribute_local_to_global() on both the left side matrix and the right side vector, as in most other tutorial programs, the matrix diagonal entries and right hand side values are chosen in such a way that the result of solving a linear system is what we want it to be, but the scaling is done automatically.
This is not good enough for us here, because we are building the right hand side independently from the matrix in different functions. Thus, for any constrained degree of freedom, we set the diagonal of the Jacobian to be one. This leaves the Jacobian matrix invertible, consistent with what the time stepper expects, and it also makes sure that if we did not make a mistake in the residual and/or in the Jacbian matrix, then asking the time stepper to check the Jacobian with a finite difference method will produce the correct result. This can be activated at run time via passing the -snes_test_jacobian option on the command line.
This is the function that actually solves the linear system with the Jacobian matrix we have previously built (in a call to the previous function during the current time step or another earlier one – time steppers are quite sophisticated in determining internally whether it is necessary to update the Jacobian matrix, or whether one can reuse it for another time step without rebuilding it; this is similar to how one can re-use the Newton matrix for several Newton steps, see for example the discussion in step-77). We could in principle not provide this function to the time stepper, and instead select a specific solver on the command line by using the -ksp_* options of PETSc. However, by providing this function, we can use a specific solver and preconditioner for the linear system, and still have the possibility to change them on the command line.
Providing a specific solver is more in line with the way we usually do things in other deal.II examples, while letting PETSc choose a generic solver, and changing it on the command line via -ksp_type is more in line with the way PETSc is usually used, and it can be a convenient approach when we are experimenting to find an optimal solver for our problem. Both options are available here, since we can still change both the solver and the preconditioner on the command line via -user_ksp_type and -user_pc_type options.
In any case, recall that the Jacobian we built in the previous function is always of the form
\[ J_\alpha = \alpha M + A \]
where \(M\) is a mass matrix and \(A\) a Laplace matrix. \(M\) is symmetric and positive definite; \(A\) is symmetric and at least positive semidefinite; \(\alpha> 0\). As a consequence, the Jacobian matrix is a symmetric and positive definite matrix, which we can efficiently solve with the Conjugate Gradient method, along with either SSOR or (if available) the algebraic multigrid implementation provided by PETSc (via the Hypre package) as preconditioner. In practice, if you wanted to solve "real" problems, one would spend some time finding which preconditioner is optimal, perhaps using PETSc's ability to read solver and preconditioner choices from the command line. But this is not the focus of this tutorial program, and so we just go with the following:
The next block of functions deals with mesh refinement. We split this process up into a "decide whether and what you want to refine" and a "please transfer these vectors from old to new mesh" phase, where the first also deals with marking cells for refinement. (The decision whether or not to refine is done in the lambda function that calls the current function.)
Breaking things into a "mark cells" function and into a "execute mesh adaptation and transfer solution vectors" function is awkward, though conceptually not difficult to understand. These two pieces of code should really be part of the same function, as they are in step-26. The issue is with what PETScWrappers::TimeStepper provides us with in these callbacks. Specifically, the "decide whether and what you want to refine" callback has access to the current solution, and so can evaluate (spatial) error estimators to decide which cells to refine. The second callback that transfers vectors from old to new mesh gets a bunch of vectors, but without the semantic information on which of these is the current solution vector. As a consequence, it cannot do the marking of cells for refinement or coarsening, and we have to do that from the first callback.
In practice, however, the problem is minor. The first of these two functions is essentially identical to the first half of the corresponding function in step-26, with the only difference that it uses the parallel::distributed::GridRefinement namespace function instead of the serial one. Once again, we make sure that we never fall below the minimum refinement level, and above the maximum one, that we can select from the parameter file.
The following function then is the second half of the correspond function in step-26. It is called by the time stepper whenever it requires to transfer the solution and any intermediate stage vectors to a new mesh. We must make sure that all input vectors are transformed into ghosted vectors before the actual transfer is executed, and that we distribute the hanging node constraints on the output vectors as soon as we have interpolated the vectors to the new mesh – i.e., that all constraints are satisfied on the vectors we transfer.
We have no way to enforce boundary conditions at this stage, however. This is because the various vectors may correspond to solutions at previous time steps if the method used here is a multistep time integrator, and so may correspond to different time points that we are not privy to.
While this could be a problem if we used the values of the solution and of the intermediate stages on the constrained degrees of freedom to compute the errors, we do not do so. Instead, we compute the errors on the differential equation, and ignore algebraic constraints; therefore we do no need to guarantee that the boundary conditions are satisfied also in the intermediate stages.
We have at our disposal the hanging node current_constraints alone, though, and therefore we enforce them on the output vectors, even if this is not really needed.
Since regenerating the constraints at each time step may be expensive, we make sure that we only do so when the time changes. We track time change by checking if the time of the boundary_values_function has changed, with respect to the time of the last call to this function. This will work most of the times, but not the very first time we call this function, since the time then may be zero and the time of the boundary_values_function is zero at construction time. We therefore also check if the number constraints in current_constraints, and if these are empty, we regenerate the constraints regardless of the time variable.
We have finally arrived at the main function of the class. At the top, it creates the mesh and sets up the variables that make up the linear system:
We then set up the time stepping object and associate the matrix we will build whenever requested for both the Jacobian matrix (see the definition above of what the "Jacobian" actually refers to) and for the matrix that will be used as a preconditioner for the Jacobian.
The real work setting up the time stepping object starts here. As discussed in the introduction, the way the PETScWrappers::TimeStepper class is used is by inverting control: At the end of this function, we will call PETScWrappers::TimeStepper::solve() which internally will run the loop over time steps, and at the appropriate places call back into user code for whatever functionality is required. What we need to do is to hook up these callback functions by assigning appropriate lambda functions to member variables of the petsc_ts object.
We start by creating lambda functions that provide information about the "implicit function" (i.e., that part of the right hand side of the ODE system that we want to treat implicitly – which in our case is the entire right hand side), a function that assembles the Jacobian matrix, and a function that solves a linear system with the Jacobian.
The next two callbacks deal with identifying and setting variables that are considered "algebraic" (rather than "differential"), i.e., for which we know what values they are supposed to have rather than letting their values be determined by the differential equation. We need to instruct the time stepper to ignore these components when computing the error in the residuals, and we do so by first providing a function that returns an IndexSet with the indices of these algebraic components of the solution vector (or rather, that subset of the locally-owned part of the vector that is algebraic, in case we are running in parallel). This first of the following two functions does that. Specifically, both nodes at which Dirichlet boundary conditions are applied, and hanging nodes are algebraically constrained. This function then returns a set of indices that is initially empty (but knows about the size of the index space) and which we then construct as the union of boundary and hanging node indices.
Following this, we then also need a function that, given a solution vector y and the current time, sets the algebraic components of that vector to their correct value. This comes down to ensuring that we have up to date constraints in the constraints variable, and then applying these constraints to the solution vector via AffineConstraints::distribute(). (It is perhaps worth noting that we could have achieved the same in solve_with_jacobian(). Whenever the time stepper solves a linear system, it follows up the call to the solver by calling the callback to set algebraic components correct. We could also have put the calls to update_current_constraints() and distribute() into the solve_with_jacobian function, but by not doing so, we can also replace the solve_with_jacobian function with a call to a PETSc solver, and we would still have the current_constraints correctly applied to the solution vector.)
The next two callbacks relate to mesh refinement. As discussed in the introduction, PETScWrappers::TimeStepper knows how to deal with the situation where we want to change the mesh. All we have to provide is a callback that returns true if we are at a point where we want to refine the mesh (and false otherwise) and that if we want to do mesh refinement does some prep work for that in the form of calling the prepare_for_coarsening_and_refinement function.
If the first callback below returns true, then PETSc TS will do some clean-up operations, and call the second of the callback functions (petsc_ts.transfer_solution_vectors_to_new_mesh) with a collection of vectors that need to be interpolated from the old to the new mesh. This may include the current solution, perhaps the current time derivative of the solution, and in the case of multistep time integrators also the solutions of some previous time steps. We hand all of these over to the interpolate() member function of this class.
The final callback is a "monitor" that is called in each time step. Here we use it to create a graphical output. Perhaps a better scheme would output the solution at fixed time intervals, rather than in every time step, but this is not the main point of this program and so we go with the easy approach:
With all of this out of the way, the rest of the function is anticlimactic: We just have to initialize the solution vector with the initial conditions and call the function that does the time stepping, and everything else will happen automatically:
The rest of the program is as it always looks. We read run-time parameters from an input file via the ParameterAcceptor class in the same way as we showed in step-60 and step-70.
When you run this program with the input file as is, you get output as follows:
We can generate output for this in the form of a video:
The solution here is driven by boundary values (the initial conditions are zero, and so is the right hand side of the equation). It takes a little bit of time for the boundary values to diffuse into the domain, and so the temperature (the solution of the heat equation) in the interior of the domain has a slight lag compared to the temperature at the boundary.
The more interesting component of this program is how easy it is to play with the details of the time stepping algorithm. Recall that the solution above is controlled by the following parameters:
Of particular interest for us here is to set the time stepping algorithm and the adaptive time step control method. The latter is set to "none" above, but there are several alternative choices for this parameter. For example, we can set parameters as follows,
What we do here is set the initial time step size to 0.025, and choose relatively large absolute and relative error tolerances of 0.01 for the time step size adaptation algorithm (for which we choose "basic"). We ask PETSc TS to use a Backward Differentiation Formula (BDF) method, and we get the following as output:
What is happening here is that apparently PETSc TS is not happy with our choice of initial time step size, and after several linear solves has reduced it to the minimum we allow it to, 0.01. The following time steps then run at a time step size of 0.01 until it decides to make it slightly larger again and (apparently) switches to a higher order method that requires more linear solves per time step but allows for a larger time step closer to our initial choice 0.025 again. It does not quite hit the final time of \(T=5\) with its time step choices, but we've got only ourselves to blame for that by setting
in the input file. If hitting the end time exactly is important to us, setting the flag to true resolves this issue.
We can even reason why PETSc eventually chooses a time step of around 0.025: The boundary values undergo a complete cosine cycle within 0.5 time units; we should expect that it takes around ten or twenty time steps to resolve each period of a cycle to reasonable accuracy, and this leads to the time step choice PETSc finds.
Not all combinations of methods, time step adaptation algorithms, and other parameters are valid, but the main messages from the experiment above that you should take away are:
The program actually runs in parallel, even though we have not used that above. Specifically, if you have configured deal.II to use MPI, then you can do mpirun -np 8 ./step-86 heat_equation.prm to run the program with 8 processes.
For the program as currently written (and in debug mode), this makes little difference: It will run about twice as fast, but take about 8 times as much CPU time. That is because the problem is just so small: Generally between 1000 and 2000 degrees of freedom. A good rule of thumb is that programs really only benefit from parallel computing if you have somewhere in the range of 50,000 to 100,000 unknowns per MPI process. But it is not difficult to adapt the program at hand here to run with a much finer mesh, or perhaps in 3d, so that one is beyond that limit and sees the benefits of parallel computing.
1.9.8