Skip to content

Working with Bifrost snapshots

Using IDL

Setting up IDL with Bifrost extension

There is lots of useful tools written in IDL for manipulating Bifrost data, thus it is worth it to setup your enviroment with IDL extensions which will help you read Bifrost data.

You should start by compiling cstagger module (for fast stagger operations), which you will find in $BIFROST/IDL/cstagger. There are folders specific for:

- `armmac` - Apple Macs with ARM cpu (M1,M2 or M3)
- `intelmac` - Apple Macs with Intel CPUs
- `linux` - Any LINUX machine

Choose one for your architeqture and simply execute make command in it.

Warning

This will work under assumplion that you have installed C compilers on your machine (for Mac OS this is done via installing xcode and on Linux, there should be a package for GNU Gcc)

Next you have to set few enviromental constants namely:

export BIFROST=/path/to/Bifrost/
export BIFROST_CSTAGGER="${BIFROST}/IDL/cstagger/linux"

Warning

Remember to select correct type of machine in definition of BIFROST_CSTAGGER. Example provide solution for Linux machines.

Best to put them in your default rc file for your shell (e.g., $HOME/.bash_profile or $HOME/.zshrc )

Now we have to expand IDL_PATH and IDL_DLM_PATH, so IDL can find IDL routines from BIFROST, you can do it from your rc file but better to do it from custom IDL startup file.

First define IDL_STARTUP in your shell profile file:

export IDL_STARTUP=$HOME/.idlstartup

Tip

If you only going to use IDL for bifrost you can use template for .idlstatup from Bifrost source folder: 

cd $HOME
ln -s $BIFROST/IDL/setup.bifrost.idl .idlstartup

And modify $HOME/.idlstartup after with following lines:

print, ''
print, ':: ... Call from .idlstartup'
print, ':: ...'
print, ':: ... setting up BIFROST/IDL & cstagger'

bifrost_dir = getenv('BIFROST')
bifrost_cstagger_dir = getenv('BIFROST_CSTAGGER')

print, ':: ... BIFROST env: ', getenv('BIFROST')
print, ':: ... BIFROST CSTAGGER env: ', getenv('BIFROST_CSTAGGER')

bifrost_idl = '+' + bifrost_dir + '/IDL/:<IDL_DEFAULT>'
bifrost_ddl =       bifrost_stagger_dir + ':<IDL_DLM_PATH>'

pref_set, 'IDL_PATH', bifrost_idl, /commit
pref_set, 'IDL_DLM_PATH', bifrost_ddl, /commit

print, ''

Warning

Note that in bifrost_ddl you have to select cstagger folder suitable for your machine (e.g., for linux)

Testing you setup

If everything is set correctly you can test it by entering any directory with Bifrost simulation and typing:

:: ... Call from .idlstartup
:: ...
:: ... BIFROST/IDL & cstagger/linux
:: ... BIFROST env: /mn/stornext/d9/data/mikolajs/Bifrost

IDL> .r cstagger
% Compiled module: XUP.
% Compiled module: YUP.
% Compiled module: ZUP.
% Compiled module: XDN.
% Compiled module: YDN.
% Compiled module: ZDN.
% Compiled module: DDXUP.
% Compiled module: DDYUP.
% Compiled module: DDZUP.
% Compiled module: DDXDN.
% Compiled module: DDYDN.
% Compiled module: DDZDN.
IDL> d = br_obj("ch012012_hion",isnap=845)
% Compiled module: BR_OBJ.
% Compiled module: BR_DATA__DEFINE.
% Compiled module: BR_HDR__DEFINE.
% Compiled module: BR_HION__DEFINE.
% Compiled module: BR_HELIUM__DEFINE.
% Compiled module: BR_AUX__DEFINE.
% BR_AUX::READTABPARAM: reading tabparam file tabparam.in
% BR_HDR::INITHDR: BR Header$Revision: 1.16 $ params file root: ch012012_hion
% BR_HDR::INITHDR: BR Data$Revision$ params file root: ch012012_hion
% Compiled module: BR_STRING3.
% BR_HDR::READDATAMESH: reading ch012012_hion.mesh
% Compiled module: INIT_STAGGER.
IDL> z = d->getz()
IDL> help,z
Z               FLOAT     = Array[1024]
IDL> z.min()
      -7.9821329
