Guiding AMR with adjoint flagging

A new approach to flagging cells for refinement was introduced in Clawpack 5.6.0 – using the solution to an adjoint problem to determine what cells in the forward solution should be refined because these cells may have an impact on the some specified functional of interest. This approach currently only works for autonomous linear problems, in which case the adjoint problem needs to be solved only once, and shifted versions of the adjoint solution can be used at any time that flagging is performed. The adjoint problem is solved first and snapshots of the adjoint are saved. These are read in at the start of the forward solution, and space-time interpolation used as needed at each regridding time.

The general approach is described in:

See Using adjoint flagging in GeoClaw for discussion of the GeoClaw version.

Adjoint flagging is appropriate when you are not interested in computing an accurate solution over the entire space-time domain, but rather are interested only in some linear functional applied to the solution at each time (or at a single time, or range of times). In one space dimension this functional has the form

\[J(t) = \int_a^b \phi(x)^T q(x,t)\,dx\]

where \(a\leq x \leq b\) is the full computational domain and \(\phi(x)\) is specified by the user as initial data for the adjoint problem that is solved backward in time. For example, if the solution is of interest only over a small range of \(x\) values, say \(x_1 \leq x \leq x_2\), then \(\phi(x)\) might be a box function with value 1 in this interval and 0 elsewhere, or \(\phi(x)\) could be a sharply peaked Gaussian about one location of interest.

In order to calculate an accurate solution near the location of interest at the final time \(T\) it may be necessary to refine the solution at other places at earlier times. The adjoint helps to identify where refinement is needed. The adjoint equation is first solved backward in time from the final time \(T\) with initial data \(\hat q(x,T) = \phi(x)\) given by the functional. The waves propagating backward from time \(T\) to some regridding time \(t_r\) in the adjoint solution identify which waves in the forward solution at time \(t_r\) will reach the location of interest at time \(T\).

Some examples for AMRClaw are available in

  • $CLAW/amrclaw/examples/acoustics_1d_adjoint

  • $CLAW/amrclaw/examples/acoustics_2d_adjoint

In each case the main directory has a subdirectory named adjoint that contains the code that must be run first in order to compute and save snapshots of the adjoint solution.

The adjoint/Makefile must point to an appropriate Riemann solver for the adjoint problem, which is a linear hyperbolic PDE with coefficient matrices that are transposes of the coefficient matrices appearing in the forward problem.

For variable-coefficient problems it is important to note that if the forward problem is in conservation form then the adjoint is not, and vice versa. For example, in one space dimension, if the forward problem is \(q_t + A(x)q_x = 0\), then the adjoint is \(\hat q_t + (A(x)^T \hat q)_x = 0\). On the other hand if the forward problem is \(q_t + (A(x)q)_x = 0\), then the adjoint is \(\hat q_t + A(x)^T \hat q_x = 0\).

Note that the eigenvalues of \(A\) are unchanged upon transforming but the left eigenvectors of \(A\) are now the right eigenvectors of \(A^T\), and these must be used in the adjoint Riemann solver. See for example $CLAW/riemann/src/rp1_acoustics_variable_adjoint.f90, used for the example in $CLAW/amrclaw/examples/acoustics_1d_adjoint/adjoint.

Boundary conditions conditions may also need to be adjusted in going from the forward to adjoint equation. The guiding principle is that boundary conditions must vanish during the integration by parts that is used to define the adjoint PDE, as described in more detail in the references.

The functional of interest is defined in the adjoint/qinit.f file that specifies “initial” conditions for the adjoint problem.

The adjoint/setrun.py file specifies the final time \(T\) (as clawdata.tfinal) and the output interval via clawdata.num_output_times, as usual. You should specify \(T\) at least as large as the final time of interest in the forward problem, and frequent enough snapshots that interpolation between them is reasonable.

You should set

clawdata.output_format = 'binary'

so that output is in binary format, since the code that reads in these snapshots in solving the forward problem assumes this format.

After solving the adjoint equation by running the code in the adjoint subdirectory in the usual manner (e.g. make .output), the code in the main directory can now be used to solve the forward problem, with the adjoint snapshots used to guide AMR.

Starting in v5.6.0 a new attribute of clawutil.data.ClawRunData is available named adjointdata. This ia an object of class amrclaw.data.AdjointData and has several attribures that should be set. For example, in $CLAW/amrclaw/examples/acoustics_1d_adjoint they are set as follows:

#------------------------------------------------------------------
# Adjoint specific data:
#------------------------------------------------------------------
# Also need to set flagging method and appropriate tolerances above

adjointdata = rundata.adjointdata
adjointdata.use_adjoint = True

# location of adjoint solution, must first be created:
adjointdata.adjoint_outdir = os.path.abspath('adjoint/_output')

# time period of interest:
adjointdata.t1 = rundata.clawdata.t0
adjointdata.t2 = rundata.clawdata.tfinal

if adjointdata.use_adjoint:
    # need an additional aux variable for inner product:
    rundata.amrdata.aux_type.append('center')
    rundata.clawdata.num_aux = len(rundata.amrdata.aux_type)
    adjointdata.innerprod_index = len(rundata.amrdata.aux_type)

The times adjointdata.t1 and adjointdata.t2 determine the time interval over which the functional is of interest, and adjointdata.adjoint_outdir specifies where the adjoint snapshots are found.

The flagging method and tolerances are set using, e.g.:

# set tolerances appropriate for adjoint flagging:

# Flag for refinement based on Richardson error estimater:
amrdata.flag_richardson = False
amrdata.flag_richardson_tol = 1e-5

# Flag for refinement using routine flag2refine:
amrdata.flag2refine = True
amrdata.flag2refine_tol = 0.01

If amrdata.flag_richardson is True then we attempt to use estimates of the one-step error generated by Richardson extrapolation together with the adjoint to perform flagging. This is still experimental. (Describe in more detail).

Otherwise it is simply inner products of the forward and adjoint solutions that are computed and a cell is flagged for refinement in cells where the magnitude of the inner project is greater than amrdata.flag2refine_tol.

Using adjoint flagging in GeoClaw

The references above contain tsunami modeling examples, as does the paper

An example can be found in

  • $CLAW/geoclaw/examples/tsunami/chile2010_adjoint

Note that GeoClaw solves the nonlinear shallow water equations while the adjoint as implemented in GeoClaw is only suitable for linear problems. To date the adjoint has only been used to guide refinement for waves propagating across the ocean as a way to identify which waves will reach a target location of interest (possibly after multiple reflections). In the deep ocean the tsunami amplitude is very small compared to the water depth and so GeoClaw is essentially solving the linear shallow water equations, linearized about the ocean at rest. Hence the adjoint problem is also solved about the ocean at rest and the adjoint equations take essentially the same form as the forward equations. The adjoint Riemann solver can be found in

  • $CLAW/riemann/src/rpn2_geoclaw_adjoint_qwave.f

  • $CLAW/riemann/src/rpt2_geoclaw_adjoint_qwave.f

Note that since in the forward problem the adjoint equation is solved using a f-wave formulation, the adjoint problem is a variable-coefficient problem in non-conservation form and is solved using the q-wave formulation in which jumps in the the solution vector are split into eigenvectors, rather than jumps in the flux. See the comments in the rpn2 solver for more details.