{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# `GiRaFFE_NRPy`: Solving the Induction Equation\n", "\n", "## Author: Patrick Nelson\n", "\n", "This notebook documents the function from the original `GiRaFFE` that calculates the flux for $A_i$ according to the method of Harten, Lax, von Leer, and Einfeldt (HLLE), assuming that we have calculated the values of the velocity and magnetic field on the cell faces according to the piecewise-parabolic method (PPM) of [Colella and Woodward (1984)](https://crd.lbl.gov/assets/pubs_presos/AMCS/ANAG/A141984.pdf), modified for the case of GRFFE. \n", "\n", "**Notebook Status:** Validated \n", "\n", "**Validation Notes:** This code has been validated by showing that it converges to the exact answer at the expected order\n", "\n", "### NRPy+ Source Code for this module: \n", "* [GiRaFFE_NRPy/Afield_flux.py](../../edit/in_progress/GiRaFFE_NRPy/Afield_flux.py)\n", "\n", "Our goal in this module is to write the code necessary to solve the induction equation \n", "$$\n", "\\partial_t A_i = \\underbrace{\\epsilon_{ijk} v^j B^k}_{\\rm Flux\\ terms} - \\underbrace{\\partial_i \\left(\\alpha \\Phi - \\beta^j A_j \\right)}_{\\rm Gauge\\ terms}.\n", "$$\n", "To properly handle the flux terms and avoiding problems with shocks, we cannot simply take a cross product of the velocity and magnetic field at the cell centers. Instead, we must solve the Riemann problem at the cell faces using the reconstructed values of the velocity and magnetic field on either side of the cell faces. The reconstruction is done using PPM (see [here](Tutorial-GiRaFFE_NRPy-PPM.ipynb)); in this module, we will assume that that step has already been done. Metric quantities are assumed to have been interpolated to cell faces, as is done in [this](Tutorial-GiRaFFE_NRPy-Metric_Face_Values.ipynb) tutorial. \n", "\n", "Tóth's [paper](https://www.sciencedirect.com/science/article/pii/S0021999100965197?via%3Dihub), Eqs. 30 and 31, are one of the first implementations of such a scheme. The original GiRaFFE used a 2D version of the algorithm from [Del Zanna, et al. (2002)](https://arxiv.org/abs/astro-ph/0210618); but since we are not using staggered grids, we can greatly simplify this algorithm with respect to the version used in the original `GiRaFFE`. Instead, we will adapt the implementations of the algorithm used in [Mewes, et al. (2020)](https://arxiv.org/abs/2002.06225) and [Giacomazzo, et al. (2011)](https://arxiv.org/abs/1009.2468), Eqs. 3-11. \n", "\n", "We first write the flux contribution to the induction equation RHS as \n", "$$\n", "\\partial_t A_i = -E_i,\n", "$$\n", "where the electric field $E_i$ is given in ideal MHD (of which FFE is a subset) as\n", "$$\n", "-E_i = \\epsilon_{ijk} v^j B^k,\n", "$$\n", "where $v^i$ is the drift velocity, $B^i$ is the magnetic field, and $\\epsilon_{ijk} = \\sqrt{\\gamma} [ijk]$ is the Levi-Civita tensor.\n", "In Cartesian coordinates, \n", "\\begin{align}\n", "-E_x &= [F^y(B^z)]_x = -[F^z(B^y)]_x \\\\\n", "-E_y &= [F^z(B^x)]_y = -[F^x(B^z)]_y \\\\\n", "-E_z &= [F^x(B^y)]_z = -[F^y(B^x)]_z, \\\\\n", "\\end{align}\n", "where \n", "$$\n", "[F^i(B^j)]_k = \\sqrt{\\gamma} (v^i B^j - v^j B^i).\n", "$$\n", "To compute the actual contribution to the RHS in some direction $i$, we average the above listed field as calculated on the $+j$, $-j$, $+k$, and $-k$ faces. That is, at some point $(i,j,k)$ on the grid,\n", "\\begin{align}\n", "-E_x(x_i,y_j,z_k) &= \\frac{1}{4} \\left( [F_{\\rm HLL}^y(B^z)]_{x(i,j+1/2,k)}+[F_{\\rm HLL}^y(B^z)]_{x(i,j-1/2,k)}-[F_{\\rm HLL}^z(B^y)]_{x(i,j,k+1/2)}-[F_{\\rm HLL}^z(B^y)]_{x(i,j,k-1/2)} \\right) \\\\\n", "-E_y(x_i,y_j,z_k) &= \\frac{1}{4} \\left( [F_{\\rm HLL}^z(B^x)]_{y(i,j,k+1/2)}+[F_{\\rm HLL}^z(B^x)]_{y(i,j,k-1/2)}-[F_{\\rm HLL}^x(B^z)]_{y(i+1/2,j,k)}-[F_{\\rm HLL}^x(B^z)]_{y(i-1/2,j,k)} \\right) \\\\\n", "-E_z(x_i,y_j,z_k) &= \\frac{1}{4} \\left( [F_{\\rm HLL}^x(B^y)]_{z(i+1/2,j,k)}+[F_{\\rm HLL}^x(B^y)]_{z(i-1/2,j,k)}-[F_{\\rm HLL}^y(B^x)]_{z(i,j+1/2,k)}-[F_{\\rm HLL}^y(B^x)]_{z(i,j-1/2,k)} \\right). \\\\\n", "\\end{align}\n", "Note the use of $F_{\\rm HLL}$ here. This change signifies that the quantity output here is from the HLLE Riemann solver. Note also the indices on the fluxes. Values of $\\pm 1/2$ indicate that these are computed on cell faces using the reconstructed values of $v^i$ and $B^i$ and the interpolated values of the metric gridfunctions. So, \n", "$$\n", "F_{\\rm HLL}^i(B^j) = \\frac{c_{\\rm min} F_{\\rm R}^i(B^j) + c_{\\rm max} F_{\\rm L}^i(B^j) - c_{\\rm min} c_{\\rm max} (B_{\\rm R}^j-B_{\\rm L}^j)}{c_{\\rm min} + c_{\\rm max}}.\n", "$$\n", "\n", "The speeds $c_\\min$ and $c_\\max$ are characteristic speeds that waves can travel through the plasma. In GRFFE, the expressions defining them reduce a function of only the metric quantities. $c_\\min$ is the negative of the minimum amongst the speeds $c_-$ and $0$ and $c_\\max$ is the maximum amongst the speeds $c_+$ and $0$. The speeds $c_\\pm = \\left. \\left(-b \\pm \\sqrt{b^2-4ac}\\right)\\middle/ \\left(2a\\right) \\right.$ must be calculated on both the left and right faces, where \n", "$$a = 1/\\alpha^2,$$ \n", "$$b = 2 \\beta^i / \\alpha^2$$\n", "and $$c = g^{ii} - (\\beta^i)^2/\\alpha^2.$$\n", "An outline of a general finite-volume method is as follows, with the current step in bold:\n", "1. The Reconstruction Step - Piecewise Parabolic Method\n", " 1. Within each cell, fit to a function that conserves the volume in that cell using information from the neighboring cells\n", " * For PPM, we will naturally use parabolas\n", " 1. Use that fit to define the state at the left and right interface of each cell\n", " 1. Apply a slope limiter to mitigate Gibbs phenomenon\n", "1. Interpolate the value of the metric gridfunctions on the cell faces\n", "1. **Solving the Riemann Problem - Harten, Lax, (This notebook, $E_i$ only)**\n", " 1. **Use the left and right reconstructed states to calculate the unique state at boundary**\n", "\n", "We will assume in this notebook that the reconstructed velocities and magnetic fields are available on cell faces as input. We will also assume that the metric gridfunctions have been interpolated on the metric faces. \n", "\n", "Solving the Riemann problem, then, consists of two substeps: First, we compute the flux through each face of the cell. Then, we add the average of these fluxes to the right-hand side of the evolution equation for the vector potential. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Table of Contents\n", "$$\\label{toc}$$\n", "\n", "This notebook is organized as follows\n", "\n", "1. [Step 1](#prelim): Preliminaries\n", "1. [Step 2](#a_i_flux): Computing the Magnetic Flux\n", " 1. [Step 2.a](#hydro_speed): GRFFE characteristic wave speeds\n", " 1. [Step 2.b](#fluxes): Compute the HLLE fluxes\n", "1. [Step 3](#code_validation): Code Validation against `GiRaFFE_NRPy.Afield_flux` NRPy+ Module\n", "1. [Step 4](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 1: Preliminaries \\[Back to [top](#toc)\\]\n", "$$\\label{prelim}$$\n", "\n", "We begin by importing the NRPy+ core functionality. We also import the Levi-Civita symbol, the GRHD module, and the GRFFE module." ] }, { "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", "from outputC import outCfunction, outputC # NRPy+: Core C code output module\n", "import sympy as sp # SymPy: The Python computer algebra package upon which NRPy+ depends\n", "import NRPy_param_funcs as par # NRPy+: Parameter interface\n", "import grid as gri # NRPy+: Functions having to do with numerical grids\n", "import indexedexp as ixp # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support\n", "import cmdline_helper as cmd # NRPy+: Multi-platform Python command-line interface\n", "\n", "thismodule = \"GiRaFFE_NRPy-Afield_flux\"\n", "\n", "import GRHD.equations as GRHD\n", "# import GRFFE.equations as GRFFE" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 2: Computing the Magnetic Flux \\[Back to [top](#toc)\\]\n", "$$\\label{a_i_flux}$$\n", "\n", "\n", "\n", "## Step 2.a: GRFFE characteristic wave speeds \\[Back to [top](#toc)\\]\n", "$$\\label{hydro_speed}$$\n", "\n", "Next, we will find the speeds at which the hydrodynamics waves propagate. We start from the speed of light (since FFE deals with very diffuse plasmas), which is $c=1.0$ in our chosen units. We then find the speeds $c_+$ and $c_-$ on each face with the function `find_cp_cm`; then, we find minimum and maximum speeds possible from among those.\n", "\n", "\n", "\n", "Below is the source code for `find_cp_cm`, edited to work with the NRPy+ version of GiRaFFE. One edit we need to make in particular is to the term `psim4*gupii` in the definition of `c`; that was written assuming the use of the conformal metric $\\tilde{g}^{ii}$. Since we are not using that here, and are instead using the ADM metric, we should not multiply by $\\psi^{-4}$.\n", "\n", "```c\n", "static inline void find_cp_cm(REAL &cplus,REAL &cminus,const REAL v02,const REAL u0,\n", " const REAL vi,const REAL lapse,const REAL shifti,\n", " const REAL gammadet,const REAL gupii) {\n", " const REAL u0_SQUARED=u0*u0;\n", " const REAL ONE_OVER_LAPSE_SQUARED = 1.0/(lapse*lapse);\n", " // sqrtgamma = psi6 -> psim4 = gammadet^(-1.0/3.0)\n", " const REAL psim4 = pow(gammadet,-1.0/3.0);\n", " //Find cplus, cminus:\n", " const REAL a = u0_SQUARED * (1.0-v02) + v02*ONE_OVER_LAPSE_SQUARED;\n", " const REAL b = 2.0* ( shifti*ONE_OVER_LAPSE_SQUARED * v02 - u0_SQUARED * vi * (1.0-v02) );\n", " const REAL c = u0_SQUARED*vi*vi * (1.0-v02) - v02 * ( gupii -\n", " shifti*shifti*ONE_OVER_LAPSE_SQUARED);\n", " REAL detm = b*b - 4.0*a*c;\n", " //ORIGINAL LINE OF CODE:\n", " //if(detm < 0.0) detm = 0.0;\n", " //New line of code (without the if() statement) has the same effect:\n", " detm = sqrt(0.5*(detm + fabs(detm))); /* Based on very nice suggestion from Roland Haas */\n", " \n", " cplus = 0.5*(detm-b)/a;\n", " cminus = -0.5*(detm+b)/a;\n", " if (cplus < cminus) {\n", " const REAL cp = cminus;\n", " cminus = cplus;\n", " cplus = cp;\n", " }\n", "}\n", "```\n", "Comments documenting this have been excised for brevity, but are reproduced in $\\LaTeX$ [below](#derive_speed).\n", "\n", "We could use this code directly, but there's substantial improvement we can make by changing the code into a NRPyfied form. Note the `if` statement; NRPy+ does not know how to handle these, so we must eliminate it if we want to leverage NRPy+'s full power. (Calls to `fabs()` are also cheaper than `if` statements.) This can be done if we rewrite this, taking inspiration from the other eliminated `if` statement documented in the above code block:\n", "```c\n", " cp = 0.5*(detm-b)/a;\n", " cm = -0.5*(detm+b)/a;\n", " cplus = 0.5*(cp+cm+fabs(cp-cm));\n", " cminus = 0.5*(cp+cm-fabs(cp-cm));\n", "```\n", "This can be simplified further, by substituting `cp` and `cm` into the below equations and eliminating terms as appropriate. First note that `cp+cm = -b/a` and that `cp-cm = detm/a`. Thus,\n", "```c\n", " cplus = 0.5*(-b/a + fabs(detm/a));\n", " cminus = 0.5*(-b/a - fabs(detm/a));\n", "```\n", "This fulfills the original purpose of the `if` statement in the original code because we have guaranteed that $c_+ \\geq c_-$.\n", "\n", "This leaves us with an expression that can be much more easily NRPyfied. So, we will rewrite the following in NRPy+, making only minimal changes to be proper Python. However, it turns out that we can make this even simpler. In GRFFE, $v_0^2$ is guaranteed to be exactly one. In GRMHD, this speed was calculated as $$v_{0}^{2} = v_{\\rm A}^{2} + c_{\\rm s}^{2}\\left(1-v_{\\rm A}^{2}\\right),$$ where the Alfvén speed $v_{\\rm A}^{2}$ $$v_{\\rm A}^{2} = \\frac{b^{2}}{\\rho_{b}h + b^{2}}.$$ So, we can see that when the density $\\rho_b$ goes to zero, $v_{0}^{2} = v_{\\rm A}^{2} = 1$. Then \n", "\\begin{align}\n", "a &= (u^0)^2 (1-v_0^2) + v_0^2/\\alpha^2 \\\\\n", "&= 1/\\alpha^2 \\\\\n", "b &= 2 \\left(\\beta^i v_0^2 / \\alpha^2 - (u^0)^2 v^i (1-v_0^2)\\right) \\\\\n", "&= 2 \\beta^i / \\alpha^2 \\\\\n", "c &= (u^0)^2 (v^i)^2 (1-v_0^2) - v_0^2 \\left(\\gamma^{ii} - (\\beta^i)^2/\\alpha^2\\right) \\\\\n", "&= -\\gamma^{ii} + (\\beta^i)^2/\\alpha^2,\n", "\\end{align}\n", "are simplifications that should save us some time; we can see that $a \\geq 0$ is guaranteed. Note that we also force `detm` to be positive. Thus, `detm/a` is guaranteed to be positive itself, rendering the calls to `nrpyAbs()` superfluous. Furthermore, we eliminate any dependence on the Valencia 3-velocity and the time compoenent of the four-velocity, $u^0$. This leaves us free to solve the quadratic in the familiar way: $$c_\\pm = \\frac{-b \\pm \\sqrt{b^2-4ac}}{2a}$$." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# We'll write this as a function so that we can calculate the expressions on-demand for any choice of i\n", "def find_cp_cm(lapse,shifti,gammaUUii):\n", " # Inputs: u0,vi,lapse,shift,gammadet,gupii\n", " # Outputs: cplus,cminus\n", "\n", " # a = 1/(alpha^2)\n", " a = sp.sympify(1)/(lapse*lapse)\n", " # b = 2 beta^i / alpha^2\n", " b = sp.sympify(2) * shifti /(lapse*lapse)\n", " # c = -g^{ii} + (beta^i)^2 / alpha^2\n", " c = - gammaUUii + shifti*shifti/(lapse*lapse)\n", "\n", " # Now, we are free to solve the quadratic equation as usual. We take care to avoid passing a\n", " # negative value to the sqrt function.\n", " detm = b*b - sp.sympify(4)*a*c\n", "\n", " import Min_Max_and_Piecewise_Expressions as noif\n", " detm = sp.sqrt(noif.max_noif(sp.sympify(0),detm))\n", " global cplus,cminus\n", " cplus = sp.Rational(1,2)*(-b/a + detm/a)\n", " cminus = sp.Rational(1,2)*(-b/a - detm/a)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In flat spacetime, where $\\alpha=1$, $\\beta^i=0$, and $\\gamma^{ij} = \\delta^{ij}$, $c_+ > 0$ and $c_- < 0$. For the HLLE solver, we will need both `cmax` and `cmin` to be positive; we also want to choose the speed that is larger in magnitude because overestimating the characteristic speeds will help damp unwanted oscillations. (However, in GRFFE, we only get one $c_+$ and one $c_-$, so we only need to fix the signs here.) Hence, the following function. \n", "\n", "We will now write a function in NRPy+ similar to the one used in the old `GiRaFFE`, allowing us to generate the expressions with less need to copy-and-paste code; the key difference is that this one will be in Python, and generate optimized C code integrated into the rest of the operations. Notice that since we eliminated the dependence on velocities, none of the input quantities are different on either side of the face. So, this function won't really do much besides guarantee that `cmax` and `cmin` are positive, but we'll leave the machinery here since it is likely to be a useful guide to somebody who wants to something similar. The only modifications we'll make are those necessary to eliminate calls to `fabs(0)` in the C code. We use the same technique as above to replace the `if` statements inherent to the `MAX()` and `MIN()` functions. " ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# We'll write this as a function, and call it within HLLE_solver, below.\n", "def find_cmax_cmin(field_comp,gamma_faceDD,beta_faceU,alpha_face):\n", " # Inputs: flux direction field_comp, Inverse metric gamma_faceUU, shift beta_faceU,\n", " # lapse alpha_face, metric determinant gammadet_face\n", " # Outputs: maximum and minimum characteristic speeds cmax and cmin\n", " # First, we need to find the characteristic speeds on each face\n", " gamma_faceUU,unusedgammaDET = ixp.generic_matrix_inverter3x3(gamma_faceDD)\n", " # Original needed for GRMHD\n", "# find_cp_cm(alpha_face,beta_faceU[field_comp],gamma_faceUU[field_comp][field_comp])\n", "# cpr = cplus\n", "# cmr = cminus\n", "# find_cp_cm(alpha_face,beta_faceU[field_comp],gamma_faceUU[field_comp][field_comp])\n", "# cpl = cplus\n", "# cml = cminus\n", " find_cp_cm(alpha_face,beta_faceU[field_comp],gamma_faceUU[field_comp][field_comp])\n", " cp = cplus\n", " cm = cminus\n", "\n", " # The following algorithms have been verified with random floats:\n", "\n", " global cmax,cmin\n", " # Now, we need to set cmax to the larger of cpr,cpl, and 0\n", "\n", " import Min_Max_and_Piecewise_Expressions as noif\n", " cmax = noif.max_noif(cp,sp.sympify(0))\n", "\n", " # And then, set cmin to the smaller of cmr,cml, and 0\n", " cmin = -noif.min_noif(cm,sp.sympify(0))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 2.b: Compute the HLLE fluxes \\[Back to [top](#toc)\\]\n", "$$\\label{fluxes}$$\n", "\n", "Here, we we calculate the flux and state vectors for the electric field. The flux vector is here given as \n", "$$\n", "[F^i(B^j)]_k = \\sqrt{\\gamma} (v^i B^j - v^j B^i).\n", "$$\n", "Here, $v^i$ is the drift velocity and $B^i$ is the magnetic field.\n", "This can be easily handled for an input flux direction $i$ with\n", "$$\n", "[F^j(B^k)]_i = \\epsilon_{ijk} v^j B^k,\n", "$$\n", "where $\\epsilon_{ijk} = \\sqrt{\\gamma} [ijk]$ and $[ijk]$ is the Levi-Civita symbol.\n", "\n", "The state vector is simply the magnetic field $B^j$." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "def calculate_flux_and_state_for_Induction(field_comp,flux_dirn, gammaDD,betaU,alpha,ValenciavU,BU):\n", " # Define Levi-Civita symbol\n", " def define_LeviCivitaSymbol_rank3(DIM=-1):\n", " if DIM == -1:\n", " DIM = par.parval_from_str(\"DIM\")\n", "\n", " LeviCivitaSymbol = ixp.zerorank3()\n", "\n", " for i in range(DIM):\n", " for j in range(DIM):\n", " for k in range(DIM):\n", " # From https://codegolf.stackexchange.com/questions/160359/levi-civita-symbol :\n", " LeviCivitaSymbol[i][j][k] = (i - j) * (j - k) * (k - i) * sp.Rational(1,2)\n", " return LeviCivitaSymbol\n", " GRHD.compute_sqrtgammaDET(gammaDD)\n", " # Here, we import the Levi-Civita tensor and compute the tensor with lower indices\n", " LeviCivitaDDD = define_LeviCivitaSymbol_rank3()\n", " for i in range(3):\n", " for j in range(3):\n", " for k in range(3):\n", " LeviCivitaDDD[i][j][k] *= GRHD.sqrtgammaDET\n", "\n", " global U,F\n", " # Flux F = \\epsilon_{ijk} v^j B^k\n", " F = sp.sympify(0)\n", " for j in range(3):\n", " for k in range(3):\n", " F += LeviCivitaDDD[field_comp][j][k] * (alpha*ValenciavU[j]-betaU[j]) * BU[k]\n", " # U = B^i\n", " U = BU[flux_dirn]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, we write a standard HLLE solver based on eq. 3.15 in [the HLLE paper](https://epubs.siam.org/doi/pdf/10.1137/1025002),\n", "$$\n", "F^{\\rm HLL} = \\frac{c_{\\rm min} F_{\\rm R} + c_{\\rm max} F_{\\rm L} - c_{\\rm min} c_{\\rm max} (U_{\\rm R}-U_{\\rm L})}{c_{\\rm min} + c_{\\rm max}}\n", "$$" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "def HLLE_solver(cmax, cmin, Fr, Fl, Ur, Ul):\n", " # This solves the Riemann problem for the flux of E_i in one direction\n", "\n", " # F^HLL = (c_\\min f_R + c_\\max f_L - c_\\min c_\\max ( st_j_r - st_j_l )) / (c_\\min + c_\\max)\n", " return (cmin*Fr + cmax*Fl - cmin*cmax*(Ur-Ul) )/(cmax + cmin)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here, we will use the function we just wrote to calculate the flux through a face. We will pass the reconstructed Valencia 3-velocity and magnetic field on either side of an interface to this function (designated as the \"left\" and \"right\" sides) along with the value of the 3-metric, shift vector, and lapse function on the interface. The parameter `flux_dirn` specifies which face through which we are calculating the flux. However, unlike when we used this method to calculate the flux term, the RHS of each component of $A_i$ does not depend on all three of the flux directions. Instead, the flux of one component of the $E_i$ field depends on flux through the faces in the other two directions. This will be handled when we generate the C function, as demonstrated in the example code after this next function.\n", "\n", "Note that we allow the user to declare their own gridfunctions if they wish, and default to declaring basic symbols if they are not provided. The default names are chosen to imply interpolation of the metric gridfunctions and reconstruction of the primitives." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "def calculate_E_i_flux(flux_dirn,alpha_face=None,gamma_faceDD=None,beta_faceU=None,\\\n", " Valenciav_rU=None,B_rU=None,Valenciav_lU=None,B_lU=None):\n", " global E_fluxD\n", " E_fluxD = ixp.zerorank1()\n", " for field_comp in range(3):\n", " find_cmax_cmin(field_comp,gamma_faceDD,beta_faceU,alpha_face)\n", " calculate_flux_and_state_for_Induction(field_comp,flux_dirn, gamma_faceDD,beta_faceU,alpha_face,\\\n", " Valenciav_rU,B_rU)\n", " Fr = F\n", " Ur = U\n", " calculate_flux_and_state_for_Induction(field_comp,flux_dirn, gamma_faceDD,beta_faceU,alpha_face,\\\n", " Valenciav_lU,B_lU)\n", " Fl = F\n", " Ul = U\n", " E_fluxD[field_comp] += HLLE_solver(cmax, cmin, Fr, Fl, Ur, Ul)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Below, we will write some example code to use the above functions to generate C code for `GiRaFFE_NRPy`. We need to write our own memory reads and writes because we need to add contributions from *both* faces in a given direction, which is expressed in the code as adding contributions from adjacent gridpoints to the RHS, which is not something `FD_outputC` can handle. The `.replace()` function calls adapt these reads and writes to the different directions. Note that, for reconstructions in a given direction, the fluxes are only added to the other two components, as can be seen in the equations we are implementing.\n", "\\begin{align}\n", "-E_x(x_i,y_j,z_k) &= \\frac{1}{4} \\left( [F_{\\rm HLL}^y(B^z)]_{x(i,j+1/2,k)}+[F_{\\rm HLL}^y(B^z)]_{x(i,j-1/2,k)}-[F_{\\rm HLL}^z(B^y)]_{x(i,j,k+1/2)}-[F_{\\rm HLL}^z(B^y)]_{x(i,j,k-1/2)} \\right) \\\\\n", "-E_y(x_i,y_j,z_k) &= \\frac{1}{4} \\left( [F_{\\rm HLL}^z(B^x)]_{y(i,j,k+1/2)}+[F_{\\rm HLL}^z(B^x)]_{y(i,j,k-1/2)}-[F_{\\rm HLL}^x(B^z)]_{y(i+1/2,j,k)}-[F_{\\rm HLL}^x(B^z)]_{y(i-1/2,j,k)} \\right) \\\\\n", "-E_z(x_i,y_j,z_k) &= \\frac{1}{4} \\left( [F_{\\rm HLL}^x(B^y)]_{z(i+1/2,j,k)}+[F_{\\rm HLL}^x(B^y)]_{z(i-1/2,j,k)}-[F_{\\rm HLL}^y(B^x)]_{z(i,j+1/2,k)}-[F_{\\rm HLL}^y(B^x)]_{z(i,j-1/2,k)} \\right). \\\\\n", "\\end{align}\n", "From this, we can see that when, for instance, we reconstruct and interpolate in the $x$-direction, we must add only to the $y$- and $z$-components of the electric field.\n", "\n", "Recall that when we reconstructed the velocity and magnetic field, we constructed to the $i-1/2$ face, so the data at $i+1/2$ is stored at $i+1$.\n" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "def generate_Afield_flux_function_files(out_dir,subdir,alpha_face,gamma_faceDD,beta_faceU,\\\n", " Valenciav_rU,B_rU,Valenciav_lU,B_lU,inputs_provided=True):\n", " if not inputs_provided:\n", " # declare all variables\n", " alpha_face = sp.symbols(alpha_face)\n", " beta_faceU = ixp.declarerank1(\"beta_faceU\")\n", " gamma_faceDD = ixp.declarerank2(\"gamma_faceDD\",\"sym01\")\n", " Valenciav_rU = ixp.declarerank1(\"Valenciav_rU\")\n", " B_rU = ixp.declarerank1(\"B_rU\")\n", " Valenciav_lU = ixp.declarerank1(\"Valenciav_lU\")\n", " B_lU = ixp.declarerank1(\"B_lU\")\n", "\n", " Memory_Read = \"\"\"const double alpha_face = auxevol_gfs[IDX4S(ALPHA_FACEGF, i0,i1,i2)];\n", " const double gamma_faceDD00 = auxevol_gfs[IDX4S(GAMMA_FACEDD00GF, i0,i1,i2)];\n", " const double gamma_faceDD01 = auxevol_gfs[IDX4S(GAMMA_FACEDD01GF, i0,i1,i2)];\n", " const double gamma_faceDD02 = auxevol_gfs[IDX4S(GAMMA_FACEDD02GF, i0,i1,i2)];\n", " const double gamma_faceDD11 = auxevol_gfs[IDX4S(GAMMA_FACEDD11GF, i0,i1,i2)];\n", " const double gamma_faceDD12 = auxevol_gfs[IDX4S(GAMMA_FACEDD12GF, i0,i1,i2)];\n", " const double gamma_faceDD22 = auxevol_gfs[IDX4S(GAMMA_FACEDD22GF, i0,i1,i2)];\n", " const double beta_faceU0 = auxevol_gfs[IDX4S(BETA_FACEU0GF, i0,i1,i2)];\n", " const double beta_faceU1 = auxevol_gfs[IDX4S(BETA_FACEU1GF, i0,i1,i2)];\n", " const double beta_faceU2 = auxevol_gfs[IDX4S(BETA_FACEU2GF, i0,i1,i2)];\n", " const double Valenciav_rU0 = auxevol_gfs[IDX4S(VALENCIAV_RU0GF, i0,i1,i2)];\n", " const double Valenciav_rU1 = auxevol_gfs[IDX4S(VALENCIAV_RU1GF, i0,i1,i2)];\n", " const double Valenciav_rU2 = auxevol_gfs[IDX4S(VALENCIAV_RU2GF, i0,i1,i2)];\n", " const double B_rU0 = auxevol_gfs[IDX4S(B_RU0GF, i0,i1,i2)];\n", " const double B_rU1 = auxevol_gfs[IDX4S(B_RU1GF, i0,i1,i2)];\n", " const double B_rU2 = auxevol_gfs[IDX4S(B_RU2GF, i0,i1,i2)];\n", " const double Valenciav_lU0 = auxevol_gfs[IDX4S(VALENCIAV_LU0GF, i0,i1,i2)];\n", " const double Valenciav_lU1 = auxevol_gfs[IDX4S(VALENCIAV_LU1GF, i0,i1,i2)];\n", " const double Valenciav_lU2 = auxevol_gfs[IDX4S(VALENCIAV_LU2GF, i0,i1,i2)];\n", " const double B_lU0 = auxevol_gfs[IDX4S(B_LU0GF, i0,i1,i2)];\n", " const double B_lU1 = auxevol_gfs[IDX4S(B_LU1GF, i0,i1,i2)];\n", " const double B_lU2 = auxevol_gfs[IDX4S(B_LU2GF, i0,i1,i2)];\n", " REAL A_rhsD0 = 0; REAL A_rhsD1 = 0; REAL A_rhsD2 = 0;\n", " \"\"\"\n", " Memory_Write = \"\"\"rhs_gfs[IDX4S(AD0GF,i0,i1,i2)] += A_rhsD0;\n", " rhs_gfs[IDX4S(AD1GF,i0,i1,i2)] += A_rhsD1;\n", " rhs_gfs[IDX4S(AD2GF,i0,i1,i2)] += A_rhsD2;\n", " \"\"\"\n", "\n", " indices = [\"i0\",\"i1\",\"i2\"]\n", " indicesp1 = [\"i0+1\",\"i1+1\",\"i2+1\"]\n", "\n", " for flux_dirn in range(3):\n", " calculate_E_i_flux(flux_dirn,alpha_face,gamma_faceDD,beta_faceU,\\\n", " Valenciav_rU,B_rU,Valenciav_lU,B_lU)\n", "\n", " E_field_to_print = [\\\n", " sp.Rational(1,4)*E_fluxD[(flux_dirn+1)%3],\n", " sp.Rational(1,4)*E_fluxD[(flux_dirn+2)%3],\n", " ]\n", "\n", " E_field_names = [\\\n", " \"A_rhsD\"+str((flux_dirn+1)%3),\n", " \"A_rhsD\"+str((flux_dirn+2)%3),\n", " ]\n", "\n", " desc = \"Calculate the electric flux on the left face in direction \" + str(flux_dirn) + \".\"\n", " name = \"calculate_E_field_D\" + str(flux_dirn) + \"_right\"\n", " outCfunction(\n", " outfile = os.path.join(out_dir,subdir,name+\".h\"), desc=desc, name=name,\n", " params =\"const paramstruct *params,const REAL *auxevol_gfs,REAL *rhs_gfs\",\n", " body = Memory_Read \\\n", " +outputC(E_field_to_print,E_field_names,\"returnstring\",params=\"outCverbose=False\")\\\n", " +Memory_Write,\n", " loopopts =\"InteriorPoints\",\n", " rel_path_to_Cparams=os.path.join(\"../\"))\n", "\n", " desc = \"Calculate the electric flux on the left face in direction \" + str(flux_dirn) + \".\"\n", " name = \"calculate_E_field_D\" + str(flux_dirn) + \"_left\"\n", " outCfunction(\n", " outfile = os.path.join(out_dir,subdir,name+\".h\"), desc=desc, name=name,\n", " params =\"const paramstruct *params,const REAL *auxevol_gfs,REAL *rhs_gfs\",\n", " body = Memory_Read.replace(indices[flux_dirn],indicesp1[flux_dirn]) \\\n", " +outputC(E_field_to_print,E_field_names,\"returnstring\",params=\"outCverbose=False\")\\\n", " +Memory_Write,\n", " loopopts =\"InteriorPoints\",\n", " rel_path_to_Cparams=os.path.join(\"../\"))\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 3: Code Validation against `GiRaFFE_NRPy.Induction_Equation` NRPy+ Module \\[Back to [top](#toc)\\]\n", "$$\\label{code_validation}$$\n", "\n", "\n", "Here, as a code validation check, we verify agreement in the SymPy expressions for the `GiRaFFE` evolution equations and auxiliary quantities we intend to use between\n", "1. this tutorial and \n", "2. the NRPy+ [GiRaFFE_NRPy.Induction_Equation](../../edit/in_progress/GiRaFFE_NRPy/Induction_Equation.py) module.\n", "\n", "Below are the gridfunction registrations we will need for testing. We will pass these to the above functions to self-validate the module that corresponds with this tutorial." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "ALL TESTS PASSED!\n" ] } ], "source": [ "all_passed=True\n", "def comp_func(expr1,expr2,basename,prefixname2=\"C2P_P2C.\"):\n", " if str(expr1-expr2)!=\"0\":\n", " print(basename+\" - \"+prefixname2+basename+\" = \"+ str(expr1-expr2))\n", " all_passed=False\n", "\n", "def gfnm(basename,idx1,idx2=None,idx3=None):\n", " if idx2 is None:\n", " return basename+\"[\"+str(idx1)+\"]\"\n", " if idx3 is None:\n", " return basename+\"[\"+str(idx1)+\"][\"+str(idx2)+\"]\"\n", " return basename+\"[\"+str(idx1)+\"][\"+str(idx2)+\"][\"+str(idx3)+\"]\"\n", "\n", "# These are the standard gridfunctions we've used before.\n", "#ValenciavU = ixp.register_gridfunctions_for_single_rank1(\"AUXEVOL\",\"ValenciavU\",DIM=3)\n", "#gammaDD = ixp.register_gridfunctions_for_single_rank2(\"AUXEVOL\",\"gammaDD\",\"sym01\")\n", "#betaU = ixp.register_gridfunctions_for_single_rank1(\"AUXEVOL\",\"betaU\")\n", "#alpha = gri.register_gridfunctions(\"AUXEVOL\",[\"alpha\"])\n", "#AD = ixp.register_gridfunctions_for_single_rank1(\"EVOL\",\"AD\",DIM=3)\n", "#BU = ixp.register_gridfunctions_for_single_rank1(\"AUXEVOL\",\"BU\",DIM=3)\n", "\n", "# We will pass values of the gridfunction on the cell faces into the function. This requires us\n", "# to declare them as C parameters in NRPy+. We will denote this with the _face infix/suffix.\n", "alpha_face = gri.register_gridfunctions(\"AUXEVOL\",\"alpha_face\")\n", "gamma_faceDD = ixp.register_gridfunctions_for_single_rank2(\"AUXEVOL\",\"gamma_faceDD\",\"sym01\")\n", "beta_faceU = ixp.register_gridfunctions_for_single_rank1(\"AUXEVOL\",\"beta_faceU\")\n", "\n", "# We'll need some more gridfunctions, now, to represent the reconstructions of BU and ValenciavU\n", "# on the right and left faces\n", "Valenciav_rU = ixp.register_gridfunctions_for_single_rank1(\"AUXEVOL\",\"Valenciav_rU\",DIM=3)\n", "B_rU = ixp.register_gridfunctions_for_single_rank1(\"AUXEVOL\",\"B_rU\",DIM=3)\n", "Valenciav_lU = ixp.register_gridfunctions_for_single_rank1(\"AUXEVOL\",\"Valenciav_lU\",DIM=3)\n", "B_lU = ixp.register_gridfunctions_for_single_rank1(\"AUXEVOL\",\"B_lU\",DIM=3)\n", "\n", "import GiRaFFE_NRPy.Afield_flux as Af\n", "\n", "expr_list = []\n", "exprcheck_list = []\n", "namecheck_list = []\n", "\n", "for flux_dirn in range(3):\n", " calculate_E_i_flux(flux_dirn,alpha_face,gamma_faceDD,beta_faceU,\\\n", " Valenciav_rU,B_rU,Valenciav_lU,B_lU)\n", " Af.calculate_E_i_flux(flux_dirn,alpha_face,gamma_faceDD,beta_faceU,\\\n", " Valenciav_rU,B_rU,Valenciav_lU,B_lU)\n", " namecheck_list.extend([gfnm(\"E_fluxD\",flux_dirn)])\n", " exprcheck_list.extend([Af.E_fluxD[flux_dirn]])\n", " expr_list.extend([E_fluxD[flux_dirn]])\n", "\n", "for mom_comp in range(len(expr_list)):\n", " comp_func(expr_list[mom_comp],exprcheck_list[mom_comp],namecheck_list[mom_comp])\n", "\n", "import sys\n", "if all_passed:\n", " print(\"ALL TESTS PASSED!\")\n", "else:\n", " print(\"ERROR: AT LEAST ONE TEST DID NOT PASS\")\n", " sys.exit(1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will also check the output C code to make sure it matches what is produced by the python module." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Output C function calculate_E_field_D0_right() to file GiRaFFE_standalone_Ccodes\\RHSs\\calculate_E_field_D0_right.h\n", "Output C function calculate_E_field_D0_left() to file GiRaFFE_standalone_Ccodes\\RHSs\\calculate_E_field_D0_left.h\n", "Output C function calculate_E_field_D1_right() to file GiRaFFE_standalone_Ccodes\\RHSs\\calculate_E_field_D1_right.h\n", "Output C function calculate_E_field_D1_left() to file GiRaFFE_standalone_Ccodes\\RHSs\\calculate_E_field_D1_left.h\n", "Output C function calculate_E_field_D2_right() to file GiRaFFE_standalone_Ccodes\\RHSs\\calculate_E_field_D2_right.h\n", "Output C function calculate_E_field_D2_left() to file GiRaFFE_standalone_Ccodes\\RHSs\\calculate_E_field_D2_left.h\n", "Output C function calculate_E_field_D0_right() to file GiRaFFE_Ccodes_validation\\RHSs\\calculate_E_field_D0_right.h\n", "Output C function calculate_E_field_D0_left() to file GiRaFFE_Ccodes_validation\\RHSs\\calculate_E_field_D0_left.h\n", "Output C function calculate_E_field_D1_right() to file GiRaFFE_Ccodes_validation\\RHSs\\calculate_E_field_D1_right.h\n", "Output C function calculate_E_field_D1_left() to file GiRaFFE_Ccodes_validation\\RHSs\\calculate_E_field_D1_left.h\n", "Output C function calculate_E_field_D2_right() to file GiRaFFE_Ccodes_validation\\RHSs\\calculate_E_field_D2_right.h\n", "Output C function calculate_E_field_D2_left() to file GiRaFFE_Ccodes_validation\\RHSs\\calculate_E_field_D2_left.h\n", "Printing difference between original C code and this code...\n", "Checking file RHSs/calculate_E_field_D0_right.h\n", "No difference. TEST PASSED!\n", "Checking file RHSs/calculate_E_field_D0_left.h\n", "No difference. TEST PASSED!\n", "Checking file RHSs/calculate_E_field_D1_right.h\n", "No difference. TEST PASSED!\n", "Checking file RHSs/calculate_E_field_D1_left.h\n", "No difference. TEST PASSED!\n", "Checking file RHSs/calculate_E_field_D2_right.h\n", "No difference. TEST PASSED!\n", "Checking file RHSs/calculate_E_field_D2_left.h\n", "No difference. TEST PASSED!\n" ] } ], "source": [ "import difflib\n", "import sys\n", "\n", "subdir = os.path.join(\"RHSs\")\n", "\n", "out_dir = os.path.join(\"GiRaFFE_standalone_Ccodes\")\n", "cmd.mkdir(out_dir)\n", "cmd.mkdir(os.path.join(out_dir,subdir))\n", "valdir = os.path.join(\"GiRaFFE_Ccodes_validation\")\n", "cmd.mkdir(valdir)\n", "cmd.mkdir(os.path.join(valdir,subdir))\n", "\n", "generate_Afield_flux_function_files(out_dir,subdir,alpha_face,gamma_faceDD,beta_faceU,\\\n", " Valenciav_rU,B_rU,Valenciav_lU,B_lU,inputs_provided=True)\n", "Af.generate_Afield_flux_function_files(valdir,subdir,alpha_face,gamma_faceDD,beta_faceU,\\\n", " Valenciav_rU,B_rU,Valenciav_lU,B_lU,inputs_provided=True)\n", "\n", "print(\"Printing difference between original C code and this code...\")\n", "# Open the files to compare\n", "files = [\"RHSs/calculate_E_field_D0_right.h\",\n", " \"RHSs/calculate_E_field_D0_left.h\",\n", " \"RHSs/calculate_E_field_D1_right.h\",\n", " \"RHSs/calculate_E_field_D1_left.h\",\n", " \"RHSs/calculate_E_field_D2_right.h\",\n", " \"RHSs/calculate_E_field_D2_left.h\"]\n", "\n", "for file in files:\n", " print(\"Checking file \" + file)\n", " with open(os.path.join(valdir,file)) as file1, open(os.path.join(out_dir,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(out_dir,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": [ "\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-Induction_Equation.pdf](Tutorial-GiRaFFE_NRPy-Induction_Equation.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": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Notebook output to PDF is only supported on Linux systems, with pdflatex installed.\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-Afield_flux\")" ] } ], "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.4" } }, "nbformat": 4, "nbformat_minor": 2 }