IDL> z.max()
       2.4980001
IDL> tg = d->getvar('tg',845)
% BR_HDR::READDATAMESH: reading ch012012_hion.mesh
% BR_DATA::FINDVAR: reading tg from ch012012_hion_845.aux isnap 845
tg.mean()IDL> tg.mean()
       31943.252

In above example we have laoded bifrost object, then z-axis of the mesh and temperature variable to calculating mean value.

Working with Bifrost data - Python

Python modules to import and parameter

d = bifrost.BifrostData()

r = np.memmap("file_name") # (1)
  1. 🙋‍♂️ this is footnote.
# importing modules for later analysis
import os
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

from matplotlib.gridspec import GridSpec

plt.rcParams['figure.figsize'] = [8, 5]

Simplified I/O routines for Bifrost raw format

# Simplified I/O routines for Bifrost raw format
def Bifrost_read_params(filename):
    ''' Reads Bifrost params file into dictionary '''

    params = {}

    # go through the file, add stuff to dictionary
    with open(filename) as f:
        for line in f:
            line = line.strip()
            if len(line) <= 1:
                continue
            if line[0] == ';':
                continue

            line = line.split(';')[0].split('=')

            key   = line[0].strip().lower() # force lowercase because IDL is case-insensitive
            value = line[1].strip()

            if (value.find('"') >= 0):
                value = value.strip('"')
            elif (value.find('[') >= 0 and value.find(']') >= 0):
                value = eval(value)
            elif (value.find("'") >= 0):
                value = value.strip("'")
            elif (value.lower() in ['.false.', '.true.']):
                value = False if value.lower() == '.false.' else True
            elif ((value.upper().find('E') >= 0) or (value.find('.') >= 0 )):
                value = float(value)
            else:
                try:
                    value = int(value)
                except:
                    print('Value exception')
                    continue
            params[key] = value

        return params

def Bifrost_load_mesh(meshfile):
    """
    Reads bifrost mesh file but return only x,y & z axis ( units [ Mm ] )
    """
    if os.path.isfile(meshfile):
        f = open(meshfile, 'r')
        for p in ['x', 'y', 'z']:
            dim = int(f.readline().strip('\n').strip())
            # quantity
            locals()[p] = np.array([float(v) for v in f.readline().strip('\n').split()])
            # quantity "down"
            locals()[p + 'dn'] = np.array([float(v) for v in f.readline().strip('\n').split()])
            # up derivative of quantity
            locals()['d%sid%sup' % (p, p)] = np.array([float(v) for v in f.readline().strip('\n').split()])
            # down derivative of quantity
            locals()['d%sid%sdn' % (p, p)] = np.array([float(v) for v in f.readline().strip('\n').split()])
        f.close()
        return eval("x,y,z")

def Bifrost_load_snapdata(snapshot):
    """
    Reads Bifrost *.snap binarry file as numpy array in dimension: [mx,my,mz,nvar]
    if "do_mhd = 1", nvar = 8, otherwise nvar = 5
        vars:
            datasnap[:,:,:,0] : r, density
            datasnap[:,:,:,1] : px, x-component of momentum
            datasnap[:,:,:,2] : py, y-component of momentum
            datasnap[:,:,:,3] : pz, z-component of momentum
            datasnap[:,:,:,4] : e, energy

            if "do_mhd == 1" :
                datasnap[:,:,:,5] : bx, x-component of magnetic field
                datasnap[:,:,:,6] : by, y-component of magnetic field
                datasnap[:,:,:,7] : bz, z-component of magnetic field
    Warning:
        variables in code units !!
    """
    snapstr = "_%03d" % snapshot
    base_name = expfolder+"/"+expname+snapstr
    idl_filename = base_name+".idl"
    snap_filename = base_name+".snap"

    params = Bifrost_read_params(idl_filename)

    if params["do_mhd"] == 1:
        nvars = 8
    else:
        nvars = 5

    snap_size = (params["mx"],params["my"],params["mz"],nvars)
    data_snap = np.memmap(snap_filename,dtype=np.float32, mode="r",order="f",shape=snap_size);
    return data_snap

