Setting up an accretion disk simulation in Astrobear 2.0

This tutorial is designed to teach you the essential knowledge needed to set up a Disk module in AstroBEAR 2.0

Assuming that you obtained a copy of AstroBEAR 2.0, every file referenced in this tutorial is contained in the folder scambler/modules/BasicDisk, namely, the directory of our disk module.

Modules and Files

These directions assume that you are able to compile a problem, in case you have never compiled AstroBEAR check the page ​How to Build AstroBEAR 2.0

In order to set up a simple accretion disk simulation proceed to select the BasicDisk problem in the modules folder, this can be done by typing the following commands from the AstroBEAR root directory:

cd modules
ln -s BasicDisk Problem

After compilation, you need to be able to run the AstroBEAR executable, a tutorial for this procedure can be found at ​How to Run AstroBEAR 2.0 page.

Warning: It is good practice to delete profile.data from the problem folder before executing AstroBEAR for the first time, failure to do so will can result in execution errors.

Getting familiar with the module

After checking that the disk module compiles and runs without any issues, you can examine the chombo files created by this first run.

If you have no experience in visualization tools, the page ​VisIt Basics is a good resource to learn the fundamentals of it.

The following image is a density plot of the disk at time t=0, a log scale has been applied to these perpendicular and parallel slices.

At his point you should start to get familiar with some of the elements in the image above:

• GxBounds: the boundaries of the plot, namely the physical boundaries of the simulation
• Radius: the radius of the disk
• Height: the height of the disk
• Ambient: the blue portion of the plot, it represents the space around the disk
• PointGravity: a point mass located at the center of the disk, not visible in this picture

Note: Numerical values reported by Visit are in computational units, scales.data contains the scales values needed to translate computational units to physical units.

Initial Conditions

Now you can proceed to examine problem.data contained in the BasicDisk folder.
The default problem.data file contains a list of initial conditions that can be changed depending on the user's needs, let's analyze some of them

ddensity=1d0	disk density (g/cc or particles/cc)
adensity=0.001d0	ambient density (g/cc or particles/cc)
dtemp=1d4	disk temperature (K)
atemp=6d6	ambient temperature (K)
radius=2e10	disk radius (cm)
height=5e9	disk height (cm)
mass=.6d-1	central particle mass (g)

Note: for this specific module, all the values reported in problem.data are in physical units.

Sometimes it can be useful to change the size of the physical boundaries of the simulation, this can be done by changing the section below contained in global.data :

GxBounds          = -10.d0,         ! Problem boundaries in computational units.
                    -10.d0,         ! format: (xlower, ylower, zlower, xupper, yupper, zupper)
                    -10.d0,         ! For 2D problems, set zlower and zupper to 0.d0.
                     10.d0,
                     10.d0,
                     10.d0

Setting pressure

By default, the disk module lets you set the temperature for a problem, however, in the case you want to set pressure there are two options:

• Use of the equation of state to calculate temperature depending on pressure and density
• Manual modification of problem.f90 to initialize the problem using pressure parameters

For the second option, it may be sufficient to examine a few lines of ProblemModuleInit contained in problem.f90:

Ambient%pressure= adensity*atemp/rScale/TempScale !cu

the line above sets the ambient pressure depending of density and temperature values in problem.data, we can change it to

Ambient%pressure= atemp/pScale

after this change the value contained in atemp will be interpreted as the pressure in physical units.
Note: If are considering to make this change in the code, you are strongly advised to change variable names to avoid misleading nomenclature.

Scales

AstroBEAR makes use of scales to adapt physical values to a numerical simulation,
Physics.data contains a section dedicated to the setting of scales:

nScale          =       0.0,    ! number density scale parameter [particles/cc] 
rScale          =       1d0,    ! density scale [g/cc], nScale is ignored !!! NOTE: if = 0 then nScale is used to find rScale
TempScale       =       1000d0, ! temperature scale parameter [Kelvin]
pScale          =       0.0,    ! pressure scale [dynes/cm^2], TempScale is ignored !!! NOTE: if = 0 then TempScale is used to find pScale
lScale          =       .5e10   ! length scale parameter [cm] (defines 1 computational unit!!!)

Density scale
There are two options to set a density scale, you can choose either one:

• nScale in particles/cc
• rScale in g/cc

