How to create photon sources/sources using a DAGMC in a way that will only create sources in materials and not void areas?

[nnelso13]

My DAGMC file contains a complex geometry that is filled by different materials and void regions throughout it:

I am trying to create photon sources within my model to calculate dose rates, but it seems to create them in the void regions as well. I think that I need to redefine my space definition or add a domain to the sources, but I have no idea on how to do this. Currently I have:

dagmc_bounds = dagmc_univ.bounding_box
min_x, min_y, min_z = dagmc_bounds.lower_left
max_x, max_y, max_z = dagmc_bounds.upper_right

source_space = openmc.stats.Box(
    lower_left=[min_x, min_y, min_z],
    upper_right=[max_x, max_y, max_z],
    only_fissionable=False
)
            
source = openmc.IndependentSource(
                space=source_space,
                energy=energy,
                particle="photon",
                strength=strength,
                #domains=[material]

But this gave a mesh plot of my tally like:

Suggesting that it was indeed creating sources in my void areas and not just my materials. I know that the source_space probably needs to be changed to only be the regions with material in my DAGMC, but I am unsure how to do this. Alternatively, I tried to set the domains to material, but this gave the error: " [ERROR: More than 95% of external source sites sampled were rejected. Please check your external source’s spatial definition]".

Thank you for taking the time to read this!

Cheers!

[Shimwell]

Hi nnelso13

Always keen to see openmc useage in fusion.

One way would be to use a mesh source.

Do you have the cad volume of your plasma geometry, if so you could do an unstructured mesh and make a really nice conformal source for the plasma.
Example tet meshing of here
Example 6 unstructured mesh source here

Also mesh sources can accept regular or cylindrical meshes so you can make a mesh source like that as well.
Exampe 7 structured mesh source here

These mesh source allow for different sources in each voxel so you can do things like get the ion density and make a stronger source where there are more ions, or get the ion temperature and make a different neutron energy distribution.

[nnelso13]

Thank you so much for the reply! I see what you are suggesting, and I think I have an idea. I am not trying to create a plasma source however, I did do this part successful using the “openmc_plasma_source” package for my neutron depletion:

my_source = TokamakSource(
    elongation=2.2,
    ion_density_centre=1.6e20,
    ion_density_peaking_factor=1.52,
    ion_density_pedestal=1.0e20,
    ion_density_separatrix=0.5e20,
    ion_temperature_centre=22,
    ion_temperature_peaking_factor=2.25,
    ion_temperature_pedestal=4.0,
    ion_temperature_separatrix=0.1,
    major_radius=610,
    minor_radius=120.0,
    pedestal_radius=114.0,
    mode="H",
    shafranov_factor=14.4,
    triangularity=0.625,
    ion_temperature_beta=7.5,
  ).make_openmc_sources()

This was the first part of my R2S script, and now I am trying to create the photon sources in my materials after the neutron activation in order to carry out a shutdown dose rate calculation. So far I have this:

results = openmc.deplete.Results("depletion_results.h5")

mesh2 = openmc.RegularMesh()
mesh2.dimension = [71, 29, 59]
mesh2.lower_left = [0, -382, -410]  # x,y,z coordinates start at 0 as this is a sector model
mesh2.upper_right = [1000, 0.0, 410]

my_gamma_settings = openmc.Settings()
my_gamma_settings.run_mode = "fixed source"
my_gamma_settings.batches = 10
my_gamma_settings.particles = 10000
my_gamma_settings.dagmc = True

energies, pSv_cm2 = openmc.data.dose_coefficients(particle="photon", geometry="AP")
dose_filter = openmc.EnergyFunctionFilter(
    energies, pSv_cm2, interpolation="cubic"
)

mesh_p_particle_filter = openmc.ParticleFilter(["photon"])
mesh_p_filter = openmc.MeshFilter(mesh2)
mesh_p_tally = openmc.Tally(name="p_flux_on_mesh")
mesh_p_tally.filters = [mesh_p_filter, dose_filter, mesh_p_particle_filter]
mesh_p_tally.scores = ["flux"]

my_photon_tallies = openmc.Tallies([mesh_p_tally])

dagmc_bounds = dagmc_univ.bounding_box
min_x, min_y, min_z = dagmc_bounds.lower_left
max_x, max_y, max_z = dagmc_bounds.upper_right

source_space = openmc.stats.Box(
    lower_left=[min_x, min_y, min_z],
    upper_right=[max_x, max_y, max_z],
    only_fissionable=False
)

statepoints_folder = Path('statepoints_folder')

for i_cool in range(1, len(timesteps)):
    photon_sources_for_timestep = []
    print(f"Making photon source for timestep {i_cool}")

    for material in mats:
        mat_id = material.id
        activated_mat = results[i_cool].get_material(str(mat_id))
        energy = activated_mat.get_decay_photon_energy(units='Bq')
        strength = energy.integral()

        if strength > 0.:
            print(f"Material ID: {mat_id}, Energy: {energy}, Strength: {strength}")

            source = openmc.IndependentSource(
                space=source_space,
                energy=energy,
                particle="photon",
                strength=strength,
                #domains=[material]
            )
            photon_sources_for_timestep.append(source)

    if not photon_sources_for_timestep:
        print(f"No valid photon sources for timestep {i_cool}")
        continue
    
    my_gamma_settings.source = photon_sources_for_timestep

    model_gamma = openmc.Model(geometry, mats, my_gamma_settings, my_photon_tallies)

    model_gamma.run(cwd=statepoints_folder / "photons" / f"photon_at_time_{i_cool}")

Which is what resulted in an output close to what I was looking for in the image in my post, but there was something amiss with the photon source creation making sources in void and material regions instead of just the materials. This came from the example (Task 11) here: link. If I understand correctly, you think it would be better to use a mesh source rather than an independent source in this scenario while using DAGMC for my geometry?

Thank you immensely!!!

[Shimwell]

Yep i think the gamma source space is not quite right.

source_space in the csg cell based example encompasses the csg cell.

I don’t think we have the cell bounding box for each cell in dagmc. Perhaps that is something pydagmc can already do or might do one day. But assuming we don’t have the cell bounding box for each dagmc cell. A mesh based approach doesn’t need that information as it gets activation on eqch cell voxel and turns each mesh element into a source.

[nnelso13]

I tried to follow the example from Task 7 to create a mesh source for my photons like this:

for i_cool in range(1, len(timesteps)):
    photon_sources_for_timestep = []
    print(f"Making photon source for timestep {i_cool}")

    for i in mesh.indices:
        mesh_index = (i[0]-1, i[1]-1, i[2]-1)
        volume = mesh.volumes[mesh_index]
        source = openmc.IndependentSource()
        
        found_activated_material = False
        for material in mats:
            mat_id = material.id
            try:
                activated_mat = results[i_cool].get_material(str(mat_id))
                energy = activated_mat.get_decay_photon_energy(units='Bq')
                source.energy = energy
                source.strength = volume
                found_activated_material = True
                break  # Break after finding the first valid activated material
            except ValueError:
                continue
        
        if not found_activated_material:
            source.strength = 0.0

        photon_sources_for_timestep.append(source)

    if not photon_sources_for_timestep:
        print(f"No valid photon sources for timestep {i_cool}")
        continue

    mesh_source = openmc.MeshSource(
        mesh=mesh,
        sources=photon_sources_for_timestep
    )

    my_gamma_settings.source = mesh_source

    model_gamma = openmc.Model(geometry, mats, my_gamma_settings, my_photon_tallies)
    model_gamma.run(cwd=statepoints_folder / "photons" / f"photon_at_time_{i_cool}")

but I am not sure if this is correct? I did also try modifying my original code:

dagmc_bounds = dagmc_univ.bounding_box
min_x, min_y, min_z = dagmc_bounds.lower_left
max_x, max_y, max_z = dagmc_bounds.upper_right

source_space = openmc.stats.Box(
    lower_left=[min_x, min_y, min_z],
    upper_right=[max_x, max_y, max_z],
    only_fissionable=False
)

statepoints_folder = Path('statepoints_folder')

for i_cool in range(1, len(timesteps)):
    photon_sources_for_timestep = []
    print(f"Making photon source for timestep {i_cool}")

    for material in mats:
        mat_id = material.id
        activated_mat = results[i_cool].get_material(str(mat_id))
        energy = activated_mat.get_decay_photon_energy(units='Bq')
        strength = energy.integral()

        if strength > 0.:
            print(f"Material ID: {mat_id}, Energy: {energy}, Strength: {strength}")
            source = openmc.IndependentSource(
                space=source_space,
                energy=energy,
                particle="photon",
                strength=strength,
                constraints={'domains': [material]},
            )
            photon_sources_for_timestep.append(source)

    if not photon_sources_for_timestep:
        print(f"No valid photon sources for timestep {i_cool}")
        continue
    
    my_gamma_settings.source = photon_sources_for_timestep

    model_gamma = openmc.Model(geometry, mats, my_gamma_settings, my_photon_tallies)

    model_gamma.run(cwd=statepoints_folder / "photons" / f"photon_at_time_{i_cool}")

after upgrading to OpenMC 0.15 due to the inclusions:

constraints (dict) – Constraints on sampled source particles. Valid keys include ‘domains’, ‘time_bounds’, ‘energy_bounds’, ‘fissionable’, and ‘rejection_strategy’. For ‘domains’, the corresponding value is an iterable ofopenmc.Cell, openmc.Material, or openmc.Universe for which sampled sites must be within. For ‘time_bounds’ and ‘energy_bounds’, the corresponding value is a sequence of floats giving the lower and upper bounds on time in [s] or energy in [eV] that the sampled particle must be within. For ‘fissionable’, the value is a bool indicating that only sites in fissionable material should be accepted. The ‘rejection_strategy’ indicates what should happen when a source particle is rejected: either ‘resample’ (pick a new particle) or ‘kill’ (accept and terminate).

this did seem to do something differently:

It literally outputted a photo with just the void regions with sources??? Confused why it would do that if the constraint is supposed to restrict the source sampling to whatever it is set to if I understand that correctly.