From 4b23c799ffd1a5d33bd8262b5a7d4aab6f12ec01 Mon Sep 17 00:00:00 2001 From: quincy-huhn98 Date: Thu, 15 May 2025 17:31:32 -0500 Subject: [PATCH 1/4] Test App Initial Commit --- .gitignore | 89 +++++++++++++++++++++++++++++++++++++++++++++++ CMakeLists.txt | 54 ++++++++++++++++++++++++++++ src/main.cc | 11 ++++++ src/test.cc | 28 +++++++++++++++ src/test.h | 15 ++++++++ src/testapp_py.cc | 12 +++++++ test.py | 8 +++++ 7 files changed, 217 insertions(+) create mode 100644 .gitignore create mode 100644 CMakeLists.txt create mode 100644 src/main.cc create mode 100644 src/test.cc create mode 100644 src/test.h create mode 100644 src/testapp_py.cc create mode 100644 test.py diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..dc23fda --- /dev/null +++ b/.gitignore @@ -0,0 +1,89 @@ +.git +cmake-build-debug/ +cmake-build-debug +.idea +bin +build*/ +doc/documentation +doc/HTMLdocs +doc/whitepages + +# root Makefile is generated by CMake automatically +Makefile + +resources/Dependencies/lua-5.3.5/ +resources/Dependencies/ncurses/ +resources/Dependencies/readline/ +resources/Dependencies/triangle/ +resources/Dependencies/glew-2.1.0/ + +# Latex compile files +*.aux +*.lof +*.log +*.lot +*.synctex.gz +*.toc +*.pdf + +# All vtk mesh files +/*.vtu +/*.pvtu +test/**/*.pvtu +test/**/*.vtu +test/modules/linear_boltzmann_solvers/dsa/SimTest_92b_DSA_PWLC.pvtu +test/modules/linear_boltzmann_solvers/dsa/SimTest_92b_DSA_PWLC_0.vtu +test/**/*.csv + +tests/BigTests/*/solutions/ + +# All exodus files +/*.e + +resources/Dependencies/.DS_Store + +resources/.DS_Store + +# python files +.DS_Store +._.DS_Store +*.pyc +*.pmesh +__pycache__ +**/__pycache__ + +*-private.sh + +#Documentation +doc/generated_files + +#visual studio code files +.vscode/ + +test/**/out/ +tutorials/**/out/ + +#Scratch directory +scratch/ + +#general data files +*.data +tests/BigTests/c5g7/Test1/solutions/ +**/*.cmesh + +# Python generated libraries +pyopensn/libopensn*.so* +pyopensn/libopensn*.dylib +pyopensn/__init__.cpython* +pyopensn.egg-info + +# Files generated by ply package +doc/scripts/parser.out +doc/scripts/parsetab.py + +doc/**/*.pvtu +doc/**/*.vtu + +tutorials/**/*.pvtu +tutorials/**/*.vtu +tutorials/**/*.py \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..8c873af --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,54 @@ +cmake_minimum_required(VERSION 3.14) +project(MyOpenSnApp LANGUAGES C CXX) + +set(CMAKE_CXX_STANDARD 17) + +find_package(MPI REQUIRED) +find_package(VTK REQUIRED COMPONENTS + CommonCore + CommonDataModel + IOLegacy + IOCore + IOXML + ParallelCore + IOParallelXML + FiltersCore + IOEnSight + IOExodus +) +find_package(OpenSn REQUIRED) + +# Optional: find pybind11 +find_package(pybind11 REQUIRED) + +find_package(caliper REQUIRED) + +find_package(Python REQUIRED COMPONENTS Interpreter Development) + +# Python module +add_library(testapp MODULE + src/test.cc + src/testapp_py.cc +) + +execute_process( + COMMAND python3 -m pybind11 --includes + OUTPUT_VARIABLE PYBIND11_INCLUDE_FLAGS + OUTPUT_STRIP_TRAILING_WHITESPACE +) +string(REPLACE "-I" "" PYBIND11_INCLUDE_DIRS "${PYBIND11_INCLUDE_FLAGS}") +separate_arguments(PYBIND11_INCLUDE_DIRS) + +target_include_directories(testapp PRIVATE ${OPENSN_INCLUDE_DIR} ${PYBIND11_INCLUDE_DIRS}) +target_link_libraries(testapp PRIVATE opensn::libopensn MPI::MPI_C caliper pybind11::module Python::Python) + +set_target_properties(testapp PROPERTIES + PREFIX "" # remove "lib" prefix + OUTPUT_NAME "testapp" # Python import name + LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/python +) + +# Optional executable +add_executable(test_app_exec src/main.cc src/test.cc) +target_include_directories(test_app_exec PRIVATE ${OPENSN_INCLUDE_DIR} ${PYBIND11_INCLUDE_DIRS}) +target_link_libraries(test_app_exec PRIVATE opensn::libopensn MPI::MPI_C caliper pybind11::module Python::Python) \ No newline at end of file diff --git a/src/main.cc b/src/main.cc new file mode 100644 index 0000000..7c59a53 --- /dev/null +++ b/src/main.cc @@ -0,0 +1,11 @@ +#include "test.h" +#include "opensn/python/lib/py_app.h" +#include "opensn/mpicpp-lite/mpicpp-lite.h" + +int main(int argc, char** argv) +{ + mpi::Environment env(argc, argv); + py::scoped_interpreter guard{}; + auto app = TestApp::Create({}); + return 0; +} diff --git a/src/test.cc b/src/test.cc new file mode 100644 index 0000000..b25f81c --- /dev/null +++ b/src/test.cc @@ -0,0 +1,28 @@ +#include "test.h" +#include "opensn/framework/logging/log.h" +#include "opensn/framework/object_factory.h" + +using namespace opensn; + +OpenSnRegisterObject(TestApp); + +InputParameters +TestApp::GetInputParameters() +{ + InputParameters params; + params.SetGeneralDescription("A dummy application for testing."); + params.AddOptionalParameter("name", "test_app", "A name for the application."); + return params; +} + +TestApp::TestApp(const InputParameters& params) : name_(params.GetParamValue("name")) +{ + opensn::log.Log() << "TestApp named " << name_ << " created."; +} + +std::shared_ptr +TestApp::Create(const ParameterBlock& params) +{ + auto& factory = opensn::ObjectFactory::GetInstance(); + return factory.Create("TestApp", params); +} diff --git a/src/test.h b/src/test.h new file mode 100644 index 0000000..2ea1caa --- /dev/null +++ b/src/test.h @@ -0,0 +1,15 @@ +#pragma once +#include "opensn/framework/object.h" + +class TestApp : public opensn::Object +{ +public: + explicit TestApp(const opensn::InputParameters& params); + +private: + const std::string name_; + +public: + static opensn::InputParameters GetInputParameters(); + static std::shared_ptr Create(const opensn::ParameterBlock& params); +}; diff --git a/src/testapp_py.cc b/src/testapp_py.cc new file mode 100644 index 0000000..e841cf1 --- /dev/null +++ b/src/testapp_py.cc @@ -0,0 +1,12 @@ +#include +#include "test.h" + +namespace py = pybind11; + +PYBIND11_MODULE(myapp, m) +{ + m.doc() = "Python bindings for my OpenSn-based TestApp"; + + py::class_>(m, "TestApp") + .def_static("Create", &TestApp::Create); +} diff --git a/test.py b/test.py new file mode 100644 index 0000000..693f29e --- /dev/null +++ b/test.py @@ -0,0 +1,8 @@ +# test.py +import sys +sys.path.insert(0, "build/python") # where CMake drops the .so + +import myapp + +app = myapp.TestApp.Create() +print("Successfully created:", app) From c935fc8e5e3ef2245c479caaacb78b2f6e0c6f32 Mon Sep 17 00:00:00 2001 From: quincy-huhn98 Date: Wed, 21 May 2025 14:46:14 -0500 Subject: [PATCH 2/4] Actions and Test --- .github/workflows/regression.yaml | 42 +++++++++++++++++++++++++++++++ test.py => tests/test.py | 0 2 files changed, 42 insertions(+) create mode 100644 .github/workflows/regression.yaml rename test.py => tests/test.py (100%) diff --git a/.github/workflows/regression.yaml b/.github/workflows/regression.yaml new file mode 100644 index 0000000..6542198 --- /dev/null +++ b/.github/workflows/regression.yaml @@ -0,0 +1,42 @@ +name: Regression Tests + +on: + schedule: + - cron: "0 0 * * *" + push: + branches: [main] + workflow_dispatch: + pull_request: + branches: + - main + +jobs: + run-tests: + runs-on: [self-hosted] + strategy: + fail-fast: false + steps: + - name: checkout testapp + uses: actions/checkout@v4 + - name: checkout opensn + uses: actions/checkout@v4 + with: + repository: Open-Sn/opensn + path: opensn + - name: install opensn + shell: bash + run: | + module load opensn/clang/17 python3/3.12.3 + cd opensn && mkdir build && mkdir install && cd build + cmake -DOPENSN_WITH_PYTHON=True -DCMAKE_INSTALL_PREFIX=../install .. && make -j && make install + - name: compile app + shell: bash + run: | + module load opensn/clang/17 python3/3.12.3 + cd OpenSnApp && mkdir build && cd build + cmake -DCMAKE_PREFIX_PATH=../../opensn/install .. && make -j + - name: test examples + shell: bash + run: | + module load opensn/clang/17 python3/3.12.3 + cd tests && ../build/test_app_exec -i test.py \ No newline at end of file diff --git a/test.py b/tests/test.py similarity index 100% rename from test.py rename to tests/test.py From 47e7a6e4c5cf659d7e02b27b4a0973671eeb1cab Mon Sep 17 00:00:00 2001 From: quincy-huhn98 Date: Fri, 31 Oct 2025 16:02:53 -0500 Subject: [PATCH 3/4] ROM App --- .github/workflows/regression.yaml | 4 +- .gitignore | 7 +- CMakeLists.txt | 87 +++- examples/2gcheckerboard/absorber_base.txt | 14 + .../2gcheckerboard/base_2gcheckerboard.py | 192 ++++++++ .../2gcheckerboard/run_rom_2gcheckerboard.py | 105 +++++ examples/2gcheckerboard/scatterer_base.txt | 14 + examples/checkerboard/base_checkerboard.py | 198 ++++++++ examples/checkerboard/run_rom_checkerboard.py | 92 ++++ examples/python/plot_errors.py | 84 ++++ examples/python/plotting.py | 116 +++++ examples/python/utils.py | 125 +++++ examples/reed/base_reed.py | 158 +++++++ examples/reed/run_rom_reed.py | 70 +++ src/main.cc | 30 +- src/py_wrappers.h | 21 + src/rom.cc | 132 ++++++ src/rom/rom_problem.cc | 439 ++++++++++++++++++ src/rom/rom_problem.h | 80 ++++ src/rom/rom_structs.h | 21 + src/rom/steady_state_rom_solver.cc | 111 +++++ src/rom/steady_state_rom_solver.h | 32 ++ src/rom_py_app.cc | 19 + src/rom_py_app.h | 20 + src/test.cc | 28 -- src/test.h | 15 - src/testapp_py.cc | 12 - tests/test.py | 8 - 28 files changed, 2147 insertions(+), 87 deletions(-) create mode 100644 examples/2gcheckerboard/absorber_base.txt create mode 100644 examples/2gcheckerboard/base_2gcheckerboard.py create mode 100644 examples/2gcheckerboard/run_rom_2gcheckerboard.py create mode 100644 examples/2gcheckerboard/scatterer_base.txt create mode 100644 examples/checkerboard/base_checkerboard.py create mode 100644 examples/checkerboard/run_rom_checkerboard.py create mode 100644 examples/python/plot_errors.py create mode 100644 examples/python/plotting.py create mode 100644 examples/python/utils.py create mode 100644 examples/reed/base_reed.py create mode 100644 examples/reed/run_rom_reed.py create mode 100644 src/py_wrappers.h create mode 100644 src/rom.cc create mode 100644 src/rom/rom_problem.cc create mode 100644 src/rom/rom_problem.h create mode 100644 src/rom/rom_structs.h create mode 100644 src/rom/steady_state_rom_solver.cc create mode 100644 src/rom/steady_state_rom_solver.h create mode 100644 src/rom_py_app.cc create mode 100644 src/rom_py_app.h delete mode 100644 src/test.cc delete mode 100644 src/test.h delete mode 100644 src/testapp_py.cc delete mode 100644 tests/test.py diff --git a/.github/workflows/regression.yaml b/.github/workflows/regression.yaml index 6542198..44076c0 100644 --- a/.github/workflows/regression.yaml +++ b/.github/workflows/regression.yaml @@ -28,7 +28,7 @@ jobs: run: | module load opensn/clang/17 python3/3.12.3 cd opensn && mkdir build && mkdir install && cd build - cmake -DOPENSN_WITH_PYTHON=True -DCMAKE_INSTALL_PREFIX=../install .. && make -j && make install + cmake -DCMAKE_INSTALL_PREFIX=../install .. && make -j && make install - name: compile app shell: bash run: | @@ -39,4 +39,4 @@ jobs: shell: bash run: | module load opensn/clang/17 python3/3.12.3 - cd tests && ../build/test_app_exec -i test.py \ No newline at end of file + cd examples/reed && python run_rom_reed.py \ No newline at end of file diff --git a/.gitignore b/.gitignore index dc23fda..dfe896d 100644 --- a/.gitignore +++ b/.gitignore @@ -86,4 +86,9 @@ doc/**/*.vtu tutorials/**/*.pvtu tutorials/**/*.vtu -tutorials/**/*.py \ No newline at end of file +tutorials/**/*.py + +examples/*/basis/* +examples/*/data/* +examples/*/output/* +examples/*/results/* \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt index 8c873af..f827cf6 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,7 +1,9 @@ cmake_minimum_required(VERSION 3.14) -project(MyOpenSnApp LANGUAGES C CXX) +project(OpenSnApp LANGUAGES C CXX) -set(CMAKE_CXX_STANDARD 17) +set(CMAKE_CXX_STANDARD 20) +set(CMAKE_CXX_STANDARD_REQUIRED ON) +set(CMAKE_CXX_EXTENSIONS OFF) find_package(MPI REQUIRED) find_package(VTK REQUIRED COMPONENTS @@ -16,19 +18,48 @@ find_package(VTK REQUIRED COMPONENTS IOEnSight IOExodus ) + +include(CMakeFindDependencyMacro) + find_package(OpenSn REQUIRED) -# Optional: find pybind11 +message(STATUS "Found OpenSn ${OpenSn_VERSION} (from: ${OpenSn_DIR})") + find_package(pybind11 REQUIRED) find_package(caliper REQUIRED) find_package(Python REQUIRED COMPONENTS Interpreter Development) +find_library(LIBROM_LIBRARY + NAMES ROM + HINTS + ENV LIBROM_DIR + ${LIBROM_DIR} + ${LIBROM_DIR}/install + PATH_SUFFIXES lib + REQUIRED) + +find_path(LIBROM_INCLUDE_DIR + NAMES librom.h + HINTS + ENV LIBROM_DIR + ${LIBROM_DIR} + ${LIBROM_DIR}/install + PATH_SUFFIXES include + REQUIRED) + +message(STATUS "libROM include dir: ${LIBROM_INCLUDE_DIR}") +message(STATUS "libROM library: ${LIBROM_LIBRARY}") + +file(GLOB_RECURSE ROM_SRCS CONFIGURE_DEPENDS ${PROJECT_SOURCE_DIR}/src/rom/*.cc) + # Python module -add_library(testapp MODULE - src/test.cc - src/testapp_py.cc +add_library(romapp MODULE + src/main.cc + src/rom_py_app.cc + src/rom.cc + ${ROM_SRCS} ) execute_process( @@ -39,16 +70,42 @@ execute_process( string(REPLACE "-I" "" PYBIND11_INCLUDE_DIRS "${PYBIND11_INCLUDE_FLAGS}") separate_arguments(PYBIND11_INCLUDE_DIRS) -target_include_directories(testapp PRIVATE ${OPENSN_INCLUDE_DIR} ${PYBIND11_INCLUDE_DIRS}) -target_link_libraries(testapp PRIVATE opensn::libopensn MPI::MPI_C caliper pybind11::module Python::Python) +target_include_directories(romapp + PRIVATE + ${PROJECT_SOURCE_DIR}/rom + ${PYBIND11_INCLUDE_DIRS} + ${OpenSn_INSTALL_PREFIX}/include/opensn + ${LIBROM_INCLUDE_DIR} +) +target_link_libraries(romapp + PRIVATE + opensn::opensn + opensn::libopensnpy + MPI::MPI_CXX + caliper + pybind11::module + Python::Python + ${LIBROM_LIBRARY}) -set_target_properties(testapp PROPERTIES - PREFIX "" # remove "lib" prefix - OUTPUT_NAME "testapp" # Python import name +set_target_properties(romapp PROPERTIES + PREFIX "" + OUTPUT_NAME "romapp" LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/python ) -# Optional executable -add_executable(test_app_exec src/main.cc src/test.cc) -target_include_directories(test_app_exec PRIVATE ${OPENSN_INCLUDE_DIR} ${PYBIND11_INCLUDE_DIRS}) -target_link_libraries(test_app_exec PRIVATE opensn::libopensn MPI::MPI_C caliper pybind11::module Python::Python) \ No newline at end of file +add_executable(rom_app_exec src/main.cc src/rom_py_app.cc src/rom.cc ${ROM_SRCS}) +target_include_directories(rom_app_exec PRIVATE + ${PROJECT_SOURCE_DIR}/rom + ${PYBIND11_INCLUDE_DIRS} + ${OpenSn_INSTALL_PREFIX}/include/opensn + ${LIBROM_INCLUDE_DIR} +) +target_link_libraries(rom_app_exec + PRIVATE + opensn::opensn + opensn::libopensnpy + MPI::MPI_C + caliper + pybind11::module + Python::Python + ${LIBROM_LIBRARY}) \ No newline at end of file diff --git a/examples/2gcheckerboard/absorber_base.txt b/examples/2gcheckerboard/absorber_base.txt new file mode 100644 index 0000000..8da84b2 --- /dev/null +++ b/examples/2gcheckerboard/absorber_base.txt @@ -0,0 +1,14 @@ +NUM_GROUPS 2 +NUM_MOMENTS 1 + +SIGMA_T_BEGIN +0 1.0 +1 1.0 +SIGMA_T_END + +TRANSFER_MOMENTS_BEGIN +M_GPRIME_G_VAL 0 0 0 0.05 +M_GPRIME_G_VAL 0 1 0 0.45 +M_GPRIME_G_VAL 0 0 1 0.45 +M_GPRIME_G_VAL 0 1 1 0.05 +TRANSFER_MOMENTS_END \ No newline at end of file diff --git a/examples/2gcheckerboard/base_2gcheckerboard.py b/examples/2gcheckerboard/base_2gcheckerboard.py new file mode 100644 index 0000000..00672ee --- /dev/null +++ b/examples/2gcheckerboard/base_2gcheckerboard.py @@ -0,0 +1,192 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +# 2D Transport test. Checkerboard https://doi.org/10.1016/j.jcp.2022.111525 + +import os +import sys +import math +import time +import numpy as np + +if "opensn_console" not in globals(): + from mpi4py import MPI + size = MPI.COMM_WORLD.size + rank = MPI.COMM_WORLD.rank + sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), "../../../../../"))) + from pyopensn.mesh import OrthogonalMeshGenerator + from pyopensn.xs import MultiGroupXS + from pyopensn.source import VolumetricSource + from pyopensn.aquad import GLCProductQuadrature2DXY + from pyopensn.solver import DiscreteOrdinatesProblem, SteadyStateSolver + from pyopensn.fieldfunc import FieldFunctionInterpolationVolume + from pyopensn.logvol import RPPLogicalVolume + +if __name__ == "__main__": + + try: + print("Absorber Sigma_a Parameter = {}".format(abs_1)) + param = abs_1 + except: + abs_1=10.0 + print("Absorber Sigma_a Nominal = {}".format(abs_1)) + + try: + print("Absorber Sigma_s Parameter = {}".format(scatt_1)) + param = scatt_1 + except: + scatt_1=10.0 + print("Absorber Sigma_s Nominal = {}".format(scatt_1)) + + + try: + print("Scatterer Sigma_a Parameter = {}".format(abs_2)) + param = abs_2 + except: + abs_2=0.0 + print("Scatterer Sigma_a Nominal = {}".format(abs_2)) + + try: + print("Scatterer Sigma_s Parameter = {}".format(scatt_2)) + param = scatt_2 + except: + scatt_2=1.0 + print("Scatterer Sigma_s Nominal = {}".format(scatt_2)) + + + try: + print("Source Parameter = {}".format(param_q)) + param = param_q + except: + param_q=1.0 + print("Source Nominal = {}".format(param_q)) + + try: + print("Parameter id = {}".format(p_id)) + except: + p_id=0 + print("Parameter id = {}".format(p_id)) + + try: + if phase == 0: + print("Offline Phase") + phase = "offline" + elif phase == 1: + print("Merge Phase") + phase = "merge" + elif phase == 2: + print("Systems Phase") + phase = "systems" + elif phase == 3: + print("Online Phase") + phase = "online" + except: + phase="offline" + print("Phase default to offline") + + # Check number of processors + num_procs = 4 + if size != num_procs: + sys.exit(f"Incorrect number of processors. Expected {num_procs} processors but got {size}.") + + # Setup mesh + nodes = [] + N = 70 + L = 7 + xmin = 0 + dx = L / N + for i in range(N + 1): + nodes.append(xmin + i * dx) + meshgen = OrthogonalMeshGenerator(node_sets=[nodes, nodes], + partitioner=KBAGraphPartitioner( + nx=2, + ny=2, + xcuts=[3.5], + ycuts=[3.5], + )) + grid = meshgen.Execute() + + # Set background (Scatterer) block ID = 0 + grid.SetUniformBlockID(0) + + # Set Source (central red square from x=3 to x=4, y=3 to y=4) block ID = 1 + logvol_src = RPPLogicalVolume(xmin=3.0, xmax=4.0, + ymin=3.0, ymax=4.0, + infz=True) + grid.SetBlockIDFromLogicalVolume(logvol_src, 1, True) + + # Set Absorbers (green 1x1 squares) block ID = 2 + absorber_centers = [ + (1,1), (3,1), (5,1), + (2,2), (4,2), + (1,3), (5,3), + (2,4), (4,4), + (1,5), (5,5) + ] + for xc, yc in absorber_centers: + vol_abs = RPPLogicalVolume( + xmin=xc+0.0, xmax=xc+1.0, + ymin=yc+0.0, ymax=yc+1.0, + infz=True + ) + grid.SetBlockIDFromLogicalVolume(vol_abs, 2, True) + + num_groups = 2 + scatterer = MultiGroupXS() + scatterer.LoadFromOpenSn("data/scatterer.xs") + + absorber = MultiGroupXS() + absorber.LoadFromOpenSn("data/absorber.xs") + + strength = [0.0 for _ in range(num_groups)] + src0 = VolumetricSource(block_ids=[0], group_strength=strength) + strength[0] = param_q + src1 = VolumetricSource(block_ids=[1], group_strength=strength) + + # Setup Physics + pquad = GLCProductQuadrature2DXY(n_polar=4, n_azimuthal=32, scattering_order=0) + + if phase == "online": + rom_options={ + "param_id":0, + "phase":phase, + "param_file":"data/params.txt", + "new_point":[scatt_1, abs_1] + } + else: + rom_options={ + "param_id":p_id, + "phase":phase + } + + phys = DiscreteOrdinatesProblem( + mesh=grid, + num_groups=num_groups, + groupsets=[ + { + "groups_from_to": [0, 1], + "angular_quadrature": pquad, + "angle_aggregation_num_subsets": 1, + "inner_linear_method": "petsc_gmres", + "l_abs_tol": 1.0e-8, + "l_max_its": 300, + "gmres_restart_interval": 100, + }, + ], + xs_map=[ + {"block_ids": [0, 1], "xs": scatterer}, + {"block_ids": [2], "xs": absorber} + ], + volumetric_sources=[src0, src1] + ) + + rom = ROMProblem(problem=phys,options=rom_options) + + ss_solver = SteadyStateROMSolver(problem=phys, rom_problem=rom) + ss_solver.Initialize() + ss_solver.Execute() + + if phase == "online": + phys.WriteFluxMoments("output/rom") + if phase == "offline": + phys.WriteFluxMoments("output/fom") \ No newline at end of file diff --git a/examples/2gcheckerboard/run_rom_2gcheckerboard.py b/examples/2gcheckerboard/run_rom_2gcheckerboard.py new file mode 100644 index 0000000..c702a6d --- /dev/null +++ b/examples/2gcheckerboard/run_rom_2gcheckerboard.py @@ -0,0 +1,105 @@ +import numpy as np +import sys, os +sys.path.insert(0, os.path.realpath("../python")) + +import plotting +import utils + +# Sampling training points +bounds = [[0.5,1.0],[7.5,12.5]] +num_params = 100 + +params = utils.sample_parameter_space(bounds, num_params) + +S_abs = [[0.0, 0.0], + [0.0, 0.0]] +sigma_t_scatt = [1.0, 1.0] + +# OFFLINE PHASE +phase = 0 + +for i in range(num_params): + S_scatt = [[1-params[i,0], params[i,0]], + [0.0, 1.0]] + utils.update_xs("scatterer_base.txt", "data/scatterer.xs", sigma_t_scatt, S_scatt) + + sigma_t_abs = [params[i,1], params[i,1]] + utils.update_xs("absorber_base.txt", "data/absorber.xs", sigma_t_abs, S_abs) + + cmd = "mpiexec -n=4 ../../build/rom_app_exec -i base_2gcheckerboard.py \ + -p phase={} -p p_id={}".format(phase,i) + utils.run_opensn(cmd) + +# MERGE PHASE +phase = 1 + +cmd = "mpiexec -n=4 ../../build/rom_app_exec -i base_2gcheckerboard.py -p phase={} -p p_id={}".format(phase, i) +utils.run_opensn(cmd) + +plotting.plot_sv(num_groups=2) + + +# SYSTEMS PHASE +phase = 2 +for i in range(num_params): + S_scatt = [[1-params[i,0], params[i,0]], + [0.0, 1.0]] + utils.update_xs("scatterer_base.txt", "data/scatterer.xs", sigma_t_scatt, S_scatt) + + sigma_t_abs = [params[i,1], params[i,1]] + utils.update_xs("absorber_base.txt", "data/absorber.xs", sigma_t_abs, S_abs) + + cmd = "mpiexec -n=4 ../../build/rom_app_exec -i base_2gcheckerboard.py \ + -p phase={} -p p_id={}".format(phase,i) + utils.run_opensn(cmd) + +np.savetxt("data/params.txt", params) + + +# Generate Test Data +test_scatt_1 = np.random.uniform(0.5,1.0,10) +test_abs_1 = np.random.uniform(7.5,12.5,10) +test = np.append(test_scatt_1[:,np.newaxis], test_abs_1[:,np.newaxis], axis=1) +np.savetxt("data/validation.txt", test) + +test = np.loadtxt("data/validation.txt") + +num_test = 10 +errors = [] +speedups = [] + +for i in range(num_test): + # ONLINE PHASE + S_scatt = [[1-test[i,0], test[i,0]], + [0.0, 1.0]] + utils.update_xs("scatterer_base.txt", "data/scatterer.xs", sigma_t_scatt, S_scatt) + + sigma_t_abs = [test[i,1], test[i,1]] + utils.update_xs("absorber_base.txt", "data/absorber.xs", sigma_t_abs, S_abs) + + phase = 3 + cmd = "mpiexec -n=4 ../../build/rom_app_exec -i base_2gcheckerboard.py \ + -p phase={} -p scatt_1={} -p abs_1={} -p p_id={}"\ + .format(phase,test[i][0],test[i][1],i) + utils.run_opensn(cmd) + rom_time = np.loadtxt("results/online_time.txt") + + # Reference FOM solution + phase = 0 + cmd = "mpiexec -n=4 ../../build/rom_app_exec -i base_2gcheckerboard.py \ + -p phase={} -p p_id={}".format(phase,i) + utils.run_opensn(cmd) + fom_time = np.loadtxt("results/offline_time.txt") + + plotting.plot_2d_flux("output/fom{}.h5", ranks=range(4), prefix="fom", pid=i) + plotting.plot_2d_flux("output/rom{}.h5", ranks=range(4), prefix="rom", pid=i) + + error = plotting.plot_2d_lineout(ranks=range(4), pid=i) + + errors.append(error) + speedups.append(fom_time/rom_time) + +print("Avg Error ", np.mean(errors)) +np.savetxt("results/errors.txt", errors) +print("Avg Speedup ", np.mean(speedups)) +np.savetxt("results/speedups.txt", speedups) \ No newline at end of file diff --git a/examples/2gcheckerboard/scatterer_base.txt b/examples/2gcheckerboard/scatterer_base.txt new file mode 100644 index 0000000..8da84b2 --- /dev/null +++ b/examples/2gcheckerboard/scatterer_base.txt @@ -0,0 +1,14 @@ +NUM_GROUPS 2 +NUM_MOMENTS 1 + +SIGMA_T_BEGIN +0 1.0 +1 1.0 +SIGMA_T_END + +TRANSFER_MOMENTS_BEGIN +M_GPRIME_G_VAL 0 0 0 0.05 +M_GPRIME_G_VAL 0 1 0 0.45 +M_GPRIME_G_VAL 0 0 1 0.45 +M_GPRIME_G_VAL 0 1 1 0.05 +TRANSFER_MOMENTS_END \ No newline at end of file diff --git a/examples/checkerboard/base_checkerboard.py b/examples/checkerboard/base_checkerboard.py new file mode 100644 index 0000000..6f34f7c --- /dev/null +++ b/examples/checkerboard/base_checkerboard.py @@ -0,0 +1,198 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +# 2D Transport test. Checkerboard https://doi.org/10.1016/j.jcp.2022.111525 + +import os +import sys +import math +import time +import numpy as np + +if "opensn_console" not in globals(): + from mpi4py import MPI + size = MPI.COMM_WORLD.size + rank = MPI.COMM_WORLD.rank + sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), "../../../../../"))) + from pyopensn.mesh import OrthogonalMeshGenerator + from pyopensn.xs import MultiGroupXS + from pyopensn.source import VolumetricSource + from pyopensn.aquad import GLCProductQuadrature2DXY + from pyopensn.solver import DiscreteOrdinatesProblem, SteadyStateSolver + from pyopensn.fieldfunc import FieldFunctionInterpolationVolume + from pyopensn.logvol import RPPLogicalVolume + +if __name__ == "__main__": + + try: + print("Absorber Sigma_a Parameter = {}".format(abs_1)) + param = abs_1 + except: + abs_1=10.0 + print("Absorber Sigma_a Nominal = {}".format(abs_1)) + + try: + print("Absorber Sigma_s Parameter = {}".format(scatt_1)) + param = scatt_1 + except: + scatt_1=10.0 + print("Absorber Sigma_s Nominal = {}".format(scatt_1)) + + + try: + print("Scatterer Sigma_a Parameter = {}".format(abs_2)) + param = abs_2 + except: + abs_2=0.0 + print("Scatterer Sigma_a Nominal = {}".format(abs_2)) + + try: + print("Scatterer Sigma_s Parameter = {}".format(scatt_2)) + param = scatt_2 + except: + scatt_2=1.0 + print("Scatterer Sigma_s Nominal = {}".format(scatt_2)) + + + try: + print("Source Parameter = {}".format(param_q)) + param = param_q + except: + param_q=1.0 + print("Source Nominal = {}".format(param_q)) + + try: + print("Parameter id = {}".format(p_id)) + except: + p_id=0 + print("Parameter id = {}".format(p_id)) + + try: + if phase == 0: + print("Offline Phase") + phase = "offline" + elif phase == 1: + print("Merge Phase") + phase = "merge" + elif phase == 2: + print("Systems Phase") + phase = "systems" + elif phase == 3: + print("Online Phase") + phase = "online" + except: + phase="offline" + print("Phase default to offline") + + # Check number of processors + num_procs = 4 + if size != num_procs: + sys.exit(f"Incorrect number of processors. Expected {num_procs} processors but got {size}.") + + # Setup mesh + nodes = [] + N = 70 + L = 7 + xmin = 0 + dx = L / N + for i in range(N + 1): + nodes.append(xmin + i * dx) + meshgen = OrthogonalMeshGenerator(node_sets=[nodes, nodes], + partitioner=KBAGraphPartitioner( + nx=2, + ny=2, + xcuts=[3.5], + ycuts=[3.5], + )) + grid = meshgen.Execute() + + # Set background (Scatterer) block ID = 0 + grid.SetUniformBlockID(0) + + # Set Source (central red square from x=3 to x=4, y=3 to y=4) block ID = 1 + logvol_src = RPPLogicalVolume(xmin=3.0, xmax=4.0, + ymin=3.0, ymax=4.0, + infz=True) + grid.SetBlockIDFromLogicalVolume(logvol_src, 1, True) + + # Set Absorbers (green 1x1 squares) block ID = 2 + absorber_centers = [ + (1,1), (3,1), (5,1), + (2,2), (4,2), + (1,3), (5,3), + (2,4), (4,4), + (1,5), (5,5) + ] + for xc, yc in absorber_centers: + vol_abs = RPPLogicalVolume( + xmin=xc+0.0, xmax=xc+1.0, + ymin=yc+0.0, ymax=yc+1.0, + infz=True + ) + grid.SetBlockIDFromLogicalVolume(vol_abs, 2, True) + + scatt_t = abs_2 + scatt_2 + num_groups = 1 + scatterer = MultiGroupXS() + scatterer.CreateSimpleOneGroup(sigma_t=scatt_t, c=scatt_2/scatt_t) + + abs_t = abs_1 + scatt_1 + + absorber = MultiGroupXS() + absorber.CreateSimpleOneGroup(sigma_t=abs_t, c=scatt_1/abs_t) + + strength = [0.0] + src0 = VolumetricSource(block_ids=[0], group_strength=strength) + strength = [param_q] + src1 = VolumetricSource(block_ids=[1], group_strength=strength) + + # Setup Physics + pquad = GLCProductQuadrature2DXY(n_polar=4, n_azimuthal=32, scattering_order=0) + #pquad = GLCProductQuadrature2DXY(n_polar=8, n_azimuthal=256, scattering_order=0) + + if phase == "online": + rom_options={ + "param_id":0, + "phase":phase, + "param_file":"data/interpolation_params.txt", + "new_point":[scatt_1, scatt_2, abs_1, abs_2, param_q] + } + else: + rom_options={ + "param_id":p_id, + "phase":phase + } + + phys = DiscreteOrdinatesProblem( + mesh=grid, + num_groups=num_groups, + groupsets=[ + { + "groups_from_to": [0, 0], + "angular_quadrature": pquad, + "angle_aggregation_num_subsets": 1, + "inner_linear_method": "petsc_gmres", + "l_abs_tol": 1.0e-8, + "l_max_its": 300, + "gmres_restart_interval": 100, + }, + ], + xs_map=[ + {"block_ids": [0, 1], "xs": scatterer}, + {"block_ids": [2], "xs": absorber} + ], + volumetric_sources= [src0, src1] + ) + + rom = ROMProblem(problem=phys,options=rom_options) + + ss_solver = SteadyStateROMSolver(problem=phys, rom_problem=rom) + ss_solver.Initialize() + ss_solver.Execute() + + phys.ComputeBalance() + + if phase == "online": + phys.WriteFluxMoments("output/rom") + if phase == "offline": + phys.WriteFluxMoments("output/fom") \ No newline at end of file diff --git a/examples/checkerboard/run_rom_checkerboard.py b/examples/checkerboard/run_rom_checkerboard.py new file mode 100644 index 0000000..9b26780 --- /dev/null +++ b/examples/checkerboard/run_rom_checkerboard.py @@ -0,0 +1,92 @@ +import numpy as np +import sys, os +sys.path.insert(0, os.path.realpath("../python")) + +import plotting +import utils + +# Sampling training points +bounds = [[0,5.0],[0.5,1.5],[7.5,12.5],[0.0,0.5],[0.1,1]] +num_params = 50 + +params = utils.sample_parameter_space(bounds, num_params) +#params = np.loadtxt("data/params.txt") + +#params = params[:num_params,:] +np.savetxt("data/interpolation_params.txt", params) +# OFFLINE PHASE +phase = 0 + +for i in range(num_params): + cmd = "mpiexec -n=4 ../../build/rom_app_exec -i base_checkerboard.py \ + -p phase={} -p scatt_1={} -p scatt_2={} -p abs_1={} -p abs_2={} -p param_q={} -p p_id={}"\ + .format(phase,params[i][0],params[i][1],params[i][2],params[i][3],params[i][4],i) + utils.run_opensn(cmd) + +# MERGE PHASE +phase = 1 + +cmd = "mpiexec -n=4 ../../build/rom_app_exec -i base_checkerboard.py -p phase={} -p p_id={}".format(phase, num_params-1) +utils.run_opensn(cmd) + +plotting.plot_sv(num_groups=1) + + +# SYSTEMS PHASE +phase = 2 +for i in range(num_params): + cmd = "mpiexec -n=4 ../../build/rom_app_exec -i base_checkerboard.py \ + -p phase={} -p scatt_1={} -p scatt_2={} -p abs_1={} -p abs_2={} -p param_q={} -p p_id={}"\ + .format(phase,params[i][0],params[i][1],params[i][2],params[i][3],params[i][4],i) + utils.run_opensn(cmd) + +np.savetxt("data/params.txt", params) + + +# Generate Test Data +test_scatt_1 = np.random.uniform(0,5.0,10) +test_scatt_2 = np.random.uniform(0.5,1.5,10) +test_abs_1 = np.random.uniform(7.5,12.5,10) +test_abs_2 = np.random.uniform(0.0,0.5,10) +test_q = np.random.uniform(0.1,1,10) +test = np.append(test_scatt_1[:,np.newaxis], test_scatt_2[:,np.newaxis], axis=1) +test = np.append(test, test_abs_1[:,np.newaxis], axis=1) +test = np.append(test, test_abs_2[:,np.newaxis], axis=1) +test = np.append(test, test_q[:,np.newaxis], axis=1) +np.savetxt("data/validation.txt", test) +# test = np.loadtxt("data/validation.txt") + +num_test = 10 +errors = [] +speedups = [] + +for i in range(num_test): + # ONLINE PHASE + phase = 3 + cmd = "mpiexec -n=4 ../../build/rom_app_exec -i base_checkerboard.py \ + -p phase={} -p scatt_1={} -p scatt_2={} -p abs_1={} -p abs_2={} -p param_q={} -p p_id={}"\ + .format(phase,test[i][0],test[i][1],test[i][2],test[i][3],test[i][4],i) + utils.run_opensn(cmd) + rom_time = np.loadtxt("results/online_time.txt") + + # Reference FOM solution + phase = 0 + cmd = "mpiexec -n=4 ../../build/rom_app_exec -i base_checkerboard.py \ + -p phase={} -p scatt_1={} -p scatt_2={} -p abs_1={} -p abs_2={} -p param_q={} -p p_id={}"\ + .format(phase,test[i][0],test[i][1],test[i][2],test[i][3],test[i][4],i) + utils.run_opensn(cmd) + fom_time = np.loadtxt("results/offline_time.txt") + + plotting.plot_2d_flux("output/fom{}.h5", ranks=range(4), prefix="fom", pid=i) + plotting.plot_2d_flux("output/rom{}.h5", ranks=range(4), prefix="rom", pid=i) + + error = plotting.plot_2d_lineout(ranks=range(4), pid=i) + + errors.append(error) + speedups.append(fom_time/rom_time) + + +print("Avg Error ", np.mean(errors)) +np.savetxt("results/errors.txt", errors) +print("Avg Speedup ", np.mean(speedups)) +np.savetxt("results/speedups.txt", speedups) \ No newline at end of file diff --git a/examples/python/plot_errors.py b/examples/python/plot_errors.py new file mode 100644 index 0000000..471b1e3 --- /dev/null +++ b/examples/python/plot_errors.py @@ -0,0 +1,84 @@ +import numpy as np +import matplotlib.pyplot as plt + +problem = "../checkerboard" + +# List of sample sizes +sample_sizes = [100, 150, 200, 250, 300, 350, 400, 450, 500] + +# Collect error data +all_errors = [] +for size in sample_sizes: + filename = f"{problem}/results/errors_{size}.txt" + data = np.loadtxt(filename) + all_errors.append(data) + +# Make the boxplot +plt.figure(figsize=(8, 5)) +plt.boxplot(all_errors, positions=sample_sizes, widths=20, showfliers=False) + +plt.xlabel("Training Set Size") +plt.ylabel("$L_2$ Error") +plt.yscale("log") +plt.grid(True, linestyle="--", alpha=0.6) + +plt.savefig(f'{problem}/results/errors_set_size.jpg') +plt.close() + +# Collect speedup data +all_speedups = [] +for size in sample_sizes: + filename = f"{problem}/results/speedups_{size}.txt" + data = np.loadtxt(filename) + all_speedups.append(data) + +# Make the boxplot +plt.figure(figsize=(8, 5)) +plt.boxplot(all_speedups, positions=sample_sizes, widths=20, showfliers=False) + +plt.xlabel("Training Set Size") +plt.ylabel("Speedup") +plt.yscale("log") +plt.grid(True, linestyle="--", alpha=0.6) + +plt.savefig(f'{problem}/results/speedups_set_size.jpg') + +# List of ranks +ranks = [5, 10, 15, 20] + +# Collect error data +all_errors = [] +for rank in ranks: + filename = f"{problem}/results/errors_{rank}.txt" + data = np.loadtxt(filename) + all_errors.append(data) + +# Make the boxplot +plt.figure(figsize=(8, 5)) +plt.boxplot(all_errors, positions=ranks, widths=2, showfliers=False) + +plt.xlabel("Rank") +plt.ylabel("$L_2$ Error") +plt.yscale("log") +plt.grid(True, linestyle="--", alpha=0.6) + +plt.savefig(f'{problem}/results/errors_rank.jpg') +plt.close() + +# Collect speedup data +all_speedups = [] +for rank in ranks: + filename = f"{problem}/results/speedups_{rank}.txt" + data = np.loadtxt(filename) + all_speedups.append(data) + +# Make the boxplot +plt.figure(figsize=(8, 5)) +plt.boxplot(all_speedups, positions=ranks, widths=2, showfliers=False) + +plt.xlabel("Rank") +plt.ylabel("Speedup") +plt.yscale("log") +plt.grid(True, linestyle="--", alpha=0.6) + +plt.savefig(f'{problem}/results/speedups_rank.jpg') \ No newline at end of file diff --git a/examples/python/plotting.py b/examples/python/plotting.py new file mode 100644 index 0000000..5f703ab --- /dev/null +++ b/examples/python/plotting.py @@ -0,0 +1,116 @@ +import numpy as np +import matplotlib.pyplot as plt +from utils import * +import scipy.interpolate +from matplotlib.colors import LogNorm + + +def plot_2d_flux(file_pattern, ranks, moment=0, prefix="fom", grid_res=200, pid=0): + """Create smooth full-color plots (not scatter) for each energy group.""" + xs, ys, vals, G = load_2d_flux(file_pattern, ranks, moment=moment) + + for g in range(G): + # Create regular grid + 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") + + vmin = max(np.nanmin(Z), 1e-10) + vmax = np.nanmax(Z) + norm = LogNorm(vmin=vmin, vmax=vmax) + + plt.figure(figsize=(6, 5)) + im = plt.imshow( + Z, + extent=[xi.min(), xi.max(), yi.min(), yi.max()], + origin="lower", + aspect="equal", + cmap="viridis", + norm=norm + ) + plt.title(f"{prefix.upper()} Group {g} (moment {moment})") + plt.xlabel("x") + plt.ylabel("y") + cbar = plt.colorbar(im) + cbar.set_label("Flux") + outpath = f"results/{prefix}_group_{g}_{pid}.png" + plt.tight_layout() + plt.savefig(outpath, dpi=200) + plt.close() + +def plot_2d_lineout(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("output/rom{}.h5", ranks, moment=moment) + xs_, ys_, vals_, G = load_2d_flux("output/fom{}.h5", ranks, moment=moment) + + for g in range(G): + # Create regular grid + 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.asarray(vals_)-np.asarray(vals))/np.linalg.norm(np.asarray(vals_)) + return error + + +def plot_sv(num_groups): + for i in range(num_groups): + S = np.loadtxt("data/singular_values_g{}.txt".format(i)) + plt.semilogy(S, 'o-') + plt.xlabel("Mode index") + plt.ylabel("Singular value") + plt.title("Singular value decay") + plt.grid(True) + plt.tight_layout() + plt.savefig("results/svd_decay_{}.jpg".format(i)) + plt.close() + + +def plot_1d_flux(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): + 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/examples/python/utils.py b/examples/python/utils.py new file mode 100644 index 0000000..8a16c49 --- /dev/null +++ b/examples/python/utils.py @@ -0,0 +1,125 @@ +import numpy as np +import subprocess +import h5py + +def run_opensn(cmd): + args = cmd.split(" ") + print(args) + process = subprocess.Popen(args, stdout=subprocess.PIPE, stderr=subprocess.PIPE, text=True) + output, errors = process.communicate() + print("Output:", output) + print("Errors:", errors) + + +def load_2d_flux(file_pattern, ranks, moment=0): + """Load (x, y, flux) grouped by energy group from HDF5 files.""" + with h5py.File(file_pattern.format(ranks[0]), "r") as f0: + num_groups = int(f0.attrs["num_groups"]) + num_moments = int(f0.attrs["num_moments"]) + + xs = [[] for _ in range(num_groups)] + ys = [[] for _ in range(num_groups)] + vals = [[] for _ in range(num_groups)] + + for r in ranks: + fp = file_pattern.format(r) + with h5py.File(fp, "r") as f: + x = f["mesh/nodes_x"][:] + y = f["mesh/nodes_y"][:] + values = f["values"][:] + num_nodes = x.size + + values = values.reshape(num_nodes, num_moments, num_groups) + for g in range(num_groups): + xs[g].append(x) + ys[g].append(y) + vals[g].append(values[:, moment, g]) + + for g in range(num_groups): + xs[g] = np.concatenate(xs[g]) + ys[g] = np.concatenate(ys[g]) + vals[g] = np.concatenate(vals[g]) + + return xs, ys, vals, num_groups + +def load_1d_flux(file_pattern, ranks, moment=0): + """Load concatenated 1-D (x, flux) data per energy group.""" + with h5py.File(file_pattern.format(ranks[0]), "r") as f0: + num_groups = int(f0.attrs["num_groups"]) + num_moments = int(f0.attrs["num_moments"]) + + xs = [[] for _ in range(num_groups)] + vals = [[] for _ in range(num_groups)] + + for r in ranks: + fp = file_pattern.format(r) + with h5py.File(fp, "r") as f: + x = f["mesh/nodes_z"][:] + values = f["values"][:] + num_nodes = x.size + + values = values.reshape(num_nodes, num_moments, num_groups) + + for g in range(num_groups): + xs[g].append(x) + vals[g].append(values[:, moment, g]) + + for g in range(num_groups): + xs[g] = np.concatenate(xs[g]) + vals[g] = np.concatenate(vals[g]) + + return xs, vals, num_groups + + +def sample_parameter_space(bounds, n_samples): + n_dim = len(bounds) + n_vertices = 2**n_dim + n_random = n_samples - n_vertices + + # Random interior samples + random_samples = np.array([ + [np.random.uniform(low, high) for (low, high) in bounds] + for _ in range(n_random) + ]) + + # Vertices of domain + vertices = np.zeros((n_vertices, n_dim)) + for i in range(n_vertices): + for d, (low, high) in enumerate(bounds): + if (i >> d) & 1: + vertices[i, d] = high + else: + vertices[i, d] = low + + samples = np.vstack([random_samples, vertices]) + return samples + + +def update_xs(in_file, out_file, sigma_t_vec, S): + with open(in_file, "r") as f: + lines = f.readlines() + + # --- 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" + + # --- 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) + + 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_GPRIME_G_VAL 0 {gprime} {g} {val:.12g}\n") + + lines[tb+1:te] = new_tm + + with open(out_file, "w") as f: + f.writelines(lines) \ No newline at end of file diff --git a/examples/reed/base_reed.py b/examples/reed/base_reed.py new file mode 100644 index 0000000..57622ce --- /dev/null +++ b/examples/reed/base_reed.py @@ -0,0 +1,158 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +# Standard Reed 1D 1-group problem + +import os +import sys + +if "opensn_console" not in globals(): + from mpi4py import MPI + size = MPI.COMM_WORLD.size + rank = MPI.COMM_WORLD.rank + sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), "../../../../../"))) + from pyopensn.mesh import OrthogonalMeshGenerator + from pyopensn.xs import MultiGroupXS + from pyopensn.source import VolumetricSource + from pyopensn.aquad import GLProductQuadrature1DSlab + from pyopensn.solver import DiscreteOrdinatesProblem, SteadyStateSolver + from pyopensn.logvol import RPPLogicalVolume + from pyopensn.rom import ROMProblem, SteadyStateROMSolver + +if __name__ == "__main__": + + try: + print("Scattering Parameter = {}".format(scatt)) + param = scatt + except: + scatt=0.0 + print("Scattering Nominal = {}".format(scatt)) + + try: + print("Cross Section Parameter = {}".format(sigma_t)) + param = sigma_t + except: + sigma_t=1.0 + print("Cross Section Nominal = {}".format(sigma_t)) + + try: + print("Source Parameter = {}".format(param_q)) + param = param_q + except: + param_q=1.0 + print("Source Nominal = {}".format(param_q)) + + try: + print("Parameter id = {}".format(p_id)) + except: + p_id=0 + print("Parameter id = {}".format(p_id)) + + try: + if phase == 0: + print("Offline Phase") + phase = "offline" + elif phase == 1: + print("Merge Phase") + phase = "merge" + elif phase == 2: + print("Systems Phase") + phase = "systems" + elif phase == 3: + print("Online Phase") + phase = "online" + except: + phase="offline" + print("Phase default to offline") + + # Create Mesh + widths = [2., 1., 2., 1., 2.] + nrefs = [200, 200, 200, 200, 200] + 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() + + # Set block IDs + z_min = 0.0 + z_max = widths[1] + for imat in range(Nmat): + z_max = z_min + widths[imat] + print("imat=", imat, ", zmin=", z_min, ", zmax=", z_max) + lv = RPPLogicalVolume(infx=True, infy=True, zmin=z_min, zmax=z_max) + grid.SetBlockIDFromLogicalVolume(lv, imat, True) + z_min = z_max + + # Add cross sections to materials + total = [50., 5., 0., sigma_t, sigma_t] + c = [0., 0., 0., scatt, scatt] + xs_map = len(total) * [None] + for imat in range(Nmat): + xs_ = MultiGroupXS() + xs_.CreateSimpleOneGroup(total[imat], c[imat]) + xs_map[imat] = { + "block_ids": [imat], "xs": xs_, + } + + # Create sources in 1st and 4th materials + src0 = VolumetricSource(block_ids=[0], group_strength=[50.]) + src1 = VolumetricSource(block_ids=[3], group_strength=[param_q]) + + # Angular Quadrature + gl_quad = GLProductQuadrature1DSlab(n_polar=128, scattering_order=0) + + # LBS block option + num_groups = 1 + + phys = DiscreteOrdinatesProblem( + mesh=grid, + num_groups=num_groups, + groupsets=[ + { + "groups_from_to": (0, num_groups - 1), + "angular_quadrature": gl_quad, + "inner_linear_method": "petsc_gmres", + "l_abs_tol": 1.0e-9, + "l_max_its": 300, + "gmres_restart_interval": 30, + }, + ], + xs_map=xs_map, + volumetric_sources= [src0, src1], + boundary_conditions= [ + {"name": "zmin", "type": "vacuum"}, + {"name": "zmax", "type": "vacuum"} + ] + ) + + if phase == "online": + rom_options = { + "param_id": 0, + "phase": phase, + "param_file": "data/params.txt", + "new_point": [scatt, param_q] + } + else: + rom_options = { + "param_id": p_id, + "phase": phase + } + + rom = ROMProblem(problem=phys,options=rom_options) + + # Initialize and execute solver + ss_solver = SteadyStateROMSolver(problem=phys, rom_problem=rom) + ss_solver.Initialize() + ss_solver.Execute() + + # compute particle balance + phys.ComputeBalance() + + if phase == "online": + phys.WriteFluxMoments("output/rom") + if phase == "offline": + phys.WriteFluxMoments("output/fom") \ No newline at end of file diff --git a/examples/reed/run_rom_reed.py b/examples/reed/run_rom_reed.py new file mode 100644 index 0000000..a1274c2 --- /dev/null +++ b/examples/reed/run_rom_reed.py @@ -0,0 +1,70 @@ +import numpy as np +import sys, os +sys.path.insert(0, os.path.realpath("../python")) + +import plotting +import utils + +# Sampling training points +bounds = [[0.0,1.0],[0.0,1.0]] +num_params = 100 + +params = utils.sample_parameter_space(bounds, num_params) + +# OFFLINE PHASE +phase = 0 + +for i, param in enumerate(params): + cmd = "mpiexec -n 2 ../../build/rom_app_exec -i base_reed.py -p phase={} -p scatt={} -p param_q={} -p p_id={}"\ + .format(phase,param[0], param[1], i) + utils.run_opensn(cmd) + +# MERGE PHASE +phase = 1 + +cmd = "mpiexec -n 2 ../../build/rom_app_exec -i base_reed.py -p phase={} -p p_id={}".format(phase, i) +utils.run_opensn(cmd) + +plotting.plot_sv(num_groups=1) + + +# SYSTEMS PHASE +phase = 2 + +for i, param in enumerate(params): + cmd = "mpiexec -n 2 ../../build/rom_app_exec -i base_reed.py -p phase={} -p scatt={} -p param_q={} -p p_id={}"\ + .format(phase,param[0], param[1], i) + utils.run_opensn(cmd) + +np.savetxt("data/params.txt", params) + +# Generate Test Data +test = np.random.uniform(0,1,[10,2]) + +errors = [] +speedups = [] + +for i, param in enumerate(test): + # ONLINE PHASE + phase = 3 + cmd = "mpiexec -n 2 ../../build/rom_app_exec -i base_reed.py -p phase={} -p scatt={} -p param_q={} -p p_id={}"\ + .format(phase,param[0],param[1], i) + utils.run_opensn(cmd) + rom_time = np.loadtxt("results/online_time.txt") + + # Reference FOM solution + phase = 0 + cmd = "mpiexec -n 2 ../../build/rom_app_exec -i base_reed.py -p phase={} -p scatt={} -p param_q={} -p p_id={}"\ + .format(phase,param[0],param[1], i) + utils.run_opensn(cmd) + fom_time = np.loadtxt("results/offline_time.txt") + + error = plotting.plot_1d_flux("output/fom{}.h5", "output/rom{}.h5", ranks=range(2), pid=i) + + errors.append(error) + speedups.append(fom_time/rom_time) + +print("Avg Error ", np.mean(errors)) +np.savetxt("results/errors.txt", errors) +print("Avg Speedup ", np.mean(speedups)) +np.savetxt("results/speedups.txt", speedups) \ No newline at end of file diff --git a/src/main.cc b/src/main.cc index 7c59a53..397544b 100644 --- a/src/main.cc +++ b/src/main.cc @@ -1,11 +1,29 @@ -#include "test.h" -#include "opensn/python/lib/py_app.h" -#include "opensn/mpicpp-lite/mpicpp-lite.h" +// SPDX-FileCopyrightText: 2025 The OpenSn Authors +// SPDX-License-Identifier: MIT + +#include "rom_py_app.h" +#include "petsc.h" int main(int argc, char** argv) { - mpi::Environment env(argc, argv); + mpi::Environment env(argc, argv); + + PetscCall(PetscInitializeNoArguments()); // NOLINT(bugprone-casting-through-void) + + int retval = EXIT_SUCCESS; + try + { py::scoped_interpreter guard{}; - auto app = TestApp::Create({}); - return 0; + rompy::ROMApp app(MPI_COMM_WORLD); // NOLINT(bugprone-casting-through-void) + retval = app.Run(argc, argv); + } + catch (...) + { + std::fprintf(stderr, "Unknown fatal error\n"); + retval = EXIT_FAILURE; + } + + PetscFinalize(); + + return retval; } diff --git a/src/py_wrappers.h b/src/py_wrappers.h new file mode 100644 index 0000000..24a1822 --- /dev/null +++ b/src/py_wrappers.h @@ -0,0 +1,21 @@ +// SPDX-FileCopyrightText: 2025 The OpenSn Authors +// SPDX-License-Identifier: MIT + +#pragma once + +#include "framework/parameters/parameter_block.h" +#include +#include + +namespace py = pybind11; +namespace opensn +{ + +// Wrap the ROM components of app +void py_rom(py::module& pyopensn); +void WrapROM(py::module& ROM); + +/// Translate a Python dictionary into a ParameterBlock. +ParameterBlock kwargs_to_param_block(const py::kwargs& params); + +} // namespace opensn \ No newline at end of file diff --git a/src/rom.cc b/src/rom.cc new file mode 100644 index 0000000..4d90902 --- /dev/null +++ b/src/rom.cc @@ -0,0 +1,132 @@ +// SPDX-FileCopyrightText: 2025 The OpenSn Authors +// SPDX-License-Identifier: MIT + +#include "py_wrappers.h" +#include "rom/rom_problem.h" +#include "rom/steady_state_rom_solver.h" +#include "modules/solver.h" +#include +#include + +namespace opensn +{ +namespace py = pybind11; + +template +static InputParameters KwargsToParams(const py::kwargs& kw) +{ + auto params = T::GetInputParameters(); + params.AssignParameters(kwargs_to_param_block(kw)); + return params; +} + + +// clang-format off +void WrapROM(py::module& m) +{ + // ROMProblem + auto rom_problem = py::class_, + Problem>( + m, "ROMProblem", + R"( + ROM controller for reduced-order modeling workflows. + + Wrapper of :cpp:class:`opensn::ROMProblem`. + )"); + + rom_problem.def( + py::init([](py::kwargs kw) + { + auto params = ROMProblem::GetInputParameters(); + params.AssignParameters(kwargs_to_param_block(kw)); + return std::make_shared(params); + }), + R"( + ROMProblem(**kwargs) + + Construct a ROMProblem directly from keyword arguments. + )"); + + rom_problem.def_static( + "Create", + [](py::kwargs kw) + { + return ROMProblem::Create(kwargs_to_param_block(kw)); + }, + R"( + Create(**kwargs) -> ROMProblem + + Factory constructor (recommended). Accepts the same kwargs as the direct constructor: + - problem : LBSProblem (required) + - options : dict (optional) + - name : str (optional) + )"); + + + // SteadyStateROMSolver + auto ss_rom_solver = + py::class_, + Solver>( + m, + "SteadyStateROMSolver", + R"( + Steady-state ROM driver. + + Wrapper of :cpp:class:`opensn::SteadyStateROMSolver`. + + 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. + )" + ); + + ss_rom_solver.def( + py::init([](py::kwargs kw) + { + auto params = KwargsToParams(kw); + + return std::make_shared(params); + }), + R"( + SteadyStateROMSolver(**kwargs) + + Construct a steady-state driver that dispatches to ROM or FOM paths + depending on the ROM options and phase. + )" + ); + + ss_rom_solver + .def("Initialize", &SteadyStateROMSolver::Initialize, + R"( + Initialize() + + Prepare the solver and ROM controller for execution. + )") + .def("Execute", &SteadyStateROMSolver::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 + )"); +} +// clang-format on + + +void py_rom(py::module& pyopensn) +{ + auto rom = pyopensn.def_submodule("rom", "ROM module."); + WrapROM(rom); +} + +} // namespace opensn \ No newline at end of file diff --git a/src/rom/rom_problem.cc b/src/rom/rom_problem.cc new file mode 100644 index 0000000..9849890 --- /dev/null +++ b/src/rom/rom_problem.cc @@ -0,0 +1,439 @@ +// SPDX-FileCopyrightText: 2024 The OpenSn Authors +// SPDX-License-Identifier: MIT + +#include "modules/linear_boltzmann_solvers/lbs_problem/lbs_problem.h" +#include "modules/linear_boltzmann_solvers/lbs_problem/iterative_methods/wgs_context.h" +#include "framework/object_factory.h" +#include "framework/logging/log.h" +#include "framework/runtime.h" +#include "rom_problem.h" +#include "linalg/BasisGenerator.h" +#include "linalg/BasisReader.h" +#include +#include + +namespace opensn +{ + +OpenSnRegisterObjectInNamespace(rom, ROMProblem); + +InputParameters +ROMProblem::GetInputParameters() +{ + InputParameters params = Problem::GetInputParameters(); + + params.SetClassName("ROMProblem"); + + params.ChangeExistingParamToOptional("name", "ROMProblem"); + + params.AddRequiredParameter>( + "problem", "An existing LBS problem to attach the ROM controller to."); + + // Optional nested ROM options block (phase, ids, files, new_point, etc.) + params.AddOptionalParameterBlock( + "options", ParameterBlock(), "ROM control options (phase, param_id, param_file, new_point, …)"); + + return params; +} + +std::shared_ptr +ROMProblem::Create(const ParameterBlock& params) +{ + auto& factory = opensn::ObjectFactory::GetInstance(); + return factory.Create("rom::ROMProblem", params); +} + +ROMProblem::ROMProblem(const InputParameters& params) + : Problem(params), + lbs_problem(params.GetSharedPtrParam("problem")) +{ + // Initialize options + if (params.IsParameterValid("options")) + { + auto options_params = ROMProblem::GetOptionsBlock(); + options_params.AssignParameters(params.GetParam("options")); + SetOptions(options_params); + } +} + +void +ROMProblem::TakeSample(int id) +{ + bool update_right_SV = false; + int max_num_snapshots = 100; + bool isIncremental = false; + const std::string basisName = "basis/snapshots_"; + + auto num_moments = lbs_problem->GetNumMoments(); + auto num_groups = lbs_problem->GetNumGroups(); + auto num_local_nodes = lbs_problem->GetLocalNodeCount(); + std::vector phi_new_local = lbs_problem->GetPhiNewLocal(); + + for (int g = 0; g < num_groups; ++g) + { + const std::string basisFileName = basisName + std::to_string(g) + "_" + std::to_string(id); + + int group_dim = num_local_nodes * num_moments; + + CAROM::Options options(group_dim, max_num_snapshots, update_right_SV); + CAROM::BasisGenerator generator(options, isIncremental, basisFileName); + + std::vector phi_group(group_dim, 0.0); + + for (int n = 0; n < num_local_nodes; ++n) + { + size_t node_base_full = n * num_moments * num_groups; + size_t node_base_group = n * num_moments; + + for (int m = 0; m < num_moments; ++m) + { + auto idx_full = node_base_full + m * num_groups + g; + auto idx_group = node_base_group + m; + + phi_group[idx_group] = phi_new_local[idx_full]; + } + } + + generator.takeSample(phi_group.data()); + generator.writeSnapshot(); + } +} + +void +ROMProblem::MergePhase(int nsnaps) +{ + bool update_right_SV = false; + int max_num_snapshots = 300; + bool isIncremental = false; + + auto num_moments = lbs_problem->GetNumMoments(); + auto num_groups = lbs_problem->GetNumGroups(); + auto num_local_nodes = lbs_problem->GetLocalNodeCount(); + auto group_dim = num_local_nodes * num_moments; + auto full_dim = num_local_nodes * num_moments * num_groups; + + CAROM::Options options(group_dim, max_num_snapshots, update_right_SV); + double tol = 1e-20; + romRank = 20; + options.setSingularValueTol(tol); + options.setMaxBasisDimension(romRank); + + for (auto g = 0; g < num_groups; ++g) + { + auto basis_prefix = "basis/basis_" + std::to_string(g); + CAROM::BasisGenerator loader(options, isIncremental, basis_prefix); + + for (auto paramID = 0; paramID < nsnaps; ++paramID) + { + auto snap_file = "basis/snapshots_" + std::to_string(g) + "_" + std::to_string(paramID) + "_snapshot"; + loader.loadSamples(snap_file, "snapshot"); + } + + loader.endSamples(); + + // Save singular values per group + if (opensn::mpi_comm.rank() == 0) + { + auto sv_file = "data/singular_values_g" + std::to_string(g) + ".txt"; + std::ofstream sv_out(sv_file); + auto S_vec = loader.getSingularValues(); + for (auto i = 0; i < S_vec->dim(); ++i) + sv_out << std::setprecision(16) << S_vec->item(i) << "\n"; + sv_out.close(); + log.Log() << "Saved singular values for group " << g << " to " << sv_file << "\n"; + } + } +} + +void +ROMProblem::ReadParamMatrix(const std::string& filename) +{ + param_points_.clear(); + + std::ifstream infile(filename); + std::string line; + + while (std::getline(infile, line)) + { + std::istringstream iss(line); + std::vector row; + double val; + while (iss >> val) + row.push_back(val); + + if (!row.empty()) + param_points_.emplace_back(row.data(), static_cast(row.size()),false,true); + } +} + +std::shared_ptr +ROMProblem::AssembleAU() +{ + 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; + + std::vector> Ugs; + 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) romRank = Ug->numColumns(); + Ugs.push_back(std::move(Ug)); + } + + auto AU = std::make_shared(num_local_dofs, romRank * num_groups, true); + + // Assuming one groupset for ROM problems + auto wgs_solvers = lbs_problem->GetWGSSolvers(); + auto raw_context = wgs_solvers.front()->GetContext(); + auto gs_context = std::dynamic_pointer_cast(raw_context); + const auto scope = gs_context->lhs_src_scope; + + 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 < romRank; ++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 * romRank + r; + for (size_t i = 0; i < static_cast(num_local_dofs); ++i) + AU->item(i, static_cast(col_idx)) = basis_local[i] - phi_new_local[i]; + } + } + return AU; +} + +std::shared_ptr +ROMProblem::AssembleRHS() +{ + 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 b = std::make_shared(num_local_dofs, true); + + // Assuming one groupset for ROM problems + auto wgs_solvers = lbs_problem->GetWGSSolvers(); + auto raw_context = wgs_solvers.front()->GetContext(); + auto gs_context_ptr = std::dynamic_pointer_cast(raw_context); + auto scope = gs_context_ptr->rhs_src_scope; + + auto& phi_old_local = lbs_problem->GetPhiOldLocal(); + auto& q_moments_local = lbs_problem->GetQMomentsLocal(); + + q_moments_local.assign(q_moments_local.size(), 0.0); + gs_context_ptr->set_source_function(gs_context_ptr->groupset, q_moments_local, phi_old_local, scope); + + // Sweep + gs_context_ptr->ApplyInverseTransportOperator(scope); + std::vector phi_new_local = lbs_problem->GetPhiNewLocal(); + + for (int i = 0; i < num_local_dofs; ++i) + (*b)(i) = phi_new_local[i]; + return b; +} + +void +ROMProblem::AssembleROM( + std::shared_ptr& AU, + std::shared_ptr& b, + const std::string& Ar_filename, + const std::string& rhs_filename) +{ + // rhs = AU^T * b + auto rhs = AU->transposeMult(*b); + + // Ar = AU^T * AU + auto Ar = AU->transposeMult(*AU); + + // Save + Ar->write(Ar_filename); + rhs->write(rhs_filename); +} + +void +ROMProblem::SolveROM( + std::shared_ptr& Ar, + std::shared_ptr& rhs) +{ + auto Ar_inv = std::make_shared(Ar->numRows(), Ar->numColumns(), false); + + Ar->inverse(*Ar_inv); + + auto c_vec = Ar_inv->mult(*rhs); + + 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; + + 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) romRank = 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); + + for (int g = 0; g < num_groups; ++g) + { + for (int r = 0; r < romRank; ++r) + { + const int cr_idx = g * romRank + r; + const double cr = (*c_vec)(cr_idx); + + 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) + { + const size_t row_phi = n * (num_moments * num_groups) + m * num_groups + static_cast(g); + phi_new_local[row_phi] += cr * col_g->item(row_g); + } + } + } +} + +void +ROMProblem::SetupInterpolator(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 + 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 + std::vector> rotations; + int rom_dim = Ar_matrices[0]->numRows(); + + for (size_t i = 0; i < Ar_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); + + Ar_interp_obj_ptr = std::make_unique( + param_points_, rotations, Ar_matrices, + ref_index, "SPD", "G", "LS", 0.999, false); + rhs_interp_obj_ptr = std::make_unique( + param_points_, rotations, rhs_vectors, + ref_index, "G", "LS", 0.999, false); +} + +void +ROMProblem::InterpolateArAndRHS( + CAROM::Vector& desired_point, + std::shared_ptr& Ar_interp, + std::shared_ptr& rhs_interp) +{ + Ar_interp = Ar_interp_obj_ptr->interpolate(desired_point); + rhs_interp = rhs_interp_obj_ptr->interpolate(desired_point); +} + +InputParameters +ROMProblem::GetOptionsBlock() +{ + InputParameters params; + + params.SetGeneralDescription("Set options from a list of parameters"); + params.AddOptionalParameter("param_id", 0, "A parameter id for parametric problems."); + params.AddOptionalParameter("phase", "offline", "The phase (offline, online, or merge) for ROM purposes."); + params.AddOptionalParameter("param_file", "", "A file containing an array of parameters for ROM."); + params.AddOptionalParameterArray("new_point", {0.0}, "New parameter point for ROM."); + + return params; +} + +void +ROMProblem::SetOptions(const InputParameters& input) +{ + auto params = ROMProblem::GetOptionsBlock(); + params.AssignParameters(input); + + for (size_t p = 0; p < params.GetNumParameters(); ++p) + { + const auto& spec = params.GetParam(p); + + if (spec.GetName() == "param_id") + options_.param_id = spec.GetValue(); + + else if (spec.GetName() == "phase") + options_.phase = spec.GetValue(); + + else if (spec.GetName() == "param_file") + options_.param_file = spec.GetValue(); + + else if (spec.GetName() == "new_point") + { + spec.RequireBlockTypeIs(ParameterBlockType::ARRAY); + + std::vector vals; + vals.reserve(spec.GetNumParameters()); + for (const auto& sub_param : spec) + vals.push_back(sub_param.GetValue()); + + options_.new_point = std::make_unique(static_cast(vals.size()), false); + for (int i = 0; i < static_cast(vals.size()); ++i) + (*(options_.new_point))(i) = vals[static_cast(i)]; + } + } // for p +} + +ROMOptions& +ROMProblem::GetOptions() +{ + return options_; +} + +} // namespace opensn diff --git a/src/rom/rom_problem.h b/src/rom/rom_problem.h new file mode 100644 index 0000000..ad80171 --- /dev/null +++ b/src/rom/rom_problem.h @@ -0,0 +1,80 @@ +// SPDX-FileCopyrightText: 2024 The OpenSn Authors +// SPDX-License-Identifier: MIT + +#pragma once + +#include "framework/parameters/input_parameters.h" +#include "modules/linear_boltzmann_solvers/lbs_problem/lbs_problem.h" +#include "linalg/Matrix.h" +#include "linalg/Vector.h" +#include "algo/manifold_interp/MatrixInterpolator.h" +#include "algo/manifold_interp/VectorInterpolator.h" +#include "rom_structs.h" +#include +#include + +namespace opensn +{ + +class ROMProblem : public Problem +{ +public: + /// Input parameters based construction. + explicit ROMProblem(const InputParameters& params); + + static InputParameters GetOptionsBlock(); + + static InputParameters GetInputParameters(); + static std::shared_ptr Create(const ParameterBlock& params); + + void SetOptions(const InputParameters& input); + /// Returns a reference to the solver options. + ROMOptions& GetOptions(); + + /// Save current Phi in the basis generator + void TakeSample(int id); + + /// Load snapshots and perform SVD + void MergePhase(int nsnaps); + + /// Load the params from file + void ReadParamMatrix(const std::string& filename); + + /// Calculate AU via sweeps + std::shared_ptr AssembleAU(); + + /// Sweep to form RHS + std::shared_ptr AssembleRHS(); + + /// Assemble the reduced system and save to file + void AssembleROM(std::shared_ptr& AU_, + std::shared_ptr& b_, + const std::string& Ar_filename, + const std::string& rhs_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); + + void InterpolateArAndRHS( + CAROM::Vector& desired_point, + std::shared_ptr& Ar_interp, + std::shared_ptr& rhs_interp); + +protected: + std::unique_ptr spatialbasis; + opensn::Vector b_; + std::unique_ptr Ar_interp_obj_ptr; + std::unique_ptr rhs_interp_obj_ptr; + std::shared_ptr lbs_problem; + ROMOptions options_; + +public: + std::vector param_points_; + int romRank; +}; + +} // namespace opensn diff --git a/src/rom/rom_structs.h b/src/rom/rom_structs.h new file mode 100644 index 0000000..064424d --- /dev/null +++ b/src/rom/rom_structs.h @@ -0,0 +1,21 @@ +// SPDX-FileCopyrightText: 2024 The OpenSn Authors +// SPDX-License-Identifier: MIT + +#pragma once + +#include "linalg/Vector.h" + +namespace opensn +{ + +struct ROMOptions +{ + int param_id = 0; + std::string phase = "offline"; + std::string param_file = ""; + std::unique_ptr new_point; + + ROMOptions() = default; +}; + +} // namespace opensn \ No newline at end of file diff --git a/src/rom/steady_state_rom_solver.cc b/src/rom/steady_state_rom_solver.cc new file mode 100644 index 0000000..decf673 --- /dev/null +++ b/src/rom/steady_state_rom_solver.cc @@ -0,0 +1,111 @@ +// SPDX-FileCopyrightText: 2025 The OpenSn Authors +// SPDX-License-Identifier: MIT + +#include "modules/linear_boltzmann_solvers/solvers/steady_state_solver.h" +#include "modules/linear_boltzmann_solvers/lbs_problem/lbs_problem.h" +#include "modules/linear_boltzmann_solvers/lbs_problem/iterative_methods/ags_linear_solver.h" +#include "framework/runtime.h" +#include "steady_state_rom_solver.h" +#include +#include +#include + +namespace opensn +{ + +InputParameters +SteadyStateROMSolver::GetInputParameters() +{ + InputParameters params = Solver::GetInputParameters(); + + params.SetGeneralDescription("Implementation of a steady state ROM solver. This solver calls the " + "across-groupset (AGS) solver offline and interfaces with libROM."); + params.ChangeExistingParamToOptional("name", "SteadyStateROMSolver"); + params.AddRequiredParameter>("problem", "An existing lbs problem"); + params.AddRequiredParameter>("rom_problem", "A ROM problem"); + + return params; +} + +SteadyStateROMSolver::SteadyStateROMSolver(const InputParameters& params) + : SteadyStateSourceSolver(params), + lbs_problem_(params.GetSharedPtrParam("problem")), + rom_problem_(params.GetSharedPtrParam("rom_problem")) +{ +} + +void +SteadyStateROMSolver::Initialize() +{ + lbs_problem_->Initialize(); +} + +void +SteadyStateROMSolver::Execute() +{ + auto& options = lbs_problem_->GetOptions(); + auto& rom_options = rom_problem_->GetOptions(); + + if (rom_options.phase == "offline") + { + auto& ags_solver = *lbs_problem_->GetAGSSolver(); + + auto start = std::chrono::high_resolution_clock::now(); + ags_solver.Solve(); + 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.txt"); + if (outfile.is_open()) { + outfile << elapsed.count() <UpdateFieldFunctions(); + + rom_problem_->TakeSample(rom_options.param_id); + } + if (rom_options.phase == "merge") + { + rom_problem_->MergePhase(rom_options.param_id); + } + if (rom_options.phase == "systems") + { + 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); + const std::string& rhs_filename = "data/rom_system_br_" + std::to_string(rom_options.param_id); + rom_problem_->AssembleROM(AU_, b_, Ar_filename, rhs_filename); + } + if (rom_options.phase == "online") + { + rom_problem_->ReadParamMatrix(rom_options.param_file); + + std::shared_ptr Ar_interp; + std::shared_ptr rhs_interp; + rom_problem_->SetupInterpolator(*rom_options.new_point); + + auto start = std::chrono::high_resolution_clock::now(); + + rom_problem_->InterpolateArAndRHS(*rom_options.new_point, Ar_interp, rhs_interp); + rom_problem_->SolveROM(Ar_interp, rhs_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.txt"); + if (outfile.is_open()) + { + outfile << elapsed.count() < +// SPDX-License-Identifier: MIT + +#pragma once + +#include "modules/linear_boltzmann_solvers/solvers/steady_state_solver.h" +#include "rom_problem.h" + +namespace opensn +{ + +class LBSProblem; + +class SteadyStateROMSolver : public SteadyStateSourceSolver +{ +protected: + std::shared_ptr lbs_problem_; + std::shared_ptr rom_problem_; + +public: + explicit SteadyStateROMSolver(const InputParameters& params); + + void Initialize(); + + void Execute(); + +public: + static InputParameters GetInputParameters(); + +}; + +} // namespace opensn \ No newline at end of file diff --git a/src/rom_py_app.cc b/src/rom_py_app.cc new file mode 100644 index 0000000..55f53c3 --- /dev/null +++ b/src/rom_py_app.cc @@ -0,0 +1,19 @@ +// SPDX-FileCopyrightText: 2025 The OpenSn Authors +// SPDX-License-Identifier: MIT + +#include "rom_py_app.h" +#include "py_wrappers.h" +#include "python/lib/console.h" + +using namespace opensn; + +namespace rompy +{ + +ROMApp::ROMApp(const mpi::Communicator& comm) + : opensnpy::PyApp(comm) +{ + opensnpy::console.BindModule(WrapROM); +} + +} // namespace rompy \ No newline at end of file diff --git a/src/rom_py_app.h b/src/rom_py_app.h new file mode 100644 index 0000000..3fce39d --- /dev/null +++ b/src/rom_py_app.h @@ -0,0 +1,20 @@ +// SPDX-FileCopyrightText: 2025 The OpenSn Authors +// SPDX-License-Identifier: MIT + +#pragma once + +#include "mpicpp-lite/mpicpp-lite.h" +#include "python/lib/py_app.h" + +namespace mpi = mpicpp_lite; + +namespace rompy +{ + +class ROMApp : public opensnpy::PyApp +{ +public: + explicit ROMApp(const mpi::Communicator& comm); +}; + +} // namespace rompy \ No newline at end of file diff --git a/src/test.cc b/src/test.cc deleted file mode 100644 index b25f81c..0000000 --- a/src/test.cc +++ /dev/null @@ -1,28 +0,0 @@ -#include "test.h" -#include "opensn/framework/logging/log.h" -#include "opensn/framework/object_factory.h" - -using namespace opensn; - -OpenSnRegisterObject(TestApp); - -InputParameters -TestApp::GetInputParameters() -{ - InputParameters params; - params.SetGeneralDescription("A dummy application for testing."); - params.AddOptionalParameter("name", "test_app", "A name for the application."); - return params; -} - -TestApp::TestApp(const InputParameters& params) : name_(params.GetParamValue("name")) -{ - opensn::log.Log() << "TestApp named " << name_ << " created."; -} - -std::shared_ptr -TestApp::Create(const ParameterBlock& params) -{ - auto& factory = opensn::ObjectFactory::GetInstance(); - return factory.Create("TestApp", params); -} diff --git a/src/test.h b/src/test.h deleted file mode 100644 index 2ea1caa..0000000 --- a/src/test.h +++ /dev/null @@ -1,15 +0,0 @@ -#pragma once -#include "opensn/framework/object.h" - -class TestApp : public opensn::Object -{ -public: - explicit TestApp(const opensn::InputParameters& params); - -private: - const std::string name_; - -public: - static opensn::InputParameters GetInputParameters(); - static std::shared_ptr Create(const opensn::ParameterBlock& params); -}; diff --git a/src/testapp_py.cc b/src/testapp_py.cc deleted file mode 100644 index e841cf1..0000000 --- a/src/testapp_py.cc +++ /dev/null @@ -1,12 +0,0 @@ -#include -#include "test.h" - -namespace py = pybind11; - -PYBIND11_MODULE(myapp, m) -{ - m.doc() = "Python bindings for my OpenSn-based TestApp"; - - py::class_>(m, "TestApp") - .def_static("Create", &TestApp::Create); -} diff --git a/tests/test.py b/tests/test.py deleted file mode 100644 index 693f29e..0000000 --- a/tests/test.py +++ /dev/null @@ -1,8 +0,0 @@ -# test.py -import sys -sys.path.insert(0, "build/python") # where CMake drops the .so - -import myapp - -app = myapp.TestApp.Create() -print("Successfully created:", app) From 9cad07399ad94db44de235ad7dda8914a9311ccb Mon Sep 17 00:00:00 2001 From: "quincy.huhn98" Date: Mon, 9 Feb 2026 09:19:11 -0600 Subject: [PATCH 4/4] Job Manager --- examples/2gcheckerboard/absorber_base.txt | 8 +- .../2gcheckerboard/base_2gcheckerboard.py | 80 ++-------- .../2gcheckerboard/checkerboard_problem_2g.py | 68 ++++++++ .../2gcheckerboard/run_rom_2gcheckerboard.py | 150 ++++++------------ examples/2gcheckerboard/scatterer_base.txt | 8 +- examples/checkerboard/base_checkerboard.py | 95 ++++------- examples/checkerboard/checkerboard_problem.py | 72 +++++++++ examples/checkerboard/run_rom_checkerboard.py | 137 ++++++---------- examples/python/plot_errors.py | 84 ---------- examples/reed/base_reed.py | 71 +++------ examples/reed/reed_problem.py | 58 +++++++ examples/reed/run_rom_reed.py | 115 ++++++-------- python/job_manager.py | 81 ++++++++++ {examples/python => python}/plotting.py | 8 +- python/rom_driver.py | 94 +++++++++++ {examples/python => python}/utils.py | 11 +- src/rom/rom_problem.cc | 2 +- src/rom/steady_state_rom_solver.cc | 4 +- 18 files changed, 588 insertions(+), 558 deletions(-) create mode 100644 examples/2gcheckerboard/checkerboard_problem_2g.py create mode 100644 examples/checkerboard/checkerboard_problem.py delete mode 100644 examples/python/plot_errors.py create mode 100644 examples/reed/reed_problem.py create mode 100644 python/job_manager.py rename {examples/python => python}/plotting.py (91%) create mode 100644 python/rom_driver.py rename {examples/python => python}/utils.py (91%) diff --git a/examples/2gcheckerboard/absorber_base.txt b/examples/2gcheckerboard/absorber_base.txt index 8da84b2..c6612e5 100644 --- a/examples/2gcheckerboard/absorber_base.txt +++ b/examples/2gcheckerboard/absorber_base.txt @@ -7,8 +7,8 @@ SIGMA_T_BEGIN SIGMA_T_END TRANSFER_MOMENTS_BEGIN -M_GPRIME_G_VAL 0 0 0 0.05 -M_GPRIME_G_VAL 0 1 0 0.45 -M_GPRIME_G_VAL 0 0 1 0.45 -M_GPRIME_G_VAL 0 1 1 0.05 +M_GFROM_GTO_VAL 0 0 0 0.05 +M_GFROM_GTO_VAL 0 1 0 0.45 +M_GFROM_GTO_VAL 0 0 1 0.45 +M_GFROM_GTO_VAL 0 1 1 0.05 TRANSFER_MOMENTS_END \ No newline at end of file diff --git a/examples/2gcheckerboard/base_2gcheckerboard.py b/examples/2gcheckerboard/base_2gcheckerboard.py index 00672ee..1ced4c8 100644 --- a/examples/2gcheckerboard/base_2gcheckerboard.py +++ b/examples/2gcheckerboard/base_2gcheckerboard.py @@ -24,65 +24,13 @@ if __name__ == "__main__": - try: - print("Absorber Sigma_a Parameter = {}".format(abs_1)) - param = abs_1 - except: - abs_1=10.0 - print("Absorber Sigma_a Nominal = {}".format(abs_1)) - - try: - print("Absorber Sigma_s Parameter = {}".format(scatt_1)) - param = scatt_1 - except: - scatt_1=10.0 - print("Absorber Sigma_s Nominal = {}".format(scatt_1)) - - - try: - print("Scatterer Sigma_a Parameter = {}".format(abs_2)) - param = abs_2 - except: - abs_2=0.0 - print("Scatterer Sigma_a Nominal = {}".format(abs_2)) - - try: - print("Scatterer Sigma_s Parameter = {}".format(scatt_2)) - param = scatt_2 - except: - scatt_2=1.0 - print("Scatterer Sigma_s Nominal = {}".format(scatt_2)) - - - try: - print("Source Parameter = {}".format(param_q)) - param = param_q - except: - param_q=1.0 - print("Source Nominal = {}".format(param_q)) - try: print("Parameter id = {}".format(p_id)) except: p_id=0 print("Parameter id = {}".format(p_id)) - try: - if phase == 0: - print("Offline Phase") - phase = "offline" - elif phase == 1: - print("Merge Phase") - phase = "merge" - elif phase == 2: - print("Systems Phase") - phase = "systems" - elif phase == 3: - print("Online Phase") - phase = "online" - except: - phase="offline" - print("Phase default to offline") + print("{} phase".format(phase)) # Check number of processors num_procs = 4 @@ -138,9 +86,7 @@ absorber = MultiGroupXS() absorber.LoadFromOpenSn("data/absorber.xs") - strength = [0.0 for _ in range(num_groups)] - src0 = VolumetricSource(block_ids=[0], group_strength=strength) - strength[0] = param_q + strength = [1.0, 0.0] src1 = VolumetricSource(block_ids=[1], group_strength=strength) # Setup Physics @@ -148,14 +94,14 @@ if phase == "online": rom_options={ - "param_id":0, + "param_id":pid, "phase":phase, "param_file":"data/params.txt", - "new_point":[scatt_1, abs_1] + "new_point":[p0, p1] } else: rom_options={ - "param_id":p_id, + "param_id":pid, "phase":phase } @@ -177,7 +123,7 @@ {"block_ids": [0, 1], "xs": scatterer}, {"block_ids": [2], "xs": absorber} ], - volumetric_sources=[src0, src1] + volumetric_sources=[src1] ) rom = ROMProblem(problem=phys,options=rom_options) @@ -186,7 +132,13 @@ ss_solver.Initialize() ss_solver.Execute() - if phase == "online": - phys.WriteFluxMoments("output/rom") - if phase == "offline": - phys.WriteFluxMoments("output/fom") \ No newline at end of file + try: + if phase == "online" and saveh5: + phys.WriteFluxMoments("output/rom_{}_".format(pid)) + if phase == "offline" and saveh5: + phys.WriteFluxMoments("output/fom_{}_".format(pid)) + except: + if phase == "online": + phys.WriteFluxMoments("output/rom") + if phase == "offline": + phys.WriteFluxMoments("output/fom") \ No newline at end of file diff --git a/examples/2gcheckerboard/checkerboard_problem_2g.py b/examples/2gcheckerboard/checkerboard_problem_2g.py new file mode 100644 index 0000000..c7c5eab --- /dev/null +++ b/examples/2gcheckerboard/checkerboard_problem_2g.py @@ -0,0 +1,68 @@ +from pathlib import Path +import utils +import plotting +import numpy as np + + +class CheckerboardProblem2G: + def __init__(self, workdir, nprocs=4, ntrain=100, ntest=10): + self.workdir = Path(workdir) + self.deck_path = self.workdir / "base_2gcheckerboard.py" + + self.nprocs = nprocs + self.ntrain = ntrain + self.ntest = ntest + + def sample_training(self): + bounds = [[0.5,1.0],[7.5,12.5]] + self.training_set = utils.sample_parameter_space(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) + + 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) + + + + def plot_results(self): + plotting.plot_sv(num_groups=2) + 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" + 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) + + 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/examples/2gcheckerboard/run_rom_2gcheckerboard.py b/examples/2gcheckerboard/run_rom_2gcheckerboard.py index c702a6d..3a6b38c 100644 --- a/examples/2gcheckerboard/run_rom_2gcheckerboard.py +++ b/examples/2gcheckerboard/run_rom_2gcheckerboard.py @@ -1,105 +1,45 @@ -import numpy as np -import sys, os -sys.path.insert(0, os.path.realpath("../python")) - -import plotting -import utils - -# Sampling training points -bounds = [[0.5,1.0],[7.5,12.5]] -num_params = 100 - -params = utils.sample_parameter_space(bounds, num_params) - -S_abs = [[0.0, 0.0], - [0.0, 0.0]] -sigma_t_scatt = [1.0, 1.0] - -# OFFLINE PHASE -phase = 0 - -for i in range(num_params): - S_scatt = [[1-params[i,0], params[i,0]], - [0.0, 1.0]] - utils.update_xs("scatterer_base.txt", "data/scatterer.xs", sigma_t_scatt, S_scatt) - - sigma_t_abs = [params[i,1], params[i,1]] - utils.update_xs("absorber_base.txt", "data/absorber.xs", sigma_t_abs, S_abs) - - cmd = "mpiexec -n=4 ../../build/rom_app_exec -i base_2gcheckerboard.py \ - -p phase={} -p p_id={}".format(phase,i) - utils.run_opensn(cmd) - -# MERGE PHASE -phase = 1 - -cmd = "mpiexec -n=4 ../../build/rom_app_exec -i base_2gcheckerboard.py -p phase={} -p p_id={}".format(phase, i) -utils.run_opensn(cmd) - -plotting.plot_sv(num_groups=2) - - -# SYSTEMS PHASE -phase = 2 -for i in range(num_params): - S_scatt = [[1-params[i,0], params[i,0]], - [0.0, 1.0]] - utils.update_xs("scatterer_base.txt", "data/scatterer.xs", sigma_t_scatt, S_scatt) - - sigma_t_abs = [params[i,1], params[i,1]] - utils.update_xs("absorber_base.txt", "data/absorber.xs", sigma_t_abs, S_abs) - - cmd = "mpiexec -n=4 ../../build/rom_app_exec -i base_2gcheckerboard.py \ - -p phase={} -p p_id={}".format(phase,i) - utils.run_opensn(cmd) - -np.savetxt("data/params.txt", params) - - -# Generate Test Data -test_scatt_1 = np.random.uniform(0.5,1.0,10) -test_abs_1 = np.random.uniform(7.5,12.5,10) -test = np.append(test_scatt_1[:,np.newaxis], test_abs_1[:,np.newaxis], axis=1) -np.savetxt("data/validation.txt", test) - -test = np.loadtxt("data/validation.txt") - -num_test = 10 -errors = [] -speedups = [] - -for i in range(num_test): - # ONLINE PHASE - S_scatt = [[1-test[i,0], test[i,0]], - [0.0, 1.0]] - utils.update_xs("scatterer_base.txt", "data/scatterer.xs", sigma_t_scatt, S_scatt) - - sigma_t_abs = [test[i,1], test[i,1]] - utils.update_xs("absorber_base.txt", "data/absorber.xs", sigma_t_abs, S_abs) - - phase = 3 - cmd = "mpiexec -n=4 ../../build/rom_app_exec -i base_2gcheckerboard.py \ - -p phase={} -p scatt_1={} -p abs_1={} -p p_id={}"\ - .format(phase,test[i][0],test[i][1],i) - utils.run_opensn(cmd) - rom_time = np.loadtxt("results/online_time.txt") - - # Reference FOM solution - phase = 0 - cmd = "mpiexec -n=4 ../../build/rom_app_exec -i base_2gcheckerboard.py \ - -p phase={} -p p_id={}".format(phase,i) - utils.run_opensn(cmd) - fom_time = np.loadtxt("results/offline_time.txt") - - plotting.plot_2d_flux("output/fom{}.h5", ranks=range(4), prefix="fom", pid=i) - plotting.plot_2d_flux("output/rom{}.h5", ranks=range(4), prefix="rom", pid=i) - - error = plotting.plot_2d_lineout(ranks=range(4), pid=i) - - errors.append(error) - speedups.append(fom_time/rom_time) - -print("Avg Error ", np.mean(errors)) -np.savetxt("results/errors.txt", errors) -print("Avg Speedup ", np.mean(speedups)) -np.savetxt("results/speedups.txt", speedups) \ No newline at end of file +# run_rom_checkerboard.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 checkerboard_problem_2g import CheckerboardProblem2G +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 = CheckerboardProblem2G(repo_root) + + run_pipeline(problem, repo_root, jm) + problem.plot_results() + + +if __name__ == "__main__": + main() diff --git a/examples/2gcheckerboard/scatterer_base.txt b/examples/2gcheckerboard/scatterer_base.txt index 8da84b2..c6612e5 100644 --- a/examples/2gcheckerboard/scatterer_base.txt +++ b/examples/2gcheckerboard/scatterer_base.txt @@ -7,8 +7,8 @@ SIGMA_T_BEGIN SIGMA_T_END TRANSFER_MOMENTS_BEGIN -M_GPRIME_G_VAL 0 0 0 0.05 -M_GPRIME_G_VAL 0 1 0 0.45 -M_GPRIME_G_VAL 0 0 1 0.45 -M_GPRIME_G_VAL 0 1 1 0.05 +M_GFROM_GTO_VAL 0 0 0 0.05 +M_GFROM_GTO_VAL 0 1 0 0.45 +M_GFROM_GTO_VAL 0 0 1 0.45 +M_GFROM_GTO_VAL 0 1 1 0.05 TRANSFER_MOMENTS_END \ No newline at end of file diff --git a/examples/checkerboard/base_checkerboard.py b/examples/checkerboard/base_checkerboard.py index 6f34f7c..a783fa6 100644 --- a/examples/checkerboard/base_checkerboard.py +++ b/examples/checkerboard/base_checkerboard.py @@ -25,64 +25,21 @@ if __name__ == "__main__": try: - print("Absorber Sigma_a Parameter = {}".format(abs_1)) - param = abs_1 + print("5 Parameter Case q = {}".format(p4)) + int_point = [p0,p1,p2,p3,p4] except: - abs_1=10.0 - print("Absorber Sigma_a Nominal = {}".format(abs_1)) + p4=1.0 + print("4 Parameter Case q nominal = {}".format(p4)) + int_point = [p0,p1,p2,p3] try: - print("Absorber Sigma_s Parameter = {}".format(scatt_1)) - param = scatt_1 + print("Parameter id = {}".format(pid)) except: - scatt_1=10.0 - print("Absorber Sigma_s Nominal = {}".format(scatt_1)) + pid=0 + print("Parameter id = {}".format(pid)) + print("{} phase".format(phase)) - try: - print("Scatterer Sigma_a Parameter = {}".format(abs_2)) - param = abs_2 - except: - abs_2=0.0 - print("Scatterer Sigma_a Nominal = {}".format(abs_2)) - - try: - print("Scatterer Sigma_s Parameter = {}".format(scatt_2)) - param = scatt_2 - except: - scatt_2=1.0 - print("Scatterer Sigma_s Nominal = {}".format(scatt_2)) - - - try: - print("Source Parameter = {}".format(param_q)) - param = param_q - except: - param_q=1.0 - print("Source Nominal = {}".format(param_q)) - - try: - print("Parameter id = {}".format(p_id)) - except: - p_id=0 - print("Parameter id = {}".format(p_id)) - - try: - if phase == 0: - print("Offline Phase") - phase = "offline" - elif phase == 1: - print("Merge Phase") - phase = "merge" - elif phase == 2: - print("Systems Phase") - phase = "systems" - elif phase == 3: - print("Online Phase") - phase = "online" - except: - phase="offline" - print("Phase default to offline") # Check number of processors num_procs = 4 @@ -131,19 +88,19 @@ ) grid.SetBlockIDFromLogicalVolume(vol_abs, 2, True) - scatt_t = abs_2 + scatt_2 + scatt_t = p3 + p1 num_groups = 1 scatterer = MultiGroupXS() - scatterer.CreateSimpleOneGroup(sigma_t=scatt_t, c=scatt_2/scatt_t) + scatterer.CreateSimpleOneGroup(sigma_t=scatt_t, c=p1/scatt_t) - abs_t = abs_1 + scatt_1 + abs_t = p2 + p0 absorber = MultiGroupXS() - absorber.CreateSimpleOneGroup(sigma_t=abs_t, c=scatt_1/abs_t) + absorber.CreateSimpleOneGroup(sigma_t=abs_t, c=p0/abs_t) strength = [0.0] src0 = VolumetricSource(block_ids=[0], group_strength=strength) - strength = [param_q] + strength = [p4] src1 = VolumetricSource(block_ids=[1], group_strength=strength) # Setup Physics @@ -152,14 +109,14 @@ if phase == "online": rom_options={ - "param_id":0, + "param_id":pid, "phase":phase, - "param_file":"data/interpolation_params.txt", - "new_point":[scatt_1, scatt_2, abs_1, abs_2, param_q] + "param_file":"data/params.txt", + "new_point":int_point } else: rom_options={ - "param_id":p_id, + "param_id":pid, "phase":phase } @@ -190,9 +147,13 @@ ss_solver.Initialize() ss_solver.Execute() - phys.ComputeBalance() - - if phase == "online": - phys.WriteFluxMoments("output/rom") - if phase == "offline": - phys.WriteFluxMoments("output/fom") \ No newline at end of file + try: + if phase == "online" and saveh5: + phys.WriteFluxMoments("output/rom_{}_".format(pid)) + if phase == "offline" and saveh5: + phys.WriteFluxMoments("output/fom_{}_".format(pid)) + except: + if phase == "online": + phys.WriteFluxMoments("output/rom") + if phase == "offline": + phys.WriteFluxMoments("output/fom") \ No newline at end of file diff --git a/examples/checkerboard/checkerboard_problem.py b/examples/checkerboard/checkerboard_problem.py new file mode 100644 index 0000000..7198df2 --- /dev/null +++ b/examples/checkerboard/checkerboard_problem.py @@ -0,0 +1,72 @@ +from pathlib import Path +import utils +import plotting +import numpy as np + + +class CheckerboardProblem: + def __init__(self, workdir, five_param=False, nprocs=4, ntrain=100, ntest=10): + self.workdir = Path(workdir) + self.deck_path = self.workdir / "base_checkerboard.py" + + self.nprocs = nprocs + self.ntrain = ntrain + self.ntest = ntest + self.five_param=five_param + + def sample_training(self): + if self.five_param: + bounds = [[0,5.0],[0.5,1.5],[7.5,12.5],[0.0,0.5],[0.1,1]] + else: + bounds = [[0,5.0],[0.5,1.5],[7.5,12.5],[0.0,0.5]] + + self.training_set = utils.sample_parameter_space(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.0,10) + test_scatt_2 = np.random.uniform(0.5,1.5,10) + test_abs_1 = np.random.uniform(7.5,12.5,10) + test_abs_2 = np.random.uniform(0.0,0.5,10) + + test = np.append(test_scatt_1[:,np.newaxis], test_scatt_2[:,np.newaxis], axis=1) + test = np.append(test, test_abs_1[:,np.newaxis], axis=1) + self.testing_set = np.append(test, test_abs_2[:,np.newaxis], axis=1) + + if self.five_param: + test_q = np.random.uniform(0.1,1,10) + self.testing_set = np.append(self.testing_set, test_q[:,np.newaxis], axis=1) + + params_path = self.workdir / "data" / "test_params.txt" + np.savetxt(str(params_path), self.testing_set) + + def update_xs(self): + print("Checkerbaord problem uses SimpleOneGroupXS, use run_pipeline_1g") + + + def plot_results(self): + plotting.plot_sv(num_groups=1) + 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" + 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) + + 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/examples/checkerboard/run_rom_checkerboard.py b/examples/checkerboard/run_rom_checkerboard.py index 9b26780..46bb3c8 100644 --- a/examples/checkerboard/run_rom_checkerboard.py +++ b/examples/checkerboard/run_rom_checkerboard.py @@ -1,92 +1,45 @@ -import numpy as np -import sys, os -sys.path.insert(0, os.path.realpath("../python")) - -import plotting -import utils - -# Sampling training points -bounds = [[0,5.0],[0.5,1.5],[7.5,12.5],[0.0,0.5],[0.1,1]] -num_params = 50 - -params = utils.sample_parameter_space(bounds, num_params) -#params = np.loadtxt("data/params.txt") - -#params = params[:num_params,:] -np.savetxt("data/interpolation_params.txt", params) -# OFFLINE PHASE -phase = 0 - -for i in range(num_params): - cmd = "mpiexec -n=4 ../../build/rom_app_exec -i base_checkerboard.py \ - -p phase={} -p scatt_1={} -p scatt_2={} -p abs_1={} -p abs_2={} -p param_q={} -p p_id={}"\ - .format(phase,params[i][0],params[i][1],params[i][2],params[i][3],params[i][4],i) - utils.run_opensn(cmd) - -# MERGE PHASE -phase = 1 - -cmd = "mpiexec -n=4 ../../build/rom_app_exec -i base_checkerboard.py -p phase={} -p p_id={}".format(phase, num_params-1) -utils.run_opensn(cmd) - -plotting.plot_sv(num_groups=1) - - -# SYSTEMS PHASE -phase = 2 -for i in range(num_params): - cmd = "mpiexec -n=4 ../../build/rom_app_exec -i base_checkerboard.py \ - -p phase={} -p scatt_1={} -p scatt_2={} -p abs_1={} -p abs_2={} -p param_q={} -p p_id={}"\ - .format(phase,params[i][0],params[i][1],params[i][2],params[i][3],params[i][4],i) - utils.run_opensn(cmd) - -np.savetxt("data/params.txt", params) - - -# Generate Test Data -test_scatt_1 = np.random.uniform(0,5.0,10) -test_scatt_2 = np.random.uniform(0.5,1.5,10) -test_abs_1 = np.random.uniform(7.5,12.5,10) -test_abs_2 = np.random.uniform(0.0,0.5,10) -test_q = np.random.uniform(0.1,1,10) -test = np.append(test_scatt_1[:,np.newaxis], test_scatt_2[:,np.newaxis], axis=1) -test = np.append(test, test_abs_1[:,np.newaxis], axis=1) -test = np.append(test, test_abs_2[:,np.newaxis], axis=1) -test = np.append(test, test_q[:,np.newaxis], axis=1) -np.savetxt("data/validation.txt", test) -# test = np.loadtxt("data/validation.txt") - -num_test = 10 -errors = [] -speedups = [] - -for i in range(num_test): - # ONLINE PHASE - phase = 3 - cmd = "mpiexec -n=4 ../../build/rom_app_exec -i base_checkerboard.py \ - -p phase={} -p scatt_1={} -p scatt_2={} -p abs_1={} -p abs_2={} -p param_q={} -p p_id={}"\ - .format(phase,test[i][0],test[i][1],test[i][2],test[i][3],test[i][4],i) - utils.run_opensn(cmd) - rom_time = np.loadtxt("results/online_time.txt") - - # Reference FOM solution - phase = 0 - cmd = "mpiexec -n=4 ../../build/rom_app_exec -i base_checkerboard.py \ - -p phase={} -p scatt_1={} -p scatt_2={} -p abs_1={} -p abs_2={} -p param_q={} -p p_id={}"\ - .format(phase,test[i][0],test[i][1],test[i][2],test[i][3],test[i][4],i) - utils.run_opensn(cmd) - fom_time = np.loadtxt("results/offline_time.txt") - - plotting.plot_2d_flux("output/fom{}.h5", ranks=range(4), prefix="fom", pid=i) - plotting.plot_2d_flux("output/rom{}.h5", ranks=range(4), prefix="rom", pid=i) - - error = plotting.plot_2d_lineout(ranks=range(4), pid=i) - - errors.append(error) - speedups.append(fom_time/rom_time) - - -print("Avg Error ", np.mean(errors)) -np.savetxt("results/errors.txt", errors) -print("Avg Speedup ", np.mean(speedups)) -np.savetxt("results/speedups.txt", speedups) \ No newline at end of file +# run_rom_checkerboard.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 checkerboard_problem import CheckerboardProblem +from rom_driver import run_pipeline_1g + +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 = CheckerboardProblem(repo_root) + + run_pipeline_1g(problem, repo_root, jm) + problem.plot_results() + + +if __name__ == "__main__": + main() diff --git a/examples/python/plot_errors.py b/examples/python/plot_errors.py deleted file mode 100644 index 471b1e3..0000000 --- a/examples/python/plot_errors.py +++ /dev/null @@ -1,84 +0,0 @@ -import numpy as np -import matplotlib.pyplot as plt - -problem = "../checkerboard" - -# List of sample sizes -sample_sizes = [100, 150, 200, 250, 300, 350, 400, 450, 500] - -# Collect error data -all_errors = [] -for size in sample_sizes: - filename = f"{problem}/results/errors_{size}.txt" - data = np.loadtxt(filename) - all_errors.append(data) - -# Make the boxplot -plt.figure(figsize=(8, 5)) -plt.boxplot(all_errors, positions=sample_sizes, widths=20, showfliers=False) - -plt.xlabel("Training Set Size") -plt.ylabel("$L_2$ Error") -plt.yscale("log") -plt.grid(True, linestyle="--", alpha=0.6) - -plt.savefig(f'{problem}/results/errors_set_size.jpg') -plt.close() - -# Collect speedup data -all_speedups = [] -for size in sample_sizes: - filename = f"{problem}/results/speedups_{size}.txt" - data = np.loadtxt(filename) - all_speedups.append(data) - -# Make the boxplot -plt.figure(figsize=(8, 5)) -plt.boxplot(all_speedups, positions=sample_sizes, widths=20, showfliers=False) - -plt.xlabel("Training Set Size") -plt.ylabel("Speedup") -plt.yscale("log") -plt.grid(True, linestyle="--", alpha=0.6) - -plt.savefig(f'{problem}/results/speedups_set_size.jpg') - -# List of ranks -ranks = [5, 10, 15, 20] - -# Collect error data -all_errors = [] -for rank in ranks: - filename = f"{problem}/results/errors_{rank}.txt" - data = np.loadtxt(filename) - all_errors.append(data) - -# Make the boxplot -plt.figure(figsize=(8, 5)) -plt.boxplot(all_errors, positions=ranks, widths=2, showfliers=False) - -plt.xlabel("Rank") -plt.ylabel("$L_2$ Error") -plt.yscale("log") -plt.grid(True, linestyle="--", alpha=0.6) - -plt.savefig(f'{problem}/results/errors_rank.jpg') -plt.close() - -# Collect speedup data -all_speedups = [] -for rank in ranks: - filename = f"{problem}/results/speedups_{rank}.txt" - data = np.loadtxt(filename) - all_speedups.append(data) - -# Make the boxplot -plt.figure(figsize=(8, 5)) -plt.boxplot(all_speedups, positions=ranks, widths=2, showfliers=False) - -plt.xlabel("Rank") -plt.ylabel("Speedup") -plt.yscale("log") -plt.grid(True, linestyle="--", alpha=0.6) - -plt.savefig(f'{problem}/results/speedups_rank.jpg') \ No newline at end of file diff --git a/examples/reed/base_reed.py b/examples/reed/base_reed.py index 57622ce..6d4fa91 100644 --- a/examples/reed/base_reed.py +++ b/examples/reed/base_reed.py @@ -22,48 +22,12 @@ if __name__ == "__main__": try: - print("Scattering Parameter = {}".format(scatt)) - param = scatt - except: - scatt=0.0 - print("Scattering Nominal = {}".format(scatt)) - - try: - print("Cross Section Parameter = {}".format(sigma_t)) - param = sigma_t - except: - sigma_t=1.0 - print("Cross Section Nominal = {}".format(sigma_t)) - - try: - print("Source Parameter = {}".format(param_q)) - param = param_q - except: - param_q=1.0 - print("Source Nominal = {}".format(param_q)) - - try: - print("Parameter id = {}".format(p_id)) + print("Parameter id = {}".format(pid)) except: p_id=0 - print("Parameter id = {}".format(p_id)) + print("Parameter id = {}".format(pid)) - try: - if phase == 0: - print("Offline Phase") - phase = "offline" - elif phase == 1: - print("Merge Phase") - phase = "merge" - elif phase == 2: - print("Systems Phase") - phase = "systems" - elif phase == 3: - print("Online Phase") - phase = "online" - except: - phase="offline" - print("Phase default to offline") + print("{} phase".format(phase)) # Create Mesh widths = [2., 1., 2., 1., 2.] @@ -88,8 +52,8 @@ z_min = z_max # Add cross sections to materials - total = [50., 5., 0., sigma_t, sigma_t] - c = [0., 0., 0., scatt, scatt] + total = [50., 5., 0., 1., 1.] + c = [0., 0., 0., p0, p0] xs_map = len(total) * [None] for imat in range(Nmat): xs_ = MultiGroupXS() @@ -100,7 +64,7 @@ # Create sources in 1st and 4th materials src0 = VolumetricSource(block_ids=[0], group_strength=[50.]) - src1 = VolumetricSource(block_ids=[3], group_strength=[param_q]) + src1 = VolumetricSource(block_ids=[3], group_strength=[p1]) # Angular Quadrature gl_quad = GLProductQuadrature1DSlab(n_polar=128, scattering_order=0) @@ -131,14 +95,14 @@ if phase == "online": rom_options = { - "param_id": 0, + "param_id": pid, "phase": phase, "param_file": "data/params.txt", - "new_point": [scatt, param_q] + "new_point": [p0, p1] } else: rom_options = { - "param_id": p_id, + "param_id": pid, "phase": phase } @@ -149,10 +113,15 @@ ss_solver.Initialize() ss_solver.Execute() - # compute particle balance - phys.ComputeBalance() + try: + if phase == "online" and saveh5: + phys.WriteFluxMoments("output/rom_{}_".format(pid)) + if phase == "offline" and saveh5: + phys.WriteFluxMoments("output/fom_{}_".format(pid)) + except: + if phase == "online": + phys.WriteFluxMoments("output/rom") + if phase == "offline": + phys.WriteFluxMoments("output/fom") + - if phase == "online": - phys.WriteFluxMoments("output/rom") - if phase == "offline": - phys.WriteFluxMoments("output/fom") \ No newline at end of file diff --git a/examples/reed/reed_problem.py b/examples/reed/reed_problem.py new file mode 100644 index 0000000..aa2b4c2 --- /dev/null +++ b/examples/reed/reed_problem.py @@ -0,0 +1,58 @@ +from pathlib import Path +import utils +import plotting +import numpy as np + + +class ReedProblem: + def __init__(self, workdir, nprocs=2, ntrain=100, ntest=10): + self.workdir = Path(workdir) + self.deck_path = self.workdir / "base_reed.py" + + self.nprocs = nprocs + self.ntrain = ntrain + self.ntest = ntest + + def sample_training(self): + bounds = [[0.0,1.0],[0.0,1.0]] + + self.training_set = utils.sample_parameter_space(bounds, self.ntrain) + + params_path = self.workdir / "data" / "params.txt" + np.savetxt(str(params_path), self.training_set) + + + def sample_testing(self): + self.testing_set = np.random.uniform(0,1,[self.ntest,2]) + + params_path = self.workdir / "data" / "test_params.txt" + np.savetxt(str(params_path), self.testing_set) + + def update_xs(self): + print("Reed problem uses SimpleOneGroupXS, use run_pipeline_1g") + + + def plot_results(self): + plotting.plot_sv(num_groups=1) + 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_flux( + str(output_dir / ("fom_{}_".format(i) + "{}.h5")), + str(output_dir / ("rom_{}_".format(i) + "{}.h5")), + ranks=range(self.nprocs), + pid=i) + + errors.append(error) + 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/examples/reed/run_rom_reed.py b/examples/reed/run_rom_reed.py index a1274c2..7839ccd 100644 --- a/examples/reed/run_rom_reed.py +++ b/examples/reed/run_rom_reed.py @@ -1,70 +1,45 @@ -import numpy as np -import sys, os -sys.path.insert(0, os.path.realpath("../python")) - -import plotting -import utils - -# Sampling training points -bounds = [[0.0,1.0],[0.0,1.0]] -num_params = 100 - -params = utils.sample_parameter_space(bounds, num_params) - -# OFFLINE PHASE -phase = 0 - -for i, param in enumerate(params): - cmd = "mpiexec -n 2 ../../build/rom_app_exec -i base_reed.py -p phase={} -p scatt={} -p param_q={} -p p_id={}"\ - .format(phase,param[0], param[1], i) - utils.run_opensn(cmd) - -# MERGE PHASE -phase = 1 - -cmd = "mpiexec -n 2 ../../build/rom_app_exec -i base_reed.py -p phase={} -p p_id={}".format(phase, i) -utils.run_opensn(cmd) - -plotting.plot_sv(num_groups=1) - - -# SYSTEMS PHASE -phase = 2 - -for i, param in enumerate(params): - cmd = "mpiexec -n 2 ../../build/rom_app_exec -i base_reed.py -p phase={} -p scatt={} -p param_q={} -p p_id={}"\ - .format(phase,param[0], param[1], i) - utils.run_opensn(cmd) - -np.savetxt("data/params.txt", params) - -# Generate Test Data -test = np.random.uniform(0,1,[10,2]) - -errors = [] -speedups = [] - -for i, param in enumerate(test): - # ONLINE PHASE - phase = 3 - cmd = "mpiexec -n 2 ../../build/rom_app_exec -i base_reed.py -p phase={} -p scatt={} -p param_q={} -p p_id={}"\ - .format(phase,param[0],param[1], i) - utils.run_opensn(cmd) - rom_time = np.loadtxt("results/online_time.txt") - - # Reference FOM solution - phase = 0 - cmd = "mpiexec -n 2 ../../build/rom_app_exec -i base_reed.py -p phase={} -p scatt={} -p param_q={} -p p_id={}"\ - .format(phase,param[0],param[1], i) - utils.run_opensn(cmd) - fom_time = np.loadtxt("results/offline_time.txt") - - error = plotting.plot_1d_flux("output/fom{}.h5", "output/rom{}.h5", ranks=range(2), pid=i) - - errors.append(error) - speedups.append(fom_time/rom_time) - -print("Avg Error ", np.mean(errors)) -np.savetxt("results/errors.txt", errors) -print("Avg Speedup ", np.mean(speedups)) -np.savetxt("results/speedups.txt", speedups) \ No newline at end of file +# run_rom_reed.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 reed_problem import ReedProblem +from rom_driver import run_pipeline_1g + +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 = ReedProblem(repo_root) + + run_pipeline_1g(problem, repo_root, jm) + problem.plot_results() + + +if __name__ == "__main__": + main() diff --git a/python/job_manager.py b/python/job_manager.py new file mode 100644 index 0000000..7c1d6c5 --- /dev/null +++ b/python/job_manager.py @@ -0,0 +1,81 @@ +import os +import subprocess + + +class JobManager: + """ + System-aware launcher for OpenSn. + Decides how to run OpenSn (srun / mpirun / direct) based on environment. + """ + + def __init__(self, system="auto", opensn_exe=None): + self.system = system + self.opensn_exe = opensn_exe + + # ------------------------- + # System detection + # ------------------------- + def detect_system(self): + if self.system != "auto": + return self.system + + # Slurm environment detection + if os.environ.get("SLURM_JOB_ID") or os.environ.get("SLURM_CLUSTER_NAME"): + return "slurm" + + return "local" + + # ------------------------- + # Command construction + # ------------------------- + def build_command(self, input_file, nprocs=1, launcher_args=None, opensn_args=None): + exe = self.opensn_exe + + tail = [exe, "-i", str(input_file)] + if opensn_args: + tail.extend(opensn_args) + + if nprocs > 1: + launcher = "mpirun" + else: + launcher = "none" + + cmd = [] + + if launcher == "none": + cmd = tail + + elif launcher == "mpirun": + cmd = ["mpirun", "-np", str(nprocs)] + if launcher_args: + cmd.extend(launcher_args) + cmd.extend(tail) + + return cmd + + # ------------------------- + # Execution + # ------------------------- + def run( + self, + input_file, + nprocs=1, + workdir=None, + launcher_args=None, + opensn_args=None, + stream_output=True, + check=False, + ): + cmd = self.build_command( + input_file=input_file, + nprocs=nprocs, + launcher_args=launcher_args, + opensn_args=opensn_args, + ) + print("OPENSN CMD:", cmd, flush=True) + if stream_output: + p = subprocess.run(cmd,cwd=workdir) + else: + p = subprocess.run(cmd,cwd=workdir, + stdout=subprocess.PIPE, stderr=subprocess.PIPE, + text=True) \ No newline at end of file diff --git a/examples/python/plotting.py b/python/plotting.py similarity index 91% rename from examples/python/plotting.py rename to python/plotting.py index 5f703ab..acfc876 100644 --- a/examples/python/plotting.py +++ b/python/plotting.py @@ -6,7 +6,7 @@ def plot_2d_flux(file_pattern, ranks, moment=0, prefix="fom", grid_res=200, pid=0): - """Create smooth full-color plots (not scatter) for each energy group.""" + """Create smooth full-color plots for each energy group.""" xs, ys, vals, G = load_2d_flux(file_pattern, ranks, moment=moment) for g in range(G): @@ -41,10 +41,10 @@ def plot_2d_flux(file_pattern, ranks, moment=0, prefix="fom", grid_res=200, pid= plt.savefig(outpath, dpi=200) plt.close() -def plot_2d_lineout(ranks, y_target=4.0, moment=0, grid_res=200, pid=0): +def plot_2d_lineout(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("output/rom{}.h5", ranks, moment=moment) - xs_, ys_, vals_, G = load_2d_flux("output/fom{}.h5", ranks, moment=moment) + 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): # Create regular grid diff --git a/python/rom_driver.py b/python/rom_driver.py new file mode 100644 index 0000000..41a68ee --- /dev/null +++ b/python/rom_driver.py @@ -0,0 +1,94 @@ +from pathlib import Path +import numpy as np + + +def ensure_problem_dirs(problem_root): + paths = { + "root": problem_root, + "data": problem_root / "data", + "basis": problem_root / "basis", + "output": problem_root / "output", + "results": problem_root / "results", + } + for p in paths.values(): + p.mkdir(parents=True, exist_ok=True) + return paths + + +def make_opensn_args(phase, pid, pvec, save_h5=False): + args = ["-p", "phase={}".format(repr(phase)),"-p", "pid={}".format(pid)] + + if pvec is not None: + for i, v in enumerate(pvec): + args.extend(["-p", "p{}={}".format(i, float(v))]) + + if save_h5: + args.extend(["-p", "saveh5=True"]) + + return args + + +def _run_one(problem, jm, workdir, phase, pid, pvec=None, save_h5=False): + opensn_args = make_opensn_args(phase=phase, pid=pid, pvec=pvec, save_h5=save_h5) + + res = jm.run( + input_file=str(problem.deck_path), + nprocs=problem.nprocs, + workdir=str(workdir), + opensn_args=opensn_args, + stream_output=True, + ) + +def _run_many(problem, jm, workdir, phase, dataset, save_h5=False): + for pid, pvec in enumerate(dataset): + problem.update_xs(pvec) + _run_one(problem, jm, workdir, phase=phase, pid=pid, pvec=pvec, save_h5=save_h5) + +def _run_many_1g(problem, jm, workdir, phase, dataset, save_h5=False): + for pid, pvec in enumerate(dataset): + _run_one(problem, jm, workdir, phase=phase, pid=pid, pvec=pvec, save_h5=save_h5) + + +def run_pipeline(problem, repo_root, jm): + paths = ensure_problem_dirs(Path(repo_root)) + + problem.sample_training() + + # OFFLINE training + _run_many(problem, jm, workdir=paths["root"], phase="offline", dataset=problem.training_set) + + # MERGE + _run_one(problem, jm, workdir=paths["root"], phase="merge", pid=problem.ntrain - 1, pvec=np.ones_like(problem.training_set[0])) + + # SYSTEMS + _run_many(problem, jm, workdir=paths["root"], phase="systems", dataset=problem.training_set) + + problem.sample_testing() + + # OFFLINE testing (save HDF5) + _run_many(problem, jm, workdir=paths["root"], phase="offline", dataset=problem.testing_set, save_h5=True) + + # ONLINE testing (save HDF5) + _run_many(problem, jm, workdir=paths["root"], phase="online", dataset=problem.testing_set, save_h5=True) + +def run_pipeline_1g(problem, repo_root, jm): + paths = ensure_problem_dirs(Path(repo_root)) + + problem.sample_training() + + # OFFLINE training + _run_many_1g(problem, jm, workdir=paths["root"], phase="offline", dataset=problem.training_set) + + # MERGE + _run_one(problem, jm, workdir=paths["root"], phase="merge", pid=problem.ntrain - 1, pvec=np.ones_like(problem.training_set[0])) + + # SYSTEMS + _run_many_1g(problem, jm, workdir=paths["root"], phase="systems", dataset=problem.training_set) + + problem.sample_testing() + + # OFFLINE testing (save HDF5) + _run_many_1g(problem, jm, workdir=paths["root"], phase="offline", dataset=problem.testing_set, save_h5=True) + + # ONLINE testing (save HDF5) + _run_many_1g(problem, jm, workdir=paths["root"], phase="online", dataset=problem.testing_set, save_h5=True) diff --git a/examples/python/utils.py b/python/utils.py similarity index 91% rename from examples/python/utils.py rename to python/utils.py index 8a16c49..c9e0c43 100644 --- a/examples/python/utils.py +++ b/python/utils.py @@ -2,15 +2,6 @@ import subprocess import h5py -def run_opensn(cmd): - args = cmd.split(" ") - print(args) - process = subprocess.Popen(args, stdout=subprocess.PIPE, stderr=subprocess.PIPE, text=True) - output, errors = process.communicate() - print("Output:", output) - print("Errors:", errors) - - def load_2d_flux(file_pattern, ranks, moment=0): """Load (x, y, flux) grouped by energy group from HDF5 files.""" with h5py.File(file_pattern.format(ranks[0]), "r") as f0: @@ -117,7 +108,7 @@ def update_xs(in_file, out_file, sigma_t_vec, S): for gprime in range(G): for g in range(G): val = float(S[gprime][g]) - new_tm.append(f"M_GPRIME_G_VAL 0 {gprime} {g} {val:.12g}\n") + new_tm.append(f"M_GFROM_GTO_VAL 0 {gprime} {g} {val:.12g}\n") lines[tb+1:te] = new_tm diff --git a/src/rom/rom_problem.cc b/src/rom/rom_problem.cc index 9849890..1a8418b 100644 --- a/src/rom/rom_problem.cc +++ b/src/rom/rom_problem.cc @@ -113,7 +113,7 @@ ROMProblem::MergePhase(int nsnaps) auto full_dim = num_local_nodes * num_moments * num_groups; CAROM::Options options(group_dim, max_num_snapshots, update_right_SV); - double tol = 1e-20; + double tol = 1e-8; romRank = 20; options.setSingularValueTol(tol); options.setMaxBasisDimension(romRank); diff --git a/src/rom/steady_state_rom_solver.cc b/src/rom/steady_state_rom_solver.cc index decf673..782861c 100644 --- a/src/rom/steady_state_rom_solver.cc +++ b/src/rom/steady_state_rom_solver.cc @@ -57,7 +57,7 @@ SteadyStateROMSolver::Execute() std::chrono::duration elapsed = end - start; if (opensn::mpi_comm.rank() == 0) { - std::ofstream outfile("results/offline_time.txt"); + std::ofstream outfile("results/offline_time_" + std::to_string(rom_options.param_id) + ".txt"); if (outfile.is_open()) { outfile << elapsed.count() <