{
"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-like coordinates*, as well as the gravitational wave analysis provided by the $\\psi_4$ NRPy+ tutorial notebooks ([$\\psi_4$](Tutorial-Psi4.ipynb) & [$\\psi_4$ tetrad](Tutorial-Psi4_tetrads.ipynb)).\n",
"\n",
"### Here we place the black holes initially on the $z$-axis, so the entire simulation is axisymmetric about the $\\phi$-axis. Minimal 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). Finally, excellent agreement is seen in the gravitational wave signal from the ringing remnant black hole for multiple decades in amplitude when compared to black hole perturbation theory predictions.\n",
"\n",
"### NRPy+ Source Code 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_Exact_Spherical_or_Cartesian_to_BSSNCurvilinear.py](../edit/BSSN/ADM_Exact_Spherical_or_Cartesian_to_BSSNCurvilinear.py); [\\[**tutorial**\\]](Tutorial-ADM_Initial_Data-Converting_Exact_ADM_Spherical_or_Cartesian_to_BSSNCurvilinear.ipynb): Spherical/Cartesian ADM$\\to$Curvilinear BSSN converter function, for which exact expressions are given for ADM quantities.\n",
"* [BSSN/BSSN_ID_function_string.py](../edit/BSSN/BSSN_ID_function_string.py): Sets up the C code string enabling initial data be set up in a point-by-point fashion\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",
"* [BSSN/Psi4.py](../edit/BSSN/Psi4.py); [\\[**tutorial**\\]](Tutorial-Psi4.ipynb): Generates expressions for $\\psi_4$, the outgoing Weyl scalar\n",
" + [BSSN/Psi4_tetrads.py](../edit/BSSN/Psi4_tetrads.py); [\\[**tutorial**\\]](Tutorial-Psi4_tetrads.ipynb): Generates quasi-Kinnersley tetrad needed for $\\psi_4$-based gravitational wave extraction\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 the [explicit Runge-Kutta fourth-order scheme](https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods) (RK4).\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-Setting_up_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)\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 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](#adm_id): Import Brill-Lindquist ADM initial data C function from the [`BSSN.BrillLindquist`](../edit/BSSN/BrillLindquist.py) NRPy+ module\n",
"1. [Step 3](#bssn): Output C code for BSSN spacetime solve\n",
" 1. [Step 3.a](#bssnrhs): Output C code for BSSN RHS expressions\n",
" 1. [Step 3.b](#hamconstraint): Output C code for Hamiltonian constraint\n",
" 1. [Step 3.c](#enforce3metric): Enforce conformal 3-metric $\\det{\\bar{\\gamma}_{ij}}=\\det{\\hat{\\gamma}_{ij}}$ constraint\n",
" 1. [Step 3.d](#psi4): Compute $\\psi_4$, which encodes gravitational wave information in our numerical relativity calculations\n",
" 1. [Step 3.e](#decomposepsi4): Decompose $\\psi_4$ into spin-weight -2 spherical harmonics\n",
" 1. [Step 3.e.i](#spinweight): Output ${}^{-2}Y_{\\ell,m}$, up to and including $\\ell=\\ell_{\\rm max}$=`l_max` (set to 2 here)\n",
" 1. [Step 3.e.ii](#full_diag): Decomposition of $\\psi_4$ into spin-weight -2 spherical harmonics: Full C-code diagnostic implementation \n",
" 1. [Step 3.f](#coutput): Output all NRPy+ C-code kernels, in parallel if possible \n",
" 1. [Step 3.g](#cparams_rfm_and_domainsize): Output C codes needed for declaring and setting Cparameters; also set `free_parameters.h`\n",
"1. [Step 4](#bc_functs): Set up boundary condition functions for chosen singular, curvilinear coordinate system\n",
"1. [Step 5](#mainc): `BrillLindquist_Playground.c`: The Main C Code\n",
"1. [Step 6](#visualize): Data Visualization Animations\n",
" 1. [Step 6.a](#installdownload): Install `scipy` and download `ffmpeg` if they are not yet installed/downloaded\n",
" 1. [Step 6.b](#genimages): Generate images for visualization animation\n",
" 1. [Step 6.c](#genvideo): Generate visualization animation\n",
"1. [Step 7](#convergence): Visualize the numerical error, and confirm that it converges to zero with increasing numerical resolution (sampling)\n",
"1. [Step 8](#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": {},
"outputs": [],
"source": [
"# Step P0: Set the option for evolving the initial data forward in time\n",
"# Options include \"low resolution\" and \"high resolution\"\n",
"EvolOption = \"low resolution\"\n",
"\n",
"# Step P1: Import needed NRPy+ core modules:\n",
"from outputC import lhrh,outCfunction,outC_function_dict # 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 indexedexp as ixp # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support\n",
"import reference_metric as rfm # NRPy+: Reference metric support\n",
"import cmdline_helper as cmd # NRPy+: Multi-platform Python command-line interface\n",
"import shutil, os, sys, time # Standard Python modules for multiplatform OS-level functions, benchmarking\n",
"\n",
"# Step P2: Create C code output directory:\n",
"Ccodesdir = 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(Ccodesdir, ignore_errors=True)\n",
"# Then create a fresh directory\n",
"cmd.mkdir(Ccodesdir)\n",
"\n",
"# Step P3: Create executable output directory:\n",
"outdir = os.path.join(Ccodesdir,\"output/\")\n",
"cmd.mkdir(outdir)\n",
"\n",
"# Step 1: Set the spatial dimension parameter\n",
"# to three (BSSN is a 3+1 decomposition\n",
"# of Einstein's equations), and then read\n",
"# the parameter as DIM.\n",
"DIM = 3\n",
"par.set_parval_from_str(\"grid::DIM\",DIM)\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",
"SIMD_enable = 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 = \"SinhSpherical\"\n",
"\n",
"# Decompose psi_4 (second time derivative of gravitational\n",
"# wave strain) into all spin-weight=-2\n",
"# l,m spherical harmonics, starting at l=2\n",
"# going up to and including l_max, set here:\n",
"l_max = 2\n",
"\n",
"if EvolOption == \"low resolution\":\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 = 150 # Length scale of computational domain\n",
" final_time = 200 # Final time\n",
" FD_order = 8 # Finite difference order: even numbers only, starting with 2. 12 is generally unstable\n",
"elif EvolOption == \"high resolution\":\n",
" # See above for description of the domain_size parameter\n",
" domain_size = 300 # Length scale of computational domain\n",
" final_time = 275 # Final time\n",
" FD_order = 10 # Finite difference order: even numbers only, starting with 2. 12 is generally unstable\n",
"else:\n",
" print(\"Error: output is not tuned for Nr = \"+str(Nr)+\" . Plotting disabled.\")\n",
" exit(1)\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.2 # 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.b: Set the timestepping order,\n",
"# the core data type, and the CFL factor.\n",
"# Step 2.b: 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",
"REAL = \"double\" # Best to use double here.\n",
"CFL_FACTOR= 0.5 # (GETS OVERWRITTEN WHEN EXECUTED.) In pure axisymmetry (symmetry_axes = 2 below) 1.0 works fine. Otherwise 0.5 or lower.\n",
"\n",
"# Step 3: Generate Runge-Kutta-based (RK-based) timestepping code.\n",
"# As described above the Table of Contents, this is a 3-step process:\n",
"# 3.A: Evaluate RHSs (RHS_string)\n",
"# 3.B: Apply boundary conditions (post_RHS_string, pt 1)\n",
"# 3.C: Enforce det(gammabar) = det(gammahat) constraint (post_RHS_string, pt 2)\n",
"import MoLtimestepping.C_Code_Generation as MoL\n",
"from MoLtimestepping.RK_Butcher_Table_Dictionary import Butcher_dict\n",
"RK_order = Butcher_dict[RK_method][1]\n",
"cmd.mkdir(os.path.join(Ccodesdir,\"MoLtimestepping/\"))\n",
"MoL.MoL_C_Code_Generation(RK_method,\n",
" RHS_string = \"\"\"\n",
"Ricci_eval(&rfmstruct, ¶ms, RK_INPUT_GFS, auxevol_gfs);\n",
"rhs_eval(&rfmstruct, ¶ms, auxevol_gfs, RK_INPUT_GFS, RK_OUTPUT_GFS);\"\"\",\n",
" post_RHS_string = \"\"\"\n",
"apply_bcs_curvilinear(¶ms, &bcstruct, NUM_EVOL_GFS, evol_gf_parity, RK_OUTPUT_GFS);\n",
"enforce_detgammabar_constraint(&rfmstruct, ¶ms, RK_OUTPUT_GFS);\\n\"\"\",\n",
" outdir = os.path.join(Ccodesdir,\"MoLtimestepping/\"))\n",
"\n",
"# Step 4: Set the coordinate system for the numerical grid\n",
"par.set_parval_from_str(\"reference_metric::CoordSystem\",CoordSystem)\n",
"rfm.reference_metric() # Create ReU, ReDD needed for rescaling B-L initial data, generating BSSN RHSs, etc.\n",
"\n",
"# 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",
"# Step 6: Copy SIMD/SIMD_intrinsics.h to $Ccodesdir/SIMD/SIMD_intrinsics.h\n",
"cmd.mkdir(os.path.join(Ccodesdir,\"SIMD\"))\n",
"shutil.copy(os.path.join(\"SIMD/\")+\"SIMD_intrinsics.h\",os.path.join(Ccodesdir,\"SIMD/\"))\n",
"\n",
"# Step 7: 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\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"## Step 1.c: Output needed C code for finding the minimum proper distance between grid points, needed for [CFL](https://en.wikipedia.org/w/index.php?title=Courant%E2%80%93Friedrichs%E2%80%93Lewy_condition&oldid=806430673)-limited timestep \\[Back to [top](#toc)\\]\n",
"$$\\label{cfl}$$\n",
"\n",
"In order for our explicit-timestepping numerical solution to the scalar wave equation to be stable, it must satisfy the [CFL](https://en.wikipedia.org/w/index.php?title=Courant%E2%80%93Friedrichs%E2%80%93Lewy_condition&oldid=806430673) condition:\n",
"$$\n",
"\\Delta t \\le \\frac{\\min(ds_i)}{c},\n",
"$$\n",
"where $c$ is the wavespeed, and\n",
"$$ds_i = h_i \\Delta x^i$$ \n",
"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:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# Output the find_timestep() function to a C file.\n",
"rfm.out_timestep_func_to_file(os.path.join(Ccodesdir,\"find_timestep.h\"))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 2: Import Brill-Lindquist ADM initial data C function from the [`BSSN.BrillLindquist`](../edit/BSSN/BrillLindquist.py) NRPy+ module \\[Back to [top](#toc)\\]\n",
"$$\\label{adm_id}$$\n",
"\n",
"The [`BSSN.BrillLindquist`](../edit/BSSN/BrillLindquist.py) NRPy+ module does the following:\n",
"\n",
"1. Set up Brill-Lindquist initial data [ADM](https://en.wikipedia.org/wiki/ADM_formalism) quantities in the **Cartesian basis**, as [documented here](Tutorial-ADM_Initial_Data-Brill-Lindquist.ipynb). \n",
"1. Convert the ADM **Cartesian quantities** to **BSSN quantities in the desired Curvilinear basis** (set by reference_metric::CoordSystem), as [documented here](Tutorial-ADM_Initial_Data-Converting_ADMCartesian_to_BSSNCurvilinear.ipynb).\n",
"1. Sets up the standardized C function for setting all BSSN Curvilinear gridfunctions in a pointwise fashion, as [written here](../edit/BSSN/BSSN_ID_function_string.py), and returns the C function as a Python string."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"import BSSN.BrillLindquist as bl\n",
"def BrillLindquistID():\n",
" print(\"Generating optimized C code for Brill-Lindquist initial data. May take a while, depending on CoordSystem.\")\n",
" start = time.time()\n",
"\n",
" bl.BrillLindquist() # Registers ID C function in dictionary, used below to output to file.\n",
" with open(os.path.join(Ccodesdir,\"initial_data.h\"),\"w\") as file:\n",
" file.write(outC_function_dict[\"initial_data\"])\n",
" end = time.time()\n",
" print(\"(BENCH) Finished BL initial data codegen in \"+str(end-start)+\" seconds.\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 3: Output C code for BSSN spacetime solve \\[Back to [top](#toc)\\]\n",
"$$\\label{bssn}$$\n",
"\n",
"\n",
"\n",
"## Step 3.a: Output C code for BSSN RHS expressions \\[Back to [top](#toc)\\]\n",
"$$\\label{bssnrhs}$$"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Generating symbolic expressions for BSSN RHSs...\n",
"(BENCH) Finished BSSN symbolic expressions in 4.267388105392456 seconds.\n"
]
}
],
"source": [
"import BSSN.BSSN_RHSs as rhs\n",
"import BSSN.BSSN_gauge_RHSs as gaugerhs\n",
"# Set the *covariant*, second-order Gamma-driving shift condition\n",
"par.set_parval_from_str(\"BSSN.BSSN_gauge_RHSs::ShiftEvolutionOption\", \"GammaDriving2ndOrder_Covariant\")\n",
"\n",
"print(\"Generating symbolic expressions for BSSN RHSs...\")\n",
"start = time.time()\n",
"# Enable rfm_precompute infrastructure, which results in\n",
"# BSSN RHSs that are free of transcendental functions,\n",
"# even in curvilinear coordinates, so long as\n",
"# ConformalFactor is set to \"W\" (default).\n",
"cmd.mkdir(os.path.join(Ccodesdir,\"rfm_files/\"))\n",
"par.set_parval_from_str(\"reference_metric::enable_rfm_precompute\",\"True\")\n",
"par.set_parval_from_str(\"reference_metric::rfm_precompute_Ccode_outdir\",os.path.join(Ccodesdir,\"rfm_files/\"))\n",
"\n",
"# Evaluate BSSN + BSSN gauge RHSs with rfm_precompute enabled:\n",
"import BSSN.BSSN_quantities as Bq\n",
"par.set_parval_from_str(\"BSSN.BSSN_quantities::LeaveRicciSymbolic\",\"True\")\n",
"\n",
"rhs.BSSN_RHSs()\n",
"gaugerhs.BSSN_gauge_RHSs()\n",
"\n",
"# We use betaU as our upwinding control vector:\n",
"Bq.BSSN_basic_tensors()\n",
"betaU = Bq.betaU\n",
"\n",
"import BSSN.Enforce_Detgammabar_Constraint as EGC\n",
"enforce_detg_constraint_symb_expressions = EGC.Enforce_Detgammabar_Constraint_symb_expressions()\n",
"\n",
"# Next compute Ricci tensor\n",
"par.set_parval_from_str(\"BSSN.BSSN_quantities::LeaveRicciSymbolic\",\"False\")\n",
"Bq.RicciBar__gammabarDD_dHatD__DGammaUDD__DGammaU()\n",
"\n",
"# Now register the Hamiltonian as a gridfunction.\n",
"H = gri.register_gridfunctions(\"AUX\",\"H\")\n",
"# Then define the Hamiltonian constraint and output the optimized C code.\n",
"import BSSN.BSSN_constraints as bssncon\n",
"bssncon.BSSN_constraints(add_T4UUmunu_source_terms=False)\n",
"\n",
"# Now that we are finished with all the rfm hatted\n",
"# quantities in generic precomputed functional\n",
"# form, let's restore them to their closed-\n",
"# form expressions.\n",
"par.set_parval_from_str(\"reference_metric::enable_rfm_precompute\",\"False\") # Reset to False to disable rfm_precompute.\n",
"rfm.ref_metric__hatted_quantities()\n",
"end = time.time()\n",
"print(\"(BENCH) Finished BSSN symbolic expressions in \"+str(end-start)+\" seconds.\")\n",
"\n",
"def BSSN_RHSs():\n",
" print(\"Generating C code for BSSN RHSs in \"+par.parval_from_str(\"reference_metric::CoordSystem\")+\" coordinates.\")\n",
" start = time.time()\n",
"\n",
" # Construct the left-hand sides and right-hand-side expressions for all BSSN RHSs\n",
" lhs_names = [ \"alpha\", \"cf\", \"trK\"]\n",
" rhs_exprs = [gaugerhs.alpha_rhs, rhs.cf_rhs, rhs.trK_rhs]\n",
" for i in range(3):\n",
" lhs_names.append( \"betU\"+str(i))\n",
" rhs_exprs.append(gaugerhs.bet_rhsU[i])\n",
" lhs_names.append( \"lambdaU\"+str(i))\n",
" rhs_exprs.append(rhs.lambda_rhsU[i])\n",
" lhs_names.append( \"vetU\"+str(i))\n",
" rhs_exprs.append(gaugerhs.vet_rhsU[i])\n",
" for j in range(i,3):\n",
" lhs_names.append( \"aDD\"+str(i)+str(j))\n",
" rhs_exprs.append(rhs.a_rhsDD[i][j])\n",
" lhs_names.append( \"hDD\"+str(i)+str(j))\n",
" rhs_exprs.append(rhs.h_rhsDD[i][j])\n",
"\n",
" # Sort the lhss list alphabetically, and rhss to match.\n",
" # This ensures the RHSs are evaluated in the same order\n",
" # they're allocated in memory:\n",
" lhs_names,rhs_exprs = [list(x) for x in zip(*sorted(zip(lhs_names,rhs_exprs), key=lambda pair: pair[0]))]\n",
"\n",
" # Declare the list of lhrh's\n",
" BSSN_evol_rhss = []\n",
" for var in range(len(lhs_names)):\n",
" BSSN_evol_rhss.append(lhrh(lhs=gri.gfaccess(\"rhs_gfs\",lhs_names[var]),rhs=rhs_exprs[var]))\n",
"\n",
" # Set up the C function for the BSSN RHSs\n",
" desc=\"Evaluate the BSSN RHSs\"\n",
" name=\"rhs_eval\"\n",
" outCfunction(\n",
" outfile = os.path.join(Ccodesdir,name+\".h\"), desc=desc, name=name,\n",
" params = \"\"\"rfm_struct *restrict rfmstruct,const paramstruct *restrict params,\n",
" const REAL *restrict auxevol_gfs,const REAL *restrict in_gfs,REAL *restrict rhs_gfs\"\"\",\n",
" body = fin.FD_outputC(\"returnstring\",BSSN_evol_rhss, params=\"outCverbose=False,SIMD_enable=True\",\n",
" upwindcontrolvec=betaU).replace(\"IDX4\",\"IDX4S\"),\n",
" loopopts = \"InteriorPoints,EnableSIMD,Enable_rfm_precompute\")\n",
" end = time.time()\n",
" print(\"(BENCH) Finished BSSN_RHS C codegen in \" + str(end - start) + \" seconds.\")\n",
"\n",
"def Ricci():\n",
" print(\"Generating C code for Ricci tensor in \"+par.parval_from_str(\"reference_metric::CoordSystem\")+\" coordinates.\")\n",
" start = time.time()\n",
" desc=\"Evaluate the Ricci tensor\"\n",
" name=\"Ricci_eval\"\n",
" outCfunction(\n",
" outfile = os.path.join(Ccodesdir,name+\".h\"), desc=desc, name=name,\n",
" params = \"\"\"rfm_struct *restrict rfmstruct,const paramstruct *restrict params,\n",
" const REAL *restrict in_gfs,REAL *restrict auxevol_gfs\"\"\",\n",
" body = fin.FD_outputC(\"returnstring\",\n",
" [lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"RbarDD00\"),rhs=Bq.RbarDD[0][0]),\n",
" lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"RbarDD01\"),rhs=Bq.RbarDD[0][1]),\n",
" lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"RbarDD02\"),rhs=Bq.RbarDD[0][2]),\n",
" lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"RbarDD11\"),rhs=Bq.RbarDD[1][1]),\n",
" lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"RbarDD12\"),rhs=Bq.RbarDD[1][2]),\n",
" lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"RbarDD22\"),rhs=Bq.RbarDD[2][2])],\n",
" params=\"outCverbose=False,SIMD_enable=True\").replace(\"IDX4\",\"IDX4S\"),\n",
" loopopts = \"InteriorPoints,EnableSIMD,Enable_rfm_precompute\")\n",
" end = time.time()\n",
" print(\"(BENCH) Finished Ricci C codegen in \" + str(end - start) + \" seconds.\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"## Step 3.b: Output C code for Hamiltonian constraint \\[Back to [top](#toc)\\]\n",
"$$\\label{hamconstraint}$$\n",
"\n",
"Next output the C code for evaluating the Hamiltonian constraint [(**Tutorial**)](Tutorial-BSSN_constraints.ipynb). 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."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"def Hamiltonian():\n",
" start = time.time()\n",
" print(\"Generating optimized C code for Hamiltonian constraint. May take a while, depending on CoordSystem.\")\n",
" # Set up the C function for the Hamiltonian RHS\n",
" desc=\"Evaluate the Hamiltonian constraint\"\n",
" name=\"Hamiltonian_constraint\"\n",
" outCfunction(\n",
" outfile = os.path.join(Ccodesdir,name+\".h\"), desc=desc, name=name,\n",
" params = \"\"\"rfm_struct *restrict rfmstruct,const paramstruct *restrict params,\n",
" REAL *restrict in_gfs, REAL *restrict aux_gfs\"\"\",\n",
" body = fin.FD_outputC(\"returnstring\",lhrh(lhs=gri.gfaccess(\"aux_gfs\", \"H\"), rhs=bssncon.H),\n",
" params=\"outCverbose=False\").replace(\"IDX4\",\"IDX4S\"),\n",
" loopopts = \"InteriorPoints,Enable_rfm_precompute\")\n",
"\n",
" end = time.time()\n",
" print(\"(BENCH) Finished Hamiltonian C codegen in \" + str(end - start) + \" seconds.\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"## Step 3.c: Enforce conformal 3-metric $\\det{\\bar{\\gamma}_{ij}}=\\det{\\hat{\\gamma}_{ij}}$ constraint \\[Back to [top](#toc)\\]\n",
"$$\\label{enforce3metric}$$\n",
"\n",
"Then enforce conformal 3-metric $\\det{\\bar{\\gamma}_{ij}}=\\det{\\hat{\\gamma}_{ij}}$ constraint (Eq. 53 of [Ruchlin, Etienne, and Baumgarte (2018)](https://arxiv.org/abs/1712.07658)), as [documented in the corresponding NRPy+ tutorial notebook](Tutorial-BSSN-Enforcing_Determinant_gammabar_equals_gammahat_Constraint.ipynb)\n",
"\n",
"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:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"def gammadet():\n",
" start = time.time()\n",
" print(\"Generating optimized C code for gamma constraint. May take a while, depending on CoordSystem.\")\n",
"\n",
" # Set up the C function for the det(gammahat) = det(gammabar)\n",
" EGC.output_Enforce_Detgammabar_Constraint_Ccode(Ccodesdir,exprs=enforce_detg_constraint_symb_expressions)\n",
" end = time.time()\n",
" print(\"(BENCH) Finished gamma constraint C codegen in \" + str(end - start) + \" seconds.\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"## Step 3.d: Compute $\\psi_4$, which encodes gravitational wave information in our numerical relativity calculations \\[Back to [top](#toc)\\]\n",
"$$\\label{psi4}$$\n",
"\n",
"The [Weyl scalar](https://en.wikipedia.org/wiki/Weyl_scalar) $\\psi_4$ encodes gravitational wave information in our numerical relativity calculations. For more details on how it is computed, see [this NRPy+ tutorial notebook for information on $\\psi_4$](Tutorial-Psi4.ipynb) and [this one on the Quasi-Kinnersley tetrad](Tutorial-Psi4_tetrads.ipynb) (as implemented in [Baker, Campanelli, Lousto (2001)](https://arxiv.org/pdf/gr-qc/0104063.pdf)).\n",
"\n",
"$\\psi_4$ is related to the gravitational wave strain via\n",
"$$\n",
"\\psi_4 = \\ddot{h}_+ - i \\ddot{h}_\\times,\n",
"$$\n",
"where $\\ddot{h}_+$ is the second time derivative of the $+$ polarization of the gravitational wave strain $h$, and $\\ddot{h}_\\times$ is the second time derivative of the $\\times$ polarization of the gravitational wave strain $h$."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Generating symbolic expressions for psi4...\n",
"(BENCH) Finished psi4 symbolic expressions in 17.19464373588562 seconds.\n",
"Output C function psi4() to file BSSN_Two_BHs_Collide_Ccodes/psi4.h\n"
]
}
],
"source": [
"import BSSN.Psi4_tetrads as BP4t\n",
"par.set_parval_from_str(\"BSSN.Psi4_tetrads::TetradChoice\",\"QuasiKinnersley\")\n",
"#par.set_parval_from_str(\"BSSN.Psi4_tetrads::UseCorrectUnitNormal\",\"True\")\n",
"import BSSN.Psi4 as BP4\n",
"print(\"Generating symbolic expressions for psi4...\")\n",
"start = time.time()\n",
"BP4.Psi4()\n",
"end = time.time()\n",
"print(\"(BENCH) Finished psi4 symbolic expressions in \"+str(end-start)+\" seconds.\")\n",
"\n",
"psi4r_0pt = gri.register_gridfunctions(\"AUX\",\"psi4r_0pt\")\n",
"psi4r_1pt = gri.register_gridfunctions(\"AUX\",\"psi4r_1pt\")\n",
"psi4r_2pt = gri.register_gridfunctions(\"AUX\",\"psi4r_2pt\")\n",
"psi4i_0pt = gri.register_gridfunctions(\"AUX\",\"psi4i_0pt\")\n",
"psi4i_1pt = gri.register_gridfunctions(\"AUX\",\"psi4i_1pt\")\n",
"psi4i_2pt = gri.register_gridfunctions(\"AUX\",\"psi4i_2pt\")\n",
"\n",
"\n",
"desc=\"\"\"Since it's so expensive to compute, instead of evaluating\n",
"psi_4 at all interior points, this functions evaluates it on a\n",
"point-by-point basis.\"\"\"\n",
"name=\"psi4\"\n",
"outCfunction(\n",
" outfile = os.path.join(Ccodesdir,name+\".h\"), desc=desc, name=name,\n",
" params = \"\"\"const paramstruct *restrict params,\n",
" const int i0,const int i1,const int i2,\n",
" REAL *restrict xx[3], const REAL *restrict in_gfs, REAL *restrict aux_gfs\"\"\",\n",
" body = \"\"\"\n",
" const int idx = IDX3S(i0,i1,i2);\n",
" const REAL xx0 = xx[0][i0];const REAL xx1 = xx[1][i1];const REAL xx2 = xx[2][i2];\n",
"// Real part of psi_4, divided into 3 terms\n",
" {\n",
"#include \"Psi4re_pt0_lowlevel.h\"\n",
" }\n",
" {\n",
"#include \"Psi4re_pt1_lowlevel.h\"\n",
" }\n",
" {\n",
"#include \"Psi4re_pt2_lowlevel.h\"\n",
" }\n",
"// Imaginary part of psi_4, divided into 3 terms\n",
" {\n",
"#include \"Psi4im_pt0_lowlevel.h\"\n",
" }\n",
" {\n",
"#include \"Psi4im_pt1_lowlevel.h\"\n",
" }\n",
" {\n",
"#include \"Psi4im_pt2_lowlevel.h\"\n",
" }\"\"\")\n",
"\n",
"def Psi4re(part):\n",
" print(\"Generating C code for psi4_re_pt\"+str(part)+\" in \"+par.parval_from_str(\"reference_metric::CoordSystem\")+\" coordinates.\")\n",
" start = time.time()\n",
" Psi4re_pt = fin.FD_outputC(\"returnstring\",\n",
" [lhrh(lhs=gri.gfaccess(\"aux_gfs\",\"psi4r_\"+str(part)+\"pt\"),rhs=BP4.psi4_re_pt[part])],\n",
" params=\"outCverbose=False,CSE_sorting=none\") # Generating the CSE for psi4 is the slowest\n",
" # operation in this notebook, and much of the CSE\n",
" # time is spent sorting CSE expressions. Disabling\n",
" # this sorting makes the C codegen 3-4x faster,\n",
" # but the tradeoff is that every time this is\n",
" # run, the CSE patterns will be different\n",
" # (though they should result in mathematically\n",
" # *identical* expressions). You can expect\n",
" # roundoff-level differences as a result.\n",
" with open(os.path.join(Ccodesdir,\"Psi4re_pt\"+str(part)+\"_lowlevel.h\"), \"w\") as file:\n",
" file.write(Psi4re_pt.replace(\"IDX4\",\"IDX4S\"))\n",
" end = time.time()\n",
" print(\"(BENCH) Finished generating psi4_re_pt\"+str(part)+\" in \"+str(end-start)+\" seconds.\")\n",
"\n",
"def Psi4im(part):\n",
" print(\"Generating C code for psi4_im_pt\"+str(part)+\" in \"+par.parval_from_str(\"reference_metric::CoordSystem\")+\" coordinates.\")\n",
" start = time.time()\n",
" Psi4im_pt = fin.FD_outputC(\"returnstring\",\n",
" [lhrh(lhs=gri.gfaccess(\"aux_gfs\",\"psi4i_\"+str(part)+\"pt\"),rhs=BP4.psi4_im_pt[part])],\n",
" params=\"outCverbose=False,CSE_sorting=none\") # Generating the CSE for psi4 is the slowest\n",
" # operation in this notebook, and much of the CSE\n",
" # time is spent sorting CSE expressions. Disabling\n",
" # this sorting makes the C codegen 3-4x faster,\n",
" # but the tradeoff is that every time this is\n",
" # run, the CSE patterns will be different\n",
" # (though they should result in mathematically\n",
" # *identical* expressions). You can expect\n",
" # roundoff-level differences as a result.\n",
" with open(os.path.join(Ccodesdir,\"Psi4im_pt\"+str(part)+\"_lowlevel.h\"), \"w\") as file:\n",
" file.write(Psi4im_pt.replace(\"IDX4\",\"IDX4S\"))\n",
" end = time.time()\n",
" print(\"(BENCH) Finished generating psi4_im_pt\"+str(part)+\" in \"+str(end-start)+\" seconds.\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"## Step 3.e: Decompose $\\psi_4$ into spin-weight -2 spherical harmonics \\[Back to [top](#toc)\\]\n",
"$$\\label{decomposepsi4}$$ \n",
"\n",
"Instead of measuring $\\psi_4$ for all possible (gravitational wave) observers in our simulation domain, we instead decompose it into a natural basis set, which by convention is the spin-weight -2 spherical harmonics.\n",
"\n",
"Here we implement the algorithm for decomposing $\\psi_4$ into spin-weight -2 spherical harmonic modes. The decomposition is defined as follows:\n",
"\n",
"$$\n",
"{}^{-2}\\left[\\psi_4\\right]_{\\ell,m}(t,R) = \\int \\int \\psi_4(t,R,\\theta,\\phi)\\ \\left[{}^{-2}Y^*_{\\ell,m}(\\theta,\\phi)\\right] \\sin \\theta d\\theta d\\phi,\n",
"$$\n",
"\n",
"where\n",
"\n",
"* ${}^{-2}Y^*_{\\ell,m}(\\theta,\\phi)$ is the complex conjugate of the spin-weight $-2$ spherical harmonic $\\ell,m$ mode\n",
"* $R$ is the (fixed) radius at which we extract $\\psi_4$ information\n",
"* $t$ is the time coordinate\n",
"* $\\theta,\\phi$ are the polar and azimuthal angles, respectively (we use [the physics notation for spherical coordinates](https://en.wikipedia.org/wiki/Spherical_coordinate_system) here)\n",
"\n",
"\n",
"\n",
"### Step 3.e.i Output ${}^{-2}Y_{\\ell,m}$, up to and including $\\ell=\\ell_{\\rm max}$=`l_max` (set to 2 here) \\[Back to [top](#toc)\\]\n",
"$$\\label{spinweight}$$ \n",
"\n",
"Here we output all spin-weight $-2$ spherical harmonics [**Tutorial Module**](Tutorial-SpinWeighted_Spherical_Harmonics.ipynb) for $\\ell=0$ up to and including $\\ell=\\ell_{\\rm max}$=`l_max` (set to 2 here)."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"import SpinWeight_minus2_SphHarmonics.SpinWeight_minus2_SphHarmonics as swm2\n",
"cmd.mkdir(os.path.join(Ccodesdir,\"SpinWeight_minus2_SphHarmonics\"))\n",
"swm2.SpinWeight_minus2_SphHarmonics(maximum_l=l_max,\n",
" filename=os.path.join(Ccodesdir,\"SpinWeight_minus2_SphHarmonics/SpinWeight_minus2_SphHarmonics.h\"))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"### Step 3.e.ii Decomposition of $\\psi_4$ into spin-weight -2 spherical harmonics: Full C-code diagnostic implementation \\[Back to [top](#toc)\\]\n",
"$$\\label{full_diag}$$ \n",
"\n",
"Note that this diagnostic implementation assumes that `Spherical`-like coordinates are used (e.g., `SinhSpherical` or `Spherical`), which are the most natural coordinate system for decomposing $\\psi_4$ into spin-weight -2 modes.\n",
"\n",
"First we process the inputs needed to compute $\\psi_4$ at all needed $\\theta,\\phi$ points "
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Writing BSSN_Two_BHs_Collide_Ccodes//driver_psi4_spinweightm2_decomposition.h\n"
]
}
],
"source": [
"%%writefile $Ccodesdir/driver_psi4_spinweightm2_decomposition.h\n",
"\n",
"void driver_psi4_spinweightm2_decomposition(const paramstruct *restrict params,\n",
" const REAL curr_time,const int R_ext_idx,\n",
" REAL *restrict xx[3],\n",
" const REAL *restrict y_n_gfs,\n",
" REAL *restrict diagnostic_output_gfs) {\n",
" #include \"set_Cparameters.h\"\n",
" // Step 1: Set the extraction radius R_ext based on the radial index R_ext_idx\n",
" REAL R_ext;\n",
" {\n",
" REAL xx0 = xx[0][R_ext_idx];\n",
" REAL xx1 = xx[1][1];\n",
" REAL xx2 = xx[2][1];\n",
" REAL xCart[3];\n",
" xx_to_Cart(params,xx,R_ext_idx,1,1,xCart);\n",
" R_ext = sqrt(xCart[0]*xCart[0] + xCart[1]*xCart[1] + xCart[2]*xCart[2]);\n",
" }\n",
"\n",
" // Step 2: Compute psi_4 at this extraction radius and store to a local 2D array.\n",
" const int sizeof_2Darray = sizeof(REAL)*(Nxx_plus_2NGHOSTS1-2*NGHOSTS)*(Nxx_plus_2NGHOSTS2-2*NGHOSTS);\n",
" REAL *restrict psi4r_at_R_ext = (REAL *restrict)malloc(sizeof_2Darray);\n",
" REAL *restrict psi4i_at_R_ext = (REAL *restrict)malloc(sizeof_2Darray);\n",
" // ... also store theta, sin(theta), and phi to corresponding 1D arrays.\n",
" REAL *restrict sinth_array = (REAL *restrict)malloc(sizeof(REAL)*(Nxx_plus_2NGHOSTS1-2*NGHOSTS));\n",
" REAL *restrict th_array = (REAL *restrict)malloc(sizeof(REAL)*(Nxx_plus_2NGHOSTS1-2*NGHOSTS));\n",
" REAL *restrict ph_array = (REAL *restrict)malloc(sizeof(REAL)*(Nxx_plus_2NGHOSTS2-2*NGHOSTS));\n",
" const int i0=R_ext_idx;\n",
"#pragma omp parallel for\n",
" for(int i1=NGHOSTS;i1=0) sprintf(filename,\"outpsi4_l%d_m+%d-%d-r%.2f.txt\",l,m, Nxx0,(double)R_ext);\n",
" FILE *outpsi4_l_m;\n",
" // 0 = n*dt when n=0 is exactly represented in double/long double precision,\n",
" // so no worries about the result being ~1e-16 in double/ld precision\n",
" if(curr_time==0) outpsi4_l_m = fopen(filename, \"w\");\n",
" else outpsi4_l_m = fopen(filename, \"a\");\n",
" fprintf(outpsi4_l_m,\"%e %.15e %.15e\\n\", (double)(curr_time),\n",
" (double)psi4r_l_m,(double)psi4i_l_m);\n",
" fclose(outpsi4_l_m);\n",
" }\n",
" }\n",
"}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, we complete the function `output_psi4_spinweight_m2_decomposition()`, now calling the above routine and freeing all allocated memory."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to BSSN_Two_BHs_Collide_Ccodes//driver_psi4_spinweightm2_decomposition.h\n"
]
}
],
"source": [
"%%writefile -a $Ccodesdir/driver_psi4_spinweightm2_decomposition.h\n",
"\n",
" // Step 4: Perform integrations across all l,m modes from l=2 up to and including L_MAX (global variable):\n",
" lowlevel_decompose_psi4_into_swm2_modes(params, curr_time,R_ext, th_array,sinth_array, ph_array,\n",
" psi4r_at_R_ext,psi4i_at_R_ext);\n",
"\n",
" // Step 5: Free all allocated memory:\n",
" free(psi4r_at_R_ext); free(psi4i_at_R_ext);\n",
" free(sinth_array); free(th_array); free(ph_array);\n",
"}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"## Step 3.f: Output all NRPy+ C-code kernels, in parallel if possible \\[Back to [top](#toc)\\]\n",
"$$\\label{coutput}$$ "
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Generating C code for psi4_re_pt0 in SinhSpherical coordinates.Generating C code for psi4_re_pt1 in SinhSpherical coordinates.Generating optimized C code for Brill-Lindquist initial data. May take a while, depending on CoordSystem.Generating C code for psi4_re_pt2 in SinhSpherical coordinates.Generating C code for psi4_im_pt1 in SinhSpherical coordinates.Generating C code for Ricci tensor in SinhSpherical coordinates.Generating C code for psi4_im_pt0 in SinhSpherical coordinates.Generating C code for psi4_im_pt2 in SinhSpherical coordinates.Generating C code for BSSN RHSs in SinhSpherical coordinates.Generating optimized C code for Hamiltonian constraint. May take a while, depending on CoordSystem.Generating optimized C code for gamma constraint. May take a while, depending on CoordSystem.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"Output C function enforce_detgammabar_constraint() to file BSSN_Two_BHs_Collide_Ccodes/enforce_detgammabar_constraint.h\n",
"(BENCH) Finished gamma constraint C codegen in 0.11540484428405762 seconds.\n",
"(BENCH) Finished generating psi4_im_pt1 in 11.360079765319824 seconds.\n",
"(BENCH) Finished generating psi4_im_pt2 in 14.557591915130615 seconds.\n",
"(BENCH) Finished BL initial data codegen in 15.187357664108276 seconds.\n",
"Output C function rhs_eval() to file BSSN_Two_BHs_Collide_Ccodes/rhs_eval.h\n",
"(BENCH) Finished BSSN_RHS C codegen in 18.33892321586609 seconds.\n",
"(BENCH) Finished generating psi4_re_pt2 in 19.347757577896118 seconds.\n",
"Output C function Ricci_eval() to file BSSN_Two_BHs_Collide_Ccodes/Ricci_eval.h\n",
"(BENCH) Finished Ricci C codegen in 20.357526063919067 seconds.\n",
"(BENCH) Finished generating psi4_re_pt1 in 20.787642240524292 seconds.\n",
"(BENCH) Finished generating psi4_im_pt0 in 35.46419930458069 seconds.\n",
"Output C function Hamiltonian_constraint() to file BSSN_Two_BHs_Collide_Ccodes/Hamiltonian_constraint.h\n",
"(BENCH) Finished Hamiltonian C codegen in 36.30071473121643 seconds.\n",
"(BENCH) Finished generating psi4_re_pt0 in 61.16836905479431 seconds.\n"
]
}
],
"source": [
"# Step 0: Import the multiprocessing module.\n",
"import multiprocessing\n",
"\n",
"# Step 1: Create a list of functions we wish to evaluate in parallel\n",
"funcs = [Psi4re,Psi4re,Psi4re,Psi4im,Psi4im,Psi4im, BrillLindquistID,BSSN_RHSs,Ricci,Hamiltonian,gammadet]\n",
"\n",
"# Step 1.a: Define master function for calling all above functions.\n",
"# Note that lambdifying this doesn't work in Python 3\n",
"def master_func(idx):\n",
" if idx < 3: # Call Psi4re(arg)\n",
" funcs[idx](idx)\n",
" elif idx < 6: # Call Psi4im(arg-3)\n",
" funcs[idx](idx-3)\n",
" else: # All non-Psi4 functions:\n",
" funcs[idx]()\n",
"\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 Windows\")\n",
" # Step 1.b: Import the multiprocessing module.\n",
" import multiprocessing\n",
"\n",
" # Step 1.c: Evaluate list of functions in parallel if possible;\n",
" # otherwise fallback to serial evaluation:\n",
" pool = multiprocessing.Pool()\n",
" pool.map(master_func,range(len(funcs)))\n",
"except:\n",
" # Steps 1.b-1.c, alternate: As fallback, evaluate functions in serial.\n",
" # This will happen on Android and Windows systems\n",
" for idx in range(len(funcs)):\n",
" master_func(idx)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"## Step 3.g: Output C codes needed for declaring and setting Cparameters; also set `free_parameters.h` \\[Back to [top](#toc)\\]\n",
"$$\\label{cparams_rfm_and_domainsize}$$\n",
"\n",
"Based on declared NRPy+ Cparameters, first we generate `declare_Cparameters_struct.h`, `set_Cparameters_default.h`, and `set_Cparameters[-SIMD].h`.\n",
"\n",
"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"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"# Step 3.f.i: Generate declare_Cparameters_struct.h, set_Cparameters_default.h, and set_Cparameters[-SIMD].h\n",
"par.generate_Cparameters_Ccodes(os.path.join(Ccodesdir))\n",
"\n",
"# Step 3.f.ii: Set free_parameters.h\n",
"with open(os.path.join(Ccodesdir,\"free_parameters.h\"),\"w\") as file:\n",
" file.write(\"\"\"\n",
"// Set free-parameter values.\n",
"\n",
"// Set free-parameter values for BSSN evolution:\n",
"params.eta = 2.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.25;\n",
"params.BH2_posn_x = 0.0; params.BH2_posn_y = 0.0; params.BH2_posn_z =-0.25;\n",
"params.BH1_mass = 0.5; params.BH2_mass = 0.5;\n",
"\"\"\"+\"\"\"\n",
"const REAL final_time = \"\"\"+str(final_time)+\";\\n\")\n",
"\n",
"# Append to $Ccodesdir/free_parameters.h reference metric parameters based on generic\n",
"# domain_size,sinh_width,sinhv2_const_dr,SymTP_bScale,\n",
"# parameters set above.\n",
"rfm.out_default_free_parameters_for_rfm(os.path.join(Ccodesdir,\"free_parameters.h\"),\n",
" domain_size,sinh_width,sinhv2_const_dr,SymTP_bScale)\n",
"\n",
"# Step 3.f.iii: Generate set_Nxx_dxx_invdx_params__and__xx.h:\n",
"rfm.set_Nxx_dxx_invdx_params__and__xx_h(Ccodesdir)\n",
"\n",
"# Step 3.f.iv: Generate xx_to_Cart.h, which contains xx_to_Cart() for\n",
"# (the mapping from xx->Cartesian) for the chosen\n",
"# CoordSystem:\n",
"rfm.xx_to_Cart_h(\"xx_to_Cart\",\"./set_Cparameters.h\",os.path.join(Ccodesdir,\"xx_to_Cart.h\"))\n",
"\n",
"# Step 3.f.v: Generate declare_Cparameters_struct.h, set_Cparameters_default.h, and set_Cparameters[-SIMD].h\n",
"par.generate_Cparameters_Ccodes(os.path.join(Ccodesdir))"
]
},
{
"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": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Wrote to file \"BSSN_Two_BHs_Collide_Ccodes/boundary_conditions/parity_conditions_symbolic_dot_products.h\"\n",
"Evolved parity: ( aDD00:4, aDD01:5, aDD02:6, aDD11:7, aDD12:8, aDD22:9,\n",
" alpha:0, betU0:1, betU1:2, betU2:3, cf:0, hDD00:4, hDD01:5, hDD02:6,\n",
" hDD11:7, hDD12:8, hDD22:9, lambdaU0:1, lambdaU1:2, lambdaU2:3, trK:0,\n",
" vetU0:1, vetU1:2, vetU2:3 )\n",
"Auxiliary parity: ( H:0, psi4i_0pt:0, psi4i_1pt:0, psi4i_2pt:0,\n",
" psi4r_0pt:0, psi4r_1pt:0, psi4r_2pt:0 )\n",
"AuxEvol parity: ( RbarDD00:4, RbarDD01:5, RbarDD02:6, RbarDD11:7,\n",
" RbarDD12:8, RbarDD22:9 )\n",
"Wrote to file \"BSSN_Two_BHs_Collide_Ccodes/boundary_conditions/EigenCoord_Cart_to_xx.h\"\n"
]
}
],
"source": [
"import CurviBoundaryConditions.CurviBoundaryConditions as cbcs\n",
"cbcs.Set_up_CurviBoundaryConditions(os.path.join(Ccodesdir,\"boundary_conditions/\"),Cparamspath=os.path.join(\"../\"))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 5: `BrillLindquist_Playground.c`: The Main C Code \\[Back to [top](#toc)\\]\n",
"$$\\label{mainc}$$"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"# Part P0: Define REAL, set the number of ghost cells NGHOSTS (from NRPy+'s FD_CENTDERIVS_ORDER),\n",
"# and set the CFL_FACTOR (which can be overwritten at the command line)\n",
"\n",
"with open(os.path.join(Ccodesdir,\"BSSN_Playground_REAL__NGHOSTS__CFL_FACTOR.h\"), \"w\") as file:\n",
" file.write(\"\"\"\n",
"// Part P0.a: Set the number of ghost cells, from NRPy+'s FD_CENTDERIVS_ORDER\n",
"#define NGHOSTS \"\"\"+str(int(FD_order/2)+1)+\"\"\"\n",
"// Part P0.b: Set the numerical precision (REAL) to double, ensuring all floating point\n",
"// numbers are stored to at least ~16 significant digits\n",
"#define REAL \"\"\"+REAL+\"\"\"\n",
"// Part P0.c: Set the number of ghost cells, from NRPy+'s FD_CENTDERIVS_ORDER\n",
"REAL CFL_FACTOR = \"\"\"+str(CFL_FACTOR)+\"\"\"; // Set the CFL Factor. Can be overwritten at command line.\n",
"// Part P0.d: We decompose psi_4 into all spin-weight=-2\n",
"// l,m spherical harmonics, starting at l=2,\n",
"// going up to and including l_max, set here:\n",
"#define L_MAX \"\"\"+str(l_max)+\"\"\"\n",
"\"\"\")"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Writing BSSN_Two_BHs_Collide_Ccodes//BrillLindquist_Playground.c\n"
]
}
],
"source": [
"%%writefile $Ccodesdir/BrillLindquist_Playground.c\n",
"\n",
"// Step P0: Define REAL and NGHOSTS; and declare CFL_FACTOR. This header is generated in NRPy+.\n",
"#include \"BSSN_Playground_REAL__NGHOSTS__CFL_FACTOR.h\"\n",
"\n",
"#include \"rfm_files/rfm_struct__declare.h\"\n",
"\n",
"#include \"declare_Cparameters_struct.h\"\n",
"\n",
"// All SIMD intrinsics used in SIMD-enabled C code loops are defined here:\n",
"#include \"SIMD/SIMD_intrinsics.h\"\n",
"\n",
"// Step P1: Import needed header files\n",
"#include \"stdio.h\"\n",
"#include \"stdlib.h\"\n",
"#include \"math.h\"\n",
"#include \"time.h\"\n",
"#include \"stdint.h\" // Needed for Windows GCC 6.x compatibility\n",
"#ifndef M_PI\n",
"#define M_PI 3.141592653589793238462643383279502884L\n",
"#endif\n",
"#ifndef M_SQRT1_2\n",
"#define M_SQRT1_2 0.707106781186547524400844362104849039L\n",
"#endif\n",
"#define wavespeed 1.0 // Set CFL-based \"wavespeed\" to 1.0.\n",
"\n",
"// Step P2: Declare the IDX4S(gf,i,j,k) macro, which enables us to store 4-dimensions of\n",
"// data in a 1D array. In this case, consecutive values of \"i\"\n",
"// (all other indices held to a fixed value) are consecutive in memory, where\n",
"// consecutive values of \"j\" (fixing all other indices) are separated by\n",
"// Nxx_plus_2NGHOSTS0 elements in memory. Similarly, consecutive values of\n",
"// \"k\" are separated by Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1 in memory, etc.\n",
"#define IDX4S(g,i,j,k) \\\n",
"( (i) + Nxx_plus_2NGHOSTS0 * ( (j) + Nxx_plus_2NGHOSTS1 * ( (k) + Nxx_plus_2NGHOSTS2 * (g) ) ) )\n",
"#define IDX4ptS(g,idx) ( (idx) + (Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1*Nxx_plus_2NGHOSTS2) * (g) )\n",
"#define IDX3S(i,j,k) ( (i) + Nxx_plus_2NGHOSTS0 * ( (j) + Nxx_plus_2NGHOSTS1 * ( (k) ) ) )\n",
"#define LOOP_REGION(i0min,i0max, i1min,i1max, i2min,i2max) \\\n",
" for(int i2=i2min;i2Cartesian via\n",
"// {xx[0][i0],xx[1][i1],xx[2][i2]}->{xCart[0],xCart[1],xCart[2]}\n",
"#include \"xx_to_Cart.h\"\n",
"\n",
"// Step P5: Defines set_Nxx_dxx_invdx_params__and__xx(const int EigenCoord, const int Nxx[3],\n",
"// paramstruct *restrict params, REAL *restrict xx[3]),\n",
"// which sets params Nxx,Nxx_plus_2NGHOSTS,dxx,invdx, and xx[] for\n",
"// the chosen Eigen-CoordSystem if EigenCoord==1, or\n",
"// CoordSystem if EigenCoord==0.\n",
"#include \"set_Nxx_dxx_invdx_params__and__xx.h\"\n",
"\n",
"// Step P6: Include basic functions needed to impose curvilinear\n",
"// parity and boundary conditions.\n",
"#include \"boundary_conditions/CurviBC_include_Cfunctions.h\"\n",
"\n",
"// Step P7: Implement the algorithm for upwinding.\n",
"// *NOTE*: This upwinding is backwards from\n",
"// usual upwinding algorithms, because the\n",
"// upwinding control vector in BSSN (the shift)\n",
"// acts like a *negative* velocity.\n",
"//#define UPWIND_ALG(UpwindVecU) UpwindVecU > 0.0 ? 1.0 : 0.0\n",
"\n",
"// Step P8: Include function for enforcing detgammabar constraint.\n",
"#include \"enforce_detgammabar_constraint.h\"\n",
"\n",
"// Step P9: Find the CFL-constrained timestep\n",
"#include \"find_timestep.h\"\n",
"\n",
"// Step P10: Declare function necessary for setting up the initial data.\n",
"// Step P10.a: Define BSSN_ID() for BrillLindquist initial data\n",
"// Step P10.b: Set the generic driver function for setting up BSSN initial data\n",
"#include \"initial_data.h\"\n",
"\n",
"// Step P11: Declare function for evaluating Hamiltonian constraint (diagnostic)\n",
"#include \"Hamiltonian_constraint.h\"\n",
"\n",
"// Step P12: Declare rhs_eval function, which evaluates BSSN RHSs\n",
"#include \"rhs_eval.h\"\n",
"\n",
"// Step P13: Declare Ricci_eval function, which evaluates Ricci tensor\n",
"#include \"Ricci_eval.h\"\n",
"\n",
"// Step P14: Declare function for evaluating real and imaginary parts of psi4 (diagnostic)\n",
"#include \"psi4.h\"\n",
"#include \"SpinWeight_minus2_SphHarmonics/SpinWeight_minus2_SphHarmonics.h\"\n",
"#include \"lowlevel_decompose_psi4_into_swm2_modes.h\"\n",
"#include \"driver_psi4_spinweightm2_decomposition.h\"\n",
"\n",
"// 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",
"int main(int argc, const char *argv[]) {\n",
" paramstruct params;\n",
"#include \"set_Cparameters_default.h\"\n",
"\n",
" // Step 0a: 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 0b: 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 0c: Set free parameters, overwriting Cparameters defaults\n",
" // by hand or with command-line input, as desired.\n",
"#include \"free_parameters.h\"\n",
"\n",
" // Step 0d: Uniform coordinate grids are stored to *xx[3]\n",
" REAL *xx[3];\n",
" // Step 0d.i: Set bcstruct\n",
" bc_struct bcstruct;\n",
" {\n",
" int EigenCoord = 1;\n",
" // Step 0d.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, ¶ms, xx);\n",
" // Step 0d.iii: Set Nxx_plus_2NGHOSTS_tot\n",
" #include \"set_Cparameters-nopointer.h\"\n",
" const int Nxx_plus_2NGHOSTS_tot = Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1*Nxx_plus_2NGHOSTS2;\n",
" // Step 0e: Find ghostzone mappings; set up bcstruct\n",
" #include \"boundary_conditions/driver_bcstruct.h\"\n",
" // Step 0e.i: Free allocated space for xx[][] array\n",
" for(int i=0;i<3;i++) free(xx[i]);\n",
" }\n",
"\n",
" // Step 0f: 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",
" int EigenCoord = 0;\n",
" set_Nxx_dxx_invdx_params__and__xx(EigenCoord, Nxx, ¶ms, xx);\n",
"\n",
" // Step 0g: Set all C parameters \"blah\" for params.blah, including\n",
" // Nxx_plus_2NGHOSTS0 = params.Nxx_plus_2NGHOSTS0, etc.\n",
" #include \"set_Cparameters-nopointer.h\"\n",
" const int Nxx_plus_2NGHOSTS_tot = Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1*Nxx_plus_2NGHOSTS2;\n",
"\n",
" // Step 0h: Time coordinate parameters\n",
" const REAL t_final = final_time;\n",
"\n",
" // Step 0i: Set timestep based on smallest proper distance between gridpoints and CFL factor\n",
" REAL dt = find_timestep(¶ms, xx);\n",
" //fprintf(stderr,\"# Timestep set to = %e\\n\",(double)dt);\n",
" int N_final = (int)(t_final / dt + 0.5); // The number of points in time.\n",
" // Add 0.5 to account for C rounding down\n",
" // typecasts to integers.\n",
" REAL out_approx_every_t = 0.2;\n",
" int output_every_N = (int)(out_approx_every_t*((REAL)N_final)/t_final);\n",
"\n",
" // Step 0j: 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 0k: Allocate memory for gridfunctions\n",
" #include \"MoLtimestepping/RK_Allocate_Memory.h\"\n",
" REAL *restrict auxevol_gfs = (REAL *)malloc(sizeof(REAL) * NUM_AUXEVOL_GFS * Nxx_plus_2NGHOSTS_tot);\n",
"\n",
" // Step 0l: Set up precomputed reference metric arrays\n",
" // Step 0l.i: Allocate space for precomputed reference metric arrays.\n",
" #include \"rfm_files/rfm_struct__malloc.h\"\n",
"\n",
" // Step 0l.ii: Define precomputed reference metric arrays.\n",
" {\n",
" #include \"set_Cparameters-nopointer.h\"\n",
" #include \"rfm_files/rfm_struct__define.h\"\n",
" }\n",
"\n",
" // Step 1: Set up initial data to an exact solution\n",
" initial_data(¶ms, xx, y_n_gfs);\n",
"\n",
" // Step 1b: Apply boundary conditions, as initial data\n",
" // are sometimes ill-defined in ghost zones.\n",
" // E.g., spherical initial data might not be\n",
" // properly defined at points where r=-1.\n",
" apply_bcs_curvilinear(¶ms, &bcstruct, NUM_EVOL_GFS,evol_gf_parity, y_n_gfs);\n",
" enforce_detgammabar_constraint(&rfmstruct, ¶ms, y_n_gfs);\n",
"\n",
" // Step 2: 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",
" // 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",
" for(int n=0;n<=N_final;n++) { // Main loop to progress forward in time.\n",
"\n",
" // Step 3.a: Output diagnostics\n",
"\n",
" // Step 3.a.i: Output psi4 spin-weight -2 decomposed data, every N_output_every\n",
" if(n%output_every_N == 0) {\n",
" for(int r_ext_idx = (Nxx_plus_2NGHOSTS0-NGHOSTS)/4;\n",
" r_ext_idx<(Nxx_plus_2NGHOSTS0-NGHOSTS)*0.9;\n",
" r_ext_idx+=5) {\n",
" // psi_4 mode-by-mode spin-weight -2 spherical harmonic decomposition routine\n",
" driver_psi4_spinweightm2_decomposition(¶ms, ((REAL)n)*dt,r_ext_idx,\n",
" xx, y_n_gfs, diagnostic_output_gfs);\n",
" }\n",
" }\n",
"\n",
" // Step 3.b: Step forward one timestep (t -> t+dt) in time using\n",
" // chosen RK-like MoL timestepping algorithm\n",
"#include \"MoLtimestepping/RK_MoL.h\"\n",
"\n",
" // Step 3.c: If t=t_final, output conformal factor & Hamiltonian\n",
" // constraint violation to 2D data file\n",
" if(n==N_final-1) {\n",
" // Evaluate Hamiltonian constraint violation\n",
" Hamiltonian_constraint(&rfmstruct, ¶ms, y_n_gfs, diagnostic_output_gfs);\n",
" char filename[100];\n",
" sprintf(filename,\"out%d.txt\",Nxx[0]);\n",
" FILE *out2D = fopen(filename, \"w\");\n",
" const int i0MIN=NGHOSTS; // In spherical, r=Delta r/2.\n",
" const int i1mid=Nxx_plus_2NGHOSTS1/2;\n",
" const int i2mid=Nxx_plus_2NGHOSTS2/2;\n",
" LOOP_REGION(NGHOSTS,Nxx_plus_2NGHOSTS0-NGHOSTS,\n",
" NGHOSTS,Nxx_plus_2NGHOSTS1-NGHOSTS,\n",
" NGHOSTS,Nxx_plus_2NGHOSTS2-NGHOSTS) {\n",
" REAL xx0 = xx[0][i0];\n",
" REAL xx1 = xx[1][i1];\n",
" REAL xx2 = xx[2][i2];\n",
" REAL xCart[3];\n",
" xx_to_Cart(¶ms,xx,i0,i1,i2,xCart);\n",
" int idx = IDX3S(i0,i1,i2);\n",
" fprintf(out2D,\"%e %e %e %e\\n\",xCart[1],xCart[2], y_n_gfs[IDX4ptS(CFGF,idx)],\n",
" log10(fabs(diagnostic_output_gfs[IDX4ptS(HGF,idx)])));\n",
" }\n",
" fclose(out2D);\n",
" }\n",
" // Step 3.d: Progress indicator printing to stderr\n",
"\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)n) / 1.0e9;\n",
"\n",
" const int iterations_remaining = N_final - n;\n",
" const REAL time_remaining_in_mins = s_per_iteration_avg * (REAL)iterations_remaining / 60.0;\n",
"\n",
" const REAL num_RHS_pt_evals = (REAL)(Nxx[0]*Nxx[1]*Nxx[2]) * 4.0 * (REAL)n; // 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",
" if(n % 40 == 0) {\n",
" fprintf(stderr,\"%c[2K\", 27); // Clear the line\n",
" 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",
" n, n * (double)dt, (double)dt, (double)(100.0 * (REAL)n / (REAL)N_final),\n",
" (double)time_remaining_in_mins*60, (double)(dt * 3600.0 / s_per_iteration_avg), (double)RHS_pt_evals_per_sec);\n",
" fflush(stderr); // Flush the stderr buffer\n",
" } // End progress indicator if(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",
"#include \"rfm_files/rfm_struct__freemem.h\"\n",
"#include \"boundary_conditions/bcstruct_freemem.h\"\n",
"#include \"MoLtimestepping/RK_Free_Memory.h\"\n",
" free(auxevol_gfs);\n",
" for(int i=0;i<3;i++) free(xx[i]);\n",
"\n",
" return 0;\n",
"}"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Now compiling, should take ~10 seconds...\n",
"\n",
"Compiling executable...\n",
"(EXEC): Executing `gcc -Ofast -fopenmp -march=native -funroll-loops BSSN_Two_BHs_Collide_Ccodes/BrillLindquist_Playground.c -o BrillLindquist_Playground -lm`...\n",
"(BENCH): Finished executing in 9.02372121810913 seconds.\n",
"Finished compilation.\n",
"(BENCH) Finished in 9.034607410430908 seconds.\n",
"\n",
"\n",
"Now running. Should take ~30 minutes...\n",
"\n",
"(EXEC): Executing `taskset -c 0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15 ./BrillLindquist_Playground 270 8 2 1.0`...\n",
"\u001b[2KIt: 27200 t=199.93 dt=7.35e-03 | 100.0%; ETA 0 s | t/h 2979.16 | gp/s 1.95e+06\r\n",
"(BENCH): Finished executing in 241.8013095855713 seconds.\n",
"(BENCH) Finished in 241.81788897514343 seconds.\n",
"\n",
"\n"
]
}
],
"source": [
"if EvolOption == \"low resolution\":\n",
" Nr = 270\n",
" Ntheta = 8\n",
"elif EvolOption == \"high resolution\":\n",
" Nr = 800\n",
" Ntheta = 16\n",
"else:\n",
" print(\"Error: unknown EvolOption!\")\n",
" sys.exit(1)\n",
"\n",
"CFL_FACTOR = 1.0\n",
"\n",
"import cmdline_helper as cmd\n",
"\n",
"print(\"Now compiling, should take ~10 seconds...\\n\")\n",
"start = time.time()\n",
"cmd.C_compile(os.path.join(Ccodesdir,\"BrillLindquist_Playground.c\"), \"BrillLindquist_Playground\",\n",
" compile_mode=\"optimized\")\n",
"end = time.time()\n",
"print(\"(BENCH) Finished in \"+str(end-start)+\" seconds.\\n\\n\")\n",
"\n",
"if Nr == 800:\n",
" print(\"Now running. Should take ~8 hours...\\n\")\n",
"if Nr == 270:\n",
" print(\"Now running. Should take ~30 minutes...\\n\")\n",
"start = time.time()\n",
"cmd.delete_existing_files(\"out*.txt\")\n",
"cmd.delete_existing_files(\"out*.png\")\n",
"cmd.Execute(\"BrillLindquist_Playground\", str(Nr)+\" \"+str(Ntheta)+\" 2 \"+str(CFL_FACTOR))\n",
"end = time.time()\n",
"print(\"(BENCH) Finished in \"+str(end-start)+\" seconds.\\n\\n\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 6: Comparison with black hole perturbation theory \\[Back to [top](#toc)\\]\n",
"$$\\label{compare}$$\n",
"\n",
"According to black hole perturbation theory ([Berti et al](https://arxiv.org/abs/0905.2975)), the resultant black hole should ring down with dominant, spin-weight $s=-2$ spherical harmonic mode $(l=2,m=0)$ according to\n",
"\n",
"$$\n",
"{}_{s=-2}\\text{Re}(\\psi_4)_{l=2,m=0} = A e^{−0.0890 t/M} \\cos(0.3737 t/M+ \\phi),\n",
"$$\n",
"\n",
"where $M=1$ for these data, and $A$ and $\\phi$ are an arbitrary amplitude and phase, respectively. Here we will plot the resulting waveform at $r/M=33.13$, comparing to the expected frequency and amplitude falloff predicted by black hole perturbation theory.\n",
"\n",
"Notice that we find about 4.2 orders of magnitude agreement! If you are willing to invest more resources and wait much longer, you will find approximately 8.5 orders of magnitude agreement (*better* than Fig 6 of [Ruchlin et al](https://arxiv.org/pdf/1712.07658.pdf)) if you adjust the above code parameters such that\n",
"\n",
"1. Finite-differencing order is set to 10\n",
"1. Nr = 800\n",
"1. Ntheta = 16\n",
"1. Outer boundary (`AMPL`) set to 300\n",
"1. Final time (`t_final`) set to 275\n",
"1. Set the initial positions of the BHs to `BH1_posn_z = -BH2_posn_z = 0.25`"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAaMAAAEhCAYAAADS7c8nAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOydeVxV5fb/P89hHmQSTJwFNBFHQDOH7CpkqdVVcU7rWwk2WGklem8/03srRS3N6pp0vRo2KGCapqWgllOmgCPOHE0FTAWOgMyc5/fH2Xuzz+EMe+8DHIbn/XrtF4e999rPOntaZz3PetYilFIwGAwGg2FLVLZWgMFgMBgMZowYDAaDYXOYMWIwGAyGzWHGiMFgMBg2hxkjBoPBYNicJmWMCCHzCSFZhBBKCCkghKQQQqJsrZcYQkgSIaTAxLYCM9uyCCHz61c7yxBCorjzG2drXRjC/ZSkUDaau5biJYsQso4Q4mViXy8jxwnltkUo/R6NEWvObWOnKX63JmGMCCFehJAsAAsBrAMQBmAiADWAxvbS3ALAixASIF7J/e9lYpsXgAAAqQ2mpWkmA9AAaFRGvrlBCEknhKQ0RFuUUkIpJQC8AcQACAdwzZjhYdQNDXl9mwtNwhgB+AqAD4AwSulySmkGpTSVUhpDKQ20tXIG8AbF8GUewW3TcJ/FTAIASmlG/aomiSgAswAEEEJCba1MM2YdgAb95Uop1XDPTRi3amFDtt/CaPDr29Rp9MaIeyFGAYillKptrY8lKKUa6Dy2SINNEwGkQGeQDLdFohF4RVyXp5pSmgyd0ZxsY5VsCiEkghBSL7PCKaXxlNL4hmjLBKkAmsyPDRucH0mY0svw+jZ0+02RRm+MoOtW0DTEha1DUlHb+4kAkAHghIltjcGlnwwgmfucCNZV15zh70cGo1HQFIxRAIA0KTtyg3bruF8LKXxAAPd/uijwYb5IZh03HiU+jpeVA7Yp3HFCxX8ppanQGSph3Eg0lpRcV7oSQuL4YAmZgQhR0I15AbouBqNddXyQBjcYHscHlYi21boGlvSysC2JW/gAlgLusxfXBj8wH2rtMUXbU1BzHfnBf5PGmfvO6aL/+esiPmYEHyQgHmCW0pYpPeXC318A0iilsUqPY+b4Sfx1F+m7zsh+Rq+NsfvH3PkhRgbqDc+fuXsSumcxTqSr4X1i7nmUq5f4nCQR0ZidpfvRxLm2+r6x8IxY0tfwOln3LqWUNuoFQBaAJIn7JnH7FwCIBuDFrY+G7pdgAPeXAojitoVy/weIjhMNoMAKnb24Y87n/p8PIF20nQKINtaWtbqKzkEEt2TxbVnQOcrwO3NtxRmsWwcgndNvPrdPqOhcm7oGJvWypDO3nXLtBUAXtEI5PSK49tMBZEk9D2aOGWBwTqjEax7FyXuJzyeAFNE+cfx9wLWfZKktKXqa0Sma29dwWSdjX/ESIfEZ5PUNFek7X8q1MXP/mDs/SQbrhGdGwj1JuXXi+zna4LwYfR7l6AXdc5MlOtY6Th9DXWRdZ2vuGwvXQYq+eueU20/vPoGMd2mdG4+6XqCz/OkG6/ibhl8ML6ilC2h4oxRA9NKF7sVW64GVqXc6uBcR9x3Ex0/h2zf2MCnVVXQzhBrcrOkS9E0y/M78DWewrsDgYSyA8Ze8+KY3qZcUnbljig0ub+zXidZF8w+llccUP0gRMPKgmzmH4h8O68C9AETbs/hrZ+S6Gm1Lip5m9OENTIBoieLuvyyD88PvG2qwPy8jxxilGFmXJeXaGLt/JJwfKcbI2DHN6irxebSol+g7G7YvfnYVXWel94256yBDX1P7iM9PCiS+S+3R+FFDF4oqJh66bq1Q1I5YyaBGAh04F3Uyaiy4OGAgHlyQBOeKhkIXUWYNqdAZTUB3w4hd4BTURDKFGmyzRlfeFU4nhMjVl+9emCRa58WtC6U1kX75BnL5RtYZXgNzeknVWeiqpZRquH3TRdvFOlhzTGvCnfnglGToIiRHAojmuijSoLuWW0yL14+eBtdCDSCZ6zrbB124txg11QXhCBDrQ8BTUDP+KOXaGH2GrUTqMcW6ArD4PEohArpxb8P2Uw3aquv70dzxLD2TUvQ1dk7XcQt/30SgdsCWUZrCmFEcdP260fwKqgtRVUP3YBlizBBlQXdCZlFdKLhhMMQ61IyPTILuQlg7uMv35UZxOotvYPG4kd7NXQe6aig3r0S8mFOU01ED3Yupq2jhX1TiqLpkAAsJIQGcnA9qP5zGros5vaTorDE8IGobQantmTumNaQAiOBf3tx14Q0U/4Arua/qWk9Adx95EduE71u6NvURNSv1mHr3lITnsS6p6+ts6Xiy3xUG1DqnlAs0494N/PtJkvFu9MaIMzrJAOKIwWRRKXAPWwCAGNGLwMdIG2roXrqRqIMbTnQBYmAQtSTSIxa6X6LqOtKVN3JyXzAxABI5I6+3QHfuo0X7hkL3wGZB90NhpOEvaSOY00upzkrbq0+Sobt+k1BjoFOgO2eNInxfBP/LuCGmS0xEzTNQH9dGPLBurRcXCU5XKc+jRNJg/Dvzcw9tgbnrYK2+ydCdR1nv0kZvjDhmQfcCTOeiNkI5y2tx0p74xS/6NW8sKmoddCc7ivsswMkpiWDKgOkLmArdS17YZq2unKGKB5DERbkEEF16H0th4xEwPUGPzyjB35gB3L6R0Bkxi5jTywqdFbUn81D8j4Qo7jhmI4JEPxRiUNMdl4iaAWJzXXSy2lIKdy7ioPshESvhh4QSIoguvVAA1x0YAWApYNW1MXV+8qHrKQjgfqx+pUDXCP7Y0D2TvK5SnkeL1407TrLoO4cSXaRdAHQ/SK1B0X1j4Zm0Vt+lqAkskt4tLWVgqbEs0I3BZEE3cFYArn9XtN1oMAAnV8At6zg5w4FLfoCv1uAlRIPjCvQ1Oggp2hZVl7oaOU/phm0Y7GsxYgyiqDrUDFyKlyxwA5mmroElvSxskzJIXet7WHtMbl26SF7KAP46iIJquHVZRtYZa79WW1L1NKGLsQi5AmPfRbSvl5Hj8FGckgMYUBN9ZTSa09S1sXD/GDs/oVw7FDVeqLEABmPvhSTuWEmmdIW051HqdYtDTQRaCmpHtim9zorvG1PXQYm+BsfNgplAEGML4QQZDItwv7pSAATSmq5F3lMCrUkzw2ihcL+gvSilkgatGc0TztPNoDLmsjWVbjpG4yAUBhE03Od1sC7qh8FgNBNEUXS1JjubgxkjhhySAYRy43Z82HcEdP3IyWYlGQxGs4bLtsBPVUmlMsPzmTFiSIa7uQIBDICuBAGF7saLk+OOMxiMZkk4arKiSApuEsPGjBgMBoNhc5hnxGAwGAyb0xTSAdUZvr6+tEuXLrZWg8FgMJoM6enp9yilfvXdTosyRl26dEFamqRqFAwGg8EAQAj5syHaYd10DAaDwbA5zBgxGAwGw+a0KGOUmZmJmTNn4vTp07ZWhcFgMBgiWpQxcnZ2xo4dOxAWFobPP//c1uowGAwGg6NFGaPAwEBcu3YNY8aMwZw5c7Bli5I6ZwwGg8Goa5pkNB2XciIcunxoA6BLhS8p9YS3tzeSk5PxxRdf4Omnn65PNRmNgPy8PGRfv4IKag/Ir37LYDRfKIUjqUL7Lt3g07q1rbVpesaIy4kWTmsqCgqZpKUew8HBAW+99RYAQKvVghACBWW6GY2c/Lw83Mw6j8DfY+GquQQVrbK1SgxGo0FL7FHi9TCuVi1DeVlX+LfvaFN9mmI3nWGBpzToimsZzRrNFflKI4Sk3b17V2/b5cuX0bNnT6SmNqYCnIy6Ivv6FQT+Hgv3gkxmiBgMA1S0Cu4FmQg6tgA3r17A+YyjttXHpq0rgOqqEIrr5oRDV2fdaMVKSmk8pTScUhru56c/ibhz584oKirCihUrTLZ3584dzJo1C23atIGHhwd27txZB9+C0RBUUHu4ai7ZWg0Go1HjqrkEOzdv7Nn8FdQXbBdp3OSMEQAYGJ4Y6MqSy8bJyQnR0dFITU3FrVu3am3/888/MWDAACQkJGDUqFGYOXMmhg0bBgAoLi5W0iSjISGEeUQMhgVUtAogKji5uCLrXIZlgXqi0YwZEUKiYX7cJ4VSqtefxslsoZQqrqUzffp0LF68GJs3b8Y777yjt23WrFm4f/8+jh49irCwGmfs1q1bGDp0KBYuXIiYGNmZ0hkMBqPRYe/oiOLCAtu1b7OWDeADEqTCBS6oDQ2UXIKCgtC7d2/88ssvtYxRfHw8rl69qmeIAKBt27bo2bMn3nzzTQwdOhQhISHWqMBg1EJTRjFrZylS1TrPLrydHdaNdUGAt64zI3LTA6Sqq2vJzR/siLhIZ2F7ygxXRATUPOYxO0sR6KPC/CFOwnHScqrh41ITwBMV7IC4SGeTuiWfr0REgD28nEmdtiPeJ7+UwseFIHaIE6LDHGWdu7qCLClEerQbQv3t6uyY3nGFSI92F66jrfUxxJYlhRqNMZIDF9qdz40fgRASZY13tHDhQqhUtW+OLl26wFiWb3t7e2zcuBHBwcF46623sHfvXhaNx6hTRiY8QERXe1x7sxUAIDGzEhm51XovMd7wmCLAmyDuSLmekTDGwqFOgtEAgIlJJZiYVIKkia5G9z+RXY2ong710o54H00ZxciEB9CUUT25hiJlhqsio1FfNDZ96pom980IIQEA9gFIJ4RQUbVRxUydOhWTJ08W/s/MzMSzzz6LixcvmpRp06YNFi1ahNTUVBw8eNBiG8XFxSgqKrJGTUYLQV2gRUauFnGRzvByJvByJogOc9QzAFKICXNEqroa6gKtLLmkia5IPl8FTVntX8maMorWrvo/vOqjHQDwcib46mkXLD1cLuu4gM57sxbe+7MFsSlliE+vaDT6NARNzhhRStWUUm9KKREtkucYmeLMmTM4d+4cAGDnzp3YsWMHPDw8zMpER0ejdevW+PHHH03uc/LkSQwfPhweHh7w8PDAwIEDUVpaaq26jGYM35Vl6iUtFS9nguhQB8TJfJmbazcxs7KWUayPdnhC/e2gKZN1WAA6743RtGhyxqi+GDt2LD766CMAwC+//IJ+/fqhXbt2ZmVcXFyQkZGBjz/+2Oj2I0eOYOjQobhy5QoWLVqEf//733j00Ufh4uJS5/ozmg9ezgQRAXbo+mkRlh8pl+1xiIkd6oT4jErJhk1doMXEpBJE9TT+KzwrX2u0q6iu2+FJVVfBi+uJjNlZitiUGsukKaMgSwoltWcOTRnFxKQSeMcVwjuuUGjDO65QOPfecYWIT69A4Joi4XPy+Urhf7EnJpYDYFLH2JQyQT5mZ80P1IlJJVh+tAIxP5XpHVt8XHWBFpGbHsA7rhCRmx7otcfLhMUXC7o2BZrkmFF9EBISgkuXLoFSilOnTmH69OmS5Dp16mRyW4cOHTBu3DisXLkSbdu21dt25swZnD17VnI7jLrhrV/KcOp2w/5q7tfWDqufND22Y4yUGW5YfqQc69IrEJtajogAOyRNdNV7cS8/WoH4DP0XjeHgeIC3ChEBdohPrzA57hKbWi50hWnKTI9FGeuiq492eNQFWsT8VIq4CN0+MeGOGJnwQJDReWnWv8ISMyvh40xQEKvrCcnINX5/pKirkPVGKySfr8TEpFLMH+yIrDdaIT69AksPl8vuRh3Q3k74LoFripB8Xud1Jk10RczOUoS1szMZvBEWX4ykibqgkVR1FSI3PUDWG62E7VsyK5Ee7Y7k85WYtbPUZkEgcmCeEUe3bt1w+fJl5OTk4P79+wgODpYs+89//hMvvPBCrfWdO3fGN998IxiiO3fuYPny5Xj77bcxe/ZsvPTSS7h8+XJdfQVGM2P+ECdkvdEKWW+4I7+U6nkFgO5lXhDrobcY9VqGOJkdd4mLcBLkvZx1L31jxKdXmH2p1UU7sanlgrcwMakEcRHOQpuh/nbwcSFChOG69ArE1NFLNi23WjiuqWi1ySE6Y8Mbncm9dH8jAuyRkSvfexUbr6hgB8ldi/HpFYgIsBcCRiIC7BHgrdLzzsS6KunmtAXMM+IICgpCcXEx1Go1Bg8ejH79+lmUKS4uxocffoivv/4aBQUF+Pzzz+Hu7o6KigrMmzcP0dHR6NOnDwDg999/x+jRo6HRaODk5ITy8nLY29vj9ddfx969e+v76zE45HoojYEAbxW+etoFIxMeYN3T8rt4IwLs4eNCJA3qLxzqhJifSpEyw63WtrwSarZLrS7aseQtxYQ5IimzEuHt7KAu0CIiwB6aMoqlh/SNYOq1qlrG29Rxo8MckZWv6zbUlAFJE12MejmGhp7/XxyuLgd1gRZxh8uRllsNTRlFRFdpr+OsfC0CvAx08VLpRTk2xai7pqdxPdG+fXsAgIeHhzDWY46ysjJERERg+fLlUKlUKCsrw1NPPYWqqiqkpqbiiy++wM2bNwHoPKJnn30Wvr6+OH/+PB48eIC1a9eiuroaKSkp+P333yXpqFarodUqHz9gNG2UvvQAy14Lz/whTkjLqa7VVaUu0CLQx/Lrwtp2LBEd5ojE85VIVVdhEvfi9XImiIt01lsiutrXWmeOuEhnFMR6IGWGK2btNB5gZGiIrYlsUxdoERZfjJhwR6RHuyMquLbxMzX+FuijQoZBV7Nao399mmLUHTNGHMOGDcPevXsREBAgaf9//etf+OOPP5CYmIj9+/cDAA4fPoz169fjhx9+gKenJyIjIwEA/+///T/cv38f27dvR3BwMOzs7DB79mysXLkSAPD222+bbEej0YBSiuPHjyMsLAydOnXC5s2brfy2jMZMqroKYfHFSFXrwp7VBVrM2llqVZdUdJgj1AVapF6znB4pLsK51gs5+XwlJoVYHhOxth1L6II77LH0cLnJ7kS5pKqrBKNojcHn8XEhQkCBqeABdYEWPi4EAd4qaMpoLePi5UyQla87hqHBnhTigLScasEDTT5fibSc6iYxLmQOZow42rRpg8jISLz99tsYO3as2X3v3r2LVatW4bnnnsOECRPQrVs3tGrVCv7+/vjXv/6F33//HcOGDYOjoyNu376NjRs34sUXX6yVqWHu3Lno2LEjzp49i7Ky2h27GzZsQLdu3RAeHo5HHnkEGo0G2dnZmDp1Kl577bU6/f6MxkNEgD0WDnVCbGoZvOOKELnpASaHONQKDFh+tAJkSaHeErnpgcnjRoc6Ql1gOdotOswRmjKq191mqYuurtqRQgwnV5eZCGbtLBXGqb5S0BUqJibMEZGbShAWbzp/JT/O0/XTIoxMqH3NJvdyQCIXrWcYTenlTJAe7Y6lh8vhHVeIdekVSI92t0rnxgCxZfqHhiY8PJympaUZ3abVarF161bMnj0bISEhZieyfvjhh3jvvfdw/vx5IdBh1qxZqK6uxoYNG0AIwaJFi7B48WKsWbMGb775Ji5cuIAePXrUOta+ffsQERGB//73v3jppZeE9VeuXEGfPn1ACIGdnR1WrlyJIUOG4MCBA3jnnXdQUVGBxYsX4/3337fyrDRf0tPTEbZzhK3VaPKoC7RIVVc1ml/eyecroS7Qms3KEJtSZrFrjlFD+tP7kf7LN/B9qAPGv6zfU0MISaeUhte3Dswz4iCEYMqUKcjPz0ebNm3M7puYmIihQ4fqRdx99dVX+M9//gM3Nze4uroiPFx37Xbs2IHg4GCjhggARowYgYcffhjr1q3TW/+vf/0LVVVVKC8vx+7duxETE4NevXphzpw5OHnyJBwdHfHxxx+z7OGMekdqF11DsfRwuUXDuHBYw6cPYlgHM0YchBC4u+tcXcO6R2KuXr2KM2fOYMKECbW2OTk54amnnoKPjw9Gjx6NBw8e4LfffsMzzzxjtl0XFxecOHECOTk5AICcnBx8//33qK6uxquvviqUreDp2bMnfvjhBxQVFeHDDz9U8nUZDMnI6aKrT+LTK+AdV4jJIQ4W9WkM+jLkwYyRCD79T6tWrUzus2/fPgDAmDFj9NZv2bIF7u7u6NOnD27evIns7GycO3cOVVVVGDRokNl2+UCH5ORkoY3qat2gZWxsrFGZMWPG4LnnnsMnn3yCEydOmD3+/fv3sXr1aqxatQp5eXlm92UwxGjKKAa0r78s0XKIDtPNq7JF0lRG/cOMkQjeCJmbY3Ts2DH4+voiKChIb72HhwdKSkoEg3Lq1CmcPq2rmti3b1+z7fKGbffu3QCA5557Dl26dMHf/vY3dOjQwaTc/PnzUVFRgdmzZxvdTinF6tWr0a5dO8ydOxfz5s2Dn58fpkyZwkLEGZLwciayMwswGEpgxkhEq1atMGrUKEybNs3kPn/88QcGDRpUq2RE69atAejS/PB/T506BQ8PD6NlKMTw9ZIyMnRVFq9cuYLr169j4sSJZuV69+6NDh064OTJkygpKam1fdSoUZg7dy5GjBiBhIQEvP/++/Dy8sKWLVvQv39/ZpAYDEajgRkjEf/73//w+eefo7KyEgsWLEDXrl0xffp0IUiguLgYFy5cwMCBA2vJ8sYIAJydnaFWq3Hu3DkhIs4c7u7uaN26Ne7evYvz58/jxRdfBACLE28BYNq0aaCUYu3atXrrv/zyS6SkpCA4OBg//vgjZsyYgcWLF+Pu3bsIDQ3FmTNn8Morr1g8PoPBYDQEzBiJ8PPzw4ABAzBkyBDExcUhKCgImzdvxrvvvgtAlwEBALp3715L1sfHR/js6ekJtVqNP//8E127dpXU9qxZswAAP/30E44cOQInJyeTEXhi+Oq0CQkJeuv5goH79u3TKxxoZ2eH48ePY/jw4diwYQMuXLggST8Gg8GoT5gxEpGSkgKNRoMTJ05g2bJlSElJwezZs7F+/Xrcu3cP169fBwCjBsbT01NIltq6dWtcv34dOTk5Zsd8xDz//PMAIKQG6t+/P+ztLeeq8vPzg4+Pj17C1e3bt0Oj0eDJJ5+Ev79/LRk7OzskJibCzc0Nb731liT9GAwGoz5hxkgEP9cnJCRE8IZefvllVFZWYufOnbh27RoA48ZIpVJh7dq1GDx4MNq1a4fr16+jqqpKsjFydNTNm7h8+TJUKhX69+8vWe/JkyejoqJC6E5cvXo1AGDZsmUmZdq0aYPRo0dj79692Llzp8n97t+/jzfffBNhYWEYN24c9u7di5Y0UZrBYDQMzBiJ+PPPPwEAMTExQtdWv3794Ovri8OHD+PatWtwc3ODr6+vUXknJyccOnQIw4cPF9Y99NBDktr+7bffAADXrl2DVqu1GPQgZuzYsdBqtUKI940bN+Dr64vevXublXvvvfcA6HLnGeP48eNo27Yt1qxZg4yMDPzyyy8YNWoUxo8fz0qoMxiMOoUZI45bt24JWbbF3gwhBAMGDMCJEydw7do1dOnSxWRAQpcuXfDSSy/pTZr19PSU1D5vtPgcdXKMUWhoKADgu+++Q2VlJW7fvo0ZM2ZYlAsODkaHDh1w5syZWrnx8vLyMHz4cJSVleHVV1/Fn3/+idu3b2PFihX48ccfERwcjPJyeWWmGdLgq47ydX0C1+gqvoqJ3PSgVokEfr3hvo0Zw6qo1kCWFMrOAm4Ndak7gxkjgc8//1zofvL29tbb1qNHD1y9ehU3btxA586dTR7jxo0b2LhxI27fvi2sk2qM+BREvA5Su/cAnSEjhODw4cPIyMhAaWmp0Yg/Y0yaNAmUUsTHx+utnzJlCsrKyrBgwQJ88cUX6NSpEzw9PfHOO+9g1KhRyM7OxujRoyXryJDHwqG6wnq6sgZuSFFX6ZWmZtQmZYZrk6zjw9DBrhx0Idvr1q3DgAEDANQuJR4YGIjS0lJkZmbqhXCbws6uZsa6VGNk2PUnpR0ePpVRbm6uMF4ktRTGvHnzAECvLMXt27dx8OBBjBo1Ch999FEtmV27dqFjx47Yv3+/2fEmRt0Q4K1C0kRXxGfIy25d38SmlJkskWCLtiMC7OstDZAtv2tLgRkjABs3boRGo8EHH3yAK1euoGPHjnrb+S6zysrKWl6TMVxcalLQ8ymGLOHmpl/xUhwqLgUfHx88ePBACNXmE7Vaon379nj44Yf1xoASExNRUVGBVatWGe2SVKlUQrYINlepYcgvZUEjjOZNizdGWq0Wn376KQYNGoQnnngCQUFBcHDQT38i9lqkGAmxMZLqGXl5eenlu/Py8pIkx9OmTRtUVVUhOzsbzs7OenOLLDF27FhcuXIFVVW6gmhxcXF46KGH9LKSG9KrVy8MGTIE2dnZ+PXXX03ud+zYMUyYMAFdu3ZFUFAQ4uLihLx7DGlk5FYjctMDJE20rs4Oj3dcIZYfKUdYfLHwWdwWv17cLegdVyhsi00pw8SkEiw/WoGYn8rgHVco1CQyHEchSwpNHoMn+XylUV1iU8qEcTOxLlLaVhdoEbnpAbzjdDWexDrxMnybljweU+0Z6i4+jqnzaEkv8fmJ2Vmqd540ZVTvfDY3pBVdb8b8+uuvuHr1qtm6QGIDZM4zGj16NHbv3g1XV1dhnbOztJoqDg4OCAkJwa5du4T/5dChQwecOHECeXl5RucWmSMoKAjl5eU4ePAgwsPDkZOTYzY/H098fDz69euH7du34/HHH6+1ffny5YiNjUWbNm3w8MMP49ChQ1iwYAH+85//ID093WRUYn3z+MbaxcwmhTjg1QGOKKmkGP1t7dRKL/RzwAv9HHGvRIuoxNpjN6+EO2JyLwfcvK/FjG362399wa3W/paITS3H0sPl0HDvorgIp1o54pYfrUB8hv6LVFMGRAZYfqzzSijSo92Rqq5C5KYSRATYI9TfDiMTHuCrp10Q1dMBE5NKkHy+Umh3YlIJ4iKchf9jdpYirJ2drDpHhscQ66J7ET8QdBnQ3k6oSRS4pkjQJWmiq8W2w+KLkTTRFREB9tx3fICsN2oSIG/JrER6tDuSz1di1s5Ss9/BXHsncqqNHsfUebSkl/j8ZORWY2TCA+EcJGZWIqpn831lt3jPaP369fD09DRaEoJHPH5jzhgtWbIEERERei9ZS6mAxJw9e1byvoZMmjQJgC4AQm4XH2+8kpKSsGPHDgDQC083Rc+ePfHss89i8+bNtfLcHThwALGxsfDw8MCVK1dw8OBBZGVlYdCgQbhx4wZCQkLw4LnRdeAAACAASURBVIHpqqQtnbgIJxTEeoC+74GsN9xxIqe6VhXX+YN1WazFS0SAtAzbk3vpjEFEgD2ietpjy7lKJJ+vRIC3SjAUC4c6YUtmjRcQO6S2QZSLsWPwuoT62yE61AFbzunaFO8XFeyAE9nSPOr49ApEBNgjgjPKfFVVsUczmavPFNXTQTD4SjB2HFPnUYpe4vMT6m8HHxeCVLWux2JdeoVVpecbO03ezBJC1lFKY5TIFhQUYOvWrXjppZf0utYMEXe1GY7tiOHrDCkt0/Dzzz8DkGfAeMTh5IblzS0REREBQJdpnJ84a64Gk5gBAwYgOTkZycnJegZxypQpAHSBEfy4WUBAAI4ePYrx48dj+/bteOaZZ4SSHA2JOU/F1YGY3e7rqjK7vaOn+e1K4AMYxN5BXTKgnR2y8rVQF+gWvly2pozqRaeFt7O+lISlYwT6qJCVX9PVFne4HGm51dCUUUR0lfa6ysrXIsBL/3d2gJcKJ7KrhXNXV1F3xo5j6jxK0cvw/MSEOSIpsxLh7eygLtAKhqw50qQ9I0JIKIBJSuW/++47lJeX65X7NoY4Os5cip5nnnkGTz75JJycGr7eSmlpTdeQ1OAFHjc3N9jb2yMnJwfp6ekghOCxxx6TJPvkk08C0FW65dm6dSvu3LmDoUOH4qmnntLbnxCCrVu34pFHHsH+/fuFTOUMywR4q+olkOFETjUCfVQI8FZhUk8HpEe7Iz3aHVlvtELKjBrDaixSTVMmTx9L0W5Z+VoE+qiEl3lMuCPSo90RFVzbAJtqO9BHhYzb+l6UWqM7rlQ9jGGsPWPHMXUelegVHeaIxPOVSFVXYVIzL+XRZI0RIUTeCL8R1q9fj/79+wuTRqVgzhjt27cPR48eRW5uriJ92rZtC0CZZyQ2mHKDHwBd5vD8/HwUFhbCz89PUl48AOjTpw+cnJyE2k1ATRoiw0ziPCqVCnv27IGPj4+QBYJhnuVHypGWU11n5b/5rrBUdRWSz1chqqcDono6IPValdAtpC7Qmh3c93ImghfDTzb1cSHCoLzUUGhxe4mc56cu0MLHhSDAWwVNGa31EjfWNs+kEAek5VQL3V/J5yuRllMta2zLEHPtGWLqPCrRy8uZICLAHksPlyMmvPl20QFN2BgBiKCUWvxZTQiJJoSkEULS7t69K6zPyMjAyZMnLXpFhkh5SfN55uQiNdjBGO3btxc+iyfdSsXX1xclJSXQarW1qthaokOHDsjLy4NWq0VVVRUyMzPh6+uLXr16mZTx9PTE9OnT8fPPP+Po0aNG96GUIiEhAb1794adnR28vb0xZ84co7WbmhtLD5cjcE2REE3GD5TX5TyawDVFmJhUgqSJLkJ3U8oMN8QdKYd3XCEmJpWY7c6a3MsBiecrEbimSDBAMWGOiNxUInRRSYEfD4vcpBv0D/BWCeMpXT8twsiE2mOLxtrm8XImSI92x9LDuu+xLr0C6dHukvWR+l3NYew8KtUrJswRmjKKUP/GUXG3viBNMeklISQCQBqlVEMIKaCUWp78AyA8PJympaUBAF577TWsX78eubm5kuYO8d5KamoqRo4caXYftVotTDqVc355eZVKJTv8OT8/Xwi0WL16Nd58801Z8q+88gq+/PJLAMDixYvNRhcaEhUVha1btyItLQ1arRYDBw7Et99+a7ZIIQAcPnwYw4YNw9ChQ3Ho0CG9bZRSvP3221i1ahX69++P7Oxs3L17F5RS+Pn5ISMjw2KWivT0dITtHCH5e7QUvOMKsW+mW7N/uTUXks9XQl2grddy6+lP70f6L9/A96EOGP/y23rbCCHplFJ5ff8KaDSeEefBxJlZIrj9vADkU0o1StsqLS3Ft99+iwkTJkgyRGKkeEZKPZwVK1YokgP0u+bMBWOYQpyJ/NatW7JkR40aBUBnXE6dOgUAeOSRRyzKDR06FN7e3jh+/Hgto/3ZZ59h1apVmDBhAtLS0vDXX3/h3r17ePLJJ3H37l307duXJWtltAiWHi63qouxqdBojBGlNJ5SGmtmSeV2jQAQQAiJIoREAfDiDJm0/DcAfvjhB9y/f192Fx1g3hjx4dBy5wjxyEmOaoh4kqsSY8RPeAUgOwBj3LhxwueEhASoVCrJc50iIiJQUVEhRBICQHl5OebPnw87OzvExcUJ383Hxwc///wzJk+ejPz8fNneH4PRlIhPr4B3XCEmhzjUW5qjxkSTixOklCaL/yeEgFIab2p/Y6xfvx5du3Y1OlHTEuJAAUPeeecd+Pj4KDIGAGRlTTCHkvb5kG6gJmmrVFq3bg13d3eo1WpcunQJ9vb2ehN/zfHGG28gKSkJGzZsEBKvLlu2DOXl5XjxxRcRGBhYS+a7776Dv78/Vq9ejddeew1hYWGy9G3pFMRKS1HFsC3RYY4twiPiaTSekVwIIV6EkPnc5/lSPaOsrCwcOHAAL774oqKXvznPaOjQofjss8/g5uaG2bNn47vvvpN17LoyRlJrKIkRB0DIzeBACIG9vT2Sk5ORl5cnq/0hQ4bA2dlZr2vwv//9Lwgh+OSTT4zKqFQqLFmyBN7e3oiOjpalK4PBaJw0WWNEKdVQSpdTSgn3Vy1FbsOGDVCpVEKJcLmYM0YLFy4UKrSuXbsWU6dOlXXsujJGcsfBAP0sE2LDJBV7e3vk5+dDq9XWynpuDkIInn76aSEcvrCwELdu3UK3bt3M5vXz8PBAly5dkJGRIRQVrAWl0JIm5/wzGA2KltgD1PZ1mZqsMVLKxo0bMWrUKFn1gsSY66b78ssvhYgvJdSVMVJyHHHXnrGy6pbw9PQUCvTJPbePPPII/vzzT+Tm5mL//v0AgLlz51qUi4uLA2C8Um1VVRWK8v9CidfDsnRhMFoaJV4PQ1tm+wSsLcoY3b9/H9nZ2YoCF3ikRNMpmbQK1J0xMswTJwWxMVISSCHOh8d7h1Lhjfe2bdtw+vRpEELw/PPPW5SLiIiAu7s7fvvtt1o/AF5++WUs/XgNzvdfjGLvEOYhMRgGaIk9ir1DcHnAB7h9Qw2qpbBTGHxVF7SoJ/TevXvw8/PD008/rfgY5jyj999/H0uWLFF87LoyRlKDB8SIu9aURAOKgx4MUwBZgi9qePz4cVy/fh1t27aVFIRBCMHw4cOxa9cuHDt2DI8++igAIDc3FwkJCXBzcwOlbyOz32LYu3kDpEX99mIwzEO10JYV4vb1LBTl5aK0uAi+D8nvoq8rWpQxun//PubOnas4QwJg3jNavHgxFi9erPjYtoymE8so0UPctdeqVSsze9aGNyJXrlzByZMnZRnT559/Hrt27UJSUpJwnDfffBOUUnz44YfoN3Aotm9YjdwbWXDz8IJKZceVWV+H4uIHmDp1Ctq1s90DyGDYGkq1eFCoga9/e/QfGmkzPVqUMaKU4sUXX7TqGFJztimBNwJKu/l4lHg24owPSoyROLw6Pz9f1riTo6Mj7OzscPv2bZSWlqJbt26SZceNG4fWrVsjPz9fWHfmzBkEBwfj9ddfh0qlwrMvvIk/9u3A7VvXUFWpy5c2fsJEbNq0CRknTyMgyHx7VRWVKCkrg7ubG1R2zLtiNC/s7R0QENwXj4x8Bq7utgv7b1HGyM3NTXZ5BUPMddNZS115RkoMprVtiz0rufWUePns7GwAkGWM7O3t8dhjj+HIkSMAgEuXLuHSpUv47LPPhO/k7OqG4U/Xjmy8B298+umn+E9ijNFw9vz8fIwbNw6HDh0CpRSOjo6YO3cuFi1apKgrlMFgmKZF/cyri8qiDeEZWYsSz8jd3bpEkjdu3BA+KzFGfLVZ/rMcnJyccPXqVRQVFeH7778HAAwaNMii3P/93/+hqqoKsbGxtbZpNBp0794dBw8ehJeXF8LDw9G7d2/ExcVh2LBhECfdZTAY1tOijJGS+TeGNFdjZG0NJrGnYK4AoSnENZjkzFMCauZI7d27V0gr1L17d4tyPXv2hJOTE3bu3Flr27PPPou8vDxERUXh7t27OHHiBNLS0rBz506cPn0avXr1QkWFtBIJDAbDMi3KGNVFF1tz7aazJqgD0K+Gq6R9sWcmpeS5mCFDhgAAjh07hqtXr8LFxUWoLmsOQghCQ0Oh0Whw7do1YX1mZiYOHjyIp556ComJiXrXfOzYsXj22Wdx584dTJ48WZaeDAbDNC3KGNUFTcEzUnIcpcldeaz1Os+ePSt8ljvpls8QfvHiRRQUFMiadPvcc88BANatWyesW7VqFVxdXbFp0yajwSSJiYnw9fXF9u3bkZmZKUtXBoNhHGaMZNIQxsgWNaasNYRKxonEiI2Z3NB03nidP38elFJJXXQ8fMqmgwcPAtBlbkhISED//v31UiSJsbOzw8aNGwFAUlqp27dvIzMzUxgTYzAYtWHGSCZNoZvOFsjN9G2I2JjJDW0nhMDd3V2ocCsnAMLb2xtubm64efMmAJ3XU1lZadGgjRkzBh06dEB6ejry8vKM7nPp0iUEBgbC398fvXr1go+PD2JjY1FaWipZPwajpdB03342ghkj4/BjPkoLC1ob6RgaGiqUI5c7ljNz5kzcv38fWq0W33zzDQBIqpW0bNkyUErx448/1tp29epV9O7dG2q1Gg8//DBiY2Mxfvx4LF++HB07dsSdO3dk6chgNHea7tvPRtSnwWjKxoiPoPv73/+uSF5cqVYJ4i41udF4ffr0QVFREbKzs3H+/HnY2dmhb9++FuWmTZuGzp07Y/v27XrrKaV4/PHHUVlZiaVLl+LixYtYtmwZNm3ahNmzZyMvLw/9+/fXK2jIYLR0mu7brxlirTEaMWJEHWkiH3d3d1y4cEEYS5GLsSJ6csjJyRE+yx3X473dn376Cbdv35Y8/kUIwZAhQ7Br1y4UFBQI65OSkpCdnY3hw4djwYIFejJr167F2LFjkZOTY3U2EAajOcGMUSPC2jRAe/bssel4RI8ePRTPV1JStkKM2JDLnefEZ3w4cuQIysvLZWXp6Ny5M7RaLTZv3iys+/rrr+Hv76+3Tsy2bdvg6emJb775Rm+yMIPRkmHGSCIjR46s9zas9Yzs7e0Vj9nYGmtDy8XzlOQaI77kxcmTJwEAr7/+umTZadOmAQB2794NALh79y727NmDmTNnom3btkZl7O3t8cUXX4BSildffdXs8SsrK/HHH3/gwIED0Gg0kvViMJoazBhJZNeuXfU+6NyUx4ysxdpzK/bI5HqYnp6eIIQIXoqcMadevXrB3t4eGRkZAICVK1eiuroaw4YNMys3ffp0PProo0hLSzM5drRlyxZ4e3tj0KBBGDFiBPz8/DBy5Eg8ePBAsn4MRlOh5b79ZOLk5AQ/P796baMlGyMpGRPMYa1H6ODggOLiYgA6b0QO/v7+gjE9cOAAAGDo0KEW5ebPn4+//voLv/zyS61tSUlJmDJlCkpKSrBy5Urs3bsXAQEB2L9/Pzp37qyXpZzBaA603LdfI8TWxig6OhpLly61SdvWGiN+npPS0HtxaHlAQIAs2cDAQF2Z86IiXLp0CW5ubnrpkUwxZswYODk51SqbXlRUhOeeew6EEOzfvx9vv/02IiMjcfHiRcycORN5eXkICwtTVNGXwWisMGPUiLC1MVq3bl2t6K+GwlpjxFeLlVN+wpg8ADz00EOyZKOjowHo0hEVFhZKNmYODg7w9fXFmTNn9AzLW2+9hYqKCsydOxePP/64sJ4Qgq+//hojRozA9evXMXfuXFl6MhiNGWaMGhG2Nka2REl1WjF8olelEYl8AIRKpZJ9jC5dugDQjSsCQL9+/STLRkZGQqvVChNnKaX47rvv4OTkZNJL3b17N1xdXbFhwwaWOZzRbGi5b79GSEs2RtYWq7ty5QoAIDc3V5H8pUuXACjr5uO75LZu3QpAWi0lnkmTJgGoMWTHjx9HWVkZ3nvvPZOZ1J2cnLB582YUFRUhISHB7PHv3r2LrVu3YvXq1bh+/bpkvRiMhqblvv0aIS3ZGPGh3VLGWozBezPWRpopSfjKd8up1WoAulpIUuGnDKSnpwMAduzYAXt7e4vh5WPHjkWfPn2wevVqo4l1KaVYsWIFOnXqhKioKMydOxddu3bF0KFDWTQeo1HSct9+jZCWbIwAYN++fTh06JAiWd6zUtpNx89NkpPxm8fZ2RkqlUrIjSdnzMnR0RE+Pj64f/8+AJ131bZtW4vpkQghGDx4MDIzM41mvXjrrbcwf/58PProo/jjjz/w2WefoXPnzjhy5Ai6dOmCoqIi6V+QwWgAWvbbr5HR0o3R6NGjhUSlcrHWGFlbdp1PQUQIkd3VN378eBQWFoJSiqtXr0qWe/vttwEA69ev11ufnp6ONWvWwNXVFd988w0GDhyI119/HdevX8eMGTNw7949vcq6DEZjoGW//RoZLd0YlZeX28wY8dF8Z86cUSTPdzMqCYDo3r078vLykJaWhurqasnpiIKCguDh4SF08fHExMQAABISEtCuXTu9bQkJCRg8eDAuX74slGhnMBoDTfbtRwiJEi+21qcuaOnGCNBPeCoH3hgpjcpr3769IjkevltNSfl2Pp/gBx98AEDahFmegQMHoqysTAjg+Ouvv5Ceno527dphwoQJRmVSU1MRGBiIBQsWsLlKjEZDk3z7EULmAwClNBlAKoCFttWobmDGCBbT6JiCT+Ejp7CeGD6QQGlZh8GDBwNQZgx5nY8fPw6gpoy6FCIjIwHURPJ9/vnnAGq8I2O4uLhgyZIlOHPmjNlJzpRS/PDDD5g6dSpmzJiB/fv3S9aLwZCLXq59Qsi7AJTWj86jlK60XiVJLKSUegMApVQDIKyB2q1XWroxKi4uVuRZADWl2pWeQ94YhYaGKpLnowCVZC0PDg4GANy7dw+AvHIa48aNQ2xsrNA1ePLkSfj7+1uMxps0aRJeeOEFxMXF4R//+EetrsWioiIMHjwY586dE9Z98803mDZtGuLj42Uno2UwLGFY+CXeJlrIgBASAUDNdc1pAIQCSKaUqk3sHw0gGpBfdK2haenGyJoXHG8MevXqpUieP/cXL15UJM8HHhiO0UiBNz5VVVUICAiQlY4oKCgIrVu3Frrp0tLS8NRTT1kMUXdwcEBkZCR+/vlnHDhwQK8WFqUUw4YNw7lz5xAYGIjVq1eDEIJt27Zhw4YNuHbtGn7++WfFYfgMhlEopU1qgc6wUNH/XgCypMiGhYXRxoxGo6EAqEqlsrUqTZK9e/fSkpISxfIJCQn04sWLimQfeeQRCoDOnDlTtqxWq6UAKAA6cOBARW0HBATQnJwcCoC+9dZbkuRSU1MpAPr000/rrd+0aRMFQIcMGULLysr0tm3ZsoUCoIGBgbL1ZDRNAKTRBni3yyuJWY9wHoy5/okUSmkqADW3ANB10xFCAgghAdSEd9RUaOmekbXw4ydKmTFjhmJZPrRbSV0mQghcXV1RUlKiKINEZWUl1Go19uzZAwCorq6WJDdixAg4Ojri4MGDwjqtVosPP/wQ/fr1w8GDB2vdk5MmTcJnn32Gw4cPY8mSJXj//fdl68tgGMMqY0QI8aCUFtaFIpRSqV2ExgxOs6g6xoxR04Uf61I6abdr167IzMxUVEAvMDAQGRkZSEpKAiC9ECQhBD169MDZs2dRVlYGZ2dn7NmzBxcvXsTy5ctN3o8///wzfHx8sHTpUvzjH/+wujAigwFIjKYjhHgYWwBMIoS8U8866sF5P8ITSwjxAqBu6l4RwIxRU4Z/IVMjqXmkwM9zUhKNx4858fONBg4cKFl25syZoJQKwRN8NB4fVGEMd3d3REdHo7y83KxndOfOHbzxxhvo27cvBg8ejC+//FLx+WG0AKT05QF4F8BVAFu4JZFbrgK40hD9iQb6BACIg278KA5AgBS5xj5mVFpaysaMmijTpk2jAGi3bt0Uyfv6+ioei/nqq68oAGpnZ0cJIVSr1UqW/f333ykA+sMPP1BKKXV3d6cODg60srLSrFxxcTFVqVS0TZs2Rrf/8ccf1MfHRxgL45f27dvTmzdvSv9yDJuDBhozkvRTnFK6Arrw6XgACyilkyilkwDEAphtvUmUB6VUTSmNpZTGc3+bvFcEMM+oKcNPVOXHjuRCOY+hVatWsmX5Gk7V1dVwcHCQlQGiX79+IIRg9erVuHfvHoqLixEUFGTxe7i5uWH8+PEoKioSJu3y3LlzB2PGjIGnpyd++uknXLp0CRcuXMDIkSORnZ2NsLAwlqyVUQvJbz9K6X1K6T4AIISM4LrpKL+OYT3MGDVd+G46pemI+GsvrjgrFXEouJwJs4AuyauLiwtOnz6NP/74AwAk562bNWsWSktLsW+f/ivgpZdeQn5+Pr7//nuMGTMG3bt3R48ePZCamop///vfuHPnDt59911ZejKaP7LffpTSa5TS/dBFvil78hhG4V9IUgegGY0HPnuCkhIUQI0RGzdunGxZvuS64Wep+Pr6oqioCIcPHwYAPPHEE5LkHn/8cTg5OWHRokXCutzcXOzatQuOjo5Gc+y99957mDdvHtauXYstW7bI1pXRfLFojAghiYSQPYSQ8YQQYWYcpfQkpXRr/arXslCpVLhw4QK2bdtma1UYMuHLR3To0EGRPP9DREkGCicnJ0E+Pz9ftnz79u2h1WqRnZ2N1q1bY9q0aZLkHB0d0aZNG5w6dUqoOMvXV4qOjjaZCX3x4sVwcHDAiy++KDkMndH8keIZbQEwiVL6A6V0PyFkAiGkS/2q1XLp0aMHS7XSBOHLRlRWViqS59MI8QX6lFJQUCBbhjegZ86cQY8ePWR1F4eGhoJSioyMDABAcnIyAGDOnDkmZVq1aoXJkyejpKQEn376qcn99u/fj8jISLRp0waBgYGYPXu2YPQZzY9adx0hZCQhpB//P6V0K6X0vvh/AIHcmBGDwUCNMTp16pQied6LUJqo1fA4cuATtZ49e1b2pNtRo0YBAH788UcUFhbi2rVr8PPzs5iw9pNPPgEArFq1yuj2jz/+GCNHjkRmZib+/ve/g1KKdevWwd/fH5cvX5alI6NpUMsYcQEJ9zkPaDy3dDGyT0TDqMhgNH6sDT7hPSpvb29F8nwpCCXRePy8JK1Wi7KyMlmyTz75JADg999/R1paGiilksqu+/n5ITg4GLdu3UJ2drbetl9//RXvvKObvrhixQrEx8cjIyMDM2fORGFhIfr27Ytbt27J0pPR+DH6BHFBClu5rrkfABDOOE0wZpwYjJYOnzRUaWg3/0JWkvXbmB5yEAc9yJXv3LkzCCEoKCjAhQsXAOjGhKTw/PPPAwDWrl0rrKOUCuuXLVuG6dOnA9DVi/r666/xwQcfoKysDEOGDJGlJ6PxI3WeEW+ctvLGCUCkyHNiXXaMFk337t0BKM88zntE1maW54v8yUGss1zPTKVSoWfPnujatSvOnTsHd3d3yZnLZ82aBUDfAO/Zswc3btxAx44djYZ///Of/8SIESNw48YNfPvtt7J0ZTRupETTTSCEbCGEvMwbHUrpNQD3RZ6T9AIsDEYzhPeIlNZjeuihh/T+KiU6Olq2jLhrT4kx69ChA7Kzs7Fnzx6UlpZKDuLw8fFBSEgIjh49Kqzjx5CWLFlisutzx44dCAoKwsqVKy2mFyovL8f169eFaD9G40WKZ+QDXeqfJwBcJ4ScIITkATjO70ApPVlP+jEYTYLWrVsD0GW1tgZru+mUjDmJPSMlnhmlFKdPn0ZOTg6cnZ1lGeTg4GDs27cP169fR2VlJY4dO4YRI0Zg8uTJZvVduHAhTp06ZXKuklarxYoVK+Dn54euXbvCzc0NEydOrJUtgtF4kGKM8qGrFzSJUuoDXT64AM4jYjAYqPGIlI5l8L/w+ag8pZw+fVq2jDgCLyoqSra8p6cnKisrUV5eLtsYhoSEoLKyEtu2bUNmZiYKCwvx8ssvw9XV1azctGnToFKpTGZymD9/PubPn4/HHnsMHTp0ACEEycnJaNu2rVCIkNG4sGiMuFBuwgctcJNd75sVYjBaGHwGhaKiIkXyvDFSmk6I59KlS7JlxJnCLRkBY7Rv3174zHuIUuGzPRw5cgS//fYbAGlZJJydndGzZ0/cunULt2/f1tu2bds2fPzxx5gyZQp27tyJmzdvori4GLNmzRKi8Qwj+Bi2R2oAw0lK6fV61oXBaLLcuXMHAJCQkKBIvq6MkZJJt+KxmbS0NNny4oAFuWNefJn4S5cuCTnupH4HPupuzZo1wjqtVisERkyfPl04n46OjoiPj8eSJUtQVlaGOXPmsHIWjQy9OFRCyLvQjREpIY9SutJ6lRiMloe1xugf//gHPvroI0XdbGKUDPT7+/sLnx977DFZsh4eHnBwcEBubi6Ki4sBwGhOO2O89NJLePfdd7F792589NFHAIDNmzcjLy8P4eHhGDt2bC2ZRYsWwcnJCQsWLEBiYqLZsSlGAyOuJwHA05qlIWpeWLM09npGjKbLjh07KAA6f/58RfIHDx6kffr0oaWlpXWsmTTA1Rtau3atbNlffvlFkD958qRseT8/P+rp6UldXFyog4ODrHpMXl5e1MXFRfh/wIABFAD99ddfTcpUV1fToKAg6unpScvLy80ef9++ffTNN9+kS5cupfn5+ZL1ak7AFvWMqK5MhOKlvg0ng9FY4dPfhIaGKpIfNmwYTp8+DWdn57pUSzZK5kn5+fkJn5WUIH/00Ufx0EMPobS0FG3atJHlHT7xxBPQarVC9oiMjAy4ublh2LBhJmVUKhUiIiJw//59fPzxx0b30Wq1eOGFFzBy5EisWbMGCxcuRI8ePXD16lXZ348hDVZAh8GoA/icctZGw9kKPqJOSW47vmQ6AOzdu1e2vJ+fn5BvTtzlJ4UnnngC5eXlyMrKwpEjR1BdXY3/+7//s5ieadmyZSCE6GV/EDNlyhR8/fXXsLOzw6BBg+Dv74/S0lKMHTsWhYWFsnRkSIMZIwajDuAn0UDZfgAAGR9JREFUi7Zt29bGmiiDD8lWUg9JnEJISToicXLYBQsWyJINDg4GoEvUyiepff/99y3KeXp6IiQkBDdv3oRGo9Hb9tNPPyEpKQleXl64dOkSjh49ilu3bmH79u24evUqRo0aJYxvMeoOq4wRSwPEYOjo2LEjKKVC+fGmRo8ePfDoo48qmiclDgdXMulWHFouDhOXAt89unXrVvzxxx9o27at5Gq5EyZMAADEx8cL60pKShAdHY2QkBCcPHkSgYG65DIqlQojRozAggULcOzYMVYAsx6QZYwIIV25YntXCSFXAGQQQq5w6YK61IuGDAaj3rG3t1dcvkJsTJRUuhWHg8vNxu3n5weVSoXs7Gz8+uuvuHv3ruSCfa+88goA/XD2VatWITc3F2vXrkWXLl1qyfz73/9G9+7dcfz4cXzxxReydGWYR65nNJLqMjEEUUq78X+hy8pgXUwpg8GwGVVVVThx4gSOHTsmW1Y8PqMkt564a1NucUBCCNzc3JCfn4979+6hdevWksftHnroIfTr1w/37+tiryilWLVqFRwcHBAeHm6yvV9//RV2dnaYN28e7t27Z7Gd3377DUFBQRg8eDCuXbsm/cu1MOQaI6NnkoukY/npGIwmCj9WU15ebtVxunbtKltGnLVBbgYHQDf+U1JSAkqp7ACIfv364fTp09BqtTh37hzy8vLQsWNHPW/PEH9/f8ybNw8VFRVC3SVjVFdX4+jRoxgzZgyys7Nx6tQpTJo0Sag9xdBHrjHyIoQs5TJ48+UjXiaErAUg/y5kMBiNAr6LTmk9Jh4liV7F0XhKovlat24tTBrmS6hLxcnJCX/99RcOHjyI77//HgAwdepUi3L//Oc/0apVKyQmJhothR4XF4eQkBA8+eST8Pf3R0hICEpLS5GWloavvvpKlo4tBVnGiOry1MUD8AYwkFu8ASynlP637tVjMBgNweuvvw7AemOkBLEBUzLPasCAAcLnzp07y5Lt27cvAGDfvn04ePAgACAyMtKinKenJ5KTk1FaWiqkMQJ0ARCDBw/GggULcOnSJQQGBuLAgQNIS0vD5s2bQQjB3Llz8eDBA1l6tgRkR9NRXaG9FZTSBdyygurqGzEYjCYKn2zUFsZIPFFWiWcl9obGjBkjS/bRRx8FAJw9e1aoVNu7d29Jso8//jjc3d3xww81BQyWLFmC33//HR4eHti3bx9OnDgh6Dd58mS8+uqrKC0txXPPPSdLz5ZAnc0zIoT0q6tjMRiMhoUfa1ESmm0tYgM4cOBA2fLiMZhBgwbJkuUj5q5evYr8/HyMGTNGckQgIQSVlZX48ccfAegSvH7++ecAgI8//hgjRoyoZdzj4uLg4uKCsrIyk8e9e/cuxowZgxdeeEEvaSylVHHEY1NAbmi3JyGkn7EFQEw96WhMj1BCSDQhJIoQMp8QEtBQbTMYzZEhQ4agQ4cOCAho+EdJ7BkpSRQrDuW2lHnBEE9PTxBCcOPGDQDyKuU6ODjAx8cHBQUFKCgowLFjx1BSUgJnZ2dMmzbNqIybmxvmzJmDvXv34s8//zS6z+LFi7F79258/fXX+O9/daMfxcXFGDRoELy9vZGamirrOzYV5HpGPgCWA5gMYIrBYjwWsn6IoJTGU0qTKaXLAcQ2YNsMRrPD3t5e8vycukZsjPLy8mTLi0ulyy2cRwiBt7e34KmYi6IzRr9+ug6hs2fPIiUlBYQQ7Ny502xdqNdeew0AMGrUKGRmZuptq6ysxHfffYfnnnsOYWFhWLp0KRwcHNCpUyccP34cWq0Wc+bMaZYReXIDGK4BiKOULhSNGS2glC4AsKx+VDRKDCHEy/JuDAZDCmq1Grm5uYrnwXz77bd6dYXkIO7KUlKPqVWrVsJnsWGSSmBgoNDu/v37ZcnyJTMOHz6MI0eOIDQ0FBEREWZlOnXqhLFjx+LSpUt62R8AIDMzExqNBoMHD8brr7+OmzdvIioqSojY8/X1xcWLF7Fjxw5ZejYFlAQw7DOxfqv16kgmDsA1rqsuGmY8I26fNEJI2t27dxtOQwajCcEPqCutpzRt2jTMmTNHkazYM1ISTWdtaLh4nIxP/yOV4cOHAwAyMjJw6tQpVFVVmR0P4pk/fz4AXTFGPiwdqMkG8e677yI4OBju7u7YunUrysvLsWTJEuTl5cHBwaGWEWsOWDRGhJCXzWzrYovABUppPICl0I1TxcBMQUCuOy+cUhouTnXPYDBqsGXWcWuj6cTGyNHRUbZ8fn6+8FlctVYKvXr1goeHBzQaDfLz83H58mVJ32Hw4MHo3LkzNBoNzp49K6w/duwYCCHo06cPwsPDMXPmTFRWVqJXr1547733sG3bNkyfPh179uxpdqXTpXhGRn8qEUImAEiGrstsvLWKcB5MnJklQrTvfErpckppGIB1AFKsbZ/BaMksWrQIgG1Cu8VtKjFG4kzjSoyR2JuS283XqlUrPPzww0IJjK5du0ryLvn5RgD0yljs27cPlFJ88MEHsLOzw7Jly/DJJ59g165dUKlUiIyMxHvvvQetVotPP/1Ulq6NHUl3HiFkL4AwAImU0le41dEAXqaUnuLKlVsF5+1I0SUCQIZYjhASSAgJpZRmmBFlMBgmKCoqAmD7eUZyo+EA/TEjuQEIgL4xUtrNd+LECQAwmlzVFC+//DLmzZuH8+fPC+tycnLg6emJv/3tbwB03403WjyBgYHo3r07Vq5ciejoaCFzeVNHypUPhG6MJhxAqsgLCgSg5j5rjAnWE/kAapXTZIaIwVBOjx49ACh7mVuLkuqwYoKCgjBw4EAsWbJEkTHljZm/v7+i3Hr8ZFkA6N69u2Q5Nzc3jB8/HpcvXwalFLm5uaioqEDv3r0telejR48GpVTInNEckGKMjlNK93GZF8RBCpRSypc8pMYE6wPO6Kj54AVCyHwAWxqqfQajOcJnLlBSdtxarPXG7O3tcfz4ccXzb/gxJz8/Pz0vSyriAIhu3brJkh09ejRu376Ns2fPCoEP48aNsyg3b948EEKwZ8+eWsUBmypSjNFAQkhnQogH5xXxwQKtCSH8lWvQMGtuflE8tyxnXhGDYR184IIt5hpZ6xnxHDp0SJEcn32ipKRE0fwdsQGTM2kWqKmMu379euTk5AAAevbsaVGuY8eOwhyn3377TW9bfHy8XlBEU0GKMVoHYCuAAugmuxYQQpZBF079D85A3a8/FRkMRkMhDjNuKOrKGCmFr6d09epVRfLiaD65Xl5oqG7E4dChQ/jf//4HQGdopDBp0iQAEBK88ri6umLbtm2y9GgMEGtuPkJIfwDhlNImkRM9PDyciqs6MhgM20MpFQIXlL6P+DEWJfKJiYmYPHmyYvmpU6di8+bNAHQZJORUu62srISTkxNatWoFOzs7FBQUQKPRCB6TOW7cuIE+ffpgxIgRmDJlCi5cuID3339ftv6WIISkU0rrPcOOZDNOCPGALoghjR8ropSeBCuqx2AwrEDpRNu6Qu44jyH9+/cXjJHceVoODg7w9PQUxn3s7e31PC1zdOrUCX379sVff/0lGFNPT09MnTpVUcVdWyMpjpIQ8iV0EXPJ0HXTsYABBoPRaPj73/8uuXvLkF69egFQlkoI0I+gM5eTzhTisutubm6yjLOHhweOHj0q/D937ly9khZNCSkZGN4BkEQpVVFKfSildgASufUMBoNhcxYtWoSNGzcqknVwcEBGRgZu3bqlSF48TqQkMjAqKkr4LNUr4jGWhWH06NGydWgMSPGMrhnmo+NCvFnQAoPBqDMefvhhxbJ79+5VnOQV0HW1KQ1rF2caV9Ll+N577wmf5XpnhvOinJ2d0alTJ9k6NAakmHFTI3ryc70zGAyGEa5duyZr4N+QBQsWoHXr1njppZfqUCtpWJsBwdHRER4eHigsLNTrspOC2Bg5ODigVatWNh+DU4qkDAxc8IIAIaQLAPklGRkMBsMIXbp0kd1FZYiSWkh1AT9PSSkrVqxAYaEuf4Bcz6hz587C5+rqakW5+RoLUjyjeAD7CSEUulQ8/M+XkfWmFYPBYDQR+OSufFkIuYjz4UkJ6RYjNl5DhgzBoUOH0KlTJxw5ckRxQIetsOgZUUrvczHmCwGkAlhGKR0gSgXEYDAYNiUlJQUHDhywSdt8Dably5crkhd7M3Kj8cQThvm0RDdv3sTFixcV6WJLJId+UEpToTNGAHRRdpTSlfWiFYPBYMjAUnXV+sTDwwN2dnaKJ+yKy2bIzUYhNmTiAI4nnngCVVVVNqlPpZRaxogQMgtAlJF99XaDrqQEM0YMBqNF4+DggGeeeQZXrlxRJC82KHJDw8WyhvnompIhAox7RoHQTW5VG9nGQwAo6yBlMBiMZsbu3bsVZx/nJ90C1nlGTR1jZ28Ll+bHLIQQFtrNYDAYAMrLy1FeXq5INiQkBB06dMCtW7esMkaOjo6oqKhQpENjoFYAgxRDJGc/BoPBYJimsrISQ4cOBWBdN52ts59bS8PXGGYwGIxmRtu2bfH0/2/v/nniuNY4jv8eJZbipFnjayndldapUq5xUqVbl+kgf16A4QVEAt0yTST8BiKTvAEE7wDSpEpiTJMqhfem970bpNtGem4xZ/AwzCws7M5zGH8/EjLM/nu8Z4Yf58zsOZ9/fq3H/vrrr2cTrc4bKNX7Ryz/sUjzLzgPADjnvffeO1updV7Vq+lu0jN6/PixpGKuu19++eVatUQijADghv7880/t7V1vMYObDLWVj71z544++aSYFOe7777Tp59+eq1aIjFMBwA39PPPP5996HRei/ickbvro48+0u+//37j9Zmi0DMCgBv67LPPzl2iPY9qGJUr3l5VGUZ///23fvzxx2vXkAPCCAAC3b9//+z7eWfcrg7xvXjxYmE1RSCMACDQ3bt3z+akm/eKuOqw3nXXY8oFYQQAge7cuaOvv/76Wo8tr767e/eufvvtt0WW1TnCCABuqQ8++ED37t3T999/r48//ji6nBvhajoAuKXeeecdTafT6DIWgp4RACBc9mFkZocN24ZmtmVm4/TvfGv1AgCyku0wnZmNJQ0lNa2ate/uj9L9jiX9IGm9w/IAYGHmvaS7j7LtGbn7kbvvSjqtbjezkaRp5X6nag4sALgVbvskp4uQbRjNMFQtoCRNzWzYdGcz2zCzYzM7fv369fKrA4A5lVMJlZ83ehtlO0w3w0rL9sbzRql3tStJq6ur/PkBIDvffvutPvzwQ3311VfRpYTpNIzMbEPFsuZtDt396JKnmepi8LQFFABk7/3339c333wTXUaoTsMo9VJuatLy3CcLeG4AQIBbd84ohc5ZTyidK7qsNwUAyFi254zSVXNjSQMz29H5IbynZrYl6UTSSNLToDIBAAuQbRilHtCJpGczbpPoFQHArXfrhukAAP1DGAEAwhFGAIBwhBEAIBxhBAAIRxgBAMIRRgCAcIQRACAcYQQACEcYAQDCEUYAgHCEEQAgHGEEAAhHGAEAwhFGAIBwhBEAIBxhBAAIRxgBAMIRRgCAcIQRACAcYQQACEcYAQDCEUYAgHCEEQAgHGEEAAhHGAEAwr0bXcBlzOzQ3Z/Uto0krUoaSHosadvdJxH1AQBuLtswMrOxpKGkcW37QNKqu+9W7nco6WHnRQIAFiLbYTp3P0qBc1q7aShpu/LzsaRhCikAwC2UbRi1cfcTSY8qm1Ylnbp7PbQkSWa2YWbHZnb8+vXrTmoEAMzn1oWRJNWCZ1PS0xn33XX3VXdfffDgwfKLAwDMrdNzRma2odnndg7d/WjO59tz94MbFwcACNNpGJUXHSxCunBhMk94AQDydCuH6dKl3dMyiMxsLbgkAMAN5Hxp90jFZd0DM9tRGsIzs6Gkn9L28u4TSQzVAcAtlW0YpavmTiQ9q22fSLoXUhQAYClu5TAdAKBfCCMAQDjCCAAQjjACAIQjjAAA4QgjAEA4wggAEI4wAgCEI4wAAOEIIwBAOMIIABCOMAIAhCOMAADhCCMAQDjCCAAQjjACAIQjjAAA4QgjAEA4wggAEI4wAgCEM3ePrqEzZvY/SX9E19HgH5L+E11ETY41SdQ1L+qaD3Vd9E93f7DsF3l32S+QmT/cfTW6iDozO86trhxrkqhrXtQ1H+qKwzAdACAcYQQACPe2hdFudAEtcqwrx5ok6poXdc2HuoK8VRcwAADy9Lb1jAAAGSKMAADhCCOcY2bPaz8PzWzLzMbp30FUbbjIzA4btrW2WVft2VLXyMw20uvum9mwy7qaaqrdHrLvt9VlZmvVr67r6py79/pL0lDSlqRx+ncQWMtI0kaqY1/SMKc6U31/1ba9rHw/kLTfYT1r1a8c3qtKG66l1w5pw/QaG8UhfOG21jZbdnu21ZVea6N2v1dd1DXrvaq1a6f7/iVtuFXu8+m1X3ZVV9RXeAFL/w9m0nCRB+Mc9Z07INPPh7X7/dVRPVkejJK2aj8/j6yr4Rdoa5t12Z4tdVX394Ekr+x3S6+r7Tmj9/2m551Ra9gxueyvXg/TmdlI0rT82d1PVYRAhKGk7crPx5KGZjbIpM6xu5/Utg0lnda2TavDK0v0L3c/kIr3w90fSVm06WbTsEgGdZVmtVlYe6Z961Fl06qk0/Q+Re5nUmb7vpmNJU3S8Fw5FFe+bvR7tTR9nw6oteHcfdJlIe5+YmaNB2PakcLqTDv/UcNNKy0PWeoYdfVgVPG+jCQdpPciuk13JP3bzMo/LMp/o+sqzWqzkPYspeApbUp6mr4Pqyu3fT8ZShqVf4yZ2bGkl5IeBte1VH0Po6waLtODcSBpWqutNG2ooa3WRcr2YHT33fSebaZNRypCKJd9bVabRbXnOWa2IWmvbF8F1ZXpvi9Jk/Qlqfi9kS5aGAbXtVR9D6MsGy6XgzEZp5rKbv4g1XekygFR1TCksWjZHoxmtuXuzyQ9S+/ToYqQzGVfa20zM2t8QAfteabs9bp7tTcStZ/luO+r5bVPZ9zWaRsuS9/DKLuGy+xgVCUQJUlmJnffrfy8Uvl+qOYhjUXL8mBMbXf2OqmX9DCdL8piX0uh09hms27rQnlerXxPzGzN3Q+i6sp035e7T8zsrLeWenCTcrg3sg2XqddhFH3w1eV2MNZqG6i4zFRmtqU352iepp9PVJy7edr+LIuR8cE4VfHX9LnXq7RnZ3WlfWms4q/5HRVXWJWvN6vNltqebXWl9+OntL28+0RSGQhLq+uS9yps37+krvW07ZWKnvd65aGdH5Nd6P3cdJUGLxtut2WMeNl1DFWc96gO5Uzc/WFOdeYivV+benMwPq+EUdh7lS6qKENnIOmoEka0IXBNvQ8jAED+ev05IwDA7UAYAQDCEUYAgHCEEQAgHGGE3uvDvF2zzPP/S5cLA9khjNBr6Rd1yOS4aZLLQzNzM9tJE14+t2Itn7XLn+HKxlcJpLYJXlM9nj7U2/bYQzN7RZhhWXr9oVdA0ra7b9Y3mtnLcibwZal82HPF3asztsvM/jKz09pMHNd9nd0UEtuX3HVD0rkF5NIHrvfSjyM1fFA3hdSKpJ3qDAXAItEzQm+lecb2G7aPVVnuYcmeqHkmhmm6bVH20qfyZ3ncMoP4iorZEB7WbygnE1VLUAGLQhihz9Zbeh5PVExw2oVx/bXszdpCe42PuIY0C0RruKXXfDHjKV6lmurKobvTrpddwduFMEIvzTg/sqViuOr+FXoSN61hqGLp8fo8aDuSHi1hEtXTGeeONiVdGGJL9z9W0TMa1m4re0Nta/4AC8M5I/TVqiozbEtn50cmKlaRvez8yiKUiwSOVcxj96WKiXLXZz/s2l6oGE5r6sEMWubJG7n7QQrJepCtpPfsiRqGO4FFomeEvhpI+m/D9rGKnsBMViwH3/p1xRqeqJgB+ijNzr4uaTWdy7o2M2sLhsZF/lIP52XLY1akNws/lv+3NKN82RuiZ4Slo2eEvmpb2O7Sv/LTZdePL7mPrtC7Guv81P9S0Wu59lV8s9ZOSpqCclOXX2mn9LzDtHxHdSZyzhdh6Qgj9NVUDVeHqQiIHak4X9L0SzYtunZQ3z6PpvNFyUi1y6trS09M09BY+TmkYdr+RMUFDzuSDs1s1HDOaaD2Rf4uDNFVzheVyvNGp/SK0DXCCH01UfPVZStp8b7yHMmy/uI/typsxVBp1dpKL+eH8jNPZraWeiZlAOy7+zMzK0NKadnzJvfrr5nOV7VdOTiufW5oImnT3avvG+eL0AnCCL2UfnE3XVm2W/Y66stOL0Jl1dDN9PNG7Rf+tqSHqYaJpC/05iKHczWl+xyW/58rvPyoYehwvf6h38oVfV+Y2aPK7S+Vem0pKL9UEaqnZnbc9RLqeLuwuB56y8yeq5iBIdvVVuvBmIJiRcUw446KFW5PUjhI0pfuvl0fpisDpiF4njfNQAHkhjBCb6We0dqMYa0spM87HenNTAhlr6m8CGNaufx6Q8Vy9QdNz1ELqDWdP/8DZIswQq+ly6iP+nw1WDkZbH3eODPbX+JnmoCF4nNG6LX0Czpk1u4OXQiipLcBjP6hZwQACEfPCAAQjjACAIQjjAAA4QgjAEA4wggAEI4wAgCEI4wAAOH+D2oSi7XAhIzHAAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%matplotlib inline\n",
"\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from matplotlib.pyplot import savefig\n",
"\n",
"from matplotlib import rc\n",
"rc('text', usetex=True)\n",
"\n",
"\n",
"if Nr == 270:\n",
" extraction_radius = \"33.13\"\n",
" Amplitude = 1.8e-2\n",
" Phase = 2.8\n",
"elif Nr == 800:\n",
" extraction_radius = \"33.64\"\n",
" Amplitude = 1.8e-2\n",
" Phase = 2.8\n",
"else:\n",
" print(\"Error: output is not tuned for Nr = \"+str(Nr)+\" . Plotting disabled.\")\n",
" exit(1)\n",
"\n",
"#Transposed for easier unpacking:\n",
"t,psi4r,psi4i = np.loadtxt(\"outpsi4_l2_m+0-\"+str(Nr)+\"-r\"+extraction_radius+\".txt\").T\n",
"\n",
"t_retarded = []\n",
"log10abspsi4r = []\n",
"bh_pert_thry = []\n",
"for i in range(len(psi4r)):\n",
" retarded_time = t[i]-np.float(extraction_radius)\n",
" t_retarded.append(retarded_time)\n",
" log10abspsi4r.append(np.log(np.float(extraction_radius)*np.abs(psi4r[i]))/np.log(10))\n",
" bh_pert_thry.append(np.log(Amplitude*np.exp(-0.0890*retarded_time)*np.abs(np.cos(0.3737*retarded_time+Phase)))/np.log(10))\n",
"\n",
"# print(bh_pert_thry)\n",
"\n",
"fig, ax = plt.subplots()\n",
"plt.title(\"Grav. Wave Agreement with BH perturbation theory\",fontsize=18)\n",
"plt.xlabel(\"$(t - R_{ext})/M$\",fontsize=16)\n",
"plt.ylabel('$\\log_{10}|\\psi_4|$',fontsize=16)\n",
"\n",
"ax.plot(t_retarded, log10abspsi4r, 'k-', label='SENR/NRPy+ simulation')\n",
"ax.plot(t_retarded, bh_pert_thry, 'k--', label='BH perturbation theory')\n",
"#ax.set_xlim([0,t_retarded[len(psi4r1)-1]])\n",
"ax.set_xlim([0,final_time - float(extraction_radius)+10])\n",
"ax.set_ylim([-13.5,-1.5])\n",
"\n",
"plt.xticks(size = 14)\n",
"plt.yticks(size = 14)\n",
"\n",
"legend = ax.legend(loc='upper right', shadow=True, fontsize='x-large')\n",
"legend.get_frame().set_facecolor('C1')\n",
"\n",
"plt.show()\n",
"\n",
"# Note that you'll need `dvipng` installed to generate the following file:\n",
"savefig(\"BHperttheorycompare.png\",dpi=150)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 7: 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-Psi4.pdf](Tutorial-Start_to_Finish-BSSNCurvilinear-Two_BHs_Collide-Psi4.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": 19,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Created Tutorial-Start_to_Finish-BSSNCurvilinear-Two_BHs_Collide-Psi4.tex,\n",
" and compiled LaTeX file to PDF file Tutorial-Start_to_Finish-\n",
" BSSNCurvilinear-Two_BHs_Collide-Psi4.pdf\n"
]
}
],
"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-Psi4\")"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}