def Bifrost_load_auxdata(snapshot):
    """
    Reads Bifrost *.aux binarry file as numpy array in dimension: [mx,my,mz,nvar]
    nvar is determined by lenght of `aux = ` list in param file.
    """
    snapstr = "_%03d" % snapshot
    base_name = expfolder+"/"+expname+snapstr
    idl_filename = base_name+".idl"
    snap_filename = base_name+".aux"

    params = Bifrost_read_params(idl_filename)
    aux_nvars = len(params["aux"].split())

    if aux_nvars == 0:
        return None
    else:
        aux_size = (params["mx"], params["my"], params["mz"], aux_nvars)
        aux_snap = np.memmap(snap_filename, dtype=np.float32,
                          mode="r", order="f", shape=aux_size)
        return aux_snap

def Bifrost_load_header(snapshot):
    """
    Reads Bifrost param file as list of lines
    """
    widthl=3
    fill='0'
    snapstr=f'{snapshot:{fill}{widthl}}'
    filename=expfolder+"/"+expname+snapstr+".idl"
    fileread = open(filename, 'r')
    lines=fileread.readlines()
    fileread.close()
    return lines

Grid Stagger

Remember that the variables are stored on a staggered grid and not all at cell centre values. It is important to remember this and account for it properly when working with Bifrost data.

Tip

Lorem ipsum dolor sit amet, consectetur adipiscing elit. Nulla et euismod nulla. Curabitur feugiat, tortor non consequat finibus, justo purus auctor massa, nec semper lorem quam in massa.

# Simplified Stagger operations (to shift values from face to cell center)
def dnup(q, shift=1, axis=0):
   i = 0 if shift==1 else -1
   if axis == 0 and q.shape[0] > 1:
      f = (q+np.roll(q,shift,axis))*0.5
      f[i,:,:] = q[i,:,:]
   elif axis == 1 and q.shape[1] > 1:
      f = (q+np.roll(q,shift,axis))*0.5
      f[:,i,:] = q[:,i,:]
   elif axis == 2 and q.shape[2] > 1:
      f = (q+np.roll(q,shift,axis))*0.5
      f[:,:,i] = q[:,:,i]
   else:
      f = q
   return f

# The actual routines you call to shift the values up or down in x, y, or z.
def xdn(f):
   return dnup(f,1,0)
def ydn(f):
   return dnup(f,1,1)
def zdn(f):
   return dnup(f,1,2)

def xup(f):
   return dnup(f,-1,0)
def yup(f):
   return dnup(f,-1,1)
def zup(f):
   return dnup(f,-1,2)

Let's look at an experiment's parameters using these

Specify which folder (i.e. which experiment) we want to look at.

# setting experiment parameters
expfolder="ExperimentA"
expname="ch16r_2d"

Read param file, to see the size information and more..

# read the initial experiment snapshot .idl" file.
# this is a text file containing parameters that will help up load in the data
params = Bifrost_read_params("ExperimentA/ch16r_2d_001.idl")

Value exception

## params now is a dictionary of all parameters for pointed *.idl file. At the very top you can find m[x,y,z] with size of the mesh
params

## getting a specific value
#tempvar=params["u_u"]
#print(tempvar)

