diff --git a/.github/workflows/nightly.yaml b/.github/workflows/nightly.yaml index fa0a34a..f23babe 100644 --- a/.github/workflows/nightly.yaml +++ b/.github/workflows/nightly.yaml @@ -49,6 +49,7 @@ jobs: cd ROM export PYTHONPATH=$PWD/python cd examples/reed && python3 run_rom_reed.py --exe ../../build/rom_app_exec + cd ../1dk && python3 run_rom_P58.py --exe ../../build/rom_app_exec - name: compile diffusion shell: bash run: | diff --git a/ROM/README.md b/ROM/README.md new file mode 100644 index 0000000..54fc309 --- /dev/null +++ b/ROM/README.md @@ -0,0 +1,43 @@ +Standalone ROM application built around OpenSn and libROM, with supporting Python tooling and example input decks. + +This repo is organized to: +- build a C++ application/binary that links against OpenSn and libROM, +- run transport problems and ROM workflows from Python drivers, +- keep reproducible example cases in a single place. + +## Repository layout + +- [src/](src/README.md) — C++ sources for the application (python bindings, app-specific code) +- [python/](python/README.md) — Python drivers/utilities (e.g., job management, ROM pipeline orchestration) +- [examples/](examples/README.md) — Example decks, runs, and reference cases + +# Quickstart + +## Dependencies + +This repository builds a standalone application that links against a libROM installation, an OpenSn installation and certain OpenSn dependencies. + +You must have the following available: + +- OpenSn +- libROM +- OpenSn dependencies (VTK, Caliper, MPI, Python with pybind11) + +--- + +## Build + +From the repository root: + +```bash +mkdir build +cd build +cmake .. \ + -DOpenSn_DIR= \ + -DlibROM_DIR= +make -j +``` +An example can be run from its respective folder with +```bash +python run_rom_*.py --exe=path/to/app/exe +``` \ No newline at end of file diff --git a/ROM/examples/1dk/H2O_mg_base.txt b/ROM/examples/1dk/H2O_mg_base.txt new file mode 100644 index 0000000..e5274c3 --- /dev/null +++ b/ROM/examples/1dk/H2O_mg_base.txt @@ -0,0 +1,14 @@ +NUM_GROUPS 2 +NUM_MOMENTS 1 + +SIGMA_T_BEGIN +0 0.88798 +1 2.9865 +SIGMA_T_END + +TRANSFER_MOMENTS_BEGIN +M_GFROM_GTO_VAL 0 0 0 0.83975 +M_GFROM_GTO_VAL 0 0 1 0.04749 +M_GFROM_GTO_VAL 0 1 0 0.000336 +M_GFROM_GTO_VAL 0 1 1 2.9676 +TRANSFER_MOMENTS_END \ No newline at end of file diff --git a/ROM/examples/1dk/P58_problem.py b/ROM/examples/1dk/P58_problem.py new file mode 100644 index 0000000..c3db200 --- /dev/null +++ b/ROM/examples/1dk/P58_problem.py @@ -0,0 +1,100 @@ +from pathlib import Path +import utils +import plotting +import numpy as np +import xs + + +class P58Problem: + def __init__(self, workdir, nprocs=2, ntrain=100, ntest=10): + self.workdir = Path(workdir) + self.deck_path = self.workdir / "base_P58.py" + + self.nprocs = nprocs + self.ntrain = ntrain + self.ntest = ntest + + self.xs = xs.CrossSections( + [ + { + "in_file": "URRb_base.txt", + "out_file": "data/URRb.xs", + "kind": "fission", + }, + ], + frac=0.2, + transfer_tol=1.0e-14, + param_mode="entrywise", + ) + + # Preserve the original curated P58 parameter domain rather than + # adopting auto-generated relative bounds from xs.py. + self.bounds = [ + [0.000836,0.001648], # sigma_f[0] + [0.029564,0.057296], # sigma_f[1] + [0.001104,0.001472], # sigma_c[0] + [0.024069,0.029244], # sigma_c[1] + [0.83807,0.83892], # S[0,0] + [0.04536,0.04635], # S[0,1] + [0.000767,0.00116], # S[1,0] + [2.8751,2.9183], # S[1,1] + ] + if len(self.bounds) != self.xs.n_params: + raise ValueError( + "P58 bounds define {} parameters, but xs.py expects {}." + .format(len(self.bounds), self.xs.n_params) + ) + + def sample_training(self): + self.training_set = utils.sample_LHS(self.bounds, self.ntrain) + + params_path = self.workdir / "data" / "params.txt" + np.savetxt(str(params_path), self.training_set) + # Also move water xs file to data + with open("H2O_mg_base.txt", "rb") as f_src, open("data/H2O_mg.xs", "wb") as f_dst: + f_dst.write(f_src.read()) + + def load_training(self): + params_path = self.workdir / "data" / "params.txt" + self.training_set = np.loadtxt(str(params_path)) + + def sample_testing(self): + self.testing_set = utils.sample_test(self.bounds, self.ntest) + + params_path = self.workdir / "data" / "test_params.txt" + np.savetxt(str(params_path), self.testing_set) + + def update_xs(self, pvec): + self.xs.write_sample(pvec) + + def plot_results(self): + plotting.plot_sv(num_groups=2) + errors = [] + k_errors = [] + speedups = [] + for i in range(self.ntest): + results_dir = self.workdir / "results" + rom_time = np.loadtxt(str(results_dir / "online_time_{}.txt".format(i))) + fom_time = np.loadtxt(str(results_dir / "offline_time_{}.txt".format(i))) + + output_dir = self.workdir / "output" + error = plotting.plot_1d_eigenvector( + str(output_dir / ("fom_{}_".format(i) + "{}.h5")), + str(output_dir / ("rom_{}_".format(i) + "{}.h5")), + ranks=range(self.nprocs), + pid=i) + k_error = np.abs( + np.loadtxt("output/fom_k_{}.txt".format(i)) + - np.loadtxt("output/rom_k_{}.txt".format(i)) + ) + + k_errors.append(k_error) + errors.append(error) + speedups.append(fom_time / rom_time) + + print("Avg Eigenvector Error ", np.mean(errors)) + np.savetxt("results/errors.txt", errors) + print("Avg k Error ", np.mean(k_errors) * 1e5, "pcm") + np.savetxt("results/k_errors.txt", k_errors) + print("Avg Speedup ", np.mean(speedups)) + np.savetxt("results/speedups.txt", speedups) diff --git a/ROM/examples/1dk/README.md b/ROM/examples/1dk/README.md new file mode 100644 index 0000000..d69ec45 --- /dev/null +++ b/ROM/examples/1dk/README.md @@ -0,0 +1,27 @@ +# 1dk/ + +**Location:** [repo root](../../README.md) / [examples/](../README.md) / `1dk/` + +1-D k-eigenvalue example and ROM driver. Based off of P58 from [LA-13511](https://doi.org/10.2172/10601) + +--- + +## Files + +- `base_P58.py` + Base OpenSn deck for problem 58. + +- `P58_problem.py` + Problem definition and parameterization for problem 58. + +- `run_rom_P58.py` + Entry-point script that runs the ROM workflow. + +--- + +## How to Run + +From this directory: + +```bash +python run_rom_P58.py --exe=path/to/app/exe diff --git a/ROM/examples/1dk/URRb_base.txt b/ROM/examples/1dk/URRb_base.txt new file mode 100644 index 0000000..54dbe37 --- /dev/null +++ b/ROM/examples/1dk/URRb_base.txt @@ -0,0 +1,29 @@ +NUM_GROUPS 2 +NUM_MOMENTS 1 + +SIGMA_T_BEGIN +0 0.88721 +1 2.9727 +SIGMA_T_END + +SIGMA_F_BEGIN +0 0.000836 +1 0.029564 +SIGMA_F_END + +NU_BEGIN +0 2.50 +1 2.50 +NU_END + +CHI_BEGIN +0 1.0 +1 0.0 +CHI_END + +TRANSFER_MOMENTS_BEGIN +M_GFROM_GTO_VAL 0 0 0 0.83892 +M_GFROM_GTO_VAL 0 0 1 0.04635 +M_GFROM_GTO_VAL 0 1 0 0.000767 +M_GFROM_GTO_VAL 0 1 1 2.9183 +TRANSFER_MOMENTS_END \ No newline at end of file diff --git a/ROM/examples/1dk/base_P58.py b/ROM/examples/1dk/base_P58.py new file mode 100644 index 0000000..7929a24 --- /dev/null +++ b/ROM/examples/1dk/base_P58.py @@ -0,0 +1,106 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +import numpy as np + +if __name__ == "__main__": + + try: + print("Parameter id = {}".format(pid)) + except: + p_id=0 + print("Parameter id = {}".format(pid)) + + print("{} phase".format(phase)) + + widths = [4.6, 1.126152] + nrefs = [50, 50] + Nmat = len(widths) + nodes = [0.] + for imat in range(Nmat): + dx = widths[imat] / nrefs[imat] + for i in range(nrefs[imat]): + nodes.append(nodes[-1] + dx) + meshgen = OrthogonalMeshGenerator(node_sets=[nodes]) + grid = meshgen.Execute() + grid.SetUniformBlockID(0) + + # Set block IDs + lv = RPPLogicalVolume(infx=True, infy=True, zmin=0.0, zmax=4.6) + grid.SetBlockIDFromLogicalVolume(lv, 1, True) + + num_groups = 2 + scatt = MultiGroupXS() + scatt.LoadFromOpenSn("data/H2O_mg.xs") + + fissile = MultiGroupXS() + fissile.LoadFromOpenSn("data/URRb.xs") + + # Angle information + n_angles = 32 # Number of discrete angles + scat_order = 0 # Scattering order + + pquad = GLProductQuadrature1DSlab(n_polar=n_angles, + scattering_order=scat_order) + + # Create and configure the discrete ordinates solver + phys = DiscreteOrdinatesProblem( + mesh=grid, + num_groups=num_groups, + groupsets=[ + { + "groups_from_to": (0, num_groups - 1), + "angular_quadrature": pquad, + "inner_linear_method": "petsc_gmres", + "l_max_its": 50, + "gmres_restart_interval": 50, + "l_abs_tol": 1.0e-10, + }, + ], + xs_map=[ + {"block_ids": [0], "xs": scatt}, + {"block_ids": [1], "xs": fissile} + ], + boundary_conditions=[ + {"name": "zmin", "type": "reflecting"}, + ], + options={ + "use_precursors": False, + "verbose_inner_iterations": False, + "verbose_outer_iterations": True, + } + ) + + if phase == "online": + new_point = [p0, p1, p2, p3, p4, p5, p6, p7] + print(new_point) + rom_options = { + "param_id": pid, + "phase": phase, + "param_file": "data/params.txt", + "new_point": new_point + } + else: + rom_options = { + "param_id": pid, + "phase": phase + } + + rom = ROMProblem(problem=phys,options=rom_options) + + # Initialize and execute solver + k_solver = PowerIterationROMSolver(problem=phys, rom_problem=rom, k_tol=1.0e-7) + k_solver.Initialize() + k_solver.Execute() + + try: + if phase == "online" and saveh5: + phys.WriteFluxMoments("output/rom_{}_".format(pid)) + np.savetxt("output/rom_k_{}.txt".format(pid), [k_solver.GetEigenvalue()]) + if phase == "offline" and saveh5: + phys.WriteFluxMoments("output/fom_{}_".format(pid)) + np.savetxt("output/fom_k_{}.txt".format(pid), [k_solver.GetEigenvalue()]) + except: + if phase == "online": + phys.WriteFluxMoments("output/rom") + if phase == "offline": + phys.WriteFluxMoments("output/fom") \ No newline at end of file diff --git a/ROM/examples/1dk/run_rom_P58.py b/ROM/examples/1dk/run_rom_P58.py new file mode 100644 index 0000000..79e553e --- /dev/null +++ b/ROM/examples/1dk/run_rom_P58.py @@ -0,0 +1,45 @@ +# run_rom_P58.py +from pathlib import Path +import argparse +import os, sys + +python_root = os.environ.get("OPENSN_PYTHON_PATH") +if python_root: + sys.path.insert(0, python_root) + +from job_manager import JobManager +from P58_problem import P58Problem +from rom_driver import run_pipeline + +def main(): + ap = argparse.ArgumentParser() + ap.add_argument( + "--exe", + type=str, + default=None, + help="OpenSn application executable (e.g. opensn, ./opensn, path/to/app)", + ) + ap.add_argument( + "--system", + type=str, + default="auto", + help="Execution system: auto, slurm, local, etc.", + ) + args = ap.parse_args() + + repo_root = Path.cwd() + + # Pass executable into the JobManager + jm = JobManager( + system=args.system, + opensn_exe=args.exe, + ) + + problem = P58Problem(repo_root, ntrain=50) + + run_pipeline(problem, repo_root, jm) + problem.plot_results() + + +if __name__ == "__main__": + main() diff --git a/ROM/examples/2gcheckerboard/README.md b/ROM/examples/2gcheckerboard/README.md new file mode 100644 index 0000000..c7349c9 --- /dev/null +++ b/ROM/examples/2gcheckerboard/README.md @@ -0,0 +1,34 @@ +# 2gcheckerboard/ + +**Location:** [repo root](../../README.md) / [examples/](../README.md) / `2gcheckerboard/` + +Two-group checkerboard/lattice ROM example. + +This case extends the checkerboard configuration to a two-group problem +with separate absorber and scatterer material definitions. + +--- + +## Files + +- `base_2gcheckerboard.py` + Base two-group OpenSn deck. + +- `checkerboard_problem_2g.py` + Problem definition and ROM configuration. + +- `run_rom_2gcheckerboard.py` + Entry-point script for running the ROM workflow. + +- `absorber_base.txt` + Baseline absorber material definition. + +- `scatterer_base.txt` + Baseline scatterer material definition. + +--- + +## How to Run + +```bash +python run_rom_2gcheckerboard.py --exe=path/to/app/exe diff --git a/ROM/examples/2gcheckerboard/checkerboard_problem_2g.py b/ROM/examples/2gcheckerboard/checkerboard_problem_2g.py index c7c5eab..bf97e67 100644 --- a/ROM/examples/2gcheckerboard/checkerboard_problem_2g.py +++ b/ROM/examples/2gcheckerboard/checkerboard_problem_2g.py @@ -1,8 +1,10 @@ from pathlib import Path -import utils -import plotting import numpy as np +import plotting +import utils +from xs import CrossSections + class CheckerboardProblem2G: def __init__(self, workdir, nprocs=4, ntrain=100, ntest=10): @@ -13,35 +15,83 @@ def __init__(self, workdir, nprocs=4, ntrain=100, ntest=10): self.ntrain = ntrain self.ntest = ntest + self.xs = CrossSections( + [ + { + "in_file": self.workdir / "scatterer_base.txt", + "out_file": self.workdir / "data" / "scatterer.xs", + "kind": "total", + }, + { + "in_file": self.workdir / "absorber_base.txt", + "out_file": self.workdir / "data" / "absorber.xs", + "kind": "total", + }, + ], + param_mode="entrywise", + ) + + self.nominal_sample = self.xs.get_nominal_sample() + self.full_bounds = list(self.xs.get_bounds()) + + # Keep the original 2-parameter interface: + # pvec[0] -> scatterer downscatter entry S(0,1) + # pvec[1] -> absorber absorption sigma_a applied to both groups + self._scatterer_downscatter_idx = self._find_total_transfer_index( + self.workdir / "scatterer_base.txt", 0, 1 + ) + self._absorber_abs_indices = [ + self._find_total_abs_index(self.workdir / "absorber_base.txt", 0), + self._find_total_abs_index(self.workdir / "absorber_base.txt", 1), + ] + + # Override the active bounds to preserve the original checkerboard + # ranges while still using CrossSections for the flattened XS layout. + self.full_bounds[self._scatterer_downscatter_idx] = (0.5, 1.0) + for idx in self._absorber_abs_indices: + self.full_bounds[idx] = (7.5, 12.5) + + self.bounds = [ + list(self.full_bounds[self._scatterer_downscatter_idx]), + list(self.full_bounds[self._absorber_abs_indices[0]]), + ] + + def _find_total_abs_index(self, in_file, group): + block = self.xs.get_block(in_file) + return block["slice"].start + group + + def _find_total_transfer_index(self, in_file, gfrom, gto): + block = self.xs.get_block(in_file) + try: + entry_idx = block["transfer_entries"].index((gfrom, gto)) + except ValueError as exc: + raise ValueError( + "Transfer entry ({}, {}) not found in {}.".format(gfrom, gto, in_file) + ) from exc + + return block["slice"].start + block["G"] + entry_idx + + def _build_full_sample(self, pvec): + full_sample = self.nominal_sample.copy() + full_sample[self._scatterer_downscatter_idx] = pvec[0] + for idx in self._absorber_abs_indices: + full_sample[idx] = pvec[1] + return full_sample + def sample_training(self): - bounds = [[0.5,1.0],[7.5,12.5]] - self.training_set = utils.sample_parameter_space(bounds, self.ntrain) + self.training_set = utils.sample_parameter_space(self.bounds, self.ntrain) params_path = self.workdir / "data" / "params.txt" np.savetxt(str(params_path), self.training_set) - def sample_testing(self): - test_scatt_1 = np.random.uniform(0.5,1.0,10) - test_abs_1 = np.random.uniform(7.5,12.5,10) - self.testing_set = np.append(test_scatt_1[:,np.newaxis], test_abs_1[:,np.newaxis], axis=1) + self.testing_set = utils.sample_test(self.bounds, self.ntest) params_path = self.workdir / "data" / "test_params.txt" np.savetxt(str(params_path), self.testing_set) def update_xs(self, pvec): - S_abs = [[0.0, 0.0], - [0.0, 0.0]] - sigma_t_scatt = [1.0, 1.0] - - S_scatt = [[1-pvec[0], pvec[0]], - [0.0, 1.0]] - utils.update_xs("scatterer_base.txt", "data/scatterer.xs", sigma_t_scatt, S_scatt) - - sigma_t_abs = [pvec[1], pvec[1]] - utils.update_xs("absorber_base.txt", "data/absorber.xs", sigma_t_abs, S_abs) - - + self.xs.write_sample(self._build_full_sample(pvec)) def plot_results(self): plotting.plot_sv(num_groups=2) @@ -53,16 +103,25 @@ def plot_results(self): fom_time = np.loadtxt(str(results_dir / "offline_time_{}.txt".format(i))) output_dir = self.workdir / "output" - plotting.plot_2d_flux(str(output_dir / ("fom_{}_".format(i) + "{}.h5")), ranks=range(4), prefix="fom", pid=i) - plotting.plot_2d_flux(str(output_dir / ("rom_{}_".format(i) + "{}.h5")), ranks=range(4), prefix="rom", pid=i) + plotting.plot_2d_flux( + str(output_dir / ("fom_{}_".format(i) + "{}.h5")), + ranks=range(4), + prefix="fom", + pid=i, + ) + plotting.plot_2d_flux( + str(output_dir / ("rom_{}_".format(i) + "{}.h5")), + ranks=range(4), + prefix="rom", + pid=i, + ) error = plotting.plot_2d_lineout(output_dir, ranks=range(4), pid=i) errors.append(error) - speedups.append(fom_time/rom_time) + speedups.append(fom_time / rom_time) print("Avg Error ", np.mean(errors)) np.savetxt(str(results_dir / "errors.txt"), errors) print("Avg Speedup ", np.mean(speedups)) np.savetxt(str(results_dir / "speedups.txt"), speedups) - diff --git a/ROM/examples/README.md b/ROM/examples/README.md new file mode 100644 index 0000000..5e399ad --- /dev/null +++ b/ROM/examples/README.md @@ -0,0 +1,35 @@ +# examples/ + +**Location:** [repo root](../README.md) / `examples/` + +Example input decks and run scripts demonstrating how to use the application and ROM pipeline. + +Each example folder typically contains: +- a base OpenSn input deck, +- a problem definition helper, +- a `run_rom_*.py` entry script you execute, +- optional reference data (e.g., base cross sections). + +## Examples + +- [`reed/`](reed/README.md) + Reed 1-D benchmark example. + +- [`checkerboard/`](checkerboard/README.md) + 2-D Checkerboard/lattice problem configuration. + +- [`2gcheckerboard/`](2gcheckerboard/README.md) + Two-group checkerboard variant with absorber/scatterer material files. + +- [`1dk/`](1dk/README.md) + 1-D k-eigenvalue test problem with H2O/U material files. + + +## Running + +Each example folder includes a `run_rom_*.py` script. Start with the folder README for details. + +Execution pattern: + +```bash +python run_rom_.py --exe=path/to/app/exe diff --git a/ROM/examples/checkerboard/README.md b/ROM/examples/checkerboard/README.md new file mode 100644 index 0000000..2db3edd --- /dev/null +++ b/ROM/examples/checkerboard/README.md @@ -0,0 +1,25 @@ +# checkerboard/ + +**Location:** [repo root](../../README.md) / [examples/](../README.md) / `checkerboard/` + +Single-group checkerboard/lattice example and ROM driver. + +--- + +## Files + +- `base_checkerboard.py` + Base deck defining geometry, materials, and solver setup. + +- `checkerboard_problem.py` + Problem definition and ROM configuration. + +- `run_rom_checkerboard.py` + Entry-point script for running the ROM workflow. + +--- + +## How to Run + +```bash +python run_rom_checkerboard.py --exe=path/to/app/exe diff --git a/ROM/examples/checkerboard/run_rom_checkerboard.py b/ROM/examples/checkerboard/run_rom_checkerboard.py index 46bb3c8..0d18683 100644 --- a/ROM/examples/checkerboard/run_rom_checkerboard.py +++ b/ROM/examples/checkerboard/run_rom_checkerboard.py @@ -42,4 +42,4 @@ def main(): if __name__ == "__main__": - main() + main() \ No newline at end of file diff --git a/ROM/examples/reed/README.md b/ROM/examples/reed/README.md new file mode 100644 index 0000000..af2da46 --- /dev/null +++ b/ROM/examples/reed/README.md @@ -0,0 +1,27 @@ +# reed/ + +**Location:** [repo root](../../README.md) / [examples/](../README.md) / `reed/` + +Reed example and ROM driver. + +--- + +## Files + +- `base_reed.py` + Base OpenSn deck for the Reed problem. + +- `reed_problem.py` + Problem definition and parameterization for the Reed example. + +- `run_rom_reed.py` + Entry-point script that runs the ROM workflow. + +--- + +## How to Run + +From this directory: + +```bash +python run_rom_reed.py --exe=path/to/app/exe diff --git a/ROM/examples/reed/base_reed.py b/ROM/examples/reed/base_reed.py index 6d4fa91..114f494 100644 --- a/ROM/examples/reed/base_reed.py +++ b/ROM/examples/reed/base_reed.py @@ -124,4 +124,3 @@ if phase == "offline": phys.WriteFluxMoments("output/fom") - diff --git a/ROM/examples/reed/run_rom_reed.py b/ROM/examples/reed/run_rom_reed.py index 7839ccd..34feb1f 100644 --- a/ROM/examples/reed/run_rom_reed.py +++ b/ROM/examples/reed/run_rom_reed.py @@ -42,4 +42,4 @@ def main(): if __name__ == "__main__": - main() + main() \ No newline at end of file diff --git a/ROM/python/README.md b/ROM/python/README.md new file mode 100644 index 0000000..f1e70a9 --- /dev/null +++ b/ROM/python/README.md @@ -0,0 +1,24 @@ +# python/ + +**Location:** [repo root](../README.md) / `python/` + +Python drivers and utilities for running example problems and ROM workflows using the application implemented in [`src/`](../src/README.md). + +This folder is intended to: +- standardize launching (local vs batch/HPC), +- keep ROM pipeline orchestration in one place, +- centralize plotting and small utilities. + +## Key modules + +- [`rom_driver.py`](rom_driver.py) + Orchestrates the ROM pipeline (e.g., offline, merge, system, online phases). + +- [`job_manager.py`](job_manager.py) + Launch abstraction for running the executable (local execution vs MPI/Slurm-style workflows). + +- [`plotting.py`](plotting.py) + Plotting helpers for example runs and ROM result summaries. + +- [`utils.py`](utils.py) + Shared utilities (loading fluxes from h5, basic sampling capabilities, and updating OpenSn .xs files). \ No newline at end of file diff --git a/ROM/python/plotting.py b/ROM/python/plotting.py index acfc876..450173c 100644 --- a/ROM/python/plotting.py +++ b/ROM/python/plotting.py @@ -79,6 +79,45 @@ def plot_2d_lineout(output_dir, ranks, y_target=4.0, moment=0, grid_res=200, pid error = np.linalg.norm(np.asarray(vals_)-np.asarray(vals))/np.linalg.norm(np.asarray(vals_)) return error +def plot_2d_lineout_eig(output_dir, ranks, y_target=4.0, moment=0, grid_res=200, pid=0): + """Plot lineout at y_target of ROM and FOM.""" + xs, ys, vals, G = load_2d_flux(str(output_dir / ("fom_{}_".format(pid) + "{}.h5")), ranks, moment=moment) + xs_, ys_, vals_, G = load_2d_flux(str(output_dir / ("rom_{}_".format(pid) + "{}.h5")), ranks, moment=moment) + + for g in range(G): + vals[g] /= np.linalg.norm(vals[g]) + vals_[g] /= np.linalg.norm(vals_[g]) + xi = np.linspace(xs[g].min(), xs[g].max(), grid_res) + yi = np.linspace(ys[g].min(), ys[g].max(), grid_res) + X, Y = np.meshgrid(xi, yi) + + # Interpolate data onto grid + Z = scipy.interpolate.griddata((xs[g], ys[g]), vals[g], (X, Y), method="linear") + Z_ = scipy.interpolate.griddata((xs[g], ys[g]), vals_[g], (X, Y), method="linear") + + # Find the row index closest to y=4 + row_idx = np.argmin(np.abs(yi - y_target)) + + # Extract data along y = 4 + rom_line = Z[row_idx, :] + fom_line = Z_[row_idx, :] + + # Plot ROM vs FOM + plt.figure(figsize=(8,5)) + plt.plot(xi, rom_line, label='ROM', color='blue') + plt.plot(xi, fom_line, label='FOM', color='orange', linestyle='--') + plt.xlabel('X') + plt.ylabel('Scalar Value') + plt.title(f'ROM vs FOM at y={y_target}') + plt.legend() + plt.grid(True) + plt.tight_layout() + plt.savefig(f"results/line_y{y_target}_rom_fom_{pid}_{g}.jpg") + plt.close() + + + error = np.linalg.norm(np.abs(np.asarray(vals_))-np.abs(np.asarray(vals)))/np.linalg.norm(np.asarray(vals_)) + return error def plot_sv(num_groups): for i in range(num_groups): @@ -112,5 +151,30 @@ def plot_1d_flux(fom_pattern, rom_pattern, ranks, moment=0, prefix="reed_ommi", plt.savefig(outpath, dpi=200) plt.close() + error = np.linalg.norm(np.array(rom_vals) - np.array(fom_vals)) / np.linalg.norm(fom_vals) + return error + +def plot_1d_eigenvector(fom_pattern, rom_pattern, ranks, moment=0, prefix="reed_ommi", pid=0): + """Compare FOM vs ROM 1-D flux.""" + fom_x, fom_vals, G = load_1d_flux(fom_pattern, ranks, moment=moment) + rom_x, rom_vals, G = load_1d_flux(rom_pattern, ranks, moment=moment) + + errors = [] + for g in range(G): + rom_vals[g] /= np.linalg.norm(rom_vals[g]) + rom_vals = np.abs(rom_vals) + fom_vals[g] /= np.linalg.norm(fom_vals[g]) + plt.figure(figsize=(6, 4)) + plt.plot(fom_x[g], fom_vals[g], "-", label="FOM") + plt.plot(rom_x[g], rom_vals[g], "--", label="ROM") + plt.xlabel("x") + plt.ylabel("Flux") + plt.grid() + plt.legend() + outpath = f"results/{prefix}_{pid}_g_{g}.png" + plt.tight_layout() + plt.savefig(outpath, dpi=200) + plt.close() + error = np.linalg.norm(np.array(rom_vals) - np.array(fom_vals)) / np.linalg.norm(fom_vals) return error \ No newline at end of file diff --git a/ROM/python/rom_driver.py b/ROM/python/rom_driver.py index bc0acef..a74a145 100644 --- a/ROM/python/rom_driver.py +++ b/ROM/python/rom_driver.py @@ -26,6 +26,8 @@ def make_opensn_args(phase, pid, pvec, save_h5=False): if save_h5: args.extend(["-p", "saveh5=True"]) + else: + args.extend(["-p", "saveh5=False"]) return args diff --git a/ROM/python/utils.py b/ROM/python/utils.py index 5b945e4..505ed4e 100644 --- a/ROM/python/utils.py +++ b/ROM/python/utils.py @@ -1,6 +1,8 @@ import numpy as np import subprocess import h5py +from scipy.stats import qmc + def load_2d_flux(file_pattern, ranks, moment=0): """Load (x, y, flux) grouped by energy group from HDF5 files.""" @@ -86,33 +88,28 @@ def sample_parameter_space(bounds, n_samples): samples = np.vstack([random_samples, vertices]) return samples +def sample_LHS(bounds, n_samples): + """Performs Latin Hypercube Sampling on a set of bounds""" + d = len(bounds) -def update_xs(in_file, out_file, sigma_t_vec, S): - """Writes a OpenSn .xs file replacing cross sections with user supplied values""" - with open(in_file, "r") as f: - lines = f.readlines() + # Create LHS sampler + sampler = qmc.LatinHypercube(d=d) - # --- SIGMA_T block --- - b = next(i for i, s in enumerate(lines) if "SIGMA_T_BEGIN" in s) - e = next(i for i, s in enumerate(lines) if "SIGMA_T_END" in s) - for i in range(b+1, e): - toks = lines[i].split() - g = int(toks[0]) - toks[1] = f"{float(sigma_t_vec[g]):.12g}" - lines[i] = " ".join(toks) + "\n" + # Generate samples in [0, 1]^d + unit_samples = sampler.random(n=n_samples) - # --- TRANSFER_MOMENTS block --- - tb = next(i for i, s in enumerate(lines) if "TRANSFER_MOMENTS_BEGIN" in s) - te = next(i for i, s in enumerate(lines) if "TRANSFER_MOMENTS_END" in s) + # Scale to physical bounds + lows = np.array([low for (low, high) in bounds]) + highs = np.array([high for (low, high) in bounds]) - G = len(sigma_t_vec) - new_tm = [] - for gprime in range(G): - for g in range(G): - val = float(S[gprime][g]) - new_tm.append(f"M_GFROM_GTO_VAL 0 {gprime} {g} {val:.12g}\n") + samples = qmc.scale(unit_samples, lows, highs) + return samples - lines[tb+1:te] = new_tm +def sample_test(bounds, n_samples): + """Performs uniform sampling on a set of bounds to generate a test set""" + samples = np.array([ + [np.random.uniform(low, high) for (low, high) in bounds] + for _ in range(n_samples) + ]) - with open(out_file, "w") as f: - f.writelines(lines) \ No newline at end of file + return samples diff --git a/ROM/python/xs.py b/ROM/python/xs.py new file mode 100644 index 0000000..18d3c47 --- /dev/null +++ b/ROM/python/xs.py @@ -0,0 +1,426 @@ +from pathlib import Path +import numpy as np + + +class CrossSections: + """ + Build a single parameter space across multiple XS files. + + Modes + ----- + param_mode="entrywise" + Original behavior: + - total XS: [sigma_a_vec, transfer_vals] + - fission XS: [sigma_f_vec, sigma_c_vec, transfer_vals] + + param_mode="block_scales" + Sample multiplicative scale factors instead of every individual entry: + - total XS: [sigma_a_scale, scatter_scale] + - fission XS: [sigma_f_scale, sigma_c_scale, scatter_scale] + + Each scale multiplies the entire corresponding block. + For example, sigma_f_scale multiplies all group values in sigma_f. + """ + + def __init__(self, xs_specs, frac=0.2, transfer_tol=0.0, param_mode="entrywise"): + self.frac = frac + self.transfer_tol = transfer_tol + self.param_mode = param_mode + + if self.param_mode not in ("entrywise", "block_scales"): + raise ValueError( + "Unknown param_mode '{}'. Expected 'entrywise' or 'block_scales'.".format( + self.param_mode + ) + ) + + self.blocks = [] + self.bounds = [] + self.n_params = 0 + + start = 0 + for spec in xs_specs: + in_file = Path(spec["in_file"]) + out_file = Path(spec["out_file"]) + kind = spec["kind"] + + if kind == "total": + data = self._read_total_xs(in_file) + if self.param_mode == "entrywise": + x_nominal = self._flatten_total_data(data) + else: + x_nominal = self._flatten_total_scales(data) + + elif kind == "fission": + data = self._read_fission_xs(in_file) + if self.param_mode == "entrywise": + x_nominal = self._flatten_fission_data(data) + else: + x_nominal = self._flatten_fission_scales(data) + + else: + raise ValueError( + "Unknown XS kind '{}'. Expected 'total' or 'fission'.".format(kind) + ) + + stop = start + len(x_nominal) + + block = { + "in_file": in_file, + "out_file": out_file, + "kind": kind, + "G": data["G"], + "transfer_entries": data["transfer_entries"], + "x_nominal": x_nominal, + "bounds": self._relative_bounds(x_nominal), + "slice": slice(start, stop), + "n_params": len(x_nominal), + "data": data, + } + + self.blocks.append(block) + self.bounds.extend(block["bounds"]) + start = stop + + self.n_params = start + + def get_bounds(self): + return self.bounds + + def get_nominal_sample(self): + return np.concatenate([block["x_nominal"] for block in self.blocks]) + + def random_sample(self, rng=None): + if rng is None: + rng = np.random.default_rng() + + vals = [rng.uniform(lo, hi) for (lo, hi) in self.bounds] + return np.asarray(vals, dtype=float) + + def get_block(self, in_file): + in_file = Path(in_file) + for block in self.blocks: + if block["in_file"] == in_file: + return block + raise KeyError("{} not found in CrossSections.".format(in_file)) + + def get_block_sample(self, sample, in_file): + block = self.get_block(in_file) + return np.asarray(sample[block["slice"]], dtype=float) + + def write_file(self, sample, in_file): + block = self.get_block(in_file) + x_block = self.get_block_sample(sample, in_file) + + if block["kind"] == "total": + if self.param_mode == "entrywise": + sigma_a_vec, S = self._unflatten_total_data(x_block, block) + else: + sigma_a_vec, S = self._unflatten_total_scales(x_block, block) + + self._update_xs( + block["in_file"], + block["out_file"], + sigma_a_vec, + S, + ) + + elif block["kind"] == "fission": + if self.param_mode == "entrywise": + sigma_f_vec, sigma_c_vec, S = self._unflatten_fission_data(x_block, block) + else: + sigma_f_vec, sigma_c_vec, S = self._unflatten_fission_scales(x_block, block) + + self._update_fission_xs( + block["in_file"], + block["out_file"], + sigma_f_vec, + sigma_c_vec, + S, + ) + else: + raise ValueError("Unknown XS kind '{}'.".format(block["kind"])) + + def write_sample(self, sample): + sample = np.asarray(sample, dtype=float) + if len(sample) != self.n_params: + raise ValueError( + "Sample has length {}, expected {}.".format(len(sample), self.n_params) + ) + + for block in self.blocks: + self.write_file(sample, block["in_file"]) + + # ------------------------------------------------------------------------- + # XS writers + # ------------------------------------------------------------------------- + + def _update_xs(self, in_file, out_file, sigma_a_vec, S): + with open(in_file, "r") as f: + lines = f.readlines() + + sigma_t_vec = np.zeros_like(sigma_a_vec, dtype=float) + + G = len(sigma_a_vec) + for gfrom in range(G): + sigma_t_vec[gfrom] = float(sigma_a_vec[gfrom]) + for gto in range(G): + sigma_t_vec[gfrom] += float(S[gfrom, gto]) + + b = next(i for i, s in enumerate(lines) if "SIGMA_T_BEGIN" in s) + e = next(i for i, s in enumerate(lines) if "SIGMA_T_END" in s) + for i in range(b + 1, e): + toks = lines[i].split() + g = int(toks[0]) + toks[1] = "{:.12g}".format(float(sigma_t_vec[g])) + lines[i] = " ".join(toks) + "\n" + + tb = next(i for i, s in enumerate(lines) if "TRANSFER_MOMENTS_BEGIN" in s) + te = next(i for i, s in enumerate(lines) if "TRANSFER_MOMENTS_END" in s) + + new_tm = [] + for gfrom in range(G): + for gto in range(G): + val = float(S[gfrom, gto]) + if abs(val) > self.transfer_tol: + new_tm.append( + "M_GFROM_GTO_VAL 0 {} {} {:.12g}\n".format(gfrom, gto, val) + ) + + lines[tb + 1:te] = new_tm + + with open(out_file, "w") as f: + f.writelines(lines) + + def _update_fission_xs(self, in_file, out_file, sigma_f_vec, sigma_c_vec, S): + with open(in_file, "r") as f: + lines = f.readlines() + + sigma_t_vec = np.zeros_like(sigma_f_vec, dtype=float) + + b = next(i for i, s in enumerate(lines) if "SIGMA_F_BEGIN" in s) + e = next(i for i, s in enumerate(lines) if "SIGMA_F_END" in s) + for i in range(b + 1, e): + toks = lines[i].split() + g = int(toks[0]) + toks[1] = "{:.12g}".format(float(sigma_f_vec[g])) + lines[i] = " ".join(toks) + "\n" + sigma_t_vec[g] = float(sigma_f_vec[g]) + float(sigma_c_vec[g]) + + tb = next(i for i, s in enumerate(lines) if "TRANSFER_MOMENTS_BEGIN" in s) + te = next(i for i, s in enumerate(lines) if "TRANSFER_MOMENTS_END" in s) + + G = len(sigma_t_vec) + new_tm = [] + for gfrom in range(G): + for gto in range(G): + val = float(S[gfrom, gto]) + if abs(val) > self.transfer_tol: + new_tm.append( + "M_GFROM_GTO_VAL 0 {} {} {:.12g}\n".format(gfrom, gto, val) + ) + sigma_t_vec[gfrom] += val + + lines[tb + 1:te] = new_tm + + b = next(i for i, s in enumerate(lines) if "SIGMA_T_BEGIN" in s) + e = next(i for i, s in enumerate(lines) if "SIGMA_T_END" in s) + for i in range(b + 1, e): + toks = lines[i].split() + g = int(toks[0]) + toks[1] = "{:.12g}".format(float(sigma_t_vec[g])) + lines[i] = " ".join(toks) + "\n" + + with open(out_file, "w") as f: + f.writelines(lines) + + # ------------------------------------------------------------------------- + # Readers + # ------------------------------------------------------------------------- + + def _read_lines(self, xs_file): + with open(xs_file, "r") as f: + return f.readlines() + + def _find_block(self, lines, begin_marker, end_marker): + b = next(i for i, s in enumerate(lines) if begin_marker in s) + e = next(i for i, s in enumerate(lines) if end_marker in s) + return b, e + + def _read_num_groups(self, lines): + line = next(s for s in lines if s.strip().startswith("NUM_GROUPS")) + return int(line.split()[1]) + + def _read_vector_block(self, lines, begin_marker, end_marker): + G = self._read_num_groups(lines) + vec = np.zeros(G, dtype=float) + + b, e = self._find_block(lines, begin_marker, end_marker) + for line in lines[b + 1:e]: + toks = line.split() + if len(toks) < 2: + continue + g = int(toks[0]) + vec[g] = float(toks[1]) + + return vec + + def _read_transfer_entries(self, lines): + entries = [] + values = [] + + b, e = self._find_block(lines, "TRANSFER_MOMENTS_BEGIN", "TRANSFER_MOMENTS_END") + for line in lines[b + 1:e]: + toks = line.split() + + if len(toks) < 5: + continue + + key = toks[0] + if key not in ("M_GFROM_GTO_VAL", "M_GPRIME_G_VAL"): + continue + + ell = int(toks[1]) + if ell != 0: + continue + + gfrom = int(toks[2]) + gto = int(toks[3]) + val = float(toks[4]) + + entries.append((gfrom, gto)) + values.append(val) + + return entries, np.array(values, dtype=float) + + def _entries_to_matrix(self, G, entries, values): + S = np.zeros((G, G), dtype=float) + for (gfrom, gto), val in zip(entries, values): + S[gfrom, gto] = val + return S + + def _read_total_xs(self, xs_file): + lines = self._read_lines(xs_file) + G = self._read_num_groups(lines) + + sigma_t_vec = self._read_vector_block(lines, "SIGMA_T_BEGIN", "SIGMA_T_END") + transfer_entries, transfer_vals = self._read_transfer_entries(lines) + + S = self._entries_to_matrix(G, transfer_entries, transfer_vals) + sigma_a_vec = sigma_t_vec - S.sum(axis=1) + + return { + "kind": "total", + "file": Path(xs_file), + "G": G, + "sigma_t_vec": sigma_t_vec, + "sigma_a_vec": sigma_a_vec, + "transfer_entries": transfer_entries, + "transfer_vals": transfer_vals, + } + + def _read_fission_xs(self, xs_file): + lines = self._read_lines(xs_file) + G = self._read_num_groups(lines) + + sigma_t_vec = self._read_vector_block(lines, "SIGMA_T_BEGIN", "SIGMA_T_END") + sigma_f_vec = self._read_vector_block(lines, "SIGMA_F_BEGIN", "SIGMA_F_END") + transfer_entries, transfer_vals = self._read_transfer_entries(lines) + + S = self._entries_to_matrix(G, transfer_entries, transfer_vals) + sigma_c_vec = sigma_t_vec - sigma_f_vec - S.sum(axis=1) + + return { + "kind": "fission", + "file": Path(xs_file), + "G": G, + "sigma_t_vec": sigma_t_vec, + "sigma_f_vec": sigma_f_vec, + "sigma_c_vec": sigma_c_vec, + "transfer_entries": transfer_entries, + "transfer_vals": transfer_vals, + } + + # ------------------------------------------------------------------------- + # Flatten / unflatten: entrywise + # ------------------------------------------------------------------------- + + def _flatten_total_data(self, data): + return np.concatenate([data["sigma_a_vec"], data["transfer_vals"]]) + + def _flatten_fission_data(self, data): + return np.concatenate([data["sigma_f_vec"], data["sigma_c_vec"], data["transfer_vals"]]) + + def _unflatten_total_data(self, x_block, block): + G = block["G"] + n_transfer = len(block["transfer_entries"]) + + sigma_a_vec = np.array(x_block[:G], dtype=float) + transfer_vals = np.array(x_block[G:G + n_transfer], dtype=float) + S = self._entries_to_matrix(G, block["transfer_entries"], transfer_vals) + + return sigma_a_vec, S + + def _unflatten_fission_data(self, x_block, block): + G = block["G"] + n_transfer = len(block["transfer_entries"]) + + i0 = 0 + i1 = i0 + G + i2 = i1 + G + i3 = i2 + n_transfer + + sigma_f_vec = np.array(x_block[i0:i1], dtype=float) + sigma_c_vec = np.array(x_block[i1:i2], dtype=float) + transfer_vals = np.array(x_block[i2:i3], dtype=float) + S = self._entries_to_matrix(G, block["transfer_entries"], transfer_vals) + + return sigma_f_vec, sigma_c_vec, S + + # ------------------------------------------------------------------------- + # Flatten / unflatten: block scales + # ------------------------------------------------------------------------- + + def _flatten_total_scales(self, data): + return np.array([1.0, 1.0], dtype=float) + + def _flatten_fission_scales(self, data): + return np.array([1.0, 1.0, 1.0], dtype=float) + + def _unflatten_total_scales(self, x_block, block): + data = block["data"] + + sigma_a_scale = float(x_block[0]) + scatter_scale = float(x_block[1]) + + sigma_a_vec = sigma_a_scale * np.asarray(data["sigma_a_vec"], dtype=float) + transfer_vals = scatter_scale * np.asarray(data["transfer_vals"], dtype=float) + S = self._entries_to_matrix(block["G"], block["transfer_entries"], transfer_vals) + + return sigma_a_vec, S + + def _unflatten_fission_scales(self, x_block, block): + data = block["data"] + + sigma_f_scale = float(x_block[0]) + sigma_c_scale = float(x_block[1]) + scatter_scale = float(x_block[2]) + + sigma_f_vec = sigma_f_scale * np.asarray(data["sigma_f_vec"], dtype=float) + sigma_c_vec = sigma_c_scale * np.asarray(data["sigma_c_vec"], dtype=float) + transfer_vals = scatter_scale * np.asarray(data["transfer_vals"], dtype=float) + S = self._entries_to_matrix(block["G"], block["transfer_entries"], transfer_vals) + + return sigma_f_vec, sigma_c_vec, S + + # ------------------------------------------------------------------------- + # Bounds + # ------------------------------------------------------------------------- + + def _relative_bounds(self, x_nominal): + x_nominal = np.asarray(x_nominal, dtype=float) + a = (1.0 - self.frac) * x_nominal + b = (1.0 + self.frac) * x_nominal + lower = np.minimum(a, b) + upper = np.maximum(a, b) + return list(zip(lower, upper)) \ No newline at end of file diff --git a/ROM/src/README.md b/ROM/src/README.md new file mode 100644 index 0000000..c0ccebe --- /dev/null +++ b/ROM/src/README.md @@ -0,0 +1,22 @@ +# src/ + +**Location:** [repo root](../README.md) / `src/` + +C++ sources for the application. This folder contains the main entry point, ROM implementation, and the Python-app bridge used to integrate with OpenSn problems. + +## Key files + +- [`main.cc`](main.cc) + Application entry point. + +- [`rom.cc`](rom.cc) + Top-level ROM plumbing/glue code. + +- [`rom_py_app.h`](rom_py_app.h), [`rom_py_app.cc`](rom_py_app.cc) + ROM “Python app” integration layer (app interface invoked from Python-side workflows). + +- [`py_wrappers.h`](py_wrappers.h) + C++ declarations for Python-facing wrappers (bindings/glue utilities). + +- [`rom/`](rom/README.md) + ROM core implementation (problem definition, solver, structs). \ No newline at end of file diff --git a/ROM/src/rom.cc b/ROM/src/rom.cc index 4d90902..b4c92d0 100644 --- a/ROM/src/rom.cc +++ b/ROM/src/rom.cc @@ -4,6 +4,7 @@ #include "py_wrappers.h" #include "rom/rom_problem.h" #include "rom/steady_state_rom_solver.h" +#include "rom/pi_keigen_rom_solver.h" #include "modules/solver.h" #include #include @@ -119,6 +120,67 @@ void WrapROM(py::module& m) - 'systems' : assemble reduced systems and write libROM files - 'online' : interpolate and solve reduced system )"); + + // PowerIterationKEigenROMSolver + auto pi_rom_solver = + py::class_, + Solver>( + m, + "PowerIterationROMSolver", + R"( + k eigen ROM driver. + + Wrapper of :cpp:class:`opensn::PowerIterationKEigenROMSolver`. + + Parameters (kwargs) + ------------------- + problem : LBSProblem + The full-order transport problem. + rom_problem : ROMProblem + The ROM controller (bases, reduced systems, interpolation). + name : str + Required solver name for logging/monitors. + )" + ); + + pi_rom_solver.def( + py::init([](py::kwargs kw) + { + auto params = KwargsToParams(kw); + + return std::make_shared(params); + }), + R"( + PowerIterationKEigenROMSolver(**kwargs) + + Construct a k-eigen driver that dispatches to ROM or FOM paths + depending on the ROM options and phase. + )" + ); + + pi_rom_solver + .def("Initialize", &PowerIterationKEigenROMSolver::Initialize, + R"( + Initialize() + + Prepare the solver and ROM controller for execution. + )") + .def("Execute", &PowerIterationKEigenROMSolver::Execute, + R"( + Execute() + + Run the solve. Behavior depends on the ROM phase: + - 'offline' : full-order solve + snapshot sample + - 'merge' : merge snapshots into bases + - 'systems' : assemble reduced systems and write libROM files + - 'online' : interpolate and solve reduced system + )") + .def("GetEigenvalue", + &PowerIterationKEigenROMSolver::GetEigenvalue, + R"( + Return the current k‑eigenvalue. + )"); } // clang-format on diff --git a/ROM/src/rom/README.md b/ROM/src/rom/README.md new file mode 100644 index 0000000..28721a5 --- /dev/null +++ b/ROM/src/rom/README.md @@ -0,0 +1,22 @@ +# src/rom/ + +**Location:** [repo root](../../README.md) / [src/](../README.md) / `rom/` + +Core reduced-order modeling (ROM) implementation: problem definition, solver(s), and shared ROM data structures. + +## Contents + +### Problem definition + +- [`rom_problem.h`](rom_problem.h), [`rom_problem.cc`](rom_problem.cc) + Defines the ROM problem object and its responsibilities (assembly, operators, state, I/O hooks, etc.). + +### Solver + +- [`steady_state_rom_solver.h`](steady_state_rom_solver.h), [`steady_state_rom_solver.cc`](steady_state_rom_solver.cc) + Steady-state ROM solver implementation. + +### Shared structs / types + +- [`rom_structs.h`](rom_structs.h) + Shared ROM data containers and configuration structs. \ No newline at end of file diff --git a/ROM/src/rom/pi_keigen_rom_solver.cc b/ROM/src/rom/pi_keigen_rom_solver.cc new file mode 100644 index 0000000..bff660f --- /dev/null +++ b/ROM/src/rom/pi_keigen_rom_solver.cc @@ -0,0 +1,139 @@ +// SPDX-FileCopyrightText: 2026 The OpenSn Authors +// SPDX-License-Identifier: MIT + +#include "modules/linear_boltzmann_solvers/lbs_problem/compute/lbs_compute.h" +#include "framework/runtime.h" +#include "pi_keigen_rom_solver.h" +#include +#include +#include + +namespace opensn +{ + +/** Returns the input-parameter schema for PowerIterationKEigenROMSolver. + * + * Extends the base power-iteration k-eigen solver schema with: + * - `rom_problem` : an existing ROMProblem instance that manages ROM workflow. + */ +InputParameters +PowerIterationKEigenROMSolver::GetInputParameters() +{ + // Start from the base PI-k-eigen solver parameters + InputParameters params = PowerIterationKEigenSolver::GetInputParameters(); + + params.SetGeneralDescription( + "Implementation of a k-eigenvalue ROM solver. Offline phase runs the " + "full-order power-iteration k-eigen solver and takes sample with libROM."); + params.ChangeExistingParamToOptional("name", "PowerIterationKEigenROMSolver"); + + params.AddRequiredParameter>( + "rom_problem", "A ROM problem"); + + return params; +} + +PowerIterationKEigenROMSolver::PowerIterationKEigenROMSolver(const InputParameters& params) + : PowerIterationKEigenSolver(params), + lbs_problem_(params.GetSharedPtrParam("problem")), + rom_problem_(params.GetSharedPtrParam("rom_problem")) +{ +} + +void +PowerIterationKEigenROMSolver::Initialize() +{ + PowerIterationKEigenSolver::Initialize(); +} + +/** Executes the requested ROM workflow phase for the k-eigenvalue solver. + * + * Supported phases are: + * - OFFLINE : runs the full-order power iteration and writes snapshots, + * - MERGE : builds the reduced bases from stored snapshots, + * - SYSTEMS : assembles and writes reduced operators, + * - ONLINE : interpolates and solves the reduced k-eigenvalue system. + */ +void +PowerIterationKEigenROMSolver::Execute() +{ + auto& rom_options = rom_problem_->GetOptions(); + auto& options = lbs_problem_->GetOptions(); + + if (rom_options.phase == Phase::OFFLINE) + { + // Run full-order k-eigen solve and time it + auto start = std::chrono::high_resolution_clock::now(); + + PowerIterationKEigenSolver::Execute(); + + auto end = std::chrono::high_resolution_clock::now(); + std::chrono::duration elapsed = end - start; + + if (opensn::mpi_comm.rank() == 0) + { + std::ofstream outfile("results/offline_time_" + std::to_string(rom_options.param_id) + ".txt"); + if (outfile.is_open()) + { + outfile << elapsed.count() << std::endl; + outfile.close(); + } + } + + rom_problem_->TakeSample(rom_options.param_id); + } + if (rom_options.phase == Phase::MERGE) + { + rom_problem_->MergePhase(rom_options.param_id); + } + if (rom_options.phase == Phase::SYSTEMS) + { + std::shared_ptr AU_ = rom_problem_->AssembleAU(); + rom_problem_->LoadUgs(); + std::shared_ptr BU_ = rom_problem_->AssembleBU(); + + const std::string Ar_filename = + "data/rom_system_Ar_" + std::to_string(rom_options.param_id); + const std::string Br_filename = + "data/rom_system_Br_" + std::to_string(rom_options.param_id); + + rom_problem_->AssembleROM(AU_, BU_, Ar_filename, Br_filename); + } + if (rom_options.phase == Phase::ONLINE) + { + rom_problem_->ReadParamMatrix(rom_options.param_file); + rom_problem_->LoadUgs(); + + std::shared_ptr Ar_interp; + std::shared_ptr Br_interp; + + rom_problem_->SetupArInterpolator(*rom_options.new_point); + rom_problem_->SetupBrInterpolator(*rom_options.new_point); + + auto start = std::chrono::high_resolution_clock::now(); + + rom_problem_->InterpolateArAndBr(*rom_options.new_point, Ar_interp, Br_interp); + + k_eff_ = rom_problem_->SolveROM(Ar_interp, Br_interp); + + auto end = std::chrono::high_resolution_clock::now(); + std::chrono::duration elapsed = end - start; + if (opensn::mpi_comm.rank() == 0) + { + std::ofstream outfile("results/online_time_" + std::to_string(rom_options.param_id) + ".txt"); + if (outfile.is_open()) + { + outfile << elapsed.count() << std::endl; + outfile.close(); + } + } + + log.Log() << "\n"; + log.Log() << " Final k-eigenvalue : " << std::setprecision(7) << k_eff_; + log.Log() << "\n\n"; + + log.Log() << "LinearBoltzmann::KEigenvalueROMSolver execution completed\n\n"; + } +} + +} // namespace opensn diff --git a/ROM/src/rom/pi_keigen_rom_solver.h b/ROM/src/rom/pi_keigen_rom_solver.h new file mode 100644 index 0000000..4f1a9e4 --- /dev/null +++ b/ROM/src/rom/pi_keigen_rom_solver.h @@ -0,0 +1,31 @@ +// SPDX-FileCopyrightText: 2026 The OpenSn Authors +// SPDX-License-Identifier: MIT + +#pragma once + +#include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/discrete_ordinates_problem.h" +#include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/solvers/pi_keigen_solver.h" +#include "rom_problem.h" + +namespace opensn +{ + +class LBSProblem; + +class PowerIterationKEigenROMSolver : public PowerIterationKEigenSolver +{ +protected: + std::shared_ptr lbs_problem_; + std::shared_ptr rom_problem_; + +public: + explicit PowerIterationKEigenROMSolver(const InputParameters& params); + + void Initialize() override; + + void Execute() override; + + static InputParameters GetInputParameters(); +}; + +} // namespace opensn diff --git a/ROM/src/rom/rom_problem.cc b/ROM/src/rom/rom_problem.cc index b961925..8ee987a 100644 --- a/ROM/src/rom/rom_problem.cc +++ b/ROM/src/rom/rom_problem.cc @@ -207,6 +207,27 @@ ROMProblem::ReadParamMatrix(const std::string& filename) } } +/** Loads the group-wise reduced bases from disk. + * + * Each group basis is read from the corresponding libROM basis file and + * stored in \c Ugs_. The reduced dimension is inferred from the first + * basis and stored in \c rom_rank. + */ +void +ROMProblem::LoadUgs() +{ + auto num_groups = lbs_problem_->GetNumGroups(); + Ugs_.reserve(num_groups); + for (auto g = 0; g < num_groups; ++g) + { + const auto basis_root = "basis/basis_" + std::to_string(g); + auto reader = std::make_unique(basis_root); + auto Ug = reader->getSpatialBasis(); + if (g == 0) rom_rank = Ug->numColumns(); + Ugs_.push_back(std::move(Ug)); + } +} + /** Assembles the global (local-DOF) matrix AU used to form reduced systems. * * For each group \c g and basis vector \c r, this method: @@ -315,6 +336,65 @@ ROMProblem::AssembleRHS() return b; } +/** Assembles the full-order operator image BU for k-eigenvalue ROM systems. + * + * For each group \c g and basis vector \c r, this method injects the basis + * column into the full-order dof layout, applies the inverse transport + * operator with fission, scattering, and fixed sources enabled, and stores + * the resulting vector as a column of BU. + * + * \return Shared pointer to BU of size (local_dofs) x (rom_rank * num_groups). + */ +std::shared_ptr +ROMProblem::AssembleBU() +{ + auto num_moments = lbs_problem_->GetNumMoments(); + auto num_groups = lbs_problem_->GetNumGroups(); + auto num_local_nodes = lbs_problem_->GetLocalNodeCount(); + const auto num_local_dofs = num_local_nodes * num_moments * num_groups; + + auto BU = std::make_shared(num_local_dofs, rom_rank * num_groups, true); + + // Assuming one groupset for ROM problems + assert(lbs_problem_->GetNumWGSSolvers() == 1); + auto raw_context = lbs_problem_->GetWGSSolver(0)->GetContext(); + auto gs_context = std::dynamic_pointer_cast(raw_context); + const auto scope = APPLY_AGS_FISSION_SOURCES | APPLY_WGS_FISSION_SOURCES | APPLY_AGS_SCATTER_SOURCES | APPLY_FIXED_SOURCES; + + auto& phi_old_local = lbs_problem_->GetPhiOldLocal(); + auto& q_moments_local = lbs_problem_->GetQMomentsLocal(); + + for (auto g = 0; g < num_groups; ++g) + { + for (auto r = 0; r < rom_rank; ++r) + { + std::vector basis_local(num_local_dofs, 0.0); + phi_old_local.assign(phi_old_local.size(), 0.0); + + auto col_g = Ugs_[g]->getColumn(r); + size_t rowg = 0; + for (size_t n = 0; n < num_local_nodes; ++n) + for (size_t m = 0; m < static_cast(num_moments); ++m, ++rowg) + { + const size_t row_phi = n * (num_moments * num_groups) + m * num_groups + g; + basis_local[row_phi] = col_g->item(rowg); + } + + q_moments_local.assign(q_moments_local.size(), 0.0); + gs_context->set_source_function(gs_context->groupset, q_moments_local, basis_local, scope); + + // Sweep + gs_context->ApplyInverseTransportOperator(scope); + std::vector phi_new_local = lbs_problem_->GetPhiNewLocal(); + + const auto col_idx = g * rom_rank + r; + for (size_t i = 0; i < static_cast(num_local_dofs); ++i) + BU->item(i, static_cast(col_idx)) = phi_new_local[i]; + } + } + return BU; +} + /** Forms and writes the reduced system (Ar, rhs) for a given AU and b. * * Computes: @@ -345,6 +425,35 @@ ROMProblem::AssembleROM( rhs->write(rhs_filename); } +/** Forms and writes the reduced matrices (Ar, Br) for a k-eigenvalue ROM. + * + * Computes: + * - Br = AU^T * BU + * - Ar = AU^T * AU + * and writes both matrices to libROM files. + * + * \param AU Full-order operator matrix assembled by AssembleAU(). + * \param BU Full-order operator image assembled by AssembleBU(). + * \param Ar_filename Output filename for Ar. + * \param Br_filename Output filename for Br. + */ +void +ROMProblem::AssembleROM( + std::shared_ptr& AU, + std::shared_ptr& BU, + const std::string& Ar_filename, + const std::string& Br_filename) +{ + // Br = AU^T * BU + auto Br = AU->transposeMult(*BU); + + // Ar = AU^T * AU + auto Ar = AU->transposeMult(*AU); + + // Save + Ar->write(Ar_filename); + Br->write(Br_filename); +} /** Solves the reduced system and reconstructs the full-order flux moments. * @@ -371,17 +480,6 @@ ROMProblem::SolveROM( auto num_local_nodes = lbs_problem_->GetLocalNodeCount(); auto num_local_dofs = num_local_nodes * num_moments * num_groups; - std::vector> Ugs; - Ugs.reserve(num_groups); - for (int g = 0; g < num_groups; ++g) - { - auto basis_root = "basis/basis_" + std::to_string(g); - auto reader_ptr = std::make_unique(basis_root); - auto Ug = reader_ptr->getSpatialBasis(); - if (g == 0) rom_rank = Ug->numColumns(); - Ugs.push_back(std::move(Ug)); - } - auto& phi_new_local = lbs_problem_->GetPhiNewLocal(); phi_new_local.assign(phi_new_local.size(), 0.0); @@ -392,7 +490,7 @@ ROMProblem::SolveROM( const int cr_idx = g * rom_rank + r; const double cr = (*c_vec)(cr_idx); - auto col_g = Ugs[g]->getColumn(r); + auto col_g = Ugs_[g]->getColumn(r); size_t row_g = 0; for (size_t n = 0; n < num_local_nodes; ++n) for (size_t m = 0; m < static_cast(num_moments); ++m, ++row_g) @@ -404,37 +502,100 @@ ROMProblem::SolveROM( } } -/** Initializes libROM interpolators for Ar and rhs over parameter space. + +/** Solves the reduced k-eigenvalue problem and reconstructs the full-order state. + * + * Forms the reduced operator Ar^{-1} Br, computes its right eigenpairs, + * selects the dominant positive real eigenvalue, and reconstructs the + * associated full-order local flux moments using the stored group-wise bases. + * + * \param Ar Reduced loss/operator matrix. + * \param Br Reduced production/operator matrix. + * \return The dominant positive real eigenvalue. + */ +double +ROMProblem::SolveROM( + std::shared_ptr& Ar, + std::shared_ptr& Br) +{ + auto Ar_inv = std::make_shared(Ar->numRows(), Ar->numColumns(), false); + auto Ar_inv_Br = std::make_shared(Ar->numRows(), Ar->numColumns(), false); + + Ar->inverse(*Ar_inv); + + Ar_inv_Br = Ar_inv->mult(*Br); + + auto eigen_pair = CAROM::NonSymmetricRightEigenSolve(*Ar_inv_Br); + + double k_eff = 0.0; + int best_col = -1; + + for (int i = 0; i < (int)eigen_pair.eigs.size(); ++i) + { + const auto& lam = eigen_pair.eigs[i]; + if (std::abs(lam.imag()) > 1.0e-10) continue; + if (lam.real() <= 0.0) continue; + if (lam.real() > k_eff) + { + k_eff = lam.real(); + best_col = i; + } + } + + auto num_moments = lbs_problem_->GetNumMoments(); + auto num_groups = lbs_problem_->GetNumGroups(); + auto num_local_nodes = lbs_problem_->GetLocalNodeCount(); + auto num_local_dofs = num_local_nodes * num_moments * num_groups; + + auto& phi_new_local = lbs_problem_->GetPhiNewLocal(); + phi_new_local.assign(phi_new_local.size(), 0.0); + + for (int g = 0; g < num_groups; ++g) + { + for (int r = 0; r < rom_rank; ++r) + { + const int cr_idx = g * rom_rank + r; + const double cr = eigen_pair.ev_real->item(cr_idx, best_col); + + auto col_g = Ugs_[g]->getColumn(r); + size_t row_g = 0; + for (size_t n = 0; n < (size_t)num_local_nodes; ++n) + for (size_t m = 0; m < (size_t)num_moments; ++m, ++row_g) + { + const size_t row_phi = + n * ((size_t)num_moments * (size_t)num_groups) + m * (size_t)num_groups + (size_t)g; + phi_new_local[row_phi] += cr * col_g->item(row_g); + } + } + } + return k_eff; +} + +/** Loads reduced matrices Ar and initializes the libROM matrix interpolator. * - * Loads all stored reduced systems from disk (one per parameter point), creates - * identity rotations, chooses a reference index (closest point), and constructs: - * - MatrixInterpolator on the SPD manifold for Ar - * - VectorInterpolator for rhs + * One reduced matrix is read for each sampled parameter point. Identity + * rotations are constructed, the closest sampled point to \p desired_point + * is used as the reference index, and the interpolator is initialized for + * subsequent online interpolation of Ar. * - * \param desired_point Parameter at which interpolation will be queried. + * \param desired_point Online parameter point used to choose the reference sample. */ void -ROMProblem::SetupInterpolator(CAROM::Vector& desired_point) +ROMProblem::SetupArInterpolator(CAROM::Vector& desired_point) { std::vector> Ar_matrices; - std::vector> rhs_vectors; // Load Ar and rhs from libROM files for (size_t i = 0; i < param_points.size(); ++i) { const std::string Ar_filename = "data/rom_system_Ar_" + std::to_string(i); - const std::string rhs_filename = "data/rom_system_br_" + std::to_string(i); // Create empty containers auto Ar = std::make_shared(); - auto rhs = std::make_shared(); - // Read matrix and vector + // Read matrix Ar->local_read(Ar_filename, opensn::mpi_comm.rank()); - rhs->local_read(rhs_filename, opensn::mpi_comm.rank()); - Ar_matrices.push_back(Ar); - rhs_vectors.push_back(rhs); } // Make Identity Rotations @@ -455,14 +616,105 @@ ROMProblem::SetupInterpolator(CAROM::Vector& desired_point) Ar_interp_obj_ptr_ = std::make_unique( param_points, rotations, Ar_matrices, ref_index, "SPD", "G", "LS", 0.999, false); +} + +/** Loads reduced RHS vectors and initializes the libROM vector interpolator. + * + * One reduced RHS vector is read for each sampled parameter point. Identity + * rotations are constructed, the closest sampled point to \p desired_point + * is used as the reference index, and the interpolator is initialized for + * subsequent online interpolation of the reduced RHS. + * + * \param desired_point Online parameter point used to choose the reference sample. + */ +void +ROMProblem::SetupRHSrInterpolator(CAROM::Vector& desired_point) +{ + std::vector> rhs_vectors; + + // Load Ar and rhs from libROM files + for (size_t i = 0; i < param_points.size(); ++i) + { + const std::string rhs_filename = "data/rom_system_br_" + std::to_string(i); + + // Create empty containers + auto rhs = std::make_shared(); + + // Read vector + rhs->local_read(rhs_filename, opensn::mpi_comm.rank()); + rhs_vectors.push_back(rhs); + } + + // Make Identity Rotations + std::vector> rotations; + int rom_dim = rhs_vectors[0]->dim(); + + for (size_t i = 0; i < rhs_vectors.size(); ++i) + { + std::shared_ptr I; + I = std::make_shared(rom_dim, rom_dim, false); + for (int j = 0; j < rom_dim; ++j) + I->item(j, j) = 1.0; + rotations.push_back(I); + } + + int ref_index = getClosestPoint(param_points, desired_point); + rhs_interp_obj_ptr_ = std::make_unique( param_points, rotations, rhs_vectors, ref_index, "G", "LS", 0.999, false); } +/** Loads reduced matrices Br and initializes the libROM matrix interpolator. + * + * One reduced matrix is read for each sampled parameter point. Identity + * rotations are constructed, the closest sampled point to \p desired_point + * is used as the reference index, and the interpolator is initialized for + * subsequent online interpolation of Br. + * + * \param desired_point Online parameter point used to choose the reference sample. + */ +void +ROMProblem::SetupBrInterpolator(CAROM::Vector& desired_point) +{ + std::vector> Br_matrices; + + // Load Ar and rhs from libROM files + for (size_t i = 0; i < param_points.size(); ++i) + { + const std::string Br_filename = "data/rom_system_Br_" + std::to_string(i); + + // Create empty containers + auto Br = std::make_shared(); + + // Read matrix + Br->local_read(Br_filename, opensn::mpi_comm.rank()); + Br_matrices.push_back(Br); + } + + // Make Identity Rotations + std::vector> rotations; + int rom_dim = Br_matrices[0]->numRows(); + + for (size_t i = 0; i < Br_matrices.size(); ++i) + { + std::shared_ptr I; + I = std::make_shared(rom_dim, rom_dim, false); + for (int j = 0; j < rom_dim; ++j) + I->item(j, j) = 1.0; + rotations.push_back(I); + } + + int ref_index = getClosestPoint(param_points, desired_point); + + Br_interp_obj_ptr_ = std::make_unique( + param_points, rotations, Br_matrices, + ref_index, "R", "G", "LS", 0.999, false); +} + /** Interpolates Ar and rhs at a desired parameter point. * - * Requires SetupInterpolator() to have been called. + * Requires SetupArInterpolator() and SetupRHSrInterpolator() to have been called. * * \param desired_point Parameter at which to interpolate. * \param Ar_interp Output interpolated reduced matrix. @@ -478,6 +730,23 @@ ROMProblem::InterpolateArAndRHSr( rhs_interp = rhs_interp_obj_ptr_->interpolate(desired_point); } +/** Interpolates Ar and Br at a desired parameter point. + * + * Requires SetupArInterpolator() and SetupBrInterpolator() to have been called. + * + * \param desired_point Parameter at which to interpolate. + * \param Ar_interp Output interpolated reduced matrix. + * \param Br_interp Output interpolated reduced production matrix. + */ +void +ROMProblem::InterpolateArAndBr( + CAROM::Vector& desired_point, + std::shared_ptr& Ar_interp, + std::shared_ptr& Br_interp) +{ + Ar_interp = Ar_interp_obj_ptr_->interpolate(desired_point); + Br_interp = Br_interp_obj_ptr_->interpolate(desired_point); +} /** Returns the schema for the nested `options` parameter block. * diff --git a/ROM/src/rom/rom_problem.h b/ROM/src/rom/rom_problem.h index 1bbd2fe..44464cc 100644 --- a/ROM/src/rom/rom_problem.h +++ b/ROM/src/rom/rom_problem.h @@ -40,11 +40,16 @@ class ROMProblem : public Problem /// Load the params from file void ReadParamMatrix(const std::string& filename); + void LoadUgs(); + /// Calculate AU via sweeps std::shared_ptr AssembleAU(); /// Sweep to form RHS std::shared_ptr AssembleRHS(); + + /// Sweep to form BU + std::shared_ptr AssembleBU(); /// Assemble the reduced system and save to file void AssembleROM(std::shared_ptr& AU, @@ -52,23 +57,47 @@ class ROMProblem : public Problem const std::string& Ar_filename, const std::string& rhs_filename); + /// Assemble the reduced system and save to file + void AssembleROM(std::shared_ptr& AU, + std::shared_ptr& BU, + const std::string& Ar_filename, + const std::string& Br_filename); + /// Solve given LHS and RHS of a ROM system void SolveROM(std::shared_ptr& Ar, std::shared_ptr& rhs); - /// Load reduced systems and initialize libROM interpolator objects - void SetupInterpolator(CAROM::Vector& desired_point); + /// Solve given LHS and RHS of a k-eigenvalue ROM system + double SolveROM(std::shared_ptr& Ar, + std::shared_ptr& Br); + + /// Load Ar and initialize libROM interpolator objects + void SetupArInterpolator(CAROM::Vector& desired_point); + + /// Load br and initialize libROM interpolator objects + void SetupRHSrInterpolator(CAROM::Vector& desired_point); + + /// Load Br and initialize libROM interpolator objects + void SetupBrInterpolator(CAROM::Vector& desired_point); void InterpolateArAndRHSr( CAROM::Vector& desired_point, std::shared_ptr& Ar_interp, std::shared_ptr& rhs_interp); + void InterpolateArAndBr( + CAROM::Vector& desired_point, + std::shared_ptr& Ar_interp, + std::shared_ptr& Br_interp); + protected: std::unique_ptr spatial_basis_; opensn::Vector b_; + std::vector> Ugs_; std::unique_ptr Ar_interp_obj_ptr_; + std::unique_ptr Br_interp_obj_ptr_; std::unique_ptr rhs_interp_obj_ptr_; + std::shared_ptr lbs_problem_; ROMOptions options_; diff --git a/ROM/src/rom/steady_state_rom_solver.cc b/ROM/src/rom/steady_state_rom_solver.cc index f0d3e21..05c66cb 100644 --- a/ROM/src/rom/steady_state_rom_solver.cc +++ b/ROM/src/rom/steady_state_rom_solver.cc @@ -13,6 +13,12 @@ namespace opensn { + +/** Returns the input-parameter schema for SteadyStateROMSolver. + * + * Extends the steady-state source solver schema with: + * - `rom_problem` : an existing ROMProblem instance that manages ROM workflow. + */ InputParameters SteadyStateROMSolver::GetInputParameters() { @@ -38,6 +44,14 @@ SteadyStateROMSolver::Initialize() { } +/** Executes the requested ROM workflow phase for the steady-state solver. + * + * Supported phases are: + * - OFFLINE : solves the full-order system and writes snapshots, + * - MERGE : builds the reduced bases from stored snapshots, + * - SYSTEMS : assembles and writes reduced operators, + * - ONLINE : interpolates and solves the reduced system. + */ void SteadyStateROMSolver::Execute() { @@ -70,6 +84,8 @@ SteadyStateROMSolver::Execute() } if (rom_options.phase == Phase::SYSTEMS) { + rom_problem_->LoadUgs(); + std::shared_ptr AU_ = rom_problem_->AssembleAU(); std::shared_ptr b_ = rom_problem_->AssembleRHS(); const std::string& Ar_filename = "data/rom_system_Ar_" + std::to_string(rom_options.param_id); @@ -79,10 +95,12 @@ SteadyStateROMSolver::Execute() if (rom_options.phase == Phase::ONLINE) { rom_problem_->ReadParamMatrix(rom_options.param_file); + rom_problem_->LoadUgs(); std::shared_ptr Ar_interp; std::shared_ptr rhs_interp; - rom_problem_->SetupInterpolator(*rom_options.new_point); + rom_problem_->SetupArInterpolator(*rom_options.new_point); + rom_problem_->SetupRHSrInterpolator(*rom_options.new_point); auto start = std::chrono::high_resolution_clock::now();