This is a slightly more advanced example in which we demonstrate the use of spatial adaptivity in time-dependent problems. We discuss the implementation of the spatially adaptive version of
oomph-lib's unsteady Newton solver,
Problem::unsteady_newton_solve(...), and explain why the assignment of initial conditions should be performed by overloading the
Problem::set_initial_condition() function. We also discuss briefly how
oomph-lib's generic dump and restart functions deal with adaptive meshes.
For this purpose consider the following problem:
in the quarter-circle domain , bounded by the coordinate axes and the unit circle, subject to Neumann boundary conditions,
along the horizontal domain boundary , and to Dirichlet boundary conditions,
where the functions and are given.
We choose the functions and so that
is the exact solution.
The solution represents the "usual" tanh profile, whose steepness is controlled by the parameter so that for the solution approaches a step. The step is oriented at an angle against the axis and its position varies periodically. The parameter controls the amplitude of the step's lateral displacement, while determines the rate at which its position changes. For , the step remains stationary for most of the period and then translates rapidly parallel to the axis, making this a very challenging problem.
The figure below shows a snapshot of the animated solution, obtained from the spatially adaptive simulation discussed below, for the parameter values
The mesh adaptation in response to the translation of the step can be seen more clearly in this contour plot, taken from another animation of the solution.
Enabling spatial adaptivity in time-dependent problems involves essentially the same steps as for steady problems:
ErrorEstimatorobject must be created and passed to the mesh.
Problem::actions_after_adapt()may be overloaded to perform any actions that are required before or after the mesh adaptation, such as the deletion or recreation of any
FaceElementsthat apply flux boundary conditions.
Once these steps have been performed, a spatially adaptive solution can be computed with a three-argument version of
oomph-lib's unsteady Newton solver
The arguments to this function are as follows:
dtspecifies the (fixed) timestep.
max_adaptspecifies the maximum number of spatial adaptations allowed.
firstindicates if the first timestep is performed. This argument is required to allow the automatic re-assignment of the initial conditions following any mesh adaptations during the computation of the first timestep.
Given these arguments, the unsteady Newton solver solves the non-linear system of spatially and temporally discretised equations to advance the solution from time to . Once the solution at time has been obtained, error estimates are computed for all elements. If any elemental error estimates are outside the target range, the solution is rejected and the mesh is adapted. In the course of mesh adaptation the existing solution (the nodal values and the history values) at time are interpolated onto the new mesh before recomputing the solution. This process is repeated until the error estimates are within the target range, or until the maximum number of adaptations, specified by the parameter
max_adapt, is exceeded, just as in the steady case.
Here is an illustration of the procedure for a 1D problem:
This procedure is the obvious generalisation of the procedure for steady problems. However, in time-dependent problems two additional issues arise:
max_adapt> 1) is therefore limited by the fact that mesh refinement cannot improve their accuracy – the history values are always given by the (possibly poor) approximations obtained by interpolation from the coarser mesh employed at the previous timestep. We therefore recommend limiting the number of spatial adaptations to
max_adapt= 1. We stress that, in practice, this is not a serious restriction because the time-integration procedure will only provide (temporally) accurate results if the timestep
dtis so small that the solution at time only differs slightly from that at time . One level of mesh adaptation per timestep should therefore be sufficient to adapt the mesh in response to these changes.
max_adaptcan (and should!) be specified when the first timestep is computed. The unsteady Newton solver
Problem::unsteady_newton_solve(...) performs the revised procedure if the boolean argument
true. In that case, the values and history values on the adapted mesh are (re-)assigned by calling the function
Problembase class. You should overload it in your derived
Problemto ensure that your specific initial conditions are assigned by the mesh adaptation procedures. [In fact, the function
Problem::set_initial_condition()is not quite empty – not re-setting the initial condition when performing mesh adaptations during the first timestep of a time-dependent simulation seems "so wrong" that the function issues a warning message. Although the overloading of this function is not strictly necessary if the initial conditions can be represented exactly by the interpolation from the coarse mesh onto the fine mesh, we consider it good practice to do so, for reasons discussed in another tutorial.]
Equipped with this background information, the driver code for our example problem is easy to understand, if somewhat lengthy. [Using an example with Dirichlet boundary conditions along the entire domain boundary would have shortened the code significantly but we deliberately chose an example with Neumann boundary conditions to demonstrate that the functions
Problem::actions_after_adapt() may be used exactly as in the steady computations.] We will not discuss the methodology for applying flux-type boundary conditions in problems with spatial adaptivity in detail, but refer to the discussion provided in the earlier steady example.
Overall, the code is a straightforward combination of the driver code for the steady Poisson problem with flux boundary conditions and spatial adaptivity and the driver code for the unsteady heat equation without spatial adaptivity.
As usual, we store the problem parameters in a namespace,
TanhSolnForUnsteadyHeat, in which we also specify the source function, the prescribed flux along the Neumann boundary and the exact solution.
As discussed elsewhere,
oomph-lib's mesh adaptation procedures require curvilinear domain boundaries to be represented by
GeomObjects which describe the object's shape via their member function
GeomObject::position(...). This function exists in two versions:
GeomObject::position(xi,r)computes the position vector,
r, to the point on/in the
GeomObject, parametrised by the vector of intrinsic coordinates,
unsigned, computes the position vector at the
t- th previous timestep.
In the current problem, the domain boundary is stationary, therefore the steady and unsteady versions of the function are identical. Here is the complete source code for the
MyUnitCircle object which we will use to represent the curvilinear domain boundary:
As before, we use command line arguments to (optionally) specify a restart file. We store the command line arguments in the namespace
CommandLineArgs and build the
Problem object, passing the pointer to the source function. Next we specify the time-interval for the simulation and set the error targets for the spatial adaptation.
We create and initialise the boolean flag that indicates if the first timestep is computed, and choose a large initial value for the number of permitted mesh adaptations. We then assign the initial conditions on the coarse initial mesh and retrieve the timestep (chosen when the initial conditions are assigned in
set_initial_condition() ) from the problem's
If the simulation has been restarted, the first timestep is not the step at which the initial condition has to be assigned, therefore we reset the
max_adapt parameters to their appropriate values. If the run is not restarted, the problem will have been built with a very coarse initial mesh (comprising just three elements). We don't need an error estimator to tell us that this is too coarse to represent the solution accurately and apply two levels of uniform refinement before solving the problem. Note that we refine the entire problem, not just the mesh to ensure that
Problem::actions_after_adapt() are executed and the equation numbering scheme is re-generated.
Problem::refine_uniformly() also interpolates the solution from the coarse initial mesh onto the refined mesh but, as discussed above, this will lead to a very poor representation of the initial condition. Therefore we re-assign the initial condition on the refined mesh and document the finite-element representation of the initial condition.
The time-stepping loop itself is very similar to that used in the example without spatial adaptivity. Here we call the three-argument version of the unsteady Newton solver
Problem::unsteady_newton_solve(...) and re-set the parameters
first to their appropriate values once the first step has been performed.
The problem constructor combines the constructors of the steady and unsteady problems. We start by creating a
DocInfo object to control the output, set the parameters for the exact solution and create the
We create the
GeomObject that describes the curvilinear domain boundary and pass it to the mesh constructor:
Next, we create the surface mesh that contains the prescribed flux elements and combine the two submeshes to the
Problem's global mesh. We create an instance of the
Z2ErrorEstimator and pass it to the bulk mesh.
We pin the nodal values on the Dirichlet boundaries and select the central node in the unrefined three-element mesh as the control node at which the solution is documented in the trace file.
Finally, we complete the build of all elements by passing the relevant function pointers to the elements, and assign the equation numbers.
The remaining member functions
are identical (or at least extremely similar) to those in previous examples, so we do not list them here. You can examine the functions in detail in the source code two_d_unsteady_heat_adapt.cc.
It is worth examining the dump and restart functions, however, as they demonstrate that the generic versions defined in the
Problem base class can also deal with adaptive problems – a non-trivial task!
Since details of their implementation are hidden from the user, we briefly comment on the various tasks performed by these functions. The main task of the
Problem::read(...) function is to read values (and history values) of all
Data objects from a file and to assign these values to the appropriate
Node) objects in the
Problem. This assumes that the
Data objects have been created, and that the
Problem's various pointer-based lookup schemes access them in the order they were in when the
Problem was dumped to the restart file. In a non-adaptive computation, the number of elements and the number of
Data objects remain constant throughout the simulation and the
Problem::read(...) function can be called as soon as the
Problem has been built – usually by its constructor. (The
Problem constructor always builds and enumerates its constituent objects in the same order.)
In a simulation with spatial adaptivity the number of elements,
Data objects varies throughout the computation. It is therefore necessary to re-generate the
Problem's refinement pattern before the
Data values can be read from the restart file. This is achieved (internally) by calling the function
RefineableMesh::dump_refinement(...) for all refineable meshes before the
Data is dumped. This function writes the
Mesh's refinement pattern to the restart file, using a format that can be read by the corresponding member function
RefineableMesh::refine(...) which adapts an unrefined mesh so that its topology and the order of its
Nodes and elements is recreated.
The plots below show the time history of various parameters.
The plots illustrate clearly how the mesh is adapted as the step moves through the domain – the peaks in the number of refined/unrefined elements per timestep coincide with the periods during which the step moves very rapidly. The increase in the error during these phases is mainly due to the temporal error – the animation shows that the computed solution lags behind the exact one. We will address this by adding adaptive time-stepping in another example.
RESLT/soln0.dat) against that obtained when the re-assignment of the initial conditions after the two calls to
mainfunction is suppressed.
mainfunction and compare the computed results against those obtained with the correct procedure.
A pdf version of this document is available.