[SPH] Add L2 error to dustysettle example#1906
Conversation
|
Thanks @tdavidcl 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: |
![gemini-code-assist[bot]](https://avatars.githubusercontent.com/in/956858?s=60&v=4){kind=small}
There was a problem hiding this comment.
Code Review
This pull request implements a periodic HCP box setup using a new pybind11 binding get_periodic_hcp_box, updates the dust settling time scale, and adds L2 error calculation and plotting. The review feedback highlights several critical improvements: first, binning the 3D SPH particle data before computing the L2 error to avoid undefined behavior with np.interp; second, wrapping the new plotting code in an MPI rank check to prevent race conditions; third, ensuring <pybind11/stl.h> is included for proper type casting in the C++ binding; and finally, addressing PEP 8 style violations and potential asymmetric bounds from negative floor division.
Important
The consumer version of Gemini Code Assist on GitHub is being sunset. Starting June 18, 2026, new organization installations will be blocked, and all code review activity will officially cease on July 17, 2026.
For more details on the timeline and next steps, please review the Help Documentation.
| field_field_sorted = field[z_field_sort_args] | ||
|
|
||
| # Interpolate field to the same grid as the reference field | ||
| field_interp = np.interp(z_ref, z_field_sorted, field_field_sorted) |
There was a problem hiding this comment.
Using np.interp directly on raw 3D SPH particle data projected onto the 1D
- SPH particles are distributed in 3D, meaning there will be multiple particles with very close or identical
$z$ coordinates but different field values (scatter).np.interpexpects a single-valued function on a strictly increasing grid. If there are duplicate$z$ values, the behavior ofnp.interpis undefined. - It does not perform any averaging or smoothing, making the L2 error highly sensitive to individual particle noise rather than the overall fluid profile.
To get a meaningful L2 error, you should first bin the particle data along the scipy.stats.binned_statistic or np.histogram) to obtain a clean 1D profile, and then compute the L2 error against the reference solution.
| # %% | ||
| # Plot the L2 errors | ||
|
|
||
| plt.figure() | ||
| for k in range(ndust): | ||
| plt.plot(t_l2, l2_errors[k], label=f"dust {k}") | ||
| print(f"L2 error for dust {k} = {l2_errors[k]}") | ||
| plt.xlabel("t") | ||
| plt.ylabel("L2 error") | ||
| plt.yscale("log") | ||
| plt.title("L2 errors") | ||
| plt.legend() | ||
| plt.tight_layout() | ||
| plt.savefig(f"{dump_folder}/plots/l2_errors.png") | ||
| plt.show() No newline at end of file |
There was a problem hiding this comment.
The new plotting block for L2 errors is executed on all MPI ranks. In a parallel run, this will cause all ranks to simultaneously attempt to plot, call plt.show() (which can block or crash in headless environments), and write to l2_errors.png (causing file corruption or race conditions).
This entire block should be wrapped in an if shamrock.sys.world_rank() == 0: check to ensure only the root process handles plotting and file saving.
| # %% | |
| # Plot the L2 errors | |
| plt.figure() | |
| for k in range(ndust): | |
| plt.plot(t_l2, l2_errors[k], label=f"dust {k}") | |
| print(f"L2 error for dust {k} = {l2_errors[k]}") | |
| plt.xlabel("t") | |
| plt.ylabel("L2 error") | |
| plt.yscale("log") | |
| plt.title("L2 errors") | |
| plt.legend() | |
| plt.tight_layout() | |
| plt.savefig(f"{dump_folder}/plots/l2_errors.png") | |
| plt.show() | |
| # %% | |
| # Plot the L2 errors | |
| if shamrock.sys.world_rank() == 0: | |
| plt.figure() | |
| for k in range(ndust): | |
| plt.plot(t_l2, l2_errors[k], label=f"dust {k}") | |
| print(f"L2 error for dust {k} = {l2_errors[k]}") | |
| plt.xlabel("t") | |
| plt.ylabel("L2 error") | |
| plt.yscale("log") | |
| plt.title("L2 errors") | |
| plt.legend() | |
| plt.tight_layout() | |
| plt.savefig(f"{dump_folder}/plots/l2_errors.png") | |
| plt.show() |
| lx = 12 | ||
| ly = 12 | ||
| lz = 512 | ||
| lmin = (-lx//2,-ly//2,-lz//2) | ||
| lmax = (lx//2,ly//2,lz//2) |
There was a problem hiding this comment.
Adhering to PEP 8 style guidelines, binary operators (like //) and commas should be separated by spaces.
Additionally, using -lx // 2 can lead to asymmetric box bounds if lx is odd (since floor division of negative numbers rounds towards negative infinity, e.g., -11 // 2 = -6 while 11 // 2 = 5). Using -(lx // 2) is more robust and ensures symmetric bounds.
| lx = 12 | |
| ly = 12 | |
| lz = 512 | |
| lmin = (-lx//2,-ly//2,-lz//2) | |
| lmax = (lx//2,ly//2,lz//2) | |
| lx = 12 | |
| ly = 12 | |
| lz = 512 | |
| lmin = (-(lx // 2), -(ly // 2), -(lz // 2)) | |
| lmax = (lx // 2, ly // 2, lz // 2) |
References
- PEP 8: Always surround binary operators with a single space on either side, and use a space after commas. (link)
| def compute_L2_error(z,field, z_ref, field_ref): | ||
|
|
||
| z_field_sort_args = np.argsort(z) | ||
| z_field_sorted = z[z_field_sort_args] | ||
| field_field_sorted = field[z_field_sort_args] | ||
|
|
||
| # Interpolate field to the same grid as the reference field | ||
| field_interp = np.interp(z_ref, z_field_sorted, field_field_sorted) | ||
|
|
||
| # Compute delta | ||
| delta = field_interp - field_ref | ||
|
|
||
| # Compute L2 func | ||
| L2_func = delta ** 2 | ||
|
|
||
| # Compute L2 integral # Use the appropriate trapezoidal integration function | ||
| if hasattr(np, "trapezoid"): | ||
| L2_integral = np.trapezoid(L2_func, z_ref) # NumPy >= 2.0 | ||
| else: | ||
| L2_integral = np.trapz(L2_func, z_ref) # NumPy < 2.0 | ||
|
|
||
|
|
||
| return np.sqrt(L2_integral) |
There was a problem hiding this comment.
Adhering to PEP 8 style guidelines, there should be a space after the comma in the function arguments (z, field), and double blank lines inside the function body should be avoided.
def compute_L2_error(z, field, z_ref, field_ref):
z_field_sort_args = np.argsort(z)
z_field_sorted = z[z_field_sort_args]
field_field_sorted = field[z_field_sort_args]
# Interpolate field to the same grid as the reference field
field_interp = np.interp(z_ref, z_field_sorted, field_field_sorted)
# Compute delta
delta = field_interp - field_ref
# Compute L2 func
L2_func = delta ** 2
# Compute L2 integral # Use the appropriate trapezoidal integration function
if hasattr(np, "trapezoid"):
L2_integral = np.trapezoid(L2_func, z_ref) # NumPy >= 2.0
else:
L2_integral = np.trapz(L2_func, z_ref) # NumPy < 2.0
return np.sqrt(L2_integral)References
- PEP 8: Use a single space after commas in function definitions, and avoid extraneous blank lines inside functions. (link)
| math_module.def( | ||
| "get_periodic_hcp_box", | ||
| [](f64 dr, std::array<i32, 3> box_min, std::array<i32, 3> box_max) { |
There was a problem hiding this comment.
The new binding get_periodic_hcp_box uses std::array<i32, 3> in its lambda signature. For pybind11 to automatically cast Python lists/tuples to std::array, the header <pybind11/stl.h> must be included. Please ensure that #include <pybind11/stl.h> is explicitly included at the top of pyshammath.cpp if it is not already transitively included, to prevent compilation or runtime type-casting errors.
for more information, see https://pre-commit.ci
Workflow reportworkflow report corresponding to commit b3fafad Pre-commit check reportSome failures were detected in base source checks checks. ❌ ruff-checkSuggested changesDetailed changes : |
1 participant
No description provided.