Refactor opacity and 1D RT#128
Conversation
Introduce dump-backed MoleculeLine and MoleculeCIA opacity modules and wire them into the radiation factory and Python bindings. The new readers resolve NetCDF dimension order from variable metadata, preserve the existing runtime units, add same-species continuum terms to molecular line opacity, and read CIA binary coefficients from the same dump format. Remove the legacy RFM-facing surface from this repo, including the old RFM sources, helper scripts, bindings, query helpers, and AMARS example tree. Update docs, stubs, sample YAML, and tests to use molecule-line and molecule-cia, and enable the replacement attenuation unit test with synthetic NetCDF coverage.
Move the NetCDF opacity helper functions from src/opacity into src/utils and flatten them into the harp namespace so both molecule readers consume the same utility surface through harp/utils/netcdf_opacity_utils.hpp. Update MoleculeLine and MoleculeCIA to store logarithmic opacity tables and logarithmic base temperatures. Both modules now interpolate ln(cross section) and ln(base temperature), exponentiate after interpolation, and keep the runtime attenuation calculation unchanged. This also adds positivity checks before taking logarithms.
Refine the NetCDF-backed molecule opacity path so sparse tables with zero entries no longer fail log-space interpolation. Centralize the positive fill helper in the shared NetCDF opacity utilities and update the attenuation tests to exercise zero-valued line and CIA cells explicitly. Add RadiationBandOptions::set_wave_lower and set_wave_upper so spectral bin edges can be configured at the band level and forwarded consistently to the active solver. Update the YAML-driven radiation-band setup and both C++ and Python example call sites to use the new API. Rewrite the molecule dump example into a one-band column-radiation workflow. The new example YAML defines species, geometry, molecular line and CIA opacity sources, and a longwave DISORT band; the C++ driver now constructs a TP column, computes per-source and total attenuation, runs RadiationBand, and writes flux/attenuation diagnostics.
Update the molecule dump column-radiation example so its state tensors match the native RadiationBand input convention directly. The example now stores conc as (ncol, nlyr, nspecies) and temp/pres as (ncol, nlyr), removes the temporary batching helpers, and passes the revised tensors through both attenuation and flux calculations.
Rename the molecule dump example source and YAML to example_jupiter_1d_rt.cpp and jupiter_1d_rt.yaml, update the examples CMake target accordingly, and keep the renamed example wired to the current Jupiter 1D RT configuration. Add --config and --output-dir command-line options to the example driver, defaulting to jupiter_1d_rt.yaml and the current directory. The example now adds the selected config directory to the resource search path and writes diagnostics under the requested output root.
Move mean molecular weight support from opacity to utils, expose the helper through the Python API, and reuse it in the Jupiter 1D RT example. Also move composition normalization and formatting helpers into compound.hpp so composition handling is shared instead of duplicated in the example. Rename the example source to jupiter_1d_rt.cpp, update the example target wiring, and replace the ad hoc output writers with reusable text-format utilities for column and spectral profiles. The generated profile tables now preserve full field names, include field/unit comments, and emit the requested flux/transmittance diagnostics.
Update the spectroscopy CLI stack to use endpoint-inclusive wavenumber ranges and align the user-facing documentation with that behavior. SpectralBandConfig now includes both the start and end wavenumber on generated grids, and the dump/plot docs and README examples now describe inclusive bounds and shared-boundary duplication for adjacent ranges. Replace the previous monolithic CLI module layout with split implementations: utils.py now owns shared CLI helpers and output naming, hitran_cia_utils.py and hitran_cia_plot.py separate CIA data logic from CIA plotting/CLI entry points, and hitran_molecule_utils.py and hitran_molecule_plot.py do the same for molecular HITRAN workflows. Remove the legacy atm_overview_cli.py, cia_plot_cli.py, molecule_plot_cli.py, shared_cli.py, hitran_lines.py, hitran_cia.py, hitran_molecule.py, and output_names.py modules in favor of the new layout. Normalize public helper naming where functions are used across module boundaries or tests. Shared helpers such as parse_composition, compute_requested_absorption_spectrum, compute_requested_transmittance, compute_overview_products, load_requested_cia_dataset, and state/selection helpers now use non-underscored names, while file-local implementation details remain private. Also update imports and tests throughout the repo to target the new modules directly, and keep the NetCDF/plot naming helpers consolidated under utils.py.
Restore working spectra CLI/runtime behavior after the refactor series. This commit fixes several issues that showed up when running pyharp-dump and pyharp-plot from an editable checkout on restricted cluster hosts: - replace the staged NetCDF temp-file write-and-move flow with direct writes to the final path - force HDF5_USE_FILE_LOCKING=FALSE during dataset writes to avoid NFS write stalls - rename the shared dataset writer to write_dataset and update dump/spectrum/ transmittance call sites and tests - break the circular import between spectrum.py and hitran_molecule_utils.py by switching the shared spectrum/transmittance helpers there to lazy imports - keep dump and plot parallel execution working when ProcessPoolExecutor cannot allocate SemLock resources by falling back to subprocess-based workers rather than crashing or degrading to serial execution - apply the subprocess fallback consistently across dump_cli.py, plot_cli.py, atm_overview.py, and hitran_molecule_plot.py - add regression coverage for the new file-locking behavior and the restricted- host parallel fallback paths - fix the atm_overview subprocess fallback signature bug discovered while exercising pyharp-plot overview These changes are limited to the Python spectra tooling and its tests; the unrelated python/csrc/pyopacity.cpp worktree change is intentionally excluded from this commit.
There was a problem hiding this comment.
Pull request overview
Refactors the spectroscopy + opacity toolchain to remove the legacy RFM workflow, switch wavenumber grids to endpoint-inclusive semantics, and introduce NetCDF dump-backed molecular opacities for 1D RT and related CLIs.
Changes:
- Update
SpectralBandConfig.grid()to include both endpoints and propagate the new range semantics across CLIs, docs, and tests. - Replace legacy RFM-based opacity handling with NetCDF dump-backed
molecule-lineandmolecule-ciaopacities (C++ core + Python bindings/stubs). - Add multiprocessing fallbacks (subprocess-based) and HDF5 file-locking mitigation for CLI workflows.
Reviewed changes
Copilot reviewed 85 out of 88 changed files in this pull request and generated 6 comments.
Show a summary per file
| File | Description |
|---|---|
| tests/spectra/test_spectrum.py | Update imports to new HITRAN CIA utils module. |
| tests/spectra/test_shared_cli.py | Point process-pool tests at spectra.utils entry points. |
| tests/spectra/test_plot_cli.py | Update plot CLI dispatch targets; adjust wn-range expectations; add fallback test. |
| tests/spectra/test_output_names.py | Move default output naming to spectra.utils. |
| tests/spectra/test_molecule_plot_cli.py | Update tests to new HITRAN molecule plot/utils split and inclusive grid. |
| tests/spectra/test_hitran_lines.py | Redirect HITRAN line helpers/tests to new molecule utils/plot modules. |
| tests/spectra/test_config.py | Add/adjust test for endpoint-inclusive spectral grids. |
| tests/spectra/test_cli.py | Update dump CLI tests for new writer + inclusive grid + fallback. |
| tests/spectra/test_cia.py | Split CIA plotting vs CIA utility functions; adjust imports/patches. |
| tests/spectra/test_atm_overview_cli.py | Move overview logic into atm_overview module; adjust parsing and fallback coverage. |
| tests/CMakeLists.txt | Re-enable test_attenuator. |
| src/utils/write_spectral_profile.hpp | New helper to write formatted spectral profiles. |
| src/utils/write_column_profile.hpp | New helper to write formatted 1D column profiles. |
| src/utils/netcdf_opacity_utils.hpp | New NetCDF reading + unit conversion + dim permutation utilities for opacities. |
| src/utils/mean_molecular_weight.hpp | Add reusable mean molecular weight utility (now shared + bound to Python). |
| src/radiation/radiation_band.hpp | Add solver wave-bound setters to options. |
| src/radiation/radiation_band.cpp | Wire new molecule-* opacities and use wave setter helpers. |
| src/opacity/rfm.cpp | Remove legacy RFM opacity implementation. |
| src/opacity/read_rayleigh.hpp | Remove unused legacy reader. |
| src/opacity/read_rayleigh.cpp | Remove unused legacy reader implementation. |
| src/opacity/read_cia_ff.hpp | Remove legacy CIA/freefree readers. |
| src/opacity/read_cia_ff.cpp | Remove legacy CIA/freefree readers implementation. |
| src/opacity/oxygen_cia.hpp | Remove legacy oxygen CIA implementation. |
| src/opacity/oxygen_cia.cpp_ | Remove legacy oxygen CIA implementation. |
| src/opacity/opacity_options.hpp | Remove RFM-centric query helpers from options header. |
| src/opacity/opacity_options.cpp | Remove RFM/multiband query helpers implementation. |
| src/opacity/nitrogen_cia.hpp | Remove legacy nitrogen CIA implementation. |
| src/opacity/nitrogen_cia.cpp_ | Remove legacy nitrogen CIA implementation. |
| src/opacity/molecule_line.hpp | Rename/repurpose interface for NetCDF dump-backed line opacity. |
| src/opacity/molecule_line.cpp | Implement NetCDF dump-backed molecular line opacity reader/interpolator. |
| src/opacity/molecule_cia.hpp | New NetCDF dump-backed CIA opacity module interface. |
| src/opacity/molecule_cia.cpp | Implement NetCDF dump-backed CIA opacity reader/interpolator. |
| src/opacity/mean_molecular_weight.hpp | Remove old opacity-local mean molecular weight utility. |
| src/compound.hpp | Add composition normalization and formatting helpers. |
| python/spectra/utils.py | Consolidate CLI utilities (naming, ranges, bands/grids) and defaults. |
| python/spectra/transmittance.py | Switch dataset writing helper to write_dataset. |
| python/spectra/spectrum.py | Switch to new HITRAN utils + new dataset writer. |
| python/spectra/plot_cli.py | Refactor dispatch targets; add subprocess fallback for parallelism; update defaults. |
| python/spectra/output_names.py | Remove module (moved into spectra.utils). |
| python/spectra/hitran_molecule_utils.py | Expand molecule workflow utilities (band/config build, CIA selection, spectrum/transmittance computation). |
| python/spectra/hitran_cia_utils.py | Keep CIA I/O + numerics; remove plotting functions. |
| python/spectra/hitran_cia_plot.py | Move CIA plotting + CLI helpers into dedicated module. |
| python/spectra/dump_cli.py | Add subprocess fallback for parallelism; switch writer; align composition parsing. |
| python/spectra/dataset_io.py | Replace tmp-write helper with env-based HDF5 locking mitigation + direct writer. |
| python/spectra/config.py | Change grid semantics to endpoint-inclusive. |
| python/spectra/atm_overview.py | Move overview CLI logic into module; add subprocess fallback; align parsing/range handling. |
| python/run_ktable_amars.py | Remove legacy script tied to RFM workflow. |
| python/run_cktable.py | Remove legacy script tied to RFM workflow. |
| python/rfmlib.py | Remove legacy RFM Python library. |
| python/pyharp.pyi | Add stubs for new bound utilities and band option setters. |
| python/opacity.pyi | Update opacity types/classes to MoleculeLine/MoleculeCIA. |
| python/klib.py | Remove legacy correlated-k helper module. |
| python/csrc/pyradiation.cpp | Bind new set_wave_lower/upper methods. |
| python/csrc/pyopacity.cpp | Bind MoleculeLine and MoleculeCIA; remove RFM binding. |
| python/csrc/pyharp.cpp | Bind new utility functions (mmw + profile writers). |
| python/init.py | Stop exporting removed rfmlib. |
| examples/test2/grey-opacity2.pt | Example artifact retained. |
| examples/test2/grey-opacity.pt | Example artifact retained. |
| examples/test2/example_jit.py | Remove legacy example script. |
| examples/test2/amars-ck.yaml | Remove legacy example config. |
| examples/jupiter_1d_rt.yaml | Add new 1D RT example using dump-backed opacities. |
| examples/example_sonora_2020_flux.py | Update to new band wave-bound setters. |
| examples/CMakeLists.txt | Replace legacy AMARS example wiring with Jupiter 1D RT example. |
| examples/2025-amars/rfm.atm | Remove legacy AMARS artifacts. |
| examples/2025-amars/pVals.txt | Remove legacy AMARS artifacts. |
| examples/2025-amars/amarsw-lbl.yaml | Remove legacy AMARS artifacts. |
| examples/2025-amars/amarsw-ck.yaml | Remove legacy AMARS artifacts. |
| examples/2025-amars/amars_sw.py | Remove legacy AMARS artifacts. |
| examples/2025-amars/amars_rt_rk3.cpp | Remove legacy AMARS artifacts. |
| examples/2025-amars/amars_rt_cli.cpp | Remove legacy AMARS artifacts. |
| examples/2025-amars/amars_rt.py | Remove legacy AMARS artifacts. |
| examples/2025-amars/amars_rt.cpp | Remove legacy AMARS artifacts. |
| examples/2025-amars/amars_lw.cpp | Remove legacy AMARS artifacts. |
| examples/2025-amars/amars-ck.yaml | Remove legacy AMARS artifacts. |
| examples/2025-amars/TVals.txt | Remove legacy AMARS artifacts. |
| examples/2025-amars/CMakeLists.txt | Remove legacy AMARS example build wiring. |
| docs/source/plot_cli.rst | Document endpoint-inclusive range semantics. |
| docs/source/opacity_classes.rst | Update documented opacity classes. |
| docs/source/opacity.rst | Update opacity type list and NetCDF schema expectations. |
| docs/source/dump_cli.rst | Document endpoint-inclusive range semantics and overlap behavior. |
| docs/source/builtin_opacities.rst | Add docs for dump-backed molecular opacities. |
| README.md | Update range semantics + CLI behavior descriptions. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| max_workers = min(len(tasks), os.cpu_count() or 1) | ||
| ctx = process_pool_context() | ||
| with ProcessPoolExecutor(max_workers=max_workers, mp_context=ctx) as executor: | ||
| return list(executor.map(worker, tasks)) | ||
| try: | ||
| with ProcessPoolExecutor(max_workers=max_workers, mp_context=ctx) as executor: | ||
| return list(executor.map(worker, tasks)) | ||
| except PermissionError: | ||
| return _parallel_plot_results_via_subprocess(tasks, worker=worker, max_workers=max_workers) |
There was a problem hiding this comment.
_parallel_plot_results() catches PermissionError for the entire ProcessPoolExecutor block, which can also swallow a real PermissionError raised by the worker itself (e.g., reading/writing files) and incorrectly trigger the subprocess fallback. Consider narrowing the try/except to only the executor construction / mp_context setup, or re-raising when the exception originates from task execution.
| def _run_task_via_subprocess(task: tuple[str, argparse.Namespace], *, worker_name: str) -> None: | ||
| code = ( | ||
| "from pyharp.spectra.plot_cli import _subprocess_worker_entry; " | ||
| "import sys; " | ||
| "_subprocess_worker_entry(sys.argv[1], sys.argv[2], sys.argv[3])" | ||
| ) | ||
| with tempfile.TemporaryDirectory(prefix="pyharp_plot_task_") as tmpdir: | ||
| task_path = Path(tmpdir) / "task.pkl" | ||
| result_path = Path(tmpdir) / "result.pkl" | ||
| with task_path.open("wb") as handle: | ||
| pickle.dump(task, handle, protocol=pickle.HIGHEST_PROTOCOL) | ||
| subprocess.run( | ||
| [sys.executable, "-c", code, str(task_path), str(result_path), worker_name], | ||
| check=True, | ||
| ) | ||
| with result_path.open("rb") as handle: | ||
| return pickle.load(handle) |
There was a problem hiding this comment.
_run_task_via_subprocess() is annotated as returning None, but it returns the unpickled result from result.pkl. This mismatch will confuse type checkers and readers; either adjust the return annotation to match what’s returned (and align the caller’s type hints), or don’t return a value and keep the function None-returning.
| def write_dataset(dataset: xr.Dataset, output_path: Path, *, engine: str = DEFAULT_NETCDF_ENGINE) -> None: | ||
| output_path.parent.mkdir(parents=True, exist_ok=True) | ||
| with hdf5_file_locking_disabled(): | ||
| dataset.to_netcdf(output_path, engine=engine) |
There was a problem hiding this comment.
write_dataset() now writes directly to the final output_path. If the process crashes or is interrupted mid-write, this can leave a partially-written/corrupt NetCDF file at the destination. Consider restoring an atomic write pattern (write to a temp file in the same directory, then replace/rename), while still using hdf5_file_locking_disabled() around the to_netcdf call.
| std::filesystem::create_directories(filename.parent_path()); | ||
| std::ofstream out(filename); | ||
| TORCH_CHECK(out, "Failed to open output file: ", filename.string()); |
There was a problem hiding this comment.
std::filesystem::create_directories(filename.parent_path()) will throw if filename has no parent directory component (e.g., "profile.txt"), because parent_path() can be empty. Guard this with a check for an empty parent path before creating directories.
| std::filesystem::create_directories(filename.parent_path()); | ||
| std::ofstream out(filename); | ||
| TORCH_CHECK(out, "Failed to open output file: ", filename.string()); |
There was a problem hiding this comment.
std::filesystem::create_directories(filename.parent_path()) can throw when filename has no directory component (empty parent_path()). Add a guard for an empty parent path before attempting to create directories so callers can pass plain filenames safely.
| inline torch::Tensor mean_molecular_weight(torch::Tensor conc) { | ||
| TORCH_CHECK(species_weights.size() == static_cast<size_t>(conc.size(-1)), | ||
| "The last dimension of 'conc' must match the number of species"); | ||
|
|
||
| auto ww = torch::tensor( | ||
| species_weights, | ||
| torch::TensorOptions().dtype(conc.dtype()).device(conc.device())); | ||
|
|
||
| int ndim = conc.dim(); | ||
| std::vector<int64_t> shape(ndim, 1); | ||
| shape.back() = static_cast<int64_t>(species_weights.size()); | ||
|
|
||
| return (conc * ww.view(shape)).sum(-1) / conc.sum(-1); | ||
| } |
There was a problem hiding this comment.
mean_molecular_weight() divides by conc.sum(-1) without guarding against the all-zero case. If a layer/column has zero total concentration, this will produce NaN/Inf and can silently propagate into outputs (e.g., write_column_profile). Consider validating the denominator is positive (or defining a fill value) before division.
No description provided.