[Stagging PR][Ramses] main PR for new AMR - #1942
Draft
Akos299 wants to merge 61 commits into
Draft
Conversation
This commit is the first in a series leading to a PR that updates the AMR algorithms. Specifically, it implements the 2:1 refinement consistency check.
Save the beautifull AMR
…stance between cell centers. Important for the AMR.
|
Important Review skippedDraft detected. Please check the settings in the CodeRabbit UI or the ⚙️ Run configurationConfiguration used: Repository UI Review profile: ASSERTIVE Plan: Pro Plus Run ID: You can disable this status message by setting the Use the checkbox below for a quick retry:
✨ Finishing Touches🧪 Generate unit tests (beta)
Comment |
Contributor
|
Thanks @Akos299 for opening this PR! You can do multiple things directly here: Once the workflow completes a message will appear displaying informations related to the run. Also the PR gets automatically reviewed by gemini, you can: |
Contributor
Workflow reportworkflow report corresponding to commit c93bc56 Pre-commit check reportSome failures were detected in base source checks checks. ❌ trailing-whitespace❌ remove-tabs❌ forbid-tabs❌ end-of-file-fixer❌ ruff-format❌ ruff-checkSuggested changesDetailed changes :diff --git a/examples/TO_MIGRATE/ramses/Jeans_instability.py b/examples/TO_MIGRATE/ramses/Jeans_instability.py
index 6a0f8221..23051264 100644
--- a/examples/TO_MIGRATE/ramses/Jeans_instability.py
+++ b/examples/TO_MIGRATE/ramses/Jeans_instability.py
@@ -99,8 +99,8 @@ def run_sim(rhog, vg, etot, cs, times, lembda=0.5, rho0=1, amp=1e-2, NJ=4):
gamma = 1.0000001
- k = 2*np.pi /lembda
- gamma_growth = np.sqrt(4*np.pi*G*rho0 - cs**2*k**2)
+ k = 2 * np.pi / lembda
+ gamma_growth = np.sqrt(4 * np.pi * G * rho0 - cs**2 * k**2)
### Gas maps
def rho_map(rmin, rmax) -> float:
@@ -179,17 +179,16 @@ def run_sim(rhog, vg, etot, cs, times, lembda=0.5, rho0=1, amp=1e-2, NJ=4):
freq = 500
dt = 0.000
t = 0
- k = 2*np.pi /lembda
+ k = 2 * np.pi / lembda
- lambdaJ = np.sqrt(np.pi*cs**2/(G*rho0))
+ lambdaJ = np.sqrt(np.pi * cs**2 / (G * rho0))
tend = 0
if lembda < lambdaJ:
- omega = np.sqrt(cs**2*k**2 - 4*np.pi*G*rho0)
- tend = 2.0 * (2*np.pi/omega) # two oscillation periods
+ omega = np.sqrt(cs**2 * k**2 - 4 * np.pi * G * rho0)
+ tend = 2.0 * (2 * np.pi / omega) # two oscillation periods
else:
-
- tend = 3.0 / gamma_growth # exponential growth
+ tend = 3.0 / gamma_growth # exponential growth
a = None
b = None
c = None
@@ -208,20 +207,20 @@ def run_sim(rhog, vg, etot, cs, times, lembda=0.5, rho0=1, amp=1e-2, NJ=4):
dist2 = xc**2 + yc**2 + zc**2
idx0 = np.argmin(dist2)
- vg_i = dic["rhovel"][idx0,0] / dic["rho"][idx0]
+ vg_i = dic["rhovel"][idx0, 0] / dic["rho"][idx0]
rg_i = dic["rho"][idx0]
e_i = dic["rhoetot"][idx0]
a = dic["rho"] - rho0
- b = (dic["xmin"] + dic["xmax"])
+ b = dic["xmin"] + dic["xmax"]
# b = dic["xmin"]
c = dic["rhovel"][:, 0] / (dic["rho"])
mask = mask = np.logical_and(
- dic["ymin"] == 0,
- dic["zmin"] == 0,
- )
+ dic["ymin"] == 0,
+ dic["zmin"] == 0,
+ )
rhog.append(rg_i - rho0)
vg.append(vg_i)
@@ -231,14 +230,13 @@ def run_sim(rhog, vg, etot, cs, times, lembda=0.5, rho0=1, amp=1e-2, NJ=4):
t += dt
dt = next_dt
- if not (lembda < lambdaJ)and np.max(dic["rho"]) > 2.0 * rho0:
+ if not (lembda < lambdaJ) and np.max(dic["rho"]) > 2.0 * rho0:
break
if tend < t + next_dt:
dt = tend - t
if t == tend:
break
-
return a, c, b, mask
@@ -289,46 +287,179 @@ def plot_analytical_solution_col_snapshots(A, k, rho0, Lambd, gam_lambd, positio
L = 2.0
lembda = L
amp = 1e-2
-rho0 = 1.
+rho0 = 1.0
G = 1.0
-cs_list_col = [0.1 ,0.15188784, 0.18302054, 0.21415324,
- 0.22453081, 0.23490838, 0.26604108, 0.27641864, 0.28679621, 0.32830648,
- 0.36981675, 0.41132702, 0.45283729, 0.50472513, 0.55661296, 0.56699053, 0.5773681 , 0.61887837, 0.62925593, 0.6396335, 0.65001107,
- 0.66038864, 0.6707662, 0.72265404, 0.73303161, 0.77454188, 0.78491944, 0.79529701, 0.83680728, 0.84718485, 0.85756242, 0.89907269,
- 0.90945025, 0.94058295, 0.97171566, 1.01322593, 1.04435863, 1.0547362, 1.0858689 ,
- 1.09624647, 1.10662403, 1.1170016, 1.12737917 ]
-
-
-cs_list_osc = [1.16837917, 1.20837917, 1.40837917, 1.56837917, 1.68837917, 1.80837917,
- 1.96837917, 2.00837917, 2.04837917, 2.08837917, 2.12837917,2.28837917, 2.44837917, 2.56837917,
- 2.72837917, 2.76837917, 2.88837917, 3.04837917, 3.16837917,
- 3.32837917, 3.48837917, 3.60837917, 3.76837917, 3.92837917, 3.96837917, 4.08837917, 4.24837917,
- 4.40837917, 4.64837917, 4.68837917, 4.80837917, 4.96837917,
- 5.00837917, 5.12837917]
-
-
-cs_list_glob = [0.1 ,0.15188784, 0.18302054, 0.21415324,
- 0.22453081, 0.23490838, 0.26604108, 0.27641864, 0.28679621, 0.32830648,
- 0.36981675, 0.41132702, 0.45283729, 0.50472513, 0.55661296, 0.56699053, 0.5773681 , 0.61887837, 0.62925593, 0.6396335, 0.65001107,
- 0.66038864, 0.6707662, 0.72265404, 0.73303161, 0.77454188, 0.78491944, 0.79529701, 0.83680728, 0.84718485, 0.85756242, 0.89907269,
- 0.90945025, 0.94058295, 0.97171566, 1.01322593, 1.04435863, 1.0547362, 1.0858689 ,
- 1.09624647, 1.10662403, 1.1170016, 1.12737917 ,
-
- 1.16837917, 1.20837917, 1.40837917, 1.56837917, 1.68837917, 1.80837917,
- 1.96837917, 2.00837917, 2.04837917, 2.08837917, 2.12837917,2.28837917, 2.44837917, 2.56837917,
- 2.72837917, 2.76837917, 2.88837917, 3.04837917, 3.16837917,
- 3.32837917, 3.48837917, 3.60837917, 3.76837917, 3.92837917, 3.96837917, 4.08837917, 4.24837917,
- 4.40837917, 4.64837917, 4.68837917, 4.80837917, 4.96837917,
- 5.00837917, 5.12837917]
+cs_list_col = [
+ 0.1,
+ 0.15188784,
+ 0.18302054,
+ 0.21415324,
+ 0.22453081,
+ 0.23490838,
+ 0.26604108,
+ 0.27641864,
+ 0.28679621,
+ 0.32830648,
+ 0.36981675,
+ 0.41132702,
+ 0.45283729,
+ 0.50472513,
+ 0.55661296,
+ 0.56699053,
+ 0.5773681,
+ 0.61887837,
+ 0.62925593,
+ 0.6396335,
+ 0.65001107,
+ 0.66038864,
+ 0.6707662,
+ 0.72265404,
+ 0.73303161,
+ 0.77454188,
+ 0.78491944,
+ 0.79529701,
+ 0.83680728,
+ 0.84718485,
+ 0.85756242,
+ 0.89907269,
+ 0.90945025,
+ 0.94058295,
+ 0.97171566,
+ 1.01322593,
+ 1.04435863,
+ 1.0547362,
+ 1.0858689,
+ 1.09624647,
+ 1.10662403,
+ 1.1170016,
+ 1.12737917,
+]
+
+
+cs_list_osc = [
+ 1.16837917,
+ 1.20837917,
+ 1.40837917,
+ 1.56837917,
+ 1.68837917,
+ 1.80837917,
+ 1.96837917,
+ 2.00837917,
+ 2.04837917,
+ 2.08837917,
+ 2.12837917,
+ 2.28837917,
+ 2.44837917,
+ 2.56837917,
+ 2.72837917,
+ 2.76837917,
+ 2.88837917,
+ 3.04837917,
+ 3.16837917,
+ 3.32837917,
+ 3.48837917,
+ 3.60837917,
+ 3.76837917,
+ 3.92837917,
+ 3.96837917,
+ 4.08837917,
+ 4.24837917,
+ 4.40837917,
+ 4.64837917,
+ 4.68837917,
+ 4.80837917,
+ 4.96837917,
+ 5.00837917,
+ 5.12837917,
+]
+
+
+cs_list_glob = [
+ 0.1,
+ 0.15188784,
+ 0.18302054,
+ 0.21415324,
+ 0.22453081,
+ 0.23490838,
+ 0.26604108,
+ 0.27641864,
+ 0.28679621,
+ 0.32830648,
+ 0.36981675,
+ 0.41132702,
+ 0.45283729,
+ 0.50472513,
+ 0.55661296,
+ 0.56699053,
+ 0.5773681,
+ 0.61887837,
+ 0.62925593,
+ 0.6396335,
+ 0.65001107,
+ 0.66038864,
+ 0.6707662,
+ 0.72265404,
+ 0.73303161,
+ 0.77454188,
+ 0.78491944,
+ 0.79529701,
+ 0.83680728,
+ 0.84718485,
+ 0.85756242,
+ 0.89907269,
+ 0.90945025,
+ 0.94058295,
+ 0.97171566,
+ 1.01322593,
+ 1.04435863,
+ 1.0547362,
+ 1.0858689,
+ 1.09624647,
+ 1.10662403,
+ 1.1170016,
+ 1.12737917,
+ 1.16837917,
+ 1.20837917,
+ 1.40837917,
+ 1.56837917,
+ 1.68837917,
+ 1.80837917,
+ 1.96837917,
+ 2.00837917,
+ 2.04837917,
+ 2.08837917,
+ 2.12837917,
+ 2.28837917,
+ 2.44837917,
+ 2.56837917,
+ 2.72837917,
+ 2.76837917,
+ 2.88837917,
+ 3.04837917,
+ 3.16837917,
+ 3.32837917,
+ 3.48837917,
+ 3.60837917,
+ 3.76837917,
+ 3.92837917,
+ 3.96837917,
+ 4.08837917,
+ 4.24837917,
+ 4.40837917,
+ 4.64837917,
+ 4.68837917,
+ 4.80837917,
+ 4.96837917,
+ 5.00837917,
+ 5.12837917,
+]
cs_list_glob = [1.11]
-
for cs in cs_list_glob:
-
times = []
rg_num = []
vg_num = []
@@ -346,7 +477,6 @@ for cs in cs_list_glob:
rg_num, vg_num, etot_num, cs, times, lembda, rho0, amp, NJ=N_J
)
-
if shamrock.sys.world_rank() == 0:
# get indexes X at (Y,Z)=(0,0)
ind = np.where(mask)[0]
@@ -356,8 +486,6 @@ for cs in cs_list_glob:
rho_last = rho_last[ind]
vel_last = vel_last[ind]
-
-
######
order = np.argsort(X)
X = X[order]
@@ -365,8 +493,6 @@ for cs in cs_list_glob:
vel_last = vel_last[order]
#####
-
-
times = np.array(times)
x0 = X[0]
t_last = times[-1]
@@ -401,17 +527,15 @@ for cs in cs_list_glob:
datas_times,
)
-
datas_spaces = np.stack((X, rho_last, dens_fix_time, vel_last, vel_fix_time)).T
np.savetxt(
f"_space_evolution_Jeans-instablity_A_{amp:.3f}_Cs_{cs:.6f}_Rho_0_{rho0}_Lambd_{lembda}_X_sz_{len(X)}_Rhog_sz_{len(rg_num)}.txt",
datas_spaces,
)
-
- #--------------------------------------
- # Plots
- #-------------------------------------
+ # --------------------------------------
+ # Plots
+ # -------------------------------------
fig, axs = plt.subplots(2, 2, figsize=(8, 8))
plt.subplots_adjust(wspace=0.25)
axs[0][0].plot(times, rg_num, "co", label="$\\rho_{num}$")
@@ -463,4 +587,3 @@ for cs in cs_list_glob:
plt.savefig(
f"Jeans_Instability_test_2_07_2026_{cs}.pdf",
)
-
diff --git a/examples/TO_MIGRATE/ramses/Kelvin_Helmholtz.py b/examples/TO_MIGRATE/ramses/Kelvin_Helmholtz.py
index 46f74424..6266b602 100644
--- a/examples/TO_MIGRATE/ramses/Kelvin_Helmholtz.py
+++ b/examples/TO_MIGRATE/ramses/Kelvin_Helmholtz.py
@@ -3,34 +3,34 @@ import os
import matplotlib.pyplot as plt
import numpy as np
-
-plt.rcParams.update({
- "font.size": 15,
- "axes.labelsize": 20,
- "axes.titlesize": 15,
- "xtick.labelsize": 15,
- "ytick.labelsize": 15,
- "legend.fontsize": 15,
-})
+plt.rcParams.update(
+ {
+ "font.size": 15,
+ "axes.labelsize": 20,
+ "axes.titlesize": 15,
+ "xtick.labelsize": 15,
+ "ytick.labelsize": 15,
+ "legend.fontsize": 15,
+ }
+)
import matplotlib as mpl
-
-mpl.rcParams.update({
- "text.usetex": True, # Use LaTeX
- "font.family": "serif",
- "font.serif": ["Latin Modern Roman"], # Match lmodern
-
- # LaTeX preamble to match your class
- "text.latex.preamble": r"""
+mpl.rcParams.update(
+ {
+ "text.usetex": True, # Use LaTeX
+ "font.family": "serif",
+ "font.serif": ["Latin Modern Roman"], # Match lmodern
+ # LaTeX preamble to match your class
+ "text.latex.preamble": r"""
\usepackage{lmodern}
\usepackage{amsmath}
\usepackage{amssymb}
""",
-
- # Optional but recommended
- "axes.unicode_minus": False
-})
+ # Optional but recommended
+ "axes.unicode_minus": False,
+ }
+)
mpl.rcParams["pgf.texsystem"] = "pdflatex"
@@ -41,12 +41,11 @@ if not shamrock.sys.is_initialized():
shamrock.sys.init("0:0")
-sim_folder = f"_to_trash/amr_athena_kelvin_helmholtz/"
+sim_folder = "_to_trash/amr_athena_kelvin_helmholtz/"
if shamrock.sys.world_rank() == 0:
os.makedirs(sim_folder, exist_ok=True)
-
# Utility for plotting
def make_cartesian_coords(nx, ny, z_val, min_x, max_x, min_y, max_y):
# Create the cylindrical coordinate grid
@@ -69,6 +68,7 @@ def make_cartesian_coords(nx, ny, z_val, min_x, max_x, min_y, max_y):
# Plot utility
################################################################
+
def plot_kh_density(ext, time, rho, nx, ny, dpi=200):
rho = np.asarray(rho).reshape(nx, ny)
plt.figure(figsize=(6, 6), dpi=dpi)
@@ -79,7 +79,7 @@ def plot_kh_density(ext, time, rho, nx, ny, dpi=200):
extent=ext,
cmap="jet",
vmin=1.0,
- vmax=2.,
+ vmax=2.0,
interpolation="nearest",
aspect="equal",
)
@@ -87,31 +87,21 @@ def plot_kh_density(ext, time, rho, nx, ny, dpi=200):
plt.xlabel("x")
plt.ylabel("y")
- plt.title(
- rf"$t = {time:.3f}$"
- )
+ plt.title(rf"$t = {time:.3f}$")
cbar = plt.colorbar(im)
cbar.set_label(r"$\rho$")
plt.tight_layout()
- plt.savefig(
- os.path.join(
- sim_folder,
- f"amr_kh_density_resx_{nx}_resy_{ny}_at_{time:.3f}.pdf"
- )
- )
+ plt.savefig(os.path.join(sim_folder, f"amr_kh_density_resx_{nx}_resy_{ny}_at_{time:.3f}.pdf"))
plt.close()
-
-
##########################################################
# Main routine for simulation
##########################################################
-
def run_simulation(output_freq, t_final, extent, base, multx, multy, multz, dens_jump):
@@ -121,11 +111,10 @@ def run_simulation(output_freq, t_final, extent, base, multx, multy, multz, dens
sz = 2 << amr_lev
# sz = 1<<1
base = base
- scale_fact = 1. / (sz * base * multx)
+ scale_fact = 1.0 / (sz * base * multx)
nx = base * sz * multx
ny = base * sz * multy
- positions = make_cartesian_coords(nx, ny, 0.2, 0, 1. - 1e-6, 0, 1. - 1e-6)
-
+ positions = make_cartesian_coords(nx, ny, 0.2, 0, 1.0 - 1e-6, 0, 1.0 - 1e-6)
####---------------
shamrock.enable_experimental_features()
@@ -145,12 +134,9 @@ def run_simulation(output_freq, t_final, extent, base, multx, multy, multz, dens
cfg.set_boundary_condition("y", "periodic")
cfg.set_boundary_condition("z", "periodic")
-
-
thre_s = 0.01
cfg.set_amr_mode_shear_based(Threshold=thre_s)
-
# err_min = 0.05
# err_max = 0.10
# cfg.set_amr_mode_pseudo_gradient_based(error_min=err_min, error_max=err_max)
@@ -160,7 +146,6 @@ def run_simulation(output_freq, t_final, extent, base, multx, multy, multz, dens
###------------
-
# P0 = 2.5
# y0_ref= 0.25
# y0 = 0.5 + y0_ref
@@ -168,18 +153,13 @@ def run_simulation(output_freq, t_final, extent, base, multx, multy, multz, dens
# sig = 0.2
# L = 0.005
-
# P0 = 10
# y1= 0.5
# y2 = 1.5
# A = 0.01
# sig = 0.2
# L = 0.05
- u_flow=1
-
-
-
-
+ u_flow = 1
P0 = 2.5
@@ -191,16 +171,14 @@ def run_simulation(output_freq, t_final, extent, base, multx, multy, multz, dens
vflow = 1.0
-
def rho_map(rmin, rmax):
_, y, _ = rmin
yloc = abs(y - 0.5)
- w = 0.5*(np.tanh((yloc - 0.25)/L) + 1.0)
-
- return w + (1.0 - w)*dens_jump
+ w = 0.5 * (np.tanh((yloc - 0.25) / L) + 1.0)
+ return w + (1.0 - w) * dens_jump
# def rho_map(rmin, rmax):
# _, y, _ = rmin
@@ -208,28 +186,21 @@ def run_simulation(output_freq, t_final, extent, base, multx, multy, multz, dens
# return 1.5 - 0.5 * np.tanh(np.abs(y - y0) / L)
# # return 1.0 + tmp_rho
-
-
def rhovel_map(rmin, rmax):
x, y, _ = rmin
- yloc = abs(y-0.5)
-
- w = 0.5*(np.tanh((yloc-0.25)/L)+1.0)
+ yloc = abs(y - 0.5)
- rho = w + (1-w)*dens_jump
+ w = 0.5 * (np.tanh((yloc - 0.25) / L) + 1.0)
- vx = (w - (1-w)*dens_jump)/rho
+ rho = w + (1 - w) * dens_jump
- vy = (
- A
- * np.cos(4*np.pi*x)
- * np.exp(-(yloc-0.25)**2/(sig*sig))
- )
+ vx = (w - (1 - w) * dens_jump) / rho
- return (rho*vx, rho*vy, 0)
+ vy = A * np.cos(4 * np.pi * x) * np.exp(-((yloc - 0.25) ** 2) / (sig * sig))
+ return (rho * vx, rho * vy, 0)
# def rhovel_map(rmin, rmax):
# # vx = 0.5 * np.tanh(np.abs(y - y0)/L)
@@ -241,10 +212,8 @@ def run_simulation(output_freq, t_final, extent, base, multx, multy, multz, dens
# vx = u_flow * tmp_vx
# vy = A*np.sin(2*np.pi*x) * tmp_vy
-
# return (rho*vx, vy*rho, 0)
-
# def rhoetot_map(rmin, rmax):
# # vx = 0.5 * np.tanh(np.abs(y - y0)/L)
# # vy = A * np.cos(4*np.pi * x)*np.exp(-(y-y0)**2/sig**2)
@@ -254,29 +223,25 @@ def run_simulation(output_freq, t_final, extent, base, multx, multy, multz, dens
# tmp_vy = np.exp(-(y-y1)**2/sig**2) + np.exp(-(y-y2)**2/sig**2)
# vx = u_flow * tmp_vx
# vy = A*np.sin(2*np.pi*x) * tmp_vy
- # P = P0
+ # P = P0
# Eint = P/(gamma - 1.0)
# Ekin = 0.5*rho*(vx**2 + vy**2)
# return Ekin + Eint
-
-
def rhoetot_map(rmin, rmax):
- rho = rho_map(rmin,rmax)
-
- mx,my,mz = rhovel_map(rmin,rmax)
+ rho = rho_map(rmin, rmax)
- kinetic = 0.5*(mx*mx+my*my+mz*mz)/rho
+ mx, my, mz = rhovel_map(rmin, rmax)
- return kinetic + P0/(gamma-1)
+ kinetic = 0.5 * (mx * mx + my * my + mz * mz) / rho
+ return kinetic + P0 / (gamma - 1)
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)
-
#####------------------------
freq = output_freq
dt = 0.0000
@@ -287,28 +252,34 @@ def run_simulation(output_freq, t_final, extent, base, multx, multy, multz, dens
model.dump_vtk("amr_athena_kelvin_helmhotz" + str(i // freq) + ".vtk")
next_dt = model.evolve_once_override_time(t, dt)
-
+
t += dt
dt = next_dt
current_time = t
-
if i % freq == 0:
- arr_rho_pos = model.render_slice("rho", "f64", positions)
- rhov_vals = model.render_slice("rhovel", "f64_3", positions)
- rhoetot_vals = model.render_slice("rhoetot", "f64", positions)
- vx = np.array(rhov_vals)[:,0] / np.array(arr_rho_pos)
- vy = np.array(rhov_vals)[:,1] / np.array(arr_rho_pos)
- P = (np.array(rhoetot_vals) - 0.5 * np.array(arr_rho_pos) * (vx**2 + vy**2))*(gamma - 1.)
-
- output = np.column_stack((np.array(arr_rho_pos), np.array(vx), np.array(vy), np.array(P)))
- filename = f"amr_kh_datas_base_{base}_lev_max_{amr_lev}_nx_{nx}_ny_{ny}_at_{current_time}.txt"
- np.savetxt(os.path.join(sim_folder,filename),
- output,
- fmt=["%.10f", "%.10f", "%.10f", "%.10f"],
- header="rho vx vy P",
- )
- plot_kh_density(extent,current_time,arr_rho_pos,nx,ny)
+ arr_rho_pos = model.render_slice("rho", "f64", positions)
+ rhov_vals = model.render_slice("rhovel", "f64_3", positions)
+ rhoetot_vals = model.render_slice("rhoetot", "f64", positions)
+ vx = np.array(rhov_vals)[:, 0] / np.array(arr_rho_pos)
+ vy = np.array(rhov_vals)[:, 1] / np.array(arr_rho_pos)
+ P = (np.array(rhoetot_vals) - 0.5 * np.array(arr_rho_pos) * (vx**2 + vy**2)) * (
+ gamma - 1.0
+ )
+
+ output = np.column_stack(
+ (np.array(arr_rho_pos), np.array(vx), np.array(vy), np.array(P))
+ )
+ filename = (
+ f"amr_kh_datas_base_{base}_lev_max_{amr_lev}_nx_{nx}_ny_{ny}_at_{current_time}.txt"
+ )
+ np.savetxt(
+ os.path.join(sim_folder, filename),
+ output,
+ fmt=["%.10f", "%.10f", "%.10f", "%.10f"],
+ header="rho vx vy P",
+ )
+ plot_kh_density(extent, current_time, arr_rho_pos, nx, ny)
if t_final < t + next_dt:
dt = t_final - t
@@ -316,7 +287,6 @@ def run_simulation(output_freq, t_final, extent, base, multx, multy, multz, dens
break
-
#####-----------------
multz = 1
@@ -332,8 +302,8 @@ multx = 1
multy = 1
-extent = [0, 1., 0, 1.]
+extent = [0, 1.0, 0, 1.0]
out_freq = 10
tend = 5
-dens_jump =1
-run_simulation(out_freq, tend, extent, base, multx, multy, multz, dens_jump)
\ No newline at end of file
+dens_jump = 1
+run_simulation(out_freq, tend, extent, base, multx, multy, multz, dens_jump)
diff --git a/examples/TO_MIGRATE/ramses/Rayleigh_Taylor.py b/examples/TO_MIGRATE/ramses/Rayleigh_Taylor.py
index 50d84106..129b1e2c 100644
--- a/examples/TO_MIGRATE/ramses/Rayleigh_Taylor.py
+++ b/examples/TO_MIGRATE/ramses/Rayleigh_Taylor.py
@@ -3,34 +3,34 @@ import os
import matplotlib.pyplot as plt
import numpy as np
-
-plt.rcParams.update({
- "font.size": 15,
- "axes.labelsize": 20,
- "axes.titlesize": 15,
- "xtick.labelsize": 15,
- "ytick.labelsize": 15,
- "legend.fontsize": 15,
-})
+plt.rcParams.update(
+ {
+ "font.size": 15,
+ "axes.labelsize": 20,
+ "axes.titlesize": 15,
+ "xtick.labelsize": 15,
+ "ytick.labelsize": 15,
+ "legend.fontsize": 15,
+ }
+)
import matplotlib as mpl
-
-mpl.rcParams.update({
- "text.usetex": True, # Use LaTeX
- "font.family": "serif",
- "font.serif": ["Latin Modern Roman"], # Match lmodern
-
- # LaTeX preamble to match your class
- "text.latex.preamble": r"""
+mpl.rcParams.update(
+ {
+ "text.usetex": True, # Use LaTeX
+ "font.family": "serif",
+ "font.serif": ["Latin Modern Roman"], # Match lmodern
+ # LaTeX preamble to match your class
+ "text.latex.preamble": r"""
\usepackage{lmodern}
\usepackage{amsmath}
\usepackage{amssymb}
""",
-
- # Optional but recommended
- "axes.unicode_minus": False
-})
+ # Optional but recommended
+ "axes.unicode_minus": False,
+ }
+)
mpl.rcParams["pgf.texsystem"] = "pdflatex"
@@ -41,12 +41,11 @@ if not shamrock.sys.is_initialized():
shamrock.sys.init("0:0")
-sim_folder = f"_to_trash/rayleigh_taylor_vl/amr/"
+sim_folder = "_to_trash/rayleigh_taylor_vl/amr/"
if shamrock.sys.world_rank() == 0:
os.makedirs(sim_folder, exist_ok=True)
-
# Utility for plotting
def make_cartesian_coords(nx, ny, z_val, min_x, max_x, min_y, max_y):
# Create the cylindrical coordinate grid
@@ -69,6 +68,7 @@ def make_cartesian_coords(nx, ny, z_val, min_x, max_x, min_y, max_y):
# Plot utility
################################################################
+
def plot_rt_density(ext, time, rho, nx, ny, dpi=200):
rho = np.asarray(rho).reshape(nx, ny)
plt.figure(figsize=(6, 12), dpi=dpi)
@@ -85,44 +85,34 @@ def plot_rt_density(ext, time, rho, nx, ny, dpi=200):
# )
im = plt.imshow(
- rho.T,
- origin="lower",
- extent=ext,
- cmap="rainbow",
- interpolation="bicubic",
- vmin=1.0,
- vmax=2.0,
- aspect="equal",
- )
+ rho.T,
+ origin="lower",
+ extent=ext,
+ cmap="rainbow",
+ interpolation="bicubic",
+ vmin=1.0,
+ vmax=2.0,
+ aspect="equal",
+ )
plt.xlabel("x")
plt.ylabel("y")
- plt.title(
- rf"$t = {time:.3f}$"
- )
+ plt.title(rf"$t = {time:.3f}$")
cbar = plt.colorbar(im)
cbar.set_label(r"$\rho$")
plt.tight_layout()
- plt.savefig(
- os.path.join(
- sim_folder,
- f"rt_density_resx_{nx}_resy_{ny}_at_{time:.3f}.pdf"
- )
- )
+ plt.savefig(os.path.join(sim_folder, f"rt_density_resx_{nx}_resy_{ny}_at_{time:.3f}.pdf"))
plt.close()
-
-
##########################################################
# Main routine for simulation
##########################################################
-
def run_simulation(output_freq, t_final, extent, base, multx, multy, multz):
@@ -137,7 +127,6 @@ def run_simulation(output_freq, t_final, extent, base, multx, multy, multz):
ny = base * sz * multy
positions = make_cartesian_coords(nx, ny, 0.2, 0, 0.5 - 1e-6, 0, 1.5 - 1e-6)
-
####---------------
shamrock.enable_experimental_features()
ctx = shamrock.Context()
@@ -153,7 +142,7 @@ def run_simulation(output_freq, t_final, extent, base, multx, multy, multz):
cfg.set_slope_lim_vanleer_sym()
cfg.set_face_time_interpolation(True)
cfg.set_status_grav_acc(True)
- cfg.set_constant_grav_acc(gx = 0., gy =-g, gz = 0.)
+ cfg.set_constant_grav_acc(gx=0.0, gy=-g, gz=0.0)
cfg.set_boundary_condition("x", "periodic")
cfg.set_boundary_condition("y", "reflective")
cfg.set_boundary_condition("z", "reflective")
@@ -165,45 +154,50 @@ def run_simulation(output_freq, t_final, extent, base, multx, multy, multz):
model.make_base_grid((0, 0, 0), (sz, sz, sz), (base * multx, base * multy, base * multz))
###------------
- y_half = 1.5 / 2.
+ y_half = 1.5 / 2.0
P0 = 2.5
-
def rho_map(rmin, rmax):
_, y, _ = rmin
yloc = y - y_half
- if yloc <= 0.:
- return 1.
+ if yloc <= 0.0:
+ return 1.0
else:
- return 2.
-
+ return 2.0
def rhovel_map(rmin, rmax):
rho = rho_map(rmin, rmax)
- x,y,_ =rmin
+ x, y, _ = rmin
xloc = x - 0.25
yloc = y - y_half
- vy = 1e-2 * (1./4.)*(1. + np.cos(4. * np.pi * xloc))*(1. + np.cos(3. * np.pi * yloc))
- return (0, vy*rho, 0)
-
+ vy = (
+ 1e-2
+ * (1.0 / 4.0)
+ * (1.0 + np.cos(4.0 * np.pi * xloc))
+ * (1.0 + np.cos(3.0 * np.pi * yloc))
+ )
+ return (0, vy * rho, 0)
def rhoetot_map(rmin, rmax):
rho = rho_map(rmin, rmax)
x, y, _ = rmin
- P = P0 -g*rho*(y-y_half)
+ P = P0 - g * rho * (y - y_half)
xloc = x - 0.25
yloc = y - y_half
- vy = 1e-2 * (1./4.)*(1. + np.cos(4. * np.pi * xloc))*(1. + np.cos(3. * np.pi * yloc))
- Eint = P/(gamma - 1.0)
- Ekin = 0.5*rho*(vy**2)
+ vy = (
+ 1e-2
+ * (1.0 / 4.0)
+ * (1.0 + np.cos(4.0 * np.pi * xloc))
+ * (1.0 + np.cos(3.0 * np.pi * yloc))
+ )
+ Eint = P / (gamma - 1.0)
+ Ekin = 0.5 * rho * (vy**2)
return Ekin + Eint
-
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)
-
#####------------------------
freq = output_freq
dt = 0.0000
@@ -214,15 +208,14 @@ def run_simulation(output_freq, t_final, extent, base, multx, multy, multz):
model.dump_vtk("Rayleigh_Taylor" + str(i // freq) + ".vtk")
next_dt = model.evolve_once_override_time(t, dt)
-
+
t += dt
dt = next_dt
current_time = t + dt
arr_rho_pos = model.render_slice("rho", "f64", positions)
if i % freq == 0:
- plot_rt_density(extent,current_time,arr_rho_pos,nx,ny)
-
+ plot_rt_density(extent, current_time, arr_rho_pos, nx, ny)
if t_final < t + next_dt:
dt = t_final - t
@@ -232,7 +225,6 @@ def run_simulation(output_freq, t_final, extent, base, multx, multy, multz):
break
-
#####-----------------
multz = 1
@@ -248,7 +240,7 @@ multx = 1
multy = 3
-extent = [0, 0.5, 0, 1.5]
+extent = [0, 0.5, 0, 1.5]
out_freq = 50
tend = 5
-run_simulation(out_freq, tend, extent, base, multx, multy, multz)
\ No newline at end of file
+run_simulation(out_freq, tend, extent, base, multx, multy, multz)
diff --git a/examples/TO_MIGRATE/ramses/collapse.py b/examples/TO_MIGRATE/ramses/collapse.py
index 0c02bd6a..d8dea7f3 100644
--- a/examples/TO_MIGRATE/ramses/collapse.py
+++ b/examples/TO_MIGRATE/ramses/collapse.py
@@ -5,8 +5,6 @@ import matplotlib.animation as animation
import matplotlib.pyplot as plt
import numpy as np
-
-
import shamrock
if not shamrock.sys.is_initialized():
@@ -61,7 +59,7 @@ min_reso = (L0 * N_J) / (lamb_J)
print(f"min reso = {min_reso}\n")
gamma = 5.0 / 3.0
-rho_c = 2.7e-11 * 1e3 # [g/cm^3 -> kg/m^3]
+rho_c = 2.7e-11 * 1e3 # [g/cm^3 -> kg/m^3]
def run_sim():
@@ -132,11 +130,11 @@ def run_sim():
def rhoe_map(rmin, rmax):
rho = rho_map(rmin, rmax)
rhov = rhovel_map(rmin, rmax)
- Ekin = 0.5 * (rhov[0]**2)/rho
+ Ekin = 0.5 * (rhov[0] ** 2) / rho
x = rho / rho_c
- P = cs_sqr * rho * (1. + x**(2./3.))
- Eint = P / (gamma - 1.0)
- return Ekin + Eint
+ P = cs_sqr * rho * (1.0 + x ** (2.0 / 3.0))
+ Eint = P / (gamma - 1.0)
+ return Ekin + Eint
model.set_field_value_lambda_f64("rho", rho_map)
model.set_field_value_lambda_f64("rhoetot", rhoe_map)
@@ -154,12 +152,12 @@ def run_sim():
dt = next_dt
if i % freq == 0:
- model.dump_vtk(f"_iso_collapse_{t/t_ff:5f}.vtk")
+ model.dump_vtk(f"_iso_collapse_{t / t_ff:5f}.vtk")
if tmax < t + next_dt:
dt = tmax - t
if t == tmax:
- model.dump_vtk(f"_iso_collapse{t/t_ff:5f}.vtk")
+ model.dump_vtk(f"_iso_collapse{t / t_ff:5f}.vtk")
break
diff --git a/examples/TO_MIGRATE/ramses/godunov_sod_test_with_analytical_solution.py b/examples/TO_MIGRATE/ramses/godunov_sod_test_with_analytical_solution.py
index ba740a5e..057f7139 100644
--- a/examples/TO_MIGRATE/ramses/godunov_sod_test_with_analytical_solution.py
+++ b/examples/TO_MIGRATE/ramses/godunov_sod_test_with_analytical_solution.py
@@ -1,4 +1,3 @@
-
import argparse
import glob
import os
@@ -7,28 +6,22 @@ import sys
import matplotlib.pyplot as plt
import numpy as np
-
import shamrock
-
-
-
######################################################
## Analytical solution object Adapted form
#####################################################
+def run_numerical_simulation(
+ max_amr_lev, multx, multy, multz, base, L, gamma, rho0, E0, alpha0, P0_out=1e-3, with_amr=0
+):
-def run_numerical_simulation(max_amr_lev, multx, multy, multz, base, L, gamma, rho0, E0, alpha0, P0_out = 1e-3, with_amr=0 ):
-
- #------------------ Context -----------------------------
+ # ------------------ Context -----------------------------
ctx = shamrock.Context()
ctx.pdata_layout_new()
model = shamrock.get_Model_Ramses(context=ctx, vector_type="f64_3", grid_repr="i64_3")
- #--------------------------------------------------------
-
-
-
+ # --------------------------------------------------------
# multx = 1
# multy = 1
@@ -38,7 +31,7 @@ def run_numerical_simulation(max_amr_lev, multx, multy, multz, base, L, gamma, r
scale_fact = L / (cell_size * base * multx)
# base = 16
- #----------------------------------- Model configuration ------------------
+ # ----------------------------------- Model configuration ------------------
cfg = model.gen_default_config()
cfg.set_scale_factor(scale_fact)
cfg.set_eos_gamma(gamma)
@@ -51,13 +44,12 @@ def run_numerical_simulation(max_amr_lev, multx, multy, multz, base, L, gamma, r
cfg.set_slope_lim_minmod()
cfg.set_face_time_interpolation(True)
-
### Radius of 1 cell (want set the explosion in a single point )
Rstart = scale_fact
# gamma = 5.0 / 3.0
### amr configuration
- if(with_amr):
+ if with_amr:
err_min = 0.30
err_max = 0.10
cfg.set_amr_mode_pseudo_gradient_based(error_min=err_min, error_max=err_max)
@@ -66,18 +58,15 @@ def run_numerical_simulation(max_amr_lev, multx, multy, multz, base, L, gamma, r
### split factors
model.init_scheduler(int(1e7), 1)
- ### set grid
+ ### set grid
model.make_base_grid(
(0, 0, 0), (cell_size, cell_size, cell_size), (base * multx, base * multy, base * multz)
)
-
def rho_map(rmin, rmax):
return 1.0
# return rho0
-
-
####--------pre-check for the Energy setup --------------
dx = scale_fact
@@ -93,13 +82,13 @@ def run_numerical_simulation(max_amr_lev, multx, multy, multz, base, L, gamma, r
Nz = cell_size * base * multz
for k in range(Nz):
- z = (k + 0.5) * dx - L/2.
+ z = (k + 0.5) * dx - L / 2.0
for j in range(Ny):
- y = (j + 0.5) * dx - L/2.
+ y = (j + 0.5) * dx - L / 2.0
for i in range(Nx):
- x = (i + 0.5) * dx - L/2.
+ x = (i + 0.5) * dx - L / 2.0
- r = np.sqrt(x*x + y*y + z*z)
+ r = np.sqrt(x * x + y * y + z * z)
if r < Rstart:
Ncells += 1
@@ -116,7 +105,6 @@ def run_numerical_simulation(max_amr_lev, multx, multy, multz, base, L, gamma, r
print(f"Injected pressure = {Pin:.6e}")
print(f"Injected energy density = {rhoe_in:.6e}")
-
# def rhoe_map(rmin, rmax):
# x_min, y_min, z_min = rmin
# x_max, y_max, z_max = rmax
@@ -134,8 +122,7 @@ def run_numerical_simulation(max_amr_lev, multx, multy, multz, base, L, gamma, r
# return rhoe_in
# else:
# return P0_out/(gamma -1.)
-
... truncated ...
|
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
No description provided.