Add wall heat load engineering module

[krystophny]

Summary

Adds standalone Python engineering module for wall heat load analysis from SIMPLE particle tracing results with energy-weighted flux calculation.

Physics features:

• Energy-weighted heat flux: q = sum(p²) * (P_alpha/N) / A_bin where p = v/v₀
• Accounts for collisional slowing down of alpha particles
• Wall area from actual 3D stellarator geometry via chartmap
• Monte Carlo error estimates: σ = q / √N_bin
• Distinguishes particle loss fraction vs energy loss fraction

Engineering capabilities:

• WallHeatMap: Loads results.nc and computes heat flux distribution (MW/m²)
• PortOptimizer: Optimizes port placement using scipy differential evolution
  • Exclusion zones (e.g., inboard side for NBI)
  • Max flux constraints
  • Multiple ports with non-overlap constraints
• compute_wall_area_from_chartmap(): Accurate 3D surface area from metric tensor
• plot_heat_flux_2d() / plot_heat_flux_3d(): Visualization with port overlays
• export_to_vtk(): VTK export for ParaView analysis

Coordinate handling:

• Maps VMEC theta [0, 2π] to heat map [-π, π]
• Handles cumulative zeta angles (applies modulo 2π)
• Proper zeta wrapping for ports near zeta = 0
• Documented conventions: theta=0 is outboard, ±π is inboard

Code quality:

• 1108 lines (slightly over soft limit but cohesive module)
• Comprehensive docstrings with physics background
• Proper error handling and input validation
• Warning when bin resolution exceeds chartmap resolution

Test coverage (38 tests):

• WallHeatMap: loading, binning, energy weighting, edge cases
• PortOptimizer: constraints, exclusion zones, max flux limits
• New features: bin_areas, flux_error, coordinate mapping, validation
• Wall area calculation: various geometries and aspect ratios
• Visualization and integration tests

Documentation:

• DOC/wall-heat-loads.md: Complete physics derivation with references
• Resolution matching guidance
• Guiding center limitations explained

Closes #284

Test plan

• All 38 Python tests pass
• Energy weighting verified against physics formulas
• Coordinate mapping tested with VMEC conventions
• Wall area matches analytic torus formula within 1.5%
• Resolution mismatch warning tested
• Reviewed by chris-architect (4 rounds, all issues resolved)

[qodo-code-review]

PR Compliance Guide 🔍

Below is a summary of compliance checks for this PR:

Security Compliance
🟢	No security concerns identified No security vulnerabilities detected by AI analysis. Human verification advised for critical code.
Ticket Compliance
🟡	🎫 #284 🟢 Provide a Python-based engineering workflow to compute and visualize alpha particle wall heat loads. Implement pysimple.engineering.WallHeatMap to bin wall hits on a (theta, zeta) grid and expose properties like total power, peak flux, and flux grid. Provide PyVista-based 3D visualization of heat flux on the wall surface and ability to overlay ports. Provide export functionality (e.g., VTK/ParaView export). Keep dependencies minimal and optional (e.g., pyvista, scipy, numpy), avoiding heavier geometry stacks. 🔴 Support loading SIMPLE output in the form described by the ticket (e.g., times_lost.dat or wall_hits.nc) containing per-particle wall hit info (including at least theta, zeta, and ideally x, y, z). Integrate chartmap-based wall tracing/geometry conversion (use chartmap file for coordinate conversion / wall surface generation; no CGAL requirement). Implement pysimple.engineering.PortOptimizer to optimize port/component placement in (theta, zeta) to minimize thermal exposure, supporting constraints such as max flux on ports and exclusion zones.
Codebase Duplication Compliance
⚪	Codebase context is not defined Follow the guide to enable codebase context checks.
Custom Compliance
🟢	Generic: Meaningful Naming and Self-Documenting Code Objective: Ensure all identifiers clearly express their purpose and intent, making code self-documenting Status: Passed Learn more about managing compliance generic rules or creating your own custom rules
	Generic: Secure Error Handling Objective: To prevent the leakage of sensitive system information through error messages while providing sufficient detail for internal debugging. Status: Passed Learn more about managing compliance generic rules or creating your own custom rules
	Generic: Secure Logging Practices Objective: To ensure logs are useful for debugging and auditing without exposing sensitive information like PII, PHI, or cardholder data. Status: Passed Learn more about managing compliance generic rules or creating your own custom rules
