{
"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` C code library: Boundary conditions\n",
"\n",
"## Author: Patrick Nelson\n",
"\n",
"## This writes and documents the C code that `GiRaFFE_NRPy` uses to apply boundary conditions to the GRFFE quantities. \n",
"\n",
"**Notebook Status:** <font color=green><b> Validated </b></font>\n",
"\n",
"**Validation Notes:** This code has been validated by confirming that it produces exactly the expected output to round-off precision for arbitrary linear inputs.\n",
"\n",
"## Introduction: \n",
"The functions and macros defined here fall into one of two categories. The [first](#linear) we will work with is the `FACE_UPDATE` family, which will act on a single face, looping over each point as defined by the parameters `i0min`, `i0max`, `i1min`, `i1max`, `i2min`, and `i2max`. The parameters `FACEX0`, `FACEX1`, and `FACEX2` define which face on which we will act; that is, two of the `FACEX` parameters must be set to `NUL` (defined as 0) while the third is set to either `MAXFACE` (defined as -1) or `MINFACE` (defined as +1). For instance, if we want to fill in a ghostzone on the +x face of our grid, we must call `FACE_UPDATE` with `FACEX0 = MAXFACE`, `FACEX1 = NUL`, and `FACEX2 = NUL`. \n",
"\n",
"Care must be taken to set `i0min`, `i0max`, `i1min`, `i1max`, `i2min`, and `i2max` in such a way as to be consistent with `FACEX0`, `FACEX1`, and `FACEX2`; failure to do so can result in bad data and out-of-bounds errors. This is handled by the function `apply_bcs`, which is a part of the [second](#apply) family of functions. Functions of this type are responsible for doing so on each face in each ghostzone.\n",
"\n",
"<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](#extrap): Extrapolation Boundary Conditions\n",
" 1. [Step 1.a](#linear): Linear Extrapolation\n",
" 1. [Step 1.b](#copy): Copy Boundary Conditions\n",
" 1. [Step 1.c](#outflow): Outflow Boundary Conditions\n",
" 1. [Step 1.d](#apply): Applying the Boundary Conditions\n",
"1. [Step 2](#exact): Exact Boundary Conditions\n",
" 1. [Step 2.a](#a_i_and_vi): Setting $A_i$ and $v^i$ exactly\n",
" 1. [Step 2.b](#apply_exact): Applying the exact Boundary Conditions to $A_i$ and $v^i$ \n",
" 1. [Step 2.c](#stilded): Setting $\\tilde{S}_i$ exactly \n",
" 1. [Step 2.d](#apply_stilded): Applying the exact Boundary Conditions to $\\tilde{S}_i$\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 # Standard Python modules for multiplatform OS-level functions\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 # NRPy+: Multi-platform Python command-line interface\n",
"outdir = os.path.join(\"GiRaFFE_NRPy\",\"GiRaFFE_Ccode_validation\",\"boundary_conditions\")\n",
"cmd.mkdir(outdir)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='extrap'></a>\n",
"\n",
"## Step 1: Extrapolation Boundary Conditions \\[Back to [top](#toc)\\]\n",
"$$\\label{extrap}$$\n",
"\n",
"<a id='linear'></a>\n",
"\n",
"### Step 1.a: Linear Extrapolation \\[Back to [top](#toc)\\]\n",
"$$\\label{linear}$$\n",
"\n",
"The first `FACE_UPDATE` macro will be basic linear extrapolation conditions. It will apply boundary conditions in the specified ghostzone for the gridfunction specified by `which_gf` in the array `gfs`. That array will be passed under that name into the functions that call `FACE_UPDATE`; by convention, `gfs` is the array of evolved gridfunctions. "
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Writing GiRaFFE_NRPy/GiRaFFE_Ccode_validation/boundary_conditions/GiRaFFE_boundary_conditions.h\n"
]
}
],
"source": [
"%%writefile $outdir/GiRaFFE_boundary_conditions.h\n",
"// Currently, we're using basic Cartesian boundary conditions, pending fixes by Zach.\n",
"// Part P8a: Declare boundary condition FACE_UPDATE macro,\n",
"// which updates a single face of the 3D grid cube\n",
"// using quadratic polynomial extrapolation.\n",
"// Basic extrapolation boundary conditions\n",
"#define FACE_UPDATE(which_gf, i0min,i0max, i1min,i1max, i2min,i2max, FACEX0,FACEX1,FACEX2) \\\n",
" for(int i2=i2min;i2<i2max;i2++) for(int i1=i1min;i1<i1max;i1++) for(int i0=i0min;i0<i0max;i0++) { \\\n",
" gfs[IDX4S(which_gf,i0,i1,i2)] = \\\n",
" +2.0*gfs[IDX4S(which_gf,i0+1*FACEX0,i1+1*FACEX1,i2+1*FACEX2)] \\\n",
" -1.0*gfs[IDX4S(which_gf,i0+2*FACEX0,i1+2*FACEX1,i2+2*FACEX2)]; \\\n",
" }\n",
"// +1.0*gfs[IDX4S(which_gf,i0+3*FACEX0,i1+3*FACEX1,i2+3*FACEX2)]; \\"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='copy'></a>\n",
"\n",
"### Step 1.b: Copy Boundary Conditions \\[Back to [top](#toc)\\]\n",
"$$\\label{copy}$$\n",
"\n",
"This macro, `FACE_UPDATE_COPY`, applies copy boundary conditions. Instead of a linear extrapolation of the data in the nearest two points in the direction of the grid interior, it simply copies the data from the nearest point in the direction of the grid interior.\n",
"\n",
"We also define `MAXFACE`, `NUL`, and `MINFACE` as constants. These should be unchanging and accessible to any function anywhere in the program. "
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to GiRaFFE_NRPy/GiRaFFE_Ccode_validation/boundary_conditions/GiRaFFE_boundary_conditions.h\n"
]
}
],
"source": [
"%%writefile -a $outdir/GiRaFFE_boundary_conditions.h\n",
"// Basic Copy boundary conditions\n",
"#define FACE_UPDATE_COPY(which_gf, i0min,i0max, i1min,i1max, i2min,i2max, FACEX0,FACEX1,FACEX2) \\\n",
" for(int i2=i2min;i2<i2max;i2++) for(int i1=i1min;i1<i1max;i1++) for(int i0=i0min;i0<i0max;i0++) { \\\n",
" gfs[IDX4S(which_gf,i0,i1,i2)] = gfs[IDX4S(which_gf,i0+1*FACEX0,i1+1*FACEX1,i2+1*FACEX2)]; \\\n",
" }\n",
"\n",
"// Part P8b: Boundary condition driver routine: Apply BCs to all six\n",
"// boundary faces of the cube, filling in the innermost\n",
"// ghost zone first, and moving outward.\n",
"const int MAXFACE = -1;\n",
"const int NUL = +0;\n",
"const int MINFACE = +1;\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='outflow'></a>\n",
"\n",
"### Step 1.c: Outflow Boundary Conditions \\[Back to [top](#toc)\\]\n",
"$$\\label{outflow}$$\n",
"\n",
"This macro, `FACE_UPDATE_OUTFLOW`, is poorly named at the moment; currently, it is a clone of the macro `FACE_UPDATE` that acts on the array `aux_gfs` instead of `gfs`. However, take note of the commented code below the macro these lines can be used to implement outflow boundary conditions with linear extrapolation. For that algorithm, the macro will accept a `which_gf_0` parameter instead of `which_gf`, and operate on the gridfunctions `which_gf_0+0`, `which_gf_0+1`, and `which_gf_0+2` (that is, the macro will act on an entire 3-vector). This must be done since the different faces and components must be handled in slightly different ways. \n",
"\n",
"In (actual) outflow boundary conditions, if a quantity is directed inwards (e.g. if $v^x < 0$ in the +x ghostzone), it is set to zero. Otherwise, we apply the standard linear extrapolation boundary condition."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to GiRaFFE_NRPy/GiRaFFE_Ccode_validation/boundary_conditions/GiRaFFE_boundary_conditions.h\n"
]
}
],
"source": [
"%%writefile -a $outdir/GiRaFFE_boundary_conditions.h\n",
"// This macro acts differently in that it acts on an entire 3-vector of gfs, instead of 1.\n",
"// which_gf_0 corresponds to the zeroth component of that vector. The if statements only\n",
"// evaluate true if the velocity is directed inwards on the face in consideration.\n",
"#define FACE_UPDATE_OUTFLOW(which_gf, i0min,i0max, i1min,i1max, i2min,i2max, FACEX0,FACEX1,FACEX2) \\\n",
" for(int i2=i2min;i2<i2max;i2++) for(int i1=i1min;i1<i1max;i1++) for(int i0=i0min;i0<i0max;i0++) { \\\n",
" aux_gfs[IDX4S(which_gf_0,i0,i1,i2)] = \\\n",
" aux_gfs[IDX4S(which_gf_0,i0+FACEX0,i1+FACEX1,i2+FACEX2)]; \\\n",
" aux_gfs[IDX4S(which_gf_0+1,i0,i1,i2)] = \\\n",
" aux_gfs[IDX4S(which_gf_0+1,i0+FACEX0,i1+FACEX1,i2+FACEX2)]; \\\n",
" aux_gfs[IDX4S(which_gf_0+2,i0,i1,i2)] = \\\n",
" aux_gfs[IDX4S(which_gf_0+2,i0+FACEX0,i1+FACEX1,i2+FACEX2)]; \\\n",
" }\n",
"/* if(FACEX0*aux_gfs[IDX4S(which_gf_0+0,i0,i1,i2)] > 0.0) { \\\n",
" aux_gfs[IDX4S(which_gf_0+0,i0,i1,i2)] = 0.0; \\\n",
" } \\\n",
" if(FACEX1*aux_gfs[IDX4S(which_gf_0+1,i0,i1,i2)] > 0.0) { \\\n",
" aux_gfs[IDX4S(which_gf_0+1,i0,i1,i2)] = 0.0; \\\n",
" } \\\n",
" if(FACEX2*aux_gfs[IDX4S(which_gf_0+2,i0,i1,i2)] > 0.0) { \\\n",
" aux_gfs[IDX4S(which_gf_0+2,i0,i1,i2)] = 0.0; \\\n",
" } \\\n",
"*/\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='apply'></a>\n",
"\n",
"### Step 1.d: Applying the Boundary Conditions \\[Back to [top](#toc)\\]\n",
"$$\\label{apply}$$\n",
"\n",
"The second category of functions we use here is responsible for applying BCs in all the ghostzones by applying the `FACE_UPDATE` macro in the correct manner for each ghostzone on each face. So, we loop over each evolved gridfunction (that is *not* a component of `StildeD`); for each gridfunction, we first define the parameters `imin` and `imax` to specify the area just outside the interior of the grid. We then call `FACE_UPDATE` on each face of the innermost ghostzone. As we go, we'll decrement each component of `imin` and increment each component of `imax`; thus, after we have done all six faces, `imin` and `imax` specify the next-innermost ghostzone. We proceed in this manner until we have covered each ghostzone. "
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to GiRaFFE_NRPy/GiRaFFE_Ccode_validation/boundary_conditions/GiRaFFE_boundary_conditions.h\n"
]
}
],
"source": [
"%%writefile -a $outdir/GiRaFFE_boundary_conditions.h\n",
"void apply_bcs_potential(const paramstruct *restrict params,REAL *gfs) {\n",
"#include \"../set_Cparameters.h\"\n",
" // First, we apply extrapolation boundary conditions to AD\n",
"#pragma omp parallel for\n",
" for(int which_gf=0;which_gf<NUM_EVOL_GFS;which_gf++) {\n",
" if(which_gf < STILDED0GF || which_gf > STILDED2GF) {\n",
" int imin[3] = { NGHOSTS, NGHOSTS, NGHOSTS };\n",
" int imax[3] = { Nxx_plus_2NGHOSTS0-NGHOSTS, Nxx_plus_2NGHOSTS1-NGHOSTS, Nxx_plus_2NGHOSTS2-NGHOSTS };\n",
" for(int which_gz = 0; which_gz < NGHOSTS; which_gz++) {\n",
" // After updating each face, adjust imin[] and imax[]\n",
" // to reflect the newly-updated face extents.\n",
" FACE_UPDATE(which_gf, imin[0]-1,imin[0], imin[1],imax[1], imin[2],imax[2], MINFACE,NUL,NUL); imin[0]--;\n",
" FACE_UPDATE(which_gf, imax[0],imax[0]+1, imin[1],imax[1], imin[2],imax[2], MAXFACE,NUL,NUL); imax[0]++;\n",
"\n",
" FACE_UPDATE(which_gf, imin[0],imax[0], imin[1]-1,imin[1], imin[2],imax[2], NUL,MINFACE,NUL); imin[1]--;\n",
" FACE_UPDATE(which_gf, imin[0],imax[0], imax[1],imax[1]+1, imin[2],imax[2], NUL,MAXFACE,NUL); imax[1]++;\n",
"\n",
" FACE_UPDATE(which_gf, imin[0],imax[0], imin[1],imax[1], imin[2]-1,imin[2], NUL,NUL,MINFACE);\n",
" imin[2]--;\n",
" FACE_UPDATE(which_gf, imin[0],imax[0], imin[1],imax[1], imax[2],imax[2]+1, NUL,NUL,MAXFACE);\n",
" imax[2]++;\n",
" }\n",
" }\n",
" }\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This code works almost identically to the above. It applies copy boundary conditions, but is currently not in use. "
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to GiRaFFE_NRPy/GiRaFFE_Ccode_validation/boundary_conditions/GiRaFFE_boundary_conditions.h\n"
]
}
],
"source": [
"%%writefile -a $outdir/GiRaFFE_boundary_conditions.h\n",
" // Then, we apply copy boundary conditions to StildeD and psi6Phi\n",
"/*#pragma omp parallel for\n",
" for(int which_gf=3;which_gf<NUM_EVOL_GFS;which_gf++) {\n",
" int imin[3] = { NGHOSTS, NGHOSTS, NGHOSTS };\n",
" int imax[3] = { Nxx_plus_2NGHOSTS0-NGHOSTS, Nxx_plus_2NGHOSTS1-NGHOSTS, Nxx_plus_2NGHOSTS2-NGHOSTS };\n",
" for(int which_gz = 0; which_gz < NGHOSTS; which_gz++) {\n",
" // After updating each face, adjust imin[] and imax[]\n",
" // to reflect the newly-updated face extents.\n",
" FACE_UPDATE_COPY(which_gf, imin[0]-1,imin[0], imin[1],imax[1], imin[2],imax[2], MINFACE,NUL,NUL); imin[0]--;\n",
" FACE_UPDATE_COPY(which_gf, imax[0],imax[0]+1, imin[1],imax[1], imin[2],imax[2], MAXFACE,NUL,NUL); imax[0]++;\n",
"\n",
" FACE_UPDATE_COPY(which_gf, imin[0],imax[0], imin[1]-1,imin[1], imin[2],imax[2], NUL,MINFACE,NUL); imin[1]--;\n",
" FACE_UPDATE_COPY(which_gf, imin[0],imax[0], imax[1],imax[1]+1, imin[2],imax[2], NUL,MAXFACE,NUL); imax[1]++;\n",
"\n",
" FACE_UPDATE_COPY(which_gf, imin[0],imax[0], imin[1],imax[1], imin[2]-1,imin[2], NUL,NUL,MINFACE); imin[2]--;\n",
" FACE_UPDATE_COPY(which_gf, imin[0],imax[0], imin[1],imax[1], imax[2],imax[2]+1, NUL,NUL,MAXFACE); imax[2]++;\n",
" }\n",
" }*/\n",
"}\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Once again, this next set of loops operates almost identically as above, but it applies BCs to the velocities instead. "
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to GiRaFFE_NRPy/GiRaFFE_Ccode_validation/boundary_conditions/GiRaFFE_boundary_conditions.h\n"
]
}
],
"source": [
"%%writefile -a $outdir/GiRaFFE_boundary_conditions.h\n",
"void apply_bcs_velocity(const paramstruct *restrict params,REAL *aux_gfs) {\n",
"#include \"../set_Cparameters.h\"\n",
" // Apply outflow/copy boundary conditions to ValenciavU by passing VALENCIAVU0 as which_gf_0\n",
"// for(int which_gf=VALENCIAVU0GF;which_gf<=VALENCIAVU2GF;which_gf++) {\n",
" const int which_gf_0 = VALENCIAVU0GF;\n",
" int imin[3] = { NGHOSTS, NGHOSTS, NGHOSTS };\n",
" int imax[3] = { Nxx_plus_2NGHOSTS0-NGHOSTS, Nxx_plus_2NGHOSTS1-NGHOSTS, Nxx_plus_2NGHOSTS2-NGHOSTS };\n",
" for(int which_gz = 0; which_gz < NGHOSTS; which_gz++) {\n",
" FACE_UPDATE_OUTFLOW(which_gf_0, imin[0]-1,imin[0], imin[1],imax[1], imin[2],imax[2], MINFACE,NUL,NUL); imin[0]--;\n",
" FACE_UPDATE_OUTFLOW(which_gf_0, imax[0],imax[0]+1, imin[1],imax[1], imin[2],imax[2], MAXFACE,NUL,NUL); imax[0]++;\n",
"\n",
" FACE_UPDATE_OUTFLOW(which_gf_0, imin[0],imax[0], imin[1]-1,imin[1], imin[2],imax[2], NUL,MINFACE,NUL); imin[1]--;\n",
" FACE_UPDATE_OUTFLOW(which_gf_0, imin[0],imax[0], imax[1],imax[1]+1, imin[2],imax[2], NUL,MAXFACE,NUL); imax[1]++;\n",
"\n",
" FACE_UPDATE_OUTFLOW(which_gf_0, imin[0],imax[0], imin[1],imax[1], imin[2]-1,imin[2], NUL,NUL,MINFACE);\n",
" imin[2]--;\n",
" FACE_UPDATE_OUTFLOW(which_gf_0, imin[0],imax[0], imin[1],imax[1], imax[2],imax[2]+1, NUL,NUL,MAXFACE);\n",
" imax[2]++;\n",
" }\n",
"// }\n",
"}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='exact'></a>\n",
"\n",
"## Step 2: Exact Boundary Conditions \\[Back to [top](#toc)\\]\n",
"$$\\label{exact}$$\n",
"\n",
"<a id='a_i_and_vi'></a>\n",
"\n",
"### Step 2.a: Setting $A_i$ and $v^i$ exactly \\[Back to [top](#toc)\\]\n",
"$$\\label{a_i_and_vi}$$\n",
"\n",
"The next algorithms we will cover are exact boundary conditions. These are a testing tool that we can use to determine if our boundary conditions are causing a problem - since we know the exact solution to the Alfvén wave at any future time, we can simply set the boundary conditions to this value. \n",
"\n",
"**IMPORTANT:** Since we have gauge freedom in specifying the vector potential $A_i$, this vector can drift in a way that has no physical effect on the system, but will cause massive inconsistencies between the ghostzones and grid interior that will propagate when taking derivatives of $A_i$.\n",
"\n",
"Note that `FACE_UPDATE_EXACT` is a function, not a macro; this is necessary because macros will not let us include header files containing the equations that we wish to use here. This forces us to pass more parameters than we did before. This function works similarly to the initial data function we use, but instead loops over a more limited portion of the grid, as determined by parameters passed from within `apply_bcs_EXACT`. Note also that when we define the x coordinate `xx0`, we shift it by the expected distance the wave should have travelled. **TODO: also multiply by the wavespeed**"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to GiRaFFE_NRPy/GiRaFFE_Ccode_validation/boundary_conditions/GiRaFFE_boundary_conditions.h\n"
]
}
],
"source": [
"%%writefile -a $outdir/GiRaFFE_boundary_conditions.h\n",
"/*// A supplement to the boundary conditions for debugging. This will overwrite data with exact conditions\n",
"void FACE_UPDATE_EXACT(const paramstruct *restrict params,REAL *restrict xx[3],\n",
" const int n, const REAL dt,REAL *out_gfs,REAL *aux_gfs,\n",
" const int i0min,const int i0max, const int i1min,const int i1max, const int i2min,const int i2max,\n",
" const int FACEX0,const int FACEX1,const int FACEX2) {\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 xx0 = xx[0][i0]-n*dt;\n",
" REAL xx1 = xx[1][i1];\n",
" REAL xx2 = xx[2][i2];\n",
" if(xx0<=lbound) {\n",
"#include \"../GiRaFFEfood_A_v_1D_tests_left.h\"\n",
" }\n",
" else if (xx0<rbound) {\n",
"#include \"../GiRaFFEfood_A_v_1D_tests_center.h\"\n",
" }\n",
" else {\n",
"#include \"../GiRaFFEfood_A_v_1D_tests_right.h\"\n",
" }\n",
" out_gfs[IDX4S(PSI6PHIGF, i0,i1,i2)] = 0.0;\n",
" }\n",
"}\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='apply_exact'></a>\n",
"\n",
"### Step 2.b: Applying the exact Boundary Conditions to $A_i$ and $v^i$ \\[Back to [top](#toc)\\]\n",
"$$\\label{apply_exact}$$\n",
"\n",
"This function is, once again, almost identical to the first portion of `apply_bcs`. The primary difference is the version of `FACE_UPDATE` it calls; furthermore, since `FACE_UPDATE_EXACT` operates on several gridfunctions simultaneously, it is not necessary to loop over gridfunctions. "
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to GiRaFFE_NRPy/GiRaFFE_Ccode_validation/boundary_conditions/GiRaFFE_boundary_conditions.h\n"
]
}
],
"source": [
"%%writefile -a $outdir/GiRaFFE_boundary_conditions.h\n",
"void apply_bcs_EXACT(const paramstruct *restrict params,REAL *restrict xx[3],\n",
" const int n, const REAL dt,\n",
" REAL *out_gfs,REAL *aux_gfs) {\n",
"#include \"../set_Cparameters.h\"\n",
" int imin[3] = { NGHOSTS, NGHOSTS, NGHOSTS };\n",
" int imax[3] = { Nxx_plus_2NGHOSTS0-NGHOSTS, Nxx_plus_2NGHOSTS1-NGHOSTS, Nxx_plus_2NGHOSTS2-NGHOSTS };\n",
" for(int which_gz = 0; which_gz < NGHOSTS; which_gz++) {\n",
" // After updating each face, adjust imin[] and imax[]\n",
" // to reflect the newly-updated face extents.\n",
" // Right now, we only want to update the xmin and xmax faces with the exact data.\n",
" FACE_UPDATE_EXACT(Nxx,Nxx_plus_2NGHOSTS,xx,n,dt,out_gfs,aux_gfs,imin[0]-1,imin[0], imin[1],imax[1], imin[2],imax[2], MINFACE,NUL,NUL);\n",
" imin[0]--;\n",
" FACE_UPDATE_EXACT(Nxx,Nxx_plus_2NGHOSTS,xx,n,dt,out_gfs,aux_gfs,imax[0],imax[0]+1, imin[1],imax[1], imin[2],imax[2], MAXFACE,NUL,NUL);\n",
" imax[0]++;\n",
"\n",
" FACE_UPDATE_EXACT(Nxx,Nxx_plus_2NGHOSTS,xx,n,dt,out_gfs,aux_gfs,imin[0],imax[0], imin[1]-1,imin[1], imin[2],imax[2], NUL,MINFACE,NUL);\n",
" imin[1]--;\n",
" FACE_UPDATE_EXACT(Nxx,Nxx_plus_2NGHOSTS,xx,n,dt,out_gfs,aux_gfs,imin[0],imax[0], imax[1],imax[1]+1, imin[2],imax[2], NUL,MAXFACE,NUL);\n",
" imax[1]++;\n",
"\n",
" FACE_UPDATE_EXACT(Nxx,Nxx_plus_2NGHOSTS,xx,n,dt,out_gfs,aux_gfs,imin[0],imax[0], imin[1],imax[1], imin[2]-1,imin[2], NUL,NUL,MINFACE);\n",
" imin[2]--;\n",
" FACE_UPDATE_EXACT(Nxx,Nxx_plus_2NGHOSTS,xx,n,dt,out_gfs,aux_gfs,imin[0],imax[0], imin[1],imax[1], imax[2],imax[2]+1, NUL,NUL,MAXFACE);\n",
" imax[2]++;\n",
" }\n",
"}\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='stilded'></a>\n",
"\n",
"### Step 2.c: Setting $\\tilde{S}_i$ exactly\n",
"\\[Back to [top](#toc)\\]\n",
"$$\\label{stilded}$$\n",
"\n",
"This function covers the gap in the above algorithm by applying exact boundary conditions to `StildeD`. There are two different options given here: one includes the header file that is used in the initial data setup to calculate `StildeD` from the 3-velocity and magnetic field the other assumes that this step was already done when `out_gfs_exact` was filled at the current timestep (the `*_exact` arrays are filled at each timestep with the exact solution to allow for convergence testing) and copies the data from there. "
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to GiRaFFE_NRPy/GiRaFFE_Ccode_validation/boundary_conditions/GiRaFFE_boundary_conditions.h\n"
]
}
],
"source": [
"%%writefile -a $outdir/GiRaFFE_boundary_conditions.h\n",
"// A supplement to the boundary conditions for debugging. This will overwrite data with exact conditions\n",
"void FACE_UPDATE_EXACT_StildeD(const paramstruct *restrict params,REAL *restrict xx[3],\n",
" REAL *out_gfs,REAL *out_gfs_exact,\n",
" const int i0min,const int i0max, const int i1min,const int i1max, const int i2min,const int i2max,\n",
" const int FACEX0,const int FACEX1,const int FACEX2) {\n",
"#include \"../set_Cparameters.h\"\n",
" // This is currently modified to calculate more exact boundary conditions for StildeD. Rename if it works.\n",
" for(int i2=i2min;i2<i2max;i2++) for(int i1=i1min;i1<i1max;i1++) for(int i0=i0min;i0<i0max;i0++) {\n",
"#include \"../GiRaFFEfood_NRPy_Stilde.h\"\n",
" }\n",
" idx = IDX3(i0,i1,i2);\n",
" out_gfs[IDX4ptS(STILDED0GF,idx)] = out_gfs_exact[IDX4ptS(STILDED0GF,idx)];\n",
" out_gfs[IDX4ptS(STILDED1GF,idx)] = out_gfs_exact[IDX4ptS(STILDED1GF,idx)];\n",
" out_gfs[IDX4ptS(STILDED2GF,idx)] = out_gfs_exact[IDX4ptS(STILDED2GF,idx)];\n",
"}\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='apply_stilded'></a>\n",
"\n",
"### Step 2.d: Applying the exact Boundary Conditions to $\\tilde{S}_i$ \\[Back to [top](#toc)\\]\n",
"$$\\label{apply_stilded}$$\n",
"\n",
"This function is nearly identical to `apply_bcs_EXACT_StildeD`, but calls `FACE_UPDATE_EXACT_StildeD` instead of `FACE_UPDATE_EXACT`."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to GiRaFFE_NRPy/GiRaFFE_Ccode_validation/boundary_conditions/GiRaFFE_boundary_conditions.h\n"
]
}
],
"source": [
"%%writefile -a $outdir/GiRaFFE_boundary_conditions.h\n",
"void apply_bcs_EXACT_StildeD(const paramstruct *restrict params,REAL *restrict xx[3],\n",
" REAL *out_gfs,REAL *out_gfs_exact) {\n",
"#include \"../set_Cparameters.h\"\n",
" int imin[3] = { NGHOSTS, NGHOSTS, NGHOSTS };\n",
" int imax[3] = { Nxx_plus_2NGHOSTS0-NGHOSTS, Nxx_plus_2NGHOSTS1-NGHOSTS, Nxx_plus_2NGHOSTS2-NGHOSTS };\n",
" for(int which_gz = 0; which_gz < NGHOSTS; which_gz++) {\n",
" // After updating each face, adjust imin[] and imax[]\n",
" // to reflect the newly-updated face extents.\n",
" // Right now, we only want to update the xmin and xmax faces with the exact data.\n",
" FACE_UPDATE_EXACT_StildeD(Nxx,Nxx_plus_2NGHOSTS,xx,out_gfs,out_gfs_exact,imin[0]-1,imin[0], imin[1],imax[1], imin[2],imax[2], MINFACE,NUL,NUL);\n",
" imin[0]--;\n",
" FACE_UPDATE_EXACT_StildeD(Nxx,Nxx_plus_2NGHOSTS,xx,out_gfs,out_gfs_exact,imax[0],imax[0]+1, imin[1],imax[1], imin[2],imax[2], MAXFACE,NUL,NUL);\n",
" imax[0]++;\n",
"\n",
" //FACE_UPDATE_EXACT_StildeD(Nxx,Nxx_plus_2NGHOSTS,xx,out_gfs,out_gfs_exact,imin[0],imax[0], imin[1]-1,imin[1], imin[2],imax[2], NUL,MINFACE,NUL);\n",
" imin[1]--;\n",
" //FACE_UPDATE_EXACT_StildeD(Nxx,Nxx_plus_2NGHOSTS,xx,out_gfs,out_gfs_exact,imin[0],imax[0], imax[1],imax[1]+1, imin[2],imax[2], NUL,MAXFACE,NUL);\n",
" imax[1]++;\n",
"\n",
" //FACE_UPDATE_EXACT_StildeD(Nxx,Nxx_plus_2NGHOSTS,xx,out_gfs,out_gfs_exact,imin[0],imax[0], imin[1],imax[1], imin[2]-1,imin[2], NUL,NUL,MINFACE);\n",
" imin[2]--;\n",
" //FACE_UPDATE_EXACT_StildeD(Nxx,Nxx_plus_2NGHOSTS,xx,out_gfs,out_gfs_exact,imin[0],imax[0], imin[1],imax[1], imax[2],imax[2]+1, NUL,NUL,MAXFACE);\n",
" imax[2]++;\n",
" }\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`\n"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Printing difference between original C code and this code...\n",
"Checking file GiRaFFE_boundary_conditions.h\n",
"No difference. TEST PASSED!\n"
]
}
],
"source": [
"import difflib\n",
"import sys\n",
"\n",
"# Define the directory that we wish to validate against:\n",
"valdir = os.path.join(\"GiRaFFE_NRPy\",\"GiRaFFE_Ccode_library\",\"boundary_conditions\")\n",
"\n",
"import GiRaFFE_NRPy.GiRaFFE_NRPy_BCs as BC\n",
"BC.GiRaFFE_NRPy_BCs(valdir)\n",
"\n",
"print(\"Printing difference between original C code and this code...\")\n",
"# Open the files to compare\n",
"files_to_check = [\"GiRaFFE_boundary_conditions.h\"]\n",
"\n",
"for file in files_to_check:\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-BCs.pdf](Tutorial-GiRaFFE_NRPy-BCs.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": 13,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Created Tutorial-GiRaFFE_NRPy-BCs.tex, and compiled LaTeX file to PDF file\n",
" Tutorial-GiRaFFE_NRPy-BCs.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-BCs\")"
]
}
],
"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.11"
}
},
"nbformat": 4,
"nbformat_minor": 4
}