libEMM: A fictious wave domain 3D CSEM modelling library bridging sequential and parallel GPU implementation

The CSEM physics is governed by the diffusive Maxwell equation (consisting of Faraday’s law and Ampère’s law), which reads in the frequency domain

where $\mathbf{E}=(E_{x},E_{y},E_{z})^{T}$ and $\mathbf{H}=(H_{x},H_{y},H_{z})^{T}$ are the electrical and magnetic fields; $\mathbf{J}=(J_{x},J_{y},J_{z})^{T}$ and $\mathbf{M}=(M_{x},M_{y},M_{z})^{T}$ are the electrical and magnetic sources. Here, we slightly abuse the same notation to denote the same quantity when switching between time domain and frequency domain, where the convention of Fourier transform $\partial_{t}\leftrightarrow-\mathrm{i}\omega$ has been adopted. The magnetic permeability is $\mu$. The conductivity is a symmetric $3\times 3$ tensor:

in which $\sigma_{ij}=\sigma_{ji}$, $i,j\in\{x,y,z\}$. The isotropic medium means only the diagonal elements of the conductivity tensor are non-zeros and the same in all directions: $\sigma_{xx}=\sigma_{yy}=\sigma_{zz}$; $\sigma_{ij}=0,i\neq j$. The vertical transverse isotropic (VTI) medium still has only the diagonal elements, that is, $\sigma_{h}:=\sigma_{xx}=\sigma_{yy},\quad\sigma_{v}=\sigma_{zz}$, where $\sigma_{h}$ and $\sigma_{v}$ stand for horizontal conductivity and vertical conductivity, respectively. The resistivity is also frequently used in CSEM system which is the inverse of the conductivity ($\rho_{ij}=1/\sigma_{ij}$).

By defining a fictitious di-electrical permittivity $\epsilon$ such that $\sigma=2\omega_{0}\varepsilon$, and multiplying the second equation with $\sqrt{-\mathrm{i}\omega/2\omega_{0}}$, (Mittet,, 2010) transformed the equation into

where we introduce a prime to define the resulting wave domain fields via

Equation (3) may be transformed into time domain as

From the electromagnetic fields in the wave domain, the frequency-domain fields can be computed on the fly during modelling using the fictitious wave transformation

where $i=\{x,y,z\}$ are component indices of the electromagnetic fields; $T_{\max}$ is the final time until the electromagnetic fields reach their steady state.

The Green’s function for point source is then obtained by normalizing the electromagnetic field with the source current. In the absence of magnetic source ($\mathbf{M}=0$), the Green’s functions due to the electrical current $J_{j}(\mathbf{x}_{s},\omega)$ are

where $G_{ij}^{EE}(\mathbf{x},\omega;\mathbf{x}_{s})$ and $G_{ij}^{HE}(\mathbf{x},\omega;\mathbf{x}_{s})$ stand for the $i$th component of the electrical and magnetic Green’s function for angular frequency $\omega$ at the spatial location $\mathbf{x}$.

Different from the attenuative formulation after (Maaø,, 2007), equation (3) is a pure wave domain equation, which can be discretized using leap-frog staggered grid finite difference method on the staggered (Yee) grid:

where $\Delta t$ is the temporal sampling, $\xi_{ijk}$ is the Levi-Civita symbol. The positioning of each field on the nonuniform staggered grid is similar to the staggered FDTD on a uniform grid. The major difference between FDTD implementations on a uniform and a nonuniform grids lies in the discretization of these spatial derivatives ($\partial_{i}E^{\prime}_{k}$ and $\partial_{k}H^{\prime}_{l}$)

In order to compute the electromagnetic field as well as its derivatives with arbitrary grid spacing, we have to do polynomial interpolation using a number of knots $x_{0},x_{1},\cdots,x_{n}$. According to Taylor expansion, we have

where $i=0,1,\cdots,n$. Let us define $a_{i}(x):=f^{(i)}(x)/i!$ and use first $n+1$ distinct nodes $x_{0},x_{1},\cdots,x_{n}$ to build a matrix system

where $\mathbf{V}^{T}(x_{0}-x,\cdots,x_{n}-x)$ is the transpose of a Vandermonde matrix determined by $x_{0}-x,\cdots,x_{n}-x$. The above expression implies that the function $f(x)$ and its derivatives up to $n$-th oder at arbitrary location $x$ can be found by inverting the Vandermonde matrix $\mathbf{a}=[\mathbf{V}^{T}]^{-1}\mathbf{f}$. When the order of the Vandermonde matrix increases, the inversion of Vandermonde matrix becomes highly ill-conditioned. Fortunately, there exists accurate and efficient algorithm (Björck and Pereyra,, 1970) to achieve this challenging job, by exploring the equivalence between Vandermonde matrix inversion and the Lagrange polynomial interpolation. The following code snippet implements it following (Golub,, 1996, Algorithm 4.6.1).

/* solve the vandermonde system V^T(x0,x1,...,xn) a = f
 * where the n+1 points are prescribed by vector x=[x0,x1,...,xn];
 * the solution is stored in vector a.
 */
void vandermonde(int n, float *x, float *a, float *f)
{
  int i,k;

  /* calculate Newton representation of the interpolating polynomial */
  for(i=0; i<=n; ++i) a[i] = f[i];

  for(k=0; k<n; k++){
    for(i=n; i>k; i--){
      a[i] = (a[i]-a[i-1])/(x[i]-x[i-k-1]);
    }
  }

  for(k=n-1; k>=0; k--){
    for(i=k; i<n; i++){
      a[i] -= a[i+1]*x[k];
    }
  }
}

