[AMR][RefinementHandler] Main PR with sod and amr tests

examples/TO_MIGRATE/ramses/interacting_blast_wave.py

@@ -0,0 +1,153 @@
+import os
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+import shamrock
+
+# If we use the shamrock executable to run this script instead of the python interpreter,
+# we should not initialize the system as the shamrock executable needs to handle specific MPI logic
+if not shamrock.sys.is_initialized():
+    shamrock.change_loglevel(1)
+    shamrock.sys.init("0:0")
+
+
+def run_1d_interacting_blast_wave_test(
+    nblocks,
+    nxbox,
+    nybox,
+):
+    timestamps = 100
+    gamma = 1.4
+
+    output_folder = "_to_trash/interacting_blast_wave/"
+    os.makedirs(output_folder, exist_ok=True)
+
+    multx = nxbox
+    multy = nybox
+    multz = nybox
+    sz = 1 << 1
+    base = nblocks
+    scale_fact = 1 / (sz * base * multx)
+
+    rez_plot = 256
+    positions_plot = [(x, 0, 0) for x in np.linspace(0, 1, rez_plot).tolist()[:-1]]
+
+    ctx = shamrock.Context()
+    ctx.pdata_layout_new()
+
+    model = shamrock.get_Model_Ramses(context=ctx, vector_type="f64_3", grid_repr="i64_3")
+
+    cfg = model.gen_default_config()
+    cfg.set_scale_factor(scale_fact)
+
+    cfg.set_eos_gamma(gamma)
+    cfg.set_Csafe(0.8)
+    cfg.set_riemann_solver_hllc()
+    cfg.set_slope_lim_vanleer_sym()
+    cfg.set_face_time_interpolation(True)
+    cfg.set_boundary_condition("x", "reflective")
+    cfg.set_boundary_condition("y", "periodic")
+    cfg.set_boundary_condition("z", "periodic")
+
+    model.set_solver_config(cfg)
+    model.init_scheduler(int(1e7), 1)
+    model.make_base_grid((0, 0, 0), (sz, sz, sz), (base * multx, base * multy, base * multz))
+
+    x0_l = 0.1
+    x0_r = 0.9
+    p0_l = 1000
+    p0_m = 0.01
+    p0_r = 100
+    tmax = 0.02
+
+    def rho_map(rmin, rmax):
+        return 1
+
+    def rhovel_map(rmin, rmax):
+        return (0, 0, 0)
+
+    def rhoetot_map(rmin, rmax):
+        x, y, z = rmin
+        if x < x0_l:
+            return p0_l / (gamma - 1.0)
+        elif x < x0_r:
+            return p0_m / (gamma - 1.0)
+        else:
+            return p0_r / (gamma - 1.0)
+
+    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)
+
+    dt_evolve = tmax / timestamps
+
+    def analysis():
+        results = []
+        for i in range(timestamps + 1):
+            model.evolve_until(dt_evolve * i)
+            rho_vals = model.render_slice("rho", "f64", positions_plot)
+            results.append(rho_vals)
+        return results
+
+    def plot_results(data, tmax):
+        arr_x = [x[0] for x in positions_plot]
+        for i, frame in enumerate(data):
+            fig, axs = plt.subplots(1, 1, figsize=(8, 8))
+            fig.suptitle(f"Interacting blast wave test at t = {dt_evolve * i}", fontsize=14)
+            axs.set_xlabel("$x$")
+            axs.set_ylabel("$\\rho$")
+            axs.grid(True, alpha=0.3)
+            axs.plot(arr_x, frame, label=f"{base * 2}", linewidth=1)
+            axs.legend()
+            plt.tight_layout()
+            plt.savefig(os.path.join(output_folder, f"Interacting_blast_wave_test_{i:04d}.png"))
+            plt.close(fig)
+
+    def gif_results(data, tmax, case_anim="inter-blast"):
+
+        arr_x = [x[0] for x in positions_plot]
+
+        import matplotlib.animation as animation
+
+        fig2, axs = plt.subplots(1, 1, figsize=(8, 8))
+        fig2.suptitle(f"{case_anim} - t = {0.0:.3f} s", fontsize=14)
+
+        # Calculate global min/max across all frames for fixed y-axis limits
+        rho_min = min(np.min(frame) for frame in data)
+        rho_max = max(np.max(frame) for frame in data)
+        # Add 5% margin to y-axis limits
+        rho_margin = (rho_max - rho_min) * 0.05
+        # Configure each axis
+        axs.set_xlabel("$x$")
+        axs.set_ylabel("$\\rho$")
+        axs.set_ylim(rho_min - rho_margin, rho_max + rho_margin)
+        axs.grid(True, alpha=0.3)
+        (line_rho,) = axs.plot(arr_x, data[0], label="$\\rho$", linewidth=2, color="C0")
+        axs.legend()
+
+        def animate(frame):
+            t = tmax * frame / timestamps
+            line_rho.set_ydata(data[frame])
+
+            fig2.suptitle(f"{case_anim} - t = {t:.3f} s", fontsize=14)
+            return (line_rho, 0)
+
+        anim = animation.FuncAnimation(
+            fig2, animate, frames=timestamps + 1, interval=50, blit=False, repeat=True
+        )
+        plt.tight_layout()
+        writer = animation.PillowWriter(fps=15, metadata=dict(artist="Me"), bitrate=1800)
+        anim.save(os.path.join(output_folder, "Interacting_blast_wave_test.gif"), writer=writer)
+        return anim
+
+    def run_and_plot():
+
+        data = analysis()
+
+        return plot_results(data, tmax), gif_results(data, tmax)
+
+    plot, anim = run_and_plot()
+
+
+run_1d_interacting_blast_wave_test(16, 2, 1)

