Start-to-Finish Example: Applying Boundary Conditions in Curvilinear Coordinates in Three Dimensions¶

Author: Zach Etienne¶

This module documents and validates basic boundary condition algorithms for curvilinear coordinate systems (e.g., Spherical, Cylindrical), based on prescription in the SENR/NRPy+ paper.¶

Introduction: Applying boundary conditions in curvilinear coordinates¶

Here is a list of four challenges we face when solving PDEs in curvilinear coordinates:

1. Challenge 1: Unlike ordinary Cartesian coordinate grids, not all boundary points on uniform curvilinear coordinate grids are outer boundary points.
2. Challenge 2: Figuring out the locations to which boundary points map requires a coordinate inversion, but some coordinate systems are not easily invertible.
3. Challenge 3: Tensors and vectors in curvilinear coordinates can change directions across "inner" boundaries, changing sign as a result.
4. Challenge 4: Coordinate singularities can appear, causing a stiffening or divergence of terms in PDEs.

Challenge 1: Inner versus outer boundary points¶

Unlike Cartesian coordinates, the boundaries of our grid in generic curvilinear coordinates are not all outer boundaries.

Consider first a computational grid in Cartesian coordinates, with $(x_0,x_1,x_2)=(x,y,z)$, that is uniform (i.e., with $\Delta x$, $\Delta y$, and $\Delta z$ all set to constant values). This grid may extend over arbitrary coordinate ranges $x_i \in [x_{i, \rm min},x_{i, \rm max}]$.

By contrast, consider now a uniform grid in spherical coordinates $(x_0,x_1,x_2)=(r,\theta,\phi)$ with constant spacing $\Delta r$, $\Delta \theta$, and $\Delta \phi$ between grid points in $r$, $\theta$, and $\phi$, respectively. Further, let's assume that these grids span all possible values of $\theta$ and $\phi$, with $r=0$ included in the domain. Then our numerical domain must satisfy the following relations

• $x_0 = r \in [0,{\rm RMAX}]$,
• $x_1 = \theta \in [0,\pi]$, and
• $x_2 = \phi \in [-\pi,\pi]$. (Notice how we do not choose $x_2= \phi \in [0,2\pi]$ so that our conversion from Cartesian to spherical coordinates is compatible with the output range from the ${\rm atan2}(y,x)$ function: $\phi={\rm atan2}(y,x)\in[-\pi,\pi]$.)

Notice that unlike Cartesian coordinates, the boundaries of this numerical grid in spherical coordinates are not all outer boundaries. For example, data stored at $\phi=-\pi$ will be identical to data at $\phi=\pi$, regardless of $r$ or $\theta$. I.e., $\phi$ satisfies periodic boundary conditions only. Further, $\theta=0$ presents a more complicated boundary condition, in which points with negative $\theta$ map to points with $|\theta|$ but at an angle of $\phi\to \phi+\pi$. Finally, negative $r$ points will map to postive $r$ points on the other side of the origin. We call these boundaries inner boundaries, as they generally map to other points in the interior (as opposed to the outer boundaries) of the grid.

As we generally cannot apply an outer boundary condition to the inner boundaries, these boundaries will need to be treated differently.

On our numerical grids, this poses some difficulty, as finite difference derivatives we compute within the numerical domain require that the grid be extended beyond the domain boundaries. In spherical coordinates, this means that we need additional grid points at, e.g., $r<0$, $\theta<0$, and $\phi>\pi$, just to name a few. Whether they be on outer or inner boundaries, we call grid points in the extended region ghost zones.

Numerical grids of $N$th order accuracy generally possess $N/2$ ghost zone points in the boundary regions (i.e., $x_i < x_{i,\rm min}$ and $x_i > x_{i, \rm max}$). While in Cartesian coordinates, these ghost zone points map to regions outside the grid domain $x_i \in [x_{i, \rm min},x_{i, \rm max}]$, in spherical coordinates, most ghost zone points map to regions inside the grid domain. For example, for some $\tilde{r}\in [0,{\rm RMAX}]$ and $\tilde{\theta}\in[0,\pi]$, the ghost zone point $(\tilde{r},\tilde{\theta},2\pi+\Delta \phi/2)$ would map to the interior point $(\tilde{r},\tilde{\theta},\Delta \phi/2)$ because the $\phi$ coordinate is periodic. Thus when given a ghost zone point in some arbitrary curvilinear coordinate system, we are faced with the problem of addressing the following two questions:

1. Does a given ghost point map to an interior point, or is it an outer boundary point (i.e., a point exterior to the domain)?
2. If the ghost zone point maps to an interior point, to which interior point does it map?

Challenge 2: Inverting coordinates¶

Coordinate systems within NRPy+ are generally Spherical-like, Cylindrical-like, SymTP-like (where SymTP is a prolate-spheroidal coordinate system), or Cartesian-like. For example, SinhSphericalv2 coordinates are exactly the same as Spherical coordinates, except we choose an odd function for the radial coordinate $r$ as a function of $x_0$:

While this coordinate choice exhibits nice properties for certain cases, the function $x_0(r)$ is not a closed-form expression. Thus finding the mapping of ghost zone points in the radial direction would require a root finder.

Is there an easier way of dealing with this problem than with a root finder?

Challenge 3: Parity: changes of direction in vectors and tensors across inner boundaries¶

When applying inner boundary conditions to vectors and tensors, we must consider how the direction or parity of vector and tensor components change across the inner boundary.

Suppose we have a vector $v^\rho$ defined at ghost zone $(-\rho,\phi,z)$ ($\rho>0$) in cylindrical coordinates. This will map to an interior point at $(\rho,\phi+\pi,z)$. At this point, the direction of the $\hat{\rho}$ unit vector flips sign. Thus we cannot simply set the value of $v^\rho$ to the value it possesses at interior point $(\rho,\phi+\pi,z)$; that would result in a sign error. Instead we have \begin{align} v^\rho(-\rho,\phi,z)&=-v^\rho(\rho,\phi+\pi,z) \\ &= \mathbf{e}^\rho\left(-\rho,\phi,z\right) \cdot \mathbf{e}^\rho\left(\rho,\phi+\pi,z\right)v^\rho(\rho,\phi+\pi,z), \end{align} where $\mathbf{e}^\rho\left(\rho,\phi,z\right)$ is the $\rho$ unit vector evaluated at point $(\rho,\phi,z)$, and $\mathbf{e}^\rho\left(-\rho,\phi,z\right) \cdot \mathbf{e}^\rho\left(\rho,\phi+\pi,z\right)$ is the dot product of the two unit vectors, which must evaluate to $\pm 1$ (i.e., the parity). Contrast this with scalars, which do not possess a sense of direction/parity.

Challenge 4: Coordinate singularities¶

Most non-Cartesian, orthogonal coordinate systems (like spherical coordinates) possess coordinate singularities.

For example, coordinate singularities in spherical coordinates lie along $\theta=0$ and $\theta=\pi$; these are points where the coordinate system focuses to a single point. For example, the coordinate singularity at the North Pole is the reason why all directions are south there. Critically, these singularities manifest as points where the reference metric or its inverse crosses through zero or diverges to $\infty$. As we derived in a previous module, the Laplacian in spherical polar coordinates takes the form $$ \nabla^2 u = \partial_r^2 u + \frac{1}{r^2} \partial_\theta^2 u + \frac{1}{r^2 \sin^2 \theta} \partial_\phi^2 u + \frac{2}{r} \partial_r u + \frac{\cos\theta}{r^2 \sin\theta} \partial_\theta u, $$ which diverges at $r=0$ and $\sin\theta=0-$precisesly at the $\theta=0$ and $\theta=\pi$ coordinate singularity.

To avoid this divergence, we simply choose that our numerical grids be cell-centered.

I.e., given the desired bounds of the grid interior to be

${\rm NGHOSTS}$ to be the number of ghost zones (assumed the same in all directions), and $\{N_0,N_1,N_2\}$ to be the desired number of points in the grid interior in the $\{x_0,x_1,x_2\}$ directions, respectively, then the numerical grid spacing in each respective direction will be given by

Given the above definitions, the complete set of indices $\{{\rm i0},{\rm i1},{\rm i2}\}$ located at $\{x_{0,{\rm i0}},x_{1,{\rm i1}},x_{2,{\rm i2}}\}$ as follows:

where

• ${\rm i0}\in [0,N_0+2\cdot{\rm NGHOSTS})$
• ${\rm i1}\in [0,N_1+2\cdot{\rm NGHOSTS})$
• ${\rm i2}\in [0,N_2+2\cdot{\rm NGHOSTS})$,

which guarantees the interior is covered by exactly $\{N_0,N_1,N_2\}$ grid points, the boundaries are covered by ${\rm NGHOSTS}$ ghost zones, and we maintain cell centering.

So for example, if we choose a numerical grid in spherical coordinates $\{r,\theta,\phi\}$, with 3 ghost zone points (needed for e.g., 6th-order-accurate centered finite differencing), and we want the grid interior to be sampled with $\{N_r,N_\theta,N_\phi\}$ grid points, then we have

where again

