Skip to content

Memory test PR - #3

Open
jhmatthews wants to merge 36 commits into
jhmatthews:memory-formattingfrom
sirocco-rt:memory
Open

Memory test PR#3
jhmatthews wants to merge 36 commits into
jhmatthews:memory-formattingfrom
sirocco-rt:memory

Conversation

@jhmatthews

Copy link
Copy Markdown
Owner

No description provided.

kslong and others added 30 commits March 19, 2026 11:06
shared/private MPI memory model

This is Phase 1 of a plan to reduce memory usage when running
Sirocco with many MPI ranks on a single node. Currently every rank
holds a complete copy of plasmamain and macromain (~700 MB per rank
for a 10K cell macro-atom model). The goal is to enable MPI-3
shared memory windows so that ranks on the same node share one
read-only copy of plasma state, while each rank maintains its own
private copy of the estimator arrays that are accumulated during
photon transport.

The key difficulty is that during photon transport, all ranks write
estimator fields (j, heat_tot, heat_lines, etc.) to any cell a
photon passes through, so the plasma struct cannot be naively placed
in shared memory. The solution is to separate fields by their MPI
communication pattern, so that in a future phase the read-only and
accumulated portions can be allocated in different memory regions.

STRUCT REORGANIZATION (source/sirocco.h):

The monolithic plasma_dummy struct (~100 fields) has been split into
three nested sub-structures, categorized by how each field is used
during the MPI-parallel photon transport phase:

  plasma_state (read-only during transport):
    Fields that are set during initialization or wind updates and
    only read during photon transport. Includes ne, rho, vol,
    xgamma, t_r, t_e, w, density[], partition[], levden[],
    recomb_simple[], kbf_use[], spectral model parameters
    (pl_alpha, exp_temp, spec_mod_type, etc.), and kappa_ff_factor.
    In a future shared-memory phase, this sub-struct will be placed
    in an MPI shared window visible to all ranks on a node.

  plasma_estimators (accumulated during transport, reduced across ranks):
    Fields that every rank increments via += as photons traverse
    cells. After transport, these are summed across ranks by
    reduce_simple_estimators() using MPI_Allreduce. Includes j,
    j_direct, j_scatt, ave_freq, ip, heat_tot, heat_lines,
    heat_ff, heat_comp, heat_photo, cool_tot, kpkt_abs, xj[],
    xave_freq[], rad_force_es[], F_vis[], F_UV[], cell_spec_flux[],
    ioniz[], heat_ion[], and photon passage counters (ntot, nioniz,
    etc.). In a future phase, each rank will hold a private copy of
    this sub-struct.

  plasma_derived (computed during wind updates, broadcast to ranks):
    Fields computed from the estimators during the wind update
    phase, then broadcast to all ranks. Includes cooling rates
    (cool_lines, cool_comp, cool_di, cool_dr, etc.), luminosities
    (lum_lines, lum_ff, lum_rr, lum_tot, etc.), ionization-band
    quantities, convergence diagnostics (converge_t_r, converge_t_e,
    gain, etc.), scatter counters, persistent radiation force
    averages, and per-ion derived arrays (scatters[], cool_rr_ion[],
    lum_rr_ion[], recomb[], etc.).

The top-level struct becomes:
  typedef struct plasma {
      int nwind, nplasma;
      struct plasma_state state;
      struct plasma_estimators est;
      struct plasma_derived derived;
  } plasma_dummy;

The same treatment was applied to macro_dummy, split into:
  macro_state     - jbar_old, gamma_old, gamma_e_old, alpha_st_old,
                    alpha_st_e_old, store_matom_matrix, matom_transition_mode
  macro_estimators - jbar, gamma, gamma_e, alpha_st, alpha_st_e,
                    recomb_sp, recomb_sp_e, matom_abs, cooling_bf,
                    cooling_bf_col, cooling_bb, and cooling totals
  macro_derived   - matom_emiss, kpkt_rates_known, matom_matrix,
                    matrix_rates_known, cooling_bb_simple_tot

Several #define constants and enum spec_mod_type_enum were moved
from inside the plasma_dummy struct body to file scope, as required
for the sub-struct definitions.

MECHANICAL REFACTORING (55 source files):