examples/TO_MIGRATE/ramses/liska_wendroff_amr.py

@@ -0,0 +1,192 @@
+import os
+
+import matplotlib
+import matplotlib.animation as animation
+import matplotlib.pyplot as plt
+import numpy as np
+
+import shamrock
+
+if not shamrock.sys.is_initialized():
+    shamrock.change_loglevel(1)
+    shamrock.sys.init("0:0")
+
+
+multx = 1
+multy = 1
+multz = 1
+max_amr_lev = 1
+cell_size = 2 << max_amr_lev  # refinement is limited to cell_size = 2
+base = 64
+scale_fact = 0.3 / (cell_size * base * multx)
+gamma = 1.4
+err_min = 0.30
+err_max = 0.10
+nx = base * 2
+ny = base * 2
+sim_folder = f"_to_trash/ramses_liska_wendroff/check_nan_reso_{nx}_hll/"
+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
+    x_vals = np.linspace(min_x, max_x, nx)
+    y_vals = np.linspace(min_y, max_y, ny)
+
+    # Create meshgrid
+    x_grid, y_grid = np.meshgrid(x_vals, y_vals, indexing="ij")
+
+    # Convert to Cartesian coordinates (z = 0 for a disc in the xy-plane)
+    z_grid = z_val * np.ones_like(x_grid)
+
+    # Flatten and stack to create list of positions
+    positions = np.column_stack([x_grid.ravel(), y_grid.ravel(), z_grid.ravel()])
+
+    return [tuple(pos) for pos in positions]
+
+
+positions = make_cartesian_coords(nx, ny, 0.2, 0, 0.3 - 1e-6, 0, 0.3 - 1e-6)
+
+
+def plot_rho_slice_cartesian(metadata, arr_rho_pos, iplot, case_name, dpi=200):
+    ext = metadata["extent"]
+
+    my_cmap = matplotlib.colormaps["viridis"].copy()  # copy the default cmap
+    my_cmap.set_bad(color="black")
+
+    arr_rho_pos = np.array(arr_rho_pos).reshape(nx, ny)
+    # bad_values_mask = ~np.isfinite(arr_rho_pos)
+    # bad_values_indices = np.argwhere(bad_values_mask)
+    # bad_values = arr_rho_pos[bad_values_mask]
+    # output = np.column_stack((bad_values_indices, bad_values))
+
+    # np.savetxt(f"bad_indices_with_values_{iplot}.txt", output,
+    #         fmt=["%d", "%d", "%.6f"],
+    #         header="row col value")
+
+    # np.savetxt(f"data_{iplot}", arr_rho_pos[:, base -2])
+
+    ampl = 1e-5
+
+    X = np.linspace(ext[0], ext[1], nx)
+    Y = np.linspace(ext[2], ext[3], ny)
+    X, Y = np.meshgrid(X, Y)
+
+    plt.figure(dpi=dpi)
+    vmin = 0
+    vmax = 2.5
+    levels = np.arange(vmin, vmax + 0.025, 0.025)
+
+    # out_of_range = (arr_rho_pos < vmin) | (arr_rho_pos > vmax)
+    # print(f"Out-of-range count for iplot = {iplot}:", np.sum(out_of_range))
+    res = plt.contourf(X, Y, arr_rho_pos, levels=levels, cmap=my_cmap)
+    cs = plt.contour(X, Y, arr_rho_pos, levels=levels, colors="black", linewidths=0.5)
+
+    plt.xlabel("x")
+    plt.ylabel("y")
+    plt.title(f"t = {metadata['time']:0.3f} [seconds]")
+    cbar = plt.colorbar(res, extend="both")
+    cbar.set_label(r"$\rho$ [code unit]")
+    plt.savefig(os.path.join(sim_folder, f"rho_{case_name}_{iplot:04d}.png"))
+    plt.close()
+
+
+from shamrock.utils.plot import show_image_sequence
+
+
+def run_case(set_bc_func, case_name):
+
+    ctx = shamrock.Context()
+    ctx.pdata_layout_new()
+
+    model = shamrock.get_Model_Ramses(context=ctx, vector_type="f64_3", grid_repr="i64_3")
+
+    cfg = model.gen_default_config()
+    cfg.set_scale_factor(scale_fact)
+    set_bc_func(cfg)
+    cfg.set_eos_gamma(gamma)
+    cfg.set_Csafe(0.5)
+    cfg.set_riemann_solver_hllc()
+    cfg.set_slope_lim_vanleer_sym()
+    cfg.set_face_time_interpolation(True)
+    # cfg.set_amr_mode_pseudo_gradient_based(error_min=err_min, error_max=err_max)
+
+    model.set_solver_config(cfg)
+    model.init_scheduler(int(1e7), 1)
+    model.make_base_grid(
+        (0, 0, 0), (cell_size, cell_size, cell_size), (base * multx, base * multy, base * multz)
+    )
+
+    def rho_map(rmin, rmax):
+
+        x, y, z = rmin
+        if (x + y) < 0.15:
+            return 0.125
+        else:
+            return 1.0
+
+    etot_L = 1.0 / (gamma - 1)
+    etot_R = 0.14 / (gamma - 1)
+
+    def rhoetot_map(rmin, rmax):
+
+        rho = rho_map(rmin, rmax)
+
+        x, y, z = rmin
+        if (x + y) < 0.15:
+            return etot_R
+        else:
+            return etot_L
+
+    def rhovel_map(rmin, rmax):
+        return (0, 0, 0)
+
+    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)
+
+    # model.evolve_once(0,0.1)
+    fact = 200
+    tmax = 0.01 * fact
+    all_t = np.linspace(0, tmax, fact)
+
+    def plot(t, iplot):
+        metadata = {"extent": [0, 0.3, 0, 0.3], "time": t}
+        arr_rho_pos = model.render_slice("rho", "f64", positions)
+        plot_rho_slice_cartesian(metadata, arr_rho_pos, iplot, case_name)
+
+    current_time = 0.0
+    for i, t in enumerate(all_t):
+        # model.dump_vtk(os.path.join(sim_folder, f"{case_name}_{i:04d}.vtk"))
+        model.evolve_until(t)
+        current_time = t
+        plot(current_time, i)
+
+    plot(current_time, len(all_t))
+
+    # If the animation is not returned only a static image will be shown in the doc
+    ani = show_image_sequence(os.path.join(sim_folder, f"rho_{case_name}_*.png"), render_gif=True)
+
+    if shamrock.sys.world_rank() == 0:
+        # To save the animation using Pillow as a gif
+        writer = animation.PillowWriter(fps=15, metadata=dict(artist="Me"), bitrate=1800)
+        ani.save(os.path.join(sim_folder, f"rho_{case_name}.gif"), writer=writer)
+
+        return ani
+    else:
+        return None
+
+
+def run_case_reflective():
+    def set_bc_func(cfg):
+        cfg.set_boundary_condition("x", "reflective")
+        cfg.set_boundary_condition("y", "reflective")
+        cfg.set_boundary_condition("z", "reflective")
+
+    return run_case(set_bc_func, "reflective")
+
+
+ani_reflective = run_case_reflective()
+# plt.show()