• ${\rm i0}\in [0,N_r+2\cdot{\rm NGHOSTS})$
• ${\rm i1}\in [0,N_\theta+2\cdot{\rm NGHOSTS})$
• ${\rm i2}\in [0,N_\phi+2\cdot{\rm NGHOSTS})$,

which guarantees the interior is covered by exactly $\{N_r,N_\theta,N_\phi\}$ grid points, the boundaries are covered by ${\rm NGHOSTS}$ ghost zones, and we maintain cell centering.

Notice that in NRPy+, we use the physics notation for spherical coordinates, where $\theta$ is the polar and $\phi$ is the azimuthal angle. Also we choose $\phi$ to range from $-\pi$ to $+\pi$, which is most useful since it is compatible with output from atan2.

Exercise to student: Given the prescription above, why do the integers $N_\theta$ and $N_\phi$ need to be even?

As Laplacians like these appear on the right-hand sides of, e.g., the scalar wave equation in curvilinear coordinates, we still have a problem of some terms becoming quite large as the coordinate singularity is approached. This issue manifests as a stiffening of the PDE, requiring that we be very careful about the precise Method of Lines timestepping algorithm used. See Cordero-Carrión & Cerdá-Durán for information on dealing with this subtlety in a second-order Runge-Kutta Method of Lines context; it was later found that the standard RK4 method maintain stable solutions to PDEs affected by this sort of stiffening.

The above discussion focuses primarily on scalar fields. However, when solving PDEs involving vectors and tensors, the vectors and tensors themselves can exhibit divergent behavior at coordinate singularities. The good news is, this singular behavior is well-understood in terms of the scale factors of the reference metric, enabling us to define rescaled version of these quantities that are well behaved (so that, e.g., they can be finite differenced).

For example, given a smooth vector in a 3D Cartesian basis $\bar{\Lambda}^{i}$, all components $\bar{\Lambda}^{x}$, $\bar{\Lambda}^{y}$, and $\bar{\Lambda}^{z}$ will be smooth (by assumption). When changing the basis to spherical coordinates (applying the appropriate Jacobian matrix transformation), we will find that since $\phi = \arctan(y/x)$, $\bar{\Lambda}^{\phi}$ is given by

Thus $\bar{\Lambda}^{\phi}$ diverges at all points where $r\sin\theta=0$ (or equivalently where $x=y=0$; i.e., the $z$-axis) due to the $\frac{1}{r\sin\theta}$ that appear in the Jacobian transformation.

This divergence might pose no problem on cell-centered grids that avoid $r \sin\theta=0$, except that the BSSN equations require that first and second derivatives of quantities like $\bar{\Lambda}^{\phi}$ be computed. Usual strategies for numerical approximation of these derivatives (e.g., finite difference methods) will "see" these divergences and errors generally will not drop to zero with increased numerical sampling of the functions at points near where the functions diverge.

However, notice that if we define $\lambda^{\phi}$ such that

then $\lambda^{\phi}$ will be smooth and non-divergent as well. The strategy when computing derivatives of $\bar{\Lambda}^{\phi}$ therefore is to perform the product rule on the above expression, computing derivatives of the scale factors analytically (i.e., exactly using a computer algebra system like SymPy), and smooth terms like $\lambda^{\phi}$ with finite-difference derivatives.

Avoiding such singularities can be generalized to arbitrary coordinate systems, so long as $\lambda^i$ is defined as:

where scalefactor[i] is the $i$th scale factor in the given coordinate system. This idea can be extended to covariant (lowered-index) vectors and arbitrary tensors, as described in the BSSN quantities tutorial notebook.

In summary, Challenge 4 is addressed via a combination of cell-centered grids, tensor rescaling, and a stable Method of Lines time stepping algorithm. This tutorial notebook will therefore focus on addressing Challenges 1 through 3, which, coincidentally, are addressed via an appropriate boundary condition algorithm.

Table of Contents¶

This notebook is organized as follows

1. Step 1: Overview of boundary condition algorithm in curvilinear coordinates
   1. Step 1.a: Addressing Challenge 1: Distinguishing inner from outer boundary points
   2. Step 1.b: Addressing Challenge 2: Eigen-Coordinate systems
   3. Step 1.c: Addressing Challenge 3: Applying parity conditions to arbitrary-rank tensors
2. Step 2: bc_struct and other C data structures used for storing boundary condition information
3. Step 3: NRPy+-based C code generation for parity conditions
   1. Step 3.a: Set up unit-vector dot products (=parity) for each of the 10 parity condition types, store to parity_conditions_symbolic_dot_products.h
   2. Step 3.b: Set parity type for each gridfunction, based on the digits at the end of its name, output to dirname+gridfunction_defines.h
