tutorial.html

<h1 id="neutron-star-merger-simulations-tutorial">Neutron Star Merger Simulations: Tutorial</h1>
<p>This tutorial will guide you through the steps necessary to perform a binary neutron star (NS) merger simulations with tabulated nuclear equation of state (EOS) and neutrino leakage.</p>
<h2 id="changes">Changes</h2>
<ul>
<li>June 9, 2020. Removed discussion on the <code>mass_conversion</code> parameter, since it is actually not needed when using the tables and scripts provided with WhiskyTHC.</li>
</ul>
<h2 id="downloading-whiskythc">Downloading WhiskyTHC</h2>
<p>Download and unpack the WhiskyTHC core distribution with the commands</p>
<pre><code>    $ wget http://personal.psu.edu/dur566/downloads/WhiskyTHC.tar</code></pre>
<p>The tarball contains scripts and microphysics tables for WhiskyTHC. It does not contain the actual source code, which is instead hosted on <a href="https://bitbucket.org">bitbucket</a>.</p>
<p>Unpack the tarball with</p>
<pre><code>    $ tar xvf WhiskyTHC.tar
    $ cd FreeTHC</code></pre>
<p>Use the GetComponents utility to checkout Cactus and WhiskyTHC</p>
<pre><code>    $ ./scripts/GetComponents --parallel thornlists/full.th</code></pre>
<p>This will download <a href="http://cactuscode.org/">Cactus</a>, <a href="http://carpetcode.org/">Carpet</a>, <a href="https://bitbucket.org/llamacode/">CTGamma</a>, and <a href="http://personal.psu.edu/dur566/whiskythc.html">WhiskyTHC</a> and place them in a folder called <code>Cactus</code>. This will also download <a href="http://bitbucket.org/dradice/batchtools">batchtools</a>, the set of scripts I use to manage simulations.</p>
<p>GetComponents has some problem with the capitalization of repository names from bitbucket, where WhiskyTHC is hosted. To workaround these issues, it is necessary to run the <code>fixsymlinks.sh</code> script after <code>GetComponents</code>.</p>
<pre><code>    $ ./scripts/fixsymlinks.sh</code></pre>
<p>NOTE: <code>GetComponents</code> will place a copy of the thornlist used for the checkout in <code>Cactus/thornlists</code>.</p>
<h2 id="compiling-whiskythc">Compiling WhiskyTHC</h2>
<p>Detailed instructions on how to compile can be found in <code>Cactus/doc/UsersGuide.pdf</code>. In brief, you need to load all necessary modules, create a configuration, and then compile it.</p>
<p>Example of module lists for some of the machines I have compiled and run WhiskyTHC can be found in the <code>Cactus/batchtools/templates/cactus</code> folder. If you are using <code>bash</code> as login shell, you can load them with the <code>source</code> command (in the <code>Cactus</code> folder):</p>
<pre><code>    $ source batchtools/templates/cactus/comet.env</code></pre>
<p>To create a configuration you will need to prepare a thornlist, that is a list of components to be built, and a config file. An example of thornlist is the <code>full.th</code> file, which contains everything I typically use in a production run. Example of optionlists for some machines where I compiled and run WhiskyTHC can be found in the <code>batchtools/templates/cactus</code> subfolder of Cactus.</p>
<p>For example, on SDSC Comet you can create a configuration</p>
<pre><code>    $ yes | make thc THORNLIST=thornlists/full.th \
                     options=batchtools/templates/cactus/comet.cfg</code></pre>
