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/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/base/reactor.py b/bluemira/base/reactor.py index 5abe97d846..865663cf59 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_properties("material", first=False)[0] + @staticmethod def _validate_cad_dim(dim: DIM_3D | DIM_2D) -> None: """ 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..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): @@ -58,8 +59,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 OpenMCResult: +class OpenMCCSGResult(OpenMCResultBase): """ Class that looks opens up the openmc universe from the statepoint file, so that the dataframes containing the relevant results @@ -104,7 +164,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 +185,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" @@ -141,10 +201,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 +276,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 +294,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""" @@ -520,6 +540,78 @@ def _tabulate( ) +@dataclass +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_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_mesh_tally = statepoint.get_tally(name="heating_on_mesh") + flux_mesh_tally = statepoint.get_tally(name="flux_on_mesh") + + 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", + 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) + + 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, + ) + + @dataclass class NeutronicsOutputParams(ParameterFrame): """ @@ -540,16 +632,41 @@ class NeutronicsOutputParams(ParameterFrame): peak_div_cu_dpa_rate: Parameter[float] @classmethod - def from_openmc_csg_result(cls, result: OpenMCResult): + 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): """ Produce output parameters from an OpenMC CSG result """ 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), @@ -559,32 +676,20 @@ def from_openmc_csg_result(cls, result: OpenMCResult): # 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 ba7cf9e609..9cbf54172c 100644 --- a/bluemira/codes/openmc/solver.py +++ b/bluemira/codes/openmc/solver.py @@ -7,19 +7,18 @@ 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 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 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 ( @@ -30,20 +29,25 @@ CodesTeardown, ) from bluemira.codes.openmc.make_csg import ( - BlanketCellArray, BluemiraNeutronicsCSG, CellStage, - DivertorCellArray, 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, + OpenMCDAGMCResult, +) 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 +from bluemira.radiation_transport.neutronics.neutronics_axisymmetric import ( + NeutronicsReactor, +) if TYPE_CHECKING: from matproplib.conditions import OperationalConditions @@ -111,92 +115,89 @@ 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""" + 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 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.pre_cell_model = pre_cell_model self.materials = materials - self.matlist = attrgetter( - "outb_sf_mat", - "outb_fw_mat", - "outb_bz_mat", - "outb_mani_mat", - "outb_vv_mat", - "divertor_mat", - "div_fw_mat", - "tf_coil_mat", - ) - @contextmanager def _base_setup(self, run_mode, *, debug: bool = False): from openmc.config import config # noqa: PLC0415 - self.files_created = set() - folder = run_mode.name.lower() - cwd = Path(self.out_path, folder) - cwd.mkdir(parents=True, exist_ok=True) - config["cross_sections"] = self.cross_section_xml - self.settings = openmc.Settings( - run_mode=run_mode.value, output={"summary": False} - ) - 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( - self, - run_mode, - tally_function: TALLY_FUNCTION_TYPE, - cell_arrays: CellStage, - material_list: list[openmc.Material], - ): - out_path = Path(self.out_path, run_mode.name.lower(), "tallies.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) + if isinstance(self.materials, MaterialsLibrary): + model.materials = self.materials.get_all_materials() + else: + model.materials = self.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(material_list, cell_arrays): + for name, scores, filters in tally_function(self.tally_mats, self.tally_geom): tally = openmc.Tally(name=name) tally.scores = [scores] - tally.filters = filters + if filters is not None: + tally.filters = filters tallies_list.append(tally) - openmc.Tallies(tallies_list).export_to_xml(out_path) - self.files_created.add(out_path) + return openmc.Tallies(tallies_list) - def run( + @abstractmethod + def plot( self, run_mode, runtime_params, @@ -205,81 +206,201 @@ def run( 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") + ): ... - def plot( + @abstractmethod + def volume( self, run_mode, runtime_params, - _eq, - _source_params, - _tally_function, + eq, + source_params, + tally_function, *, 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) + ): ... - def volume( + def run( self, run_mode, runtime_params, - _source_params, - _tally_function, + 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 = self.pre_cell_model.bounding_box + 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 - 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") + +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.mat_list = attrgetter( + "outb_sf_mat", + "outb_fw_mat", + "outb_bz_mat", + "outb_mani_mat", + "outb_vv_mat", + "divertor_mat", + "div_fw_mat", + "tf_coil_mat", + ) + + @property + def tally_mats(self) -> list[openmc.Material]: + return self.mat_list(self.materials) + + @property + def tally_geom(self) -> CellStage: + return self.cell_arrays + + def _create_geometry(self): + universe = openmc.Universe(cells=self.cell_arrays.cells) + geometry = openmc.Geometry(universe) + return universe, geometry + + 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 + ) + + def volume( + self, + run_mode, + 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, + ): + super().__init__(codes_name, cross_section_xml, eq, source, materials) + self.dag_model_path = dag_model_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 + + @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, + *_args, + debug: bool = False, + ) -> tuple[openmc.Model, PlotConfig]: + return self._plot( + run_mode, runtime_params, self.universe.bounding_box, debug=debug + ) + + def volume( + self, + run_mode, + runtime_params, + *_args, + debug: bool = False, + ) -> tuple[openmc.Model, None]: + return self._volume( + run_mode, + runtime_params, + self.universe, + self.universe.bounding_box, + debug=debug, ) @@ -291,43 +412,78 @@ 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, - 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""" - 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): + def volume( + self, run_mode, model: openmc.Model, _config: None, *, 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): +class OpenMCCSGTeardown(CodesTeardown): """Teardown task for OpenMC solver""" def __init__( @@ -341,85 +497,89 @@ 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 - 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, - files_created, - source_rate: float, - source_triton_rate: float, + source_info: SourceInfo, statepoint_file, - *, - 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, - statepoint_file, + source_info.rate, + source_info.triton_rate, + 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, - 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 fig.axis def volume( self, _universe, - files_created, _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" ) - for cell in self.cells + for cell in self.cell_arrays.cells } +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, + ): + """Run stage for Teardown task""" + result = OpenMCDAGMCResult.from_run( + universe, source_info.rate, source_info.triton_rate, Path(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 fig.axis + + def volume( + self, + _universe, + _source_params, + _statepoint_file, + ) -> dict[int, float]: + """Stochastic volume stage for teardown task""" + raise NotImplementedError + + TALLY_FUNCTION_TYPE = Callable[ - [list[openmc.Material], BlanketCellArray, DivertorCellArray], + [list[openmc.Material], CellStage | openmc.Geometry], tuple[ str, str, @@ -428,26 +588,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 @@ -457,17 +615,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""" @@ -486,9 +633,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,47 +659,134 @@ def _single_run( runtime_params: OpenMCSimulationRuntimeParameters, *, debug=False, - ) -> OpenMCResult | dict[int, float]: + ) -> OpenMCCSGResult | dict[int, float]: + result = None + if setup := self._get_execution_method(self._setup, run_mode): + model, config = setup( + run_mode, + runtime_params, + self.eq, + source_params, + self.tally_function, + debug=debug, + ) + if run := self._get_execution_method(self._run, run_mode): + 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 + + +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.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.cell_arrays, + self.neutronics_model, ) + self._run = self.run_cls(self.out_path, self.name) self._teardown = self.teardown_cls( self.cell_arrays, - self.pre_cell_model, + self.neutronics_model, self.out_path, self.name, ) - result = None - if setup := self._get_execution_method(self._setup, run_mode): - result = setup( - run_mode, - runtime_params, - self.eq, - source_params, - self.tally_function, - debug=debug, - ) - 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", - ), - ) - return result + 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..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, ): @@ -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/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("-", "_") 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 # ===================================================== diff --git a/eudemo/config/build_config.json b/eudemo/config/build_config.json index d71201aa1f..5bc670fb9d 100644 --- a/eudemo/config/build_config.json +++ b/eudemo/config/build_config.json @@ -118,23 +118,44 @@ "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": true, + "dagmc_export_dir": "$path_expand:./dagmc_export", + "cross_section_xml": "$path_expand:./cross_section_data/cross_sections.xml", + "particles": 10000, + "batches": 2, + "converter_config": { + "converter_type": "fast_ctd", + "imprint_geometry": true, + "imprint_per_compound": false, + "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": { "run_mode": "run", 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..5135575cc4 100644 --- a/eudemo/eudemo/neutronics/run.py +++ b/eudemo/eudemo/neutronics/run.py @@ -7,8 +7,12 @@ from __future__ import annotations +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 from bluemira.radiation_transport.neutronics.blanket_data import ( @@ -27,7 +31,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 @@ -54,8 +58,8 @@ def _get_wires_from_components( ) -def run_neutronics( - params: dict | ParameterFrame, +def run_csg_neutronics( + params: ParameterFrame, build_config: dict, blanket: ComponentManager, vacuum_vessel: ComponentManager, @@ -64,7 +68,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 +114,73 @@ 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: # noqa: PLR2004 + res = outputs[0] + params.update_from_frame(outputs[1]) + else: + 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, + # reactor.blanket, + reactor.coil_structures, + ], + "group_by_materials": True, + }, + converter_config=build_config.get("converter_config", {}), + ) + + +def run_dagmc_neutronics( + reactor, + params, + build_config, + eq: Equilibrium, + 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=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, + ) + + return solver.execute(build_config.get("run_mode", "run")) diff --git a/eudemo/eudemo/reactor.py b/eudemo/eudemo/reactor.py index e170fe1a10..315fc0830e 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 @@ -712,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_neutronics( - reactor_config.params_for("Neutronics").global_params, - reactor_config.config_for("Neutronics"), + *run_csg_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, @@ -728,7 +706,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: @@ -882,14 +860,14 @@ def save_reactor(reactor, reactor_config, folder_name): n_TF=reactor_config.global_params.n_TF.value, ) - export_dagmc_model( + # TODO @je-cook: Put the results somewhere + dagmc_out = run_dagmc_neutronics( reactor, - reactor_config.config_for("CAD_Neutronics"), + reactor_config.params_for("Neutronics", "DAGMC"), + reactor_config.config_for("Neutronics", "DAGMC"), + reference_eq, ) - # 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( diff --git a/examples/radiation_transport/run_cad_neutronics.ex.py b/examples/radiation_transport/run_cad_neutronics.ex.py index 6a7e733e51..deeeb97d23 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,263 @@ # ## 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. +""" + +# %% +from pathlib import Path + +import matplotlib.pyplot as plt +import numpy as np + +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 +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 + + +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 diff --git a/pyproject.toml b/pyproject.toml index 30bcd271e6..6460c9808f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -91,13 +91,14 @@ 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"] dagmc = [ "cgal>=6.0.1", "fast_ctd @ git+https://github.com/Fusion-Power-Plant-Framework/fast_ctd@0.1.1", + "bluemira[openmc]" ] [build-system] 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