{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<script async src=\"https://www.googletagmanager.com/gtag/js?id=UA-59152712-8\"></script>\n",
"<script>\n",
" window.dataLayer = window.dataLayer || [];\n",
" function gtag(){dataLayer.push(arguments);}\n",
" gtag('js', new Date());\n",
"\n",
" gtag('config', 'UA-59152712-8');\n",
"</script>\n",
"\n",
"# `GiRaFFE_NRPy`: A-to-B code\n",
"\n",
"## Author: Patrick Nelson\n",
"\n",
"<a id='intro'></a>\n",
"\n",
"**Notebook Status:** <font color=green><b> Validated </b></font>\n",
"\n",
"**Validation Notes:** This code produces the expected magnetic fields for arbitrary vector potentials.\n",
"\n",
"# This module presents the functionality of [GiRaFFE_NRPy_A2B.py](../../edit/in_progress/GiRaFFE_NRPy/GiRaFFE_NRPy_A2B.py).\n",
"\n",
"## Introduction: \n",
"This writes and documents the C code that `GiRaFFE_NRPy` uses to compute the magnetic fields from the vector potential. This is a relatively straightforward calculation, but requires care to be taken in the ghost zones. \n",
"\n",
"We will need to compute $B^i$ everywhere in order to evolve $\\tilde{S}_i$. However, $B^i = \\epsilon^{ijk} \\partial_j A_k$ requires derivatives of $A_i$, so getting $B^i$ in the ghostzones (and not just on the interior) will require some finesse. To do so, we will first compute the derivatives on the interior normally. Then, in the ghost zones, we will check if the point is on any face. If it is, we will use an appropriately-shifted stencil to avoid accessing points that are not on our grid. This will let us compute the derivative at the outermost gridpoints.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='toc'></a>\n",
"\n",
"# Table of Contents\n",
"$$\\label{toc}$$\n",
"\n",
"This notebook is organized as follows\n",
"\n",
"1. [Step 1](#magnetic_field): The Magnetic Field in terms of the Vector Potential\n",
"1. [Step 2](#ccode): Write out C code\n",
" 1. [Step 2.a](#ghostzones): Compute the magnetic field at a point in a ghostzone\n",
" 1. [Step 2.b](#driver): Looping over all ghostzones\n",
"1. [Step 3](#code_validation): Code Validation against original C code\n",
"1. [Step 4](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"# Step 0: Add NRPy's directory to the path\n",
"# https://stackoverflow.com/questions/16780014/import-file-from-parent-directory\n",
"import os,sys\n",
"nrpy_dir_path = os.path.join(\"..\")\n",
"if nrpy_dir_path not in sys.path:\n",
" sys.path.append(nrpy_dir_path)\n",
"\n",
"import cmdline_helper as cmd\n",
"outdir = \"GiRaFFE_NRPy/GiRaFFE_Ccode_validation/A2B/\"\n",
"cmd.mkdir(outdir)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='magnetic_field'></a>\n",
"\n",
"## Step 1: The Magnetic Field in terms of the Vector Potential \\[Back to [top](#toc)\\]\n",
"$$\\label{magnetic_field}$$\n",
"\n",
"We start in the usual way - import the modules we need. We will also import the Levi-Civita symbol from `indexedexp.py` and use it to set the Levi-Civita tensor $\\epsilon^{ijk} = [ijk]/\\sqrt{\\gamma}$."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# Step 1: The A-to-B driver\n",
"from outputC import outCfunction, lhrh # NRPy+: Core C code output module\n",
"import finite_difference as fin # NRPy+: Finite difference C code generation module\n",
"import NRPy_param_funcs as par # NRPy+: Parameter interface\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",
"\n",
"thismodule = \"GiRaFFE_NRPy_A2B\"\n",
"\n",
"CoordSystem = \"Cartesian\"\n",
"\n",
"# Set the finite-differencing order to 2\n",
"par.set_parval_from_str(\"finite_difference::FD_CENTDERIVS_ORDER\", 2)\n",
"\n",
"# Declare the three metric and compute the sqrt of its determinant.\n",
"gammaDD = ixp.register_gridfunctions_for_single_rank2(\"AUXEVOL\",\"gammaDD\",\"sym01\")\n",
"import GRHD.equations as gh\n",
"gh.compute_sqrtgammaDET(gammaDD)\n",
"\n",
"# Import the Levi-Civita symbol and build the corresponding tensor.\n",
"# We already have a handy function to define the Levi-Civita symbol in indexedexp.py\n",
"LeviCivitaUUU = ixp.LeviCivitaTensorUUU_dim3_rank3(gh.sqrtgammaDET)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We will now create the structures we need: register the gridfunctions `AD` and `BU`, declare `AD_dD` for $A_{i,j}$ and zero `BU` for $B^i$. Then, we'll build the standard form of `BU` that we will use: $$B^i = \\epsilon^{ijk} A_{k,j}.$$"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"AD = ixp.register_gridfunctions_for_single_rank1(\"EVOL\",\"AD\")\n",
"BU = ixp.register_gridfunctions_for_single_rank1(\"AUXEVOL\",\"BU\")\n",
"AD_dD = ixp.declarerank2(\"AD_dD\",\"nosym\")\n",
"BU = ixp.zerorank1() # BU is already registered as a gridfunction, but we need to zero its values and declare it in this scope.\n",
"for i in range(3):\n",
" for j in range(3):\n",
" for k in range(3):\n",
" BU[i] += LeviCivitaUUU[i][j][k] * AD_dD[k][j]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='ccode'></a>\n",
"\n",
"## Step 2: Write out C code \\[Back to [top](#toc)\\]\n",
"$$\\label{ccode}$$\n",
"\n",
"<a id='ghostzones'></a>\n",
"\n",
"### Step 2.a: Compute the magnetic field at a point in a ghostzone \\[Back to [top](#toc)\\]\n",
"$$\\label{ghostzones}$$\n",
"\n",
"The rest of this file will consist of two functions. The first, `compute_A2B_in_ghostzones()`, will loop over the region specified by `i0min`, `i0max`, `i1min`, `i1max`, `i2min`, and `i2max`, passed to the function from the driver function. At each point, we consider the derivatives in each of the three directions that we will need to compute the curl; if we are on the edge of the grid, we will shift the stencil we use one point so that it no longer extends to outside our computational domain. Once we have all the derivatives we will need, we simply calculate the magnetic field as we otherwise would.\n",
"\n",
"With the header files created to actually calculate $B^i$ from $A_k$, we will write the C code to direct how those are used. We will start by including needed header files. We will then define `REAL` as type `double` (so that we can easily change it everywhere if we want to use a different precision) and the `IDX` macros that we use to access specific points in the gridfunction arrays. "
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Overwriting GiRaFFE_NRPy/GiRaFFE_Ccode_validation/A2B//driver_AtoB.h\n"
]
}
],
"source": [
"%%writefile $outdir/driver_AtoB.h\n",
"#include <math.h>\n",
"#include <stdio.h>\n",
"#include <stdlib.h>\n",
"\n",
"#ifndef REAL\n",
"#define REAL double\n",
"#include \"NGHOSTS.h\" // A NRPy+-generated file, which is set based on FD_CENTDERIVS_ORDER.\n",
"#include \"../CurviBoundaryConditions/gridfunction_defines.h\"\n",
"#define IDX4S(g,i,j,k) \\\n",
"( (i) + Nxx_plus_2NGHOSTS0 * ( (j) + Nxx_plus_2NGHOSTS1 * ( (k) + Nxx_plus_2NGHOSTS2 * (g) ) ) )\n",
"#endif\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note the order of the `if` statements here; we first check if a point is *not* on an outermost-face because that is most likely. This will save time during execution by reducing the number of times the `if` statements are evaluated."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Overwriting GiRaFFE_NRPy/GiRaFFE_Ccode_validation/A2B//driver_AtoB.h\n"
]
}
],
"source": [
"%%writefile $outdir/driver_AtoB.h\n",
"void compute_A2B_in_ghostzones(const paramstruct *restrict params,REAL *restrict in_gfs,REAL *restrict auxevol_gfs,\n",
" const int i0min,const int i0max,\n",
" const int i1min,const int i1max,\n",
" const int i2min,const int i2max) {\n",
"#include \"../set_Cparameters.h\"\n",
" for(int i2=i2min;i2<i2max;i2++) for(int i1=i1min;i1<i1max;i1++) for(int i0=i0min;i0<i0max;i0++) {\n",
" REAL dx_Ay,dx_Az,dy_Ax,dy_Az,dz_Ax,dz_Ay;\n",
" // Check to see if we're on the +x or -x face. If so, use a downwinded- or upwinded-stencil, respectively.\n",
" // Otherwise, use a centered stencil.\n",
" if (i0 > 0 && i0 < Nxx_plus_2NGHOSTS0-1) {\n",
" dx_Ay = invdx0*(-1.0/2.0*in_gfs[IDX4S(AD1GF, i0-1,i1,i2)] + (1.0/2.0)*in_gfs[IDX4S(AD1GF, i0+1,i1,i2)]);\n",
" dx_Az = invdx0*(-1.0/2.0*in_gfs[IDX4S(AD2GF, i0-1,i1,i2)] + (1.0/2.0)*in_gfs[IDX4S(AD2GF, i0+1,i1,i2)]);\n",
" }\n",
" else if (i0==0) {\n",
" dx_Ay = invdx0*(-3.0/2.0*in_gfs[IDX4S(AD1GF, i0,i1,i2)] + 2*in_gfs[IDX4S(AD1GF, i0+1,i1,i2)] - 1.0/2.0*in_gfs[IDX4S(AD1GF, i0+2,i1,i2)]);\n",
" dx_Az = invdx0*(-3.0/2.0*in_gfs[IDX4S(AD2GF, i0,i1,i2)] + 2*in_gfs[IDX4S(AD2GF, i0+1,i1,i2)] - 1.0/2.0*in_gfs[IDX4S(AD2GF, i0+2,i1,i2)]);\n",
" }\n",
" else {\n",
" dx_Ay = invdx0*((3.0/2.0)*in_gfs[IDX4S(AD1GF, i0,i1,i2)] - 2*in_gfs[IDX4S(AD1GF, i0-1,i1,i2)] + (1.0/2.0)*in_gfs[IDX4S(AD1GF, i0-2,i1,i2)]);\n",
" dx_Az = invdx0*((3.0/2.0)*in_gfs[IDX4S(AD2GF, i0,i1,i2)] - 2*in_gfs[IDX4S(AD2GF, i0-1,i1,i2)] + (1.0/2.0)*in_gfs[IDX4S(AD2GF, i0-2,i1,i2)]);\n",
" }\n",
" // As above, but in the y direction.\n",
" if (i1 > 0 && i1 < Nxx_plus_2NGHOSTS1-1) {\n",
" dy_Ax = invdx1*(-1.0/2.0*in_gfs[IDX4S(AD0GF, i0,i1-1,i2)] + (1.0/2.0)*in_gfs[IDX4S(AD0GF, i0,i1+1,i2)]);\n",
" dy_Az = invdx1*(-1.0/2.0*in_gfs[IDX4S(AD2GF, i0,i1-1,i2)] + (1.0/2.0)*in_gfs[IDX4S(AD2GF, i0,i1+1,i2)]);\n",
" }\n",
" else if (i1==0) {\n",
" dy_Ax = invdx1*(-3.0/2.0*in_gfs[IDX4S(AD0GF, i0,i1,i2)] + 2*in_gfs[IDX4S(AD0GF, i0,i1+1,i2)] - 1.0/2.0*in_gfs[IDX4S(AD0GF, i0,i1+2,i2)]);\n",
" dy_Az = invdx1*(-3.0/2.0*in_gfs[IDX4S(AD2GF, i0,i1,i2)] + 2*in_gfs[IDX4S(AD2GF, i0,i1+1,i2)] - 1.0/2.0*in_gfs[IDX4S(AD2GF, i0,i1+2,i2)]);\n",
" }\n",
" else {\n",
" dy_Ax = invdx1*((3.0/2.0)*in_gfs[IDX4S(AD0GF, i0,i1,i2)] - 2*in_gfs[IDX4S(AD0GF, i0,i1-1,i2)] + (1.0/2.0)*in_gfs[IDX4S(AD0GF, i0,i1-2,i2)]);\n",
" dy_Az = invdx1*((3.0/2.0)*in_gfs[IDX4S(AD2GF, i0,i1,i2)] - 2*in_gfs[IDX4S(AD2GF, i0,i1-1,i2)] + (1.0/2.0)*in_gfs[IDX4S(AD2GF, i0,i1-2,i2)]);\n",
" }\n",
" // As above, but in the z direction.\n",
" if (i2 > 0 && i2 < Nxx_plus_2NGHOSTS2-1) {\n",
" dz_Ax = invdx2*(-1.0/2.0*in_gfs[IDX4S(AD0GF, i0,i1,i2-1)] + (1.0/2.0)*in_gfs[IDX4S(AD0GF, i0,i1,i2+1)]);\n",
" dz_Ay = invdx2*(-1.0/2.0*in_gfs[IDX4S(AD1GF, i0,i1,i2-1)] + (1.0/2.0)*in_gfs[IDX4S(AD1GF, i0,i1,i2+1)]);\n",
" }\n",
" else if (i2==0) {\n",
" dz_Ax = invdx2*(-3.0/2.0*in_gfs[IDX4S(AD0GF, i0,i1,i2)] + 2*in_gfs[IDX4S(AD0GF, i0,i1,i2+1)] - 1.0/2.0*in_gfs[IDX4S(AD0GF, i0,i1,i2+2)]);\n",
" dz_Ay = invdx2*(-3.0/2.0*in_gfs[IDX4S(AD1GF, i0,i1,i2)] + 2*in_gfs[IDX4S(AD1GF, i0,i1,i2+1)] - 1.0/2.0*in_gfs[IDX4S(AD1GF, i0,i1,i2+2)]);\n",
" }\n",
" else {\n",
" dz_Ax = invdx2*((3.0/2.0)*in_gfs[IDX4S(AD0GF, i0,i1,i2)] - 2*in_gfs[IDX4S(AD0GF, i0,i1,i2-1)] + (1.0/2.0)*in_gfs[IDX4S(AD0GF, i0,i1,i2-2)]);\n",
" dz_Ay = invdx2*((3.0/2.0)*in_gfs[IDX4S(AD1GF, i0,i1,i2)] - 2*in_gfs[IDX4S(AD1GF, i0,i1,i2-1)] + (1.0/2.0)*in_gfs[IDX4S(AD1GF, i0,i1,i2-2)]);\n",
" }\n",
" // Compute the magnetic field in the normal way, using the previously calculated derivatives.\n",
" const double gammaDD00 = auxevol_gfs[IDX4S(GAMMADD00GF, i0,i1,i2)];\n",
" const double gammaDD01 = auxevol_gfs[IDX4S(GAMMADD01GF, i0,i1,i2)];\n",
" const double gammaDD02 = auxevol_gfs[IDX4S(GAMMADD02GF, i0,i1,i2)];\n",
" const double gammaDD11 = auxevol_gfs[IDX4S(GAMMADD11GF, i0,i1,i2)];\n",
" const double gammaDD12 = auxevol_gfs[IDX4S(GAMMADD12GF, i0,i1,i2)];\n",
" const double gammaDD22 = auxevol_gfs[IDX4S(GAMMADD22GF, i0,i1,i2)];\n",
" /*\n",
" * NRPy+ Finite Difference Code Generation, Step 2 of 1: Evaluate SymPy expressions and write to main memory:\n",
" */\n",
" const double invsqrtg = 1.0/sqrt(gammaDD00*gammaDD11*gammaDD22\n",
" - gammaDD00*gammaDD12*gammaDD12\n",
" + 2*gammaDD01*gammaDD02*gammaDD12\n",
" - gammaDD11*gammaDD02*gammaDD02\n",
" - gammaDD22*gammaDD01*gammaDD01);\n",
" auxevol_gfs[IDX4S(BU0GF, i0,i1,i2)] = (dy_Az-dz_Ay)*invsqrtg;\n",
" auxevol_gfs[IDX4S(BU1GF, i0,i1,i2)] = (dz_Ax-dx_Az)*invsqrtg;\n",
" auxevol_gfs[IDX4S(BU2GF, i0,i1,i2)] = (dx_Ay-dy_Ax)*invsqrtg;\n",
" }\n",
"}\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='driver'></a>\n",
"\n",
"### Step 2.b: Multiple-order algorithm for interior \\[Back to [top](#toc)\\]\n",
"$$\\label{driver}$$\n",
"\n",
"To help avoid Gibbs' phenomenon, we will implement an algorithm to lower the finite-differencing order near shocks. This will be done by calculating the magnetic field at multiple orders and using the lower order if the different orders disagree. What we write here will take the place of a NRPy+-generated 'body', as the parameter of the NRPy+ function `outCfunction()`. It will act on each point and each component individually. This will be much easier with some auxiliary functions, which we will write first. One will be a basic relative error function, another will compute a component of the magnetic field based on the input vector potential, and the last will return the best-choice order. "
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to GiRaFFE_NRPy/GiRaFFE_Ccode_validation/A2B//driver_AtoB.h\n"
]
}
],
"source": [
"%%writefile -a $outdir/driver_AtoB.h\n",
"REAL relative_error(REAL a, REAL b) {\n",
" if((a+b)!=0.0) {\n",
" return 2.0*fabs(a-b)/fabs(a+b);\n",
" }\n",
" else {\n",
" return 0.0;\n",
" }\n",
"}\n",
"\n",
"#define M2 0\n",
"#define M1 1\n",
"#define P0 2\n",
"#define P1 3\n",
"#define P2 4\n",
"#define CN4 0\n",
"#define CN2 1\n",
"#define UP2 2\n",
"#define DN2 3\n",
"#define UP1 4\n",
"#define DN1 5\n",
"void compute_Bx_pointwise(REAL *Bx, const REAL invdy, const REAL *Ay, const REAL invdz, const REAL *Az) {\n",
" REAL dz_Ay,dy_Az;\n",
" dz_Ay = invdz*((Ay[P1]-Ay[M1])*2.0/3.0 - (Ay[P2]-Ay[M2])/12.0);\n",
" dy_Az = invdy*((Az[P1]-Az[M1])*2.0/3.0 - (Az[P2]-Az[M2])/12.0);\n",
" Bx[CN4] = dy_Az - dz_Ay;\n",
"\n",
" dz_Ay = invdz*(Ay[P1]-Ay[M1])/2.0;\n",
" dy_Az = invdy*(Az[P1]-Az[M1])/2.0;\n",
" Bx[CN2] = dy_Az - dz_Ay;\n",
"\n",
" dz_Ay = invdz*(-1.5*Ay[P0]+2.0*Ay[P1]-0.5*Ay[P2]);\n",
" dy_Az = invdy*(-1.5*Az[P0]+2.0*Az[P1]-0.5*Az[P2]);\n",
" Bx[UP2] = dy_Az - dz_Ay;\n",
"\n",
" dz_Ay = invdz*(1.5*Ay[P0]-2.0*Ay[M1]+0.5*Ay[M2]);\n",
" dy_Az = invdy*(1.5*Az[P0]-2.0*Az[M1]+0.5*Az[M2]);\n",
" Bx[DN2] = dy_Az - dz_Ay;\n",
"\n",
" dz_Ay = invdz*(Ay[P1]-Ay[P0]);\n",
" dy_Az = invdy*(Az[P1]-Az[P0]);\n",
" Bx[UP1] = dy_Az - dz_Ay;\n",
"\n",
" dz_Ay = invdz*(Ay[P0]-Ay[M1]);\n",
" dy_Az = invdy*(Az[P0]-Az[M1]);\n",
" Bx[DN1] = dy_Az - dz_Ay;\n",
"}\n",
"\n",
"#define TOLERANCE_A2B 1.0e-4\n",
"REAL find_accepted_Bx_order(REAL *Bx) {\n",
" REAL accepted_val = Bx[CN4];\n",
" REAL Rel_error_o2_vs_o4 = relative_error(Bx[CN2],Bx[CN4]);\n",
" REAL Rel_error_oCN2_vs_oDN2 = relative_error(Bx[CN2],Bx[DN2]);\n",
" REAL Rel_error_oCN2_vs_oUP2 = relative_error(Bx[CN2],Bx[UP2]);\n",
"\n",
" if(Rel_error_o2_vs_o4 > TOLERANCE_A2B) {\n",
" accepted_val = Bx[CN2];\n",
" if(Rel_error_o2_vs_o4 > Rel_error_oCN2_vs_oDN2 || Rel_error_o2_vs_o4 > Rel_error_oCN2_vs_oUP2) {\n",
" // Should we use AND or OR in if statement?\n",
" if(relative_error(Bx[UP2],Bx[UP1]) < relative_error(Bx[DN2],Bx[DN1])) {\n",
" accepted_val = Bx[UP2];\n",
" }\n",
" else {\n",
" accepted_val = Bx[DN2];\n",
" }\n",
" }\n",
" }\n",
" return accepted_val;\n",
"}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The algorithm will start by reading in all the quantities it will need and storing them as temporary variables. This includes all three components of the vector potential in two different directions each, reflecting the derivatives we will need to take. Five points are needed in each direction in order to fill the templates for fourth- and second-order centered-differencing and first-order forward- and backward-differencing. We calculate the inverse of the square root of the metric determinant and then compute all three components of the magnetic field at all orders. Then, if the relative differences between orders are found to be high, the order is dropped. The final result is written to memory.\n",
"\n",
"***This is currently not used, it does not reliably improve shock handling***"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"# order_lowering_body = \"\"\"REAL AD0_1[5],AD0_2[5],AD1_2[5],AD1_0[5],AD2_0[5],AD2_1[5];\n",
"# const double gammaDD00 = auxevol_gfs[IDX4S(GAMMADD00GF, i0,i1,i2)];\n",
"# const double gammaDD01 = auxevol_gfs[IDX4S(GAMMADD01GF, i0,i1,i2)];\n",
"# const double gammaDD02 = auxevol_gfs[IDX4S(GAMMADD02GF, i0,i1,i2)];\n",
"# const double gammaDD11 = auxevol_gfs[IDX4S(GAMMADD11GF, i0,i1,i2)];\n",
"# const double gammaDD12 = auxevol_gfs[IDX4S(GAMMADD12GF, i0,i1,i2)];\n",
"# const double gammaDD22 = auxevol_gfs[IDX4S(GAMMADD22GF, i0,i1,i2)];\n",
"# AD0_2[M2] = in_gfs[IDX4S(AD0GF, i0,i1,i2-2)];\n",
"# AD0_2[M1] = in_gfs[IDX4S(AD0GF, i0,i1,i2-1)];\n",
"# AD0_1[M2] = in_gfs[IDX4S(AD0GF, i0,i1-2,i2)];\n",
"# AD0_1[M1] = in_gfs[IDX4S(AD0GF, i0,i1-1,i2)];\n",
"# AD0_1[P0] = AD0_2[P0] = in_gfs[IDX4S(AD0GF, i0,i1,i2)];\n",
"# AD0_1[P1] = in_gfs[IDX4S(AD0GF, i0,i1+1,i2)];\n",
"# AD0_1[P2] = in_gfs[IDX4S(AD0GF, i0,i1+2,i2)];\n",
"# AD0_2[P1] = in_gfs[IDX4S(AD0GF, i0,i1,i2+1)];\n",
"# AD0_2[P2] = in_gfs[IDX4S(AD0GF, i0,i1,i2+2)];\n",
"# AD1_2[M2] = in_gfs[IDX4S(AD1GF, i0,i1,i2-2)];\n",
"# AD1_2[M1] = in_gfs[IDX4S(AD1GF, i0,i1,i2-1)];\n",
"# AD1_0[M2] = in_gfs[IDX4S(AD1GF, i0-2,i1,i2)];\n",
"# AD1_0[M1] = in_gfs[IDX4S(AD1GF, i0-1,i1,i2)];\n",
"# AD1_2[P0] = AD1_0[P0] = in_gfs[IDX4S(AD1GF, i0,i1,i2)];\n",
"# AD1_0[P1] = in_gfs[IDX4S(AD1GF, i0+1,i1,i2)];\n",
"# AD1_0[P2] = in_gfs[IDX4S(AD1GF, i0+2,i1,i2)];\n",
"# AD1_2[P1] = in_gfs[IDX4S(AD1GF, i0,i1,i2+1)];\n",
"# AD1_2[P2] = in_gfs[IDX4S(AD1GF, i0,i1,i2+2)];\n",
"# AD2_1[M2] = in_gfs[IDX4S(AD2GF, i0,i1-2,i2)];\n",
"# AD2_1[M1] = in_gfs[IDX4S(AD2GF, i0,i1-1,i2)];\n",
"# AD2_0[M2] = in_gfs[IDX4S(AD2GF, i0-2,i1,i2)];\n",
"# AD2_0[M1] = in_gfs[IDX4S(AD2GF, i0-1,i1,i2)];\n",
"# AD2_0[P0] = AD2_1[P0] = in_gfs[IDX4S(AD2GF, i0,i1,i2)];\n",
"# AD2_0[P1] = in_gfs[IDX4S(AD2GF, i0+1,i1,i2)];\n",
"# AD2_0[P2] = in_gfs[IDX4S(AD2GF, i0+2,i1,i2)];\n",
"# AD2_1[P1] = in_gfs[IDX4S(AD2GF, i0,i1+1,i2)];\n",
"# AD2_1[P2] = in_gfs[IDX4S(AD2GF, i0,i1+2,i2)];\n",
"# const double invsqrtg = 1.0/sqrt(gammaDD00*gammaDD11*gammaDD22\n",
"# - gammaDD00*gammaDD12*gammaDD12\n",
"# + 2*gammaDD01*gammaDD02*gammaDD12\n",
"# - gammaDD11*gammaDD02*gammaDD02\n",
"# - gammaDD22*gammaDD01*gammaDD01);\n",
"\n",
"# REAL BU0[4],BU1[4],BU2[4];\n",
"# compute_Bx_pointwise(BU0,invdx2,AD1_2,invdx1,AD2_1);\n",
"# compute_Bx_pointwise(BU1,invdx0,AD2_0,invdx2,AD0_2);\n",
"# compute_Bx_pointwise(BU2,invdx1,AD0_1,invdx0,AD1_0);\n",
"\n",
"# auxevol_gfs[IDX4S(BU0GF, i0,i1,i2)] = find_accepted_Bx_order(BU0)*invsqrtg;\n",
"# auxevol_gfs[IDX4S(BU1GF, i0,i1,i2)] = find_accepted_Bx_order(BU1)*invsqrtg;\n",
"# auxevol_gfs[IDX4S(BU2GF, i0,i1,i2)] = find_accepted_Bx_order(BU2)*invsqrtg;\n",
"# \"\"\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='driver'></a>\n",
"\n",
"### Step 2.c: Looping over all ghostzones \\[Back to [top](#toc)\\]\n",
"$$\\label{driver}$$\n",
"\n",
"This function is responsible for driving the logic needed to compute $B^i$ from $A_i$ whenever it is needed in the main code (typically done along with applying BCs after a timestep in the main code). We first loop over the interior and calculate $B^i$ there at the FD order currently set in NRPy+. We then set `imin` and `imax` such that they define the interior of the grid. \n",
"\n",
"With the Interior handled, we call `compute_A2B_in_ghostzones` several times. Note that `imin` and `imax` are passed such that each call loops over one face in one ghostzone (e.g., the innermost ghostzone on the $+x$ face). On the first pass through the loop, we calculate $B^i$ in the innermost ghostzone. We also decrement `imin` and increment `imax` on each pass through the loop. This serves two purposes: the next time through the loop, `compute_A2B_in_ghostzones()` will act on the next-outer ghost zone. on the next-outer ghost zone, and the size of each face will increase, eventually allowing us to act on each point. "
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"desc=\"Compute the magnetic field from the vector potential everywhere, including ghostzones\"\n",
"# body = order_lowering_body,\n",
"name=\"driver_A_to_B\"\n",
"driver_Ccode = outCfunction(\n",
" outfile = \"returnstring\", desc=desc, name=name,\n",
" params = \"const paramstruct *restrict params,REAL *restrict in_gfs,REAL *restrict auxevol_gfs\",\n",
" body = fin.FD_outputC(\"returnstring\",[lhrh(lhs=gri.gfaccess(\"out_gfs\",\"BU0\"),rhs=BU[0]),\\\n",
" lhrh(lhs=gri.gfaccess(\"out_gfs\",\"BU1\"),rhs=BU[1]),\\\n",
" lhrh(lhs=gri.gfaccess(\"out_gfs\",\"BU2\"),rhs=BU[2])]),\n",
" postloop = \"\"\"\n",
" int imin[3] = { NGHOSTS_A2B, NGHOSTS_A2B, NGHOSTS_A2B };\n",
" int imax[3] = { NGHOSTS+Nxx0, NGHOSTS+Nxx1, NGHOSTS+Nxx2 };\n",
" // Now, we loop over the ghostzones to calculate the magnetic field there.\n",
" for(int which_gz = 0; which_gz < NGHOSTS_A2B; which_gz++) {\n",
" // After updating each face, adjust imin[] and imax[]\n",
" // to reflect the newly-updated face extents.\n",
" compute_A2B_in_ghostzones(params,in_gfs,auxevol_gfs,imin[0]-1,imin[0], imin[1],imax[1], imin[2],imax[2]); imin[0]--;\n",
" compute_A2B_in_ghostzones(params,in_gfs,auxevol_gfs,imax[0],imax[0]+1, imin[1],imax[1], imin[2],imax[2]); imax[0]++;\n",
"\n",
" compute_A2B_in_ghostzones(params,in_gfs,auxevol_gfs,imin[0],imax[0], imin[1]-1,imin[1], imin[2],imax[2]); imin[1]--;\n",
" compute_A2B_in_ghostzones(params,in_gfs,auxevol_gfs,imin[0],imax[0], imax[1],imax[1]+1, imin[2],imax[2]); imax[1]++;\n",
"\n",
" compute_A2B_in_ghostzones(params,in_gfs,auxevol_gfs,imin[0],imax[0], imin[1],imax[1], imin[2]-1,imin[2]); imin[2]--;\n",
" compute_A2B_in_ghostzones(params,in_gfs,auxevol_gfs,imin[0],imax[0], imin[1],imax[1], imax[2],imax[2]+1); imax[2]++;\n",
" }\n",
"\"\"\",\n",
" loopopts=\"InteriorPoints\",\n",
" rel_path_to_Cparams=os.path.join(\"../\")).replace(\"= NGHOSTS\",\"= NGHOSTS_A2B\").replace(\"NGHOSTS+Nxx0\",\"Nxx_plus_2NGHOSTS0-NGHOSTS_A2B\").replace(\"NGHOSTS+Nxx1\",\"Nxx_plus_2NGHOSTS1-NGHOSTS_A2B\").replace(\"NGHOSTS+Nxx2\",\"Nxx_plus_2NGHOSTS2-NGHOSTS_A2B\")\n",
"\n",
"with open(os.path.join(outdir,\"driver_AtoB.h\"),\"a\") as file:\n",
" file.write(driver_Ccode)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='code_validation'></a>\n",
"\n",
"# Step 3: Code Validation against original C code \\[Back to [top](#toc)\\]\n",
"$$\\label{code_validation}$$\n",
"\n",
"To validate the code in this tutorial we check for agreement between the files\n",
"\n",
"1. that were written in this tutorial and\n",
"1. those that are stored in `GiRaFFE_NRPy/GiRaFFE_Ccode_library` or generated by `GiRaFFE_NRPy_A2B.py`\n"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Printing difference between original C code and this code...\n",
"Checking file driver_AtoB.h\n",
"No difference. TEST PASSED!\n"
]
}
],
"source": [
"# Define the directory that we wish to validate against:\n",
"valdir = \"GiRaFFE_NRPy/GiRaFFE_Ccode_library/A2B/\"\n",
"\n",
"import GiRaFFE_NRPy.GiRaFFE_NRPy_A2B as A2B\n",
"A2B.GiRaFFE_NRPy_A2B(valdir,gammaDD,AD,BU)\n",
"\n",
"import difflib\n",
"import sys\n",
"\n",
"print(\"Printing difference between original C code and this code...\")\n",
"# Open the files to compare\n",
"file = \"driver_AtoB.h\"\n",
"\n",
"print(\"Checking file \" + file)\n",
"with open(os.path.join(valdir,file)) as file1, open(os.path.join(outdir,file)) as file2:\n",
" # Read the lines of each file\n",
" file1_lines = file1.readlines()\n",
" file2_lines = file2.readlines()\n",
" num_diffs = 0\n",
" for line in difflib.unified_diff(file1_lines, file2_lines, fromfile=os.path.join(valdir+file), tofile=os.path.join(outdir+file)):\n",
" sys.stdout.writelines(line)\n",
" num_diffs = num_diffs + 1\n",
" if num_diffs == 0:\n",
" print(\"No difference. TEST PASSED!\")\n",
" else:\n",
" print(\"ERROR: Disagreement found with .py file. See differences above.\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='latex_pdf_output'></a>\n",
"\n",
"# Step 4: 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-GiRaFFE_NRPy_C_code_library-A2B](TTutorial-GiRaFFE_NRPy_C_code_library-A2B.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": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Created Tutorial-GiRaFFE_NRPy-A2B.tex, and compiled LaTeX file to PDF file\n",
" Tutorial-GiRaFFE_NRPy-A2B.pdf\n"
]
}
],
"source": [
"import cmdline_helper as cmd # NRPy+: Multi-platform Python command-line interface\n",
"cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-GiRaFFE_NRPy-A2B\",location_of_template_file=os.path.join(\"..\"))"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"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.8.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}