All field accesses across the codebase were updated:
  xplasma->ne          -> xplasma->state.ne
  xplasma->j           -> xplasma->est.j
  xplasma->lum_tot     -> xplasma->derived.lum_tot
  plasmamain[n].density -> plasmamain[n].state.density
  macromain[n].jbar    -> macromain[n].est.jbar
  (and so on for all fields)

Eight field names conflicted with fields in other structs (WindPtr,
PhotPtr, ions, etc.): w, vol, f1, f2, ip, xi, nbands, xgamma.
For these, only accesses through known PlasmaPtr variable names
(xplasma, plasmamain[...]) were modified, to avoid false positives.

MPI COMMUNICATION (communicate_plasma.c, communicate_macro.c):

The MPI_Pack/MPI_Unpack sequences in broadcast_plasma_grid(),
reduce_simple_estimators(), broadcast_updated_plasma_properties(),
broadcast_wind_luminosity(), broadcast_wind_cooling(), and the
macro equivalents were all updated to use the new field paths.
The pack/unpack order and buffer size constants (N_BASIC_DOUBLES=73,
N_BASIC_INTS=22) are unchanged — this is a purely structural
refactoring with no change to the communication protocol.

DYNAMIC MEMORY (source/gridwind.c):

calloc_dyn_plasma() and calloc_estimators() updated to allocate
arrays within the appropriate sub-structs (e.g., plasmamain[n].
state.density, macromain[n].est.jbar).

WIND SAVE I/O (source/windsave.c):

fwrite/fread calls for dynamic arrays updated to new field paths.
Note: the binary wind save format changes because sizeof(plasma_dummy)
and sizeof(macro_dummy) changed due to sub-struct padding. Old wind
save files are incompatible with this version.

VERSION BUMP (source/Makefile):

VERSION changed from 1.2 to 2.0 to reflect the struct layout change
and wind save incompatibility.

BUG FIXES:

source/signal.c - Fix race condition in xsignal_rm():
  Replaced fopen() existence check followed by system("rm ...") with
  a single remove() call that ignores ENOENT. The old code had a
  TOCTOU race when multiple processes called xsignal_rm concurrently.

py_progs/run_check.py - Fix xwindsave2table() version fallback:
  Fixed three bugs that prevented the fallback from working:
  (1) Header check now accepts 'Sirocco' in addition to 'Python'.
  (2) Binary name now includes hyphen: windsave2table-VERSION.
  (3) Fixed syntax error in print statement on line 158.

py_progs/regression_check.py - Fall back to .spec for model discovery:
  When .out.pf files are missing, the comparison now discovers models
  from .spec files instead of failing with zero models.

py_progs/regression_plot.py - Same .spec fallback for plot generation:
  Applied the identical fallback so comparison plots are generated
  even when .out.pf files are absent.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
…totypes

Convert ~300 K&R-style function definitions to ANSI C style across 102
source files, moving parameter type declarations into the function
signature. Remove ~90 redundant local function declarations with empty
parentheses (e.g., `double dot ();`) that are already prototyped in
headers. For standalone files not including sirocco.h, replace removed
declarations with properly typed forward declarations.

Additional fixes:
- Remove -Wno-deprecated-non-prototype from Makefile (no longer needed)
- Add #include <sys/stat.h> for mkdir() in parse.c and inspect_wind.c
- Fix type mismatch in swind.c: char choice -> int choice (matches templates.h)
- Initialize fptr=NULL in swind_sub.c to silence uninitialized warning
- Remove unused f_base/f_cat variables in reverb.c

This eliminates all Clang -Wdeprecated-non-prototype warnings and
brings the codebase closer to C99/C23 compliance. The build now
produces zero warnings across all targets (sirocco, swind, windsave2table,
windsave2fits, inspect_wind, modify_wind, etc.).

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
Document the three-way split of plasma_dummy and macro_dummy into
state/est/derived sub-structures in programmer_notes.rst and
mpi_comms.rst, including field access patterns, communication
function mappings, and guidance for adding new variables.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
This commit implements Phase 2 of the shared-memory optimization for
SIROCCO.  When running with multiple MPI ranks on the same node,
variable-length plasma and macro-atom arrays now use MPI-3 shared
memory windows (MPI_Win_allocate_shared) so that only one physical
copy exists per node, rather than one per rank.  This significantly
reduces memory consumption for large models on multi-core nodes.

Phase 2a — contiguous block allocation:

All dynamic plasma arrays (density, partition, levden, ioniz, etc.)
are now allocated as single contiguous blocks in calloc_dyn_plasma()
and calloc_estimators(), with each cell's pointer set to the correct
offset within the block.  This replaces the previous pattern of
separate calloc() calls per cell and is a prerequisite for shared
memory, which requires contiguous regions.  Block base pointers and
MPI_Win handles are tracked in new global structs plasma_block_ptrs
(type plasma_blocks) and macro_block_ptrs (type macro_blocks),
declared in sirocco.h.

Phase 2b — MPI-3 shared memory windows:

A node-local communicator (node_comm) is created during MPI init
via MPI_Comm_split_type(MPI_COMM_TYPE_SHARED).  The helper functions
alloc_block_double() and alloc_block_int() accept a use_shared flag:
when TRUE and np_mpi_global > 1, only the node leader (node_rank==0)
allocates memory via MPI_Win_allocate_shared; other ranks obtain a
pointer to the same physical memory via MPI_Win_shared_query.

The allocation strategy follows the sub-structure split:

  - State arrays (density, partition, levden, recomb_simple, etc.)
    are shared — read-only during photon transport.
  - Estimator arrays (ioniz, heat_ion, heat_inner_ion, inner_ioniz)
    are always private — each rank accumulates independently.
  - Derived arrays (recomb, cool_rr_ion, lum_rr_ion, cool_dr_ion,
    inner_recomb) are shared — computed during wind updates.
  - scatters and xscatters are private despite being derived,
    because they are incremented during photon transport.

The same pattern applies to macro-atom arrays in calloc_estimators().

Race condition fix in sobolev():

The sobolev() function in resonate.c previously modified
xplasma->state.density[nion] temporarily during photon transport
to pass an interpolated density to two_level_atom().  With shared
memory, this created a race condition where other ranks could read
the temporarily corrupted value.  Fixed by adding a density_override
parameter to two_level_atom(): when >= 0, it overrides the density
read from xplasma->state.density, eliminating the need to modify
shared state.  All callers pass -1.0 for normal behaviour.

Cleanup and synchronisation:

MPI_Barrier(node_comm) is called after every broadcast that writes
to shared dynamic arrays, ensuring visibility to all node-local
ranks.  Cleanup at program exit in janitor.c frees contiguous blocks
via the block pointer structs; shared blocks are simply NULLed (the
MPI runtime frees them at MPI_Finalize), while private blocks are
freed with free().

Documentation:

Updated mpi_comms.rst to replace the "Future shared-memory model"
placeholder with a complete description of the implementation.
Updated programmer_notes.rst to reflect the current shared-memory
model and document thread-safety constraints on state arrays.
Updated doxygen headers for two_level_atom(), sobolev(), and
calloc_dyn_plasma().

Files modified:
  source/sirocco.c          — node_comm setup via MPI_Comm_split_type
  source/sirocco.h          — plasma_blocks/macro_blocks structs, MPI_Win handles
  source/sirocco_extern_init.c — node_comm, node_rank, node_size globals
  source/gridwind.c         — contiguous block allocation with shared memory
  source/janitor.c          — block-aware cleanup for shared/private memory
  source/lines.c            — density_override parameter for two_level_atom()
  source/resonate.c         — sobolev() no longer modifies shared state
  source/templates.h        — updated two_level_atom() prototype
  source/communicate_plasma.c — MPI_Barrier(node_comm) after broadcasts
  source/communicate_macro.c  — MPI_Barrier(node_comm) after broadcasts
  source/estimators_macro.c — updated two_level_atom() call
  source/macro_accelerate.c — updated two_level_atom() call
  source/swind_ion.c        — updated two_level_atom() call
  source/wind_updates2d.c   — removed debug log message
  docs/sphinx/source/developer/mpi_comms.rst        — full shared-memory docs
  docs/sphinx/source/developer/programmer_notes.rst  — thread-safety notes

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
Convert 18 fixed-size arrays in plasma_state (spectral model parameters)
and plasma_derived (persistent radiation field averages) from inline
arrays to pointers into shared contiguous blocks, reducing
sizeof(plasma_dummy) by ~2.3 KB/cell.

State arrays moved: f1, f2, spec_mod_type, pl_alpha, pl_log_w,
exp_temp, exp_w, fmin_mod, fmax_mod (4 combined blocks).
Derived arrays moved: F_vis_persistent, F_UV_persistent,
F_Xray_persistent, rad_force_*_persist, F_UV_ang_*_persist.

