diff --git a/.github/workflows/regression.yaml b/.github/workflows/regression.yaml new file mode 100644 index 0000000..44076c0 --- /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 -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 examples/reed && python run_rom_reed.py \ No newline at end of file diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..dfe896d --- /dev/null +++ b/.gitignore @@ -0,0 +1,94 @@ +.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 + +examples/*/basis/* +examples/*/data/* +examples/*/output/* +examples/*/results/* \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..f827cf6 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,111 @@ +cmake_minimum_required(VERSION 3.14) +project(OpenSnApp LANGUAGES C CXX) + +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 + CommonCore + CommonDataModel + IOLegacy + IOCore + IOXML + ParallelCore + IOParallelXML + FiltersCore + IOEnSight + IOExodus +) + +include(CMakeFindDependencyMacro) + +find_package(OpenSn REQUIRED) + +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(romapp MODULE + src/main.cc + src/rom_py_app.cc + src/rom.cc + ${ROM_SRCS} +) + +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(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(romapp PROPERTIES + PREFIX "" + OUTPUT_NAME "romapp" + LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/python +) + +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 new file mode 100644 index 0000000..397544b --- /dev/null +++ b/src/main.cc @@ -0,0 +1,29 @@ +// 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); + + PetscCall(PetscInitializeNoArguments()); // NOLINT(bugprone-casting-through-void) + + int retval = EXIT_SUCCESS; + try + { + py::scoped_interpreter guard{}; + 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