next up previous [pdf]

Next: Single Shot Examples Up: IWAVE Structure and Basic Previous: IWaveInfo

Case: Constant-density Acoustics

A simple example illustrating the framework described above is the IWAVE implementation of the constant-density acoustic wave equation with Dirichlet (pressure-free) boundary conditions, connecting the acoustic potential field $ u({\bf x},t)$ and a right-hand side $ f({\bf x},t)$ representing a source of mechanical energy, defined in a spatial domain $ \Omega$ over a suitable time interval,

$\displaystyle \frac{\partial^2 u}{\partial t^2} - c^2 \nabla^2 u = f,
$

$\displaystyle u({\bf x},0) = u_0({\bf x}),  \frac{\partial u}{\partial t}({\bf x},0) =
v({\bf x})   {\bf x}\in \Omega,
$

$\displaystyle u({\bf x},t) = 0,  {\bf x}\in \partial \Omega.
$

The examples to be discussed use the centered difference approximation (Kelly et al., 1976)

$\displaystyle u^{n+1} = 2 u^{n} - u^{n-1} + \Delta t^2 c^2 Lu^n + \Delta t^2 f^n
$

in which $ L$ is a regular grid difference approximation to the Laplacian, and $ u^n$ represents the array of acoustic potential samples for time $ n\Delta t$ . The choice of $ L$ used below is a sum of centered coordinate second difference operators of order $ 2k$ , $ k=1, 2, 4,...$ resulting in a scheme of formal order $ 2$ in time and $ 2k$ in space. Lax-Wendroff extension to higher order time approximation fits this pattern also.

Since each array element in $ u^{n-1}$ appears exactly once in a loop through the array, it is possible to store only the two arrays for time indices $ n-1$ and $ n$ , represented by RARRs up and uc respectively, and store $ c^2$ in the RARR csq. The the three-level scheme above becomes

up = 2 * uc - up + dt2 * csq .* L uc + dt2 * f^n
[swap up, uc]
With the type of discrete finite difference Laplacian described above, the grids for uc, up, and csq are all primal, and commensurable. Since csq must exist along with all of its metadata (its grid information, basically) in the scope of the simulation, it is natural to read the primal grid geometry from it. Thus an appropriate iwave_fields array for acoustic constant density modeling is
FIELD IWaveInfo::iwave_fields[]
= {
  {"csq",    0,    0,  {0, 0, 0}},
  {"uc",     1,    0,  {0, 0, 0}},
  {"up",     1,    0,  {0, 0, 0}},
  {"",       0,    0,  {0, 0, 0}}
};
The last line functions the same way as the traling NUL for C strings, that is, to signal the end of the structure.

Inspection of the pseudo-code above reveals that up and csq need to be available at precisely the same gridpoints, whereas uc must store additional gridpoints around the boundary (``halo'' or ``ghost'' points) in order that the Laplacian can be built at the csq gridpoints (the ``physical'' grid). So the memory allocations for uc and up differ, and the ``swap'' mentioned in the algorithm exchanges only values of the two fields at physical grid points. The algorithm must be completed with a boundary loop which updates the non-physical gridpoints of uc, for example with odd reflection implementing a Dirichlet condition.

The design described in the preceding paragraphs is realized in the IWAVE acoustic constant density package, RSFSRC/trip/iwave/acd. The standalone executable implementing the various options provided by IWAVE is also called acd. It can be built as part of a Madagascar top-down build, in which case it shows up as RSFROOT/bin/sfacd and can be referenced as acd in Madagascar Flows, or standalone via invocation of scons in either RSFSRC/trip or RSFSRC/trip/iwave. In the latter case, the dependency RSFSRC/trip/rvl must be built first. The standalone build has the virtue of permitting local control of build environment. The RSFSRC/trip/admin directory includes a number of example configuration scripts for build options - to use, copy one of these, or create a similar file, as config.py in the root build directory.

The standalone build results in the executable command

  RSFSRC/trip/iwave/acd/main/acd.x.
