tutorial from Nvidia GTC

My talk is attached in *pdf. Overall the conference felt like a rock’n'roll show, my first non-academic meeting of this nature. Thoroughly enjoyable although a little difficult to tell what people are doing since it’s mostly proprietary work. The announcement of fermi was certainly exciting; ~8x times the double prec performance, ECC and many more cores/memory…

Talk outline:

• Company Background

• CUDA accelerates Geophysics:
  • Data Processing w/ Linear Algebra

  • Sparse Matrices
  • Seismic Imaging:Kirchhoff Time Migration
    • Frequency Filtering

    • Traveltime/Weight Calculations
    • GPU Migration
• Summary

nvidiaSummit0909

Split-Step Algorithm

This is a fast method for the solution of parabolic PDE’s, relying on the FFT implementation of the Fourier Transform to speed things up. Here it’s applied to an acoustic/seismic problem, following the development of Kuperman/Jackson in “Imaging of Complex Media with Acoustic and Seismic Waves”. Consider first the Helmholtz equation for a point source at r’,z’:

where G is the Green’s function, K the wavenumber (function of the frequency omega and sound speed c). Assuming azimuthal symmetry, G may be expressed as a product of two functions:

and similarly K, now a product of (constant) K0 and index of refraction n. Substituting into the Helmholtz equation gives two PDE’s:

The first PDE has Bessel functions as solutions; taking the assymptotic outgoing Hankel function solution and substituting it into the second, with the narrow angle approximation (second derivative of psi with respect to r much smaller than first derivative wrt r), one finds for the second (parabolic) PDE:

where chi is the fourier transform of psi. Assuming the variation of n is insignificant, in the wave-space domain the PDE and solution are:

Finally, the inverse FT gives the field as a function of depth (Delta r = r-r0, r0 the boundary value):

The following is a little hack in octave, to demonstrate the solution method:

function [Psi,s]=test_ssa(N,Dr,P0,n,K0)

% to test split step algo on parabolic PDE:
%
% \frac{\partial^2 \psi}{\partial r^2}
% +2iK0\frac{\partial \psi}{\partial r}
% +K0^2(n^2-1)\psi = 0
%
% K0*n^2 = (const.) wave number
% N=marching steps in field
% P0=boundary field data
% Dr=field step size

% wjb 0609

%initialize
m=length(P0); Psi=zeros(N,m); Psi(1,:)=P0;

%Nyquist thm, s (wave space) range

s=fftshift(-m/2:m/2-1)./Dr;

for k=2:N

Psi(k,:)=exp((i*K0/2)*(n^2-1)*Dr).*ifft(exp(-(i*Dr*s.^2)./(2*K0)).*fft(Psi(k-1,:)));

end

Custom Electric Potential

I would like an electric potential with specific qualities, namely, quadratic in two dimensions. On the one hand, I could try and do a little thought experiment, as to which boundary gives the desired properties. Realistically though on the length scale I need it, the odds of producing the boundary are zilch. Much easier in terms of the multipole expansion. With this in mind, and knowing that copper wire along a cylinder makes for an easy setup, and that at least four terms are required for a quadratic potential, the following should work: V(x=0,z=a)=k, V(x=0,z=-a)=k, V(x=-a,z=0)=-k, V(x=a,z=0)=-k, with each wire along y. Assuming long wires, the potential for all four in the x/z plane is:

To confirm this has the desired behavior near the origin (and hopefully over some reasonable distance), I turn to the Taylor multi-variable expansion. Example partial derivatives:

Substituting into the expansion:

does indeed reveal that V ~ (2/{a*a})[z*z-x*x], as desired:

MQMAS sim paper in press

Abstract: The majority of nuclei available for study in solid state Nuclear Magnetic Resonance have half-integer spin I > 1/2, with corresponding electric quadrupole moment. As such, they may couple with a surrounding electric field gradient. This effect introduces anisotropic line broadening to spectra, arising from distinct chemical species within polycrystalline solids. In Multiple Quantum Magic Angle Spinning (MQMAS) experiments, a second frequency dimension is created, devoid of quadrupolar anisotropy. As a result, the center of gravity of peaks in the high resolution dimension are functions of isotropic quadrupole and chemical shifts alone. However, for complex materials, these parameters take on a stochastic nature due in turn to structural and chemical disorder. Lineshapes may still overlap in the isotropic dimension, complicating the task of assignment and interpretation. A distributed computational approach is presented here which permits simulation of the MQMAS spectrum, generated by random variates from model distributions of isotropic chemical and quadrupole shifts. Owing to the non-convex nature of the least squared cost function between experimental and simulated spectra, simulated annealing is used to optimize the simulation parameters. In this manner, local chemical environments for disordered materials may be characterized, and via a re-sampling approach, error estimates for parameters produced.

Key words: Nuclear Magnetic Resonance, Multiple Quantum Magic Angle Spinning, OpenMP, Sobol sequence, quasi-random numbers, simulated annealing, distribution functions, quadrupole interaction.

Padé Approximants

You may have a Taylor series which is slowly convergent or downright divergent, in which case you might want to try a Padé approximant. This is essentially a method which approximates a function in terms of the ratio of two polynomials, of orders N,M. A particularly useful representation is in continued fractions; the (N+1)th member of the Padé sequence:

for J >= 0 is given by:

Bender & Orszag in their book give an algorithm for constants c, I’ve worked out the first few by hand for you. There’s no such thing as a free lunch; you might find a converged value for your function, but it’s at the expense of nasty analysis in finding progressively higher orders for c…