Here you will find a step by step guide to downloading, configuring, and running the Einstein Toolkit. You may use this tutorial on a workstation or laptop, or on a supported cluster. Configuring the Einstein Toolkit on an unsupported cluster is beyond the scope of this tutorial. If you find something that does not work, please feel free to mail [email protected]
On Mac, please first
sudo xcodebuild -license
sudo xcode-select --install
sudo port -N install pkgconfig gcc9 openmpi-gcc9 fftw-3 gsl jpeg zlib hdf5 +fortran +gfortran openssl ld64
sudo port select mpi openmpi-gcc9-fortran
sudo port select gcc mp-gcc9
On Debian/Ubuntu/Mint use this command (the strange syntax is to suport all three of them):
$(sudo -l sudo) su -c 'apt-get install -y subversion gcc git numactl libgsl-dev libpapi-dev python libhwloc-dev make libopenmpi-dev libhdf5-openmpi-dev libfftw3-dev libssl-dev liblapack-dev g++ curl gfortran patch pkg-config libhdf5-dev libjpeg-turbo?-dev'
On Fedora use this command:
sudo dnf install -y libjpeg-turbo-devel gcc git lapack-devel make subversion gcc-c++ which papi-devel python hwloc-devel openmpi-devel hdf5-openmpi-devel openssl-devel libtool-ltdl-devel numactl-devel gcc-gfortran findutils hdf5-devel fftw-devel patch gsl-devel pkgconfig
module load mpi/openmpi-x86_64
You will have to repeat the module load command once in each new shell each time you would like to compile or run the code. You may have to log out and back in for the module command to become visible.
On Centos use this command:
su -c 'yum install -y epel-release'
su -c 'yum install -y libjpeg-turbo-devel gcc git lapack-devel make subversion gcc-c++ which papi-devel hwloc-devel openmpi-devel hdf5-openmpi-devel openssl-devel libtool-ltdl-devel numactl-devel gcc-gfortran hdf5-devel fftw-devel patch gsl-devel'
module load mpi/openmpi-x86_64
You will have to repeat the module load command once in each new shell each time you would like to compile or run the code. You may have to log out and back in for the module command to become visible.
On OpenSuse use this command:
sudo zypper install -y curl gcc git lapack-devel make subversion gcc-c++ which papi-devel hwloc-devel openmpi-devel libopenssl-devel libnuma-devel gcc-fortran hdf5-devel libfftw3-3 patch gsl-devel pkg-config
sudo mpi-selector --set $(mpi-selector --list | head -n1)
You will only have to execute the mpi-selector once, after that log out and back in to make the mpirun and mpicc commands visible without which Cactus will compile very slowly and fail to run.
This notebook is intended to be used online on the Einstein Toolkit tutorial server, offline as a read-only document, as a jupyter notebook that you can download and also in your own docker container using ndslabs/jupyter-et. To make all of these work some setting need to be tweaked, which we do in the next cell.
# this allows you to use "cd" in cells to change directories instead of requiring "%cd"
%automagic on
# override IPython's default %%bash to not buffer all output
from IPython.core.magic import register_cell_magic
@register_cell_magic
def bash(line, cell): get_ipython().system(cell)
# Some versions of OpenMPI prevent oversubscribing cpus, which may happen if simfactory's
# number of cores detection is imperfect.
# OpenMPI by default pins MPI ranks to specific cores, which causes issues on shared
# system like the tutorial cluster.
import os
os.environ["OMPI_MCA_rmaps_base_oversubscribe"] = "true"
os.environ["OMPI_MCA_hwloc_base_binding_policy"] = "none"
A script called GetComponents is used to fetch the components of the Einstein Toolkit. GetComponents serves as convenient wrapper around lower level tools like git and svn to download the codes that make up the Einstein toolkit from their individual repositories. You may download and make it executable as follows:
Note: By default, the cells in this notebook are Python 3 commands. However, cells that start with %%bash are executed in a bash shell. If you wish to run these commands outside the notebook and in a bash shell, cut and paste only the part after the initial %%bash. The next cell tweaks some setting to make it easier to copy and paste commands from this notebook into a shell if you would like to follow along in a terminal instead.
cd ~/
%%bash
curl -kLO https://raw.githubusercontent.com/gridaphobe/CRL/ET_2019_03/GetComponents
chmod a+x GetComponents
GetComponents accepts a thorn list as an argument. To check out the needed components:
%%bash
./GetComponents --parallel https://bitbucket.org/einsteintoolkit/manifest/raw/ET_2019_03/einsteintoolkit.th
cd ~/Cactus
The recommended way to compile the Einstein Toolkit is to use the Simulation Factory ("SimFactory").
The ET depends on various libraries, and needs to interact with machine-specific queueing systems and MPI implementations. As such, it needs to be configured for a given machine. For this, it uses SimFactory. Generally, configuring SimFactory means providing an optionlist, for specifying library locations and build options, a submit script for using the batch queueing system, and a runscript, for specifying how Cactus should be run, e.g. which mpirun command to use.
%%bash
./simfactory/bin/sim setup-silent
After this step is complete you will find your machine's default setup under ./simfactory/mdb/machines/<hostname >.ini You can edit some of these settings freely, such as "description", "basedir" etc. Some entry edits could result in simulation start-up warnings and/or errors such as "ppn" (processor-per-node meaning number of cores on your machine), "num-threads" (number of threads per core) so such edits must be done with some care.
Assuming that SimFactory has been successfully set up on your machine, you should be able to build the Einstein Toolkit with the command below. The option "--mdbkey make 'make -j2'" sets the make command that will be used by the script. The number used is the number of processors used when building. Even in parallel, this step may take a while, as it compiles all the thorns specified in manifest/einsteintoolkit.th.
Note that the "cat" command on the end of the line is to prevent problems with the display in Jupyter when output comes from multiple threads.
Note: Using too many threads to compile on the test machine may result in compiler failures.
%%bash
./simfactory/bin/sim build --mdbkey make 'make -j2' --thornlist ../einsteintoolkit.th
You can now run the Einstein Toolkit with a simple test parameter file.
%%bash
./simfactory/bin/sim create-run helloworld \
--parfile arrangements/CactusExamples/HelloWorld/par/HelloWorld.par
The above command will run the simulation naming it "helloworld" and display its log output to screen.
If you see
INFO (HelloWorld): Hello World!anywhere in the above output, then congratulations, you have successfully downloaded, compiled and run the Einstein Toolkit! You may now want to try some of the other tutorials to explore some interesting physics examples.
What follows is the much more computationally intensive example of simulating a static TOV star. Just below this cell you can see the contents of a Cactus parameter file to simulate a single, spherical symmetric star using the Einstein Toolkit. The parameter file has been set up to run to completion in about 10 minutes, making it useful for a tutorial but too coarsely resolved to do science with it.
%%writefile par/tov_ET.par
# Example parameter file for a static TOV star. Everything is evolved, but
# because this is a solution to the GR and hydro equations, nothing changes
# much. What can be seen is the initial perturbation (due to numerical errors)
# ringing down (look at the density maximum), and later numerical errors
# governing the solution. Try higher resolutions to decrease this error.
# Some basic stuff
ActiveThorns = "Time MoL"
ActiveThorns = "Coordbase CartGrid3d Boundary StaticConformal"
ActiveThorns = "SymBase ADMBase TmunuBase HydroBase InitBase ADMCoupling ADMMacros"
ActiveThorns = "IOUtil"
ActiveThorns = "Formaline"
ActiveThorns = "SpaceMask CoordGauge Constants LocalReduce aeilocalinterp LoopControl"
ActiveThorns = "Carpet CarpetLib CarpetReduce CarpetRegrid2 CarpetInterp"
ActiveThorns = "CarpetIOASCII CarpetIOScalar CarpetIOHDF5 CarpetIOBasic"
# Finalize
Cactus::terminate = "time"
Cactus::cctk_final_time = 400 #800 # divide by ~203 to get ms
# Termination Trigger
ActiveThorns = "TerminationTrigger"
TerminationTrigger::max_walltime = 24 # hours
TerminationTrigger::on_remaining_walltime = 0 # minutes
TerminationTrigger::check_file_every = 512
TerminationTrigger::termination_file = "TerminationTrigger.txt"
TerminationTrigger::termination_from_file = "yes"
TerminationTrigger::create_termination_file = "yes"
# grid parameters
Carpet::domain_from_coordbase = "yes"
CartGrid3D::type = "coordbase"
CartGrid3D::domain = "full"
CartGrid3D::avoid_origin = "no"
CoordBase::xmin = 0.0
CoordBase::ymin = 0.0
CoordBase::zmin = 0.0
CoordBase::xmax = 24.0
CoordBase::ymax = 24.0
CoordBase::zmax = 24.0
# Change these parameters to change resolution. The ?max settings above
# have to be multiples of these. 'dx' is the size of one cell in x-direction.
# Making this smaller means using higher resolution, because more points will
# be used to cover the same space.
CoordBase::dx = 2.0
CoordBase::dy = 2.0
CoordBase::dz = 2.0
CarpetRegrid2::regrid_every = 0
CarpetRegrid2::num_centres = 1
CarpetRegrid2::num_levels_1 = 2
CarpetRegrid2::radius_1[1] = 12.0
CoordBase::boundary_size_x_lower = 3
CoordBase::boundary_size_y_lower = 3
CoordBase::boundary_size_z_lower = 3
CoordBase::boundary_size_x_upper = 3
CoordBase::boundary_size_y_upper = 3
CoordBase::boundary_size_z_upper = 3
CoordBase::boundary_shiftout_x_lower = 1
CoordBase::boundary_shiftout_y_lower = 1
CoordBase::boundary_shiftout_z_lower = 1
CoordBase::boundary_shiftout_x_upper = 0
CoordBase::boundary_shiftout_y_upper = 0
CoordBase::boundary_shiftout_z_upper = 0
ActiveThorns = "ReflectionSymmetry"
ReflectionSymmetry::reflection_x = "yes"
ReflectionSymmetry::reflection_y = "yes"
ReflectionSymmetry::reflection_z = "yes"
ReflectionSymmetry::avoid_origin_x = "no"
ReflectionSymmetry::avoid_origin_y = "no"
ReflectionSymmetry::avoid_origin_z = "no"
# storage and coupling
TmunuBase::stress_energy_storage = yes
TmunuBase::stress_energy_at_RHS = yes
TmunuBase::timelevels = 1
TmunuBase::prolongation_type = none
HydroBase::timelevels = 3
ADMMacros::spatial_order = 4
SpaceMask::use_mask = "yes"
Carpet::enable_all_storage = no
Carpet::use_buffer_zones = "yes"
Carpet::poison_new_timelevels = "yes"
Carpet::check_for_poison = "no"
Carpet::init_3_timelevels = no
Carpet::init_fill_timelevels = "yes"
CarpetLib::poison_new_memory = "yes"
CarpetLib::poison_value = 114
# system specific Carpet paramters
Carpet::max_refinement_levels = 10
driver::ghost_size = 3
Carpet::prolongation_order_space = 3
Carpet::prolongation_order_time = 2
# Time integration
time::dtfac = 0.25
MoL::ODE_Method = "rk4"
MoL::MoL_Intermediate_Steps = 4
MoL::MoL_Num_Scratch_Levels = 1
# check all physical variables for NaNs
# This can save you computing time, so it's not a bad idea to do this
# once in a whioe.
ActiveThorns = "NaNChecker"
NaNChecker::check_every = 16384
NaNChecker::action_if_found = "terminate" #"terminate", "just warn", "abort"
NaNChecker::check_vars = "ADMBase::metric ADMBase::lapse ADMBase::shift HydroBase::rho HydroBase::eps HydroBase::press HydroBase::vel"
# Hydro paramters
ActiveThorns = "EOS_Omni GRHydro"
HydroBase::evolution_method = "GRHydro"
GRHydro::riemann_solver = "Marquina"
GRHydro::GRHydro_eos_type = "Polytype"
GRHydro::GRHydro_eos_table = "2D_Polytrope"
GRHydro::recon_method = "ppm"
GRHydro::GRHydro_stencil = 3
GRHydro::bound = "none"
GRHydro::rho_abs_min = 1.e-10
# Parameter controlling finite difference order of the Christoffel symbols in GRHydro
GRHydro::sources_spatial_order = 4
# Curvature evolution parameters
ActiveThorns = "GenericFD NewRad"
ActiveThorns = "ML_BSSN ML_BSSN_Helper"
ADMBase::evolution_method = "ML_BSSN"
ADMBase::lapse_evolution_method = "ML_BSSN"
ADMBase::shift_evolution_method = "ML_BSSN"
ADMBase::dtlapse_evolution_method= "ML_BSSN"
ADMBase::dtshift_evolution_method= "ML_BSSN"
ML_BSSN::timelevels = 3
ML_BSSN::harmonicN = 1 # 1+log
ML_BSSN::harmonicF = 2.0 # 1+log
ML_BSSN::evolveA = 1
ML_BSSN::evolveB = 1
# NOTE: The following parameters select geodesic slicing. This choice only enables you to evolve stationary spacetimes.
# They will not allow you to simulate a collapsing TOV star.
ML_BSSN::ShiftGammaCoeff = 0.0
ML_BSSN::AlphaDriver = 0.0
ML_BSSN::BetaDriver = 0.0
ML_BSSN::advectLapse = 0
ML_BSSN::advectShift = 0
ML_BSSN::MinimumLapse = 1.0e-8
ML_BSSN::my_initial_boundary_condition = "extrapolate-gammas"
ML_BSSN::my_rhs_boundary_condition = "NewRad"
# Some dissipation to get rid of high-frequency noise
ActiveThorns = "SphericalSurface Dissipation"
Dissipation::verbose = "no"
Dissipation::epsdis = 0.01
Dissipation::vars = "
ML_BSSN::ML_log_confac
ML_BSSN::ML_metric
ML_BSSN::ML_curv
ML_BSSN::ML_trace_curv
ML_BSSN::ML_Gamma
ML_BSSN::ML_lapse
ML_BSSN::ML_shift
"
# init parameters
InitBase::initial_data_setup_method = "init_some_levels"
# Use TOV as initial data
ActiveThorns = "TOVSolver"
HydroBase::initial_hydro = "tov"
ADMBase::initial_data = "tov"
ADMBase::initial_lapse = "tov"
ADMBase::initial_shift = "tov"
ADMBase::initial_dtlapse = "zero"
ADMBase::initial_dtshift = "zero"
# Parameters for initial star
TOVSolver::TOV_Rho_Central[0] = 1.28e-3
TOVSolver::TOV_Gamma = 2
TOVSolver::TOV_K = 100
# Set equation of state for evolution
EOS_Omni::poly_gamma = 2
EOS_Omni::poly_k = 100
EOS_Omni::gl_gamma = 2
EOS_Omni::gl_k = 100
# I/O
# Use (create if necessary) an output directory named like the
# parameter file (minus the .par)
IO::out_dir = ${parfile}
# Write one file overall per output (variable/group)
# In production runs, comment this or set to "proc" to get one file
# per MPI process
IO::out_mode = "onefile"
# Some screen output
IOBasic::outInfo_every = 512
IOBasic::outInfo_vars = "Carpet::physical_time_per_hour HydroBase::rho{reductions='maximum'}"
# Scalar output
IOScalar::outScalar_every = 512
IOScalar::one_file_per_group = "yes"
IOScalar::outScalar_reductions = "norm1 norm2 norm_inf sum maximum minimum"
IOScalar::outScalar_vars = "
HydroBase::rho{reductions='maximum'}
HydroBase::press{reductions='maximum'}
HydroBase::eps{reductions='minimum maximum'}
HydroBase::vel{reductions='minimum maximum'}
HydroBase::w_lorentz{reductions='minimum maximum'}
ADMBase::lapse{reductions='minimum maximum'}
ADMBase::shift{reductions='minimum maximum'}
ML_BSSN::ML_Ham{reductions='norm1 norm2 maximum minimum norm_inf'}
ML_BSSN::ML_mom{reductions='norm1 norm2 maximum minimum norm_inf'}
GRHydro::dens{reductions='minimum maximum sum'}
Carpet::timing{reductions='average'}
"
# 1D ASCII output. Disable for production runs!
IOASCII::out1D_every = 2048
IOASCII::one_file_per_group = yes
IOASCII::output_symmetry_points = no
IOASCII::out1D_vars = "
HydroBase::rho
HydroBase::press
HydroBase::eps
HydroBase::vel
ADMBase::lapse
ADMBase::metric
ADMBase::curv
ML_BSSN::ML_Ham
ML_BSSN::ML_mom
"
# 2D HDF5 output
CarpetIOHDF5::output_buffer_points = "no"
CarpetIOHDF5::out2D_every = 2048
CarpetIOHDF5::out2D_vars = "
HydroBase::rho
HydroBase::eps
HydroBase::vel
HydroBase::w_lorentz
ADMBase::lapse
ADMBase::shift
ADMBase::metric
ML_BSSN::ML_Ham
ML_BSSN::ML_mom
"
%%bash
# start simulation
./simfactory/bin/sim create-submit tov_ET \
--parfile=par/tov_ET.par --procs=2 --num-threads=1 --walltime=0:20:0
The submit command submitted the simulation to the queueing system to run in the background asking for a maximum runtime of 20 minutes. On your laptop it will start right away, on a cluster the queuing system will wait until a sufficient number of nodes is able to start your simulation.
You can check the status of the simulation with the command below. You can run this command repeatedly until the job shows
[ACTIVE (FINISHED)...as it's state. Prior to that, it may show up as QUEUED or RUNNING.
%%bash
simfactory/bin/sim list-simulations tov_ET
To watch its output use the show-output command of simfactory. Interrupt the kernel (or press CTRL-C if copying & pasting these commands to a terminal) once you wish to stop watching or once the simulation has ended.
%%bash
# watch log output, following along as new output is produced
./simfactory/bin/sim show-output --follow tov_ET
You can leave out the --follow option if you would like to see all output up to this point.
Since the submit command was used to start the simulation, it is running in the background and you have to use simfactory commands to interact with it. The next cell shows how to list simulations.
Remember that you have to interrupt the kernel to stop show-output and be able to execute the cells below.
%%bash
simfactory/bin/sim list-simulations
Simfactory offers a stop command to abort a running simulation. The next cell has the command intentionally commented out to prevent accidental stopping of your very first simulation.
%%bash
#simfactory/bin/sim stop tov_ET
after this the simulation changes to the "FINISHED" state indicating it is no longer running.
The following commands will generate a simple line plot of the data. They will work in a python script as easily as they do in the notebook (just remove the "%matplotlib inline" directive).
# This cell enables inline plotting in the notebook
%matplotlib inline
import matplotlib
import numpy as np
import matplotlib.pyplot as plt
%%bash
# find out the directory that Cactus stored its output in
simfactory/bin/sim get-output-dir tov_ET
# store output directory location for later use
# in ipython you can also use this:
# outdir = ! simfactory/bin/sim get-output-dir tov_ET
# outdir = outdir[0]
import os
outdir = os.popen("simfactory/bin/sim get-output-dir tov_ET").read().rstrip("\n")
print ("Output is written to", outdir)
Numpy has a routine called genfromtxt() which is an extremely efficient reader of textual arrays of floating point numbers. This is well-suited to Cactus .asc files.
# format of the path: SIMULATION-NAME/output-NNNN/PARFILE-NAME
lin_data = np.genfromtxt(outdir+"/tov_ET/hydrobase-rho.maximum.asc")
This is all you need to do to plot the data once you've loaded it. Note, this uses Python array notation to grab columns 1 and 2 of the data file.
plt.plot(lin_data[:,1],lin_data[:,2]/lin_data[0,2], label="central density")
plt.xlabel(r'$t$ [$M_{\odot}$]');
plt.ylabel(r'$\rho_c / \rho_c(0)$');
plt.legend();
# create small dataset to show what plot should look like
def sparsify(lin_data, sparsity):
# drop unwanted datapoint
sparse_data = lin_data[::sparsity,:]
# compute min, max of dataset then difference to minimum and quantize to 8 bit precisison
minval = np.amin(sparse_data[:,2])
maxval = np.amax(sparse_data[:,2])
diff = sparse_data[:,2] - minval
quant_diff = np.minimum(np.maximum(np.round(diff / (maxval - minval) * 255.5), 0), 255).astype('int')
# timesteps are equidistant and start at 0 so we only need the stepsize
delta_t = sparse_data[1,1] - sparse_data[0,1]
# string rep of 8bit differences
quant_diff_s = ""
for i in quant_diff: quant_diff_s += "%02x" % i
print ('"%s", %.3e, %.3e, %.3e' % (quant_diff_s, minval, maxval, delta_t))
# create a low fidelity representation of every 10th datapoint and output all data a string
sparsify(lin_data, 10)
# this cell shows the expected plot using previously stored data
# reconstruct plot data from saved strings
(quant_diff_s, minval, maxval, delta_t) = \
("ff8baee2e5d2ac70320c0007182c404f5b656f7b8897a8bbcddde8eeede8ddcfc0b0a29589817b777473757a8189929ca6b0bac4cbd0d3d4d4d2cfcbc7c2bdb8b4b0adaaa9a8a9abaeb3b8bcc1c5c8cccf",
1.235e-03, 1.280e-03, 5.000e+00)
quant_diff = np.array(bytearray.fromhex(quant_diff_s))
rec_vals = quant_diff / 255. * (maxval- minval) + minval
rec_time = np.arange(0,len(quant_diff)) * delta_t
# plot them, including your results if you have them
plt.plot(rec_time, rec_vals/rec_vals[0],
label="central density (stored values)")
try: plt.plot(lin_data[:,1],lin_data[:,2]/lin_data[0,2], label="central density (your results)")
except: pass
plt.xlabel(r'$t$ [$M_{\odot}$]');
plt.ylabel(r'$\rho_c / \rho_c(0)$');
plt.legend(loc='lower right');