<p>If this command is successful you can then compile with</p>
<pre><code>    $ make thc -j4</code></pre>
<h2 id="microphysics-tables">Microphysics Tables</h2>
<p>WhiskyTHC uses tabulated properties of NS matter. To run a simulation you will need three tables:</p>
<ol type="1">
<li><p>A zero-temperature beta-equilibrated slice of the EOS to be used with <a href="http://www.lorene.obspm.fr/">Lorene</a> to generate the binary initial data. You will also need a version of the same table in “Pizza” format which is used by WhiskyTHC to initialize the electron fraction from the Lorene data.</p></li>
<li><p>An equation of state table in HDF5 format to be used for the evolution.</p></li>
<li><p>A table with chemical potantials used to compute weak rates.</p></li>
</ol>
<p>The format of these tables is documented in the code, see the <code>THCExtra/EOS_Thermal_Table3d</code> and <code>THCExtra/WeakRates</code> thorns. <code>THCExtra/tools/eos_table</code> contains scripts that can convert tables from <a href="https://sntheory.org/equationofstate">stellarcollapse.org</a> to the format used by WhiskyTHC and can extract zero-temperature slices in Lorene and Pizza format.</p>
<p>For this tutorial we will use the LS220 tables in the <code>Data/EOS/LS220</code> folder.</p>
<h2 id="initial-data">Initial Data</h2>
<p>Please refer to the <a href="http://www.lorene.obspm.fr/">Lorene</a> documentation for how to create new initial data. In this tutorial, we will use the initial data distributed together with WhiskyTHC and available in the <code>Data/LoreneData</code> folder. In particular, the instructions refer to the binary <code>LS220_T001_13641364_45km</code> (LS220 EOS; initial temperature 0.01 MeV; 1.364 Msun vs 1.364 Msun; initial separation 45 km) that can be found in the <code>Data/LoreneData/GW170817/R01</code> folder. The same procedure can be used to run other binaries.</p>
<h2 id="parameter-file">Parameter File</h2>
<p>The parameter file (parfile in brief) contains all necessary inputs to setup a binary simulation. In this tutorial we will use the example parfile <code>parfiles/bns_gw170817_ls220_q1.par</code> distributed with WhiskyTHC. The parfile starts with a list of modules (thorns in Cactus jargon) used for the simulation:</p>
<pre><code>    ActiveThorns = &quot;
        ADMBase
        ...
    &quot;</code></pre>
<p>This is followed by the parameters for each thorn. You can find more details about the parfile in the Cactus documentation <code>Cactus/doc/UserGuide.pdf</code> and in the documentation of each thorn, which can be compiled with the command (in the Cactus folder)</p>
<pre><code>    $ make [ThornName]-ThornDoc</code></pre>
<p>For example</p>
<pre><code>    $ make THC_Core-ThornDoc</code></pre>
<p>generates the file <code>doc/ThornDoc/THCCore/THC_Core/documentation.pdf</code>.</p>
<p>Here, I only give a summary of the first parameters you may want to play with.</p>
<p>This sets the maximum wallclock time:</p>
<pre><code>    TerminationTrigger::max_walltime = @WALLTIME_HOURS@</code></pre>
<p>At the end of this time the simulation will checkpoint and terminate. If you use batchtools to manage your simulations <code>@WALLTIME_HOURS@</code> will be automatically replaced when you create a segment according to the settings you have chosen for batchtools.</p>
<p>The initial data file and EOS are set by</p>
<pre><code>    PizzaIDBase::eos_file     = &quot;@HOME@/FreeTHC/Data/EOS/LS220/0.01MeV/LS_220_25-Sept-2017.pizza&quot;
    LoreneID::lorene_bns_file = &quot;@HOME@/FreeTHC/Data/LoreneData/GW170817/R01/LS220_T001_13641364_45km/resu.d&quot;</code></pre>
<p>Again note that if you are using <code>batchtools</code> <code>@HOME@</code> will be replaced by the path to your home directory and I am assuming that you checked out WhiskyTHC there. You might need to adjust these paths if you have placed WhiskyTHC and/or the Data distribution in another location.</p>
<p>The tables used for the evolution are set with the following instructions in the parfile:</p>
<pre><code>    EOS_Thermal_Table3d::eos_db_loc   = &quot;@HOME@/FreeTHC/Data/EOS&quot;
    EOS_Thermal_Table3d::eos_folder   = &quot;LS220&quot;
    EOS_Thermal_Table3d::eos_filename = &quot;LS_220_hydro_27-Sep-2014.h5&quot;
    WeakRates::table_filename         = &quot;@HOME@/FreeTHC/Data/EOS/LS220/LS_220_weak_27-Sep-2014.h5&quot;</code></pre>
