[SPH][Doc] add Bonnor-Ebert sphere collapse setup#1472
[SPH][Doc] add Bonnor-Ebert sphere collapse setup#1472tdavidcl wants to merge 2 commits intoShamrock-code:mainfrom
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: |
Summary of ChangesHello @tdavidcl, 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 adds a new example script to the documentation, showcasing a complete workflow for simulating the collapse of a Bonnor-Ebert sphere using the 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. You can also get AI-powered code generation, chat, as well as code reviews directly in the IDE at no cost with the Gemini Code Assist IDE Extension. Footnotes
|
Workflow reportworkflow report corresponding to commit 669e0a0 Pre-commit check reportSome failures were detected in base source checks checks. ❌ ruff-checkSuggested changesDetailed changes : |
![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 introduces an example script for simulating a Bonnor-Ebert sphere collapse. The script is a valuable addition to the documentation. My review provides feedback to enhance the script's performance, readability, and maintainability. Key suggestions include vectorizing loops using NumPy/SciPy for better performance, refactoring duplicated code into helper functions, improving the structure of the main simulation loop for clarity, and addressing potential issues that could block non-interactive execution.
|
|
||
| plt.grid() | ||
|
|
||
| plt.show() |
There was a problem hiding this comment.
The plt.show() call will block script execution until the plot window is manually closed. For a documentation example script that might be run non-interactively, this will cause it to hang. Consider saving the figure to a file instead, or removing this call if the plot is not essential for the script's logic that follows. The later part of the script saves figures, which is a better practice for non-interactive runs.
| plt.show() | |
| # plt.show() |
| plt.title("t = {:0.3f} [year]".format(model.get_time())) | ||
|
|
||
| plt.colorbar() | ||
| plt.show() |
There was a problem hiding this comment.
The plt.show() call will block script execution until the plot window is manually closed. For a documentation example script that might be run non-interactively, this will cause it to hang. Consider saving the figure to a file instead, or removing this call if the plot is not essential for the script's logic that follows.
| plt.show() | |
| # plt.show() |
| rho_c_rho_0 = [] | ||
| dimless_arr = [] | ||
| for i in range(len(t)): | ||
|
|
||
| dimless_mass = xidudxi[i] * (4 * np.pi * dens[i] ** -1) ** (-1 / 2) | ||
| # print("rho_c/rho_0",dens[i]**-1) | ||
| # print("xidudxi[-1]",xidudxi[i]) | ||
| # print("(4*np.pi*dens[-1]**-1)**(1/2)",(4*np.pi*dens[i]**-1)**(-1/2)) | ||
| # print(xidudxi[i]*(4*np.pi*dens[i]**-1)**(-1/2)) | ||
|
|
||
| rho_c_rho_0.append(dens[i] ** -1) | ||
| dimless_arr.append(dimless_mass) |
There was a problem hiding this comment.
This loop can be vectorized using NumPy for improved performance and readability. This avoids manual iteration and is more idiomatic for NumPy arrays.
| rho_c_rho_0 = [] | |
| dimless_arr = [] | |
| for i in range(len(t)): | |
| dimless_mass = xidudxi[i] * (4 * np.pi * dens[i] ** -1) ** (-1 / 2) | |
| # print("rho_c/rho_0",dens[i]**-1) | |
| # print("xidudxi[-1]",xidudxi[i]) | |
| # print("(4*np.pi*dens[-1]**-1)**(1/2)",(4*np.pi*dens[i]**-1)**(-1/2)) | |
| # print(xidudxi[i]*(4*np.pi*dens[i]**-1)**(-1/2)) | |
| rho_c_rho_0.append(dens[i] ** -1) | |
| dimless_arr.append(dimless_mass) | |
| rho_c_rho_0 = dens**-1 | |
| dimless_arr = xidudxi * (4 * np.pi * dens**-1) ** (-1 / 2) |
| r_ = r_ | ||
| rho_ = rho_ |
There was a problem hiding this comment.
These assignments are redundant as they assign variables to themselves. They can be removed to improve code clarity.
| primitive = [] | ||
| accum = 0 | ||
| dx = np.diff(r_) | ||
| dx = np.concatenate([[dx[0]], dx]) | ||
| for r_val, rho_val, dx_val in zip(r_, rho_, dx): | ||
| r_val = r_val | ||
| primitive.append(accum) | ||
| # print(r_val, rho_val, accum, dx_val) | ||
| accum += rho_val * dx_val | ||
|
|
||
| primitive = np.array(primitive) |
There was a problem hiding this comment.
The manual calculation of the primitive (cumulative integral) is verbose and can be simplified. Using scipy.integrate.cumulative_trapezoid is more efficient, readable, and less error-prone. The line r_val = r_val is also a no-op and should be removed.
| primitive = [] | |
| accum = 0 | |
| dx = np.diff(r_) | |
| dx = np.concatenate([[dx[0]], dx]) | |
| for r_val, rho_val, dx_val in zip(r_, rho_, dx): | |
| r_val = r_val | |
| primitive.append(accum) | |
| # print(r_val, rho_val, accum, dx_val) | |
| accum += rho_val * dx_val | |
| primitive = np.array(primitive) | |
| primitive = scipy.integrate.cumulative_trapezoid(rho_, r_, initial=0) |
| m_s_codeu = codeu.get("m") * codeu.get("s", power=-1) | ||
| kg_m3_codeu = codeu.get("kg") * codeu.get("m", power=-3) | ||
|
|
||
| print(f"m_s_codeu: {kg_m3_codeu}") |
There was a problem hiding this comment.
There appears to be a copy-paste error in this print statement. The f-string variable name is m_s_codeu, but the value being printed is kg_m3_codeu.
| print(f"m_s_codeu: {kg_m3_codeu}") | |
| print(f"m_s_codeu: {m_s_codeu}") |
| # lattice volume | ||
| part_vol_lattice = 0.74 * part_vol | ||
|
|
||
| dr = (part_vol_lattice / ((4.0 / 3.0) * 3.1416)) ** (1.0 / 3.0) |
There was a problem hiding this comment.
Avoid using magic numbers like 3.1416 for pi. Use np.pi for better precision and code readability.
| dr = (part_vol_lattice / ((4.0 / 3.0) * 3.1416)) ** (1.0 / 3.0) | |
| dr = (part_vol_lattice / ((4.0 / 3.0) * np.pi)) ** (1.0 / 3.0) |
|
|
||
| dpi = 200 | ||
| plt.figure(num=1, clear=True, dpi=dpi) | ||
| import copy |
There was a problem hiding this comment.
import copy should be at the top of the file. This import is also repeated multiple times inside the analysis function (e.g., lines 460, 484, 509, 534). Imports should be at the top of the file for clarity and to avoid performance issues from re-importing within a loop. Please move this and other import copy statements to the top of the script.
| plt.show() | ||
|
|
||
|
|
||
| def analysis(i): |
There was a problem hiding this comment.
The analysis function contains a lot of duplicated code for plotting. The blocks for generating plots for arr_rho, arr_rho2, arr_cs, and arr_vxyz are nearly identical. This can be refactored into a helper function to reduce code duplication and improve maintainability. Also, the my_cmap colormap is created from scratch multiple times within the function; it could be created once and reused.
Example refactoring for plotting:
def plot_slice(data, cmap, title, filename, extent, dpi=200):
plt.figure(num=1, clear=True, dpi=dpi)
res = plt.imshow(
data,
cmap=cmap,
origin="lower",
extent=extent,
)
plt.xlabel("x")
plt.ylabel("y")
plt.title(title)
plt.colorbar()
plt.savefig(filename)
plt.close()
# In analysis function:
# ...
# my_cmap = ... (created once)
# title = f"t = {model.get_time():0.3f} [year]"
# plot_slice(arr_rho, my_cmap, title, f"collapse/rho_slice_{i:04}.png", [-ext, ext, -ext, ext])
# ...| i = 0 | ||
| t = 0 | ||
|
|
||
| steps = 10 | ||
| while True: | ||
| rho_max = analysis(i) | ||
|
|
||
| for j in range(steps): | ||
| model.timestep() | ||
|
|
||
| i += 1 | ||
| t += 2.0e-3 | ||
|
|
||
| if i > 1: | ||
| analysis(i) | ||
| break |
There was a problem hiding this comment.
The main simulation loop is confusing and contains an unused variable.
- The
while Trueloop with anif i > 1: breakcondition is hard to follow. Aforloop would be more explicit about the number of iterations. - The variable
tis initialized and incremented in the loop but never used. It should be removed.
Consider refactoring the loop for clarity.
| i = 0 | |
| t = 0 | |
| steps = 10 | |
| while True: | |
| rho_max = analysis(i) | |
| for j in range(steps): | |
| model.timestep() | |
| i += 1 | |
| t += 2.0e-3 | |
| if i > 1: | |
| analysis(i) | |
| break | |
| steps = 10 | |
| num_analysis_steps = 3 | |
| for i in range(num_analysis_steps): | |
| rho_max = analysis(i) | |
| if i < num_analysis_steps - 1: | |
| for _ in range(steps): | |
| model.timestep() |
1 participant
No description provided.