From de63529ddcb7586ec9e5a84eb007f085593aac51 Mon Sep 17 00:00:00 2001 From: james <81617086+je-cook@users.noreply.github.com> Date: Mon, 3 Nov 2025 09:42:09 +0000 Subject: [PATCH 01/17] =?UTF-8?q?=F0=9F=8E=A8=20First=20commit=20of=20dagm?= =?UTF-8?q?c=20solver?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- bluemira/codes/openmc/__init__.py | 2 +- bluemira/codes/openmc/output.py | 2 +- bluemira/codes/openmc/solver.py | 488 +++++++++++++++--- bluemira/codes/openmc/tallying.py | 21 + bluemira/codes/wrapper.py | 2 +- eudemo/eudemo/model_managers.py | 4 +- eudemo/eudemo/neutronics/run.py | 14 +- eudemo/eudemo/reactor.py | 12 + .../run_cad_neutronics.ex.py | 299 ++++++++++- 9 files changed, 736 insertions(+), 108 deletions(-) diff --git a/bluemira/codes/openmc/__init__.py b/bluemira/codes/openmc/__init__.py index 2e89083aed..13fa1d4b40 100644 --- a/bluemira/codes/openmc/__init__.py +++ b/bluemira/codes/openmc/__init__.py @@ -5,6 +5,6 @@ # SPDX-License-Identifier: LGPL-2.1-or-later """OpenMC interface""" -from bluemira.codes.openmc.solver import OpenMCNeutronicsSolver as Solver +from bluemira.codes.openmc.solver import OpenMCCSGNeutronicsSolver as Solver __all__ = ["Solver"] diff --git a/bluemira/codes/openmc/output.py b/bluemira/codes/openmc/output.py index 3722330e11..7a20042897 100644 --- a/bluemira/codes/openmc/output.py +++ b/bluemira/codes/openmc/output.py @@ -540,7 +540,7 @@ class NeutronicsOutputParams(ParameterFrame): peak_div_cu_dpa_rate: Parameter[float] @classmethod - def from_openmc_csg_result(cls, result: OpenMCResult): + def from_openmc_csg_result(cls, result: OpenMCCSGResult): """ Produce output parameters from an OpenMC CSG result """ diff --git a/bluemira/codes/openmc/solver.py b/bluemira/codes/openmc/solver.py index ba7cf9e609..0f24edda6a 100644 --- a/bluemira/codes/openmc/solver.py +++ b/bluemira/codes/openmc/solver.py @@ -13,7 +13,7 @@ from enum import auto from operator import attrgetter from pathlib import Path -from typing import TYPE_CHECKING, Literal, TypeAlias +from typing import TYPE_CHECKING, Any, Literal, TypeAlias import numpy as np import openmc @@ -37,12 +37,15 @@ make_cell_arrays, ) from bluemira.codes.openmc.material import MaterialsLibrary -from bluemira.codes.openmc.output import NeutronicsOutputParams, OpenMCResult +from bluemira.codes.openmc.output import ( + NeutronicsOutputParams, + OpenMCCSGResult, +) from bluemira.codes.openmc.params import ( OpenMCNeutronicsSolverParams, PlasmaSourceParameters, ) -from bluemira.codes.openmc.tallying import filter_cells +from bluemira.codes.openmc.tallying import csg_filter_cells, dagmc_tallys from bluemira.equilibria.equilibrium import Equilibrium if TYPE_CHECKING: @@ -111,7 +114,28 @@ class OpenMCSimulationRuntimeParameters: ] -class OpenMCSetup(CodesSetup): +@dataclass +class PlotConfig: + width: list[float] + pixels: list[int] + basis: str = "xz" + colour_by: str = "cell" + show_overlaps: bool = True + + +@dataclass +class SourceInfo: + rate: float + triton_rate: float + + +@dataclass +class FigureData: + axis: Any + path: Path + + +class OpenMCBaseSetup(CodesSetup, ABC): """Setup task for OpenMC solver""" def __init__( @@ -135,8 +159,142 @@ def __init__( self.source = source self._source_rate = 1.0 self._source_triton_rate = 1.0 - self.blanket_cell_array = cell_arrays.blanket - self.divertor_cell_array = cell_arrays.divertor + self.materials = materials + + def _base_setup(self, run_mode, *, debug: bool = False): + from openmc.config import config # noqa: PLC0415 + + config["cross_sections"] = self.cross_section_xml + + settings = openmc.Settings(run_mode=run_mode.value, output={"summary": False}) + + self.universe, self.geometry = self._create_geometry() + + settings.verbosity = 10 if debug else 6 + + return settings + + def _create_model( + self, settings: openmc.Settings, tallies: openmc.Tallies | None = None + ) -> openmc.Model: + model = openmc.Model(geometry=self.geometry, tallies=tallies, settings=settings) + model.materials = self.materials # .get_all_materials() + return model + + @abstractmethod + def _create_geometry(self) -> tuple[openmc.Universe, openmc.geometry]: ... + + def _create_tallies(self, tally_function: TALLY_FUNCTION_TYPE) -> openmc.Tallies: + tallies_list = [] + for name, scores, filters in tally_function(self.tally_mats, self.tally_geom): + tally = openmc.Tally(name=name) + tally.scores = [scores] + if filters is not None: + tally.filters = filters + tallies_list.append(tally) + + return openmc.Tallies(tallies_list) + + @abstractmethod + def plot( + self, + run_mode, + runtime_params, + eq, + source_params, + tally_function, + *, + debug: bool = False, + ): ... + + @abstractmethod + def volume( + self, + run_mode, + runtime_params, + eq, + source_params, + tally_function, + *, + debug: bool = False, + ): ... + + def run( + self, + run_mode, + runtime_params, + eq, + source_params, + tally_function, + *, + debug: bool = False, + ) -> tuple[openmc.Model, SourceInfo]: + """Run stage for setup openmc""" + settings = self._base_setup(run_mode, debug=debug) + settings.particles = runtime_params.particles + settings.source, source_rate, source_t_rate = self.source(eq, source_params) + settings.batches = int(runtime_params.batches) + settings.photon_transport = runtime_params.photon_transport + settings.electron_treatment = runtime_params.electron_treatment + + tallies = self._create_tallies(tally_function) + + return ( + self._create_model(settings, tallies), + SourceInfo(source_rate, source_t_rate), + ) + + def _plot( + self, run_mode, runtime_params, bounding_box, *, debug: bool = False + ) -> tuple[openmc.Model, PlotConfig]: + """Plot stage for setup openmc""" + settings = self._base_setup(run_mode, debug=debug) + z_max, _z_min, r_max, _r_min = bounding_box + plot_width_0 = r_max * 2.1 + plot_width_1 = z_max * 3.1 + basis = runtime_params.plot_axis + pixels = [ + int(plot_width_0 * runtime_params.plot_pixel_per_metre), + int(plot_width_1 * runtime_params.plot_pixel_per_metre), + ] + width = raw_uc([plot_width_0, plot_width_1], "m", "cm") + + return ( + self._create_model(settings), + PlotConfig(width, pixels, basis), + ) + + def _volume( + self, run_mode, runtime_params, domain, bounding_box, *, debug: bool = False + ) -> tuple[openmc.Model, None]: + """Stochastic volume stage for setup openmc""" + z_max, z_min, r_max, r_min = bounding_box + + min_xyz = (r_min, r_min, z_min) + max_xyz = (r_max, r_max, z_max) + settings = self._base_setup(run_mode, debug=debug) + settings.volume_calculations = openmc.VolumeCalculation( + domain, + runtime_params.particles, + raw_uc(min_xyz, "m", "cm"), + raw_uc(max_xyz, "m", "cm"), + ) + return self._create_model(settings), None + + +class OpenMCCSGSetup(OpenMCBaseSetup): + def __init__( + self, + codes_name: str, + cross_section_xml: str, + eq: Equilibrium, + source, + materials, + cell_arrays, + pre_cell_model, + ): + super().__init__(codes_name, cross_section_xml, eq, source, materials) + self.cell_arrays = cell_arrays self.pre_cell_model = pre_cell_model self.materials = materials self.matlist = attrgetter( @@ -181,9 +339,28 @@ def _base_setup(self, run_mode, *, debug: bool = False): def _set_tallies( self, run_mode, - tally_function: TALLY_FUNCTION_TYPE, - cell_arrays: CellStage, - material_list: list[openmc.Material], + runtime_params, + *_args, + debug: bool = False, + ) -> tuple[openmc.Model, None]: + return self._volume( + run_mode, + runtime_params, + self.cell_arrays.cells, + self.pre_cell_model.bounding_box, + debug=debug, + ) + + +class OpenMCDAGSetup(OpenMCBaseSetup): + def __init__( + self, + codes_name: str, + cross_section_xml: str, + eq: Equilibrium, + source, + materials, + dag_model_path: Path, ): out_path = Path(self.out_path, run_mode.name.lower(), "tallies.xml") tallies_list = [] @@ -258,28 +435,15 @@ def volume( self, run_mode, runtime_params, - _source_params, - _tally_function, - *, + *_args, debug: bool = False, - ): - """Stochastic volume stage for setup openmc""" - z_max, z_min, r_max, r_min = self.pre_cell_model.bounding_box - - min_xyz = (r_min, r_min, z_min) - max_xyz = (r_max, r_max, z_max) - - with self._base_setup(run_mode, debug=debug): - self.settings.volume_calculations = openmc.VolumeCalculation( - self.cells, - runtime_params.particles, - raw_uc(min_xyz, "m", "cm"), - raw_uc(max_xyz, "m", "cm"), - ) - self.files_created.add(Path(self.out_path, run_mode.name.lower(), "summary.h5")) - # single batch - self.files_created.add( - Path(self.out_path, run_mode.name.lower(), "statepoint.1.h5") + ) -> tuple[openmc.Model, None]: + return self._volume( + run_mode, + runtime_params, + self.universe, + self.universe.bounding_box, + debug=debug, ) @@ -291,19 +455,32 @@ def __init__(self, out_path: Path, codes_name: str): self.out_path = out_path - def _run(self, run_mode, *, debug: bool = False): + def _run(self, run_mode, function, **kwargs): """Run openmc""" folder = run_mode.name.lower() - cwd = Path(self.out_path, folder) - cwd.mkdir(parents=True, exist_ok=True) - _timing( - openmc.run, + return _timing( + function, "Executed in", f"Running OpenMC in {folder} mode", debug_info_str=False, )( + openmc_exec="openmc", + **kwargs, + ) + + def run( + self, run_mode, model: openmc.Model, _config: SourceInfo, *, debug: bool = False + ): + """Run stage for run task""" + folder = run_mode.name.lower() + cwd = Path(self.out_path, folder) + cwd.mkdir(parents=True, exist_ok=True) + return self._run( + run_mode, + model.run, + cwd=cwd, + path=cwd, output=debug, - threads=None, geometry_debug=False, restart_file=None, tracks=False, @@ -320,11 +497,35 @@ def run(self, run_mode, *, debug: bool = False): def plot(self, run_mode, *, debug: bool = False): """Plot stage for run task""" - self._run(run_mode, debug=debug) + folder = run_mode.name.lower() + cwd = Path(self.out_path, folder) + cwd.mkdir(parents=True, exist_ok=True) + return FigureData( + self._run( + run_mode, + model.plot, + width=config.width, + pixels=config.pixels, + basis=config.basis, + color_by=config.colour_by, + show_overlaps=config.show_overlaps, + ), + cwd / "geometry.png", + ) def volume(self, run_mode, *, debug: bool = False): """Stochastic volume stage for run task""" - self._run(run_mode, debug=debug) + folder = run_mode.name.lower() + cwd = Path(self.out_path, folder) + cwd.mkdir(parents=True, exist_ok=True) + return self._run( + run_mode, + model.calculate_volumes, + cwd=cwd, + path=cwd, + output=debug, + mpi_args=None, + ) class OpenMCTeardown(CodesTeardown): @@ -366,9 +567,7 @@ def delete_files(files_created): def run( self, universe, - files_created, - source_rate: float, - source_triton_rate: float, + source_info: SourceInfo, statepoint_file, *, delete_files: bool = False, @@ -387,21 +586,14 @@ def run( self.delete_files(files_created) return result, output_params - def plot( - self, - _universe, - files_created, - *_args, - delete_files: bool = False, - ): + def plot(self, _universe, _source_info, fig: FigureData, **kwargs): """Plot stage for Teardown task""" - if delete_files: - self.delete_files(files_created) + fig.axis.get_figure().savefig(fig.path) + return axis def volume( self, _universe, - files_created, _source_params, _statepoint_file, *, @@ -418,6 +610,50 @@ def volume( } +class OpenMCDAGTeardown(CodesTeardown): + def __init__(self, out_path: str, codes_name: str): + super().__init__(None, codes_name) + + self.out_path = out_path + + def run( + self, + universe, + source_info: SourceInfo, + statepoint_file, + *, + delete_files: bool = False, + ): + """Run stage for Teardown task""" + result = OpenMCDAGResult.from_run( + universe, source_info.rate, source_info.triton_rate, statepoint_file + ) + output_params = NeutronicsOutputParams.from_openmc_dag_result(result) + + return result, output_params + + def plot( + self, + _universe, + _source_info, + fig: FigureData, + **_kwargs, + ): + """Plot stage for Teardown task""" + fig.axis.get_figure().save_fig(fig.path) + return axis + + def volume( + self, + _universe, + _source_params, + _statepoint_file, + *, + delete_files: bool = False, + ) -> dict[int, float]: + """Stochastic volume stage for teardown task""" + + TALLY_FUNCTION_TYPE = Callable[ [list[openmc.Material], BlanketCellArray, DivertorCellArray], tuple[ @@ -486,9 +722,8 @@ def tally_function(self) -> TALLY_FUNCTION_TYPE: def tally_function(self, value: TALLY_FUNCTION_TYPE): self._tally_function = value - def execute(self, *, debug=False) -> OpenMCResult | dict[int, float]: + def execute(self, run_mode, *, debug=False) -> OpenMCCSGResult | dict[int, float]: """Execute the setup, run, and teardown tasks, in order.""" - run_mode = self.build_config.get("run_mode", self.run_mode_cls.RUN) if isinstance(run_mode, str): run_mode = self.run_mode_cls.from_string(run_mode) @@ -513,25 +748,7 @@ def _single_run( runtime_params: OpenMCSimulationRuntimeParameters, *, debug=False, - ) -> OpenMCResult | dict[int, float]: - self._setup = self.setup_cls( - self.out_path, - self.name, - str(self.build_config["cross_section_xml"]), - self.eq, - self.source, - self.cell_arrays, - self.pre_cell_model, - self.materials, - ) - self._run = self.run_cls(self.out_path, self.name) - self._teardown = self.teardown_cls( - self.cell_arrays, - self.pre_cell_model, - self.out_path, - self.name, - ) - + ) -> OpenMCCSGResult | dict[int, float]: result = None if setup := self._get_execution_method(self._setup, run_mode): result = setup( @@ -545,15 +762,120 @@ def _single_run( if run := self._get_execution_method(self._run, run_mode): result = run(run_mode, debug=debug) if teardown := self._get_execution_method(self._teardown, run_mode): - result = teardown( - self._setup.universe, - self._setup.files_created, - self._setup._source_rate, - self._setup._source_triton_rate, - Path( - self.out_path, - run_mode.name.lower(), - f"statepoint.{runtime_params.batches}.h5", - ), - ) + result = teardown(self._setup.universe, config, result) return result + + +class OpenMCCSGNeutronicsSolver(OpenMCNeutronicsSolver): + setup_cls: type[OpenMCCSGSetup] = OpenMCCSGSetup + teardown_cls: type[CodesTeardown] = OpenMCCSGTeardown + + def __init__( + self, + params: dict | ParameterFrame, + build_config: dict, + eq: Equilibrium, + source: NeutronSourceCreator, + neutronics_model: NeutronicsReactor, + op_cond: OperationalConditions, + tally_function: TALLY_FUNCTION_TYPE | None = None, + ): + super().__init__( + params, + build_config, + eq, + source, + ) + self.neutronics_model = neutronics_model + self.materials = MaterialsLibrary.from_neutronics_materials( + self.neutronics_model.material_library, op_cond + ) + + self.cell_arrays = make_cell_arrays( + self.neutronics_model, + BluemiraNeutronicsCSG(), + self.materials, + control_id=True, + ) + + self.tally_function = ( + csg_filter_cells if tally_function is None else tally_function + ) + + def _single_run( + self, + run_mode: OpenMCRunModes, + source_params: PlasmaSourceParameters, + runtime_params: OpenMCSimulationRuntimeParameters, + *, + debug=False, + ) -> OpenMCCSGResult | dict[int, float]: + self._setup = self.setup_cls( + self.name, + str(self.build_config["cross_section_xml"]), + self.eq, + self.source, + self.materials, + self.cell_arrays, + self.neutronics_model, + ) + + self._run = self.run_cls(self.out_path, self.name) + self._teardown = self.teardown_cls( + self.cell_arrays, + self.neutronics_model, + self.out_path, + self.name, + ) + + return super()._single_run(run_mode, source_params, runtime_params, debug=debug) + + +class OpenMCDAGMCNeutronicsSolver(OpenMCNeutronicsSolver): + setup_cls: type[OpenMCDAGSetup] = OpenMCDAGSetup + teardown_cls: type[CodesTeardown] = OpenMCDAGTeardown + + def __init__( + self, + params: dict | ParameterFrame, + build_config: dict, + eq: Equilibrium, + source: NeutronSourceCreator, + dagmc_model_path: Path, + materials, + tally_function: TALLY_FUNCTION_TYPE | None = None, + ): + super().__init__( + params, + build_config, + eq, + source, + ) + self.dagmc_model_path = dagmc_model_path + self.materials = materials + + self.tally_function = dagmc_tallys if tally_function is None else tally_function + + def _single_run( + self, + run_mode: OpenMCRunModes, + source_params: PlasmaSourceParameters, + runtime_params: OpenMCSimulationRuntimeParameters, + *, + debug=False, + ) -> OpenMCCSGResult | dict[int, float]: + self._setup = self.setup_cls( + self.name, + str(self.build_config["cross_section_xml"]), + self.eq, + self.source, + self.materials, + self.dagmc_model_path, + ) + + self._run = self.run_cls(self.out_path, self.name) + self._teardown = self.teardown_cls( + self.out_path, + self.name, + ) + return super()._single_run(run_mode, source_params, runtime_params, debug=debug) diff --git a/bluemira/codes/openmc/tallying.py b/bluemira/codes/openmc/tallying.py index e194aceb25..78d902abbe 100644 --- a/bluemira/codes/openmc/tallying.py +++ b/bluemira/codes/openmc/tallying.py @@ -94,3 +94,24 @@ def filter_cells( ("damage", "damage-energy", [cell_filter]), # used to get the EUROFER OBMP ) + + +def dagmc_tallys(material_list, model): + heating_cell_tally = openmc.Tally(name="heating") + heating_cell_tally.scores = ["heating"] + + # record the total TBR + tbr_cell_tally = openmc.Tally(name="tbr") + tbr_cell_tally.scores = ["(n,Xt)"] + + # mesh that covers the geometry + mesh = openmc.RegularMesh.from_domain(model, dimension=(100, 100, 100)) + mesh_filter = openmc.MeshFilter(mesh) + + return [ + ("heating", "heating", None), + ("heating_on_mesh", "heating", [mesh_filter]), + ("tbr", "(n,Xt)", None), + ("tbr_on_mesh", "(n,Xt)", [mesh_filter]), + ("flux_on_mesh", "flux", [mesh_filter]), + ] diff --git a/bluemira/codes/wrapper.py b/bluemira/codes/wrapper.py index 263407e79c..0f90b0318b 100644 --- a/bluemira/codes/wrapper.py +++ b/bluemira/codes/wrapper.py @@ -120,5 +120,5 @@ def neutronics_code_solver( """ neutron = get_code_interface(module) return neutron.Solver( - params, build_config, neutronics_model, eq, source, op_cond, tally_function + params, build_config, eq, source, neutronics_model, op_cond, tally_function ) diff --git a/eudemo/eudemo/model_managers.py b/eudemo/eudemo/model_managers.py index ed6cfdae29..38c515a04f 100644 --- a/eudemo/eudemo/model_managers.py +++ b/eudemo/eudemo/model_managers.py @@ -21,7 +21,7 @@ from bluemira.base.look_and_feel import bluemira_warn if TYPE_CHECKING: - from bluemira.codes.openmc.output import OpenMCResult + from bluemira.codes.openmc.output import OpenMCCSGResult from bluemira.equilibria.run import Snapshot from eudemo.eudemo.neutronics.run import EUDEMONeutronicsCSGReactor @@ -128,7 +128,7 @@ class NeutronicsManager: def __init__( self, csg_reactor: EUDEMONeutronicsCSGReactor, - results: OpenMCResult | dict[int, float], + results: OpenMCCSGResult | dict[int, float], ): self.csg_reactor = csg_reactor self.results = results diff --git a/eudemo/eudemo/neutronics/run.py b/eudemo/eudemo/neutronics/run.py index afd7842e2e..41b79190aa 100644 --- a/eudemo/eudemo/neutronics/run.py +++ b/eudemo/eudemo/neutronics/run.py @@ -27,7 +27,7 @@ from bluemira.base.parameter_frame import ParameterFrame from bluemira.base.reactor import ComponentManager - from bluemira.codes.openmc.output import OpenMCResult + from bluemira.codes.openmc.output import OpenMCCSGResult from bluemira.codes.openmc.solver import NeutronSourceCreator from bluemira.equilibria.equilibrium import Equilibrium from bluemira.geometry.wire import BluemiraWire @@ -55,7 +55,7 @@ def _get_wires_from_components( def run_neutronics( - params: dict | ParameterFrame, + params: ParameterFrame, build_config: dict, blanket: ComponentManager, vacuum_vessel: ComponentManager, @@ -64,7 +64,7 @@ def run_neutronics( op_cond: OperationalConditions, source: NeutronSourceCreator | None = None, tally_function=None, -) -> tuple[EUDEMONeutronicsCSGReactor, OpenMCResult | dict[int, float]]: +) -> tuple[EUDEMONeutronicsCSGReactor, OpenMCCSGResult | dict[int, float]]: """Runs the neutronics model Returns @@ -110,8 +110,12 @@ def run_neutronics( tally_function=tally_function, ) - res, new_params = solver.execute() + outputs = solver.execute(build_config.get("run_mode", "run")) - params.update_from_frame(new_params) + if len(outputs) == 2: + res = outputs[0] + params.update_from_frame(outputs[1]) + else: + res = outputs return neutronics_csg, res diff --git a/eudemo/eudemo/reactor.py b/eudemo/eudemo/reactor.py index e170fe1a10..11d5bf241c 100644 --- a/eudemo/eudemo/reactor.py +++ b/eudemo/eudemo/reactor.py @@ -887,6 +887,18 @@ def save_reactor(reactor, reactor_config, folder_name): reactor_config.config_for("CAD_Neutronics"), ) + from bluemira.codes.openmc.sovler import OpenMCDAGMCNeutronicsSolver + + OpenMCDAGMCNeutronicsSolver( + params, build_config, eq, source, neutronics_model, op_cond, tally_function + ) + + debug = [upper_port_koz_xz, eq_port_koz_xz, lower_port_koz_xz] + debug.extend(reactor.pf_coils.xz_boundary) + # I know there are clashes, I need to put in dynamic bounds on position opt to + # include coil XS. + show_cad(debug) + # reactor.plot("xz") # reactor.show_cad(n_sectors=2) diff --git a/examples/radiation_transport/run_cad_neutronics.ex.py b/examples/radiation_transport/run_cad_neutronics.ex.py index 6a7e733e51..3dad725ed2 100644 --- a/examples/radiation_transport/run_cad_neutronics.ex.py +++ b/examples/radiation_transport/run_cad_neutronics.ex.py @@ -25,15 +25,14 @@ """ # %% -import json from pathlib import Path import openmc import vtk from vtkmodules.util import numpy_support -from bluemira.base.file import get_bluemira_root from bluemira.base.look_and_feel import bluemira_print +from bluemira.codes.openmc.params import OpenMCNeutronicsSolverParams from bluemira.codes.openmc.sources import create_ring_source from bluemira.materials.cache import establish_material_cache, get_cached_material @@ -50,11 +49,11 @@ dag_model_path = par / "OptimisedReactor.h5m" meta_data_path = par / "OptimisedReactor.meta.json" -# %% -establish_material_cache([ - Path(get_bluemira_root()) / "data" / "materials" / "materials.json", - Path(get_bluemira_root()) / "data" / "materials" / "mixtures.json", -]) +# # %% +# establish_material_cache([ +# Path(get_bluemira_root()) / "data" / "materials" / "materials.json", +# Path(get_bluemira_root()) / "data" / "materials" / "mixtures.json", +# ]) omc_output_path = par / "omc" # Ensure OpenMC output directory exists omc_output_path.mkdir(parents=True, exist_ok=True) @@ -63,23 +62,293 @@ # ## Running the DAGMC model in OpenMC # %% +# --- +# jupyter: +# jupytext: +# cell_metadata_filter: tags,-all +# notebook_metadata_filter: -jupytext.text_representation.jupytext_version +# text_representation: +# extension: .py +# format_name: percent +# format_version: '1.3' +# kernelspec: +# display_name: Python 3 (ipykernel) +# language: python +# name: python3 +# --- + +# %% tags=["remove-cell"] +# SPDX-FileCopyrightText: 2021-present M. Coleman, J. Cook, F. Franza +# SPDX-FileCopyrightText: 2021-present I.A. Maione, S. McIntosh +# SPDX-FileCopyrightText: 2021-present J. Morris, D. Short +# +# SPDX-License-Identifier: LGPL-2.1-or-later + +""" +Optimised reactor example, showing how to export a reactor +to a DAGMC model and run it in OpenMC. +""" + +# %% +import math +from dataclasses import dataclass +from pathlib import Path + +import matplotlib.pyplot as plt +import numpy as np + +from bluemira.base.builder import Builder +from bluemira.base.components import Component, PhysicalComponent +from bluemira.base.designer import Designer +from bluemira.base.file import get_bluemira_root +from bluemira.base.logs import set_log_level +from bluemira.base.look_and_feel import bluemira_print +from bluemira.base.parameter_frame import Parameter, ParameterFrame +from bluemira.base.reactor import ComponentManager, Reactor +from bluemira.builders.plasma import Plasma, PlasmaBuilder +from bluemira.builders.tools import apply_component_display_options +from bluemira.display import show_cad +from bluemira.display.palettes import BLUE_PALETTE +from bluemira.display.plotter import PlotOptions, plot_2d +from bluemira.equilibria.coils import Coil, CoilSet +from bluemira.equilibria.diagnostics import PicardDiagnostic, PicardDiagnosticOptions +from bluemira.equilibria.equilibrium import Equilibrium +from bluemira.equilibria.grid import Grid +from bluemira.equilibria.optimisation.constraints import ( + FieldNullConstraint, + IsofluxConstraint, + MagneticConstraintSet, +) +from bluemira.equilibria.optimisation.problem import ( + UnconstrainedTikhonovCurrentGradientCOP, +) +from bluemira.equilibria.profiles import CustomProfile +from bluemira.equilibria.shapes import JohnerLCFS +from bluemira.equilibria.solve import PicardIterator +from bluemira.geometry.face import BluemiraFace +from bluemira.geometry.optimisation import optimise_geometry +from bluemira.geometry.parameterisations import PrincetonD +from bluemira.geometry.tools import ( + distance_to, + interpolate_bspline, + make_polygon, + offset_wire, + revolve_shape, + sweep_shape, +) +from bluemira.geometry.wire import BluemiraWire +from bluemira.materials.cache import establish_material_cache +from bluemira.optimisation import Algorithm + + +from bluemira.codes.openmc.solver import OpenMCDAGMCNeutronicsSolver + + +def _lcfs_parameterisation(R_0, A): + return JohnerLCFS({ + "r_0": {"value": R_0}, + "z_0": {"value": 0.0}, + "a": {"value": R_0 / A}, + "kappa_u": {"value": 1.6}, + "kappa_l": {"value": 1.9}, + "delta_u": {"value": 0.4}, + "delta_l": {"value": 0.4}, + "phi_u_neg": {"value": 0.0}, + "phi_u_pos": {"value": 0.0}, + "phi_l_neg": {"value": 45.0}, + "phi_l_pos": {"value": 30.0}, + }) + + +def ref_eq(R_0, A) -> Equilibrium: # noqa: D103 + x = [5.4, 14.0, 17.75, 17.75, 14.0, 7.0, 2.77, 2.77, 2.77, 2.77, 2.77] + z = [9.26, 7.9, 2.5, -2.5, -7.9, -10.5, 7.07, 4.08, -0.4, -4.88, -7.86] + dx = [0.6, 0.7, 0.5, 0.5, 0.7, 1.0, 0.4, 0.4, 0.4, 0.4, 0.4] + dz = [0.6, 0.7, 0.5, 0.5, 0.7, 1.0, 2.99 / 2, 2.99 / 2, 5.97 / 2, 2.99 / 2, 2.99 / 2] + + coils = [] + j = 1 + for i, (xi, zi, dxi, dzi) in enumerate(zip(x, z, dx, dz, strict=False)): + if j > 6: # noqa: PLR2004 + j = 1 + ctype = "PF" if i < 6 else "CS" # noqa: PLR2004 + coil = Coil( + xi, + zi, + current=0, + dx=dxi, + dz=dzi, + ctype=ctype, + name=f"{ctype}_{j}", + ) + coils.append(coil) + j += 1 + + coilset = CoilSet(*coils) + coilset.assign_material("CS", j_max=16.5e6, b_max=12.5) + coilset.assign_material("PF", j_max=12.5e6, b_max=11.0) + + cs = coilset.get_coiltype("CS") + cs.fix_sizes() + cs.discretisation = 0.3 + + B_0 = 4.8901 # T + I_p = 19.07e6 # A + + grid = Grid(3.0, 13.0, -10.0, 10.0, 65, 65) + + profiles = CustomProfile( + np.array([ + 86856, + 86506, + 84731, + 80784, + 74159, + 64576, + 52030, + 36918, + 20314, + 4807, + 0.0, + ]), + -np.array([ + 0.125, + 0.124, + 0.122, + 0.116, + 0.106, + 0.093, + 0.074, + 0.053, + 0.029, + 0.007, + 0.0, + ]), + R_0=R_0, + B_0=B_0, + I_p=I_p, + ) + + eq = Equilibrium(coilset, grid, profiles, psi=None) + + lcfs = ( + _lcfs_parameterisation(R_0, A).create_shape().discretise(byedges=True, ndiscr=50) + ) + + x_bdry, z_bdry = lcfs.x, lcfs.z + arg_inner = np.argmin(x_bdry) + + isoflux = IsofluxConstraint( + x_bdry, + z_bdry, + x_bdry[arg_inner], + z_bdry[arg_inner], + tolerance=0.5, # Difficult to choose... + constraint_value=0.0, # Difficult to choose... + ) + + xp_idx = np.argmin(z_bdry) + x_point = FieldNullConstraint( + x_bdry[xp_idx], + z_bdry[xp_idx], + tolerance=1e-4, # [T] + ) + current_opt_problem = UnconstrainedTikhonovCurrentGradientCOP( + coilset, eq, MagneticConstraintSet([isoflux, x_point]), gamma=1e-7 + ) + diagnostic_plotting = PicardDiagnosticOptions(plot=PicardDiagnostic.EQ) + program = PicardIterator( + eq, + current_opt_problem, + fixed_coils=True, + relaxation=0.2, + diagnostic_plotting=diagnostic_plotting, + ) + program() + + plt.show() + + return eq + + +eq = ref_eq(9, 3) + +"Toroidal_Field_Coil_2015" +"Homogenised_HCPB_2015_v3_BZ" + +establish_material_cache(["eurofusion_materials.library", "matproplib"]) + +openmc_mats = [ + get_cached_material(mat_name).convert( + "openmc", op_cond={"temperature": 294, "pressure": 101325} + ) + for mat_name in ["Toroidal_Field_Coil_2015", "Homogenised_HCPB_2015_v3_BZ"] +] + +params = OpenMCNeutronicsSolverParams( + R_0=Parameter(name="R_0", value=9, unit="m"), + profile_rho_ped=Parameter( + name="profile_rho_ped", + value=0.94, + unit="dimensionless", + source="Input", + long_name="Pedestal location in normalized radius", + ), + P_fus=Parameter( + name="P_fus", + value=2000, + unit="megawatt", + source="Input", + long_name="Total fusion power", + ), + n_profile_alpha=Parameter( + name="n_profile_alpha", value=1.0, unit="dimensionless", source="Input" + ), + n_e_core=Parameter(name="n_e_core", value=1.5e20, unit="1/m^3"), + n_e_ped=Parameter(name="n_e_ped", value=8e19, unit="1/m^3"), + n_e_sep=Parameter(name="n_e_sep", value=3e19, unit="1/m^3"), + T_profile_alpha=Parameter(name="T_profile_alpha", value=1.45, unit="dimensionless"), + T_profile_beta=Parameter(name="T_profile_beta", value=2.0, unit="dimensionless"), + T_e_core=Parameter(name="T_e_core", value=20, unit="kiloelectron_volt"), + T_e_ped=Parameter(name="T_e_ped", value=5.5, unit="kiloelectron_volt"), + T_e_sep=Parameter(name="T_e_sep", value=0.1, unit="kiloelectron_volt"), + T_ie_ratio=Parameter(name="T_ie_ratio", value=1, unit="dimensionless"), + n_i_fuel=Parameter(name="n_i_fuel", value=8e19, unit="1/m^3"), + n_e=Parameter(name="n_e", value=1e20, unit="1/m^3"), + shaf_shift=Parameter(name="shaf_shift", value=0.5, unit="meter"), +) + +solver = OpenMCDAGMCNeutronicsSolver( + params, + { + "particles": 10000, + "cross_section_xml": Path( + "eudemo/config/cross_section_data/cross_sections.xml" + ).resolve(), + }, + eq, + lambda *args: (create_ring_source(ring_source_major_radius, 0), 1, 1), + dag_model_path, + openmc_mats, +) +solver.execute("run") +import pdb + +pdb.set_trace() dagmc_univ = openmc.DAGMCUniverse( filename=dag_model_path.as_posix(), auto_geom_ids=True, ).bounded_universe() -# load model materials -with open(meta_data_path) as meta_file: - bom = json.load(meta_file)["bom"] -openmc_mats = [ - get_cached_material(mat_name).convert("openmc", op_cond={"temperature": 294}) - for mat_name in bom -] + +# # load model materials +# with open(meta_data_path) as meta_file: +# bom = json.load(meta_file)["bom"] # load DAG model geometry = openmc.Geometry(dagmc_univ) -my_source = create_ring_source(ring_source_major_radius, 0) +# my_source = settings = openmc.Settings() settings.batches = n_batches From f78539a03e79f01f98555bcad2df5841bfc0f2cf Mon Sep 17 00:00:00 2001 From: Ocean Date: Mon, 20 Oct 2025 16:58:27 +0100 Subject: [PATCH 02/17] Created the correct tally for total power. (#4037) * fix upper port height * P_fus_DT and e_mult * ruff stuff * fix total heating and pass in pfus * Created the correct tally for total power. * Fixed a duplicate name * typo correction * Total power is now consistent with what you get when adding up vacuum vessel power + blanket power + divertor power. * Added the appropriate NWL tallies (#4040) * Renamed the Setup, Run, and Teardown classes to something more specific. * Changed the tallies to reflect the NWL properly. * Add comments/renamed functions/added arguments to make the code more readable. * Includes the previously-missed cells in the total power tallies now. * Bug fix * Added comment to clarify flux tally's 'mean' is actually volume integrated. * Added tallies and calculations to get the damage values in DPA/FPY as well. * Removed all traces of the word neutron wall load/NWL because its definition is too ambiguous/not useful. * removed the calculation about fw heating as well as we aren't using it. (Can be undone if necessary). * removed unused function --------- Co-authored-by: matti * rebase pain * omg no deliberate nans * clean up diagnostic printing * remove duplicate function * ruff * avoid pandas warnings * ruff * change DPACoefficients into function * use * typing --------- Co-authored-by: Matti --- bluemira/radiation_transport/neutronics/constants.py | 1 + 1 file changed, 1 insertion(+) diff --git a/bluemira/radiation_transport/neutronics/constants.py b/bluemira/radiation_transport/neutronics/constants.py index 52ea07477c..e744c49545 100644 --- a/bluemira/radiation_transport/neutronics/constants.py +++ b/bluemira/radiation_transport/neutronics/constants.py @@ -7,6 +7,7 @@ import numpy.typing as npt from periodictable import elements +import numpy.typing as npt from bluemira.base.constants import ( ELECTRON_MOLAR_MASS, From 177d818e65f1fdd983a075c692be9b78815603aa Mon Sep 17 00:00:00 2001 From: james <81617086+je-cook@users.noreply.github.com> Date: Wed, 12 Nov 2025 21:44:37 +0000 Subject: [PATCH 03/17] get to running dagmc --- bluemira/base/builder.py | 2 +- bluemira/base/reactor.py | 5 ++ bluemira/codes/openmc/solver.py | 7 ++- .../neutronics/constants.py | 1 - eudemo/eudemo/neutronics/run.py | 51 +++++++++++++++++- eudemo/eudemo/reactor.py | 53 ++----------------- .../run_cad_neutronics.ex.py | 34 +----------- 7 files changed, 68 insertions(+), 85 deletions(-) diff --git a/bluemira/base/builder.py b/bluemira/base/builder.py index 738f7e50e4..7bfb6cbb5a 100644 --- a/bluemira/base/builder.py +++ b/bluemira/base/builder.py @@ -136,7 +136,7 @@ def get_material(self, component_name: str | None = None) -> Material | None: f"No corresponding material found for {component_name} in {mats}" ) - return mat + return mat() if callable(mat) else mat def component_tree( self, diff --git a/bluemira/base/reactor.py b/bluemira/base/reactor.py index 5abe97d846..61c7db164d 100644 --- a/bluemira/base/reactor.py +++ b/bluemira/base/reactor.py @@ -133,6 +133,11 @@ def tree(self) -> str: """ return self.component().tree() + @property + def materials(self): + """Get all materials in the Manager""" + return self.component().get_component_propertues("material", first=False) + @staticmethod def _validate_cad_dim(dim: DIM_3D | DIM_2D) -> None: """ diff --git a/bluemira/codes/openmc/solver.py b/bluemira/codes/openmc/solver.py index 0f24edda6a..5a96d9b972 100644 --- a/bluemira/codes/openmc/solver.py +++ b/bluemira/codes/openmc/solver.py @@ -178,7 +178,10 @@ def _create_model( self, settings: openmc.Settings, tallies: openmc.Tallies | None = None ) -> openmc.Model: model = openmc.Model(geometry=self.geometry, tallies=tallies, settings=settings) - model.materials = self.materials # .get_all_materials() + if isinstance(self.materials, MaterialsLibrary): + model.materials = self.materials.get_all_materials() + else: + model.materials = self.materials return model @abstractmethod @@ -589,7 +592,7 @@ def run( def plot(self, _universe, _source_info, fig: FigureData, **kwargs): """Plot stage for Teardown task""" fig.axis.get_figure().savefig(fig.path) - return axis + return fig.axis def volume( self, diff --git a/bluemira/radiation_transport/neutronics/constants.py b/bluemira/radiation_transport/neutronics/constants.py index e744c49545..52ea07477c 100644 --- a/bluemira/radiation_transport/neutronics/constants.py +++ b/bluemira/radiation_transport/neutronics/constants.py @@ -7,7 +7,6 @@ import numpy.typing as npt from periodictable import elements -import numpy.typing as npt from bluemira.base.constants import ( ELECTRON_MOLAR_MASS, diff --git a/eudemo/eudemo/neutronics/run.py b/eudemo/eudemo/neutronics/run.py index 41b79190aa..7f38ec105e 100644 --- a/eudemo/eudemo/neutronics/run.py +++ b/eudemo/eudemo/neutronics/run.py @@ -7,8 +7,10 @@ from __future__ import annotations +from pathlib import Path from typing import TYPE_CHECKING +from bluemira.codes.openmc.solver import OpenMCDAGMCNeutronicsSolver from bluemira.codes.openmc.sources import make_tokamak_source from bluemira.codes.wrapper import neutronics_code_solver from bluemira.radiation_transport.neutronics.blanket_data import ( @@ -54,7 +56,7 @@ def _get_wires_from_components( ) -def run_neutronics( +def run_csg_neutronics( params: ParameterFrame, build_config: dict, blanket: ComponentManager, @@ -119,3 +121,50 @@ def run_neutronics( res = outputs return neutronics_csg, res + + +def export_dagmc_model(reactor, build_config): + """ + Export the reactor model to a DAGMC model. + + Parameters + ---------- + reactor : EUDEMO + The reactor instance to export. + build_config : dict + The build configuration parameters. + """ + if build_config.get("export_dagmc_model", False): + reactor.save_cad( + directory=build_config.get("dagmc_export_dir", None), + cad_format="dagmc", + construction_params={ + "without_components": [reactor.plasma], + "group_by_materials": True, + }, + ) + + +def run_dagmc_neutronics( + reactor, + params, + build_config, + eq: Equilibrium, + source: NeutronSourceCreator | None = None, + tally_function=None, +): + export_dagmc_model(reactor, build_config) + + solver = OpenMCDAGMCNeutronicsSolver( + params, + build_config, + eq, + source=source or make_tokamak_source, + dagmc_model_path=build_config.get("dagmc_export_dir", Path.cwd()), + materials=reactor.materials, + tally_function=tally_function, + ) + + outputs = solver.execute(build_config.get("run_mode", "run")) + + return outputs diff --git a/eudemo/eudemo/reactor.py b/eudemo/eudemo/reactor.py index 11d5bf241c..f54e4e3fde 100644 --- a/eudemo/eudemo/reactor.py +++ b/eudemo/eudemo/reactor.py @@ -96,7 +96,7 @@ ) from eudemo.maintenance.upper_port import UpperPortKOZDesigner from eudemo.model_managers import EquilibriumManager, NeutronicsManager -from eudemo.neutronics.run import run_neutronics +from eudemo.neutronics.run import run_csg_neutronics, run_dagmc_neutronics from eudemo.params import EUDEMOReactorParams from eudemo.pf_coils import PFCoil, PFCoilsDesigner, build_pf_coils_component from eudemo.power_cycle import SteadyStatePowerCycleSolver @@ -534,28 +534,6 @@ def build_radiation_plugs( return builder.build() -def export_dagmc_model(reactor: EUDEMO, build_config): - """ - Export the reactor model to a DAGMC model. - - Parameters - ---------- - reactor : EUDEMO - The reactor instance to export. - build_config : dict - The build configuration parameters. - """ - if build_config.get("export_dagmc_model", False): - reactor.save_cad( - directory=build_config.get("dagmc_export_dir", None), - cad_format="dagmc", - construction_params={ - "without_components": [reactor.plasma], - "group_by_materials": True, - }, - ) - - def save_reactor(reactor, reactor_config, folder_name): """ Save a reactor to a folder data-structure @@ -715,7 +693,7 @@ def save_reactor(reactor, reactor_config, folder_name): if reactor_config.config_for("Neutronics").get("enabled", False): neutronics_start = time.time() reactor.neutronics = NeutronicsManager( - *run_neutronics( + *run_csg_neutronics( reactor_config.params_for("Neutronics").global_params, reactor_config.config_for("Neutronics"), blanket=reactor.blanket, @@ -882,34 +860,13 @@ def save_reactor(reactor, reactor_config, folder_name): n_TF=reactor_config.global_params.n_TF.value, ) - export_dagmc_model( + run_dagmc_neutronics( reactor, + reactor_config.params_for("CAD_Neutronics"), reactor_config.config_for("CAD_Neutronics"), + reference_eq, ) - from bluemira.codes.openmc.sovler import OpenMCDAGMCNeutronicsSolver - - OpenMCDAGMCNeutronicsSolver( - params, build_config, eq, source, neutronics_model, op_cond, tally_function - ) - - debug = [upper_port_koz_xz, eq_port_koz_xz, lower_port_koz_xz] - debug.extend(reactor.pf_coils.xz_boundary) - # I know there are clashes, I need to put in dynamic bounds on position opt to - # include coil XS. - show_cad(debug) - - # reactor.plot("xz") - # reactor.show_cad(n_sectors=2) - - sspc_solver = SteadyStatePowerCycleSolver(reactor_config.global_params) - sspc_result = sspc_solver.execute() - reactor_config.global_params.P_el_net.set_value( - sspc_result["P_el_net"], "BLUEMIRA" - ) - - lcfs = ClosedFluxSurface(reference_eq.get_LCFS()) - reactor_config.global_params.V_p.set_value(lcfs.volume, "BLUEMIRA") end = time.time() diff --git a/examples/radiation_transport/run_cad_neutronics.ex.py b/examples/radiation_transport/run_cad_neutronics.ex.py index 3dad725ed2..deeeb97d23 100644 --- a/examples/radiation_transport/run_cad_neutronics.ex.py +++ b/examples/radiation_transport/run_cad_neutronics.ex.py @@ -90,26 +90,13 @@ """ # %% -import math -from dataclasses import dataclass from pathlib import Path import matplotlib.pyplot as plt import numpy as np -from bluemira.base.builder import Builder -from bluemira.base.components import Component, PhysicalComponent -from bluemira.base.designer import Designer -from bluemira.base.file import get_bluemira_root -from bluemira.base.logs import set_log_level -from bluemira.base.look_and_feel import bluemira_print -from bluemira.base.parameter_frame import Parameter, ParameterFrame -from bluemira.base.reactor import ComponentManager, Reactor -from bluemira.builders.plasma import Plasma, PlasmaBuilder -from bluemira.builders.tools import apply_component_display_options -from bluemira.display import show_cad -from bluemira.display.palettes import BLUE_PALETTE -from bluemira.display.plotter import PlotOptions, plot_2d +from bluemira.base.parameter_frame import Parameter +from bluemira.codes.openmc.solver import OpenMCDAGMCNeutronicsSolver from bluemira.equilibria.coils import Coil, CoilSet from bluemira.equilibria.diagnostics import PicardDiagnostic, PicardDiagnosticOptions from bluemira.equilibria.equilibrium import Equilibrium @@ -125,23 +112,6 @@ from bluemira.equilibria.profiles import CustomProfile from bluemira.equilibria.shapes import JohnerLCFS from bluemira.equilibria.solve import PicardIterator -from bluemira.geometry.face import BluemiraFace -from bluemira.geometry.optimisation import optimise_geometry -from bluemira.geometry.parameterisations import PrincetonD -from bluemira.geometry.tools import ( - distance_to, - interpolate_bspline, - make_polygon, - offset_wire, - revolve_shape, - sweep_shape, -) -from bluemira.geometry.wire import BluemiraWire -from bluemira.materials.cache import establish_material_cache -from bluemira.optimisation import Algorithm - - -from bluemira.codes.openmc.solver import OpenMCDAGMCNeutronicsSolver def _lcfs_parameterisation(R_0, A): From 25d2a7c0c67c6e6560cb7ed289a059921987ad04 Mon Sep 17 00:00:00 2001 From: james <81617086+je-cook@users.noreply.github.com> Date: Thu, 13 Nov 2025 11:04:21 +0000 Subject: [PATCH 04/17] =?UTF-8?q?=F0=9F=90=9B=20Fix=20materials=20getter?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- bluemira/base/reactor.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/bluemira/base/reactor.py b/bluemira/base/reactor.py index 61c7db164d..865663cf59 100644 --- a/bluemira/base/reactor.py +++ b/bluemira/base/reactor.py @@ -136,7 +136,7 @@ def tree(self) -> str: @property def materials(self): """Get all materials in the Manager""" - return self.component().get_component_propertues("material", first=False) + return self.component().get_component_properties("material", first=False)[0] @staticmethod def _validate_cad_dim(dim: DIM_3D | DIM_2D) -> None: From 8ed595737f7fa1c1897c94fdc0a3d097c8094e75 Mon Sep 17 00:00:00 2001 From: james <81617086+je-cook@users.noreply.github.com> Date: Thu, 13 Nov 2025 11:20:46 +0000 Subject: [PATCH 05/17] =?UTF-8?q?=E2=9C=A8=20Dagmc=20running?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- eudemo/config/build_config.json | 39 +++++++++++++++++++-------------- eudemo/eudemo/neutronics/run.py | 28 +++++++++++++++++------ eudemo/eudemo/reactor.py | 13 ++++++++--- 3 files changed, 53 insertions(+), 27 deletions(-) diff --git a/eudemo/config/build_config.json b/eudemo/config/build_config.json index d71201aa1f..a1ebd1162d 100644 --- a/eudemo/config/build_config.json +++ b/eudemo/config/build_config.json @@ -118,23 +118,28 @@ "material": "Homogenised_Divertor_2015" }, "Neutronics": { - "enabled": true, - "cross_section_xml": "$path_expand:./cross_section_data/cross_sections.xml", - "particles": 10000, - "batches": 2, - "photon_transport": true, - "electron_treatment": "ttb", - "run_mode": "run_and_plot", - "openmc_write_summary": false, - "blanket_type": "HCPB", - "plot_axis": "xz", - "plot_pixel_per_metre": 100, - "neutronics_output_path": "$path_expand:./neutronics", - "show_data": false - }, - "CAD_Neutronics": { - "export_dagmc_model": false, - "dagmc_export_dir": "$path_expand:./dagmc_export" + "CSG": { + "enabled": false, + "cross_section_xml": "$path_expand:./cross_section_data/cross_sections.xml", + "particles": 10000, + "batches": 2, + "photon_transport": true, + "electron_treatment": "ttb", + "run_mode": "run_and_plot", + "openmc_write_summary": false, + "blanket_type": "HCPB", + "plot_axis": "xz", + "plot_pixel_per_metre": 100, + "neutronics_output_path": "$path_expand:./neutronics", + "show_data": false + }, + "DAGMC": { + "export_dagmc_model": false, + "dagmc_export_dir": "$path_expand:./dagmc_export", + "cross_section_xml": "$path_expand:./cross_section_data/cross_sections.xml", + "particles": 10000, + "batches": 2 + } }, "TF coils": { "run_mode": "run", diff --git a/eudemo/eudemo/neutronics/run.py b/eudemo/eudemo/neutronics/run.py index 7f38ec105e..6538344877 100644 --- a/eudemo/eudemo/neutronics/run.py +++ b/eudemo/eudemo/neutronics/run.py @@ -10,6 +10,8 @@ from pathlib import Path from typing import TYPE_CHECKING +from matproplib.library.fluids import Void + from bluemira.codes.openmc.solver import OpenMCDAGMCNeutronicsSolver from bluemira.codes.openmc.sources import make_tokamak_source from bluemira.codes.wrapper import neutronics_code_solver @@ -114,7 +116,7 @@ def run_csg_neutronics( outputs = solver.execute(build_config.get("run_mode", "run")) - if len(outputs) == 2: + if len(outputs) == 2: # noqa: PLR2004 res = outputs[0] params.update_from_frame(outputs[1]) else: @@ -139,7 +141,10 @@ def export_dagmc_model(reactor, build_config): directory=build_config.get("dagmc_export_dir", None), cad_format="dagmc", construction_params={ - "without_components": [reactor.plasma], + "without_components": [ + reactor.plasma, + reactor.blanket, + ], "group_by_materials": True, }, ) @@ -153,18 +158,27 @@ def run_dagmc_neutronics( source: NeutronSourceCreator | None = None, tally_function=None, ): + """Creates and runs the DAGMC neutronics model""" # noqa: DOC201 export_dagmc_model(reactor, build_config) + mats = {"undef_material": Void(name="undef_material")} + for m in reactor.materials: + if m is not None and m.name not in mats: + mats[m.name] = m + solver = OpenMCDAGMCNeutronicsSolver( params, build_config, eq, source=source or make_tokamak_source, - dagmc_model_path=build_config.get("dagmc_export_dir", Path.cwd()), - materials=reactor.materials, + dagmc_model_path=Path( + build_config.get("dagmc_export_dir", Path.cwd()), f"{reactor.name}.h5m" + ), + materials=[ + m.convert("openmc", {"temperature": 298, "pressure": 101325}) + for m in mats.values() + ], tally_function=tally_function, ) - outputs = solver.execute(build_config.get("run_mode", "run")) - - return outputs + return solver.execute(build_config.get("run_mode", "run")) diff --git a/eudemo/eudemo/reactor.py b/eudemo/eudemo/reactor.py index f54e4e3fde..e2c7aebace 100644 --- a/eudemo/eudemo/reactor.py +++ b/eudemo/eudemo/reactor.py @@ -690,12 +690,12 @@ def save_reactor(reactor, reactor_config, folder_name): cut_angle, ) - if reactor_config.config_for("Neutronics").get("enabled", False): + if reactor_config.config_for("Neutronics", "CSG").get("enabled", False): neutronics_start = time.time() reactor.neutronics = NeutronicsManager( *run_csg_neutronics( - reactor_config.params_for("Neutronics").global_params, - reactor_config.config_for("Neutronics"), + reactor_config.params_for("Neutronics", "CSG").global_params, + reactor_config.config_for("Neutronics", "CSG"), blanket=reactor.blanket, vacuum_vessel=reactor.vacuum_vessel, ivc_shapes=ivc_shapes, @@ -875,6 +875,13 @@ def save_reactor(reactor, reactor_config, folder_name): a_string = f"{reactor_config.global_params.A.value:.2f}".replace(".", "_") folder_name = f"results_v02/A_{a_string}" + run_dagmc_neutronics( + reactor, + reactor_config.params_for("Neutronics", "DAGMC"), + reactor_config.config_for("Neutronics", "DAGMC"), + reference_eq, + ) + Path(folder_name).mkdir(exist_ok=True, parents=True) filename = f"{folder_name}/run_time.json" with open(filename, "w") as f: From a11a9e037c147fd519c56f00b9626ecf2d72a780 Mon Sep 17 00:00:00 2001 From: james <81617086+je-cook@users.noreply.github.com> Date: Mon, 17 Nov 2025 12:15:19 +0000 Subject: [PATCH 06/17] =?UTF-8?q?=F0=9F=9A=A7=20Bad=20rebase=20fixes?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- bluemira/base/constants.py | 9 +- bluemira/codes/openmc/output.py | 2 +- bluemira/codes/openmc/solver.py | 204 +++++++++++------------------- bluemira/codes/openmc/tallying.py | 2 +- eudemo/eudemo/reactor.py | 8 +- 5 files changed, 82 insertions(+), 143 deletions(-) diff --git a/bluemira/base/constants.py b/bluemira/base/constants.py index 6bf328d165..7023d3f103 100644 --- a/bluemira/base/constants.py +++ b/bluemira/base/constants.py @@ -11,7 +11,7 @@ from __future__ import annotations from enum import Enum, auto -from typing import TYPE_CHECKING +from typing import TYPE_CHECKING, TypeVar import numpy as np import numpy.typing as npt @@ -427,11 +427,14 @@ def units_compatible(unit_1: str, unit_2: str) -> bool: return True +ArrayLike = TypeVar("ArrayLike") + + def raw_uc( - value: npt.ArrayLike, + value: ArrayLike, unit_from: str | ureg.Unit, unit_to: str | ureg.Unit, -) -> float | np.ndarray: +) -> ArrayLike: """ Raw unit converter diff --git a/bluemira/codes/openmc/output.py b/bluemira/codes/openmc/output.py index 7a20042897..dcc1651174 100644 --- a/bluemira/codes/openmc/output.py +++ b/bluemira/codes/openmc/output.py @@ -59,7 +59,7 @@ def get_percent_err(row): @dataclass -class OpenMCResult: +class OpenMCCSGResult: """ Class that looks opens up the openmc universe from the statepoint file, so that the dataframes containing the relevant results diff --git a/bluemira/codes/openmc/solver.py b/bluemira/codes/openmc/solver.py index 5a96d9b972..aa891b785f 100644 --- a/bluemira/codes/openmc/solver.py +++ b/bluemira/codes/openmc/solver.py @@ -7,8 +7,8 @@ from __future__ import annotations +from abc import ABC, abstractmethod from collections.abc import Callable -from contextlib import contextmanager from dataclasses import dataclass, fields from enum import auto from operator import attrgetter @@ -30,10 +30,8 @@ CodesTeardown, ) from bluemira.codes.openmc.make_csg import ( - BlanketCellArray, BluemiraNeutronicsCSG, CellStage, - DivertorCellArray, make_cell_arrays, ) from bluemira.codes.openmc.material import MaterialsLibrary @@ -47,6 +45,9 @@ ) from bluemira.codes.openmc.tallying import csg_filter_cells, dagmc_tallys from bluemira.equilibria.equilibrium import Equilibrium +from bluemira.radiation_transport.neutronics.neutronics_axisymmetric import ( + NeutronicsReactor, +) if TYPE_CHECKING: from matproplib.conditions import OperationalConditions @@ -138,22 +139,19 @@ class FigureData: class OpenMCBaseSetup(CodesSetup, ABC): """Setup task for OpenMC solver""" + tally_mats: openmc.Materials + tally_geom: openmc.Geometry | CellStage + def __init__( self, - out_path: str, codes_name: str, cross_section_xml: str, eq: Equilibrium, source, - cell_arrays, - pre_cell_model, materials, ): super().__init__(None, codes_name) - self.out_path = out_path - self.cell_arrays = cell_arrays - self.cells = cell_arrays.cells self.cross_section_xml = cross_section_xml self.eq = eq self.source = source @@ -299,8 +297,7 @@ def __init__( super().__init__(codes_name, cross_section_xml, eq, source, materials) self.cell_arrays = cell_arrays self.pre_cell_model = pre_cell_model - self.materials = materials - self.matlist = attrgetter( + self.mat_list = attrgetter( "outb_sf_mat", "outb_fw_mat", "outb_bz_mat", @@ -311,35 +308,31 @@ def __init__( "tf_coil_mat", ) - @contextmanager - def _base_setup(self, run_mode, *, debug: bool = False): - from openmc.config import config # noqa: PLC0415 + @property + def tally_mats(self) -> list[openmc.Material]: + return self.mat_list(self.materials) - self.files_created = set() - folder = run_mode.name.lower() - cwd = Path(self.out_path, folder) - cwd.mkdir(parents=True, exist_ok=True) + @property + def tally_geom(self) -> CellStage: + return self.cell_arrays - config["cross_sections"] = self.cross_section_xml + def _create_geometry(self): + universe = openmc.Universe(cells=self.cell_arrays.cells) + geometry = openmc.Geometry(universe) + return universe, geometry - self.settings = openmc.Settings( - run_mode=run_mode.value, output={"summary": False} + def plot( + self, + run_mode, + runtime_params, + *_args, + debug: bool = False, + ) -> tuple[openmc.Model, PlotConfig]: + return self._plot( + run_mode, runtime_params, self.pre_cell_model.bounding_box, debug=debug ) - self.universe = openmc.Universe(cells=self.cells) - self.geometry = openmc.Geometry(self.universe) - self.settings.verbosity = 10 if debug else 6 - try: - yield - finally: - for obj, pth in ( - (self.settings, Path(self.out_path, folder, "settings.xml")), - (self.geometry, Path(self.out_path, folder, "geometry.xml")), - (self.materials, Path(self.out_path, folder, "materials.xml")), - ): - obj.export_to_xml(pth) - self.files_created.add(pth) - - def _set_tallies( + + def volume( self, run_mode, runtime_params, @@ -365,74 +358,35 @@ def __init__( materials, dag_model_path: Path, ): - out_path = Path(self.out_path, run_mode.name.lower(), "tallies.xml") - tallies_list = [] - for name, scores, filters in tally_function(material_list, cell_arrays): - tally = openmc.Tally(name=name) - tally.scores = [scores] - tally.filters = filters - tallies_list.append(tally) + super().__init__(codes_name, cross_section_xml, eq, source, materials) + self.dag_model_path = dag_model_path - openmc.Tallies(tallies_list).export_to_xml(out_path) - self.files_created.add(out_path) + def _create_geometry(self): + universe = openmc.DAGMCUniverse( + filename=self.dag_model_path.as_posix(), + auto_geom_ids=True, + ).bounded_universe() + geometry = openmc.Geometry(universe) + return universe, geometry - def run( - self, - run_mode, - runtime_params, - eq, - source_params, - tally_function, - *, - debug: bool = False, - ): - """Run stage for setup openmc""" - with self._base_setup(run_mode, debug=debug): - self.settings.particles = runtime_params.particles - self.settings.source, source_rate, source_t_rate = self.source( - eq, source_params - ) - self.settings.batches = int(runtime_params.batches) - self.settings.photon_transport = runtime_params.photon_transport - self.settings.electron_treatment = runtime_params.electron_treatment - self._set_tallies( - run_mode, - tally_function, - self.cell_arrays, - material_list=self.materials.get_all_materials(), - ) - self._source_rate = source_rate - self._source_triton_rate = source_t_rate - self.files_created.add(f"statepoint.{runtime_params.batches}.h5") - self.files_created.add("tallies.out") + @property + def tally_mats(self) -> list[openmc.Material]: + return self.materials + + @property + def tally_geom(self) -> openmc.Geometry: + return self.geometry def plot( self, run_mode, runtime_params, - _eq, - _source_params, - _tally_function, - *, + *_args, debug: bool = False, - ): - """Plot stage for setup openmc""" - with self._base_setup(run_mode, debug=debug): - z_max, _z_min, r_max, _r_min = self.pre_cell_model.bounding_box - plot_width_0 = r_max * 2.1 - plot_width_1 = z_max * 3.1 - plot = openmc.Plot() - plot.basis = runtime_params.plot_axis - plot.pixels = [ - int(plot_width_0 * runtime_params.plot_pixel_per_metre), - int(plot_width_1 * runtime_params.plot_pixel_per_metre), - ] - plot.width = raw_uc([plot_width_0, plot_width_1], "m", "cm") - plot.show_overlaps = True - - plot_pth = Path(self.out_path, run_mode.name.lower(), "plots.xml") - openmc.Plots([plot]).export_to_xml(plot_pth) - self.files_created.add(plot_pth) + ) -> tuple[openmc.Model, PlotConfig]: + return self._plot( + run_mode, runtime_params, self.universe.bounding_box, debug=debug + ) def volume( self, @@ -487,18 +441,14 @@ def run( geometry_debug=False, restart_file=None, tracks=False, - cwd=cwd, - openmc_exec="openmc", mpi_args=None, event_based=False, - path_input=None, + threads=None, ) - def run(self, run_mode, *, debug: bool = False): - """Run stage for run task""" - self._run(run_mode, debug=debug) - - def plot(self, run_mode, *, debug: bool = False): + def plot( + self, run_mode, model: openmc.Model, config: PlotConfig, *, debug: bool = False + ): """Plot stage for run task""" folder = run_mode.name.lower() cwd = Path(self.out_path, folder) @@ -516,7 +466,9 @@ def plot(self, run_mode, *, debug: bool = False): cwd / "geometry.png", ) - def volume(self, run_mode, *, debug: bool = False): + def volume( + self, run_mode, model: openmc.Model, _config: None, *, debug: bool = False + ): """Stochastic volume stage for run task""" folder = run_mode.name.lower() cwd = Path(self.out_path, folder) @@ -531,7 +483,7 @@ def volume(self, run_mode, *, debug: bool = False): ) -class OpenMCTeardown(CodesTeardown): +class OpenMCCSGTeardown(CodesTeardown): """Teardown task for OpenMC solver""" def __init__( @@ -545,7 +497,6 @@ def __init__( self.out_path = out_path self.cell_arrays = cell_arrays - self.cells = cell_arrays.cells # list[openmc.Cell] self.pre_cell_model = pre_cell_model @staticmethod @@ -576,11 +527,11 @@ def run( delete_files: bool = False, ): """Run stage for Teardown task""" - result = OpenMCResult.from_run( + result = OpenMCCSGResult.from_run( universe, self.cell_arrays, - source_rate, - source_triton_rate, + source_info.rate, + source_info.triton_rate, statepoint_file, ) output_params = NeutronicsOutputParams.from_openmc_csg_result(result) @@ -609,7 +560,7 @@ def volume( cell.id: raw_uc( np.nan if cell.volume is None else cell.volume, "cm^3", "m^3" ) - for cell in self.cells + for cell in self.cell_arrays.cells } @@ -658,7 +609,7 @@ def volume( TALLY_FUNCTION_TYPE = Callable[ - [list[openmc.Material], BlanketCellArray, DivertorCellArray], + [list[openmc.Material], CellStage | openmc.Geometry], tuple[ str, str, @@ -667,26 +618,24 @@ def volume( ] -class OpenMCNeutronicsSolver(CodesSolver): +class OpenMCNeutronicsSolver(CodesSolver, ABC): """OpenMC 2D neutronics solver""" + params: OpenMCNeutronicsSolverParams + setup_cls: type[OpenMCBaseSetup] + teardown_cls: type[CodesTeardown] + name: str = OPENMC_NAME param_cls: type[OpenMCNeutronicsSolverParams] = OpenMCNeutronicsSolverParams - params: OpenMCNeutronicsSolverParams run_mode_cls: type[OpenMCRunModes] = OpenMCRunModes - setup_cls: type[CodesSetup] = OpenMCSetup - run_cls: type[CodesTask] = OpenMCRun - teardown_cls: type[CodesTeardown] = OpenMCTeardown + run_cls: type[OpenMCRun] = OpenMCRun def __init__( self, params: dict | ParameterFrame, build_config: dict, - neutronics_pre_cell_model, eq: Equilibrium, source: NeutronSourceCreator, - op_cond: OperationalConditions, - tally_function: TALLY_FUNCTION_TYPE | None = None, ): self.params = make_parameter_frame(params, self.param_cls) self.build_config = build_config @@ -696,17 +645,6 @@ def __init__( self.eq = eq self.source = source - self.pre_cell_model = neutronics_pre_cell_model - self.materials = MaterialsLibrary.from_neutronics_materials( - self.pre_cell_model.material_library, op_cond - ) - - self.cell_arrays = make_cell_arrays( - self.pre_cell_model, BluemiraNeutronicsCSG(), self.materials, control_id=True - ) - - self.tally_function = filter_cells if tally_function is None else tally_function - @property def source(self) -> NeutronSourceCreator: """Source term for OpenMC""" @@ -754,7 +692,7 @@ def _single_run( ) -> OpenMCCSGResult | dict[int, float]: result = None if setup := self._get_execution_method(self._setup, run_mode): - result = setup( + model, config = setup( run_mode, runtime_params, self.eq, @@ -763,7 +701,7 @@ def _single_run( debug=debug, ) if run := self._get_execution_method(self._run, run_mode): - result = run(run_mode, debug=debug) + result = run(run_mode, model, config, debug=debug) if teardown := self._get_execution_method(self._teardown, run_mode): result = teardown(self._setup.universe, config, result) return result diff --git a/bluemira/codes/openmc/tallying.py b/bluemira/codes/openmc/tallying.py index 78d902abbe..08330f1126 100644 --- a/bluemira/codes/openmc/tallying.py +++ b/bluemira/codes/openmc/tallying.py @@ -12,7 +12,7 @@ from bluemira.codes.openmc.make_csg import CellStage -def filter_cells( +def csg_filter_cells( material_list, csg_model: CellStage, ): diff --git a/eudemo/eudemo/reactor.py b/eudemo/eudemo/reactor.py index e2c7aebace..ea811cdf30 100644 --- a/eudemo/eudemo/reactor.py +++ b/eudemo/eudemo/reactor.py @@ -46,7 +46,6 @@ from bluemira.builders.radiation_shield import RadiationShieldBuilder from bluemira.builders.thermal_shield import CryostatTSBuilder, VVTSBuilder from bluemira.equilibria.equilibrium import Equilibrium -from bluemira.equilibria.flux_surfaces import ClosedFluxSurface from bluemira.equilibria.profiles import Profile from bluemira.equilibria.run import Snapshot from bluemira.geometry.coordinates import Coordinates @@ -99,7 +98,6 @@ from eudemo.neutronics.run import run_csg_neutronics, run_dagmc_neutronics from eudemo.params import EUDEMOReactorParams from eudemo.pf_coils import PFCoil, PFCoilsDesigner, build_pf_coils_component -from eudemo.power_cycle import SteadyStatePowerCycleSolver from eudemo.radial_build import radial_build from eudemo.tf_coils import TFCoil, TFCoilBuilder, TFCoilDesigner from eudemo.vacuum_vessel import VacuumVessel, VacuumVesselBuilder @@ -706,7 +704,7 @@ def save_reactor(reactor, reactor_config, folder_name): neutronics_end = time.time() run_time_track["CSG neutronics"] = neutronics_end - neutronics_start - if reactor_config.config_for("Neutronics")["show_data"]: + if reactor_config.config_for("Neutronics", "CSG")["show_data"]: reactor.neutronics.plot() bluemira_print_clean(f"{reactor.neutronics}") else: @@ -862,8 +860,8 @@ def save_reactor(reactor, reactor_config, folder_name): run_dagmc_neutronics( reactor, - reactor_config.params_for("CAD_Neutronics"), - reactor_config.config_for("CAD_Neutronics"), + reactor_config.params_for("Neutronics", "DAGMC"), + reactor_config.config_for("Neutronics", "DAGMC"), reference_eq, ) From cd613d363025a00ffec768cf7ff45abdf89baf1d Mon Sep 17 00:00:00 2001 From: james <81617086+je-cook@users.noreply.github.com> Date: Mon, 17 Nov 2025 13:48:46 +0000 Subject: [PATCH 07/17] =?UTF-8?q?=F0=9F=9A=A7=20Initial=20openmcdagmc=20re?= =?UTF-8?q?sults?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- bluemira/codes/openmc/output.py | 21 ++++++++++++++++-- bluemira/codes/openmc/solver.py | 39 ++++----------------------------- eudemo/eudemo/neutronics/run.py | 1 + 3 files changed, 24 insertions(+), 37 deletions(-) diff --git a/bluemira/codes/openmc/output.py b/bluemira/codes/openmc/output.py index dcc1651174..384c17a1b3 100644 --- a/bluemira/codes/openmc/output.py +++ b/bluemira/codes/openmc/output.py @@ -104,7 +104,7 @@ def from_run( cell_arrays: CellStage, src_rate: float, src_triton_rate: float, - statepoint_file: str = "", + statepoint_file: Path, ): """Create results class from run statepoint""" # Create cell and material name dictionaries to allow easy mapping to dataframe @@ -125,7 +125,7 @@ def from_run( } # Loads up the output file from the simulation - statepoint = openmc.StatePoint(statepoint_file) + statepoint = openmc.StatePoint(statepoint_file.as_posix()) tbr, tbr_err = cls._load_tbr(statepoint, src_rate, src_triton_rate) blanket_power, blanket_power_err = cls._load_filter_power_err( statepoint, src_rate, "breeding blanket power" @@ -520,6 +520,19 @@ def _tabulate( ) +@dataclass +class OpenMCDAGMCResult: + @classmethod + def from_file( + cls, + universe: openmc.Universe, + src_rate: float, + src_triton_rate: float, + statepoint_file: Path, + ): + statepoint = openmc.StatePoint(statepoint_file.as_posix()) + + @dataclass class NeutronicsOutputParams(ParameterFrame): """ @@ -539,6 +552,10 @@ class NeutronicsOutputParams(ParameterFrame): peak_vv_iron_dpa_rate: Parameter[float] peak_div_cu_dpa_rate: Parameter[float] + @classmethod + def from_openmc_dagmc_result(cls, result: OpenMCDAGMCResult): + return cls() + @classmethod def from_openmc_csg_result(cls, result: OpenMCCSGResult): """ diff --git a/bluemira/codes/openmc/solver.py b/bluemira/codes/openmc/solver.py index aa891b785f..0648591695 100644 --- a/bluemira/codes/openmc/solver.py +++ b/bluemira/codes/openmc/solver.py @@ -19,7 +19,6 @@ import openmc from bluemira.base.constants import raw_uc -from bluemira.base.look_and_feel import bluemira_debug from bluemira.base.parameter_frame import ParameterFrame, make_parameter_frame from bluemira.base.tools import _timing from bluemira.codes.interface import ( @@ -38,6 +37,7 @@ from bluemira.codes.openmc.output import ( NeutronicsOutputParams, OpenMCCSGResult, + OpenMCDAGMCResult, ) from bluemira.codes.openmc.params import ( OpenMCNeutronicsSolverParams, @@ -499,32 +499,11 @@ def __init__( self.cell_arrays = cell_arrays self.pre_cell_model = pre_cell_model - @staticmethod - def delete_files(files_created): - """Remove files generated during the run (mainly .xml files.)""" - removed_files, failed_to_remove_files = [], [] - for file_name in files_created: - if (f := file_name).exists(): - f.unlink() - removed_files.append(file_name) - else: - failed_to_remove_files.append(file_name) - - if removed_files: - bluemira_debug(f"Removed files {removed_files}") - if failed_to_remove_files: - bluemira_debug( - f"Attempted to remove files {failed_to_remove_files} but " - "they don't exists." - ) - def run( self, universe, source_info: SourceInfo, statepoint_file, - *, - delete_files: bool = False, ): """Run stage for Teardown task""" result = OpenMCCSGResult.from_run( @@ -532,12 +511,10 @@ def run( self.cell_arrays, source_info.rate, source_info.triton_rate, - statepoint_file, + Path(statepoint_file), ) output_params = NeutronicsOutputParams.from_openmc_csg_result(result) - if delete_files: - self.delete_files(files_created) return result, output_params def plot(self, _universe, _source_info, fig: FigureData, **kwargs): @@ -550,12 +527,8 @@ def volume( _universe, _source_params, _statepoint_file, - *, - delete_files: bool = False, ) -> dict[int, float]: """Stochastic volume stage for teardown task""" - if delete_files: - self.delete_files(files_created) return { cell.id: raw_uc( np.nan if cell.volume is None else cell.volume, "cm^3", "m^3" @@ -575,12 +548,10 @@ def run( universe, source_info: SourceInfo, statepoint_file, - *, - delete_files: bool = False, ): """Run stage for Teardown task""" - result = OpenMCDAGResult.from_run( - universe, source_info.rate, source_info.triton_rate, statepoint_file + result = OpenMCDAGMCResult.from_run( + universe, source_info.rate, source_info.triton_rate, Path(statepoint_file) ) output_params = NeutronicsOutputParams.from_openmc_dag_result(result) @@ -602,8 +573,6 @@ def volume( _universe, _source_params, _statepoint_file, - *, - delete_files: bool = False, ) -> dict[int, float]: """Stochastic volume stage for teardown task""" diff --git a/eudemo/eudemo/neutronics/run.py b/eudemo/eudemo/neutronics/run.py index 6538344877..22bc119a97 100644 --- a/eudemo/eudemo/neutronics/run.py +++ b/eudemo/eudemo/neutronics/run.py @@ -144,6 +144,7 @@ def export_dagmc_model(reactor, build_config): "without_components": [ reactor.plasma, reactor.blanket, + reactor.coil_structures, ], "group_by_materials": True, }, From 739c2f84d9c2b36349fd502be32a868ab46ffffb Mon Sep 17 00:00:00 2001 From: james <81617086+je-cook@users.noreply.github.com> Date: Mon, 17 Nov 2025 13:54:17 +0000 Subject: [PATCH 08/17] =?UTF-8?q?=F0=9F=90=9B=20More=20rebase=20misses?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- eudemo/eudemo/reactor.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/eudemo/eudemo/reactor.py b/eudemo/eudemo/reactor.py index ea811cdf30..422b35d7d4 100644 --- a/eudemo/eudemo/reactor.py +++ b/eudemo/eudemo/reactor.py @@ -46,6 +46,7 @@ from bluemira.builders.radiation_shield import RadiationShieldBuilder from bluemira.builders.thermal_shield import CryostatTSBuilder, VVTSBuilder from bluemira.equilibria.equilibrium import Equilibrium +from bluemira.equilibria.flux_surfaces import ClosedFluxSurface from bluemira.equilibria.profiles import Profile from bluemira.equilibria.run import Snapshot from bluemira.geometry.coordinates import Coordinates @@ -98,6 +99,7 @@ from eudemo.neutronics.run import run_csg_neutronics, run_dagmc_neutronics from eudemo.params import EUDEMOReactorParams from eudemo.pf_coils import PFCoil, PFCoilsDesigner, build_pf_coils_component +from eudemo.power_cycle import SteadyStatePowerCycleSolver from eudemo.radial_build import radial_build from eudemo.tf_coils import TFCoil, TFCoilBuilder, TFCoilDesigner from eudemo.vacuum_vessel import VacuumVessel, VacuumVesselBuilder @@ -865,6 +867,14 @@ def save_reactor(reactor, reactor_config, folder_name): reference_eq, ) + sspc_solver = SteadyStatePowerCycleSolver(reactor_config.global_params) + sspc_result = sspc_solver.execute() + reactor_config.global_params.P_el_net.set_value( + sspc_result["P_el_net"], "BLUEMIRA" + ) + + lcfs = ClosedFluxSurface(reference_eq.get_LCFS()) + reactor_config.global_params.V_p.set_value(lcfs.volume, "BLUEMIRA") end = time.time() From f480ff431a11cc3951679ce4768650240aa95cc3 Mon Sep 17 00:00:00 2001 From: james <81617086+je-cook@users.noreply.github.com> Date: Wed, 19 Nov 2025 09:17:51 +0000 Subject: [PATCH 09/17] =?UTF-8?q?=E2=AC=86=EF=B8=8F=20Enforce=20recent=20t?= =?UTF-8?q?okamak=20source?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 30bcd271e6..4ce418cc50 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -91,7 +91,7 @@ process = ["process @ git+https://github.com/ukaea/process@v3.2.2"] openmc = [ "openmc>=0.14.0", "openmc_data==2.3.2", - "tokamak-neutron-source >=0.1.1", + "tokamak-neutron-source >=0.1.4", ] polyscope = ["polyscope"] radiation = ["cherab"] From 548206b919591992357746b9f426409815b75cbe Mon Sep 17 00:00:00 2001 From: james <81617086+je-cook@users.noreply.github.com> Date: Wed, 19 Nov 2025 13:55:34 +0000 Subject: [PATCH 10/17] =?UTF-8?q?=F0=9F=9A=A7=20Tmp=20copy=20from=20exampl?= =?UTF-8?q?e?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- bluemira/codes/openmc/output.py | 71 ++++++++++++++++++++++++++++++++- 1 file changed, 70 insertions(+), 1 deletion(-) diff --git a/bluemira/codes/openmc/output.py b/bluemira/codes/openmc/output.py index 384c17a1b3..dcff7e0702 100644 --- a/bluemira/codes/openmc/output.py +++ b/bluemira/codes/openmc/output.py @@ -16,7 +16,7 @@ from tabulate import tabulate from bluemira.base.constants import raw_uc -from bluemira.base.look_and_feel import bluemira_debug +from bluemira.base.look_and_feel import bluemira_debug, bluemira_print from bluemira.base.parameter_frame._frame import ParameterFrame from bluemira.base.parameter_frame._parameter import Parameter from bluemira.codes.openmc.make_csg import CellStage @@ -532,6 +532,75 @@ def from_file( ): statepoint = openmc.StatePoint(statepoint_file.as_posix()) + tbr_cell_tally = statepoint.get_tally(name="tbr") + tbr_mesh_tally = statepoint.get_tally(name="tbr_on_mesh") + heating_cell_tally = statepoint.get_tally(name="heating") + heating_mesh_tally = statepoint.get_tally(name="heating_on_mesh") + flux_mesh_tally = statepoint.get_tally(name="flux_on_mesh") + + bluemira_print(f"The reactor has a TBR of {tbr_cell_tally.mean.sum()}") + bluemira_print( + f"Standard deviation on the TBR is {tbr_cell_tally.std_dev.sum()}" + ) + + bluemira_print( + f"The heating of {heating_cell_tally.mean.sum() / 1e6} MeV " + "per source particle is deposited" + ) + bluemira_print( + f"Standard deviation on the heating tally is {heating_cell_tally.std_dev.sum()}" + ) + + mesh = tbr_mesh_tally.find_filter(openmc.MeshFilter).mesh + mesh.write_data_to_vtk( + filename="tbr_mesh_mean.vtk", + datasets={"mean": tbr_mesh_tally.mean}, + ) + + model_w = universe.bounding_box.width + + heating_mesh_mean = heating_mesh_tally.mean.reshape(100, 100, 100) + flux_mesh_mean = flux_mesh_tally.mean.reshape(100, 100, 100) + scaling = tuple( + round(t / c) for c, t in zip(heating_mesh_mean.shape, model_w, strict=False) + ) + # If you want to view only the +ve x half of the tally, uncomment the next line + # heating_mesh_mean[:, :, heating_mesh_mean.shape[2] // 2 : -1] = 0 + numpy_to_vtk(heating_mesh_mean, "heating_mesh_mean", scaling) + numpy_to_vtk(flux_mesh_mean, "flux_mesh_mean", scaling) + + +import vtk +from vtkmodules.util import numpy_support + + +def numpy_to_vtk(data, output_name, scaling=(1, 1, 1)): + """Convert a numpy array to a VTK image data file.""" + data_type = vtk.VTK_FLOAT + shape = data.shape + + flat_data_array = data.flatten() + vtk_data = numpy_support.numpy_to_vtk( + num_array=flat_data_array, deep=True, array_type=data_type + ) + vtk_data.SetName(output_name) + + half_x = int(0.5 * scaling[0] * (shape[0] - 1)) + half_y = int(0.5 * scaling[1] * (shape[1] - 1)) + half_z = int(0.5 * scaling[2] * (shape[2] - 1)) + + img = vtk.vtkImageData() + img.GetPointData().SetScalars(vtk_data) + img.SetSpacing(scaling[0], scaling[1], scaling[2]) + img.SetDimensions(shape[0], shape[1], shape[2]) + img.SetOrigin(-half_x, -half_y, -half_z) + + # Save the VTK file + writer = vtk.vtkXMLImageDataWriter() + writer.SetFileName(f"{output_name}.vti") + writer.SetInputData(img) + writer.Write() + @dataclass class NeutronicsOutputParams(ParameterFrame): From 251961bcea257b8d44df0fba372adc1b8a309a96 Mon Sep 17 00:00:00 2001 From: james <81617086+je-cook@users.noreply.github.com> Date: Wed, 19 Nov 2025 16:02:38 +0000 Subject: [PATCH 11/17] =?UTF-8?q?=E2=9C=A8=20Initial=20'working'=20version?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- bluemira/codes/openmc/output.py | 225 ++++++++++++++++++++------------ bluemira/codes/openmc/solver.py | 3 +- eudemo/eudemo/reactor.py | 10 +- 3 files changed, 142 insertions(+), 96 deletions(-) diff --git a/bluemira/codes/openmc/output.py b/bluemira/codes/openmc/output.py index dcff7e0702..3d8bb980a5 100644 --- a/bluemira/codes/openmc/output.py +++ b/bluemira/codes/openmc/output.py @@ -16,7 +16,7 @@ from tabulate import tabulate from bluemira.base.constants import raw_uc -from bluemira.base.look_and_feel import bluemira_debug, bluemira_print +from bluemira.base.look_and_feel import bluemira_debug from bluemira.base.parameter_frame._frame import ParameterFrame from bluemira.base.parameter_frame._parameter import Parameter from bluemira.codes.openmc.make_csg import CellStage @@ -58,8 +58,67 @@ def get_percent_err(row): return row["std. dev."] / row["mean"] * 100.0 +class OpenMCResultBase: + """Base class for openmc results""" + + @staticmethod + def _load_dataframe_from_statepoint( + statepoint: openmc.StatePoint, tally_name: str + ) -> pd.DataFrame: + return statepoint.get_tally(name=tally_name).get_pandas_dataframe() + + @staticmethod + def _convert_dict_contents(dataset: dict[str, dict[int, list[str | float]]]): + for k, v in dataset.items(): + vals = list(v.values()) if isinstance(v, dict) else v + dataset[k] = vals if isinstance(vals[0], str) else np.array(vals) + return dataset + + @staticmethod + def dt_neuton_power(src_triton_rate): + """Calculate dt neutron power""" + return 0.8 * E_DT_fusion() * src_triton_rate + + @staticmethod + def energy_multiplication(dt_neutron_power, total_power): + """Calculate energy multiplication""" + return total_power / dt_neutron_power + + @staticmethod + def multiplication_power(energy_multiplication, dt_neutron_power): + """Calculate multiplication power""" + return (energy_multiplication - 1.0) * dt_neutron_power + + @classmethod + def _load_filter_power_err( + cls, statepoint, src_rate: float, filter_name: str + ) -> tuple[float, float]: + """ + Power is initially loaded as eV/source particle. To convert to Watt, we need the + source particle rate. + + Parameters + ---------- + filter_name: + the literal name that was used in tallying.py to refer to this tally. + src_rate: + source particle rate. + + Returns + ------- + power: + The total power [W]. + errors: + The absolute error on the total power [W]. RMS of errors from each cell. + """ + df = cls._load_dataframe_from_statepoint(statepoint, filter_name) + powers = raw_uc(df["mean"].to_numpy() * src_rate, "eV/s", "W") + errors = raw_uc(df["std. dev."].to_numpy() * src_rate, "eV/s", "W") + return powers.sum(), np.sqrt((errors**2).sum()) + + @dataclass -class OpenMCCSGResult: +class OpenMCCSGResult(OpenMCResultBase): """ Class that looks opens up the openmc universe from the statepoint file, so that the dataframes containing the relevant results @@ -141,10 +200,10 @@ def from_run( statepoint, src_rate, "total power" ) - dt_neuton_power = 0.8 * E_DT_fusion() * src_triton_rate - e_mult = total_power / dt_neuton_power - e_mult_err = total_power_err / dt_neuton_power - mult_power = (e_mult - 1.0) * dt_neuton_power + dt_n_power = cls.dt_neuton_power(src_triton_rate) + e_mult = cls.energy_multiplication(dt_n_power, total_power) + e_mult_err = cls.energy_multiplication(dt_n_power, total_power_err) + mult_power = cls.multiplication_power(e_mult, dt_n_power) all_fluxes = cls._load_fluxes(statepoint, cell_names, cell_vols, src_rate) damage = cls._load_damage( statepoint, cell_names, cell_vols, cell_arrays, src_rate @@ -216,19 +275,6 @@ def _load_volume_calculation_from_file( return vol_results, cell_volumes - @staticmethod - def _load_dataframe_from_statepoint( - statepoint: openmc.StatePoint, tally_name: str - ) -> pd.DataFrame: - return statepoint.get_tally(name=tally_name).get_pandas_dataframe() - - @staticmethod - def _convert_dict_contents(dataset: dict[str, dict[int, list[str | float]]]): - for k, v in dataset.items(): - vals = list(v.values()) if isinstance(v, dict) else v - dataset[k] = vals if isinstance(vals[0], str) else np.array(vals) - return dataset - @classmethod def _load_tbr(cls, statepoint, source_rate: float, source_triton_rate: float): """ @@ -247,33 +293,6 @@ def _load_tbr(cls, statepoint, source_rate: float, source_triton_rate: float): # Single tally, so std dev scales linearly return scale * tbr_df["mean"].iloc[0], scale * tbr_df["std. dev."].iloc[0] - @classmethod - def _load_filter_power_err( - cls, statepoint, src_rate: float, filter_name: str - ) -> tuple[float, float]: - """ - Power is initially loaded as eV/source particle. To convert to Watt, we need the - source particle rate. - - Parameters - ---------- - filter_name: - the literal name that was used in tallying.py to refer to this tally. - src_rate: - source particle rate. - - Returns - ------- - power: - The total power [W]. - errors: - The absolute error on the total power [W]. RMS of errors from each cell. - """ - df = cls._load_dataframe_from_statepoint(statepoint, filter_name) - powers = raw_uc(df["mean"].to_numpy() * src_rate, "eV/s", "W") - errors = raw_uc(df["std. dev."].to_numpy() * src_rate, "eV/s", "W") - return powers.sum(), np.sqrt((errors**2).sum()) - @classmethod def _load_heating(cls, statepoint, mat_names, src_rate): """Load the heating (sorted by material) dataframe""" @@ -521,36 +540,43 @@ def _tabulate( @dataclass -class OpenMCDAGMCResult: +class OpenMCDAGMCResult(OpenMCResultBase): + """ + Open up the openmc universe from the statepoint file, + so that the dataframes containing the relevant results + can be generated and reformatted by its methods. + """ + + tbr: float + tbr_err: float + e_mult: float + e_mult_err: float + mult_power: float + src_rate: float + + statepoint: openmc.StatePoint + statepoint_file: Path + statepoint_file: openmc.StatePoint + @classmethod - def from_file( + def from_run( cls, universe: openmc.Universe, src_rate: float, src_triton_rate: float, statepoint_file: Path, ): + """Create results class from run statepoint""" statepoint = openmc.StatePoint(statepoint_file.as_posix()) tbr_cell_tally = statepoint.get_tally(name="tbr") tbr_mesh_tally = statepoint.get_tally(name="tbr_on_mesh") - heating_cell_tally = statepoint.get_tally(name="heating") heating_mesh_tally = statepoint.get_tally(name="heating_on_mesh") flux_mesh_tally = statepoint.get_tally(name="flux_on_mesh") - bluemira_print(f"The reactor has a TBR of {tbr_cell_tally.mean.sum()}") - bluemira_print( - f"Standard deviation on the TBR is {tbr_cell_tally.std_dev.sum()}" - ) - - bluemira_print( - f"The heating of {heating_cell_tally.mean.sum() / 1e6} MeV " - "per source particle is deposited" - ) - bluemira_print( - f"Standard deviation on the heating tally is {heating_cell_tally.std_dev.sum()}" + total_power, total_power_err = cls._load_filter_power_err( + statepoint, src_rate, "heating" ) - mesh = tbr_mesh_tally.find_filter(openmc.MeshFilter).mesh mesh.write_data_to_vtk( filename="tbr_mesh_mean.vtk", @@ -569,9 +595,25 @@ def from_file( numpy_to_vtk(heating_mesh_mean, "heating_mesh_mean", scaling) numpy_to_vtk(flux_mesh_mean, "flux_mesh_mean", scaling) + dt_n_power = cls.dt_neuton_power(src_triton_rate) + e_mult = cls.energy_multiplication(dt_n_power, total_power) + e_mult_err = cls.energy_multiplication(dt_n_power, total_power_err) + mult_power = cls.multiplication_power(e_mult, dt_n_power) + return cls( + tbr=tbr_cell_tally.mean.sum(), + tbr_err=tbr_cell_tally.std_dev.sum(), + e_mult=e_mult, + e_mult_err=e_mult_err, + src_rate=src_rate, + mult_power=mult_power, + statepoint_file=statepoint_file, + statepoint=statepoint, + ) + -import vtk -from vtkmodules.util import numpy_support +# TODO @je-cook: move this function somewhere sensible # noqa: TD003 +import vtk # noqa: E402 +from vtkmodules.util import numpy_support # noqa: E402 def numpy_to_vtk(data, output_name, scaling=(1, 1, 1)): @@ -622,8 +664,29 @@ class NeutronicsOutputParams(ParameterFrame): peak_div_cu_dpa_rate: Parameter[float] @classmethod - def from_openmc_dagmc_result(cls, result: OpenMCDAGMCResult): - return cls() + def from_openmc_dag_result(cls, result: OpenMCDAGMCResult): + """ + Produce output parameters from an OpenMC DAGMC result + """ + source = "OpenMC DAGMC" + + return cls( + Parameter("e_mult", result.e_mult, unit="", source=source), + Parameter("TBR", result.tbr, unit="", source=source), + Parameter("P_n_blanket", None, unit="W", source=source), + Parameter("P_n_divertor", None, unit="W", source=source), + Parameter("P_n_vessel", None, unit="W", source=source), + # No auxiliaries (e.g. HCD, port plugs modelled here) + Parameter("P_n_aux", 0.0, unit="W", source=source), + Parameter("P_n_e_mult", result.mult_power, unit="W", source=source), + Parameter( # can't get this without coupling to D1S/R2S/involving fispact + "P_n_decay", 0.0, unit="W", source=source + ), + Parameter("peak_eurofer_dpa_rate", None, unit="dpa/fpy", source=source), + Parameter("peak_bb_iron_dpa_rate", None, unit="dpa/fpy", source=source), + Parameter("peak_vv_iron_dpa_rate", None, unit="dpa/fpy", source=source), + Parameter("peak_div_cu_dpa_rate", None, unit="dpa/fpy", source=source), + ) @classmethod def from_openmc_csg_result(cls, result: OpenMCCSGResult): @@ -632,10 +695,10 @@ def from_openmc_csg_result(cls, result: OpenMCCSGResult): """ source = "OpenMC CSG" - peak_eurofer_dpa_rate = result.damage["eurofer damage"]["dpa/fpy"].max() - peak_bb_iron_dpa_rate = result.damage["blanket damage"]["dpa/fpy"].max() - peak_vv_iron_dpa_rate = result.damage["VV damage"]["dpa/fpy"].max() - peak_div_cu_dpa_rate = result.damage["divertor damage"]["dpa/fpy"].max() + peak_eurofer_dpa = result.damage["eurofer damage"]["dpa/fpy"].max() + peak_bb_iron_dpa = result.damage["blanket damage"]["dpa/fpy"].max() + peak_vv_iron_dpa = result.damage["VV damage"]["dpa/fpy"].max() + peak_div_cu_dpa = result.damage["divertor damage"]["dpa/fpy"].max() return cls( Parameter("e_mult", result.e_mult, unit="", source=source), Parameter("TBR", result.tbr, unit="", source=source), @@ -645,32 +708,20 @@ def from_openmc_csg_result(cls, result: OpenMCCSGResult): # No auxiliaries (e.g. HCD, port plugs modelled here) Parameter("P_n_aux", 0.0, unit="W", source=source), Parameter("P_n_e_mult", result.mult_power, unit="W", source=source), - Parameter( + Parameter( # can't get this without coupling to D1S/R2S/involving fispact "P_n_decay", 0.0, unit="W", source=source - ), # can't get this without coupling to D1S/R2S/involving fispact + ), Parameter( - "peak_eurofer_dpa_rate", - peak_eurofer_dpa_rate, - unit="dpa/fpy", - source=source, + "peak_eurofer_dpa_rate", peak_eurofer_dpa, unit="dpa/fpy", source=source ), Parameter( - "peak_bb_iron_dpa_rate", - peak_bb_iron_dpa_rate, - unit="dpa/fpy", - source=source, + "peak_bb_iron_dpa_rate", peak_bb_iron_dpa, unit="dpa/fpy", source=source ), Parameter( - "peak_vv_iron_dpa_rate", - peak_vv_iron_dpa_rate, - unit="dpa/fpy", - source=source, + "peak_vv_iron_dpa_rate", peak_vv_iron_dpa, unit="dpa/fpy", source=source ), Parameter( - "peak_div_cu_dpa_rate", - peak_div_cu_dpa_rate, - unit="dpa/fpy", - source=source, + "peak_div_cu_dpa_rate", peak_div_cu_dpa, unit="dpa/fpy", source=source ), ) diff --git a/bluemira/codes/openmc/solver.py b/bluemira/codes/openmc/solver.py index 0648591695..9cbf54172c 100644 --- a/bluemira/codes/openmc/solver.py +++ b/bluemira/codes/openmc/solver.py @@ -566,7 +566,7 @@ def plot( ): """Plot stage for Teardown task""" fig.axis.get_figure().save_fig(fig.path) - return axis + return fig.axis def volume( self, @@ -575,6 +575,7 @@ def volume( _statepoint_file, ) -> dict[int, float]: """Stochastic volume stage for teardown task""" + raise NotImplementedError TALLY_FUNCTION_TYPE = Callable[ diff --git a/eudemo/eudemo/reactor.py b/eudemo/eudemo/reactor.py index 422b35d7d4..315fc0830e 100644 --- a/eudemo/eudemo/reactor.py +++ b/eudemo/eudemo/reactor.py @@ -860,7 +860,8 @@ def save_reactor(reactor, reactor_config, folder_name): n_TF=reactor_config.global_params.n_TF.value, ) - run_dagmc_neutronics( + # TODO @je-cook: Put the results somewhere + dagmc_out = run_dagmc_neutronics( reactor, reactor_config.params_for("Neutronics", "DAGMC"), reactor_config.config_for("Neutronics", "DAGMC"), @@ -883,13 +884,6 @@ def save_reactor(reactor, reactor_config, folder_name): a_string = f"{reactor_config.global_params.A.value:.2f}".replace(".", "_") folder_name = f"results_v02/A_{a_string}" - run_dagmc_neutronics( - reactor, - reactor_config.params_for("Neutronics", "DAGMC"), - reactor_config.config_for("Neutronics", "DAGMC"), - reference_eq, - ) - Path(folder_name).mkdir(exist_ok=True, parents=True) filename = f"{folder_name}/run_time.json" with open(filename, "w") as f: From 1dec4f62b0ee37bd19c9980ec74bcca43b1524be Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Wed, 19 Nov 2025 16:04:35 +0000 Subject: [PATCH 12/17] Updated dependencies in requirements files [skip dependency] --- requirements/uv/all.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements/uv/all.txt b/requirements/uv/all.txt index 06bbe6c34b..de450c69a7 100644 --- a/requirements/uv/all.txt +++ b/requirements/uv/all.txt @@ -20,7 +20,7 @@ asteval==1.0.6 # via # bluemira (pyproject.toml) # neutronics-material-maker -astroid>=3.3.11 +astroid==3.3.11 # via sphinx-autoapi asttokens==3.0.0 # via From fa499e90545de945c37df4744b2ec9f676e5d4e8 Mon Sep 17 00:00:00 2001 From: matti Date: Thu, 20 Nov 2025 14:19:50 +0100 Subject: [PATCH 13/17] add converter args --- eudemo/config/build_config.json | 20 ++++++++++++++++++-- eudemo/eudemo/neutronics/run.py | 3 ++- 2 files changed, 20 insertions(+), 3 deletions(-) diff --git a/eudemo/config/build_config.json b/eudemo/config/build_config.json index a1ebd1162d..20b2ca7765 100644 --- a/eudemo/config/build_config.json +++ b/eudemo/config/build_config.json @@ -134,11 +134,27 @@ "show_data": false }, "DAGMC": { - "export_dagmc_model": false, + "export_dagmc_model": true, "dagmc_export_dir": "$path_expand:./dagmc_export", "cross_section_xml": "$path_expand:./cross_section_data/cross_sections.xml", "particles": 10000, - "batches": 2 + "batches": 2, + "converter_config": { + "converter_type": "fast_ctd", + "imprint_geometry": true, + "imprint_per_compound": true, + "minimum_include_volume": 1.0, + "fix_step_to_brep_geometry": false, + "merge_dist_tolerance": 0.001, + "lin_deflection_tol": 0.001, + "lin_deflection_is_absolute": false, + "angular_deflection_tol": 0.5, + "run_make_watertight": true, + "save_vtk_model": true, + "enable_ext_debug_logging": false, + "use_cached_files": false, + "clean_up_cached": true + } } }, "TF coils": { diff --git a/eudemo/eudemo/neutronics/run.py b/eudemo/eudemo/neutronics/run.py index 22bc119a97..875d5353fc 100644 --- a/eudemo/eudemo/neutronics/run.py +++ b/eudemo/eudemo/neutronics/run.py @@ -143,11 +143,12 @@ def export_dagmc_model(reactor, build_config): construction_params={ "without_components": [ reactor.plasma, - reactor.blanket, + # reactor.blanket, reactor.coil_structures, ], "group_by_materials": True, }, + **build_config.get("converter_config", {}), ) From 6118f6efadde7d071b0f6ced8d13b1660a95c3af Mon Sep 17 00:00:00 2001 From: matti Date: Thu, 20 Nov 2025 14:45:04 +0100 Subject: [PATCH 14/17] use config correctly --- eudemo/config/build_config.json | 2 +- eudemo/eudemo/neutronics/run.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/eudemo/config/build_config.json b/eudemo/config/build_config.json index 20b2ca7765..5bc670fb9d 100644 --- a/eudemo/config/build_config.json +++ b/eudemo/config/build_config.json @@ -142,7 +142,7 @@ "converter_config": { "converter_type": "fast_ctd", "imprint_geometry": true, - "imprint_per_compound": true, + "imprint_per_compound": false, "minimum_include_volume": 1.0, "fix_step_to_brep_geometry": false, "merge_dist_tolerance": 0.001, diff --git a/eudemo/eudemo/neutronics/run.py b/eudemo/eudemo/neutronics/run.py index 875d5353fc..5135575cc4 100644 --- a/eudemo/eudemo/neutronics/run.py +++ b/eudemo/eudemo/neutronics/run.py @@ -148,7 +148,7 @@ def export_dagmc_model(reactor, build_config): ], "group_by_materials": True, }, - **build_config.get("converter_config", {}), + converter_config=build_config.get("converter_config", {}), ) From 1d308934f67839d4baba6d8c6eb6ebf3738b5272 Mon Sep 17 00:00:00 2001 From: james <81617086+je-cook@users.noreply.github.com> Date: Fri, 21 Nov 2025 11:50:48 +0000 Subject: [PATCH 15/17] =?UTF-8?q?=F0=9F=93=8C=20Add=20openmc=20deps=20to?= =?UTF-8?q?=20dagmc=20extras?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- pyproject.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/pyproject.toml b/pyproject.toml index 4ce418cc50..6460c9808f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -98,6 +98,7 @@ radiation = ["cherab"] dagmc = [ "cgal>=6.0.1", "fast_ctd @ git+https://github.com/Fusion-Power-Plant-Framework/fast_ctd@0.1.1", + "bluemira[openmc]" ] [build-system] From a1f890cb2984f249ad355e1fdd3abc92cf724596 Mon Sep 17 00:00:00 2001 From: james <81617086+je-cook@users.noreply.github.com> Date: Fri, 21 Nov 2025 12:04:00 +0000 Subject: [PATCH 16/17] =?UTF-8?q?=F0=9F=90=9B=20Warn=20instead=20of=20cras?= =?UTF-8?q?h=20without=20eurofusion=20materials?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- bluemira/materials/cache.py | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/bluemira/materials/cache.py b/bluemira/materials/cache.py index c9b4422e00..de61961335 100644 --- a/bluemira/materials/cache.py +++ b/bluemira/materials/cache.py @@ -16,6 +16,7 @@ from matproplib.library.fluids import Void from matproplib.material import Material +from bluemira.base.look_and_feel import bluemira_warn from bluemira.materials.error import MaterialsError from bluemira.utilities.tools import get_module @@ -63,13 +64,18 @@ def __getattr__(self, value: str): try: super().__getattribute__(value) except AttributeError: - if value in self._material_package: - return self._material_package[value] + if value in self._material_packages: + return self._material_packages[value] raise def load_from_package(self, package): """Load material package""" - self._material_packages = [get_module(p) for p in package] + self._material_packages = [] + for p in package: + try: + self._material_packages.append(get_module(p)) + except ImportError as ie: # noqa: PERF203 + bluemira_warn(f"Unable to find {p}: {ie.msg}") def _get_material(self, name): name = name.replace("-", "_") From 4c402128c7e601cb92cc3620473859bb4733d3d1 Mon Sep 17 00:00:00 2001 From: james <81617086+je-cook@users.noreply.github.com> Date: Fri, 21 Nov 2025 12:09:25 +0000 Subject: [PATCH 17/17] =?UTF-8?q?=F0=9F=8E=A8=20Move=20vtk=20utility?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- bluemira/codes/openmc/output.py | 34 +------------------------------- bluemira/utilities/tools.py | 35 +++++++++++++++++++++++++++++++++ 2 files changed, 36 insertions(+), 33 deletions(-) diff --git a/bluemira/codes/openmc/output.py b/bluemira/codes/openmc/output.py index 3d8bb980a5..a22efb0a40 100644 --- a/bluemira/codes/openmc/output.py +++ b/bluemira/codes/openmc/output.py @@ -28,6 +28,7 @@ from bluemira.radiation_transport.neutronics.zero_d_neutronics import ( ZeroDNeutronicsResult, ) +from bluemira.utilities.tools import numpy_to_vtk def get_percent_err(row): @@ -611,39 +612,6 @@ def from_run( ) -# TODO @je-cook: move this function somewhere sensible # noqa: TD003 -import vtk # noqa: E402 -from vtkmodules.util import numpy_support # noqa: E402 - - -def numpy_to_vtk(data, output_name, scaling=(1, 1, 1)): - """Convert a numpy array to a VTK image data file.""" - data_type = vtk.VTK_FLOAT - shape = data.shape - - flat_data_array = data.flatten() - vtk_data = numpy_support.numpy_to_vtk( - num_array=flat_data_array, deep=True, array_type=data_type - ) - vtk_data.SetName(output_name) - - half_x = int(0.5 * scaling[0] * (shape[0] - 1)) - half_y = int(0.5 * scaling[1] * (shape[1] - 1)) - half_z = int(0.5 * scaling[2] * (shape[2] - 1)) - - img = vtk.vtkImageData() - img.GetPointData().SetScalars(vtk_data) - img.SetSpacing(scaling[0], scaling[1], scaling[2]) - img.SetDimensions(shape[0], shape[1], shape[2]) - img.SetOrigin(-half_x, -half_y, -half_z) - - # Save the VTK file - writer = vtk.vtkXMLImageDataWriter() - writer.SetFileName(f"{output_name}.vti") - writer.SetInputData(img) - writer.Write() - - @dataclass class NeutronicsOutputParams(ParameterFrame): """ diff --git a/bluemira/utilities/tools.py b/bluemira/utilities/tools.py index d825790bcd..a9e2d4221e 100644 --- a/bluemira/utilities/tools.py +++ b/bluemira/utilities/tools.py @@ -28,8 +28,10 @@ import nlopt import numpy as np import numpy.typing as npt +import vtk from PySide6.QtWidgets import QApplication from matplotlib import colors +from vtkmodules.util import numpy_support from bluemira.base.constants import E_I, E_IJ, E_IJK from bluemira.base.file import force_file_extension @@ -121,6 +123,39 @@ def json_writer( return None +# ===================================================== +# vtk utilities +# ===================================================== + + +def numpy_to_vtk(data, output_name, scaling=(1, 1, 1)): + """Convert a numpy array to a VTK image data file.""" + data_type = vtk.VTK_FLOAT + shape = data.shape + + flat_data_array = data.flatten() + vtk_data = numpy_support.numpy_to_vtk( + num_array=flat_data_array, deep=True, array_type=data_type + ) + vtk_data.SetName(output_name) + + half_x = int(0.5 * scaling[0] * (shape[0] - 1)) + half_y = int(0.5 * scaling[1] * (shape[1] - 1)) + half_z = int(0.5 * scaling[2] * (shape[2] - 1)) + + img = vtk.vtkImageData() + img.GetPointData().SetScalars(vtk_data) + img.SetSpacing(scaling[0], scaling[1], scaling[2]) + img.SetDimensions(shape[0], shape[1], shape[2]) + img.SetOrigin(-half_x, -half_y, -half_z) + + # Save the VTK file + writer = vtk.vtkXMLImageDataWriter() + writer.SetFileName(f"{output_name}.vti") + writer.SetInputData(img) + writer.Write() + + # ===================================================== # csv writer utilities # =====================================================