[SPH] Add Ormel variant of COALA source term#1830
Conversation
|
Thanks @mlombart 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 adds a new source file reference to the CMake configuration and includes a macro undefinition in the NodeEvolveDustCOALASourceTerm.hpp header. A critical issue was identified: the newly referenced source file, NodeEvolveDustCOALASourceTermOrmel.cpp, is missing from the pull request, which will result in a build failure.
|
/gemini review |
for more information, see https://pre-commit.ci
There was a problem hiding this comment.
Code Review
This pull request introduces the NodeEvolveDustCOALASourceTermOrmel module, implementing the Ormel model for differential velocities within the dust coagulation source term. Feedback identifies a critical correctness bug where stopping times are accessed without particle offsets, alongside performance concerns regarding redundant calculations in the kernel and inefficient buffer management. Additionally, the reviewer highlighted a potential division-by-zero risk and requested a refactor of the LaTeX documentation logic to address inaccuracies and code duplication.
| Tscal ts_i = ts[i]; | ||
| Tscal ts_j = ts[j]; |
There was a problem hiding this comment.
The stopping time ts is an IFieldSpan, which contains per-particle per-bin data. Accessing it as ts[i] and ts[j] without the particle offset id_a_d (which is id_a * nbins) will incorrectly use the stopping times of the first few particles in the patch for every particle in the simulation. This is a critical correctness bug.
Tscal ts_i = ts[id_a_d + i];
Tscal ts_j = ts[id_a_d + j];| Tscal t_L = cs_a / a_ext_a; | ||
| Tscal nh = rho_gas_a / (mu * mh); | ||
| Tscal Re = 62e6 * sycl::sqrt(nh / 1e5) * sycl::sqrt(T_a / 10.); | ||
| Tscal t_eta = t_L / sycl::sqrt(Re); |
There was a problem hiding this comment.
These variables (t_L, nh, Re, t_eta) are independent of the bin indices i and j. Calculating them inside the dv lambda results in
|
|
||
| Tscal eps_diff = 1e-16; | ||
|
|
||
| Tscal t_L = cs_a / a_ext_a; |
There was a problem hiding this comment.
If a_ext_a (the magnitude of external acceleration) is zero, this will result in a division by zero, setting t_L to infinity and potentially causing NaNs in subsequent calculations (e.g., in St or x_St). A safety guard or a small epsilon should be added.
Tscal t_L = cs_a / (a_ext_a + 1e-20);| sham::DeviceBuffer<Tscal> massgrid_buf(nbins + 1, dev_sched); | ||
| massgrid_buf.copy_from_stdvec(massgrid); | ||
|
|
||
| sham::DeviceBuffer<Tscal> tensor_tabflux_coag_buf(nbins * nbins * nbins, dev_sched); | ||
| tensor_tabflux_coag_buf.copy_from_stdvec(tensor_tabflux_coag); |
There was a problem hiding this comment.
Allocating and copying massgrid_buf and tensor_tabflux_coag_buf every time _impl_evaluate_internal is called (i.e., every timestep) is highly inefficient, especially for the
| auto rhodust_eps = get_ro_edge_base(0).get_tex_symbol(); | ||
| auto massgrid = get_ro_edge_base(1).get_tex_symbol(); | ||
| auto tensor_tabflux_coag = get_ro_edge_base(2).get_tex_symbol(); | ||
| auto part_counts = get_ro_edge_base(3).get_tex_symbol(); | ||
| auto s_j = get_ro_edge_base(5).get_tex_symbol(); | ||
| auto delta_v_j = get_ro_edge_base(6).get_tex_symbol(); | ||
| auto S_coag = get_rw_edge_base(0).get_tex_symbol(); | ||
|
|
||
| std::string tex = R"tex( | ||
| COALA dust coagulation source term, DG $k=0$ (Lombart et al., 2021) | ||
|
|
||
| Per gas particle $a$ and mass bin $j$ (monofluid: $\rho_{{\rm d},j,a}} = {s_j}_{j,a}^2$): | ||
|
|
||
| \begin{align} | ||
| \rho_{{\rm d},j,a} &= {s_j}_{j,a}^2 \\ | ||
| \Delta m_j &= {massgrid}_{j+1} - {massgrid}_j \\ | ||
| g_{j,a} &= \begin{cases} | ||
| \rho_{{\rm d},j,a} / \Delta m_j & \rho_{{\rm d},j,a} > \rho_{\rm eps} \\ | ||
| 0 & \text{otherwise} | ||
| \end{cases} \\ | ||
| \mathrm{dv}_{l,m,a} &= \left| {delta_v_j}_{m,a} - {delta_v_j}_{l,a} \right| \\ | ||
| \mathrm{flux}_{j,a} &= \sum_{l,m} | ||
| {tensor_tabflux_coag}_{j,l,m}\, | ||
| \mathrm{dv}_{l,m,a}\, g_{l,a}\, g_{m,a} \\ | ||
| {S_coag}_{0,a} &= -\mathrm{flux}_{0,a}, \quad | ||
| {S_coag}_{j,a} = \mathrm{flux}_{j-1,a} - \mathrm{flux}_{j,a} | ||
| \quad (j \ge 1) \\ | ||
| a &\in [0, {part_counts}), \quad j,l,m \in [0, N_{\rm bins}) \\ | ||
| \rho_{\rm eps} &= {rhodust_eps}, \quad N_{\rm bins} = {nbins} | ||
| \end{align} | ||
| )tex"; |
There was a problem hiding this comment.
The _impl_get_tex implementation appears to be a copy-paste from NodeEvolveDustCOALASourceTerm and contains several errors (incorrect indices, missing edges, outdated formulas). Refactor this duplicated logic into a helper function or lambda to improve readability and maintainability, as per repository guidelines.
References
- Refactor duplicated logic into a helper function or lambda to improve readability and maintainability.
Workflow reportworkflow report corresponding to commit 7017bb7 Pre-commit check report |
2 participants
No description provided.