{'mx': 512,
 'my': 1,
 'mz': 496,
 'mb': 5,
 'nstep': 200000,
 'nstepstart': 1,
 'debug': 0,
 'time_lim': -1.0,
 'periodic_x': 1,
 'periodic_y': 1,
 'periodic_z': 0,
 'ndim': 3,
 'reorder': 1,
 'u_l': 100000000.0,
 'u_t': 100.0,
 'u_r': 1e-07,
 'u_p': 100000.0,
 'u_u': 1000000.0,
 'u_kr': 0.1,
 'u_ee': 1000000000000.0,
 'u_e': 100000.0,
 'u_te': 100000000000.0,
 'u_tg': 9692.0,
 'u_b': 1121.0,
 'meshfile': 'ch16r_2d.mesh',
 'dx': 0.03134,
 'dy': 1.0,
 'dz': 0.06683,
 'cdt': 0.3,
 'dt': 0.0001131,
 't': 0.0,
 'timestepdebug': 0,
 'nu1': 0.05,
 'nu2': 0.1,
 'nu3': 0.15,
 'nu_r': 0.02,
 'nu_r_mz': 0.1,
 'nu_r_k': -1,
 'nu_r_z': -0.5,
 'nu_r_w': 0.072,
 'nu_ee': 0.05,
 'nu_ee_mz': 0.1,
 'nu_ee_k': 0,
 'nu_ee_z': -0.5,
 'nu_ee_w': 0.072,
 'symmetric_e': 1,
 'symmetric_b': 0,
 'grav': -2.74,
 'eta3': 0.3,
 'ca_max': 0.0,
 'mhddebug': 0,
 'do_mhd': 1,
 'mhdclean': -1,
 'mhdclean_ub': 0,
 'mhdclean_lb': 0,
 'snapname': 'ch16r_2d_nobdrchk',
 'isnap': 1,
 'large_memory': 1,
 'nsnap': 100000000,
 'nscr': 100,
 'aux': ' qpdv qjoule qvisc qtot qgenrad qspitz florentz1 florentz2 florentz3 fpress1 fpress2 fpress3 fstress1 fstress2 fstress3 frdiff1 frdiff2 frdiff3 fpdiff1 fpdiff2 fpdiff3 qediff   ',
 'dtsnap': 0.01,
 'newaux': 0,
 'rereadpars': 100000000,
 'dtscr': 1000.0,
 'tsnap': 0.0,
 'tscr': 1.0,
 'boundarychk': 0,
 'max_r': 5,
 'smooth_r': 3,
 'divhc_niter': 1000,
 'divhc_cfl': 0.4,
 'divhc_r': 0.18,
 'divhc_vxr': 0.0,
 'divhc_vyr': 0.0,
 'divhc_vzr': 0.95,
 'divhc_tol': 1e-05,
 'qmax': 8.0,
 'noneq': 0,
 'expieos': 0,
 'do_hion': 0,
 'do_helium': 0,
 'gamma': 1.667,
 'debug_square': 0,
 'maxlnelim': 32.3,
 'minlnrlim': -34.0,
 'tabinputfile': ' tabparam.in  ',
 'do_rad': 1,
 'dtrad': 0.003,
 'quadrature': 3,
 'zrefine': 3,
 'maxiter': 300,
 'taustream': 1e-05,
 'accuracy': 0.001,
 'strictint': 1,
 'linear': 1,
 'monotonic': 1,
 'minbin': 1,
 'maxbin': 4,
 'dualsweep': 0,
 'teff': 5832.0,
 'timing': 0,
 'spitzer': 'impl',
 'debug_spitzer': 0,
 'info_spitzer': 1,
 'spitzer_amp': 1.0,
 'theta_mg': 0.95,
 'dtgerr': 1e-06,
 'ntest_mg': 10,
 'tgb0': 500000.0,
 'tgb1': 10000.0,
 'tau_tg': 0.005,
 'fix_grad_tg': 0,
 'niter_mg': [2, 20, 10, 10, 30],
 'bmin': 1e-05,
 'do_genrad': 1,
 'genradfile': 'radtab.dat',
 'debug_genrad': 0,
 'incrad_detail': 0,
 'incrad_quad': 3,
 'dtincrad': 0.001,
 'dtincrad_lya': 0.0001,
 'debug_incrad': 0,
 'ck_nrcorks': 253952,
 'ck_xstart': 1,
 'ck_xend': 512,
 'ck_xstep': 1,
 'ck_ystart': 1,
 'ck_yend': 1,
 'ck_ystep': 1,
 'ck_zstart': 1,
 'ck_zend': 496,
 'ck_zstep': 1,
 'ck_placement': 'grid           ',
 'ck_relocate': 0,
 'ck_inject': 1,
 'ck_max_density': 2,
 'ck_max_voidsize': 1.0,
 'ck_inject_delay': -1,
 'ck_highest_id': 253952,
 'bctypelower': 'mfgc',
 'bcldebug': 0,
 'tau_bcl': 1.0,
 'tau_ee_bcl': 100.0,
 'tau_d2_bcl': 100.0,
 'tau_d5_bcl': 100000.0,
 'tau_d6_bcl': 100000.0,
 'tau_d7_bcl': 100000.0,
 'tau_d8_bcl': 100000.0,
 's0': 15.3791008,
 'e0': 0.0700427964,
 'r0': 0.0186682008,
 'cs0': 0.654462993,
 'p0': 0.00481020007,
 'nsmooth_bcl': 50,
 'nextrap_bcl': 1,
 'naver_bcl': 0,
 'nsmoothbxy_bcl': 0,
 'nsmoothbz_bcl': 0,
 'nclean_bcl': 50,
 'nclean_lbl': -4,
 'nclean_ubl': 5,
 'bctypeupper': 'dcs',
 't_bdry': 1.0,
 'rbot': 264.23999,
 'ebot': 2900.0,
 'nclean_bcu': -1,
 'nsmoothbxy_bcu': 0,
 'pb_spot': 0.0,
 'spot_upper': 0.0,
 'beta_lim': 100.0,
 'bx0': 0.0,
 'by0': 0.0,
 'bz0': 0.0,
 'x0_bcu': 0.0,
 'x1_bcu': 16.0,
 'y0_bcu': 0.0,
 'y1_bcu': 16.0,
 'uz_bcu': 0.00200000009,
 'strtb': 'n'}

