ROM k-eigenvalue implementation

.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: |

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=<path-to-opensn-install> \
+  -DlibROM_DIR=<path-to-librom-install>
+make -j
+```
+An example can be run from its respective folder with
+```bash
+python run_rom_*.py --exe=path/to/app/exe
+```

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

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)

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

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

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")

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()