Also fixes a pre-existing bug in broadcast_plasma_cell() where
spec_mod_type was packed with &cell->state.spec_mod_type (address of
array/pointer) and count 1, instead of cell->state.spec_mod_type with
count NXBANDS. The buffer size calculation is updated to match.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
…OC=500

Remove the compile-time N_PHOT_PROC=500 upper bound and allocate n_bf_in/n_bf_out
as private contiguous blocks sized to nphot_total (e.g. 290 for the standard
macro-atom dataset). This saves ~1680 bytes/cell (~20 MB/rank for 12K cells).

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
Move wmain allocation to MPI_Win_allocate_shared so all ranks on the same
node share a single copy of the wind geometry array (~25 MB savings/rank
for a 300x300 grid).  The reverb path-tracking fields (paths, line_paths)
are moved out of wind_dummy into a separate per-rank wind_paths_main array
since they are written during photon transport.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
cell_spec_flux was a fixed 1000-element array in plasma_estimators,
consuming ~98 MB/rank for the big model (12,271 cells). Convert to a
dynamically allocated pointer with runtime sizing via geo.nbins_in_cell_spec
(default 100, max 1000). Users can restore original behavior with
-cell_spec_dim 1000 on the command line.

geo.cell_freq (frequency bin boundaries) also converted from fixed array
to dynamic allocation in bands_init(). Both are properly freed on exit.

cell_spec_flux uses a private contiguous block (cell_spec_flux_block)
since it is an estimator accumulated during transport.

Backward compatibility: geo.cell_freq pointer nulled after fread in
wind_read(); nbins_in_cell_spec validated after reading old windsave
files. Command-line -cell_spec_dim overrides windsave values on restart.

Saves ~88 MB/rank at default (100 bins) for 12K-cell models.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
The transition probability matrix (matom_matrix, nrows×nrows doubles per
plasma cell) was previously allocated per-cell with separate calloc calls,
leaving each MPI rank holding a full independent copy.  Since the matrix
is written only during wind updates (each rank fills its own slice in
calc_all_matom_matrices) and is strictly read-only during photon transport,
it is a natural fit for MPI-3 shared memory.

Changes:
- calloc_matom_matrix() now allocates a single contiguous shared block
  (macro_block_ptrs.matom_matrix_block) via alloc_block_double with
  use_shared=TRUE, replacing the per-cell allocate_macro_matrix() loop.
  A private per-rank row-pointer array (matom_matrix_rowptrs) preserves
  the double** interface without duplicating the data.
- broadcast_macro_atom_state_matrix() required no changes: it already
  packs/unpacks via matom_matrix[0] which now points into the shared
  block, and the existing MPI_Barrier(node_comm) ensures visibility to
  all node-local ranks.
- janitor.c: replaced per-cell free loop with free_plasma_block on the
  shared block plus free() of the row-pointer array.
- unit_test_model.c: replaced free_and_null(matom_matrix) with a NULL
  assignment since the data is now owned by the block.
- Memory reporting updated: matom_matrix_block included in shared total.
- mpi_comms.rst updated to document the new layout and savings.

For the h20_hetop_standard80 dataset (85 macro levels) on a 300×300 grid
(~12,270 plasma cells) with 24 ranks on one node, this eliminates
23 × 726 MB ≈ 16.7 GB of duplicated physical memory.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
Two low-complexity optimisations in communicate_macro.c:

1. Skip broadcast_macro_atom_state_matrix for single-node runs.
   matom_matrix lives in MPI-3 shared memory, so a MPI_Barrier is
   sufficient when all ranks share the same node (num_nodes==1).
   This avoids the pack/Bcast/unpack cycle and its comm buffer.

2. Replace the monolithic 2x NPLASMA*nlines cooling_bb Allreduce
   (2x585 MB for big.pf) with a chunked MPI_IN_PLACE Allreduce
   targeting ~50 MB per chunk (~13 calls for big.pf). Peak transient
   memory drops from ~1.17 GB to ~50 MB per rank.

Documentation updated in docs/sphinx/source/developer/mpi_comms.rst.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
On macOS with OpenMPI 5.0.9, MPI_Win_allocate_shared maps the shared
window with invalid (read-only) permissions for non-allocating ranks,
causing a SEGV_ACCERR (signal code 2) the first time any non-rank-0
process writes to wmain.  Additionally, kern.sysv.shmmax is only 4 MB
on macOS, far too small for a full wind grid.