Let the $i$-th row, $j$-th column of the inverse matrix $[\mathbf{V}^{T}]^{-1}$ be $w_{ij}$, i.e. $([\mathbf{V}^{T}]^{-1})_{ij}=w_{ij},i,j=0,\cdots,n$. It follows that

There are a number of situations that requires the use of interpolation weights repeatedly at every time step:

1. 1.

   The same finite difference coefficients will be used at the same grid location to do numerical modelling;

2. 2.

   The same interpolation weights may be used to record fields at the same location in different time steps.

3. 3.

   The injection of the source time function needs the same interpolation weights to be applied.

The last two cases are due to the fact that the sources and the receivers may reside in the arbitrary location, not necessarily on the grid. As a consequence, it is necessary to tabulate these weights instead of repeating the computation according to the inversion algorithm. Since applying the algorithm in operator form does not need the matrix coefficients explicitly, it is not obvious how to obtain the interpolation weights $w_{ij}$. We present an effective way here to make it. The idea is the following: feeding the algorithm with the vector $f$ with $j$th element 1 and other elements 0, the output of the Vandermonde inversion algorithm will be the $j$th column of the matrix $[\mathbf{V}^{T}]^{-1}$:

By looping over index $j$ from 0 to $n$, we obtain all columns of $[\mathbf{V}^{T}]^{-1}$. Over the rectilinear non-uniform grid, we construct multi-dimensional finite difference stencil using the tensor product of three 1D interpolation operator.

Consider the staggered finite difference approximation of the first derivatives in $x$ direction using $2r$ non-equidistant nodes centered at $x=x_{i+1/2}$ for the magnetic component and $x=x_{i}$ for the electric component. The forward operator $D_{x}^{+}$ to discretize 1st derivative of electrical field is given by

and the backward operator $D_{x}^{-}$ to discretize 1st derivative of magnetic field is

where we have paired the differencing terms to highlight the similarity to the uniform grid staggered finite differencing. The coefficients $c_{j}^{+}$ and $c_{j}^{-}$, $j=-r+1,\cdots,r$ are the 2nd row of the inverse of the matrix $\mathbf{V}^{T}(x_{i+r-1/2}-x_{i},\cdots,x_{i-r+1/2}-x_{i})$ and $\mathbf{V}^{T}(x_{i+r}-x_{i},\cdots,x_{i-r+1}-x_{i})$ in (11). Using $2r$ nodes, we achieve higher accuracy up to $2r$-th order in space.

Denote the radius $r$ as emf->rd. Using Vandermonde matrix inversion, we compute the finite difference coefficients of $2r$th order and store it in emf->v3s (here 3 denotes z direction and s means grid staggering) to compute $\partial_{z}H_{x}$ and $\partial_{z}H_{y}$, yielding

//i1min,i1max=lower and upper bounds for i1 - x direction
//i2min,i2max=lower and upper bounds for i2 - y direction
//i3min,i3max=lower and upper bounds for i3 - z direction
for(i3=i3min; i3<=i3max; ++i3){
  for(i2=i2min; i2<=i2max; ++i2){
    for(i1=i1min; i1<=i1max; ++i1){
     ...
     if(emf->rd==1){
       D3H2 = emf->v3s[i3][0]*emf->H2[i3-1][i2][i1]
         + emf->v3s[i3][1]*emf->H2[i3][i2][i1];
       D3H1 = emf->v3s[i3][0]*emf->H1[i3-1][i2][i1]
         + emf->v3s[i3][1]*emf->H1[i3][i2][i1];
     }else if(emf->rd==2){
       D3H2 = emf->v3s[i3][0]*emf->H2[i3-2][i2][i1]
         + emf->v3s[i3][1]*emf->H2[i3-1][i2][i1]
         + emf->v3s[i3][2]*emf->H2[i3][i2][i1]
         + emf->v3s[i3][3]*emf->H2[i3+1][i2][i1];
       D3H1 = emf->v3s[i3][0]*emf->H1[i3-2][i2][i1]
         + emf->v3s[i3][1]*emf->H1[i3-1][i2][i1]
         + emf->v3s[i3][2]*emf->H1[i3][i2][i1]
         + emf->v3s[i3][3]*emf->H1[i3+1][i2][i1];
     }else if(emf->rd==3){
       D3H2 = emf->v3s[i3][0]*emf->H2[i3-3][i2][i1]
         + emf->v3s[i3][1]*emf->H2[i3-2][i2][i1]
         + emf->v3s[i3][2]*emf->H2[i3-1][i2][i1]
         + emf->v3s[i3][3]*emf->H2[i3][i2][i1]
         + emf->v3s[i3][4]*emf->H2[i3+1][i2][i1]
         + emf->v3s[i3][5]*emf->H2[i3+2][i2][i1];
       D3H1 = emf->v3s[i3][0]*emf->H1[i3-3][i2][i1]
         + emf->v3s[i3][1]*emf->H1[i3-2][i2][i1]
         + emf->v3s[i3][2]*emf->H1[i3-1][i2][i1]
         + emf->v3s[i3][3]*emf->H1[i3][i2][i1]
         + emf->v3s[i3][4]*emf->H1[i3+1][i2][i1]
         + emf->v3s[i3][5]*emf->H1[i3+2][i2][i1];
     }
     ...
    }
  }
}

The above procedure allows us to use high order finite difference scheme to accurately compute the electromagnetic fields and their derivatives, typically with arbitrary grid spacing in the rectilinear grid. This opens the door for CSEM modelling using high order FDTD on a nonuniform grid in a consistent framework.

