NRPy+ Tutorial: Solving the Scalar Wave Equation with NumPy¶

Authors: Zach Etienne, Patrick Nelson, Terrence Pierre Jacques, Thiago Assumpção, Leo Werneck, Brandon Clark¶

(Note on Authors: Zach wrote the NRPy+ infrastructure, as well as this notebook and its Python code; Patrick wrote the first version of the NRPy+-based Einstein Toolkit thorns for solving the scalar wave equation in 2018; Terrence rewrote these thorns to the latest NRPy+ standards in 2020, along the lines of the BaikalETK thorns; Thiago extended the scalar wave initial data infrastructure and contributed fixes to the original NRPy+ scalar wave notebooks; Leo created the boundary condition animation below; and Brandon established NRPy+ notebook formatting standards.)

This notebook was first written as a tutorial to introduce NRPy+ during the 2020 Einstein Toolkit Workshop.

This notebook constructs and validates a code that solves the scalar wave equation $\partial_t^2 u = c^2 \nabla^2 u$ using NumPy, to both review the basic components of a numerical PDE solver and motivate the use of NRPy+.¶

This notebook was written to explicitly outline the basic algorithms for numerically solving hyperbolic PDEs, including Einstein's equations of general relativity in e.g., the BSSN formalism.

While the codes here are written by hand, the objective of the notebook is to motivate the use of NRPy+, which can be used to generate much of this code either automatically or from validated templates, greatly reducing the possibility of human error.

Notebook Status: Validated

Validation Notes: The code developed in this tutorial notebook has been confirmed to converge to the exact solution at the expected rate, as documented below.

Table of Contents¶

$$\label{toc}$$

1. Step 0: Prerequisites & Additional Reading
2. Step 1: Introduction: The scalar wave equation
   1. Step 1.a: Mathematical problem statement
   2. Step 1.b: Chosen solution to the scalar wave equation
   3. Step 1.c: Initial condition
3. Step 2: The Method of Lines (MoL)
4. Step 3: NumPy Implementation: Basic Algorithm
   1. Step 3.a: Set up the numerical grid and free parameters
   2. Step 3.b: Set up the initial data
   3. Step 3.c: Allocate memory for the gridfunctions storing $u$ and $v$, and define the indexing macro function
   4. Step 3.d: Define the right-hand sides of the PDEs
   5. Step 3.e: Boundary Conditions
   6. Step 3.f: The Method of Lines
   7. Step 3.g: The main driver function
5. Step 4: Argh, the code is SLOW! Why use NRPy+ instead?
6. Step 5: Error analysis & code validation: Confirming numerical errors converge to zero at the expected rate
7. Step 6: Exercises for students
8. Step 7: Output this notebook to $\LaTeX$-formatted PDF file

Step 0: Prerequisites & Additional Reading [Back to top]¶

$$\label{prereqs}$$

This tutorial assumes basic familiarity with computer programming, undergraduate mathematics, and computational physics or numerical analysis.

For additional reading, please consult the following links:

• Online Resources for Numerical Analysis & Basic Mathematics
• Numerical Recipes in C Feel free to use the free "obsolete" (but not really!) versions. Note that Numerical Recipes were written by numerical relativists!

Step 1: Introduction: The scalar wave equation [Back to top]¶

$$\label{intro}$$

We will write a Python code (based in NumPy) to numerically solve the scalar wave equation

$$\partial_t^2 u(t,x,y,z) = c^2 \nabla^2 u(t,x,y,z),$$

given initial conditions $u(t=t_0,x,y,z)$.

Step 1.a: Mathematical problem statement [Back to top]¶

$$\label{problem_statement}$$

We will numerically solve the scalar wave equation as an initial value problem in Cartesian coordinates:

$$\partial_t^2 u = c^2 \nabla^2 u,$$ 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 (we'll stick with simple extrapolation boundary conditions at first).

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 in time

$$\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.

Step 1.b: Chosen solution to the scalar wave equation [Back to top]¶

$$\label{solution}$$

Here we will implement the spherical Gaussian solution to the scalar wave equation, which consists of ingoing and outgoing wavefronts:

$$\begin{align} u(r,t) &= u_{\rm out}(r,t) + u_{\rm in}(r,t) + u_{\rm offset},\ \ \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, $u_{\rm offset}$ is a freely specifiable constant offset, and $\sigma$ is the width of the Gaussian (i.e., the "standard deviation"). Next, we'll demonstrate using SymPy (a computer algebra system written in Python, on which NRPy+ depends) that the above indeed is a solution to the scalar wave equation.

In Cartesian coordinates we have

$$\partial_t^2 u(t,x,y,z) = c^2 \nabla^2 u(t,x,y,z),$$ and we know that $$r = \sqrt{x^2+y^2+z^2},$$ which we implement below using SymPy, the Python computer algebra package that NRPy+ uses.

Here is a visualization of this solution over time $u_{\rm exact}(t)$ in the $x-z$ plane, generated using matplotlib:

Now let's confirm the solution solves the PDE:

$$\partial_t^2 u(t,x,y,z) = c^2 \nabla^2 u(t,x,y,z),$$ by confirming that $$\partial_t^2 u(t,x,y,z) - c^2 \nabla^2 u(t,x,y,z) = 0,$$ using SymPy to compute the second derivatives of $u$ symbolically. To make it easier for SymPy to simplify the resulting expression, we will split up the above equation into:

$$\partial_t^2 (u_{\rm in} + u_{\rm out} + u_{\rm offset}) - c^2 \nabla^2 (u_{\rm in} + u_{\rm out} + u_{\rm offset}) = 0,$$

and confirm that

$$\begin{align} \partial_t^2 u_{\rm in} - c^2 \nabla^2 u_{\rm in} &= 0 \\ \partial_t^2 u_{\rm out} - c^2 \nabla^2 u_{\rm out} &= 0 \\ \partial_t^2 u_{\rm offset} - c^2 \nabla^2 u_{\rm offset} &= 0, \end{align}$$ which must be the case since the scalar wave equation is a linear PDE, in which each of the waves (ingoing, outgoing, and constant) must satisfy the wave equation separately:

Step 1.c: Initial Condition [Back to top]¶

$$\label{initial_condition}$$

We will choose the above solution at $t=0$, $u(t=0,x,y,z)$, as our initial data and adopt the Method of Lines (described next) to advance the solution forward in time (i.e., solve the initial value problem).

Step 2: The Method of Lines (MoL) [Back to top]¶

$$\label{mol}$$

Once we have set our initial conditions (usually referred to as our "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.

Here 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").

Here's how MoL works.

The RK4 method is usually presented for solving the ODE

$$y'(t) = f(y,t)$$ as follows. Given initial data $y(t_0)=y_0$, one can construct the solution at any later time via the algorithm: $$\begin{align} k_1 &= f(y_n, t_n), \\ k_2 &= f(y_n + \frac{1}{2}\Delta tk_1, t_n + \frac{\Delta t}{2}), \\ k_3 &= f(y_n + \frac{1}{2}\Delta tk_2, t_n + \frac{\Delta t}{2}), \\ k_4 &= f(y_n + \Delta tk_3, t_n + \Delta t), \\ y_{n+1} &= y_n + \frac{1}{6}\Delta t(k_1 + 2k_2 + 2k_3 + k_4) + \mathcal{O}\big((\Delta t)^5\big). \end{align}$$

Our PDE involves two variables $u$ and $v$, and the algorithm generalizes in exactly the same manner as it would if we were solving a system of coupled ODEs with two variables. Further, our PDE does not contain explicit time dependence, which simplifies the algorithm a bit:

$$\begin{align} k_{1,u} &= f_u(u_n,v_n) = f_u(v_n) = v_n, \\ k_{1,v} &= f_v(u_n,v_n) = f_v(u_n) = c^2\nabla^2 u_n, \\ k_{2,u} &= f_u\left(v_n + \frac{1}{2}\Delta tk_{1,v}\right) = v_n + \frac{1}{2}\Delta tk_{1,v}\\ k_{2,v} &= f_v\left(u_n + \frac{1}{2}\Delta tk_{1,u}\right) = c^2\nabla^2 \left(u_n + \frac{1}{2}\Delta tk_{1,u}\right), \\ k_{3,u} &= f_u\left(v_n + \frac{1}{2}\Delta tk_{2,v}\right) = v_n + \frac{1}{2}\Delta tk_{2,v}\\ k_{3,v} &= f_v\left(u_n + \frac{1}{2}\Delta tk_{2,u}\right) = c^2\nabla^2 \left(u_n + \frac{1}{2}\Delta tk_{2,u}\right), \\ k_{4,u} &= f_u(v_n + \Delta tk_{3,v}) = v_n + \Delta tk_{3,v}\\ k_{4,v} &= f_v(u_n + \Delta tk_{3,u}) = c^2\nabla^2 \left(u_n + \Delta tk_{3,u}\right), \\ u_{n+1} &= u_n + \frac{1}{6}\Delta t(k_{1,u} + 2k_{2,u} + 2k_{3,u} + k_{4,u}) + \mathcal{O}\big((\Delta t)^5\big)\\ v_{n+1} &= v_n + \frac{1}{6}\Delta t(k_{1,v} + 2k_{2,v} + 2k_{3,v} + k_{4,v}) + \mathcal{O}\big((\Delta t)^5\big). \end{align}$$

Thus, given initial data $u_0$ and $v_0$, we can use the above algorithm to advance the solution forward in time by one timestep, to $u_1$ and $v_1$. Recall the $\nabla^2 u$ terms in the above expressions are computed using finite-difference derivatives. Since finite-difference derivatives require neighboring points to be evaluated, we only evaluate the $k_i$'s in the interior of the grid; at each step, we apply boundary conditions to fill in the outermost neighboring points (called ghost zones).

Step 3: NumPy Implementation: Basic Algorithm [Back to top]¶

$$\label{basicalg}$$

We will store the numerical solution $u$ and its time derivative $v$, at a given instant in time on a three-dimensional numerical grid. Since these variables are defined at each point on the numerical grid, we call them gridfunctions.

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.

Armed with these definitions, the basic algorithm for solving the scalar wave equation initial value problem, based on the Method of Lines (see section above) is outlined below.

1. Set up the numerical grid and free parameters
2. Allocate memory for gridfunctions, including temporary storage needed for the RK4 time integration.
3. Set gridfunction values to initial data.
4. 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.

In the following sections, we will implement this algorithm to solve the scalar wave equation in 3D by hand using NumPy, and to motivate the use of NRPy+.

Step 3.a: NumPy Implementation: Set up the numerical grid and free parameters [Back to top]¶

$$\label{numgrid_freeparams}$$

We will solve the scalar wave equation on a uniform Cartesian grid with Nx by Ny by Nz coordinate points in the $x$, $y$, and $z$ directions respectively. Since the grid is uniform, we can describe the $x$ coordinate of any gridpoint with a single integer $i$, and the same holds true for $y$ and $z$, for which we will use integers $j$ and $k$. Thus we will label each gridpoint $(x_i,y_j,z_k)$.

Let's choose a "cell-centered" grid, which will store the solution at points

$$x_i \in \{..., -\frac{3}{2} \Delta x, -\frac{1}{2} \Delta x, +\frac{1}{2} \Delta x, +\frac{3}{2} \Delta x ...\};$$ and we will allow for two additional ghost zones on the outer boundary to account for the fourth-order finite differencing we will implement to numerically compute $\nabla^2 u$. Thus the expression for computing $x_i$ will be

$$x_i = x_{\rm min} + \left( (i-\text{NGHOSTS}) + \frac{1}{2} \right) \Delta x,$$

where $\Delta x$ is the spacing between gridpoints, and $x_{\rm min}$ denotes the minimum grid value in $x$. We will solve this equation on a cube centered on the origin with the domain_size=10, where $x_{\rm min}=$ -domain_size, $y_{\rm min}=$ -domain_size, $z_{\rm min}=$ -domain_size, $x_{\rm max}=$ +domain_size, and so forth. We'll also choose Nx=Ny=Nz, so that

$$\Delta x = \Delta y = \Delta z = \frac{x_{\rm max}-x_{\rm min}}{\text{Nx}}.$$

Next we set the free parameters for the scalar wave solution:

Then we set the timestep, which is governed by the CFL condition, and the final time t_final, relative to the chosen start time $t_0$ (usually $t_0=0$), so that the points closest to origin aren't affected by the approximate boundary condition:

Step 3.b: NumPy Implementation: Set up the initial data [Back to top]¶

$$\label{numpy_id}$$

Now we'll set up exact_solution_all_points(time, u, v), which numerically evaluates the solution for $u$ and $v$ at all gridpoints at a given numerical time time.

Recall the exact solution is given by

$$\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}$$

Exercise for students: Prove that at $t=0$, $v=\partial_t u \equiv 0$.

The problem is, SymPy expressions need to be converted to NumPy expressions; otherwise using functions like sp.N() will be incredibly slow. So we attempt to fix this by some simple string manipulations, some for $v$ were done by hand using the below Python function.

To store the solution $u$ and $v$ at all gridpoints on our numerical grid cube requires

2*Nx*Ny*Nz*double

bytes of memory, where double is the amount of memory storage (in bytes) needed to store one double-precision number (this is 8, by the way).

Step 3.c: NumPy Implementation: Allocate memory for the gridfunctions storing $u$ and $v$, and define the indexing macro function [Back to top]¶

$$\label{numpy_gfs}$$

As is done in the Einstein Toolkit and native NRPy+ codes, instead of declaring multi-dimensional arrays (e.g., a 3D array), we will instead declare $u$ and $v$ as one-dimensional arrays u[ijk] and v[ijk], each with (Nx+2*NGHOSTS)*(Ny+2*NGHOSTS)*(Nz+2*NGHOSTS) gridpoints. To access data an arbitrary point $(x_i,y_j,z_k)$, we need only call a simple function to find the correct index ijk given the grid indices i, j, and k, which label the point $(x_i,y_j,z_k)$:

$$\verb| (i,j,k) = i + (Nx+2*NGHOSTS)*j + (Nx+2*NGHOSTS)*(Ny+2*NGHOSTS)*k = i + (Nx+2*NGHOSTS)*(j + (Ny+2*NGHOSTS)*k)|$$

Step 3.d: NumPy Implementation: Define the right-hand sides of the PDEs [Back to top]¶

$$\label{numpy_rhss}$$

Next, we define the right-hand sides of the $u$ and $v$ equations:

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

Again we'll approximate the $\nabla^2 u$ using fourth-order finite-difference derivatives (also see the NRPy+ tutorial on how to compute these expressions automatically or by hand using simple matrix methods).

Here we'll just use the Wikipedia article on finite-difference coefficients to construct the expressions for

$$(\nabla u)_{i,j,k} = (\partial_x^2 u)_{i,j,k} + (\partial_y^2 u)_{i,j,k} + (\partial_z^2 u)_{i,j,k}$$

by hand:

The fourth-order finite difference stencil for $(\partial_x^2 u)_{i,j,k}$ is written

$$\begin{align} (\partial_x^2 u)_{i,j,k} &= \left[-\frac{1}{12} u_{i-2,j,k} + \frac{4}{3} u_{i-1,j,k} - \frac{5}{2} u_{i,j,k} + \frac{4}{3} u_{i+1,j,k} - \frac{1}{12} u_{i+2,j,k}\right]\frac{1}{(\Delta x)^2} \\ &= \left[-\frac{1}{12} \left(u_{i-2,j,k} + u_{i+2,j,k}\right) + \frac{4}{3} \left(u_{i-1,j,k}+u_{i+1,j,k}\right) - \frac{5}{2} u_{i,j,k}\right]\frac{1}{(\Delta x)^2}, \end{align}$$ and the expressions can be written for $(\partial_y^2 u)_{i,j,k}$ and $(\partial_z^2 u)_{i,j,k}$ can be immediately written based on this pattern:

Step 3.e: NumPy Implementation: Boundary Conditions [Back to top]¶

$$\label{numpy_bcs}$$

Notice that the above code does not fill the input gridfunctions $u$ and $v$ in the ghost zones, which will be updated at each Runge-Kutta substep (as outlined next). We will need to apply our spatial boundary conditions to fill in these points. For simplicity let's choose quadratic extrapolation boundary conditions.

For example, suppose we are on the lower boundary point in $x$: $u_{1,j,k}$. Then this boundary condition will be written as the quadratic polynomial extrapolation taking data from the interior:

$$u_{1,j,k} = 3 u_{2,j,k} - 3 u_{3,j,k} + u_{4,j,k}.$$

Similarly, for the upper boundary point in $x$, the condition becomes:

$$u_{\text{Nx}-2,j,k} = 3 u_{\text{Nx}-3,j,k} - 3 u_{\text{Nx}-4,j,k} + u_{\text{Nx}-5,j,k}.$$

We'll apply this algorithm from the innermost boundary point outward, using the approach of filling in the (green-colored) ghost zones as illustrated here in 2 dimensions (courtesy Leo Werneck). Extension to 3 dimensions is straightforward.

Step 3.f: NumPy Implementation: The Method of Lines [Back to top]¶

$$\label{numpy_mol}$$

Next, we'll set up the Method of Lines (MoL) routine for Runge-Kutta fourth order (RK4), which takes the solution at a given iteration in time $n$, and enables us to advance the solution forward to iteration $n+1$, as outlined above:

$$\begin{align} k_{1,u} &= f_u(u_n,v_n) = f_u(v_n) = v_n, \\ k_{1,v} &= f_v(u_n,v_n) = f_v(u_n) = c^2\nabla^2 u_n, \\ k_{2,u} &= f_u\left(v_n + \frac{1}{2}\Delta tk_{1,v}\right) = v_n + \frac{1}{2}\Delta tk_{1,v}\\ k_{2,v} &= f_v\left(u_n + \frac{1}{2}\Delta tk_{1,u}\right) = c^2\nabla^2 \left(u_n + \frac{1}{2}\Delta tk_{1,u}\right), \\ k_{3,u} &= f_u\left(v_n + \frac{1}{2}\Delta tk_{2,v}\right) = v_n + \frac{1}{2}\Delta tk_{2,v}\\ k_{3,v} &= f_v\left(u_n + \frac{1}{2}\Delta tk_{2,u}\right) = c^2\nabla^2 \left(u_n + \frac{1}{2}\Delta tk_{2,u}\right), \\ k_{4,u} &= f_u(v_n + \Delta tk_{3,v}) = v_n + \Delta tk_{3,v}\\ k_{4,v} &= f_v(u_n + \Delta tk_{3,u}) = c^2\nabla^2 \left(u_n + \Delta tk_{3,u}\right), \\ u_{n+1} &= u_n + \frac{1}{6}\Delta t(k_{1,u} + 2k_{2,u} + 2k_{3,u} + k_{4,u}) + \mathcal{O}\big((\Delta t)^5\big)\\ v_{n+1} &= v_n + \frac{1}{6}\Delta t(k_{1,v} + 2k_{2,v} + 2k_{3,v} + k_{4,v}) + \mathcal{O}\big((\Delta t)^5\big). \end{align}$$

We will store $k_1$ through $k_4$ as additional gridfunctions, one each for $u$ and $v$, and another gridfunction for $u$ and $v$ (u_tmp and v_tmp, respectively) for the input into $f_u()$ and $f_v()$ functions:

... then implement a single timestep by calling the eval_rhs_all_interior_points() function above with appropriate inputs. Recall that the RK4 algorithm is given by

$$\begin{align} k_1 &= f(y_n, t_n), \\ k_2 &= f(y_n + \frac{1}{2}\Delta tk_1, t_n + \frac{\Delta t}{2}), \\ k_3 &= f(y_n + \frac{1}{2}\Delta tk_2, t_n + \frac{\Delta t}{2}), \\ k_4 &= f(y_n + \Delta tk_3, t_n + \Delta t), \\ y_{n+1} &= y_n + \frac{1}{6}\Delta t(k_1 + 2k_2 + 2k_3 + k_4) + \mathcal{O}\big((\Delta t)^5\big). \end{align}$$

Step 3.g: NumPy Implementation: The main driver function [Back to top]¶

$$\label{numpy_driver}$$

Finally, we'll write the main driver function, which as a diagnostic outputs the relative error between numerical and exact solutions at the closest point to the center of the numerical grid.

Step 4: Argh, the code is SLOW! Why use NRPy+ instead? [Back to top]¶

$$\label{too_slow}$$

By default the above code outputs data on the Nx=Ny=Nz=24 = $24^3$ numerical grid, and it takes around 16 seconds to complete (on mybinder).

If we were to double the resolution to $48^3$ (keeping the domain size fixed), the number of gridpoints we need to update increases by a factor of 8, and the timestep reduces by a factor of 2; hence the total cost is about 16x higher. Thus it will take roughly 16 seconds, times 16, or roughly 4 minutes to complete a $48^3$ resolution simulation on the same CPU. Similarly, it should take a bit over an hour to complete a simulation with $96^3$ resolution!

One reason we wish to use NRPy+ to convert human-friendly Python expressions (written in SymPy) to highly optimized C code is the speed. As we'll see, a C code implementing exactly the same algorithm for solving the scalar wave equation can generate results roughly 10,000x faster!

Given how slowly the above Python code solves the scalar wave equation, solving Einstein's equations of general relativity in 3 dimensions with a Python code would be a futile effort. However, speed of execution isn't the only reason to use NRPy+. Here are some more reasons:

1. NRPy+ contains a rigid syntax for SymPy symbols that
   1. enables you to specify derivatives (e.g., $\partial_j u$= u_dD[j]) and output C code at arbitrary finite differencing (FD) order
      1. tutorial on FD derivatives;
      2. tutorial on computing FD coefficients;
      3. sample C code tutorial
   2. allows for tensorial expressions to be input unambiguously in Einstein-like notation (e.g., $\gamma_{ij}=$ gammaDD[i][j])
      1. tutorial on indexed (e.g., tensorial) expressions
2. NRPy+ automatically implements 15 different RK-like timestepping methods for MoL
   1. tutorial on NRPy+ dictionary of RK methods;
   2. tutorial validating the NRPy+ RK dictionary;
   3. tutorial on C code MoL implementation
3. NRPy+ supports boundary conditions in Cartesian and singular curvilinear coordinates
   1. tutorial on general boundary condition driver;
4. NRPy+ provides support for solving the scalar wave equation stably in curvilinear coordinates, including Cartesian, spherical-like, cylindrical-like, and prolate-spheroidal-like coordinates
   1. tutorial on Cartesian scalar wave equation in NRPy+;
   2. tutorial on C code scalar wave implementation;
   3. tutorial on NRPy+ curvilinear coordinates (reference metric) support;
   4. tutorial on scalar wave equation in curvilinear coordinates
   5. Einstein Toolkit thorns
      1. tutorial on WaveToyNRPy, for solving the scalar wave equation
      2. tutorial on IDScalarWaveNRPy, for setting up scalar wave initial data
5. NRPy+ implements a covariant BSSN formulation that supports Cartesian, spherical-like, and cylindrical-like coordinates, and its boundary condition driver automatically sets up correct boundary conditions for any tensor in any orthogonal coordinate system!
   1. BSSN overview tutorial; contains links to several other tutorials
   2. Einstein Toolkit thorns
      1. tutorial on Baikal & BaikalVacuum, for solving Einstein's equations in Cartesian coordinates
6. NRPy+ contains multiple initial data sets for Einstein's equations of general relativity, and a means to quickly validate they satisfy the Einstein constraint equations on numerical grids
   1. BSSN initial data validation;
   2. static trumpet initial data;
   3. "UIUC" spinning black hole initial data;
   4. shifted Kerr-Schild initial data;
   5. Brill-Lindquist initial data;
   6. Fishbone-Moncrief black hole accretion disk initial data;
   7. piecewise-polytrope TOV initial data
7. NRPy+ contains multiple diagnostics for spacetime evolutions
   1. computing $\psi_4$ in Cartesian coordinates
   2. computing $\psi_4$ in curvilinear coordinates
   3. $\psi_4$ tetrads in curvilinear coordinates
   4. Einstein Toolkit thorn
      1. WeylScal4NRPy, a WeylScal4 clone written in NRPy+

Step 5: Error analysis & code validation: Confirming numerical errors converge to zero at the expected rate [Back to top]¶

$$\label{error_analysis}$$

The above code is really slow, so to bypass the long wait, results at $24^3$, $48^3$, and $96^3$ are precomputed and output to files in the next cell:

We say that our scheme is fourth-order-accurate in the truncation error, so the numerical solution at a given point $(t,x,y,z)$, $u_{\rm num}$, should satisfy the equation

$$u_{\rm num} = u_{\rm exact} + \mathcal{O}(\Delta x^4) + \mathcal{O}(\Delta t^4),$$

where $u_{\rm exact}$ is the exact solution, and $\mathcal{O}(\Delta x^4)$ and $\mathcal{O}(\Delta t^4)$ are terms proportional to $\Delta x^4$ and $\Delta t^4$, respectively. However note that the CFL condition for this PDE requires that $\Delta x \propto \Delta t$, so we can simplify the above expression to

$$u_{\rm num} = u_{\rm exact} + \mathcal{O}(\Delta x^4).$$

Therefore, the relative error between the numerical and the exact solution should be given to good approximation by

$$\begin{align} E_{\rm Rel} &= \left| \frac{u_{\rm num} - u_{\rm exact}}{u_{\rm exact}}\right| \\ &\propto \Delta x^4 \\ \implies E_{\rm Rel} &= k \Delta x^4, \end{align}$$

where $k$ is the proportionality constant, divided by $u_{\rm exact}$.

Therefore, taking the logarithm of both sides of the equation, we get:

$$\begin{align} \log_{10} E_{\rm Rel} &= \log_{10} (k [\Delta x]^4) \\ &= \log_{10} ([\Delta x]^4) + \log_{10} (k) \\ &= 4 \log_{10} (\Delta x) + \log_{10} (k) \end{align}$$

$\Delta x$ is proportional to 1/Nx, so if we perform the simulation at twice the resolution (i.e., $\Delta x\to \Delta x/2$), $\log_{10} E_{\rm Rel}$ should drop by

$$\begin{align} 4 \log_{10} (\Delta x) - 4 \log_{10} (\Delta x/2) &= 4 \log_{10} \frac{\Delta x}{\Delta x/2} \\ &= 4 \log_{10} 2 \approx 1.20. \end{align}$$

In the below plot we show that when the logarithmic relative error $\log_{10} E_{\rm Rel}$ versus time in the $48^3$ case is shifted upward by 1.2 (or for the $96^3$ case by 2.4), we observe perfect overlap with $\log_{10} E_{\rm Rel}$ in the $24^3$ case (except at $t=0$ when the numerical solution is set to the exact solution). This is a common way in numerical relativity to present the convergence of numerical errors to zero, demonstrating that our code is working as expected.

Step 6: Exercises for students [Back to top]¶

$$\label{student_exercises}$$

1. Adjust the above code to make it run twice as fast on the same numerical grid while generating exactly the same results (stored in the files above). Bonus: Can you make it run any faster than twice as fast (while still being written in "pure" Python, using NumPy)?
2. How much should the (absolute value of the) relative error |rel_error| drop if we were to double the resolution, assuming instead 6th-order convergence? How will this adjust the log10(|rel_error|)?
3. Why did we add a nonzero constant offset to the exact solution? (Hint: Start from the definition of relative error.)
4. What will happen to the convergence order if we continue the simulation for a much longer time, say to $t=2$? Why? (Hint: Convergence order will plummet!)

Step 7: 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-Solving_the_Scalar_Wave_Equation_with_NumPy.pdf. (Note that clicking on this link may not work; you may need to open the PDF file through another means.)