print(f"Experiment '{expfolder}' has following size: mx = {params['mx']}, my = {params['my']}, mz = {params['mz']}")

if (params["do_mhd"] == 1):
   print("This is MHD experiment with 8 variables in snap file.")
else:
   print("This is HD (Hydro) experiment with 5 variables in snap file (no magnetic field)")

Experiment 'ExperimentA' has following size: mx = 512, my = 1, mz = 496
This is MHD experiment with 8 variables in snap file.

Bifrost Snapshot file

Bifrost *.snap files are FORTRAN unformatted binarry files with 8 or 5 records of size [mx, my, mz]. These are easy to ready with numpy (see implementation of Bifrost_load_snapdata() )

First lets read snap 002

snap = Bifrost_load_snapdata(2) # read snap files {expname}_001.snap from {expfolder}
print(snap.shape)
test=np.squeeze(snap)
print(test.shape)

(512, 1, 496, 8)
(512, 496, 8)

Note that snap is 4 dimensional array, where 4th dimension is for different mhd variables. Note also how numpy's squeeze routine will remove redundant dimensions

Let's save them under separate names

rho = snap[:,:,:,0] # first is density
px  = snap[:,:,:,1] # momentum, x component i.e., ux * rho
py  = snap[:,:,:,2] # momentum, y component i.e., uy * rho
pz  = snap[:,:,:,3] # momentum, z component i.e., uz * rho
e   = snap[:,:,:,4] # energy

if (params["do_mhd"] == 1): # then we have also magnetic field components
   bx = snap[:,:,:,5] # x component of mag field
   by = snap[:,:,:,6] # y component of mag field
   bz = snap[:,:,:,7] # z component of mag field

All these variables are in code units and need to be multiply (scaled) by some constant to represent meanigful data. We store all these constant in *.idl (param) file, thus we can for instance quickly define density in CGS units:

rho_cgs = rho * params["u_r"] # rho now in g/cm^3
px_cgs  = px * params["u_r"] * params["u_u"] # g/cm^2/s
py_cgs  = py * params["u_r"] * params["u_u"] # g/cm^2/s
pz_cgs  = pz * params["u_r"] * params["u_u"] # g/cm^2/s
e_cgs = e * params["u_e"] # erg/cm^3
if (params["do_mhd"] == 1):
    bx_cgs = bx * params["u_b"] # G
    by_cgs = by * params["u_b"] # G
    bz_cgs = bz * params["u_b"] # G

Danger zone!

Not all variables are co-alined on the experiment mesh as explained in the introduction, thus we can not simply define velocity from momentum with ux = px / rho. Staggering need to be taken into account. All scallars are cell-centered, but vector quatities are define on faces or edges, thus proper interpolation need to be applied. For the puprpose of this tutorial we are going to use simplified formulation of stagger operations.

