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.