[cpinte]

Type of PR:
Bug fix / modification to existing code

Description:
Modified Random Walk (MRW) for the temperature calculation

Adds an optional Modified Random Walk (Min et al. 2009, Robitaille 2010) to the
Monte Carlo temperature calculation: in optically thick cells, packets move by
diffusion steps to the cell wall instead of individual interactions. Energy is
deposited with a dB/dT-weighted absorption opacity, consistent with the
Bjorkman & Wood re-emission.

• Options: -mrw [gamma] or -gamma_mrw gamma (default gamma = 5), -no_mrw.
  Off by default.
• Works on cylindrical, spherical and Voronoi grids; LTE grains only (switched
  off with non-equilibrium grains or PAHs).
• Pinte et al. 2009 benchmark, temperature step: 2.6-3.3x faster at
  tau = 1e5-1e6, 2.4x at tau = 1e9. Temperatures agree with standard MC
  within noise, except the last midplane cells before the dark zone (2-4 % low
  at gamma = 5). With -diff_approx, this offset propagates into the dark zone
  (-0.7 to -1.8 % on average): use gamma = 10 or no MRW for benchmark-quality
  dark-zone temperatures.
• Unit tests (random-walk sampling, wall distances on all grids, mean
  opacities) and integration tests against standard MC.

Testing:
Updated unit and integration tests

Did you run the botscheck that the code and comments follow the code of conduct? yes

Did you update relevant documentation in the docs directory? yes

[github-actions]

Code Review by Gemini

Summary of Review

The pull request introduces the Modified Random Walk (MRW) algorithm (Min et al. 2009; Robitaille 2010) to accelerate Monte Carlo temperature calculations in optically thick grid regions across cylindrical, spherical, and 3D Voronoi geometries.

Overall, this is an exceptionally high-quality PR:

• Physical Correctness: Accurately implements the $dB/dT$-weighted Rosseland transport mean extinction ($\chi_{\text{trans}}$) and absorption opacity ($\kappa_{\text{abs, dB/dT}}$) consistent with the Bjorkman & Wood (2001) re-emission scheme.
• Numerical Stability: Uses the Jacobi $\theta_4$ Poisson summation identity for $y > 0.5$ in evaluating the exit probability $P(t)$, eliminating numerical cancellation as $y \to 1$.
• Parallel & Memory Safety: OpenMP loops and live energy deposition into xKJ_abs and xN_abs strictly preserve thread safety using thread-indexed data structures (id).
• Bug Fixes: Fixes pre-existing bugs in distance_to_closest_wall_sph (which previously misindexed cos_phi_lim with theta indices) and distance_to_closest_wall_Voronoi (which previously omitted the vector norm in the bisector distance denominator).
• Verification: Includes extensive unit tests (test_mrw.f90) and Python integration tests (test_mrw.py).

---

Detailed Findings

1. Numerical Stability & Zeta Sampling

• Function: cumulative_zeta_from_y & initialize_cumulative_zeta (src/MRW.f90)
• Assessment: The dual-series representation (direct series for $y \le 0.5$ and Jacobi $\theta_4$ transformation for $y > 0.5$) ensures full double-precision accuracy without loss of precision as $y \to 1$. The $O(1)$ inverse lookup table interp_y_from_zeta clamps inputs away from 0 and 1, ensuring -log(y) remains finite and strictly positive.

2. Grid Geometry & Nearest-Wall Distance

• Files: src/cylindrical_grid.f90, src/spherical_grid.f90, src/Voronoi.f90
• Assessment:
  • distance_to_closest_wall_sph now correctly uses cos_theta_lim(ja) with $1 \le j_a \le n_z$, resolving array out-of-bound errors on 3D or fine-theta spherical grids.
  • distance_to_closest_wall_Voronoi correctly calculates exact Euclidean distance dot_product(n, p - r) / norm_n to bisector planes and handles stellar surface boundaries.
  • All coordinate distance functions safely return 0.0_dp for boundary/ghost cells, disabling MRW at domain edges where path diffusion approximations are invalid.

3. OpenMP & Thread Safety

• Files: src/thermal_emission.f90, src/dust_transfer.f90
• Assessment:
  • init_MRW_opacities uses OpenMP collapse(2) with explicit private loop variables and shared read-only arrays.
  • Photon propagation during MRW steps deposits energy via xKJ_abs(icell, id) and updates packet counters xN_abs(icell, 1, id), where id is the thread index, eliminating race conditions.

---

Conclusion

The PR is mathematically sound, numerically stable, thread-safe, and thoroughly tested. It is safe to merge.

[github-actions]

Code Review by Gemini

Physical / Numerical Purpose

The changes introduce a Modified Random Walk (MRW) acceleration scheme (Min et al. 2009; Robitaille 2010) to the Monte Carlo thermal equilibrium calculation. In optically thick cells ($\tau \gg 1$), packet transport by individual scattering steps is replaced by macroscopic diffusion steps inside spheres bounded by the nearest cell wall. Deposited energy is accumulated continuously using $dB/dT$-weighted Planck mean absorption opacities ($\kappa_{\text{abs,MRW}}$), and exiting packets are resampled at the local updated dust temperature with an anisotropic emergent-flux exit distribution ($p(\mu)$) and Chandrasekhar extrapolation length ($z_{\text{extrap}} = 0.7104$).

