Geometric axial symmetry Boundary Conditions

User Support
1

Hello to everyone,

I am studying a fusion reactor, that for simplicity has been designed as cylindrical and I want to take advantage of the axial symmetry of the cylinder to reduce the computational cost of the calculation. The design includes an isotropic neutronic source at the center of the cylinder, and as far as I run the calculation in full geometry, no problem detected. I added the axial symmetry as follow

surf1 = openmc.ZCylinder(surface_id=1, x0=0, y0=0, r=radius[0], name=‘surf 1’)
surf2 = openmc.ZCylinder(surface_id=2, x0=0, y0=0, r=radius[1], name=‘surf 2’)
surf3 = openmc.ZCylinder(surface_id=3, x0=0, y0=0, r=radius[2], name=‘surf 3’)
surf3.boundary_type = ‘vacuum’
xplane = openmc.XPlane(x0=0, boundary_type=‘periodic’)
yplane = openmc.YPlane(y0=0, boundary_type=‘periodic’)
xplane.periodic_surface = yplane
min_z = openmc.ZPlane(z0=-5, boundary_type=‘reflective’)
max_z = openmc.ZPlane(z0=5, boundary_type=‘reflective’)
cell_1 = openmc.Cell(cell_id=1, name=‘cell 1’)
cell_2 = openmc.Cell(cell_id=2, name=‘cell 2’)
cell_3 = openmc.Cell(cell_id=3, name=‘cell 3’)
cell_1.region = -surf1 & -max_z & +min_z & -xplane & +yplane
cell_2.region = -surf2 & +surf1 & -max_z & +min_z & -xplane & +yplane
cell_3.region = -surf3 & +surf2 & -max_z & +min_z & -xplane & +yplane

with a source thought as follow:
source = openmc.Source(strength=1.0e+19)
source.space = openmc.stats.Box((0, 0, -5), (0, 0, 5))
source.angle = openmc.stats.Isotropic()

And I get the error:

ERROR: More than 95% of external source sites sampled were rejected. Please
** check your external source definition.**

As far as I understood from the OpenMC guide, the lines concerning the creation of both the x=0 and y=0 planes should be implemented as I reported here.
And yet, if I run a calculation, I get the afore mentioned error.

I want to thank in advance anyone that can give me some tips or further instructions and help me to solve my issue.

Best regards,
Tony

2

I suspect the issue is that your source is exactly coincident with the planes x=0 and y=0. If you give the source x and y coordinates some very small positive value, that may help to avoid the issue, e.g.

source.space = openmc.stats.Box((1e-6, 1e-6, -5), (1e-6, 1e-6, 5))
3

Firstly, I want to thank you for answering.

Back to the business, I have already thought to something like that, in fact, when I read it, I doubted it would have worked. I still tried, and it surprisingly worked for a couple of times (I copied and pasted your line). Then again, I get the following output:

===============> FIXED SOURCE TRANSPORT SIMULATION <===============

Simulating batch 1
Traceback (most recent call last):
File “/path/to/file/axial_symmetry_study/axial_SY_script.py”, line 85, in
openmc.run(threads=8)
File “/opt/anaconda3/envs/myenv/lib/python3.9/site-packages/openmc/executor.py”, line 218, in run
_run(args, output, cwd)
File “/opt/anaconda3/envs/myenv/lib/python3.9/site-packages/openmc/executor.py”, line 28, in _run
raise subprocess.CalledProcessError(p.returncode, ’ '.join(args),
subprocess.CalledProcessError: Command ‘openmc -s 8’ died with <Signals.SIGSEGV: 11>.
ERROR: More than 95% of external source sites sampled were rejected. Please
check your external source definition.

As before, I thought it could be something not related to the source position, but as soon as I run with the whole cylinder, the problem disappear.
Do you think it might help you to comprehend better the issue if I would upload the script?

Best regards,
Tony

4

@tony_emme Yes, if you’re able to upload the script, I can try to diagnose what’s going wrong with the source.

1 Like
5

Thanks a lot, here is the script (I tried to upload it, but I can’t):

import openmc
from openmc import stats
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

# Copper - Pure
copper = openmc.Material(name='Copper')
copper.add_element('Cu', 1.0, 'ao')
copper.set_density('g/cm3', 8.96)
# Silver - pure
silver = openmc.Material(name='Silver')
silver.add_element('Ag', 1.0, 'ao')
silver.set_density('g/cm3', 10.49)
materials = openmc.Materials((copper, silver))
materials.export_to_xml()

radius = np.array([100])
wall_thickness = np.array([25, 15])
for ii in range(len(wall_thickness)):
    radius = np.append(radius, radius[-1] + wall_thickness[ii])
numb_layers = len(radius)
print(radius)

surf1 = openmc.ZCylinder(surface_id=1, x0=0, y0=0, r=radius[0], name='surf 1')
surf2 = openmc.ZCylinder(surface_id=2, x0=0, y0=0, r=radius[1], name='surf 2')  # copper
surf3 = openmc.ZCylinder(surface_id=3, x0=0, y0=0, r=radius[2], name='surf 3')  # silver

surf3.boundary_type = 'vacuum'
xplane = openmc.XPlane(x0=0, boundary_type='periodic')
yplane = openmc.YPlane(y0=0, boundary_type='periodic')
# plane = openmc.Plane(a=1.0, b=1.0, d=0, boundary_type='periodic', periodic_surface=yplane)
yplane.periodic_surface = xplane
min_z = openmc.ZPlane(z0=-5, boundary_type='reflective')
max_z = openmc.ZPlane(z0=5, boundary_type='reflective')
# Cells
cell_1 = openmc.Cell(cell_id=1, name='cell 1')
cell_2 = openmc.Cell(cell_id=2, name='cell 2')
cell_3 = openmc.Cell(cell_id=3, name='cell 3')
# Cell regions
cell_1.region = -surf1 & -max_z & +min_z & -xplane & +yplane  # plasma
cell_2.region = -surf2 & +surf1 & -max_z & +min_z & -xplane & +yplane  # copper
cell_3.region = -surf3 & +surf2 & -max_z & +min_z & -xplane & +yplane  # silver

# cell_1.region = -surf1 & -max_z & +min_z  # plasma
# cell_2.region = -surf2 & +surf1 & -max_z & +min_z  # copper
# cell_3.region = -surf3 & +surf2 & -max_z & +min_z  # silver
# Register Materials with Cells
cell_2.fill = copper  # scrape off layer
cell_3.fill = silver
# Universe
universe = openmc.Universe(name='universe', cells=[cell_1, cell_2, cell_3])
geometry = openmc.Geometry(universe)
geometry.export_to_xml()

# plot geometry
# universe.plot(width=[-700, 700], pixels=[100, 100])
# plt.show()

# Tallies
tallies_file = openmc.Tallies()
cell2_filter = openmc.CellFilter(cell_2)
energy_filter = openmc.EnergyFilter([0., 20.e6])
scores = ['flux']
c2_tally = openmc.Tally(tally_id=2, name='tally c2')
c2_tally.filters = [cell2_filter, energy_filter]
c2_tally.scores = scores
tallies_file.append(c2_tally)
# Instantiate a Tallies collection and export to XML
tallies_file.export_to_xml()
# settings
settings_file = openmc.Settings()
settings_file.run_mode = 'fixed source'
settings_file._survival_biasing = 'true'
# source definition
source = openmc.Source(strength=1.0e+19)
source.space = openmc.stats.Box((1e-6, 1e-6, -5), (1e-6, 1e-6, 5))
source.angle = openmc.stats.Isotropic()
source.energy = openmc.stats.Discrete([14.1e6], [1.0])
settings_file.source = source
settings_file.batches = 10
settings_file.generations_per_batch = 100
settings_file.particles = 500
settings_file.export_to_xml()
# run
openmc.run(threads=8)
# import section
sp = openmc.StatePoint('statepoint.10.h5')
c2_output_tally = sp.get_tally(scores=['flux'], filters=[cell2_filter])
df2 = c2_output_tally.get_pandas_dataframe()
pd.options.display.float_format = '{:.7e}'.format
print(df2)

Much appreciated the help.
Best regards,
Tony

6

When you use a rotationally-periodic boundary condition, the normal vectors for the periodic surfaces need to face inward to the geometry. This means if you are pairing an x-plane and y-plane, your geometry should be defined in the first quadrant (x>0 and y>0). This means that you should be using +xplane instead of -xplane in the following:

Once you make that change, your model should run appropriately.

7

Thanks a lot @paulromano, it works. It makes sense, although I have a concern: on the website, (https://docs.openmc.org/en/stable/usersguide/geometry.html) when boundary conditions are explained, there is a yellow box which invites people to pay particular attention when dealing with periodic B.C. It suggests to define the geometry in the first quadrant, i.e., above the y-plane and to the right of the x-plane. Would that imply to define the first quadrant as the one on the top right hand side rather than the top left hand side? If I remember correctly, I think I defined my geometry coherently with both how the cartesian plane is reproduced and the statement in the yellow box. Am I missing something?

Finally, the same approach works for a portion of the geometry smaller than a quadrant, let’s say, an octant, am I right?

Anyway, these are miscellaneous, so thanks a lot for the help you gave me so far, it would definitely help me to reduce the computational cost of my model, and I am truly grateful for that.

8

Yes, the usual convention is that the first quadrant is the top right hand side. In your original definition, you were using the region -xplane, and since xplane is equivalent to x=0, -xplane is the region for which x<0. +xplane will give you x>0 which is what is needed for the first quadrant.

One of our developers just implemented support for rotationally periodic boundary conditions other than x- and y-planes, so in the latest developmental version it is possible to do a 1/8-core model, or a 1/6 model (useful for hexagonal geometries). This will be part of the next release. The advice I gave before (make sure the normals for the surfaces point inward to the geometry) still stands for this new capability.

9

and that is exactly my concern: if I define the geometry on the top rhs, as we agree, and I plot it, I get a quarter of circumference there. But running the calculation I do get the error for which I opened this query as well. On the contrary, if I define the geometry on the top left hand side (which is fine for my purposes), and therefore with both the norm vectors of xplane and plane >0, it works. But it is on the top lhs, not the top rhs.

Now, mine is just a concern of how it is defined the cartesian plane on OpenMC, you helped me to solve my issue. I just want to understand why the norm vectors result positive on the top lhs quadrant rather than the top rhs. Simply this

10

The x-plane surface is defined in OpenMC as x - x0 = 0, and as such the surface normal is always pointing to the right (because the leading coefficient on the ‘x’ term is just 1). Note that the surface normal is defined for the surface, not the half-space that you’re using, so whether you have -xplane or +xplane, the normal is pointing to the right in both cases. If you really wanted an x-plane with a normal pointing to the left, you’d have to use a general plane (Ax + By + Cz - D = 0) with A=-1, B=0, C=0.

1 Like
11

Dear Paul,

Thanks a lot for the clarification and the patience,
Best regards,

Tony