[Ramses] Implement Ramses like small r fix

[tdavidcl]

• Reproducer
• small r fix test

[tdavidcl]

import matplotlib.pyplot as plt
import numpy as np

import shamrock

from scipy.interpolate import griddata
from turbustat.simulator import make_3dfield

kb_cgs = 1.38064852e-16 # cm2 g s-2 K-1   # Boltzman constant
mH_cgs = 1.6737236e-24 # g                # hydrogen mass

T_init=30.0 # K
mu=2.37
units_time=3.1556926e13
units_length=3.0857e18
scale_v=units_length/units_time

# Helper function to obtain ...?
def convert_to_cell_coords(dic):
    cmin = dic["cell_min"]
    cmax = dic["cell_max"]

    xmin = []
    ymin = []
    zmin = []
    xmax = []
    ymax = []
    zmax = []

    for i in range(len(cmin)):
        m, M = cmin[i], cmax[i]

        mx, my, mz = m
        Mx, My, Mz = M

        for j in range(8):
            a, b = model.get_cell_coords(((mx, my, mz), (Mx, My, Mz)), j)

            x, y, z = a
            xmin.append(x)
            ymin.append(y)
            zmin.append(z)

            x, y, z = b
            xmax.append(x)
            ymax.append(y)
            zmax.append(z)

    dic["xmin"] = np.array(xmin)
    dic["ymin"] = np.array(ymin)
    dic["zmin"] = np.array(zmin)
    dic["xmax"] = np.array(xmax)
    dic["ymax"] = np.array(ymax)
    dic["zmax"] = np.array(zmax)

    return dic


def generate_turb_field(ncells, mach3d=1, powerlaw=3.5):
    # Mach number
    mach3d = 1

    # powerlaw index
    # empirically Gaussian fluctuations power law (modelled by the fBm) is about 3.4 - 3.6 (J.-F. Robitaille)
    powerlaw=3.5

    # random number seeds
    seed1 = 824329
    seed2 = 129862
    seed3 = 786276

    # generate fields for the three components of the velocity
    # the average velocity will be zero (no bulk motion)
    # the velocity dispersion will be 1 (before rescaling to Mach)
    vx = make_3dfield(ncells, powerlaw=powerlaw, randomseed=seed1)
    vy = make_3dfield(ncells, powerlaw=powerlaw, randomseed=seed2)
    vz = make_3dfield(ncells, powerlaw=powerlaw, randomseed=seed3)

    # calculate sound speed in code units from temperature
    sound_speed = (kb_cgs*T_init/(mu*mH_cgs))**0.5 / scale_v # code units

    # rescale velocity to Mach requested Mach number
    sigma_new = (mach3d/np.sqrt(3))*sound_speed # code units
    vx = vx * sigma_new
    vy = vy * sigma_new
    vz = vz * sigma_new

    return vx, vy, vz

''' --- Setup Shamrock with RAMSES model --- '''

ctx = shamrock.Context()
ctx.pdata_layout_new()

model = shamrock.get_Model_Ramses(context=ctx, vector_type="f64_3", grid_repr="i64_3")


''' --- Configure the solver --- '''

cfg = model.gen_default_config()

# Set grid parameters
levelmin = 7             # minimal refinement level
levelmax = 7             # maximum refinement level
boxlen = 4               # size of the box
num_cells = 2**levelmin  # total number of cells in 1 dimension
base = 2**(levelmin-1)   # number of blocks (ramses grids) in coarse grid
sz = 2                   # number of cells per block, here ye have octs like ramses
dx = boxlen / num_cells  # cell size
cfg.set_scale_factor(dx)

# EOS
gamma = 1.4
cfg.set_eos_gamma(gamma)

# Solver options
cfg.set_riemann_solver_hllc()
cfg.set_slope_lim_minmod()

# Boundary conditions
cfg.set_boundary_condition("x", "periodic")
cfg.set_boundary_condition("y", "periodic")
cfg.set_boundary_condition("z", "periodic")

# Apply solver configuration
model.set_solver_config(cfg)


''' --- Set up the work scheduler --- '''

# initialize scheduler
crit_split = int(1e7)  # threshold number of cells for dividing a patch in two
crit_merge = 1         # threshold number of cells for merging patches
model.init_scheduler(crit_split, crit_merge) # we have only 1 patch