---

Code Review Findings

1. High-Performance & Parallel Computing

• Thread Safety (dust_transfer.f90, MRW.f90, thermal_emission.f90):
  • Random number generation inside make_MRW_step and sample_exit_direction uses thread-indexed SPRNG streams (stream(id)), avoiding RNG race conditions across OpenMP threads.
  • Energy deposition (xKJ_abs(icell, id)) and cell hit counters (xN_abs(icell, 1, id)) are accumulated per thread id, maintaining OpenMP race-freedom.
• Cache Locality (thermal_emission.f90, mem.f90):
  • Arrays chi_Rosseland_MRW(n_T, p_Nc) and kappa_abs_MRW(n_T, p_Nc) are allocated with n_T as the fast (first) dimension. Inner loops over temperature (do T = 1, n_T) access contiguous memory columns, adhering to Fortran column-major order.
• OpenMP Parallelization (thermal_emission.f90):
  • The init_MRW_opacities loop is cleanly parallelized over (icell, T) via !$omp parallel do collapse(2) with all thread-local scalars declared private.

2. Numerical Precision & Stability

• Poisson Summation Identity (MRW.f90: cumulative_zeta_from_y):
  • Replacing the truncated power series for $y > 0.5$ with the Jacobi $\vartheta_4$ Poisson summation identity ($\sum_k \exp(-\pi^2 (k+1/2)^2 / (-\ln y))$) eliminates catastrophic numerical cancellation as $y \to 1$.
• Floating-Point Literals:
  • Double-precision constants and literals are consistently specified using _dp (e.g. 5.0_dp, 0.7104_dp, 1.0e-5_dp), preventing accidental single-precision implicit conversions.
• Division by Zero & Boundary Protections:
  • Proper guards are in place for zero total opacity (chi_trans_wl > tiny_dp), zero cell volume, zero normal norms (norm_n <= tiny_dp), and edge cases where $y \to 0$ or $y \to 1$.

3. Grid Geometry & Radiative Transfer Physics

• Exact Nearest-Wall Distances (cylindrical_grid.f90, spherical_grid.f90, Voronoi.f90):
  • Cylindrical/Spherical: Bounds distances accurately against radial, elevation cone, and azimuthal plane boundaries. Correctly handles 2D midplane cell symmetry (s4 = huge(1.0_dp) for $j=1$).
  • 3D Voronoi: Distance to bisector planes is computed via signed projection $(p - r) \cdot \hat{n}$. Cells on outer domain boundaries (id_n <= 0), cut cells, or star-containing cells safely return 0.0_dp to disable MRW near boundaries.
  • Coordinate Fix: Fixed a bug where sin_phi_lim evaluated to 1.0d300 instead of $\pm 1$ at $\phi = \pi/2$ or $3\pi/2$.
• Hemisphere Transition Fix (spherical_grid.f90):
  • Midplane crossings at $z \approx 0$ assign hemisphere index thetaj1 based on direction cosine sign (w < 0), avoiding grid mapping lockup near boundary cones.

4. Testing & Verification

• Extensive Fortran unit test suite (test_suite/unit/test_mrw.f90) verifies:
  • Analytical diffusion kernel properties and moments ($\langle -\ln y \rangle = \pi^2 / 6$).
  • Nearest-wall distances against 1000-ray sweeps and sphere containment on 2D/3D Cylindrical, Spherical, and Voronoi grids.
  • $dB/dT$-weighted mean opacities against analytical Riemann-zeta solutions ($\beta = 2$ power-law opacity).
• Pytest suite (test_suite/test_mrw.py) validates CLI flags, temperature accuracy, and SED consistency against standard Monte Carlo.

---

Recommendations / Minor Suggestions

1. Thread-Safety Guard for Lookup Table Initialization (src/MRW.f90):
   In initialize_cumulative_zeta(), lzeta_initialized is a saved variable:

   logical, save :: lzeta_initialized = .false.

   init_MRW_cells() currently calls this sequentially during pre-computation before transport starts. To guard against potential future invocation inside a parallel region, consider adding an $omp single or if (lzeta_initialized) return short-circuit before any table writes.

---

Summary & Conclusion

• Overall Assessment: Excellent, high-quality PR. The mathematical derivation, parallel execution, numerical safeguards, and grid geometry logic are rigorous and complete.
• Key Strengths: Poisson summation identity for $\zeta(y)$, comprehensive distance-to-wall calculations across all grid geometries, and exemplary Fortran unit test coverage.
• Merge Recommendation: Safe to merge.

[github-actions]

ERROR: # Code Review by Gemini
An unexpected error occurred in Client call: 503 UNAVAILABLE. {'error': {'code': 503, 'message': 'This model is currently experiencing high demand. Spikes in demand are usually temporary. Please try again later.', 'status': 'UNAVAILABLE'}}

[github-actions]

Code Review by Gemini

Here is a code review for the Modified Random Walk (MRW) pull request.