Following the standard von Neumann analysis, the author have proved (Yang and Mittet,, 2023) that the following condition must be satisfied in the generic rectilinear nonuniform grid for stable numerical modelling

where $v_{\max}=\max\{v\}$ in which $v:=1/\sqrt{\mu\epsilon}$ is the velocity of the propagating EM field; $D_{x}^{\max}$, $D_{y}^{\max}$ and $D_{z}^{\max}$ are the maximum value of the the discretized first derivatives along $x$, $y$ and $z$ directions. In particular, we have

and similar estimations for $D_{y}^{\max}$ and $D_{z}^{\max}$. To minimize the number of time steps $N_{t}$, we would like to use the largest possible $\Delta t$ without violation of the above stability condition. The temporal sampling $\Delta t$ is therefore automatically determined based on the pre-determined grid spacing, i.e., $\Delta t=0.99/\eta$.

As shown in the Appendix, the dispersion error on the uniform grid can be easily measured. However, over the non-uniform grid, the discritized spatial derivative operators in (13) and (14) immediately make the dispersion analysis difficult (the wavenumbers are also space dependent). In view of the difficulties to practically measure the dispersion error on nonuniform grid, there is no control of the dispersion error in wave domain. Another solid reason for not doing so is that the final electromagnetic field we compute is in diffusion domain instead of wave domain. Physically speaking, the diffusive EM field is dispersive in itself, thus the physical dispersion will domainate the computed field after conversion from wave domain into diffusion domain. This means the numeric dispersion error in wave domain will eventually be annihilated by the physical dispersion. Note that in our method, the physical dispersion has been achieved analytically according to equation (4), thus free of numeric error.

At each time step, the top boundary (air-water interface in marine CSEM or air-formation interface in land CSEM) must be properly implemented with high order FDTD (Wang and Hohmann,, 1993; Oristaglio and Hohmann,, 1984). In the air, the conductivity is zero. Without electrical source, the Ampère’s law then reduces to

which leads to the wave-number domain relation

where $k_{x}$, $k_{y}$ and $k_{z}$ are wave-numbers in x-, y- and z- directions. Since $\nabla\cdot\mathbf{H}=0$, from the relation $\nabla\times\nabla\times\mathbf{H}=\nabla(\nabla\cdot\mathbf{H})-\Delta\mathbf{H}$ we obtain

which implies $k_{x}^{2}+k_{y}^{2}+k_{z}^{2}=0$, hence $k_{z}=\pm\mathrm{i}\sqrt{k_{x}^{2}+k_{y}^{2}}$. It is critical to choose a correct sign here such that the field vanishes at infinite far distance ($\mathbf{H}\rightarrow 0$ when $z\rightarrow 0$), yielding $k_{z}=-\mathrm{i}\sqrt{k_{x}^{2}+k_{y}^{2}}$. This leads to the following implementation

Note we have made a sign correction to (Mittet,, 2010, equation 46). In order to use high-order FDTD scheme, we also need to extrapolate electrical fields by forcing (Oristaglio and Hohmann,, 1984),

at the air interface.

The above boundary conditions result in the extrapolation of the fields to ghost points using Fourier transform. Because the use of fast Fourier transform (FFT) assumes equal grid spacing, the air-water boundary condition requires a uniform grid in both x- and y- directions. The use of non-uniform grid in x- and y- directions is feasible by interpolating between non-uniform grid and uniform grid. Our numerical result (Yang and Mittet,, 2023) shows that it leads to degraded modelling accuracy and thus not implemented in libEMM.

In order to mimic the wave propagation to infinity, we apply perfectly matched layer (PML) (Chew and Weedon,, 1994) boundary condition in other directions except the top part of the model. It helps to attenuate the potential reflection using truncated computing mesh. The application of PML is equivalent to perform coordinate stretching $\partial_{\tilde{x}}=s_{x}^{-1}\partial_{x}$ with a complex factor $s_{x}$ ($|s_{x}|>1$). The convolutional PML (CPML) (Roden and Gedney,, 2000) uses $s_{x}=\kappa_{x}+\frac{d_{x}(x)}{\alpha_{x}+i\omega}(\kappa_{x}\geq 1)$ such that

where the memory variable $\psi_{x}$ can then be computed in a recursive manner during time stepping

where

The damping profile $d_{x}(x)$ and the constant $\kappa_{x}$ are chosen according to (Komatitsch and Martin,, 2007). We implement the CPML for $x$, $y$ and $z$ directions in the same way in FDTD. Assume the resistivity/conductivity model is cube of size n1*n2*n3. The model will be padded with nb CPML layers and ne number of extra buffers on each side:

    emf->nbe = emf->nb + emf->ne;
    emf->n1pad = emf->n1 + 2*emf->nbe;//total number of gridpoints in x
    emf->n2pad = emf->n2 + 2*emf->nbe;//total number of gridpoints in y
    emf->n3pad = emf->n3 + 2*emf->nbe;//total number of gridpoints in z

