{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"# Start-to-Finish Example: Applying Boundary Conditions in Curvilinear Coordinates in Three Dimensions\n",
"\n",
"## Authors: Zach Etienne & Terrence Pierre Jacques\n",
"\n",
"## 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](https://arxiv.org/abs/1712.07658).\n",
"\n",
"\n",
"\n",
"## Introduction: Applying boundary conditions in curvilinear coordinates\n",
"$$\\label{intro}$$\n",
"\n",
"We face four challenges when solving PDEs in curvilinear coordinates:\n",
"\n",
"1. [Challenge 1](#innerouterbcs): Unlike ordinary Cartesian coordinate grids, not all boundary points on uniform curvilinear coordinate grids are outer boundary points.\n",
"1. [Challenge 2](#coordinversion): Figuring out the locations to which boundary points map requires a coordinate inversion, but some coordinate systems are not easily invertible.\n",
"1. [Challenge 3](#parity): Tensors and vectors in curvilinear coordinates can change directions across \"inner\" boundaries, changing sign as a result. \n",
"1. [Challenge 4](#coordsingularities): Coordinate singularities can appear, causing a stiffening or divergence of terms in PDEs."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"### Challenge 1: Inner versus outer boundary points\n",
"$$\\label{innerouterbcs}$$\n",
"\n",
"Unlike Cartesian coordinates, the boundaries of our grid in generic curvilinear coordinates are not all outer boundaries. \n",
"\n",
"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}]$.\n",
"\n",
"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\n",
"\n",
"+ $x_0 = r \\in [0,{\\rm RMAX}]$,\n",
"+ $x_1 = \\theta \\in [0,\\pi]$, and\n",
"+ $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]$.)\n",
"\n",
"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. \n",
"\n",
"As we generally cannot apply an outer boundary condition to the inner boundaries, these boundaries will need to be treated differently.\n",
"\n",
"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*.\n",
"\n",
"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:\n",
"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)?\n",
"1. If the ghost zone point maps to an interior point, to which interior point does it map?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"### Challenge 2: Inverting coordinates\n",
"$$\\label{coordinversion}$$\n",
"\n",
"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$:\n",
"\n",
"$$\n",
"r(x_0) = {\\rm AMPL} \\left[ {\\rm const\\_dr} x_0 + \\sinh\\left(\\frac{x_0}{\\rm SINHW}\\right) / \\sinh\\left(\\frac{1}{\\rm SINHW}\\right) \\right].\n",
"$$\n",
"\n",
"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. \n",
"\n",
"*Is there an easier way of dealing with this problem than with a root finder?*"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"### Challenge 3: Parity: changes of direction in vectors and tensors across inner boundaries\n",
"$$\\label{parity}$$\n",
"\n",
"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.\n",
"\n",
"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\n",
"\\begin{align}\n",
"v^\\rho(-\\rho,\\phi,z)&=-v^\\rho(\\rho,\\phi+\\pi,z) \\\\\n",
"&= \\mathbf{e}^\\rho\\left(-\\rho,\\phi,z\\right) \\cdot \\mathbf{e}^\\rho\\left(\\rho,\\phi+\\pi,z\\right)v^\\rho(\\rho,\\phi+\\pi,z),\n",
"\\end{align}\n",
"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."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"### Challenge 4: Coordinate singularities\n",
"$$\\label{coordsingularities}$$\n",
"\n",
"Most non-Cartesian, orthogonal coordinate systems (like spherical coordinates) possess *coordinate singularities*. \n",
"\n",
"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](Tutorial-ScalarWaveCurvilinear.ipynb), the Laplacian in spherical polar coordinates takes the form\n",
"$$\n",
"\\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,\n",
"$$\n",
"which diverges at $r=0$ and $\\sin\\theta=0-$precisesly at the $\\theta=0$ and $\\theta=\\pi$ coordinate singularity. \n",
"\n",
"To avoid this divergence, we simply choose that our numerical grids be **cell-centered**. \n",
"\n",
"I.e., given the desired bounds of the grid interior to be \n",
"\n",
"\\begin{align}\n",
"x_0 &\\in [x_{0,\\ \\rm min},x_{0,\\ \\rm max}]\\\\\n",
"x_1 &\\in [x_{1,\\ \\rm min},x_{1,\\ \\rm max}]\\\\\n",
"x_2 &\\in [x_{2,\\ \\rm min},x_{2,\\ \\rm max}],\n",
"\\end{align}\n",
"\n",
"${\\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\n",
"\n",
"\\begin{align}\n",
"dx_0 &= \\frac{x_{0,\\ \\rm max} - x_{0,\\ \\rm min}}{N_0} \\\\\n",
"dx_1 &= \\frac{x_{1,\\ \\rm max} - x_{1,\\ \\rm min}}{N_1} \\\\\n",
"dx_2 &= \\frac{x_{2,\\ \\rm max} - x_{2,\\ \\rm min}}{N_2}.\n",
"\\end{align}\n",
"\n",
"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:\n",
"\n",
"\\begin{align}\n",
"x_{0,{\\rm i0}} &= x_{0,\\ \\rm min} + \\left[({\\rm i0}-{\\rm NGHOSTS}) + \\frac{1}{2}\\right] dx_0 \\\\\n",
"x_{1,{\\rm i1}} &= x_{1,\\ \\rm min} + \\left[({\\rm i1}-{\\rm NGHOSTS}) + \\frac{1}{2}\\right] dx_1 \\\\\n",
"x_{2,{\\rm i2}} &= x_{2,\\ \\rm min} + \\left[({\\rm i2}-{\\rm NGHOSTS}) + \\frac{1}{2}\\right] dx_2 \\\\\n",
"\\end{align}\n",
"where\n",
"\n",
"* ${\\rm i0}\\in [0,N_0+2\\cdot{\\rm NGHOSTS})$\n",
"* ${\\rm i1}\\in [0,N_1+2\\cdot{\\rm NGHOSTS})$\n",
"* ${\\rm i2}\\in [0,N_2+2\\cdot{\\rm NGHOSTS})$,\n",
"\n",
"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."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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\n",
"\n",
"\\begin{align}\n",
"{\\rm NGHOSTS} &= 3 \\\\\n",
"dr &= \\frac{r_{\\rm max} - 0}{N_r} \\\\\n",
"d\\theta &= \\frac{\\pi - 0}{N_\\theta} \\\\\n",
"d\\phi &= \\frac{\\pi - (-\\pi)}{N_\\phi} \\\\\n",
"r_{{\\rm i0}} &= 0 + \\left[({\\rm i0}-{\\rm NGHOSTS}) + \\frac{1}{2}\\right] dx_0 \\\\\n",
"&= \\left[({\\rm i0}-3) + \\frac{1}{2}\\right] dx_0 \\\\\n",
"\\theta_{{\\rm i1}} &= 0 + \\left[({\\rm i1}-{\\rm NGHOSTS}) + \\frac{1}{2}\\right] dx_1 \\\\\n",
"&= \\left[({\\rm i1}-3) + \\frac{1}{2}\\right] dx_1 \\\\\n",
"\\phi_{{\\rm i2}} &= -\\pi + \\left[({\\rm i2}-{\\rm NGHOSTS}) + \\frac{1}{2}\\right] dx_2 \\\\\n",
"&= -\\pi + \\left[({\\rm i2}-3) + \\frac{1}{2}\\right] dx_2, \\\\\n",
"\\end{align}\n",
"\n",
"where again\n",
"* ${\\rm i0}\\in [0,N_r+2\\cdot{\\rm NGHOSTS})$\n",
"* ${\\rm i1}\\in [0,N_\\theta+2\\cdot{\\rm NGHOSTS})$\n",
"* ${\\rm i2}\\in [0,N_\\phi+2\\cdot{\\rm NGHOSTS})$,\n",
"\n",
"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.\n",
"\n",
"Notice that in NRPy+, we use the [physics](https://en.wikipedia.org/wiki/Spherical_coordinate_system) 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`](https://en.wikipedia.org/wiki/Atan2).\n",
"\n",
"**Exercise to student**: Given the prescription above, why do the integers $N_\\theta$ and $N_\\phi$ need to be even?\n",
"\n",
"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](Tutorial-Method_of_Lines-C_Code_Generation.ipynb) timestepping algorithm used. See [Cordero-Carrión & Cerdá-Durán](https://arxiv.org/abs/1211.5930) 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.\n",
"\n",
"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).\n",
"\n",
"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\n",
"\n",
"\\begin{align}\n",
"\\bar{\\Lambda}^{\\phi} &= \\frac{\\partial \\phi}{\\partial x} \\bar{\\Lambda}^{x} + \n",
"\\frac{\\partial \\phi}{\\partial y} \\bar{\\Lambda}^{y} + \n",
"\\frac{\\partial \\phi}{\\partial z} \\bar{\\Lambda}^{z} \\\\\n",
"&= -\\frac{y}{x^2+y^2} \\bar{\\Lambda}^{x} + \n",
"\\frac{x}{x^2+y^2} \\bar{\\Lambda}^{y} \\\\\n",
"&= -\\frac{y}{(r \\sin\\theta)^2} \\bar{\\Lambda}^{x} + \n",
"\\frac{x}{(r \\sin\\theta)^2} \\bar{\\Lambda}^{y} \\\\\n",
"&= -\\frac{r \\sin\\theta \\sin\\phi}{(r \\sin\\theta)^2} \\bar{\\Lambda}^{x} + \n",
"\\frac{r \\sin\\theta \\cos\\phi}{(r \\sin\\theta)^2} \\bar{\\Lambda}^{y}\\\\\n",
"&= -\\frac{\\sin\\phi}{r \\sin\\theta} \\bar{\\Lambda}^{x} + \n",
"\\frac{\\cos\\phi}{r \\sin\\theta} \\bar{\\Lambda}^{y}\\\\\n",
"\\end{align}\n",
"\n",
"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.\n",
"\n",
"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.\n",
"\n",
"However, notice that if we define $\\lambda^{\\phi}$ such that\n",
"\n",
"$$\\bar{\\Lambda}^{\\phi} = \\frac{1}{r\\sin\\theta} \\lambda^{\\phi},$$\n",
"\n",
"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.\n",
"\n",
"Avoiding such singularities can be generalized to arbitrary coordinate systems, so long as $\\lambda^i$ is defined as:\n",
"\n",
"$$\\bar{\\Lambda}^{i} = \\frac{\\lambda^i}{\\text{scalefactor[i]}} ,$$\n",
"\n",
"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](Tutorial-BSSN_quantities.ipynb#rescaling_tensors).\n",
"\n",
"**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.**"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Table of Contents\n",
"$$\\label{toc}$$\n",
"\n",
"This notebook is organized as follows\n",
"\n",
"1. [Step 1](#basic_algorithm): Overview of boundary condition algorithm in curvilinear coordinates\n",
" 1. [Step 1.a](#challenge1): Addressing Challenge 1: Distinguishing inner from outer boundary points\n",
" 1. [Step 1.b](#challenge2): Addressing Challenge 2: Eigen-Coordinate systems\n",
" 1. [Step 1.c](#challenge3): Addressing Challenge 3: Applying parity conditions to arbitrary-rank tensors\n",
"1. [Step 2](#ccode_bc_struct): `bc_struct` and other C data structures used for storing boundary condition information\n",
"1. [Step 3](#nrpycodegen): NRPy+-based C code generation for parity conditions\n",
" 1. [Step 3.a](#dotproducts): Set up unit-vector dot products (=parity) for each of the 10 parity condition types, store to `parity_conditions_symbolic_dot_products()`\n",
" 1. [Step 3.b](#set_parity_type): Set parity type for each gridfunction, based on the digits at the end of its name, output to `dirname+gridfunction_defines.h`\n",
"1. [Step 4](#set_up__bc_gz_map_and_parity_condns): `set_up__bc_gz_map_and_parity_condns()`: C function for distinguishing inner from outer boundary points, and setting parity conditions\n",
"1. [Step 5](#set_bc_struct): `set_bcstruct()`: Using information from `set_up__bc_gz_map_and_parity_condns()` as input, set `bcstruct`\n",
"1. [Step 6](#bcstruct_c_code_driver): `driver_bcstruct()`: 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`\n",
"1. [Step 7](#extrap_bcs_curvilinear): `\"extrapolation\"` outer boundary conditions: apply quadratic polynomial extrapolation BCs\n",
"1. [Step 8](#radiation_bcs_curvilinear): `\"radiation\"` outer boundary conditions: apply `NewRad`-style BCs\n",
"1. [Step 9](#apply_bcs_curvilinear): Main driver function: `apply_bcs_curvilinear()`: quickly apply boundary and parity conditions with information from `bcstruct`\n",
"1. [Step 10](#start2finish): `CurviBC_Playground.c`: Start-to-Finish C code module for testing & validating curvilinear boundary conditions\n",
" 1. [Step 10.a](#register_gfs): Register gridfunctions of all 10 parity types in NRPy+; output gridfunction aliases to `CurviBoundaryConditions/gridfunction_defines.h`\n",
" 1. [Step 10.b](#validate): Set up test data for Curvilinear Boundary Conditions code validation\n",
" 1. [Step 10.c](#mainc): `CurviBC_Playground`'s `main.c` Code\n",
" 1. [Step 10.d](#curvibc_setupall): Add all CurviBC C codes to C function dictionary, and add CurviBC definitions to `NRPy_basic_defines.h`\n",
"1. [Step 10](#senr_compare): Validation: Compare with original SENR results\n",
"1. [Step 11](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 1: Overview of boundary condition algorithm in curvilinear coordinates \\[Back to [top](#toc)\\]\n",
"$$\\label{basic_algorithm}$$\n",
"\n",
"Here we **review** the basic algorithm for addressing Challenges [1](#innerouterbcs) [2](#coordinversion), and [3](#parity) discussed in the [Introduction](#intro) above. \n",
"\n",
"The algorithm itself is **implemented** as C code in Steps [**2**](#bc_struct) (data structures), [**3**](#set_up__bc_gz_map_and_parity_condns) (searching entire grid for inner and outer boundary points, and setting parities), [**4**](#set_bcstruct) (setting data structures for quick and efficient implementation of outer boundaries), and [**5**](#apply_bcs_curvilinear) (function to apply inner & outer boundary conditions).\n",
"\n",
"\n",
"\n",
"## Step 1.a: Addressing Challenge 1: Distinguishing inner from outer boundary points \\[Back to [top](#toc)\\]\n",
"$$\\label{challenge1}$$\n",
"\n",
"At each ghost zone grid point $\\mathbf{d}_{\\rm gz}=(x_0,x_1,x_2)$, we will do the following:\n",
"\n",
"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)$. \n",
" 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.\n",
" 1. 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.\n",
"\n",
"In detail, the algorithm is as follows:\n",
"\n",
"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\n",
" + $x(r,\\theta,\\phi) = r \\sin(\\theta) \\cos(\\phi) = x_0 \\sin(x_1) \\cos(x_2)$\n",
" + $y(r,\\theta,\\phi) = r \\sin(\\theta) \\sin(\\phi) = x_0 \\sin(x_1) \\sin(x_2)$\n",
" + $z(r,\\theta,\\phi) = r \\cos(\\theta) = x_0 \\cos(x_1)$\n",
"1. 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:\n",
" + $r_{\\rm in/OB} = x_{0, \\rm in/OB} = \\sqrt{x^2+y^2+z^2} \\in [0,\\infty)$\n",
" + $\\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]$\n",
" + $\\phi_{\\rm in/OB} = x_{2, \\rm in/OB} = {\\rm atan2}(y,x) \\in [-\\pi,\\pi]$ [Wikipedia article on atan2()](https://en.wikipedia.org/w/index.php?title=Atan2&oldid=859313982)\n",
"\n",
"1. 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:\n",
" + $\\rm{i0}_{\\rm in/OB}=\\frac{r_{\\rm in/OB} - r_{\\rm min}}{\\Delta r} - \\frac{1}{2}$\n",
" + $\\rm{i1}_{\\rm in/OB}=\\frac{\\theta_{\\rm in/OB} - \\theta_{\\rm min}}{\\Delta \\theta} - \\frac{1}{2}$\n",
" + $\\rm{i2}_{\\rm in/OB}=\\frac{\\phi_{\\rm in/OB} - \\phi_{\\rm min}}{\\Delta \\phi} - \\frac{1}{2}$\n",
"\n",
"1. 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)`.\n",
"\n",
"Following the prescription in the [SENR/NRPy+ paper](https://arxiv.org/abs/1712.07658), 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.\n",
"\n",
"\n",
"\n",
"## Step 1.b: Addressing Challenge 2: Eigen-Coordinate Systems \\[Back to [top](#toc)\\]\n",
"$$\\label{challenge2}$$\n",
"\n",
"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. \n",
"\n",
"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.\n",
"\n",
"\n",
"\n",
"## Step 1.c: Addressing Challenge 3: Applying parity conditions to arbitrary-rank tensors \\[Back to [top](#toc)\\]\n",
"$$\\label{challenge3}$$\n",
"\n",
"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.\n",
"\n",
"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$).\n",
"\n",
"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\"):\n",
"$$\n",
"P_1 = \\mathbf{e}^i\\left(\\mathbf{d}_{\\rm gz}\\right) \\cdot \\mathbf{e}^i\\left(\\mathbf{d}_{\\rm new}\\right).\n",
"$$\n",
"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}$\n",
"1. Repeat the above for the remaining components of the rank-$N$ tensor $j\\in \\{2,3,...,N\\}$, storing each $P_j$.\n",
"1. The tensor mapping from $\\mathbf{d}_{\\rm gz}$ to $\\mathbf{d}_{\\rm new}$ for this tensor $T^{ijk...}_{mnp...}$ will be given by\n",
"$$\n",
"T^{ijk...}_{lmn...}(x_0,x_1,x_2)_{\\rm gz} = \\prod_{\\ell=1}^N P_\\ell T^{ijk...}_{mnp...}(x_0,x_1,x_2)_{\\rm new}.\n",
"$$\n",
"\n",
"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 \n",
"+ 1 parity condition (the trivial +1) for scalars (rank-0 tensors)\n",
"+ 3 parity conditions for all rank-1 tensors (corresponding to each direction)\n",
"+ 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)\n",
"\n",
"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:\n",
"\n",
"The appropriate dot products determining parity condition are assigned to each gridfunction based on the following numbering:\n",
"\n",
"Tensor type | Parity type | Dot product(s) determining parity condition (see equation above)\n",
"--- | --- | ---\n",
"Scalar (Rank-0 tensor) | 0 | (*none*)\n",
"Rank-1 tensor in **i0** direction | 1 | $\\mathbf{e}^0\\left(\\mathbf{d}_{\\rm gz}\\right) \\cdot \\mathbf{e}^0\\left(\\mathbf{d}_{\\rm new}\\right)$\n",
"Rank-1 tensor in **i1** direction | 2 | $\\mathbf{e}^1\\left(\\mathbf{d}_{\\rm gz}\\right) \\cdot \\mathbf{e}^1\\left(\\mathbf{d}_{\\rm new}\\right)$\n",
"Rank-1 tensor in **i2** direction | 3 | $\\mathbf{e}^2\\left(\\mathbf{d}_{\\rm gz}\\right) \\cdot \\mathbf{e}^2\\left(\\mathbf{d}_{\\rm new}\\right)$\n",
"Rank-2 tensor in **i0-i0** direction | 4 | $\\left[\\mathbf{e}^0\\left(\\mathbf{d}_{\\rm gz}\\right) \\cdot \\mathbf{e}^0\\left(\\mathbf{d}_{\\rm new}\\right)\\right]^2 = 1$\n",
"Rank-2 tensor in **i0-i1** direction | 5 | $\\left[\\mathbf{e}^0\\left(\\mathbf{d}_{\\rm gz}\\right) \\cdot \\mathbf{e}^0\\left(\\mathbf{d}_{\\rm new}\\right)\\right]\\left[\\mathbf{e}^1\\left(\\mathbf{d}_{\\rm gz}\\right) \\cdot \\mathbf{e}^1\\left(\\mathbf{d}_{\\rm new}\\right)\\right]$\n",
"Rank-2 tensor in **i0-i2** direction | 6 | $\\left[\\mathbf{e}^0\\left(\\mathbf{d}_{\\rm gz}\\right) \\cdot \\mathbf{e}^0\\left(\\mathbf{d}_{\\rm new}\\right)\\right]\\left[\\mathbf{e}^2\\left(\\mathbf{d}_{\\rm gz}\\right) \\cdot \\mathbf{e}^2\\left(\\mathbf{d}_{\\rm new}\\right)\\right]$\n",
"Rank-2 tensor in **i1-i1** direction | 7 | $\\left[\\mathbf{e}^1\\left(\\mathbf{d}_{\\rm gz}\\right) \\cdot \\mathbf{e}^1\\left(\\mathbf{d}_{\\rm new}\\right)\\right]^2 = 1$\n",
"Rank-2 tensor in **i1-i2** direction | 8 | $\\left[\\mathbf{e}^1\\left(\\mathbf{d}_{\\rm gz}\\right) \\cdot \\mathbf{e}^1\\left(\\mathbf{d}_{\\rm new}\\right)\\right]\\left[\\mathbf{e}^2\\left(\\mathbf{d}_{\\rm gz}\\right) \\cdot \\mathbf{e}^2\\left(\\mathbf{d}_{\\rm new}\\right)\\right]$\n",
"Rank-2 tensor in **i2-i2** direction | 9 | $\\left[\\mathbf{e}^2\\left(\\mathbf{d}_{\\rm gz}\\right) \\cdot \\mathbf{e}^2\\left(\\mathbf{d}_{\\rm new}\\right)\\right]^2 = 1$\n",
"\n",
"\n",
"In the following few steps, we will document the data structures for and implementation of this boundary condition algorithm."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 2: `bc_struct`: Data structure for storing boundary condition information \\[Back to [top](#toc)\\]\n",
"$$\\label{ccode_bc_struct}$$\n",
"\n",
"Here we define `bc_struct`, a C data structure that stores all information needed to apply boundary conditions at all boundary points. \n",
"\n",
"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:\n",
"\n",
"```C\n",
"typedef struct __bcstruct__ {\n",
" outer_bc **outer; // Array of 1D arrays, of length \n",
" // [NGHOSTS][num_ob_gz_pts[which_outer_ghostzone_point]]\n",
"\n",
" inner_bc **inner; // Array of 1D arrays, of length \n",
" // [NGHOSTS][num_ib_gz_pts[which_inner_ghostzone_point]]\n",
"\n",
" // Arrays storing number of outer/inner boundary ghostzone points at each ghostzone,\n",
" // of length NGHOSTS:\n",
" int *num_ob_gz_pts; \n",
" int *num_ib_gz_pts;\n",
"} bc_struct;\n",
"\n",
"```\n",
"\n",
"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\". \n",
"\n",
"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.\n",
"\n",
"The same reasoning holds when considering ghost zones that are not outer boundary points. \n",
"\n",
"\n",
"Next we summarize all basic data structures that appear within `bc_struct`. \n",
"\n",
"* `inner_bc`:\n",
"```c\n",
"typedef struct __inner_bc__ {\n",
" gz_map inner_bc_dest_pt;\n",
" gz_map inner_bc_src_pt;\n",
" int8_t parity[10]; // We store the 10 parity conditions in 10 int8_t integers, \n",
" // one for each condition. Note that these conditions can \n",
" // only take one of two values: +1 or -1, hence the use of \n",
" // int8_t, the smallest C data type.\n",
"} inner_bc;\n",
"```\n",
" * `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.\n",
" * `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.\n",
" * `parity[10]` Parity information ($\\pm 1$) at the given inner ghost zone, for all 10 gridfunction parity types.\n",
"* `outer_bc`: \n",
"```c \n",
"typedef struct __outer_bc__ {\n",
" gz_map outer_bc_dest_pt;\n",
" int8_t FACEi0,FACEi1,FACEi2; // FACEi* takes values of -1, 0, and +1 only,\n",
" // corresponding to MAXFACE, NUL, and MINFACE\n",
" // respectively.\n",
" // Thus int8_t (one byte each, the smallest C \n",
" // type) is sufficient.\n",
"} outer_bc;\n",
"```\n",
" * 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 \n",
" * `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.\n",
"* `num_ib/ob_gz_pts[which_gz]`: The number of inner/outer boundary points on ghost zone layer `which_gz`"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"execution": {
"iopub.execute_input": "2021-09-08T15:47:14.965489Z",
"iopub.status.busy": "2021-09-08T15:47:14.964975Z",
"iopub.status.idle": "2021-09-08T15:47:15.183581Z",
"shell.execute_reply": "2021-09-08T15:47:15.183117Z"
}
},
"outputs": [],
"source": [
"# Step P1: Import needed NRPy+ core modules:\n",
"from outputC import outputC # NRPy+: Core C code output module\n",
"from outputC import outC_NRPy_basic_defines_h_dict # NRPy+: Core C code output module\n",
"from outputC import add_to_Cfunction_dict # NRPy+: Core C code output module\n",
"from outputC import Cfunction # NRPy+: Core C code output module\n",
"from outputC import indent_Ccode # NRPy+: Core C code output module\n",
"import NRPy_param_funcs as par # NRPy+: Parameter interface\n",
"import sympy as sp # SymPy: The Python computer algebra package upon which NRPy+ depends\n",
"import grid as gri # NRPy+: Functions having to do with numerical grids\n",
"import indexedexp as ixp # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support\n",
"import reference_metric as rfm # NRPy+: Reference metric support\n",
"import cmdline_helper as cmd # NRPy+: Multi-platform Python command-line interface\n",
"import finite_difference as fin # NRPy+: Finite-difference module\n",
"import shutil, os, sys # Standard Python modules for multiplatform OS-level functions\n",
"from UnitTesting.assert_equal import check_zero # NRPy+: Checks whether an expression evaluates to zero.\n",
"\n",
"# Step P2: Create C code output directory:\n",
"Ccodesrootdir = os.path.join(\"CurviBoundaryConditions_Ccodes_new_way/\")\n",
"# P2.a: First remove C code output directory if it exists\n",
"# Courtesy https://stackoverflow.com/questions/303200/how-do-i-remove-delete-a-folder-that-is-not-empty\n",
"shutil.rmtree(Ccodesrootdir, ignore_errors=True)\n",
"# P2.b: Then create a fresh directory\n",
"# Then create a fresh directory\n",
"cmd.mkdir(Ccodesrootdir)\n",
"\n",
"# CoordSystem = \"Cylindrical\"\n",
"CoordSystem = \"Spherical\"\n",
"outer_bcs_type = \"extrapolation\"\n",
"# outer_bcs_type = \"radiation\"\n",
"outer_bcs_FDORDER = 4 # Set to -1 to adopt the same value as finite_difference::FD_CENTDERIVS_ORDER\n",
"\n",
"# Set the finite differencing order to 4; although this module doesn't compute\n",
"# finite difference derivatives, it does need to set the number of ghost zone\n",
"# cells, which is generally based on NGHOSTS, which itself depends\n",
"# on finite_difference::FD_CENTDERIVS_ORDER.\n",
"par.set_parval_from_str(\"finite_difference::FD_CENTDERIVS_ORDER\", 4)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Define/declare core data structures for curvilinear BCs within `NRPy_basic_defines.h`:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"execution": {
"iopub.execute_input": "2021-09-08T15:47:15.187904Z",
"iopub.status.busy": "2021-09-08T15:47:15.187482Z",
"iopub.status.idle": "2021-09-08T15:47:15.189166Z",
"shell.execute_reply": "2021-09-08T15:47:15.188788Z"
}
},
"outputs": [],
"source": [
"# First register basic C data structures/macros inside NRPy_basic_defines.h\n",
"def NRPy_basic_defines_CurviBC_data_structures():\n",
" return r\"\"\"\n",
"// NRPy+ Curvilinear Boundary Conditions: Core data structures\n",
"// Documented in: Tutorial-Start_to_Finish-Curvilinear_BCs.ipynb\n",
"typedef struct __ghostzone_map__ {\n",
" short i0,i1,i2; // i0,i1,i2 stores values from -1 (used to indicate outer boundary)\n",
" // to Nxx_plus_2NGHOSTS*. We assume that grid extents beyond the\n",
" // limits of short (i.e., beyond about 32,000) are unlikely. This\n",
" // can be easily extended if needed, though.\n",
"} gz_map;\n",
"\n",
"typedef struct __parity__ {\n",
" int8_t parity[10]; // We store the 10 parity conditions in 10 int8_t integers,\n",
" // one for each condition. Note that these conditions can\n",
" // only take one of two values: +1 or -1, hence the use of\n",
" // int8_t, the smallest C data type.\n",
"} parity_condition;\n",
"\n",
"typedef struct __inner_bc__ {\n",
" gz_map inner_bc_dest_pt;\n",
" gz_map inner_bc_src_pt;\n",
" int8_t parity[10]; // We store the 10 parity conditions in 10 int8_t integers,\n",
" // one for each condition. Note that these conditions can\n",
" // only take one of two values: +1 or -1, hence the use of\n",
" // int8_t, the smallest C data type.\n",
"} inner_bc;\n",
"\n",
"typedef struct __outer_bc__ {\n",
" gz_map outer_bc_dest_pt;\n",
" int8_t FACEi0,FACEi1,FACEi2; // FACEi* takes values of -1, 0, and +1 only,\n",
" // corresponding to MAXFACE, NUL, and MINFACE\n",
" // respectively.\n",
" // Thus int8_t (one byte each, the smallest C\n",
" // type) is sufficient.\n",
"} outer_bc;\n",
"\n",
"typedef struct __bcstruct__ {\n",
" outer_bc **outer; // Array of 1D arrays, of length\n",
" // [NGHOSTS][num_ob_gz_pts[which_outer_ghostzone_point]]\n",
"\n",
" inner_bc **inner; // Array of 1D arrays, of length\n",
" // [NGHOSTS][num_ib_gz_pts[which_inner_ghostzone_point]]\n",
"\n",
" // Arrays storing number of outer/inner boundary ghostzone points at each ghostzone,\n",
" // of length NGHOSTS:\n",
" int *num_ob_gz_pts;\n",
" int *num_ib_gz_pts;\n",
"} bc_struct;\n",
"\"\"\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 3: NRPy+-based C code generation for parity conditions \\[Back to [top](#toc)\\]\n",
"$$\\label{nrpycodegen}$$\n",
"\n",
"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.\n",
"\n",
"However, there are two parts to the overall algorithm that must be generated by NRPy+, namely\n",
"\n",
"1. [Step 3.a](#dotproducts): `parity_conditions_symbolic_dot_products()`: Based on the chosen reference metric, sets up the needed unit-vector dot products for each of the 10 parity condition types.\n",
"1. [Step 3.b](#set_parity_type): 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`\n",
"\n",
"\n",
"\n",
"## 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()` \\[Back to [top](#toc)\\]\n",
"$$\\label{dotproducts}$$\n",
"\n",
"Next we generate the C code necessary to perform needed dot products for filling in the parity condition arrays inside `bcstruct`. \n",
"\n",
"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](#challenge2). We paste the table here again, for quick reference:\n",
"\n",
"The appropriate symbolic dot products determining parity condition are assigned to each gridfunction based on the following numbering:\n",
"\n",
"Tensor type | Parity type | Dot product(s) determining parity condition (\n",
"--- | --- | ---\n",
"Scalar (Rank-0 tensor) | 0 | (*none*)\n",
"Rank-1 tensor in **i0** direction | 1 | $\\mathbf{e}^0\\left(\\mathbf{d}_{\\rm gz}\\right) \\cdot \\mathbf{e}^0\\left(\\mathbf{d}_{\\rm new}\\right)$\n",
"Rank-1 tensor in **i1** direction | 2 | $\\mathbf{e}^1\\left(\\mathbf{d}_{\\rm gz}\\right) \\cdot \\mathbf{e}^1\\left(\\mathbf{d}_{\\rm new}\\right)$\n",
"Rank-1 tensor in **i2** direction | 3 | $\\mathbf{e}^2\\left(\\mathbf{d}_{\\rm gz}\\right) \\cdot \\mathbf{e}^2\\left(\\mathbf{d}_{\\rm new}\\right)$\n",
"Rank-2 tensor in **i0-i0** direction | 4 | $\\left[\\mathbf{e}^0\\left(\\mathbf{d}_{\\rm gz}\\right) \\cdot \\mathbf{e}^0\\left(\\mathbf{d}_{\\rm new}\\right)\\right]^2 = 1$\n",
"Rank-2 tensor in **i0-i1** direction | 5 | $\\left[\\mathbf{e}^0\\left(\\mathbf{d}_{\\rm gz}\\right) \\cdot \\mathbf{e}^0\\left(\\mathbf{d}_{\\rm new}\\right)\\right]\\left[\\mathbf{e}^1\\left(\\mathbf{d}_{\\rm gz}\\right) \\cdot \\mathbf{e}^1\\left(\\mathbf{d}_{\\rm new}\\right)\\right]$\n",
"Rank-2 tensor in **i0-i2** direction | 6 | $\\left[\\mathbf{e}^0\\left(\\mathbf{d}_{\\rm gz}\\right) \\cdot \\mathbf{e}^0\\left(\\mathbf{d}_{\\rm new}\\right)\\right]\\left[\\mathbf{e}^2\\left(\\mathbf{d}_{\\rm gz}\\right) \\cdot \\mathbf{e}^2\\left(\\mathbf{d}_{\\rm new}\\right)\\right]$\n",
"Rank-2 tensor in **i1-i1** direction | 7 | $\\left[\\mathbf{e}^1\\left(\\mathbf{d}_{\\rm gz}\\right) \\cdot \\mathbf{e}^1\\left(\\mathbf{d}_{\\rm new}\\right)\\right]^2 = 1$\n",
"Rank-2 tensor in **i1-i2** direction | 8 | $\\left[\\mathbf{e}^1\\left(\\mathbf{d}_{\\rm gz}\\right) \\cdot \\mathbf{e}^1\\left(\\mathbf{d}_{\\rm new}\\right)\\right]\\left[\\mathbf{e}^2\\left(\\mathbf{d}_{\\rm gz}\\right) \\cdot \\mathbf{e}^2\\left(\\mathbf{d}_{\\rm new}\\right)\\right]$\n",
"Rank-2 tensor in **i2-i2** direction | 9 | $\\left[\\mathbf{e}^2\\left(\\mathbf{d}_{\\rm gz}\\right) \\cdot \\mathbf{e}^2\\left(\\mathbf{d}_{\\rm new}\\right)\\right]^2 = 1$\n",
"\n",
"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: \n",
"\n",
"```\n",
"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);\n",
"```\n",
"\n",
"To wit (as described above), there are 10 parity types for BSSN evolved variables, which include tensors up to and including rank-2."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"execution": {
"iopub.execute_input": "2021-09-08T15:47:15.196318Z",
"iopub.status.busy": "2021-09-08T15:47:15.195711Z",
"iopub.status.idle": "2021-09-08T15:47:15.198183Z",
"shell.execute_reply": "2021-09-08T15:47:15.197651Z"
}
},
"outputs": [],
"source": [
"# Set unit-vector dot products (=parity) for each of the 10 parity condition types\n",
"def parity_conditions_symbolic_dot_products():\n",
" parity = ixp.zerorank1(DIM=10)\n",
" UnitVectors_inner = ixp.zerorank2()\n",
" xx0_inbounds,xx1_inbounds,xx2_inbounds = sp.symbols(\"xx0_inbounds xx1_inbounds xx2_inbounds\", real=True)\n",
" for i in range(3):\n",
" for j in range(3):\n",
" UnitVectors_inner[i][j] = rfm.UnitVectors[i][j].subs(rfm.xx[0],xx0_inbounds).subs(rfm.xx[1],xx1_inbounds).subs(rfm.xx[2],xx2_inbounds)\n",
" # Type 0: scalar\n",
" parity[0] = sp.sympify(1)\n",
" # Type 1: i0-direction vector or one-form\n",
" # Type 2: i1-direction vector or one-form\n",
" # Type 3: i2-direction vector or one-form\n",
" for i in range(3):\n",
" for Type in range(1,4):\n",
" parity[Type] += rfm.UnitVectors[Type-1][i]*UnitVectors_inner[Type-1][i]\n",
" # Type 4: i0i0-direction rank-2 tensor\n",
" # parity[4] = parity[1]*parity[1]\n",
" # Type 5: i0i1-direction rank-2 tensor\n",
" # Type 6: i0i2-direction rank-2 tensor\n",
" # Type 7: i1i1-direction rank-2 tensor\n",
" # Type 8: i1i2-direction rank-2 tensor\n",
" # Type 9: i2i2-direction rank-2 tensor\n",
" count = 4\n",
" for i in range(3):\n",
" for j in range(i,3):\n",
" parity[count] = parity[i+1]*parity[j+1]\n",
" count = count + 1\n",
"\n",
" lhs_strings = []\n",
" for i in range(10):\n",
" lhs_strings.append(\"REAL_parity_array[\"+str(i)+\"]\")\n",
" outstr = \"\"\"\n",
" // NRPy+ Curvilinear Boundary Conditions: Unit vector dot products for all\n",
" // ten parity conditions, in given coordinate system.\n",
" // Needed for automatically determining sign of tensor across coordinate boundary.\n",
" // Documented in: Tutorial-Start_to_Finish-Curvilinear_BCs.ipynb\n",
"\"\"\"\n",
" return outstr + outputC(parity,lhs_strings, filename=\"returnstring\", params=\"preindent=4\")\n",
"\n",
"# print(\"\\n\\nExample: parity type 1's dot product is given by: \\n\"+lhs_strings[1]+\" = \"+str(parity[1]))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"## 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](#toc)\\]\n",
"$$\\label{set_parity_type}$$\n",
"\n",
"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 rather robust algorithm, because `gri.register_gridfunctions()` in [grid.py](../edit/grid.py) will throw an error if a gridfunction's **base** name ends in an integer. This strict syntax was added with the express purpose of making it easier to set parity types based solely on the gridfunction name.\n",
"\n",
"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 `NRPy_basic_defines.h`."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"execution": {
"iopub.execute_input": "2021-09-08T15:47:15.214520Z",
"iopub.status.busy": "2021-09-08T15:47:15.211805Z",
"iopub.status.idle": "2021-09-08T15:47:15.215728Z",
"shell.execute_reply": "2021-09-08T15:47:15.216134Z"
},
"scrolled": false
},
"outputs": [],
"source": [
"def NRPy_basic_defines_set_gridfunction_defines_with_parity_types(verbose=True):\n",
" # First add human-readable gridfunction aliases (grid.py) to NRPy_basic_defines dictionary,\n",
" evolved_variables_list, auxiliary_variables_list, auxevol_variables_list = gri.gridfunction_lists()\n",
"\n",
" # Step 3.b: set the parity conditions on all gridfunctions in gf_list,\n",
" # based on how many digits are at the end of their names\n",
" def set_parity_types(list_of_gf_names):\n",
" parity_type = []\n",
" for name in list_of_gf_names:\n",
" for gf in gri.glb_gridfcs_list:\n",
" if gf.name == name:\n",
" parity_type__orig_len = len(parity_type)\n",
" if gf.DIM < 3 or gf.DIM > 4:\n",
" print(\"Error: Cannot currently specify parity conditions on gridfunctions with DIM<3 or >4.\")\n",
" sys.exit(1)\n",
" if gf.rank == 0:\n",
" parity_type.append(0)\n",
" elif gf.rank == 1:\n",
" if gf.DIM == 3:\n",
" parity_type.append(int(gf.name[-1]) + 1) # = 1 for e.g., beta^0; = 2 for e.g., beta^1, etc.\n",
" elif gf.DIM == 4:\n",
" parity_type.append(int(gf.name[-1])) # = 0 for e.g., b4^0; = 1 for e.g., beta^1, etc.\n",
" elif gf.rank == 2:\n",
" if gf.DIM == 3:\n",
" # element of a list; a[-2] the\n",
" # second-to-last element, etc.\n",
" idx0 = gf.name[-2]\n",
" idx1 = gf.name[-1]\n",
" if idx0 == \"0\" and idx1 == \"0\":\n",
" parity_type.append(4)\n",
" elif (idx0 == \"0\" and idx1 == \"1\") or (idx0 == \"1\" and idx1 == \"0\"):\n",
" parity_type.append(5)\n",
" elif (idx0 == \"0\" and idx1 == \"2\") or (idx0 == \"2\" and idx1 == \"0\"):\n",
" parity_type.append(6)\n",
" elif idx0 == \"1\" and idx1 == \"1\":\n",
" parity_type.append(7)\n",
" elif (idx0 == \"1\" and idx1 == \"2\") or (idx0 == \"2\" and idx1 == \"1\"):\n",
" parity_type.append(8)\n",
" elif idx0 == \"2\" and idx1 == \"2\":\n",
" parity_type.append(9)\n",
" elif gf.DIM == 4:\n",
" idx0 = gf.name[-2]\n",
" idx1 = gf.name[-1]\n",
" # g4DD00 = g_{tt} : parity type = 0\n",
" # g4DD01 = g_{tx} : parity type = 1\n",
" # g4DD02 = g_{ty} : parity type = 2\n",
" # g4DD0a = g_{ta} : parity type = a\n",
" if idx0 == \"0\":\n",
" parity_type.append(int(idx1))\n",
" elif idx1 == \"0\":\n",
" parity_type.append(int(idx0))\n",
" if idx0 == \"1\" and idx1 == \"1\":\n",
" parity_type.append(4)\n",
" elif (idx0 == \"1\" and idx1 == \"2\") or (idx0 == \"2\" and idx1 == \"1\"):\n",
" parity_type.append(5)\n",
" elif (idx0 == \"1\" and idx1 == \"3\") or (idx0 == \"3\" and idx1 == \"1\"):\n",
" parity_type.append(6)\n",
" elif idx0 == \"2\" and idx1 == \"2\":\n",
" parity_type.append(7)\n",
" elif (idx0 == \"2\" and idx1 == \"3\") or (idx0 == \"3\" and idx1 == \"2\"):\n",
" parity_type.append(8)\n",
" elif idx0 == \"3\" and idx1 == \"3\":\n",
" parity_type.append(9)\n",
" if len(parity_type) == parity_type__orig_len:\n",
" print(\"Error: Could not figure out parity type for \"+gf.gftype+\" gridfunction: \" + gf.name,gf.DIM,gf.name[-2],gf.name[-1],gf.rank)\n",
" sys.exit(1)\n",
" if len(parity_type) != len(list_of_gf_names):\n",
" print(\"Error: For some reason the length of the parity types list did not match the length of the gf list.\")\n",
" sys.exit(1)\n",
" return parity_type\n",
"\n",
" evol_parity_type = set_parity_types(evolved_variables_list)\n",
" aux_parity_type = set_parity_types(auxiliary_variables_list)\n",
" auxevol_parity_type = set_parity_types(auxevol_variables_list)\n",
"\n",
" # Output all gridfunctions to Ccodesrootdir/gridfunction_defines.h\n",
" # ... then append to the file the parity type for each gridfunction.\n",
" outstr = \"\"\"\n",
"/* PARITY TYPES FOR ALL GRIDFUNCTIONS.\n",
" * SEE \\\"Tutorial-Start_to_Finish-Curvilinear_BCs.ipynb\\\" FOR DEFINITIONS. */\n",
"\"\"\"\n",
" if len(evolved_variables_list) > 0:\n",
" outstr += \"static const int8_t evol_gf_parity[\" + str(len(evolved_variables_list)) + \"] = { \"\n",
" for i in range(len(evolved_variables_list) - 1):\n",
" outstr += str(evol_parity_type[i]) + \", \"\n",
" outstr += str(evol_parity_type[len(evolved_variables_list) - 1]) + \" };\\n\"\n",
"\n",
" if len(auxiliary_variables_list) > 0:\n",
" outstr += \"static const int8_t aux_gf_parity[\" + str(len(auxiliary_variables_list)) + \"] = { \"\n",
" for i in range(len(auxiliary_variables_list) - 1):\n",
" outstr += str(aux_parity_type[i]) + \", \"\n",
" outstr += str(aux_parity_type[len(auxiliary_variables_list) - 1]) + \" };\\n\"\n",
"\n",
" if len(auxevol_variables_list) > 0:\n",
" outstr += \"static const int8_t auxevol_gf_parity[\" + str(len(auxevol_variables_list)) + \"] = { \"\n",
" for i in range(len(auxevol_variables_list) - 1):\n",
" outstr += str(auxevol_parity_type[i]) + \", \"\n",
" outstr += str(auxevol_parity_type[len(auxevol_variables_list) - 1]) + \" };\\n\"\n",
"\n",
" if verbose == True:\n",
" for i in range(len(evolved_variables_list)):\n",
" print(\"Evolved gridfunction \\\"\" + evolved_variables_list[i] + \"\\\" has parity type \" + str(\n",
" evol_parity_type[i]) + \".\")\n",
" for i in range(len(auxiliary_variables_list)):\n",
" print(\"Auxiliary gridfunction \\\"\" + auxiliary_variables_list[i] + \"\\\" has parity type \" + str(\n",
" aux_parity_type[i]) + \".\")\n",
" for i in range(len(auxevol_variables_list)):\n",
" print(\"AuxEvol gridfunction \\\"\" + auxevol_variables_list[i] + \"\\\" has parity type \" + str(\n",
" auxevol_parity_type[i]) + \".\")\n",
" return outstr"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# 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](#toc)\\]\n",
"$$\\label{set_up__bc_gz_map_and_parity_condns}$$\n",
"\n",
"This step is performed using only the [Eigen-Coordinate](#challenge2) corresponding to the chosen `reference_metric::CoordSystem`.\n",
"\n",
"First we generate the C code needed for applying boundary conditions in generic coordinate systems, using the [Eigen-Coordinate](#challenge2) approach described above:\n",
"1. $\\left(x(x_0,x_1,x_2),y(x_0,x_1,x_2),z(x_0,x_1,x_2)\\right)$, (`EigenCoord_xx_to_Cart()`)\n",
"1. $\\left(x_0(x,y,z),x_1(x,y,z),x_2(x,y,z)\\right)$, (`EigenCoord_Cart_to_xx()`):"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"execution": {
"iopub.execute_input": "2021-09-08T15:47:15.220481Z",
"iopub.status.busy": "2021-09-08T15:47:15.219875Z",
"iopub.status.idle": "2021-09-08T15:47:15.221957Z",
"shell.execute_reply": "2021-09-08T15:47:15.221580Z"
}
},
"outputs": [],
"source": [
"def EigenCoord_xx_to_Cart(i012suffix = \"\"):\n",
" # Step 1: Find the Eigen-Coordinate and set up the Eigen-Coordinate's reference metric:\n",
" CoordSystem_orig = par.parval_from_str(\"reference_metric::CoordSystem\")\n",
" par.set_parval_from_str(\"reference_metric::CoordSystem\",rfm.get_EigenCoord())\n",
" rfm.reference_metric()\n",
"\n",
" # Step 2: Output C code for the Eigen-Coordinate mapping from xx->Cartesian:\n",
" outstr = \"\"\" {\n",
" // xx_to_Cart for EigenCoordinate \"\"\"+rfm.get_EigenCoord()+r\"\"\" (orig coord = \"\"\"+CoordSystem_orig+r\"\"\"):\n",
" REAL xx0 = xx[0][i0\"\"\"+i012suffix+\"\"\"];\n",
" REAL xx1 = xx[1][i1\"\"\"+i012suffix+\"\"\"];\n",
" REAL xx2 = xx[2][i2\"\"\"+i012suffix+\"\"\"];\\n\"\"\"+ \\\n",
" outputC([rfm.xx_to_Cart[0],rfm.xx_to_Cart[1],rfm.xx_to_Cart[2]],\n",
" [\"xCart[0]\",\"xCart[1]\",\"xCart[2]\"],\n",
" \"returnstring\", params=\"preindent=3\")+\" }\\n\"\n",
"\n",
" # Step 3: Restore reference_metric::CoordSystem back to the original CoordSystem\n",
" par.set_parval_from_str(\"reference_metric::CoordSystem\",CoordSystem_orig)\n",
" rfm.reference_metric()\n",
"\n",
" # Step 4: Return EigenCoord xx_to_Cart C code\n",
" return outstr"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"execution": {
"iopub.execute_input": "2021-09-08T15:47:15.227418Z",
"iopub.status.busy": "2021-09-08T15:47:15.227017Z",
"iopub.status.idle": "2021-09-08T15:47:15.228957Z",
"shell.execute_reply": "2021-09-08T15:47:15.228562Z"
}
},
"outputs": [],
"source": [
"def EigenCoord_Cart_to_xx():\n",
" # Step 1: Find the Eigen-Coordinate and set up the Eigen-Coordinate's reference metric:\n",
" CoordSystem_orig = par.parval_from_str(\"reference_metric::CoordSystem\")\n",
" par.set_parval_from_str(\"reference_metric::CoordSystem\",rfm.get_EigenCoord())\n",
" rfm.reference_metric()\n",
"\n",
" # Step 2: Output the Eigen-Coordinate mapping from Cartesian->xx:\n",
" # Step 2.a: Sanity check: First make sure that rfm.Cart_to_xx has been set. Error out if not!\n",
" if rfm.Cart_to_xx[0] == 0 or rfm.Cart_to_xx[1] == 0 or rfm.Cart_to_xx[2] == 0:\n",
" print(\"ERROR: rfm.Cart_to_xx[], which maps Cartesian -> xx, has not been set for\")\n",
" print(\" reference_metric::CoordSystem = \"+par.parval_from_str(\"reference_metric::CoordSystem\"))\n",
" print(\" Boundary conditions in curvilinear coordinates REQUiRE this be set.\")\n",
" sys.exit(1)\n",
" # Step 2.b: Output C code for the Eigen-Coordinate mapping from Cartesian->xx:\n",
" outstr = \"\"\" // Cart_to_xx for EigenCoordinate \"\"\"+rfm.get_EigenCoord()+r\"\"\" (orig coord = \"\"\"+CoordSystem_orig+\");\\n\"\n",
" outstr += outputC([rfm.Cart_to_xx[0],rfm.Cart_to_xx[1],rfm.Cart_to_xx[2]],\n",
" [\"Cart_to_xx0_inbounds\",\"Cart_to_xx1_inbounds\",\"Cart_to_xx2_inbounds\"],\n",
" filename=\"returnstring\", params=\"preindent=2\")\n",
"\n",
" # Step 3: Restore reference_metric::CoordSystem back to the original CoordSystem\n",
" par.set_parval_from_str(\"reference_metric::CoordSystem\",CoordSystem_orig)\n",
" rfm.reference_metric()\n",
"\n",
" # Step 4: Return EigenCoord Cart_to_xx C code\n",
" return outstr"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next we set up `set_up__bc_gz_map_and_parity_condns()`, which 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](#basic_algorithm). \n",
"\n",
"That is to say, at each coordinate point $(x_0,x_1,x_2)$ situated at grid index `(i0,i1,i2)`:\n",
"\n",
"1. Convert the curvilinear coordinate $(x_0,x_1,x_2)$ to the corresponding Cartesian coordinate $(x,y,z)$\n",
"1. 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)$.\n",
"1. 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`.\n",
"1. 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](#challenge2). The C code for computing the needed symbolic dot products is generated above in [Step 2](#dotproducts)."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"execution": {
"iopub.execute_input": "2021-09-08T15:47:15.236052Z",
"iopub.status.busy": "2021-09-08T15:47:15.235636Z",
"iopub.status.idle": "2021-09-08T15:47:15.237789Z",
"shell.execute_reply": "2021-09-08T15:47:15.237434Z"
}
},
"outputs": [],
"source": [
"def add_to_Cfunction_dict_set_up__bc_gz_map_and_parity_condns(rel_path_to_Cparams=os.path.join(\".\")):\n",
" includes = [os.path.join(rel_path_to_Cparams, \"NRPy_basic_defines.h\"),\n",
" os.path.join(rel_path_to_Cparams, \"NRPy_function_prototypes.h\")]\n",
" desc = \"\"\n",
" c_type = \"void\"\n",
" name = \"set_up__bc_gz_map_and_parity_condns\"\n",
" params = \"\"\"const paramstruct *restrict params,\n",
" REAL *restrict xx[3], gz_map *restrict bc_gz_map,parity_condition *restrict bc_parity_conditions\"\"\"\n",
" body = r\"\"\"\n",
" // xx[0][j] = xxmin[0] + ((REAL)(j-NGHOSTS) + (1.0/2.0))*dxx0;\n",
" // -> xxmin[0] = xx[0][0] - ((REAL)(0-NGHOSTS) + (1.0/2.0))*dxx0\n",
" const REAL xxmin[3] = { xx[0][0] - ((REAL)(0-NGHOSTS) + (1.0/2.0))*dxx0,\n",
" xx[1][0] - ((REAL)(0-NGHOSTS) + (1.0/2.0))*dxx1,\n",
" xx[2][0] - ((REAL)(0-NGHOSTS) + (1.0/2.0))*dxx2 };\n",
" //fprintf(stderr,\"hey inside setbc: %e %e %e | %e %e\\n\",xxmin[0],xxmin[1],xxmin[2],xx[0][0],dxx0);\n",
" LOOP_REGION(0,Nxx_plus_2NGHOSTS0,0,Nxx_plus_2NGHOSTS1,0,Nxx_plus_2NGHOSTS2) {\n",
" // Step 1: Convert the (curvilinear) coordinate (x0,x1,x2) to Cartesian coordinates\n",
" REAL xCart[3];\n",
"\"\"\" + EigenCoord_xx_to_Cart(i012suffix=\"\") + r\"\"\"\n",
" //EigenCoord_xx_to_Cart(params, xx, i0,i1,i2, xCart);\n",
" REAL Cartx = xCart[0];\n",
" REAL Carty = xCart[1];\n",
" REAL Cartz = xCart[2];\n",
"\n",
" // Step 2: Find the (i0_inbounds,i1_inbounds,i2_inbounds) corresponding to the above Cartesian coordinate.\n",
" // If (i0_inbounds,i1_inbounds,i2_inbounds) is in a ghost zone, then it must equal (i0,i1,i2), and\n",
" // the point is an outer boundary point.\n",
" // Otherwise (i0_inbounds,i1_inbounds,i2_inbounds) is in the grid interior, and data at (i0,i1,i2)\n",
" // must be replaced with data at (i0_inbounds,i1_inbounds,i2_inbounds), but multiplied by the\n",
" // appropriate parity condition (+/- 1).\n",
" REAL Cart_to_xx0_inbounds,Cart_to_xx1_inbounds,Cart_to_xx2_inbounds;\n",
"\"\"\" + EigenCoord_Cart_to_xx() + r\"\"\"\n",
" int i0_inbounds = (int)( (Cart_to_xx0_inbounds - xxmin[0] - (1.0/2.0)*dxx0 + ((REAL)NGHOSTS)*dxx0)/dxx0 + 0.5 );\n",
" int i1_inbounds = (int)( (Cart_to_xx1_inbounds - xxmin[1] - (1.0/2.0)*dxx1 + ((REAL)NGHOSTS)*dxx1)/dxx1 + 0.5 );\n",
" int i2_inbounds = (int)( (Cart_to_xx2_inbounds - xxmin[2] - (1.0/2.0)*dxx2 + ((REAL)NGHOSTS)*dxx2)/dxx2 + 0.5 );\n",
"\n",
" // Step 2.a: (Sanity/validation check) Convert the interior point\n",
" // x0(i0_inbounds),x1(i1_inbounds),x2(i2_inbounds) to Cartesian coordinates,\n",
" // make sure that the Cartesian coordinate matches the Cartesian coordinate of\n",
" // x0(i0),x1(i1),x2(i2). If not, error out!\n",
" REAL xCart_orig[3]; for(int ii=0;ii<3;ii++) xCart_orig[ii] = xCart[ii];\n",
" //EigenCoord_xx_to_Cart(params, xx, i0_inbounds,i1_inbounds,i2_inbounds, xCart);\n",
"\"\"\" + EigenCoord_xx_to_Cart(i012suffix=\"_inbounds\") + r\"\"\"\n",
"\n",
"//fprintf(stderr,\"Cartesian agreement: ( %.15e %.15e %.15e ) ?= ( %.15e %.15e %.15e )\\n\",\n",
"// (double)xCart_orig[0],(double)xCart_orig[1],(double)xCart_orig[2],\n",
"// (double)xCart[0],(double)xCart[1],(double)xCart[2]);\n",
"\n",
"#define EPS_ABS 1e-8\n",
" if(fabs( (double)(xCart_orig[0] - xCart[0]) ) > EPS_ABS ||\n",
" fabs( (double)(xCart_orig[1] - xCart[1]) ) > EPS_ABS ||\n",
" fabs( (double)(xCart_orig[2] - xCart[2]) ) > EPS_ABS) {\n",
" fprintf(stderr,\"Error. Cartesian disagreement: ( %.15e %.15e %.15e ) != ( %.15e %.15e %.15e )\\n\",\n",
" (double)xCart_orig[0],(double)xCart_orig[1],(double)xCart_orig[2],\n",
" (double)xCart[0],(double)xCart[1],(double)xCart[2]);\n",
" exit(1);\n",
" }\n",
"\n",
" // Step 3: Set bc_gz_map and bc_parity_conditions.\n",
" if(i0_inbounds-i0 == 0 && i1_inbounds-i1 == 0 && i2_inbounds-i2 == 0) {\n",
" // Step 3.a: Iff we are on an outer boundary point or in the grid\n",
" // interior, i0_inbounds==i0, i1_inbounds==i1, and\n",
" // i2_inbounds==i2, and inner boundary conditions do not\n",
" // apply: set bc_gz_map to -1, and parity=1.\n",
" bc_gz_map[IDX3S(i0,i1,i2)].i0=-1;\n",
" bc_gz_map[IDX3S(i0,i1,i2)].i1=-1;\n",
" bc_gz_map[IDX3S(i0,i1,i2)].i2=-1;\n",
" for(int which_parity=0; which_parity<10; which_parity++) {\n",
" bc_parity_conditions[IDX3S(i0,i1,i2)].parity[which_parity] = 1;\n",
" }\n",
" } else {\n",
" // Step 3.b: If we are on an *inner* boundary point:\n",
" // 1. Set bc_gz_map at (i0,i1,i2) to the point\n",
" // in the interior to which this boundary\n",
" // point maps, and\n",
" // 2. Perform the unit vector dot products\n",
" // necessary to set all 10 possible parity\n",
" // conditions, calling function\n",
" // set_parity_from_unit_vector_dot_product()\n",
" bc_gz_map[IDX3S(i0,i1,i2)].i0=i0_inbounds;\n",
" bc_gz_map[IDX3S(i0,i1,i2)].i1=i1_inbounds;\n",
" bc_gz_map[IDX3S(i0,i1,i2)].i2=i2_inbounds;\n",
" const REAL xx0 = xx[0][i0];\n",
" const REAL xx1 = xx[1][i1];\n",
" const REAL xx2 = xx[2][i2];\n",
" const REAL xx0_inbounds = xx[0][i0_inbounds];\n",
" const REAL xx1_inbounds = xx[1][i1_inbounds];\n",
" const REAL xx2_inbounds = xx[2][i2_inbounds];\n",
" REAL REAL_parity_array[10];\n",
" {\n",
" // Evaluate dot products needed for setting parity\n",
" // conditions at a given point (xx0,xx1,xx2),\n",
" // using C code generated by NRPy+\n",
"\"\"\"+parity_conditions_symbolic_dot_products()+r\"\"\"\n",
" }\n",
" //eval_symbolic_dot_products_to_set_parity_conditions(params, REAL_parity_array, xx0,xx1,xx2,\n",
" // xx0_inbounds,xx1_inbounds,xx2_inbounds);\n",
" for(int whichparity=0;whichparity<10;whichparity++) {\n",
" //printf(\"Good? Parity %d evaluated to %e\\n\",whichparity,(double)REAL_parity_array[whichparity]);\n",
" // Perform sanity check on parity array output: should be +1 or -1 to within 8 significant digits:\n",
" if( (REAL_parity_array[whichparity] > 0 && fabs(REAL_parity_array[whichparity] - (+1)) > 1e-8) ||\n",
" (REAL_parity_array[whichparity] <= 0 && fabs(REAL_parity_array[whichparity] - (-1)) > 1e-8) ) {\n",
" fprintf(stderr,\"Error at point (%d %d %d); (%e %e %e); maps to (%e %e %e).\\n\",\n",
" i0,i1,i2, xx0,xx1,xx2, xx0_inbounds,xx1_inbounds,xx2_inbounds);\n",
" fprintf(stderr,\"Parity evaluated to %e , which is not within 8 significant digits of +1 or -1.\\n\",\n",
" REAL_parity_array[whichparity]);\n",
" exit(1);\n",
" }\n",
" if(REAL_parity_array[whichparity] < 0.0) bc_parity_conditions[IDX3S(i0,i1,i2)].parity[whichparity] = -1;\n",
" if(REAL_parity_array[whichparity] > 0.0) bc_parity_conditions[IDX3S(i0,i1,i2)].parity[whichparity] = +1;\n",
" } // END for(int whichparity=0;whichparity<10;whichparity++)\n",
" } // END if(i0_inbounds-i0 == 0 && i1_inbounds-i1 == 0 && i2_inbounds-i2 == 0)\n",
" } // END LOOP_REGION(0,Nxx_plus_2NGHOSTS0,0,Nxx_plus_2NGHOSTS1,0,Nxx_plus_2NGHOSTS2)\n",
"\"\"\"\n",
" add_to_Cfunction_dict(\n",
" includes=includes,\n",
" desc=desc,\n",
" c_type=c_type, name=name, params=params,\n",
" body=body,\n",
" rel_path_to_Cparams=rel_path_to_Cparams)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 5: `set_bcstruct()`: Using information from `set_up__bc_gz_map_and_parity_condns()` as input, set `bc_struct` \\[Back to [top](#toc)\\]\n",
"$$\\label{set_bc_struct}$$\n",
"\n",
"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`.\n",
"\n",
"The algorithm follows the `bc_struct` data type defined above in [Step 2](#bc_struct), 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. \n",
"\n",
"1. Count the number of outer and inner boundary points, store to `num_ob_pts` and `num_ib_pts`, respectively.\n",
"1. 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.\n",
" 1. At all outer boundary ghost zones, allocate memory a single member of the `outer_bc` data type.\n",
" 1. At all inner boundary ghost zones, allocate memory a single member of the `inner_bc` data type.\n",
"1. 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.\n",
"1. Store information needed for outer boundary conditions, to `outer_bc_dest_pt` and `outer_bc_face` arrays.\n",
"1. Store information needed for inner boundary conditions, including interior point to which inner ghost zone maps, and parity conditions for all 10 gridfunction types."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"execution": {
"iopub.execute_input": "2021-09-08T15:47:15.245220Z",
"iopub.status.busy": "2021-09-08T15:47:15.244606Z",
"iopub.status.idle": "2021-09-08T15:47:15.246704Z",
"shell.execute_reply": "2021-09-08T15:47:15.246202Z"
}
},
"outputs": [],
"source": [
"def add_to_Cfunction_dict_set_bcstruct():\n",
" includes = [\"NRPy_basic_defines.h\", \"NRPy_function_prototypes.h\"]\n",
"\n",
" prefunc = \"\"\"\n",
"static const int8_t MAXFACE = -1;\n",
"static const int8_t NUL = +0;\n",
"static const int8_t MINFACE = +1;\n",
"\"\"\"\n",
" desc = \"\"\"\n",
"set_bcstruct() loops from the innermost boundary\n",
" ghostzones on the cube (\"which_gz==0\",\n",
" corresponding to the single layer of ghostzones\n",
" closest to the interior data), and at each\n",
" ghostzone layer, we apply the following 5-step\n",
" algorithm:\n",
"Step 1: Count the number of outer and inner\n",
" boundary points, store to\n",
" num_ob_pts and num_ib_pts, respectively.\n",
"Step 2: Now that we know the number of outer\n",
" boundary points on this ghostzone layer,\n",
" allocate memory needed for storing the\n",
" outer and inner boundary condition data.\n",
"Step 2.a: At all outer boundary ghost zones, allocate\n",
" memory for a single member of the outer_bc\n",
" data type.\n",
"Step 2.b: At all inner boundary ghost zones, allocate\n",
" memory for a single member of the inner_bc\n",
" data type.\n",
"Step 3: Store the number of outer and inner boundary\n",
" points on each ghostzone layer, where e.g.,\n",
" which_gz==0 corresponds to the innermost\n",
" ghostzones on the numerical domain.\n",
"Step 4: Store information needed for outer boundary\n",
" conditions, to outer_bc_dest_pt and\n",
" outer_bc_face arrays.\n",
"Step 5: Store information needed for inner boundary\n",
" conditions, including interior point to which\n",
" inner ghost zone maps, and parity conditions\n",
" for all 10 gridfunction types.\n",
"\"\"\"\n",
" c_type = \"void\"\n",
"\n",
" name = \"set_bcstruct\"\n",
" params = \"\"\"const paramstruct *restrict params,\n",
" gz_map *restrict bc_gz_map,\n",
" parity_condition *bc_parity_conditions,\n",
" bc_struct *restrict bcstruct\"\"\"\n",
" body = r\"\"\"\n",
" int imin[3] = { NGHOSTS, NGHOSTS, NGHOSTS };\n",
" int imax[3] = { Nxx_plus_2NGHOSTS0-NGHOSTS, Nxx_plus_2NGHOSTS1-NGHOSTS, Nxx_plus_2NGHOSTS2-NGHOSTS };\n",
"\n",
" // Loop from the innermost ghostzone on the cube (which_gz==0) and work outward.\n",
" // This ordering is necessary, as ghostzones at which_gz==1 will generally\n",
" // depend on ghostzones at which_gz==0 being already set.\n",
" for(int which_gz = 0; which_gz < NGHOSTS; which_gz++) {\n",
"\n",
" // Step 1: Count the number of outer and inner\n",
" // boundary points, store to\n",
" // num_ob_pts and num_ib_pts, respectively.\n",
"\"\"\"\n",
" body += r\"\"\"\n",
"#define COUNT_INNER_OR_OUTER if(bc_gz_map[IDX3S(i0,i1,i2)].i0==-1) { num_ob_pts++;} else { num_ib_pts++; }\n",
"\"\"\"\n",
" body += r\"\"\"\n",
" int num_ob_pts = 0;\n",
" int num_ib_pts = 0;\n",
" LOOP_REGION(imin[0]-1,imin[0], imin[1],imax[1], imin[2],imax[2]) { COUNT_INNER_OR_OUTER } imin[0]--;\n",
" LOOP_REGION(imax[0],imax[0]+1, imin[1],imax[1], imin[2],imax[2]) { COUNT_INNER_OR_OUTER } imax[0]++;\n",
" LOOP_REGION(imin[0],imax[0], imin[1]-1,imin[1], imin[2],imax[2]) { COUNT_INNER_OR_OUTER } imin[1]--;\n",
" LOOP_REGION(imin[0],imax[0], imax[1],imax[1]+1, imin[2],imax[2]) { COUNT_INNER_OR_OUTER } imax[1]++;\n",
" LOOP_REGION(imin[0],imax[0], imin[1],imax[1], imin[2]-1,imin[2]) { COUNT_INNER_OR_OUTER } imin[2]--;\n",
" LOOP_REGION(imin[0],imax[0], imin[1],imax[1], imax[2],imax[2]+1) { COUNT_INNER_OR_OUTER } imax[2]++;\n",
"\n",
" // Step 2: Now that we know the number of outer boundary points on this ghostzone\n",
" // layer, we allocate memory needed for storing the outer and inner boundary\n",
" // condition data.\n",
"\n",
" // Step 2.a: At all outer boundary ghost zones, allocate memory for a single member of the outer_bc\n",
" // data type.\n",
" bcstruct->outer[which_gz] = (outer_bc *)malloc(sizeof(outer_bc)*num_ob_pts);\n",
" // Step 2.b: At all inner boundary ghost zones, allocate memory for a single member of the inner_bc\n",
" // data type.\n",
" bcstruct->inner[which_gz] = (inner_bc *)malloc(sizeof(inner_bc)*num_ib_pts);\n",
"\n",
" // Step 3: Store the number of outer and inner boundary points on each ghostzone layer, where e.g.,\n",
" // which_gz==0 corresponds to the innermost ghostzones on the numerical domain.\n",
" bcstruct->num_ob_gz_pts[which_gz] = num_ob_pts;\n",
" bcstruct->num_ib_gz_pts[which_gz] = num_ib_pts;\n",
"\n",
" // Reset imin[] and imax[], to prepare for the next step.\n",
" for(int ii=0;ii<3;ii++) {imin[ii]++; imax[ii]--;}\n",
"\n",
" // Step 4: Store information needed for outer boundary conditions, to outer_bc_dest_pt[which_gz][]\n",
"// and outer_bc_face[which_gz][] arrays:\n",
"\"\"\"\n",
" body += \"\"\"#define OB_SET(facei0,facei1,facei2) if(bc_gz_map[IDX3S(i0,i1,i2)].i0==-1) { \\\\\"\"\"\n",
" body += r\"\"\"\n",
" bcstruct->outer[which_gz][pt].outer_bc_dest_pt.i0 = i0; \\\n",
" bcstruct->outer[which_gz][pt].outer_bc_dest_pt.i1 = i1; \\\n",
" bcstruct->outer[which_gz][pt].outer_bc_dest_pt.i2 = i2; \\\n",
" bcstruct->outer[which_gz][pt].FACEi0= facei0; \\\n",
" bcstruct->outer[which_gz][pt].FACEi1= facei1; \\\n",
" bcstruct->outer[which_gz][pt].FACEi2= facei2; \\\n",
" pt++; }\n",
"\n",
" int pt = 0;\n",
" LOOP_REGION(imin[0]-1,imin[0], imin[1],imax[1], imin[2],imax[2]) {OB_SET(MINFACE,NUL,NUL)} imin[0]--;\n",
" LOOP_REGION(imax[0],imax[0]+1, imin[1],imax[1], imin[2],imax[2]) {OB_SET(MAXFACE,NUL,NUL)} imax[0]++;\n",
" LOOP_REGION(imin[0],imax[0], imin[1]-1,imin[1], imin[2],imax[2]) {OB_SET(NUL,MINFACE,NUL)} imin[1]--;\n",
" LOOP_REGION(imin[0],imax[0], imax[1],imax[1]+1, imin[2],imax[2]) {OB_SET(NUL,MAXFACE,NUL)} imax[1]++;\n",
" LOOP_REGION(imin[0],imax[0], imin[1],imax[1], imin[2]-1,imin[2]) {OB_SET(NUL,NUL,MINFACE)} imin[2]--;\n",
" LOOP_REGION(imin[0],imax[0], imin[1],imax[1], imax[2],imax[2]+1) {OB_SET(NUL,NUL,MAXFACE)} imax[2]++;\n",
" // fprintf(stderr,\"num OB points with which_gz = %d: %d | should be: %d\\n\",which_gz,pt,num_ob_gz_pts[which_gz]);\n",
"\n",
" // Reset imin[] and imax[], to prepare for the next step.\n",
" for(int ii=0;ii<3;ii++) {imin[ii]++; imax[ii]--;}\n",
"\n",
" // Step 5: Store information needed for inner boundary conditions, including interior point to which\n",
" // inner ghost zone maps, and parity conditions for all 10 gridfunction types.\n",
"#define IB_SET if(bc_gz_map[IDX3S(i0,i1,i2)].i0!=-1) { \\\n",
" bcstruct->inner[which_gz][pt].inner_bc_dest_pt.i0=i0; \\\n",
" bcstruct->inner[which_gz][pt].inner_bc_dest_pt.i1=i1; \\\n",
" bcstruct->inner[which_gz][pt].inner_bc_dest_pt.i2=i2; \\\n",
" bcstruct->inner[which_gz][pt].inner_bc_src_pt.i0 =bc_gz_map[IDX3S(i0,i1,i2)].i0; \\\n",
" bcstruct->inner[which_gz][pt].inner_bc_src_pt.i1 =bc_gz_map[IDX3S(i0,i1,i2)].i1; \\\n",
" bcstruct->inner[which_gz][pt].inner_bc_src_pt.i2 =bc_gz_map[IDX3S(i0,i1,i2)].i2; \\\n",
" for(int ii=0;ii<10;ii++) { \\\n",
" bcstruct->inner[which_gz][pt].parity[ii] = \\\n",
" (int8_t)bc_parity_conditions[IDX3S(i0,i1,i2)].parity[ii]; } \\\n",
" pt++; }\n",
"\n",
" pt = 0;\n",
" LOOP_REGION(imin[0]-1,imin[0], imin[1],imax[1], imin[2],imax[2]) {IB_SET} imin[0]--;\n",
" LOOP_REGION(imax[0],imax[0]+1, imin[1],imax[1], imin[2],imax[2]) {IB_SET} imax[0]++;\n",
" LOOP_REGION(imin[0],imax[0], imin[1]-1,imin[1], imin[2],imax[2]) {IB_SET} imin[1]--;\n",
" LOOP_REGION(imin[0],imax[0], imax[1],imax[1]+1, imin[2],imax[2]) {IB_SET} imax[1]++;\n",
" LOOP_REGION(imin[0],imax[0], imin[1],imax[1], imin[2]-1,imin[2]) {IB_SET} imin[2]--;\n",
" LOOP_REGION(imin[0],imax[0], imin[1],imax[1], imax[2],imax[2]+1) {IB_SET} imax[2]++;\n",
"\n",
" } // END for(int which_gz = 0; which_gz < NGHOSTS; which_gz++)\n",
"\"\"\"\n",
" add_to_Cfunction_dict(\n",
" includes=includes, prefunc=prefunc,\n",
" desc=desc,\n",
" c_type=c_type, name=name, params=params,\n",
" body=indent_Ccode(body),\n",
" rel_path_to_Cparams=os.path.join(\".\"))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 6: `driver_bcstruct()`: 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](#toc)\\]\n",
"$$\\label{bcstruct_c_code_driver}$$\n",
"\n",
"`driver_bcstruct()` 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.\n",
"\n",
"For completeness, we also output `bcstruct_freemem()`, which frees memory allocated within `bcstruct`."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"execution": {
"iopub.execute_input": "2021-09-08T15:47:15.251820Z",
"iopub.status.busy": "2021-09-08T15:47:15.251421Z",
"iopub.status.idle": "2021-09-08T15:47:15.253335Z",
"shell.execute_reply": "2021-09-08T15:47:15.252944Z"
}
},
"outputs": [],
"source": [
"# driver_bcstruct() allocates memory for boundary condition pointers\n",
"# and calls functions referenced above, including\n",
"# set_up__bc_gz_map_and_parity_condns() and set_bcstruct(). It then\n",
"# frees unneeded memory after bcstruct has been fully initialized.\n",
"def add_to_Cfunction_dict_driver_bcstruct():\n",
" includes = [\"NRPy_basic_defines.h\", \"NRPy_function_prototypes.h\"]\n",
" desc = \"\"\"driver_bcstruct(): Set up bcstruct.\n",
"WARNING: Needs Nxx_plus_2NGHOSTS_tot * [sizeof(gz_map) + sizeof(parity_condition) ] bytes of temporary memory!\n",
" (This memory will be deallocated at end of function, in freemem_bcstruct().)\n",
" Note that gz_map consists of 3*sizeof(short = 2 bytes), and parity_condition consists of 10*sizeof(int8_t = 1 byte)\n",
" Thus the total cost is about 16 bytes at all gridpoints, or 2 double-precision gridfunctions.\n",
" To avoid memory overload, be sure to set this up prior to allocation of other gridfunctions.\n",
"STEPS:\n",
"1. Allocate memory for bc_gz_map, which maps inner gzs to the\n",
" appropriate interior gridpoint. In the case of outer gzs, the\n",
" gz point maps to itself.\n",
"2. Allocate storage for parity_condition, which\n",
"\"\"\"\n",
" c_type = \"void\"\n",
" name = \"driver_bcstruct\"\n",
" params = \"\"\"const paramstruct *restrict params, bc_struct *restrict bcstruct, REAL *restrict xx[3]\"\"\"\n",
" set_bcstruct_extra_arg = \"\"\n",
" body = \"\"\"\n",
"const int Nxx_plus_2NGHOSTS_tot = Nxx_plus_2NGHOSTS0 * Nxx_plus_2NGHOSTS1 * Nxx_plus_2NGHOSTS2;\n",
"\n",
"// Step 1: Allocate memory storage for bc_gz_map, which\n",
"// in the case a boundary point is a *parity*\n",
"// boundary, is set to the interior, non-\n",
"// boundary point corresponding to the same\n",
"// Cartesian gridpoint. Otherwise bc_gz_map\n",
"// is set to (i0,i1,i2) = (-1,-1,-1).\n",
"gz_map *restrict bc_gz_map = (gz_map *restrict)malloc(sizeof(gz_map)*Nxx_plus_2NGHOSTS_tot);\n",
"\n",
"// Step 2: Allocate memory storage for bc_parity_conditions,\n",
"// which store parity conditions for all 10\n",
"// gridfunction types at all grid points.\n",
"parity_condition *restrict bc_parity_conditions = (parity_condition *restrict)malloc(sizeof(parity_condition)*Nxx_plus_2NGHOSTS_tot);\n",
"\n",
"// Step 3: Set bc_gz_map and bc_parity_conditions at *all*\n",
"// points; on the boundary and otherwise.\n",
"set_up__bc_gz_map_and_parity_condns(params, xx, bc_gz_map,\n",
" bc_parity_conditions);\n",
"\n",
"// Step 4: Declare and allocate memory for bcstruct,\n",
"// which will store all information needed for\n",
"// applying the boundary conditions.\n",
"bcstruct->outer = (outer_bc **)malloc(sizeof(outer_bc *)*NGHOSTS);\n",
"bcstruct->inner = (inner_bc **)malloc(sizeof(inner_bc *)*NGHOSTS);\n",
"bcstruct->num_ob_gz_pts = ( int *)malloc(sizeof(int)*NGHOSTS);\n",
"bcstruct->num_ib_gz_pts = ( int *)malloc(sizeof(int)*NGHOSTS);\n",
"\n",
"// Step 5: Store all information needed to quickly and\n",
"// efficiently apply boundary conditions. This\n",
"// function transfers all information from\n",
"// bc_gz_map (defined at *all gridpoints*) into\n",
"// bcstruct (defined only at boundary points).\n",
"// Thus when this function has finished,\n",
"// bc_gz_map is no longer needed.\n",
"set_bcstruct(params,bc_gz_map,\n",
" bc_parity_conditions,\n",
" bcstruct\"\"\"+set_bcstruct_extra_arg+r\"\"\");\n",
"// Do not apply outer boundary conditions on grids that do not possess the outer boundary!\n",
"//if(params->has_outer_boundary == 0) {\n",
"// for(int i=0;inum_ob_gz_pts[i] = 0;\n",
"//}\n",
"\n",
"// Step 6: As described in Step 4, bc_gz_map is no\n",
"// longer needed at this point, so we free its\n",
"// memory. Farewell, friend!\n",
"free(bc_gz_map);\n",
"free(bc_parity_conditions);\n",
"\"\"\"\n",
" add_to_Cfunction_dict(\n",
" includes=includes,\n",
" desc=desc,\n",
" c_type=c_type, name=name, params=params,\n",
" body=indent_Ccode(body, \" \"),\n",
" rel_path_to_Cparams=os.path.join(\".\"))"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"execution": {
"iopub.execute_input": "2021-09-08T15:47:15.257125Z",
"iopub.status.busy": "2021-09-08T15:47:15.256577Z",
"iopub.status.idle": "2021-09-08T15:47:15.259040Z",
"shell.execute_reply": "2021-09-08T15:47:15.258473Z"
}
},
"outputs": [],
"source": [
"def add_to_Cfunction_dict_freemem_bcstruct():\n",
" includes = [\"NRPy_basic_defines.h\", \"NRPy_function_prototypes.h\"]\n",
" desc = \"Free memory allocated within bcstruct\"\n",
" c_type = \"void\"\n",
" name = \"freemem_bcstruct\"\n",
" params = \"const paramstruct *restrict params, const bc_struct *restrict bcstruct\"\n",
" body = r\"\"\" for(int i=0;iouter[i]); free(bcstruct->inner[i]); }\n",
" free(bcstruct->outer); free(bcstruct->inner);\n",
" free(bcstruct->num_ob_gz_pts); free(bcstruct->num_ib_gz_pts);\n",
"\"\"\"\n",
" add_to_Cfunction_dict(\n",
" includes=includes,\n",
" desc=desc,\n",
" c_type=c_type, name=name, params=params,\n",
" body=body,\n",
" rel_path_to_Cparams=os.path.join(\".\"))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 7: `\"extrapolation\"` outer boundary conditions: apply quadratic polynomial extrapolation \\[Back to [top](#toc)\\]\n",
"$$\\label{extrap_bcs_curvilinear}$$\n",
"\n",
"As an option, quadratic extrapolation may be applied to each outer boundary point `(i0,i1,i2)`, as follows.\n",
"\n",
"Suppose the outer boundary point is at the `i0=max(i0)` face. Then we fit *known* data at `i0-3`, `i0-2`, and `i0-1` [i.e., $f_{-3}=f(x_{i0-3}=x_{-3})$, $f_{-2}=f(x_{i0-2}=x_{-2})$, and $f_{-1}=f(x_{i0-1}=x_{-1})$] to the unique quadratic polynomial:\n",
"\n",
"\\begin{align}\n",
"f_{-3} &= c_2 x_{-3}^2 + c_1 x_{-3} + c_0 \\\\\n",
"f_{-2} &= c_2 x_{-2}^2 + c_1 x_{-2} + c_0 \\\\\n",
"f_{-1} &= c_2 x_{-1}^2 + c_1 x_{-1} + c_0 \\\\\n",
"\\end{align}\n",
"\n",
"We wish to extrapolate to $f_0=f(x_0)$. Since our grid has uniform spacing, \n",
"\n",
"* $x_{-3}=x_0-3\\Delta x$, \n",
"* $x_{-2}=x_0-2\\Delta x$, and\n",
"* $x_{-1}=x_0-\\Delta x$.\n",
"\n",
"The extrapolated value $f_0$ cannot depend on the choice of the fiducial $x_0$ (i.e., it will hold for *any* choice of $x_0$), so without loss of generality we will set $x_0=0$:\n",
"\n",
"$$\n",
"\\mathbf{A c} =\n",
"\\left[\n",
"\\begin{array}{ccc}\n",
" 1 & x_{-3} & x_{-3}^2 \\\\\n",
" 1 & x_{-2} & x_{-2}^2 \\\\\n",
" 1 & x_{-1} & x_{-1}^2 \\\\\n",
"\\end{array}\n",
"\\right]\n",
"\\left[\n",
"\\begin{array}{c}\n",
" c_0 \\\\\n",
" c_1 \\\\\n",
" c_2 \\\\\n",
"\\end{array}\n",
"\\right]\n",
"=\n",
"\\left[\n",
"\\begin{array}{ccc}\n",
" 1 & -3 \\Delta x & 9 \\Delta x^2 \\\\\n",
" 1 & -2 \\Delta x & 4 \\Delta x^2 \\\\\n",
" 1 & -1 \\Delta x & \\Delta x^2 \\\\\n",
"\\end{array}\n",
"\\right]\n",
"\\left[\n",
"\\begin{array}{c}\n",
" c_0 \\\\\n",
" c_1 \\\\\n",
" c_2 \\\\\n",
"\\end{array}\n",
"\\right]\n",
"=\n",
"\\left[\n",
"\\begin{array}{ccc}\n",
" 1 & -3 & 9 \\\\\n",
" 1 & -2 & 4 \\\\\n",
" 1 & -1 & 1 \\\\\n",
"\\end{array}\n",
"\\right]\n",
"\\left[\n",
"\\begin{array}{c}\n",
" c_0 \\\\\n",
" c_1 \\Delta x \\\\\n",
" c_2 \\Delta x^2 \\\\\n",
"\\end{array}\n",
"\\right]\n",
"=\n",
"\\left[\n",
"\\begin{array}{c}\n",
" f_{-3} \\\\\n",
" f_{-2} \\\\\n",
" f_{-1} \\\\\n",
"\\end{array}\n",
"\\right]\n",
"= \\mathbf{f}.\n",
"$$\n",
"\n",
"This is known as the [Vandermonde matrix](https://en.wikipedia.org/wiki/Vandermonde_matrix) for the quadratic polynomial, and its solution for $c_0$, $c_1$, and $c_2$ will be unique. But before we invert the matrix, we note that as we wish to solve for $f(x_0)=f(0)$, our quadratic polynomial simplifies to:\n",
"\n",
"$$\n",
"f_{0} = c_2 x_{0}^2 + c_1 x_{0} + c_0 = c_0.\n",
"$$\n",
"\n",
"Thus we need only extract the value of $c_0$."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"⎡ 1 -3 3 ⎤\n",
"⎢ ⎥\n",
"⎢3/2 -4 5/2⎥\n",
"⎢ ⎥\n",
"⎣1/2 -1 1/2⎦\n"
]
}
],
"source": [
"from sympy import symbols,Matrix,factor,pretty_print,simplify,latex\n",
"MM = Matrix([[1,-3,9],\n",
" [1,-2,4],\n",
" [1,-1,1]])\n",
"# print(latex(factor(MM.inv())))\n",
"pretty_print(factor(MM.inv()))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Thus we get \n",
"\n",
"$$\n",
"f_0 = c_0 = f_{-3} - 3 f_{-2} + 3 f_{-1},\n",
"$$\n",
"To determine the coefficient at the `i0=min(i0)` face, the above analysis can be repeated:"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"⎡ 1 -3 3 ⎤\n",
"⎢ ⎥\n",
"⎢-3/2 4 -5/2⎥\n",
"⎢ ⎥\n",
"⎣1/2 -1 1/2 ⎦\n"
]
}
],
"source": [
"from sympy import symbols,Matrix,factor,pretty_print,simplify,latex\n",
"MM = Matrix([[1,+3,9],\n",
" [1,+2,4],\n",
" [1,+1,1]])\n",
"# print(latex(factor(MM.inv())))\n",
"pretty_print(factor(MM.inv()))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"but the $c_0$ coefficient is basically the same, just replace $f_{-3}\\to f_{+3}$, etc:\n",
"\n",
"$$\n",
"f_0 = c_0 = f_{+3} - 3 f_{+2} + 3 f_{+1}.\n",
"$$\n",
"\n",
"The resulting extrapolation algorithm appears in the `BC_UPDATE_OUTER()` macro below. I.e.,\n",
"\n",
"```c\n",
"#define EXTRAP_BC_UPDATE_OUTER(which_gf, i0,i1,i2, FACEX0,FACEX1,FACEX2) { \\\n",
" const int idx3 = IDX3S(i0,i1,i2); \\\n",
" gfs[IDX4S(which_gf,i0,i1,i2)] = \\\n",
" +3.0*gfs[IDX4S(which_gf,i0+1*FACEX0,i1+1*FACEX1,i2+1*FACEX2)] \\\n",
" -3.0*gfs[IDX4S(which_gf,i0+2*FACEX0,i1+2*FACEX1,i2+2*FACEX2)] \\\n",
" +1.0*gfs[IDX4S(which_gf,i0+3*FACEX0,i1+3*FACEX1,i2+3*FACEX2)]; \\\n",
" }\n",
"```\n",
"\n",
"This function is meant to be called within a [Method of Lines timestepping algorithm](Tutorial-Method_of_Lines-C_Code_Generation.ipynb), at or near the end of each substep."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 8: `\"radiation\"` outer boundary conditions: apply `NewRad`-style BCs \\[Back to [top](#toc)\\]\n",
"$$\\label{radiation_bcs_curvilinear}$$\n",
"\n",
"Here we describe outgoing radiation (a.k.a., Sommerfeld) boundary conditions.\n",
"\n",
"Our implementation follows that of the [Einstein Toolkit](https://einsteintoolkit.org/)'s [`NewRad` thorn](https://einsteintoolkit.org/thornguide/EinsteinEvolve/NewRad/documentation.html), but extends to arbitrary finite-difference order.\n",
"\n",
"The basic *ansatz* is that at the outer boundary, an arbitrary field $f$ behaves as an outgoing spherical wave:\n",
"\n",
"$$\n",
"f(r,t) = f_{r\\to\\infty} + \\frac{w(r - ct)}{r} + \\frac{K}{r^2},\n",
"$$\n",
"where\n",
"\n",
"* $f_{r\\to\\infty}$ is the (constant) value the field $f$ tends to as $r\\to\\infty$.\n",
"* $w(r - ct)/r$ is a solution to the spherically symmetric wave equation ($\\partial_t^2 (r w) - c^2 \\partial_r^2 (r w) = 0$) for an *outgoing* spherical wave.\n",
"* $K/r^2$ corrects for the next-leading-order radial falloff, where $K$ is a constant that is determined by analyzing the behavior of $f(r)$ just inside the outer boundary of the numerical domain.\n",
"\n",
"This boundary condition approach thus contains two free parameters: $c$ (the wavespeed for field $f$) and $f_{r\\to\\infty}$. Both of these should be known for any given field prior to attempting a solution.\n",
"\n",
"For convenience we apply this boundary condition not to $f$ itself but to $\\partial_t f$, which will always be computed in a [Method of Lines](Tutorial-Method_of_Lines-C_Code_Generation.ipynb) approach within `NRPy+`.\n",
"\n",
"Taking the first time derivative we get:\n",
"\n",
"$$\n",
"\\partial_t f = -c \\frac{w'(r - ct)}{r}.\n",
"$$\n",
"\n",
"As $w$ represents an outgoing wave, in which temporal and spatial derivatives are directly related, we will find it convenient to also compute the radial derivative $\\partial_r f$ as well:\n",
"\n",
"$$\n",
"\\partial_r f = \\frac{w'(r - ct)}{r} - \\frac{w(r - ct)}{r^2} - 2\\frac{K}{r^3}.\n",
"$$\n",
"\n",
"Combining these two equations we get:\n",
"\n",
"$$\n",
"\\partial_t f = -c \\left(\\partial_r f + \\frac{w(r - ct)}{r^2} + 2\\frac{K}{r^3}\\right).\n",
"$$\n",
"\n",
"Next we use the *ansatz* to compute $w(r - ct)/r^2$:\n",
"\n",
"\\begin{align}\n",
"f(r,t) &= f_{r\\to\\infty} + \\frac{w(r - ct)}{r} + \\frac{K}{r^2} \\\\\n",
"\\implies \\frac{w(r - ct)}{r^2} &= \\frac{f - f_{r\\to\\infty}}{r} - \\frac{K}{r^3}.\n",
"\\end{align}\n",
"\n",
"This enables us to rewrite $\\partial_t f$ as\n",
"\n",
"\\begin{align}\n",
"\\partial_t f &= -c \\left(\\partial_r f + \\frac{f - f_{r\\to\\infty}}{r} - \\frac{K}{r^3} + 2\\frac{K}{r^3}\\right) \\\\\n",
"&= -c \\left(\\partial_r f + \\frac{f - f_{r\\to\\infty}}{r} + \\frac{K}{r^3}\\right) \\\\\n",
"&= -\\frac{c}{r} \\left[r \\partial_r f + \\left(f - f_{r\\to\\infty}\\right)\\right] + \\frac{k}{r^3},\n",
"\\end{align}\n",
"where $k=-Kc$ just re-expresses the constant."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Numerical Implementation\n",
"\n",
"We start with the equation derived in the previous section:\n",
"$$\n",
"\\partial_t f = -\\frac{c}{r} \\left[r \\partial_r f + \\left(f - f_{r\\to\\infty}\\right)\\right] + \\frac{k}{r^3}.\n",
"$$\n",
"\n",
"First note that the right-hand-side of this equation is applied *at the current time* (i.e., the same time at which $\\partial_t f$ is evaluated in the Method of Lines timestepping). I.e., it is applied prior to the MoL update. Further at the current time, $f$ is known at *all* gridpoints (interior and boundary).\n",
"\n",
"### $\\partial_r f$ term:\n",
"\n",
"As this boundary condition must be applicable to *any* curvilinear coordinate system (and not just spherical coordinates), $\\partial_r f$ must be computed via\n",
"\n",
"$$\n",
"\\partial_r f = \\frac{\\partial x^i}{\\partial r} \\partial_{i} f,\n",
"$$\n",
"\n",
"where $\\partial_i f$ is evaluated using finite-difference derivatives, suitably upwinded to avoid extending beyond the domain of the grid.\n",
"\n",
"The term $\\partial x^i/\\partial r$ cannot be immediately computed, as $x^i$ is never written as a function of $r$. However, in all curvilinear coordinates $r$, $\\theta$, and $\\phi$ are constructed from $x^i$. Thus we have exact expressions for the Jacobian:\n",
"\n",
"$$\n",
"J_i^j = \\frac{\\partial x^j_{\\rm Sph}}{\\partial x^i_{\\rm Curv}},\n",
"$$\n",
"\n",
"and $\\partial x^i/\\partial r$ can be computed from the inverse of this matrix (computed within `NRPy+`):\n",
"\n",
"$$\n",
"(J^{-1})_j^i = \\frac{\\partial x^j_{\\rm Curv}}{\\partial x^i_{\\rm Sph}}.\n",
"$$\n",
"\n",
"For example let's set the coordinate system to Cylindrical, as exact expressions do exist for cylindrical coordinates in terms of spherical coordinates for this case:\n",
"\n",
"* $\\rho = r \\sin\\theta$\n",
"* $\\phi = \\phi$\n",
"* $z = r \\cos\\theta$\n",
"\n",
"Thus\n",
"\n",
"$$\n",
"\\partial_r \\rho = \\sin\\theta = \\frac{\\rho}{\\sqrt{\\rho^2 + z^2}}.\n",
"$$\n",
"Note that the negative root is not considered, as $\\theta \\in [0,\\pi]$.\n",
"\n",
"Again these expressions for cylindrical coordinates in terms of spherical *do not exist* within `NRPy+` (as they are not generally easy to describe). Instead we'll use the inverse Jacobian trick to confirm that the `NRPy+` expression is correct:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Arbitrary-offset finite-difference derivatives\n",
"\n",
"See this [Wikipedia article](https://en.wikipedia.org/wiki/Finite_difference_coefficient) for validating coefficients."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"scrolled": false
},
"outputs": [],
"source": [
"# import finite_difference as fin\n",
"# import sympy as sp\n",
"\n",
"def get_arb_offset_FD_coeffs_indices(FDORDER, offset, deriv):\n",
" # deriv = 1 <-- 1st derivative\n",
" Minv = fin.setup_FD_matrix__return_inverse(FDORDER+1, offset)\n",
" indices = []\n",
" coeffs = []\n",
" for i in range(FDORDER+1):\n",
" indices.append(i-int(FDORDER/2) + offset)\n",
" coeffs.append(Minv[i, deriv])\n",
" return coeffs, indices\n",
"\n",
"# FDORDER=4\n",
"# for offset in range(-int(FDORDER/2), int(FDORDER/2)+1):\n",
"# print(get_arb_offset_FD_coeffs_indices(FDORDER, offset, 1))"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"def setup_Cfunction_FD1_arbitrary_upwind(dirn, BC_FDORDER=-1):\n",
" default_FDORDER = par.parval_from_str(\"finite_difference::FD_CENTDERIVS_ORDER\")\n",
" if BC_FDORDER == -1:\n",
" BC_FDORDER = default_FDORDER\n",
"\n",
" par.set_parval_from_str(\"finite_difference::FD_CENTDERIVS_ORDER\", BC_FDORDER)\n",
"\n",
" includes = []\n",
" desc = \"Compute 1st derivative finite-difference derivative with arbitrary upwind\"\n",
" c_type = \"static inline REAL\"\n",
" name = \"FD1_arbitrary_upwind_x\"+str(dirn)+\"_dirn\"\n",
" params = \"\"\"const paramstruct *restrict params, const REAL *restrict gf,\n",
" const int i0,const int i1,const int i2, const int offset\"\"\"\n",
" body = r\"\"\"switch(offset) {\n",
"\"\"\"\n",
" tmp_list = []\n",
" for offset in range(0, int(BC_FDORDER / 2) + 1):\n",
" tmp_list.append(offset)\n",
" if offset > 0:\n",
" tmp_list.append(-offset)\n",
"\n",
" for offset in tmp_list:\n",
" body += \"case \" + str(offset) + \":\\n\"\n",
" body += \" return (\"\n",
" coeffs, indices = get_arb_offset_FD_coeffs_indices(BC_FDORDER, offset, 1)\n",
" for i, coeff in enumerate(coeffs):\n",
" if coeff == 0:\n",
" continue # skip this iteration if coeff=0\n",
" offset = str(indices[i])\n",
" if i > 0:\n",
" body += \" \"\n",
" if offset == \"0\":\n",
" body += \"+\"+str(sp.ccode(coeff))+\"*gf[IDX3S(i0,i1,i2)]\\n\"\n",
" else:\n",
" if dirn == 0:\n",
" body += \"+\"+str(sp.ccode(coeff))+\"*gf[IDX3S(i0+\"+offset+\",i1,i2)]\\n\"\n",
" elif dirn == 1:\n",
" body += \"+\"+str(sp.ccode(coeff))+\"*gf[IDX3S(i0,i1+\"+offset+\",i2)]\\n\"\n",
" elif dirn == 2:\n",
" body += \"+\"+str(sp.ccode(coeff))+\"*gf[IDX3S(i0,i1,i2+\"+offset+\")]\\n\"\n",
" body = body[:-1].replace(\"+-\", \"-\") + \") * invdx\"+str(dirn)+\";\\n\"\n",
" body += \"\"\"}\n",
"return 0.0 / 0.0; // poison output if offset computed incorrectly\n",
"\"\"\"\n",
" rel_path_to_Cparams = os.path.join(\".\")\n",
"\n",
" _prototype, func = Cfunction(includes=includes,\n",
" desc=desc,\n",
" c_type=c_type, name=name, params=params,\n",
" body=indent_Ccode(body, \" \"),\n",
" rel_path_to_Cparams=rel_path_to_Cparams)\n",
"\n",
" par.set_parval_from_str(\"finite_difference::FD_CENTDERIVS_ORDER\", default_FDORDER)\n",
"\n",
" return func"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### $\\partial_r f$ term:\n",
"\n",
"As this boundary condition must be applicable to *any* curvilinear coordinate system (and not just spherical coordinates), $\\partial_r f$ must be computed via\n",
"\n",
"$$\n",
"\\partial_r f = \\frac{\\partial x^i}{\\partial r} \\partial_{i} f,\n",
"$$\n",
"\n",
"where $\\partial_i f$ is evaluated using finite-difference derivatives, suitably upwinded to avoid extending beyond the domain of the grid.\n",
"\n",
"The term $\\partial x^i/\\partial r$ cannot be immediately computed, as $x^i$ is never written as a function of $r$. However, in all curvilinear coordinates $r$, $\\theta$, and $\\phi$ are constructed from $x^i$. Thus we have exact expressions for the Jacobian:\n",
"\n",
"$$\n",
"J_i^j = \\frac{\\partial x^j_{\\rm Sph}}{\\partial x^i_{\\rm Curv}},\n",
"$$\n",
"\n",
"and $\\partial x^i/\\partial r$ can be computed from the inverse of this matrix (computed within `NRPy+`):\n",
"\n",
"$$\n",
"(J^{-1})_j^i = \\frac{\\partial x^j_{\\rm Curv}}{\\partial x^i_{\\rm Sph}}.\n",
"$$\n",
"\n",
"For example let's set the coordinate system to Cylindrical, as exact expressions do exist for cylindrical coordinates in terms of spherical coordinates for this case:\n",
"\n",
"* $\\rho = r \\sin\\theta$\n",
"* $\\phi = \\phi$\n",
"* $z = r \\cos\\theta$\n",
"\n",
"Thus\n",
"\n",
"$$\n",
"\\partial_r \\rho = \\sin\\theta = \\frac{\\rho}{\\sqrt{\\rho^2 + z^2}}.\n",
"$$\n",
"Note that the negative root is not considered, as $\\theta \\in [0,\\pi]$.\n",
"\n",
"Again these expressions for cylindrical coordinates in terms of spherical *do not exist* within `NRPy+` (as they are not generally easy to describe). Instead we'll use the inverse Jacobian trick to confirm that the `NRPy+` expression is correct:"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"xx0/sqrt(xx0**2 + xx2**2)\n"
]
}
],
"source": [
"par.set_parval_from_str(\"reference_metric::CoordSystem\", \"Cylindrical\")\n",
"rfm.reference_metric()\n",
"\n",
"def compute_Jacobian_and_inverseJacobian_tofrom_Spherical():\n",
" # Step 2.a: First construct Jacobian matrix:\n",
" Jac_dUSph_dDrfmUD = ixp.zerorank2()\n",
" for i in range(3):\n",
" for j in range(3):\n",
" Jac_dUSph_dDrfmUD[i][j] = sp.diff(rfm.xxSph[i], rfm.xx[j])\n",
" Jac_dUrfm_dDSphUD, dummyDET = ixp.generic_matrix_inverter3x3(Jac_dUSph_dDrfmUD)\n",
" return Jac_dUSph_dDrfmUD, Jac_dUrfm_dDSphUD\n",
"\n",
"# Jac_dUrfm_dDSphUD is the inverse Jacobian\n",
"Jac_dUSph_dDrfmUD, Jac_dUrfm_dDSphUD = compute_Jacobian_and_inverseJacobian_tofrom_Spherical()\n",
"\n",
"print(sp.simplify(Jac_dUrfm_dDSphUD[0][0]))"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"def setup_Cfunction_compute_partial_r_f(BC_FDORDER=-1):\n",
" desc = \"Compute \\partial_r f\"\n",
" c_type = \"static inline REAL\"\n",
" name = \"compute_partial_r_f\"\n",
" params = \"\"\"const paramstruct *restrict params, REAL *restrict xx[3], const REAL *restrict gfs,\n",
" const int which_gf, const int dest_i0,const int dest_i1,const int dest_i2,\n",
" const int FACEi0,const int FACEi1,const int FACEi2,\n",
" const REAL partial_x0_r, const REAL partial_x1_r, const REAL partial_x2_r\"\"\"\n",
" Jac_dUSph_dDrfmUD, Jac_dUrfm_dDSphUD = compute_Jacobian_and_inverseJacobian_tofrom_Spherical()\n",
"\n",
" default_FDORDER = par.parval_from_str(\"finite_difference::FD_CENTDERIVS_ORDER\")\n",
" if BC_FDORDER == -1:\n",
" BC_FDORDER = default_FDORDER\n",
"\n",
" body = r\"\"\" ///////////////////////////////////////////////////////////\n",
"\n",
" // FD1_stencil_radius = BC_FDORDER/2 = \"\"\" + str(int(BC_FDORDER/2)) + r\"\"\"\n",
" const int FD1_stencil_radius = \"\"\" + str(int(BC_FDORDER/2)) + r\"\"\";\n",
"\n",
" const int ntot = Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1*Nxx_plus_2NGHOSTS2;\n",
"\n",
" ///////////////////////////////////////////////////////////\n",
" // Next we'll compute partial_xi f, using a maximally-centered stencil.\n",
" // The {i0,i1,i2}_offset parameters set the offset of the maximally-centered\n",
" // stencil, such that an offset=0 implies a centered stencil.\n",
"\n",
" // CHECK: Nxx_plus_2NGHOSTS0=10; FD1_stencil_radius=2. Then Nxx_plus_2NGHOSTS0-FD1_stencil_radius-1 = 7\n",
" // if dest_i0 = 9, we get i0_offset=7-9=-2, so the (4th order) deriv\n",
" // stencil is: -4,-3,-2,-1,0\n",
"\n",
" // CHECK: if FD1_stencil_radius=2 and dest_i0 = 1, we get i0_offset = FD1_stencil_radius-dest_i0 = 1,\n",
" // so the (4th order) deriv stencil is: -1,0,1,2,3\n",
"\n",
" // CHECK: if FD1_stencil_radius=2 and dest_i0 = 0, we get i0_offset = FD1_stencil_radius-1 = 2,\n",
" // so the (4th order) deriv stencil is: 0,1,2,3,4\n",
"\"\"\"\n",
" for i in range(3):\n",
" si = str(i)\n",
" if check_zero(Jac_dUrfm_dDSphUD[i][0]):\n",
" body += \" const REAL partial_x\"+si+\"_f=0.0;\\n\"\n",
" else:\n",
" body += \" int i\"+si+\"_offset = FACEi\"+si+\"; // up/downwind on the faces. This offset should never go out of bounds.\\n\"\n",
" body += \" if(dest_i\"+si+\" < FD1_stencil_radius) i\"+si+\"_offset = FD1_stencil_radius-dest_i\"+si+\";\\n\"\n",
" body += \" else if(dest_i\"+si+\" > (Nxx_plus_2NGHOSTS\"+si+\"-FD1_stencil_radius-1)) i\"+si+\"_offset = (Nxx_plus_2NGHOSTS\"+si+\"-FD1_stencil_radius-1) - dest_i\"+si+\";\\n\"\n",
" body += \" const REAL partial_x\"+si+\"_f=FD1_arbitrary_upwind_x\"+si+\"_dirn(params,&gfs[which_gf*ntot],dest_i0,dest_i1,dest_i2,i\"+si+\"_offset);\\n\\n\"\n",
" body += r\"\"\" return partial_x0_r*partial_x0_f + partial_x1_r*partial_x1_f + partial_x2_r*partial_x2_f;\n",
"\"\"\"\n",
" rel_path_to_Cparams = os.path.join(\".\")\n",
"\n",
" _prototype, func = Cfunction(\n",
" includes=[], desc=desc, c_type=c_type, name=name, params=params, body=body,\n",
" rel_path_to_Cparams=rel_path_to_Cparams)\n",
" return func"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"def setup_Cfunction_r_and_partial_xi_r_derivs():\n",
" desc = \"Compute r(xx0,xx1,xx2).\"\n",
" c_type = \"static inline void\"\n",
" name = \"r_and_partial_xi_r_derivs\"\n",
" params = \"\"\"const paramstruct *restrict params,const REAL xx0,const REAL xx1,const REAL xx2,\n",
" REAL *r, REAL *partial_x0_r,REAL *partial_x1_r,REAL *partial_x2_r\"\"\"\n",
" Jac_dUSph_dDrfmUD, Jac_dUrfm_dDSphUD = compute_Jacobian_and_inverseJacobian_tofrom_Spherical()\n",
" body = outputC([rfm.xxSph[0],\n",
" sp.simplify(Jac_dUrfm_dDSphUD[0][0]),\n",
" sp.simplify(Jac_dUrfm_dDSphUD[1][0]),\n",
" sp.simplify(Jac_dUrfm_dDSphUD[2][0])],\n",
" [\"*r\", \"*partial_x0_r\", \"*partial_x1_r\", \"*partial_x2_r\"], filename=\"returnstring\",\n",
" params=\"preindent=1,outCverbose=False,includebraces=False\")\n",
" rel_path_to_Cparams = os.path.join(\".\")\n",
"\n",
" _prototype, func = Cfunction(\n",
" includes=[], desc=desc, c_type=c_type, name=name, params=params, body=body,\n",
" rel_path_to_Cparams=rel_path_to_Cparams)\n",
" return func"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Evaluating the constant $k$ in the $k/r^3$ term\n",
"\n",
"First we note that in the *ansatz*, if $f(r,t)$ perfectly satisfied the outgoing wave equation, then $k=0$ and\n",
"\n",
"$$\n",
"\\left[\\partial_t f\\right]_{\\rm Outgoing\\ Wave} = -\\frac{c}{r} \\left[r \\partial_r f + \\left(f - f_{r\\to\\infty}\\right)\\right].\n",
"$$\n",
"\n",
"However, note that $\\partial_t f$ is evaluated *at all points in the interior* of the grid, so the difference between perfect satisfaction of the radially outgoing wave equation and the actual solution is known at these points. Let's call this difference $\\xi$:\n",
"\n",
"$$\n",
"\\xi = \\partial_t f - \\left[\\partial_t f\\right]_{\\rm Outgoing\\ Wave} \\equiv \\frac{k}{r^3}.\n",
"$$\n",
"\n",
"We compute $\\xi$ at a neighboring interior point with $r=r_{\\rm int}$, $\\xi_{\\rm int}$, so that\n",
"\n",
"$$\n",
"\\xi_{\\rm int} = [\\partial_t f]_{\\rm int} - \\left[\\partial_t f\\right]_{\\rm Outgoing\\ Wave,\\ int} \\equiv \\frac{k}{r_{\\rm int}^3}.\n",
"$$\n",
"\n",
"In this way we obtain $k$:\n",
"\n",
"$$\n",
"k = r_{\\rm int}^3 \\left([\\partial_t f]_{\\rm int} - \\left[\\partial_t f\\right]_{\\rm Outgoing\\ Wave,\\ int}\\right)\n",
"$$\n",
"\n",
"To determine the appropriate interior point, we simply keep track of the current boundary face being updated and choose the nearest neighbor in the direction of the grid interior."
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"def setup_Cfunction_radiation_bcs_curvilinear(BC_FDORDER=-1):\n",
" includes = []\n",
" prefunc = \"\"\n",
" Jac_dUSph_dDrfmUD, Jac_dUrfm_dDSphUD = compute_Jacobian_and_inverseJacobian_tofrom_Spherical()\n",
" for i in range(3):\n",
" # Do not generate FD1_arbitrary_upwind_xj_dirn() if the symbolic expression for dxj/dr == 0!\n",
" if not check_zero(Jac_dUrfm_dDSphUD[i][0]):\n",
" prefunc += setup_Cfunction_FD1_arbitrary_upwind(dirn=i, BC_FDORDER=BC_FDORDER)\n",
" prefunc += setup_Cfunction_r_and_partial_xi_r_derivs()\n",
" prefunc += setup_Cfunction_compute_partial_r_f(BC_FDORDER=BC_FDORDER)\n",
" desc = r\"\"\"*** Apply radiation BCs to all outer boundaries. ***\n",
"\"\"\"\n",
" c_type = \"static inline void\"\n",
" name = \"radiation_bcs_curvilinear\"\n",
" params = \"\"\"const paramstruct *restrict params, const bc_struct *restrict bcstruct,REAL *restrict xx[3],\n",
" const REAL *restrict gfs, REAL *restrict gfs_rhss,\n",
" const int which_gf, const int dest_i0,const int dest_i1,const int dest_i2,\n",
" const int FACEi0,const int FACEi1,const int FACEi2\"\"\"\n",
" body = r\"\"\"// Nearest \"interior\" neighbor of this gridpoint, based on current face\n",
"const int dest_i0_int=dest_i0+1*FACEi0, dest_i1_int=dest_i1+1*FACEi1, dest_i2_int=dest_i2+1*FACEi2;\n",
"REAL r, partial_x0_r,partial_x1_r,partial_x2_r;\n",
"REAL r_int, partial_x0_r_int,partial_x1_r_int,partial_x2_r_int;\n",
"r_and_partial_xi_r_derivs(params,xx[0][dest_i0], xx[1][dest_i1], xx[2][dest_i2], &r, &partial_x0_r, &partial_x1_r, &partial_x2_r);\n",
"r_and_partial_xi_r_derivs(params,xx[0][dest_i0_int],xx[1][dest_i1_int],xx[2][dest_i2_int],&r_int, &partial_x0_r_int,&partial_x1_r_int,&partial_x2_r_int);\n",
"const REAL partial_r_f = compute_partial_r_f(params,xx,gfs, which_gf,dest_i0, dest_i1, dest_i2,\n",
" FACEi0,FACEi1,FACEi2,\n",
" partial_x0_r, partial_x1_r, partial_x2_r);\n",
"const REAL partial_r_f_int = compute_partial_r_f(params,xx,gfs, which_gf,dest_i0_int,dest_i1_int,dest_i2_int,\n",
" FACEi0,FACEi1,FACEi2,\n",
" partial_x0_r_int,partial_x1_r_int,partial_x2_r_int);\n",
"\n",
"const int idx3 = IDX3S(dest_i0,dest_i1,dest_i2);\n",
"const int idx3_int = IDX3S(dest_i0_int,dest_i1_int,dest_i2_int);\n",
"\n",
"const REAL partial_t_f_int = gfs_rhss[IDX4ptS(which_gf, idx3_int)];\n",
"\n",
"const REAL c = gridfunctions_wavespeed[which_gf];\n",
"const REAL f_infinity = gridfunctions_f_infinity[which_gf];\n",
"const REAL f = gfs[IDX4ptS(which_gf, idx3)];\n",
"const REAL f_int = gfs[IDX4ptS(which_gf, idx3_int)];\n",
"const REAL partial_t_f_int_outgoing_wave = -c * (partial_r_f_int + (f_int - f_infinity) / r_int);\n",
"\n",
"const REAL k = r_int*r_int*r_int * (partial_t_f_int - partial_t_f_int_outgoing_wave);\n",
"\n",
"const REAL rinv = 1.0 / r;\n",
"const REAL partial_t_f_outgoing_wave = -c * (partial_r_f + (f - f_infinity) * rinv);\n",
"\n",
"gfs_rhss[IDX4ptS(which_gf, idx3)] = partial_t_f_outgoing_wave + k * rinv*rinv*rinv;\n",
"\"\"\"\n",
" rel_path_to_Cparams = os.path.join(\".\")\n",
"\n",
" _prototype, func = Cfunction(\n",
" includes=includes, prefunc=prefunc,\n",
" desc=desc,\n",
" c_type=c_type, name=name, params=params,\n",
" body=indent_Ccode(body, \" \"),\n",
" rel_path_to_Cparams=rel_path_to_Cparams)\n",
" return func"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 9: Main driver function: `apply_bcs_curvilinear()`: quickly apply boundary and parity conditions with information from `bcstruct` \\[Back to [top](#toc)\\]\n",
"$$\\label{apply_bcs_curvilinear}$$\n",
"\n",
"`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. Outer boundary grid points are filled using either quadratic polynomial extrapolation (`\"extrapolation\"`) or radiation (`\"radiation\"`) boundary conditions, described previously."
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"execution": {
"iopub.execute_input": "2021-09-08T15:47:15.264277Z",
"iopub.status.busy": "2021-09-08T15:47:15.263692Z",
"iopub.status.idle": "2021-09-08T15:47:15.265449Z",
"shell.execute_reply": "2021-09-08T15:47:15.265065Z"
}
},
"outputs": [],
"source": [
"def add_to_Cfunction_dict_apply_bcs_curvilinear(outer_bcs_type=\"extrapolation\", BC_FDORDER=-1):\n",
" includes = [\"NRPy_basic_defines.h\", \"NRPy_function_prototypes.h\"]\n",
" prefunc = \"\"\n",
" if outer_bcs_type == \"extrapolation\":\n",
" prefunc = r\"\"\"// Declare boundary condition EXTRAP_BC_UPDATE_OUTER macro,\n",
"// which updates a single outer boundary point\n",
"// of the 3D grid cube using quadratic polynomial\n",
"// extrapolation.\n",
"#define EXTRAP_BC_UPDATE_OUTER(which_gf, i0,i1,i2, FACEX0,FACEX1,FACEX2) { \\\n",
" const int idx3 = IDX3S(i0,i1,i2); \\\n",
" gfs[IDX4S(which_gf,i0,i1,i2)] = \\\n",
" +3.0*gfs[IDX4S(which_gf,i0+1*FACEX0,i1+1*FACEX1,i2+1*FACEX2)] \\\n",
" -3.0*gfs[IDX4S(which_gf,i0+2*FACEX0,i1+2*FACEX1,i2+2*FACEX2)] \\\n",
" +1.0*gfs[IDX4S(which_gf,i0+3*FACEX0,i1+3*FACEX1,i2+3*FACEX2)]; \\\n",
" }\n",
"\"\"\"\n",
" elif outer_bcs_type == \"radiation\":\n",
" prefunc = setup_Cfunction_radiation_bcs_curvilinear(BC_FDORDER=BC_FDORDER)\n",
" else:\n",
" print(\"outer_bcs_type == \" + outer_bcs_type + \" NOT SUPPORTED.\")\n",
" sys.exit(1)\n",
" desc = r\"\"\"Curvilinear boundary condition driver routine: Apply BCs to all six\n",
" boundary faces of the 3D numerical domain, filling in the\n",
" innermost ghost zone layer first, and moving outward.\n",
"\"\"\"\n",
" c_type = \"void\"\n",
" name = \"apply_bcs_curvilinear\"\n",
" params = \"\"\"const paramstruct *restrict params, const bc_struct *restrict bcstruct,\n",
" const int NUM_GFS, const int8_t *restrict gfs_parity, REAL *restrict xx[3],\n",
" REAL *restrict gfs, REAL *restrict gfs_rhss\"\"\"\n",
" body = r\"\"\"#pragma omp parallel for\n",
" for(int which_gf=0;which_gfnum_ob_gz_pts[which_gz];pt++) {\n",
"\"\"\"\n",
" if outer_bcs_type == \"radiation\":\n",
" body += r\"\"\" // *** Apply radiation BCs to all outer boundary points. ***\n",
" radiation_bcs_curvilinear(params, bcstruct, xx, gfs, gfs_rhss, which_gf,\n",
" bcstruct->outer[which_gz][pt].outer_bc_dest_pt.i0,\n",
" bcstruct->outer[which_gz][pt].outer_bc_dest_pt.i1,\n",
" bcstruct->outer[which_gz][pt].outer_bc_dest_pt.i2,\n",
" bcstruct->outer[which_gz][pt].FACEi0,\n",
" bcstruct->outer[which_gz][pt].FACEi1,\n",
" bcstruct->outer[which_gz][pt].FACEi2);\n",
"\"\"\"\n",
" elif outer_bcs_type == \"extrapolation\":\n",
" body += r\"\"\" // *** Apply 2nd-order polynomial extrapolation BCs to all outer boundary points. ***\n",
" EXTRAP_BC_UPDATE_OUTER(which_gf,\n",
" bcstruct->outer[which_gz][pt].outer_bc_dest_pt.i0,\n",
" bcstruct->outer[which_gz][pt].outer_bc_dest_pt.i1,\n",
" bcstruct->outer[which_gz][pt].outer_bc_dest_pt.i2,\n",
" bcstruct->outer[which_gz][pt].FACEi0,\n",
" bcstruct->outer[which_gz][pt].FACEi1,\n",
" bcstruct->outer[which_gz][pt].FACEi2);\n",
"\"\"\"\n",
" body += r\"\"\" }\n",
"\n",
" // Then apply INNER (parity) boundary conditions:\n",
" for(int pt=0;ptnum_ib_gz_pts[which_gz];pt++) {\n",
" const int i0dest = bcstruct->inner[which_gz][pt].inner_bc_dest_pt.i0;\n",
" const int i1dest = bcstruct->inner[which_gz][pt].inner_bc_dest_pt.i1;\n",
" const int i2dest = bcstruct->inner[which_gz][pt].inner_bc_dest_pt.i2;\n",
" const int i0src = bcstruct->inner[which_gz][pt].inner_bc_src_pt.i0;\n",
" const int i1src = bcstruct->inner[which_gz][pt].inner_bc_src_pt.i1;\n",
" const int i2src = bcstruct->inner[which_gz][pt].inner_bc_src_pt.i2;\n",
"\"\"\"\n",
" inner_bc_str = \"\"\" gfs[IDX4S(which_gf,i0dest,i1dest,i2dest)] =\n",
" bcstruct->inner[which_gz][pt].parity[gfs_parity[which_gf]] * gfs[IDX4S(which_gf, i0src,i1src,i2src)];\"\"\"\n",
" if outer_bcs_type == \"radiation\":\n",
" body += inner_bc_str.replace(\"gfs[IDX\", \"gfs_rhss[IDX\")\n",
" else:\n",
" body += inner_bc_str\n",
" body += r\"\"\"\n",
" } // END for(int pt=0;pt\n",
"\n",
"# Step 10: `CurviBC_Playground.c`: Start-to-Finish C code module for testing & validating curvilinear boundary conditions \\[Back to [top](#toc)\\]\n",
"$$\\label{start2finish}$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"## Step 10.a: Register gridfunctions of all 10 parity types; output gridfunction aliases to `CurviBoundaryConditions/gridfunction_defines.h` \\[Back to [top](#toc)\\]\n",
"$$\\label{register_gfs}$$ \n",
"\n",
"Here we \n",
"\n",
"1. Register within NRPy+ one gridfunction per each of the 10 parity conditions, and then \n",
"1. output to `CurviBoundaryConditions/gridfunction_defines.h` the corresponding gridfunction aliases and parity info, so the C code can\n",
" 1. Access each gridfunction by its human-friendly alias (e.g., `test_gfs[RANKONEU0GF][idx]` instead of `test_gfs[6][idx]`), and\n",
" 1. Access each gridfunction parity by the same alias (e.g., `bcstruct->inner_bc_parity[which_gz][pt].parity[gf_parity[RANKONEU0GF]]`)"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"execution": {
"iopub.execute_input": "2021-09-08T15:47:15.269753Z",
"iopub.status.busy": "2021-09-08T15:47:15.269360Z",
"iopub.status.idle": "2021-09-08T15:47:15.271299Z",
"shell.execute_reply": "2021-09-08T15:47:15.270905Z"
}
},
"outputs": [],
"source": [
"# Step 10.a: Register gridfunctions of all 10 parity types\n",
"\n",
"# 6 gridfunctions, corresponding to all unique rank-2 tensor components:\n",
"ranktwosymmDD = ixp.register_gridfunctions_for_single_rank2(\"EVOL\",\"ranktwosymmDD\", \"sym01\", f_infinity=0.0, wavespeed=1.0)\n",
"# 3 gridfunctions, corresponding to all unique rank-1 tensor components:\n",
"rankoneU = ixp.register_gridfunctions_for_single_rank1(\"EVOL\",\"rankoneU\", f_infinity=0.0, wavespeed=1.0)\n",
"# 1 rank-0 (scalar) gridfunction\n",
"rankzero = ixp.gri.register_gridfunctions(\"EVOL\",\"rankzero\", f_infinity=1.0, wavespeed=sp.sqrt(2.0))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"## Step 10.b: Set up test data for Curvilinear Boundary Conditions code validation \\[Back to [top](#toc)\\]\n",
"$$\\label{validate}$$\n",
"\n",
"We will validate this curvilinear boundary condition module by comparing its results with the original (trusted) SENR code, as follows:\n",
"\n",
"* **Discrete data test**:\n",
" 1. Fill all 10 gridfunctions at each gridpoint with the unique gridpoint integer index `IDX3S(i0,i1,i2)`\n",
" 1. Apply curvilinear boundary conditions\n",
" 1. Compare output data at all gridpoints with those from the original SENR code. Agreement should be perfect.\n",
"\n",
"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:\n",
"\n",
"* **Smooth data test** (TODO):\n",
" 1. Fill all 10 gridfunctions with data that are smooth in the Cartesian basis.\n",
" 1. Apply Jacobian transformation to all data points, to convert to curvilinear basis\n",
" 1. Apply curvilinear boundary conditions\n",
" 1. Apply Jacobian transformation to all data points, to convert back to Cartesian basis\n",
" 1. 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."
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"execution": {
"iopub.execute_input": "2021-09-08T15:47:15.274546Z",
"iopub.status.busy": "2021-09-08T15:47:15.274158Z",
"iopub.status.idle": "2021-09-08T15:47:15.276042Z",
"shell.execute_reply": "2021-09-08T15:47:15.275654Z"
}
},
"outputs": [],
"source": [
"def add_to_Cfunction_dict_CurviBC_Discrete_initial_data():\n",
" add_to_Cfunction_dict(\n",
" includes = [\"NRPy_basic_defines.h\", \"NRPy_function_prototypes.h\"],\n",
" desc =\"Populate in_gfs[] with discrete, regular data\",\n",
" c_type =\"void\",\n",
" name =\"CurviBC_Discrete_initial_data\",\n",
" params =\"const paramstruct *restrict params,REAL *restrict in_gfs\",\n",
" body =\"\"\"\n",
"#pragma omp parallel for\n",
" for(int i=0;i\n",
"\n",
"## Step 10.c: `CurviBC_Playground`'s `main.c` Code \\[Back to [top](#toc)\\]\n",
"$$\\label{mainc}$$"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {
"execution": {
"iopub.execute_input": "2021-09-08T15:47:15.281826Z",
"iopub.status.busy": "2021-09-08T15:47:15.281400Z",
"iopub.status.idle": "2021-09-08T15:47:15.283297Z",
"shell.execute_reply": "2021-09-08T15:47:15.282912Z"
}
},
"outputs": [],
"source": [
"def add_to_Cfunction_dict_main__CurviBC_Playground():\n",
" includes = [\"NRPy_basic_defines.h\", \"NRPy_function_prototypes.h\"]\n",
" desc = \"\"\"// main() function:\n",
"// Step 0: Read command-line input, set up grid structure, allocate memory for gridfunctions, set up coordinates\n",
"// Step 1: Write test data to gridfunctions\n",
"// Step 2: Overwrite all data in ghost zones with NaNs\n",
"// Step 3: Apply curvilinear boundary conditions\n",
"// Step 4: Print gridfunction data after curvilinear boundary conditions have been applied\n",
"// Step 5: Free all allocated memory\n",
"\"\"\"\n",
" c_type = \"int\"\n",
" name = \"main\"\n",
" params = \"int argc, const char *argv[]\"\n",
" body = r\"\"\" // Step 0a: Read command-line input, error out if nonconformant\n",
" if(argc != 2) {\n",
" fprintf(stderr,\"Error: Expected one command-line argument: ./CurviBC_Playground [test type: Smooth or Discrete],\\n\");\n",
" exit(1);\n",
" }\n",
"\n",
" griddata_struct griddata;\n",
" set_Cparameters_to_default(&griddata.params);\n",
"\n",
" char CoordSystem_name[100];\n",
" snprintf(CoordSystem_name, 100, \"\"\"\n",
" body += \"\\\"\"+CoordSystem+\"\\\"\"\n",
" body += r\"\"\");\n",
"\n",
" // Step 0b: Set number of gridpoints...\n",
" const int Nxx[3] = { 4, 4, 4 };\n",
"\n",
" // Step 0c: Set test type to Smooth or Discrete\n",
" char test_type[100];\n",
" snprintf(test_type, 100, \"%s\", argv[1]);\n",
" if(strncmp(\"Smooth\", test_type, 100) != 0 &&\n",
" strncmp(\"Discrete\",test_type, 100) != 0) {\n",
" fprintf(stderr,\"Error: test type = %s not supported. Choose Smooth or Discrete (CASE SENSITIVE).\\n\",test_type);\n",
" exit(1);\n",
" }\n",
"\n",
" // Step 0d: Uniform coordinate grids are stored to *xx[3]\n",
" // Step 0d.i: Set bcstruct\n",
" {\n",
" int EigenCoord;\n",
" EigenCoord = 1;\n",
" // Step 0d.ii: Call set_Nxx_dxx_invdx_params__and__xx(), which sets\n",
" // params Nxx,Nxx_plus_2NGHOSTS,dxx,invdx, and xx[] for the\n",
" // chosen Eigen-CoordSystem.\n",
" set_Nxx_dxx_invdx_params__and__xx(EigenCoord, Nxx, &griddata.params, griddata.xx);\n",
" // Step 0e: Find ghostzone mappings; set up bcstruct\n",
" driver_bcstruct(&griddata.params, &griddata.bcstruct, griddata.xx);\n",
" // Step 0e.i: Free allocated space for xx[][] array\n",
" for(int i=0;i<3;i++) free(griddata.xx[i]);\n",
"\n",
" // Step 0f: Call set_Nxx_dxx_invdx_params__and__xx(), which sets\n",
" // params Nxx,Nxx_plus_2NGHOSTS,dxx,invdx, and xx[] for the\n",
" // chosen (non-Eigen) CoordSystem.\n",
" EigenCoord = 0;\n",
" set_Nxx_dxx_invdx_params__and__xx(EigenCoord, Nxx, &griddata.params, griddata.xx);\n",
" }\n",
"\n",
" // Step 0g: Set all C parameters \"blah\" for params.blah, including\n",
" // Nxx_plus_2NGHOSTS0 = params.Nxx_plus_2NGHOSTS0, etc.\n",
" const int Nxx_plus_2NGHOSTS0 = griddata.params.Nxx_plus_2NGHOSTS0;\n",
" const int Nxx_plus_2NGHOSTS1 = griddata.params.Nxx_plus_2NGHOSTS1;\n",
" const int Nxx_plus_2NGHOSTS2 = griddata.params.Nxx_plus_2NGHOSTS2;\n",
"\n",
" // Step 0h: Allocate memory for gridfunctions\n",
" griddata.test_gfs = (REAL *)malloc(sizeof(REAL) * NUM_EVOL_GFS * Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1*Nxx_plus_2NGHOSTS2);\n",
" REAL *restrict test_gfs = griddata.test_gfs;\n",
"\"\"\"\n",
" if outer_bcs_type == \"extrapolation\":\n",
" body += \" REAL *restrict test_gfs_rhss = test_gfs;\"\n",
" elif outer_bcs_type == \"radiation\":\n",
" body += \"\"\" griddata.test_gfs_rhss = (REAL *)malloc(sizeof(REAL) * NUM_EVOL_GFS * Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1*Nxx_plus_2NGHOSTS2);\n",
" REAL *restrict test_gfs_rhss = griddata.test_gfs_rhss;\"\"\"\n",
"\n",
" body += r\"\"\"\n",
" // Step 1: Write test data to gridfunctions\n",
" if(strncmp(\"Discrete\", test_type, 100)==0) {\n",
" CurviBC_Discrete_initial_data(&griddata.params, test_gfs);\n",
" CurviBC_Discrete_initial_data(&griddata.params, test_gfs_rhss);\n",
" } else {\n",
" fprintf(stderr, \"Sorry, curvilinear boundary conditions test = %s not yet supported. Feel free to contribute!\\n\",test_type);\n",
" }\n",
"\n",
" // Step 2: Overwrite all data in ghost zones with NaNs\n",
" LOOP_REGION(0,Nxx_plus_2NGHOSTS0, 0,Nxx_plus_2NGHOSTS1, 0,Nxx_plus_2NGHOSTS2) {\n",
" for(int gf=0;gf= Nxx_plus_2NGHOSTS0-NGHOSTS) test_gfs_rhss[idx4] = +(0.0 / 0.0);\n",
" if(i1 < NGHOSTS || i1 >= Nxx_plus_2NGHOSTS1-NGHOSTS) test_gfs_rhss[idx4] = +(0.0 / 0.0);\n",
" if(i2 < NGHOSTS || i2 >= Nxx_plus_2NGHOSTS2-NGHOSTS) test_gfs_rhss[idx4] = +(0.0 / 0.0);\n",
" }\n",
" }\n",
"\n",
" // Step 3: Apply curvilinear boundary conditions\n",
" apply_bcs_curvilinear(&griddata.params, &griddata.bcstruct, NUM_EVOL_GFS, evol_gf_parity,\n",
" griddata.xx, test_gfs, test_gfs_rhss);\n",
"\n",
" // Step 4: Print gridfunction data after curvilinear boundary conditions have been applied:\n",
" char filename[120]; sprintf(filename,\"out4x4x4-%s-NGHOSTS4oFD.txt\", CoordSystem_name);\n",
" FILE *outfile = fopen(filename, \"w\");\n",
"\n",
" LOOP_REGION(0,Nxx_plus_2NGHOSTS0, 0,Nxx_plus_2NGHOSTS1, 0,Nxx_plus_2NGHOSTS2) {\n",
" fprintf(outfile, \"%d %d %d | \", i0,i1,i2);\n",
" for(int gf=0;gf\n",
"\n",
"## Step 10.d: Add all CurviBC C codes to C function dictionary, and add CurviBC definitions to `NRPy_basic_defines.h` \\[Back to [top](#toc)\\]\n",
"$$\\label{curvibc_setupall}$$"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {
"execution": {
"iopub.execute_input": "2021-09-08T15:47:15.286703Z",
"iopub.status.busy": "2021-09-08T15:47:15.286318Z",
"iopub.status.idle": "2021-09-08T15:47:15.288180Z",
"shell.execute_reply": "2021-09-08T15:47:15.287794Z"
}
},
"outputs": [],
"source": [
"# Always call this after ALL gridfunctions have been registered!\n",
"def CurviBoundaryConditions_register_C_functions_and_NRPy_basic_defines(verbose=True):\n",
" # First register C functions needed by CurviBCs\n",
" add_to_Cfunction_dict_set_bcstruct()\n",
" add_to_Cfunction_dict_driver_bcstruct()\n",
"# add_to_Cfunction_dict_apply_extrap_bcs_curvilinear()\n",
" add_to_Cfunction_dict_freemem_bcstruct()\n",
" # Then set up the dictionary entry for CurviBC in NRPy_basic_defines\n",
" Nbd_str = NRPy_basic_defines_CurviBC_data_structures()\n",
" Nbd_str += NRPy_basic_defines_set_gridfunction_defines_with_parity_types(verbose=verbose)\n",
" outC_NRPy_basic_defines_h_dict[\"CurviBoundaryConditions\"] = Nbd_str"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 11: Generate C code for generating numerical grids \\[Back to [top](#toc)\\]\n",
"$$\\label{ccode_numgridgen}$$"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {
"execution": {
"iopub.execute_input": "2021-09-08T15:47:15.295762Z",
"iopub.status.busy": "2021-09-08T15:47:15.294840Z",
"iopub.status.idle": "2021-09-08T15:47:18.611539Z",
"shell.execute_reply": "2021-09-08T15:47:18.611295Z"
},
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Evolved gridfunction \"rankoneU0\" has parity type 1.\n",
"Evolved gridfunction \"rankoneU1\" has parity type 2.\n",
"Evolved gridfunction \"rankoneU2\" has parity type 3.\n",
"Evolved gridfunction \"ranktwosymmDD00\" has parity type 4.\n",
"Evolved gridfunction \"ranktwosymmDD01\" has parity type 5.\n",
"Evolved gridfunction \"ranktwosymmDD02\" has parity type 6.\n",
"Evolved gridfunction \"ranktwosymmDD11\" has parity type 7.\n",
"Evolved gridfunction \"ranktwosymmDD12\" has parity type 8.\n",
"Evolved gridfunction \"ranktwosymmDD22\" has parity type 9.\n",
"Evolved gridfunction \"rankzero\" has parity type 0.\n"
]
}
],
"source": [
"import outputC as outC\n",
"# Then we set the coordinate system for the numerical grid\n",
"par.set_parval_from_str(\"reference_metric::CoordSystem\", CoordSystem)\n",
"rfm.reference_metric()\n",
"# Step 5: Generate & register set_Nxx_dxx_invdx_params__and__xx(), which sets\n",
"# params Nxx,Nxx_plus_2NGHOSTS,dxx,invdx, and xx[] for the\n",
"# chosen (not necessarily Eigen-) CoordSystem.\n",
"rfm.add_to_Cfunc_dict_set_Nxx_dxx_invdx_params__and__xx()\n",
"\n",
"add_to_Cfunction_dict_apply_bcs_curvilinear(outer_bcs_type=outer_bcs_type,\n",
" BC_FDORDER=outer_bcs_FDORDER)\n",
"add_to_Cfunction_dict_set_up__bc_gz_map_and_parity_condns()\n",
"\n",
"\n",
"outC.outputC_register_C_functions_and_NRPy_basic_defines() # #define M_PI, etc.\n",
"# Declare paramstruct, register set_Cparameters_to_default(),\n",
"# and output declare_Cparameters_struct.h and set_Cparameters[].h:\n",
"outC.NRPy_param_funcs_register_C_functions_and_NRPy_basic_defines(os.path.join(Ccodesrootdir))\n",
"\n",
"gri.register_C_functions_and_NRPy_basic_defines(enable_griddata_struct=True,\n",
" enable_bcstruct_in_griddata_struct=True,\n",
" enable_rfmstruct=False,\n",
" enable_MoL_gfs_struct=False,\n",
" extras_in_griddata_struct=[\"REAL *restrict test_gfs\",\n",
" \"REAL *restrict test_gfs_rhss\", ]) # #define IDX3S(), etc.\n",
"fin.register_C_functions_and_NRPy_basic_defines(NGHOSTS_account_for_onezone_upwind=True) # #define NGHOSTS, etc.\n",
"rfm.register_NRPy_basic_defines()\n",
"CurviBoundaryConditions_register_C_functions_and_NRPy_basic_defines()\n",
"# Comment out above line and uncomment below line to confirm independent Python module agrees.\n",
"# import CurviBoundaryConditions.CurviBoundaryConditions_new_way as CBC\n",
"# CBC.CurviBoundaryConditions_register_C_functions_and_NRPy_basic_defines()\n",
"\n",
"# initial data function:\n",
"add_to_Cfunction_dict_CurviBC_Discrete_initial_data()\n",
"\n",
"# main function:\n",
"add_to_Cfunction_dict_main__CurviBC_Playground()\n",
"\n",
"outC.construct_NRPy_basic_defines_h(Ccodesrootdir)\n",
"outC.construct_NRPy_function_prototypes_h(Ccodesrootdir)\n",
"# The following is run from inside cmd.new_C_compile() in the next code cell.\n",
"# outC.construct_Makefile_from_outC_function_dict(Ccodesrootdir, \"CurviBC_Playground\")\n",
"# print(outC.outC_function_dict[\"set_Nxx_dxx_invdx_params__and__xx\"])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 12: Validation: Compile & compare with original (trusted) SENR results \\[Back to [top](#toc)\\]\n",
"$$\\label{senr_compare}$$"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {
"execution": {
"iopub.execute_input": "2021-09-08T15:47:18.614645Z",
"iopub.status.busy": "2021-09-08T15:47:18.614349Z",
"iopub.status.idle": "2021-09-08T15:47:19.240208Z",
"shell.execute_reply": "2021-09-08T15:47:19.239773Z"
},
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(EXEC): Executing `make -j34`...\n",
"(BENCH): Finished executing in 0.2019340991973877 seconds.\n",
"Finished compilation.\n",
"(EXEC): Executing `taskset -c 0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15 ./CurviBC_Playground Discrete`...\n",
"(BENCH): Finished executing in 0.20223331451416016 seconds.\n"
]
}
],
"source": [
"import cmdline_helper as cmd\n",
"# from outputC import construct_Makefile_from_outC_function_dict\n",
"# construct_Makefile_from_outC_function_dict(Ccodesrootdir, \"CurviBC_Playground\", compiler_opt_option=\"fast\")\n",
"cmd.new_C_compile(Ccodesrootdir, \"CurviBC_Playground\", compiler_opt_option=\"fast\") # fastdebug or debug also supported\n",
"os.chdir(Ccodesrootdir)\n",
"cmd.Execute(\"CurviBC_Playground\", \"Discrete\")\n",
"#cmd.Execute(\"CurviBC_Playground\", \"4 4 4 Discrete\", \"out4x4x4-Spherical-NGHOSTS4oFD.txt\")\n",
"os.chdir(\"..\")"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {
"execution": {
"iopub.execute_input": "2021-09-08T15:47:19.244765Z",
"iopub.status.busy": "2021-09-08T15:47:19.244251Z",
"iopub.status.idle": "2021-09-08T15:47:19.246131Z",
"shell.execute_reply": "2021-09-08T15:47:19.246437Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Spherical boundary condition comparison test between this tutorial notebook & trusted original SENR code: PASSED\n"
]
}
],
"source": [
"import filecmp\n",
"if outer_bcs_type == \"extrapolation\":\n",
" if \"Cylindrical\" in CoordSystem:\n",
" if filecmp.cmp(os.path.join(Ccodesrootdir, \"out4x4x4-\"+CoordSystem+\"-NGHOSTS4oFD.txt\"),\n",
" os.path.join(\"CurviBoundaryConditions\", \"SENRout4x4x4-Cylindrical_NGHOSTS4oFD.txt\")) == False:\n",
" print(\"ERROR: \"+CoordSystem+\" boundary conditions malfunction!\")\n",
" sys.exit(1)\n",
" elif \"Spherical\" in CoordSystem:\n",
" if filecmp.cmp(os.path.join(Ccodesrootdir, \"out4x4x4-\"+CoordSystem+\"-NGHOSTS4oFD.txt\"),\n",
" os.path.join(\"CurviBoundaryConditions\", \"SENRout4x4x4-Spherical_NGHOSTS4oFD.txt\")) == False:\n",
" print(\"ERROR: \"+CoordSystem+\" boundary conditions malfunction!\")\n",
" sys.exit(1)\n",
" else:\n",
" print(\"ERROR: \"+CoordSystem+\" coordinate system comparison unavailable!\")\n",
" sys.exit(1)\n",
" print(CoordSystem + \" boundary condition comparison test between this tutorial notebook & trusted original SENR code: PASSED\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 12: Output this notebook to $\\LaTeX$-formatted PDF file \\[Back to [top](#toc)\\]\n",
"$$\\label{latex_pdf_output}$$\n",
"\n",
"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\n",
"[Tutorial-Start_to_Finish-Curvilinear_BCs_new_way.pdf](Tutorial-Start_to_Finish-Curvilinear_BCs_new_way.pdf) (Note that clicking on this link may not work; you may need to open the PDF file through another means.)"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {
"execution": {
"iopub.execute_input": "2021-09-08T15:47:19.249219Z",
"iopub.status.busy": "2021-09-08T15:47:19.248835Z",
"iopub.status.idle": "2021-09-08T15:47:22.888138Z",
"shell.execute_reply": "2021-09-08T15:47:22.888544Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Created Tutorial-Start_to_Finish-Curvilinear_BCs_new_way.tex, and compiled\n",
" LaTeX file to PDF file Tutorial-Start_to_Finish-\n",
" Curvilinear_BCs_new_way.pdf\n"
]
}
],
"source": [
"import cmdline_helper as cmd # NRPy+: Multi-platform Python command-line interface\n",
"cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-Start_to_Finish-Curvilinear_BCs_new_way\")"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.0"
}
},
"nbformat": 4,
"nbformat_minor": 2
}