Problem with dagmc simulation

[gglamzn]

OneDrive_1_16-07-2024.zip](https://github.com/user-attachments/files/16256217/OneDrive_1_16-07-2024.zip)`
Hello,

I wanted to perform a DAGMC simulation of a geometry as follows:

It is a very simple geometry with two wall thicknesses and a floor, all made of concrete.

Here is how I proceed to convert and simulate the geometry (.step file) created in a CAD software:

First, I use GMSH to create the mesh, and then I use Python:

First function:

convert_file.py

`# Libraries

from cad_to_dagmc import MeshToDagmc
import cadquery as cq
import os
from cad_to_dagmc import CadToDagmc
import time

def output_path(output_dir, output_file_name):
    """
    Useful for having the output_path of the h5m file in the other functions.

    Args:
        output_dir (str): The directory where the output h5m file will be saved.
        output_file_name (str): The name of h5m output file.
    
    Returns:
        output_path (str): The path of the h5m file
    """
    return os.path.join(output_dir, output_file_name)

def mesh_geometry_to_h5m(material_tags_list, input_file_path, output_dir, output_file_name):
    """
    Converts a MESH file to a DAGMC h5m file usable for simulations on OpenMC.

    Args:
        material_tags_list (list): The materials tags that will be applied to the volumes.
        input_file_path (str): The path of the mesh file to convert 
        output_dir (str): The directory where the output h5m file will be saved.
        output_file_name (str): The name of h5m output file.
    
    Returns:
        None
    """
    # Initialize MeshToDagmc object
    mesh_converter = MeshToDagmc(filename=input_file_path)
    file_path = output_path(output_dir, output_file_name)
    # Export DAGMC h5m file
    mesh_converter.export_dagmc_h5m_file(material_tags=material_tags_list,
                                         filename=file_path)

material_tags_list = ['mat_fw', 'mat_sw', 'mat_ceiling']
input_file_path = 'path/to/parametric_model.msh'
output_dir = 'dir'
output_file_name = 'parametric_gmsh.h5m'

start_time = time.time()
mesh_geometry_to_h5m(material_tags_list, input_file_path, output_dir, output_file_name)
end_time = time.time()

simulation_time = end_time - start_time 
print(f'Meshing by gmsh took {simulation_time:.2f} `seconds')`

Conversion from .msh to .h5m, applying the different materials to the 3 volumes (WALL1, WALL2, floor).

Second function:

gmsh_sim.py

CAD Geometry Simulation Script

from cad_to_dagmc import CadToDagmc
import openmc
import os
import random
import time

def simulation_cad_geometry(material_tags_list, output_xml_path, input_dir, input_file_name):
    """
    Run an OpenMC simulation.

    Args:
        material_tags_list (list): The materials tags that will be applied to the volumes.
        output_xml_path (str): The path of the XML result file.
        input_dir (str): The directory where the output h5m file is located.
        input_file_name (str): The name of the input (.h5m) file.

    Returns:
        tuple: A tuple containing the result file and mesh object.
              result_file (str): The path to the file containing the results.
              mesh (openmc.RegularMesh): The mesh object used in the simulation.
    """
    try:
        # Define materials based on tags
        materials = openmc.Materials()

        for tag in material_tags_list:
            material = openmc.Material(name=tag)
            if tag == 'mat_fw':
                material.add_nuclide('H1', 0.010, percent_type='wo')
                material.add_nuclide('C0', 0.001416, percent_type='wo')
                material.add_nuclide('O16', 0.5267, percent_type='wo')
                material.add_nuclide('Na23', 0.016, percent_type='wo')
                material.add_nuclide('Mg24', 0.002, percent_type='wo')
                material.add_nuclide('Al27', 0.034, percent_type='wo')
                material.add_element('Si', 0.337, percent_type='wo')
                material.add_nuclide('B10', 0.00033, percent_type='wo')
                material.add_nuclide('K39', 0.013, percent_type='wo')
                material.add_nuclide('Ca40', 0.044, percent_type='wo')
                material.add_nuclide('Fe56', 0.014, percent_type='wo')
                material.set_density('g/cm3', 2.3)
            elif tag == 'mat_sw':
                material.add_nuclide('H1', 0.006, percent_type='wo')
                material.add_nuclide('C0', 0.0024, percent_type='wo')
                material.add_nuclide('O16', 0.5973, percent_type='wo')
                material.add_nuclide('Na23', 0.016, percent_type='wo')
                material.add_nuclide('Mg24', 0.002, percent_type='wo')
                material.add_nuclide('Al27', 0.038, percent_type='wo')
                material.add_element('Si', 0.337, percent_type='wo')
                material.add_nuclide('K39', 0.013, percent_type='wo')
                material.add_nuclide('Ca40', 0.044, percent_type='wo')
                material.add_nuclide('Fe56', 0.014, percent_type='wo')
                material.set_density('g/cm3', 2.3)
            elif tag == 'mat_ceiling':
                material.add_nuclide('H1', 0.010, percent_type='wo')
                material.add_nuclide('C0', 0.001416, percent_type='wo')
                material.add_nuclide('O16', 0.5267, percent_type='wo')
                material.add_nuclide('Na23', 0.016, percent_type='wo')
                material.add_nuclide('Mg24', 0.002, percent_type='wo')
                material.add_nuclide('Al27', 0.034, percent_type='wo')
                material.add_element('Si', 0.337, percent_type='wo')
                material.add_nuclide('B10', 0.00033, percent_type='wo')
                material.add_nuclide('K39', 0.013, percent_type='wo')
                material.add_nuclide('Ca40', 0.044, percent_type='wo')
                material.add_nuclide('Fe56', 0.014, percent_type='wo')
                material.set_density('g/cm3', 2.3)
            materials.append(material)

        air = openmc.Material(5, "air")
        air.add_nuclide('N14', 0.79)
        air.add_nuclide('O16', 0.21)
        air.set_density('g/cm3', 0.0012)
        materials.append(air)

        materials.export_to_xml()
        materials.cross_sections = 'path/to/crosssections'

        h5mfile_path = os.path.join(input_dir, input_file_name)
        dag_univ = openmc.DAGMCUniverse(h5mfile_path)

        L_square_sw = 1400
        height_root = 100
        height_ground = 100
        height_walls = 100
        height_ceiling = 100

        dag_univ = openmc.DAGMCUniverse(h5mfile_path, auto_geom_ids=True)
        bounding_region = dag_univ.bounding_region(bounded_type='box', boundary_type='vacuum', starting_id=10000, padding_distance=0.0)
        containing_cell = openmc.Cell(region=bounding_region, fill=dag_univ)

        air_cell = openmc.Cell(fill=air)
        air_universe = openmc.Universe(cells=[air_cell])
        containing_cell.fill = air_universe

        root_universe = openmc.Universe(cells=[containing_cell])
        geometry = openmc.Geometry(root_universe)
        geometry.export_to_xml()

        btch = 10
        part = 1 * 10**6
        inbatch = 0

        my_settings = openmc.Settings()
        my_settings.batches = btch
        my_settings.inactive = inbatch
        my_settings.particles = part
        my_settings.seed = int(random.random() * 100) + 1
        my_settings.particle = "neutron"
        my_settings.run_mode = 'fixed source'

        my_source = openmc.Source()
        my_source.space = openmc.stats.Box((-25, -10, -10), (25, 10, 10), only_fissionable=False)
        my_source.angle = openmc.stats.Isotropic()
        my_source.energy = openmc.stats.Discrete([2.5e6], [1])
        my_settings.source = my_source
        my_settings.survival_biasing = True
        my_settings.photon_transport = True
        my_settings.export_to_xml()

        mesh = openmc.RegularMesh()
        num_divisions_x_y = int(L_square_sw / 10)
        num_divisions_z = int((height_ground + height_ceiling + height_walls) / 10)
        mesh.dimension = [num_divisions_x_y, num_divisions_x_y, num_divisions_z]

        mesh.lower_left = [-L_square_sw / 2, -L_square_sw / 2, -height_ground]
        mesh.upper_right = [+L_square_sw / 2, +L_square_sw / 2, height_walls]

        mesh_filter = openmc.MeshFilter(mesh)
        mesh_tally = openmc.Tally(name='tallies_on_mesh')
        mesh_tally.filters = [mesh_filter]
        mesh_tally.scores = ['flux']

        my_tallies = openmc.Tallies([mesh_tally])
        my_model = openmc.Model(materials=materials, geometry=geometry, settings=my_settings, tallies=my_tallies)

        start_time = time.time()
        output_filename = my_model.run()
        end_time = time.time()
        simulation_time = end_time - start_time 
        print(f'Simulation time is {simulation_time:.2f} seconds')

        return 1

    except Exception as e:
        print(f"An error occurred: {e}")

output_xml_path = 'building_model_gmsh.xml'
material_tags_list = ['mat_fw', 'mat_sw', 'mat_ceiling']
input_dir = ''
input_file_name = 'parametric_python.h5m'

simulation_cad_geometry(material_tags_list, output_xml_path, input_dir, input_file_name)

The goal of my import is to have a box surrounding my CAD geometry by about 1 meter on all sides. I also want to fill the box with air, knowing that the concrete materials have already been assigned to the volumes.

I perform my simulation and I do not get anything that resembles what might be expected. There is no attenuation in the flux at the walls.

In my opinion, the function that collects data from the output file tallies.out is correct as I have used it in other simulations and it displayed expected results. I think that the way I coded the geometry import in gmsh_sim.py and the setup of the cells and universes might be incorrect, and this could be the source of the wrong results . Finally , to compare , I added a plot of openmc simulation with an openmc-made geometry wich is similar to the cad model so you can see what I would want to except orignally from the malfunctionnin simulation

Could you please help me?

Have a good day.

Best regards.

[shimwell]

Thanks for raising the issue and adding so much detail.

I wonder if it is possible to include the CAD file in the attached zip file.

One thing I would check is the units, often cad produces geometry in mm and openmc/dagmc assume accept units in cm so that might be the issue. I think there is a scaling option on cad-to-dagmc so that is easy to fix.

[gglamzn]

Hello,

I have a zip file containing the codes as well as the CAD geometry in .step format. For your information, the geometry is based on a square with dimensions of 14x14 meters. It consists of three volumes: the two walls and the floor, all made of different varieties of concrete.

Have a good day,
Best regards

simulations_files.zip

[shimwell]

I loaded the step file and I see the coordinates of the inner square are much bigger than the plotted region

the mouse icon shows this area is 1500mm which openmc / dagmc will assume is in cm so it will be 1500cm.

perhaps you can scale your gmesh or we can add a scaling argument to the MeshToDagmc class

so the mesh would need scaling to put it into the plotted region

[gglamzn]

Hello,

Indeed, the point you're showing me is at the correct coordinates, at least the ones I intended on the CAD model when I did the model .
I understand the issue you're pointing out. How can I modify the functions to incorporate the scaling you suggest?

It's true that when I created the geometry in the CAD software, I used a scale in meters. Perhaps for future modeling, I should use a scale in centimeters to match the base scale used in OpenMC/DAGMC ?

Have a nice day.

Best regards,

[gglamzn]

Hello,

I am sorry to bother you again, but despite your advice, I have not been able to solve my problem.

I did not find the scaling button on Gmsh.
I am having difficulty understanding how I can implement a scaling argument in the MeshToDagmc class function.
Could you please provide me with some assistance?

Best regards.

[shimwell]

The values of the coordiates could be divided by a scaling factor

[gglamzn]

Hi ,

I get that but how could I get acess to the step file coordinates?

[shimwell]

As your scripts pass in a gmsh mesh file then I think you might want to scale the mesh. I would recommend looking in the gmsh docs for a suitable function. Perhaps Mesh.MeshSizeFactor does the job. If you find a way to do it in gmsh then we can add an argument to the MeshToDagmc class that makes use of the gmsh capability.

[gglamzn]

Hello,

Since the last time, I searched the gmsh library for a function to scale the geometry. Unfortunately, I got a bit lost in the documentation and decided to proceed differently.

I found a method that seems to work:

I create my CAD geometry in centimeters, with a scale factor of 1/10. When I import the geometry into gmsh, a factor of 10 seems to be applied to the geometry.

Since these two factors cancel each other out, I then convert the geometry to .h5m format.

I then verify the consistency of the distances with openmc_plot, and with this method, it seems that I achieve the desired importation with the correct dimensions. I'm able to conduct simulations that more closely resemble what I should expect.

I have two questions:

Is this method, despite being somewhat unconventional but seemingly functional, actually a mirage that I have misinterpreted?
How could I create a box around my geometry and fill all the empty space, including the space between the walls of the CAD geometry, with air if possible?

I have attached several files within this message: the different basement files at various stages of conversion, the functions to convert to .h5m and display the geometry, and some images.
docs.zip

Best regards

[shimwell]

Is this method, despite being somewhat unconventional but seemingly functional, actually a mirage that I have misinterpreted?

units in CAD can be a bit of a mess, the often (but not always default to mm) sometimes not included in the step file, sometimes they are. Most DAGMC workflows ignore the units in CAD files and just assume the numbers provided are in cm. So I guess what you have is the normal way of going from CAD to neutronics.

How could I create a box around my geometry and fill all the empty space, including the space between the walls of the CAD geometry, with air if possible?

I think you want the implicit material option. DAGMC adds a void material to any space not covered by cad geometry. this is called an implicit void. The material assignment can be changed and you can ask DAGMC to use air instead of a void.

take a peak at this section of code

cad_to_dagmc/src/cad_to_dagmc/core.py

Lines 337 to 341 in 5cd38e4

	def export_dagmc_h5m_file(
	self,
	material_tags: typing.Iterable[str],
	implicit_complement_material_tag: typing.Optional[str] = None,
	filename: str = "dagmc.h5m",

see it has a argument called implicit_complement_material_tag . you would assign that to a string (perhaps 'air') and make a material named 'air' on the openmc side and that should work 🤞

[gglamzn]

Hi ,
I want to thank you for your precious insinghts throughout all the problems I faced .

I've been able to fill with air all the empty space as the png that I attached to this mail.

Best regards.
cad_filled.zip