4. Step 4: set_up__bc_gz_map_and_parity_condns(): C function for distinguishing inner from outer boundary points, and setting parity conditions
5. Step 5: set_bcstruct(): Using information from set_up__bc_gz_map_and_parity_condns() as input, set bcstruct
6. Step 6: driver_bcstruct.h: C code driver for declaring bc_struct data type, the bcstruct instance of said data type, and calling set_up__bc_gz_map_and_parity_condns() and set_bcstruct() to fill bcstruct
7. Step 7: apply_bcs_curvilinear() C function: quickly apply boundary and parity conditions with information from bcstruct
8. Step 8: CurviBC_Playground.c: Start-to-Finish C code module for testing & validating curvilinear boundary conditions
   1. Step 8.a: Generate C code for generating numerical grids
   2. Step 8.b: Register gridfunctions of all 10 parity types in NRPy+; output gridfunction aliases to CurviBoundaryConditions/gridfunction_defines.h
   3. Step 8.c: Set up test data for Curvilinear Boundary Conditions code validation
   4. Step 8.d: CurviBC_Playground.c : The Main C Code
9. Step 9: Validation: Compare with original SENR results
10. Step 10: Output this notebook to $\LaTeX$-formatted PDF file

Step 1: Overview of boundary condition algorithm in curvilinear coordinates [Back to top]¶

Here we review the basic algorithm for addressing Challenges 1 2, and 3 discussed in the Introduction above.

The algorithm itself is implemented as C code in Steps 2 (data structures), 3 (searching entire grid for inner and outer boundary points, and setting parities), 4 (setting data structures for quick and efficient implementation of outer boundaries), and 5 (function to apply inner & outer boundary conditions).

Step 1.a: Addressing Challenge 1: Distinguishing inner from outer boundary points [Back to top]¶

At each ghost zone grid point $\mathbf{d}_{\rm gz}=(x_0,x_1,x_2)$, we will do the following:

1. Evaluate the Cartesian coordinate $\left(x(x_0,x_1,x_2),y(x_0,x_1,x_2),z(x_0,x_1,x_2)\right)$, corresponding to this grid point. Then evaluate the inverse $\mathbf{d}_{\rm new}=\left(x_0(x,y,z),x_1(x,y,z),x_2(x,y,z)\right)$.
   1. If $\mathbf{d}_{\rm new} \ne \mathbf{d}_{\rm gz}$, then the ghost zone grid point maps to a point in the grid interior, which is exactly the case described in the above section. To distinguish this case from an "outer boundary condition", we shall henceforth refer to it variously as an application of an "interior", "inner", or "parity" boundary condition.
   2. Ghost zone points for which $\mathbf{d}_{\rm new} \equiv \mathbf{d}_{\rm gz}$ are on the outer boundary of the grid, and standard outer boundary conditions should be applied.

In detail, the algorithm is as follows:

1. Convert the coordinate $(x_0,x_1,x_2)$ for the ghost zone point to Cartesian coordinates $\left(x(x_0,x_1,x_2),y(x_0,x_1,x_2),z(x_0,x_1,x_2)\right)$. For example, if we choose ordinary spherical coordinates $(x_0,x_1,x_2)=(r,\theta,\phi)$, then

   • $x(r,\theta,\phi) = r \sin(\theta) \cos(\phi) = x_0 \sin(x_1) \cos(x_2)$
   • $y(r,\theta,\phi) = r \sin(\theta) \sin(\phi) = x_0 \sin(x_1) \sin(x_2)$
   • $z(r,\theta,\phi) = r \cos(\theta) = x_0 \cos(x_1)$

2. Once we have $(x,y,z)$, we then find the corresponding value $(x_0,x_1,x_2)_{\rm in/OB}=(r,\theta,\phi)_{\rm in/OB}$ in the grid interior or outer boundary, via the simple inverse formula:

   • $r_{\rm in/OB} = x_{0, \rm in/OB} = \sqrt{x^2+y^2+z^2} \in [0,\infty)$
   • $\theta_{\rm in/OB} = x_{1, \rm in/OB} = {\rm acos}\left(\frac{z}{\sqrt{x^2+y^2+z^2}}\right) \in [0,\pi]$
   • $\phi_{\rm in/OB} = x_{2, \rm in/OB} = {\rm atan2}(y,x) \in [-\pi,\pi]$ Wikipedia article on atan2()

3. If $(x_0,x_1,x_2)_{\rm in/OB}$ is the same as the original $(x_0,x_1,x_2)$, then we know $(x_0,x_1,x_2)$ is an outer boundary point (in spherical coordinates, at $r>{\rm RMAX}$), and we store (i0,i1,i2)$_{\rm in/OB} = (-1,-1,-1)$. Otherwise, we know that $(x_0,x_1,x_2)$ maps to some interior point at index (i0,i1,i2), which we store:

   • $\rm{i0}_{\rm in/OB}=\frac{r_{\rm in/OB} - r_{\rm min}}{\Delta r} - \frac{1}{2}$
   • $\rm{i1}_{\rm in/OB}=\frac{\theta_{\rm in/OB} - \theta_{\rm min}}{\Delta \theta} - \frac{1}{2}$
   • $\rm{i2}_{\rm in/OB}=\frac{\phi_{\rm in/OB} - \phi_{\rm min}}{\Delta \phi} - \frac{1}{2}$

4. When updating a ghost zone point (i0,i1,i2) in the domain exterior, if the corresponding (i0,i1,i2)$_{\rm in/OB}$ was set to $(-1,-1,-1)$, then we apply outer boundary conditions. Otherwise, we simply copy the data from the interior point at (i0,i1,i2)$_{\rm in/OB}$ to (i0,i1,i2).

Following the prescription in the SENR/NRPy+ paper, we will implement curvilinear boundary conditions for rank-0, rank-1, and symmetric rank-2 tensors in three dimensions; as this is the same dimension and highest rank needed for BSSN.

Step 1.b: Addressing Challenge 2: Eigen-Coordinate Systems [Back to top]¶

Suppose we were to rewrite the spherical coordinate $r$ as an arbitrary odd function of $x_0$ instead of $x_0$ itself. In that case, $r(-x_0)=-r(x_0)$, and all parity conditions remain unchanged. However the inverse function, $x_0(r)$, may not be writable as a closed-form expression, requiring a Newton-Raphson root finder to find the appropriate boundary mappings.

To greatly simplify the algorithm in the case of arbitrary $r(x_0)$ in Spherical-like coordinates, or $\rho(x_0)$ or $z(x_2)$ in Cylindrical-like coordinates, we note that the coordinate mappings and parities for all Spherical-like coordinate systems are identical to the mappings and parities for ordinary Spherical coordinates. The same holds true for Cylindrical-like and SymTP-like coordinate systems. Thus so long as we know the correct "Eigen-Coordinate system" (i.e., Spherical in the case of SinhSpherical or SinhSphericalv2; Cylindrical in the case of SinhCylindrical; SymTP in the case of SinhSymTP; etc.) there is no need for a Newton-Raphson root finder to set up the boundary conditions.

Step 1.c: Addressing Challenge 3: Applying parity conditions to arbitrary-rank tensors [Back to top]¶

Above we presented the strategy for applying parity boundary conditions to a single component of a vector. Here we outline the generic algorithm for arbitrary-rank tensors.

Continuing the discussion from the previous section, we assume $\mathbf{d}_{\rm new} \ne \mathbf{d}_{\rm gz}$ (otherwise we would apply the outer boundary condition algorithm). Next suppose we are given a generic rank-$N$ tensor ($N>0$).

1. The first component of the rank-$N$ tensor corresponds to some direction with unit vector $\mathbf{e}^i$; e.g., $v^r$ corresponds to the $\mathbf{e}^r$ direction. Compute the dot product of the unit vector $\mathbf{e}^i$ evaluated at points $\mathbf{d}_{\rm gz}$ and $\mathbf{d}_{\rm new}$. Define this dot product as $P_1$ ("$P$" for "parity"):

1. $P_1$ will take the value of $\pm 1$, depending on the unit-vector direction and the points $\mathbf{d}_{\rm gz}$ and $\mathbf{d}_{\rm new}$
2. Repeat the above for the remaining components of the rank-$N$ tensor $j\in \{2,3,...,N\}$, storing each $P_j$.
3. The tensor mapping from $\mathbf{d}_{\rm gz}$ to $\mathbf{d}_{\rm new}$ for this tensor $T^{ijk...}_{mnp...}$ will be given by

In this formulation of BSSN, we only need to deal with rank-0, rank-1, and symmetric rank-2 tensors. Further, our basis consists of 3 directions, so there are a total of

• 1 parity condition (the trivial +1) for scalars (rank-0 tensors)
• 3 parity conditions for all rank-1 tensors (corresponding to each direction)
• 6 parity conditions for all symmetric rank-2 tensors (corresponding to the number of elements in the lower or upper triangle of a $3\times3$ matrix, including the diagonal)

Thus we must keep track of the behavior of 10 separate parity conditions, which can be evaluated once the numerical grid has been set up, for all time. The following Table outlines the correct conditions for each:

The appropriate dot products determining parity condition are assigned to each gridfunction based on the following numbering:

In the following few steps, we will document the data structures for and implementation of this boundary condition algorithm.

Step 2: bc_struct: Data structure for storing boundary condition information [Back to top]¶

Here we define bc_struct, a C data structure that stores all information needed to apply boundary conditions at all boundary points.

The information needed to fill a ghost zone depends on whether it exists on an inner or an outer boundary, and the bc_struct data structure is constructed to account for this fact:

typedef struct __bcstruct__ {
    outer_bc **outer; // Array of 1D arrays, of length 
                      //   [NGHOSTS][num_ob_gz_pts[which_outer_ghostzone_point]]

    inner_bc **inner; // Array of 1D arrays, of length 
                      //   [NGHOSTS][num_ib_gz_pts[which_inner_ghostzone_point]]

    // Arrays storing number of outer/inner boundary ghostzone points at each ghostzone,
    //   of length NGHOSTS:
    int     *num_ob_gz_pts; 
    int     *num_ib_gz_pts;
} bc_struct;

Ghost zones must be filled in from the inside outward, as e.g., the outermost ghost zones may depend on ghost zones closer to the grid interior being set. We thus store ghost zones in arrays that point to a particular layer of ghost zones. This explains the fact that we declare inner and outer in bc_struct with the ** prefix, denoting that these are "pointers to pointers".

For example outer[0] points to the set of ghost zone cells on the outer boundary that are in the layer of ghost zones just adjacent to the grid interior. Further, outer[0][0] points to a outer_bc struct containing all information needed to fill "outer boundary point zero" with valid ata. Note that the numbering of the boundary points (the j index in outer[i][j]) is rather arbitrary, but each point has a unique label, and there are no duplicates. This ensures efficiency and locality in memory.

The same reasoning holds when considering ghost zones that are not outer boundary points.

Next we summarize all basic data structures that appear within bc_struct.

• inner_bc:

typedef struct __inner_bc__ {
    gz_map inner_bc_dest_pt;
    gz_map inner_bc_src_pt;
    int8_t parity[10]; // We store the 10 parity conditions in 10 int8_t integers, 
                       // one for each condition. Note that these conditions can 
                       // only take one of two values: +1 or -1, hence the use of 
                       // int8_t, the smallest C data type.
} inner_bc;

* `inner_bc_dest_pt.{i0,i1,i2}`: Location `(i0,i1,i2)` of inner ghost zone point `which_pt` on the `which_gz` ghost zone layer.
* `inner_bc_src_pt.{i0,i1,i2}`: Location of the interior grid point to which the `inner_bc_dest_pt.{i0,i1,i2}` inner ghost zone maps.
* `parity[10]` Parity information ($\pm 1$) at the given inner ghost zone, for all 10 gridfunction parity types.

• outer_bc:

typedef struct __outer_bc__ {
    gz_map outer_bc_dest_pt;
    int8_t FACEi0,FACEi1,FACEi2; // FACEi* takes values of -1, 0, and +1 only,
                                 // corresponding to MAXFACE, NUL, and MINFACE
                                 // respectively.
                                 // Thus int8_t (one byte each, the smallest C 
                                 // type) is sufficient.
} outer_bc;