Length scale
In this tutorial we set the disk radius to 2*10¹⁰cm, with that, the length scale parameter (lScale) is set to .5*10¹⁰cm. The resulting ratio R_disk/lScale is 4.
For convenience purposes, a rule of thumb is to keep R_disk/lScale within one order of magnitude from the value above:
for example, an ideal length scale for R_disk = 2*10¹⁵cm would be .5*10¹⁵cm or 1*10¹⁵cm.

Warning: Always check physics.data for

MinTemp         =       1e-3,      ! minimum allowed temperature for the system, without cooling, this should be zero.
MinDensity      =       1e-3       ! Minimum computational density

and set them to appropriate values depending on the chosen scales.

Softening

Softening is a numerical trick commonly used in simulations to prevent divergences when a particle comes close to another and the force of gravity goes to infinity.
In our case, softening approximates the force of gravity to a well defined value as the disk radius gets smaller.

Problem.data provides a field in which the user can set the preferred softening parameter:

soft_radius=.333d0 !diks radii
soft_function=2  !NOSOFT=0, SPLINESOFT=1, PLUMMERSOFT=2

The plot below shows 3 different scenarios:

• No Softening (soft_function=0) : the force of gravity goes to infinity as the radius decreases
• Softening enabled (soft_function=2, soft_radius=.1d0) : the force of gravity has a well defined value at r=0
• Increased softening radius, (soft_function=2, soft_radius=.1d3): softening starts at a greater radial distance, the force of gravity has a lower value at r=0 compared to the previous case

Astrobear provides different softening techniques such as SplineSoft and PlummerSoft, the user can set the softening function depending on the simulation needs.

Warning: It is good practice to keep the soft_radius between .1 and .3. Failure to do so may result in undesired outcomes.

AMR

The resolution of the simulation can be changed by modifying the following fields contained in ​global.data:

MaxLevel          = 2            ! Maximum refinement level for this simulation.
GmX               = 32,          ! number of cells in x direction
                    32,          ! number of cells in y direction
                    32,          ! number of cells in z direction (set to 1 for 2D problems).

MaxLevel regulates the maximum amount of AMR levels, while GmX sets the number of cells in a fixed grid array.

Warning: global.data contains another array entry, Domain%mGlobal, which should be set to (1, 1, 1, GmX(1), GmX(2), GmX(3)), for the GmX values seen above, Domain%mGlobal should look like:

Domain%mGlobal=1,1,1,32,32,32

In this specific module refining is accomplished by defining a region around the disk that will always be marked for refinement. For more information on how this is done please visit ​Controlling Refinement in AstroBEAR 2.0.

Below is an image of our mesh in the default module:

A subroutine named ProblemSetErrFlag in problem.f90 is responsible for marking the regions for refinement and giving the pattern seen above.

Depending on the initial conditions of a specific run it may be the case that a disk undergoes gravitational collapse and its boundaries can expand within reasonable limits. This event can cause part of the disk to outflow the refined area.

ProblemSetErrFlag can be altered to address this issue. By changing the following line

if (sqrt((x-xloc(1))**2+(y-xloc(2))**2).le.radius/lScale+2d0*dx) Info%ErrFlag(i,j,k)=1

to

if (sqrt((x-xloc(1))**2+(y-xloc(2))**2).le.radius*1.2d0/lScale+2d0*dx) Info%ErrFlag(i,j,k)=1

we can increase the radius in which additional refinement in enabled.

Here is a comparison of our mesh before and after this change.

Notes:

• You need to recompile AstroBEAR after making any changes to problem.f90
• Keep in mind that increasing the refined region also increases computational time

Mach Number

When setting up a new simulation it is good practice to check that the Mach values are below a certain threshold, in most of numerical simulations this value is ~60.
Values above this threshold may create difficulties in evolving a model and therefore create computational errors.

To plot the Mach numbers in Visit you can check ​Chombo expressions and locate Mach in the expression list.

The plot below represents a disk simulation that respects the limits stated above:

Additional parameters

Now that you have acquired the essential knowledge to set up your disk module, you should be capable to make other changes such as:

• Equation of state
• Solvers
• MHD
• Boundary conditions
• Self gravity
• Accretion

You are advised to analyze every .data file to familiarize with other aspects of this module that were not discussed in this tutorial.