fusion-energy / openmc-dagmc-wrapper Goto Github PK

Currently, if users set a tally name on instanciation, it gets overwritten:

my_tally = odw.CellTally(
    "(n,Xa)", target="copper", materials=materials, name="my_custom_name"
)

I think it can be set afterwards, but that's confusing.

Currently the method NeutronicsModel.export_xml is like 450 line long which makes it hard to read and also hard to override when needed.

Maybe it could be split into at least three parts as the tet_mesh_filename, mesh_tally_3D and mesh_tally_2D seem independent enough

As we require dagmc and openmc for this package to work the pip install will not be enough to install everything that is needed.

A conda install will be needed to get everything deployed on the users computer.

currently the openmc conda install does not include dagmc so we can't fully complete this issue.

I've made a branch to get the work started https://github.com/fusion-energy/openmc-dagmc-wrapper/tree/conda_install

Will finish this once the conda install for openmc includes dagmc

The 2d and 3d mesh resolution can be tricky for users to find when the mesh is not of equal sizes in all directions.

A feature should be added that allows users to specify the number of mesh elements required and it makes the mesh with cubes / squares accordingly

The resulting 2D or 3D mesh will have elements with equal length sides

pshriwise wrote this code snippet a while back and I have not got around to incorporating it quite yet.

I shall get this done next week

#! /bin/env python

​

from argparse import ArgumentParser

​

from pymoab import core, types

​

​

_VALID_GEOM_TYPES =  ("Group", "Volume", "Surface", "Curve", "Vertex")

​

​

def ids_by_type(mb, geom_type):

    """

    Returns the ids for a given geometry entity type.

​

    Parameters

    ----------

​

    mb : pymoab.core

        A PyMOAB instance with the file of interest already loaded

    geom_type : string

        The type of geometry entity to provide IDs for.

        One of ("Group", "Volume", "Surface", "Curve", "Vertex")

​

    Returns

    -------

    list[int] : a sorted list of IDs for the geometry type

​

    """

    if geom_type not in _VALID_GEOM_TYPES:

        raise


    # retrieve the category tag on the instance

    try:

        cat_tag = mb.tag_get_handle(types.CATEGORY_TAG_NAME)

    except types.MB_ENTITY_NOT_FOUND:

        raise RuntimeError("The category tag could not be found in the PyMOAB instance."

                           "Please check that the DAGMC file has been loaded.")

    # get the id tag

    gid_tag = mb.tag_get_handle(types.GLOBAL_ID_TAG_NAME)

    # get the set of entities using the provided category tag name

    # (0 means search on the instance's root set)

    ents = mb.get_entities_by_type_and_tag(0, types.MBENTITYSET, [cat_tag], [geom_type])

    # retrieve the IDs of the entities

    ids = mb.tag_get_data(gid_tag, ents).flatten()

    return sorted(list(ids))

​

if __name__ == "__main__":

    ap = ArgumentParser(description="Program to report DAGMC IDs of various entities")

    ap.add_argument('filename', type=str, help='Filename of the DAGMC model')

    ap.add_argument('geom_type', type=str, help='Type of geometry entity to report IDs for')

    args = ap.parse_args()

    # create a new PyMOAB instance and load the specified DAGMC file

    mb = core.Core()

    mb.load_file(args.filename)


    # get the IDs for the requested geometry type

    ids = ids_by_type(mb, args.geom_type)

    print("{} IDs in {}: {}".format(args.geom_type, args.filename, ids))

It looks like the examples submersion_reactor.ipynb and text_example.ipynb have got out of date and should perhaps be updated or moved to the neutronics_workflow repo

Currently mesh tallies can accept bounding_box in the form of a list of two coordinates or a h5m file.

In the case of using a h5m file it would be useful to be able to also increase or decrease the mesh bounding set found.

This will allow us to mesh beyond the geometry limits

As a work around one can also place a small cube in the far corner of the geometry