<p>The last section of the parfile is dedicated to selecting the output variables and frequency and is self-explanatory.</p>
<h2 id="adaptive-mesh-refinement">Adaptive Mesh Refinement</h2>
<p>WhiskyTHC uses a simple adaptive mesh refinement (AMR) algorithm built on top of the <a href="http://carpetcode.org/">Carpet</a> infrastructure. In the Carpet AMR model, the grid is constructed starting from a uniform, low-resolution, Cartesian box. Where higher resolution is required, the user can create hierarchies of refinement levels made of nested boxes, with each level having twice the resolution of the parent level (the level number 0 is the global uniform grid). Carpet supports the creation of a number of these hierarchies that can have varying number of levels and position.</p>
<p>For NS merger simulations we use three hierarchies or regions:</p>
<ul>
<li><p>Regions 1 and 2, are centered on the position of the NS.</p></li>
<li><p>Regions 3 is centered at the coordinate origin and includes both binaries;</p></li>
</ul>
<p>The thorn <code>BNSTrackerGen</code> is responsible for moving the regions 1 and 2. During the inspiral, region 3 only has few active refinement level boxes and is used to resolve the gravitational wave (GW) “wave zone”. Once the two NSs merger, we switch off regions 1 and 2 and increase the number of levels in region 3 to create high-resolution boxes that contain the whole merger remnant.</p>
<p>The size of the computational domain in code units (G = c = Msun = 1) and the basic resolution are specified by</p>
<pre><code>    CoordBase::xmin     = -1024
    CoordBase::xmax     =  1024
    CoordBase::ymin     = -1024
    CoordBase::ymax     =  1024
    CoordBase::zmin     =  0
    CoordBase::zmax     =  1024

    CoordBase::spacing  = &quot;numcells&quot;
    CoordBase::ncells_x = 128
    CoordBase::ncells_y = 128
    CoordBase::ncells_z = 64</code></pre>
<p>For this particular run the resolution on the coarsest grid is 8 Msun = 11.81 km. For merger simulations we use 7 refinement levels (including level 0), so the resolution in the region containing the NSs is 16 / 2^6 = 0.250 Msun = 369.2 m. Note that this is half of the typical resolution we use in production simulations.</p>
<p>The AMR setup is specified by the <code>Carpet</code>, <code>CarpetRegrid2</code>, and <code>BNSTrackerGen</code> options. The only parameters you might want to experiment with, before becoming familiar with the inner working of the code, are the initial positions of the two NS on the grid:</p>
<pre><code>    CarpetRegrid2::position_x_1 = 15.2365094091
    CarpetRegrid2::position_x_2 = 15.2365094091</code></pre>
<p>and the sizes of the innermost boxes, which should be large enough to cover the entire NSs during the inspiral:</p>
<pre><code>    CarpetRegrid2::radius_1[6] = 10.0
    CarpetRegrid2::radius_2[6] = 10.0</code></pre>
<p>and the postmerger</p>
<pre><code>    CarpetRegrid2::radius_3[6] = 10.0</code></pre>
<p>Note that the leakage scheme also needs to know the position of the two NSs at the initial time. This is because the neutrino optical depth is initially estimated using two auxiliary spherical grids centered at the locations of the NSs:</p>
<pre><code>    THC_LeakageBase::center_grid1[0] = 15.2365094091
    THC_LeakageBase::center_grid2[0] = -15.2365094091</code></pre>
<h2 id="running-the-simulation">Running the Simulation</h2>
<p>On most supercomputers, simulations should be run in a special high-performance file system that we assume to be called <code>$SCRATCH</code> in this tutorial.</p>
<p>First we create a folder for the simulation:</p>
<pre><code>    $ cd $SCRATCH
    $ mkdir -p simulations/GW170817/LS220_M13641364</code></pre>
<p>The simulation folder needs to be in a nested folder several levels below <code>$SCRATCH</code>, because the Lorene initial data reader will look for the EOS table in the folder <code>../../../Eos_tables</code> (relative to the working directory of the simulation). For this reason, we also need to create a symbolic link in the <code>$SCRATCH</code> folder:</p>
<pre><code>    $ cd $SCRATCH
    $ ln -s $HOME/FreeTHC/Data/EOS/Lorene Eos_tables</code></pre>