Fix: add #ifdef __APPLE__ guards in calloc_wind() (gridwind.c) and
free_wind_grid() (janitor.c) so that macOS always uses private calloc
instead of MPI_Win_allocate_shared.  The Linux shared-memory path is
unchanged.  broadcast_wind_grid() already synchronises wmain across
ranks, so correctness is preserved on macOS.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
Fix all 44 -Wdeprecated-non-prototype warnings on macOS (Apple clang).
All affected functions used old-style K&R parameter declarations; these
have been converted to standard C prototypes. No logic changes.

Also update the local get_models() forward declarations in sirocco.c
and setup.c from empty-parameter style to full prototypes, eliminating
the two remaining "passing arguments without prototype" warnings.

Build is now warning-clean on macOS with Apple clang.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
- communicate_wind.c: replace VLA block_offsets[count] with fixed size
  block_offsets[2]; remove now-unused 'count' variable. Apple clang
  warns on VLAs with const-int size under -Wgnu-folding-constant.

- swind_sub.c, swind_ion.c, rad_hydro_files.c: remove old-style
  FILE *fopen() forward declarations. fopen() is already declared by
  <stdio.h>; the empty-parameter declaration conflicts with the system
  prototype under -Wdeprecated-non-prototype.

Build is now warning-clean on macOS for both make sirocco and make all.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
- Update test_define_wind.c to access plasma fields via ->state sub-struct
  (rho, ne, t_e, t_r, density) following the plasma struct refactor
- Fix cleanup_model to use free_domains/free_wind_grid/free_plasma_grid/
  free_macro_grid instead of per-cell frees, since arrays are now offsets
  into contiguous blocks managed by plasma_block_ptrs
- Fix free() calls on non-heap CYLVAR pointers in unit_test_model.c
- Initialize node_comm/leader_comm in unit_test_main.c to match sirocco.c,
  fixing MPI_ERR_COMM crashes on MPI_Barrier(node_comm)
- Add missing NULL assignments after free() in janitor.c for zdom, wmain,
  plasmamain and macromain, preventing double-free on successive test runs
- Update .pf test parameter files to include new ionization modes
  (LTE_iterate, matrix_multishot) and geometry option (iso) added on
  this branch

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
Resolves conflict in py_progs/run_check.py: memory's version correctly
returns True (skip versioned windsave2table) when the sig file doesn't
name a known code, and adds the missing else-branch for no sig files.
Dev's change (FALSE -> False) was superseded by memory's broader fix.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
Race conditions (Linux MPI multi-rank runs):

- sirocco.c: Move MPI_Finalize() to after clean_on_exit() and Log_close().
  free_wind_grid() calls MPI_Win_free(&wmain_win), which must run before
  MPI_Finalize(); the previous ordering caused SEGV_MAPERR at exit when
  running agn_macro.pf with multiple MPI ranks.

- gridwind.c (calloc_wind): Add MPI_Barrier(node_comm) after the node
  leader's memset of the shared wmain block, so all node-local ranks see
  the zeroed memory before any rank begins writing wind-cell fields.

