RAMSES Fundamentals
1. Prerequisites
You will need a fortran compiler to compile RAMSES, along with some version of MPI if you have access to multiple cores and want to run RAMSES in parallel mode. The tutorials require running python3 in a Jupyter notebook, and will use different python packages and analysis tools (see the prerequisites section of each tutorial). This ‘Fundamentals’ tutorial specifically requires the Osyris package for visualization of RAMSES results.
For getting all this set up, we recommend that you start by following the instructions in general setup requirements.
2. First steps with RAMSES
RAMSES (Teyssier 2002) is a popular code used to produce astrophysical simulations, using the adaptive mesh refinement (AMR) numerical technique, together with state-of-the art star formation and feedback models. Each directory contains a set of files with a given common purpose. For example, amr/ contains all Fortran 90 routines dealing with the AMR grid management and MPI communications, while hydro/ contains all Fortran 90 routines dealing with hydrodynamics. RAMSES is written in Fortran 90, which can be used on massively parallel architectures thanks to the MPI library (MPI stands for Message Passing Interface). Therefore, before producing RAMSES simulations, we need to compile the RAMSES code, as explained below.
In order to obtain and compile the code, you can follow the first steps of the documentation at this link, summarised (and updated) below.
2.1 Downloading the code
To get the RAMSES code, open a terminal, and clone the code into the Fundamentals sub-directory:
[ ]:
%%bash
git clone https://github.com/ramses-organisation/ramses
2.2 Compiling the code: the Makefile
The first directory you are interested in is the bin/ directory, in which the code will be compiled thanks to the set of instructions written in the Makefile.
[ ]:
%%bash
cd ramses/bin
cat Makefile
The first lines of the Makefile correspond to the so-called Compilation Time Parameters. They are used to e.g.:
set the type of compiler to be used, the name of the RAMSES executable
define the number of dimensions and the total number of variables needed for the simulation
enable or disable the compilation of some of the physical modules, and of the MPI routines
Other preprocessor directives are then defined below, in the variable DEFINES.
Based on the choice of compiler (COMPILER=GNU or INTEL), Fortran compiler options and directives are defined. The F90 variable invokes the Fortran compiler (also depends on the use of MPI), while FFLAGS contains other compilation flags and preprocessor directives (also depends on the use of debug options). While both F90 and FFLAGS can be set up using default values, one might need or want to edit them according to the compiler, the machine, or the simulation. In this
case, you can manually edit F90 and FFLAGS in the Makefile.
Finally, the rest of the Makefile contains the list and path of the files to be compiled.
Once the Makefile is ready, you can (re)compile the code:
[ ]:
%%bash
cd ramses/bin
make clean
make
2.3 Setting up the parameters: the namelist
To run a simulation, RAMSES needs a parameter file, which is a Fortran namelist (brief intro here, and more details about the RAMSES namelist here). Namelists can be read and written as text files, usually ending with .nml. They are used to assign values to RAMSES variables, to be used instead of the default values initialised in the code.
Namelist files are organised in namelist blocks, or sections (that may be provided in any order) that define parameters related to different parts of the code. Each block starts with &BLOCK_NAME and ends with the character /. Values are then declared in the corresponding block using standard Fortran namelist syntax, i.e. as parameter_name = value (where values for character variables must be enclosed in quotes, and logical values must be specified using .true. or .false.).
The following 4 parameter blocks are mandatory and must always be present in the namelist file:
RUN_PARAMS: defines some of the most basic properties of the runAMR_PARAMS: defines the global mesh properties and memory footprint of the codeINIT_PARAMS: tells RAMSES where to look for initial condition filesOUTPUT_PARAMS: defines the output strategy
The other blocks are optional. They must be present in the file only if they are relevant to the current run:
POISSON_PARAMS: specifies runtime parameters for the Poisson solverBOUNDARY_PARAMS: sets up boundary conditions for the current simulationUNIT_PARAMS: contains the parameters related to units that are set “by hand” for non cosmological runsCOOLING_PARAMS: contains parameters related to cooling / basic chemistrySF_PARAMS: contains the parameters related to star formationFEEDBACK_PARAMS: contains the parameters related to stellar feedbackHYDRO_PARAMS: specifies runtime parameters for the Godunov solverREFINE_PARAMS: defines the refinement strategySINK_PARAMS: contains the parameters related to the sink particle implementation, for star formation or black holesCLUMPFIND_PARAMS: contains the parameters related to the built-in RAMSES clump finder (PHEW)RT_PARAMSandRT_GROUPS: sets radiative properties for RAMSES RHD runsMOVIE_PARAMS: contains the parameters related to the movies, for creating on the fly visualizations of a simulation
The variables you can set or adjust in each namelist block are described in more details in the the RAMSES namelist documentation and in the different tutorials.
3. Sod shock tube
This standard one-dimensional hydro experiment is part of the automatic RAMSES test suite and is often used in the literature to demonstrate the accuracy of hydrodynamics codes (see e.g. the original RAMSES paper by Romain Teyssier). It starts with a box split in two uniform halves with a large difference in density and pressure and follows the flow of gas from the high-pressure half into the low-pressure one.
The namelist blocks are:
RUN_PARAMSstates that this is a pure-hydro run and allows sub-cycling except on the three coarsest levels,AMR_PARAMSallows refinement levels 3 to 10, i.e. effetive resolutions of \(2^3=8\) up to \(2^{10}=1024\) cells ,BOUNDARY_PARAMS, sets reflective boundaries on both sides,INIT_PARAMSsets up a high-pressure high-density half and a low-pressure low-density half,OUTPUT_PARAMSsets one output for the run, in addition to the initial conditions output, and sets at what time the simulation makes that output and stops,HYDRO_PARAMSsets the adiabatic index and some details for the hydro solverREFINE_PARAMSsets the refinement strategy, in this case on gradients in density, velocity, and pressure
Now let’s run the test. First compile a 1d RAMSES executable:
[ ]:
%%bash
cd ramses/bin
make clean
make NDIM=1 PATCH= SOLVER=hydro EXEC=ramses_sod
Then run the test, using the sod-tube.nml namelist:
[ ]:
%%bash
cd sod-tube
../ramses/bin/ramses_sod1d sod-tube.nml
You should see RAMSES start and spit out a log of its progress in the terminal. Once it is finished (this should be less than a second), you’ll see two output directories, named output_00001 and output_00002.
Feel free to look at what is in those output directories. They contain .out files, for the amr grid structure and the hydro quantities in each cell, amr_XXXXX.outYYYYY and hydro_XXXXX.outYYYYY respectively. Here, the XXXXX stands for the number of the output, and the YYYYY for the number of the processor making the output (since you are running in non-parallel mode, there is only one). There are also several text files containing various kinds of information about the run, such
as the amount of particles of various types (in header, but there are zero particles in this run), the hydro variables included in the run (in hydro_file_descriptor), various run properties and units (in info), and a list over how the computation is split time-wise into the different components of RAMSES (in timer).
After the run is done, you can use the two cells below to plot the final result and compare to the analytic solution.
[ ]:
# Code for plotting analytic Sod solution, taken from Thomas Chuna: https://github.com/chunatho/SodShockSolution
# Just run this cell to compile and then proceed to the next cell to read the RAMSES output and compare to the analytic solution.
import numpy as np
from scipy.optimize import newton
def f(P, pL, pR, cL, cR, gg):
"""
Function to find the roots of f!
"""
a = (gg-1)*(cR/cL)*(P-1)
b = np.sqrt( 2*gg*(2*gg + (gg+1)*(P-1) ) )
return P - pL/pR*( 1 - a/b )**(2.*gg/(gg-1.))
def SodShockAnalytic(rL, uL, pL, rR, uR, pR, xs, x0, T, gg):
"""
Analtyic solution to Sod Shock
Parameters
----------
rL, uL, pL, rR, uR, pR : Initial conditions of the Reimann problem
xs: position array (e.g. xs = [0,dx,2*dx,...,(Nx-1)*dx])
x0: THIS IS AN INDEX! the array index where the interface sits.
T: the desired solution time
gg: adiabatic constant 1.4=7/5 for a 3D diatomic gas
"""
dx = xs[1]
Nx = len(xs)
v_analytic = np.zeros((3,Nx),dtype='float64')
# compute speed of sound
cL = np.sqrt(gg*pL/rL)
cR = np.sqrt(gg*pR/rR)
# compute P
P = newton(f, 0.5, args=(pL, pR, cL, cR, gg), tol=1e-12)
# compute region positions right to left
# region R
c_shock = uR + cR*np.sqrt( (gg-1+P*(gg+1)) / (2*gg) )
x_shock = x0 + int(np.floor(c_shock*T/dx))
v_analytic[0,x_shock-1:] = rR
v_analytic[1,x_shock-1:] = uR
v_analytic[2,x_shock-1:] = pR
# region 2
alpha = (gg+1)/(gg-1)
c_contact = uL + 2*cL/(gg-1)*( 1-(P*pR/pL)**((gg-1.)/2/gg) )
x_contact = x0 + int(np.floor(c_contact*T/dx))
v_analytic[0,x_contact:x_shock-1] = (1 + alpha*P)/(alpha+P)*rR
v_analytic[1,x_contact:x_shock-1] = c_contact
v_analytic[2,x_contact:x_shock-1] = P*pR
# region 3
r3 = rL*(P*pR/pL)**(1/gg)
p3 = P*pR
c_fanright = c_contact - np.sqrt(gg*p3/r3)
x_fanright = x0 + int(np.ceil(c_fanright*T/dx))
v_analytic[0,x_fanright:x_contact] = r3
v_analytic[1,x_fanright:x_contact] = c_contact
v_analytic[2,x_fanright:x_contact] = P*pR
# region 4
c_fanleft = -cL
x_fanleft = x0 + int(np.ceil(c_fanleft*T/dx))
u4 = 2 / (gg+1) * (cL + (xs[x_fanleft:x_fanright]-xs[x0])/T )
v_analytic[0,x_fanleft:x_fanright] = rL*(1 - (gg-1)/2.*u4/cL)**(2/(gg-1))
v_analytic[1,x_fanleft:x_fanright] = u4
v_analytic[2,x_fanleft:x_fanright] = pL*(1 - (gg-1)/2.*u4/cL)**(2*gg/(gg-1))
# region L
v_analytic[0,:x_fanleft] = rL
v_analytic[1,:x_fanleft] = uL
v_analytic[2,:x_fanleft] = pL
return v_analytic
[ ]:
# Load RAMSES output with Osyris, plot the results, and compare to the analytic solution
import osyris
import matplotlib.pyplot as plt
snapshot = 2 # Final output of the simulation -- change to 1 to see initial conditions
data = osyris.RamsesDataset(snapshot, path="./sod-tube/").load()
time_sec = data.meta["time"].to("sec").magnitude
gamma = data.meta["gamma"]
x = (data["mesh"]["position_x"]).to("cm").values
order = x.argsort()
x = x[order]
rho = data["mesh"]["density"][order]
vx = data["mesh"]["velocity_x"][order]
p = data["mesh"]["pressure"][order]
amrlev = data["mesh"]["level"][order]
# We can plot directly with Osyris, but here we use our own plots for more flexibility
fig, ax = plt.subplots(nrows=1, ncols=3, figsize=(12, 3), constrained_layout=True)
ax[0].plot(x,rho,'o',color='black',markerfacecolor='none')
ax[0].set_xlabel('Distance (cm)')
ax[0].set_ylabel(r'$\rm Density \ (g/cm^3)$')
ax[1].plot(x,vx,'o',color='black',markerfacecolor='none')
ax[1].set_xlabel('Distance (cm)')
ax[1].set_ylabel('Velocity (cm/s)')
ax[2].plot(x,p,'o',color='black',markerfacecolor='none',label="RAMSES at t=%.3f s"%(time_sec))
ax[2].set_xlabel('Distance (cm)')
ax[2].set_ylabel(r'$\rm Pressure \ (erg/cm^3)$')
axAMR = ax[2].twinx()
axAMR.plot(x,amrlev,color='black',ls='dotted')
axAMR.set_ylabel('AMR Level')
ax[2].plot([1],[1],color='black',ls='dotted',label="AMR level") # just for the legend
# Overplot analytic result
rL, uL, pL = 1.0, 0.0, 1
rR, uR, pR = 0.125, 0.0, .1
Nx = 1000
xs = np.linspace(0,1,Nx)
analytic = SodShockAnalytic(rL, uL, pL, rR, uR, pR, xs, Nx//2, time_sec, gamma)
ax[0].plot(xs,analytic[0])
ax[1].plot(xs,analytic[1])
ax[2].plot(xs,analytic[2],label="Analytic result")
ax[2].legend()
Before going on to the next section, feel free to play a bit with some of the namelist parameters, e.g. changing the hydro solver in HYDRO_PARAMS or the refinement strategy in REFINE_PARAMS, and see how your solution is affected.
4. A Sedov-Taylor blast in 3D
We now turn to 3d experiments, starting with a very simple setup of a thermal blast corresponding to one type-II supernova, i.e. \(E = 10^{51}\) ergs, in a homogeneous medium.
First compile a 3d version of the code, using MPI:
[ ]:
%%bash
cd ramses/bin
make clean
make NDIM=3 PATCH= SOLVER=hydro MPI=1 EXEC=ramses_sedov
Then run the test on four cores, using the sedov3d.nml namelist. You can run this right here in the notebook, but it is much better to copy-paste the lines below to the terminal, so you don’t have to wait for the run to finish before starting to look at the outputs in the cells further down.
[ ]:
%%bash
cd sedov
mpirun -n 4 ../ramses/bin/ramses_sedov3d sedov3d.nml
After the run is done, or even while it is running, you can use the cells below to plot the final result and compare to an analytic solution for the temperature versus density.
For this experiment setup, i.e. an initially homogeneous environment and no radiative cooling, an analytic solution exists for the evolution of the blast profile. The next cell computes this analytic profile and needs to be run only once to compare to the RAMSES run.
[ ]:
# This cell contains a routine to plot analytic Sedov profiles in N dimensions.
# Just run the cell to compile it and then go to the cell below to use the routine.
##########################################################################
# The routine returns 4 arrays of unknown size (determined
# automatically by the routine). The output variables are
# the following dimensionless quantities: r (position from the
# point like explosion, d (density), u (velocity) and p
# (pressure). To recover the true value, you have to rescale
# these dimensionless values to the true values, defining first
# the total energy E_0, the initial mass density rho_0 and the
# time t you consider and finally computing the true values
# using the following scaling laws:
#
# r = r * (E_0/rho_0)^(1./(dim+2.)) * t^(2./(dim+2.))
# d = d * rho_0
# u = u * (E_0/rho_0)^(1./(dim+2.)) * t^(-dim/(dim+2.))
# p = p * (E_0/rho_0)^(2./(dim+2.)) * t^(-2.*dim/(dim+2.)) * rho_0
#
# This routine is adapted from an IDL routine by R Teyssier
def sedovana(gamma=1.4,dim=1):
import numpy as np
#---------------------------------------------------------------------
n = int(dim)
print('gamma=',gamma)
g=float(gamma)
n=float(n)
n1=1000 ; n2=1000
vmax=4.e0/(n+2.)/(g+1.)
vmin=2.e0/(n+2.)/g
temp = np.linspace(0,n1-1,n1)
v=vmin+10.e0**(-10.e0*(1.0e0-(temp+1)/float(n1)))*(vmax-vmin)
a2 = (1.-g)/(2.*(g-1.)+n)
a1 = (n+2.)*g/(2.+n*(g-1.)) * ( 2.*n*(2.-g)/g/(n+2.)**2 - a2 )
a3 = n/(2.*(g-1)+n)
a4 = a1*(n+2.)/(2.-g)
a5 = 2./(g-2.)
a6 = g/(2.*(g-1.)+n)
a7 = a1*(2.+n*(g-1.))/(n*(2.-g))
r1 = ((n+2.)*(g+1.)/4.*v)**(-2./(2.+n)) \
*((g+1.)/(g-1.) *( (n+2.)*g/2.*v-1.) )**(-a2) \
*( (n+2.)*(g+1.) / ( (n+2)*(g+1)-2.*(2.+n*(g-1.)) ) * (1.-(2.+n*(g-1.))/2.*v) )**(-a1)
u1 = (n+2.)*(g+1.)/4.*v*r1
d1 = ((g+1.)/(g-1.)*((n+2.)*g/2.*v-1.))**(a3) \
*((g+1.)/(g-1.)*(1.-(n+2.)/2.*v ))**(a5) \
*((n+2.)*(g+1.)/( (n+2)*(g+1)-2.*(2.+n*(g-1.))) *(1.-(2.+n*(g-1.))/2.*v) )**(a4)
p1 = ((n+2.)*(g+1.)/4.*v)**(2.*n/(2.+n)) \
*((g+1.)/(g-1.)*(1.-(n+2.)/2.*v ))**(a5+1.) \
*((n+2.)*(g+1.)/( (n+2)*(g+1)-2.*(2.+n*(g-1.))) *(1.-(2.+n*(g-1.))/2.*v) )**(a4-2.*a1)
temp=np.linspace(0,n2-1,n2)
r2=r1[0]*(temp+0.5e0)/float(n2)
u2=u1[0]*r2/r1[0]
d2=d1[0]*(r2/r1[0])**(n/(g-1.0e0))
p2=p1[0]*(r2/r2)
r=np.zeros(len(r1)+len(r2)+2)
r[0:len(r2)]=r2
r[len(r2):len(r1)+len(r2)]=r1
r[len(r1)+len(r2)]=r1.max()
r[len(r1)+len(r2)+1]=r1.max()+1000.
d=np.zeros(len(r))
d[:]=r[:]
d[0:len(r2)]=d2
d[len(r2):len(r1)+len(r2)]=d1
d[len(r1)+len(r2)]=1./((g+1.)/(g-1.))
d[len(r1)+len(r2)+1]=1./((g+1.)/(g-1.))
u=np.zeros(len(r))
u[:]=r[:]
u[0:len(r2)]=u2
u[len(r2):len(r1)+len(r2)]=u1
u[len(r1)+len(r2)]=0.
u[len(r1)+len(r2)+1]=0.
p=np.zeros(len(r))
p[:]=r[:]
p[0:len(r2)]=p2
p[len(r2):len(r1)+len(r2)]=p1
p[len(r1)+len(r2)]=0.
p[len(r1)+len(r2)+1]=0.
d=d*(g+1.)/(g-1.)
u=u*4./(n+2.)/(g+1.)
p=p*8./(n+2.)**2./(g+1.)
nn=len(r)
vol=np.zeros(nn)
vol[:]=r[:]
for i in range(1,nn):
vol[i]=r[i]**n-r[i-1]**n
vol[0]=r[0]**n
const=1.
if n == 1.: const=2.0e0
if n == 2.: const=3.1415926535897931
if n == 3.: const=4.*3.1415926535897931/3.
vol=vol*const
int1=(d*u*u/2.e0)*vol
int2=p/(g-1.e0)*vol
sum1=sum(int1)
sum2=sum(int2)
summ=sum1+sum2
print('chi0=',summ**(-1./(2.+n)))
chi0=summ**(-1./(2.+n))
r=r*chi0
u=u*chi0
p=p*chi0**2
return r,d,u,p
Now make a map of the simulated Sedov explosion, using the Osyris package:
[ ]:
import osyris
import matplotlib.pyplot as plt
import numpy as np
snap=10 # You can change this!
mydata = osyris.RamsesDataset(snap, path="./sedov/").load()
print(mydata)
kpc = osyris.units("kpc")
center=osyris.Vector(0.5,y=0.5,z=0.5,unit=kpc)
# Map of density and refinement level
osyris.map(mydata["mesh"].layer("density"), norm="log", direction="z", cmap="RdGy_r",dx=1*kpc,origin=center)
time_myr = mydata.meta["time"].to("Myr").magnitude
time_sec = mydata.meta["time"].to("sec").magnitude
plt.annotate('%.2f Myr'%(time_myr), xy=(0.15, 0.8),xycoords='figure fraction')
egy_tot = np.sum(mydata['mesh']['pressure']*mydata['mesh']['dx']**3)/0.4
print('The total energy is %.3e %s'%(egy_tot.values, egy_tot.label))
# Map of refinement level
osyris.map(mydata["mesh"].layer("level"), cmap="viridis" ,dx=1*kpc,origin=center)
[ ]:
# Compute analytic profile corresponding to the simulation time, time_sec
import matplotlib.pyplot as plt
dim=3.
r,d,u,p=sedovana(dim=dim)
n = len(r)
r=r[0:n-1]
d=d[0:n-1]
u=u[0:n-1]
p=p[0:n-1]
E_0=1e51
rho_0 = 0.1 * 1.66e-24
r_cgs = r * (E_0/rho_0)**(1./(dim+2.)) * time_sec**(2./(dim+2.))
d_cgs = d * rho_0
u_cgs = u * (E_0/rho_0)**(1./(dim+2.)) * time_sec**(-dim/(dim+2.))
p_cgs = p * (E_0/rho_0)**(2./(dim+2.)) * time_sec**(-2.*dim/(dim+2.)) * rho_0
plt.figure()
plt.yscale('log')
plt.xlabel(r'${\rm r \ [kpc]}$')
plt.ylabel(r'${\rho \ [g \ cm^{-3}]}$')
plt.plot(r_cgs/3.08e21,d_cgs)
plt.show()
[ ]:
# Compare analytic and simulation profiles
import matplotlib.pyplot as plt
plt.figure()
plt.yscale('log')
plt.xlabel(r'${\rm r \ [kpc]}$')
plt.ylabel(r'${\rho \ [g \ cm^{-3}]}$')
plt.ylim(1e-30,3e-24)
plt.plot(r_cgs/3.08e21,d_cgs, color='blue', label='Analytic profile')
rad = (mydata["mesh"]["position"] - center).to("kpc").norm.values
plt.scatter(rad, mydata['mesh']['density'].values,s=0.1, color='red', label='Simulation cells')
plt.legend()
plt.show()
You can now try to change the two cells above to compare the simulated pressure and velocity to the analytic solution (instead of density).
Generally, the comparison you should see is OK in terms of the position of the expanding shock blast, but the simulated profile does not really fit the analytic solution perfectly. What do you think are the reasons for that and what could be done in the simulation to improve the comparison?
4.1 Diverging from the analytic solution
In reality, even the position of the blast in the Sedov-Taylor solution is not a very accurate description of the late evolution of SN blasts. In part this is because the environment is usually far from being homogeneous, but perhaps more importantly there is the effect of radiative cooling.
Try to run a new Sedov-Taylor blast experiment where you turn on radiative cooling in the namelist. What happens to your blast and how does it now compare to the analytic solution?
You can then try to change to Solar metallicity by setting z_ave = 1 in the namelist, i.e. assume a Solar metallicity for all gas. What happens now and why?
4.2 Code units (optional)
This exercise is for understanding how to convert between the RAMSES code units and physical units. The exercise is optional but probably very useful, since if you go on to seriously use RAMSES, you’ll find that a large part of your time, frustration, and confusion will be due to unit conversions.
The namelist defines two regions, in INIT_PARAMS; a background ‘square’ low-pressure region of 0.1 hydrogen atom per cubic centimetre, defined via the d_region parameter, and a ‘point’ region where energy is injected via the p_region parameter. Note that for the ‘square’ region, p_region is interpreted as a pressure, i.e. energy per unit volume, but for the ‘point’ region it is just an energy (i.e. not per volume). The d_region and p_region values are given in code
units, and the exercise is to see how they convert to more understandable units like CGS.
In RAMSES this conversion is made with the three base units, which are also given in the namelist. Each of these converts a RAMSES variable, in code units, to cgs units. In other words, you multiply a density given by RAMSES, with units_density to convert to \(\rm{ g /cm^3}\), you multiply a length in RAMSES with units_length to convert to \(\rm cm\) and a time with units_time to convert to seconds. For example, you can find that the box length of 1 given in code units in the
namelist corresponds to 1 kpc and that the initial homogeneous density given by d_region corresponds to \(1.66 \times 10^{-25} \ {\rm g / cm^3}\). Any unit conversion can be derived from these three basic units. For example, speed is converted from code units to \(\rm cm/s\) with
I now leave it to you to derive how the value of p_region\(= 8.6 \times 10^{-7}\) in the Sedov namelist corresponds to \(10^{51}\) ergs. Keep in mind that p_region indicates a pressure \(p\) and you want to derive an energy \(e\), and the two are related via
where \(\gamma= 1.4\) is the ratio of specific heats, also given in the namelist.
4.3 Exploration
You can now play with the Sedov blast run if you like, e.g. change the resolution, refinement strategy, or the hydro solver. You can also experiment with running the same experiment in 1d or 2d (for which you can also compare to analytic solutions).
5. A Strömgren circle in 2D
Now we go for a bit of radiative transfer (RT) and radiation hydrodynamics (RHD), using the M1 RT method in RAMSES-RT.
The Strömgren sphere is one of the few radiative transfer tests which has a simple analytic result to compare to. In two dimensions, in the limit of a constant radiation source emitting \(\dot N\) photons per second into a stationary and homogeneous pure hydrogen gas with number density \(\rm n_H\), a so-called Strömgren circle (sphere in 3d) of ionised hydrogen (HII) forms, initially expands, and eventually settles at the Strömgren radius, \(r_S\), where the rate of recombinations inside the circle equals the source luminosity, i.e.
where \(\alpha^{\rm B}\) is the case-B recombination rate of the gas. Solving for the 2d Strömgren radius, we get
The propagation of the ionization front (I-front) with time t can easily be derived to have the form
where \(t_{\rm rec} = \left( n_{\rm H} \alpha^{\rm B} \right)^{-1}\) is the recombination time of the gas. (Pawlik & Schaye 2008, showed that the analytic expression above for \(r_{\rm S}\) is actually too small by 5%, due to the approximation of a step function in the ionised fraction at the Strömgren radius.)
5.1 Compare to the analytic Strömgren radius, with a reduced light speed
First, compile RAMSES-RT:
[ ]:
%%bash
cd ramses/bin
make clean
make NDIM=2 SOLVER=hydro RT=1 NGROUPS=1 NIONS=3 MPI=1 EXEC=ramses_rhd
Then run the Strömgren circle experiment:
[ ]:
%%bash
cd stromgren
mpirun -np 4 ../ramses/bin/ramses_rhd2d stromgren2d.nml
The namelist (stromgren2d.nml) sets up an isotropic source of radiation, with a luminosity of \(\dot N = 2 \times 10^{28}\) ionizing photons per second, in the corner of a 1 kpc wide cubic box. The homogeneous gas density is set to \(n_{\rm H} = 10^{-1} \ {\rm cm}^{-2}\). In the namelist, the rt_TConst=1d4 parameter forces the temperature everywhere to \(10^4\) K, which sets a constant recombination rate \(\alpha^{\rm B} = 2.6 \times 10^{-13} \ \rm s^{-1} \rm cm^{2}\). From
the equation above, we then get \(r_{\rm S} = 0.51\) kpc, the recombination time is \(t_{\rm rec} = 1.2\) Myr, and the Strömgren circle crossing time is \(t_{\rm cross} \equiv r_{\rm S}/c = 1650\) years.
Now use the cell below to plot the results and compare the ionization front (or I-front) to the analytic expressions for \(r_{\rm S}\) and \(r_{\rm I}\), shown in red and yellow quarter-circles respectively. The I-front should converge towards the red circle and should do so following the time-evolution of the yellow quarter-circle.
[ ]:
import osyris
import numpy as np
import matplotlib.pyplot as plt
import math
from ipywidgets import interactive
kb = 1.3807e-16 # Boltzmann constant in cgs
mp = 1.6749e-24 # Proton mass in grams
def show_stromgren_maps(snapshot=1):
fig, ax = plt.subplots(nrows=2, ncols=2, figsize=(9, 7), constrained_layout=False)
mydata = osyris.RamsesDataset(snapshot, path="./stromgren").load()
kpc = osyris.units("kpc")
center=osyris.Vector(0.5,y=0.5,unit=kpc)
# Map of gas density
osyris.map(mydata["mesh"].layer("density"), cmap="viridis", norm="log", vmin=1e-2, vmax=1e0
,dx=1*kpc,origin=center,ax=ax[0][0])
# Map of hydrogen ionization fraction
mydata["mesh"]["scalar_00"].name = r'${\rm x_{HII}}$'
osyris.map(mydata["mesh"].layer("scalar_00"), cmap="summer", norm="log", vmin=1e-4, vmax=1e0
,dx=1*kpc,origin=center,ax=ax[1][0])
# Map of temperature
mydata["mesh"]["T"] = mydata["mesh"]["pressure"]/mydata["mesh"]["density"] * mp/kb # Create new var for temperature
mydata["mesh"]["T"].name = r'${\rm T/\mu \ [K]}$'
mydata["mesh"]["T"].unit = ''
osyris.map(mydata["mesh"].layer("T"), cmap="coolwarm", norm="log", vmin=3e3, vmax=3e5
,dx=1*kpc,origin=center,ax=ax[0][1])
# Map of photon flux
mydata["mesh"]["photon_flux_01"] = mydata["mesh"]["photon_flux_01"]*mydata.meta["unit_l"]/mydata.meta["unit_t"]
mydata["mesh"]["photon_flux_01"].name = r'${\rm cN_{\gamma} \ [cm^{-1} \ s^{-1}]}$'
osyris.map(mydata["mesh"].layer("photon_flux_01"), cmap="bone", norm="log", vmin=1e5, vmax=1e8
,dx=1*kpc,origin=center,ax=ax[1][1])
time_myr = mydata.meta["time"].to("Myr").magnitude
plt.annotate('%.2f Myr'%(time_myr), xy=(0.1, 0.8),xycoords='figure fraction',color="grey")
show_stromgren_radius=True
if show_stromgren_radius:
Ndot = 2e28 # Photon injection rate in #/s
nH = 1e-1 # Initial gas density in #/cm3
rec_rate = 2.6e-13 # Case B rec. rate for T~10^4 K
myr2s = 60.*60.*24.*365.*1e6 # Myr to second conversion
r_S_cm = math.sqrt(Ndot/np.pi/nH**2/rec_rate) # in Centimeters
r_S = r_S_cm / 3.08e21 # in kpc
r_S = r_S * 1.05 # Pawlik&Schaye 2008 correction
t_rec_s = 1./nH/rec_rate # recomb. time in s
t_rec = t_rec_s/myr2s # --> Myr
c_cgs = 3e10 # Light speed (cgs)
crs1=plt.Circle((-0.5,-0.5),radius=r_S,color='r',alpha=0.5,fill=False,linewidth=3)
ax[1][0].add_patch(crs1) # Stromgren (final) radius
r_I = r_S*(1.-math.exp(-time_myr/t_rec))**0.5
cri1=plt.Circle((-0.5,-0.5),radius=r_I,color='y',alpha=0.5,fill=False,linewidth=3)
ax[1][0].add_patch(cri1) # Predicted current radius
interactive_plot = interactive(show_stromgren_maps, snapshot=(1, 46))
output = interactive_plot.children[-1]
output.layout.height = '700px'
interactive_plot
# In the lower left map below, the red circle shows the final Stromgren radius
# and the yellow circle shows the expected radius at the time corresponding to the snapshot
In the setup given to you, the light speed has been reduced to one-thousandth of the actual light speed (with the namelist parameter rt_c_fraction=0.001). This is good for performance – the lower the light speed the larger the timestep and the faster the simulation – but as you can see in the output, it is not very good for accuracy, as the actual front is lagging well behind the analytic result shown with the yellow circle.
In the RAMSES-RT code paper (Rosdahl et al., 2013), we derive a framework for setting a minimum value for the reduced light speed, while still retaining accuracy. The light speed fraction should be set to \(\min\left(1;\sim 10 \times t_{\rm cross}/\tau_{\rm sim}\right)\) where \(\tau_{\rm sim}\) is some simulation timescale within which we want to get an accurate I-front position. Use this expression to set a new (higher) value for the reduced light speed. Then verify that the I-front is accurately simulated within the timescale you choose (by comparing to the yellow expanding circle).
5.2 Turn on the photoionization heating of the gas
Now turn off the constant \(T = 10^4\) K (comment out rt_TConst=1d4 in the namelist), thus allowing the gas to cool or be heated by the photons. What changes in the evolution of the I-front, and why? (If you need a hint, look into appendix C, figure C.1 of the thesis of J Rosdahl.)
5.3 Turn on the hydrodynamic response of the gas
Now let’s make it more interesting: ‘unfreeze’ the gas (static=.false.) and allow it to react to the photoionization heating. Look at the I-front expansion in this new run. What happens, and why?
5.4 Play with RHD parameters
It is interesting to play with some parameters and consider the effects they have. You’re free to play with whatever parameters you like, but some ideas are:
Photon properties
Change the photoionization cross section and photon energy. The current values are averages for a \(T = 10^5\) K blackbody radiation. Increasing the energy and decreasing the cross sections corresponds to hotter stellar atmospheres, i.e. bigger stars. What does the cross section do to the I-front? What about the photon energy?
Adaptive mesh refinement
Turn on AMR refinement. Refine on e.g. mass, and gradients in density or ionization fractions. Can you get RAMSES-RT to focus the resolution on the I-front, but leave other regions coarsely resolved?
Spectral resolution
Try out three photon groups. This allows you to approximate some spectral resolution compared to the single photon group (Grey approximation) we’ve been working with. For the photon group properties (in the RAMSES-RT namelist), try a \(T = 10^5\) K blackbody spectrum:
rt_nsource=3 !--Idealized radiation sources------------
rt_src_group=1,2,3 ! Photon groups emitted into
rt_source_type=3*’point’ ! Source geometry (point or square)
rt_n_source=0.89d28, 0.99d28, 0.12d28 ! Total of 2d28 photons/s
&RT_GROUPS ! Blackbody at T=1d5 Kelvin
group_csn(1,:)= 3.0d-18, 0.,0. ! group 1-> HI, HeI, HeII
group_cse(1,:)= 2.8d-18, 0.,0. ! group 1-> HI, HeI, HeII
group_csn(2,:)= 5.7d-19, 0.,0. ! group 2-> HI, HeI, HeII
group_cse(2,:)= 5.0d-19, 0.,0. ! group 2-> HI, HeI, HeII
group_csn(3,:)= 7.9d-20, 0.,0. ! group 3-> HI, HeI, HeII
group_cse(3,:)= 7.5d-20, 0.,0. ! group 3-> HI, HeI, HeII
group_egy = 18.8, 35.1, 65.6 ! energy in units of eV
/
What happens to the I-front? You can also compare with a harder source of radiation, here a blackbody at \(T = 3 \times 10^5\) K:
rt_nsource=3 !--Idealized radiation sources------------
rt_src_group=1,2,3 ! Photon groups emitted into
rt_source_type=3*’point’ ! Source geometry (point or square)
rt_n_source=0.19d28, 0.63d28, 1.2d28 ! Total of 2d28 photons/s
&RT_GROUPS ! Blackbody at T=3d5 Kelvin
group_csn(1,:)= 2.8d-18, 0.,0.
group_cse(1,:)= 2.6d-18, 0.,0.
group_csn(2,:)= 4.1d-19, 0.,0.
group_cse(2,:)= 3.6d-19, 0.,0.
group_csn(3,:)= 3.9d-20, 0.,0.
group_cse(3,:)= 2.8d-20, 0.,0.
group_egy = 19.4, 39.6, 99.0
/
Helium
If you still have time (remember there are more problems to do), you can add helium to the mix. If you’ve already set up the three photon groups in the previous problem, all you need to do is specify in the namelist the hydrogen and helium mass fractions (X=0.76 and Y=0.24), and set the cross sections for helium ionization. For a \(T = 10^5\) K blackbody, separated into HI-ionizing, HeI-ionizing and HeII-ionizing photon groups, these are:
rt_nsource=3 !--Idealized radiation sources------------
rt_src_group=1,2,3 ! Photon groups emitted into
rt_source_type=3*’point’ ! Source geometry (point or square)
rt_n_source=0.89d28, 0.99d28, 0.12d28 ! Total of 2d28 photons/s
&RT_GROUPS ! Blackbody at T=1d5 Kelvin
group_csn(1,:)= 3.0d-18, 0., 0.
group_cse(1,:)= 2.8d-18, 0., 0.
group_csn(2,:)= 5.7d-19, 4.5d-18, 0.
group_cse(2,:)= 5.0d-19, 4.1d-18, 0.
group_csn(3,:)= 7.9d-20, 1.2d-18, 1.1d-18
group_cse(3,:)= 7.5d-20, 1.1d-18, 1.0d-18
group_egy = 18.8, 35.1, 65.6
/
What happens to the I-front?
Other
You could try inserting multiple radiation sources into the frame, or beams of radiation.
[ ]: