convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.C
¶This module is currently under development
We will now use the cmdline_helper.py NRPy+ module to create the source directory within the IllinoisGRMHD
NRPy+ directory, if it does not exist yet.
# Step 0: Creation of the IllinoisGRMHD source directory
# Step 0a: Add NRPy's directory to the path
# https://stackoverflow.com/questions/16780014/import-file-from-parent-directory
import os,sys
nrpy_dir_path = os.path.join("..","..")
if nrpy_dir_path not in sys.path:
sys.path.append(nrpy_dir_path)
# Step 0b: Load up cmdline_helper and create the directory
import cmdline_helper as cmd
IGM_src_dir_path = os.path.join("..","src")
cmd.mkdir(IGM_src_dir_path)
# Step 0c: Create the output file path
outfile_path__ADM_to_BSSN__det_gammabar_eq_1__C = os.path.join(IGM_src_dir_path,"convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.C")
In this module we will explain the procedure used within IllinoisGRMHD
to compute the conformal metric, $\bar\gamma_{ij}$, and its inverse, $\bar\gamma^{ij}$, from the physical metric $\gamma_{ij}$. We will also describe how to compute $\phi$, the conformal factor, and $\psi\equiv e^{\phi}$. Finally, we also explain the procedure used to enforce the constraint $\bar\gamma = \det\left(\bar\gamma_{ij}\right) = 1$.
A note on notation: The notation used throughout the NRPy tutorial notebooks for the conformal metric is $\bar\gamma_{ij}$. However, in the literature the notation $\tilde\gamma_{ij}$ is also used to represent the conformal metric. It is important to note that in IllinoisGRMHD
we refer to the conformal metric using the latter notation, i.e. ${\rm gtij} := \tilde\gamma_{ij}$ and ${\rm gtupij} := \tilde\gamma^{ij}$. To keep the discussion consistent with the other notebooks, however, we will still present the discussion using the notation $\tilde\gamma_{ij} \to \bar\gamma_{ij}$. Bottom line, here $\tilde\gamma_{ij}$ and $\bar\gamma_{ij}$ represent exactly the same quantity, the conformal metric.
convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.C
[Back to top]¶We start by reading in the physical metric ${\rm gij\_physL} := \gamma_{ij}$ and computing its determinant
$$ \boxed{ \begin{align} \gamma = \det\left(\gamma_{ij}\right) &= \gamma_{xx}\gamma_{yy}\gamma_{zz} + \gamma_{xy}\gamma_{yz}\gamma_{xz} + \gamma_{xz}\gamma_{xy}\gamma_{yz}\\ &- \gamma_{xz}\gamma_{yy}\gamma_{xz} - \gamma_{xy}\gamma_{xy}\gamma_{zz} - \gamma_{xx}\gamma_{yz}\gamma_{yz} \end{align} }\ . $$Notice that we have used the fact that $\gamma_{ij}$ is symmetric above, i.e. $\gamma_{ij}=\gamma_{ji}$.
%%writefile $outfile_path__ADM_to_BSSN__det_gammabar_eq_1__C
#include <cmath>
#include <cstdio>
#include <cstdlib>
#ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER
#include "cctk.h"
#include "cctk_Parameters.h"
#endif
void IllinoisGRMHD_convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij
(const cGH *cctkGH,const int *cctk_lsh,
CCTK_REAL *gxx,CCTK_REAL *gxy,CCTK_REAL *gxz,CCTK_REAL *gyy,CCTK_REAL *gyz,CCTK_REAL *gzz,CCTK_REAL *alp,
CCTK_REAL *gtxx,CCTK_REAL *gtxy,CCTK_REAL *gtxz,CCTK_REAL *gtyy,CCTK_REAL *gtyz,CCTK_REAL *gtzz,
CCTK_REAL *gtupxx,CCTK_REAL *gtupxy,CCTK_REAL *gtupxz,CCTK_REAL *gtupyy,CCTK_REAL *gtupyz,CCTK_REAL *gtupzz,
CCTK_REAL *phi,CCTK_REAL *psi,CCTK_REAL *lapm1) {
#ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER
DECLARE_CCTK_PARAMETERS;
#endif
#pragma omp parallel for
for(int k=0;k<cctk_lsh[2];k++) for(int j=0;j<cctk_lsh[1];j++) for(int i=0;i<cctk_lsh[0];i++) {
int index=CCTK_GFINDEX3D(cctkGH,i,j,k);
CCTK_REAL gxx_physL=gxx[index];
CCTK_REAL gxy_physL=gxy[index];
CCTK_REAL gxz_physL=gxz[index];
CCTK_REAL gyy_physL=gyy[index];
CCTK_REAL gyz_physL=gyz[index];
CCTK_REAL gzz_physL=gzz[index];
/**********************************************************************
* Compute \tilde{\gamma_{ij}}, phi, and psi (BSSN) from g_{ij} (ADM) *
**********************************************************************/
CCTK_REAL gijdet = gxx_physL * gyy_physL * gzz_physL + gxy_physL * gyz_physL * gxz_physL + gxz_physL * gxy_physL * gyz_physL
- gxz_physL * gyy_physL * gxz_physL - gxy_physL * gxy_physL * gzz_physL - gxx_physL * gyz_physL * gyz_physL;
gijdet = fabs(gijdet);
Writing ../src/convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.C
The BSSN formalism relates the physical metric, $\gamma_{ij}$, to a conformal metric, $\bar\gamma{ij}$, via the relation
$$ \boxed{\bar\gamma_{ij} = e^{-4\phi}\gamma_{ij}}\ , $$with the constraint that $\bar\gamma \equiv \det\left(\bar\gamma_{ij}\right) = 1$. This immediately implies that
$$ e^{-12\phi} \gamma = \bar\gamma = 1 \implies -12\phi = \log\left(\frac{1}{\gamma}\right) = - \log\gamma \implies \boxed{\phi = \frac{1}{12}\log\gamma}\ . $$Useful quantities to compute are
$$ \boxed{ \begin{align} \psi &\equiv e^{\phi}\\ \psi^{-4} &\equiv e^{-4\phi} \end{align} }\ . $$All the boxed quantities above are evaluated below.
%%writefile -a $outfile_path__ADM_to_BSSN__det_gammabar_eq_1__C
CCTK_REAL phiL = (1.0/12.0) * log(gijdet);
CCTK_REAL psiL = exp(phiL);
CCTK_REAL Psim4 = 1.0/(psiL*psiL*psiL*psiL);
CCTK_REAL gtxxL = gxx_physL*Psim4;
CCTK_REAL gtxyL = gxy_physL*Psim4;
CCTK_REAL gtxzL = gxz_physL*Psim4;
CCTK_REAL gtyyL = gyy_physL*Psim4;
CCTK_REAL gtyzL = gyz_physL*Psim4;
CCTK_REAL gtzzL = gzz_physL*Psim4;
Appending to ../src/convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.C
The BSSN formalism evolves the determinant of the conformal metric through $\partial_{t}\bar\gamma = 0$. Since initially $\bar\gamma=1$, one would expect that $\bar\gamma=1$ at all times. However, numerical errors can cause the value of the determinant of the conformal metric to drift away from unity. To remedy this, we enforce the constraint $\bar\gamma=1$ by transforming the conformal metric according to the algebraic condition (cf. eq. 53 of Ruchlin, Etienne, and Baumgarte (2018) with $\hat\gamma = 1$).
$$ \boxed{\bar{\gamma}_{ij} \to \left(\frac{1}{\bar{\gamma}}\right)^{1/3} \bar{\gamma}_{ij}}\ . $$We then start by evaluating
$$ \boxed{ \begin{align} \bar\gamma = \det\left(\bar\gamma_{ij}\right) &= \bar\gamma_{xx}\bar\gamma_{yy}\bar\gamma_{zz} + \bar\gamma_{xy}\bar\gamma_{yz}\bar\gamma_{xz} + \bar\gamma_{xz}\bar\gamma_{xy}\bar\gamma_{yz}\\ &- \bar\gamma_{xz}\bar\gamma_{yy}\bar\gamma_{xz} - \bar\gamma_{xy}\bar\gamma_{xy}\bar\gamma_{zz} - \bar\gamma_{xx}\bar\gamma_{yz}\bar\gamma_{yz} \end{align} }\ . $$In the code below we also define
$$ \boxed{{\rm gtijdet\_Fm1o3} = \left(\frac{1}{\bar{\gamma}}\right)^{1/3}}\ . $$%%writefile -a $outfile_path__ADM_to_BSSN__det_gammabar_eq_1__C
/*********************************
* Apply det gtij = 1 constraint *
*********************************/
CCTK_REAL gtijdet = gtxxL * gtyyL * gtzzL + gtxyL * gtyzL * gtxzL + gtxzL * gtxyL * gtyzL -
gtxzL * gtyyL * gtxzL - gtxyL * gtxyL * gtzzL - gtxxL * gtyzL * gtyzL;
CCTK_REAL gtijdet_Fm1o3 = fabs(1.0/cbrt(gtijdet));
gtxxL = gtxxL * gtijdet_Fm1o3;
gtxyL = gtxyL * gtijdet_Fm1o3;
gtxzL = gtxzL * gtijdet_Fm1o3;
gtyyL = gtyyL * gtijdet_Fm1o3;
gtyzL = gtyzL * gtijdet_Fm1o3;
gtzzL = gtzzL * gtijdet_Fm1o3;
if(gtijdet<0.0) {
#ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER
CCTK_VWarn(CCTK_WARN_ALERT,__LINE__, __FILE__, CCTK_THORNSTRING,
#else
fprintf(stderr,
#endif
"WARNING: det[3-metric]<0.0 at point %d %d %d | cctk_lsh: %d %d %d. Hopefully this is occurring in gz's! gtij_phys = %.2e %.2e %.2e %.2e %.2e %.2e gtij_new = %.2e %.2e %.2e %.2e %.2e %.2e | gijdet = %.2e | gtijdet = %.2e\n",
i,j,k,cctk_lsh[0],cctk_lsh[1],cctk_lsh[2],gxx_physL,gxy_physL,gxz_physL,gyy_physL,gyz_physL,gzz_physL,gtxxL,gtxyL,gtxzL,gtyyL,gtyzL,gtzzL,-gijdet,gtijdet);
}
Appending to ../src/convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.C
Finally, we update the gridfunctions for $\phi$, $\psi$, $\bar\gamma_{ij}$, $\gamma_{ij}$, and $\bar\gamma^{ij}$, respectively, in the code below.
If $\mathcal{M}$ is a $3\times3$ symmetric matrix and $\det(\mathcal{M})\neq0$, recall that we have
$$ \mathcal{M}^{-1} = \begin{pmatrix} a & b & c\\ b & d & e\\ c & e & f \end{pmatrix}^{-1} = \frac{1}{\det(\mathcal{M})} \begin{pmatrix} A & B & C\\ B & D & E\\ C & E & F \end{pmatrix}\ , $$where
$$ \boxed{ \begin{align} A &= \ \ \left(df - ee\right)\\ B &= -\left(bf - ce\right)\\ C &= \ \ \left(be - cd\right)\\ D &= \ \ \left(af - cc\right)\\ E &= -\left(ae - bc\right)\\ F &= \ \ \left(ad - bb\right) \end{align} }\ . $$Notice that if we replace $\mathcal{M}\to\bar\gamma_{ij}$, then we have also $\det(\mathcal{M})\to1$.
%%writefile -a $outfile_path__ADM_to_BSSN__det_gammabar_eq_1__C
CCTK_REAL Psi4 = psiL*psiL*psiL*psiL;
/*****************************************
* Set all the needed BSSN gridfunctions *
*****************************************/
phi[index] = phiL;
psi[index] = psiL;
lapm1[index] = alp[index] - 1.0;
gtxx[index] = gtxxL;
gtxy[index] = gtxyL;
gtxz[index] = gtxzL;
gtyy[index] = gtyyL;
gtyz[index] = gtyzL;
gtzz[index] = gtzzL;
gxx[index] = gtxxL*Psi4;
gxy[index] = gtxyL*Psi4;
gxz[index] = gtxzL*Psi4;
gyy[index] = gtyyL*Psi4;
gyz[index] = gtyzL*Psi4;
gzz[index] = gtzzL*Psi4;
gtupxx[index] = ( gtyyL * gtzzL - gtyzL * gtyzL );
gtupxy[index] = - ( gtxyL * gtzzL - gtyzL * gtxzL );
gtupxz[index] = ( gtxyL * gtyzL - gtyyL * gtxzL );
gtupyy[index] = ( gtxxL * gtzzL - gtxzL * gtxzL );
gtupyz[index] = - ( gtxxL * gtyzL - gtxyL * gtxzL );
gtupzz[index] = ( gtxxL * gtyyL - gtxyL * gtxyL );
}
}
Appending to ../src/convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.C
# Verify if the code generated by this tutorial module
# matches the original IllinoisGRMHD source code
# First download the original IllinoisGRMHD source code
import urllib
from os import path
original_IGM_file_url = "https://bitbucket.org/zach_etienne/wvuthorns/raw/5611b2f0b17135538c9d9d17c7da062abe0401b6/IllinoisGRMHD/src/convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.C"
original_IGM_file_name = "convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij-original.C"
original_IGM_file_path = os.path.join(IGM_src_dir_path,original_IGM_file_name)
# Then download the original IllinoisGRMHD source code
# We try it here in a couple of ways in an attempt to keep
# the code more portable
try:
original_IGM_file_code = urllib.request.urlopen(original_IGM_file_url).read().decode("utf-8")
# Write down the file the original IllinoisGRMHD source code
with open(original_IGM_file_path,"w") as file:
file.write(original_IGM_file_code)
except:
try:
original_IGM_file_code = urllib.urlopen(original_IGM_file_url).read().decode("utf-8")
# Write down the file the original IllinoisGRMHD source code
with open(original_IGM_file_path,"w") as file:
file.write(original_IGM_file_code)
except:
# If all else fails, hope wget does the job
!wget -O $original_IGM_file_path $original_IGM_file_url
# Perform validation
Validation__ADM_to_BSSN__det_gammabar_eq_1__C = !diff $original_IGM_file_path $outfile_path__ADM_to_BSSN__det_gammabar_eq_1__C
if Validation__ADM_to_BSSN__det_gammabar_eq_1__C == []:
# If the validation passes, we do not need to store the original IGM source code file
!rm $original_IGM_file_path
print("Validation test for convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.C: PASSED!")
else:
# If the validation fails, we keep the original IGM source code file
print("Validation test for convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.C: FAILED!")
# We also print out the difference between the code generated
# in this tutorial module and the original IGM source code
print("Diff:")
for diff_line in Validation__ADM_to_BSSN__det_gammabar_eq_1__C:
print(diff_line)
Validation test for convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.C: FAILED! Diff: 1d0 < #include "cctk.h" 4a4,5 > #ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER > #include "cctk.h" 5a7 > #endif 13c15,20 < DECLARE_CCTK_PARAMETERS; --- > > #ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER > DECLARE_CCTK_PARAMETERS; > #endif > > 31a39 > 42c50,51 < --- > > 46c55 < CCTK_REAL gtijdet = gtxxL * gtyyL * gtzzL + gtxyL * gtyzL * gtxzL + gtxzL * gtxyL * gtyzL - --- > CCTK_REAL gtijdet = gtxxL * gtyyL * gtzzL + gtxyL * gtyzL * gtxzL + gtxzL * gtxyL * gtyzL - 52,55c61,64 < gtxyL = gtxyL * gtijdet_Fm1o3; < gtxzL = gtxzL * gtijdet_Fm1o3; < gtyyL = gtyyL * gtijdet_Fm1o3; < gtyzL = gtyzL * gtijdet_Fm1o3; --- > gtxyL = gtxyL * gtijdet_Fm1o3; > gtxzL = gtxzL * gtijdet_Fm1o3; > gtyyL = gtyyL * gtijdet_Fm1o3; > gtyzL = gtyzL * gtijdet_Fm1o3; 58,60c67,76 < if(gtijdet<0.0) { CCTK_VWarn(CCTK_WARN_ALERT,__LINE__, __FILE__, CCTK_THORNSTRING, < "WARNING: det[3-metric]<0.0 at point %d %d %d | cctk_lsh: %d %d %d. Hopefully this is occurring in gz's! gtij_phys = %.2e %.2e %.2e %.2e %.2e %.2e gtij_new = %.2e %.2e %.2e %.2e %.2e %.2e | gijdet = %.2e | gtijdet = %.2e", < i,j,k,cctk_lsh[0],cctk_lsh[1],cctk_lsh[2],gxx_physL,gxy_physL,gxz_physL,gyy_physL,gyz_physL,gzz_physL,gtxxL,gtxyL,gtxzL,gtyyL,gtyzL,gtzzL,-gijdet,gtijdet); } --- > if(gtijdet<0.0) { > #ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER > CCTK_VWarn(CCTK_WARN_ALERT,__LINE__, __FILE__, CCTK_THORNSTRING, > #else > fprintf(stderr, > #endif > "WARNING: det[3-metric]<0.0 at point %d %d %d | cctk_lsh: %d %d %d. Hopefully this is occurring in gz's! gtij_phys = %.2e %.2e %.2e %.2e %.2e %.2e gtij_new = %.2e %.2e %.2e %.2e %.2e %.2e | gijdet = %.2e | gtijdet = %.2e\n", > i,j,k,cctk_lsh[0],cctk_lsh[1],cctk_lsh[2],gxx_physL,gxy_physL,gxz_physL,gyy_physL,gyz_physL,gzz_physL,gtxxL,gtxyL,gtxzL,gtyyL,gtyzL,gtzzL,-gijdet,gtijdet); > } > 92a109 >
The following code cell converts this Jupyter notebook into a proper, clickable $\LaTeX$-formatted PDF file. After the cell is successfully run, the generated PDF may be found in the root NRPy+ tutorial directory, with filename Tutorial-IllinoisGRMHD__convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.pdf (Note that clicking on this link may not work; you may need to open the PDF file through another means).
latex_nrpy_style_path = os.path.join(nrpy_dir_path,"latex_nrpy_style.tplx")
#!jupyter nbconvert --to latex --template $latex_nrpy_style_path --log-level='WARN' Tutorial-IllinoisGRMHD__convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.ipynb
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.tex
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.tex
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.tex
!rm -f Tut*.out Tut*.aux Tut*.log