[AMR][RefinementHandler] Main PR with sod and amr tests#1719
Draft
Akos299 wants to merge 55 commits into
Draft
Conversation
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
Summary of ChangesHello, I'm Gemini Code Assist1! I'm currently reviewing this pull request and will post my feedback shortly. In the meantime, here's a summary to help you and other reviewers quickly get up to speed! This pull request significantly enhances the Adaptive Mesh Refinement (AMR) capabilities of the RAMSES solver by introducing a new pseudo-gradient based refinement criterion. It includes comprehensive C++ backend modifications to support this new mode, alongside new Python examples that validate its functionality for common hydrodynamic test cases. The changes also refine the grid management process by enforcing the 2:1 refinement rule and updating internal data structures to better track and utilize AMR levels. Highlights
Using Gemini Code AssistThe full guide for Gemini Code Assist can be found on our documentation page, here are some quick tips. Invoking Gemini You can request assistance from Gemini at any point by creating a comment using either
Customization To customize Gemini Code Assist for GitHub experience, repository maintainers can create a configuration file and/or provide a custom code review style guide (such as PEP-8 for Python) by creating and adding files to a Limitations & Feedback Gemini Code Assist may make mistakes. Please leave feedback on any instances where its feedback is incorrect or counter productive. You can react with 👍 and 👎 on @gemini-code-assist comments. If you're interested in giving your feedback about your experience with Gemini Code Assist for GitHub and other Google products, sign up here. Footnotes
|
![gemini-code-assist[bot]](https://avatars.githubusercontent.com/in/956858?s=60&v=4){kind=small}
Contributor
There was a problem hiding this comment.
Code Review
This pull request introduces a new pseudo-gradient-based AMR refinement mode, including new PseudoGradientBased structs and related methods in SolverConfig.hpp and AMRGridRefinementHandler.hpp, along with corresponding Python bindings. It also adds new utility functions for pseudo-gradient calculation in SlopeLimitedGradientUtilities.hpp and a rho_primitive field to SolverStorage.hpp and Solver.cpp for use in the new AMR logic. The changes also include updates to ConsToPrimGas to handle this new field. Review comments highlight a critical bug where refine_list is moved multiple times, leading to undefined behavior, and several improvement opportunities such as removing unused parameters and redundant computations, clarifying confusing comments, addressing potential logical inconsistencies in err_min/err_max values, and suggesting more robust convergence checks for iteration loops in the AMR refinement process.
| enforce_two_to_one_refinement(std::move(refine_list)); | ||
|
|
||
| /////// enforce 2:1 for derefinement ////// | ||
| enforce_two_to_one_derefinement(std::move(derefine_list), std::move(refine_list)); |
Contributor
There was a problem hiding this comment.
The refine_list is moved to enforce_two_to_one_refinement on line 1246 and then moved again to enforce_two_to_one_derefinement here. Moving an already moved object results in undefined behavior. If the content of refine_list is still needed after the first move, it should be passed by reference or copied.
| using namespace sham; | ||
| using namespace sham::details; | ||
|
|
||
| auto get_avg_neigh = [&](auto &graph_links, u32 dir) -> T { |
Contributor
There was a problem hiding this comment.
The u32 dir parameter in the get_avg_neigh lambda is unused. Please remove it to improve code clarity and avoid potential confusion.
Suggested change
| auto get_avg_neigh = [&](auto &graph_links, u32 dir) -> T { | |
| auto get_avg_neigh = [&](auto &graph_links) -> T { |
|
|
||
|
|
||
| def rhovel_map(rmin, rmax): | ||
| rho = rho_map(rmin, rmax) |
Contributor
There was a problem hiding this comment.
The rho variable is calculated here but is not used within the rhovel_map function. This is a redundant computation and can be removed.
| next_dt = model.evolve_once_override_time(t, dt) | ||
| if i == 0: | ||
| dic0 = convert_to_cell_coords(ctx.collect_data()) | ||
| dX0.append(dic0["ymax"][i] - dic0["ymin"][i]) |
Contributor
There was a problem hiding this comment.
This line appends to dX0 inside a loop, but the condition if i == 0 means it will only execute once, causing dX0 to contain a single element repeated many times. If dX0 is meant to store initial cell sizes, it should be populated only once outside the loop, or its purpose needs clarification.
Comment on lines
+282
to
+293
| auto get_coord = [](u32 i) -> std::array<u32, dim> { | ||
| constexpr u32 NsideBlockPow = 1; | ||
| constexpr u32 Nside = 1U << NsideBlockPow; | ||
| constexpr u32 side_size = Nside; | ||
| constexpr u32 block_size = shambase::pow_constexpr<dim>(Nside); | ||
|
|
||
| if constexpr (dim == 3) { | ||
| const u32 tmp = i >> NsideBlockPow; | ||
| // This line is why derefinement never happens | ||
| // return {i % Nside, (tmp) % Nside, (tmp ) >> NsideBlockPow}; | ||
| return {(tmp) >> NsideBlockPow, (tmp) % Nside, i % Nside}; | ||
| } |
Contributor
There was a problem hiding this comment.
The commented-out line // This line is why derefinement never happens is confusing. If this was a bug that has been fixed, the comment should be removed. If it's still a potential issue or a historical note, please clarify its meaning and current relevance.
| // // // enforce 2:1 at parent level | ||
| // /////////////////////////////////////////////////////////////////////////////////// | ||
|
|
||
| for (int it = 0; it < 3; it++) { |
Contributor
There was a problem hiding this comment.
The for (int it = 0; it < 3; it++) loop uses a fixed number of iterations for convergence. Similar to the refinement loop, consider adding a convergence check or a more dynamic stopping condition for robustness.
| u64 id_patch, | ||
| shamrock::patch::Patch p, | ||
| shamrock::patch::PatchDataLayer &pdat, | ||
| Tscal err_min, |
Contributor
There was a problem hiding this comment.
The p parameter in the RefineCritPseudoGradientAccessor constructor is unused. Please remove it.
Suggested change
| Tscal err_min, | |
| u64 id_patch, |
| u64 id_patch, | ||
| shamrock::patch::Patch p, | ||
| shamrock::patch::PatchDataLayer &pdat, | ||
| Tscal err_min, |
Contributor
There was a problem hiding this comment.
The p parameter in the RefineCritPseudoGradientAccessor::finalize method is unused. Please remove it.
Suggested change
| Tscal err_min, | |
| u64 id_patch, |
Comment on lines
+37
to
+38
| err_min = 0.25 | ||
| err_max = 0.15 |
Contributor
There was a problem hiding this comment.
Similar to sedo_amr.py, err_min (0.25) is greater than err_max (0.15). This configuration is unusual for refinement/derefinement thresholds and could lead to unexpected behavior. Please confirm if this is the intended logic.
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.
Akos299
force-pushed
the
features/all-amr-dev
branch
from
0f59064 to
cc83d20
Compare
Akos299
force-pushed
the
features/all-amr-dev
branch
from
cc83d20 to
9f73d96
Compare
Save the beautifull AMR
tdavidcl
reviewed
Apr 17, 2026
tdavidcl
reviewed
Apr 17, 2026
| internal_refine_grid(shambase::DistributedData<sham::DeviceBuffer<u32>> &&refine_list) { | ||
|
|
||
| using namespace shamrock::patch; | ||
| internal_refine_grid(shambase::DistributedData<sham::DeviceBuffer<u32>> &&dd_refine_flags) { |
Member
There was a problem hiding this comment.
Those renaming could be a sub PR also
Member
|
i convert the PR to draft since it is not intended to be merged in one block and will be used mostly for stagging |
tdavidcl
requested changes
May 2, 2026
|
|
||
| struct Edges { | ||
| const shamrock::solvergraph::PatchDataFieldDDShared<T> &ghost_fields; | ||
| shamrock::solvergraph::PatchDataFieldDDShared<T> &ghost_fields; |
Member
There was a problem hiding this comment.
Suggested change
| shamrock::solvergraph::PatchDataFieldDDShared<T> &ghost_fields; | |
| const shamrock::solvergraph::PatchDataFieldDDShared<T> &ghost_fields; |
Comment on lines
+41
to
+42
| __internal_set_ro_edges({}); | ||
| __internal_set_rw_edges({ghost_fields, fields}); |
Member
There was a problem hiding this comment.
Suggested change
| __internal_set_ro_edges({}); | |
| __internal_set_rw_edges({ghost_fields, fields}); | |
| __internal_set_ro_edges({ghost_fields}); | |
| __internal_set_rw_edges({ fields}); |
tdavidcl
reviewed
May 2, 2026
| fields.get_refs().for_each([&](u32 id_patch, PatchDataField<T> &field) { | ||
| // TODO: currently we guess the GZ size using the input, we should use in place ghost | ||
| // zones really ... | ||
| PatchDataField<T> temp_pdat(field.get_name(), field.get_nvar()); |
Member
There was a problem hiding this comment.
unused
Suggested change
| PatchDataField<T> temp_pdat(field.get_name(), field.get_nvar()); |
…stance between cell centers. Important for the AMR.
Contributor
Workflow reportworkflow report corresponding to commit a7c9015 Light CI is enabled. This will only run the basic tests and not the full tests. Pre-commit check reportSome failures were detected in base source checks checks. ❌ ruff-check❌ forbid-tabs❌ end-of-file-fixer❌ trailing-whitespace❌ ruff-format❌ remove-tabsSuggested changesDetailed changes :diff --git a/examples/TO_MIGRATE/ramses/collapse.py b/examples/TO_MIGRATE/ramses/collapse.py
index 1aa9f7f9..57dc812f 100644
--- a/examples/TO_MIGRATE/ramses/collapse.py
+++ b/examples/TO_MIGRATE/ramses/collapse.py
@@ -58,7 +58,7 @@ min_reso = (L0 * N_J) / (lamb_J)
print(f"min reso = {min_reso}\n")
gamma = 5.0 / 3.0
-rho_c = 3.7e-13*1e3 #[g/cm^3 ==> kg/m^3]
+rho_c = 3.7e-13 * 1e3 # [g/cm^3 ==> kg/m^3]
def run_sim():
@@ -131,13 +131,12 @@ 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.))
+ P = cs_sqr * rho * (1.0 + x ** (2.0 / 3.0))
Eint = P / (gamma - 1.0)
- return Ekin + Eint
-
+ return Ekin + Eint
model.set_field_value_lambda_f64("rho", rho_map)
model.set_field_value_lambda_f64("rhoetot", rhoe_map)
@@ -164,7 +163,6 @@ def run_sim():
dt = tmax - t
if t == tmax:
if shamrock.sys.world_rank() == 0:
-
model.dump_vtk(f"base32_NJ16_Max_8_bar_collapse{i:04d}.vtk")
times.append(t)
break
diff --git a/examples/TO_MIGRATE/ramses/interacting_blast_wave.py b/examples/TO_MIGRATE/ramses/interacting_blast_wave.py
index 4ca5a6a7..ae5796b1 100644
--- a/examples/TO_MIGRATE/ramses/interacting_blast_wave.py
+++ b/examples/TO_MIGRATE/ramses/interacting_blast_wave.py
@@ -61,7 +61,6 @@ def run_1d_interacting_blast_wave_test(
p0_m = 0.01
p0_r = 100
-
def rho_map(rmin, rmax):
return 1
@@ -81,71 +80,76 @@ def run_1d_interacting_blast_wave_test(
model.set_field_value_lambda_f64("rhoetot", rhoetot_map)
model.set_field_value_lambda_f64_3("rhovel", rhovel_map)
-
-
def analysis():
results = []
for dt in timestep_lists:
model.evolve_until(dt)
- rho_vals = model.render_slice("rho", "f64", positions_plot)
+ rho_vals = model.render_slice("rho", "f64", positions_plot)
rhov_vals = model.render_slice("rhovel", "f64_3", positions_plot)
rhoetot_vals = model.render_slice("rhoetot", "f64", positions_plot)
vx = np.array(rhov_vals)[:, 0] / np.array(rho_vals)
P = (np.array(rhoetot_vals) - 0.5 * np.array(rho_vals) * vx**2) * (gamma - 1)
- e_int = P/((gamma - 1.) * np.array(rho_vals))
+ e_int = P / ((gamma - 1.0) * np.array(rho_vals))
results_dic = {
- "rho": np.array(rho_vals),
- "vx": np.array(vx),
- "P": np.array(P),
- "e_int":np.array(e_int)
- }
+ "rho": np.array(rho_vals),
+ "vx": np.array(vx),
+ "P": np.array(P),
+ "e_int": np.array(e_int),
+ }
results.append(results_dic)
-
- output = np.column_stack((np.array(rho_vals), np.array(vx), np.array(P), np.array(e_int)))
- filename= f"data_interacting_blast_wave_reso_{2*base*multx}_at_{dt}.txt"
- np.savetxt(os.path.join(output_folder,filename),
- output,
- fmt=["%.10f", "%.10f", "%.10f", "%.10f"],
- header="rho vx P e_int",
- )
+
+ output = np.column_stack(
+ (np.array(rho_vals), np.array(vx), np.array(P), np.array(e_int))
+ )
+ filename = f"data_interacting_blast_wave_reso_{2 * base * multx}_at_{dt}.txt"
+ np.savetxt(
+ os.path.join(output_folder, filename),
+ output,
+ fmt=["%.10f", "%.10f", "%.10f", "%.10f"],
+ header="rho vx P e_int",
+ )
return results
def plot_results(data):
arr_x = [x[0] for x in positions_plot]
for i, frame in enumerate(data):
fig, axs = plt.subplots(2, 2, figsize=(8, 8), dpi=100, constrained_layout=True)
-
- ## density
- axs[0,0].set_xlabel("$x$")
- axs[0,0].set_ylabel("$\\rho$")
- axs[0,0].grid(True, alpha=0.3)
- axs[0,0].scatter(arr_x, frame["rho"], label=f"{2*base*multx}", s=10, marker = "*")
- axs[0,0].legend(loc=0)
+
+ ## density
+ axs[0, 0].set_xlabel("$x$")
+ axs[0, 0].set_ylabel("$\\rho$")
+ axs[0, 0].grid(True, alpha=0.3)
+ axs[0, 0].scatter(arr_x, frame["rho"], label=f"{2 * base * multx}", s=10, marker="*")
+ axs[0, 0].legend(loc=0)
## velocity
- axs[0,1].set_xlabel("$x$")
- axs[0,1].set_ylabel("$v_\\mathrm{x}$")
- axs[0,1].grid(True, alpha=0.3)
- axs[0,1].scatter(arr_x, frame["vx"], label=f"{2*base*multx}", s=10, marker ="*" )
- axs[0,1].legend(loc=0)
-
-
+ axs[0, 1].set_xlabel("$x$")
+ axs[0, 1].set_ylabel("$v_\\mathrm{x}$")
+ axs[0, 1].grid(True, alpha=0.3)
+ axs[0, 1].scatter(arr_x, frame["vx"], label=f"{2 * base * multx}", s=10, marker="*")
+ axs[0, 1].legend(loc=0)
+
## pressure
- axs[1,1].set_xlabel("$x$")
- axs[1,1].set_ylabel("$\\mathrm{P}$")
- axs[1,1].grid(True, alpha=0.3)
- axs[1,1].scatter(arr_x, frame["P"], label=f"{2*base*multx}", s=10, marker ="*" )
- axs[1,1].legend(loc=0)
+ axs[1, 1].set_xlabel("$x$")
+ axs[1, 1].set_ylabel("$\\mathrm{P}$")
+ axs[1, 1].grid(True, alpha=0.3)
+ axs[1, 1].scatter(arr_x, frame["P"], label=f"{2 * base * multx}", s=10, marker="*")
+ axs[1, 1].legend(loc=0)
## internal energy
- axs[1,0].set_xlabel("$x$")
- axs[1,0].set_ylabel("$\\mathrm{e}_\\mathrm{int}$")
- axs[1,0].grid(True, alpha=0.3)
- axs[1,0].scatter(arr_x, frame["e_int"], label=f"{2*base*multx}", s=10, marker ="*" )
- axs[1,0].legend(loc=0)
-
+ axs[1, 0].set_xlabel("$x$")
+ axs[1, 0].set_ylabel("$\\mathrm{e}_\\mathrm{int}$")
+ axs[1, 0].grid(True, alpha=0.3)
+ axs[1, 0].scatter(arr_x, frame["e_int"], label=f"{2 * base * multx}", s=10, marker="*")
+ axs[1, 0].legend(loc=0)
+
plt.tight_layout()
- plt.savefig(os.path.join(output_folder, f"Interacting_blast_wave_test_at_{timestep_lists[i]}_nx_reso_{base*2*multx}.png"))
+ plt.savefig(
+ os.path.join(
+ output_folder,
+ f"Interacting_blast_wave_test_at_{timestep_lists[i]}_nx_reso_{base * 2 * multx}.png",
+ )
+ )
plt.close(fig)
# def gif_results(data, tmax, case_anim="inter-blast"):
@@ -185,7 +189,6 @@ def run_1d_interacting_blast_wave_test(
# axs[0,1].set_ylim(vx_min - vx_margin, vx_max + vx_margin)
# axs[0,1].grid(True, alpha=0.3)
-
# (line_rho,) = axs[0,0].plot(arr_x, data[0], label="$\\rho$", linewidth=2, color="C0")
# axs[0,0].legend()
@@ -213,4 +216,4 @@ def run_1d_interacting_blast_wave_test(
plot = run_and_plot()
-run_1d_interacting_blast_wave_test(8, 2<<2, 1)
+run_1d_interacting_blast_wave_test(8, 2 << 2, 1)
diff --git a/examples/TO_MIGRATE/ramses/shu_osher.py b/examples/TO_MIGRATE/ramses/shu_osher.py
index 2fbed524..d81d23f5 100644
--- a/examples/TO_MIGRATE/ramses/shu_osher.py
+++ b/examples/TO_MIGRATE/ramses/shu_osher.py
@@ -1,5 +1,3 @@
-
-
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
import numpy as np
@@ -12,9 +10,6 @@ ctx.pdata_layout_new()
model = shamrock.get_Model_Ramses(context=ctx, vector_type="f64_3", grid_repr="i64_3")
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,
@@ -24,7 +19,6 @@ if not shamrock.sys.is_initialized():
shamrock.sys.init("0:0")
-
def run_1d_Shu_Osher_test(
nblocks,
nxbox,
@@ -71,37 +65,33 @@ def run_1d_Shu_Osher_test(
model.init_scheduler(int(1e7), 1)
model.make_base_grid((0, 0, 0), (sz, sz, sz), (base * multx, base * multy, base * multz))
-
x_s = 0.2
A0 = 0.2
- fp = 5.
-
-
+ fp = 5.0
def rho_map(rmin, rmax):
- x,y,z = rmin
- rho = rho=1 + A0 * np.sin(fp * np.pi* x)
+ x, y, z = rmin
+ rho = rho = 1 + A0 * np.sin(fp * np.pi * x)
if x < x_s:
- rho=3.857143
+ rho = 3.857143
return rho
-
def rhovel_map(rmin, rmax):
- rho = rho_map(rmin,rmax)
- x,y,z = rmin
- rhovec=(0,0,0)
+ rho = rho_map(rmin, rmax)
+ x, y, z = rmin
+ rhovec = (0, 0, 0)
if x < x_s:
rhovec = (rho * 2.629369, 0, 0)
return rhovec
def rhoetot_map(rmin, rmax):
- rho = rho_map(rmin,rmax)
+ rho = rho_map(rmin, rmax)
rhov = rhovel_map(rmin, rmax)
- rhoe_kin = 0.5 * (rhov[0]**2)/rho
- rhoe_int = 1.0/(gamma - 1.0)
+ rhoe_kin = 0.5 * (rhov[0] ** 2) / rho
+ rhoe_int = 1.0 / (gamma - 1.0)
x, y, z = rmin
if x < x_s:
- rhoe_int = 10.33333/(gamma - 1.0)
+ rhoe_int = 10.33333 / (gamma - 1.0)
return rhoe_int + rhoe_kin
model.set_field_value_lambda_f64("rho", rho_map)
@@ -109,73 +99,80 @@ def run_1d_Shu_Osher_test(
model.set_field_value_lambda_f64_3("rhovel", rhovel_map)
def analysis():
-
+
results = []
# for i in range(timestamps+1):
- for dt in timestep_list:
+ for dt in timestep_list:
# model.evolve_until(dt_evolve * i)
model.evolve_until(dt)
- rho_vals = model.render_slice("rho", "f64", positions_plot)
+ rho_vals = model.render_slice("rho", "f64", positions_plot)
rhov_vals = model.render_slice("rhovel", "f64_3", positions_plot)
rhoetot_vals = model.render_slice("rhoetot", "f64", positions_plot)
vx = np.array(rhov_vals)[:, 0] / np.array(rho_vals)
P = (np.array(rhoetot_vals) - 0.5 * np.array(rho_vals) * vx**2) * (gamma - 1)
- e_int = P/((gamma - 1.) * np.array(rho_vals))
+ e_int = P / ((gamma - 1.0) * np.array(rho_vals))
results_dic = {
- "rho": np.array(rho_vals),
- "vx": np.array(vx),
- "P": np.array(P),
- "e_int":np.array(e_int)
- }
+ "rho": np.array(rho_vals),
+ "vx": np.array(vx),
+ "P": np.array(P),
+ "e_int": np.array(e_int),
+ }
results.append(results_dic)
-
- output = np.column_stack((np.array(rho_vals), np.array(vx), np.array(P), np.array(e_int)))
- filename= f"data_Shu_Osher_reso_{2*base*multx}_at_{dt}.txt"
- np.savetxt(os.path.join(output_folder,filename),
- output,
- fmt=["%.10f", "%.10f", "%.10f", "%.10f"],
- header="rho vx P e_int",
- )
- return results
+ output = np.column_stack(
+ (np.array(rho_vals), np.array(vx), np.array(P), np.array(e_int))
+ )
+ filename = f"data_Shu_Osher_reso_{2 * base * multx}_at_{dt}.txt"
+ np.savetxt(
+ os.path.join(output_folder, filename),
+ output,
+ fmt=["%.10f", "%.10f", "%.10f", "%.10f"],
+ header="rho vx P e_int",
+ )
+ return results
def plot_results(data):
arr_x = [x[0] for x in positions_plot]
for i, frame in enumerate(data):
fig, axs = plt.subplots(2, 2, figsize=(8, 8), dpi=100, constrained_layout=True)
-
- ## density
- axs[0,0].set_xlabel("$x$")
- axs[0,0].set_ylabel("$\\rho$")
- axs[0,0].grid(True, alpha=0.3)
- axs[0,0].scatter(arr_x, frame["rho"], label=f"{2*base*multx}", s=10, marker = "*")
- axs[0,0].legend(loc=0)
+
+ ## density
+ axs[0, 0].set_xlabel("$x$")
+ axs[0, 0].set_ylabel("$\\rho$")
+ axs[0, 0].grid(True, alpha=0.3)
+ axs[0, 0].scatter(arr_x, frame["rho"], label=f"{2 * base * multx}", s=10, marker="*")
+ axs[0, 0].legend(loc=0)
## velocity
- axs[0,1].set_xlabel("$x$")
- axs[0,1].set_ylabel("$v_\\mathrm{x}$")
- axs[0,1].grid(True, alpha=0.3)
- axs[0,1].scatter(arr_x, frame["vx"], label=f"{2*base*multx}", s=10, marker ="*" )
- axs[0,1].legend(loc=0)
-
-
+ axs[0, 1].set_xlabel("$x$")
+ axs[0, 1].set_ylabel("$v_\\mathrm{x}$")
+ axs[0, 1].grid(True, alpha=0.3)
+ axs[0, 1].scatter(arr_x, frame["vx"], label=f"{2 * base * multx}", s=10, marker="*")
+ axs[0, 1].legend(loc=0)
+
## pressure
- axs[1,1].set_xlabel("$x$")
- axs[1,1].set_ylabel("$\\mathrm{P}$")
- axs[1,1].grid(True, alpha=0.3)
- axs[1,1].scatter(arr_x, frame["P"], label=f"{2*base*multx}", s=10, marker ="*" )
- axs[1,1].legend(loc=0)
+ axs[1, 1].set_xlabel("$x$")
+ axs[1, 1].set_ylabel("$\\mathrm{P}$")
+ axs[1, 1].grid(True, alpha=0.3)
+ axs[1, 1].scatter(arr_x, frame["P"], label=f"{2 * base * multx}", s=10, marker="*")
+ axs[1, 1].legend(loc=0)
## internal energy
- axs[1,0].set_xlabel("$x$")
- axs[1,0].set_ylabel("$\\mathrm{e}_\\mathrm{int}$")
- axs[1,0].grid(True, alpha=0.3)
- axs[1,0].scatter(arr_x, frame["e_int"], label=f"{2*base*multx}", s=10, marker ="*" )
- axs[1,0].legend(loc=0)
-
+ axs[1, 0].set_xlabel("$x$")
+ axs[1, 0].set_ylabel("$\\mathrm{e}_\\mathrm{int}$")
+ axs[1, 0].grid(True, alpha=0.3)
+ axs[1, 0].scatter(arr_x, frame["e_int"], label=f"{2 * base * multx}", s=10, marker="*")
+ axs[1, 0].legend(loc=0)
+
plt.tight_layout()
- plt.savefig(os.path.join(output_folder, f"Shu_Osher_test_at_{timestep_list[i]}_nx_reso_{base*2*multx}.png"))
+ plt.savefig(
+ os.path.join(
+ output_folder,
+ f"Shu_Osher_test_at_{timestep_list[i]}_nx_reso_{base * 2 * multx}.png",
+ )
+ )
plt.close(fig)
+
def run_and_plot():
data = analysis()
@@ -185,4 +182,4 @@ def run_1d_Shu_Osher_test(
plot = run_and_plot()
-run_1d_Shu_Osher_test(8, 2<<5, 1)
\ No newline at end of file
+run_1d_Shu_Osher_test(8, 2 << 5, 1)
diff --git a/examples/TO_MIGRATE/ramses/sod_amr_thesis.py b/examples/TO_MIGRATE/ramses/sod_amr_thesis.py
index 80c0d79b..b9850f2a 100644
--- a/examples/TO_MIGRATE/ramses/sod_amr_thesis.py
+++ b/examples/TO_MIGRATE/ramses/sod_amr_thesis.py
@@ -11,6 +11,7 @@ ctx.pdata_layout_new()
model = shamrock.get_Model_Ramses(context=ctx, vector_type="f64_3", grid_repr="i64_3")
+
def run_case(nb_blocks, amr_lev, case="amr"):
multx = 1
multy = 1
@@ -31,11 +32,10 @@ def run_case(nb_blocks, amr_lev, case="amr"):
cfg.set_boundary_condition("z", "reflective")
cfg.set_riemann_solver_hllc()
-
cfg.set_slope_lim_minmod()
cfg.set_face_time_interpolation(True)
- if(case == "amr"):
+ if case == "amr":
err_min = 0.25
err_max = 0.15
cfg.set_amr_mode_pseudo_gradient_based(error_min=err_min, error_max=err_max)
@@ -47,7 +47,6 @@ def run_case(nb_blocks, amr_lev, case="amr"):
(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 < 0.5:
@@ -55,11 +54,9 @@ def run_case(nb_blocks, amr_lev, case="amr"):
else:
return 0.125
-
etot_L = 1.0 / (gamma - 1)
etot_R = 0.1 / (gamma - 1)
-
def rhoetot_map(rmin, rmax):
x, y, z = rmin
if x < 0.5:
@@ -67,16 +64,13 @@ def run_case(nb_blocks, amr_lev, case="amr"):
else:
return etot_R
-
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)
-
def convert_to_cell_coords(dic):
cmin = dic["cell_min"]
@@ -117,7 +111,6 @@ def run_case(nb_blocks, amr_lev, case="amr"):
return dic
-
t_target = 0.245
dt = 0
@@ -198,20 +191,22 @@ def run_case(nb_blocks, amr_lev, case="amr"):
## density
ax00_1 = axs[0, 0]
- if(case == "amr"):
+ if case == "amr":
ax00_2 = ax00_1.twinx()
ax00_2.set_ylim(np.floor(l.min()), np.ceil(l.max())) # optional but important
ax00_2.yaxis.set_major_locator(ticker.MultipleLocator(1))
ax00_2.yaxis.set_major_formatter(ticker.FormatStrFormatter("%d"))
- ax00_1.scatter(X, rho, rasterized=True, s=12 * np.ones(X.shape), color="red", label="numeric")
+ ax00_1.scatter(
+ X, rho, rasterized=True, s=12 * np.ones(X.shape), color="red", label="numeric"
+ )
ax00_1.plot(arr_x, arr_rho, ls="--", lw=2.0, color="black", label="exact")
ax00_1.set_xlabel("$X$")
ax00_1.set_ylabel("density")
ax00_1.legend(loc=0)
ax00_1.grid()
- if(case == "amr"):
+ if case == "amr":
idx = np.argsort(X)
ax00_2.plot(X[idx], l[idx], color="purple", marker="*", linewidth=1.0, ls="-.")
ax00_2.set_ylabel("AMR level")
@@ -219,7 +214,7 @@ def run_case(nb_blocks, amr_lev, case="amr"):
## velocity
ax01_1 = axs[0, 1]
- if(case=="amr"):
+ if case == "amr":
ax01_2 = ax01_1.twinx()
ax01_2.set_ylim(np.floor(l.min()), np.ceil(l.max())) # optional but important
ax01_2.yaxis.set_major_locator(ticker.MultipleLocator(1))
@@ -232,7 +227,7 @@ def run_case(nb_blocks, amr_lev, case="amr"):
ax01_1.legend(loc=0)
ax01_1.grid()
- if(case=="amr"):
+ if case == "amr":
idx = np.argsort(X)
ax01_2.plot(X[idx], l[idx], color="purple", marker="*", linewidth=1.0, ls="-.")
ax01_2.set_ylabel("AMR level")
@@ -240,7 +235,7 @@ def run_case(nb_blocks, amr_lev, case="amr"):
## pressure
ax11_1 = axs[1, 1]
- if(case =="amr"):
+ if case == "amr":
ax11_2 = ax11_1.twinx()
ax11_2.set_ylim(np.floor(l.min()), np.ceil(l.max())) # optional but important
ax11_2.yaxis.set_major_locator(ticker.MultipleLocator(1))
@@ -254,8 +249,8 @@ def run_case(nb_blocks, amr_lev, case="amr"):
ax11_1.set_ylabel("pressure")
ax11_1.legend(loc=0)
ax11_1.grid()
-
- if(case == "amr"):
+
+ if case == "amr":
idx = np.argsort(X)
ax11_2.plot(X[idx], l[idx], color="purple", marker="*", linewidth=1.0, ls="-.")
ax11_2.set_ylabel("AMR level")
@@ -263,14 +258,19 @@ def run_case(nb_blocks, amr_lev, case="amr"):
## internal energy
ax10_1 = axs[1, 0]
- if(case=="amr"):
+ if case == "amr":
ax10_2 = ax10_1.twinx()
ax10_2.set_ylim(np.floor(l.min()), np.ceil(l.max())) # optional but important
ax10_2.yaxis.set_major_locator(ticker.MultipleLocator(1))
ax10_2.yaxis.set_major_formatter(ticker.FormatStrFormatter("%d"))
ax10_1.scatter(
- X, internal_energy, rasterized=True, s=12 * np.ones(X.shape), color="red", label="numeric"
+ X,
+ internal_energy,
+ rasterized=True,
+ s=12 * np.ones(X.shape),
+ color="red",
+ label="numeric",
)
ax10_1.plot(
arr_x, arr_P / (arr_rho * (gamma - 1.0)), ls="--", lw=2.0, color="black", label="exact"
@@ -280,7 +280,7 @@ def run_case(nb_blocks, amr_lev, case="amr"):
ax10_1.legend(loc=0)
ax10_1.grid()
- if(case=="amr"):
+ if case == "amr":
idx = np.argsort(X)
ax10_2.plot(X[idx], l[idx], color="purple", marker="*", linewidth=1.0, ls="-.")
ax10_2.set_ylabel("AMR level")
@@ -289,4 +289,4 @@ def run_case(nb_blocks, amr_lev, case="amr"):
plt.savefig(f"sod-{case}-new-bb.pdf")
-run_case(16, 3,"amr")
+run_case(16, 3, "amr")
diff --git a/src/shammath/include/shammath/riemann.hpp b/src/shammath/include/shammath/riemann.hpp
index 7d59f633..38e30318 100644
--- a/src/shammath/include/shammath/riemann.hpp
+++ b/src/shammath/include/shammath/riemann.hpp
@@ -146,15 +146,15 @@ namespace shammath {
/** This just for testing purpose. Will be remove*/
- auto m_H = 1.67262192e-27; //[kg]
- auto kb = 1.380649e-23;
- auto mu = 2.3; // molecular gas
- auto T = 10.;
- auto rho_c = 3.7e-13 * 1e3; // [g/cm^3 ===> kg/m^3]
+ auto m_H = 1.67262192e-27; //[kg]
+ auto kb = 1.380649e-23;
+ auto mu = 2.3; // molecular gas
+ auto T = 10.;
+ auto rho_c = 3.7e-13 * 1e3; // [g/cm^3 ===> kg/m^3]
auto cs0_sqr = (kb * T) / (mu * m_H);
- prim.press = cons.rho * cs0_sqr * (1.0 + sycl::pow(cons.rho / rho_c, 2./3.));
+ prim.press = cons.rho * cs0_sqr * (1.0 + sycl::pow(cons.rho / rho_c, 2. / 3.));
return prim;
}
@@ -180,19 +180,17 @@ namespace shammath {
template<class Tvec>
inline constexpr shambase::VecComponent<Tvec> sound_speed(
PrimState<Tvec> prim, shambase::VecComponent<Tvec> gamma) {
- //return sycl::sqrt(gamma * prim.press / prim.rho);
- //
- auto rho_c = 3.7e-13 * 1e3; // [g/cm^3 ===> kg/m^3]
+ // return sycl::sqrt(gamma * prim.press / prim.rho);
+ //
+ auto rho_c = 3.7e-13 * 1e3; // [g/cm^3 ===> kg/m^3]
auto m_H = 1.67262192e-27; // [kg]
- auto kb = 1.380649e-23; // []
- auto mu = 2.3; // molecular gas
-
-
- auto T = 10.;
- auto cs0_sqr = (kb * T) / (mu * m_H);
- auto cs_sqr = cs0_sqr * (1. + (5.0/3.0) * sycl::pow(prim.rho/rho_c,2./3.));
+ auto kb = 1.380649e-23; // []
+ auto mu = 2.3; // molecular gas
+ auto T = 10.;
+ auto cs0_sqr = (kb * T) / (mu * m_H);
+ auto cs_sqr = cs0_sqr * (1. + (5.0 / 3.0) * sycl::pow(prim.rho / rho_c, 2. / 3.));
// return sycl::sqrt(gamma * prim.press / prim.rho);
return sycl::sqrt(cs_sqr);
@@ -408,28 +406,27 @@ namespace shammath {
const auto csL = sound_speed(primL, gamma);
const auto csR = sound_speed(primR, gamma);
-
-/*
- if (sycl::isnan(csL) || sycl::isnan(csR)) {
- logger::raw_ln(
- "Nan in HLLC solver \t csL \t",
- csL,
- "\t pL \t",
- primL.press,
- "\t rhoL \t",
- primL.rho,
- "\t csR \t",
- csR,
- "\t pR \t",
- primR.press,
- "\t rhoR \t",
- primR.rho,
- "\t gamma \t",
- gamma,
- "\t\n");
- }
-
-*/
+ /*
+ if (sycl::isnan(csL) || sycl::isnan(csR)) {
+ logger::raw_ln(
+ "Nan in HLLC solver \t csL \t",
+ csL,
+ "\t pL \t",
+ primL.press,
+ "\t rhoL \t",
+ primL.rho,
+ "\t csR \t",
+ csR,
+ "\t pR \t",
+ primR.press,
+ "\t rhoR \t",
+ primR.rho,
+ "\t gamma \t",
+ gamma,
+ "\t\n");
+ }
+
+ */
// Left and right state fluxes
const auto FL = hydro_flux_x(cL, gamma);
diff --git a/src/shammodels/ramses/include/shammodels/ramses/modules/ComputeFluxUtilities.hpp b/src/shammodels/ramses/include/shammodels/ramses/modules/ComputeFluxUtilities.hpp
index c87d0c85..179307e9 100644
--- a/src/shammodels/ramses/include/shammodels/ramses/modules/ComputeFluxUtilities.hpp
+++ b/src/shammodels/ramses/include/shammodels/ramses/modules/ComputeFluxUtilities.hpp
@@ -201,15 +201,13 @@ namespace shammodels::basegodunov::modules {
return "";
};
- /*
- std::string cur_direction = get_dir_name();
- std::string kernel_name = (std::string) "compute " + flux_name + cur_direction;
- const char *_kernel_name = kernel_name.c_str();
- */
-
- constexpr const char* _kernel_name = "compute_flux";
-
+ /*
+ std::string cur_direction = get_dir_name();
+ std::string kernel_name = (std::string) "compute " + flux_name + cur_direction;
+ const char *_kernel_name = kernel_name.c_str();
+ */
+ constexpr const char *_kernel_name = "compute_flux";
sham::EventList depends_list;
auto rho = rho_face_dir.get_read_access(depends_list);
@@ -259,12 +257,9 @@ namespace shammodels::basegodunov::modules {
using d_Flux = DustFluxCompute<Tvec, mode, dir>;
-
-
std::string flux_name
= (mode == DustRiemannSolverMode::DHLL) ? "dust hll flux " : "dust huang-bai flux ";
-
auto get_dir_name = [&]() {
if constexpr (dir == Direction::xp) {
return "xp";
@@ -284,14 +279,14 @@ namespace shammodels::basegodunov::modules {
return "";
};
- /*
- std::string cur_direction = get_dir_name();
- std::string kernel_name = (std::string) "compute " + flux_name + cur_direction;
- const char *_kernel_name = kernel_name.c_str();
+ /*
+ std::string cur_direction = get_dir_name();
+ std::string kernel_name = (std::string) "compute " + flux_name + cur_direction;
+ const char *_kernel_name = kernel_name.c_str();
- */
+ */
- constexpr const char* _kernel_name = "dust_compute_flux";
+ constexpr const char *_kernel_name = "dust_compute_flux";
sham::EventList depends_list;
diff --git a/src/shammodels/ramses/src/modules/AMRGridRefinementHandler.cpp b/src/shammodels/ramses/src/modules/AMRGridRefinementHandler.cpp
index 614eaa83..52502b9b 100644
--- a/src/shammodels/ramses/src/modules/AMRGridRefinementHandler.cpp
+++ b/src/shammodels/ramses/src/modules/AMRGridRefinementHandler.cpp
@@ -1319,7 +1319,7 @@ void shammodels::basegodunov::modules::AMRGridRefinementHandler<Tvec, TgridVec>:
auto mu = 2.3; // molecular gas
auto gamma = 5. / 3.;
- auto rho_c = 3.7e-13 * 1e3; // [g/cm^3 ===> kg/m^3]
+ auto rho_c = 3.7e-13 * 1e3; // [g/cm^3 ===> kg/m^3]
//--------------------------------
Tscal block_max_jeans_length = shambase::VectorProperties<Tscal>::get_zero();
@@ -1330,20 +1330,20 @@ void shammodels::basegodunov::modules::AMRGridRefinementHandler<Tvec, TgridVec>:
block_rho_max = sham::details::g_sycl_max(block_rho_max, block_rho[idx]);
}
- //Tscal block_min_jeans_length = sycl::sqrt(
- // (shamunits::pi<Tscal> * kb * T_0)
- // / (N_J * N_J * ctes.G() * block_rho_max * mu * m_H));
+ // Tscal block_min_jeans_length = sycl::sqrt(
+ // (shamunits::pi<Tscal> * kb * T_0)
+ // / (N_J * N_J * ctes.G() * block_rho_max * mu * m_H));
// auto cs = T_0;
// auto G = 1.;
// Tscal block_min_jeansength = sycl::sqrt(shamunits::pi<Tscal> * cs*cs / ( N_J * N_J
// * G * block_rho_max));
- //
- auto cs0_sqr = (kb * T_0) / (mu * m_H);
- auto cs_sqr = cs0_sqr * (1. + (5.0/3.0) * sycl::pow(block_rho_max/rho_c,2./3.));
+ //
+ auto cs0_sqr = (kb * T_0) / (mu * m_H);
+ auto cs_sqr = cs0_sqr * (1. + (5.0 / 3.0) * sycl::pow(block_rho_max / rho_c, 2. / 3.));
-
- Tscal block_min_jeans_length = sycl::sqrt(shamunits::pi<Tscal> * cs_sqr / ( N_J * N_J * ctes.G() * block_rho_max));
+ Tscal block_min_jeans_length = sycl::sqrt(
+ shamunits::pi<Tscal> * cs_sqr / (N_J * N_J * ctes.G() * block_rho_max));
should_refine = false;
should_derefine = false;
diff --git a/src/shammodels/ramses/src/modules/BICGSTABInit.cpp b/src/shammodels/ramses/src/modules/BICGSTABInit.cpp
index ba7eea9f..5a8d65c1 100644
--- a/src/shammodels/ramses/src/modules/BICGSTABInit.cpp
+++ b/src/shammodels/ramses/src/modules/BICGSTABInit.cpp
@@ -113,23 +113,21 @@ namespace {
: rho[cell_global_id] - mean_rho)
* dV;
- /*
- *
- if (sycl::isnan(rho[cell_global_id]) || sycl::isnan(phi[cell_global_id])) {
- logger::raw_ln(
- "nan in INIT @ \t",
- cell_global_id,
- "\t rho = \t ",
- rho[cell_global_id],
- "\t phi = \t",
- phi[cell_global_id],
- "\t \t");
- }
-
-
- */
-
-
+ /*
+ *
+ if (sycl::isnan(rho[cell_global_id]) || sycl::isnan(phi[cell_global_id])) {
+ logger::raw_ln(
+ "nan in INIT @ \t",
+ cell_global_id,
+ "\t rho = \t ",
+ rho[cell_global_id],
+ "\t phi = \t",
+ phi[cell_global_id],
+ "\t \t");
+ }
+
+
+ */
phi_res[cell_global_id] = b_rhs - Aphi;
phi_res_bis[cell_global_id] = 0.5 * (b_rhs - Aphi);
diff --git a/src/shammodels/ramses/src/modules/ComputeCFL.cpp b/src/shammodels/ramses/src/modules/ComputeCFL.cpp
index 7a6e36a9..77aa425e 100644
--- a/src/shammodels/ramses/src/modules/ComputeCFL.cpp
+++ b/src/shammodels/ramses/src/modules/ComputeCFL.cpp
@@ -97,18 +97,18 @@ auto shammodels::basegodunov::modules::ComputeCFL<Tvec, TgridVec>::compute_cfl()
constexpr Tscal div = 1. / 3.;
- Tscal cs = sound_speed(prim_state, gamma);
+ Tscal cs = sound_speed(prim_state, gamma);
/** Will be remove later. Only for testing the the spherical collapse */
- /*
- auto m_H = 1.67262192e-27; //[kg]
- auto kb = 1.380649e-23;
- auto mu = 2.3; // molecular gas
- auto T = 10;
- auto cs0_sqr = (kb * T) / (mu * m_H);
- Tscal cs = cs0_sqr;
-
- */
+ /*
+ auto m_H = 1.67262192e-27; //[kg]
+ auto kb = 1.380649e-23;
+ auto mu = 2.3; // molecular gas
+ auto T = 10;
+ auto cs0_sqr = (kb * T) / (mu * m_H);
+ Tscal cs = cs0_sqr;
+
+ */
/** ------------------------ **/
Tscal vnorm = sycl::length(prim_state.vel);
Tscal dt = C_safe * dx * div / (cs + vnorm);
diff --git a/src/shammodels/ramses/src/modules/ConsToPrimGas.cpp b/src/shammodels/ramses/src/modules/ConsToPrimGas.cpp
index 0b359d09..a4bda8e2 100644
--- a/src/shammodels/ramses/src/modules/ConsToPrimGas.cpp
+++ b/src/shammodels/ramses/src/modules/ConsToPrimGas.cpp
@@ -74,22 +74,22 @@ namespace {
// P[i] = (prim_state.press > 0.0) ? prim_state.press : 1e-6;
// rho_field[i] = prim_state.rho;
-
... 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.
This PR is the main from which modular ones will be extract from.
sedov_amr_clean_v1.mp4