It would make more sense if the order in mesh_resolution was the same as in plane ('xz', 'xy', or 'yz')

The newest image with python3.8 as default is currently used. However we could update to the newer miniconda image and create a python3.8 environment in the following way

FROM continuumio/miniconda3:4.10.3 as dependencies

# base image is python 3.9 but python 3.8 is needed for Cubit
# this creates and activates a new python 3.8 conda enviroment
RUN conda create --name py38 python=3.8
RUN echo "conda activate py38" >> ~/.bashrc
SHELL ["/bin/bash", "--login", "-c"]

TetMeshTally and TetMeshTallies should be imported in the init file.

we currently do

my_model = odw.Model(
    materials=materials, geometry=geometry, settings=settings, tallies=tallies
)

and openmc does this

my_model = openmc.model.Model(
    materials=materials, geometry=geometry, settings=settings, tallies=tallies
)

perhaps we should be more like openmc

I recall adding the remove_tag_from_h5m_file function with a plan to replace it with the tested https://github.com/svalinn/remove_dagmc_tags but have not gotten arround to it quite yet. Perhaps now is a good time to remove the whole function as the next release of DAGMC will be able to add graveyards automatically svalinn/DAGMC#714

On the latest fusion-neutronics-workflow image, running this is not possible:

import openmc_dagmc_wrapper as odw

geometry = odw.Geometry("geometry.h5m")

materials = odw.Materials(
    h5m_filename=geometry.h5m_filename,
    correspondence_dict={"water": "H2O"}
)