<p>We will use batchtools to create and run the simulation, to this aim we need to add batchtools to the system <code>PATH</code>. If you are using bash, you can achieve this with</p>
<pre><code>    $ export PATH=$PATH:$HOME/FreeTHC/Cactus/batchtools/bin</code></pre>
<p>This line can be added to your <code>$HOME/.bashrc</code> file to make batchtools permanently available.</p>
<p>Now we are ready to setup the simulation</p>
<pre><code>    $ cd $SCRATCH/simulations/GW170817/LS220_M13641364
    $ batchtools init \
        --parfile ~/FreeTHC/parfiles/bns_gw170817_ls220_q1.par \
        --exe ~/FreeTHC/Cactus/exe/cactus_thc \
        --batch ~/FreeTHC/Cactus/batchtools/templates/cactus/comet.sub</code></pre>
<p>This will create a <code>BATCH</code> folder in the current directory. Please edit the <code>BATCH/CONFIG</code> file to select the allocation you plan to use and other options. See the <a href="http://bitbucket.org/dradice/batchtools">batchtools documentation</a> for more details. This example simulation require about 35 Gb of memory and should be ideally run using between 64 and 256 CPU cores.</p>
<p>The last step is to create the first simulation segment and submit it to the queue:</p>
<pre><code>    $ batchtools makesegment
    $ batchtools submit</code></pre>
<p>This will create a folder <code>output-0000</code> which will be populated with the output of the simulation once the job is run by the queue manager of the cluster (e.g., SLURM). Subsequent restarts will be named as <code>output-0001</code>, <code>output-0002</code>, etc. The separation of the run in multiple chunks is done for safety reasons, because we do not want a bad file system write to corrupt existing data and we want to be able to restart the simulation at intermediate points if it is necessary to make adjustments.</p>
<p>Depending on how many cores you have allocated for this simulation, a full run might require several restarts to reach the time of merger.</p>
<h2 id="analyze-and-visualize-results">Analyze and Visualize Results</h2>
<p>Note that WhiskyTHC uses geometrized units with G = c = Msun = 1 (but temperatures are given in MeV).</p>
<p>Some of the key diagnostics outputs from WhiskyTHC are:</p>
<ul>
<li>rho.maximum.asc : maximum density (geometrized units)</li>
<li>dens.norm1.asc : total rest-mass density (excluded regions inside event horizons) divided by the domain volume</li>
<li>H.norm2.asc : norm2 of the Hamiltonian constrain violation</li>
<li>outflow_det_X.asc : flux of material crossing the detector sphere number X (specified in the parfile)</li>
<li>mp_Psi4_l2_m2_r400.00.asc: real and imaginary part of the l = m = 2 multipole of the Weyl scalar Psi4 extracted at a radius r = 400 Msun.</li>
<li>bnstrackergen-bns_positions: output from the BNS tracker</li>
</ul>
<p>For example here we plot the curvature gravitational wave using gnuplot:</p>
<pre><code>    gnuplot&gt; plot &quot;&lt;cat output-00??/data/mp_Psi4_l2_m2_r400.00.asc&quot; u 1:($2*400) w l</code></pre>
<p>This yields the following figure:</p>
<figure>
<img src="./tutorial.png" alt="" /><figcaption>l=m=2 multipole of the Weyl scalar Psi4</figcaption>
</figure>
<p>For more advanced analysis and visualization of WhiskyTHC results, I developed two python packages:</p>
<ul>
<li><a href="https://bitbucket.org/dradice/scidata">scidata</a>: utilities to read Cactus/Carpet data into python for further analysis.</li>
<li><a href="https://bitbucket.org/dradice/scivis">scivis</a>: library of 2D visualization scripts (uses scidata).</li>
</ul>
<p>Please refer to the <code>README</code> files of those packages for instructions on how to install and use them.</p>
<p>For complex visualizations, including 3D volume renderings, you can use <a href="https://wci.llnl.gov/simulation/computer-codes/visit">VisIt</a>. VisIt already includes a <a href="http://cactuscode.org/documentation/visualization/VisIt/">plugin to read Cactus/Carpet data</a>. Note that this plugin can only handle correctly vertex-centered AMR data. It is still possible to use it to visualize cell-centered AMR data with VisIt and with this plugin, however this produces small, but visible, artifacts at refinement level boundaries.</p>