{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"# Start-to-Finish Example: Head-On Black Hole Collision\n",
"\n",
"## Author: Zach Etienne\n",
"### Formatting improvements courtesy Brandon Clark\n",
"\n",
"## This module implements a basic numerical relativity code to merge two black holes in *spherical coordinates*. The method consists of generating the necessary C source code for initial data, progressing it with a Runge-Kutta fourth-order scheme, and visualizing the black hole evolution by plotting the conformal factor. The validation confirms numerical error convergence to zero at higher resolutions.\n",
"\n",
"### Here we place the black holes initially on the $z$-axis, so the entire simulation is axisymmetric about the $\\phi$-axis. Not sampling in the $\\phi$ direction greatly speeds up the simulation.\n",
"\n",
"**Notebook Status:** Validated \n",
"\n",
"**Validation Notes:** This module has been validated to exhibit convergence to zero of the Hamiltonian constraint violation at the expected order to the exact solution *after a short numerical evolution of the initial data* (see [plots at bottom](#convergence)), and all quantities have been validated against the [original SENR code](https://bitbucket.org/zach_etienne/nrpy).\n",
"\n",
"### NRPy+ Source Code & Tutorials for this module: \n",
"* [BSSN/BrillLindquist.py](../edit/BSSN/BrillLindquist.py); [\\[**tutorial**\\]](Tutorial-ADM_Initial_Data-Brill-Lindquist.ipynb): Brill-Lindquist initial data; sets all ADM variables in Cartesian basis: \n",
"* [BSSN/ADM_Initial_Data_Reader__BSSN_Converter.py](../edit/BSSN/ADM_Initial_Data_Reader__BSSN_Converter.py); [\\[**tutorial**\\]](Tutorial-ADM_Initial_Data_Reader__BSSN_Converter.ipynb): Registers the C function for our \"universal\" initial data reader/converter initial_data_reader__convert_ADM_Cartesian_to_BSSN().\n",
"* [MoLtimestepping/MoL.py](../edit/MoLtimestepping/MoL.py); ([**NRPy+ Tutorial documentation**](Tutorial-Method_of_Lines-C_Code_Generation.ipynb)): Registers C functions for Method of Lines time integration, to push our initial data forward in time.\n",
"* [CurviBoundaryConditions/CurviBoundaryConditions.py](../edit/CurviBoundaryConditions/CurviBoundaryConditions.py); [\\[**tutorial**\\]](Tutorial-Start_to_Finish-Curvilinear_BCs.ipynb): Registers the C function for our \"universal\" initial data reader/converter initial_data_reader__convert_ADM_Cartesian_to_BSSN(). This applies boundary conditions to BSSN quantities (including $\\lambda^i$, which is computed via finite difference derivatives and thus only defined in grid interior)\n",
"* [BSSN/BSSN_Ccodegen_library.py](../edit/BSSN/BSSN_Ccodegen_library.py); [\\[**tutorial**\\]](Tutorial-BSSN_time_evolution-C_codegen_library.ipynb): Implements a number of helper functions for generating C codes from symbolic expressions generated in the following modules/tutorials:\n",
" * [BSSN/BSSN_constraints.py](../edit/BSSN/BSSN_constraints.py); [\\[**tutorial**\\]](Tutorial-BSSN_constraints.ipynb): Hamiltonian constraint in BSSN curvilinear basis/coordinates\n",
" * [BSSN/BSSN_RHSs.py](../edit/BSSN/BSSN_RHSs.py); [\\[**tutorial**\\]](Tutorial-BSSN_time_evolution-BSSN_RHSs.ipynb): Generates the right-hand sides for the BSSN evolution equations in singular, curvilinear coordinates\n",
" * [BSSN/BSSN_gauge_RHSs.py](../edit/BSSN/BSSN_gauge_RHSs.py); [\\[**tutorial**\\]](Tutorial-BSSN_time_evolution-BSSN_gauge_RHSs.ipynb): Generates the right-hand sides for the BSSN gauge evolution equations in singular, curvilinear coordinates\n",
"\n",
"## Introduction:\n",
"Here we use NRPy+ to generate the C source code necessary to set up initial data for two black holes (Brill-Lindquist, [Brill & Lindquist, Phys. Rev. 131, 471, 1963](https://journals.aps.org/pr/abstract/10.1103/PhysRev.131.471); see also Eq. 1 of [Brandt & Brügmann, arXiv:gr-qc/9711015v1](https://arxiv.org/pdf/gr-qc/9711015v1.pdf)). Then we use it to generate the RHS expressions for [Method of Lines](https://reference.wolfram.com/language/tutorial/NDSolveMethodOfLines.html) time integration based on an [explicit Runge-Kutta fourth-order scheme](https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods) (RK4 is chosen below, but multiple options exist). \n",
"\n",
"The entire algorithm is outlined as follows, with links to the relevant NRPy+ tutorial notebooks listed at each step:\n",
"\n",
"1. Allocate memory for gridfunctions, including temporary storage for the Method of Lines time integration\n",
" * [**NRPy+ tutorial on Method of Lines algorithm**](Tutorial-Method_of_Lines-C_Code_Generation.ipynb).\n",
"1. Set gridfunction values to initial data \n",
" * [**NRPy+ tutorial on Brill-Lindquist initial data**](Tutorial-ADM_Initial_Data-Brill-Lindquist.ipynb)\n",
" * [**NRPy+ tutorial on validating Brill-Lindquist initial data**](Tutorial-Start_to_Finish-BSSNCurvilinear-Exact_Initial_Data.ipynb).\n",
"1. Next, integrate the initial data forward in time using the Method of Lines coupled to a Runge-Kutta explicit timestepping algorithm:\n",
" 1. At the start of each iteration in time, output the Hamiltonian constraint violation \n",
" * [**NRPy+ tutorial on BSSN constraints**](Tutorial-BSSN_constraints.ipynb).\n",
" 1. At each RK time substep, do the following:\n",
" 1. Evaluate BSSN RHS expressions \n",
" * [**NRPy+ tutorial on BSSN right-hand sides**](Tutorial-BSSN_time_evolution-BSSN_RHSs.ipynb) ([**BSSN Introduction Notebook**](Tutorial-BSSN_formulation.ipynb))\n",
" * [**NRPy+ tutorial on BSSN gauge condition right-hand sides**](Tutorial-BSSN_time_evolution-BSSN_gauge_RHSs.ipynb) \n",
" 1. Apply singular, curvilinear coordinate boundary conditions [*a la* the SENR/NRPy+ paper](https://arxiv.org/abs/1712.07658)\n",
" * [**NRPy+ tutorial on setting up singular, curvilinear boundary conditions**](Tutorial-Start_to_Finish-Curvilinear_BCs.ipynb)\n",
" 1. Enforce constraint on conformal 3-metric: $\\det{\\bar{\\gamma}_{ij}}=\\det{\\hat{\\gamma}_{ij}}$ \n",
" * [**NRPy+ tutorial on enforcing $\\det{\\bar{\\gamma}_{ij}}=\\det{\\hat{\\gamma}_{ij}}$ constraint**](Tutorial-BSSN_enforcing_determinant_gammabar_equals_gammahat_constraint.ipynb)\n",
"1. Repeat the above steps at two numerical resolutions to confirm convergence to zero."
]
},
{
"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](#initializenrpy): Set core NRPy+ parameters for numerical grids and reference metric\n",
"1. [Step 2](#ccodegen): Generate C code kernels for BSSN expressions, in parallel if possible\n",
" 1. [Step 2.a](#rfm_ccodegen): Generate C code kernels for reference metric\n",
"1. [Step 3](#cparams_rfm_and_domainsize): Set `free_parameters.h`; also output C codes needed for declaring and setting Cparameters\n",
"1. [Step 4](#bc_functs): Set up boundary condition functions for chosen singular, curvilinear coordinate system\n",
"1. [Step 5](#mainc): `BrillLindquist_Playground`: The C code `main()` function\n",
"1. [Step 6](#compileexec): Compile generated C codes & perform the black hole collision calculation\n",
"1. [Step 7](#visualize): Visualize the output!\n",
" 1. [Step 7.a](#installdownload): Install `scipy` and download `ffmpeg` if they are not yet installed/downloaded\n",
" 1. [Step 7.b](#genimages): Generate images for visualization animation\n",
" 1. [Step 7.c](#genvideo): Generate visualization animation\n",
"1. [Step 8](#convergence): Plot the numerical error at the end of the simulation, and confirm that it converges to zero with increasing numerical resolution (sampling)\n",
"1. [Step 9](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 1: Set core NRPy+ parameters for numerical grids and reference metric \\[Back to [top](#toc)\\]\n",
"$$\\label{initializenrpy}$$"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"execution": {
"iopub.execute_input": "2021-09-23T22:08:25.778737Z",
"iopub.status.busy": "2021-09-23T22:08:25.774943Z",
"iopub.status.idle": "2021-09-23T22:08:26.199113Z",
"shell.execute_reply": "2021-09-23T22:08:26.198629Z"
}
},
"outputs": [
{
"ename": "ModuleNotFoundError",
"evalue": "No module named 'sympy'",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mModuleNotFoundError\u001b[0m Traceback (most recent call last)",
"Cell \u001b[0;32mIn[1], line 2\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[38;5;66;03m# Step P1: Import needed NRPy+ core modules:\u001b[39;00m\n\u001b[0;32m----> 2\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01moutputC\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m add_to_Cfunction_dict, outC_function_master_list \u001b[38;5;66;03m# NRPy+: Core C code output module\u001b[39;00m\n\u001b[1;32m 3\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mfinite_difference\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01mfin\u001b[39;00m \u001b[38;5;66;03m# NRPy+: Finite difference C code generation module\u001b[39;00m\n\u001b[1;32m 4\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mNRPy_param_funcs\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01mpar\u001b[39;00m \u001b[38;5;66;03m# NRPy+: Parameter interface\u001b[39;00m\n",
"File \u001b[0;32m~/Documents/VS_Code/Etienne/forked/nrpytutorial/outputC.py:18\u001b[0m\n\u001b[1;32m 12\u001b[0m __all__ \u001b[38;5;241m=\u001b[39m [\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mlhrh\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124moutCparams\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mnrpyAbs\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124msuperfast_uniq\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mcheck_if_string__error_if_not\u001b[39m\u001b[38;5;124m'\u001b[39m,\n\u001b[1;32m 13\u001b[0m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124moutputC\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mparse_outCparams_string\u001b[39m\u001b[38;5;124m'\u001b[39m,\n\u001b[1;32m 14\u001b[0m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124moutC_NRPy_basic_defines_h_dict\u001b[39m\u001b[38;5;124m'\u001b[39m,\n\u001b[1;32m 15\u001b[0m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124moutC_function_prototype_dict\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124moutC_function_dict\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mCfunction\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124madd_to_Cfunction_dict\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124moutCfunction\u001b[39m\u001b[38;5;124m'\u001b[39m]\n\u001b[1;32m 17\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mloop\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01mlp\u001b[39;00m \u001b[38;5;66;03m# NRPy+: C code loop interface\u001b[39;00m\n\u001b[0;32m---> 18\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mNRPy_param_funcs\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01mpar\u001b[39;00m \u001b[38;5;66;03m# NRPy+: parameter interface\u001b[39;00m\n\u001b[1;32m 19\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mSIMD\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m expr_convert_to_SIMD_intrins \u001b[38;5;66;03m# NRPy+: SymPy expression => SIMD intrinsics interface\u001b[39;00m\n\u001b[1;32m 20\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mcse_helpers\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m cse_preprocess,cse_postprocess \u001b[38;5;66;03m# NRPy+: CSE preprocessing and postprocessing\u001b[39;00m\n",
"File \u001b[0;32m~/Documents/VS_Code/Etienne/forked/nrpytutorial/NRPy_param_funcs.py:10\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[38;5;66;03m# As documented in the NRPy+ tutorial module\u001b[39;00m\n\u001b[1;32m 2\u001b[0m \u001b[38;5;66;03m# Tutorial-Coutput__Parameter_Interface.ipynb\u001b[39;00m\n\u001b[1;32m 3\u001b[0m \u001b[38;5;66;03m# this core NRPy+ module is used for\u001b[39;00m\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 7\u001b[0m \u001b[38;5;66;03m# Author: Zachariah B. Etienne\u001b[39;00m\n\u001b[1;32m 8\u001b[0m \u001b[38;5;66;03m# zachetie **at** gmail **dot* com\u001b[39;00m\n\u001b[0;32m---> 10\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01msympy\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01msp\u001b[39;00m \u001b[38;5;66;03m# Import SymPy\u001b[39;00m\n\u001b[1;32m 11\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01msys\u001b[39;00m \u001b[38;5;66;03m# Standard Python: OS-independent system functions\u001b[39;00m\n\u001b[1;32m 12\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mcollections\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m namedtuple \u001b[38;5;66;03m# Standard Python: Enable namedtuple data type\u001b[39;00m\n",
"\u001b[0;31mModuleNotFoundError\u001b[0m: No module named 'sympy'"
]
}
],
"source": [
"# Step P1: Import needed NRPy+ core modules:\n",
"from outputC import add_to_Cfunction_dict, outC_function_master_list # NRPy+: Core C code output module\n",
"import finite_difference as fin # NRPy+: Finite difference C code generation module\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 reference_metric as rfm # NRPy+: Reference metric support\n",
"import cmdline_helper as cmd # NRPy+: Multi-platform Python command-line interface\n",
"from pickling import unpickle_NRPy_env # NRPy+: Pickle/unpickle NRPy+ environment, for parallel codegen\n",
"import shutil, os, sys # Standard Python modules for multiplatform OS-level functions, benchmarking\n",
"\n",
"# Step P2: Create C code output directory:\n",
"Ccodesrootdir = os.path.join(\"BSSN_Two_BHs_Collide_Ccodes\")\n",
"# First remove C code output directory if it exists\n",
"# Courtesy https://stackoverflow.com/questions/303200/how-do-i-remove-delete-a-folder-that-is-not-empty\n",
"shutil.rmtree(Ccodesrootdir, ignore_errors=True)\n",
"# Then create a fresh directory\n",
"cmd.mkdir(Ccodesrootdir)\n",
"\n",
"# Step P3: Create executable output directory:\n",
"outdir = os.path.join(Ccodesrootdir, \"output\")\n",
"cmd.mkdir(outdir)\n",
"\n",
"# Step 1.a: Enable SIMD-optimized code?\n",
"# I.e., generate BSSN and Ricci C code kernels using SIMD-vectorized\n",
"# compiler intrinsics, which *greatly improve the code's performance*,\n",
"# though at the expense of making the C-code kernels less\n",
"# human-readable.\n",
"# * Important note in case you wish to modify the BSSN/Ricci kernels\n",
"# here by adding expressions containing transcendental functions\n",
"# (e.g., certain scalar fields):\n",
"# Note that SIMD-based transcendental function intrinsics are not\n",
"# supported by the default installation of gcc or clang (you will\n",
"# need to use e.g., the SLEEF library from sleef.org, for this\n",
"# purpose). The Intel compiler suite does support these intrinsics\n",
"# however without the need for external libraries.\n",
"enable_SIMD = True\n",
"\n",
"# Step 1.b: Enable reference metric precomputation.\n",
"enable_rfm_precompute = True\n",
"\n",
"if enable_SIMD and not enable_rfm_precompute:\n",
" print(\"ERROR: SIMD does not currently handle transcendental functions,\\n\")\n",
" print(\" like those found in rfmstruct (rfm_precompute).\\n\")\n",
" print(\" Therefore, enable_SIMD==True and enable_rfm_precompute==False\\n\")\n",
" print(\" is not supported.\\n\")\n",
" sys.exit(1)\n",
"\n",
"# Step 1.c: Enable \"FD functions\". In other words, all finite-difference stencils\n",
"# will be output as inlined static functions. This is essential for\n",
"# compiling highly complex FD kernels with using certain versions of GCC;\n",
"# GCC 10-ish will choke on BSSN FD kernels at high FD order, sometimes\n",
"# taking *hours* to compile. Unaffected GCC versions compile these kernels\n",
"# in seconds. FD functions do not slow the code performance, but do add\n",
"# another header file to the C source tree.\n",
"# With gcc 7.5.0, enable_FD_functions=True decreases performance by 10%\n",
"enable_FD_functions = True\n",
"\n",
"# Step 2: Set some core parameters, including CoordSystem MoL timestepping algorithm,\n",
"# FD order, floating point precision, and CFL factor:\n",
"# Choices are: Spherical, SinhSpherical, SinhSphericalv2, Cylindrical, SinhCylindrical,\n",
"# SymTP, SinhSymTP\n",
"CoordSystem = \"Spherical\"\n",
"par.set_parval_from_str(\"reference_metric::CoordSystem\", CoordSystem)\n",
"rfm.reference_metric()\n",
"\n",
"# Step 2.a: Set outer boundary condition\n",
"outer_bc_type = \"RADIATION_OUTER_BCS\" # can be EXTRAPOLATION_OUTER_BCS or RADIATION_OUTER_BCS\n",
"radiation_BC_FD_order = 2\n",
"\n",
"# Step 2.b: Set defaults for Coordinate system parameters.\n",
"# These are perhaps the most commonly adjusted parameters,\n",
"# so we enable modifications at this high level.\n",
"\n",
"# domain_size sets the default value for:\n",
"# * Spherical's params.RMAX\n",
"# * SinhSpherical*'s params.AMAX\n",
"# * Cartesians*'s -params.{x,y,z}min & .{x,y,z}max\n",
"# * Cylindrical's -params.ZMIN & .{Z,RHO}MAX\n",
"# * SinhCylindrical's params.AMPL{RHO,Z}\n",
"# * *SymTP's params.AMAX\n",
"domain_size = 7.5 # Needed for all coordinate systems.\n",
"\n",
"# sinh_width sets the default value for:\n",
"# * SinhSpherical's params.SINHW\n",
"# * SinhCylindrical's params.SINHW{RHO,Z}\n",
"# * SinhSymTP's params.SINHWAA\n",
"sinh_width = 0.4 # If Sinh* coordinates chosen\n",
"\n",
"# sinhv2_const_dr sets the default value for:\n",
"# * SinhSphericalv2's params.const_dr\n",
"# * SinhCylindricalv2's params.const_d{rho,z}\n",
"sinhv2_const_dr = 0.05 # If Sinh*v2 coordinates chosen\n",
"\n",
"# SymTP_bScale sets the default value for:\n",
"# * SinhSymTP's params.bScale\n",
"SymTP_bScale = 0.5 # If SymTP chosen\n",
"\n",
"# Step 2.c: Set the order of spatial and temporal derivatives;\n",
"# the core data type, and the CFL factor.\n",
"# RK_method choices include: Euler, \"RK2 Heun\", \"RK2 MP\", \"RK2 Ralston\", RK3, \"RK3 Heun\", \"RK3 Ralston\",\n",
"# SSPRK3, RK4, DP5, DP5alt, CK5, DP6, L6, DP8\n",
"RK_method = \"RK4\"\n",
"FD_order = 4 # Finite difference order: even numbers only, starting with 2. 12 is generally unstable\n",
"REAL = \"double\" # Best to use double here.\n",
"default_CFL_FACTOR= 0.5 # (GETS OVERWRITTEN IF SPECIFIED AT COMMAND LINE.)\n",
" # In pure axisymmetry (symmetry_axes = 2 below) 1.0 works fine. Otherwise 0.5 or lower."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {
"iopub.execute_input": "2021-09-23T22:08:26.203895Z",
"iopub.status.busy": "2021-09-23T22:08:26.203369Z",
"iopub.status.idle": "2021-09-23T22:08:26.205606Z",
"shell.execute_reply": "2021-09-23T22:08:26.205163Z"
}
},
"outputs": [],
"source": [
"# Step 5: Set the finite differencing order to FD_order (set above).\n",
"par.set_parval_from_str(\"finite_difference::FD_CENTDERIVS_ORDER\", FD_order)\n",
"\n",
"# Directory for reference_metric precomputation header files:\n",
"rfm_precompute_Ccode_outdir = os.path.join(Ccodesrootdir, \"rfm_files/\")\n",
"if enable_rfm_precompute:\n",
" cmd.mkdir(os.path.join(Ccodesrootdir, \"rfm_files/\"))\n",
" par.set_parval_from_str(\"reference_metric::rfm_precompute_Ccode_outdir\", rfm_precompute_Ccode_outdir)\n",
"\n",
"# Step 6: Copy SIMD/SIMD_intrinsics.h to $Ccodesrootdir/SIMD/SIMD_intrinsics.h\n",
"if enable_SIMD:\n",
" cmd.mkdir(os.path.join(Ccodesrootdir,\"SIMD\"))\n",
" shutil.copy(os.path.join(\"SIMD/\")+\"SIMD_intrinsics.h\",os.path.join(Ccodesrootdir,\"SIMD/\"))\n",
"\n",
"# Step 7: Set finite_difference::enable_FD_functions appropriately. Defaults to False\n",
"if enable_FD_functions:\n",
" par.set_parval_from_str(\"finite_difference::enable_FD_functions\", enable_FD_functions)\n",
"\n",
"# Step 8: If enable_SIMD, then copy SIMD/SIMD_intrinsics.h to $Ccodesrootdir/SIMD/SIMD_intrinsics.h\n",
"cmd.mkdir(os.path.join(Ccodesrootdir,\"SIMD\"))\n",
"if enable_SIMD:\n",
" shutil.copy(os.path.join(\"SIMD\", \"SIMD_intrinsics.h\"), os.path.join(Ccodesrootdir, \"SIMD\"))\n",
"\n",
"# Step 9: Set the direction=2 (phi) axis to be the symmetry axis; i.e.,\n",
"# axis \"2\", corresponding to the i2 direction.\n",
"# This sets all spatial derivatives in the phi direction to zero.\n",
"par.set_parval_from_str(\"indexedexp::symmetry_axes\",\"2\")\n",
"OMP_pragma_on = \"i1\" # structure OpenMP loops to parallelize, not over i2 (phi direction), but i1 (theta direction)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Step 10: Generate Runge-Kutta-based (RK-based) timestepping code.\n",
"# As described above the Table of Contents, this is a 3-step process:\n",
"# 10.A: Evaluate RHSs (RHS_string)\n",
"# 10.B: Apply boundary conditions (post_RHS_string, pt 1)\n",
"# 10.C: Enforce det(gammabar) = det(gammahat) constraint (post_RHS_string, pt 2)\n",
"import MoLtimestepping.MoL as MoL\n",
"# from MoLtimestepping.RK_Butcher_Table_Dictionary import Butcher_dict\n",
"# RK_order = Butcher_dict[RK_method][1]\n",
"RHS_string = \"\"\"\n",
"Ricci_eval(params, rfmstruct, RK_INPUT_GFS, auxevol_gfs);\n",
"rhs_eval( params, rfmstruct, auxevol_gfs, RK_INPUT_GFS, RK_OUTPUT_GFS);\"\"\"\n",
"RHS_string += \"\"\"\n",
"if(params->outer_bc_type == RADIATION_OUTER_BCS)\n",
" apply_bcs_outerradiation_and_inner(params, bcstruct, griddata->xx,\n",
" gridfunctions_wavespeed,gridfunctions_f_infinity,\n",
" RK_INPUT_GFS, RK_OUTPUT_GFS);\"\"\"\n",
"\n",
"# Extrapolation BCs are applied to the evolved gridfunctions themselves after the MoL update\n",
"post_RHS_string = \"\"\"if(params->outer_bc_type == EXTRAPOLATION_OUTER_BCS)\n",
" apply_bcs_outerextrap_and_inner(params, bcstruct, RK_OUTPUT_GFS);\n",
"\"\"\"\n",
"post_RHS_string += \"enforce_detgammahat_constraint(params, rfmstruct, RK_OUTPUT_GFS);\\n\"\n",
"\n",
"if not enable_rfm_precompute:\n",
" RHS_string = RHS_string.replace(\"rfmstruct\", \"xx\")\n",
" post_RHS_string = post_RHS_string.replace(\"rfmstruct\", \"xx\")\n",
"\n",
"MoL.register_C_functions_and_NRPy_basic_defines(RK_method,\n",
" RHS_string=RHS_string, post_RHS_string=post_RHS_string,\n",
" enable_rfm=enable_rfm_precompute, enable_curviBCs=True, enable_SIMD=False)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 2: Generate C code kernels for BSSN expressions, in parallel if possible \\[Back to [top](#toc)\\]\n",
"$$\\label{ccodegen}$$\n",
"\n",
"In the following code cell, we create a list of Python functions, which each registers a single C code function in `outputC`'s `outC_function_dict` dictionary. These Python functions are defined in \n",
"\n",
"1. [`BSSN.Initial_Data_Reader__BSSN_Converter`](../edit/BSSN.Initial_Data_Reader__BSSN_Converter.py); [\\[**tutorial**\\]](Tutorial-ADM_Initial_Data_Reader__BSSN_Converter.ipynb), which registers the C function for our \"universal\" initial data reader/converter `initial_data_reader__convert_ADM_Cartesian_to_BSSN()`.\n",
"1. the [`BSSN.BSSN_Ccodegen_library`](../edit/BSSN/BSSN_Ccodegen_library.py) NRPy+ module [\\[**tutorial**\\]](Tutorial-BSSN_time_evolution-C_codegen_library.ipynb), which contains Python functions for generating C code from symbolic expressions constructed within the following NRPy+ modules/tutorials:\n",
" 1. [BSSN/BSSN_constraints.py](../edit/BSSN/BSSN_constraints.py); [\\[**tutorial**\\]](Tutorial-BSSN_constraints.ipynb): Hamiltonian constraint in BSSN curvilinear basis/coordinates\n",
" 1. [BSSN/BSSN_RHSs.py](../edit/BSSN/BSSN_RHSs.py); [\\[**tutorial**\\]](Tutorial-BSSN_time_evolution-BSSN_RHSs.ipynb): Generates the right-hand sides for the BSSN evolution equations in singular, curvilinear coordinates\n",
" 1. [BSSN/BSSN_gauge_RHSs.py](../edit/BSSN/BSSN_gauge_RHSs.py); [\\[**tutorial**\\]](Tutorial-BSSN_time_evolution-BSSN_gauge_RHSs.ipynb): Generates the right-hand sides for the BSSN gauge evolution equations in singular, curvilinear coordinates\n",
" 1. [BSSN/Enforce_Detgammahat_Constraint.py](../edit/BSSN/Enforce_Detgammahat_Constraint.py); [**tutorial**](Tutorial-BSSN_enforcing_determinant_gammabar_equals_gammahat_constraint.ipynb): Generates symbolic expressions for enforcing the $\\det{\\bar{\\gamma}}=\\det{\\hat{\\gamma}}$ constraint\n",
"\n",
"\n",
"Next, from within a `multiprocessing` environment, we then call all the Python C-code generation functions in this list in parallel (if `multiprocessing` is supported). This is quite useful, as these functions take several seconds to complete.\n",
"\n",
"Within each `multiprocessing` process, the current NRPy+ environment is cloned, and a new function is registered to the `outC_function_dict` dictionary. Thus when each process completes, it contains a unique NRPy+ environment, with only its function registered. We address this by saving each process' NRPy+ environment and sending it back in a common binary format known as a `pickle`, using NRPy+'s [`pickling`](../edit/pickling.py) module. The environments are combined in an unpickling such that all functions exist within the same `outC_function_dict` dictionary.\n",
"\n",
"To make the current environment fully consistent, we call `reference_metric.py` to register all its associated C functions (stored in globals) and contributions to `NRPy_basic_defines.h`."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# First set up the symbolic expressions for Brill-Lindquist initial data\n",
"import BSSN.BrillLindquist as BLID\n",
"BLID.BrillLindquist()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next we register the `ID_function()` needed for setting ADM spacetime quantities within our \"universal\" initial data reader/converter, as described in [this tutorial](Tutorial-ADM_Initial_Data_Reader__BSSN_Converter.ipynb). Appropriately, we refer to this `ID_function()` by the name `BrillLindquist()`."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import BSSN.ADM_Initial_Data_Reader__BSSN_Converter as IDread\n",
"\n",
"IDread.add_to_Cfunction_dict_exact_ADM_ID_function(\"BrillLindquist\", \"Cartesian\",\n",
" BLID.alpha, BLID.betaU, BLID.BU, BLID.gammaDD, BLID.KDD)\n",
"IDread.register_NRPy_basic_defines(include_T4UU=False)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {
"iopub.execute_input": "2021-09-23T22:08:26.232538Z",
"iopub.status.busy": "2021-09-23T22:08:26.231967Z",
"iopub.status.idle": "2021-09-23T22:08:41.510761Z",
"shell.execute_reply": "2021-09-23T22:08:41.510225Z"
}
},
"outputs": [],
"source": [
"# Step 2: Generate C code kernels for BSSN expressions, in parallel if possible;\n",
"import BSSN.BSSN_Ccodegen_library as BCL\n",
"# Step 2.a: Create a list of functions we wish to evaluate in parallel (if possible)\n",
"# Create lists for all BSSN functions\n",
"codegen_funcs = [IDread.add_to_Cfunction_dict_initial_data_reader__convert_ADM_Sph_or_Cart_to_BSSN]\n",
"codegen_funcs.append(BCL.add_rhs_eval_to_Cfunction_dict)\n",
"codegen_funcs.append(BCL.add_Ricci_eval_to_Cfunction_dict)\n",
"codegen_funcs.append(BCL.add_BSSN_constraints_to_Cfunction_dict)\n",
"codegen_funcs.append(BCL.add_enforce_detgammahat_constraint_to_Cfunction_dict)\n",
"\n",
"# Step 2.b: Define master functions for parallelization.\n",
"# Note that lambdifying this doesn't work in Python 3\n",
"def master_func(arg):\n",
" if codegen_funcs[arg] == IDread.add_to_Cfunction_dict_initial_data_reader__convert_ADM_Sph_or_Cart_to_BSSN:\n",
" ret = codegen_funcs[arg](addl_includes=None, rel_path_to_Cparams=os.path.join(\".\"),\n",
" input_Coord=\"Cartesian\", include_T4UU=False)\n",
" else:\n",
" if enable_rfm_precompute:\n",
" # We use rfm_precompute for all BSSN functions:\n",
" par.set_parval_from_str(\"reference_metric::enable_rfm_precompute\", \"True\")\n",
" rfm.reference_metric()\n",
"\n",
" if codegen_funcs[arg].__name__ == \"add_BSSN_constraints_to_Cfunction_dict\":\n",
" ret = codegen_funcs[arg](includes=[\"NRPy_basic_defines.h\"],\n",
" rel_path_to_Cparams=os.path.join(\".\"), output_H_only=True,\n",
" enable_rfm_precompute=enable_rfm_precompute, enable_SIMD=enable_SIMD,\n",
" OMP_pragma_on=OMP_pragma_on)\n",
" elif codegen_funcs[arg].__name__ == \"add_rhs_eval_to_Cfunction_dict\" or \\\n",
" codegen_funcs[arg].__name__ == \"add_Ricci_eval_to_Cfunction_dict\":\n",
" ret = codegen_funcs[arg](includes=[\"NRPy_basic_defines.h\"],\n",
" rel_path_to_Cparams=os.path.join(\".\"),\n",
" enable_rfm_precompute=enable_rfm_precompute, enable_SIMD=enable_SIMD,\n",
" OMP_pragma_on=OMP_pragma_on)\n",
" elif codegen_funcs[arg].__name__ == \"add_enforce_detgammahat_constraint_to_Cfunction_dict\":\n",
" ret = codegen_funcs[arg](includes=[\"NRPy_basic_defines.h\"],\n",
" rel_path_to_Cparams=os.path.join(\".\"),\n",
" enable_rfm_precompute=enable_rfm_precompute, OMP_pragma_on=OMP_pragma_on)\n",
" else:\n",
" print(\"ERROR: DID NOT RECOGNIZE FUNCTION \" + codegen_funcs[arg].__name__ + \"\\n\")\n",
" sys.exit(1)\n",
" if enable_rfm_precompute:\n",
" par.set_parval_from_str(\"reference_metric::enable_rfm_precompute\", \"False\")\n",
" rfm.ref_metric__hatted_quantities()\n",
" return ret\n",
"\n",
"\n",
"NRPyEnvVars = []\n",
"raised_exception = False\n",
"try:\n",
" if os.name == 'nt':\n",
" # It's a mess to get working in Windows, so we don't bother. :/\n",
" # https://medium.com/@grvsinghal/speed-up-your-python-code-using-multiprocessing-on-windows-and-jupyter-or-ipython-2714b49d6fac\n",
" raise Exception(\"Parallel codegen currently not available in certain environments, e.g., Windows\")\n",
"\n",
" # Step 2.d: Import the multiprocessing module.\n",
" import multiprocessing\n",
"\n",
" # Step 2.e: Evaluate list of functions in parallel if possible;\n",
" # otherwise fallback to serial evaluation:\n",
" pool = multiprocessing.Pool()\n",
" NRPyEnvVars.append(pool.map(master_func, range(len(codegen_funcs))))\n",
" pool.terminate()\n",
" pool.join()\n",
"except:\n",
" print(\"FAILED PARALLEL CODEGEN!\")\n",
" NRPyEnvVars = [] # Reset, as pickling/unpickling unnecessary for serial codegen (see next line)\n",
"\n",
" # Steps 2.d-e, alternate: As fallback, evaluate functions in serial.\n",
" # This will happen on Android and Windows systems\n",
" for i, func in enumerate(codegen_funcs):\n",
" master_func(i)\n",
" raised_exception = True\n",
"\n",
"outCfunc_master_list = outC_function_master_list\n",
"if not raised_exception:\n",
" outCfunc_master_list = unpickle_NRPy_env(NRPyEnvVars)\n",
" for el in outCfunc_master_list:\n",
" if el not in outC_function_master_list: # in case there are duplicate funcs, which can happen\n",
" # if finite_difference_functions = True\n",
" outC_function_master_list += [el]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"## Step 2.a: Generate C code kernels for reference metric \\[Back to [top](#toc)\\]\n",
"$$\\label{rfm_ccodegen}$$\n",
"\n",
"In the [reference_metric](../edit/reference_metric.py) NRPy+ module, `register_C_functions_and_NRPy_basic_defines()` registers the following C functions to `outC_Cfunction_dict`:\n",
"\n",
"1. `find_timestep()`: Finds the minimum spacing between adjacent gridpoints on our numerical grid $\\min(ds_i)$, and sets the timestep according to the [CFL](https://en.wikipedia.org/w/index.php?title=Courant%E2%80%93Friedrichs%E2%80%93Lewy_condition&oldid=806430673) 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.\n",
"1. `xx_to_Cart()`: Input = uniformly sampled coordinate xx0,xx1,xx2 (e.g., r,theta,phi in Spherical coordinates). Output = Cartesian coordinate (x,y,z).\n",
"1. `set_Nxx_dxx_invdx_params__and__xx()`: Sets `Nxx{0,1,2}`, `Nxx_plus_2NGHOSTS{0,1,2}`, `dxx{0,1,2}`, and `invdx{0,1,2}`; and defines `xx[3][]`.\n",
"1. `Cart_to_xx_and_nearest_i0i1i2()`: Input = Cartesian coordinate (x,y,z). Output = uniformly sampled coordinate xx0,xx1,xx2 (e.g., r,theta,phi in Spherical coordinates), as well as corresponding grid index `i0,i1,i2`."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {
"iopub.execute_input": "2021-09-23T22:08:41.574056Z",
"iopub.status.busy": "2021-09-23T22:08:41.548577Z",
"iopub.status.idle": "2021-09-23T22:08:41.938615Z",
"shell.execute_reply": "2021-09-23T22:08:41.938185Z"
}
},
"outputs": [],
"source": [
"# Generate & register C function set_Nxx_dxx_invdx_params__and__xx()\n",
"# Generate & register C function xx_to_Cart() for\n",
"# (the mapping from xx->Cartesian) for the chosen\n",
"# CoordSystem:\n",
"# Generate & register the find_timestep() function\n",
"\n",
"# Sets reference_metric globals: NRPy_basic_defines_str, rfm_struct__malloc, rfm_struct__define, rfm_struct__freemem\n",
"if enable_rfm_precompute:\n",
" par.set_parval_from_str(\"reference_metric::rfm_precompute_Ccode_outdir\", rfm_precompute_Ccode_outdir)\n",
" par.set_parval_from_str(\"reference_metric::enable_rfm_precompute\", \"True\")\n",
" par.set_parval_from_str(\"reference_metric::rfm_precompute_to_Cfunctions_and_NRPy_basic_defines\", \"True\")\n",
" rfm.reference_metric()\n",
"\n",
"rfm.register_C_functions(enable_rfm_precompute=enable_rfm_precompute, use_unit_wavespeed_for_find_timestep=True)\n",
"rfm.register_NRPy_basic_defines(enable_rfm_precompute=enable_rfm_precompute)\n",
"\n",
"if enable_rfm_precompute:\n",
" par.set_parval_from_str(\"reference_metric::enable_rfm_precompute\", \"False\")\n",
" rfm.ref_metric__hatted_quantities()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 3: Set `free_parameters.h`; also output C codes needed for declaring and setting Cparameters \\[Back to [top](#toc)\\]\n",
"$$\\label{cparams_rfm_and_domainsize}$$\n",
"\n",
"First 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."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {
"iopub.execute_input": "2021-09-23T22:08:41.942239Z",
"iopub.status.busy": "2021-09-23T22:08:41.941715Z",
"iopub.status.idle": "2021-09-23T22:08:41.943573Z",
"shell.execute_reply": "2021-09-23T22:08:41.943985Z"
}
},
"outputs": [],
"source": [
"# Step 3.e.i: Set free_parameters.h\n",
"outstr = r\"\"\"// Set free-parameter values.\n",
"\n",
"// Outer boundary condition choice:\n",
"params.outer_bc_type = \"\"\"+outer_bc_type+r\"\"\";\n",
"\n",
"// Set the default CFL Factor. Can be overwritten at command line.\n",
"REAL CFL_FACTOR = \"\"\"+str(default_CFL_FACTOR)+r\"\"\";\n",
"\n",
"// Set free-parameter values for BSSN evolution:\n",
"params.eta = 1.0;\n",
"\n",
"// Set free parameters for the (Brill-Lindquist) initial data\n",
"params.BH1_posn_x = 0.0; params.BH1_posn_y = 0.0; params.BH1_posn_z =+0.5;\n",
"params.BH2_posn_x = 0.0; params.BH2_posn_y = 0.0; params.BH2_posn_z =-0.5;\n",
"params.BH1_mass = 0.5; params.BH2_mass = 0.5;\n",
"\"\"\"\n",
"\n",
"# Append to $Ccodesrootdir/free_parameters.h reference metric parameters based on generic\n",
"# domain_size,sinh_width,sinhv2_const_dr,SymTP_bScale,\n",
"# parameters set above.\n",
"outstr += rfm.out_default_free_parameters_for_rfm(\"returnstring\",\n",
" domain_size,sinh_width,sinhv2_const_dr,SymTP_bScale)\n",
"with open(os.path.join(Ccodesrootdir,\"free_parameters.h\"),\"w\") as file:\n",
" file.write(outstr.replace(\"params.\", \"griddata.params.\"))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 4: Set up boundary condition functions for chosen singular, curvilinear coordinate system \\[Back to [top](#toc)\\]\n",
"$$\\label{bc_functs}$$\n",
"\n",
"Next apply singular, curvilinear coordinate boundary conditions [as documented in the corresponding NRPy+ tutorial notebook](Tutorial-Start_to_Finish-Curvilinear_BCs.ipynb)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {
"iopub.execute_input": "2021-09-23T22:08:41.952190Z",
"iopub.status.busy": "2021-09-23T22:08:41.951777Z",
"iopub.status.idle": "2021-09-23T22:08:42.181582Z",
"shell.execute_reply": "2021-09-23T22:08:42.181905Z"
}
},
"outputs": [],
"source": [
"import CurviBoundaryConditions.CurviBoundaryConditions as CBC\n",
"CBC.CurviBoundaryConditions_register_NRPy_basic_defines()\n",
"CBC.CurviBoundaryConditions_register_C_functions(radiation_BC_FD_order=radiation_BC_FD_order)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 5: The C code `main()` function for `BrillLindquist_Playground` \\[Back to [top](#toc)\\]\n",
"$$\\label{mainc}$$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {
"iopub.execute_input": "2021-09-23T22:08:42.191346Z",
"iopub.status.busy": "2021-09-23T22:08:42.183557Z",
"iopub.status.idle": "2021-09-23T22:08:42.193048Z",
"shell.execute_reply": "2021-09-23T22:08:42.192711Z"
}
},
"outputs": [],
"source": [
"def add_to_Cfunction_dict_main__BrillLindquist_Playground():\n",
" includes = [\"NRPy_basic_defines.h\", \"NRPy_function_prototypes.h\", \"time.h\"]\n",
" desc = \"\"\"// main() function:\n",
"// Step 0: Read command-line input, set up grid structure, allocate memory for gridfunctions, set up coordinates\n",
"// Step 1: Set up initial data to an exact solution\n",
"// Step 2: Start the timer, for keeping track of how fast the simulation is progressing.\n",
"// Step 3: Integrate the initial data forward in time using the chosen RK-like Method of\n",
"// Lines timestepping algorithm, and output periodic simulation diagnostics\n",
"// Step 3.a: Output 2D data file periodically, for visualization\n",
"// Step 3.b: Step forward one timestep (t -> t+dt) in time using\n",
"// chosen RK-like MoL timestepping algorithm\n",
"// Step 3.c: If t=t_final, output conformal factor & Hamiltonian\n",
"// constraint violation to 2D data file\n",
"// Step 3.d: Progress indicator printing to stderr\n",
"// Step 4: Free all allocated memory\n",
"\"\"\"\n",
" c_type = \"int\"\n",
" name = \"main\"\n",
" params = \"int argc, const char *argv[]\"\n",
" body = r\"\"\" griddata_struct griddata;\n",
" set_Cparameters_to_default(&griddata.params);\n",
"\n",
" // Step 0.a: Set free parameters, overwriting Cparameters defaults\n",
" // by hand or with command-line input, as desired.\n",
"#include \"free_parameters.h\"\n",
"\n",
" // Step 0.b: Read command-line input, error out if nonconformant\n",
" if((argc != 4 && argc != 5) || atoi(argv[1]) < NGHOSTS || atoi(argv[2]) < NGHOSTS || atoi(argv[3]) < 2 /* FIXME; allow for axisymmetric sims */) {\n",
" fprintf(stderr,\"Error: Expected three command-line arguments: ./BrillLindquist_Playground Nx0 Nx1 Nx2,\\n\");\n",
" fprintf(stderr,\"where Nx[0,1,2] is the number of grid points in the 0, 1, and 2 directions.\\n\");\n",
" fprintf(stderr,\"Nx[] MUST BE larger than NGHOSTS (= %d)\\n\",NGHOSTS);\n",
" exit(1);\n",
" }\n",
" if(argc == 5) {\n",
" CFL_FACTOR = strtod(argv[4],NULL);\n",
" if(CFL_FACTOR > 0.5 && atoi(argv[3])!=2) {\n",
" fprintf(stderr,\"WARNING: CFL_FACTOR was set to %e, which is > 0.5.\\n\",CFL_FACTOR);\n",
" fprintf(stderr,\" This will generally only be stable if the simulation is purely axisymmetric\\n\");\n",
" fprintf(stderr,\" However, Nx2 was set to %d>2, which implies a non-axisymmetric simulation\\n\",atoi(argv[3]));\n",
" }\n",
" }\n",
" // Step 0.c: Set up numerical grid structure, first in space...\n",
" const int Nxx[3] = { atoi(argv[1]), atoi(argv[2]), atoi(argv[3]) };\n",
" if(Nxx[0]%2 != 0 || Nxx[1]%2 != 0 || Nxx[2]%2 != 0) {\n",
" fprintf(stderr,\"Error: Cannot guarantee a proper cell-centered grid if number of grid cells not set to even number.\\n\");\n",
" fprintf(stderr,\" For example, in case of angular directions, proper symmetry zones will not exist.\\n\");\n",
" exit(1);\n",
" }\n",
"\n",
" // Step 0.d: Uniform coordinate grids are stored to *xx[3]\n",
" // Step 0.d.i: Set bcstruct\n",
" {\n",
" int EigenCoord;\n",
" EigenCoord = 1;\n",
" // Step 0.d.ii: Call set_Nxx_dxx_invdx_params__and__xx(), which sets\n",
" // params Nxx,Nxx_plus_2NGHOSTS,dxx,invdx, and xx[] for the\n",
" // chosen Eigen-CoordSystem.\n",
" set_Nxx_dxx_invdx_params__and__xx(EigenCoord, Nxx, &griddata.params, griddata.xx);\n",
" // Step 0.e: Find ghostzone mappings; set up bcstruct\n",
" bcstruct_set_up(&griddata.params, griddata.xx, &griddata.bcstruct);\n",
" // Step 0.e.i: Free allocated space for xx[][] array\n",
" for(int i=0;i<3;i++) free(griddata.xx[i]);\n",
"\n",
" // Step 0.f: Call set_Nxx_dxx_invdx_params__and__xx(), which sets\n",
" // params Nxx,Nxx_plus_2NGHOSTS,dxx,invdx, and xx[] for the\n",
" // chosen (non-Eigen) CoordSystem.\n",
" EigenCoord = 0;\n",
" set_Nxx_dxx_invdx_params__and__xx(EigenCoord, Nxx, &griddata.params, griddata.xx);\n",
" }\n",
"\n",
" // Step 0.g: Time coordinate parameters\n",
" griddata.params.t_final = domain_size; /* Final time is set so that at t=t_final,\n",
" * data at the origin have not been corrupted\n",
" * by the approximate outer boundary condition */\n",
"\n",
" // Step 0.h: Set timestep based on smallest proper distance between gridpoints and CFL factor\n",
" griddata.params.dt = find_timestep(&griddata.params, griddata.xx, CFL_FACTOR);\n",
" //fprintf(stderr,\"# Timestep set to = %e\\n\",(double)dt);\n",
" // Ntot = N_outputs * output_every_N -> output_every_N = Ntot/N_outputs; set N_outputs=800\n",
" int output_every_N = (int)((griddata.params.t_final/griddata.params.dt + 0.5) / 800.0);\n",
" if(output_every_N == 0) output_every_N = 1;\n",
"\n",
" // Step 0.i: Error out if the number of auxiliary gridfunctions outnumber evolved gridfunctions.\n",
" // This is a limitation of the RK method. You are always welcome to declare & allocate\n",
" // additional gridfunctions by hand.\n",
" if(NUM_AUX_GFS > NUM_EVOL_GFS) {\n",
" fprintf(stderr,\"Error: NUM_AUX_GFS > NUM_EVOL_GFS. Either reduce the number of auxiliary gridfunctions,\\n\");\n",
" fprintf(stderr,\" or allocate (malloc) by hand storage for *diagnostic_output_gfs. \\n\");\n",
" exit(1);\n",
" }\n",
"\n",
" // Step 0.j: Declare struct for gridfunctions and allocate memory for y_n_gfs gridfunctions\n",
" MoL_malloc_y_n_gfs(&griddata.params, &griddata.gridfuncs);\n",
"\"\"\"\n",
" if enable_rfm_precompute:\n",
" body += \"\"\"\n",
" // Step 0.k: Set up precomputed reference metric arrays\n",
" // Step 0.k.i: Allocate space for precomputed reference metric arrays.\n",
" rfm_precompute_rfmstruct_malloc(&griddata.params, &griddata.rfmstruct);\n",
"\n",
" // Step 0.k.ii: Define precomputed reference metric arrays.\n",
" rfm_precompute_rfmstruct_define(&griddata.params, griddata.xx, &griddata.rfmstruct);\"\"\"\n",
" body += r\"\"\"\n",
" // Step 0.l: Set up initial data to an exact solution (Brill-Lindquist)\n",
" ID_persist_struct ID_persist;\n",
" initial_data_reader__convert_ADM_Cartesian_to_BSSN(&griddata, &ID_persist, BrillLindquist);\n",
"\n",
" // Step 0.m: Allocate memory for non-y_n_gfs gridfunctions\n",
" MoL_malloc_non_y_n_gfs(&griddata.params, &griddata.gridfuncs);\n",
"\n",
" // Step 0.n: Apply boundary conditions, as the BSSN\n",
" // quantity lambda^i, defined using finite-\n",
" // difference derivatives, is undefined in\n",
" // ghost zones.\n",
" apply_bcs_outerextrap_and_inner(&griddata.params, &griddata.bcstruct, griddata.gridfuncs.y_n_gfs);\n",
" enforce_detgammahat_constraint(&griddata.params, &griddata.rfmstruct, griddata.gridfuncs.y_n_gfs);\n",
"\n",
" // Step 1: Start the timer, for keeping track of how fast the simulation is progressing.\n",
"#ifdef __linux__ // Use high-precision timer in Linux.\n",
" struct timespec start, end;\n",
" clock_gettime(CLOCK_REALTIME, &start);\n",
"#else // Resort to low-resolution, standards-compliant timer in non-Linux OSs\n",
" // http://www.cplusplus.com/reference/ctime/time/\n",
" time_t start_timer,end_timer;\n",
" time(&start_timer); // Resolution of one second...\n",
"#endif\n",
"\n",
" const int Nxx_plus_2NGHOSTS0 = griddata.params.Nxx_plus_2NGHOSTS0;\n",
" const int Nxx_plus_2NGHOSTS1 = griddata.params.Nxx_plus_2NGHOSTS1;\n",
" const int Nxx_plus_2NGHOSTS2 = griddata.params.Nxx_plus_2NGHOSTS2;\n",
"\n",
" // Step 3: Integrate the initial data forward in time using the chosen RK-like Method of\n",
" // Lines timestepping algorithm, and output periodic simulation diagnostics\n",
" while(griddata.params.time < griddata.params.t_final)\n",
" { // Main loop to progress forward in time.\n",
"\n",
" // Step 3.a: Output 2D data file periodically, for visualization\n",
" if(griddata.params.n%100 == 0) {\n",
" // Evaluate BSSN constraints (currently only Hamiltonian constraint violation computed).\n",
" Ricci_eval(&griddata.params, &griddata.rfmstruct, griddata.gridfuncs.y_n_gfs, griddata.gridfuncs.auxevol_gfs); // <- depends on Ricci.\n",
" BSSN_constraints(&griddata.params, &griddata.rfmstruct, griddata.gridfuncs.y_n_gfs, griddata.gridfuncs.auxevol_gfs, griddata.gridfuncs.diagnostic_output_gfs);\n",
"\n",
" char filename[100]; sprintf(filename,\"out%d-%08d.txt\",Nxx[0],griddata.params.n);\n",
" FILE *out2D = fopen(filename, \"w\");\n",
" LOOP_REGION(NGHOSTS,Nxx_plus_2NGHOSTS0-NGHOSTS,\n",
" NGHOSTS,Nxx_plus_2NGHOSTS1-NGHOSTS,\n",
" NGHOSTS,Nxx_plus_2NGHOSTS2-NGHOSTS) {\n",
" const int idx = IDX3S(i0,i1,i2);\n",
" REAL xCart[3]; xx_to_Cart(&griddata.params,griddata.xx,i0,i1,i2, xCart);\n",
" fprintf(out2D,\"%e %e %e %e\\n\",\n",
" xCart[1],xCart[2],\n",
" griddata.gridfuncs.y_n_gfs[IDX4ptS(CFGF,idx)],log10(fabs(griddata.gridfuncs.diagnostic_output_gfs[IDX4ptS(HGF,idx)])));\n",
" }\n",
" fclose(out2D);\n",
" }\n",
"\n",
" // Step 3.b: Step forward one timestep (t -> t+dt) in time using\n",
" // chosen RK-like MoL timestepping algorithm\n",
" MoL_step_forward_in_time(&griddata);\n",
"\n",
" // Step 3.c: If t=t_final, output conformal factor & Hamiltonian\n",
" // constraint violation to 2D data file\n",
" if(griddata.params.time > (griddata.params.t_final - griddata.params.dt)) {\n",
" // Evaluate BSSN constraints (currently only Hamiltonian constraint violation computed).\n",
" Ricci_eval(&griddata.params, &griddata.rfmstruct, griddata.gridfuncs.y_n_gfs, griddata.gridfuncs.auxevol_gfs); // <- depends on Ricci.\n",
" BSSN_constraints(&griddata.params, &griddata.rfmstruct, griddata.gridfuncs.y_n_gfs, griddata.gridfuncs.auxevol_gfs, griddata.gridfuncs.diagnostic_output_gfs);\n",
"\n",
" char filename[100]; sprintf(filename,\"out%d.txt\",Nxx[0]);\n",
" FILE *out2D = fopen(filename, \"w\");\n",
" LOOP_REGION(NGHOSTS,Nxx_plus_2NGHOSTS0-NGHOSTS,\n",
" NGHOSTS,Nxx_plus_2NGHOSTS1-NGHOSTS,\n",
" NGHOSTS,Nxx_plus_2NGHOSTS2-NGHOSTS) {\n",
" int idx = IDX3S(i0,i1,i2);\n",
" REAL xCart[3]; xx_to_Cart(&griddata.params,griddata.xx,i0,i1,i2, xCart);\n",
" fprintf(out2D,\"%e %e %e %e\\n\",\n",
" xCart[1],xCart[2],\n",
" griddata.gridfuncs.y_n_gfs[IDX4ptS(CFGF,idx)],log10(fabs(griddata.gridfuncs.diagnostic_output_gfs[IDX4ptS(HGF,idx)])));\n",
" }\n",
" fclose(out2D);\n",
" }\n",
" // Step 3.d: Progress indicator printing to stderr.\n",
" // Note that we're at the end of an iteration, so n=0\n",
" // indicates we are actually at the start of iteration 1.\n",
" if((griddata.params.n + 1) % 10 == 0 ||\n",
" (griddata.params.time > griddata.params.t_final + griddata.params.dt) ) {\n",
" // Step 3.d.i: Measure average time per iteration\n",
"#ifdef __linux__ // Use high-precision timer in Linux.\n",
" clock_gettime(CLOCK_REALTIME, &end);\n",
" const long long unsigned int time_in_ns = 1000000000L * (end.tv_sec - start.tv_sec) + end.tv_nsec - start.tv_nsec;\n",
"#else // Resort to low-resolution, standards-compliant timer in non-Linux OSs\n",
" time(&end_timer); // Resolution of one second...\n",
" REAL time_in_ns = difftime(end_timer,start_timer)*1.0e9+0.5; // Round up to avoid divide-by-zero.\n",
"#endif\n",
" const REAL s_per_iteration_avg = ((REAL)time_in_ns / (REAL)(griddata.params.n + 1)) / 1.0e9;\n",
" // We are at the n+1st iteration now ---------------------v -------^\n",
" const int iterations_remaining = (int)(griddata.params.t_final/griddata.params.dt + 0.5) - (griddata.params.n + 1);\n",
" const REAL seconds_remaining = (int)(s_per_iteration_avg * (REAL)iterations_remaining);\n",
" const int time_remaining__hrs = (int)(seconds_remaining / 3600.0);\n",
" const int time_remaining__mins = (int)(seconds_remaining / 60.0) % 60;\n",
" const int time_remaining__secs = (int)(seconds_remaining) - time_remaining__mins*60 - time_remaining__hrs*3600;\n",
" //const REAL num_RHS_pt_evals = (REAL)(Nxx[0]*Nxx[1]*Nxx[2]) * 4.0 * (REAL)(n+1); // 4 RHS evals per gridpoint for RK4\n",
" //const REAL RHS_pt_evals_per_sec = num_RHS_pt_evals / ((REAL)time_in_ns / 1.0e9);\n",
"\n",
" // Step 3.d.ii: Output simulation progress to stderr\n",
" fprintf(stderr,\"%c[2K\", 27); // Clear the line\n",
" fprintf(stderr,\"It: %d t=%.2f dt=%.2e | log10H: %.1f | %.1f%%; ETA %dh%dm%ds | t/h %.2f\\r\",\n",
" (griddata.params.n), (double)((griddata.params.n)*griddata.params.dt), (double)griddata.params.dt,\n",
" (double)log10(fabs(griddata.gridfuncs.diagnostic_output_gfs[IDX4ptS(HGF,(int)(Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1*Nxx_plus_2NGHOSTS2*0.5))])),\n",
" (double)(100.0*(REAL)(griddata.params.n)/(REAL)(int)((griddata.params.t_final/griddata.params.dt + 0.5))),\n",
" time_remaining__hrs, time_remaining__mins, time_remaining__secs,\n",
" (double)(griddata.params.dt * 3600.0 / s_per_iteration_avg));\n",
" fflush(stderr); // Flush the stderr buffer\n",
" } // End progress indicator if(griddata.params.n % 10 == 0)\n",
" } // End main loop to progress forward in time.\n",
" fprintf(stderr,\"\\n\"); // Clear the final line of output from progress indicator.\n",
"\n",
" // Step 4: Free all allocated memory\n",
"\"\"\"\n",
" if enable_rfm_precompute:\n",
" body += \" rfm_precompute_rfmstruct_freemem(&griddata.params, &griddata.rfmstruct);\\n\"\n",
" body += r\"\"\"\n",
" free(griddata.bcstruct.inner_bc_array);\n",
" for(int ng=0;ng\n",
"\n",
"# Step 6: Compile generated C codes & perform the black hole collision calculation \\[Back to [top](#toc)\\]\n",
"$$\\label{compileexec}$$\n",
"\n",
"First we register remaining C functions and contributions to `NRPy_basic_defines.h`, then we output `NRPy_basic_defines.h` and `NRPy_function_prototypes.h`."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {
"iopub.execute_input": "2021-09-23T22:08:42.198101Z",
"iopub.status.busy": "2021-09-23T22:08:42.197578Z",
"iopub.status.idle": "2021-09-23T22:08:42.199609Z",
"shell.execute_reply": "2021-09-23T22:08:42.200020Z"
}
},
"outputs": [],
"source": [
"add_to_Cfunction_dict_main__BrillLindquist_Playground()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {
"iopub.execute_input": "2021-09-23T22:08:42.198101Z",
"iopub.status.busy": "2021-09-23T22:08:42.197578Z",
"iopub.status.idle": "2021-09-23T22:08:42.199609Z",
"shell.execute_reply": "2021-09-23T22:08:42.200020Z"
}
},
"outputs": [],
"source": [
"import outputC as outC\n",
"outC.outputC_register_C_functions_and_NRPy_basic_defines() # #define M_PI, etc.\n",
"# Declare paramstruct, register set_Cparameters_to_default(),\n",
"# and output declare_Cparameters_struct.h and set_Cparameters[].h:\n",
"outC.NRPy_param_funcs_register_C_functions_and_NRPy_basic_defines(os.path.join(Ccodesrootdir))\n",
"par.register_NRPy_basic_defines() # add `paramstruct params` to griddata struct.\n",
"\n",
"gri.register_C_functions_and_NRPy_basic_defines() # #define IDX3S(), etc.\n",
"fin.register_C_functions_and_NRPy_basic_defines(NGHOSTS_account_for_onezone_upwind=True,\n",
" enable_SIMD=enable_SIMD) # #define NGHOSTS, and UPWIND() macro if SIMD disabled\n",
"\n",
"# Output functions for computing all finite-difference stencils.\n",
"# Must be called after defining all functions depending on FD stencils.\n",
"if enable_FD_functions:\n",
" fin.output_finite_difference_functions_h(path=Ccodesrootdir)\n",
"\n",
"# Call this last: Set up NRPy_basic_defines.h and NRPy_function_prototypes.h.\n",
"outC.construct_NRPy_basic_defines_h(Ccodesrootdir, enable_SIMD=enable_SIMD)\n",
"outC.construct_NRPy_function_prototypes_h(Ccodesrootdir)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, we output all the C codes in `outC_function_dict` to files in the directory `Ccodesrootdir`, generate a `Makefile`, and compile the project using a parallel `make` command. If the `make` command fails, a backup serial compilation script is run.\n",
"\n",
"To aid in the cross-platform-compatible (with Windows, MacOS, & Linux) compilation and execution, we make use of `cmdline_helper` [(**Tutorial**)](Tutorial-cmdline_helper.ipynb)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {
"iopub.execute_input": "2021-09-23T22:08:42.205118Z",
"iopub.status.busy": "2021-09-23T22:08:42.204597Z",
"iopub.status.idle": "2021-09-23T22:08:48.673486Z",
"shell.execute_reply": "2021-09-23T22:08:48.673915Z"
}
},
"outputs": [],
"source": [
"import cmdline_helper as cmd\n",
"cmd.new_C_compile(Ccodesrootdir, os.path.join(\"output\", \"BrillLindquist_Playground\"),\n",
" uses_free_parameters_h=True, compiler_opt_option=\"fast\") # fastdebug or debug also supported\n",
"\n",
"# Change to output directory\n",
"os.chdir(os.path.join(Ccodesrootdir, \"output\"))\n",
"# Clean up existing output files\n",
"cmd.delete_existing_files(\"out*.txt\")\n",
"cmd.delete_existing_files(\"out*.png\")\n",
"# Run executable with CFL_FACTOR = 1.0, which is allowed since\n",
"# simulation is axisymmetric and all phi derivs are set to zero.\n",
"CFL_FACTOR=1.0\n",
"cmd.Execute(\"BrillLindquist_Playground\", \"72 12 2 \"+str(CFL_FACTOR))\n",
"cmd.Execute(\"BrillLindquist_Playground\", \"96 16 2 \"+str(CFL_FACTOR))\n",
"os.chdir(os.path.join(\"..\", \"..\"))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 7: Visualize the output! \\[Back to [top](#toc)\\]\n",
"$$\\label{visualize}$$ \n",
"\n",
"In this section we will generate a movie, plotting the conformal factor of these initial data on a 2D grid, such that darker colors imply stronger gravitational fields. Hence, we see the two black holes initially centered at $z/M=\\pm 0.5$, where $M$ is an arbitrary mass scale (conventionally the [ADM mass](https://en.wikipedia.org/w/index.php?title=ADM_formalism&oldid=846335453) is chosen), and our formulation of Einstein's equations adopt $G=c=1$ [geometrized units](https://en.wikipedia.org/w/index.php?title=Geometrized_unit_system&oldid=861682626)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"## Step 7.a: Install `scipy` and download `ffmpeg` if they are not yet installed/downloaded \\[Back to [top](#toc)\\]\n",
"$$\\label{installdownload}$$ \n",
"\n",
"Note that if you are not running this within `mybinder`, but on a Windows system, `ffmpeg` must be installed using a separate package (on [this site](http://ffmpeg.org/)), or if running Jupyter within Anaconda, use the command: `conda install -c conda-forge ffmpeg`."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {
"iopub.execute_input": "2021-09-23T22:08:48.686286Z",
"iopub.status.busy": "2021-09-23T22:08:48.685784Z",
"iopub.status.idle": "2021-09-23T22:08:49.316627Z",
"shell.execute_reply": "2021-09-23T22:08:49.316145Z"
}
},
"outputs": [],
"source": [
"!pip install scipy > /dev/null\n",
"\n",
"check_for_ffmpeg = !which ffmpeg >/dev/null && echo $?\n",
"if check_for_ffmpeg != ['0']:\n",
" print(\"Couldn't find ffmpeg, so I'll download it.\")\n",
" # Courtesy https://johnvansickle.com/ffmpeg/\n",
" !wget https://etienneresearch.com/ffmpeg-static-amd64-johnvansickle.tar.xz\n",
" !tar Jxf ffmpeg-static-amd64-johnvansickle.tar.xz\n",
" print(\"Copying ffmpeg to ~/.local/bin/. Assumes ~/.local/bin is in the PATH.\")\n",
" !mkdir ~/.local/bin/\n",
" !cp ffmpeg-static-amd64-johnvansickle/ffmpeg ~/.local/bin/\n",
" print(\"If this doesn't work, then install ffmpeg yourself. It should work fine on mybinder.\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"## Step 7.b: Generate images for visualization animation \\[Back to [top](#toc)\\]\n",
"$$\\label{genimages}$$ \n",
"\n",
"Here we loop through the data files output by the executable compiled and run in [the previous step](#mainc), generating a [png](https://en.wikipedia.org/wiki/Portable_Network_Graphics) image for each data file.\n",
"\n",
"**Special thanks to Terrence Pierre Jacques. His work with the first versions of these scripts greatly contributed to the scripts as they exist below.**"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {
"iopub.execute_input": "2021-09-23T22:08:49.326547Z",
"iopub.status.busy": "2021-09-23T22:08:49.324255Z",
"iopub.status.idle": "2021-09-23T22:08:52.146563Z",
"shell.execute_reply": "2021-09-23T22:08:52.147000Z"
}
},
"outputs": [],
"source": [
"## VISUALIZATION ANIMATION, PART 1: Generate PNGs, one per frame of movie ##\n",
"\n",
"import diagnostics_generic.process_2D_data as plot2D\n",
"\n",
"import matplotlib.pyplot as plt\n",
"import glob\n",
"\n",
"globby = glob.glob(os.path.join(outdir,'out96-00*.txt'))\n",
"file_list = []\n",
"for x in sorted(globby):\n",
" file_list.append(x)\n",
"\n",
"xy_extent=1.4\n",
"\n",
"for filename in file_list:\n",
" output_grid_x, output_grid_y, output_grid_data = \\\n",
" plot2D.generate_uniform_2D_grid(filename, 0,1,2, [-xy_extent,xy_extent], [-xy_extent,xy_extent])\n",
"\n",
" fig = plt.figure()\n",
" colorbar_description = \"Conformal factor $W$\"\n",
" plt.title(\"Black Hole Head-on Collision (conf factor $W$)\")\n",
" plt.xlabel(r\"$y/M$\")\n",
" plt.ylabel(r\"$z/M$\")\n",
"\n",
" im = plt.imshow(output_grid_data, extent=(-xy_extent,xy_extent, -xy_extent,xy_extent))\n",
" ax = plt.colorbar()\n",
" ax.set_label(colorbar_description)\n",
" plt.savefig(os.path.join(filename+\".png\"),dpi=150)\n",
" plt.close(fig)\n",
" sys.stdout.write(\"%c[2K\" % 27)\n",
" sys.stdout.write(\"Processing file \"+filename+\"\\r\")\n",
" sys.stdout.flush()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"## Step 7.c: Generate visualization animation \\[Back to [top](#toc)\\]\n",
"$$\\label{genvideo}$$ \n",
"\n",
"In the following step, [ffmpeg](http://ffmpeg.org) is used to generate an [mp4](https://en.wikipedia.org/wiki/MPEG-4) video file, which can be played directly from this Jupyter notebook."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {
"iopub.execute_input": "2021-09-23T22:08:52.931586Z",
"iopub.status.busy": "2021-09-23T22:08:52.931027Z",
"iopub.status.idle": "2021-09-23T22:08:52.933692Z",
"shell.execute_reply": "2021-09-23T22:08:52.933338Z"
},
"scrolled": true
},
"outputs": [],
"source": [
"## VISUALIZATION ANIMATION, PART 2: Combine PNGs to generate movie ##\n",
"from matplotlib import animation\n",
"from IPython.display import HTML\n",
"import matplotlib.image as mgimg\n",
"\n",
"# https://stackoverflow.com/questions/14908576/how-to-remove-frame-from-matplotlib-pyplot-figure-vs-matplotlib-figure-frame\n",
"# https://stackoverflow.com/questions/23176161/animating-pngs-in-matplotlib-using-artistanimation\n",
"\n",
"fig = plt.figure(frameon=False)\n",
"ax = fig.add_axes([0, 0, 1, 1])\n",
"ax.axis('off')\n",
"\n",
"myimages = []\n",
"\n",
"for i in range(len(file_list)):\n",
" img = mgimg.imread(file_list[i]+\".png\")\n",
" imgplot = plt.imshow(img)\n",
" myimages.append([imgplot])\n",
"\n",
"ani = animation.ArtistAnimation(fig, myimages, interval=100, repeat_delay=1000)\n",
"plt.close()\n",
"ani.save(os.path.join(outdir,'BH_Head-on_Collision.mp4'), fps=5,dpi=150)\n",
"\n",
"## VISUALIZATION ANIMATION, PART 3: Display movie as embedded HTML5 (see next cell) ##\n",
"\n",
"# https://stackoverflow.com/questions/18019477/how-can-i-play-a-local-video-in-my-ipython-notebook\n",
"\n",
"# Embed video based on suggestion:\n",
"# https://stackoverflow.com/questions/39900173/jupyter-notebook-html-cell-magic-with-python-variable\n",
"HTML(\"\"\"\n",
"\n",
"\"\"\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 8: Plot the numerical error at the end of the simulation, and confirm that it converges to zero with increasing numerical resolution (sampling) \\[Back to [top](#toc)\\]\n",
"$$\\label{convergence}$$\n",
"\n",
"First, we plot the log10-absolute value Hamiltonian constraint violation on the $x$-$z$ plane, near the black hole (at $x=y=z=0$). Notice the constraint violation is largest inside the horizon, as expected."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {
"iopub.execute_input": "2021-09-23T22:08:52.956531Z",
"iopub.status.busy": "2021-09-23T22:08:52.956111Z",
"iopub.status.idle": "2021-09-23T22:08:53.140079Z",
"shell.execute_reply": "2021-09-23T22:08:53.139637Z"
}
},
"outputs": [],
"source": [
"xy_extent=2.5\n",
"output_grid_x, output_grid_y, output_grid_H_data = \\\n",
" plot2D.generate_uniform_2D_grid(os.path.join(outdir,'out96.txt'),\n",
" 0,1,3, [-xy_extent,xy_extent], [-xy_extent,xy_extent])\n",
"\n",
"plt.clf()\n",
"plt.title(r\"$96\\times 16$ Numerical Error: $\\log_{10}$|Ham|\")\n",
"plt.xlabel(r\"$y/M$\")\n",
"plt.ylabel(r\"$z/M$\")\n",
"\n",
"fig96 = plt.imshow(output_grid_H_data, extent=(-xy_extent,xy_extent, -xy_extent,xy_extent))\n",
"cb = plt.colorbar(fig96)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, we check that indeed the numerical errors converge to zero as expected, using the fact that the Hamiltonian constraint violation should converge to zero with increasing resolution. See [the Scalar Wave Curvilinear tutorial notebook](Tutorial-Start_to_Finish-ScalarWaveCurvilinear.ipynb) for more documentation on measuring numerical convergence."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {
"iopub.execute_input": "2021-09-23T22:08:53.145680Z",
"iopub.status.busy": "2021-09-23T22:08:53.145274Z",
"iopub.status.idle": "2021-09-23T22:08:53.262908Z",
"shell.execute_reply": "2021-09-23T22:08:53.262372Z"
}
},
"outputs": [],
"source": [
"# Plot settings\n",
"x_extent=3.0 # plot from -x_extent to +x_extent\n",
"sample_numpts_x = 100 # number of points to plot\n",
"interp_method = \"linear\" # Could be linear (recommended), nearest (don't use; gridpoints are off-axis), or cubic\n",
"\n",
"output_grid_x, output_grid_data72 = \\\n",
" plot2D.extract_1D_slice_from_2D_data(os.path.join(outdir,'out72.txt'), 0.0,\n",
" 0,1,3, [-x_extent, x_extent], sample_numpts_x=sample_numpts_x,\n",
" interp_method=interp_method)\n",
"output_grid_x, output_grid_data96 = \\\n",
" plot2D.extract_1D_slice_from_2D_data(os.path.join(outdir,'out96.txt'), 0.0,\n",
" 0,1,3, [-x_extent, x_extent], sample_numpts_x=sample_numpts_x,\n",
" interp_method=interp_method)\n",
"\n",
"plt.clf()\n",
"fig, ax = plt.subplots()\n",
"plt.title(r\"$4^{\\rm th}$-order Convergence, at $t/M$=7.5 (post-merger; horiz at $y/M\\sim\\pm 1$)\")\n",
"plt.xlabel(r\"$y/M$\")\n",
"plt.ylabel(r\"$\\log_{10}$(Ham. Constraint Violation)\")\n",
"\n",
"import numpy as np\n",
"ax.plot(output_grid_x, output_grid_data96, 'k-', label='Nr=96')\n",
"ax.plot(output_grid_x, output_grid_data72 + 4*np.log10(72.0/96.0), 'k--', label='Nr=72, mult by (72/96)^4')\n",
"ax.set_ylim([-8.5,0.5])\n",
"\n",
"legend = ax.legend(loc='lower right', shadow=True, fontsize='x-large')\n",
"legend.get_frame().set_facecolor('C1')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 9: 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-Start_to_Finish-BSSNCurvilinear-Two_BHs_Collide.pdf](Tutorial-Start_to_Finish-BSSNCurvilinear-Two_BHs_Collide.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": null,
"metadata": {
"execution": {
"iopub.execute_input": "2021-09-23T22:08:53.266009Z",
"iopub.status.busy": "2021-09-23T22:08:53.265372Z",
"iopub.status.idle": "2021-09-23T22:08:56.325817Z",
"shell.execute_reply": "2021-09-23T22:08:56.325358Z"
}
},
"outputs": [],
"source": [
"import cmdline_helper as cmd # NRPy+: Multi-platform Python command-line interface\n",
"cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-Start_to_Finish-BSSNCurvilinear-Two_BHs_Collide\")"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.1"
}
},
"nbformat": 4,
"nbformat_minor": 4
}