GiRaFFE_NRPy
: Induction Equation¶GiRaFFE
.¶Notebook Status: Validated
Validation Notes: This module validates the routines in Tutorial-GiRaFFE_NRPy-Afield_flux_handwritten.
This notebook validates our algorithm to compute the flux of $\epsilon_{ijk} v^j B^k$ through cell faces, which contributes to the right-hand side of the evolution equation for $A_i$, for use in GiRaFFE_NRPy
. Because the original GiRaFFE
used staggered grids and we do not, we can not trivially do a direct comparison to the old code. Instead, we will compare the numerical results with the expected analytic results.
It is, in general, good coding practice to unit test functions individually to verify that they produce the expected and intended output. Here, we expect our functions to produce the correct cross product between the velocity and magnetic field in an arbitrary spacetime. To that end, we will choose functions with relatively simple analytic forms.
TERRENCE SAYS: wording is confusing here
In this test, we will generate analytic forms for the magnetic field and three-velocity. We will need to do this in a $7 \times 7 \times 7$ cube in order to run the PPM routine on the data to generate the inputs compute the non-gauge terms of the RHS of the induction equation, since it would be impractical to usefully spoof the left- and right-face values for the HLLE solver. We care about this here, because we are comparing against an analytic expression and not the old code.
When this notebook is run, the difference between the approximate and exact right-hand sides will be output to text files that can be found in the same directory as this notebook. These will be read in in Step 3, and used there to confirm second order convergence of the algorithm.
This notebook is organized as follows
GiRaFFE_NRPy
FilesAfield_flux_unit_test.c
: The Main C CodeWe'll start by appending the relevant paths to sys.path
so that we can access sympy modules in other places. Then, we'll import NRPy+ core functionality and set up a directory in which to carry out our test. We must also set the desired finite differencing order.
import os, sys, shutil # Standard Python modules for multiplatform OS-level functions
# First, we'll add the parent directory to the list of directories Python will check for modules.
nrpy_dir_path = os.path.join("..")
if nrpy_dir_path not in sys.path:
sys.path.append(nrpy_dir_path)
nrpy_dir_path = os.path.join("..","..")
if nrpy_dir_path not in sys.path:
sys.path.append(nrpy_dir_path)
from outputC import outCfunction, outputC, lhrh # NRPy+: Core C code output module
import sympy as sp # SymPy: The Python computer algebra package upon which NRPy+ depends
import finite_difference as fin # NRPy+: Finite difference C code generation module
import NRPy_param_funcs as par # NRPy+: Parameter interface
import grid as gri # NRPy+: Functions having to do with numerical grids
import indexedexp as ixp # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support
import reference_metric as rfm # NRPy+: Reference metric support
import cmdline_helper as cmd # NRPy+: Multi-platform Python command-line interface
Ccodesdir = "Start-to-Finish-UnitTests/Afield_flux_UnitTest/"
# First remove C code output directory if it exists
# Courtesy https://stackoverflow.com/questions/303200/how-do-i-remove-delete-a-folder-that-is-not-empty
shutil.rmtree(Ccodesdir, ignore_errors=True)
# Then create a fresh directory
cmd.mkdir(Ccodesdir)
outdir = os.path.join(Ccodesdir,"output/")
cmd.mkdir(outdir)
thismodule = "Start_to_Finish_UnitTest-GiRaFFE_NRPy-Afield_flux"
Use_Shock_Data = True
First, we'll set some arbitrary functions for the three-metric $\gamma_{ij}$, the shift vector $\beta^i$, and the lapse function $\alpha$. We will need these to compute the characteristic speeds in each direction. In addition, we will need the metric determinant to exactly compute the electric field at the test point.
We will use this basic form for $\gamma_{ij}$. We will thus write a function with the following arbitrary equations.
\begin{align} \gamma_{xx} &= ax^3 + by^3 + cz^3 + dy^2 + ez^2 + 1 \\ \gamma_{yy} &= gx^3 + hy^3 + lz^3 + mx^2 + nz^2 + 1 \\ \gamma_{zz} &= px^3 + qy^3 + rz^3 + sx^2 + ty^2 + 1. \\ \gamma_{xy} &= \frac{1}{10} \exp\left(-\left((x-b)^2+(y-c)^2+(z-d)^2\right)\right) \\ \gamma_{xz} &= \frac{1}{10} \exp\left(-\left((x-g)^2+(y-h)^2+(z-l)^2\right)\right) \\ \gamma_{yz} &= \frac{1}{10} \exp\left(-\left((x-n)^2+(y-o)^2+(z-p)^2\right)\right), \\ \end{align}a,b,c,d,e,f,g,h,l,m,n,o,p,q,r,s,t = par.Cparameters("REAL",thismodule,["a","b","c","d","e","f","g","h","l","m","n","o","p","q","r","s","t"],1e300)
# Note that the above default value allows us to set these directly in the C code
M_PI = par.Cparameters("#define",thismodule,["M_PI"], "")
par.set_parval_from_str("reference_metric::CoordSystem","Cartesian")
rfm.reference_metric()
x = rfm.xx_to_Cart[0]
y = rfm.xx_to_Cart[1]
z = rfm.xx_to_Cart[2]
gammaDD = ixp.register_gridfunctions_for_single_rank2("AUXEVOL","gammaDD","sym01",DIM=3)
betaU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","betaU",DIM=3)
alpha = gri.register_gridfunctions("AUXEVOL","alpha")
gammaDD[0][0] = a*x**3 + b*y**3 + c*z**3 + d*y**2 + e*z**2 + sp.sympify(1)
gammaDD[1][1] = g*x**3 + h*y**3 + l*z**3 + m*x**2 + n*z**2 + sp.sympify(1)
gammaDD[2][2] = p*x**3 + q*y**3 + r*z**3 + s*x**2 + t*y**2 + sp.sympify(1)
gammaDD[0][1] = sp.Rational(1,10) * sp.exp(-((x-b)**2 + (y-c)**2 + (z-d)**2))
gammaDD[0][2] = sp.Rational(1,10) * sp.exp(-((x-g)**2 + (y-h)**2 + (z-l)**2))
gammaDD[1][2] = sp.Rational(1,10) * sp.exp(-((x-n)**2 + (y-o)**2 + (z-p)**2))
betaU[0] = sp.sympify(0)#(sp.sympify(2)/M_PI) * sp.atan(a*x + b*y + c*z)
betaU[1] = sp.sympify(0)#(sp.sympify(2)/M_PI) * sp.atan(b*x + c*y + a*z)
betaU[2] = sp.sympify(0)#(sp.sympify(2)/M_PI) * sp.atan(c*x + a*y + b*z)
alpha = sp.sympify(1)#sp.sympify(1) - sp.sympify(1) / (sp.sympify(2) + x**2 + y**2 + z**2)
metric_gfs_to_print = [
lhrh(lhs=gri.gfaccess("aux_gfs","gammaDD00"),rhs=gammaDD[0][0]),
lhrh(lhs=gri.gfaccess("aux_gfs","gammaDD01"),rhs=gammaDD[0][1]),
lhrh(lhs=gri.gfaccess("aux_gfs","gammaDD02"),rhs=gammaDD[0][2]),
lhrh(lhs=gri.gfaccess("aux_gfs","gammaDD11"),rhs=gammaDD[1][1]),
lhrh(lhs=gri.gfaccess("aux_gfs","gammaDD12"),rhs=gammaDD[1][2]),
lhrh(lhs=gri.gfaccess("aux_gfs","gammaDD22"),rhs=gammaDD[2][2]),
lhrh(lhs=gri.gfaccess("aux_gfs","betaU0"),rhs=betaU[0]),
lhrh(lhs=gri.gfaccess("aux_gfs","betaU1"),rhs=betaU[1]),
lhrh(lhs=gri.gfaccess("aux_gfs","betaU2"),rhs=betaU[2]),
lhrh(lhs=gri.gfaccess("aux_gfs","alpha"),rhs=alpha),
]
desc = "Calculate the metric gridfunctions"
name = "calculate_metric_gfs"
outCfunction(
outfile = os.path.join(Ccodesdir,name+".h"), desc=desc, name=name,
params ="const paramstruct *restrict params,REAL *restrict xx[3],REAL *restrict auxevol_gfs",
body = fin.FD_outputC("returnstring",metric_gfs_to_print,params="outCverbose=False"),
loopopts="AllPoints,Read_xxs")
Output C function calculate_metric_gfs() to file Start-to-Finish-UnitTests/Afield_flux_UnitTest/calculate_metric_gfs.h
Here, we'll generate some functions for the velocity. Let's choose arctangents, since those have asymptotes that can be easily manipulated to prevent accidentally setting superluminal speeds.
\begin{align}
\bar{v}^x &= \frac{1}{5\pi} \arctan(a(x-k) + by + cz) \\
\bar{v}^y &= \frac{1}{5\pi} \arctan(bx + c(y-k) + az) \\
\bar{v}^z &= \frac{1}{5\pi} \arctan(cx + ay + b(z-k)) \\
\end{align}
We offset by $k$ so we have smooth data in the region of interest when Use_Shock_Data
is false.
If we want to add a jump at the origin, we can simply add $\max(0,5x)$ to the argument of the arctangent. This will add a shock in the $x$-direction. The maximum will be described without the use of if statements using the Min_Max_and_Piecewise_Expressions
module.
import Min_Max_and_Piecewise_Expressions as noif
offset = sp.sympify(5)
args = ixp.zerorank1()
args[0] = a*(x-offset) + b*y + c*z
args[1] = b*x + c*(y-offset) + a*z
args[2] = c*x + a*y + b*(z-offset)
if Use_Shock_Data:
for i in range(3):
args[i] += noif.max_noif(0,5*x)
ValenciavU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","ValenciavU")
for i in range(3):
ValenciavU[i] = (sp.Rational(1,5)/M_PI)*sp.atan(args[i])
We'll also need some functions for the magnetic field. Exponentials sound fun. \begin{align} B^x &= \exp\left(\frac{1}{10}(ey+fz)\right) \\ B^y &= \exp\left(\frac{1}{10}(fz+dx)\right) \\ B^z &= \exp\left(\frac{1}{10}(dx+ey)\right) \\ \end{align} In this case, we'll add $\max{0,5x}$ to the field to add the jump.
BU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","BU")
BU[0] = sp.exp(sp.Rational(1,10)*(e*y+f*z))
BU[1] = sp.exp(sp.Rational(1,10)*(f*z+d*x))
BU[2] = sp.exp(sp.Rational(1,10)*(d*x+e*y))
if Use_Shock_Data:
for i in range(3):
BU[i] += noif.max_noif(0,5*x)
BU_to_print = [
lhrh(lhs=gri.gfaccess("out_gfs","BU0"),rhs=BU[0]),
lhrh(lhs=gri.gfaccess("out_gfs","BU1"),rhs=BU[1]),
lhrh(lhs=gri.gfaccess("out_gfs","BU2"),rhs=BU[2]),
]
desc = "Calculate sample magnetic field data"
name = "calculate_BU"
outCfunction(
outfile = os.path.join(Ccodesdir,name+".h"), desc=desc, name=name,
params ="const paramstruct *params,REAL *xx[3],REAL *auxevol_gfs",
body = fin.FD_outputC("returnstring",BU_to_print,params="outCverbose=False"),
loopopts="AllPoints,Read_xxs")
ValenciavU_to_print = [
lhrh(lhs=gri.gfaccess("out_gfs","ValenciavU0"),rhs=ValenciavU[0]),
lhrh(lhs=gri.gfaccess("out_gfs","ValenciavU1"),rhs=ValenciavU[1]),
lhrh(lhs=gri.gfaccess("out_gfs","ValenciavU2"),rhs=ValenciavU[2]),
]
desc = "Calculate sample velocity data"
name = "calculate_ValenciavU"
outCfunction(
outfile = os.path.join(Ccodesdir,name+".h"), desc=desc, name=name,
params ="const paramstruct *params,REAL *xx[3],REAL *auxevol_gfs",
body = fin.FD_outputC("returnstring",ValenciavU_to_print,params="outCverbose=False"),
loopopts="AllPoints,Read_xxs")
Output C function calculate_BU() to file Start-to-Finish-UnitTests/Afield_flux_UnitTest/calculate_BU.h Output C function calculate_ValenciavU() to file Start-to-Finish-UnitTests/Afield_flux_UnitTest/calculate_ValenciavU.h
# Step 3.d.i: Generate declare_Cparameters_struct.h, set_Cparameters_default.h, and set_Cparameters[-SIMD].h
NGHOSTS = par.Cparameters("int",thismodule,["NGHOSTS"], 3)
# Step 3.d.ii: Set free_parameters.h
with open(os.path.join(Ccodesdir,"free_parameters.h"),"w") as file:
file.write("""
// Override parameter defaults with values based on command line arguments and NGHOSTS.
// We'll use this grid. It has one point and one ghost zone.
params.NGHOSTS = 3;
params.Nxx0 = atoi(argv[1]);
params.Nxx1 = atoi(argv[2]);
params.Nxx2 = atoi(argv[3]);
params.Nxx_plus_2NGHOSTS0 = params.Nxx0 + 2*params.NGHOSTS;
params.Nxx_plus_2NGHOSTS1 = params.Nxx1 + 2*params.NGHOSTS;
params.Nxx_plus_2NGHOSTS2 = params.Nxx2 + 2*params.NGHOSTS;
// Step 0d: Set up space and time coordinates
// Step 0d.i: Declare \Delta x^i=dxx{0,1,2} and invdxx{0,1,2}, as well as xxmin[3] and xxmax[3]:
const REAL xxmin[3] = {-0.1,-0.1,-0.1};
const REAL xxmax[3] = { 0.1, 0.1, 0.1};
params.dxx0 = (xxmax[0] - xxmin[0]) / ((REAL)params.Nxx0);
params.dxx1 = (xxmax[1] - xxmin[1]) / ((REAL)params.Nxx1);
params.dxx2 = (xxmax[2] - xxmin[2]) / ((REAL)params.Nxx2);
//printf("dxx0,dxx1,dxx2 = %.5e,%.5e,%.5e\\n",params.dxx0,params.dxx1,params.dxx2);
params.invdx0 = 1.0 / params.dxx0;
params.invdx1 = 1.0 / params.dxx1;
params.invdx2 = 1.0 / params.dxx2;
\n""")
# Generates declare_Cparameters_struct.h, set_Cparameters_default.h, and set_Cparameters[-SIMD].h
par.generate_Cparameters_Ccodes(os.path.join(Ccodesdir))
GiRaFFE_NRPy
Files [Back to top]¶Here, we generate the functions we want to test by calling the function found here and documented in this tutorial.
subdir = "PPM"
cmd.mkdir(os.path.join(Ccodesdir,subdir))
import GiRaFFE_NRPy.GiRaFFE_NRPy_PPM as PPM
PPM.GiRaFFE_NRPy_PPM(os.path.join(Ccodesdir,subdir))
import GiRaFFE_NRPy.GiRaFFE_NRPy_Afield_flux_handwritten as Af
# We will pass values of the gridfunction on the cell faces into the function. This requires us
# to declare them as C parameters in NRPy+. We will denote this with the _face infix/suffix.
alpha_face = gri.register_gridfunctions("AUXEVOL","alpha_face")
gamma_faceDD = ixp.register_gridfunctions_for_single_rank2("AUXEVOL","gamma_faceDD","sym01")
beta_faceU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","beta_faceU")
# We'll need some more gridfunctions, now, to represent the reconstructions of BU and ValenciavU
# on the right and left faces
Valenciav_rU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","Valenciav_rU",DIM=3)
B_rU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","B_rU",DIM=3)
Valenciav_lU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","Valenciav_lU",DIM=3)
B_lU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","B_lU",DIM=3)
AD = ixp.register_gridfunctions_for_single_rank1("EVOL","AD",DIM=3)
subdir = "RHSs"
cmd.mkdir(os.path.join(Ccodesdir,subdir))
Af.GiRaFFE_NRPy_Afield_flux(os.path.join(Ccodesdir,subdir))
subdir = "FCVAL"
cmd.mkdir(os.path.join(Ccodesdir,subdir))
# We also need to interpolate the metric values to cell faces, which is easily done with this module:
import GiRaFFE_NRPy.GiRaFFE_NRPy_Metric_Face_Values as FCVAL
FCVAL.GiRaFFE_NRPy_FCVAL(os.path.join(Ccodesdir,subdir))
import GRHD.equations as gh
gh.compute_sqrtgammaDET(gammaDD)
LeviCivitaDDD = ixp.LeviCivitaTensorDDD_dim3_rank3(gh.sqrtgammaDET)
driftvU = ixp.zerorank1()
for i in range(3):
# TODO: don't use _face metric gridfunctions once we're in curved space!
driftvU[i] = alpha*ValenciavU[i]-betaU[i]
A_rhsD = ixp.zerorank1()
for i in range(3):
for j in range(3):
for k in range(3):
A_rhsD[i] += LeviCivitaDDD[i][j][k]*driftvU[j]*BU[k]
A_rhsD_to_print = [
lhrh(lhs=gri.gfaccess("rhs_gfs","AD0"),rhs=A_rhsD[0]),
lhrh(lhs=gri.gfaccess("rhs_gfs","AD1"),rhs=A_rhsD[1]),
lhrh(lhs=gri.gfaccess("rhs_gfs","AD2"),rhs=A_rhsD[2]),
]
desc = "Calculate analytic electric field, part of the right-hand side of AD"
name = "calculate_E_exactD"
outCfunction(
outfile = os.path.join(Ccodesdir,name+".h"), desc=desc, name=name,
params ="const paramstruct *params,REAL *xx[3],const REAL *auxevol_gfs,REAL *rhs_gfs",
body = fin.FD_outputC("returnstring",A_rhsD_to_print,params="outCverbose=False"),
loopopts ="AllPoints,Read_xxs")
Output C function calculate_E_exactD() to file Start-to-Finish-UnitTests/Afield_flux_UnitTest/calculate_E_exactD.h
Afield_flux_unit_test.c
: The Main C Code [Back to top]¶Now that we have our vector potential and analytic magnetic field to compare against, we will start writing our unit test. We'll also import common C functionality, define REAL
, the number of ghost zones, and the faces, and set the standard macros for NRPy+ style memory access. After the preliminaries and setting data on our grid, we compare the output between calculate_E_exactD()
and calculate_E_field_flat_all_in_one()
.
%%writefile $Ccodesdir/Afield_flux_unit_test.c
// These are common packages that we are likely to need.
#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include "string.h" // Needed for strncmp, etc.
#include "stdint.h" // Needed for Windows GCC 6.x compatibility
#include <time.h> // Needed to set a random seed.
#define REAL double
#include "declare_Cparameters_struct.h"
REAL a,b,c,d,e,f,g,h,l,m,n,o,p,q,r,s,t;
// Standard NRPy+ memory access:
#define IDX3S(i,j,k) ( (i) + Nxx_plus_2NGHOSTS0 * ( (j) + Nxx_plus_2NGHOSTS1 * (k) ) )
#define IDX4S(g,i,j,k) \
( (i) + Nxx_plus_2NGHOSTS0 * ( (j) + Nxx_plus_2NGHOSTS1 * ( (k) + Nxx_plus_2NGHOSTS2 * (g) ) ) )
#define IDX4ptS(g,idx) ( (idx) + (Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1*Nxx_plus_2NGHOSTS2) * (g) )
#define LOOP_REGION(i0min,i0max, i1min,i1max, i2min,i2max) \
for(int i2=i2min;i2<i2max;i2++) for(int i1=i1min;i1<i1max;i1++) for(int i0=i0min;i0<i0max;i0++)
Writing Start-to-Finish-UnitTests/Afield_flux_UnitTest//Afield_flux_unit_test.c
We'll now define the gridfunction names.
%%writefile -a $Ccodesdir/Afield_flux_unit_test.c
// Let's also #define the NRPy+ gridfunctions
#define VALENCIAVU0GF 0
#define VALENCIAVU1GF 1
#define VALENCIAVU2GF 2
#define BU0GF 3
#define BU1GF 4
#define BU2GF 5
#define VALENCIAV_RU0GF 6
#define VALENCIAV_RU1GF 7
#define VALENCIAV_RU2GF 8
#define B_RU0GF 9
#define B_RU1GF 10
#define B_RU2GF 11
#define VALENCIAV_LU0GF 12
#define VALENCIAV_LU1GF 13
#define VALENCIAV_LU2GF 14
#define B_LU0GF 15
#define B_LU1GF 16
#define B_LU2GF 17
#define GAMMA_FACEDD00GF 18
#define GAMMA_FACEDD01GF 19
#define GAMMA_FACEDD02GF 20
#define GAMMA_FACEDD11GF 21
#define GAMMA_FACEDD12GF 22
#define GAMMA_FACEDD22GF 23
#define BETA_FACEU0GF 24
#define BETA_FACEU1GF 25
#define BETA_FACEU2GF 26
#define ALPHA_FACEGF 27
#define GAMMADD00GF 28
#define GAMMADD01GF 29
#define GAMMADD02GF 30
#define GAMMADD11GF 31
#define GAMMADD12GF 32
#define GAMMADD22GF 33
#define BETAU0GF 34
#define BETAU1GF 35
#define BETAU2GF 36
#define ALPHAGF 37
#define GAMMA_FACEUU00GF 38
#define GAMMA_FACEUU01GF 39
#define GAMMA_FACEUU02GF 40
#define GAMMA_FACEUU11GF 41
#define GAMMA_FACEUU12GF 42
#define GAMMA_FACEUU22GF 43
#define GAMMAUU00GF 44
#define GAMMAUU01GF 45
#define GAMMAUU02GF 46
#define GAMMAUU11GF 47
#define GAMMAUU12GF 48
#define GAMMAUU22GF 49
#define NUM_AUXEVOL_GFS 50
#define AD0GF 0
#define AD1GF 1
#define AD2GF 2
#define NUM_EVOL_GFS 3
Appending to Start-to-Finish-UnitTests/Afield_flux_UnitTest//Afield_flux_unit_test.c
Now, we'll handle the different A2B codes. There are several things to do here. First, we'll add #include
s to the C code so that we have access to the functions we want to test. We must also create a directory and copy the files to that directory. We will choose to do this in the subfolder A2B
relative to this tutorial.
%%writefile -a $Ccodesdir/Afield_flux_unit_test.c
// Some specific definitions needed for this file
typedef struct __gf_and_gz_struct__ {
REAL *gf;
int gz_lo[4],gz_hi[4];
} gf_and_gz_struct;
const int VX=0,VY=1,VZ=2,BX=3,BY=4,BZ=5;
const int NUM_RECONSTRUCT_GFS = 6;
#include "PPM/reconstruct_set_of_prims_PPM_GRFFE_NRPy.c"
#include "PPM/loop_defines_reconstruction_NRPy.h"
#include "FCVAL/interpolate_metric_gfs_to_cell_faces.h"
#include "calculate_metric_gfs.h"
#include "calculate_BU.h"
#include "calculate_ValenciavU.h"
#include "calculate_E_exactD.h"
// These are the functions we want to test.
#include "RHSs/calculate_E_field_flat_all_in_one.h"
Appending to Start-to-Finish-UnitTests/Afield_flux_UnitTest//Afield_flux_unit_test.c
Now, we'll write the main method. First, we'll set up the grid. In this test, we cannot use only one point. To properly perform reconstruction and have a point at the center of the grid that's ideal for convergence testing, we'll need at least a $9 \times 9 \times 9$ grid; for the higher resolution, we'll double the interior points in one direction, yielding a $11 \times 11 \times 11$ grid. Then, we'll write the magnetic fields and drift velocity. After that, we'll calculate the electric field contribution to the $A_i$ RHS in two ways.
%%writefile -a $Ccodesdir/Afield_flux_unit_test.c
int main(int argc, const char *argv[]) {
paramstruct params;
#include "set_Cparameters_default.h"
// Step 0c: Set free parameters, overwriting Cparameters defaults
// by hand or with command-line input, as desired.
#include "free_parameters.h"
#include "set_Cparameters-nopointer.h"
// We'll define our grid slightly different from how we normally would. We let our outermost
// ghostzones coincide with xxmin and xxmax instead of the interior of the grid. This means
// that the ghostzone points will have identical positions so we can do convergence tests of them. // Step 0d.ii: Set up uniform coordinate grids
REAL *xx[3];
xx[0] = (REAL *)malloc(sizeof(REAL)*Nxx_plus_2NGHOSTS0);
xx[1] = (REAL *)malloc(sizeof(REAL)*Nxx_plus_2NGHOSTS1);
xx[2] = (REAL *)malloc(sizeof(REAL)*Nxx_plus_2NGHOSTS2);
for(int j=0;j<Nxx_plus_2NGHOSTS0;j++) xx[0][j] = xxmin[0] + (j-NGHOSTS)*dxx0;
for(int j=0;j<Nxx_plus_2NGHOSTS1;j++) xx[1][j] = xxmin[1] + (j-NGHOSTS)*dxx1;
for(int j=0;j<Nxx_plus_2NGHOSTS2;j++) xx[2][j] = xxmin[2] + (j-NGHOSTS)*dxx2;
//for(int j=0;j<Nxx_plus_2NGHOSTS0;j++) printf("x[%d] = %.5e\n",j,xx[0][j]);
// This is the array to which we'll write the NRPy+ variables.
REAL *auxevol_gfs = (REAL *)malloc(sizeof(REAL) * NUM_AUXEVOL_GFS * Nxx_plus_2NGHOSTS2 * Nxx_plus_2NGHOSTS1 * Nxx_plus_2NGHOSTS0);
REAL *rhs_gfs = (REAL *)malloc(sizeof(REAL) * NUM_EVOL_GFS * Nxx_plus_2NGHOSTS2 * Nxx_plus_2NGHOSTS1 * Nxx_plus_2NGHOSTS0);
REAL *rhs_exact_gfs = (REAL *)malloc(sizeof(REAL) * NUM_EVOL_GFS * Nxx_plus_2NGHOSTS2 * Nxx_plus_2NGHOSTS1 * Nxx_plus_2NGHOSTS0);
a = (double)(rand()%20-10);
b = (double)(rand()%20-10);
c = (double)(rand()%20-10);
d = (double)(rand()%20-10);
e = (double)(rand()%20-10);
f = (double)(rand()%20-10);
g = (double)(rand()%20-10);
h = (double)(rand()%20-10);
l = (double)(rand()%20-10);
m = (double)(rand()%20-10);
n = (double)(rand()%20-10);
o = (double)(rand()%20-10);
p = (double)(rand()%20-10);
q = (double)(rand()%20-10);
r = (double)(rand()%20-10);
s = (double)(rand()%20-10);
t = (double)(rand()%20-10);
//printf("a,b,c,d,e,f = %f,%f,%f,%f,%f,%f\n",a,b,c,d,e,f);
// First, calculate the test metric on our grid:
calculate_metric_gfs(¶ms,xx,auxevol_gfs);
// Calculate the initial data and the exact solution.
calculate_BU(¶ms,xx,auxevol_gfs);
calculate_ValenciavU(¶ms,xx,auxevol_gfs);
//for(int j=0;j<Nxx_plus_2NGHOSTS0;j++) xx[0][j] -= 0.5*dxx0;
calculate_E_exactD(¶ms,xx,auxevol_gfs,rhs_exact_gfs);
// Set up the pointers and constants for the reconstruction procedure.
gf_and_gz_struct in_prims[NUM_RECONSTRUCT_GFS], out_prims_r[NUM_RECONSTRUCT_GFS], out_prims_l[NUM_RECONSTRUCT_GFS];
const int Nxxp2NG012 = Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1*Nxx_plus_2NGHOSTS2;
REAL temporary[Nxxp2NG012];
int ww=0;
in_prims[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAVU0GF;
out_prims_r[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_RU0GF;
out_prims_l[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_LU0GF;
ww++;
in_prims[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAVU1GF;
out_prims_r[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_RU1GF;
out_prims_l[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_LU1GF;
ww++;
in_prims[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAVU2GF;
out_prims_r[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_RU2GF;
out_prims_l[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_LU2GF;
ww++;
in_prims[ww].gf = auxevol_gfs + Nxxp2NG012*BU0GF;
out_prims_r[ww].gf = auxevol_gfs + Nxxp2NG012*B_RU0GF;
out_prims_l[ww].gf = auxevol_gfs + Nxxp2NG012*B_LU0GF;
ww++;
in_prims[ww].gf = auxevol_gfs + Nxxp2NG012*BU1GF;
out_prims_r[ww].gf = auxevol_gfs + Nxxp2NG012*B_RU1GF;
out_prims_l[ww].gf = auxevol_gfs + Nxxp2NG012*B_LU1GF;
ww++;
in_prims[ww].gf = auxevol_gfs + Nxxp2NG012*BU2GF;
out_prims_r[ww].gf = auxevol_gfs + Nxxp2NG012*B_RU2GF;
out_prims_l[ww].gf = auxevol_gfs + Nxxp2NG012*B_LU2GF;
ww++;
// Prims are defined AT ALL GRIDPOINTS, so we set the # of ghostzones to zero:
for(int i=0;i<NUM_RECONSTRUCT_GFS;i++) for(int j=1;j<=3;j++) { in_prims[i].gz_lo[j]=0; in_prims[i].gz_hi[j]=0; }
// Left/right variables are not yet defined, yet we set the # of gz's to zero by default:
for(int i=0;i<NUM_RECONSTRUCT_GFS;i++) for(int j=1;j<=3;j++) { out_prims_r[i].gz_lo[j]=0; out_prims_r[i].gz_hi[j]=0; }
for(int i=0;i<NUM_RECONSTRUCT_GFS;i++) for(int j=1;j<=3;j++) { out_prims_l[i].gz_lo[j]=0; out_prims_l[i].gz_hi[j]=0; }
int which_prims_to_reconstruct[6] = {VX, VY, VZ, BX, BY, BZ};
int num_prims_to_reconstruct = 6;
// In each direction, perform the PPM reconstruction procedure.
// Then, add the fluxes to the RHS as appropriate.
for(int flux_dirn=0;flux_dirn<3;flux_dirn++) {
// This function is housed in the file: "reconstruct_set_of_prims_PPM_GRFFE_NRPy.c"
interpolate_metric_gfs_to_cell_faces(¶ms,auxevol_gfs,flux_dirn+1);
reconstruct_set_of_prims_PPM_GRFFE_NRPy(¶ms, auxevol_gfs, flux_dirn+1, num_prims_to_reconstruct,
which_prims_to_reconstruct, in_prims, out_prims_r, out_prims_l, temporary);
for(int count=0;count<=1;count++) {
// This function is written to be general, using notation that matches the forward permutation added to AD2,
// i.e., [F_HLL^x(B^y)]_z corresponding to flux_dirn=0, count=1.
// The SIGN parameter is necessary because
// -E_z(x_i,y_j,z_k) = 0.25 ( [F_HLL^x(B^y)]_z(i+1/2,j,k)+[F_HLL^x(B^y)]_z(i-1/2,j,k)
// -[F_HLL^y(B^x)]_z(i,j+1/2,k)-[F_HLL^y(B^x)]_z(i,j-1/2,k) )
// Note the negative signs on the reversed permutation terms!
// By cyclically permuting with flux_dirn, we
// get contributions to the other components, and by incrementing count, we get the backward permutations:
// Let's suppose flux_dirn = 0. Then we will need to update Ay (count=0) and Az (count=1):
// flux_dirn=count=0 -> AD0GF+(flux_dirn+1+count)%3 = AD0GF + (0+1+0)%3=AD1GF <- Updating Ay!
// (flux_dirn)%3 = (0)%3 = 0 Vx
// (flux_dirn-count+2)%3 = (0-0+2)%3 = 2 Vz . Inputs Vx, Vz -> SIGN = -1 ; 2.0*((REAL)count)-1.0=-1 check!
// flux_dirn=0,count=1 -> AD0GF+(flux_dirn+1+count)%3 = AD0GF + (0+1+1)%3=AD2GF <- Updating Az!
// (flux_dirn)%3 = (0)%3 = 0 Vx
// (flux_dirn-count+2)%3 = (0-1+2)%3 = 1 Vy . Inputs Vx, Vy -> SIGN = +1 ; 2.0*((REAL)count)-1.0=2-1=+1 check!
// Let's suppose flux_dirn = 1. Then we will need to update Az (count=0) and Ax (count=1):
// flux_dirn=1,count=0 -> AD0GF+(flux_dirn+1+count)%3 = AD0GF + (1+1+0)%3=AD2GF <- Updating Az!
// (flux_dirn)%3 = (1)%3 = 1 Vy
// (flux_dirn-count+2)%3 = (1-0+2)%3 = 0 Vx . Inputs Vy, Vx -> SIGN = -1 ; 2.0*((REAL)count)-1.0=-1 check!
// flux_dirn=count=1 -> AD0GF+(flux_dirn+1+count)%3 = AD0GF + (1+1+1)%3=AD0GF <- Updating Ax!
// (flux_dirn)%3 = (1)%3 = 1 Vy
// (flux_dirn-count+2)%3 = (1-1+2)%3 = 2 Vz . Inputs Vy, Vz -> SIGN = +1 ; 2.0*((REAL)count)-1.0=2-1=+1 check!
// Let's suppose flux_dirn = 2. Then we will need to update Ax (count=0) and Ay (count=1):
// flux_dirn=2,count=0 -> AD0GF+(flux_dirn+1+count)%3 = AD0GF + (2+1+0)%3=AD0GF <- Updating Ax!
// (flux_dirn)%3 = (2)%3 = 2 Vz
// (flux_dirn-count+2)%3 = (2-0+2)%3 = 1 Vy . Inputs Vz, Vy -> SIGN = -1 ; 2.0*((REAL)count)-1.0=-1 check!
// flux_dirn=2,count=1 -> AD0GF+(flux_dirn+1+count)%3 = AD0GF + (2+1+1)%3=AD1GF <- Updating Ay!
// (flux_dirn)%3 = (2)%3 = 2 Vz
// (flux_dirn-count+2)%3 = (2-1+2)%3 = 0 Vx . Inputs Vz, Vx -> SIGN = +1 ; 2.0*((REAL)count)-1.0=2-1=+1 check!
calculate_E_field_flat_all_in_one(¶ms,
&auxevol_gfs[IDX4ptS(VALENCIAV_RU0GF+(flux_dirn)%3, 0)],&auxevol_gfs[IDX4ptS(VALENCIAV_RU0GF+(flux_dirn-count+2)%3, 0)],
&auxevol_gfs[IDX4ptS(VALENCIAV_LU0GF+(flux_dirn)%3, 0)],&auxevol_gfs[IDX4ptS(VALENCIAV_LU0GF+(flux_dirn-count+2)%3, 0)],
&auxevol_gfs[IDX4ptS(B_RU0GF +(flux_dirn)%3, 0)],&auxevol_gfs[IDX4ptS(B_RU0GF +(flux_dirn-count+2)%3, 0)],
&auxevol_gfs[IDX4ptS(B_LU0GF +(flux_dirn)%3, 0)],&auxevol_gfs[IDX4ptS(B_LU0GF +(flux_dirn-count+2)%3, 0)],
&auxevol_gfs[IDX4ptS(B_RU0GF +(flux_dirn-count+2)%3, 0)],
&auxevol_gfs[IDX4ptS(B_LU0GF +(flux_dirn-count+2)%3, 0)],
&auxevol_gfs[IDX4ptS(GAMMA_FACEDD00GF,0)],&auxevol_gfs[IDX4ptS(GAMMA_FACEDD01GF,0)],&auxevol_gfs[IDX4ptS(GAMMA_FACEDD02GF,0)],
&auxevol_gfs[IDX4ptS(GAMMA_FACEDD11GF,0)],&auxevol_gfs[IDX4ptS(GAMMA_FACEDD12GF,0)],&auxevol_gfs[IDX4ptS(GAMMA_FACEDD22GF,0)],
&auxevol_gfs[IDX4ptS(BETA_FACEU0GF+flux_dirn,0)],&auxevol_gfs[IDX4ptS(BETA_FACEU0GF+(flux_dirn-count+2)%3,0)],&auxevol_gfs[IDX4ptS(ALPHA_FACEGF,0)],
&rhs_gfs[IDX4ptS(AD0GF+(flux_dirn+1+count)%3,0)], 2.0*((REAL)count)-1.0, flux_dirn);
}
}
char filename[100];
sprintf(filename,"out%d-numer.txt",Nxx0);
FILE *out2D = fopen(filename, "w");
// We print the difference between approximate and exact numbers.
int i0 = Nxx_plus_2NGHOSTS0/2;
int i1 = Nxx_plus_2NGHOSTS1/2;
int i2 = Nxx_plus_2NGHOSTS2/2;
printf("Numerical: %.15e\n",rhs_gfs[IDX4S(AD0GF,i0,i1,i2)]);
printf("Analytic: %.15e\n",rhs_exact_gfs[IDX4S(AD0GF,i0,i1,i2)]);
printf("Numerical: %.15e\n",rhs_gfs[IDX4S(AD1GF,i0,i1,i2)]);
printf("Analytic: %.15e\n",rhs_exact_gfs[IDX4S(AD1GF,i0,i1,i2)]);
printf("Numerical: %.15e\n",rhs_gfs[IDX4S(AD2GF,i0,i1,i2)]);
printf("Analytic: %.15e\n\n",rhs_exact_gfs[IDX4S(AD2GF,i0,i1,i2)]);
//printf("i0,i1,i2 = %d,%d,%d\n",i0,i1,i2);
fprintf(out2D,"%.16e\t%.16e\t%.16e\t %e %e %e\n",
rhs_exact_gfs[IDX4S(AD0GF,i0,i1,i2)]-rhs_gfs[IDX4S(AD0GF,i0,i1,i2)],
rhs_exact_gfs[IDX4S(AD1GF,i0,i1,i2)]-rhs_gfs[IDX4S(AD1GF,i0,i1,i2)],
rhs_exact_gfs[IDX4S(AD2GF,i0,i1,i2)]-rhs_gfs[IDX4S(AD2GF,i0,i1,i2)],
xx[0][i0],xx[1][i1],xx[2][i2]
);
fclose(out2D);
}
Appending to Start-to-Finish-UnitTests/Afield_flux_UnitTest//Afield_flux_unit_test.c
Now that we have our file, we can compile it and run the executable.
cmd.C_compile(os.path.join(Ccodesdir,"Afield_flux_unit_test.c"), os.path.join(outdir,"Afield_flux_unit_test"))
# Change to output directory
os.chdir(outdir)
cmd.Execute(os.path.join("Afield_flux_unit_test"),"2 2 2","out2.txt")
# To do a convergence test, we'll also need a second grid with twice the resolution.
cmd.Execute(os.path.join("Afield_flux_unit_test"),"4 4 4","out4.txt")
Compiling executable... (EXEC): Executing `gcc -Ofast -fopenmp -march=native -funroll-loops Start-to-Finish-UnitTests/Afield_flux_UnitTest/Afield_flux_unit_test.c -o Start-to-Finish-UnitTests/Afield_flux_UnitTest/output/Afield_flux_unit_test -lm`... (BENCH): Finished executing in 0.6064903736114502 seconds. Finished compilation. (EXEC): Executing `taskset -c 0,1,2,3,4,5 ./Afield_flux_unit_test 2 2 2`... (BENCH): Finished executing in 0.20341801643371582 seconds. (EXEC): Executing `taskset -c 0,1,2,3,4,5 ./Afield_flux_unit_test 4 4 4`... (BENCH): Finished executing in 0.20840764045715332 seconds.
For testing purposes, we have only checked these algorithms on a small grid. By construction, we have only guaranteed ourselves output from the functions we are testing at a point, so we will simply print the convergence order at that point after processing our outputs below.
import numpy as np
import matplotlib.pyplot as plt
Data1 = np.loadtxt("out2-numer.txt")
Data2 = np.loadtxt("out4-numer.txt")
""
def conv_order(lowres,highres):
order = np.log2(np.abs(lowres/highres))
if Use_Shock_Data:
if order < 0.7 or order > 2.3:
sys.exit(1)
else:
if order < 1.7 or order > 2.3:
sys.exit(1)
return order
print("The following quantities converge at the listed order (should be ~1 for Shock data, ~2 otherwise):")
print("A_rhsD0: "+str(conv_order(Data1[0],Data2[0])))
print("A_rhsD1: "+str(conv_order(Data1[1],Data2[1])))
print("A_rhsD2: "+str(conv_order(Data1[2],Data2[2])))
The following quantities converge at the listed order (should be ~1 for Shock data, ~2 otherwise): A_rhsD0: 2.0234555321265675 A_rhsD1: 1.0081097780531285 A_rhsD2: 1.1291244848636108
The following code cell converts this Jupyter notebook into a proper, clickable $\LaTeX$-formatted PDF file. After the cell is successfully run, the generated PDF may be found in the root NRPy+ tutorial directory, with filename Tutorial-Start_to_Finish_UnitTest-GiRaFFE_NRPy-A2B.pdf (Note that clicking on this link may not work; you may need to open the PDF file through another means.)
import cmdline_helper as cmd # NRPy+: Multi-platform Python command-line interface
# Change to NRPy directory
os.chdir("../../../")
cmd.output_Jupyter_notebook_to_LaTeXed_PDF("Tutorial-Start_to_Finish_UnitTest-GiRaFFE_NRPy-Afield_flux",location_of_template_file=os.path.join(".."))
Created Tutorial-Start_to_Finish_UnitTest-GiRaFFE_NRPy-Afield_flux.tex, and compiled LaTeX file to PDF file Tutorial-Start_to_Finish_UnitTest- GiRaFFE_NRPy-Afield_flux.pdf