Calculating Velocity

ux_cgs = xup(px) / rho * params["u_u"] # xup moves values half grid in x-direction, params["u_u"] scales velocity to CGS
uy_cgs = yup(py) / rho * params["u_u"]
uz_cgs = zup(pz) / rho * params["u_u"]

Auxiliary variables

Additionally to *.snap binary files you might have optional *.aux files which stores derived quanties such as: temperature, pressure, varius heating terms and many many more. To find out what has been saved in given *.aux file we have to explore params variable aux which provide full list of extra arrays.

print(f"Experiment '{expfolder}' has n = {len(params['aux'].split())} auxiary variables. Here is a full list:")

aux_list = params['aux'].split()

for i, var in enumerate(aux_list):
   print(f"{i:02} : {var}")

Experiment 'ExperimentA' has n = 22 auxiary variables. Here is a full list:
00 : qpdv
01 : qjoule
02 : qvisc
03 : qtot
04 : qgenrad
05 : qspitz
06 : florentz1
07 : florentz2
08 : florentz3
09 : fpress1
10 : fpress2
11 : fpress3
12 : fstress1
13 : fstress2
14 : fstress3
15 : frdiff1
16 : frdiff2
17 : frdiff3
18 : fpdiff1
19 : fpdiff2
20 : fpdiff3
21 : qediff

Aux variable names

The first letter generally indicates the type of term you are dealing with

q: heating
f: force

To get an acceleration term one can convert by dividing the force by density (remembering to adjust for the staggered values of course).

Directions

1=x, 2=y, 3=z.

see here for more details

Lets read aux file in similar fashion as snap file:

aux = Bifrost_load_auxdata(2)
aux.shape # note that last dimension should be == to number of aux variables in list above
(512, 1, 496, 21)

Lets extract qjoule variable (Joule heating) and convert units to CGS i.e., [erg/cm^3/s]

qjoule = aux[:,:,:,1] * params["u_e"]/params["u_t"] # 01 on the list, `u_e` scaling factor # since this is an energy term, `u_t` for time scaling

Displaying a Bifrost atmosphere on the correct grid

Z-axis in pixel space is not very useful. We have to read mesh data to find physical units for z-axis.

x,y,z = Bifrost_load_mesh("common/ch16r_2d.mesh") # axis values in Mm units

2D Plots

plt.imshow(np.log10(np.squeeze(rho_cgs.T)),cmap='jet') # .T is for transpose because of the python array ordering
plt.colorbar()
plt.title(r"$\log_{10}(rho)$ [g.cm^3]")
plt.xlabel("x-axis [pixels no]")
plt.ylabel("z-axis [pixels no]")
png

2D plot taking into account mesh cell sizes (non-uniform z axis)

from matplotlib.image import NonUniformImage


fig, ax = plt.subplots()
fig.suptitle('Rho as NonUniformImage')
im = NonUniformImage(ax, interpolation=None, extent=(x.min(), x.max(), z.min(), z.max()),cmap='jet')
im.set_data(x, z, np.log10(np.squeeze(rho_cgs.T)))     # squeeze removes the redundant y-axis from the data
ax.images.append(im)
ax.set_xlim(x.min(), x.max())
ax.set_ylim(z.max(), z.min())
ax.set_xlabel("x-axis [Mm]")
ax.set_ylabel("z-axis [Mm]")
ax.axes.set_aspect('equal')
png

1D Plots (average profiles) 

def plot_1d(var,title,log=True,zaxis=None):
   if zaxis is None:
      if log:
         plt.semilogy(var)
      else:
         plt.plot(var)
      plt.xlabel('z-axis [pixels]')
   else:
      if log:
         plt.semilogy(zaxis,var)
      else:
         plt.plot(zaxis,var)
      plt.xlabel('z-axis [Mm]')
   plt.title(title)

Lets calculate and plot average density profile from rho_cgs

rho_avrg = np.average(rho_cgs,axis=(0,1))     # take average over x (and y) values
#print(np.shape(rho_cgs))                     # uncomment to see why averaged over axis 0 and axis 1

plt.semilogy(rho_avrg)
plt.xlabel('z-axis')
plt.ylabel('rho (density) [g/cm^3]')
plt.title("Average density profile")
png

now insert z grid axis values

plt.semilogy(z,rho_avrg)
plt.xlabel('z-axis [Mm]')
plt.ylabel('rho (density) [g/cm^3]')
plt.title("Average density profile")
png

Inspecting non-uniformity of z axis

# Note that Bifrost mesh is non-uniform in z-direction
# (smaller cells in transition region and chromosphere)
dz = np.ediff1d(z)
plt.plot(z[1:],dz)
plt.xlabel("z-axis [Mm]")
plt.ylabel("dz [Mm]")

# photosphere ~ 0.0 Mm
# positive values below Photosphere
# negative values above
png

Extension: 2D plotting with interpolation to regular grid

from scipy.interpolate import interp1d

new_z = np.linspace(z.min(),z.max(),512) # define new regular z-axis
fint = interp1d(z,  rho_cgs, axis=-1, kind='cubic') # define interpolation object for rho_cgs
rho_cgs_on_new_z = fint(new_z)
rho_cgs_on_new_z.shape # now rho is on regular grid "new_z"

plt.imshow(np.log10(np.squeeze(rho_cgs_on_new_z.T)), cmap='jet', extent=(x.min(),x.max(),z.max(),z.min()))
plt.colorbar()
plt.title(r"$\log_{10}(rho)$ [g.cm^3]")
plt.xlabel("x-axis [Mm]")
plt.ylabel("z-axis [Mm]")

png

Streamlines

Now let's get some modulus B-field and Velocity streamlines.

# get bx,bz magnetic filed components

# extract data

snap = Bifrost_load_snapdata(10)
params = Bifrost_read_params("ExperimentA/ch16r_2d_010.idl")
bx = np.squeeze(snap[:,:,:,5]) # x component of mag field
by = np.squeeze(snap[:,:,:,6]) # y component of mag field
bz = np.squeeze(snap[:,:,:,7]) # z component of mag field
rho = np.log10(np.squeeze(snap[:,:,:,0])).T

bs = np.sqrt(bx*bx + by*by + bz*bz)

# plot
imgplot = plt.imshow(rho, extent=[0,params['mx'],params['mz'],0], cmap='jet')

# title and labels
plt.title("Rho with magnetic field stream lines")
plt.ylabel(r"z (Cell)")
plt.xlabel(r"x (Cell)")
plt.colorbar()

# add stream plots
x = np.arange(0, params['mx'])
y = np.arange(0, params['mz'])
Xgrid, Ygrid = np.meshgrid(x, y)
plt.streamplot(Xgrid, Ygrid, bx.T, bz.T, color="White", integration_direction="both")

png

# Control stream line width with field strength

imgplot = plt.imshow(rho, extent=[0,params['mx'],params['mz'],0], cmap='jet')

# title and labels
plt.title("Rho with magnetic field stream lines")
plt.ylabel(r"z (Cell)")
plt.xlabel(r"x (Cell)")
plt.colorbar()

# add stream plots
x = np.arange(0, params['mx'])
y = np.arange(0, params['mz'])
Xgrid, Ygrid = np.meshgrid(x, y)
plt.streamplot(Xgrid, Ygrid, bx.T, bz.T, color="White", integration_direction="both",linewidth=bs.T*3000.0,density=1.8)
png

How-To

How to use newaux, i.e., how to expand aux variable list in the existing simulation.

It is possible to expand aux variable array for the range of snapshots without re-running the entire simulation. To do so: Edit *.idl param file of the first snapshot in the range of interest. In this file, edit aux param in io section, providing a new set of variables. Change the value of newaux variable from 0 to the last snapshot number in your range. If you intend to expand aux only for one snapshot, set this number to isnap value from *.idl file you are currently editing.

To expand aux variable, Bifrost will have to execute one iteration that can be costly for big snapshots, making sure you will execute this procedure on a big enough machine.

Important:

New aux variables will be evaluated for time dt away from basic MHD variables stored in *.snap files i.e., after this procedure, aux variables will not be co-temporal with basic MHD variables, and the time difference will be exactly dt. This temporal difference is recorded in each param file as the value of dtaux. If dtaux == 0, MHD and aux variables are cotemporal.

Note that dt values for individual snapshots can be different, but you can optionally set dtaux in first snap in the range to some constant (small) value which will be use as new dt value in production of new aux list.

How to get conversion factors for magnetic field, current density and electric field?

Here we aim at finding conversion factors u_b, u_i and u_el, from Bifrost units to CGS, for the magnetic field, current density and electric field, respectively (see also here).

First note that Maxwell's equations differ in Bifrost units, such that 

dB'dt' = - rot_stagger(E') ; Faraday's
rot_stagger(B) = J' = I ; Ampere's
and that the Ohm's law is then express as E' = eta - v' x B', where ' values represent quantities in Bifrost units and eta is the non-ideal part of the electric field in such units, where the hyper-diffusive resistivity is coded in Bifrost.

Let's now recap that in CGS we have instead 

dBdt = - c.rot(E) ; Faraday's
rot(B) = 4pi/c * J ; Ampere's
E = 4pi*diff*J/c^2 - (v x B)/c ; Ohm's

Getting u_b

Bifrost units give by definition that I = J' = rot_stagger(B') = J *u_b /c_cgs = rot(B * u_b) ie. J *sqrt(4pi)/c = sqrt(4pi) * rot(B) , which is fine in CGS and also means that J is rescaled by sqrt(4pi)/c in the Bifrost system (see below).

Getting u_i

Let's say we want to reconstruct rot(B) in CGS, we will then have to do rot_stagger(B') *u_b/u_l, (nothing more as B=B'*u_b by definition).

Using Ampere's law, we want J=rot(B)*c/4pi in CGS, then we will have to do the same ie. J = rot_stagger(B') *u_b/u_l *c_cgs/4pi

however, we also have J'= I = rot_stagger(B') in the Bifrost system of unit, ie. I and rot_stagger(B') are identical in Bifrost output

so to get back to J in cgs from J'=I=rot_stagger(B'), we take back the previous formula J = rot_stagger(B') *u_b/u_l *c_cgs/4pi = u_i * I , which gives us finally u_i = u_b/u_l *c_cgs/4pi by definition

Getting u_e

from the Faraday's law 

dB'dt' = - rot_stagger(E') ; Bifrost
dBdt = -c_cgs *rot(E)` ; CGS

Knowing dBdt = dB'dt' * u_b /u_t this yields -c_cgs * rot(E) = -rot_stagger(E') * u_b/u_t then rot() = rot_stagger() / u_l so -c_cgs * rot(E) = -rot(E') * u_b*u_l/u_t

and finally by definition rot(E) = rot(E')*u_el So u_el = u_b*u_l/u_t/c_cgs by definition.

How to get conversion factors for magnetic diffusivities?

eta_hall and eta_amb are equivalent to eta_ohm in terms of unit (Martinez-Sykora et al. 2012). MHD equations give that fundamentally (if you omit Hall and Ambipol. effect. for the moment) eta_ohm = (the diffusive part of the electric field = E_ohm) / (rotational of B) so this leads you to (length)^2/(time) units. (Let's note that E_ohm = etax,y and z in Bifrost output)

If you are in Bifrost units, then this is etax_ohm = E_ohm / J = 1/sigma (= 'etax'/'Ix' from Bifrost output). BUT If you want to get (length)^2/(time) in CGS, then you have from MHD eqs: E + (u x B)/c = J/sigma = 4pi*eta_ohm*J/c^2 hence eta_ohm = ( E + (v x B)/c ) * c^2/4pi/J then from Bifrost you need to take (etax/Ix) * u_el/u_i * c^2/4pi then the conversion term to CGS for a diffusivity is u_el/u_i * c^2/4pi = u_l^2/u_t

Same exercice in SI: you have then E + u x B = mu0*eta_ohm*J hence eta_ohm = ( E + v x B ) /J /mu0 then from Bifrost you need to take (etax/Ix) * u_el/u_i/mu0 then the conversion term to SI for a diffusivity is u_el/u_i/mu0 *1.e-4 (1.e-4 for converting from cm^2 to m^2).