GiRaFFE_NRPy
: $A_k$ to $B^i$¶GiRaFFE
.¶Notebook Status: Validated
Validation Notes: This module will validate the routines in Tutorial-GiRaFFE_NRPy-A2B.
This notebook validates our A-to-B solver 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 curl operator in an arbitrary spacetime. To that end, we will choose functions that are easy to differentiate, but lack the symmetries that would trivialize the finite-difference algorithm. Higher-order polynomials are one such type of function.
When this notebook is run, if Use_Gaussian_Data
is True
, the difference between the approximate and exact magnetic field 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. Otherwise, is Use_Gaussian_Data
is False
, polynomial data will be used and the significant digits of agreement between the approximate and exact magnetic field will be printed to the screen right after the code is run here.
This notebook is organized as follows
A2B_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)
from outputC import outCfunction, 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/A2B_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-A2B"
# Set the finite-differencing order to 2
par.set_parval_from_str("finite_difference::FD_CENTDERIVS_ORDER", 2)
Use_Gaussian_Data = True
a,b,c,d,e,f,g,h,l,m,n,o,p,q,r,s,t,u = par.Cparameters("REAL",thismodule,["a","b","c","d","e","f","g","h","l","m","n","o","p","q","r","s","t","u"],1e300)
gammaDD = ixp.register_gridfunctions_for_single_rank2("AUXEVOL","gammaDD","sym01")
AD = ixp.register_gridfunctions_for_single_rank1("EVOL","AD")
BU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","BU")
We will start with the simplest case - testing the second-order solver. In second-order finite-differencing, we use a three-point stencil that can exactly differentiate polynomials up to quadratic. So, we will use cubic functions three variables. For instance,
\begin{align} A_x &= ax^3 + by^3 + cz^3 + dy^2 + ez^2 + f \\ A_y &= gx^3 + hy^3 + lz^3 + mx^2 + nz^2 + p \\ A_z &= px^3 + qy^3 + rz^3 + sx^2 + ty^2 + u. \\ \end{align}It will be much simpler to let NRPy+ handle most of this work. So, we will import the core functionality of NRPy+, build the expressions, and then output them using outputC()
.
if not Use_Gaussian_Data:
is_gaussian = par.Cparameters("int",thismodule,"is_gaussian",0)
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]
AD[0] = a*x**3 + b*y**3 + c*z**3 + d*y**2 + e*z**2 + f
AD[1] = g*x**3 + h*y**3 + l*z**3 + m*x**2 + n*z**2 + o
AD[2] = p*x**3 + q*y**3 + r*z**3 + s*x**2 + t*y**2 + u
Alternatively, we might want to use different functions for the vector potential. Here, we'll give some 3D Gaussians: \begin{align} A_x &= a e^{-((x-b)^2+(y-c)^2+(z-d)^2)} \\ A_y &= f e^{-((x-g)^2+(y-h)^2+(z-l)^2)} \\ A_z &= m e^{-((x-n)^2+(y-o)^2+(z-p)^2)}, \\ \end{align} where $e$ is the natural number.
if Use_Gaussian_Data:
is_gaussian = par.Cparameters("int",thismodule,"is_gaussian",1)
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]
AD[0] = a * sp.exp(-((x-b)**2 + (y-c)**2 + (z-d)**2))
AD[1] = f * sp.exp(-((x-g)**2 + (y-h)**2 + (z-l)**2))
AD[2] = m * sp.exp(-((x-n)**2 + (y-o)**2 + (z-p)**2))
Next, we'll let NRPy+ compute derivatives analytically according to $$B^i = \frac{[ijk]}{\sqrt{\gamma}} \partial_j A_k.$$ Then we can carry out two separate tests to verify the numerical derivatives. First, we will verify that when we let the cubic terms be zero, the two calculations of $B^i$ agree to roundoff error. Second, we will verify that when we set the cubic terms, our error is dominated by trunction error that converges to zero at the expected rate.
We will need a sample metric $\gamma_{ij}$ for $\sqrt{\gamma}$. 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}
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[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))
import GRHD.equations as gh
gh.compute_sqrtgammaDET(gammaDD)
LeviCivitaUUU = ixp.LeviCivitaTensorUUU_dim3_rank3(gh.sqrtgammaDET)
B_analyticU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","B_analyticU")
for i in range(3):
B_analyticU[i] = 0
for j in range(3):
for k in range(3):
B_analyticU[i] += LeviCivitaUUU[i][j][k] * sp.diff(AD[k],rfm.xx_to_Cart[j])
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]),\
]
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/A2B_UnitTest/calculate_metric_gfs.h
We also should write a function that will use the analytic formulae for $B^i$.
B_analyticU_to_print = [\
lhrh(lhs=gri.gfaccess("out_gfs","B_analyticU0"),rhs=B_analyticU[0]),\
lhrh(lhs=gri.gfaccess("out_gfs","B_analyticU1"),rhs=B_analyticU[1]),\
lhrh(lhs=gri.gfaccess("out_gfs","B_analyticU2"),rhs=B_analyticU[2]),\
]
desc = "Calculate the exact magnetic field"
name = "calculate_exact_BU"
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",B_analyticU_to_print,params="outCverbose=False"),
loopopts="AllPoints,Read_xxs")
Output C function calculate_exact_BU() to file Start-to-Finish-UnitTests/A2B_UnitTest/calculate_exact_BU.h
We'll now write a function to set the vector potential $A_k$. This simply uses NRPy+ to generate most of the code from the expressions we wrote at the beginning. Then, we'll need to call the function from the module GiRaFFE_NRPy_A2B
to generate the code we need. Also, we will declare the parameters for the vector potential functions.
AD_to_print = [\
lhrh(lhs=gri.gfaccess("out_gfs","AD0"),rhs=AD[0]),\
lhrh(lhs=gri.gfaccess("out_gfs","AD1"),rhs=AD[1]),\
lhrh(lhs=gri.gfaccess("out_gfs","AD2"),rhs=AD[2]),\
]
desc = "Calculate the vector potential"
name = "calculate_AD"
outCfunction(
outfile = os.path.join(Ccodesdir,name+".h"), desc=desc, name=name,
params ="const paramstruct *restrict params,REAL *restrict xx[3],REAL *restrict out_gfs",
body = fin.FD_outputC("returnstring",AD_to_print,params="outCverbose=False"),
loopopts="AllPoints,Read_xxs")
# cmd.mkdir(os.path.join(out_dir))
import GiRaFFE_NRPy.GiRaFFE_NRPy_A2B as A2B
# We'll generate these into the A2B subdirectory since that's where the functions
# we're testing expect them to be.
AD = ixp.declarerank1("AD") # Make sure these aren't analytic expressions
gammaDD = ixp.declarerank2("gammaDD","sym01")
A2B.GiRaFFE_NRPy_A2B(os.path.join(Ccodesdir,"A2B"),gammaDD,AD,BU)
Output C function calculate_AD() to file Start-to-Finish-UnitTests/A2B_UnitTest/calculate_AD.h
# Step 3.d.i: Generate declare_Cparameters_struct.h, set_Cparameters_default.h, and set_Cparameters[-SIMD].h
# par.generate_Cparameters_Ccodes(os.path.join(out_dir))
# 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.Nxx0 = atoi(argv[1]);
params.Nxx1 = atoi(argv[2]);
params.Nxx2 = atoi(argv[3]);
params.Nxx_plus_2NGHOSTS0 = params.Nxx0 + 2*NGHOSTS;
params.Nxx_plus_2NGHOSTS1 = params.Nxx1 + 2*NGHOSTS;
params.Nxx_plus_2NGHOSTS2 = params.Nxx2 + 2*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.01,-0.01,-0.01};
const REAL xxmax[3] = { 0.01, 0.01, 0.01};
params.dxx0 = (xxmax[0] - xxmin[0]) / ((REAL)params.Nxx_plus_2NGHOSTS0-1.0);
params.dxx1 = (xxmax[1] - xxmin[1]) / ((REAL)params.Nxx_plus_2NGHOSTS1-1.0);
params.dxx2 = (xxmax[2] - xxmin[2]) / ((REAL)params.Nxx_plus_2NGHOSTS2-1.0);
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))
A2B_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. Once we are done with preliminaries and defining data on our grid, we will compare the output from calculate_exact_BU()
against driver_A_to_B
to confirm that our error converges to zero.
%%writefile $Ccodesdir/A2B_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"
const int MAXFACE = -1;
const int NUL = +0;
const int MINFACE = +1;
const int NGHOSTS = 3;
const int NGHOSTS_A2B = 3;
REAL a,b,c,d,e,f,g,h,l,m,n,o,p,q,r,s,t,u;
// Standard NRPy+ memory access:
#define IDX4S(g,i,j,k) \
( (i) + Nxx_plus_2NGHOSTS0 * ( (j) + Nxx_plus_2NGHOSTS1 * ( (k) + Nxx_plus_2NGHOSTS2 * (g) ) ) )
Writing Start-to-Finish-UnitTests/A2B_UnitTest//A2B_unit_test.c
We'll now define the gridfunction names.
%%writefile -a $Ccodesdir/A2B_unit_test.c
// Let's also #define the NRPy+ gridfunctions
#define AD0GF 0
#define AD1GF 1
#define AD2GF 2
#define NUM_EVOL_GFS 3
#define GAMMADD00GF 0
#define GAMMADD01GF 1
#define GAMMADD02GF 2
#define GAMMADD11GF 3
#define GAMMADD12GF 4
#define GAMMADD22GF 5
#define B_ANALYTICU0GF 6
#define B_ANALYTICU1GF 7
#define B_ANALYTICU2GF 8
#define BU0GF 9
#define BU1GF 10
#define BU2GF 11
#define NUM_AUXEVOL_GFS 12
Appending to Start-to-Finish-UnitTests/A2B_UnitTest//A2B_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, as generated above. We will choose to do this in the subfolder A2B
relative to this tutorial.
%%writefile -a $Ccodesdir/A2B_unit_test.c
#include "A2B/driver_AtoB.h" // This file contains both functions we need.
#include "calculate_exact_BU.h"
#include "calculate_AD.h"
#include "calculate_metric_gfs.h"
Appending to Start-to-Finish-UnitTests/A2B_UnitTest//A2B_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. As we are testing a three-point stencil, we can get away with a minimal $3 \times 3 \times 3$ grid. Then, we'll write the A fields. After that, we'll calculate the magnetic field two ways.
%%writefile -a $Ccodesdir/A2B_unit_test.c
int main(int argc, const char *argv[]) {
paramstruct params;
#include "set_Cparameters_default.h"
// Let the last argument be the test we're doing. 1 = coarser grid, 0 = finer grid.
int do_quadratic_test = atoi(argv[4]);
// 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] + ((REAL)(j))*dxx0;
for(int j=0;j<Nxx_plus_2NGHOSTS1;j++) xx[1][j] = xxmin[1] + ((REAL)(j))*dxx1;
for(int j=0;j<Nxx_plus_2NGHOSTS2;j++) xx[2][j] = xxmin[2] + ((REAL)(j))*dxx2;
//for(int j=0;j<Nxx_plus_2NGHOSTS0;j++) printf("x[%d] = %.5e\n",j,xx[0][j]);
//for(int i=0;i<Nxx_plus_2NGHOSTS0;i++) printf("xx[0][%d] = %.15e\\n",i,xx[0][i]);
// 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 *evol_gfs = (REAL *)malloc(sizeof(REAL) * NUM_EVOL_GFS * Nxx_plus_2NGHOSTS2 * Nxx_plus_2NGHOSTS1 * Nxx_plus_2NGHOSTS0);
for(int i=0;i<Nxx_plus_2NGHOSTS0;i++) for(int j=0;j<Nxx_plus_2NGHOSTS1;j++) for(int k=0;k<Nxx_plus_2NGHOSTS1;k++) {
auxevol_gfs[IDX4S(BU0GF,i,j,k)] = 0.0;
auxevol_gfs[IDX4S(BU1GF,i,j,k)] = 0.0;
auxevol_gfs[IDX4S(BU2GF,i,j,k)] = 0.0;
}
// We now want to set up the vector potential. First, we must set the coefficients.
if(is_gaussian) {
// Gaussian coefficients:
// Magnitudes:
a = (double)(rand()%20)/5.0;
f = (double)(rand()%20)/5.0;
m = (double)(rand()%20)/5.0;
// Offsets:
b = (double)(rand()%10-5)/1000.0;
c = (double)(rand()%10-5)/1000.0;
d = (double)(rand()%10-5)/1000.0;
g = (double)(rand()%10-5)/1000.0;
h = (double)(rand()%10-5)/1000.0;
l = (double)(rand()%10-5)/1000.0;
n = (double)(rand()%10-5)/1000.0;
o = (double)(rand()%10-5)/1000.0;
p = (double)(rand()%10-5)/1000.0;
/*printf("Offsets: b,c,d = %f,%f,%f\n",b,c,d);
printf("Offsets: g,h,l = %f,%f,%f\n",g,h,l);
printf("Offsets: n,o,p = %f,%f,%f\n",n,o,p);*/
// First, calculate the test data on our grid:
}
else {
// Polynomial coefficients
// We will use random integers between -10 and 10. For the first test, we let the
// Cubic coefficients remain zero. Those are a,b,c,g,h,l,p,q, and r.
d = (double)(rand()%20-10);
e = (double)(rand()%20-10);
f = (double)(rand()%20-10);
m = (double)(rand()%20-10);
n = (double)(rand()%20-10);
o = (double)(rand()%20-10);
s = (double)(rand()%20-10);
t = (double)(rand()%20-10);
u = (double)(rand()%20-10);
}
calculate_metric_gfs(¶ms,xx,auxevol_gfs);
if(do_quadratic_test && !is_gaussian) {
calculate_AD(¶ms,xx,evol_gfs);
// We'll also calculate the exact solution for B^i
calculate_exact_BU(¶ms,xx,auxevol_gfs);
// And now for the numerical derivatives:
driver_A_to_B(¶ms,evol_gfs,auxevol_gfs);
printf("This test uses quadratic vector potentials, so the magnetic fields should agree to roundoff error.\n");
printf("Below, each row represents one point. Each column represents a component of the magnetic field.\n");
printf("Shown is the number of Significant Digits of Agreement, at least 13 is good, higher is better:\n\n");
for(int i2=0;i2<Nxx_plus_2NGHOSTS2;i2++) for(int i1=0;i1<Nxx_plus_2NGHOSTS1;i1++) for(int i0=0;i0<Nxx_plus_2NGHOSTS0;i0++) {
printf("i0,i1,i2=%d,%d,%d; SDA: %.3f, %.3f, %.3f\n",i0,i1,i2,
1.0-log10(2.0*fabs(auxevol_gfs[IDX4S(B_ANALYTICU0GF,i0,i1,i2)]-auxevol_gfs[IDX4S(BU0GF,i0,i1,i2)])/(fabs(auxevol_gfs[IDX4S(B_ANALYTICU0GF,i0,i1,i2)])+fabs(auxevol_gfs[IDX4S(BU0GF,i0,i1,i2)])+1.e-15)),
1.0-log10(2.0*fabs(auxevol_gfs[IDX4S(B_ANALYTICU1GF,i0,i1,i2)]-auxevol_gfs[IDX4S(BU1GF,i0,i1,i2)])/(fabs(auxevol_gfs[IDX4S(B_ANALYTICU1GF,i0,i1,i2)])+fabs(auxevol_gfs[IDX4S(BU1GF,i0,i1,i2)])+1.e-15)),
1.0-log10(2.0*fabs(auxevol_gfs[IDX4S(B_ANALYTICU2GF,i0,i1,i2)]-auxevol_gfs[IDX4S(BU2GF,i0,i1,i2)])/(fabs(auxevol_gfs[IDX4S(B_ANALYTICU2GF,i0,i1,i2)])+fabs(auxevol_gfs[IDX4S(BU2GF,i0,i1,i2)])+1.e-15))
);
/*printf("%.3f, %.3f, %.3f\n",
auxevol_gfs[IDX4S(BU0GF,i0,i1,i2)],
auxevol_gfs[IDX4S(BU1GF,i0,i1,i2)],
auxevol_gfs[IDX4S(BU2GF,i0,i1,i2)]
);*/
}
}
if(!is_gaussian) {
// Now, we'll set the cubic coefficients:
a = (double)(rand()%20-10);
b = (double)(rand()%20-10);
c = (double)(rand()%20-10);
g = (double)(rand()%20-10);
h = (double)(rand()%20-10);
l = (double)(rand()%20-10);
p = (double)(rand()%20-10);
q = (double)(rand()%20-10);
r = (double)(rand()%20-10);
// First, calculate the test data on our grid:
calculate_metric_gfs(¶ms,xx,auxevol_gfs);
}
// And recalculate on our initial grid:
calculate_AD(¶ms,xx,evol_gfs);
// We'll also calculate the exact solution for B^i
calculate_exact_BU(¶ms,xx,auxevol_gfs);
// And now for the numerical derivatives:
driver_A_to_B(¶ms,evol_gfs,auxevol_gfs);
char filename[100];
sprintf(filename,"out%d-numer.txt",Nxx0);
FILE *out2D = fopen(filename, "w");
for(int i2=0;i2<Nxx_plus_2NGHOSTS2;i2++) for(int i1=0;i1<Nxx_plus_2NGHOSTS1;i1++) for(int i0=0;i0<Nxx_plus_2NGHOSTS0;i0++) {
// We print the difference between approximate and exact numbers.
fprintf(out2D,"%.16e\t%.16e\t%.16e %e %e %e\n",
auxevol_gfs[IDX4S(B_ANALYTICU0GF,i0,i1,i2)]-auxevol_gfs[IDX4S(BU0GF,i0,i1,i2)],
auxevol_gfs[IDX4S(B_ANALYTICU1GF,i0,i1,i2)]-auxevol_gfs[IDX4S(BU1GF,i0,i1,i2)],
auxevol_gfs[IDX4S(B_ANALYTICU2GF,i0,i1,i2)]-auxevol_gfs[IDX4S(BU2GF,i0,i1,i2)],
xx[0][i0],xx[1][i1],xx[2][i2]
);
}
fclose(out2D);
}
Appending to Start-to-Finish-UnitTests/A2B_UnitTest//A2B_unit_test.c
import time
print("Now compiling, should take ~2 seconds...\n")
start = time.time()
cmd.C_compile(os.path.join(Ccodesdir,"A2B_unit_test.c"), os.path.join(outdir,"A2B_unit_test"))
end = time.time()
print("Finished in "+str(end-start)+" seconds.\n\n")
# Change to output directory
os.chdir(outdir)
print("Now running...\n")
start = time.time()
cmd.Execute(os.path.join("A2B_unit_test"),"1 1 1 1")
if Use_Gaussian_Data:
# To do a convergence test, we'll also need a second grid with twice the resolution.
cmd.Execute(os.path.join("A2B_unit_test"),"7 7 7 1")
end = time.time()
print("Finished in "+str(end-start)+" seconds.\n\n")
Now compiling, should take ~2 seconds... Compiling executable... (EXEC): Executing `gcc -Ofast -fopenmp -march=native -funroll-loops Start-to-Finish-UnitTests/A2B_UnitTest/A2B_unit_test.c -o Start-to-Finish-UnitTests/A2B_UnitTest/output/A2B_unit_test -lm`... (BENCH): Finished executing in 0.6097359657287598 seconds. Finished compilation. Finished in 0.6173965930938721 seconds. Now running... (EXEC): Executing `taskset -c 0,1,2,3,4,5 ./A2B_unit_test 1 1 1 1`... (BENCH): Finished executing in 0.21168088912963867 seconds. (EXEC): Executing `taskset -c 0,1,2,3,4,5 ./A2B_unit_test 7 7 7 1`... (BENCH): Finished executing in 0.20528268814086914 seconds. Finished in 0.4351034164428711 seconds.
Now that we have shown that when we use a quadratic vector potential, we get roundoff-level agreement (which is to be expected, since the finite-differencing used approximates the underlying function with a quadratic), we will use do a convergence test to show that when we can't exactly model the function, the truncation error dominates and converges to zero at the expected rate. For this, we use cubic functions for the vector potential. In the code above, we output the difference beteween the numeric and exact magnetic fields at the overlapping, non-edge, non-vertex points of two separate grids. Here, we import that data and calculate the convergence in the usual way, $$ k = \log_2 \left( \frac{F - F_1}{F - F_2} \right), $$ where $k$ is the convergence order, $F$ is the exact solution, $F_1$ is the approximate solution on the coarser grid with resolution $\Delta x$, and $F_2$ is the approximate solution on the finer grid with resolution $\Delta x/2$.
Here, we will calculate the convergence of the L2 Norm over the points in each region: $$ | B^i_{\rm approx} - B^i_{\rm exact}| = \sqrt{\frac{1}{N} \sum_{ijk} \left( B^i_{\rm approx} - B^i_{\rm exact} \right)^2} $$
import numpy as np
import matplotlib.pyplot as plt
Data1 = np.loadtxt("out1-numer.txt")
Data2 = np.loadtxt("out7-numer.txt")
"A"
# print("Convergence test: All should be approximately 2\n")
# convergence = np.log(np.divide(np.abs(Data1),np.abs(Data2)))/np.log(2)
# for i in range(len(convergence[:,0])):
# print(convergence[i,:])
def IDX4(i,j,k,Nxx_plus_2NGHOSTS0,Nxx_plus_2NGHOSTS1,Nxx_plus_2NGHOSTS2):
return (i) + Nxx_plus_2NGHOSTS0 * ( (j) + Nxx_plus_2NGHOSTS1 * ( (k) + Nxx_plus_2NGHOSTS2 * (0) ) )
comp = 0 # 0->Bx, 1->By, 2->Bz
# First, let's do this over the interior
N = 7 # This is the number of total gridpoints
nface = 0 # This is the number of points we are taking the norm of.
nint = 0 # This is the number of points we are taking the norm of.
L2_1 = 0
L2_1_xm = 0 # We declare one L2 norm for each face.
L2_1_xp = 0
L2_1_ym = 0
L2_1_yp = 0
L2_1_zm = 0
L2_1_zp = 0
for k in range(N):
for j in range(N):
for i in range(N):
if i==0:
L2_1_xm += Data1[IDX4(i,j,k,N,N,N),comp]**2
nface += 1
if i==N-1:
L2_1_xp += Data1[IDX4(i,j,k,N,N,N),comp]**2
if j==0:
L2_1_ym += Data1[IDX4(i,j,k,N,N,N),comp]**2
if j==N-1:
L2_1_yp += Data1[IDX4(i,j,k,N,N,N),comp]**2
if k==0:
L2_1_zm += Data1[IDX4(i,j,k,N,N,N),comp]**2
if k==N-1:
L2_1_zp += Data1[IDX4(i,j,k,N,N,N),comp]**2
if not (i%(N-1)==0 or j%(N-1)==0 or k%(N-1)==0):
L2_1 += Data1[IDX4(i,j,k,N,N,N),comp]**2
nint += 1
L2_1 = np.sqrt(L2_1/(nint))
L2_1_xm = np.sqrt(L2_1_xm/(nface))
L2_1_xp = np.sqrt(L2_1_xp/(nface))
L2_1_ym = np.sqrt(L2_1_ym/(nface))
L2_1_yp = np.sqrt(L2_1_yp/(nface))
L2_1_zm = np.sqrt(L2_1_zm/(nface))
L2_1_zp = np.sqrt(L2_1_zp/(nface))
N = 13 # This is the number of total gridpoints
nface = 0 # This is the number of points we are taking the norm of.
nint = 0 # This is the number of points we are taking the norm of.
L2_2 = 0
L2_2_xm = 0
L2_2_xp = 0
L2_2_ym = 0
L2_2_yp = 0
L2_2_zm = 0
L2_2_zp = 0
for k in range(N):
for j in range(N):
for i in range(N):
if i==0:
L2_2_xm += Data2[IDX4(i,j,k,N,N,N),comp]**2
nface += 1
if i==N-1:
L2_2_xp += Data2[IDX4(i,j,k,N,N,N),comp]**2
if j==0:
L2_2_ym += Data2[IDX4(i,j,k,N,N,N),comp]**2
if j==N-1:
L2_2_yp += Data2[IDX4(i,j,k,N,N,N),comp]**2
if k==0:
L2_2_zm += Data2[IDX4(i,j,k,N,N,N),comp]**2
if k==N-1:
L2_2_zp += Data2[IDX4(i,j,k,N,N,N),comp]**2
if not (i%(N-1)==0 or j%(N-1)==0 or k%(N-1)==0):
L2_2 += Data2[IDX4(i,j,k,N,N,N),comp]**2
nint += 1
L2_2 = np.sqrt(L2_2/(nint))
L2_2_xm = np.sqrt(L2_2_xm/(nface))
L2_2_xp = np.sqrt(L2_2_xp/(nface))
L2_2_ym = np.sqrt(L2_2_ym/(nface))
L2_2_yp = np.sqrt(L2_2_yp/(nface))
L2_2_zm = np.sqrt(L2_2_zm/(nface))
L2_2_zp = np.sqrt(L2_2_zp/(nface))
def conv_order(low,high):
order = np.log2(low/high)
if order < 1.6:
sys.exit(1)
return order
print("Face | Res | L2 norm | Conv. Order")
print(" Int | Dx | " + "{:.7f}".format(L2_1) + " | -- ")
print(" -- | Dx/2 | " + "{:.7f}".format(L2_2) + " | " + "{:.5f}".format(conv_order(L2_1,L2_2)))
print(" -x | Dx | " + "{:.7f}".format(L2_1_xm) + " | -- ")
print(" -- | Dx/2 | " + "{:.7f}".format(L2_2_xm) + " | " + "{:.5f}".format(conv_order(L2_1_xm,L2_2_xm)))
print(" +x | Dx | " + "{:.7f}".format(L2_1_xp) + " | -- ")
print(" -- | Dx/2 | " + "{:.7f}".format(L2_2_xp) + " | " + "{:.5f}".format(conv_order(L2_1_xp,L2_2_xp)))
print(" -y | Dx | " + "{:.7f}".format(L2_1_ym) + " | -- ")
print(" -- | Dx/2 | " + "{:.7f}".format(L2_2_ym) + " | " + "{:.5f}".format(conv_order(L2_1_ym,L2_2_ym)))
print(" +y | Dx | " + "{:.7f}".format(L2_1_yp) + " | -- ")
print(" -- | Dx/2 | " + "{:.7f}".format(L2_2_yp) + " | " + "{:.5f}".format(conv_order(L2_1_yp,L2_2_yp)))
print(" -z | Dx | " + "{:.7f}".format(L2_1_zm) + " | -- ")
print(" -- | Dx/2 | " + "{:.7f}".format(L2_2_zm) + " | " + "{:.5f}".format(conv_order(L2_1_zm,L2_2_zm)))
print(" +z | Dx | " + "{:.7f}".format(L2_1_zp) + " | -- ")
print(" -- | Dx/2 | " + "{:.7f}".format(L2_2_zp) + " | " + "{:.5f}".format(conv_order(L2_1_zp,L2_2_zp)))
Face | Res | L2 norm | Conv. Order Int | Dx | 0.0000005 | -- -- | Dx/2 | 0.0000001 | 1.90520 -x | Dx | 0.0000008 | -- -- | Dx/2 | 0.0000002 | 2.08857 +x | Dx | 0.0000008 | -- -- | Dx/2 | 0.0000002 | 2.08857 -y | Dx | 0.0000008 | -- -- | Dx/2 | 0.0000002 | 1.64224 +y | Dx | 0.0000016 | -- -- | Dx/2 | 0.0000004 | 1.87830 -z | Dx | 0.0000010 | -- -- | Dx/2 | 0.0000002 | 2.09971 +z | Dx | 0.0000008 | -- -- | Dx/2 | 0.0000002 | 1.99483
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-A2B",location_of_template_file=os.path.join(".."))
Created Tutorial-Start_to_Finish_UnitTest-GiRaFFE_NRPy-A2B.tex, and compiled LaTeX file to PDF file Tutorial-Start_to_Finish_UnitTest- GiRaFFE_NRPy-A2B.pdf