An example for the memory variables associated with $\partial_{z}H_{x}$ and $\partial_{z}H_{y}$ are given as follows:

    ...
    if(i3<emf->nb){
      emf->memD3H2[i3][i2][i1] = emf->b3[i3]*emf->memD3H2[i3][i2][i1]
                               + emf->a3[i3]*D3H2;
      emf->memD3H1[i3][i2][i1] = emf->b3[i3]*emf->memD3H1[i3][i2][i1]
                               + emf->a3[i3]*D3H1;
      D3H2 += emf->memD3H2[i3][i2][i1];
      D3H1 += emf->memD3H1[i3][i2][i1];
    }else if(i3>emf->n3pad-1-emf->nb){
      j3 = emf->n3pad-1-i3;
      k3 = j3+emf->nb;
      emf->memD3H2[k3][i2][i1] = emf->b3[j3]*emf->memD3H2[k3][i2][i1]
                               + emf->a3[j3]*D3H2;
      emf->memD3H1[k3][i2][i1] = emf->b3[j3]*emf->memD3H1[k3][i2][i1]
                               + emf->a3[j3]*D3H1;
      D3H2 += emf->memD3H2[k3][i2][i1];
      D3H1 += emf->memD3H1[k3][i2][i1];
    }
    ...

After this step, the curl operator of the fields can then be obtained. For the curl of magnetic fields, the formula

will translate into

    ...
    emf->curlH1[i3][i2][i1] = D2H3-D3H2;
    emf->curlH2[i3][i2][i1] = D3H1-D1H3;
    emf->curlH3[i3][i2][i1] = D1H2-D2H1;
    ...

The above code implements the CPML regions and the interior regions in an elegant and compatible frame.

Equation (6) shows that many time steps may be required to compute the frequency domain electromagnetic fields. The frequency domain data should be integrated on the fly during forward modelling process thanks to the discrete time Fourier transform (DTFT).

where the time has been discretized as $t_{n}=n\Delta t$; $u\in\{E_{x},E_{y},E_{z},H_{x},H_{y},H_{z}\}$ is one component of the electromagnetic fields; $N_{t}$ is the total number of time steps needed. The following code snippet examplifies the DTFT for $E_{x}$ component at it time step:

void dtft_emf(emf_t *emf, int it, float ***E1, float _Complex ****fwd_E1)
/*<DTFT of the electromagnetic fields (emf) >*/
{
  int i1,i2,i3,ifreq;
  float _Complex factor;

  int i1min=emf->nb;//lower bound for index i1
  int i2min=emf->nb;//lower bound for index i2
  int i3min=emf->nb;//lower bound for index i3
  int i1max=emf->n1pad-1-emf->nb;//upper bound for index i1
  int i2max=emf->n2pad-1-emf->nb;//upper bound for index i2
  int i3max=emf->n3pad-1-emf->nb;//upper bound for index i3

  for(ifreq=0; ifreq<emf->nfreq; ++ifreq){
    factor = emf->expfactor[it][ifreq];

    for(i3=i3min; i3<i3max; ++i3){
      for(i2=i2min; i2<i2max; ++i2){
        for(i1=i1min; i1<i1max; ++i1){
          fwd_E1[ifreq][i3][i2][i1] += E1[i3][i2][i1]*factor;
        }
      }
    }
  }
}

where E1 and fwd_E1 stand for time-domain and frequency-domain $E_{x}$ field. Note that we only perform the computation in the region of interest without absorbing boundaries by specifying the index bounds. We do the same for all other field components.

There have been a significant amount of effort in the literature to estimate the total modelling time $T_{\max}$ (Wang and Hohmann,, 1993; Mittet,, 2010). The key for this estimation is to ensure that the frequency domain EM field reaches its steady state. In practice, the final number of time steps is usually over estimated for safety. In our implementation, we therefore regularly check the convergence of the frequency domain field until no contribution added after more number of time steps. This can be computationally expensive if the convergence check is very frequent. This is particularly true if one wishes to output the Green’s function at every grid point for every frequency of interest. Here, we take a short cut based on judicious physical intuition. Physically, the lowest frequency component of the CSEM has the largest penetration depth, according to the relation between the frequency and the skin depth $\delta=\sqrt{2/(\omega\mu\sigma)}$. Thus, using lowest frequency at the boundaries of the computing volume should be sufficient to monitor the evolution of the fields. The simulation will be automatically terminated once the field is converged, thus avoiding extra computation.

When an interface is present in the medium, additional effort is required to handle the high medium contrast in order to achieve precise modelling. This can be achieved by averaging the medium (Davydycheva et al.,, 2003). libEMM adapts it for VTI medium, even though this method is capable to handle full anisotropy. The key idea of this method is to compute the effective medium parameters for the horizontal conductivity and vertical resistivity by averaging over the cell size

The above procedure is based on the fact that the tangential field components see the model as a combination of resistors in parallel, while the normal field component sees the model as a combination of resistors in series. In the convention of Soviet literature, the method is often coined homogenization, as the same formula holds for all grid points, regardless of whether the point is in a homogeneous region or in the neighbourhood of an interface.

Our finite difference modelling is carried out on a rectilinear mesh, which is simply the tensor product of 1D non-equispaced meshes. To generate these 1D meshes, the rule of geometrical progression is used (Mulder,, 2006, Appendix C). Assume we have the total grid length $L$ divided into $n$ intervals ($n+1$ nodes) with a common ratio $q>1$. Denote the smallest interval $\Delta x=x_{1}-x_{0}$. Thus, the relation between $L$ and $\Delta x$ is

Due to the stability requirement and the resulting computational cost in the modelling, we are restricted to the smallest interval $\Delta x$ and a given number of nodes to discretize over certain distance. Since equation (27) does not yield an explicit expression for the stretching factor $q$, the question boils down to finding the optimal $q$ from fixed $n$, $\Delta x$ and $L$. The relation in (27) is equivalent to

which inspires us to carry out a number of fixed point iterations:

Since $|g^{\prime}(q)|<1$ holds for all $q>1$, the scheme proposed here is guaranteed to converge. This idea has been implemented in the following code snippet.

