{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "# Tutorial-IllinoisGRMHD: `driver_evaluate_MHD_rhs.C`\n", "\n", "## Authors: Leo Werneck & Zach Etienne\n", "\n", "**This module is currently under development**\n", "\n", "## In this tutorial module we explain the driver functions that compute the right-hand side (RHS) of the MHD equations within `IllinoisGRMHD`\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](#header_files): **Load up necessary ETK and `IllinoisGRMHD` header files**\n", "1. [Step 3](#driver_mhd_rhs_function): **The ` IllinoisGRMHD_driver_evaluate_MHD_rhs()` function**\n", " 1. [Step 3.1](#eos_parameters): *Equation of state (EOS) parameters*\n", " 1. [Step 3.2](#set_pointers_grmhd_gfs): *Set pointers to GRMHD gridfunctions*\n", " 1. [Step 3.3](#admbase_to_bssnbase): *Convert ADM variables to BSSN variables*\n", " 1. [Step 3.4](#pointers_metric_tmunu_gfs): *Setting up pointers to the metric and stress-energy tensor gridfunctions*\n", " 1. [Step 3.5](#initialize_rhss): *Initialization of the RHS variables*\n", " 1. [Step 3.6](#tau_rhs_ext_curv_and_tupmunu): *Compute extrinsic curvature terms from the RHS of $\\partial_{t}\\tilde\\tau$ and $T^{\\mu\\nu}$*\n", " 1. [Step 3.7](#computing_ftilde): *Computing ${\\rm ftilde}$*\n", " 1. [Step 3.8](#rhs_mhd_and_a_i): *The RHSs of $\\rho_{\\star}$, $\\tilde\\tau$, $\\tilde{S}_{i}$, and $A_{i}$*\n", " 1. [Step 3.8.1](#reconstructing_vx_vy_by_along_x): Reconstructing $\\left\\{v^{x}, v^{y}, B^{y}_{\\rm stagger}\\right\\}$ along the $x$-direction\n", " 1. [Step 3.8.2](#fluxes_x_dirn): Evaluating $\\partial_{x}\\boldsymbol{F}$\n", " 1. [Step 3.8.3](#reconstructing_vx_vy_by_along_y): Reconstructing $\\left\\{v^{x}, v^{y}, B^{y}_{\\rm stagger}\\right\\}$ along the $y$-direction\n", " 1. [Step 3.8.4](#fluxes_y_dirn): Evaluating $\\partial_{y}\\boldsymbol{F}$\n", " 1. [Step 3.8.5](#rhs_az_no_gauge_terms): Evaluating $\\left[\\partial_{t}A_{z}\\right]_{\\rm no\\ gauge\\ terms}$\n", " 1. [Step 3.8.6](#multiple_reconstructions): Multiple reconstructions\n", " 1. [Step 3.8.7](#fluxes_z_dirn): Evaluating $\\partial_{z}\\boldsymbol{F}$\n", " 1. [Step 3.8.8](#rhs_ax_no_gauge_terms): Evaluating $\\left[\\partial_{t}A_{x}\\right]_{\\rm no\\ gauge\\ terms}$\n", " 1. [Step 3.8.9](#rhs_ay_no_gauge_terms): Evaluating $\\left[\\partial_{t}A_{y}\\right]_{\\rm no\\ gauge\\ terms}$\n", " 1. [Step 3.8.10](#rhs_psi6phi_and_ai_gauge_terms): Evaluating $\\partial_{t}\\left[\\psi^{6}\\Phi\\right]$ and $\\left[\\partial_{t}A_{i}\\right]_{\\rm gauge\\ terms}$\n", "1. [Step 4](#driver_evaluate_MHD_rhs__h): **The `driver_evaluate_MHD_rhs.h` header file**\n", "1. [Step 5](#code_validation): **Code validation**\n", " 1. [Step 5.a](#code_validation_driver_evaluate_MHD_rhs__c): *`driver_evaluate_MHD_rhs.C`*\n", " 1. [Step 5.b](#code_validation_driver_evaluate_MHD_rhs__h): *`driver_evaluate_MHD_rhs.h`*\n", "1. [Step 6](#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__driver_evaluate_MHD_rhs__C = os.path.join(IGM_src_dir_path,\"driver_evaluate_MHD_rhs.C\")\n", "outfile_path__driver_evaluate_MHD_rhs__h = os.path.join(IGM_src_dir_path,\"driver_evaluate_MHD_rhs.h\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 1: Introduction \\[Back to [top](#toc)\\]\n", "$$\\label{introduction}$$\n", "\n", "We will start by creating the file `driver_evaluate_MHD_rhs.C` and writing dwn the preamble of the file, which contains useful references and information to the user.\n", "\n", "We remind the reader of the \"[generalized Lorenz gauge condition](https://arxiv.org/pdf/1207.3354.pdf)\",\n", "\n", "$$\n", "\\nabla_{\\mu}\\mathcal{A^{\\mu}} = \\xi n_{\\mu}\\mathcal{A^{\\mu}}\\ ,\n", "$$\n", "\n", "where $n_{\\mu} = \\left(\\alpha,0,0,0\\right)$ is the normal unit vector, $\\mathcal{A}_{\\mu}$ is the magnetic 4-vector potential, and $xi$ is a parameter with dimensions 1/Length, just like the $\\eta$ parameter in the gamma-driving shift condition, so its value is chosen so that the CFL condition remains satisfied." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Writing ../src/driver_evaluate_MHD_rhs.C\n" ] } ], "source": [ "%%writefile $outfile_path__driver_evaluate_MHD_rhs__C\n", "/*********************************************\n", " * Evaluate RHS of GRMHD & induction equations\n", " * (vector potential prescription), using the\n", " * generalized Lorenz gauge condition for the\n", " * EM gauge.\n", " *\n", " * Based originally on the Illinois GRMHD code,\n", " * written by Matt Duez, Yuk Tung Liu, and Branson\n", " * Stephens (original version), and then developed\n", " * primarily by Zachariah Etienne, Yuk Tung Liu,\n", " * and Vasileios Paschalidis.\n", " *\n", " * Rewritten for public release in 2013\n", " * by Zachariah B. Etienne\n", " *\n", " * References:\n", " * Original unigrid GRMHD evolution prescription:\n", " * http://arxiv.org/abs/astro-ph/0503420\n", " * Vector potential formulation in full GR:\n", " * http://arxiv.org/abs/1007.2848\n", " * Improved EM gauge conditions for AMR grids:\n", " * http://arxiv.org/abs/1110.4633\n", " * Generalized Lorenz gauge prescription:\n", " * http://arxiv.org/abs/1207.3354\n", " *\n", " * Note that the Generalized Lorenz gauge strength\n", " * parameter has units of 1/M, just like the \\eta\n", " * parameter in the gamma-driving shift condition,\n", " * so setting it too large will result in violation\n", " * of the CFL condition.\n", " *\n", " * This version of PPM implements the standard\n", " * Colella & Woodward PPM, though modified as in GRHydro\n", " * to have 3 ghostzones instead of 4.\n", " *********************************************/" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 2: Load up necessary ETK and `IllinoisGRMHD` header files \\[Back to [top](#toc)\\]\n", "$$\\label{header_files}$$\n", "\n", "Here we load all necessary ETK and `IllinoisGRMHD` file, as well as some standard C++ libraries." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/driver_evaluate_MHD_rhs.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__driver_evaluate_MHD_rhs__C\n", "\n", "\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\" /* Generic #define's and function prototypes */\n", "#include \"driver_evaluate_MHD_rhs.h\" /* Function prototypes for this file only */\n", "#include \"IllinoisGRMHD_EoS_lowlevel_functs.C\"\n", "#include \"inlined_functions.C\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 3: The ` IllinoisGRMHD_driver_evaluate_MHD_rhs()` function \\[Back to [top](#toc)\\]\n", "$$\\label{driver_mhd_rhs_function}$$\n", "\n", "This is the basic function declaration. We set up basic ETK parameters and verify double precision is being used." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/driver_evaluate_MHD_rhs.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__driver_evaluate_MHD_rhs__C\n", "\n", "\n", "extern \"C\" void IllinoisGRMHD_driver_evaluate_MHD_rhs(CCTK_ARGUMENTS) {\n", " DECLARE_CCTK_ARGUMENTS;\n", " DECLARE_CCTK_PARAMETERS;\n", " int levelnumber = GetRefinementLevel(cctkGH);\n", "\n", " if(CCTK_Equals(verbose, \"essential+iteration output\")) {\n", " CCTK_VInfo(CCTK_THORNSTRING,\"***** Iter. # %d, Lev: %d, Integrating to time: %e *****\",cctk_iteration,levelnumber,cctk_delta_time/cctk_levfac[0]+cctk_time);\n", " }\n", "\n", " if( sizeof(CCTK_REAL) < 8 ) CCTK_VError(VERR_DEF_PARAMS,\"Error: IllinoisGRMHD assumes that CCTK_REAL is a double precision number. Setting otherwise will likely cause havoc with the conserv_to_prims solver.\");\n", "\n", " if(cctk_nghostzones[0]<3 || cctk_nghostzones[1]<3 || cctk_nghostzones[2]<3) { CCTK_VError(VERR_DEF_PARAMS,\"ERROR. Need at least 3 ghostzones for IllinoisGRMHD evolutions.\"); }\n", "\n", " CCTK_REAL dX[3] = { CCTK_DELTA_SPACE(0), CCTK_DELTA_SPACE(1), CCTK_DELTA_SPACE(2) };" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 3.1: Equation of state (EOS) parameters \\[Back to [top](#toc)\\]\n", "$$\\label{eos_parameters}$$\n", "\n", "Next we set up the EOS struct, which is defined in the `IllinoisGRMHD_headers.h` header file. We set the following parameters:\n", "\n", "* $\\rm neos$: number of EOS (currently set to 1, which is a $\\Gamma$-law EOS)\n", "* $\\rm K\\_poly$: this is the constant $\\kappa$ from the polytropic EOS $P = \\kappa \\rho_{0}^{\\Gamma}$\n", "* $\\rm rho\\_poly$: $\\rho_{0}$, fluid rest-mass\n", "* $\\rm P\\_poly$: $P$, pressure\n", "* $\\rm gamma\\_th$: $\\Gamma_{\\rm th}$, the constant parameter which determines the conversion efficiency of kinetic to thermal energy at shocks\n", "* $\\rm eps\\_poly$: $\\epsilon$, specific internal energy\n", "* $\\rm k\\_poly$: $\\kappa$, polytropic EOS constant\n", "* $\\rm gamma\\_poly$: $\\Gamma$, the $\\Gamma$-law EOS constant" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/driver_evaluate_MHD_rhs.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__driver_evaluate_MHD_rhs__C\n", "\n", "\n", " /**********************************\n", " * Piecewise Polytropic EOS Patch *\n", " * Setting up the EOS struct *\n", " **********************************/\n", " /*\n", " * The short piece of code below takes care\n", " * of initializing the EOS parameters.\n", " * Please refer to the \"inlined_functions.C\"\n", " * source file for the documentation on the\n", " * function.\n", " */\n", " eos_struct eos;\n", " initialize_EOS_struct_from_input(eos);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 3.2: Set pointers to GRMHD gridfunctions \\[Back to [top](#toc)\\]\n", "$$\\label{set_pointers_grmhd_gfs}$$" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/driver_evaluate_MHD_rhs.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__driver_evaluate_MHD_rhs__C\n", "\n", "\n", " // in_prims,out_prims_r, and out_prims_l are arrays of pointers to the actual gridfunctions.\n", " gf_and_gz_struct in_prims[MAXNUMVARS],out_prims_r[MAXNUMVARS],out_prims_l[MAXNUMVARS];\n", " int which_prims_to_reconstruct[MAXNUMVARS],num_prims_to_reconstruct;\n", "\n", " /* SET POINTERS TO GRMHD GRIDFUNCTIONS */\n", " // The order here MATTERS, and must be consistent with the global variable declarations in\n", " // evaluate_MHD_rhs_headers.h (look for RHOB=0, etc.)\n", " // For example, in_prims[0] _must_ be rho_b.\n", " int ww=0;\n", " in_prims[ww].gf=rho_b; out_prims_r[ww].gf=rho_br; out_prims_l[ww].gf=rho_bl; ww++;\n", " in_prims[ww].gf=P; out_prims_r[ww].gf=Pr; out_prims_l[ww].gf=Pl; ww++;\n", " in_prims[ww].gf=vx; out_prims_r[ww].gf=vxr; out_prims_l[ww].gf=vxl; ww++;\n", " in_prims[ww].gf=vy; out_prims_r[ww].gf=vyr; out_prims_l[ww].gf=vyl; ww++;\n", " in_prims[ww].gf=vz; out_prims_r[ww].gf=vzr; out_prims_l[ww].gf=vzl; ww++;\n", " in_prims[ww].gf=Bx; out_prims_r[ww].gf=Bxr; out_prims_l[ww].gf=Bxl; ww++;\n", " in_prims[ww].gf=By; out_prims_r[ww].gf=Byr; out_prims_l[ww].gf=Byl; ww++;\n", " in_prims[ww].gf=Bz; out_prims_r[ww].gf=Bzr; out_prims_l[ww].gf=Bzl; ww++;\n", " in_prims[ww].gf=Bx_stagger; out_prims_r[ww].gf=Bx_staggerr; out_prims_l[ww].gf=Bx_staggerl; ww++;\n", " in_prims[ww].gf=By_stagger; out_prims_r[ww].gf=By_staggerr; out_prims_l[ww].gf=By_staggerl; ww++;\n", " in_prims[ww].gf=Bz_stagger; out_prims_r[ww].gf=Bz_staggerr; out_prims_l[ww].gf=Bz_staggerl; ww++;\n", " in_prims[ww].gf=vxr; out_prims_r[ww].gf=vxrr; out_prims_l[ww].gf=vxrl; ww++;\n", " in_prims[ww].gf=vyr; out_prims_r[ww].gf=vyrr; out_prims_l[ww].gf=vyrl; ww++;\n", " in_prims[ww].gf=vzr; out_prims_r[ww].gf=vzrr; out_prims_l[ww].gf=vzrl; ww++;\n", " in_prims[ww].gf=vxl; out_prims_r[ww].gf=vxlr; out_prims_l[ww].gf=vxll; ww++;\n", " in_prims[ww].gf=vyl; out_prims_r[ww].gf=vylr; out_prims_l[ww].gf=vyll; ww++;\n", " in_prims[ww].gf=vzl; out_prims_r[ww].gf=vzlr; out_prims_l[ww].gf=vzll; ww++;\n", "\n", " // Prims are defined AT ALL GRIDPOINTS, so we set the # of ghostzones to zero:\n", " for(int i=0;i\n", "\n", "## Step 3.3: Convert ADM variables to BSSN variables \\[Back to [top](#toc)\\]\n", "$$\\label{admbase_to_bssnbase}$$\n", "\n", "We summarize here the algorithm of the `IllinoisGRMHD_convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij()` function, which is explained in detail in [this tutorial module]() (**Link not available yet - TODO**):\n", "\n", "* First, $\\gamma\\equiv \\det\\left(\\gamma_{ij}\\right)$, where $\\gamma_{ij}$ is the physical spatial metric, is evaluated.\n", "* Then, $\\phi$, the conformal factor, is computed via the relation $\\phi = \\frac{1}{12}\\log\\gamma$.\n", "* Next, we compute $e^{-4\\phi}$.\n", "* Then, the conformal metric, $\\bar{\\gamma}_{ij} = e^{-4\\phi}\\gamma_{ij}$, is computed.\n", "* Next, the condition $\\bar\\gamma = 1$ is enforced, by first computing $\\gamma$ then performing $\\bar\\gamma_{ij}\\to\\left(\\frac{1}{\\bar\\gamma}\\right)^{1/3}\\bar\\gamma_{ij}$.\n", "* Then, $\\gamma_{ij}$ is computed from $\\bar\\gamma_{ij}$ *after* the condition $\\bar\\gamma = 1$ is enforced via the inverse relation $\\gamma_{ij} = e^{4\\phi}\\bar\\gamma_{ij}$.\n", "* Finally, we compute the inverse conformal metric $\\bar\\gamma^{ij}$.\n", "\n", "**A note on notation:** in the C code, we have the following identifications between the quantities described above and the C variables:\n", "\n", "*Input and temporary variables:*\n", "* $\\gamma_{ij} := {\\rm gij\\_physL}$\n", "* $\\det(\\gamma_{ij}) := {\\rm gijdet}$\n", "* $\\phi := {\\rm phiL}$\n", "* $\\psi \\equiv e^{\\phi} := {\\rm psiL}$\n", "* $\\bar\\gamma_{ij} := {\\rm gtijL}$ (the \"t\" stands for the notation where the conformal metric is written as $\\tilde\\gamma_{ij}$ instead of $\\bar\\gamma_{ij}$. In our discussion we use the latter to keep our notation consistent with other NRPy notebooks).\n", "* $\\det(\\bar\\gamma_{ij}) := {\\rm gtijdet}$\n", "* $\\left(\\frac{1}{\\bar\\gamma}\\right)^{1/3} := {\\rm gtijdet\\_Fm1o3}$\n", "\n", "*Output/gridfunction variables:*\n", "* $\\gamma_{ij} := {\\rm gij}$ (Physical metric)\n", "* $\\bar\\gamma_{ij} := {\\rm gtij}$ (Conformal metric)\n", "* $\\bar\\gamma^{ij} := {\\rm gtupij}$ (Inverse conformal metric)\n", "* $\\phi := {\\rm phi}$\n", "* $\\psi := {\\rm psi}$" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/driver_evaluate_MHD_rhs.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__driver_evaluate_MHD_rhs__C\n", "\n", "\n", " // Convert ADM variables (from ADMBase) to the BSSN-based variables expected by this routine.\n", " IllinoisGRMHD_convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij(cctkGH,cctk_lsh, gxx,gxy,gxz,gyy,gyz,gzz,alp,\n", " gtxx,gtxy,gtxz,gtyy,gtyz,gtzz,\n", " gtupxx,gtupxy,gtupxz,gtupyy,gtupyz,gtupzz,\n", " phi_bssn,psi_bssn,lapm1);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 3.4: Setting up pointers to the metric and stress-energy tensor gridfunctions \\[Back to [top](#toc)\\]\n", "$$\\label{pointers_metric_tmunu_gfs}$$" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/driver_evaluate_MHD_rhs.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__driver_evaluate_MHD_rhs__C\n", "\n", "\n", " /* SET POINTERS TO METRIC GRIDFUNCTIONS */\n", " CCTK_REAL *metric[NUMVARS_FOR_METRIC_FACEVALS]; // \"metric\" here is array of pointers to the actual gridfunctions.\n", " ww=0;\n", " metric[ww]=phi_bssn;ww++;\n", " metric[ww]=psi_bssn;ww++;\n", " metric[ww]=gtxx; ww++;\n", " metric[ww]=gtxy; ww++;\n", " metric[ww]=gtxz; ww++;\n", " metric[ww]=gtyy; ww++;\n", " metric[ww]=gtyz; ww++;\n", " metric[ww]=gtzz; ww++;\n", " metric[ww]=lapm1; ww++;\n", " metric[ww]=betax; ww++;\n", " metric[ww]=betay; ww++;\n", " metric[ww]=betaz; ww++;\n", " metric[ww]=gtupxx; ww++;\n", " metric[ww]=gtupyy; ww++;\n", " metric[ww]=gtupzz; ww++;\n", "\n", " /* SET POINTERS TO STRESS-ENERGY TENSOR GRIDFUNCTIONS */\n", " CCTK_REAL *TUPmunu[10];// \"TUPmunu\" here is array of pointers to the actual gridfunctions.\n", " ww=0;\n", " TUPmunu[ww]=TUPtt; ww++;\n", " TUPmunu[ww]=TUPtx; ww++;\n", " TUPmunu[ww]=TUPty; ww++;\n", " TUPmunu[ww]=TUPtz; ww++;\n", " TUPmunu[ww]=TUPxx; ww++;\n", " TUPmunu[ww]=TUPxy; ww++;\n", " TUPmunu[ww]=TUPxz; ww++;\n", " TUPmunu[ww]=TUPyy; ww++;\n", " TUPmunu[ww]=TUPyz; ww++;\n", " TUPmunu[ww]=TUPzz; ww++;" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 3.5: Initialization of the RHS variables \\[Back to [top](#toc)\\]\n", "$$\\label{initialize_rhss}$$" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/driver_evaluate_MHD_rhs.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__driver_evaluate_MHD_rhs__C\n", "\n", "\n", " // 1) First initialize {rho_star_rhs,tau_rhs,st_x_rhs,st_y_rhs,st_z_rhs} to zero\n", "#pragma omp parallel for\n", " for(int k=0;k\n", "\n", "## Step 3.6: Compute extrinsic curvature terms from the RHS of $\\partial_{t}\\tilde\\tau$ and $T^{\\mu\\nu}$ \\[Back to [top](#toc)\\]\n", "$$\\label{tau_rhs_ext_curv_and_tupmunu}$$" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/driver_evaluate_MHD_rhs.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__driver_evaluate_MHD_rhs__C\n", "\n", " // Here, we:\n", " // 1) Compute tau_rhs extrinsic curvature terms, and\n", " // 2) Compute TUPmunu.\n", " // This function is housed in the file: \"compute_tau_rhs_extrinsic_curvature_terms_and_TUPmunu.C\"\n", " compute_tau_rhs_extrinsic_curvature_terms_and_TUPmunu(cctkGH,cctk_lsh,cctk_nghostzones,dX,metric,in_prims,TUPmunu,\n", " eos, Gamma_th,\n", " gtupxy,gtupxz,gtupyz,\n", " kxx,kxy,kxz,kyy,kyz,kzz,\n", " tau_rhs);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 3.7: Computing ${\\rm ftilde}$ \\[Back to [top](#toc)\\]\n", "$$\\label{computing_ftilde}$$\n", "\n", "This is part of the flattening scheme of the PPM algorithm. The main reference to look at is [Colella & Woodward (1983)](https://crd.lbl.gov/assets/pubs_presos/AMCS/ANAG/A141984.pdf). The equations implemented can be found in Appendix A (particularly eqs. (A.1) and (A.2)), while the flattening method is introduced and discussed in section 4. More will follow when we talk about the `reconstruct_set_of_prims_PPM.C` file of `IllinoisGRMHD`." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/driver_evaluate_MHD_rhs.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__driver_evaluate_MHD_rhs__C\n", "\n", " int flux_dirn;\n", " flux_dirn=1;\n", " // First compute ftilde, which is used for flattening left and right face values\n", " // This function is housed in the file: \"reconstruct_set_of_prims_PPM.C\"\n", " ftilde_gf_compute(cctkGH,cctk_lsh,flux_dirn, in_prims, ftilde_gf);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 3.8: The RHSs of $\\rho_{\\star}$, $\\tilde\\tau$, $\\tilde{S}_{i}$, and $A_{i}$ \\[Back to [top](#toc)\\]\n", "$$\\label{rhs_mhd_and_a_i}$$\n", "\n", "This part of the code evaluates the RHSs of $\\rho_{\\star}$, $\\tilde\\tau$, and $\\tilde{S}_{i}$, i.e.\n", "\n", "$$\n", "\\partial_{t}\n", "\\begin{bmatrix}\n", "\\rho_{\\star}\\\\\n", "\\tilde\\tau\\\\\n", "\\tilde{S}_{i}\n", "\\end{bmatrix}\n", "=\n", "-\\partial_{j}\n", "\\underbrace{\\begin{bmatrix}\n", "\\rho_{\\star}v^{j}\\\\\n", "\\alpha^{2}\\sqrt{\\gamma}T^{0j} - \\rho_{\\star}v^{j}\\\\\n", "\\alpha\\sqrt{\\gamma}T^{j}_{\\ i}\n", "\\end{bmatrix}}_{\\rm Flux\\ terms}\n", "+\n", "\\underbrace{\\begin{bmatrix}\n", "0\\\\\n", "s\\\\\n", "\\frac{1}{2}\\alpha\\sqrt{\\gamma}T^{\\alpha\\beta}\\partial_{i}g_{\\alpha\\beta}\n", "\\end{bmatrix}}_{\\rm Source\\ terms}\\ .\n", "$$\n", "\n", "At the same time, we are also interested in evaluating the RHS of the evolution equation for $A_{i}$, namely\n", "\n", "$$\n", "\\partial_{t}A_{i} = \\epsilon_{ijk}v^{j}\\tilde{B}^{k} - \\partial_{i}\\left(\\alpha\\Phi - \\beta^{j}A_{j}\\right) = \\psi^{6}\\epsilon_{ijk}v^{j}B^{k} - \\underbrace{\\partial_{i}\\left(\\alpha\\Phi - \\beta^{j}A_{j}\\right)}_{\\rm Gauge\\ terms}\n", "$$\n", "\n", "The following summary greatly oversimplifies what the code below does, but it is enough for the user to understand the purpose of the algorithm:\n", "\n", "1. Compute $\\partial_{x}\\boldsymbol{F}$, then $\\partial_{y}\\boldsymbol{F}$, and finally $\\left[\\partial_{t}A_{z}\\right]_{\\rm no\\ gauge\\ terms}$\n", "2. Compute $\\partial_{y}\\boldsymbol{F}$, then $\\left[\\partial_{t}A_{x}\\right]_{\\rm no\\ gauge\\ terms}$\n", "3. Compute $\\left[\\partial_{t}A_{y}\\right]_{\\rm no\\ gauge\\ terms}$\n", "4. Add gauge terms to $\\partial_{t}A_{i}$\n", "\n", "Now, in between every step of the summary above, care must be taken to evaluate the gridfunctions at the appropriate gridpoints (see table below for the location of each variable in the computational grid).\n", "\n", "\n", "\n", "| Variable(s) | Gridpoint location in the computational grid |\n", "|------------------------------------------------------------------|-----------------------------------------------|\n", "| Metric terms, $\\vec{P}$, $\\rho_*$, $\\tilde{S}_i$, $\\tilde{\\tau}$ | $(i,j,k)$ |\n", "| $B^x$, $\\tilde{B}^x$ | $(i+\\frac{1}{2},j,k)$ |\n", "| $B^y$, $\\tilde{B}^y$ | $(i,j+\\frac{1}{2},k)$ |\n", "| $B^z$, $\\tilde{B}^z$ | $(i,j,k+\\frac{1}{2})$ |\n", "| $A_x$ | $(i,j+\\frac{1}{2},k+\\frac{1}{2})$ |\n", "| $A_y$ | $(i+\\frac{1}{2},j,k+\\frac{1}{2})$ |\n", "| $A_z$ | $(i+\\frac{1}{2},j+\\frac{1}{2},k)$ |\n", "| $\\sqrt{\\gamma}\\Phi$ | $(i+\\frac{1}{2},j+\\frac{1}{2},k+\\frac{1}{2})$ |\n", "\n", "$$\\label{table_staggerings}$$\n", "\n", "For example, we know that\n", "\n", "$$\n", "\\left[\\partial_{t}A_{z}\\right]_{\\rm no\\ gauge\\ terms} \\equiv \\psi^{6}\\left(v^{x}B^{y} - v^{y}B^{x}\\right)\\ .\n", "$$\n", "\n", "But once we evaluate $v^{x}$ and $v^{y}$, we know them at the point $(i,j,k)$. Similarly, the gridfunction $B^{x}$ is known at $(i+\\frac{1}{2},j,k)$, while $B^{y}$ is known at $(i,j+\\frac{1}{2},k)$. This means that we are not able to immediately evaluate the equation above, since determining $A_{z}$ at $(i+\\frac{1}{2},j+\\frac{1}{2},k)$ requires knowning $\\left\\{v^{x},v^{y},B^{x},B^{y}\\right\\}$ at $(i+\\frac{1}{2},j+\\frac{1}{2},k)$ as well. To this end, we reconstruct the variables $\\left\\{v^{x},v^{y},B^{x},B^{y}\\right\\}$ using the PPM method at the desired staggered point. An analogous procedure is required in order to determine the RHS of $\\partial_{t}A_{x}$ and $\\partial_{t}A_{y}$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "### Step 3.8.1: Reconstructing $\\left\\{v^{x}, v^{y}, B^{y}_{\\rm stagger}\\right\\}$ along the $x$-direction \\[Back to [top](#toc)\\]\n", "$$\\label{reconstructing_vx_vy_by_along_x}$$\n", "\n", "We want to evaluate $\\partial_{x}\\boldsymbol{F}$. It is important that we keep in the back of our minds our intention of evaluating the RHS of $\\left[\\partial_{t}A_{z}\\right]_{\\rm no\\ gauge\\ terms}$ as well, since then we can reconstruct $\\left\\{v^{x},v^{y},B^{x},B^{y}\\right\\}$ cleverly, as we need them at the same gridpoint as $A_{z}$ (see the table in the end of [step 3.8](#table_staggerings)).\n", "\n", "We start by reconstructing $\\left\\{\\rho_{0},P,v^{i},B^{i}, B^{y}_{\\rm stagger}\\right\\}$ in the $x$-direction, keeping in mind that after the reconstruction we will know:\n", "\n", "1. The flux variables at $\\left(i-\\frac{1}{2},j,k\\right)$, so that we can evaluate $\\partial_{x}\\boldsymbol{F}_{i,j,k}=dx^{-1}\\left(\\boldsymbol{F}_{i+1/2,j,k}-\\boldsymbol{F}_{i-1/2,j,k}\\right)$\n", "1. The velocities $\\left\\{v^{x},v^{y}\\right\\}$ at $\\left(i-\\frac{1}{2},j,k\\right)$\n", "1. The staggered value $B^{y}_{\\rm stagger}$ at $\\left(i-\\frac{1}{2},j+\\frac{1}{2},k\\right)$" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/driver_evaluate_MHD_rhs.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__driver_evaluate_MHD_rhs__C\n", "\n", "\n", "\n", " /* There are two stories going on here:\n", " * 1) Computation of \\partial_x on RHS of \\partial_t {rho_star,tau,mhd_st_{x,y,z}},\n", " * via PPM reconstruction onto (i-1/2,j,k), so that\n", " * \\partial_x F = [ F(i+1/2,j,k) - F(i-1/2,j,k) ] / dx\n", " * 2) Computation of \\partial_t A_i, where A_i are *staggered* gridfunctions,\n", " * where A_x is defined at (i,j+1/2,k+1/2), A_y at (i+1/2,j,k+1/2), etc.\n", " * Ai_rhs = \\partial_t A_i = \\epsilon_{ijk} \\psi^{6} v^j B^k,\n", " * where \\epsilon_{ijk} is the flat-space antisymmetric operator.\n", " * 2A) Az_rhs is defined at (i+1/2,j+1/2,k), and it depends on {Bx,By,vx,vy},\n", " * so the trick is to reconstruct {Bx,By,vx,vy} cleverly to get to these\n", " * staggered points. For example:\n", " * 2Aa) vx and vy are at (i,j,k), and we reconstruct them to (i-1/2,j,k) below. After\n", " * this, we'll reconstruct again in the y-dir'n to get {vx,vy} at (i-1/2,j-1/2,k)\n", " * 2Ab) By_stagger is at (i,j+1/2,k), and we reconstruct below to (i-1/2,j+1/2,k). */\n", " ww=0;\n", " which_prims_to_reconstruct[ww]=RHOB; ww++;\n", " which_prims_to_reconstruct[ww]=PRESSURE; ww++;\n", " which_prims_to_reconstruct[ww]=VX; ww++;\n", " which_prims_to_reconstruct[ww]=VY; ww++;\n", " which_prims_to_reconstruct[ww]=VZ; ww++;\n", " //which_prims_to_reconstruct[ww]=BX_CENTER; ww++;\n", " which_prims_to_reconstruct[ww]=BY_CENTER; ww++;\n", " which_prims_to_reconstruct[ww]=BZ_CENTER; ww++;\n", " which_prims_to_reconstruct[ww]=BY_STAGGER;ww++;\n", " num_prims_to_reconstruct=ww;\n", " // This function is housed in the file: \"reconstruct_set_of_prims_PPM.C\"\n", " reconstruct_set_of_prims_PPM(cctkGH,cctk_lsh,flux_dirn,num_prims_to_reconstruct,which_prims_to_reconstruct,\n", " eos,in_prims,out_prims_r,out_prims_l,ftilde_gf,temporary);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "### Step 3.8.2: Evaluating $\\partial_{x}\\boldsymbol{F}$ \\[Back to [top](#toc)\\]\n", "$$\\label{fluxes_x_dirn}$$\n", "\n", "Next we set the face values of $B^{x}$ (which are needed for the computation of MHD flux terms) by making them consistent with $B^{x}_{\\rm stagger}$.\n", "\n", "After that, we evaluate $\\partial_{x}\\boldsymbol{F}$ and add it to the RHS of $\\partial_{t}\\left(\\rho_{\\star},\\tilde\\tau,\\tilde{S}_{i}\\right)$. It is important to notice that, as we mentioned, $A_{z}$ is defined at $\\left(i+\\frac{1}{2},j+\\frac{1}{2},k\\right)$, but other functions, like $v^{x}$ and $v^{y}$, are now known only at $\\left(i-\\frac{1}{2},j-\\frac{1}{2},k\\right)$. The function `add_fluxes_and_source_terms_to_hydro_rhss()` below takes care of this, and we will study the process in more detail when we look at the `add_fluxes_and_source_terms_to_hydro_rhss.C` file of `IllinoisGRMHD`." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/driver_evaluate_MHD_rhs.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__driver_evaluate_MHD_rhs__C\n", "\n", " //Right and left face values of BI_CENTER are used in mhdflux computation (first to compute b^a).\n", " // Instead of reconstructing, we simply set B^x face values to be consistent with BX_STAGGER.\n", "#pragma omp parallel for\n", " for(int k=0;k\n", "\n", "### Step 3.8.3: Reconstructing $\\left\\{v^{x}, v^{y}, B^{y}_{\\rm stagger}\\right\\}$ along the $y$-direction \\[Back to [top](#toc)\\]\n", "$$\\label{reconstructing_vx_vy_by_along_y}$$\n", "\n", "We want to evaluate $\\partial_{y}\\boldsymbol{F}$. At this point we must remember that $v^{x}$ and $v^{y}$ have already been reconstruct along the $x$-direction and are now known at $\\left(i-\\frac{1}{2},j,k\\right)$. Our goal is to reconstruct these quantities at $\\left(i+\\frac{1}{2},j+\\frac{1}{2},k\\right)$.\n", "\n", "We then reconstruct $\\left\\{\\rho_{0},P,v^{i},B^{i}, B^{i}_{\\rm stagger}\\right\\}$ in the $y$-direction, keeping in mind that after the reconstruction we will know:\n", "\n", "1. The flux variables at $\\left(i-\\frac{1}{2},j-\\frac{1}{2},k\\right)$, so that we can evaluate $\\partial_{y}\\boldsymbol{F}_{i,j,k}=dy^{-1}\\left(\\boldsymbol{F}_{i,j+1/2,k}-\\boldsymbol{F}_{i,j-1/2,k}\\right)$\n", "1. The velocities $\\left\\{v^{x},v^{y}\\right\\}$ at $\\left(i-\\frac{1}{2},j-\\frac{1}{2},k\\right)$\n", "1. The staggered value $B^{x}_{\\rm stagger}$ at $\\left(i+\\frac{1}{2},j-\\frac{1}{2},k\\right)$\n", "1. The staggered value $B^{y}_{\\rm stagger}$ at $\\left(i-\\frac{1}{2},j+\\frac{1}{2},k\\right)$\n", "1. The staggered value $B^{z}_{\\rm stagger}$ at $\\left(i,j-\\frac{1}{2},k+\\frac{1}{2}\\right)$" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/driver_evaluate_MHD_rhs.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__driver_evaluate_MHD_rhs__C\n", "\n", "\n", " // Note that we have already reconstructed vx and vy along the x-direction,\n", " // at (i-1/2,j,k). That result is stored in v{x,y}{r,l}. Bx_stagger data\n", " // are defined at (i+1/2,j,k).\n", " // Next goal: reconstruct Bx, vx and vy at (i+1/2,j+1/2,k).\n", " flux_dirn=2;\n", " // First compute ftilde, which is used for flattening left and right face values\n", " // This function is housed in the file: \"reconstruct_set_of_prims_PPM.C\"\n", " ftilde_gf_compute(cctkGH,cctk_lsh,flux_dirn, in_prims, ftilde_gf);\n", "\n", " // in_prims[{VXR,VXL,VYR,VYL}].gz_{lo,hi} ghostzones are set to all zeros, which\n", " // is incorrect. We fix this below.\n", " // [Note that this is a cheap operation, copying only 8 integers and a pointer.]\n", " in_prims[VXR]=out_prims_r[VX];\n", " in_prims[VXL]=out_prims_l[VX];\n", " in_prims[VYR]=out_prims_r[VY];\n", " in_prims[VYL]=out_prims_l[VY];\n", "\n", " /* There are two stories going on here:\n", " * 1) Computation of \\partial_y on RHS of \\partial_t {rho_star,tau,mhd_st_{x,y,z}},\n", " * via PPM reconstruction onto (i,j-1/2,k), so that\n", " * \\partial_y F = [ F(i,j+1/2,k) - F(i,j-1/2,k) ] / dy\n", " * 2) Computation of \\partial_t A_i, where A_i are *staggered* gridfunctions,\n", " * where A_x is defined at (i,j+1/2,k+1/2), A_y at (i+1/2,j,k+1/2), etc.\n", " * Ai_rhs = \\partial_t A_i = \\epsilon_{ijk} \\psi^{6} v^j B^k,\n", " * where \\epsilon_{ijk} is the flat-space antisymmetric operator.\n", " * 2A) Az_rhs is defined at (i+1/2,j+1/2,k), and it depends on {Bx,By,vx,vy},\n", " * so the trick is to reconstruct {Bx,By,vx,vy} cleverly to get to these\n", " * staggered points. For example:\n", " * 2Aa) VXR = [right-face of vx reconstructed along x-direction above] is at (i-1/2,j,k),\n", " * and we reconstruct it to (i-1/2,j-1/2,k) below. Similarly for {VXL,VYR,VYL}\n", " * 2Ab) Bx_stagger is at (i+1/2,j,k), and we reconstruct to (i+1/2,j-1/2,k) below\n", " * 2Ac) By_stagger is at (i-1/2,j+1/2,k) already for Az_rhs, from the previous step.\n", " * 2B) Ax_rhs is defined at (i,j+1/2,k+1/2), and it depends on {By,Bz,vy,vz}.\n", " * Again the trick is to reconstruct these onto these staggered points.\n", " * 2Ba) Bz_stagger is at (i,j,k+1/2), and we reconstruct to (i,j-1/2,k+1/2) below */\n", " ww=0;\n", " // NOTE! The order of variable reconstruction is important here,\n", " // as we don't want to overwrite {vxr,vxl,vyr,vyl}!\n", " which_prims_to_reconstruct[ww]=VXR; ww++;\n", " which_prims_to_reconstruct[ww]=VYR; ww++;\n", " which_prims_to_reconstruct[ww]=VXL; ww++;\n", " which_prims_to_reconstruct[ww]=VYL; ww++;\n", " num_prims_to_reconstruct=ww;\n", " // This function is housed in the file: \"reconstruct_set_of_prims_PPM.C\"\n", " reconstruct_set_of_prims_PPM(cctkGH,cctk_lsh,flux_dirn,num_prims_to_reconstruct,which_prims_to_reconstruct,\n", " eos,in_prims,out_prims_r,out_prims_l,ftilde_gf,temporary);\n", " ww=0;\n", " // Reconstruct other primitives last!\n", " which_prims_to_reconstruct[ww]=RHOB; ww++;\n", " which_prims_to_reconstruct[ww]=PRESSURE; ww++;\n", " which_prims_to_reconstruct[ww]=VX; ww++;\n", " which_prims_to_reconstruct[ww]=VY; ww++;\n", " which_prims_to_reconstruct[ww]=VZ; ww++;\n", " which_prims_to_reconstruct[ww]=BX_CENTER; ww++;\n", " //which_prims_to_reconstruct[ww]=BY_CENTER; ww++;\n", " which_prims_to_reconstruct[ww]=BZ_CENTER; ww++;\n", " which_prims_to_reconstruct[ww]=BX_STAGGER;ww++;\n", " which_prims_to_reconstruct[ww]=BZ_STAGGER;ww++;\n", " num_prims_to_reconstruct=ww;\n", " // This function is housed in the file: \"reconstruct_set_of_prims_PPM.C\"\n", " reconstruct_set_of_prims_PPM(cctkGH,cctk_lsh,flux_dirn,num_prims_to_reconstruct,which_prims_to_reconstruct,\n", " eos,in_prims,out_prims_r,out_prims_l,ftilde_gf,temporary);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "### Step 3.8.4: Evaluating $\\partial_{y}\\boldsymbol{F}$ \\[Back to [top](#toc)\\]\n", "$$\\label{fluxes_y_dirn}$$" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/driver_evaluate_MHD_rhs.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__driver_evaluate_MHD_rhs__C\n", "\n", " //Right and left face values of BI_CENTER are used in mhdflux computation (first to compute b^a).\n", " // Instead of reconstructing, we simply set B^y face values to be consistent with BY_STAGGER.\n", "#pragma omp parallel for\n", " for(int k=0;k\n", "\n", "### Step 3.8.5: Evaluating $\\left[\\partial_{t}A_{z}\\right]_{\\rm no\\ gauge\\ terms}$ \\[Back to [top](#toc)\\]\n", "$$\\label{rhs_az_no_gauge_terms}$$\n", "\n", "As a friendly reminder, we summarize the known gridpoint location of the needed gridfunctions here:\n", "\n", "| Staggered and unstaggered variables | Gridpoint location at which the variable is known |\n", "|-----------------------------------------------------|----------------------------------------------------|\n", "| $\\left(v^{x}\\right)_{r,l},\\left(v^{y}\\right)_{r,l}$ | $(i-\\frac{1}{2},j-\\frac{1}{2},k)$ |\n", "| $\\left(B^{z}_{\\rm stagger}\\right)_{r,l}$ | $(i+\\frac{1}{2},j-\\frac{1}{2},k)$ |\n", "| $\\left(B^{y}_{\\rm stagger}\\right)_{r,l}$ | $(i-\\frac{1}{2},j+\\frac{1}{2},k)$ |\n", "| $\\phi$ | $(i,j,k)$ |\n", "\n", "We start by interpolating $\\phi$ to $\\left(i+\\frac{1}{2},j,k\\right)$, followed by a second interpolation so that $\\phi$ is known at $\\left(i+\\frac{1}{2},j+\\frac{1}{2},k\\right)$." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/driver_evaluate_MHD_rhs.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__driver_evaluate_MHD_rhs__C\n", "\n", "\n", " /*****************************************\n", " * COMPUTING RHS OF A_z, BOOKKEEPING NOTE:\n", " * We want to compute\n", " * \\partial_t A_z - [gauge terms] = \\psi^{6} (v^x B^y - v^y B^x).\n", " * A_z is defined at (i+1/2,j+1/2,k).\n", " * ==========================\n", " * Where defined | Variables\n", " * (i-1/2,j-1/2,k)| {vxrr,vxrl,vxlr,vxll,vyrr,vyrl,vylr,vyll}\n", " * (i+1/2,j-1/2,k)| {Bx_stagger_r,Bx_stagger_l} (see Table 1 in arXiv:1007.2848)\n", " * (i-1/2,j+1/2,k)| {By_stagger_r,By_stagger_l} (see Table 1 in arXiv:1007.2848)\n", " * (i,j,k) | {phi}\n", " * ==========================\n", " ******************************************/\n", " // Interpolates to i+1/2\n", "#define IPH(METRICm1,METRICp0,METRICp1,METRICp2) (-0.0625*((METRICm1) + (METRICp2)) + 0.5625*((METRICp0) + (METRICp1)))\n", " // Next compute phi at (i+1/2,j+1/2,k):\n", "#pragma omp parallel for\n", " for(int k=0;k\n", "\n", "### Step 3.8.6: Multiple reconstructions \\[Back to [top](#toc)\\]\n", "$$\\label{multiple_reconstructions}$$\n", "\n", "1. $\\left\\{\\rho_{\\star},P,v^{i}, B^{i}\\right\\}$ from $\\left(i,j,k\\right)$ to $\\left(i,j,k-\\frac{1}{2}\\right)$\n", "1. $\\left[\\partial_{t}A_{x}\\right]_{\\rm no\\ gauge\\ terms}$ is defined at $\\left(i,j+\\frac{1}{2},k+\\frac{1}{2}\\right)$\n", " 1. $\\left(v^{y}\\right)_{r,l}$ and $\\left(v^{z}\\right)_{r,l}$ are at $\\left(i,j-\\frac{1}{2},k\\right)$, so we reconstruct them to $\\left(i,j-\\frac{1}{2},k-\\frac{1}{2}\\right)$\n", " 1. $\\left(B^{z}_{\\rm stagger}\\right)_{r,l}$ is already known at $\\left(i,j-\\frac{1}{2},k+\\frac{1}{2}\\right)$\n", " 1. $B^{y}_{\\rm stagger}$ is at $\\left(i,j+\\frac{1}{2},k\\right)$, so we reconstruct it to $\\left(i,j+\\frac{1}{2},k-\\frac{1}{2}\\right)$\n", "1. $\\left[\\partial_{t}A_{y}\\right]_{\\rm no\\ gauge\\ terms}$ is defined at $\\left(i+\\frac{1}{2},j,k+\\frac{1}{2}\\right)$\n", " 1. $v^{x}$ and $v^{z}$ are reconstructed to $\\left(i,j,k-\\frac{1}{2}\\right)$. We'll reconstruct them to $\\left(i,j-\\frac{1}{2},k-\\frac{1}{2}\\right)$ later\n", " 1. $B^{z}_{\\rm stagger}$ is already known at $\\left(i,j,k+\\frac{1}{2}\\right)$. We'll reconstruct them to $\\left(i-\\frac{1}{2},j-\\frac{1}{2},k+\\frac{1}{2}\\right)$ later\n", " 1. $B^{x}_{\\rm stagger}$ is at $\\left(i,j+\\frac{1}{2},k\\right)$, so we reconstruct it to $\\left(i,j+\\frac{1}{2},k-\\frac{1}{2}\\right)$" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/driver_evaluate_MHD_rhs.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__driver_evaluate_MHD_rhs__C\n", "\n", "\n", " // in_prims[{VYR,VYL,VZR,VZL}].gz_{lo,hi} ghostzones are not correct, so we fix\n", " // this below.\n", " // [Note that this is a cheap operation, copying only 8 integers and a pointer.]\n", " in_prims[VYR]=out_prims_r[VY];\n", " in_prims[VYL]=out_prims_l[VY];\n", " in_prims[VZR]=out_prims_r[VZ];\n", " in_prims[VZL]=out_prims_l[VZ];\n", "\n", " flux_dirn=3;\n", " // First compute ftilde, which is used for flattening left and right face values\n", " // This function is housed in the file: \"reconstruct_set_of_prims_PPM.C\"\n", " ftilde_gf_compute(cctkGH,cctk_lsh,flux_dirn, in_prims, ftilde_gf);\n", "\n", " /* There are two stories going on here:\n", " * 1) Single reconstruction to (i,j,k-1/2) for {rho,P,vx,vy,vz,Bx,By,Bz} to compute\n", " * z-dir'n advection terms in \\partial_t {rho_star,tau,mhd_st_{x,y,z}} at (i,j,k)\n", " * 2) Multiple reconstructions for *staggered* gridfunctions A_i:\n", " * Ai_rhs = \\partial_t A_i = \\epsilon_{ijk} \\psi^{6} v^j B^k,\n", " * where \\epsilon_{ijk} is the flat-space antisymmetric operator.\n", " * 2A) Ax_rhs is defined at (i,j+1/2,k+1/2), depends on v{y,z} and B{y,z}\n", " * 2Aa) v{y,z}{r,l} are at (i,j-1/2,k), so we reconstruct here to (i,j-1/2,k-1/2)\n", " * 2Ab) Bz_stagger{r,l} are at (i,j-1/2,k+1/2) already.\n", " * 2Ac) By_stagger is at (i,j+1/2,k), and below we reconstruct its value at (i,j+1/2,k-1/2)\n", " * 2B) Ay_rhs is defined at (i+1/2,j,k+1/2), depends on v{z,x} and B{z,x}.\n", " * 2Ba) v{x,z} are reconstructed to (i,j,k-1/2). Later we'll reconstruct again to (i-1/2,j,k-1/2).\n", " * 2Bb) Bz_stagger is at (i,j,k+1/2). Later we will reconstruct to (i-1/2,j,k+1/2).\n", " * 2Bc) Bx_stagger is at (i+1/2,j,k), and below we reconstruct its value at (i+1/2,j,k-1/2)\n", " */\n", " ww=0;\n", " // NOTE! The order of variable reconstruction is important here,\n", " // as we don't want to overwrite {vxr,vxl,vyr,vyl}!\n", " which_prims_to_reconstruct[ww]=VYR; ww++;\n", " which_prims_to_reconstruct[ww]=VZR; ww++;\n", " which_prims_to_reconstruct[ww]=VYL; ww++;\n", " which_prims_to_reconstruct[ww]=VZL; ww++;\n", " num_prims_to_reconstruct=ww;\n", " // This function is housed in the file: \"reconstruct_set_of_prims_PPM.C\"\n", " reconstruct_set_of_prims_PPM(cctkGH,cctk_lsh,flux_dirn,num_prims_to_reconstruct,which_prims_to_reconstruct,\n", " eos,in_prims,out_prims_r,out_prims_l,ftilde_gf,temporary);\n", " // Reconstruct other primitives last!\n", " ww=0;\n", " which_prims_to_reconstruct[ww]=RHOB; ww++;\n", " which_prims_to_reconstruct[ww]=PRESSURE; ww++;\n", " which_prims_to_reconstruct[ww]=VX; ww++;\n", " which_prims_to_reconstruct[ww]=VY; ww++;\n", " which_prims_to_reconstruct[ww]=VZ; ww++;\n", " which_prims_to_reconstruct[ww]=BX_CENTER; ww++;\n", " which_prims_to_reconstruct[ww]=BY_CENTER; ww++;\n", " //which_prims_to_reconstruct[ww]=BZ_CENTER; ww++;\n", " which_prims_to_reconstruct[ww]=BX_STAGGER; ww++;\n", " which_prims_to_reconstruct[ww]=BY_STAGGER; ww++;\n", " num_prims_to_reconstruct=ww;\n", " // This function is housed in the file: \"reconstruct_set_of_prims_PPM.C\"\n", " reconstruct_set_of_prims_PPM(cctkGH,cctk_lsh,flux_dirn,num_prims_to_reconstruct,which_prims_to_reconstruct,\n", " eos,in_prims,out_prims_r,out_prims_l,ftilde_gf,temporary);\n", " //Right and left face values of BI_CENTER are used in mhdflux computation (first to compute b^a).\n", " // Instead of reconstructing, we simply set B^z face values to be consistent with BZ_STAGGER.\n", "#pragma omp parallel for\n", " for(int k=0;k\n", "\n", "### Step 3.8.7: Evaluating $\\partial_{z}\\boldsymbol{F}$ \\[Back to [top](#toc)\\]\n", "$$\\label{fluxes_z_dirn}$$" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/driver_evaluate_MHD_rhs.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__driver_evaluate_MHD_rhs__C\n", "\n", " // Then add fluxes to RHS for hydro variables {rho_b,P,vx,vy,vz}:\n", " // This function is housed in the file: \"add_fluxes_and_source_terms_to_hydro_rhss.C\"\n", " add_fluxes_and_source_terms_to_hydro_rhss(flux_dirn,cctkGH,cctk_lsh,cctk_nghostzones,dX, metric,in_prims,TUPmunu,\n", " num_prims_to_reconstruct,out_prims_r,out_prims_l,eos,\n", " cmax_z,cmin_z,\n", " rho_star_flux,tau_flux,st_x_flux,st_y_flux,st_z_flux,\n", " rho_star_rhs,tau_rhs,st_x_rhs,st_y_rhs,st_z_rhs);\n", "\n", " // in_prims[{VYR,VYL,VZR,VZL}].gz_{lo,hi} ghostzones are not set correcty.\n", " // We fix this below.\n", " // [Note that this is a cheap operation, copying only 8 integers and a pointer.]\n", " in_prims[VXR]=out_prims_r[VX];\n", " in_prims[VZR]=out_prims_r[VZ];\n", " in_prims[VXL]=out_prims_l[VX];\n", " in_prims[VZL]=out_prims_l[VZ];\n", " // FIXME: lines above seem to be inconsistent with lines below.... Possible bug, not major enough to affect evolutions though.\n", " in_prims[VZR].gz_lo[1]=in_prims[VZR].gz_hi[1]=0;\n", " in_prims[VXR].gz_lo[1]=in_prims[VXR].gz_hi[1]=0;\n", " in_prims[VZL].gz_lo[1]=in_prims[VZL].gz_hi[1]=0;\n", " in_prims[VXL].gz_lo[1]=in_prims[VXL].gz_hi[1]=0;" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "### Step 3.8.8: Evaluating $\\left[\\partial_{t}A_{x}\\right]_{\\rm no\\ gauge\\ terms}$ \\[Back to [top](#toc)\\]\n", "$$\\label{rhs_ax_no_gauge_terms}$$\n", "\n", "As a friendly reminder, we summarize the known gridpoint location of the needed gridfunctions here:\n", "\n", "| Staggered and unstaggered variables | Gridpoint location at which the variable is known |\n", "|-----------------------------------------------------|----------------------------------------------------|\n", "| $\\left(v^{y}\\right)_{r,l},\\left(v^{z}\\right)_{r,l}$ | $(i,j-\\frac{1}{2},k-\\frac{1}{2})$ |\n", "| $\\left(B^{y}_{\\rm stagger}\\right)_{r,l}$ | $(i,j+\\frac{1}{2},k-\\frac{1}{2})$ |\n", "| $\\left(B^{z}_{\\rm stagger}\\right)_{r,l}$ | $(i,j-\\frac{1}{2},k+\\frac{1}{2})$ |\n", "| $\\phi$ | $(i,j,k)$ |\n", "\n", "We start by interpolating $\\phi$ to $\\left(i,j+\\frac{1}{2},k+\\frac{1}{2}\\right)$." ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/driver_evaluate_MHD_rhs.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__driver_evaluate_MHD_rhs__C\n", "\n", " /*****************************************\n", " * COMPUTING RHS OF A_x, BOOKKEEPING NOTE:\n", " * We want to compute\n", " * \\partial_t A_x - [gauge terms] = \\psi^{6} (v^y B^z - v^z B^y).\n", " * A_x is defined at (i,j+1/2,k+1/2).\n", " * ==========================\n", " * Where defined | Variables\n", " * (i,j-1/2,k-1/2)| {vyrr,vyrl,vylr,vyll,vzrr,vzrl,vzlr,vzll}\n", " * (i,j+1/2,k-1/2)| {By_stagger_r,By_stagger_l} (see Table 1 in arXiv:1007.2848)\n", " * (i,j-1/2,k+1/2)| {Bz_stagger_r,Bz_stagger_l} (see Table 1 in arXiv:1007.2848)\n", " * (i,j,k) | {phi}\n", " * ==========================\n", " ******************************************/\n", " // Next compute phi at (i,j+1/2,k+1/2):\n", "#pragma omp parallel for\n", " for(int k=1;k\n", "\n", "### Step 3.8.9: Evaluating $\\left[\\partial_{t}A_{y}\\right]_{\\rm no\\ gauge\\ terms}$ \\[Back to [top](#toc)\\]\n", "$$\\label{rhs_ay_no_gauge_terms}$$\n", "\n", "As a friendly reminder, we summarize the known gridpoint location of the needed gridfunctions here:\n", "\n", "| Staggered and unstaggered variables | Gridpoint location at which the variable is known |\n", "|-----------------------------------------------------|----------------------------------------------------|\n", "| $\\left(v^{x}\\right)_{r,l},\\left(v^{z}\\right)_{r,l}$ | $(i-\\frac{1}{2},j,k-\\frac{1}{2})$ |\n", "| $\\left(B^{x}_{\\rm stagger}\\right)_{r,l}$ | $(i+\\frac{1}{2},j,k-\\frac{1}{2})$ |\n", "| $\\left(B^{z}_{\\rm stagger}\\right)_{r,l}$ | $(i-\\frac{1}{2},j,k+\\frac{1}{2})$ |\n", "| $\\phi$ | $(i,j,k)$ |\n", "\n", "We start by interpolating $\\phi$ to $\\left(i+\\frac{1}{2},j,k+\\frac{1}{2}\\right)$." ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/driver_evaluate_MHD_rhs.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__driver_evaluate_MHD_rhs__C\n", "\n", "\n", " // We reprise flux_dirn=1 to finish up computations of Ai_rhs's!\n", " flux_dirn=1;\n", " // First compute ftilde, which is used for flattening left and right face values\n", " // This function is housed in the file: \"reconstruct_set_of_prims_PPM.C\"\n", " ftilde_gf_compute(cctkGH,cctk_lsh,flux_dirn, in_prims, ftilde_gf);\n", "\n", " ww=0;\n", " // NOTE! The order of variable reconstruction is important here,\n", " // as we don't want to overwrite {vxr,vxl,vyr,vyl}!\n", " which_prims_to_reconstruct[ww]=VXR; ww++;\n", " which_prims_to_reconstruct[ww]=VZR; ww++;\n", " which_prims_to_reconstruct[ww]=VXL; ww++;\n", " which_prims_to_reconstruct[ww]=VZL; ww++;\n", " which_prims_to_reconstruct[ww]=BZ_STAGGER;ww++;\n", " num_prims_to_reconstruct=ww;\n", " // This function is housed in the file: \"reconstruct_set_of_prims_PPM.C\"\n", " reconstruct_set_of_prims_PPM(cctkGH,cctk_lsh,flux_dirn,num_prims_to_reconstruct,which_prims_to_reconstruct,\n", " eos,in_prims,out_prims_r,out_prims_l,ftilde_gf,temporary);\n", "\n", " /*****************************************\n", " * COMPUTING RHS OF A_y, BOOKKEEPING NOTE:\n", " * We want to compute\n", " * \\partial_t A_y - [gauge terms] = \\psi^{6} (v^z B^x - v^x B^z).\n", " * A_y is defined at (i+1/2,j,k+1/2).\n", " * ==========================\n", " * Where defined | Variables\n", " * (i-1/2,j,k-1/2)| {vxrr,vxrl,vxlr,vxll,vzrr,vzrl,vzlr,vzll}\n", " * (i+1/2,j,k-1/2)| {Bx_stagger_r,Bx_stagger_l} (see Table 1 in arXiv:1007.2848)\n", " * (i-1/2,j,k+1/2)| {Bz_stagger_r,Bz_stagger_l} (see Table 1 in arXiv:1007.2848)\n", " * (i,j,k) | {phi}\n", " * ==========================\n", " ******************************************/\n", " // Next compute phi at (i+1/2,j,k+1/2):\n", "#pragma omp parallel for\n", " for(int k=1;k\n", "\n", "### Step 3.8.10: Evaluating $\\partial_{t}\\left[\\psi^{6}\\Phi\\right]$ and $\\left[\\partial_{t}A_{i}\\right]_{\\rm gauge\\ terms}$ \\[Back to [top](#toc)\\]\n", "$$\\label{rhs_psi6phi_and_ai_gauge_terms}$$\n", "\n", "Finally, we compute\n", "\n", "$$\n", "\\partial_{t}\\left[\\sqrt{\\gamma}\\Phi\\right] =\n", "\\partial_{t}\\left[\\psi^{6}\\Phi\\right] = \n", "-\\partial_{j}\\left(\\alpha\\sqrt{\\gamma}A^{j} - \\beta^{j}\\left[\\sqrt{\\gamma}\\Phi\\right]\\right)\n", "-\\xi\\alpha\\left[\\sqrt{\\gamma}\\Phi\\right]\\ ,\n", "$$\n", "\n", "and\n", "\n", "$$\n", "\\left[\\partial_{t}A_{i}\\right]_{\\rm gauge\\ terms} = -\\partial_{i}\\left(\\alpha\\Phi - \\beta^{j}A_{j}\\right)\\ .\n", "$$\n", "\n", "Notice that we will need $A^{i}$ to compute $\\partial_{t}\\left[\\psi^{6}\\Phi\\right]$, but we only have $A_{i}$, so we need to determine $\\bar\\gamma^{ij}$ (${\\rm gtupij}$)." ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/driver_evaluate_MHD_rhs.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__driver_evaluate_MHD_rhs__C\n", "\n", "\n", " // Next compute psi6phi_rhs, and add gauge terms to A_i_rhs terms!\n", " // Note that in the following function, we don't bother with reconstruction, instead interpolating.\n", " // We need A^i, but only have A_i. So we add gtupij to the list of input variables.\n", " CCTK_REAL *interp_vars[MAXNUMINTERP];\n", " ww=0;\n", " interp_vars[ww]=betax; ww++;\n", " interp_vars[ww]=betay; ww++;\n", " interp_vars[ww]=betaz; ww++;\n", " interp_vars[ww]=gtupxx; ww++;\n", " interp_vars[ww]=gtupxy; ww++;\n", " interp_vars[ww]=gtupxz; ww++;\n", " interp_vars[ww]=gtupyy; ww++;\n", " interp_vars[ww]=gtupyz; ww++;\n", " interp_vars[ww]=gtupzz; ww++;\n", " interp_vars[ww]=psi_bssn;ww++;\n", " interp_vars[ww]=lapm1; ww++;\n", " interp_vars[ww]=Ax; ww++;\n", " interp_vars[ww]=Ay; ww++;\n", " interp_vars[ww]=Az; ww++;\n", " int max_num_interp_variables=ww;\n", " if(max_num_interp_variables>MAXNUMINTERP) {CCTK_VError(VERR_DEF_PARAMS,\"Error: Didn't allocate enough space for interp_vars[].\"); }\n", " // We are FINISHED with v{x,y,z}{r,l} and P{r,l} so we use these 8 gridfunctions' worth of space as temp storage.\n", " Lorenz_psi6phi_rhs__add_gauge_terms_to_A_i_rhs(cctkGH,cctk_lsh,cctk_nghostzones,dX,interp_vars,psi6phi,\n", " vxr,vyr,vzr,vxl,vyl,vzl,Pr,Pl,\n", " psi6phi_rhs,Ax_rhs,Ay_rhs,Az_rhs);\n", "\n", "\n", " return;\n", " /*\n", " // FUN DEBUGGING TOOL (trust me!):\n", " #pragma omp parallel for\n", " for(int k=0;k\n", "\n", "# Step 4: The `driver_evaluate_MHD_rhs.h` header file \\[Back to [top](#toc)\\]\n", "$$\\label{driver_evaluate_MHD_rhs__h}$$\n", "\n", "Now we generate the header file for the `driver_evaluate_MHD_rhs.C` file." ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Writing ../src/driver_evaluate_MHD_rhs.h\n" ] } ], "source": [ "%%writefile $outfile_path__driver_evaluate_MHD_rhs__h\n", "#ifndef DRIVER_EVALUATE_MHD_RHS_H_\n", "#define DRIVER_EVALUATE_MHD_RHS_H_\n", "\n", "/* PRIVATE FUNCTIONS, Called within driver_evaluate_MHD_rhs.C ONLY */\n", "static void ftilde_gf_compute(const cGH *cctkGH,const int *cctk_lsh,const int flux_dirn,gf_and_gz_struct *input,CCTK_REAL *ftilde_gf);\n", "static void reconstruct_set_of_prims_PPM(const cGH *cctkGH,const int *cctk_lsh,const int flux_dirn,const int num_prims_to_reconstruct,const int *which_prims_to_reconstruct,\n", " eos_struct &eosi,gf_and_gz_struct *in_prims,gf_and_gz_struct *out_prims_r,gf_and_gz_struct *out_prims_l,\n", " CCTK_REAL *ftilde_gf,CCTK_REAL *temporary);\n", "\n", "static void compute_tau_rhs_extrinsic_curvature_terms_and_TUPmunu\n", "(const cGH *cctkGH,const int *cctk_lsh,const int *cctk_nghostzones,CCTK_REAL *dX,CCTK_REAL **metric,gf_and_gz_struct *prims,\n", " CCTK_REAL **TUPmunu,eos_struct &eos, CCTK_REAL Gamma_th,\n", " CCTK_REAL *gupxy,CCTK_REAL *gupxz,CCTK_REAL *gupyz,\n", " CCTK_REAL *kxx,CCTK_REAL *kxy,CCTK_REAL *kxz,CCTK_REAL *kyy,CCTK_REAL *kyz,CCTK_REAL *kzz,\n", " CCTK_REAL *tau_rhs);\n", "static void A_i_rhs_no_gauge_terms(const int A_dirn, const cGH *cctkGH,const int *cctk_lsh,const int *cctk_nghostzones,gf_and_gz_struct *out_prims_r,gf_and_gz_struct *out_prims_l,\n", " CCTK_REAL *phi_interped,CCTK_REAL *cmax_1,CCTK_REAL *cmin_1,CCTK_REAL *cmax_2,CCTK_REAL *cmin_2, CCTK_REAL *A3_rhs);\n", "\n", "static void Lorenz_psi6phi_rhs__add_gauge_terms_to_A_i_rhs(const cGH *cctkGH,const int *cctk_lsh,const int *cctk_nghostzones,CCTK_REAL *dX,CCTK_REAL **interp_vars,CCTK_REAL *psi6phi,\n", " CCTK_REAL *shiftx_iphjphkph,CCTK_REAL *shifty_iphjphkph,CCTK_REAL *shiftz_iphjphkph,\n", " CCTK_REAL *alpha_iphjphkph,CCTK_REAL *alpha_Phi_minus_betaj_A_j_iphjphkph,CCTK_REAL *alpha_sqrtg_Ax_interp,\n", " CCTK_REAL *alpha_sqrtg_Ay_interp,CCTK_REAL *alpha_sqrtg_Az_interp,\n", " CCTK_REAL *psi6phi_rhs,CCTK_REAL *Ax_rhs,CCTK_REAL *Ay_rhs,CCTK_REAL *Az_rhs);\n", "\n", "static void add_fluxes_and_source_terms_to_hydro_rhss(const int flux_dirn,const cGH *cctkGH,const int *cctk_lsh,const int *cctk_nghostzones,CCTK_REAL *dX,\n", " CCTK_REAL **metric,gf_and_gz_struct *in_prims,CCTK_REAL **TUPmunu,\n", " int numvars_reconstructed,gf_and_gz_struct *out_prims_r,gf_and_gz_struct *out_prims_l,eos_struct &eos,\n", " CCTK_REAL *cmax,CCTK_REAL *cmin,\n", " CCTK_REAL *rho_star_flux,CCTK_REAL *tau_flux,CCTK_REAL *st_x_flux,CCTK_REAL *st_y_flux,CCTK_REAL *st_z_flux,\n", " CCTK_REAL *rho_star_rhs,CCTK_REAL *tau_rhs,CCTK_REAL *st_x_rhs,CCTK_REAL *st_y_rhs,CCTK_REAL *st_z_rhs);\n", "\n", "#include \"harm_primitives_headers.h\"\n", "\n", "#endif /* DRIVER_EVALUATE_MHD_RHS_H_ */\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 5: Code validation \\[Back to [top](#toc)\\]\n", "$$\\label{code_validation}$$\n", "\n", "\n", "\n", "## Step 5.a: `driver_evaluate_MHD_rhs.C` \\[Back to [top](#toc)\\]\n", "$$\\label{code_validation_driver_evaluate_MHD_rhs__c}$$\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": 26, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Validation test for driver_evaluate_MHD_rhs.C: FAILED!\n", "Diff:\n", "3,4c3,4\n", "< * (vector potential prescription), using the \n", "< * generalized Lorenz gauge condition for the \n", "---\n", "> * (vector potential prescription), using the\n", "> * generalized Lorenz gauge condition for the\n", "7c7\n", "< * Based originally on the Illinois GRMHD code, \n", "---\n", "> * Based originally on the Illinois GRMHD code,\n", "10,11c10,11\n", "< * primarily by Zachariah Etienne, Yuk Tung Liu, \n", "< * and Vasileios Paschalidis. \n", "---\n", "> * primarily by Zachariah Etienne, Yuk Tung Liu,\n", "> * and Vasileios Paschalidis.\n", "13c13\n", "< * Rewritten for public release in 2013 \n", "---\n", "> * Rewritten for public release in 2013\n", "17c17\n", "< * Original unigrid GRMHD evolution prescription: \n", "---\n", "> * Original unigrid GRMHD evolution prescription:\n", "32c32\n", "< * This version of PPM implements the standard \n", "---\n", "> * This version of PPM implements the standard\n", "34c34\n", "< * to have 3 ghostzones instead of 4. \n", "---\n", "> * to have 3 ghostzones instead of 4.\n", "36a37\n", "> \n", "45a47\n", "> #include \"IllinoisGRMHD_EoS_lowlevel_functs.C\"\n", "47a50\n", "> \n", "56c59\n", "< \n", "---\n", "> \n", "63c66,77\n", "< // FIXME: only for single gamma-law EOS.\n", "---\n", "> \n", "> /**********************************\n", "> * Piecewise Polytropic EOS Patch *\n", "> * Setting up the EOS struct *\n", "> **********************************/\n", "> /*\n", "> * The short piece of code below takes care\n", "> * of initializing the EOS parameters.\n", "> * Please refer to the \"inlined_functions.C\"\n", "> * source file for the documentation on the\n", "> * function.\n", "> */\n", "65,72c79,80\n", "< eos.neos=neos;\n", "< eos.K_poly=K_poly;\n", "< eos.rho_tab[0]=rho_tab[0];\n", "< eos.P_tab[0]=P_tab[0];\n", "< eos.gamma_th=gamma_th;\n", "< eos.eps_tab[0]=eps_tab[0];\n", "< eos.k_tab[0]=k_tab[0]; eos.k_tab[1]=k_tab[1];\n", "< eos.gamma_tab[0]=gamma_tab[0]; eos.gamma_tab[1]=gamma_tab[1];\n", "---\n", "> initialize_EOS_struct_from_input(eos);\n", "> \n", "106a115\n", "> \n", "112a122\n", "> \n", "145a156\n", "> \n", "162a174\n", "> \n", "167c179,180\n", "< compute_tau_rhs_extrinsic_curvature_terms_and_TUPmunu(cctkGH,cctk_lsh,cctk_nghostzones,dX,metric,in_prims,TUPmunu,eos,\n", "---\n", "> compute_tau_rhs_extrinsic_curvature_terms_and_TUPmunu(cctkGH,cctk_lsh,cctk_nghostzones,dX,metric,in_prims,TUPmunu,\n", "> eos, Gamma_th,\n", "170a184\n", "> \n", "177a192\n", "> \n", "179,180c194,195\n", "< * 1) Computation of \\partial_x on RHS of \\partial_t {rho_star,tau,mhd_st_{x,y,z}}, \n", "< * via PPM reconstruction onto (i-1/2,j,k), so that \n", "---\n", "> * 1) Computation of \\partial_x on RHS of \\partial_t {rho_star,tau,mhd_st_{x,y,z}},\n", "> * via PPM reconstruction onto (i-1/2,j,k), so that\n", "205a221\n", "> \n", "220c236,237\n", "< // Note that we have already reconstructed vx and vy along the x-direction, \n", "---\n", "> \n", "> // Note that we have already reconstructed vx and vy along the x-direction,\n", "237,239c254,256\n", "< /* There are two stories going on here: \n", "< * 1) Computation of \\partial_y on RHS of \\partial_t {rho_star,tau,mhd_st_{x,y,z}}, \n", "< * via PPM reconstruction onto (i,j-1/2,k), so that \n", "---\n", "> /* There are two stories going on here:\n", "> * 1) Computation of \\partial_y on RHS of \\partial_t {rho_star,tau,mhd_st_{x,y,z}},\n", "> * via PPM reconstruction onto (i,j-1/2,k), so that\n", "264c281\n", "< reconstruct_set_of_prims_PPM(cctkGH,cctk_lsh,flux_dirn,num_prims_to_reconstruct,which_prims_to_reconstruct, \n", "---\n", "> reconstruct_set_of_prims_PPM(cctkGH,cctk_lsh,flux_dirn,num_prims_to_reconstruct,which_prims_to_reconstruct,\n", "280c297\n", "< reconstruct_set_of_prims_PPM(cctkGH,cctk_lsh,flux_dirn,num_prims_to_reconstruct,which_prims_to_reconstruct, \n", "---\n", "> reconstruct_set_of_prims_PPM(cctkGH,cctk_lsh,flux_dirn,num_prims_to_reconstruct,which_prims_to_reconstruct,\n", "281a299\n", "> \n", "295a314\n", "> \n", "298c317\n", "< * We want to compute \n", "---\n", "> * We want to compute\n", "314c333\n", "< temporary[CCTK_GFINDEX3D(cctkGH,i,j,k)]= \n", "---\n", "> temporary[CCTK_GFINDEX3D(cctkGH,i,j,k)]=\n", "320a340\n", "> \n", "324c344,345\n", "< // in_prims[{VYR,VYL,VZR,VZL}].gz_{lo,hi} ghostzones are not correct, so we fix \n", "---\n", "> \n", "> // in_prims[{VYR,VYL,VZR,VZL}].gz_{lo,hi} ghostzones are not correct, so we fix\n", "337c358\n", "< /* There are two stories going on here: \n", "---\n", "> /* There are two stories going on here:\n", "361c382\n", "< reconstruct_set_of_prims_PPM(cctkGH,cctk_lsh,flux_dirn,num_prims_to_reconstruct,which_prims_to_reconstruct, \n", "---\n", "> reconstruct_set_of_prims_PPM(cctkGH,cctk_lsh,flux_dirn,num_prims_to_reconstruct,which_prims_to_reconstruct,\n", "377c398\n", "< reconstruct_set_of_prims_PPM(cctkGH,cctk_lsh,flux_dirn,num_prims_to_reconstruct,which_prims_to_reconstruct, \n", "---\n", "> reconstruct_set_of_prims_PPM(cctkGH,cctk_lsh,flux_dirn,num_prims_to_reconstruct,which_prims_to_reconstruct,\n", "384a406\n", "> \n", "396,397c418,419\n", "< in_prims[VXR]=out_prims_r[VX]; \n", "< in_prims[VZR]=out_prims_r[VZ]; \n", "---\n", "> in_prims[VXR]=out_prims_r[VX];\n", "> in_prims[VZR]=out_prims_r[VZ];\n", "404a427\n", "> \n", "407c430\n", "< * We want to compute \n", "---\n", "> * We want to compute\n", "421c444\n", "< temporary[CCTK_GFINDEX3D(cctkGH,i,j,k)]= \n", "---\n", "> temporary[CCTK_GFINDEX3D(cctkGH,i,j,k)]=\n", "427a451\n", "> \n", "430a455\n", "> \n", "447c472\n", "< reconstruct_set_of_prims_PPM(cctkGH,cctk_lsh,flux_dirn,num_prims_to_reconstruct,which_prims_to_reconstruct, \n", "---\n", "> reconstruct_set_of_prims_PPM(cctkGH,cctk_lsh,flux_dirn,num_prims_to_reconstruct,which_prims_to_reconstruct,\n", "452c477\n", "< * We want to compute \n", "---\n", "> * We want to compute\n", "466c491\n", "< temporary[CCTK_GFINDEX3D(cctkGH,i,j,k)]= \n", "---\n", "> temporary[CCTK_GFINDEX3D(cctkGH,i,j,k)]=\n", "472c497,498\n", "< \n", "---\n", "> \n", "> \n", "475a502\n", "> \n", "498c525\n", "< Lorenz_psi6phi_rhs__add_gauge_terms_to_A_i_rhs(cctkGH,cctk_lsh,cctk_nghostzones,dX,interp_vars,psi6phi, \n", "---\n", "> Lorenz_psi6phi_rhs__add_gauge_terms_to_A_i_rhs(cctkGH,cctk_lsh,cctk_nghostzones,dX,interp_vars,psi6phi,\n", "505c532\n", "< // FUN DEBUGGING TOOL (trust me!): \n", "---\n", "> // FUN DEBUGGING TOOL (trust me!):\n", "532a560\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/driver_evaluate_MHD_rhs.C\"\n", "original_IGM_file_name = \"driver_evaluate_MHD_rhs-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__driver_evaluate_MHD_rhs__C = !diff $original_IGM_file_path $outfile_path__driver_evaluate_MHD_rhs__C\n", "\n", "if Validation__driver_evaluate_MHD_rhs__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 driver_evaluate_MHD_rhs.C: PASSED!\")\n", "else:\n", " # If the validation fails, we keep the original IGM source code file\n", " print(\"Validation test for driver_evaluate_MHD_rhs.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__driver_evaluate_MHD_rhs__C:\n", " print(diff_line)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 5.b: `driver_evaluate_MHD_rhs.h` \\[Back to [top](#toc)\\]\n", "$$\\label{code_validation_driver_evaluate_MHD_rhs__h}$$\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": 27, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Validation test for driver_evaluate_MHD_rhs.h: FAILED!\n", "Diff:\n", "6c6\n", "< static void reconstruct_set_of_prims_PPM(const cGH *cctkGH,const int *cctk_lsh,const int flux_dirn,const int num_prims_to_reconstruct,const int *which_prims_to_reconstruct, \n", "---\n", "> static void reconstruct_set_of_prims_PPM(const cGH *cctkGH,const int *cctk_lsh,const int flux_dirn,const int num_prims_to_reconstruct,const int *which_prims_to_reconstruct,\n", "12c12\n", "< CCTK_REAL **TUPmunu,eos_struct &eos,\n", "---\n", "> CCTK_REAL **TUPmunu,eos_struct &eos, CCTK_REAL Gamma_th,\n", "34a35\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/driver_evaluate_MHD_rhs.h\"\n", "original_IGM_file_name = \"driver_evaluate_MHD_rhs-original.h\"\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__driver_evaluate_MHD_rhs__h = !diff $original_IGM_file_path $outfile_path__driver_evaluate_MHD_rhs__h\n", "\n", "if Validation__driver_evaluate_MHD_rhs__h == []:\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 driver_evaluate_MHD_rhs.h: PASSED!\")\n", "else:\n", " # If the validation fails, we keep the original IGM source code file\n", " print(\"Validation test for driver_evaluate_MHD_rhs.h: 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__driver_evaluate_MHD_rhs__h:\n", " print(diff_line)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 6: 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__driver_evaluate_MHD_rhs.pdf](Tutorial-IllinoisGRMHD__driver_evaluate_MHD_rhs.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": 28, "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__driver_evaluate_MHD_rhs.ipynb\n", "#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__driver_evaluate_MHD_rhs.tex\n", "#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__driver_evaluate_MHD_rhs.tex\n", "#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__driver_evaluate_MHD_rhs.tex\n", "!rm -f Tut*.out Tut*.aux Tut*.log" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.3" } }, "nbformat": 4, "nbformat_minor": 2 }