# Create the base grid and distribute the work
bottom_corner = (0, 0, 0)
num_blocks = (base, base, base)
num_cells_per_grid = (sz, sz, sz)
model.make_base_grid(bottom_corner, num_cells_per_grid, num_blocks)


''' --- Initial conditions --- '''

# Uniform density and pressure field, so that T=~30K everywhere
rho_init = 0.28954719470909174   # density
#P_init = 0.030361026190591216    # pressure
P_init = rho_init * (kb_cgs*T_init/(mu*mH_cgs)) / scale_v**2  # pressure
print('DEBUG: P_init=', P_init)
etherm_init = P_init / (gamma - 1) # energy

# generate turbulent velocity field
vx, vy, vz = generate_turb_field(num_cells)

# Function for the initial density in the entire domain
def rho_map(rmin, rmax):
    return rho_init

# Function for the initial energy in the entire domain
def rhoetot_map(rmin, rmax):
    rho = rho_map(rmin, rmax)
    vxloc, vyloc, vzloc = rhovel_map(rmin, rmax)
    ekin = 0.5*(vxloc**2 + vyloc**2 + vzloc**2)
    return etherm_init*rho + ekin

# Function for the initial velocity in the entire domain
def rhovel_map(rmin, rmax):
    rho = rho_map(rmin, rmax)
    x, y, z = rmin
    ix = int((x / dx))
    iy = int((y / dx))
    iz = int((z / dx))
    vxloc = vx[ix, iy, iz] *rho
    vyloc = vy[ix, iy, iz] *rho
    vzloc = vz[ix, iy, iz] *rho
    return (vxloc, vyloc, vzloc)

# Set initial field values in the entire computational domain
# using the above functions which are functions of the 3D cell edge coordinates
model.set_field_value_lambda_f64("rho", rho_map)
model.set_field_value_lambda_f64("rhoetot", rhoetot_map)
model.set_field_value_lambda_f64_3("rhovel", rhovel_map)


''' --- Evolve the model until specified time --- '''

freq = 50
dt = 0.00001
t = 0
t_end = 0.3

#next_dt = model.evolve_once_override_time(t, dt)
#'''
#for i in range(1):
i = 0
while True:
    next_dt = model.evolve_once_override_time(t, dt)

    t += dt
    dt = next_dt

    if t_end < t + next_dt:
        dt = t_end - t
    if t == t_end:
        break
    i += 1
#'''

# write output at end time


''' --- Plot results with python --- '''

# Retrieve simulation results
dic = convert_to_cell_coords(ctx.collect_data())
print(dic)
print(dic.keys())
print('DEBUG', dic['cell_min'].shape)
print('DEBUG', dic['xmin'].shape)
print('DEBUG', dic['rho'].shape)
print(np.min(dic['rho']), np.max(dic['rho']), np.isnan(dic['rho']).any())


def plot_results(dic, what='projections'):

    # order cell coordinates into a regular 3D grid
    xmin, ymin, zmin = 0,0,0
    xmax, ymax, zmax = boxlen, boxlen, boxlen
    xpx = np.linspace(xmin,xmax-dx,num_cells)
    ypx = np.linspace(ymin,ymax-dx,num_cells)
    zpx = np.linspace(zmin,zmax-dx,num_cells)
    grid_x, grid_y, grid_z = np.meshgrid(xpx,ypx,zpx)
    points = np.transpose([dic["xmin"],dic["ymin"],dic["zmin"]])

    # convert 1D quantity arrays to regular 3D
    rho = griddata(points, dic["rho"], (grid_x,grid_y, grid_z),method='nearest')
    vel = [griddata(points, dic["rhovel"][:,0] / dic["rho"] * scale_v / 1e5 ,(grid_x,grid_y, grid_z),method='nearest'),
           griddata(points, dic["rhovel"][:,1] / dic["rho"] * scale_v / 1e5 ,(grid_x,grid_y, grid_z),method='nearest'),
           griddata(points, dic["rhovel"][:,2] / dic["rho"] * scale_v / 1e5 ,(grid_x,grid_y, grid_z),method='nearest')]
    rhoekin = 0.5*(dic["rhovel"][:,0]**2 + dic["rhovel"][:,1]**2 + dic["rhovel"][:,2]**2)
    p   = (dic["rhoetot"] - rhoekin)/dic["rho"] * (gamma-1)
    temperature = p / dic["rho"] * (mu*mH_cgs)/kb_cgs * scale_v**2
    p = griddata(points, p, (grid_x,grid_y, grid_z),method='nearest')
    temperature = griddata(points, temperature, (grid_x,grid_y, grid_z),method='nearest')

    # make figure
    fig, ax = plt.subplots(nrows=3, ncols=3, figsize=(12, 8))
    cb = []

    if what == 'projections':
        # plot projections along each axis
        for i in range(3):
            # plot density projections
            im = ax[0,i].imshow(np.sum(rho, axis=i), origin="lower", aspect='equal',
                                extent=[0, boxlen, 0, boxlen])
            cb.append(plt.colorbar(im, ax=ax[0,i], label='density'))
            # plot velocity projections
            im = ax[1,i].imshow(np.sum(vel[i], axis=i), origin="lower", aspect='equal',
                                extent=[0, boxlen, 0, boxlen], cmap='magma')
            cb.append(plt.colorbar(im, ax=ax[1,i], label='velocity'))
            # plot pressure or temperature projections
            im = ax[2,i].imshow(np.sum(temperature, axis=i), origin="lower", aspect='equal',
                                extent=[0, boxlen, 0, boxlen])
            cb.append(plt.colorbar(im, ax=ax[2,i], label='temperature'))

    else:
        # plot slices along each axis
        ind = int(num_cells / 2)
        im = ax[0,0].imshow(rho[ind,:,:], origin="lower", aspect='equal',
                            extent=[0, boxlen, 0, boxlen])
        cb.append(plt.colorbar(im, ax=ax[0,0], label='density [c.u.]'))
        im = ax[0,1].imshow(rho[:,ind,:], origin="lower", aspect='equal',
                            extent=[0, boxlen, 0, boxlen])
        cb.append(plt.colorbar(im, ax=ax[0,1], label='density [c.u.]'))
        im = ax[0,2].imshow(rho[:,:,ind], origin="lower", aspect='equal',
                            extent=[0, boxlen, 0, boxlen])
        cb.append(plt.colorbar(im, ax=ax[0,2], label='density [c.u.]'))
        im = ax[1,0].imshow(vel[0][ind,:,:], origin="lower", aspect='equal',
                            extent=[0, boxlen, 0, boxlen], cmap='coolwarm')
        cb.append(plt.colorbar(im, ax=ax[1,0], label='velocity [km/s]'))
        im = ax[1,1].imshow(vel[1][:,ind,:], origin="lower", aspect='equal',
                            extent=[0, boxlen, 0, boxlen], cmap='coolwarm')
        cb.append(plt.colorbar(im, ax=ax[1,1], label='velocity [km/s]'))
        im = ax[1,2].imshow(vel[2][:,:,ind], origin="lower", aspect='equal',
                            extent=[0, boxlen, 0, boxlen],  cmap='coolwarm')
        cb.append(plt.colorbar(im, ax=ax[1,2], label='velocity [km/s]'))
        im = ax[2,0].imshow(temperature[ind,:,:], origin="lower", aspect='equal',
                            extent=[0, boxlen, 0, boxlen], cmap='magma')
        cb.append(plt.colorbar(im, ax=ax[2,0], label='temperature [K]'))
        im = ax[2,1].imshow(temperature[:,ind,:], origin="lower", aspect='equal',
                            extent=[0, boxlen, 0, boxlen], cmap='magma')
        cb.append(plt.colorbar(im, ax=ax[2,1], label='temperature [K]'))
        im = ax[2,2].imshow(temperature[:,:,ind], origin="lower", aspect='equal',
                            extent=[0, boxlen, 0, boxlen], cmap='magma')
        cb.append(plt.colorbar(im, ax=ax[2,2], label='temperature [K]'))

    for i in range(3):
        ax[i,0].set_xlabel('x')
        ax[i,0].set_ylabel('z')
        ax[i,1].set_xlabel('y')
        ax[i,1].set_ylabel('z')
        ax[i,2].set_xlabel('x')
        ax[i,2].set_ylabel('y')

    for c in cb:
        c.ax.yaxis.set_label_coords(-1.1, 0.5)

    plt.savefig("decaying_turbulence.png")
    plt.show()

plot_results(dic, what='slices')

[tinecolman]

This setup works for lvl 6 but not lvl 7 or when you integrate for longer times.