float create_nugrid(int n, float len, float dx, float *x)
/*< generate power-law nonuniform grid by geometric progression >*/
{
  int i;
  float q, qq;

  if(fabs(n*dx-len)<1e-15) {
    for(i=0; i<=n; i++) x[i] = i*dx;
    return 1;
  }

  q = 1.1;//initial guess for the solution
  while(1){
    qq = pow(len*(q-1.)/dx + 1., 1./n);
    if(fabs(qq-q)<1e-15) break;
    q = qq;
  }
  for(x[0]=0,i=1; i<=n; i++) x[i] = (pow(q,i) - 1.)*dx/(q-1.);

  return q;
}

Since programming 3D CSEM modelling involves sophisticated memory allocation and initialization using pointers, libEMM implements every computational module following the same naming convention to ensure a memory-safe code. Throughout the software, all routines with names xxxx_init() and xxxx_clos() serve as the constructor and the destructor. This allows the C routine xxxx.c resembling C++ class and Fortran modules. The following details the major steps of FDTD based CSEM modelling in fictitious wave domain, as sketched in Algorithm 1:

• •

  The pointer emf pointing to the data structure of electromagnetic fields will be initialized by emf_init(emf). The relevant parameters requiring computing memory (for example, emf->rho11, emf->rho22 and emf->rho33) will also be allocated for reading input model in emf_init(emf). Another routine named emf_close(emf) for destroying the variables allocated before will be placed in the same source file emf_init_close.c.

• •

  Following the same convention, the pointer acqui pointing to the data structure of survey acquisition will be initialized by acqui_init(). This initialization allocates variables and reads the input files by different MPI processor. The routine name acqui_close() will deallocate variables after the modelling. These are enclosed in the source file acqui_init_close.c.

• •

  The interpolation operator will be initialized by interpolation_init() prior to modelling and destroyed by interpolation_close() following the completion of the modelling jobs.

• •

  The computing model must be extended with CPML layers and some extra buffer layers. This is achieved in extend_model.c, with initialization by extend_model_init() and variable deallocation by extend_model_close().

• •

  The frequency domain EM fields are initialized by dtft_emf_init() and deallocated by dtft_emf_close(). The DTFT of EM fields will be performed at every time step by dtft_emf() in the same source file dtft_emf.c.

• •

  Air-wave boundary condition are implemented in airwave_bc.c, initialized by airwave_bc_init() and destroyed in airwave_bc_close(). At each time step, the airwave treatment for electrical and magnetic fields will be carried out via airwave_bc_update_E() and airwave_bc_update_H(), respectively.

• •

  The time-domain electromagnetic fields for FDTD modelling is initialized in fdtd_init() and destroyed in fdtd_close(). The leap-frog time stepping at each time step will call fdtd_curlH to compute $\nabla\times\mathbf{H}$, fdtd_update_E() to update electrical field $\mathbf{E}$, fdtd_curlE to compute $\nabla\times\mathbf{E}$, fdtd_update_H() to update magnetic field $\mathbf{H}$. These routines are bundled in fdtd.c.

• •

  The injection of the electromagnetic source for forward simulation, done by inject_electric_src_fwd() and inject_magnetic_src_fwd(), resides in inject_src_fwd.c.

• •

  The convergence is checked every 100 time steps by check_convergence().

• •

  After time domain modelling, the Green’s function is computed via compute_green_function().

• •

  The CSEM data is finally extracted via extract_emf() and written out by write_data().

To allocate memory for the variables used in CSEM modelling, we dynamically allocate arrays of pointers of arrays of pointers … using contiguous memory chunk, to construct multidimensional arrays such that these arrays can be referenced in the same way as static arrays in C. That is, the $(i,j,k)$-th element can be referenced simply by array[k][j][i]. The rightmost index i is the fastest index while the leftmost index k is the slowest one. The array of size n1*n2*n3 starts from the first element array[0][0][0], and ends with the last one array[n3-1][n2-1][n1-1].

The GPU acceleration technology has been mature after more than one decade of the development effort. Besides the successful applications in seismic community such as reverse time migration and full waveform inversion (Yang et al.,, 2014, 2015), the use of GPU for 3D CSEM is scarcely reported. Here we highlight some key strategies when implementing GPU-based 3D CSEM modeling.

1. 1.

   Perform computation intensive part on GPU without frequent CPU-GPU data traffic. Note that CPU emphasis on low latency, while GPU emphasis on high throughput. In order to avoid substantial amount of communication latency, we design the code to perform all necessary pre-computation as much as possible before porting to GPU-based time stepping. The pre-computation includes validation of the input parameter for sanity check, preparation of the medium on extended domain, tabulating the interpolating weights at the source and receiver locations in a lookup table. The air-water boundary condition is directly computed on device with CUFFT library.

2. 2.

   Use shared memory with tiled algorithm (Micikevicius,, 2009). The key idea is to explore the fast L1 cache of the shared memory, which is often quite small memory size but 2 order of magnitude faster than accessing the global memory. Tiling forces multiple threads to jointly focus on a subset of the input data at each phase of execution. Here we map the global memory onto 2D shared memory blocks and slide it along the third dimension, as illustrated in Figure 1. The reading and writing of the data on GPU will be performed block by block manner rather than element by element, thus much more efficient. Since we have to repeatedly reuse the same quantity at different location in the finite difference stencil, the tiling technique with shared memory enhances locality of data access.