---

Overall Assessment

The PR introduces a Modified Random Walk (MRW) implementation to accelerate Monte Carlo radiative transfer in optically thick regions, along with updates to wall-distance algorithms across cylindrical, spherical, and Voronoi grids. The core Fortran MRW logic and numerical transport routines are well-designed and covered by comprehensive unit/integration tests.

However, there are a few issues—most notably a truncated CMakeLists.txt file, missing/misnamed C++ header files, and a loop-index bug when reallocating star arrays in hydro readers that overwrites existing star accretion properties.

---

Key Findings & Action Items

1. Truncated CMake File

• Severity: CRITICAL
• File: src/CMakeLists.txt (line 40)
• Problem: src/CMakeLists.txt is incomplete and ends abruptly on line 40 with set(CMAKE.
• Impact: Any build attempting to configure using CMake will fail immediately with a syntax error.
• Suggested Fix: Complete or restore the CMakeLists.txt file before merging.

---

2. Missing C++ Header (mcfost_interface.hpp) and Misnamed Source File

• Severity: HIGH
• Files: src/interface_template.cpp (line 1), src/mcfost2hydro.cpp (lines 1–2)
• Problem:
  1. interface_template.cpp includes "mcfost_interface.hpp", but this header file is missing from the repository diff.
  2. src/mcfost2hydro.cpp starts with header guards (#ifndef MCFOST_INTERFACE_HPP ... #define MCFOST_INTERFACE_HPP) and contains standard header definitions for namespace MCFOST, but has a .cpp file extension.
• Impact: C++ builds or interface consumers will fail to compile due to missing header files.
• Suggested Fix: Rename src/mcfost2hydro.cpp to src/mcfost_interface.hpp (or add the missing mcfost_interface.hpp header) and adjust the include path accordingly.

---

3. Overwriting Existing Star Properties During Reallocation in Hydro Readers

• Severity: HIGH (Correctness / Data Corruption)

• Files: src/read_fargo3d.f90, src/read_athena++.f90, src/read_pluto.f90

• Problem: When reallocating the star array for new sink particles, the loop index for initializing new stars starts at n_stars_old instead of n_stars_old + 1:

  ! Saving if the accretion rate was forced
  allocate(star_old(n_stars_old))
  if (allocated(star)) then
     star_old(:) = star(:)
     deallocate(star)
  endif
  allocate(star(n_stars))
  do i=1, min(n_stars, n_stars_old)
     star(i)%force_Mdot = star_old(i)%force_Mdot
     star(i)%Mdot = star_old(i)%Mdot
  enddo
  ! If we have new stars
  do i=n_stars_old,n_stars   ! <--- BUG: starts at n_stars_old instead of n_stars_old + 1
     star(i)%force_Mdot = .false.
     star(i)%Mdot = 0.
  enddo
  deallocate(star_old)

• Impact: For i = n_stars_old, the second loop immediately overwrites the force_Mdot and Mdot properties preserved in the first loop. (Note that src/density.f90 correctly uses do i=n_stars_old+1, n_stars).

• Suggested Fix: Change the loop lower bound in read_fargo3d.f90, read_athena++.f90, and read_pluto.f90:

  do i = n_stars_old + 1, n_stars
     star(i)%force_Mdot = .false.
     star(i)%Mdot = 0.0
  enddo

---

4. Unintended Planet Initial Velocity in Athena++ Reader

• Severity: MEDIUM (Physical Validity)

• File: src/read_athena++.f90

• Problem: In read_athena++.f90, the vertical velocity component vz for the hardcoded planet properties was changed to 1.0:

  - real, parameter :: vx=0.0, vy=1.0, vz=0.0 ! circular orbit in the midplane (vz was 1)
  + real, parameter :: vx=0.0, vy=1.0, vz=1.0

• Impact: Giving a midplane planet ($z=0$) a non-zero vertical velocity $v_z = 1.0 \times u_{\text{velocity}}$ ($\approx 30\text{ km/s}$) launches the planet out of the midplane, contradicting the inline comment and standard midplane circular orbit assumption.

• Suggested Fix: Revert vz to 0.0_dp:

  real, parameter :: vx=0.0_dp, vy=1.0_dp, vz=0.0_dp

---

5. OpenMP Thread Safety in initialize_cumulative_zeta

• Severity: LOW
• File: src/MRW.f90 (lines 74–105)
• Problem: initialize_cumulative_zeta uses logical, save :: lzeta_initialized = .false. and modifies shared module variables y_MRW and zeta.
• Impact: While initialize_cumulative_zeta is currently invoked during serial setup (init_MRW_cells), if it is ever called from inside an OpenMP parallel region, simultaneous initialization attempts by multiple threads will result in a data race on y_MRW and zeta.
• Suggested Fix: Guard initialization with an OpenMP critical section or single construct, or ensure initialization occurs during serial startup before parallel loops execute.

---

Summary

The PR introduces a valuable performance feature (Modified Random Walk) with strong unit and integration test coverage. Once src/CMakeLists.txt, the C++ header naming, and the hydro reader star loop indices are corrected, the PR should be ready to merge.