🔴	Generic: Robust Error Handling and Edge Case Management Objective: Ensure comprehensive error handling that provides meaningful context and graceful degradation Status: Missing input validation: The new code does not validate key boundary/edge inputs (e.g., trace_time could be zero/negative causing division-by-zero, and port widths can be non-positive causing invalid bounds/regions), leading to crashes or incorrect behavior without actionable error context. Referred Code wall_area = 4 * np.pi**2 * major_radius * minor_radius * 1e-4 bin_area = wall_area / (n_theta * n_zeta) eV_to_J = 1.602e-19 particle_energy_J = particle_energy_eV * eV_to_J total_energy_J = n_lost * particle_energy_J total_power_W = total_energy_J / trace_time total_power_MW = total_power_W * 1e-6 power_per_particle_MW = total_power_MW / n_lost if n_lost > 0 else 0 flux_grid = hit_count * power_per_particle_MW / bin_area return cls( theta_edges=theta_edges, zeta_edges=zeta_edges, flux_grid=flux_grid, hit_count=hit_count, total_power=total_power_MW, particle_energy_eV=particle_energy_eV, n_lost=n_lost, n_total=n_total, ... (clipped 248 lines) Learn more about managing compliance generic rules or creating your own custom rules
	Generic: Security-First Input Validation and Data Handling Objective: Ensure all data inputs are validated, sanitized, and handled securely to prevent vulnerabilities Status: Unvalidated external inputs: External inputs and user-controlled parameters (NetCDF contents, trace_time, n_theta/n_zeta, and port width/initial values) are not validated/sanitized, enabling invalid shapes/values (e.g., non-positive widths or zero trace_time) to propagate into unsafe computations and optimizer bounds. Referred Code def from_netcdf( cls, results_file: str | Path, n_theta: int = 64, n_zeta: int = 128, particle_energy_eV: float = 3.5e6, major_radius: float = 1000.0, minor_radius: float = 100.0, ) -> "WallHeatMap": """ Load results from NetCDF and compute heat flux distribution. Args: results_file: Path to results.nc from SIMPLE n_theta: Number of poloidal bins n_zeta: Number of toroidal bins particle_energy_eV: Alpha particle energy (default 3.5 MeV) major_radius: Major radius in cm (for area calculation) minor_radius: Minor radius in cm (for area calculation) Returns: ... (clipped 320 lines) Learn more about managing compliance generic rules or creating your own custom rules
⚪
Update

Compliance status legend

🟢 - Fully Compliant
🟡 - Partial Compliant
🔴 - Not Compliant
⚪ - Requires Further Human Verification
🏷️ - Compliance label

[qodo-code-review]

PR Code Suggestions ✨

Explore these optional code suggestions:

Category	Suggestion	Impact
Possible issue	Remove unsuitable gradient-based optimizer Remove the L-BFGS-B optimization option in _run_optimizer as it is unsuitable for the binned data, and exclusively use the more appropriate differential_evolution method. python/pysimple/engineering.py [413-436] def _run_optimizer( self, method: str, bounds: list, x0: list ) -> tuple[np.ndarray, bool, str]: """Run the selected optimization method.""" - from scipy.optimize import differential_evolution, minimize + from scipy.optimize import differential_evolution - if method == "differential_evolution": - result = differential_evolution( - self._objective, - bounds=bounds, - seed=42, - maxiter=100, - tol=1e-4, - polish=True, + if method != "differential_evolution": + raise ValueError( + "Only 'differential_evolution' method is supported for this problem." ) - else: - result = minimize( - self._objective, - x0=np.array(x0), - bounds=bounds, - method="L-BFGS-B", - ) + + result = differential_evolution( + self._objective, + bounds=bounds, + seed=42, + maxiter=100, + tol=1e-4, + polish=True, + ) return result.x, result.success, str(result.message) Apply / Chat Suggestion importance[1-10]: 9 __ Why: The suggestion correctly identifies that the L-BFGS-B optimizer is unsuitable for the piecewise constant objective function, preventing potential optimization failures and ensuring correct results.	High
	✅ Enable a disabled unit test Suggestion Impact: The commit removed the unconditional @pytest.mark.skipif(True, ...) that disabled test_solve_single_port and replaced it with a runtime conditional skip using pytest.importorskip("scipy"), so the test now runs when scipy is installed and is skipped otherwise. code diff: - @pytest.mark.skipif( - True, # Skip by default - requires scipy - reason="Optimization test requires scipy" - ) def test_solve_single_port(self, mock_results_file: Path): - """Test optimization with single port.""" - heat_map = WallHeatMap.from_netcdf(mock_results_file) + """Test optimization with single port finds low-flux region.""" + pytest.importorskip("scipy", reason="scipy required for optimization") + + heat_map = WallHeatMap.from_netcdf( + mock_results_file, total_alpha_power_MW=TEST_ALPHA_POWER_MW + ) opt = PortOptimizer(heat_map) opt.add_port("test_port", theta_width=0.3, zeta_width=0.3) @@ -227,10 +279,55 @@ assert result.ports[0].name == "test_port" assert -np.pi <= result.ports[0].theta_center <= np.pi assert 0 <= result.ports[0].zeta_center <= 2 * np.pi + # Optimized position should have less flux than center of heat map + center_flux = heat_map.integrated_flux_in_region(-0.15, 0.15, np.pi - 0.15, np.pi + 0.15) + assert result.total_flux_on_ports <= center_flux or result.total_flux_on_ports < 0.1 + Enable the test_solve_single_port test by changing the pytest.mark.skipif condition to be dependent on the availability of the scipy module. test/python/test_engineering.py [213-229] @pytest.mark.skipif( - True, # Skip by default - requires scipy + "scipy" not in sys.modules, reason="Optimization test requires scipy" ) def test_solve_single_port(self, mock_results_file: Path): """Test optimization with single port.""" heat_map = WallHeatMap.from_netcdf(mock_results_file) opt = PortOptimizer(heat_map) opt.add_port("test_port", theta_width=0.3, zeta_width=0.3) result = opt.solve() assert result.success assert len(result.ports) == 1 assert result.ports[0].name == "test_port" assert -np.pi <= result.ports[0].theta_center <= np.pi assert 0 <= result.ports[0].zeta_center <= 2 * np.pi [Suggestion processed] Suggestion importance[1-10]: 9 __ Why: The suggestion correctly identifies and fixes an unconditionally skipped test, ensuring that the core optimization functionality is actually validated by the test suite.	High
	✅ Implement unused chartmap file argument Suggestion Impact: The commit introduces substantial chartmap-based geometry support: it adds functions to compute wall surface area (and per-bin areas) from a chartmap NetCDF and wires a new `chartmap_file` argument into `WallHeatMap.from_netcdf` to use real wall area instead of a torus approximation. However, it does not implement the specific suggested change of using `chartmap_file` inside `plot_heat_flux_3d` to load X/Y/Z for the 3D mesh; the chartmap usage is applied to heat-flux scaling/area computation rather than the 3D plotting mesh construction. code diff: +def compute_wall_area_from_chartmap(chartmap_file: str | Path) -> float: + """ + Compute actual wall surface area from chartmap NetCDF file. + + The chartmap contains Cartesian coordinates (x, y, z) on a grid of + (rho, theta, zeta). The wall surface is at rho=1 (last index). + + Surface area is computed as: + A = n_fp * integral |dr/dtheta x dr/dzeta| dtheta dzeta + + where the cross product magnitude gives the local surface element. + + Args: + chartmap_file: Path to chartmap NetCDF file + + Returns: + Total wall surface area in m^2 + + Note: + This gives the actual 3D shaped stellarator wall area, which can be + 1.5-2x larger than the simple torus approximation 4*pi^2*R*a. + """ + if not HAS_NETCDF4: + raise ImportError("netCDF4 required: pip install netCDF4") + + with nc.Dataset(chartmap_file, 'r') as ds: + # Coordinates (theta and zeta are 1D, x/y/z are 3D) + theta = ds.variables['theta'][:] + zeta = ds.variables['zeta'][:] + nfp = int(ds.variables['num_field_periods'][...]) + + # Get boundary surface: x, y, z at rho=-1 (last index) + # File dims are (zeta, theta, rho), boundary is rho=-1 + x_wall = ds.variables['x'][:, :, -1] # shape (nzeta, ntheta) + y_wall = ds.variables['y'][:, :, -1] + z_wall = ds.variables['z'][:, :, -1] + + # Convert from cm to m + x_wall = x_wall * 0.01 + y_wall = y_wall * 0.01 + z_wall = z_wall * 0.01 + + nzeta, ntheta = x_wall.shape + dtheta = theta[1] - theta[0] if len(theta) > 1 else 2 * np.pi + dzeta = zeta[1] - zeta[0] if len(zeta) > 1 else 2 * np.pi / nfp + + # Compute surface element |dr/dtheta x dr/dzeta| at each grid point + # Use central differences for interior, one-sided at boundaries + # Note: theta is periodic (wraps), but zeta may not wrap if data is for + # one field period only (Cartesian x,y are not periodic in zeta within + # a single field period - they rotate by 2*pi/nfp between periods) + area = 0.0 + for iz in range(nzeta): + for it in range(ntheta): + # Theta: periodic, always use central difference with wrapping + it_p = (it + 1) % ntheta + it_m = (it - 1) % ntheta + dr_dtheta = np.array([ + (x_wall[iz, it_p] - x_wall[iz, it_m]) / (2 * dtheta), + (y_wall[iz, it_p] - y_wall[iz, it_m]) / (2 * dtheta), + (z_wall[iz, it_p] - z_wall[iz, it_m]) / (2 * dtheta), + ]) + + # Zeta: use one-sided differences at boundaries + # (Cartesian coords don't wrap within one field period) + if iz == 0: + # Forward difference at left boundary + dr_dzeta = np.array([ + (x_wall[1, it] - x_wall[0, it]) / dzeta, + (y_wall[1, it] - y_wall[0, it]) / dzeta, + (z_wall[1, it] - z_wall[0, it]) / dzeta, + ]) + elif iz == nzeta - 1: + # Backward difference at right boundary + dr_dzeta = np.array([ + (x_wall[iz, it] - x_wall[iz - 1, it]) / dzeta, + (y_wall[iz, it] - y_wall[iz - 1, it]) / dzeta, + (z_wall[iz, it] - z_wall[iz - 1, it]) / dzeta, + ]) + else: + # Central difference for interior + dr_dzeta = np.array([ + (x_wall[iz + 1, it] - x_wall[iz - 1, it]) / (2 * dzeta), + (y_wall[iz + 1, it] - y_wall[iz - 1, it]) / (2 * dzeta), + (z_wall[iz + 1, it] - z_wall[iz - 1, it]) / (2 * dzeta), + ]) + + # Cross product magnitude = surface element + cross = np.cross(dr_dtheta, dr_dzeta) + area += np.sqrt(np.sum(cross**2)) * dtheta * dzeta + + # Multiply by number of field periods to get full torus + return float(area * nfp) + + +def compute_bin_areas_from_chartmap( + chartmap_file: str | Path, + theta_edges: np.ndarray, + zeta_edges: np.ndarray, +) -> np.ndarray: + """ + Compute actual area of each (theta, zeta) bin from chartmap geometry. + + This properly accounts for the 3D stellarator shape rather than using + a uniform area approximation. + + Args: + chartmap_file: Path to chartmap NetCDF file + theta_edges: Bin edges for poloidal angle (n_theta + 1) + zeta_edges: Bin edges for toroidal angle (n_zeta + 1) + + Returns: + bin_areas: Array of shape (n_theta, n_zeta) with area in m^2 per bin + """ + if not HAS_NETCDF4: + raise ImportError("netCDF4 required: pip install netCDF4") + + with nc.Dataset(chartmap_file, 'r') as ds: + theta = ds.variables['theta'][:] + zeta = ds.variables['zeta'][:] + nfp = int(ds.variables['num_field_periods'][...]) + + x_wall = ds.variables['x'][:, :, -1] * 0.01 # cm -> m + y_wall = ds.variables['y'][:, :, -1] * 0.01 + z_wall = ds.variables['z'][:, :, -1] * 0.01 + + nzeta_cm, ntheta_cm = x_wall.shape + dtheta_cm = theta[1] - theta[0] if len(theta) > 1 else 2 * np.pi + dzeta_cm = zeta[1] - zeta[0] if len(zeta) > 1 else 2 * np.pi / nfp + + # Compute surface element at each chartmap grid point + surface_element = np.zeros((nzeta_cm, ntheta_cm)) + for iz in range(nzeta_cm): + iz_p = (iz + 1) % nzeta_cm + iz_m = (iz - 1) % nzeta_cm + for it in range(ntheta_cm): + it_p = (it + 1) % ntheta_cm + it_m = (it - 1) % ntheta_cm + + dr_dtheta = np.array([ + (x_wall[iz, it_p] - x_wall[iz, it_m]) / (2 * dtheta_cm), + (y_wall[iz, it_p] - y_wall[iz, it_m]) / (2 * dtheta_cm), + (z_wall[iz, it_p] - z_wall[iz, it_m]) / (2 * dtheta_cm), + ]) + dr_dzeta = np.array([ + (x_wall[iz_p, it] - x_wall[iz_m, it]) / (2 * dzeta_cm), + (y_wall[iz_p, it] - y_wall[iz_m, it]) / (2 * dzeta_cm), + (z_wall[iz_p, it] - z_wall[iz_m, it]) / (2 * dzeta_cm), + ]) + cross = np.cross(dr_dtheta, dr_dzeta) + surface_element[iz, it] = np.sqrt(np.sum(cross**2)) + + # Integrate surface element over each output bin + n_theta = len(theta_edges) - 1 + n_zeta = len(zeta_edges) - 1 + bin_areas = np.zeros((n_theta, n_zeta)) + + for i_theta in range(n_theta): + th_min, th_max = theta_edges[i_theta], theta_edges[i_theta + 1] + for i_zeta in range(n_zeta): + ze_min, ze_max = zeta_edges[i_zeta], zeta_edges[i_zeta + 1] + + # Find chartmap grid points in this bin and sum their contributions + area_sum = 0.0 + for iz in range(nzeta_cm): + zeta_val = zeta[iz] + # Handle periodicity: zeta in [0, 2pi/nfp) + if not (ze_min <= zeta_val < ze_max): + continue + for it in range(ntheta_cm): + theta_val = theta[it] + # Handle periodicity: theta in [-pi, pi) or [0, 2pi) + if not (th_min <= theta_val < th_max): + continue + area_sum += surface_element[iz, it] * dtheta_cm * dzeta_cm + + bin_areas[i_theta, i_zeta] = area_sum * nfp + + return bin_areas + + @dataclass class WallHeatMap: """ Wall heat flux distribution from SIMPLE particle tracing results. Bins lost particles by their wall hit location (theta, zeta) and computes - heat flux in MW/m^2 based on particle energy and wall area. + heat flux in MW/m^2. Accounts for collisional slowing down by using + energy-weighted binning based on final velocity p = v/v0. + + The heat flux calculation is: + q_bin = sum(p_i^2 for i in bin) * (P_alpha / N_total) / A_bin + + where: + - p_i = v/v0: normalized final velocity (from zend[3]) + - P_alpha: total alpha power (MW), user-provided + - N_total: total particles traced + - A_bin: area of this bin (m^2) Attributes: theta_edges: Bin edges for poloidal angle (rad) zeta_edges: Bin edges for toroidal angle (rad) flux_grid: Heat flux in MW/m^2, shape (n_theta, n_zeta) + energy_grid: Energy weight sum per bin (sum of p^2) hit_count: Number of particles per bin - total_power: Total deposited power in MW - particle_energy_eV: Energy per particle in eV + lost_power: Power lost to wall (MW) = energy_loss_fraction * total_alpha_power + total_alpha_power: Total alpha power from fusion (MW), user input n_lost: Number of lost particles n_total: Total number of particles traced trace_time: Tracing time in seconds + wall_area: Total wall area in m^2 + energy_loss_fraction: Fraction of energy lost = sum(p^2) / N_total """ theta_edges: np.ndarray zeta_edges: np.ndarray flux_grid: np.ndarray + energy_grid: np.ndarray hit_count: np.ndarray - total_power: float - particle_energy_eV: float + lost_power: float + total_alpha_power: float + energy_loss_fraction: float n_lost: int n_total: int trace_time: float @@ -66,30 +314,55 @@ minor_radius: float = 0.0 _wall_positions: Optional[np.ndarray] = field(default=None, repr=False) _loss_times: Optional[np.ndarray] = field(default=None, repr=False) + _final_p: Optional[np.ndarray] = field(default=None, repr=False) @classmethod def from_netcdf( cls, results_file: str | Path, + total_alpha_power_MW: float, n_theta: int = 64, n_zeta: int = 128, - particle_energy_eV: float = 3.5e6, major_radius: float = 1000.0, minor_radius: float = 100.0, + chartmap_file: Optional[str | Path] = None, ) -> "WallHeatMap": """ Load results from NetCDF and compute heat flux distribution. + The heat flux is computed using the Monte Carlo approach: + q_bin = sum(p_i^2 for i in bin) * (P_alpha / N_total) / A_bin + + This requires the user to specify total_alpha_power_MW based on the + plasma scenario. For D-T fusion: P_alpha = P_fusion / 5. + Args: results_file: Path to results.nc from SIMPLE + total_alpha_power_MW: Total alpha power from fusion reactions (MW). + For D-T: P_alpha = P_fusion / 5. E.g., 3 GW fusion -> 600 MW. n_theta: Number of poloidal bins n_zeta: Number of toroidal bins - particle_energy_eV: Alpha particle energy (default 3.5 MeV) - major_radius: Major radius in cm (for area calculation) - minor_radius: Minor radius in cm (for area calculation) + major_radius: Major radius in cm (for torus wall area approximation) + minor_radius: Minor radius in cm (for torus wall area approximation) + chartmap_file: Optional path to chartmap NetCDF for actual wall geometry. + If provided, uses real 3D surface area from metric tensor instead + of simple torus approximation. Returns: - WallHeatMap with binned heat flux data + WallHeatMap with physically-scaled heat flux in MW/m^2 + + Example: + # With actual wall geometry from chartmap + heat_map = WallHeatMap.from_netcdf( + "results.nc", + total_alpha_power_MW=600.0, + chartmap_file="wout.chartmap.nc", + ) + + Note: + Guiding center limitation: The reported wall positions are guiding + center positions, not actual wall impact points. For 3.5 MeV alphas, + the gyroradius is ~5-10 cm, which smears the heat load pattern. """ if not HAS_NETCDF4: raise ImportError("netCDF4 required: pip install netCDF4") @@ -109,64 +382,89 @@ lost_mask = class_lost n_lost = int(np.sum(lost_mask)) + theta_edges = np.linspace(-np.pi, np.pi, n_theta + 1) + zeta_edges = np.linspace(0, 2 * np.pi, n_zeta + 1) + + # Compute wall area: use chartmap if provided, else torus approximation + if chartmap_file is not None: + wall_area_m2 = compute_wall_area_from_chartmap(chartmap_file) + else: + # Simple torus: A = 4 * pi^2 * R0 * a, convert cm to m + wall_area_m2 = 4 * np.pi**2 * major_radius * minor_radius * 1e-4 + Implement the logic to use the chartmap_file argument in plot_heat_flux_3d to load exact wall geometry for more accurate 3D visualizations. python/pysimple/engineering.py [568-623] def plot_heat_flux_3d( heat_map: WallHeatMap, chartmap_file: Optional[str | Path] = None, cmap: str = "hot", clim: Optional[tuple[float, float]] = None, show: bool = True, ): """ Plot 3D heat flux on wall surface using PyVista. Args: heat_map: WallHeatMap to plot chartmap_file: Optional chartmap NetCDF for exact wall geometry cmap: Colormap name clim: Color limits (min, max) show: Whether to display the plot Returns: PyVista plotter object """ try: import pyvista as pv except ImportError: raise ImportError("pyvista required for 3D plots: pip install pyvista") - R0 = heat_map.major_radius * 1e-2 - a = heat_map.minor_radius * 1e-2 + if chartmap_file and Path(chartmap_file).exists(): + if not HAS_NETCDF4: + raise ImportError("netCDF4 is required to read chartmap_file") + with nc.Dataset(chartmap_file, 'r') as ds: + X = ds.variables['X'][:] + Y = ds.variables['Y'][:] + Z = ds.variables['Z'][:] + else: + R0 = heat_map.major_radius * 1e-2 + a = heat_map.minor_radius * 1e-2 - theta = heat_map.theta_centers - zeta = heat_map.zeta_centers - theta_grid, zeta_grid = np.meshgrid(theta, zeta, indexing='ij') + theta = heat_map.theta_centers + zeta = heat_map.zeta_centers + theta_grid, zeta_grid = np.meshgrid(theta, zeta, indexing='ij') - R = R0 + a * np.cos(theta_grid) - X = R * np.cos(zeta_grid) - Y = R * np.sin(zeta_grid) - Z = a * np.sin(theta_grid) + R = R0 + a * np.cos(theta_grid) + X = R * np.cos(zeta_grid) + Y = R * np.sin(zeta_grid) + Z = a * np.sin(theta_grid) grid = pv.StructuredGrid(X, Y, Z) grid["heat_flux"] = heat_map.flux_grid.flatten(order='F') plotter = pv.Plotter() plotter.add_mesh( grid, scalars="heat_flux", cmap=cmap, clim=clim, show_scalar_bar=True, scalar_bar_args={"title": "Heat Flux (MW/m^2)"}, ) plotter.add_axes() plotter.set_background("white") if show: plotter.show() return plotter [Suggestion processed] Suggestion importance[1-10]: 7 __ Why: The suggestion correctly identifies that the chartmap_file argument is unused and provides a reasonable implementation to add this functionality, making the 3D plot more accurate.	Medium
	Improve flux integration accuracy Refactor integrated_flux_in_region to use np.searchsorted on bin edges instead of checking bin centers, improving the accuracy of the power calculation. python/pysimple/engineering.py [206-220] def integrated_flux_in_region( self, theta_min: float, theta_max: float, zeta_min: float, zeta_max: float, ) -> float: """Compute total power deposited in a region (MW).""" - theta_mask = (self.theta_centers >= theta_min) & \ - (self.theta_centers <= theta_max) - zeta_mask = (self.zeta_centers >= zeta_min) & \ - (self.zeta_centers <= zeta_max) - region_hits = self.hit_count[np.ix_(theta_mask, zeta_mask)] + i_theta_min = np.searchsorted(self.theta_edges, theta_min, side='left') + i_theta_max = np.searchsorted(self.theta_edges, theta_max, side='right') + i_zeta_min = np.searchsorted(self.zeta_edges, zeta_min, side='left') + i_zeta_max = np.searchsorted(self.zeta_edges, zeta_max, side='right') + + region_hits = self.hit_count[i_theta_min:i_theta_max, i_zeta_min:i_zeta_max] fraction = np.sum(region_hits) / self.n_lost if self.n_lost > 0 else 0 return self.total_power * fraction Apply / Chat Suggestion importance[1-10]: 7 __ Why: The suggestion improves the accuracy of the flux integration by using bin edges (np.searchsorted) for slicing, rather than an approximation based on bin centers.	Medium
	✅ Check positive trace_time Suggestion Impact: The commit removed the division by trace_time entirely by changing the heat-flux/power computation to be based on user-supplied total_alpha_power_MW and energy-weighted binning, so the potential divide-by-zero risk from trace_time is implicitly eliminated rather than handled via an explicit trace_time > 0 check. code diff: - wall_area = 4 * np.pi**2 * major_radius * minor_radius * 1e-4 - bin_area = wall_area / (n_theta * n_zeta) - - eV_to_J = 1.602e-19 - particle_energy_J = particle_energy_eV * eV_to_J - total_energy_J = n_lost * particle_energy_J - total_power_W = total_energy_J / trace_time - total_power_MW = total_power_W * 1e-6 - - power_per_particle_MW = total_power_MW / n_lost if n_lost > 0 else 0 - flux_grid = hit_count * power_per_particle_MW / bin_area + # Energy-weighted histogram: sum of p^2 per bin + energy_grid, _, _ = np.histogram2d( + theta_lost, zeta_lost, + bins=[theta_edges, zeta_edges], + weights=energy_weight + ) + + # Energy-weighted heat flux calculation: + # q_bin = sum(p_i^2 for i in bin) * (P_alpha / n_total) / A_bin + bin_area = wall_area_m2 / (n_theta * n_zeta) + power_density = total_alpha_power_MW / n_total # MW per MC particle + flux_grid = energy_grid * power_density / bin_area # MW/m^2 + + # Energy loss fraction = sum(p^2 for all lost) / n_total + total_energy_lost = np.sum(energy_weight) + energy_loss_fraction = total_energy_lost / n_total + + # Lost power based on energy, not particle count + lost_power = energy_loss_fraction * total_alpha_power_MW In WallHeatMap.from_netcdf, add a check to ensure trace_time is positive before it is used in a division to prevent potential errors. python/pysimple/engineering.py [106-149] trace_time = float(ds.trace_time) + if trace_time <= 0: + raise ValueError(f"trace_time must be positive, got {trace_time}") ... total_energy_J = n_lost * particle_energy_J total_power_W = total_energy_J / trace_time [To ensure code accuracy, apply this suggestion manually] Suggestion importance[1-10]: 7 __ Why: The suggestion correctly identifies a potential division-by-zero error if trace_time is not positive and adds a necessary validation check to prevent it.	Medium
General	✅ Penalize exclusion zones Suggestion Impact: The commit updates PortOptimizer._objective to penalize ports whose (theta, zeta) center lies in an exclusion zone by adding a large 1e6 penalty to the objective value (total + penalty). This enforces exclusion zones through the objective, matching the suggestion’s intent, though implemented as an additive penalty rather than an immediate early return. code diff: def _objective(self, x: np.ndarray) -> float: - """Objective: minimize total flux on all ports.""" + """Objective: minimize total flux on all ports, penalize exclusion zones.""" total = 0.0 + penalty = 0.0 for i, port_spec in enumerate(self._ports): theta_c = x[2 * i] zeta_c = x[2 * i + 1] @@ -356,7 +655,10 @@ zeta_c - zeta_w / 2, zeta_c + zeta_w / 2, ) total += flux - return total + # Add large penalty for ports in exclusion zones + if self._in_exclusion_zone(theta_c, zeta_c): + penalty += 1e6 + return total + penalty In the _objective function, add a penalty to reject port positions that fall within defined exclusion zones by returning a large value. python/pysimple/engineering.py [346-348] def _objective(self, x: np.ndarray) -> float: + # penalize exclusion zone violations + for i in range(len(self._ports)): + theta_c, zeta_c = x[2*i], x[2*i+1] + if self._in_exclusion_zone(theta_c, zeta_c): + return 1e6 total = 0.0 [Suggestion processed] Suggestion importance[1-10]: 6 __ Why: The suggestion correctly proposes to enforce exclusion zones within the optimizer's objective function, which is a more robust approach than relying solely on optimizer bounds.	Low
	✅ Enforce max flux constraint Suggestion Impact: The commit adds enforcement of _max_flux_constraint, but not via a post-solve success/message update. Instead, it computes the peak flux density over port regions during optimization and applies a penalty when the peak exceeds _max_flux_constraint, effectively pushing the optimizer away from violating solutions. code diff: def _objective(self, x: np.ndarray) -> float: - """Objective: minimize total flux on all ports.""" + """Objective: minimize total flux on all ports, penalize constraints.""" total = 0.0 + penalty = 0.0 + max_flux_density = 0.0 + for i, port_spec in enumerate(self._ports): theta_c = x[2 * i] zeta_c = x[2 * i + 1] theta_w = port_spec["theta_width"] zeta_w = port_spec["zeta_width"] + + theta_min = theta_c - theta_w / 2 + theta_max = theta_c + theta_w / 2 + zeta_min = zeta_c - zeta_w / 2 + zeta_max = zeta_c + zeta_w / 2 + flux = self.heat_map.integrated_flux_in_region( - theta_c - theta_w / 2, theta_c + theta_w / 2, - zeta_c - zeta_w / 2, zeta_c + zeta_w / 2, + theta_min, theta_max, zeta_min, zeta_max ) total += flux - return total + + # Check peak flux density on this port for constraint + if self._max_flux_constraint is not None: + port_flux_density = self._get_peak_flux_in_region( + theta_min, theta_max, zeta_min, zeta_max + ) + max_flux_density = max(max_flux_density, port_flux_density) + + # Penalty for port overlapping exclusion zones (check all corners + center) + if self._port_overlaps_exclusion(theta_c, zeta_c, theta_w, zeta_w): + penalty += 1e6 + + # Penalty for exceeding max flux constraint + if self._max_flux_constraint is not None: + if max_flux_density > self._max_flux_constraint: + excess = max_flux_density - self._max_flux_constraint + penalty += 1e4 * excess # Proportional penalty + + return total + penalty After optimization in the solve method, check if result.max_flux_on_ports violates the _max_flux_constraint and update the result's success status accordingly. python/pysimple/engineering.py [368-394] def solve(self, method: str = "differential_evolution") -> OptimizationResult: ... x_opt, success, message = self._run_optimizer(method, bounds, x0) -- return self._build_result(x_opt, success, message) + result = self._build_result(x_opt, success, message) + if self._max_flux_constraint is not None and result.max_flux_on_ports > self._max_flux_constraint: + result.success = False + result.message += " Max flux constraint violated." + return result [To ensure code accuracy, apply this suggestion manually] Suggestion importance[1-10]: 6 __ Why: The suggestion correctly implements the check for the _max_flux_constraint, which was set but never enforced, making the feature complete and functional.	Low
Update