3. 3.

   A linear increasing index with consecutive memory access. Each CUDA warp contains 32 processing cores executed in a single instruction multiple thread fashion. The CUDA kernel will serialize the operations when the elements reside in different warps, which leads to dramatic performance penalty. The linear increasing index with consecutive memory access is therefore crucial to design efficient the CUDA kernels.

4. 4.

   Improved scheme for convergence check of the fields. The modelling can be terminated only if the frequency domain EM fields converges at every grid point. This naturally leads to a reduction operation at all steps of convergence check. To do this efficiently on GPU, the author implemented a parallel reduction scheme in (Yang,, 2021). In this paper, we do it in a more judicious manner by only checking 8 corners of interior computing cube without absorbing boundaries, as sketched in Figure 2. This eliminates the need for parallel reduction and dramatically reduces the required number of floating point operations. The scheme can equally be applied in the sequential implementation.

5. 5.

   Use of mapped pinned memory. Due to convergence check, timestepping modelling involves the communication between the device (GPU) and host (CPU), which violates the first principle aforementioned. To avoid this, we invoke mapped pinned memory which is desirable for the CPU and GPU to share a buffer without explicit buffer allocations and copies. It is well-known that extensive use of mapped pinned memory may hit the performance. Here, only 8 corners of the complex-valued EM field has to be backed up using zero-copy pinned memory, which minimizes the potential performance deterioration.

The above techniques imply that parallel GPU implementation requires a significant adaptation of each subroutine in Algorithm 1 into different CUDA kernel functions, while the computational workflow remains the same but runs in on GPU device instead of CPU host. In GPU mode, the device must be initialized first, and then performing dedicated operations by these CUDA kernel functions, using different number of multithreaded block and grid sizes. This has been done by cuda_modelling.cu, as a replacement of the sequential modelling conducted by do_modeling.c. Some simple calculations may be sophisticated to achieve but can easily performed on CPU host only once. These operations, including the extension of the model for absorbing boundaries, the setup of CPML layers and DTFD coefficients, the computation of interpolation weights to inject the source and extract data at receiver loations, have been carried out priori to cuda_modeling.cu: GPU will directly copy them into device memory for direct utilization. This highlights the importance of exploring the distinct advantages of CPU and GPU computation.

All relevant CUDA kernels are placed in cuda_fdtd.cuh. In particular, the computation of the spatial derivatives in FDTD are $\nabla\times\mathbf{E}$ and $\nabla\times\mathbf{H}$, which have been parallelized by GPU kernels cuda_fdtd_curlE_nugrid and cuda_fdtd_curlH_nugrid, as the alternatives to CPU routines fdtd_curlE and fdtd_curlH. According to tiled algorithm, the first two dimenions are parallelized using shared memory, while the 3rd dimension strides sequentially, yielding the following code

__global__ void cuda_fdtd_curlH_nugrid(...)ΨΨ
{
  const int i1 = threadIdx.x + blockDim.x*blockIdx.x;
  const int i2 = threadIdx.y + blockDim.y*blockIdx.y;
  const int t1 = threadIdx.x + 2;//t1>=2 && t1<BlockSize1+2
  const int t2 = threadIdx.y + 2;//t2>=2 && t2<BlockSize2+2

  __shared__ float sH1[BlockSize2+3][BlockSize1+3];
  __shared__ float sH2[BlockSize2+3][BlockSize1+3];
  __shared__ float sH3[BlockSize2+3][BlockSize1+3];

  int in_idx = i1+n1pad*i2;
  int out_idx = 0;
  int stride = n1pad*n2pad;

  for(i3=i3min; i3<=i3max; i3++){//i3-2>=0 && i3+1<n3pad
    ...
    __syncthreads();

    if(validrw){
      ...
      curlH1[out_idx] = D2H3 - D3H2;
      curlH2[out_idx] = D3H1 - D1H3;
      curlH3[out_idx] = D1H2 - D2H1;
    }
  }//end for i3
}

Note that it is important to synchronize GPU threads before proceeding to next stage. For efficiency, the timestepping of $\mathbf{E}$ and $\mathbf{H}$ fields also conforms the tiled mapping, i.e.,

__global__ void cuda_fdtd_update_E(...)
{
  const int i1 = threadIdx.x + blockDim.x*blockIdx.x;
  const int i2 = threadIdx.y + blockDim.y*blockIdx.y;
  int i3, id;

  for(i3=0; i3<n3pad; i3++){
    if(i1<n1pad && i2<n2pad){
      id = i1+n1pad*(i2+ n2pad*i3);
      E1[id] += dt*inveps11[id]*curlH1[id];
      E2[id] += dt*inveps22[id]*curlH2[id];
      E3[id] += dt*inveps33[id]*curlH3[id];
    }
  }
}

The convergence of the field is checked directly on GPU by cuda_check_convergence. Since the result of convergence check must be visible on host, we use cudaHostAlloc instead of commonly used routine cudaMalloc, to allocate memory to record the number of converged corners in 3D cube. Correspondingly, the pinned memory must be freed using cudaFreeHost rather than cudaFree. To be concrete, t, these correpond to the following lines of code:

  cudaHostAlloc(&h_ncorner, sizeof(int), cudaHostAllocMapped);

and

  cudaFreeHost(h_ncorner);

libEMM is a light-weight yet powerful 3D CSEM modelling toolkit which can run using single CPU processor, multiple MPI parallel process and also CUDA-enabled GPU mode. This implies that the software has the following dependencies:

• •

  FFTW for fast Fourier transform;

• •

  MPICH (or other types of MPI) for multi-processor parallelization over distributed memory architecture;

• •

  Nvidia CUDA package if one wish to benefit from GPU parallelization on device shared memory.

