Bill Spotz Scott Collis Sandia National Laboratories HIGH

Bill Spotz, Scott Collis Sandia National Laboratories HIGH ORDER FINITE DIFFERENCE METHODS FORWARD AND INVERSE ACOUSTIC PROBLEMS Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the United States Department of Energy’s National Nuclear Security Administration under contract DE-AC 04 -94 AL 85000.

Background n n Sandia Computational Research Center has a framework for solving PDEs: DGM (Discontinuous Galerkin Method) • Low and high-order unstructured meshes, local p-refinement • Physics independent: compressible and incompressible Euler, Navier. Stokes, advection-diffusion, Burgers, Darcy, Helmholtz, Poisson, wave equations, acoustic, etc. • Adjoint-based optimization • Object-oriented, component-based software design; massive parallelism Much of DGM is independent of discontinuous Galerkin: • Time integration, optimization, inversion, etc. We would like to compare DG results and performance against traditional finite difference methods To support this, we are developing the Discretization Tool Kit (DTK) by • Extracting concrete components: time integrators, optimizers, etc. • Abstracting general components: fields, meshes, adjoints, controls, etc.

Motivation n A device (e. g. a speaker) is used to generate waves in an acoustic medium (e. g. air) Receivers are set up to measure the resulting wave field Sandia is developing an application built on top of DGM (DTK), to solve inverse problems such as: • • Determine the signature of the initial disturbance (source inversion) Determine the material properties of the media

Acoustic Equations Symbol Definition p Pressure v Velocity t Time r Density Forward Problem Assume r, k, x, f, f are known, solve for p(t), v(t). k Bulk modulus x Spatial kernel Inverse Problem Measure p(t), v(t) at discrete locations, solve for (one or more of) r, k, f, f. f Source f External forcings High-level algorithm: • Does computed solution “match” recorded data? • If not, adjust control parameters. Repeat using adjoint-based optimization algorithms

DTK (Planned)

DGM & FDM Forward Solver Infrastructure

Finite Difference Grid 3 D Stencil where v = (u, v, w) Stagger? Variable x y z p Y Y Y u N Y Y v Y N Y w Y Y N

High Order, Staggered, Finite Differences We will derive approximations for similarly , but all derivatives will behave Taylor series about grid point i + ½: Subtract (eliminating ui+1/2 and even m summation terms):

High Order, Staggered, Finite Differences Generalize for symmetric points about i+1/2 at distance n grid lengths: A high-order formula will be a weighted sum of these formulae: Here, N is the (half) size of the symmetric finite difference stencil, and , n = 1, 2, … N, are the stencil coefficients.

High Order, Staggered, Finite Differences Let’s isolate the first derivative (m = 1) and truncate the inner summation at m = N (because we know this will give us a square linear system): To obtain an O(Dx 2 N) approximation to for m = 2, 3, … N. at i+1/2, set

High Order, Staggered, Finite Differences This reduces to simple linear systems for solving the stencil coefficients. For example, N = 1: N = 2: N = 3:

High Order, Staggered, Finite Differences N Order 1 2 4 9 8 -1 24 3 6 75 64 -25 384 3 640 4 8 1, 225 1, 024 -245 3, 072 49 5, 120 -5 7, 168 5 10 19, 854 16, 384 -735 8, 192 567 40, 960 -405 229, 376 35 294, 912 6 12 160, 083 131, 072 -5, 445 1, 835, 008 847 2, 359, 296 -12, 705 22, 869 131, 072 1, 310, 720 -63 2, 883, 584

Source Inversion Problem Verification Procedure: • Pick a mathematical form for the source wavelet • Run the forward model and record the results at the receiver locations • Using only the receiver results, run the inverse problem to compute an approximate source wavelet • Compare the computed source to the known source. Source of choice: Ricker wavelet where fp is the peak frequency and t 0 is the time delay

Source Inversion Formulation PDE-constrained optimization problem: Subject to acoustic wave equations in conservative form: Symbol Definition Measured pressure F Space of admissible sources T Time horizon W Computational domain Nr Number of receivers f Source (speaker) xr Spatial receiver kernel

DTK Optimization Problem

Implementing Acoustic Inverse

Inverting for the Source It Obj Fn 0 0. 0166021 1 0. 000628603 2 7. 28457 e-05 3 1. 47980 e-05 4 3. 80999 e-06 5 1. 13207 e-06

Summary n n n The Discontinuous Galerkin Method (DGM) framework has been generalized to a Discretization Tool Kit (DTK) • DGM still exists as a discretization method under DTK • New discretization methods can be added … starting with the Finite Difference Method (FDM) FDM has been added to the DTK suite • Automated parallel decomposition and halo updates • Flexible support for staggered fields • Support for 2, 4, 6, … 12 orders of accuracy • Forward solutions handled by DTK: : Time. Integrate FDM Source Inversion has been demonstrated • Results are indistinguishable from DGM (simple domain and mesh) • Framework in place for more complex inversions: w w Media inversion, source+media inversion Multiphysics