* outer_bc_dest_pt.{i0,i1,i2}`: Location `(i0,i1,i2)` of outer ghost zone point `which_pt` on the `which_gz` ghost zone layer 
* `int8_t FACEi0,FACEi1,FACEi2`: Specifies to which face of the numerical domain outer boundary point `which_pt` on the `which_gz` ghost zone layer corresponds. Many outer boundary conditions depend on some extrapolation from inner points. Thus knowing on which face a given outer boundary point lies provides needed direction for extrapolation.

• num_ib/ob_gz_pts[which_gz]: The number of inner/outer boundary points on ghost zone layer which_gz

Step 3: NRPy+-based C code generation for parity conditions [Back to top]¶

Much of the algorithm needed for setting up bcstruct requires a loop over all gridpoints on the numerical grid. As the precise numerical grids are chosen at C runtime, that part of the algorithm must be run entirely within a static C code.

However, there are two parts to the overall algorithm that must be generated by NRPy+, namely

1. Step 3.a: parity_conditions_symbolic_dot_products.h: Based on the chosen reference metric, sets up the needed unit-vector dot products for each of the 10 parity condition types.
2. Step 3.b: Set parity type for each gridfunction registered within NRPy+, based on the digits at the end of each gridfunction name, append result to dirname+gridfunction_defines.h

Step 3.a: Set up unit-vector dot products (=parity) for each of the 10 parity condition types, store to parity_conditions_symbolic_dot_products.h [Back to top]¶

Next we generate the C code necessary to perform needed dot products for filling in the parity condition arrays inside bcstruct.

Using the unit vectors defined in rfm.UnitVectors[][] (in reference_metric.py), each unit vector takes as input either $\mathbf{d}_{\rm gz} = (x_0,x_1,x_2)_{\rm IB}$=(xx0,xx1,xx2) or $\mathbf{d}_{\rm new} = (x_0,x_1,x_2)_{\rm in}$=(xx0_inbounds,xx1_inbounds,xx2_inbounds) as summarized in the table above in Step 1. We paste the table here again, for quick reference:

The appropriate symbolic dot products determining parity condition are assigned to each gridfunction based on the following numbering:

Looping over all 10 parity types, the corresponding symbolic expressions for dot product(s) is output to C code. For example, in Spherical coordinates, parity type 1's dot product output as C code is given by:

parity[1] = sin(xx1)*sin(xx1_inbounds)*sin(xx2)*sin(xx2_inbounds) + sin(xx1)*sin(xx1_inbounds)*cos(xx2)*cos(xx2_inbounds) + cos(xx1)*cos(xx1_inbounds);

To wit (as described above), there are 10 parity types for BSSN evolved variables, which include tensors up to and including rank-2.

Step 3.b: Set parity type for each gridfunction, based on the digits at the end of its name, output to dirname+gridfunction_defines.h [Back to top]¶

For example, if the gridfunction name ends with "01", then (based on the table above) the set_parity_types() function below will set the parity_type of that gridfunction to 5. We can be assured this is a robust algorithm, because gri.register_gridfunctions() in grid.py will throw an error if a gridfunction base name ends in an integer. This strict syntax was added with the express purpose of making it easier to set parity types based on the gridfunction name.

After each parity type is found, we store the parity type of each gridfunction to const int8_t arrays evol_gf_parity and aux_gf_parity, appended to the end of gridfunction_defines.h.

Step 4: set_up__bc_gz_map_and_parity_condns(): C function for distinguishing inner from outer boundary points, and setting parity conditions [Back to top]¶

This step is performed using only the Eigen-Coordinate corresponding to the chosen reference_metric::CoordSystem.

set_up__bc_gz_map_and_parity_condns() loops over all points on the numerical grid (interior and ghost zone points included), with the aim of filling gz_map *bc_gz_map and parity_condition *bc_parity_conditions at each point. To do so, the function implements the algorithm described above in Step 1.

That is to say, at each coordinate point $(x_0,x_1,x_2)$ situated at grid index (i0,i1,i2):

1. Convert the curvilinear coordinate $(x_0,x_1,x_2)$ to the corresponding Cartesian coordinate $(x,y,z)$
2. Find the Cartesian grid point (i0_inbounds,i1_inbounds,i2_inbounds) in the grid interior or outer boundary corresponding to this Cartesian coordinate $(x,y,z)$.
3. If and only if we are on an outer boundary ghost zone or in the grid interior, i0_inbounds==i0, i1_inbounds==i1, and i2_inbounds==i2, and inner boundary conditions do not apply: set bc_gz_map to $(-1,-1,-1)$, and for all 10 gridfunction parities, set parity=1.
4. If i0_inbounds==i0, i1_inbounds==i1, and i2_inbounds==i2, does not hold true, then (i0,i1,i2) is an inner boundary point: set bc_gz_map to (i0_inbounds,i1_inbounds,i2_inbounds), and for all 10 gridfunction parities, evaluate all dot products of the unit vectors evaluated at (i0_inbounds,i1_inbounds,i2_inbounds) and (i0,i1,i2), as prescribed above in Step 1. The C code for computing the needed symbolic dot products is generated above in Step 2, and is imported via an #include statement in below C code function eval_symbolic_dot_products_to_set_parity_conditions().

Step 5: set_bcstruct(): Using information from set_up__bc_gz_map_and_parity_condns() as input, set bc_struct [Back to top]¶

As described above, set_up__bc_gz_map_and_parity_condns() sets gz_map *bc_gz_map and parity_condition *bc_parity_conditions at all grid points. While this information could be used directly to apply boundary conditions, it is more efficient (both in memory and CPU time) to instead store only the needed information to the bcstruct array, so that when applying boundary conditions, we can simply loop over bcstruct.

The algorithm follows the bc_struct data type defined above in Step 2, and loops over all boundary ghost zones, from the innermost layer (which_gz=0) outward (to which_gz=NGHOSTS-1). This is necessary, as, for example, some ghost zone points on the which_gz=1 layer depend on ghost zones being set on the which_gz=0 layer.

1. Count the number of outer and inner boundary points, store to num_ob_pts and num_ib_pts, respectively.
2. Now that we know the number of outer boundary points on this ghostzone layer, allocate memory needed for storing the outer and inner boundary condition data.
   1. At all outer boundary ghost zones, allocate memory a single member of the outer_bc data type.
   2. At all inner boundary ghost zones, allocate memory a single member of the inner_bc data type.
3. Store the number of outer and inner boundary points on each ghostzone layer, where e.g.,which_gz==0 corresponds to the innermost ghostzones on the numerical domain.
4. Store information needed for outer boundary conditions, to outer_bc_dest_pt and outer_bc_face arrays.
5. Store information needed for inner boundary conditions, including interior point to which inner ghost zone maps, and parity conditions for all 10 gridfunction types.

Step 6: driver_bcstruct.h: C code driver for declaring bc_struct data type, the bcstruct instance of said data type, and calling set_up__bc_gz_map_and_parity_condns() and set_bcstruct() to fill bcstruct [Back to top]¶

CurviBoundaryConditions/driver_bcstruct.h is meant to be #include'd in a C main() function. It allocates memory for boundary condition pointers and calls functions referenced above, including set_up__bc_gz_map_and_parity_condns() and set_bcstruct(). It then frees unneeded memory after bcstruct has been fully initialized.

For completeness, we also output

• CurviBoundaryConditions/bcstruct_freemem.h, which frees memory allocated within bcstruct. It is meant to be #include'd at the end of a C main() function, and
• CurviBoundaryConditions/CurviBCs_Cfunctions.h, which #include's all the needed C source functions in order.

Step 7: apply_bcs_curvilinear() C function: quickly apply boundary and parity conditions with information from bcstruct [Back to top]¶

apply_bcs_curvilinear() loops over all NUM_GFS gridfunctions in the gfs IDX4 array and, using a bcstruct filled in by the set_bcstruct() function above, applies boundary conditions to each ghost zone layer, starting with the innermost layer and working outward.

This function is meant to be called within a Method of Lines timestepping algorithm, at or near the end of each substep.

Step 8: CurviBC_Playground.c: Start-to-Finish C code module for testing & validating curvilinear boundary conditions [Back to top]¶

Step 8.a: Generate C code for generating numerical grids [Back to top]¶

Next we generate the C code needed for applying boundary conditions in generic coordinate systems, using the Eigen-Coordinate approach described above:

1. $\left(x(x_0,x_1,x_2),y(x_0,x_1,x_2),z(x_0,x_1,x_2)\right)$, (CurviBoundaryConditions/EigenCoord_xx_to_Cart.h)
2. $\left(x_0(x,y,z),x_1(x,y,z),x_2(x,y,z)\right)$, (CurviBoundaryConditions/EigenCoord_Cart_to_xx.h):

Step 8.b: Register gridfunctions of all 10 parity types; output gridfunction aliases to CurviBoundaryConditions/gridfunction_defines.h [Back to top]¶

Here we

1. Register within NRPy+ one gridfunction per each of the 10 parity conditions, and then
2. output to CurviBoundaryConditions/gridfunction_defines.h the corresponding gridfunction aliases and parity info, so the C code can
   1. Access each gridfunction by its human-friendly alias (e.g., test_gfs[RANKONEU0GF][idx] instead of test_gfs[6][idx]), and
   2. Access each gridfunction parity by the same alias (e.g., bcstruct->inner_bc_parity[which_gz][pt].parity[gf_parity[RANKONEU0GF]])

Step 8.c: Set up test data for Curvilinear Boundary Conditions code validation [Back to top]¶

We will validate this curvilinear boundary condition module by comparing its results with the original (trusted) SENR code, as follows:

• Discrete data test:
  1. Fill all 10 gridfunctions at each gridpoint with the unique gridpoint integer index IDX3S(i0,i1,i2)
  2. Apply curvilinear boundary conditions
  3. Compare output data at all gridpoints with those from the original SENR code. Agreement should be perfect.

Another (future, to-be-implemented) test, which will enable us to validate coordinate systems that do not exist within the original SENR code, is described below:

• Smooth data test (TODO):
  1. Fill all 10 gridfunctions with data that are smooth in the Cartesian basis.
  2. Apply Jacobian transformation to all data points, to convert to curvilinear basis
  3. Apply curvilinear boundary conditions
  4. Apply Jacobian transformation to all data points, to convert back to Cartesian basis
  5. Compute difference between original Cartesian data and transformed data. Difference should be zero (to within roundoff) at all points except those that are influenced by outer boundary conditions.

Step 8.d: CurviBC_Playground.c: The Main C Code [Back to top]¶