The installation of CUDA package is optional while the first two are mandatory. libEMM uses the C programming language to achieve high performance computing. To compile the code, one will use the command make to compile the code following the given Makefile, which may be edited according to the path of the user’s installation for the above packages. The compiled executable will be named as fdtd and placed in the folder bin.

libEMM is capable to automatically parse the argument list based on the parameter specified in the token. The value for each keyword will be picked up after a equal sign =. Multiple values for the same keyword can be specified via comma-separated parameters. Different keywords will be distinguished via the space.

To specify the survey layout, three acquisition files in ASCII format must be given to prescribe the source and receiver locations as well as the their connection table. These files will be assigned to the arguments as

fsrc=sources.txt frec=receivers.txt fsrcrec=src_rec_table.txt

where the file sources.txt gives the source locations; the file receivers.txt gives the receiver locations; the file src_rec_table.txt establishes the connection between sources and receivers. When the code runs with MPI, each process will read these files and then carry out independent modelling jobs according to their process index (which uniquely determines a source). The source/receiver file is in the following form:

      x     y      z       azimuth    dip        iTx/iRx
-8000.0   0.0    275.0        0        0           1
-7975.0   0.0    275.0        0        0           2
-7950.0   0.0    275.0        0        0           3
-7925.0   0.0    275.0        0        0           4
    ...

The above 6 columns corresponds to $(x,y,z)$ coordinates (in meters), the azimuth and dip angles in rad and the index of the source/receivers. Thanks to the interpolation operator computed by inverting Vandermonde matrix, the sources and receivers can be accurately injected/extracted at arbitrary locations without the need to sit on grid. The connection table reads

 iTx iRx
 1    1
 1    2
 1    3
 ...

 2    1
 2    2
 2    3
 ...

which connects each source/transmitter indexed by iTx are associated with a number of receivers index by iRx.

libEMM carries out CSEM modelling over a 3D cube with physical dimensions specified by the bounds. Modelling over a domain $X=X_{1}\times X_{2}\times X_{3}$ with $X_{1}=X_{2}=[-10000,10000]$m, $X_{3}=[0,3000]$ reads

x1min=-10000 x1max=10000 x2min=-10000 x2max=10000 x3min=0 x3max=3000

The parameters n1, n2 and n3 specify the dimension of the model in x, y and z directions, while d1, d2 and d3 specify the smallest grid spacing (in meters) along each coordinate. To do modelling in a resistivity cube of $101^{3}$ grid points, $\Delta x=\Delta y=200$ m, $\Delta z=50$ m on uniform grid, one has to write in the argument list

n1=101 n2=101 n3=101 d1=200 d2=200 d3=50

These parameters will be checked by libEMM in order to match the bounds of the domain.

Because libEMM allows the CSEM modelling using both uniform grid and non-uniform grid, the grid coordinate in binary format, for argument fx3nu corresponding to the gridding in z direction, must be provided in order to model on non-uniform grid along z. Note that the file must store exactly n3 floating points. Non-uniform gridding in x- and y- directions has been forbidden due to degraded modelling accuracy (Yang and Mittet,, 2023). When the modelling is performed on uniform grid, there is no need to provide input for fx3nu.

Three resistivity files in binary format are expected to feed the argument list for numerical modelling, as specified in the following

frho11=rho11 frho22=rho22 frho33=rho33

where the files rho11, rho22 and rho33 are the resistivities of size n1*n2*n3 after homogenization, implying a half grid shift in x, y and z directions, respectively. These files should be built using model building tools outside the modelling kernel.

libEMM specifies the active source channel and receiver channels via the keywords - chsrc and chrec. For example, recording CSEM data of components $E_{x}$ and $E_{y}$ from an electrical source $J_{x}$ reads

chsrc=Ex chrec=Ex,Ey

Similarly, one can specify the simulation frequencies as

freqs=0.25,0.75,1.25

in which 0.25 Hz, 0.75 Hz and 1.25 Hz will be the frequencies extracted from time-domain CSEM modelling engine.

The half length of the finite-difference operator is determined by parameter rd. It is possible to choose rd=1, 2 or 3 in CPU mode, but only rd=2 is supported on GPU. This is because dramatic accuracy improvement has been observed in modelling when switching from the 2nd to the 4th order scheme, but increasing from 4th order to even higher order requires more computational effort without significant accuracy gain (Yang and Mittet,, 2023).

We pad the model with the same number of CPML layers on each side of the cube, specified by the parameter nb. In addition, we add ne number of buffers as a smooth transition zone between CPML and interior domain. A general principle for choosing ne is to ensure the number of buffer layers larger than the half length of finite-difference stencil, for example, $\mathtt{ne=rd+5}$.

The outcome after running a 3D CSEM modelling will be output as ASCII files named as emf_XXXX.txt, where XXXX is the index of the sources/transmitters. Running the simulation using 4 MPI process leads to

emf_0001.txt emf_0002.txt emf_0003.txt emf_0004.txt

Each file will follow the same convention to print out the transmitter index iTx, the receiver index iRx, the receiver channel Ex/Ey/Ez/Hx/Hy/Hz, the frequency index ifreq and the real and imaginary part of the complex-valed EM field in frequency domain:

