Start-to-Finish Example: Scalar Field Collapse with regrids¶
Authors: Leonardo Werneck & Zachariah B. Etienne¶
This module sets up spherically symmetric, time-symmetric initial data for a scalar field collapse in Spherical coordinates, as documented in this NRPy+ module (the initial data is shown to satisfy the Hamiltonian constraint in this tutorial module), which is then evolved forward in time. The tutorial notebook is very similar to the tutorial notebook on gravitational collapse of a massless scalar field, except that we will use regridding in order to adaptively increase the resolution of the run as we follow the critical solution¶
WARNING: The results in this tutorial notebook (see plot at the bottom) *do not* reproduce the ones in Werneck et al. (2021). The fine-tuning of the critical solution is extremely sensitive to any changes in the code, including round-off level disagreements which can be caused by something as simple as using different compilers. The code below chooses a different evolved conformal factor than the original NRPyCritCol code and while this choice of conformal factor results in a formulation of the BSSN equations which is equivalent to the one used by the original, it leads to round-off level disagreements. This ultimately changes the value of the critical amplitude, $\eta_{*}$, and, consequently, $\eta_{\rm weak}$ and $\eta_{\rm strong}$. If you want to reproduce Fig. 3 from the paper, then we suggest using an older version of NRPyCritCol code, which is the one that was used to generate that plot.
Exercise to the reader: Obtain a fine-tuning, $\delta\eta\equiv 1 - \eta_{\rm weak}/\eta_{\rm strong}$, of order $10^{-10}$ or better using the code below. In order to do this you will need to add new regrids and/or adjust the existing regridding parameters. You can use the regridding times/parameters used by the original NRPyCritCol code as a guide.
Results from this tutorial notebook have been used in the paper Werneck et al. (2021)¶
The entire algorithm is outlined below, with NRPy+-based components highlighted in green.
- Allocate memory for gridfunctions, including temporary storage for the RK4 time integration.
- Set gridfunction values to initial data.
- Evolve the system forward in time using RK4 time integration. At each RK4 substep, do the following:
- Evaluate BSSN RHS expressions.
- Apply singular, curvilinear coordinate boundary conditions a la the SENR/NRPy+ paper
- Apply constraints on conformal 3-metric: $\det{\bar{\gamma}_{ij}}=\det{\hat{\gamma}_{ij}}$
- At the end of each iteration in time, output the Hamiltonian constraint violation.
- Repeat above steps at two numerical resolutions to confirm convergence to zero.
References¶
- Akbarian & Choptuik (2015): Useful to understand the theoretical framework
- Baumgarte (2018): Useful to understand the theoretical framework
- Baumgarte & Shapiro's Numerical Relativity: Section 6.2.2 - Useful to understand how to solve the Hamiltonian constraint
- Werneck et al. (2021): Mathematical formulation and numerical implementation of NRPyCritCol
Table of Contents¶
$$\label{toc}$$- Step 1 Set core NRPy+ parameters for numerical grids and reference metric
- Step 2 Set up ADM initial data for the Scalar Field
- Step 3 Convert ADM initial data to BSSN-in-curvilinear coordinates
- Step 4 Output C code for BSSN spacetime evolution
- Step 4.a Set up the BSSN and ScalarField right-hand-side (RHS) expressions, and add the rescaled $T^{\mu\nu}$ source terms
- Step 4.b Output the Hamiltonian constraint
- Step 4.c Enforce conformal 3-metric $\det{\bar{\gamma}_{ij}}=\det{\hat{\gamma}_{ij}}$ constraint
- Step 4.d Generate C code kernels for BSSN expressions, in parallel if possible
- Step 4.e Output C codes needed for declaring and setting Cparameters; also set
free_parameters.h
- Step 5 Set up boundary condition functions for chosen singular, curvilinear coordinate system
- Step 6 The main C code:
ScalarFieldCollapse_Playground.c - Step 7 Visualization: self-similar behavior of the lapse function
- Step 8 Output this module as $\LaTeX$-formatted PDF file
In [1]:
# Step P1: Import needed NRPy+ core modules:
import shutil, os, sys # Standard Python modules for multiplatform OS-level functions
sys.path.append("..")
from outputC import lhrh,outCfunction # NRPy+: Core C code output module
import NRPy_param_funcs as par # NRPy+: Parameter interface
import finite_difference as fin # NRPy+: Finite difference C code generation module
import grid as gri # NRPy+: Functions having to do with numerical grids
import indexedexp as ixp # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support
import reference_metric as rfm # NRPy+: Reference metric support
import cmdline_helper as cmd # NRPy+: Multi-platform Python command-line interface
# Step P2: Create C code output directory:
Ccodesdir = os.path.join("BSSN_ScalarFieldCollapse_Ccodes")
# First remove C code output directory if it exists
# Courtesy https://stackoverflow.com/questions/303200/how-do-i-remove-delete-a-folder-that-is-not-empty
# !rm -r ScalarWaveCurvilinear_Playground_Ccodes
shutil.rmtree(Ccodesdir, ignore_errors=True)
# Then create a fresh directory
cmd.mkdir(Ccodesdir)
# Step P3: Create executable output directory:
outdir = os.path.join(Ccodesdir,"output")
cmd.mkdir(outdir)
# Step 1: Set the spatial dimension parameter
# to three this time, and then read
# the parameter as DIM.
par.set_parval_from_str("grid::DIM",3)
DIM = par.parval_from_str("grid::DIM")
# Step 2: Set some core parameters, including CoordSystem MoL timestepping algorithm,
# FD order, floating point precision, and CFL factor:
# Choices are: Spherical, SinhSpherical, SinhSphericalv2, Cylindrical, SinhCylindrical,
# SymTP, SinhSymTP
CoordSystem = "SinhSpherical"
# Step 2.a: Set defaults for Coordinate system parameters.
# These are perhaps the most commonly adjusted parameters,
# so we enable modifications at this high level.
domain_size = 64
# sinh_width sets the default value for:
# * SinhSpherical's params.SINHW
# * SinhCylindrical's params.SINHW{RHO,Z}
# * SinhSymTP's params.SINHWAA
sinh_width = 0.2 # If Sinh* coordinates chosen
# sinhv2_const_dr sets the default value for:
# * SinhSphericalv2's params.const_dr
# * SinhCylindricalv2's params.const_d{rho,z}
sinhv2_const_dr = 0.05# If Sinh*v2 coordinates chosen
# SymTP_bScale sets the default value for:
# * SinhSymTP's params.bScale
SymTP_bScale = 0.5 # If SymTP chosen
# Step 2.b: Set the order of spatial and temporal derivatives;
# the core data type, and the CFL factor.
# RK_method choices include: Euler, "RK2 Heun", "RK2 MP", "RK2 Ralston", RK3, "RK3 Heun", "RK3 Ralston",
# SSPRK3, RK4, DP5, DP5alt, CK5, DP6, L6, DP8
RK_method = "RK4"
FD_order = 4 # Finite difference order: even numbers only, starting with 2. 12 is generally unstable
REAL = "double" # Best to use double here.
CFL_FACTOR= 0.5
# Set the lapse & shift conditions
LapseCondition = "OnePlusLog"
ShiftCondition = "GammaDriving2ndOrder_Covariant"
# Step 3: Generate Runge-Kutta-based (RK-based) timestepping code.
# As described above the Table of Contents, this is a 3-step process:
# 3.A: Evaluate RHSs (RHS_string)
# 3.B: Apply boundary conditions (post_RHS_string, pt 1)
# 3.C: Enforce det(gammabar) = det(gammahat) constraint (post_RHS_string, pt 2)
import MoLtimestepping.C_Code_Generation as MoL
from MoLtimestepping.RK_Butcher_Table_Dictionary import Butcher_dict
RK_order = Butcher_dict[RK_method][1]
cmd.mkdir(os.path.join(Ccodesdir,"MoLtimestepping/"))
MoL.MoL_C_Code_Generation(RK_method,
RHS_string = """
Ricci_eval(&rfmstruct, ¶ms, RK_INPUT_GFS, auxevol_gfs);
rhs_eval(&rfmstruct, ¶ms, auxevol_gfs, RK_INPUT_GFS, RK_OUTPUT_GFS);""",
post_RHS_string = """
apply_bcs_curvilinear(¶ms, &bcstruct, NUM_EVOL_GFS, evol_gf_parity, RK_OUTPUT_GFS);
enforce_detgammahat_constraint(&rfmstruct, ¶ms, RK_OUTPUT_GFS);\n""",
outdir = os.path.join(Ccodesdir,"MoLtimestepping/"))
# Step 4: Set the coordinate system for the numerical grid
par.set_parval_from_str("reference_metric::CoordSystem",CoordSystem)
rfm.reference_metric() # Create ReU, ReDD needed for rescaling B-L initial data, generating BSSN RHSs, etc.
# Step 5: Set the finite differencing order to FD_order (set above).
par.set_parval_from_str("finite_difference::FD_CENTDERIVS_ORDER", FD_order)
# Step 6: Copy SIMD/SIMD_intrinsics.h to $Ccodesdir/SIMD/SIMD_intrinsics.h
cmd.mkdir(os.path.join(Ccodesdir,"SIMD"))
shutil.copy(os.path.join("..","SIMD","SIMD_intrinsics.h"),os.path.join(Ccodesdir,"SIMD"))
# Step 7: Impose spherical symmetry by demanding that all
# derivatives in the angular directions vanish
par.set_parval_from_str("indexedexp::symmetry_axes","12")
Step 1.a: Writing the regridding C code [Back to top]¶
$$\label{regridding_ccode}$$NRPyCritCol was designed to study critical phenomena in the context of gravitational collapse. If the scalar field's initial profile is given by, for example,
$$ \varphi(r) = \eta\exp\left(\frac{r^{2}}{\sigma^{2}}\right), $$one finds that there are values of $\eta$ that result in complete dispersal of the scalar field and values of $\eta$ that result in gravitational collapse of the scalar field. Subcritical The result of the first case is flat spacetime, while the result of the second case is the formation of a black hole. It is then conjectured that there is a critical value of $\eta$, denoted $\eta_{*}$, which separates subcritical ($\eta<\eta_{*}$) runs from supercritical ($\eta>\eta_{*}$) runs.
As we increase the level of fine-tuning, i.e. the closer we get to $\eta_{*}$, the more structure we observe in our solution. The critical solution, for which $\eta=\eta_{*}$, is characterized by the self-similar behavior that can be observed in spacetime and scalar fields quantity. The same solution repeats itself on ever smaller spatiotemporal scales, requiring more and more resolution in order to be resolved.
Starting a run with very high grid resolution, which implies very small time steps (remember the CFL condition), should allow us to resolve the structures that appear in the solution, but increases the runtime of the code to a point where studying critical phenomena becomes impractical. To remedy this, we adopt a different strategy: regridding.
We begin our run with a grid of modest resolution, enough to resolve the first structures appearing during the time evolution. After some time, we define a new grid which contains higher resolution than the original one. We then interpolate our solution from the original grid onto the new, higher resolution grid, and continue the evolution using the higher resolution grid. This characterizes 1 regrid.
Our numerical grids are defined by our coordinate system $(\mathtt{xx0},\mathtt{xx1},\mathtt{xx2})$, which is related to spherical coordinates $(r,\theta,\phi)$ via
$$ \begin{aligned} r &= \mathcal{A}\frac{\sinh\left(\mathtt{xx0}/w\right)}{\sinh\left(1/w\right)},\\ \theta &= \mathtt{xx1},\\ \phi &= \mathtt{xx2}. \end{aligned} $$We increase the resolution around $r\sim0$ by doing one of two things:
- Decreasing the value of $w$,
- Decreasing the value of $\mathcal{A}$.
Doing 1. increases the sampling density along the radial direction, while doing 2. brings the outer boundary $r_{\rm max} = \mathcal{A}$ closer to the center of the simulation while also increasing the resolution near the origin.
The code below is written to interface with the regridding function and struct defined in ScalarField/NRPyCritCol_regridding.h. We define:
- The times at which regrids will take place,
- Whether we want to change the value of $w$ ("sinhW") or $\mathcal{A}$ ("sinhA"),
- The new value of $w$ or $\mathcal{A}$ on the new grid.
In [2]:
# Define a named tuple to store the regrid information
from collections import namedtuple
regrid_tuple = namedtuple("regrid_tuple","time type new_parameter_value")
# Set up the regridding times, types, and the new parameter value
regrid_info_list = [regrid_tuple(5.100,"sinhW",0.19),
regrid_tuple(5.125,"sinhW",0.18),
regrid_tuple(5.150,"sinhW",0.17),
regrid_tuple(5.175,"sinhW",0.16),
regrid_tuple(5.200,"sinhW",0.15),
regrid_tuple(5.225,"sinhW",0.14),
regrid_tuple(5.250,"sinhW",0.13),
regrid_tuple(5.275,"sinhW",0.12),
regrid_tuple(5.300,"sinhW",0.11),
regrid_tuple(5.325,"sinhW",0.10),
regrid_tuple(6.460,"sinhA",56.0),
regrid_tuple(6.465,"sinhA",48.0),
regrid_tuple(6.470,"sinhA",40.0),
regrid_tuple(6.475,"sinhA",32.0),
regrid_tuple(6.480,"sinhA",24.0),
regrid_tuple(6.485,"sinhA",16.0)]
# Write the C code for it
with open(os.path.join(Ccodesdir,"set_regrid_struct.h"),"w") as file:
file.write(""" {
int num_regrids = 0;
""")
counter = 0
for regrid_params in regrid_info_list:
counter += 1
file.write("""
// Regrid #%d
if( num_regrids < max_number_of_regrids ) {
regrid_params.regrid_time[num_regrids] = %.15e;
regrid_params.new_ampl_or_sinhW[num_regrids] = %.15e;
regrid_params.regrid_key[num_regrids] = %s;
regrid_params.regrid_counter[num_regrids] = num_regrids;
num_regrids++;
}
"""%(counter,regrid_params.time,regrid_params.new_parameter_value,
"sinhw_regrid" if regrid_params.type == "sinhW" else "ampl_regrid"))
file.write(r"""
// Sanity check
if( num_regrids != max_number_of_regrids ) {
fprintf(stderr,"ERROR: Expected %d regrids but got %d\n",max_number_of_regrids,num_regrids);
exit(1);
}
// Print regridding information
for(int i=0;i<max_number_of_regrids;i++) {
printf("Regrid #%02d: %s, %.4lf, %.4lf\n",
i+1,
regrid_params.regrid_key[i] == sinhw_regrid ? "sinhW" : "sinhA",
regrid_params.regrid_time[i],
regrid_params.new_ampl_or_sinhW[i]);
}
}
""")
print("Wrote to file '%s'"%os.path.join(Ccodesdir,"set_regrid_struct.h"))
Wrote to file 'BSSN_ScalarFieldCollapse_Ccodes/set_regrid_struct.h'
Step 1.b: Output needed C code for finding the minimum proper distance between grid points, needed for CFL-limited timestep [Back to top]¶
$$\label{cfl}$$In order for our explicit-timestepping numerical solution to the scalar wave equation to be stable, it must satisfy the CFL condition: $$ \Delta t \le \frac{\min(ds_i)}{c}, $$ where $c$ is the wavespeed, and $$ds_i = h_i \Delta x^i$$ is the proper distance between neighboring gridpoints in the $i$th direction (in 3D, there are 3 directions), $h_i$ is the $i$th reference metric scale factor, and $\Delta x^i$ is the uniform grid spacing in the $i$th direction:
In [3]:
# Output the find_timestep() function to a C file.
rfm.out_timestep_func_to_file(os.path.join(Ccodesdir,"find_timestep.h"))
Step 2: Set up ADM initial data for the Scalar Field [Back to top]¶
$$\label{initial_data}$$As documented in the scalar field Gaussian pulse initial data NRPy+ tutorial notebook, we will now set up the scalar field initial data, storing the densely-sampled result to file.
The initial data function ScalarField_InitialData requires SciPy, so let's make sure it's installed.
In [4]:
!pip install scipy numpy > /dev/null
Next call the ScalarField_InitialData() function from the ScalarField/ScalarField_InitialData.py NRPy+ module (see the tutorial notebook).
In [5]:
# Step 2.a: Import necessary Python and NRPy+ modules
import ScalarField.ScalarField_InitialData as sfid
# Step 2.b: Set the initial data parameters
ID_Family = "Gaussian_pulse"
pulse_center = 0
pulse_width = 1
Nr = 30000
rmax = domain_size*1.1
# Step 2.c: Generate the initial data
eta_weak = 0.303326061
sfid.ScalarField_InitialData(os.path.join(outdir,"SFID_weak.txt"),ID_Family,
eta_weak,pulse_center,pulse_width,Nr,rmax)
eta_strong = 0.303326062
sfid.ScalarField_InitialData(os.path.join(outdir,"SFID_strong.txt"),ID_Family,
eta_strong,pulse_center,pulse_width,Nr,rmax)
sfid.NRPy_param_funcs_register_C_functions_and_NRPy_basic_defines(Ccodesdir=Ccodesdir)
Generated the ADM initial data for the gravitational collapse
of a massless scalar field in Spherical coordinates.
Type of initial condition: Scalar field: "Gaussian" Shell
ADM quantities: Time-symmetric
Lapse condition: Pre-collapsed
Parameters: amplitude = 0.303326061,
center = 0,
width = 1,
domain size = 70.4,
number of points = 30000,
Initial data file = BSSN_ScalarFieldCollapse_Ccodes/output/SFID_weak.txt.
Generated the ADM initial data for the gravitational collapse
of a massless scalar field in Spherical coordinates.
Type of initial condition: Scalar field: "Gaussian" Shell
ADM quantities: Time-symmetric
Lapse condition: Pre-collapsed
Parameters: amplitude = 0.303326062,
center = 0,
width = 1,
domain size = 70.4,
number of points = 30000,
Initial data file = BSSN_ScalarFieldCollapse_Ccodes/output/SFID_strong.txt.
Output C function ID_scalarfield_ADM_quantities() to file BSSN_ScalarFieldCollapse_Ccodes/ID_scalarfield_ADM_quantities.h
Output C function ID_scalarfield_spherical() to file BSSN_ScalarFieldCollapse_Ccodes/ID_scalarfield_spherical.h
Output C function ID_scalarfield_xx0xx1xx2_to_BSSN_xx0xx1xx2() to file BSSN_ScalarFieldCollapse_Ccodes/ID_scalarfield_xx0xx1xx2_to_BSSN_xx0xx1xx2.h
Output C function ID_scalarfield() to file BSSN_ScalarFieldCollapse_Ccodes/ID_scalarfield.h
Step 3: Convert ADM initial data to BSSN-in-curvilinear coordinates [Back to top]¶
$$\label{adm_id_spacetime}$$This is an automated process, taken care of by BSSN.ADM_Numerical_Spherical_or_Cartesian_to_BSSNCurvilinear, and documented in this tutorial notebook.
In [6]:
import BSSN.ADM_Numerical_Spherical_or_Cartesian_to_BSSNCurvilinear as AtoBnum
AtoBnum.Convert_Spherical_or_Cartesian_ADM_to_BSSN_curvilinear("Spherical","ID_scalarfield_ADM_quantities",
Ccodesdir=Ccodesdir,loopopts="")
Output C function ID_BSSN_lambdas() to file BSSN_ScalarFieldCollapse_Ccodes/ID_BSSN_lambdas.h Output C function ID_ADM_xx0xx1xx2_to_BSSN_xx0xx1xx2__ALL_BUT_LAMBDAs() to file BSSN_ScalarFieldCollapse_Ccodes/ID_ADM_xx0xx1xx2_to_BSSN_xx0xx1xx2__ALL_BUT_LAMBDAs.h Output C function ID_BSSN__ALL_BUT_LAMBDAs() to file BSSN_ScalarFieldCollapse_Ccodes/ID_BSSN__ALL_BUT_LAMBDAs.h
Step 4: Output C code for BSSN spacetime evolution [Back to top]¶
$$\label{bssn}$$Step 4.a: Set up the BSSN and ScalarField right-hand-side (RHS) expressions, and add the rescaled $T^{\mu\nu}$ source terms [Back to top]¶
$$\label{bssnrhs}$$BSSN.BSSN_RHSs() sets up the RHSs assuming a spacetime vacuum: $T^{\mu\nu}=0$. (This might seem weird, but remember that, for example, spacetimes containing only single or binary black holes are vacuum spacetimes.) Here, using the BSSN.BSSN_stress_energy_source_terms (tutorial) NRPy+ module, we add the $T^{\mu\nu}$ source terms to these equations.
In [7]:
import time
import BSSN.BSSN_RHSs as rhs
import BSSN.BSSN_gauge_RHSs as gaugerhs
par.set_parval_from_str("BSSN.BSSN_gauge_RHSs::LapseEvolutionOption", LapseCondition)
par.set_parval_from_str("BSSN.BSSN_gauge_RHSs::ShiftEvolutionOption", ShiftCondition)
print("Generating symbolic expressions for BSSN RHSs...")
start = time.time()
# Enable rfm_precompute infrastructure, which results in
# BSSN RHSs that are free of transcendental functions,
# even in curvilinear coordinates, so long as
# ConformalFactor is set to "W" (default).
cmd.mkdir(os.path.join(Ccodesdir,"rfm_files/"))
par.set_parval_from_str("reference_metric::enable_rfm_precompute","True")
par.set_parval_from_str("reference_metric::rfm_precompute_Ccode_outdir",os.path.join(Ccodesdir,"rfm_files/"))
# Evaluate BSSN + BSSN gauge RHSs with rfm_precompute enabled:
import BSSN.BSSN_quantities as Bq
par.set_parval_from_str("BSSN.BSSN_quantities::LeaveRicciSymbolic","True")
rhs.BSSN_RHSs()
# Evaluate the Scalar Field RHSs
import ScalarField.ScalarField_RHSs as sfrhs
sfrhs.ScalarField_RHSs()
# Compute ScalarField T^{\mu\nu}
# Compute the scalar field energy-momentum tensor
import ScalarField.ScalarField_Tmunu as sfTmunu
sfTmunu.ScalarField_Tmunu()
T4UU = sfTmunu.T4UU
import BSSN.BSSN_stress_energy_source_terms as Bsest
Bsest.BSSN_source_terms_for_BSSN_RHSs(T4UU)
rhs.trK_rhs += Bsest.sourceterm_trK_rhs
for i in range(DIM):
# Needed for Gamma-driving shift RHSs:
rhs.Lambdabar_rhsU[i] += Bsest.sourceterm_Lambdabar_rhsU[i]
# Needed for BSSN RHSs:
rhs.lambda_rhsU[i] += Bsest.sourceterm_lambda_rhsU[i]
for j in range(DIM):
rhs.a_rhsDD[i][j] += Bsest.sourceterm_a_rhsDD[i][j]
gaugerhs.BSSN_gauge_RHSs()
# We use betaU as our upwinding control vector:
Bq.BSSN_basic_tensors()
betaU = Bq.betaU
import BSSN.Enforce_Detgammahat_Constraint as EGC
enforce_detg_constraint_symb_expressions = EGC.Enforce_Detgammahat_Constraint_symb_expressions()
# Next compute Ricci tensor
par.set_parval_from_str("BSSN.BSSN_quantities::LeaveRicciSymbolic","False")
Bq.RicciBar__gammabarDD_dHatD__DGammaUDD__DGammaU()
# Now register the Hamiltonian as a gridfunction.
H = gri.register_gridfunctions("AUX","H")
# Then define the Hamiltonian constraint and output the optimized C code.
import BSSN.BSSN_constraints as bssncon
bssncon.BSSN_constraints(add_T4UUmunu_source_terms=False)
Bsest.BSSN_source_terms_for_BSSN_constraints(T4UU)
bssncon.H += Bsest.sourceterm_H
# Add Kreiss-Oliger dissipation
diss_strength = par.Cparameters("REAL","ScalarFieldCollapse",["diss_strength"],0.1)
alpha_dKOD = ixp.declarerank1("alpha_dKOD")
cf_dKOD = ixp.declarerank1("cf_dKOD")
trK_dKOD = ixp.declarerank1("trK_dKOD")
sf_dKOD = ixp.declarerank1("sf_dKOD")
sfM_dKOD = ixp.declarerank1("sfM_dKOD")
betU_dKOD = ixp.declarerank2("betU_dKOD","nosym")
vetU_dKOD = ixp.declarerank2("vetU_dKOD","nosym")
lambdaU_dKOD = ixp.declarerank2("lambdaU_dKOD","nosym")
aDD_dKOD = ixp.declarerank3("aDD_dKOD","sym01")
hDD_dKOD = ixp.declarerank3("hDD_dKOD","sym01")
for k in range(3):
gaugerhs.alpha_rhs += diss_strength*alpha_dKOD[k]*rfm.ReU[k]
rhs.cf_rhs += diss_strength* cf_dKOD[k]*rfm.ReU[k]
rhs.trK_rhs += diss_strength* trK_dKOD[k]*rfm.ReU[k]
sfrhs.sf_rhs += diss_strength* sf_dKOD[k]*rfm.ReU[k]
sfrhs.sfM_rhs += diss_strength* sfM_dKOD[k]*rfm.ReU[k]
for i in range(3):
if "2ndOrder" in ShiftCondition:
gaugerhs.bet_rhsU[i] += diss_strength* betU_dKOD[i][k]*rfm.ReU[k]
gaugerhs.vet_rhsU[i] += diss_strength* vetU_dKOD[i][k]*rfm.ReU[k]
rhs.lambda_rhsU[i] += diss_strength*lambdaU_dKOD[i][k]*rfm.ReU[k]
for j in range(3):
rhs.a_rhsDD[i][j] += diss_strength*aDD_dKOD[i][j][k]*rfm.ReU[k]
rhs.h_rhsDD[i][j] += diss_strength*hDD_dKOD[i][j][k]*rfm.ReU[k]
# Now that we are finished with all the rfm hatted
# quantities in generic precomputed functional
# form, let's restore them to their closed-
# form expressions.
par.set_parval_from_str("reference_metric::enable_rfm_precompute","False") # Reset to False to disable rfm_precompute.
rfm.ref_metric__hatted_quantities()
end = time.time()
print("(BENCH) Finished BSSN symbolic expressions in "+str(end-start)+" seconds.")
def BSSN_plus_ScalarField_RHSs():
print("Generating C code for BSSN RHSs in "+par.parval_from_str("reference_metric::CoordSystem")+" coordinates.")
start = time.time()
# Construct the left-hand sides and right-hand-side expressions for all BSSN RHSs
lhs_names = [ "alpha", "cf", "trK", "sf", "sfM" ]
rhs_exprs = [gaugerhs.alpha_rhs, rhs.cf_rhs, rhs.trK_rhs, sfrhs.sf_rhs, sfrhs.sfM_rhs]
for i in range(3):
lhs_names.append( "betU"+str(i))
rhs_exprs.append(gaugerhs.bet_rhsU[i])
lhs_names.append( "lambdaU"+str(i))
rhs_exprs.append(rhs.lambda_rhsU[i])
lhs_names.append( "vetU"+str(i))
rhs_exprs.append(gaugerhs.vet_rhsU[i])
for j in range(i,3):
lhs_names.append( "aDD"+str(i)+str(j))
rhs_exprs.append(rhs.a_rhsDD[i][j])
lhs_names.append( "hDD"+str(i)+str(j))
rhs_exprs.append(rhs.h_rhsDD[i][j])
# Sort the lhss list alphabetically, and rhss to match.
# This ensures the RHSs are evaluated in the same order
# they're allocated in memory:
lhs_names,rhs_exprs = [list(x) for x in zip(*sorted(zip(lhs_names,rhs_exprs), key=lambda pair: pair[0]))]
# Declare the list of lhrh's
BSSN_evol_rhss = []
for var in range(len(lhs_names)):
BSSN_evol_rhss.append(lhrh(lhs=gri.gfaccess("rhs_gfs",lhs_names[var]),rhs=rhs_exprs[var]))
# Set up the C function for the BSSN RHSs
desc="Evaluate the BSSN RHSs"
name="rhs_eval"
outCfunction(
outfile = os.path.join(Ccodesdir,name+".h"), desc=desc, name=name,
params = """rfm_struct *restrict rfmstruct,const paramstruct *restrict params,
const REAL *restrict auxevol_gfs,const REAL *restrict in_gfs,REAL *restrict rhs_gfs""",
body = fin.FD_outputC("returnstring",BSSN_evol_rhss, params="outCverbose=False,enable_SIMD=True",
upwindcontrolvec=betaU),
loopopts = "InteriorPoints,enable_SIMD,enable_rfm_precompute,pragma_on_i0")
end = time.time()
print("(BENCH) Finished BSSN_RHS C codegen in " + str(end - start) + " seconds.")
def Ricci():
print("Generating C code for Ricci tensor in "+par.parval_from_str("reference_metric::CoordSystem")+" coordinates.")
start = time.time()
desc="Evaluate the Ricci tensor"
name="Ricci_eval"
outCfunction(
outfile = os.path.join(Ccodesdir,name+".h"), desc=desc, name=name,
params = """rfm_struct *restrict rfmstruct,const paramstruct *restrict params,
const REAL *restrict in_gfs,REAL *restrict auxevol_gfs""",
body = fin.FD_outputC("returnstring",
[lhrh(lhs=gri.gfaccess("auxevol_gfs","RbarDD00"),rhs=Bq.RbarDD[0][0]),
lhrh(lhs=gri.gfaccess("auxevol_gfs","RbarDD01"),rhs=Bq.RbarDD[0][1]),
lhrh(lhs=gri.gfaccess("auxevol_gfs","RbarDD02"),rhs=Bq.RbarDD[0][2]),
lhrh(lhs=gri.gfaccess("auxevol_gfs","RbarDD11"),rhs=Bq.RbarDD[1][1]),
lhrh(lhs=gri.gfaccess("auxevol_gfs","RbarDD12"),rhs=Bq.RbarDD[1][2]),
lhrh(lhs=gri.gfaccess("auxevol_gfs","RbarDD22"),rhs=Bq.RbarDD[2][2])],
params="outCverbose=False,enable_SIMD=True"),
loopopts = "InteriorPoints,enable_SIMD,enable_rfm_precompute,pragma_on_i0")
end = time.time()
print("(BENCH) Finished Ricci C codegen in " + str(end - start) + " seconds.")
Generating symbolic expressions for BSSN RHSs... (BENCH) Finished BSSN symbolic expressions in 11.03122878074646 seconds.
Step 4.b: Output the Hamiltonian constraint [Back to top]¶
$$\label{hamconstraint}$$Next output the C code for evaluating the Hamiltonian constraint (Tutorial). In the absence of numerical error, this constraint should evaluate to zero. However it does not due to numerical (typically truncation and roundoff) error. We will therefore measure the Hamiltonian constraint violation to gauge the accuracy of our simulation, and, ultimately determine whether errors are dominated by numerical finite differencing (truncation) error as expected.
In [8]:
def Hamiltonian():
start = time.time()
print("Generating optimized C code for Hamiltonian constraint. May take a while, depending on CoordSystem.")
# Set up the C function for the Hamiltonian RHS
desc="Evaluate the Hamiltonian constraint"
name="Hamiltonian_constraint"
outCfunction(
outfile = os.path.join(Ccodesdir,name+".h"), desc=desc, name=name,
params = """rfm_struct *restrict rfmstruct,const paramstruct *restrict params,
REAL *restrict in_gfs, REAL *restrict auxevol_gfs, REAL *restrict aux_gfs""",
body = fin.FD_outputC("returnstring",lhrh(lhs=gri.gfaccess("aux_gfs", "H"), rhs=bssncon.H),
params="outCverbose=False"),
loopopts = "InteriorPoints,enable_rfm_precompute,pragma_on_i0")
end = time.time()
print("(BENCH) Finished Hamiltonian C codegen in " + str(end - start) + " seconds.")
Step 4.c: Enforce conformal 3-metric $\det{\bar{\gamma}_{ij}}=\det{\hat{\gamma}_{ij}}$ constraint [Back to top]¶
$$\label{enforce3metric}$$Then enforce conformal 3-metric $\det{\bar{\gamma}_{ij}}=\det{\hat{\gamma}_{ij}}$ constraint (Eq. 53 of Ruchlin, Etienne, and Baumgarte (2018)), as documented in the corresponding NRPy+ tutorial notebook
Applying curvilinear boundary conditions should affect the initial data at the outer boundary, and will in general cause the $\det{\bar{\gamma}_{ij}}=\det{\hat{\gamma}_{ij}}$ constraint to be violated there. Thus after we apply these boundary conditions, we must always call the routine for enforcing the $\det{\bar{\gamma}_{ij}}=\det{\hat{\gamma}_{ij}}$ constraint:
In [9]:
def gammadet():
start = time.time()
print("Generating optimized C code for gamma constraint. May take a while, depending on CoordSystem.")
# Set up the C function for the det(gammahat) = det(gammabar)
EGC.output_Enforce_Detgammahat_Constraint_Ccode(Ccodesdir,exprs=enforce_detg_constraint_symb_expressions)
end = time.time()
print("(BENCH) Finished gamma constraint C codegen in " + str(end - start) + " seconds.")
In [10]:
# Step 4.d: C code kernel generation
# Step 4.d.i: Create a list of functions we wish to evaluate in parallel
funcs = [BSSN_plus_ScalarField_RHSs,Ricci,Hamiltonian,gammadet]
try:
if os.name == 'nt':
# It's a mess to get working in Windows, so we don't bother. :/
# https://medium.com/@grvsinghal/speed-up-your-python-code-using-multiprocessing-on-windows-and-jupyter-or-ipython-2714b49d6fac
raise Exception("Parallel codegen currently not available in Windows")
# Step 4.d.ii: Import the multiprocessing module.
import multiprocess as multiprocessing
# Step 4.d.iii: Define master function for parallelization.
# Note that lambdifying this doesn't work in Python 3
def master_func(arg):
funcs[arg]()
# Step 4.d.iv: Evaluate list of functions in parallel if possible;
# otherwise fallback to serial evaluation:
pool = multiprocessing.Pool()
pool.map(master_func,range(len(funcs)))
except:
# Steps 4.d.iii-4.d.v, alternate: As fallback, evaluate functions in serial.
for func in funcs:
func()
Generating optimized C code for Hamiltonian constraint. May take a while, depending on CoordSystem.Generating C code for Ricci tensor in SinhSpherical coordinates.Generating C code for BSSN RHSs in SinhSpherical coordinates.Generating optimized C code for gamma constraint. May take a while, depending on CoordSystem. Output C function enforce_detgammahat_constraint() to file BSSN_ScalarFieldCollapse_Ccodes/enforce_detgammahat_constraint.h (BENCH) Finished gamma constraint C codegen in 0.2685582637786865 seconds. Output C function Ricci_eval() to file BSSN_ScalarFieldCollapse_Ccodes/Ricci_eval.h (BENCH) Finished Ricci C codegen in 24.411128044128418 seconds. Output C function rhs_eval() to file BSSN_ScalarFieldCollapse_Ccodes/rhs_eval.h (BENCH) Finished BSSN_RHS C codegen in 25.715099811553955 seconds. Output C function Hamiltonian_constraint() to file BSSN_ScalarFieldCollapse_Ccodes/Hamiltonian_constraint.h (BENCH) Finished Hamiltonian C codegen in 43.34551215171814 seconds.
Step 4.e: Output C codes needed for declaring and setting Cparameters; also set free_parameters.h [Back to top]¶
$$\label{cparams_rfm_and_domainsize}$$Based on declared NRPy+ Cparameters, first we generate declare_Cparameters_struct.h, set_Cparameters_default.h, and set_Cparameters[-SIMD].h.
Then we output free_parameters.h, which sets initial data parameters, as well as grid domain & reference metric parameters, applying domain_size and sinh_width/SymTP_bScale (if applicable) as set above
In [11]:
# Step 4.e.i: Generate declare_Cparameters_struct.h, set_Cparameters_default.h, and set_Cparameters[-SIMD].h
par.generate_Cparameters_Ccodes(os.path.join(Ccodesdir))
# Step 4.e.ii: Set free_parameters.h
# Output to $Ccodesdir/free_parameters.h reference metric parameters based on generic
# domain_size,sinh_width,sinhv2_const_dr,SymTP_bScale,
# parameters set above.
rfm.out_default_free_parameters_for_rfm(os.path.join(Ccodesdir,"free_parameters.h"),
domain_size,sinh_width,sinhv2_const_dr,SymTP_bScale)
# Step 4.e.iii: Generate set_Nxx_dxx_invdx_params__and__xx.h:
rfm.set_Nxx_dxx_invdx_params__and__xx_h(Ccodesdir)
# Step 4.e.iv: Generate xx_to_Cart.h, which contains xx_to_Cart() for
# (the mapping from xx->Cartesian) for the chosen
# CoordSystem:
rfm.xx_to_Cart_h("xx_to_Cart","./set_Cparameters.h",os.path.join(Ccodesdir,"xx_to_Cart.h"))
# Step 4.e.v: Generate declare_Cparameters_struct.h, set_Cparameters_default.h, and set_Cparameters[-SIMD].h
par.generate_Cparameters_Ccodes(os.path.join(Ccodesdir))
Step 5: Set up boundary condition functions for chosen singular, curvilinear coordinate system [Back to top]¶
$$\label{bc_functs}$$Next apply singular, curvilinear coordinate boundary conditions as documented in the corresponding NRPy+ tutorial notebook.
In [12]:
import CurviBoundaryConditions.CurviBoundaryConditions as cbcs
cmd.mkdir(os.path.join(Ccodesdir,"boundary_conditions"))
cbcs.Set_up_CurviBoundaryConditions(os.path.join(Ccodesdir,"boundary_conditions/"),Cparamspath=os.path.join("../"),enable_copy_of_static_Ccodes=False)
# Manually copy the static files required by boundary conditions
for file in ["apply_bcs_curvilinear.h","BCs_data_structs.h","bcstruct_freemem.h","CurviBC_include_Cfunctions.h",
"driver_bcstruct.h","set_bcstruct.h","set_up__bc_gz_map_and_parity_condns.h"]:
shutil.copy(os.path.join("..","CurviBoundaryConditions","boundary_conditions",file),
os.path.join(Ccodesdir,"boundary_conditions"))
with open(os.path.join(Ccodesdir,"boundary_conditions","CurviBC_include_Cfunctions.h"),"a") as file:
file.write("\n#include \"apply_bcs_curvilinear.h\"")
Wrote to file "BSSN_ScalarFieldCollapse_Ccodes/boundary_conditions/parity_conditions_symbolic_dot_products.h"
Evolved parity: ( aDD00:4, aDD01:5, aDD02:6, aDD11:7, aDD12:8, aDD22:9,
alpha:0, betU0:1, betU1:2, betU2:3, cf:0, hDD00:4, hDD01:5, hDD02:6,
hDD11:7, hDD12:8, hDD22:9, lambdaU0:1, lambdaU1:2, lambdaU2:3, sf:0,
sfM:0, trK:0, vetU0:1, vetU1:2, vetU2:3 )
Auxiliary parity: ( H:0 )
AuxEvol parity: ( RbarDD00:4, RbarDD01:5, RbarDD02:6, RbarDD11:7,
RbarDD12:8, RbarDD22:9 )
Wrote to file "BSSN_ScalarFieldCollapse_Ccodes/boundary_conditions/EigenCoord_Cart_to_xx.h"
In [13]:
# Part P0: Define REAL, set the number of ghost cells NGHOSTS (from NRPy+'s FD_CENTDERIVS_ORDER),
# and set the CFL_FACTOR (which can be overwritten at the command line)
with open(os.path.join(Ccodesdir,"ScalarFieldCollapse_Playground_REAL__NGHOSTS__CFL_FACTOR.h"), "w") as file:
file.write("""
// Part P0.a: Set the number of ghost cells, from NRPy+'s FD_CENTDERIVS_ORDER
#define NGHOSTS """+str(int(FD_order/2)+1)+"""
// Part P0.b: Set the numerical precision (REAL) to double, ensuring all floating point
// numbers are stored to at least ~16 significant digits
#define REAL """+REAL+"""
// Part P0.c: Set the number of ghost cells, from NRPy+'s FD_CENTDERIVS_ORDER
REAL CFL_FACTOR = """+str(CFL_FACTOR)+"""; // Set the CFL Factor. Can be overwritten at command line.\n""")
files = ["NRPyCritCol_regridding.h","ScalarField_output_central_values.h"]
for file in files:
shutil.copyfile(os.path.join(file),os.path.join(Ccodesdir,file))
outfile = os.path.join(Ccodesdir,"rfm_files","rfm_struct__define-pointer.h")
shutil.copyfile(os.path.join(Ccodesdir,"rfm_files","rfm_struct__define.h"),
outfile)
with open(outfile,"r") as file:
file_contents = file.read()
with open(outfile,"w") as file:
file.write(file_contents.replace("rfmstruct.","rfmstruct->"))
In [14]:
%%writefile $Ccodesdir/ScalarFieldCollapse_Playground.c
// Step P0: Define REAL and NGHOSTS; and declare CFL_FACTOR. This header is generated in NRPy+.
#include "ScalarFieldCollapse_Playground_REAL__NGHOSTS__CFL_FACTOR.h"
#include "rfm_files/rfm_struct__declare.h"
#include "declare_Cparameters_struct.h"
// All SIMD intrinsics used in SIMD-enabled C code loops are defined here:
#include "SIMD/SIMD_intrinsics.h"
// Step P1: Import needed header files
#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include "time.h"
#include "stdint.h" // Needed for Windows GCC 6.x compatibility
#ifndef M_PI
#define M_PI 3.141592653589793238462643383279502884L
#endif
#ifndef M_SQRT1_2
#define M_SQRT1_2 0.707106781186547524400844362104849039L
#endif
#define wavespeed 1.0 // Set CFL-based "wavespeed" to 1.0.
#define alpha_threshold (2e-3) // Value below which we rule gravitational collapse has happened
// Step P2: Declare the IDX4S(gf,i,j,k) macro, which enables us to store 4-dimensions of
// data in a 1D array. In this case, consecutive values of "i"
// (all other indices held to a fixed value) are consecutive in memory, where
// consecutive values of "j" (fixing all other indices) are separated by
// Nxx_plus_2NGHOSTS0 elements in memory. Similarly, consecutive values of
// "k" are separated by Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1 in memory, etc.
#define IDX4S(g,i,j,k) \
( (i) + Nxx_plus_2NGHOSTS0 * ( (j) + Nxx_plus_2NGHOSTS1 * ( (k) + Nxx_plus_2NGHOSTS2 * (g) ) ) )
#define IDX4ptS(g,idx) ( (idx) + (Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1*Nxx_plus_2NGHOSTS2) * (g) )
#define IDX3S(i,j,k) ( (i) + Nxx_plus_2NGHOSTS0 * ( (j) + Nxx_plus_2NGHOSTS1 * ( (k) ) ) )
#define LOOP_REGION(i0min,i0max, i1min,i1max, i2min,i2max) \
for(int i2=i2min;i2<i2max;i2++) for(int i1=i1min;i1<i1max;i1++) for(int i0=i0min;i0<i0max;i0++)
#define LOOP_ALL_GFS_GPS(ii) _Pragma("omp parallel for") \
for(int (ii)=0;(ii)<Nxx_plus_2NGHOSTS_tot*NUM_EVOL_GFS;(ii)++)
// Step P3: Set UUGF and VVGF macros, as well as xx_to_Cart()
#include "boundary_conditions/gridfunction_defines.h"
// Step P4: Set xx_to_Cart(const paramstruct *restrict params,
// REAL *restrict xx[3],
// const int i0,const int i1,const int i2,
// REAL xCart[3]),
// which maps xx->Cartesian via
// {xx[0][i0],xx[1][i1],xx[2][i2]}->{xCart[0],xCart[1],xCart[2]}
#include "xx_to_Cart.h"
// Step P5: Defines set_Nxx_dxx_invdx_params__and__xx(const int EigenCoord, const int Nxx[3],
// paramstruct *restrict params, REAL *restrict xx[3]),
// which sets params Nxx,Nxx_plus_2NGHOSTS,dxx,invdx, and xx[] for
// the chosen Eigen-CoordSystem if EigenCoord==1, or
// CoordSystem if EigenCoord==0.
#include "set_Nxx_dxx_invdx_params__and__xx.h"
// Step P6: Include basic functions needed to impose curvilinear
// parity and boundary conditions.
#include "boundary_conditions/CurviBC_include_Cfunctions.h"
// Step P7: Implement the algorithm for upwinding.
// *NOTE*: This upwinding is backwards from
// usual upwinding algorithms, because the
// upwinding control vector in BSSN (the shift)
// acts like a *negative* velocity.
//#define UPWIND_ALG(UpwindVecU) UpwindVecU > 0.0 ? 1.0 : 0.0
// Step P8: Include function for enforcing detgammabar constraint.
#include "enforce_detgammahat_constraint.h"
// Step P9: Find the CFL-constrained timestep
#include "find_timestep.h"
// Step P10: Declare initial data input struct:
// stores data from initial data solver,
// so they can be put on the numerical grid.
typedef struct __ID_inputs {
int interp_stencil_size;
int numlines_in_file;
REAL *r_arr,*sf_arr,*psi4_arr,*alpha_arr;
} ID_inputs;
// Part P11: Declare all functions for setting up ScalarField initial data.
/* Routines to interpolate the ScalarField solution and convert to ADM & T^{munu}: */
#include "../ScalarField_interp.h"
#include "ID_scalarfield_ADM_quantities.h"
#include "ID_scalarfield_spherical.h"
#include "ID_scalarfield_xx0xx1xx2_to_BSSN_xx0xx1xx2.h"
#include "ID_scalarfield.h"
/* Next perform the basis conversion and compute all needed BSSN quantities */
#include "ID_ADM_xx0xx1xx2_to_BSSN_xx0xx1xx2__ALL_BUT_LAMBDAs.h"
#include "ID_BSSN__ALL_BUT_LAMBDAs.h"
#include "ID_BSSN_lambdas.h"
// Step P12: Set the generic driver function for setting up BSSN initial data
void initial_data(const paramstruct *restrict params,const bc_struct *restrict bcstruct,
const rfm_struct *restrict rfmstruct,
REAL *restrict xx[3], REAL *restrict auxevol_gfs, REAL *restrict in_gfs) {
#include "set_Cparameters.h"
// Step 1: Set up ScalarField initial data
// Step 1.a: Read ScalarField initial data from data file
// Open the data file:
char filename[100];
sprintf(filename,"./SFID.txt");
FILE *fp = fopen(filename, "r");
if (fp == NULL) {
fprintf(stderr,"ERROR: could not open file %s\n",filename);
exit(1);
}
// Count the number of lines in the data file:
int numlines_in_file = count_num_lines_in_file(fp);
// Allocate space for all data arrays:
REAL *r_arr = (REAL *)malloc(sizeof(REAL)*numlines_in_file);
REAL *sf_arr = (REAL *)malloc(sizeof(REAL)*numlines_in_file);
REAL *psi4_arr = (REAL *)malloc(sizeof(REAL)*numlines_in_file);
REAL *alpha_arr = (REAL *)malloc(sizeof(REAL)*numlines_in_file);
// Read from the data file, filling in arrays
// read_datafile__set_arrays() may be found in ScalarField/ScalarField_interp.h
if(read_datafile__set_arrays(fp,r_arr,sf_arr,psi4_arr,alpha_arr) == 1) {
fprintf(stderr,"ERROR WHEN READING FILE %s!\n",filename);
exit(1);
}
fclose(fp);
const int interp_stencil_size = 12;
ID_inputs SF_in;
SF_in.interp_stencil_size = interp_stencil_size;
SF_in.numlines_in_file = numlines_in_file;
SF_in.r_arr = r_arr;
SF_in.sf_arr = sf_arr;
SF_in.psi4_arr = psi4_arr;
SF_in.alpha_arr = alpha_arr;
// Step 1.b: Interpolate data from data file to set BSSN gridfunctions
ID_scalarfield(params,xx,SF_in, in_gfs);
ID_BSSN__ALL_BUT_LAMBDAs(params,xx,SF_in, in_gfs);
apply_bcs_curvilinear(params, bcstruct, NUM_EVOL_GFS, evol_gf_parity, in_gfs);
enforce_detgammahat_constraint(rfmstruct, params, in_gfs);
ID_BSSN_lambdas(params, xx, in_gfs);
apply_bcs_curvilinear(params, bcstruct, NUM_EVOL_GFS, evol_gf_parity, in_gfs);
enforce_detgammahat_constraint(rfmstruct, params, in_gfs);
free(r_arr);
free(sf_arr);
free(psi4_arr);
free(alpha_arr);
}
// Step P11: Declare function for evaluating Hamiltonian constraint (diagnostic)
#include "Hamiltonian_constraint.h"
// Step P12: Declare rhs_eval function, which evaluates BSSN RHSs
#include "rhs_eval.h"
// Step P13: Declare Ricci_eval function, which evaluates Ricci tensor
#include "Ricci_eval.h"
#define max_number_of_regrids (16)
#include "ScalarField_output_central_values.h"
#include "NRPyCritCol_regridding.h"
REAL rho_max = 0.0;
// main() function:
// Step 0: Read command-line input, set up grid structure, allocate memory for gridfunctions, set up coordinates
// Step 1: Set up initial data to an exact solution
// Step 2: Start the timer, for keeping track of how fast the simulation is progressing.
// Step 3: Integrate the initial data forward in time using the chosen RK-like Method of
// Lines timestepping algorithm, and output periodic simulation diagnostics
// Step 3.a: Output 2D data file periodically, for visualization
// Step 3.b: Step forward one timestep (t -> t+dt) in time using
// chosen RK-like MoL timestepping algorithm
// Step 3.c: If t=t_final, output conformal factor & Hamiltonian
// constraint violation to 1D data file
// Step 3.d: Progress indicator printing to stderr
// Step 4: Free all allocated memory
int main(int argc, const char *argv[]) {
paramstruct params;
#include "set_Cparameters_default.h"
// Step 0a: Read command-line input, error out if nonconformant
if((argc != 4 && argc != 5) || atoi(argv[1]) < NGHOSTS || atoi(argv[2]) < 2 || atoi(argv[3]) < 2 /* FIXME; allow for axisymmetric sims */) {
fprintf(stderr,"Error: Expected three command-line arguments: ./ScalarFieldCollapse_Playground Nx0 Nx1 Nx2,\n");
fprintf(stderr,"where Nx[0,1,2] is the number of grid points in the 0, 1, and 2 directions.\n");
fprintf(stderr,"Nx[] MUST BE larger than NGHOSTS (= %d)\n",NGHOSTS);
exit(1);
}
if(argc == 5) {
CFL_FACTOR = strtod(argv[4],NULL);
if(CFL_FACTOR > 0.5 && atoi(argv[3])!=2) {
fprintf(stderr,"WARNING: CFL_FACTOR was set to %e, which is > 0.5.\n",CFL_FACTOR);
fprintf(stderr," This will generally only be stable if the simulation is purely axisymmetric\n");
fprintf(stderr," However, Nx2 was set to %d>2, which implies a non-axisymmetric simulation\n",atoi(argv[3]));
}
}
// Step 0b: Set up numerical grid structure, first in space...
const int Nxx[3] = { atoi(argv[1]), atoi(argv[2]), atoi(argv[3]) };
if(Nxx[0]%2 != 0 || Nxx[1]%2 != 0 || Nxx[2]%2 != 0) {
fprintf(stderr,"Error: Cannot guarantee a proper cell-centered grid if number of grid cells not set to even number.\n");
fprintf(stderr," For example, in case of angular directions, proper symmetry zones will not exist.\n");
exit(1);
}
// Step 0c: Set free parameters, overwriting Cparameters defaults
// by hand or with command-line input, as desired.
#include "free_parameters.h"
// Start with no KO dissipation
params.diss_strength = 0.0;
params.eta = 0.5/0.218; // 0.5/M_ADM
// Set up regridding struct
NRPyCritCol_regrid_params_struct regrid_params;
#include "set_regrid_struct.h"
// Step 0d: Uniform coordinate grids are stored to *xx[3]
REAL *xx[3];
// Step 0d.i: Set bcstruct
bc_struct bcstruct;
{
int EigenCoord = 1;
// Step 0d.ii: Call set_Nxx_dxx_invdx_params__and__xx(), which sets
// params Nxx,Nxx_plus_2NGHOSTS,dxx,invdx, and xx[] for the
// chosen Eigen-CoordSystem.
set_Nxx_dxx_invdx_params__and__xx(EigenCoord, Nxx, ¶ms, xx);
// Step 0d.iii: Set Nxx_plus_2NGHOSTS_tot
#include "set_Cparameters-nopointer.h"
const int Nxx_plus_2NGHOSTS_tot = Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1*Nxx_plus_2NGHOSTS2;
// Step 0e: Find ghostzone mappings; set up bcstruct
#include "boundary_conditions/driver_bcstruct.h"
// Step 0e.i: Free allocated space for xx[][] array
for(int i=0;i<3;i++) free(xx[i]);
}
// Step 0f: Call set_Nxx_dxx_invdx_params__and__xx(), which sets
// params Nxx,Nxx_plus_2NGHOSTS,dxx,invdx, and xx[] for the
// chosen (non-Eigen) CoordSystem.
int EigenCoord = 0;
set_Nxx_dxx_invdx_params__and__xx(EigenCoord, Nxx, ¶ms, xx);
// Step 0g: Set all C parameters "blah" for params.blah, including
// Nxx_plus_2NGHOSTS0 = params.Nxx_plus_2NGHOSTS0, etc.
#include "set_Cparameters-nopointer.h"
const int Nxx_plus_2NGHOSTS_tot = Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1*Nxx_plus_2NGHOSTS2;
// Step 0h: Time coordinate parameters
REAL t_final = 6.6; /* Final time is set so that at t=t_final,
* data at the origin have not been corrupted
* by the approximate outer boundary condition */
// Step 0i: Set timestep based on smallest proper distance between gridpoints and CFL factor
REAL dt = find_timestep(¶ms, xx);
//fprintf(stderr,"# Timestep set to = %e\n",(double)dt);
int N_final = (int)(t_final / dt + 0.5); // The number of points in time.
// Add 0.5 to account for C rounding down
// typecasts to integers.
int output_every_N = 20;//(int)((REAL)N_final/800.0);
if(output_every_N == 0) output_every_N = 1;
// Step 0j: Error out if the number of auxiliary gridfunctions outnumber evolved gridfunctions.
// This is a limitation of the RK method. You are always welcome to declare & allocate
// additional gridfunctions by hand.
if(NUM_AUX_GFS > NUM_EVOL_GFS) {
fprintf(stderr,"Error: NUM_AUX_GFS > NUM_EVOL_GFS. Either reduce the number of auxiliary gridfunctions,\n");
fprintf(stderr," or allocate (malloc) by hand storage for *diagnostic_output_gfs. \n");
exit(1);
}
// Step 0k: Allocate memory for gridfunctions
#include "MoLtimestepping/RK_Allocate_Memory.h"
REAL *restrict auxevol_gfs = (REAL *)malloc(sizeof(REAL) * NUM_AUXEVOL_GFS * Nxx_plus_2NGHOSTS_tot);
// Step 0l: Set up precomputed reference metric arrays
// Step 0l.i: Allocate space for precomputed reference metric arrays.
#include "rfm_files/rfm_struct__malloc.h"
// Step 0l.ii: Define precomputed reference metric arrays.
{
#include "set_Cparameters-nopointer.h"
#include "rfm_files/rfm_struct__define.h"
}
// Step 1: Set up initial data to an exact solution
initial_data(¶ms,&bcstruct, &rfmstruct, xx, auxevol_gfs, y_n_gfs);
// Step 1b: Apply boundary conditions, as initial data
// are sometimes ill-defined in ghost zones.
// E.g., spherical initial data might not be
// properly defined at points where r=-1.
apply_bcs_curvilinear(¶ms, &bcstruct, NUM_EVOL_GFS,evol_gf_parity, y_n_gfs);
enforce_detgammahat_constraint(&rfmstruct, ¶ms, y_n_gfs);
// Step 2: Start the timer, for keeping track of how fast the simulation is progressing.
#ifdef __linux__ // Use high-precision timer in Linux.
struct timespec start, end;
clock_gettime(CLOCK_REALTIME, &start);
#else // Resort to low-resolution, standards-compliant timer in non-Linux OSs
// http://www.cplusplus.com/reference/ctime/time/
time_t start_timer,end_timer;
time(&start_timer); // Resolution of one second...
#endif
int number_of_regrids_performed = 0;
REAL t = 0.0;
// Step 3: Integrate the initial data forward in time using the chosen RK-like Method of
// Lines timestepping algorithm, and output periodic simulation diagnostics
for(int n=0;n<=N_final;n++) { // Main loop to progress forward in time.
// Step 3.a: Step forward one timestep (t -> t+dt) in time using
// chosen RK-like MoL timestepping algorithm
#include "MoLtimestepping/RK_MoL.h"
t += dt;
// Step 3.b: Output central values
int lapse_collapsed = output_central_values(t,¶ms,y_n_gfs);
if( lapse_collapsed ) break;
// Step 3.c: Check if it's time to regrid
if( (t > regrid_params.regrid_time[number_of_regrids_performed]) &&
(number_of_regrids_performed == regrid_params.regrid_counter[number_of_regrids_performed]) &&
(number_of_regrids_performed < max_number_of_regrids) ) {
regrid(regrid_params.regrid_key[number_of_regrids_performed],
n,t,
regrid_params.new_ampl_or_sinhW[number_of_regrids_performed],
&bcstruct,&rfmstruct,¶ms,
&N_final,&t_final,&dt,
xx,diagnostic_output_gfs,y_n_gfs);
number_of_regrids_performed++;
// Turn on KO dissipation
params.diss_strength = 1.0;
}
// Step 3.d: Progress indicator printing to stderr
// Step 3.d.i: Measure average time per iteration
#ifdef __linux__ // Use high-precision timer in Linux.
clock_gettime(CLOCK_REALTIME, &end);
const long long unsigned int time_in_ns = 1000000000L * (end.tv_sec - start.tv_sec) + end.tv_nsec - start.tv_nsec;
#else // Resort to low-resolution, standards-compliant timer in non-Linux OSs
time(&end_timer); // Resolution of one second...
REAL time_in_ns = difftime(end_timer,start_timer)*1.0e9+0.5; // Round up to avoid divide-by-zero.
#endif
const REAL s_per_iteration_avg = ((REAL)time_in_ns / (REAL)n) / 1.0e9;
const int iterations_remaining = N_final - n;
const REAL time_remaining_in_mins = s_per_iteration_avg * (REAL)iterations_remaining / 60.0;
const REAL num_RHS_pt_evals = (REAL)(Nxx[0]*Nxx[1]*Nxx[2]) * 4.0 * (REAL)n; // 4 RHS evals per gridpoint for RK4
const REAL RHS_pt_evals_per_sec = num_RHS_pt_evals / ((REAL)time_in_ns / 1.0e9);
// Step 3.d.ii: Output simulation progress to stderr
if(n%10 == 0) {
fprintf(stderr,"%c[2K", 27); // Clear the line
fprintf(stderr,"It: %d t=%.2f dt=%.2e | %.1f%%; ETA %.0f s | t/h %.2f | gp/s %.2e\r", // \r is carriage return, move cursor to the beginning of the line
n, t, (double)dt, (double)(100.0 * (REAL)n / (REAL)N_final),
(double)time_remaining_in_mins*60, (double)(dt * 3600.0 / s_per_iteration_avg), (double)RHS_pt_evals_per_sec);
fflush(stderr); // Flush the stderr buffer
} // End progress indicator if(n % 10 == 0)
} // End main loop to progress forward in time.
fprintf(stderr,"\n"); // Clear the final line of output from progress indicator.
// Step 4: Free all allocated memory
#include "rfm_files/rfm_struct__freemem.h"
#include "boundary_conditions/bcstruct_freemem.h"
#include "MoLtimestepping/RK_Free_Memory.h"
free(auxevol_gfs);
for(int i=0;i<3;i++) free(xx[i]);
return 0;
}
Writing BSSN_ScalarFieldCollapse_Ccodes/ScalarFieldCollapse_Playground.c
In [15]:
import os
import time
import cmdline_helper as cmd
print("Now compiling, should take ~10 seconds...\n")
start = time.time()
cmd.C_compile(os.path.join(Ccodesdir,"ScalarFieldCollapse_Playground.c"),
os.path.join(outdir,"ScalarFieldCollapse_Playground"),compile_mode="optimized")
end = time.time()
print("(BENCH) Finished in "+str(end-start)+" seconds.\n")
# Change to output directory
os.chdir(outdir)
# Clean up existing output files
cmd.delete_existing_files("out*.txt")
cmd.delete_existing_files("out*.png")
# Run executable
print(os.getcwd())
start = time.time()
shutil.copyfile("SFID_weak.txt","SFID.txt")
cmd.Execute("ScalarFieldCollapse_Playground", "320 2 2 "+str(CFL_FACTOR),"out320.txt")
shutil.copyfile("out_central_values.dat","out_weak.dat")
os.remove("out_central_values.dat")
shutil.copyfile("SFID_strong.txt","SFID.txt")
cmd.Execute("ScalarFieldCollapse_Playground", "320 2 2 "+str(CFL_FACTOR),"out320.txt")
shutil.copyfile("out_central_values.dat","out_strong.dat")
os.remove("out_central_values.dat")
end = time.time()
print("(BENCH) Finished in "+str(end-start)+" seconds.\n")
# Return to root directory
os.chdir(os.path.join("../../"))
Now compiling, should take ~10 seconds... Compiling executable... (EXEC): Executing `gcc -std=gnu99 -Ofast -fopenmp -march=native -funroll-loops BSSN_ScalarFieldCollapse_Ccodes/ScalarFieldCollapse_Playground.c -o BSSN_ScalarFieldCollapse_Ccodes/output/ScalarFieldCollapse_Playground -lm`... (BENCH): Finished executing in 6.10278582572937 seconds. Finished compilation. (BENCH) Finished in 6.122713088989258 seconds. /Users/werneck/Codes/nrpytutorial/ScalarField/BSSN_ScalarFieldCollapse_Ccodes/output (EXEC): Executing `./ScalarFieldCollapse_Playground 320 2 2 0.5`... It: 24300 t=6.60 dt=1.78e-05 | 100.0%; ETA 0 s | t/h 5.98 | gp/s 4.77e+05557e+14 (BENCH): Finished executing in 261.22367572784424 seconds. (EXEC): Executing `./ScalarFieldCollapse_Playground 320 2 2 0.5`... It: 23280 t=6.58 dt=1.78e-05 | 95.8%; ETA 11 s | t/h 6.05 | gp/s 4.83e+05514e+14 (BENCH): Finished executing in 247.43084621429443 seconds. (BENCH) Finished in 508.7111668586731 seconds.
In [21]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import ConnectionPatch
from matplotlib.patches import Rectangle
from IPython.display import Image
def draw_rectangle_and_connecting_lines(ax1_in,ax2_in,rec_left,rec_right,rec_bot,rec_top,connect):
# Draw rectangle
con = Rectangle((rec_left,rec_bot),rec_right-rec_left,rec_top-rec_bot,edgecolor='black',facecolor='none',lw=0.5,ls='--')
ax1_in.add_artist(con)
# Draw connecting lines from plot 1 to plot 2
if connect == 'bottom_to_top':
xys_1 = [(rec_left,rec_top),(rec_right,rec_top)]
xys_2 = [(rec_left,rec_bot),(rec_right,rec_bot)]
for i in range(2):
con = ConnectionPatch(xyA=xys_1[i], xyB=xys_2[i], coordsA="data", coordsB="data",
axesA=ax2_in, axesB=ax1_in,color="black",ls='--',lw=0.5)
ax1_in.add_artist(con)
elif connect == 'left_to_right':
xys_1 = [(rec_left,rec_top),(rec_left,rec_bot)]
xys_2 = [(rec_right,rec_top),(rec_right,rec_bot)]
for i in range(2):
con = ConnectionPatch(xyA=xys_1[i], xyB=xys_2[i], coordsA="data", coordsB="data",
axesA=ax1_in, axesB=ax2_in,color="black",ls='--',lw=0.5)
ax1_in.add_artist(con)
# Output file name
outfile = "lapse_self_similarity.png"
# Load weak data
t_w,alp_w,sf_w = np.loadtxt(os.path.join(Ccodesdir,"output","out_weak.dat")).T
# Load strong data
t_s,alp_s,sf_s = np.loadtxt(os.path.join(Ccodesdir,"output","out_strong.dat")).T
fig = plt.figure()
linewidth = 1
color_w = 'blue'
color_s = 'orange'
# Top panel
ax1 = plt.subplot(2, 1, 1)
plt.grid(ls=':')
plt.plot(t_w,alp_w,lw=linewidth,c=color_w,label=r"$\eta_{\rm weak} = 0.303326061$")
plt.plot(t_s,alp_s,lw=linewidth,c=color_s,ls='--',label=r"$\eta_{\rm strong} = 0.303326062$")
plt.legend(loc=2,markerfirst=False)
plt.xlim(-0.4,7.4)
plt.ylim(-0.2,1.2)
plt.xticks([0,1,2,3,4,5,6,7],['0','1','2','3','4','5','6','7'])
plt.yticks([0,0.2,0.4,0.6,0.8,1],['0.0','0.2','0.4','0.6','0.8','1.0'])
# Bottom right panel
xl2 = [+5.6,+6.8]
yl2 = [-0.1,+0.9]
ax2 = plt.subplot(2, 2, 4)
plt.grid(ls=':')
plt.plot(t_w,alp_w,lw=linewidth,c=color_w)
plt.plot(t_s,alp_s,lw=linewidth,c=color_s,ls='--')
plt.xlim(xl2[0],xl2[1])
plt.ylim(yl2[0],yl2[1])
plt.xticks([5.8,6.2,6.6],['5.9','6.2','6.5'])
plt.yticks([0,0.4,0.8],['0.0','0.4','0.8'])
# Bottom left panel
xl3 = [+6.51,+6.61]
yl3 = [-0.05,+0.55]
ax3 = plt.subplot(2, 2, 3)
plt.grid(ls=':')
plt.plot(t_w,alp_w,lw=linewidth,c=color_w)
plt.plot(t_s,alp_s,lw=linewidth,c=color_s,ls='--')
plt.xlim(xl3[0],xl3[1])
plt.ylim(yl3[0],yl3[1])
plt.xticks([6.52,6.56,6.6],['6.54','6.56','6.58'])
plt.yticks([0,0.25,0.5],['0.00','0.25','0.50'])
draw_rectangle_and_connecting_lines(ax1,ax2,xl2[0],xl2[1],yl2[0],yl2[1],'bottom_to_top')
draw_rectangle_and_connecting_lines(ax2,ax3,xl3[0],xl3[1],yl3[0],yl3[1],'left_to_right')
# Set labels
fig.text(0.5, 0.03, r"$t$", ha='center', va='center')
fig.text(0.03, 0.5, r"$\alpha_{\rm central}$", ha='center', va='center', rotation='vertical')
plt.savefig(outfile,dpi=300,bbox_inches='tight',facecolor='white')
plt.close(fig)
Image(outfile)
Out[21]:
Step 8: Output this module as $\LaTeX$-formatted PDF file [Back to top]¶
$$\label{output_to_pdf}$$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-Start_to_Finish-BSSNCurvilinear-ScalarField_Collapse.pdf (Note that clicking on this link may not work; you may need to open the PDF file through another means.)
In [17]:
import cmdline_helper as cmd # NRPy+: Multi-platform Python command-line interface
shutil.copy(os.path.join("..","latex_nrpy_style.tplx"),"latex_nrpy_style.tplx")
cmd.output_Jupyter_notebook_to_LaTeXed_PDF("Tutorial-Start_to_Finish-BSSNCurvilinear-ScalarField_Collapse_with_regrids")
os.remove("latex_nrpy_style.tplx")
Created Tutorial-Start_to_Finish-BSSNCurvilinear-
ScalarField_Collapse_with_regrids.tex, and compiled LaTeX file to PDF
file Tutorial-Start_to_Finish-BSSNCurvilinear-
ScalarField_Collapse_with_regrids.pdf