In the remainder of this paper, I will refer to this command as acd.x. The SConstruct file in the project subdirectory of the paper directory is also configured to use this standalone-built command.

The IWAVE acoustic constant density implementation includes code for the acoustic simulator and its derivatives (with respect to velocity-squared) of orders 1 and 2, and their adjoints, built from a simple numerical kernel (or set of kernels) of truncation orders in space 2, 4, and 8, and truncation order 2 in time. My research group has used the Tapenade (Hascoët and Pascual, 2013) automatic differentiation package to produce the code for derivatives and adjoints. For example, the signature of an implementation of the $ (2,2k)$ scheme for 3D is

acd_3d_[2k](float *** uc3, 
            float *** up3, 
            float *** csq3, 
            int * s, 
            int * e,
            ...,
            int * lbc,
            int * rbc);
in which uc3 is the 3-dimensional array view of the array uc, and so on; s and e are the vectors of start and end indices for the loop over gridpoints; and ... stands in for a list of difference formula coefficients, the the number and value of which depend on the order ($ 2k$ ). The integer arrays lbc and rbc flag whether the left and right boundaries of the computational domain, delimited by s and e, are external (physical) boundaries or internal boundaries. In the former case, phsical boundary conditions must be applied; these are also part of the code.

Tapenade produces similar code for the first derivative of this stencil (with respect to the uc, up, and csq arguments, with signature

void acd_3d_[2k]_d(float *** uc, 
                   float *** ucd, 
                   float *** up, 
                   float *** upd, 
                   float *** csq,
                   float *** csqd, 
                   int * s, 
                   int * e, 
                   ...,
                   int *lbc,
                   int *rbc);
in which ucd, upd, and csqd are the perturbations of the arrays without the d's.

These kernels can be folded into an obvious case list, switched by the inputs to the timestep interface described above.

Tapenade output is not entirely suitable for immediate use: some minor cleanup is necessary, and any serious optimizations (vectorization, for instance) will need to be applied in a tuning phase. However the code as it comes from the package is correct and reasonably readable, and can serve as a baseline with which to verify tuned versions.

Definition of a command based on the fields and functions described above requires one more piece of information: the connection of fields to external data sources and sinks, intermediated by i/o keywords. Many choices are possible; one reasonable choice for the standalone command option and constant density acoustics (acd.x) is:

IOKEY IWaveInfo::iwave_iokeys[]
= {
  {"csq",    0, true,  true },
  {"data",   1, false, true },
  {"source", 1, true,  false},
  {"movie",  1, false, false},
  {"initc",  1, true,  false},
  {"initp",  2, true,  false},
  {"",       0, false, false}
};
Clearly the velocity (or rather velocity-squared) must be made available. Two outputs from uc are identified, "data" and "movie": while nothing about the specs demands this usage, the first is intended for trace output, the second for time slices, as the keywords choices are intended to suggest. Since the precise mechanism of I/O is inherent in the data unit (file structure, for instance) rather than the directed by the code, in fact these mnemonic suggestions could be ignored, and "data" used to store a movie, for example. However it is an intended use case that movies might be generated at a byproduct of trace generation, so two output slots are provided. Similarly, several input keywords suggest a right-hand side input (time dependent force divergence traces) ("source") and Cauchy data ("initc", "initp").

Note that the discrete Cauchy data represent pressure at two successive time levels, whereas the natural Cauchy data for the wave equation would provide presure and its time derivative. An application accepting this natural Cauchy data would need to pre-process it into discrete Cauchy data as indicated above. I have elected to ``un-bundle'' this type of pre-processing, that is, it is not included in the IWAVE code itself. Similarly, the natural SEGY representation of the RHS source traces needs to be pre-processed to code the source positions as receiver coordinates, as reviewed below.



Subsections
next up previous [pdf]

Next: Single Shot Examples Up: IWAVE Structure and Basic Previous: IWaveInfo

2015-04-20