Traceback (most recent call last):
  File "/home/fusion-neutronics-workflow/examples/run_neutronics.py", line 8, in <module>
    materials = odw.Materials(
  File "/opt/conda/lib/python3.9/site-packages/openmc_dagmc_wrapper/Materials.py", line 25, in __init__
    self.set_openmc_materials()
  File "/opt/conda/lib/python3.9/site-packages/openmc_dagmc_wrapper/Materials.py", line 41, in set_openmc_materials
    openmc_material = create_material(material_tag, material_entry)
  File "/opt/conda/lib/python3.9/site-packages/openmc_dagmc_wrapper/utils.py", line 10, in create_material
    openmc_material = nmm.Material.from_library(
  File "/opt/conda/lib/python3.9/site-packages/neutronics_material_maker/material.py", line 865, in from_library
    return Material(name=name, **entry)
  File "/opt/conda/lib/python3.9/site-packages/neutronics_material_maker/material.py", line 199, in __init__
    self._make_openmc_material()
  File "/opt/conda/lib/python3.9/site-packages/neutronics_material_maker/material.py", line 645, in _make_openmc_material
    openmc_material = self._add_density(openmc_material)
  File "/opt/conda/lib/python3.9/site-packages/neutronics_material_maker/material.py", line 758, in _add_density
    from CoolProp.CoolProp import PropsSI
ModuleNotFoundError: No module named 'CoolProp'

after installing CoolProp with

pip install Coolprop

Traceback (most recent call last):
  File "/home/fusion-neutronics-workflow/examples/run_neutronics.py", line 8, in <module>
    materials = odw.Materials(
  File "/opt/conda/lib/python3.9/site-packages/openmc_dagmc_wrapper/Materials.py", line 25, in __init__
    self.set_openmc_materials()
  File "/opt/conda/lib/python3.9/site-packages/openmc_dagmc_wrapper/Materials.py", line 41, in set_openmc_materials
    openmc_material = create_material(material_tag, material_entry)
  File "/opt/conda/lib/python3.9/site-packages/openmc_dagmc_wrapper/utils.py", line 10, in create_material
    openmc_material = nmm.Material.from_library(
  File "/opt/conda/lib/python3.9/site-packages/neutronics_material_maker/material.py", line 865, in from_library
    return Material(name=name, **entry)
  File "/opt/conda/lib/python3.9/site-packages/neutronics_material_maker/material.py", line 199, in __init__
    self._make_openmc_material()
  File "/opt/conda/lib/python3.9/site-packages/neutronics_material_maker/material.py", line 645, in _make_openmc_material
    openmc_material = self._add_density(openmc_material)
  File "/opt/conda/lib/python3.9/site-packages/neutronics_material_maker/material.py", line 768, in _add_density
    raise ValueError(
ValueError: Material.temperature is needed to calculate the density

As I understand it, we need to provide a temperature to work out the density.

Problem: with this layout, we can't.

The workaround:

import openmc_dagmc_wrapper as odw
import neutronics_material_maker as nmm

geometry = odw.Geometry("geometry.h5m")

water_material = nmm.Material(
            "H2O",
            density="PropsSI('D', 'T', temperature, 'P', pressure, 'Water')/1000.",
            density_unit="g/cm3",
            percent_type="ao",
            temperature=120 + 273.15,  # K
            pressure=3.3e6,  # Pa
        )

materials = odw.Materials(
    h5m_filename=geometry.h5m_filename,
    correspondence_dict={"water": water_material }
)

Surely that's fine as is, just wondered if this use case could be simplified. Happy with the workaround anyway

Looks like setings.source defautls to a particular source. Perhaps this could be disabled and the ops package could be advertised via a ValueError if the source is missing

perhaps the units of heating when process_results(fusion_energy_per_pulse=1) is used should be joules per pulse instead of just joules

when making a 3d mesh tally like so

tally1 = odw.MeshTally3D(
    mesh_resolution=(200, 200, 200),
    bounding_box=None,
    tally_type="heating",
)

I am getting the following error

  File "/opt/conda/lib/python3.8/site-packages/openmc_dagmc_wrapper/Tally.py", line 254, in __init__
    self.add_mesh_filter(bounding_box)
  File "/opt/conda/lib/python3.8/site-packages/openmc_dagmc_wrapper/Tally.py", line 266, in add_mesh_filter
    mesh.lower_left = self.bounding_box[0]
TypeError: 'NoneType' object is not subscriptable

Effective dose should be added as an option on cells and 2d / 3d regular mesh tallies

This is example is for a surface but might be helpful
https://github.com/ukaea/openmc_workshop/blob/main/tasks/task_09_CSG_surface_tally_dose/1_surface_dose_from_gamma_source.ipynb

OpenMC docs on dose coefficients https://docs.openmc.org/en/stable/pythonapi/generated/openmc.data.dose_coefficients.html

ICRP publication https://journals.sagepub.com/doi/pdf/10.1016/j.icrp.2011.10.001

The default value of bounding_box in MeshTally2D is None.
However, it has to either be a string (h5m filename) or a list of tuples.
It can never be None, therefore, this argument shouldn't have a default value.

Running,

tally = odw.MeshTally2D('(n,Xa)', plane="xz")

produces (on v0.2.2):

  File "stage_3_run_neutronics_simulation.py", line 75, in <module>
    tally = odw.MeshTally2D('(n,Xa)', plane="xz")#, bounding_box=[(800-700, -1, -700), (800+700, 1, 700)])
  File "/opt/conda/lib/python3.8/site-packages/openmc_dagmc_wrapper/Tally.py", line 318, in __init__
    self.create_mesh(bounding_box)
  File "/opt/conda/lib/python3.8/site-packages/openmc_dagmc_wrapper/Tally.py", line 394, in create_mesh
    mesh.lower_left = self.bounding_box[0]
TypeError: 'NoneType' object is not subscriptable

Would the pip install of paramak-neutronics work without installing paramak first?

I don't see it in the requirements.

The plotly module is missing from the dependencies container. Perhaps because it's missing from requirements-test.txt

I tried to run a sim using correspondence_dict with FLiBe but I get an error asking me to give the temperature:

import openmc_dagmc_wrapper as odw

material_tag_to_material_dict = {
    "lead": "Lead",
    "flibe": "FLiBe",
    "inner_tank_wall": "eurofer",
}

materials = odw.Materials(
    h5m_filename='dagmc_not_merged.h5m',
    correspondence_dict=material_tag_to_material_dict,
)

  File "mwe.py", line 9, in <module>
    materials = odw.Materials(
  File "/opt/conda/lib/python3.8/site-packages/openmc_dagmc_wrapper/Materials.py", line 25, in __init__
    self.set_openmc_materials()
  File "/opt/conda/lib/python3.8/site-packages/openmc_dagmc_wrapper/Materials.py", line 41, in set_openmc_materials
    openmc_material = create_material(material_tag, material_entry)
  File "/opt/conda/lib/python3.8/site-packages/openmc_dagmc_wrapper/utils.py", line 10, in create_material
    openmc_material = nmm.Material.from_library(
  File "/opt/conda/lib/python3.8/site-packages/neutronics_material_maker/material.py", line 875, in from_library
    return Material(name=name, **entry)
  File "/opt/conda/lib/python3.8/site-packages/neutronics_material_maker/material.py", line 197, in __init__
    self._make_openmc_material()
  File "/opt/conda/lib/python3.8/site-packages/neutronics_material_maker/material.py", line 655, in _make_openmc_material
    openmc_material = self._add_density(openmc_material)
  File "/opt/conda/lib/python3.8/site-packages/neutronics_material_maker/material.py", line 765, in _add_density
    raise ValueError(
ValueError: Material.temperature is needed to calculate the density

A user would want to have multiple targets on a CellTally. For instance, to compute the TBR on volumes 1, 2 and 3:

my_tally = odw.CellTally(tally_type='TBR', target=[1, 2, 3])

This is not possible currently. We would have to modify CellTally.set_filters...

Or it might be easier to just let the users do

my_tally = openmc.Tally()
my_tally.scores = ['(n,Xt)']  # yes it's not TBR
my_tally.filters = [openmc.MaterialFilter(i) for i in [1, 2, 3]]

Version 1: abstraction (loss of flexibility) + maintainance/development to make it better
Version 2: 2 additional lines 😄

Looks like we need a test for the utils function remove_tag_from_h5m_file

There might be one in the paramak that got lost during the seperation of the two codes

the post processing could be easier if tally names follow a pattern.

perhaps something like this would allow the volume ids and material_tags to be found afterwards.

f"volume_{volume_number}_material_{material_tag}_{tally_type}"

Alternatively volume id and material id can be found from the filters. but only if both filters are added.

@RemDelaporteMathurin

We have a lot of places where we parse in the h5m geometry to check the bounding box or materials present. We are passing information between the class of settings, materials, tallies.

How about we just make the whole model in one class. Then all the methods have access to the h5m geometry and other aspects of the model. The example here uses the current approach is not much shorter than using native openmc. Combining into a single object could reduce that example to about half the line count

import openmc_dagmc_wrapper as odw

my_model = odw.model(
    h5m_filename='dagmc.h5m',
    materials = {
        "mat_inboard_tf_coils": "copper",
        "mat_center_column_shield": "tungsten",
        "mat_firstwall": "tungsten",
        "mat_blanket": "Li4SiO4",
        "mat_blanket_rear_wall": "eurofer",
        "mat_divertor_upper": "tungsten",
        "mat_divertor_lower": "tungsten",
    },
    tallies= [
        odw.MeshTally3D(voxels=1000, tally_type="(n,Xt)"),
        odw.MeshTally2D(tally_type="(n,Xt)", plane="xy")
    ],
    batches=4,
    particles=1000,
    source= odw.FusionRingSource(fuel="DT", radius=350, angles=(0, 3.14))
)

# starts the simulation
statepoint_file = my_model.run()

When running the following code I got the following error

  settings = odw.FusionSettings()
  settings.batches = 1
  settings.particles = 100

processing TBR
/opt/openmc/openmc/tallies.py:278: RuntimeWarning: invalid value encountered in true_divide
  self._std_dev[nonzero] = np.sqrt((self.sum_sq[nonzero]/n -

currently just 709 group structure is supported, ideally it would be good to offer options to the other energy spectra

def create_material(material_tag: str, material_entry):
    if isinstance(material_entry, str):
        openmc_material = nmm.Material.from_library(
            name=material_entry, material_id=None
        ).openmc_material
    elif isinstance(material_entry, openmc.Material):
        # sets the material name in the event that it had not been set
        openmc_material = material_entry
    elif isinstance(material_entry, (nmm.Material)):
        # sets the material tag in the event that it had not been set
        openmc_material = material_entry.openmc_material
    else:
        raise TypeError(
            "materials must be either a str, \
            openmc.Material, nmm.MultiMaterial or nmm.Material object \
            not a ",
            type(material_entry),
            material_entry,
        )
    openmc_material.name = material_tag
    return openmc_material

I find it more intuitive if it's just

def create_material(material_entry):
    if isinstance(material_entry, str):
        openmc_material = nmm.Material.from_library(
            name=material_entry, material_id=None
        ).openmc_material
    elif isinstance(material_entry, openmc.Material):
        # sets the material name in the event that it had not been set
        openmc_material = material_entry
    elif isinstance(material_entry, (nmm.Material)):
        # sets the material tag in the event that it had not been set
        openmc_material = material_entry.openmc_material
    else:
        raise TypeError(
            "materials must be either a str, \
            openmc.Material, nmm.MultiMaterial or nmm.Material object \
            not a ",
            type(material_entry),
            material_entry,
        )
    return openmc_material

and the name of the material is assigned in odw.Materials.set_openmc_materials

CylindricalMesh would be a very useful feature for fusion applications.

https://docs.openmc.org/en/stable/pythonapi/generated/openmc.CylindricalMesh.html#openmc.CylindricalMesh

Any objections to changing the repo name from openmc-dagmc-wrapper to openmc_dagmc_wrapper.

This would help as then we would have one name everywhere.

Python packages need _ in the import statement and - doesn't work

I think the largest_radius = max(max(bbox[0]), max(bbox[1])) logic when setting the radius of the vacuum boundary condition is not quite correct and cuts of the corners of the geometry within the sphere. Might have to add some trigonometry here

The default behaviour when running a tet mesh is to process the results and produce a vtk file.

However this can be disabled with UnstructuredMesh.output = False

Then we can call the write_data_to_vtk function in the wrapper https://github.com/openmc-dev/openmc/blob/eacd8f6bba414cac12f65059a3c17226fd0a201f/openmc/mesh.py#L722

We could do with adding a scaling factor to the write_data_to_vtk so that the results can be written out in whatever units required (e.g Watts would need scaling from eV to joules and by fusion power)

The compiling required to install all the dependencies takes around 90 mins.

Work could be done on the conda so they all work together which would speed up the install time

Release of MOAB v5.3 -> conda install moab gets updated and so does pymoab
New release of DAGMC can be issued
OpenMC conda can be changed to include DAGMC

It would also be nice to include:
Embree which has a conda install already
Double Down which I don't think has a conda install

If reflective surfaces are required, the variable geom in NeutronicsModel.export_xml() needs to be modified accordingly.

Currently, the only way of achieving this is to completely override the NeutronicsModel.export_xml() method.

It would be more convenient if an openmc.Geometry object could be passed.

To get the mesh tally results it would be useful to have a script that can load a vtk and take slices

DAGMC viz could be a solution but I believe this is python 2 and Visit code.

Perhaps a solution based on the well supported vtk python package would be ideal. vtk is already in the dependency stack and has examples here https://kitware.github.io/vtk-examples/site/Python/

I am keen to :

• modernize the setup.py

• add a conda install

• merge develop into main

• perhaps rename package fusion-openmc-dagmc-wrapper ?

• make a release

What do you think @RemDelaporteMathurin

Perhaps angles in degrees or radians should be added to this doc string

Elsewhere in the code we use radians