iTx  iRx  chrec  ifreq Ψ  emf_real Ψ   emf_imag
1    1 Ψ   Ex Ψ   1 Ψ -1.086630e-14 Ψ 2.026699e-16
1    2 Ψ   Ex Ψ   1 Ψ -1.148639e-14 Ψ 3.976523e-16
1    3 Ψ   Ex Ψ   1 Ψ -1.213745e-14 Ψ 6.254443e-16
1    4 Ψ   Ex Ψ   1 Ψ -1.282039e-14 Ψ 8.894743e-16
1    5 Ψ   Ex Ψ   1 Ψ -1.353619e-14 Ψ 1.192894e-15
1    6 Ψ   Ex Ψ   1 Ψ -1.428555e-14 Ψ 1.539961e-15
1    7 Ψ   Ex Ψ   1 Ψ -1.506933e-14 Ψ 1.934350e-15
1    8 Ψ   Ex Ψ   1 Ψ -1.588794e-14 Ψ 2.380965e-15
1    9 Ψ   Ex Ψ   1 Ψ -1.674208e-14 Ψ 2.884079e-15
1    10    Ex Ψ   1 Ψ -1.763176e-14 Ψ 3.449250e-15
1    11    Ex Ψ   1 Ψ -1.855739e-14 Ψ 4.081327e-15
...

To run the 3D CSEM modelling, one may write the parameters in a shell script run.sh and launch the code using bash run.sh. The following is an example shell script running with 25 MPI process, each of the MPI process parallelized by 4 OpenMP threads:

#!/usr/bin/bash

n1=101
n2=101
n3=101
d1=200
d2=200
d3=50

export OMP_NUM_THREADS=4
mpirun -n 25 ../bin/fdtd \
       mode=0 \
       fsrc=sources.txt \
       frec=receivers.txt \
       fsrcrec=src_rec_table.txt \
       frho11=rho11 \
       frho22=rho22 \
       frho33=rho33 \
       chsrc=Ex \
       chrec=Ex \
       x1min=-10000 x1max=10000 \
       x2min=-10000 x2max=10000 \
       x3min=0 x3max=5000 \
       n1=$n1 n2=$n2 n3=$n3 \
       d1=$d1 d2=$d2 d3=$d3 \
       nb=12 ne=6 \
       freqs=0.25,0.75,1.25 \
       rd=2

where there are 25 sources in the acquisition. It is also possible to do modelling for some specific sources, for example,

mpirun -n 2 ../bin/fdtd  shots=12,14 ...

will only simulate CSEM data for transmitter-12 and transmitter-14. The output CSEM response will be named as emf_0012.txt and emf_0014.txt. The code will terminate if the number of process launched for mpirun exceeds the total number of sources given in the file sources.txt.

The code can be compiled with nvcc compiler for GPU modelling using command ‘make GPU=1’. Assume there is only one GPU card available on the laptop. Once the executable is created after compilation, one can directly run the code by executing the command with proper parameters:

../bin/fdtd ...

Adding nvprof in the script provides us a convenient way to profile the computing time spent on different CUDA kernels

nvprof --log-file profiling.txt ../bin/fdtd ...

where the profiling log will be recorded in profiling.txt.

The first test is a layered resistivity model shown in Figure 3: the model has 5 layers, the air of $10^{12}\;\Omega\cdot\mbox{m}$, the water of $0.3125\;\Omega\cdot\mbox{m}$, and two layers of formations. Between two layers of formations, there exists a resistive layer of $100\;\Omega\cdot\mbox{m}$ mimicking hydrocarbon bearing formations. The 3D FDTD code is now cross-validated against the semi-analytic solution calculated by Dipole1D program (Key,, 2009). We discretize the model of dimension $10\;\mbox{km}\times 10\;\mbox{km}\times 5\;\mbox{km}$ including a resistor at the depth 1250-1350 m, with grid spacing $\Delta x=\Delta y=100\;\mbox{m}$, $\Delta z=50\;\mbox{m}$. This results in a large 3D grid of size $101^{3}$. A point source is placed at 50 m above the seabed. Figures 4a and 4b show very good agreement between the FDTD method and the semi-analytic solution, in terms of the amplitude and the phase of the $E_{x}$ component. Figures 4c and 4d show less than $1.5\%$ of the amplitude error and less than $1^{\circ}$ of the phase error at most of the offsets. The error at the near offset are very large. Fortunately the near offset data are not relevant for practical CSEM imaging. This example confirms the good accuracy of the FDTD method.

To mimic practical marine CSEM environment, we perform CSEM model using a 3D resistivity model of 2D structure including seafloor bathymetry, as shown in Figure 5. To catch the impact of seafloor bathymetry, we mesh the model densely around the seabed, while gradually stretching the grid along the depth below. The x-z section of the mesh is shown in Figure 6. With such a model, we compute a reference solution using MARE2DEM (Key,, 2016) using finite element method. It extends the model tens of kilometers in each direction, to mimic that EM fields propagate to very far distance while avoiding possible edge/boundary effect. The regions of interest has been densely gridded using Delaunay triangulation, while big cells are employed at far distance away from interior part of the model. In our finite difference modelling, the PML boundary condition attenuates the artificial reflections in the computational domain within tens of grids to achieve the same behavior. Both the amplitude and phase match very well between FDTD modelling and MARE2DEM solution, as can be seen in Figure 7. Figure 7c shows that the amplitude error are bounded within 3% for all frequencies at relevant offsets; the amplitude error at far offset and higher frequencies becomes larger when approaching to noise level ($10^{-15}\;\mbox{V/m}^{2}$).

When libEMM is launched using multiple MPI process, we can simultaneously model CSEM data associated with different transmitters while all receivers are active/alive. This mimics the practical situations in real marine CSEM acquisitions. For a survey layout in Figure 8, libEMM helps us to obtain both inline and broadside data, as shown in Figure 9.