- communicate_wind.c (broadcast_wind_grid): Add MPI_Barrier(node_comm)
  at the end of the broadcast, matching the pattern used in all plasma
  broadcast functions. Without this, a fast rank could proceed to read
  stale zero-initialised wmain cells (e.g. producing "zero volume but
  flagged inwind" errors in 1d_sn.pf).

- define_wind.c (create_wind_grid): Add two MPI_Barrier(node_comm) calls
  bracketing make_coordinate_grid() and wind_complete(). The first
  prevents a fast rank from writing import-derived inwind values to shared
  wmain while a slow rank's init loop is still writing W_NOT_ASSIGNED to
  the same cells (triggered cv_standard_import.pf failure). The second
  ensures all ranks have finished coordinate grid setup before entering
  the parallel volume/velocity loop.

Documentation:

- docs/sphinx/source/developer/mpi_comms.rst: Document the new barriers
  and add a new section "Platform differences: macOS vs Linux" explaining
  why these races are invisible on macOS (wmain uses private calloc per
  rank via #ifdef __APPLE__) but fatal on Linux (genuinely shared pages).
  Includes a rule of thumb for when Linux testing is required.

Unit test warning fixes:

- source/tests/Makefile: Make -Wno-deprecated-non-prototype Darwin-only
  (Linux GCC does not recognise the flag); add -Wno-unused-result since
  test utility functions intentionally ignore fscanf return values on
  known-good test data files.

- source/tests/tests/test_matrix.c: Initialise matrix/inverse/vector
  pointers and size variables to NULL/0 to silence -Wmaybe-uninitialized.

- source/tests/tests/test_define_wind.c, source/tests/unit_test_model.c:
  Widen path buffers from LINELENGTH to 2*LINELENGTH and update snprintf
  size arguments accordingly, eliminating -Wformat-truncation warnings.

All regression tests pass with mpirun -n 4 and all 11 unit tests pass
with make check.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
Replaces the old MPI communication page with a new page that leads with
an explanation of the v2.0 mixed memory model, contrasting it with the
pure-MPI approach of v1.2. Adds explicit descriptions of both
communication modes (broadcast for state/derived, reduce for
estimators), a decision guide for adding new variables, and coverage of
spectra accumulation via communicate_spectra.c.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
Added node communicator setup (node_comm, leader_comm, node_rank,
node_size, num_nodes) after MPI_Init, matching sirocco.c, so that
MPI_Win_allocate_shared in gridwind.c has a valid communicator.
Also added MPI_Comm_free + MPI_Finalize before all exit() calls.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
…ables docs

- Add -fno-common to all CFLAGS variants in source/Makefile to eliminate
  the ld alignment warning on macOS caused by the large __DATA,__common
  section (~149 MB of uninitialized global arrays).

- Fix windsave.c: on macOS, wmain is allocated as a private calloc per rank
  (not MPI shared memory), so every rank must read wmain from the wind save
  file independently. Guard the shared-memory read path with #ifndef __APPLE__
  to match the allocation strategy in gridwind.c. Without this fix,
  non-rank-0 processes had zeroed wmain, causing dvds_ave=0 everywhere and
  1M+ p_escape warnings that aborted MPI runs.

- Add docs/sphinx/source/c_executables.rst: new page documenting all
  standalone C helper programs (swind, windsave2table, sirocco_optd,
  windsave2fits, rad_hydro_files, modify_wind, inspect_wind), including
  full option descriptions and output file listings for rad_hydro_files.

- Wire c_executables into index.rst toctree and update output/model.rst
  to cross-reference the new page.

- Fix RST errors in mixed_memory_model.rst: replace ^ title underlines
  with - so subsections of the MPI-3 shared memory section are correctly
  at level 3 rather than skipping from level 2 to level 4.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
* Fix MPI_Finalize abort on MPICH by explicitly freeing shared-memory windows and communicators

Shared-memory MPI_Win handles for plasma and macro blocks were not freed
before MPI_Finalize, and node_comm/leader_comm were never freed at all.
OpenMPI tolerates this; MPICH OFI transport aborts with OFI poll failed
when it encounters unreleased resources during send-queue flushing.

Fix: explicitly call MPI_Win_free on all plasma (15) and macro (6+1)
shared windows in free_plasma_grid/free_macro_grid, add MPI_Comm_free
for node_comm and leader_comm, and add MPI_Barrier(MPI_COMM_WORLD)
before clean_on_exit to synchronize all ranks before releasing shared
resources. Applied to all three MPI_Finalize sites: normal exit,
grid-only early exit, and the max-time timeout path in signal.c.


* Fix unit_test.c MPI cleanup: add node/leader comm init and MPI_Finalize

unit_test.c called MPI_Init but never set up node_comm or leader_comm,
and returned without calling MPI_Finalize. On MPICH this produces the
same OFI abort seen in sirocco. Fix mirrors unit_test_main.c and
sirocco.c: add MPI_Comm_split_type/MPI_Comm_split for node_comm and
leader_comm after MPI_Init, MPI_Barrier + MPI_Comm_free + MPI_Finalize
before return. Also replace exit(1) in zparse with Exit(1) so that if
LUM_TEST is ever re-enabled, argument errors go through MPI_Abort
rather than silently exiting one rank.
* Fix MPI_Finalize abort on MPICH by explicitly freeing shared-memory windows and communicators

Shared-memory MPI_Win handles for plasma and macro blocks were not freed
before MPI_Finalize, and node_comm/leader_comm were never freed at all.
OpenMPI tolerates this; MPICH OFI transport aborts with OFI poll failed
when it encounters unreleased resources during send-queue flushing.

Fix: explicitly call MPI_Win_free on all plasma (15) and macro (6+1)
shared windows in free_plasma_grid/free_macro_grid, add MPI_Comm_free
for node_comm and leader_comm, and add MPI_Barrier(MPI_COMM_WORLD)
before clean_on_exit to synchronize all ranks before releasing shared
resources. Applied to all three MPI_Finalize sites: normal exit,
grid-only early exit, and the max-time timeout path in signal.c.



* Fix unit_test.c MPI cleanup: add node/leader comm init and MPI_Finalize

unit_test.c called MPI_Init but never set up node_comm or leader_comm,
and returned without calling MPI_Finalize. On MPICH this produces the
same OFI abort seen in sirocco. Fix mirrors unit_test_main.c and
sirocco.c: add MPI_Comm_split_type/MPI_Comm_split for node_comm and
leader_comm after MPI_Init, MPI_Barrier + MPI_Comm_free + MPI_Finalize
before return. Also replace exit(1) in zparse with Exit(1) so that if
LUM_TEST is ever re-enabled, argument errors go through MPI_Abort
rather than silently exiting one rank.



* rad_hydro_files: add t_r to hc output and output fill-corrected density to pcon

Output t_r alongside t_e in the heating/cooling file headers and data rows.
In the pcon output, multiply rho and ne by the domain fill factor to recover
the cell-averaged density originally supplied by the hydro code (Sirocco
stores rho/fill internally). Simplify old_ne to ne*fill directly.
When a photon lands at exactly r = rmax (the outer boundary of a spherical
cell) with an outgoing or tangential direction, ds_to_sphere returns VERY_BIG
for both the inner sphere (impact parameter > rmin for a tangential path) and
the outer sphere (photon is on or outside it, moving away).  The previous fix
returned DFUDGE, but a tangential nudge of DFUDGE does not change r in
floating-point arithmetic — sqrt(rmax^2 + DFUDGE^2) == rmax at double
precision for rmax ~ 1e11 cm and DFUDGE ~ 0.01 cm — so the photon stays
trapped at r = rmax and the error fires on every subsequent call to
calculate_ds.  Each iteration logs a dfreq=0 Error in resonate.c; in dense,
optically thick models the accumulated count reaches the max_errors = 1e6
abort threshold and kills the run.

Fix: return ds_to_sphere(wmain[n].rmax + DFUDGE, p) instead of DFUDGE.
For any photon classified in cell n (r <= rmax), the expanded sphere at
rmax + DFUDGE always yields a finite positive root, and the move places the
photon at r = rmax + DFUDGE so that where_in_grid correctly assigns it to the
next cell on the following call.  This reduces the VERY_BIG error count from
~900k to ~800 and the dfreq=0 count from 1,000,001 (abort) to ~16 for a dense
T_init=1000K shell-wind test case run over 20 ionization cycles.

The bug cannot trigger at the inner cell boundary because the photon is always
inside the outer sphere there, guaranteeing ds_to_sphere(rmax) has a finite
forward root.  For multi-cell winds the problem is self-terminating (the
photon exits to the next cell), but the outermost wind cell is vulnerable in
the same way as the shell-wind case here.

Note: the same bug exists in any branch derived from main that has not
received this fix.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
Bug 1 — dvwind_ds_cmf: finite-difference step below floating-point precision

The step ds = 1e-6 * half_cell_width is used to evaluate the velocity
gradient numerically.  For thin shells at large radii (e.g. the 10 cm
test shell at r = 1e11 cm), ds = 5e-6 cm falls below ULP(r) ≈ 2.2e-5 cm,
so rmin + ds == rmin in double precision.  The test position does not move,
v2 == v1 for every direction tried, and dvds ≈ 0 for all 10 000 trials in
randwind_thermal_trapping.  The rejection loop never accepts a direction,
exhausts NSCAT_MAX, and logs an Error on every call.  In models with many
wind photons near the inner boundary the error count reaches max_errors
(1 000 000) and aborts.

Fix: clamp ds to max(1e-6 * half_cell, min(100*DBL_EPSILON*r, 0.1*half_cell)).
This guarantees the step is representable while remaining local to the cell.

Bug 2 — get_dvds_max: missing frac[] weights in interpolation

coord_fraction returns corner indices nnn[] and bilinear weights frac[]
that sum to 1.  Every other interpolation in the codebase writes
  value += frac[nn] * quantity[nnn[nn]]
but get_dvds_max accumulated the corner dvds_max values without the frac
weights, returning a sum instead of a weighted average.  For a photon
midway between two wind cells each with dvds_max ≈ 0.9 the function
returned 1.8, over-estimating p_norm in randwind_thermal_trapping and
making the acceptance criterion unnecessarily strict.

Fix: add frac[nn] * to the accumulation, consistent with all other uses.

Also removed a dead vsub() call (line 88 in the original) that was
immediately overwritten by the following vsub().

Note: both bugs exist on every branch derived from main that has not
received this fix.  Cherry-pick commit to port.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
Two changes to the extra-diagnostics file (files.extra / *.ext.txt):

1. parse.c: build the path inside files.diagfolder so ext.txt files
   appear alongside the .diag files rather than cluttering the working
   directory.

2. diag.c: guard the fopen() in init_extra_diagnostics on
   (modes.save_extract_photons || modes.track_resonant_scatters) instead
   of the broader modes.extra_diagnostics flag.  Previously the file was
   created (one per MPI rank) whenever any @diag option was active, even
   though only those two modes ever write to epltptr.  The files were
   always empty unless one of those two modes was explicitly enabled.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
The fix cherry-picked from x3d used wmain[n].rmax, which does not exist
in the memory branch wind_dummy struct.  On this branch the outer cell
boundary is zdom[ndom].wind_x[ix + 1], consistent with the existing
ds_to_sphere calls above it.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
…rmat

Reads a Sirocco macro-atom master data file and produces a simple-atom
equivalent by converting:
  IonM      -> IonV        (elem_ions files; nlte preserved for level tracking)
  LevMacro  -> LevTop      (levels files; islp/qqnum set to -1 as placeholders)
  LinMacro  -> Line        (lines files; filtered by lower level index and
                            oscillator strength to keep only the most
                            important transitions, default ll <= 1)
  PhotMacS  -> PhotTopS    (phot files; levl used as ilv, islp=-1 to match
  PhotMac   -> PhotTop      converted LevTop entries)

Files that need no conversion (topbase phot, collision data, atomic/ files)
are referenced in-place in the new master file.  All output is written
relative to the input master file location so test runs stay out of the
main xdata directory.

Usage:
  macro2simple.py data/h20_hetop_standard80.dat [--max-lower-level N]
                  [--min-osc-strength F] [--outdir DIR]

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
kslong and others added 6 commits June 17, 2026 10:53
macOS case-insensitive filesystems collide macro2simple.py with Macro2Simple.py.
ConvertMacro2Simple.py provides the same functionality without the name collision.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
macro2simple.py collides with Macro2Simple.py on macOS case-insensitive
filesystems. Functionality is preserved in ConvertMacro2Simple.py.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
…compat

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
Matches x3d branch content: build commands, architecture overview, coding
conventions, and both Claude Code rules (version number + macOS filenames).

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
Two related bugs introduced when cell_spec_flux was converted from a
fixed array to a dynamic pointer (9362f1f):

1. windsave.c: cell_spec_flux was never written to or read from the
   wind_save file, leaving the Jnu extension in windsave2fits output
   all zeros. Add fwrite/fread in the per-cell dynamic-array block.

2. sirocco.c: -cell_spec_dim N was silently ignored on new runs because
   init_geo() resets geo.nbins_in_cell_spec to 100 after parse_command_line
   captures the flag. The restart and previous-run paths already
   re-applied cmd_nbins_in_cell_spec after wind_read; add the same
   restore+validate block for RUN_TYPE_NEW.

Note: the windsave.c change alters the wind_save binary format; files
written before this fix cannot be read correctly by the updated code.

memory, polar, and biconic branches have both bugs and need
cherry-pick of this commit.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
# Conflicts:
#	source/communicate_macro.c
#	source/communicate_plasma.c
#	source/cooling.c
#	source/emission.c
#	source/estimators_macro.c
#	source/ionization.c
#	source/sirocco.h
#	source/test_cooling.c
#	source/wind_updates2d.c
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants