Axial flux plot

[jujutsukaisen]

I am simulating a subcritical pile where I want to plot the Flux vs Z-axis. My code is as follows:

my_tallies = openmc.Tallies()

# Define the 1D z-axis mesh tally with a RegularMesh
mesh1Dz = openmc.RegularMesh(name='1D_mesh')
mesh1Dz.lower_left = (30.0, 107.0, -88.0)
mesh1Dz.upper_right = (64.0, 120.6, 88.0)
mesh1Dz.dimension = [1, 1, 200] 
meshfilter1Dz = openmc.MeshFilter(mesh1Dz)

# Define a tally for various parameters
tally1Dz = openmc.Tally(name='1D_tally')
tally1Dz.scores = ['flux'] 
tally1Dz.filters = [meshfilter1Dz]
my_tallies.append(tally1Dz)
my_tallies.export_to_xml()

# tally 1D
tally1Dz = sp.get_tally(name='1D_tally')

# flux
tally1Dzflux = tally1Dz.get_slice(scores=['flux']).mean.flatten()/mesh1Dz.volumes[0][0][0]

# z axis
zaxis = np.linspace(mesh1Dz.lower_left[2], mesh1Dz.upper_right[2], mesh1Dz.dimension[2])

# plot
plt.figure()
plt.plot(zaxis, tally1Dzflux, color='red')
plt.xlabel('z position')
plt.ylabel('Flux [neutrons/cm$^2$-s]')
plt.legend()
plt.grid(True)
plt.tight_layout()

My questioon: Is the normalization ok. Do I divide by the volume of mesh whether it’s 1D, 2D or 3D mesh. [Note: I have defined source strength (n/s) in the source definition. So, didn’t multiply the tally by source strength]

[jujutsukaisen]

I am asking this as I plot a 2D mesh and normalize it with the 2D mesh volume, I am getting an order of 10 higher than what I got with normilizing 1D mesh flux.

# Create 2D mesh using RegularMesh and define tally
mesh2Dxz = openmc.RegularMesh().from_domain(
    my_geometry, 
    dimension=[200, 1, 200]
)

meshfilter2Dxz = openmc.MeshFilter(mesh2Dxz)
tally2Dxz = openmc.Tally(name='2D_tally')
tally2Dxz.filters = [meshfilter2Dxz]
tally2Dxz.scores = ['flux'] 
my_tallies.append(tally2Dxz)

# Accessing the flux tally and setting up the slice with the right shape
tally1Dz = sp.get_tally(name="2D_tally")
slice_tally = tally1Dz.get_slice(scores=['flux'])/mesh2Dxz.volumes[0][0][0]
slice_tally.mean.shape = (mesh2Dxz.dimension[0], mesh2Dxz.dimension[2])  # setting the resolution to the mesh dimensions

# Plotting the 2D data 
fig, ax = plt.subplots()
cax = ax.imshow(slice_tally.mean, 
                extent=mesh2Dxz.bounding_box.extent['xz'], 
                origin='lower', 
                norm=LogNorm())  # Apply logarithmic normalization

# Adding colorbar to represent the flux values with a log scale
fig.colorbar(cax, ax=ax, label='Flux (n/cm²-s)')

# Adding labels for axes
ax.set_xlabel("X Position (cm)")
ax.set_ylabel("Z Position (cm)")
ax.set_title("Flux Distribution in XZ Plane")

plt.show()

Should not both the answers be same or I am missing something?

The results you can see above.
If anyone can provide me with some reasoning, I’d be very grateful

[wahidluthfi]

Hi Amanda, welcome back to the community.
Just to make sure, your 2D mesh uses a different domain than your 1D mesh, right? I see that you use your my_geometry for your 2D mesh.domain size, while in 1D, you use (30.0, 107.0,- 88.0) to (64.0, 120.6, 88.0). because this domain size will affect the number of particles being scored to your tally bin. the bigger domain size will give you more particles within it, and that is where you need to know your mesh volume to calculate the average flux in that mesh.
Have you checked your calculated tally with some metric that you already know? i.e. Have you checked whether the total power of your model is equal to the power if you already knew the power? or using your source strength or total number of particles? Sometimes we just need to do some checks just to make sure that we have proper input for our case. maybe this discussion could help you

Photon flux with a layer of water and air

Hi Hedi, since you want to find the 10,000 particles that you declare as source strength, I am making a simple case that I think could help us to understand what openmc did in photon transport. I am making a simple 3D geometry consisting of a small cube with a size 1x1x1 cm in front of a shield cell, also a cube 1x1x1 on its back. A photon source was placed at -5cm (point source monodirectional monoenergetic) with a source strength of 10,000 as you want. [image] when all material being used…

[jujutsukaisen]

Thanks wahid for your response and pointing out that I didn’t use the same domain size which is actually from another test run and I totally overlooked before posting the 2D example.

However, I can assure you I have got the posted figures for the same domain.

my_tallies = openmc.Tallies()

# Define the 1D z-axis mesh tally with a RegularMesh and define tally
mesh1Dz = openmc.RegularMesh().from_domain(
    my_geometry, 
    dimension=[1, 1, 200]
)

meshfilter1Dz = openmc.MeshFilter(mesh1Dz)
tally1Dz = openmc.Tally(name='1D_tally')
tally1Dz.scores = ['flux'] 
tally1Dz.filters = [meshfilter1Dz]
my_tallies.append(tally1Dz)

# Create 2D mesh using RegularMesh and define tally
mesh2Dxz = openmc.RegularMesh().from_domain(
    my_geometry, 
    dimension=[200, 1, 200]
)

meshfilter2Dxz = openmc.MeshFilter(mesh2Dxz)
tally2Dxz = openmc.Tally(name='2D_tally')
tally2Dxz.filters = [meshfilter2Dxz]
tally2Dxz.scores = ['flux'] 
my_tallies.append(tally2Dxz)

my_tallies.export_to_xml()

# tally 1D
tally1Dz = sp.get_tally(name='1D_tally')
tally1Dzflux = tally1Dz.get_slice(scores=['flux']).mean.flatten()/mesh1Dz.volumes[0][0][0]
zaxis = np.linspace(mesh1Dz.lower_left[2], mesh1Dz.upper_right[2], mesh1Dz.dimension[2])

# Plotting the 1D data 
plt.figure()
plt.plot(zaxis, tally1Dzflux, color='red')
plt.xlabel('z position (cm)')
plt.ylabel('Flux [n/cm$^2$-s]')
plt.title("Flux Along Z-axis") 
plt.grid(True)

plt.gca().yaxis.set_major_formatter(ScalarFormatter(useMathText=True))
plt.ticklabel_format(axis="y", style="sci", scilimits=(0,0))  

plt.tight_layout()
plt.savefig("flux_along_z_axis.png")  
plt.show()

# Accessing the flux tally and setting up the slice with the right shape
tally1Dz = sp.get_tally(name="2D_tally")
slice_tally = tally1Dz.get_slice(scores=['flux'])/mesh2Dxz.volumes[0][0][0]
slice_tally.mean.shape = (mesh2Dxz.dimension[0], mesh2Dxz.dimension[2])  # setting the resolution to the mesh dimensions

# Plotting the 2D data 
fig, ax = plt.subplots()
cax = ax.imshow(slice_tally.mean, 
                extent=mesh2Dxz.bounding_box.extent['xz'], 
                origin='lower', 
                norm=LogNorm())  # Apply logarithmic normalization

fig.colorbar(cax, ax=ax, label='Flux (n/cm²-s)')
ax.set_xlabel("X Position (cm)")
ax.set_ylabel("Z Position (cm)")
ax.set_title("Flux Distribution in XZ Plane")
plt.savefig("flux_distribution_xz_plane.png", dpi=300)  

plt.show()

Regarding power, the value is unknown and I have set the source strength in my settings as

