Setting up an accretion disk simulation in Astrobear 2.0

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 in 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 likely 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 the first run of this module.

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

The following image is a density plot of the disk at time t=0, a log scale has been applied to the 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 area around the disk
• PointGravity: a point mass located at the center of the disk

Initial Conditions

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

ddensity=1d0	disk density
adensity=0.001d0	ambient density
dtemp=1d4	disk temperature
atemp=6d6	ambient temperature
radius=2e10	disk radius
height=5e9	disk height
mass=.6d-1	central particle mass

Sometimes it can be useful to change the size of the physical boundaries of the simulation, this can be done by modifying 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

Scales

Scales are parameters that can actively affect the outcome of a 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!!!)

It is good practice to adapt the scales to the specific needs of a simulation, let's see an example:
In the code listed in the Initial Conditions section we set the disk radius to be 2e10cm, with that, we set our length scale parameter (lScale) to be .5e10cm
this will result in R_disk/lScale = 4.

Softening

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

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: the force of gravity goes to infinity as the radius decreases
• Softening enabled: the force of gravity has a well defined value at r=0
• Increased softening radius: 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.

Fixed grid to AMR

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

MaxLevel          = 2            ! Maximum 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 usually be set to (1, 1, 1, GmX(1), GmX(2), GmX(3)), changing Gmx without modifying Domain%mGlobal as well could lead to instabilities.

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

Below is an image of what the mesh looks like in the default module:

A subroutine named ProblemSetErrFlag in problem.f90 is responsible for marking the regions for refinement in our Disk module and giving the Mesh 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 a significant part of the disk to outflow the refined area, and a consequential loss of information.

ProblemSetErrFlag can be altered to address this issue, changing

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

the refined radius can be increased.

Mach Number

MHD and Solvers