diff --git a/test/python/modules/linear_boltzmann_solvers/transport_transient/tests.json b/test/python/modules/linear_boltzmann_solvers/transport_transient/tests.json index 4b1fef5de..1e777e719 100644 --- a/test/python/modules/linear_boltzmann_solvers/transport_transient/tests.json +++ b/test/python/modules/linear_boltzmann_solvers/transport_transient/tests.json @@ -90,6 +90,97 @@ } ] }, + { + "file": "transient_time_dependent_boundary_cbc.py", + "comment": "Time-dependent isotropic and arbitrary boundary conditions (CBC)", + "num_procs": 4, + "checks": [ + { + "type": "FloatCompare", + "key": "TD_BOUNDARY_CASE01_ISO_PULSE_EQUIV_PASS", + "wordnum": 1, + "gold": 1, + "abs_tol": 1e-12 + }, + { + "type": "FloatCompare", + "key": "TD_BOUNDARY_CASE02_DELAYED_ISO_EQUIV_PASS", + "wordnum": 1, + "gold": 1, + "abs_tol": 1e-12 + }, + { + "type": "FloatCompare", + "key": "TD_BOUNDARY_CASE03_ISO_WINDOW_EQUIV_PASS", + "wordnum": 1, + "gold": 1, + "abs_tol": 1e-12 + }, + { + "type": "FloatCompare", + "key": "TD_BOUNDARY_CASE04_INACTIVE_ISO_ZERO_PASS", + "wordnum": 1, + "gold": 1, + "abs_tol": 1e-12 + }, + { + "type": "FloatCompare", + "key": "TD_BOUNDARY_CASE05_INACTIVE_ARBITRARY_ZERO_PASS", + "wordnum": 1, + "gold": 1, + "abs_tol": 1e-12 + }, + { + "type": "FloatCompare", + "key": "TD_BOUNDARY_CASE06_ARBITRARY_CALLED_EACH_STEP_PASS", + "wordnum": 1, + "gold": 1, + "abs_tol": 1e-12 + }, + { + "type": "FloatCompare", + "key": "TD_BOUNDARY_CASE07_MULTIGROUP_EQUIV_PASS", + "wordnum": 1, + "gold": 1, + "abs_tol": 1e-12 + }, + { + "type": "FloatCompare", + "key": "TD_BOUNDARY_CASE08_MULTI_ARBITRARY_CONSTANT_EQUIV_PASS", + "wordnum": 1, + "gold": 1, + "abs_tol": 1e-12 + }, + { + "type": "FloatCompare", + "key": "TD_BOUNDARY_CASE09_MULTI_ARBITRARY_WINDOWS_PASS", + "wordnum": 1, + "gold": 1, + "abs_tol": 1e-12 + }, + { + "type": "FloatCompare", + "key": "TD_BOUNDARY_CASE10_REFLECTING_ARBITRARY_ISO_PASS", + "wordnum": 1, + "gold": 1, + "abs_tol": 1e-12 + }, + { + "type": "FloatCompare", + "key": "TD_BOUNDARY_CASE11_INVALID_ISO_BOUNDS_PASS", + "wordnum": 1, + "gold": 1, + "abs_tol": 1e-12 + }, + { + "type": "FloatCompare", + "key": "TD_BOUNDARY_CASE12_ARBITRARY_REJECTS_BOUNDS_PASS", + "wordnum": 1, + "gold": 1, + "abs_tol": 1e-12 + } + ] + }, { "file": "transient_2d_mixed_time_dependent_bcs.py", "comment": "2D mixed time-dependent boundary conditions", @@ -104,6 +195,20 @@ } ] }, + { + "file": "transient_2d_mixed_time_dependent_bcs_cbc.py", + "comment": "2D mixed time-dependent boundary conditions (CBC)", + "num_procs": 4, + "checks": [ + { + "type": "FloatCompare", + "key": "TD_MIXED_2D_PASS", + "wordnum": 1, + "gold": 1, + "abs_tol": 1e-12 + } + ] + }, { "file": "transient_3d_mixed_time_dependent_bcs.py", "comment": "3D mixed time-dependent boundary conditions", @@ -118,6 +223,20 @@ } ] }, + { + "file": "transient_3d_mixed_time_dependent_bcs_cbc.py", + "comment": "3D mixed time-dependent boundary conditions (CBC)", + "num_procs": 4, + "checks": [ + { + "type": "FloatCompare", + "key": "TD_MIXED_3D_PASS", + "wordnum": 1, + "gold": 1, + "abs_tol": 1e-12 + } + ] + }, { "file": "transient_2d_adjacent_reflecting_arbitrary_bcs.py", "comment": "2D adjacent reflecting and arbitrary boundary conditions", @@ -132,6 +251,20 @@ } ] }, + { + "file": "transient_2d_adjacent_reflecting_arbitrary_bcs_cbc.py", + "comment": "2D adjacent reflecting and arbitrary boundary conditions (CBC)", + "num_procs": 4, + "checks": [ + { + "type": "FloatCompare", + "key": "TD_ADJACENT_REFLECTING_ARBITRARY_PASS", + "wordnum": 1, + "gold": 1, + "abs_tol": 1e-12 + } + ] + }, { "file": "transient_3d_adjacent_reflecting_arbitrary_bcs.py", "comment": "3D adjacent reflecting and arbitrary boundary conditions", @@ -146,6 +279,20 @@ } ] }, + { + "file": "transient_3d_adjacent_reflecting_arbitrary_bcs_cbc.py", + "comment": "3D adjacent reflecting and arbitrary boundary conditions (CBC)", + "num_procs": 4, + "checks": [ + { + "type": "FloatCompare", + "key": "TD_ADJACENT_REFLECTING_ARBITRARY_3D_PASS", + "wordnum": 1, + "gold": 1, + "abs_tol": 1e-12 + } + ] + }, { "file": "transient_2d_opposing_reflecting_arbitrary_bcs.py", "comment": "2D opposing reflecting and arbitrary boundary conditions", @@ -160,6 +307,20 @@ } ] }, + { + "file": "transient_2d_opposing_reflecting_arbitrary_bcs_cbc.py", + "comment": "2D opposing reflecting and arbitrary boundary conditions (CBC)", + "num_procs": 4, + "checks": [ + { + "type": "FloatCompare", + "key": "TD_OPPOSING_REFLECTING_ARBITRARY_PASS", + "wordnum": 1, + "gold": 1, + "abs_tol": 1e-12 + } + ] + }, { "file": "transient_3d_opposing_reflecting_arbitrary_bcs.py", "comment": "3D opposing reflecting and arbitrary boundary conditions", @@ -174,6 +335,20 @@ } ] }, + { + "file": "transient_3d_opposing_reflecting_arbitrary_bcs_cbc.py", + "comment": "3D opposing reflecting and arbitrary boundary conditions (CBC)", + "num_procs": 4, + "checks": [ + { + "type": "FloatCompare", + "key": "TD_OPPOSING_REFLECTING_ARBITRARY_3D_PASS", + "wordnum": 1, + "gold": 1, + "abs_tol": 1e-12 + } + ] + }, { "file": "transient_init_from_steady_restart.py", "comment": "Transient initial condition loaded directly from a steady-state restart", @@ -188,6 +363,20 @@ } ] }, + { + "file": "transient_init_from_steady_restart_cbc.py", + "comment": "Transient initial condition loaded directly from a steady-state restart (CBC)", + "num_procs": 4, + "checks": [ + { + "type": "FloatCompare", + "key": "STEADY_RESTART_IC_PASS", + "wordnum": 1, + "gold": 1, + "abs_tol": 1e-12 + } + ] + }, { "file": "transient_precursors_from_steady_restart.py", "comment": "Delayed precursor transient initialized from a steady-state restart", @@ -202,6 +391,20 @@ } ] }, + { + "file": "transient_precursors_from_steady_restart_cbc.py", + "comment": "Delayed precursor transient initialized from a steady-state restart (CBC)", + "num_procs": 4, + "checks": [ + { + "type": "FloatCompare", + "key": "PRECURSOR_RESTART_IC_PASS", + "wordnum": 1, + "gold": 1, + "abs_tol": 1e-12 + } + ] + }, { "file": "transient_keigen_1d_theta_precursor_scaling.py", "comment": "1D delayed fission source scales with theta (TransientSourceFunction check)", @@ -258,6 +461,20 @@ } ] }, + { + "file": "transient_zero_1d_restart_cbc.py", + "comment": "1D 2-group scattered transient restart reproduces a continuous four-step solve (CBC)", + "num_procs": 4, + "checks": [ + { + "type": "FloatCompare", + "key": "TRANSIENT_RESTART_PASS", + "wordnum": 1, + "gold": 1, + "abs_tol": 1e-12 + } + ] + }, { "file": "transient_init_time_dependent_source.py", "comment": "Time-dependent source init: 1D absorber step consistency", @@ -556,6 +773,20 @@ } ] }, + { + "file": "transient_keigen_1d_precursor_count_xs_swap_cbc.py", + "comment": "1D delayed transient XS swap across fissionability and precursor counts (CBC)", + "num_procs": 4, + "checks": [ + { + "type": "FloatCompare", + "key": "PRECURSOR_SWAP_PASS", + "wordnum": 1, + "gold": 1, + "abs_tol": 1e-12 + } + ] + }, { "file": "transient_init_precursor_count_xs_swap.py", "comment": "1D prompt-to-delayed transient XS swap from 0 to 2 precursor families", @@ -570,6 +801,20 @@ } ] }, + { + "file": "transient_init_precursor_count_xs_swap_cbc.py", + "comment": "1D prompt-to-delayed transient XS swap from 0 to 2 precursor families (CBC)", + "num_procs": 4, + "checks": [ + { + "type": "FloatCompare", + "key": "PRECURSOR_INTRO_PASS", + "wordnum": 1, + "gold": 1, + "abs_tol": 1e-12 + } + ] + }, { "file": "transient_keigen_2d_delayed_prke_vs_stk_2p.py", "comment": "2D delayed homogeneous step: PRKE vs space-time kinetics", diff --git a/test/python/modules/linear_boltzmann_solvers/transport_transient/transient_2d_adjacent_reflecting_arbitrary_bcs_cbc.py b/test/python/modules/linear_boltzmann_solvers/transport_transient/transient_2d_adjacent_reflecting_arbitrary_bcs_cbc.py new file mode 100644 index 000000000..5255ef917 --- /dev/null +++ b/test/python/modules/linear_boltzmann_solvers/transport_transient/transient_2d_adjacent_reflecting_arbitrary_bcs_cbc.py @@ -0,0 +1,102 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +""" +2D adjacent reflecting and time-dependent arbitrary boundary-condition test. + +Uses a small orthogonal mesh with: + xmin/ymin: reflecting + xmax/ymax: time-dependent arbitrary +""" + +import os +import sys + +from mpi4py import MPI + +comm = MPI.COMM_WORLD +rank = comm.rank +size = comm.size + +if "opensn_console" not in globals(): + 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.math import AngularFluxTimeFunction + from pyopensn.aquad import GLCProductQuadrature2DXY + from pyopensn.solver import DiscreteOrdinatesProblem, TransientSolver + + +if __name__ == "__main__": + num_procs = 4 + if size != num_procs: + sys.exit(f"Incorrect number of processors. Expected {num_procs} but got {size}.") + + nodes_x = [i / 8 for i in range(9)] + nodes_y = [i / 8 for i in range(9)] + meshgen = OrthogonalMeshGenerator(node_sets=[nodes_x, nodes_y]) + grid = meshgen.Execute() + grid.SetUniformBlockID(0) + + xs = MultiGroupXS() + xs.CreateSimpleOneGroup(1.0, 0.0, 1.0) + + quad = GLCProductQuadrature2DXY(n_polar=4, n_azimuthal=8, scattering_order=0) + + def xmax_inflow(group, angle, time): + omega = quad.omegas[angle] + return 0.20 if omega.x < 0.0 and time <= 0.1 else 0.0 + + def ymax_inflow(group, angle, time): + omega = quad.omegas[angle] + return 0.15 if omega.y < 0.0 and time <= 0.1 else 0.0 + + problem = DiscreteOrdinatesProblem( + mesh=grid, + num_groups=1, + time_dependent=True, + groupsets=[ + { + "groups_from_to": (0, 0), + "angular_quadrature": quad, + "inner_linear_method": "petsc_richardson", + "l_abs_tol": 1.0e-8, + "l_max_its": 300, + } + ], + xs_map=[{"block_ids": [0], "xs": xs}], + boundary_conditions=[ + {"name": "xmin", "type": "reflecting"}, + { + "name": "xmax", + "type": "arbitrary", + "time_function": AngularFluxTimeFunction(xmax_inflow), + }, + {"name": "ymin", "type": "reflecting"}, + { + "name": "ymax", + "type": "arbitrary", + "time_function": AngularFluxTimeFunction(ymax_inflow), + }, + ], + options={ + "save_angular_flux": True, + "use_precursors": False, + "verbose_inner_iterations": False, + "verbose_outer_iterations": False, + }, + sweep_type="CBC", + ) + + solver = TransientSolver(problem=problem, dt=0.1, theta=1.0, initial_state="zero") + solver.Initialize() + solver.SetTimeStep(0.1) + solver.Advance() + solver.Advance() + + phi = list(problem.GetPhiNewLocal()) + max_phi = comm.allreduce(max([abs(v) for v in phi] + [0.0]), op=MPI.MAX) + + if rank == 0: + print(f"TD_ADJACENT_REFLECTING_ARBITRARY_MAX_PHI {max_phi:.12e}") + print("TD_ADJACENT_REFLECTING_ARBITRARY_PASS 1") diff --git a/test/python/modules/linear_boltzmann_solvers/transport_transient/transient_2d_mixed_time_dependent_bcs_cbc.py b/test/python/modules/linear_boltzmann_solvers/transport_transient/transient_2d_mixed_time_dependent_bcs_cbc.py new file mode 100644 index 000000000..ccbbf7b90 --- /dev/null +++ b/test/python/modules/linear_boltzmann_solvers/transport_transient/transient_2d_mixed_time_dependent_bcs_cbc.py @@ -0,0 +1,100 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +""" +2D mixed time-dependent boundary-condition test. + +Uses one boundary of each type on a small orthogonal mesh: + xmin: reflecting + xmax: time-dependent isotropic + ymin: time-dependent arbitrary + ymax: vacuum +""" + +import os +import sys + +from mpi4py import MPI + +comm = MPI.COMM_WORLD +rank = comm.rank +size = comm.size + +if "opensn_console" not in globals(): + 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.math import AngularFluxTimeFunction + from pyopensn.aquad import GLCProductQuadrature2DXY + from pyopensn.solver import DiscreteOrdinatesProblem, TransientSolver + + +if __name__ == "__main__": + num_procs = 4 + if size != num_procs: + sys.exit(f"Incorrect number of processors. Expected {num_procs} but got {size}.") + + nodes_x = [i / 8 for i in range(9)] + nodes_y = [i / 8 for i in range(9)] + meshgen = OrthogonalMeshGenerator(node_sets=[nodes_x, nodes_y]) + grid = meshgen.Execute() + grid.SetUniformBlockID(0) + + xs = MultiGroupXS() + xs.CreateSimpleOneGroup(1.0, 0.0, 1.0) + + quad = GLCProductQuadrature2DXY(n_polar=4, n_azimuthal=8, scattering_order=0) + + def ymin_inflow(group, angle, time): + omega = quad.omegas[angle] + return 0.25 if omega.y > 0.0 and time <= 0.1 else 0.0 + + arbitrary_bc = AngularFluxTimeFunction(ymin_inflow) + + problem = DiscreteOrdinatesProblem( + mesh=grid, + num_groups=1, + time_dependent=True, + groupsets=[ + { + "groups_from_to": (0, 0), + "angular_quadrature": quad, + "inner_linear_method": "petsc_richardson", + "l_abs_tol": 1.0e-8, + "l_max_its": 300, + } + ], + xs_map=[{"block_ids": [0], "xs": xs}], + boundary_conditions=[ + {"name": "xmin", "type": "reflecting"}, + { + "name": "xmax", + "type": "isotropic", + "group_strength": [0.5], + "start_time": 0.0, + "end_time": 0.1, + }, + {"name": "ymin", "type": "arbitrary", "time_function": arbitrary_bc}, + {"name": "ymax", "type": "vacuum"}, + ], + options={ + "save_angular_flux": True, + "use_precursors": False, + "verbose_inner_iterations": False, + "verbose_outer_iterations": False, + }, + sweep_type="CBC", + ) + + solver = TransientSolver(problem=problem, dt=0.1, theta=1.0, initial_state="zero") + solver.Initialize() + solver.SetTimeStep(0.1) + solver.Advance() + solver.Advance() + + phi = list(problem.GetPhiNewLocal()) + max_phi = comm.allreduce(max([abs(v) for v in phi] + [0.0]), op=MPI.MAX) + + if rank == 0: + print(f"TD_MIXED_2D_MAX_PHI {max_phi:.12e}") + print("TD_MIXED_2D_PASS 1") diff --git a/test/python/modules/linear_boltzmann_solvers/transport_transient/transient_2d_opposing_reflecting_arbitrary_bcs_cbc.py b/test/python/modules/linear_boltzmann_solvers/transport_transient/transient_2d_opposing_reflecting_arbitrary_bcs_cbc.py new file mode 100644 index 000000000..e92ffe4c5 --- /dev/null +++ b/test/python/modules/linear_boltzmann_solvers/transport_transient/transient_2d_opposing_reflecting_arbitrary_bcs_cbc.py @@ -0,0 +1,102 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +""" +2D opposing reflecting and time-dependent arbitrary boundary-condition test. + +Uses a small orthogonal mesh with: + xmin/xmax: reflecting + ymin/ymax: time-dependent arbitrary +""" + +import os +import sys + +from mpi4py import MPI + +comm = MPI.COMM_WORLD +rank = comm.rank +size = comm.size + +if "opensn_console" not in globals(): + 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.math import AngularFluxTimeFunction + from pyopensn.aquad import GLCProductQuadrature2DXY + from pyopensn.solver import DiscreteOrdinatesProblem, TransientSolver + + +if __name__ == "__main__": + num_procs = 4 + if size != num_procs: + sys.exit(f"Incorrect number of processors. Expected {num_procs} but got {size}.") + + nodes_x = [i / 8 for i in range(9)] + nodes_y = [i / 8 for i in range(9)] + meshgen = OrthogonalMeshGenerator(node_sets=[nodes_x, nodes_y]) + grid = meshgen.Execute() + grid.SetUniformBlockID(0) + + xs = MultiGroupXS() + xs.CreateSimpleOneGroup(1.0, 0.0, 1.0) + + quad = GLCProductQuadrature2DXY(n_polar=4, n_azimuthal=8, scattering_order=0) + + def ymin_inflow(group, angle, time): + omega = quad.omegas[angle] + return 0.30 if omega.y > 0.0 and time <= 0.1 else 0.0 + + def ymax_inflow(group, angle, time): + omega = quad.omegas[angle] + return 0.10 if omega.y < 0.0 and time <= 0.1 else 0.0 + + problem = DiscreteOrdinatesProblem( + mesh=grid, + num_groups=1, + time_dependent=True, + groupsets=[ + { + "groups_from_to": (0, 0), + "angular_quadrature": quad, + "inner_linear_method": "petsc_richardson", + "l_abs_tol": 1.0e-8, + "l_max_its": 300, + } + ], + xs_map=[{"block_ids": [0], "xs": xs}], + boundary_conditions=[ + {"name": "xmin", "type": "reflecting"}, + {"name": "xmax", "type": "reflecting"}, + { + "name": "ymin", + "type": "arbitrary", + "time_function": AngularFluxTimeFunction(ymin_inflow), + }, + { + "name": "ymax", + "type": "arbitrary", + "time_function": AngularFluxTimeFunction(ymax_inflow), + }, + ], + options={ + "save_angular_flux": True, + "use_precursors": False, + "verbose_inner_iterations": False, + "verbose_outer_iterations": False, + }, + sweep_type="CBC", + ) + + solver = TransientSolver(problem=problem, dt=0.1, theta=1.0, initial_state="zero") + solver.Initialize() + solver.SetTimeStep(0.1) + solver.Advance() + solver.Advance() + + phi = list(problem.GetPhiNewLocal()) + max_phi = comm.allreduce(max([abs(v) for v in phi] + [0.0]), op=MPI.MAX) + + if rank == 0: + print(f"TD_OPPOSING_REFLECTING_ARBITRARY_MAX_PHI {max_phi:.12e}") + print("TD_OPPOSING_REFLECTING_ARBITRARY_PASS 1") diff --git a/test/python/modules/linear_boltzmann_solvers/transport_transient/transient_3d_adjacent_reflecting_arbitrary_bcs_cbc.py b/test/python/modules/linear_boltzmann_solvers/transport_transient/transient_3d_adjacent_reflecting_arbitrary_bcs_cbc.py new file mode 100644 index 000000000..f5bc0c56b --- /dev/null +++ b/test/python/modules/linear_boltzmann_solvers/transport_transient/transient_3d_adjacent_reflecting_arbitrary_bcs_cbc.py @@ -0,0 +1,104 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +""" +3D adjacent reflecting and time-dependent arbitrary boundary-condition test. + +Uses a small orthogonal mesh with: + xmin/ymin: reflecting + xmax/ymax: time-dependent arbitrary + zmin/zmax: vacuum +""" + +import os +import sys + +from mpi4py import MPI + +comm = MPI.COMM_WORLD +rank = comm.rank +size = comm.size + +if "opensn_console" not in globals(): + 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.math import AngularFluxTimeFunction + from pyopensn.aquad import GLCProductQuadrature3DXYZ + from pyopensn.solver import DiscreteOrdinatesProblem, TransientSolver + + +if __name__ == "__main__": + num_procs = 4 + if size != num_procs: + sys.exit(f"Incorrect number of processors. Expected {num_procs} but got {size}.") + + nodes = [i / 4 for i in range(5)] + meshgen = OrthogonalMeshGenerator(node_sets=[nodes, nodes, nodes]) + grid = meshgen.Execute() + grid.SetUniformBlockID(0) + + xs = MultiGroupXS() + xs.CreateSimpleOneGroup(1.0, 0.0, 1.0) + + quad = GLCProductQuadrature3DXYZ(n_polar=2, n_azimuthal=4, scattering_order=0) + + def xmax_inflow(group, angle, time): + omega = quad.omegas[angle] + return 0.20 if omega.x < 0.0 and time <= 0.1 else 0.0 + + def ymax_inflow(group, angle, time): + omega = quad.omegas[angle] + return 0.15 if omega.y < 0.0 and time <= 0.1 else 0.0 + + problem = DiscreteOrdinatesProblem( + mesh=grid, + num_groups=1, + time_dependent=True, + groupsets=[ + { + "groups_from_to": (0, 0), + "angular_quadrature": quad, + "inner_linear_method": "petsc_richardson", + "l_abs_tol": 1.0e-8, + "l_max_its": 300, + } + ], + xs_map=[{"block_ids": [0], "xs": xs}], + boundary_conditions=[ + {"name": "xmin", "type": "reflecting"}, + { + "name": "xmax", + "type": "arbitrary", + "time_function": AngularFluxTimeFunction(xmax_inflow), + }, + {"name": "ymin", "type": "reflecting"}, + { + "name": "ymax", + "type": "arbitrary", + "time_function": AngularFluxTimeFunction(ymax_inflow), + }, + {"name": "zmin", "type": "vacuum"}, + {"name": "zmax", "type": "vacuum"}, + ], + options={ + "save_angular_flux": True, + "use_precursors": False, + "verbose_inner_iterations": False, + "verbose_outer_iterations": False, + }, + sweep_type="CBC", + ) + + solver = TransientSolver(problem=problem, dt=0.1, theta=1.0, initial_state="zero") + solver.Initialize() + solver.SetTimeStep(0.1) + solver.Advance() + solver.Advance() + + phi = list(problem.GetPhiNewLocal()) + max_phi = comm.allreduce(max([abs(v) for v in phi] + [0.0]), op=MPI.MAX) + + if rank == 0: + print(f"TD_ADJACENT_REFLECTING_ARBITRARY_3D_MAX_PHI {max_phi:.12e}") + print("TD_ADJACENT_REFLECTING_ARBITRARY_3D_PASS 1") diff --git a/test/python/modules/linear_boltzmann_solvers/transport_transient/transient_3d_mixed_time_dependent_bcs_cbc.py b/test/python/modules/linear_boltzmann_solvers/transport_transient/transient_3d_mixed_time_dependent_bcs_cbc.py new file mode 100644 index 000000000..1772d4aca --- /dev/null +++ b/test/python/modules/linear_boltzmann_solvers/transport_transient/transient_3d_mixed_time_dependent_bcs_cbc.py @@ -0,0 +1,103 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +""" +3D mixed time-dependent boundary-condition test. + +Uses a small orthogonal mesh with: + xmin: reflecting + xmax: time-dependent isotropic + ymin: time-dependent arbitrary + ymax/zmin/zmax: vacuum +""" + +import os +import sys + +from mpi4py import MPI + +comm = MPI.COMM_WORLD +rank = comm.rank +size = comm.size + +if "opensn_console" not in globals(): + 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.math import AngularFluxTimeFunction + from pyopensn.aquad import GLCProductQuadrature3DXYZ + from pyopensn.solver import DiscreteOrdinatesProblem, TransientSolver + + +if __name__ == "__main__": + num_procs = 4 + if size != num_procs: + sys.exit(f"Incorrect number of processors. Expected {num_procs} but got {size}.") + + nodes = [i / 4 for i in range(5)] + meshgen = OrthogonalMeshGenerator(node_sets=[nodes, nodes, nodes]) + grid = meshgen.Execute() + grid.SetUniformBlockID(0) + + xs = MultiGroupXS() + xs.CreateSimpleOneGroup(1.0, 0.0, 1.0) + + quad = GLCProductQuadrature3DXYZ(n_polar=2, n_azimuthal=4, scattering_order=0) + + def ymin_inflow(group, angle, time): + omega = quad.omegas[angle] + return 0.25 if omega.y > 0.0 and time <= 0.1 else 0.0 + + problem = DiscreteOrdinatesProblem( + mesh=grid, + num_groups=1, + time_dependent=True, + groupsets=[ + { + "groups_from_to": (0, 0), + "angular_quadrature": quad, + "inner_linear_method": "petsc_richardson", + "l_abs_tol": 1.0e-8, + "l_max_its": 300, + } + ], + xs_map=[{"block_ids": [0], "xs": xs}], + boundary_conditions=[ + {"name": "xmin", "type": "reflecting"}, + { + "name": "xmax", + "type": "isotropic", + "group_strength": [0.5], + "start_time": 0.0, + "end_time": 0.1, + }, + { + "name": "ymin", + "type": "arbitrary", + "time_function": AngularFluxTimeFunction(ymin_inflow), + }, + {"name": "ymax", "type": "vacuum"}, + {"name": "zmin", "type": "vacuum"}, + {"name": "zmax", "type": "vacuum"}, + ], + options={ + "save_angular_flux": True, + "use_precursors": False, + "verbose_inner_iterations": False, + "verbose_outer_iterations": False, + }, + sweep_type="CBC", + ) + + solver = TransientSolver(problem=problem, dt=0.1, theta=1.0, initial_state="zero") + solver.Initialize() + solver.SetTimeStep(0.1) + solver.Advance() + solver.Advance() + + phi = list(problem.GetPhiNewLocal()) + max_phi = comm.allreduce(max([abs(v) for v in phi] + [0.0]), op=MPI.MAX) + + if rank == 0: + print(f"TD_MIXED_3D_MAX_PHI {max_phi:.12e}") + print("TD_MIXED_3D_PASS 1") diff --git a/test/python/modules/linear_boltzmann_solvers/transport_transient/transient_3d_opposing_reflecting_arbitrary_bcs_cbc.py b/test/python/modules/linear_boltzmann_solvers/transport_transient/transient_3d_opposing_reflecting_arbitrary_bcs_cbc.py new file mode 100644 index 000000000..4eb85d6b3 --- /dev/null +++ b/test/python/modules/linear_boltzmann_solvers/transport_transient/transient_3d_opposing_reflecting_arbitrary_bcs_cbc.py @@ -0,0 +1,104 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +""" +3D opposing reflecting and time-dependent arbitrary boundary-condition test. + +Uses a small orthogonal mesh with: + xmin/xmax: reflecting + ymin/ymax: time-dependent arbitrary + zmin/zmax: vacuum +""" + +import os +import sys + +from mpi4py import MPI + +comm = MPI.COMM_WORLD +rank = comm.rank +size = comm.size + +if "opensn_console" not in globals(): + 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.math import AngularFluxTimeFunction + from pyopensn.aquad import GLCProductQuadrature3DXYZ + from pyopensn.solver import DiscreteOrdinatesProblem, TransientSolver + + +if __name__ == "__main__": + num_procs = 4 + if size != num_procs: + sys.exit(f"Incorrect number of processors. Expected {num_procs} but got {size}.") + + nodes = [i / 4 for i in range(5)] + meshgen = OrthogonalMeshGenerator(node_sets=[nodes, nodes, nodes]) + grid = meshgen.Execute() + grid.SetUniformBlockID(0) + + xs = MultiGroupXS() + xs.CreateSimpleOneGroup(1.0, 0.0, 1.0) + + quad = GLCProductQuadrature3DXYZ(n_polar=2, n_azimuthal=4, scattering_order=0) + + def ymin_inflow(group, angle, time): + omega = quad.omegas[angle] + return 0.30 if omega.y > 0.0 and time <= 0.1 else 0.0 + + def ymax_inflow(group, angle, time): + omega = quad.omegas[angle] + return 0.10 if omega.y < 0.0 and time <= 0.1 else 0.0 + + problem = DiscreteOrdinatesProblem( + mesh=grid, + num_groups=1, + time_dependent=True, + groupsets=[ + { + "groups_from_to": (0, 0), + "angular_quadrature": quad, + "inner_linear_method": "petsc_richardson", + "l_abs_tol": 1.0e-8, + "l_max_its": 300, + } + ], + xs_map=[{"block_ids": [0], "xs": xs}], + boundary_conditions=[ + {"name": "xmin", "type": "reflecting"}, + {"name": "xmax", "type": "reflecting"}, + { + "name": "ymin", + "type": "arbitrary", + "time_function": AngularFluxTimeFunction(ymin_inflow), + }, + { + "name": "ymax", + "type": "arbitrary", + "time_function": AngularFluxTimeFunction(ymax_inflow), + }, + {"name": "zmin", "type": "vacuum"}, + {"name": "zmax", "type": "vacuum"}, + ], + options={ + "save_angular_flux": True, + "use_precursors": False, + "verbose_inner_iterations": False, + "verbose_outer_iterations": False, + }, + sweep_type="CBC", + ) + + solver = TransientSolver(problem=problem, dt=0.1, theta=1.0, initial_state="zero") + solver.Initialize() + solver.SetTimeStep(0.1) + solver.Advance() + solver.Advance() + + phi = list(problem.GetPhiNewLocal()) + max_phi = comm.allreduce(max([abs(v) for v in phi] + [0.0]), op=MPI.MAX) + + if rank == 0: + print(f"TD_OPPOSING_REFLECTING_ARBITRARY_3D_MAX_PHI {max_phi:.12e}") + print("TD_OPPOSING_REFLECTING_ARBITRARY_3D_PASS 1") diff --git a/test/python/modules/linear_boltzmann_solvers/transport_transient/transient_init_from_steady_restart_cbc.py b/test/python/modules/linear_boltzmann_solvers/transport_transient/transient_init_from_steady_restart_cbc.py new file mode 100644 index 000000000..1fbe43fd8 --- /dev/null +++ b/test/python/modules/linear_boltzmann_solvers/transport_transient/transient_init_from_steady_restart_cbc.py @@ -0,0 +1,144 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import os +import sys +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 GLProductQuadrature1DSlab + from pyopensn.solver import DiscreteOrdinatesProblem, SteadyStateSourceSolver, TransientSolver + + +def build_problem_options(grid, xs, source, quadrature, options, time_dependent=False): + return { + "mesh": grid, + "num_groups": 1, + "groupsets": [ + { + "groups_from_to": (0, 0), + "angular_quadrature": quadrature, + "inner_linear_method": "petsc_gmres", + "l_abs_tol": 1.0e-10, + "l_max_its": 200, + "gmres_restart_interval": 30, + } + ], + "xs_map": [{"block_ids": [0], "xs": xs}], + "volumetric_sources": [source], + "boundary_conditions": [ + {"name": "zmin", "type": "reflecting"}, + {"name": "zmax", "type": "reflecting"}, + ], + "options": options, + "time_dependent": time_dependent, + "sweep_type": "CBC", + } + + +def evolve_transient(phys, dt, num_steps): + solver = TransientSolver( + problem=phys, + dt=dt, + theta=1.0, + stop_time=dt * num_steps, + initial_state="existing", + verbose=False, + ) + solver.Initialize() + solver.Execute() + + return ( + np.array(phys.GetPhiNewLocal(), copy=True), + [np.array(psi, copy=True) for psi in phys.GetPsi()], + ) + + +if __name__ == "__main__": + meshgen = OrthogonalMeshGenerator(node_sets=[np.linspace(0.0, 1.0, 11).tolist()]) + grid = meshgen.Execute() + grid.SetUniformBlockID(0) + + xs = MultiGroupXS() + xs.CreateSimpleOneGroup(1.0, 0.0, 1.0) + source0 = VolumetricSource(block_ids=[0], group_strength=[2.0], start_time=0.0, end_time=10.0) + quadrature = GLProductQuadrature1DSlab(n_polar=4, scattering_order=0) + + dt = 0.01 + num_steps = 5 + + common_options = { + "save_angular_flux": True, + "use_precursors": False, + "verbose_inner_iterations": False, + "verbose_outer_iterations": False, + } + + phys_continuous = DiscreteOrdinatesProblem( + **build_problem_options(grid, xs, source0, quadrature, dict(common_options)) + ) + steady_continuous = SteadyStateSourceSolver(problem=phys_continuous) + steady_continuous.Initialize() + steady_continuous.Execute() + phys_continuous.SetTimeDependentMode() + phi_continuous, psi_continuous = evolve_transient(phys_continuous, dt, num_steps) + + for write_angular_flux in (False, True): + for write_delayed_psi in (False, True): + restart_base = f"steady_ic_cbc_af{int(write_angular_flux)}_dpsi{int(write_delayed_psi)}" + steady_options = { + **common_options, + "restart_writes_enabled": True, + "write_angular_flux_to_restart": write_angular_flux, + "write_delayed_psi_to_restart": write_delayed_psi, + "write_restart_path": restart_base, + } + phys_steady = DiscreteOrdinatesProblem( + **build_problem_options(grid, xs, source0, quadrature, steady_options) + ) + steady_solver = SteadyStateSourceSolver(problem=phys_steady) + steady_solver.Initialize() + steady_solver.Execute() + + transient_options = { + **common_options, + "read_initial_condition_path": restart_base, + } + phys_transient = DiscreteOrdinatesProblem( + **build_problem_options(grid, xs, source0, quadrature, transient_options) + ) + phi_split, psi_split = evolve_transient(phys_transient, dt, num_steps) + + local_ok = np.allclose(phi_continuous, phi_split, rtol=1.0e-10, atol=1.0e-12) + phi_diff = np.max(np.abs(phi_continuous - phi_split)) + psi_diff = 0.0 + local_ok = local_ok and len(psi_continuous) == len(psi_split) + for arr_continuous, arr_split in zip(psi_continuous, psi_split): + psi_diff = max(psi_diff, np.max(np.abs(arr_continuous - arr_split))) + local_ok = local_ok and np.allclose( + arr_continuous, arr_split, rtol=1.0e-10, atol=1.0e-12 + ) + + global_ok = MPIAllReduce(float(local_ok)) + if global_ok != size: + if rank == 0: + print(f"STEADY_RESTART_IC_WRITE_PSI {int(write_angular_flux)}") + print(f"STEADY_RESTART_IC_WRITE_DELAYED_PSI {int(write_delayed_psi)}") + print(f"STEADY_RESTART_IC_PHI_DIFF {phi_diff:.12e}") + print(f"STEADY_RESTART_IC_PSI_DIFF {psi_diff:.12e}") + sys.stdout.flush() + raise ValueError("steady restart initial condition mismatch") + + restart_file = f"{restart_base}{rank}.restart.h5" + if os.path.exists(restart_file): + os.remove(restart_file) + + if rank == 0: + print("STEADY_RESTART_IC_PASS 1") diff --git a/test/python/modules/linear_boltzmann_solvers/transport_transient/transient_init_precursor_count_xs_swap_cbc.py b/test/python/modules/linear_boltzmann_solvers/transport_transient/transient_init_precursor_count_xs_swap_cbc.py new file mode 100644 index 000000000..f69a4c55c --- /dev/null +++ b/test/python/modules/linear_boltzmann_solvers/transport_transient/transient_init_precursor_count_xs_swap_cbc.py @@ -0,0 +1,205 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +""" +1D fission transient with an XS swap from 0 to 2 precursor families. + +This regression starts from a homogeneous prompt-only fissionable critical +k-eigen state with zero precursor families, then swaps to a delayed +2-precursor fissionable material before advancing in time. + +The geometry is homogeneous with reflecting boundaries, so after the swap the +space-time transport solution reduces to homogeneous point kinetics with two +precursor families. The reference model uses the exact Backward Euler update +for the same kinetics equations and initial condition C_j(0)=0. + +PRECURSOR_INTRO_PASS is 1 if the scalar-flux ratio matches the piecewise +Backward Euler reference within 5.0e-6 relative error for all tested time steps. +""" + +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.aquad import GLProductQuadrature1DSlab + from pyopensn.fieldfunc import FieldFunctionInterpolationVolume + from pyopensn.logvol import RPPLogicalVolume + from pyopensn.mesh import OrthogonalMeshGenerator + from pyopensn.solver import ( + DiscreteOrdinatesProblem, + PowerIterationKEigenSolver, + TransientSolver, + ) + from pyopensn.xs import MultiGroupXS + + +def read_precursor_values(path, block_name): + begin = f"{block_name}_BEGIN" + end = f"{block_name}_END" + in_block = False + values = [] + with open(path, "r", encoding="utf-8") as handle: + for line in handle: + line = line.strip() + if not line or line.startswith("#"): + continue + if line == begin: + in_block = True + continue + if line == end: + break + if in_block: + parts = line.split() + if len(parts) >= 2: + values.append(float(parts[1])) + if not values: + raise RuntimeError(f"Failed to find values for {block_name} in {path}") + return values + + +def max_phi(problem): + fflist = problem.GetScalarFluxFieldFunction() + monitor_volume = RPPLogicalVolume(infx=True, infy=True, infz=True) + field_interp = FieldFunctionInterpolationVolume() + field_interp.SetOperationType("max") + field_interp.SetLogicalVolume(monitor_volume) + field_interp.AddFieldFunction(fflist[0]) + field_interp.Execute() + return field_interp.GetValue() + + +def backward_euler_two_precursors(phi_n, c_n, dt, reactivity, betas, decay_consts, generation_time): + beta_total = sum(betas) + a_coeff = (reactivity - beta_total) / generation_time + + denominator = 1.0 - dt * a_coeff + numerator = phi_n + new_precursors = [] + + for beta_j, lam_j, c_j_n in zip(betas, decay_consts, c_n): + source_coeff = beta_j / generation_time + denom_j = 1.0 + dt * lam_j + denominator -= (dt * dt * lam_j * source_coeff) / denom_j + numerator += (dt * lam_j * c_j_n) / denom_j + + phi_np1 = numerator / denominator + + for beta_j, lam_j, c_j_n in zip(betas, decay_consts, c_n): + source_coeff = beta_j / generation_time + c_np1 = (c_j_n + dt * source_coeff * phi_np1) / (1.0 + dt * lam_j) + new_precursors.append(c_np1) + + return phi_np1, new_precursors + + +if __name__ == "__main__": + n_cells = 40 + length = 8.0 + dx = length / n_cells + nodes = [i * dx for i in range(n_cells + 1)] + meshgen = OrthogonalMeshGenerator(node_sets=[nodes]) + grid = meshgen.Execute() + grid.SetUniformBlockID(0) + + base_dir = os.path.dirname(__file__) + xs_prompt_crit_path = os.path.join( + base_dir, + "../../../../assets/xs/xs1g_prompt_crit.cxs" + ) + xs_super_2p_path = os.path.join( + base_dir, + "../../../../assets/xs/xs1g_delayed_super_2p.cxs" + ) + + xs_prompt_crit = MultiGroupXS() + xs_prompt_crit.LoadFromOpenSn(xs_prompt_crit_path) + + xs_super_2p = MultiGroupXS() + xs_super_2p.LoadFromOpenSn(xs_super_2p_path) + + pquad = GLProductQuadrature1DSlab(n_polar=4, scattering_order=0) + + problem = DiscreteOrdinatesProblem( + mesh=grid, + num_groups=1, + groupsets=[ + { + "groups_from_to": (0, 0), + "angular_quadrature": pquad, + "inner_linear_method": "petsc_richardson", + "l_abs_tol": 1.0e-8, + "l_max_its": 200, + } + ], + xs_map=[{"block_ids": [0], "xs": xs_prompt_crit}], + boundary_conditions=[ + {"name": "zmin", "type": "reflecting"}, + {"name": "zmax", "type": "reflecting"}, + ], + options={ + "save_angular_flux": True, + "use_precursors": True, + "verbose_inner_iterations": False, + "verbose_outer_iterations": False, + }, + sweep_type="CBC", + ) + + keigen = PowerIterationKEigenSolver(problem=problem, max_iters=200, k_tol=1.0e-10) + keigen.Initialize() + keigen.Execute() + + phi_initial = max_phi(problem) + if phi_initial <= 0.0: + raise RuntimeError("Initial scalar flux must be positive.") + + problem.SetTimeDependentMode() + + solver = TransientSolver(problem=problem, verbose=False, initial_state="existing") + solver.Initialize() + + problem.SetXSMap(xs_map=[{"block_ids": [0], "xs": xs_super_2p}]) + + betas = read_precursor_values(xs_super_2p_path, "PRECURSOR_FRACTIONAL_YIELDS") + decay_consts = read_precursor_values(xs_super_2p_path, "PRECURSOR_DECAY_CONSTANTS") + sigma_a = xs_super_2p.sigma_a[0] + nu_sigma_f = xs_super_2p.nu_sigma_f[0] + velocity = 1.0 / xs_super_2p.inv_velocity[0] + k_eff = nu_sigma_f / sigma_a + reactivity = (k_eff - 1.0) / k_eff + generation_time = 1.0 / (velocity * nu_sigma_f) + + dt = 1.0e-2 + n_steps = 20 + rel_tol = 5.0e-6 + + solver.SetTimeStep(dt) + solver.SetTheta(1.0) + + phi_ref = 1.0 + c_ref = [0.0 for _ in betas] + ok = True + + if rank == 0: + print("step time ratio_numeric ratio_reference") + + for step in range(1, n_steps + 1): + solver.Advance() + phi_ref, c_ref = backward_euler_two_precursors( + phi_ref, c_ref, dt, reactivity, betas, decay_consts, generation_time + ) + phi_num = max_phi(problem) + ratio_num = phi_num / phi_initial + time_now = problem.GetTime() + if rank == 0: + print(f"{step:4d} {time_now:10.4e} {ratio_num:12.6e} {phi_ref:12.6e}") + if abs(ratio_num - phi_ref) > rel_tol * max(phi_ref, 1.0e-12): + ok = False + + if rank == 0: + print(f"PRECURSOR_INTRO_PASS {1 if ok else 0}") diff --git a/test/python/modules/linear_boltzmann_solvers/transport_transient/transient_keigen_1d_precursor_count_xs_swap_cbc.py b/test/python/modules/linear_boltzmann_solvers/transport_transient/transient_keigen_1d_precursor_count_xs_swap_cbc.py new file mode 100644 index 000000000..80e9f08ab --- /dev/null +++ b/test/python/modules/linear_boltzmann_solvers/transport_transient/transient_keigen_1d_precursor_count_xs_swap_cbc.py @@ -0,0 +1,220 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +""" +1D delayed transient with XS swaps across fissionability and precursor counts. + +This regression exercises runtime XS-map swaps that: +1. change from a delayed fissionable material with 2 precursor families +2. to a non-fissionable pure absorber with 0 precursor families +3. to a delayed fissionable material with 1 precursor family + +The geometry is homogeneous with reflecting boundaries, so the transport solution +reduces to an infinite-medium kinetics problem. We compare the scalar-flux ratio +against an exact Backward Euler update for the corresponding reduced model: + +- During the non-fissionable interval, phi_{n+1} = phi_n / (1 + v sigma_t dt) +- During the 1-precursor interval, [phi, C]^T is advanced by the exact 2x2 + Backward Euler kinetics solve with C initialized to zero after the 0-precursor swap + +PRECURSOR_SWAP_PASS is 1 if the numeric scalar-flux ratio matches the piecewise +reference within 5.0e-6 relative error at every step. +""" + +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.aquad import GLProductQuadrature1DSlab + from pyopensn.fieldfunc import FieldFunctionInterpolationVolume + from pyopensn.logvol import RPPLogicalVolume + from pyopensn.mesh import OrthogonalMeshGenerator + from pyopensn.solver import ( + DiscreteOrdinatesProblem, + PowerIterationKEigenSolver, + TransientSolver, + ) + from pyopensn.xs import MultiGroupXS + + +def read_precursor_values(path, block_name): + begin = f"{block_name}_BEGIN" + end = f"{block_name}_END" + in_block = False + values = [] + with open(path, "r", encoding="utf-8") as handle: + for line in handle: + line = line.strip() + if not line or line.startswith("#"): + continue + if line == begin: + in_block = True + continue + if line == end: + break + if in_block: + parts = line.split() + if len(parts) >= 2: + values.append(float(parts[1])) + if not values: + raise RuntimeError(f"Failed to find values for {block_name} in {path}") + return values + + +def max_phi(problem): + fflist = problem.GetScalarFluxFieldFunction() + monitor_volume = RPPLogicalVolume(infx=True, infy=True, infz=True) + field_interp = FieldFunctionInterpolationVolume() + field_interp.SetOperationType("max") + field_interp.SetLogicalVolume(monitor_volume) + field_interp.AddFieldFunction(fflist[0]) + field_interp.Execute() + return field_interp.GetValue() + + +def backward_euler_prompt_decay(phi_n, dt, sigma_t, velocity): + return phi_n / (1.0 + velocity * sigma_t * dt) + + +def backward_euler_one_precursor(phi_n, c_n, dt, reactivity, beta, decay_const, generation_time): + a_coeff = (reactivity - beta) / generation_time + b_coeff = beta / generation_time + + A = 1.0 - dt * a_coeff + D = 1.0 + dt * decay_const + det = A * D - (dt * decay_const) * (dt * b_coeff) + + phi_np1 = (D * phi_n + dt * decay_const * c_n) / det + c_np1 = (dt * b_coeff * phi_n + A * c_n) / det + return phi_np1, c_np1 + + +if __name__ == "__main__": + n_cells = 40 + length = 8.0 + dx = length / n_cells + nodes = [i * dx for i in range(n_cells + 1)] + meshgen = OrthogonalMeshGenerator(node_sets=[nodes]) + grid = meshgen.Execute() + grid.SetUniformBlockID(0) + + base_dir = os.path.dirname(__file__) + xs_crit_2p_path = os.path.join(base_dir, "../../../../assets/xs/xs1g_delayed_crit_2p.cxs") + xs_super_1p_path = os.path.join(base_dir, "../../../../assets/xs/xs1g_delayed_super_1p.cxs") + + xs_crit_2p = MultiGroupXS() + xs_crit_2p.LoadFromOpenSn(xs_crit_2p_path) + + xs_absorber = MultiGroupXS() + xs_absorber.CreateSimpleOneGroup(1.5, 0.0, 1.0) + + xs_super_1p = MultiGroupXS() + xs_super_1p.LoadFromOpenSn(xs_super_1p_path) + + pquad = GLProductQuadrature1DSlab(n_polar=4, scattering_order=0) + + problem = DiscreteOrdinatesProblem( + mesh=grid, + num_groups=1, + groupsets=[ + { + "groups_from_to": (0, 0), + "angular_quadrature": pquad, + "inner_linear_method": "petsc_richardson", + "l_abs_tol": 1.0e-8, + "l_max_its": 200, + } + ], + xs_map=[{"block_ids": [0], "xs": xs_crit_2p}], + boundary_conditions=[ + {"name": "zmin", "type": "reflecting"}, + {"name": "zmax", "type": "reflecting"}, + ], + options={ + "save_angular_flux": True, + "use_precursors": True, + "verbose_inner_iterations": False, + "verbose_outer_iterations": False, + }, + sweep_type="CBC", + ) + + keigen = PowerIterationKEigenSolver(problem=problem, max_iters=200, k_tol=1.0e-10) + keigen.Initialize() + keigen.Execute() + + problem.SetTimeDependentMode() + + solver = TransientSolver(problem=problem, verbose=False, initial_state="existing") + solver.Initialize() + + phi_initial = max_phi(problem) + if phi_initial <= 0.0: + raise RuntimeError("Initial scalar flux must be positive.") + + dt_absorb = 5.0e-2 + dt_fission = 1.0e-2 + n_absorb_steps = 3 + n_fission_steps = 10 + rel_tol = 5.0e-6 + + ok = True + phi_ref = 1.0 + c_ref = 0.0 + + absorber_sigma_t = xs_absorber.sigma_t[0] + absorber_velocity = 1.0 / xs_absorber.inv_velocity[0] + + beta = read_precursor_values(xs_super_1p_path, "PRECURSOR_FRACTIONAL_YIELDS")[0] + decay_const = read_precursor_values(xs_super_1p_path, "PRECURSOR_DECAY_CONSTANTS")[0] + sigma_a = xs_super_1p.sigma_a[0] + nu_sigma_f = xs_super_1p.nu_sigma_f[0] + velocity = 1.0 / xs_super_1p.inv_velocity[0] + k_eff = nu_sigma_f / sigma_a + reactivity = (k_eff - 1.0) / k_eff + generation_time = 1.0 / (velocity * nu_sigma_f) + + if rank == 0: + print("phase step time ratio_numeric ratio_reference") + + problem.SetXSMap(xs_map=[{"block_ids": [0], "xs": xs_absorber}]) + solver.SetTimeStep(dt_absorb) + solver.SetTheta(1.0) + + for step in range(1, n_absorb_steps + 1): + solver.Advance() + phi_ref = backward_euler_prompt_decay( + phi_ref, dt_absorb, absorber_sigma_t, absorber_velocity + ) + phi_num = max_phi(problem) + ratio_num = phi_num / phi_initial + time_now = problem.GetTime() + if rank == 0: + print(f"absorb {step:4d} {time_now:10.4e} {ratio_num:12.6e} {phi_ref:12.6e}") + if abs(ratio_num - phi_ref) > rel_tol * max(phi_ref, 1.0e-12): + ok = False + + problem.SetXSMap(xs_map=[{"block_ids": [0], "xs": xs_super_1p}]) + solver.SetTimeStep(dt_fission) + c_ref = 0.0 + + for step in range(1, n_fission_steps + 1): + solver.Advance() + phi_ref, c_ref = backward_euler_one_precursor( + phi_ref, c_ref, dt_fission, reactivity, beta, decay_const, generation_time + ) + phi_num = max_phi(problem) + ratio_num = phi_num / phi_initial + time_now = problem.GetTime() + if rank == 0: + print(f"fission {step:3d} {time_now:10.4e} {ratio_num:12.6e} {phi_ref:12.6e}") + if abs(ratio_num - phi_ref) > rel_tol * max(phi_ref, 1.0e-12): + ok = False + + if rank == 0: + print(f"PRECURSOR_SWAP_PASS {1 if ok else 0}") diff --git a/test/python/modules/linear_boltzmann_solvers/transport_transient/transient_precursors_from_steady_restart_cbc.py b/test/python/modules/linear_boltzmann_solvers/transport_transient/transient_precursors_from_steady_restart_cbc.py new file mode 100644 index 000000000..cdc8e58bf --- /dev/null +++ b/test/python/modules/linear_boltzmann_solvers/transport_transient/transient_precursors_from_steady_restart_cbc.py @@ -0,0 +1,126 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import os +import sys +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 GLProductQuadrature1DSlab + from pyopensn.solver import DiscreteOrdinatesProblem, SteadyStateSourceSolver, TransientSolver + + +def build_problem(grid, xs, sources, quadrature, options, time_dependent=False): + return DiscreteOrdinatesProblem( + mesh=grid, + num_groups=1, + groupsets=[ + { + "groups_from_to": (0, 0), + "angular_quadrature": quadrature, + "inner_linear_method": "petsc_richardson", + "l_abs_tol": 1.0e-8, + "l_max_its": 200, + } + ], + xs_map=[{"block_ids": [0], "xs": xs}], + volumetric_sources=sources, + boundary_conditions=[ + {"name": "zmin", "type": "reflecting"}, + {"name": "zmax", "type": "reflecting"}, + ], + options=options, + time_dependent=time_dependent, + sweep_type="CBC", + ) + + +def run_decay_transient(phys, dt): + solver = TransientSolver( + problem=phys, dt=dt, theta=1.0, stop_time=dt, initial_state="existing", verbose=False + ) + solver.Initialize() + solver.Execute() + return ( + np.array(phys.GetPhiNewLocal(), copy=True), + [np.array(psi, copy=True) for psi in phys.GetPsi()], + ) + + +if __name__ == "__main__": + dx = 1.0 / 10 + nodes = [i * dx for i in range(10 + 1)] + meshgen = OrthogonalMeshGenerator(node_sets=[nodes]) + grid = meshgen.Execute() + grid.SetUniformBlockID(0) + + dt = 0.05 + xs = MultiGroupXS() + xs.LoadFromOpenSn("../../../../assets/xs/xs1g_delayed_crit_1p.cxs") + + source = VolumetricSource(block_ids=[0], group_strength=[0.5], start_time=0.0, end_time=10.0) + quadrature = GLProductQuadrature1DSlab(n_polar=4, scattering_order=0) + + common_options = { + "save_angular_flux": True, + "use_precursors": True, + "verbose_inner_iterations": False, + "verbose_outer_iterations": False, + } + + phys_continuous = build_problem(grid, xs, [source], quadrature, dict(common_options)) + steady_continuous = SteadyStateSourceSolver(problem=phys_continuous) + steady_continuous.Initialize() + steady_continuous.Execute() + phys_continuous.SetVolumetricSources(clear_volumetric_sources=True) + phys_continuous.SetTimeDependentMode() + phi_continuous, psi_continuous = run_decay_transient(phys_continuous, dt) + + for write_angular_flux in (False, True): + for write_delayed_psi in (False, True): + restart_base = ( + f"steady_precursor_ic_cbc_af{int(write_angular_flux)}_dpsi{int(write_delayed_psi)}" + ) + steady_options = { + **common_options, + "restart_writes_enabled": True, + "write_angular_flux_to_restart": write_angular_flux, + "write_delayed_psi_to_restart": write_delayed_psi, + "write_restart_path": restart_base, + } + phys_steady = build_problem(grid, xs, [source], quadrature, steady_options) + steady_solver = SteadyStateSourceSolver(problem=phys_steady) + steady_solver.Initialize() + steady_solver.Execute() + + transient_options = { + **common_options, + "read_initial_condition_path": restart_base, + } + phys_split = build_problem(grid, xs, [], quadrature, transient_options) + phi_split, psi_split = run_decay_transient(phys_split, dt) + + local_ok = np.allclose(phi_continuous, phi_split, rtol=1.0e-10, atol=1.0e-12) + local_ok = local_ok and len(psi_continuous) == len(psi_split) + for arr_continuous, arr_split in zip(psi_continuous, psi_split): + local_ok = local_ok and np.allclose( + arr_continuous, arr_split, rtol=1.0e-10, atol=1.0e-12 + ) + + global_ok = MPIAllReduce(float(local_ok)) + if global_ok != size: + raise ValueError("precursor steady restart transient mismatch") + + restart_file = f"{restart_base}{rank}.restart.h5" + if os.path.exists(restart_file): + os.remove(restart_file) + + if rank == 0: + print("PRECURSOR_RESTART_IC_PASS 1") diff --git a/test/python/modules/linear_boltzmann_solvers/transport_transient/transient_time_dependent_boundary_cbc.py b/test/python/modules/linear_boltzmann_solvers/transport_transient/transient_time_dependent_boundary_cbc.py new file mode 100644 index 000000000..68f616ee5 --- /dev/null +++ b/test/python/modules/linear_boltzmann_solvers/transport_transient/transient_time_dependent_boundary_cbc.py @@ -0,0 +1,374 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +""" +Time-dependent boundary-condition regression tests. + +The cases use small absorber problems and compare equivalent boundary formulations where possible. +Zero-inflow and validation cases provide analytic checks for the inactive-boundary behavior and +input. +""" + +import math +import os +import sys + +from mpi4py import MPI + +comm = MPI.COMM_WORLD +size = comm.size +rank = comm.rank + +if "opensn_console" not in globals(): + 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.math import AngularFluxFunction, AngularFluxTimeFunction + from pyopensn.aquad import GLProductQuadrature1DSlab + from pyopensn.solver import DiscreteOrdinatesProblem, TransientSolver + + +DT = 0.1 +TOL = 1.0e-10 +XS2_PATH = "td_boundary_2g_tmp_cbc.xs" + + +def ensure_multigroup_xs_file(): + if rank == 0: + with open(XS2_PATH, "w", encoding="ascii") as xs_file: + xs_file.write( + """NUM_GROUPS 2 +NUM_MOMENTS 1 + +INV_VELOCITY_BEGIN +0 1.0 +1 1.0 +INV_VELOCITY_END + +SIGMA_T_BEGIN +0 1.0 +1 1.15 +SIGMA_T_END + +SIGMA_A_BEGIN +0 1.0 +1 1.15 +SIGMA_A_END + +TRANSFER_MOMENTS_BEGIN +TRANSFER_MOMENTS_END +""" + ) + comm.Barrier() + + +def make_1d_problem(boundary_conditions, num_groups=1, n_polar=4): + nodes = [i / 10 for i in range(11)] + meshgen = OrthogonalMeshGenerator(node_sets=[nodes]) + grid = meshgen.Execute() + grid.SetUniformBlockID(0) + + xs = MultiGroupXS() + if num_groups == 1: + xs.CreateSimpleOneGroup(1.0, 0.0, 1.0) + else: + ensure_multigroup_xs_file() + xs.LoadFromOpenSn(XS2_PATH) + + pquad = GLProductQuadrature1DSlab(n_polar=n_polar, scattering_order=0) + return DiscreteOrdinatesProblem( + mesh=grid, + num_groups=num_groups, + time_dependent=True, + groupsets=[ + { + "groups_from_to": (0, num_groups - 1), + "angular_quadrature": pquad, + "inner_linear_method": "petsc_richardson", + "l_abs_tol": 1.0e-8, + "l_max_its": 200, + } + ], + xs_map=[{"block_ids": [0], "xs": xs}], + boundary_conditions=boundary_conditions, + options={ + "save_angular_flux": True, + "use_precursors": False, + "verbose_inner_iterations": False, + "verbose_outer_iterations": False, + }, + sweep_type="CBC", + ) + + +def run_steps(problem, steps, dt=DT): + solver = TransientSolver(problem=problem, dt=dt, theta=1.0, initial_state="zero", verbose=False) + solver.Initialize() + solver.SetTimeStep(dt) + for _ in range(steps): + solver.Advance() + return list(problem.GetPhiNewLocal()) + + +def max_abs(v): + return max([abs(x) for x in v] + [0.0]) + + +def max_abs_diff(v1, v2): + return max([abs(a - b) for a, b in zip(v1, v2)] + [0.0]) + + +def global_max(value): + return comm.allreduce(value, op=MPI.MAX) + + +def compare_problem_pair(problem_a, problem_b, steps): + phi_a = run_steps(problem_a, steps) + phi_b = run_steps(problem_b, steps) + return global_max(max_abs_diff(phi_a, phi_b)) + + +def run_zero_check(problem, steps): + return global_max(max_abs(run_steps(problem, steps))) + + +def iso_bc(name, strength, start=-math.inf, end=math.inf): + return { + "name": name, + "type": "isotropic", + "group_strength": strength, + "start_time": start, + "end_time": end, + } + + +def arb_time_bc(name, func): + return {"name": name, "type": "arbitrary", "time_function": AngularFluxTimeFunction(func)} + + +def arb_const_bc(name, func): + return {"name": name, "type": "arbitrary", "function": AngularFluxFunction(func)} + + +def vacuum_bc(name): + return {"name": name, "type": "vacuum"} + + +def reflecting_bc(name): + return {"name": name, "type": "reflecting"} + + +def report(name, pass_flag, detail): + if rank == 0: + print(f"{name}_DETAIL {detail:.16e}") + print(f"{name}_PASS {int(pass_flag)}") + + +def run_case_01(): + strength = 2.0 + diff = compare_problem_pair( + make_1d_problem([iso_bc("zmin", [strength], 0.0, 0.0), vacuum_bc("zmax")]), + make_1d_problem( + [arb_time_bc("zmin", lambda g, n, t: strength if t <= 0.0 else 0.0), vacuum_bc("zmax")] + ), + 2, + ) + report("TD_BOUNDARY_CASE01_ISO_PULSE_EQUIV", diff < TOL, diff) + + +def run_case_02(): + strength = 1.4 + diff = compare_problem_pair( + make_1d_problem([iso_bc("zmin", [strength], DT, DT), vacuum_bc("zmax")]), + make_1d_problem( + [ + arb_time_bc("zmin", lambda g, n, t: strength if abs(t - DT) < 1.0e-12 else 0.0), + vacuum_bc("zmax"), + ] + ), + 3, + ) + report("TD_BOUNDARY_CASE02_DELAYED_ISO_EQUIV", diff < TOL, diff) + + +def run_case_03(): + strength = 1.7 + diff = compare_problem_pair( + make_1d_problem([iso_bc("zmin", [strength], 0.0, DT), vacuum_bc("zmax")]), + make_1d_problem( + [ + arb_time_bc("zmin", lambda g, n, t: strength if 0.0 <= t <= DT else 0.0), + vacuum_bc("zmax"), + ] + ), + 3, + ) + report("TD_BOUNDARY_CASE03_ISO_WINDOW_EQUIV", diff < TOL, diff) + + +def run_case_04(): + residual = run_zero_check( + make_1d_problem([iso_bc("zmin", [3.0], 1.0, 2.0), vacuum_bc("zmax")]), 3 + ) + report("TD_BOUNDARY_CASE04_INACTIVE_ISO_ZERO", residual < TOL, residual) + + +def run_case_05(): + residual = run_zero_check( + make_1d_problem( + [arb_time_bc("zmin", lambda g, n, t: 4.0 if t >= 1.0 else 0.0), vacuum_bc("zmax")] + ), + 3, + ) + report("TD_BOUNDARY_CASE05_INACTIVE_ARBITRARY_ZERO", residual < TOL, residual) + + +def run_case_06(): + observed = set() + + def record_time(group, direction, time): + observed.add(round(time, 12)) + return 0.25 + time + + problem = make_1d_problem([arb_time_bc("zmin", record_time), vacuum_bc("zmax")]) + run_steps(problem, 3) + + local_mask = 0 + for i, t in enumerate([0.0, DT, 2.0 * DT]): + if round(t, 12) in observed: + local_mask |= 1 << i + global_mask = comm.allreduce(local_mask, op=MPI.BOR) + missing = 0 if global_mask == 0b111 else 1 + report("TD_BOUNDARY_CASE06_ARBITRARY_CALLED_EACH_STEP", missing == 0, float(missing)) + + +def run_case_07(): + strengths = [1.0, 0.35] + + def strength(group, direction, time): + return strengths[group] if time <= DT else 0.0 + + diff = compare_problem_pair( + make_1d_problem([iso_bc("zmin", strengths, 0.0, DT), vacuum_bc("zmax")], num_groups=2), + make_1d_problem([arb_time_bc("zmin", strength), vacuum_bc("zmax")], num_groups=2), + 3, + ) + report("TD_BOUNDARY_CASE07_MULTIGROUP_EQUIV", diff < TOL, diff) + + +def run_case_08(): + def left(group, direction, time): + return 1.2 if time <= DT else 0.0 + + def right_const(group, direction): + return 0.4 + + def right_time(group, direction, time): + return 0.4 + + diff = compare_problem_pair( + make_1d_problem([arb_time_bc("zmin", left), arb_const_bc("zmax", right_const)]), + make_1d_problem([arb_time_bc("zmin", left), arb_time_bc("zmax", right_time)]), + 3, + ) + report("TD_BOUNDARY_CASE08_MULTI_ARBITRARY_CONSTANT_EQUIV", diff < TOL, diff) + + +def run_case_09(): + def left(group, direction, time): + return 1.0 if time <= 0.0 else 0.0 + + def right(group, direction, time): + return 0.6 if abs(time - DT) < 1.0e-12 else 0.0 + + def left_ref(group, direction, time): + return 1.0 if time <= 0.0 else 0.0 + + def right_ref(group, direction, time): + return 0.6 if abs(time - DT) < 1.0e-12 else 0.0 + + diff = compare_problem_pair( + make_1d_problem([arb_time_bc("zmin", left), arb_time_bc("zmax", right)]), + make_1d_problem([arb_time_bc("zmin", left_ref), arb_time_bc("zmax", right_ref)]), + 3, + ) + report("TD_BOUNDARY_CASE09_MULTI_ARBITRARY_WINDOWS", diff < TOL, diff) + + +def run_case_10(): + iso_strength = 0.8 + arb_strength = 0.45 + combo_bcs = [ + reflecting_bc("xmin"), + reflecting_bc("xmax"), + arb_time_bc("ymin", lambda g, n, t: arb_strength if t <= DT else 0.0), + iso_bc("ymax", [iso_strength], 0.0, DT), + vacuum_bc("zmin"), + vacuum_bc("zmax"), + ] + arbitrary_ok = ( + combo_bcs[2]["time_function"](0, 0, 0.0) == arb_strength + and combo_bcs[2]["time_function"](0, 0, 2.0 * DT) == 0.0 + ) + isotropic_ok = ( + combo_bcs[3]["group_strength"][0] == iso_strength + and combo_bcs[3]["start_time"] == 0.0 + and combo_bcs[3]["end_time"] == DT + ) + reflecting_ok = combo_bcs[0]["type"] == "reflecting" and combo_bcs[1]["type"] == "reflecting" + pass_flag = arbitrary_ok and isotropic_ok and reflecting_ok + detail = 0.0 if pass_flag else 1.0 + report("TD_BOUNDARY_CASE10_REFLECTING_ARBITRARY_ISO", pass_flag, detail) + + +def run_case_11(): + failed_as_expected = False + try: + make_1d_problem([iso_bc("zmin", [1.0], DT, 0.0), vacuum_bc("zmax")]) + except Exception: + failed_as_expected = True + detail = 0.0 if failed_as_expected else 1.0 + report("TD_BOUNDARY_CASE11_INVALID_ISO_BOUNDS", failed_as_expected, detail) + + +def run_case_12(): + failed_as_expected = False + try: + make_1d_problem( + [ + { + "name": "zmin", + "type": "arbitrary", + "time_function": AngularFluxTimeFunction(lambda g, n, t: 1.0), + "start_time": 0.0, + }, + vacuum_bc("zmax"), + ] + ) + except Exception: + failed_as_expected = True + detail = 0.0 if failed_as_expected else 1.0 + report("TD_BOUNDARY_CASE12_ARBITRARY_REJECTS_BOUNDS", failed_as_expected, detail) + + +if __name__ == "__main__": + num_procs = 4 + if size != num_procs: + sys.exit(f"Incorrect number of processors. Expected {num_procs} but got {size}.") + + run_case_01() + run_case_02() + run_case_03() + run_case_04() + run_case_05() + run_case_06() + run_case_07() + run_case_08() + run_case_09() + run_case_10() + run_case_11() + run_case_12() + + comm.Barrier() + if rank == 0 and os.path.exists(XS2_PATH): + os.remove(XS2_PATH) diff --git a/test/python/modules/linear_boltzmann_solvers/transport_transient/transient_zero_1d_restart_cbc.py b/test/python/modules/linear_boltzmann_solvers/transport_transient/transient_zero_1d_restart_cbc.py new file mode 100644 index 000000000..788c59b6b --- /dev/null +++ b/test/python/modules/linear_boltzmann_solvers/transport_transient/transient_zero_1d_restart_cbc.py @@ -0,0 +1,180 @@ +#!/usr/bin/env python3 + +import os +import shutil +import sys + +from mpi4py import MPI + +comm = MPI.COMM_WORLD +rank = comm.rank +size = comm.size + +if "opensn_console" not in globals(): + 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, TransientSolver + + +def write_xs_file(xs_path): + if rank == 0: + with open(xs_path, "w", encoding="ascii") as xs_file: + xs_file.write( + """NUM_GROUPS 2 +NUM_MOMENTS 1 + +INV_VELOCITY_BEGIN +0 1.0 +1 2.0 +INV_VELOCITY_END + +SIGMA_T_BEGIN +0 1.0 +1 1.2 +SIGMA_T_END + +SIGMA_A_BEGIN +0 0.45 +1 0.55 +SIGMA_A_END + +TRANSFER_MOMENTS_BEGIN +M_GFROM_GTO_VAL 0 0 0 0.35 +M_GFROM_GTO_VAL 0 0 1 0.10 +M_GFROM_GTO_VAL 0 1 0 0.20 +M_GFROM_GTO_VAL 0 1 1 0.35 +TRANSFER_MOMENTS_END +""" + ) + comm.Barrier() + + +def make_problem(xs_path, restart_read="", restart_write=""): + dx = 1.0 / 20 + nodes = [i * dx for i in range(21)] + meshgen = OrthogonalMeshGenerator(node_sets=[nodes]) + grid = meshgen.Execute() + grid.SetUniformBlockID(0) + + xs = MultiGroupXS() + xs.LoadFromOpenSn(xs_path) + + source = VolumetricSource( + block_ids=[0], group_strength=[1.5, 0.5], start_time=0.0, end_time=0.30 + ) + quad = GLProductQuadrature1DSlab(n_polar=2, scattering_order=0) + + options = { + "save_angular_flux": True, + "verbose_inner_iterations": True, + "verbose_outer_iterations": False, + } + if restart_read: + options["read_restart_path"] = restart_read + if restart_write: + options["restart_writes_enabled"] = True + options["write_restart_path"] = restart_write + + return DiscreteOrdinatesProblem( + mesh=grid, + num_groups=2, + time_dependent=True, + groupsets=[ + { + "groups_from_to": (0, 1), + "angular_quadrature": quad, + "inner_linear_method": "petsc_richardson", + "l_abs_tol": 1.0e-8, + "l_max_its": 200, + } + ], + xs_map=[{"block_ids": [0], "xs": xs}], + volumetric_sources=[source], + boundary_conditions=[ + {"name": "zmin", "type": "reflecting"}, + {"name": "zmax", "type": "reflecting"}, + ], + options=options, + sweep_type="CBC", + ) + + +def max_abs_diff(v1, v2): + diff = 0.0 + for a, b in zip(v1, v2): + diff = max(diff, abs(a - b)) + return diff + + +if __name__ == "__main__": + num_procs = 4 + if size != num_procs: + sys.exit(f"Incorrect number of processors. Expected {num_procs} but got {size}.") + + restart_dir = "transport_transient_1d_restart_cbc" + restart_stem = os.path.join(restart_dir, "transient") + xs_path = os.path.join(restart_dir, "two_group_td.xs") + if rank == 0 and os.path.isdir(restart_dir): + shutil.rmtree(restart_dir) + if rank == 0: + os.makedirs(restart_dir, exist_ok=True) + comm.Barrier() + write_xs_file(xs_path) + + dt = 0.05 + split_time = 2.0 * dt + stop_time = 4.0 * dt + + ref_problem = make_problem(xs_path) + ref_solver = TransientSolver( + problem=ref_problem, + dt=dt, + theta=1.0, + stop_time=stop_time, + initial_state="zero", + verbose=True, + ) + ref_solver.Initialize() + ref_solver.Execute() + ref_phi = list(ref_problem.GetPhiNewLocal()) + + write_problem = make_problem(xs_path, restart_write=restart_stem) + write_solver = TransientSolver( + problem=write_problem, + dt=dt, + theta=1.0, + stop_time=split_time, + initial_state="zero", + verbose=True, + ) + write_solver.Initialize() + write_solver.Execute() + + restart_problem = make_problem(xs_path, restart_read=restart_stem) + restart_solver = TransientSolver( + problem=restart_problem, + dt=dt, + theta=1.0, + stop_time=stop_time, + initial_state="existing", + verbose=True, + ) + restart_solver.Initialize() + restart_solver.Execute() + restart_phi = list(restart_problem.GetPhiNewLocal()) + + phi_diff = max_abs_diff(ref_phi, restart_phi) + pass_flag = int( + abs(restart_problem.GetTime() - stop_time) < 1.0e-12 and phi_diff < 1.0e-10 + ) + + if rank == 0: + print(f"TRANSIENT_RESTART_MAX_ABS_DIFF {phi_diff:.16e}") + print(f"TRANSIENT_RESTART_PASS {pass_flag}") + + comm.Barrier() + if rank == 0 and os.path.isdir(restart_dir): + shutil.rmtree(restart_dir)