Generating C Code for the Scalar Wave Equation in Cartesian Coordinates¶

Authors: Zach Etienne & Thiago Assumpção¶

Formatting improvements courtesy Brandon Clark¶

This module generates the C Code for the Scalarwave in Cartesian coordinates and sets up either a monochromatic plane wave or spherical Gaussian Initial Data. The module emphasizes the use of the Method of Lines, thus transforming the partial differential equation problem into a more manageable ordinary differential equation problem.¶

Notebook Status: Validated

Validation Notes: This tutorial notebook has been confirmed to be self-consistent with its corresponding NRPy+ module, as documented below (right-hand-side expressions; initial data expressions). In addition, all expressions have been validated against a trusted code (the original SENR/NRPy+ code).

NRPy+ Source Code for this module:¶

• ScalarWave/ScalarWave_RHSs.py
• ScalarWave/InitialData.py

Introduction:¶

Problem Statement¶

We wish to numerically solve the scalar wave equation as an initial value problem in Cartesian coordinates: $$\partial_t^2 u = c^2 \nabla^2 u \text{,}$$ where $u$ (the amplitude of the wave) is a function of time and space: $u = u(t,x,y,...)$ (spatial dimension as-yet unspecified) and $c$ is the wave speed, subject to some initial condition

$$u(0,x,y,...) = f(x,y,...)$$

and suitable spatial boundary conditions.

As described in the next section, we will find it quite useful to define $$v(t,x,y,...) = \partial_t u(t,x,y,...).$$

In this way, the second-order PDE is reduced to a set of two coupled first-order PDEs

\begin{align} \partial_t u &= v \\ \partial_t v &= c^2 \nabla^2 u. \end{align}

We will use NRPy+ to generate efficient C codes capable of generating both initial data $u(0,x,y,...) = f(x,y,...)$; $v(0,x,y,...)=g(x,y,...)$, as well as finite-difference expressions for the right-hand sides of the above expressions. These expressions are needed within the Method of Lines to "integrate" the solution forward in time.

The Method of Lines¶

Once we have initial data, we "evolve it forward in time", using the Method of Lines. In short, the Method of Lines enables us to handle

1. the spatial derivatives of an initial value problem PDE using standard finite difference approaches, and
2. the temporal derivatives of an initial value problem PDE using standard strategies for solving ordinary differential equations (ODEs), so long as the initial value problem PDE can be written in the form

$$\partial_t \vec{f} = \mathbf{M}\ \vec{f},$$

where $\mathbf{M}$ is an $N\times N$ matrix filled with differential operators that act on the $N$-element column vector $\vec{f}$. $\mathbf{M}$ may not contain $t$ or time derivatives explicitly; only spatial partial derivatives are allowed to appear inside $\mathbf{M}$. The scalar wave equation as written in the previous module, \begin{equation} \partial_t \begin{bmatrix} u \\ v \end{bmatrix}= \begin{bmatrix} 0 & 1 \\ c^2 \nabla^2 & 0 \end{bmatrix} \begin{bmatrix} u \\ v \end{bmatrix}, \end{equation} satisfies this requirement.

Thus we can treat the spatial derivatives $\nabla^2 u$ of the scalar wave equation using standard finite-difference approaches, and the temporal derivatives $\partial_t u$ and $\partial_t v$ using standard approaches for solving ODEs. In the next module, we will apply the highly robust explicit Runge-Kutta fourth-order scheme (RK4), used widely for numerically solving ODEs, to "march" (integrate) the solution vector $\vec{f}$ forward in time from its initial value ("initial data").

Basic Algorithm¶

The basic algorithm for solving the scalar wave equation initial value problem, based on the Method of Lines (see section above), is outlined below, with NRPy+-based components highlighted in green. We will review how NRPy+ generates these core components in this module.

1. Allocate memory for gridfunctions, including temporary storage for the RK4 time integration.
2. Set gridfunction values to initial data.
3. Evolve the system forward in time using RK4 time integration. At each RK4 substep, do the following:
   1. Evaluate scalar wave RHS expressions.
   2. Apply boundary conditions.

We refer to the right-hand side of the equation $\partial_t \vec{f} = \mathbf{M}\ \vec{f}$ as the RHS. In this case, we refer to the $\mathbf{M}\ \vec{f}$ as the "scalar wave RHSs". In the following sections we will

1. use NRPy+ to cast the scalar wave RHS expressions -- in finite difference form -- into highly efficient C code,
   1. first in one spatial dimension with fourth-order finite differences,
   2. and then in three spatial dimensions with tenth-order finite differences; we will
2. use NRPy+ to generate monochromatic plane-wave initial data for the scalar wave equation, where the wave propagates in an arbitrary direction.

As for the $\nabla^2 u$ term, spatial derivatives are handled in NRPy+ via finite differencing.

We will sample the solution $\{u,v\}$ at discrete, uniformly-sampled points in space and time. For simplicity, let's assume that we consider the wave equation in one spatial dimension. Then the solution at any sampled point in space and time is given by $$u^n_i = u(t_n,x_i) = u(t_0 + n \Delta t, x_0 + i \Delta x),$$ where $\Delta t$ and $\Delta x$ represent the temporal and spatial resolution, respectively. $v^n_i$ is sampled at the same points in space and time.

Table of Contents¶

$$\label{toc}$$

1. Step 1: Initialize core NRPy+ modules
2. Step 2: Scalar Wave RHSs in One Spatial Dimension, Fourth-Order Finite Differencing
   1. Step 2.a: C-code output example: Scalar wave RHSs with 4th order finite difference stencils
3. Step 3: Scalar Wave RHSs in Three Spatial Dimensions
   1. Step 3.a: Code Validation against ScalarWave.ScalarWave_RHSs NRPy+ module
   2. Step 3.b: C-code output example: Scalar wave RHSs with 10th order finite difference stencils and SIMD enabled
4. Step 4: Setting up Initial Data for the Scalar Wave Equation
   1. Step 4.a: The Monochromatic Plane-Wave Solution
   2. Step 4.b: The Spherical Gaussian Solution (Courtesy Thiago Assumpção)
5. Step 5: Code Validation against ScalarWave.InitialData NRPy+ module
6. Step 6: Output this notebook to $\LaTeX$-formatted PDF file

Step 1: Initialize core NRPy+ modules [Back to top]¶

$$\label{initializenrpy}$$

Let's start by importing all the needed modules from NRPy+:

Step 2: Scalar Wave RHSs in One Spatial Dimension [Back to top]¶

$$\label{rhss1d}$$

To minimize complication, we will first restrict ourselves to solving the wave equation in one spatial dimension, so \begin{align} \partial_t u &= v \\ \partial_t v &= c^2 \nabla^2 u \\ &= c^2 \partial_x^2 u. \end{align}

We will construct SymPy expressions of the right-hand sides of $u$ and $v$ using NRPy+ finite-difference notation to represent the derivative, so that finite-difference C-code kernels can be easily constructed.

Extension of this operator to higher spatial dimensions when using NRPy+ is straightforward, as we will see below.

Step 2.a: C-code output example: Scalar wave RHSs with 4th order finite difference stencils [Back to top]¶

$$\label{ccode1d}$$

As was discussed in the finite difference section of the tutorial, NRPy+ approximates derivatives using finite-difference methods, the second-order derivative $\partial_x^2$ accurate to fourth-order in uniform grid spacing $\Delta x$ (from fitting the unique 4th-degree polynomial to 5 sample points of $u$) is given by \begin{equation} \left[\partial_x^2 u(t,x)\right]j = \frac{1}{(\Delta x)^2} \left( -\frac{1}{12} \left(u{j+2} + u_{j-2}\right)

• \frac{4}{3} \left(u_{j+1} + u_{j-1}\right)

• \frac{5}{2} u_j \right)

• \mathcal{O}\left((\Delta x)^4\right).

\end{equation}

Success! Notice that indeed NRPy+ was able to compute the spatial derivative operator, \begin{equation} \left[\partial_x^2 u(t,x)\right]j \approx \frac{1}{(\Delta x)^2} \left( -\frac{1}{12} \left(u{j+2} + u_{j-2}\right)

• \frac{4}{3} \left(u_{j+1} + u_{j-1}\right)

• \frac{5}{2} u_j \right),

\end{equation} correctly (it's easier to read in the "Original SymPy expressions" comment block at the top of the C output).

As NRPy+ is designed to generate codes in arbitrary coordinate systems, instead of sticking with Cartesian notation for 3D coordinates, $x,y,z$, we instead adopt $x_0,x_1,x_2$ for our coordinate labels. Thus you will notice the appearance of invdx0$=1/\Delta x_0$, where $\Delta x_0$ is the (uniform) grid spacing in the zeroth or $x_0$ direction. In this case, $x_0$ represents the $x$ direction.

Step 3: Scalar Wave RHSs in Three Spatial Dimensions [Back to top]¶

$$\label{rhss3d}$$

Let's next repeat the same process, only this time for the scalar wave equation in 3 spatial dimensions (3D).

Step 3.a: Validate SymPy expressions against ScalarWave.ScalarWave_RHSs NRPy+ module [Back to top]¶

$$\label{code_validation1}$$

Here, as a code validation check, we verify agreement in the SymPy expressions for the RHSs of the three-spatial-dimension Scalar Wave equation (i.e., uu_rhs and vv_rhs) between

1. this tutorial and
2. the NRPy+ ScalarWave.ScalarWave_RHSs module.

Step 3.b: C-code output example: Scalar wave RHSs with 10th order finite difference stencils and SIMD enabled [Back to top]¶

$$\label{ccode3d}$$

Next, we'll output the above expressions as Ccode, using the NRPy+ finite-differencing C code kernel generation infrastructure. This code will represent spatial derivatives as 10th-order finite differences and output the C code with SIMD enabled. (Common-subexpression elimination is enabled by default.)

Step 4: Setting up Initial Data for the Scalar Wave Equation [Back to top]¶

$$\label{id}$$

Step 4.a: The Monochromatic Plane-Wave Solution [Back to top]¶

$$\label{planewave}$$

The solution to the scalar wave equation for a monochromatic (single-wavelength) wave traveling in the $\hat{k}$ direction is $$u(\vec{x},t) = f(\hat{k}\cdot\vec{x} - c t),$$ where $\hat{k}$ is a unit vector. We choose $f(\hat{k}\cdot\vec{x} - c t)$ to take the form $$ f(\hat{k}\cdot\vec{x} - c t) = \sin\left(\hat{k}\cdot\vec{x} - c t\right) + 2, $$ where we add the $+2$ to ensure that the exact solution never crosses through zero. In places where the exact solution passes through zero, the relative error (i.e., the most common error measure used to check that the numerical solution converges to the exact solution) is undefined. Also, $f(\hat{k}\cdot\vec{x} - c t)$ plus a constant is still a solution to the wave equation.

Next we verify that $f(\hat{k}\cdot\vec{x} - c t)$ satisfies the wave equation, by computing $$\left(c^2 \nabla^2 - \partial_t^2 \right)\ f\left(\hat{k}\cdot\vec{x} - c t\right),$$ and confirming the result is exactly zero.

Step 4.b: The Spherical Gaussian Solution [Back to top]¶

$$\label{sphericalgaussian}$$

Here we will implement the spherical Gaussian solution, which consists of ingoing and outgoing wavefronts: \begin{align} u(r,t) &= u_{\rm out}(r,t) + u_{\rm in}(r,t) + 1,\ \ \text{where}\\ u_{\rm out}(r,t) &=\frac{r-ct}{r} \exp\left[\frac{-(r-ct)^2}{2 \sigma^2}\right] \\ u_{\rm in}(r,t) &=\frac{r+ct}{r} \exp\left[\frac{-(r+ct)^2}{2 \sigma^2}\right] \\ \end{align} where $c$ is the wavespeed, and $\sigma$ is the width of the Gaussian (i.e., the "standard deviation").

Since the wave equation is linear, both the leftgoing and rightgoing waves must satisfy the wave equation, which implies that their sum also satisfies the wave equation.

Next we verify that $u(r,t)$ satisfies the wave equation, by computing $$\left(c^2 \nabla^2 - \partial_t^2 \right)\left\{u_{\rm R}(r,t)\right\},$$

and

$$\left(c^2 \nabla^2 - \partial_t^2 \right)\left\{u_{\rm L}(r,t)\right\},$$

are separately zero. We do this because SymPy has difficulty simplifying the combined expression.

Step 5: Code Validation against ScalarWave.InitialData NRPy+ module [Back to top]¶

$$\label{code_validation2}$$

As a code validation check, we will verify agreement in the SymPy expressions for plane-wave initial data for the Scalar Wave equation between

1. this tutorial and
2. the NRPy+ ScalarWave.InitialData module.

Step 6: Output this notebook to $\LaTeX$-formatted PDF file [Back to top]¶

$$\label{latex_pdf_output}$$

The following code cell converts this Jupyter notebook into a proper, clickable $\LaTeX$-formatted PDF file. After the cell is successfully run, the generated PDF may be found in the root NRPy+ tutorial directory, with filename Tutorial-ScalarWave.pdf. (Note that clicking on this link may not work; you may need to open the PDF file through another means.)