{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"# Tutorial-IllinoisGRMHD: compute_B_and_Bstagger_from_A.C\n",
"\n",
"## Authors: Leo Werneck & Zach Etienne\n",
"\n",
"**This module is currently under development**\n",
"\n",
"## This notebook outlines the computation of $B^{i}$ and $B^{i}_{\\rm stagger}$ from $A_{i}$.\n",
"\n",
"### Required and recommended citations:\n",
"\n",
"* **(Required)** Etienne, Z. B., Paschalidis, V., Haas R., Mösta P., and Shapiro, S. L. IllinoisGRMHD: an open-source, user-friendly GRMHD code for dynamical spacetimes. Class. Quantum Grav. 32 (2015) 175009. ([arxiv:1501.07276](http://arxiv.org/abs/1501.07276)).\n",
"* **(Required)** Noble, S. C., Gammie, C. F., McKinney, J. C., Del Zanna, L. Primitive Variable Solvers for Conservative General Relativistic Magnetohydrodynamics. Astrophysical Journal, 641, 626 (2006) ([astro-ph/0512420](https://arxiv.org/abs/astro-ph/0512420)).\n",
"* **(Recommended)** Del Zanna, L., Bucciantini N., Londrillo, P. An efficient shock-capturing central-type scheme for multidimensional relativistic flows - II. Magnetohydrodynamics. A&A 400 (2) 397-413 (2003). DOI: 10.1051/0004-6361:20021641 ([astro-ph/0210618](https://arxiv.org/abs/astro-ph/0210618))."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Table of Contents\n",
"$$\\label{toc}$$\n",
"\n",
"This module is organized as follows\n",
"\n",
"0. [Step 0](#src_dir): **Source directory creation**\n",
"1. [Step 1](#introduction): **Introduction**\n",
"1. [Step 2](#compute_b_and_bstagger_from_a__c): **`compute_B_and_Bstagger_from_A.C`**\n",
"1. [Step n-1](#code_validation): **Code validation**\n",
"1. [Step n](#latex_pdf_output): **Output this notebook to $\\LaTeX$-formatted PDF file**"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 0: Source directory creation \\[Back to [top](#toc)\\]\n",
"$$\\label{src_dir}$$\n",
"\n",
"We will now use the [cmdline_helper.py NRPy+ module](Tutorial-Tutorial-cmdline_helper.ipynb) to create the source directory within the `IllinoisGRMHD` NRPy+ directory, if it does not exist yet."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"# Step 0: Creation of the IllinoisGRMHD source directory\n",
"# Step 0a: 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",
"# Step 0b: Load up cmdline_helper and create the directory\n",
"import cmdline_helper as cmd\n",
"IGM_src_dir_path = os.path.join(\"..\",\"src\")\n",
"cmd.mkdir(IGM_src_dir_path)\n",
"\n",
"# Step 0c: Create the output file path\n",
"outfile_path__compute_B_and_Bstagger_from_A__C = os.path.join(IGM_src_dir_path,\"compute_B_and_Bstagger_from_A.C\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 1: Introduction \\[Back to [top](#toc)\\]\n",
"$$\\label{introduction}$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 2: `compute_B_and_Bstagger_from_A.C` \\[Back to [top](#toc)\\]\n",
"$$\\label{compute_b_and_bstagger_from_a__c}$$"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Writing ../src/compute_B_and_Bstagger_from_A.C\n"
]
}
],
"source": [
"%%writefile $outfile_path__compute_B_and_Bstagger_from_A__C\n",
"#include \"cctk.h\"\n",
"#include \n",
"#include \n",
"#include \n",
"#include \n",
"#include \"cctk_Arguments.h\"\n",
"#include \"cctk_Parameters.h\"\n",
"#include \"IllinoisGRMHD_headers.h\"\n",
"\n",
"#define LOOP_DEFINE_SIMPLE \\\n",
" _Pragma(\"omp parallel for\") \\\n",
" for(int k=0;k 0) - (0 > imax_minus_i) ),j,k);\n",
" CCTK_REAL Psi_ip1 = psi_bssn[indexip1jk];\n",
" Bx_stagger[actual_index] *= Psim3/(Psi_ip1*Psi_ip1*Psi_ip1);\n",
"\n",
" /**************/\n",
" /* By_stagger */\n",
" /**************/\n",
"\n",
" index = CCTK_GFINDEX3D(cctkGH,shiftedi,j,shiftedk);\n",
" indexim1 = CCTK_GFINDEX3D(cctkGH,shiftedim1,j,shiftedk);\n",
" indexkm1 = CCTK_GFINDEX3D(cctkGH,shiftedi,j,shiftedkm1);\n",
" // Set By_stagger = \\partial_z A_x - \\partial_x A_z\n",
" By_stagger[actual_index] = (Ax[index]-Ax[indexkm1])*dzi - (Az[index]-Az[indexim1])*dxi;\n",
"\n",
" // Now multiply By and By_stagger by 1/sqrt(gamma(i,j+1/2,k)]) = 1/sqrt(1/2 [gamma + gamma_jp1]) = exp(-6 x 1/2 [phi + phi_jp1] )\n",
" int jmax_minus_j = (cctk_lsh[1]-1)-j;\n",
" int indexijp1k = CCTK_GFINDEX3D(cctkGH,i,j + ( (jmax_minus_j > 0) - (0 > jmax_minus_j) ),k);\n",
" CCTK_REAL Psi_jp1 = psi_bssn[indexijp1k];\n",
" By_stagger[actual_index] *= Psim3/(Psi_jp1*Psi_jp1*Psi_jp1);\n",
"\n",
" /**************/\n",
" /* Bz_stagger */\n",
" /**************/\n",
"\n",
" index = CCTK_GFINDEX3D(cctkGH,shiftedi,shiftedj,k);\n",
" indexim1 = CCTK_GFINDEX3D(cctkGH,shiftedim1,shiftedj,k);\n",
" indexjm1 = CCTK_GFINDEX3D(cctkGH,shiftedi,shiftedjm1,k);\n",
" // Set Bz_stagger = \\partial_x A_y - \\partial_y A_x\n",
" Bz_stagger[actual_index] = (Ay[index]-Ay[indexim1])*dxi - (Ax[index]-Ax[indexjm1])*dyi;\n",
"\n",
" // Now multiply Bz_stagger by 1/sqrt(gamma(i,j,k+1/2)]) = 1/sqrt(1/2 [gamma + gamma_kp1]) = exp(-6 x 1/2 [phi + phi_kp1] )\n",
" int kmax_minus_k = (cctk_lsh[2]-1)-k;\n",
" int indexijkp1 = CCTK_GFINDEX3D(cctkGH,i,j,k + ( (kmax_minus_k > 0) - (0 > kmax_minus_k) ));\n",
" CCTK_REAL Psi_kp1 = psi_bssn[indexijkp1];\n",
" Bz_stagger[actual_index] *= Psim3/(Psi_kp1*Psi_kp1*Psi_kp1);\n",
"\n",
"\n",
" }\n",
"\n",
" LOOP_DEFINE_SIMPLE {\n",
" // Look Mom, no if() statements!\n",
" int shiftedim1 = (i-1)*(i!=0); // This way, i=0 yields shiftedim1=0 and shiftedi=1, used below for our COPY boundary condition.\n",
" int shiftedi = shiftedim1+1;\n",
"\n",
" int shiftedjm1 = (j-1)*(j!=0);\n",
" int shiftedj = shiftedjm1+1;\n",
"\n",
" int shiftedkm1 = (k-1)*(k!=0);\n",
" int shiftedk = shiftedkm1+1;\n",
"\n",
" int index,indexim1,indexjm1,indexkm1;\n",
"\n",
" int actual_index = CCTK_GFINDEX3D(cctkGH,i,j,k);\n",
"\n",
" // For the lower boundaries, the following applies a \"copy\"\n",
" // boundary condition on Bi and Bi_stagger where needed.\n",
" // E.g., Bx(imin,j,k) = Bx(imin+1,j,k)\n",
" // We find the copy BC works better than extrapolation.\n",
" /******/\n",
" /* Bx */\n",
" /******/\n",
" index = CCTK_GFINDEX3D(cctkGH,shiftedi,j,k);\n",
" indexim1 = CCTK_GFINDEX3D(cctkGH,shiftedim1,j,k);\n",
" // Set Bx = 0.5 ( Bx_stagger + Bx_stagger_im1 )\n",
" // \"Grid\" Bx_stagger(i,j,k) is actually Bx_stagger(i+1/2,j,k)\n",
" Bx[actual_index] = 0.5 * ( Bx_stagger[index] + Bx_stagger[indexim1] );\n",
"\n",
" /******/\n",
" /* By */\n",
" /******/\n",
" index = CCTK_GFINDEX3D(cctkGH,i,shiftedj,k);\n",
" indexjm1 = CCTK_GFINDEX3D(cctkGH,i,shiftedjm1,k);\n",
" // Set By = 0.5 ( By_stagger + By_stagger_im1 )\n",
" // \"Grid\" By_stagger(i,j,k) is actually By_stagger(i,j+1/2,k)\n",
" By[actual_index] = 0.5 * ( By_stagger[index] + By_stagger[indexjm1] );\n",
"\n",
" /******/\n",
" /* Bz */\n",
" /******/\n",
" index = CCTK_GFINDEX3D(cctkGH,i,j,shiftedk);\n",
" indexkm1 = CCTK_GFINDEX3D(cctkGH,i,j,shiftedkm1);\n",
" // Set Bz = 0.5 ( Bz_stagger + Bz_stagger_im1 )\n",
" // \"Grid\" Bz_stagger(i,j,k) is actually Bz_stagger(i,j+1/2,k)\n",
" Bz[actual_index] = 0.5 * ( Bz_stagger[index] + Bz_stagger[indexkm1] );\n",
" }\n",
"\n",
" // Finish up by setting symmetry ghostzones on Bx, By, Bz, and their staggered variants.\n",
" CCTK_REAL gridfunc_syms_Bx[3] = {-1, 1,-Sym_Bz};\n",
" IllinoisGRMHD_set_symmetry_gzs_staggered(cctkGH,cctk_lsh,x,y,z, Bx , gridfunc_syms_Bx,0,0,0);\n",
" IllinoisGRMHD_set_symmetry_gzs_staggered(cctkGH,cctk_lsh,x,y,z, Bx_stagger, gridfunc_syms_Bx,1,0,0);\n",
" CCTK_REAL gridfunc_syms_By[3] = { 1,-1,-Sym_Bz};\n",
" IllinoisGRMHD_set_symmetry_gzs_staggered(cctkGH,cctk_lsh,x,y,z, By , gridfunc_syms_By,0,0,0);\n",
" IllinoisGRMHD_set_symmetry_gzs_staggered(cctkGH,cctk_lsh,x,y,z, By_stagger, gridfunc_syms_By,0,1,0);\n",
" CCTK_REAL gridfunc_syms_Bz[3] = { 1, 1, Sym_Bz};\n",
" IllinoisGRMHD_set_symmetry_gzs_staggered(cctkGH,cctk_lsh,x,y,z, Bz , gridfunc_syms_Bz,0,0,0);\n",
" IllinoisGRMHD_set_symmetry_gzs_staggered(cctkGH,cctk_lsh,x,y,z, Bz_stagger, gridfunc_syms_Bz,0,0,1);\n",
"}\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step n-1: Code validation \\[Back to [top](#toc)\\]\n",
"$$\\label{code_validation}$$\n",
"\n",
"First we download the original `IllinoisGRMHD` source code and then compare it to the source code generated by this tutorial notebook."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Validation test for compute_B_and_Bstagger_from_A.C: FAILED!\n",
"Diff:\n",
"56c56\n",
"< // For the lower boundaries, the following applies a \"copy\" \n",
"---\n",
"> // For the lower boundaries, the following applies a \"copy\"\n",
"114c114\n",
"< \n",
"---\n",
"> \n",
"133c133\n",
"< // For the lower boundaries, the following applies a \"copy\" \n",
"---\n",
"> // For the lower boundaries, the following applies a \"copy\"\n",
"175a176\n",
"> \n"
]
}
],
"source": [
"# Verify if the code generated by this tutorial module\n",
"# matches the original IllinoisGRMHD source code\n",
"\n",
"# First download the original IllinoisGRMHD source code\n",
"import urllib\n",
"from os import path\n",
"\n",
"original_IGM_file_url = \"https://bitbucket.org/zach_etienne/wvuthorns/raw/5611b2f0b17135538c9d9d17c7da062abe0401b6/IllinoisGRMHD/src/compute_B_and_Bstagger_from_A.C\"\n",
"original_IGM_file_name = \"compute_B_and_Bstagger_from_A-original.C\"\n",
"original_IGM_file_path = os.path.join(IGM_src_dir_path,original_IGM_file_name)\n",
"\n",
"# Then download the original IllinoisGRMHD source code\n",
"# We try it here in a couple of ways in an attempt to keep\n",
"# the code more portable\n",
"try:\n",
" original_IGM_file_code = urllib.request.urlopen(original_IGM_file_url).read().decode(\"utf-8\")\n",
" # Write down the file the original IllinoisGRMHD source code\n",
" with open(original_IGM_file_path,\"w\") as file:\n",
" file.write(original_IGM_file_code)\n",
"except:\n",
" try:\n",
" original_IGM_file_code = urllib.urlopen(original_IGM_file_url).read().decode(\"utf-8\")\n",
" # Write down the file the original IllinoisGRMHD source code\n",
" with open(original_IGM_file_path,\"w\") as file:\n",
" file.write(original_IGM_file_code)\n",
" except:\n",
" # If all else fails, hope wget does the job\n",
" !wget -O $original_IGM_file_path $original_IGM_file_url\n",
"\n",
"# Perform validation\n",
"Validation__compute_B_and_Bstagger_from_A__C = !diff $original_IGM_file_path $outfile_path__compute_B_and_Bstagger_from_A__C\n",
"\n",
"if Validation__compute_B_and_Bstagger_from_A__C == []:\n",
" # If the validation passes, we do not need to store the original IGM source code file\n",
" !rm $original_IGM_file_path\n",
" print(\"Validation test for compute_B_and_Bstagger_from_A.C: PASSED!\")\n",
"else:\n",
" # If the validation fails, we keep the original IGM source code file\n",
" print(\"Validation test for compute_B_and_Bstagger_from_A.C: FAILED!\")\n",
" # We also print out the difference between the code generated\n",
" # in this tutorial module and the original IGM source code\n",
" print(\"Diff:\")\n",
" for diff_line in Validation__compute_B_and_Bstagger_from_A__C:\n",
" print(diff_line)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step n: 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-IllinoisGRMHD__compute_B_and_Bstagger_from_A.pdf](Tutorial-IllinoisGRMHD__compute_B_and_Bstagger_from_A.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": 4,
"metadata": {},
"outputs": [],
"source": [
"latex_nrpy_style_path = os.path.join(nrpy_dir_path,\"latex_nrpy_style.tplx\")\n",
"#!jupyter nbconvert --to latex --template $latex_nrpy_style_path --log-level='WARN' Tutorial-IllinoisGRMHD__compute_B_and_Bstagger_from_A.ipynb\n",
"#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__compute_B_and_Bstagger_from_A.tex\n",
"#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__compute_B_and_Bstagger_from_A.tex\n",
"#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__compute_B_and_Bstagger_from_A.tex\n",
"!rm -f Tut*.out Tut*.aux Tut*.log"
]
}
],
"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.11.1"
}
},
"nbformat": 4,
"nbformat_minor": 4
}