# Source definition
am_be_source = openmc.Source()
am_be_source.space = openmc.stats.Box((44.2, 93.2, -31.5), (44.7, 93.7, -30))
am_be_source.angle = openmc.stats.Isotropic()
am_be_source.energy = openmc.stats.Discrete([4.3e6], [1])
am_be_source.time = openmc.stats.Uniform(0, 1)
am_be_source.strength = 10 * 2 * 10**6


pu_be_source_1 = openmc.Source()
pu_be_source_1.space = openmc.stats.Box((37.2, 128.2, -31.5), (37.7, 128.7, -30))
pu_be_source_1.angle = openmc.stats.Isotropic()
pu_be_source_1.energy = openmc.stats.Discrete([4.1e6], [1])
pu_be_source_1.time = openmc.stats.Uniform(0, 1)
pu_be_source_1.strength = 1.8 * 10**6


pu_be_source_2 = openmc.Source()
pu_be_source_2.space = openmc.stats.Box((51.2, 128.2, -31.5), (51.7, 128.7, -30))
pu_be_source_2.angle = openmc.stats.Isotropic()
pu_be_source_2.energy = openmc.stats.Discrete([4.1e6], [1])
pu_be_source_2.time = openmc.stats.Uniform(0, 1)
pu_be_source_2.strength = 1.8 * 10**6


pu_be_source_3 = openmc.Source()
pu_be_source_3.space = openmc.stats.Box((37.2, 93.2, -31.5), (37.7, 93.7, -30))
pu_be_source_3.angle = openmc.stats.Isotropic()
pu_be_source_3.energy = openmc.stats.Discrete([4.1e6], [1])
pu_be_source_3.time = openmc.stats.Uniform(0, 1)
pu_be_source_3.strength = 1.8 * 10**6


pu_be_source_4 = openmc.Source()
pu_be_source_4.space = openmc.stats.Box((51.2, 93.2, -31.5), (51.7, 93.7, -30))
pu_be_source_4.angle = openmc.stats.Isotropic()
pu_be_source_4.energy = openmc.stats.Discrete([4.1e6], [1])
pu_be_source_4.time = openmc.stats.Uniform(0, 1)
pu_be_source_4.strength = 1.8 * 10**6


# Combine all sources into one list
sources = [am_be_source, pu_be_source_1, pu_be_source_2, pu_be_source_3, pu_be_source_4]

I have messeged you the script thinking the script might help you to understand you better what I am doing and where I might be mistaking. So, if you can spare some of your time to look at that, I’d be grateful and TIA

By the way, thanks for referring to that post but I had already gone through that before posting my question as I was still uncleared

[wahidluthfi]

Hi Amanda, sorry for the late response.
I have checked your input script, and I think your mesh tally definition is good, but the interpretations are quite complicated.
in your case, when you use 1D mesh while using the same domain (space) definition, then you are tallying a larger mesh which is then averaged with the larger volume. that’s why your flux value is lower in 1D z-axis plot in comparison to 2D xz plot because of the mesh region (volume) in your 2D xz is 96.8748 cc while in your 1D z is 19374.96 cc, larger but most of it has a low neutron population.
If you use the same volume by indicating a smaller domain, i.e. near your neutron source but with a proper x width to make sure that it covers the same region as your 2D XZ plot, then you will get the proper average neutron flux.
Also, I hope you noticed that your 2D xz plot itself is just reporting the average neutron flux along the y-axis from -14 to 165 when you see it from the xz plane. so the actual neutron flux locally at some point near your source could be higher.
I am troubleshooting your problem with eigenvalue mode but I think it still works with fixed source. here is the 2D xz plot in fixed source

and here is the 1D z-axis plot, sorry that I can’t get the x position as close as the peak neutron flux shown in xz plot (on the left). but you can change it later as you want.

while this is the result when I use the eigenvalue problem with 100k particles, 50 inactive, and 150 total batches. 2D xz plot

and here is the 1D z-axis plot,

I hope you can get my point since I am not good at English, sorry.