Dagmc.h5m set boundary

qianni · November 18, 2022, 2:35am

hi,I I wanted to partition the geometry by a tetrahedral grid and then perform the fuel depletion calculation, so I created a geometry and drew the tetrahedral grid, exported its grid nodes, and converted it into a h5m file by the vertices to h5m plug-in. Each small volume has defined the geometry but how to give the boundaries of each volume grid。when i run the model,it appare error:No boundary conditions were applied to any surfaces!

import numpy as np
import re
import os
from vertices_to_h5m import vertices_to_h5m
from itertools import combinations
import linecache
from math import log10
import openmc
import math
from matplotlib import pyplot
import dagmc_h5m_file_inspector as di
import openmc.lib
import openmc.lib
import copy
from time import sleep
from PySide2 import QtCore, QtGui
import openmc.model

def cell():
model = openmc.model.Model()

with open('twist.txt', "r", encoding="UTF-8") as r:
	 alltxt=r.read()	
	 alltxt=alltxt.replace('%', '')    
a = alltxt.find('Coordinates')
b = alltxt.find('Elements (tetrahedra)')
c=alltxt[a:b]
d=alltxt[b:]   


f1 = open("./vertices.txt", "w+")
f1.write(c)
f1.close()
f2 = open("./triangle.txt", "w+")
f2.write(d)
f2.close()

vertices = np.loadtxt('vertices.txt', skiprows=1, dtype=float)
#print(vertices)

triangles1=np.loadtxt('triangle.txt', skiprows=1, dtype=int)

#print(triangles1)
triangles2=[]

count = len(open("triangle.txt").readlines())

for data in triangles1:
    c = list(combinations(data,3))# pai lie  zu he 
    #print(c)
    cc=np.array(c)-1
    triangles2.append(cc)  
#print(triangles2)    


count = len(open("triangle.txt").readlines())-1
print("count=",count)

mat_list=[]
fuel_material  = openmc.Material(material_id=0)
fuel_material.set_density('g/cm3', 9.7012)
fuel_material.add_nuclide('Zr90', 0.255449, 'wo')
fuel_material.add_nuclide('Zr91', 0.055707, 'wo')
fuel_material.add_nuclide('Zr92', 0.085150, 'wo')
fuel_material.add_nuclide('Zr94', 0.086292, 'wo')
fuel_material.add_nuclide('Zr96', 0.013902, 'wo')
fuel_material.add_nuclide('U235', 0.097811, 'wo')
fuel_material.temperature=800
fuel_material.name= "mat0"
mat_list.append(fuel_material)

for i in range(1,count):
	mat = fuel_material.clone()
	mat.name= "mat"+str(i)
	mat_list.append(mat)



# Instantiate a Materials collection and export to XML
model.material = openmc.Materials(mat_list)



mats=[]
for i in range(0,count ):
	mats.append( "mat"+str(i))
vertices_to_h5m(
    vertices= vertices,
    triangles=triangles2,
    material_tags=mats,
    h5m_filename="dagmc.h5m",
 )


dagunv=openmc.DAGMCUniverse(filename="dagmc.h5m")
fuel_pin = openmc.Cell(name='Fuel', fill=dagunv)

model.geometry.root_universe = openmc.DAGMCUniverse(filename="dagmc.h5m")
model.geometry.root_universe.add_cells(fuel_pin)


print(di.get_volumes_from_h5m("dagmc.h5m"))

print(di.get_materials_from_h5m("dagmc.h5m"))










###############################################################################
#                   Exporting to OpenMC plots.xml file
###############################################################################

plot_xy = openmc.Plot(plot_id=1)
plot_xy.filename = 'plot_xy'
plot_xy.origin = [0, 0, 0]
plot_xy.basis = "xy"
plot_xy.width = [1.3, 1.3]
plot_xy.pixels = [400, 400]
plot_xy.color_by = 'material'

plot_xz = openmc.Plot(plot_id=2)
plot_xz.filename = 'plot_xz'
plot_xz.origin = [0, 0, 0]
plot_xz.basis = "xz"
plot_xz.width = [1.3, 0.2]
plot_xz.pixels = [400, 400]
plot_xz.color_by = 'material'


# Instantiate a Plots object and export to XML
model.plots= openmc.Plots([plot_xy, plot_xz])



###############################################################################
#                   Exporting to OpenMC settings.xml file
###############################################################################

# Instantiate a Settings object, set all runtime parameters, and export to XML

model.settings.particles = 100
model.settings.batches = 10
model.settings.inactive = 2
model.settings.output = {'summary' : False}


return model

model=openmc.cell()
surfaces = model.geometry.get_all_surfaces()
# modify the boundary condition of the
# planar surfaces bounding the assembly
for surface in surfaces.values():
	surface.boundary_type = 'reflective'

model.run()
                                                                                         Sincerely hope to get your answer

this is the plot of dagmc.h5m

Jiankai · November 18, 2022, 3:30pm

I would like to add a cell with boundary condition that can contain the dagmc universe.
For example:

container = openmc.ZCylinder(r=1.0, boundary_type="reflective")
dagunv=openmc.DAGMCUniverse(filename="dagmc.h5m")
fuel_pin = openmc.Cell(region=-container, fill=dagunv)

Shimwell · November 18, 2022, 8:53pm

Great to see a few package being useful. I never know if anyone uses these things I throw out on github but great to see vertices-to-h5m and dagmc-file-inspector being useful.

I think Jiankai has solved it but I just wanted to add that there is a auto bounding option as well.

Traditionally DAGMC geometry was made with a hollow cube surrounding the geometry which had the material tag graveyard assigned to it. This acted as the vacuum boundary and killed particles when they got to this outer volume. This functionality is still supported so you could make another set of triangles to encompass your geometry.

More recently openmc was able to contain a graveyard free dagmc geometry in a cage surface like Jiankai example above. This makes visualization easier as there is no large encompassing dagmc volume.

Really recently we added this pull request which simplifies things further

github.com/openmc-dev/openmc

added DAGMCUniverse.bounding_region

openmc-dev:develop ← shimwell:adding_dagmcuniverse_bounding_region

opened 06:16PM - 13 Jul 22 UTC

shimwell

+107 -9

As discussed with @pshriwise in issue #2104 and building on the work from PR #21…11 this is another small PR that adds the ability to automatically generate a bounding region around a DAGMC geometry. This adds the ```DAGMCUniverse.bounding_region()``` method This would be useful for creating a spherical or box shaped vacuum surface to bound the DAGMC geometry or other uses (boundary type is user controllable) This is based on some code from another micro package I made a while back :recycle: [here](https://github.com/fusion-energy/openmc-dagmc-wrapper/blob/9c4485d4f40e7df2bac5333be68fe6d4e0caa570/openmc_dagmc_wrapper/Geometry.py#L114-L153) and [here](https://github.com/fusion-energy/openmc-dagmc-wrapper/blob/9c4485d4f40e7df2bac5333be68fe6d4e0caa570/openmc_dagmc_wrapper/Geometry.py#L50-L55) I've added some simple tests to the DAGMC universe test and type checking Current method of bringing a graveyard free DAGMC geometry file into to OpenMC ```python # makes use of the dagmc geometry dag_univ = openmc.DAGMCUniverse("dagmc.h5m") # creates an edge of universe boundary at a large radius vac_surf = openmc.Sphere(r=10000, surface_id=9999, boundary_type="vacuum") # specifies the region as below the universe boundary region = -vac_surf # creates a cell from the region and fills the cell with the dagmc geometry containing_cell = openmc.Cell(cell_id=9999, region=region, fill=dag_univ) model = openmc.Model() model.geometry = openmc.Geometry(root=[containing_cell]) ``` method possible with this PR of bringing a DAGMC geometry file into to OpenMC ```python # loads up dagmc file and automatically bounds the geometry with a vacuum surface bound_dag_univ= openmc.DAGMCUniverse("dagmc.h5m").bounded_universe() model = openmc.Model() model.geometry = openmc.Geometry(root=bound_dag_univ) ```

This allows combined universe to be made in one line

dagunv = openmc.DAGMCUniverse("dagmc.h5m").bounded_universe()

qianni · November 27, 2022, 1:21pm

Thank you very much for your answer. I added the boundary to the geometry according to your method.

I changed a cylindrical rod bundle to carry out comprehensive calculation. vertice_to_h5m method was still adopted, and the geometry was divided into 5508 regions, each region was provided with material and volume, which was exported as dagmc file. tally set by materialfilter I wanted to see the calculation in three dimensions, but I didn’t know how to do Post processing.

from itertools import combinations
import numpy as np
import re
import os
from vertices_to_h5m import vertices_to_h5m
import linecache
from math import log10
import openmc
import math
from matplotlib import pyplot
import dagmc_h5m_file_inspector as di
import openmc.lib
import copy
from time import sleep
from PySide2 import QtCore, QtGui

def volumecalculate(k):
npoint_a=elementsList[k][0]
npoint_b = elementsList2[k][1]
npoint_c = elementsList2[k][2]
npoint_d = elementsList2[k][3]
coor_a = verticesList2[npoint_a]
coor_b = verticesList2[npoint_b]
coor_c = verticesList2[npoint_c]
coor_d = verticesList2[npoint_d]
vec_ad = coor_a - coor_d
vec_bd = coor_b - coor_d
vec_cd = coor_c - coor_d
volume=np.abs(np.dot(vec_cd,np.cross(vec_ad,vec_bd)))/6
return volume

numberOfVertices = ‘number of mesh vertices’
coordOfmesh = ‘Mesh vertex coordinates’
numberOfElements = ‘number of elements’
connectOfElement = ‘Elements’
indicesOfGeometry = ‘Geometric entity indices’
totalNumberOfVertices = 0
lineCount = 0
verticesList =
elementsList =
materialList =
trianglesList =
lineCountOfCoord = 1000
lineCountOfElement = 1000000
lineCountOfGeom = 2000000

for line in open (‘./Untitled.mphtxt’, ‘rt’):
lineCount = lineCount + 1
if line.find(numberOfVertices) > 0:
comment = line.split(’ # ‘)
totalNumberOfVertices = int(comment[0])
elif line.find(coordOfmesh) > 0:
lineCountOfCoord = lineCount
elif (lineCountOfCoord < lineCount <= lineCountOfCoord+totalNumberOfVertices):
line = line.strip().split(’ ‘)
verticesList.append([float(x) for x in line]) # vertices’ coordinates obtained
elif line.find(numberOfElements) > 0:
comment = line.split(’ # ‘)
totalNumberOfElements = int(comment[0])
elif line.find(connectOfElement) > 0:
lineCountOfElement = lineCount
elif (lineCountOfElement < lineCount <= lineCountOfElement+totalNumberOfElements):
line = line.strip().split(’ ‘)
elementsList.append([int(x) for x in line]) # element connnectivity obtained
elif line.find(indicesOfGeometry) > 0:
lineCountOfGeom = lineCount
elif (lineCountOfGeom < lineCount <= lineCountOfGeom+totalNumberOfElements):
line = line.strip().split(’ ')
materialList.append([int(x) for x in line]) # domain which elements belong to obtained
else:
continue
#print(verticesList)
#print(elementsList)
#print(materialList)
#print(len(verticesList))
#print(len(elementsList))
#print(len(materialList))

Each element contains four triangles

#dets1 = np.array(verticesList)
#np.savetxt(“./vertices.txt”, dets1,fmt=‘%f’,delimiter=’ ')

#dets2 = np.array(elementsList)
#np.savetxt(“./triangle.txt”, dets2,fmt=‘%d’,delimiter=’ ')

#vertices = np.loadtxt(‘vertices.txt’, dtype=float)
#print(vertices)
#triangles1=np.loadtxt(‘triangle.txt’,dtype=int)

#triangles2=
#count = len(open(“triangle.txt”).readlines())

for data in elementsList:
c = list(combinations(data,3))
cc=np.array(c)
trianglesList.append(cc)
#print(trianglesList)

verticesList2=np.array(verticesList)
elementsList2=np.array(elementsList)

#print(len(trianglesList))
count=len(trianglesList)
print(count)

import random
doman1_id=
doman2_id=
doman3_id=
doman4_id=
doman5_id=
doman6_id=
doman7_id=
doman8_id=
doman9_id=
doman10_id=
doman11_id=
doman12_id=
doman13_id=
doman14_id=
doman15_id=
doman16_id=
for i in range(len(materialList)):
if(materialList[i][0]==1):
doman1_id.append(i)
if(materialList[i][0]==2):
doman2_id.append(i)
if(materialList[i][0]==3):
doman3_id.append(i)
if(materialList[i][0]==4):
doman4_id.append(i)
if(materialList[i][0]==5):
doman5_id.append(i)
if(materialList[i][0]==6):
doman6_id.append(i)
if(materialList[i][0]==7):
doman7_id.append(i)
if(materialList[i][0]==8):
doman8_id.append(i)
if(materialList[i][0]==9):
doman9_id.append(i)
if(materialList[i][0]==10):
doman10_id.append(i)
if(materialList[i][0]==11):
doman11_id.append(i)
if(materialList[i][0]==12):
doman12_id.append(i)
if(materialList[i][0]==13):
doman13_id.append(i)
if(materialList[i][0]==14):
doman14_id.append(i)
if(materialList[i][0]==15):
doman15_id.append(i)
if(materialList[i][0]==16):
doman16_id.append(i)
#print(doman1_id)

#print(len(doman16_id))

uo2 = openmc.Material(name=‘UO2 fuel at 2.4% wt enrichment’)
uo2.set_density(‘g/cm3’, 10.29769)
uo2.add_element(‘U’, 1., enrichment=4.4)
uo2.add_element(‘O’, 2.)

mat_list =
mats=
for i in range(len(doman7_id)):
mat = uo2.clone()
mat.name= “mat”+str(doman7_id[i])
mat.volume= volumecalculate(i)
mat_list.append(mat)
mats.append(“mat”+str(doman7_id[i]))

for i in range(len(doman8_id)):
mat = uo2.clone()
mat.name= “mat”+str(doman8_id[i])
mat.volume= volumecalculate(i)
mat_list.append(mat)
mats.append( “mat”+str(doman8_id[i]))

for i in range(len(doman12_id)):
mat = uo2.clone()
mat.name= “mat”+str(doman12_id[i])
mat.volume= volumecalculate(i)
mat_list.append(mat)
mats.append( “mat”+str(doman12_id[i]))

for i in range(len(doman13_id)):
mat = uo2.clone()
mat.name= “mat”+str(doman13_id[i])
mat.volume= volumecalculate(i)
mat_list.append(mat)
mats.append( “mat”+str(doman13_id[i]))

helium = openmc.Material(name=‘Helium for gap’)
helium.set_density(‘g/cm3’, 0.001598)
helium.add_element(‘He’, 2.4044e-4)

for i in range(len(doman5_id)):
mat = helium .clone()
mat.name= “mat”+str(doman5_id[i])
mat.volume= volumecalculate(i)
#print(mat.volume)
mat_list.append(mat)
mats.append( “mat”+str(doman5_id[i]))

for i in range(len(doman6_id)):
mat = helium .clone()
mat.name= “mat”+str(doman6_id[i])
mat.volume= volumecalculate(i)
mat_list.append(mat)
mats.append( “mat”+str(doman6_id[i]))

for i in range(len(doman11_id)):
mat = helium .clone()
mat.name= “mat”+str(doman11_id[i])
mat.volume= volumecalculate(i)
mat_list.append(mat)
mats.append( “mat”+str(doman11_id[i]))

for i in range(len(doman14_id)):
mat = helium .clone()
mat.name= “mat”+str(doman14_id[i])
mat.volume= volumecalculate(i)
mat_list.append(mat)
mats.append( “mat”+str(doman14_id[i]))

zircaloy = openmc.Material(name=‘Zircaloy 4’)
zircaloy.set_density(‘g/cm3’, 6.55)
zircaloy.add_element(‘Sn’, 0.014, ‘wo’)
zircaloy.add_element(‘Fe’, 0.00165, ‘wo’)
zircaloy.add_element(‘Cr’, 0.001, ‘wo’)
zircaloy.add_element(‘Zr’, 0.98335, ‘wo’)

for i in range(len(doman3_id)):
mat = zircaloy.clone()
mat.name= “mat”+str(doman3_id[i])
mat.volume= volumecalculate(i)
mat_list.append(mat)
mats.append( “mat”+str(doman3_id[i]))

for i in range(len(doman4_id)):
mat = zircaloy.clone()
mat.name= “mat”+str(doman4_id[i])
mat.volume= volumecalculate(i)
mat_list.append(mat)
mats.append( “mat”+str(doman4_id[i]))

for i in range(len(doman10_id)):
mat = zircaloy.clone()
mat.name= “mat”+str(doman10_id[i])
mat.volume= volumecalculate(i)
mat_list.append(mat)
mats.append( “mat”+str(doman10_id[i]))

for i in range(len(doman15_id)):
mat = zircaloy.clone()
mat.name= “mat”+str(doman15_id[i])
mat.volume= volumecalculate(i)
mat_list.append(mat)
mats.append( “mat”+str(doman15_id[i]))

borated_water = openmc.Material(name=‘Borated water’)
borated_water.set_density(‘g/cm3’, openmc.data.water_density(temperature=600, pressure=15.5132))
borated_water.add_element(‘H’, 2.0)
borated_water.add_element(‘O’, 1.0)
borated_water.temperature = 600
borated_water.add_s_alpha_beta(‘c_H_in_H2O’)

elementsList2=np.array(elementsList)
verticesList2=np.array(verticesList)

for i in range(len(doman1_id)):
mat = borated_water.clone()
mat.name= “mat”+str(doman1_id[i])
mat.volume= volumecalculate(i)
mat_list.append(mat)
mats.append( “mat”+str(doman1_id[i]))

for i in range(len(doman2_id)):
mat =borated_water.clone()
mat.name= “mat”+str(doman2_id[i])
mat.volume= volumecalculate(i)
mat_list.append(mat)
mats.append( “mat”+str(doman2_id[i]))

for i in range(len(doman9_id)):
mat =borated_water.clone()
mat.name= “mat”+str(doman9_id[i])
mat.volume= volumecalculate(i)
mat_list.append(mat)
mats.append( “mat”+str(doman9_id[i]))

for i in range(len(doman16_id)):
mat = borated_water.clone()
mat.name= “mat”+str(doman16_id[i])
mat.volume= volumecalculate(i)
mat_list.append(mat)
mats.append( “mat”+str(doman16_id[i]))

materials = openmc.Materials(mat_list)
materials.export_to_xml()
#print(len(mats))

vertices_to_h5m(
vertices= verticesList,
triangles=trianglesList,
material_tags=mats,
h5m_filename=“dagmc.h5m”,
)

top = openmc.ZPlane( z0=+0.25, name=‘top’, boundary_type = ‘reflective’)
bottom = openmc.ZPlane( z0=-0.25, name=‘bottom’, boundary_type = ‘reflective’)
left = openmc.XPlane( x0=-6.3, name=‘left’, boundary_type = ‘reflective’)
right = openmc.XPlane( x0=+6.3, name=‘right’, boundary_type = ‘reflective’)
up = openmc.YPlane( y0=+6.3, name=‘up’, boundary_type = ‘reflective’)
down = openmc.YPlane( y0=-6.3, name=‘down’, boundary_type = ‘reflective’)
container = -top& +bottom & -up & +down & +left & -right

dagunv=openmc.DAGMCUniverse(filename=“dagmc.h5m”,auto_geom_ids=True)
fuel_pin = openmc.Cell(region=container , fill=dagunv)

Geometry = openmc.Geometry([fuel_pin])
Geometry.export_to_xml()

mat_filter = openmc.MaterialFilter(mat_list)
tallydep = openmc.Tally(tally_id=1, name=‘Spatial Tally 1’)
tallydep.filters = [mat_filter]
tallydep.scores = [‘fission-q-recoverable’, ‘heating’, ‘heating-local’, ‘fission’,‘flux’]
tallies = openmc.Tallies([tallydep])

tallies.export_to_xml()

Question 2: I used openmc-plotter to establish the geometry as shown in the following figure, some areas are white, is there a problem with the geometry? But it can be calculated normally.

qianni · November 27, 2022, 1:22pm

Thank you very much for your answer. I learned a lot. Very good

qianni · November 27, 2022, 1:31pm

As I had already divided regions when building the geometry, I thought it was not very good to set up libmesh Unstructured Mesh for statistical results, so I gave up the cellfilter and adopted the materialfilter, which may not be good. Maybe you have a better suggestion

Shimwell · November 28, 2022, 10:01am

Just to check my understanding. It sounds like you have changed from unstructured mesh to structured meshes? I think the unstructured libmesh route is perhaps the most flexible as it allows the mesh to follow volume boundaries and can also accept mixed (hex, tet etc) mesh elements. One disadvantage of the unstructured mesh is that simulations tend to be a bit slower than structured meshes.

We do have a nice selection of structured meshes (cylindrical, spherical, regular, rectilinear) and it could be possible that one of these suits your geometry well. Then you would have a nice fast simulation and a nice overlay of mesh voxels on the geometry.

Shimwell · November 28, 2022, 10:35am

@pshriwise is the best person to comment on this plot of unstructured mesh using the openmc-plotter.

‘’’ Question 2: I used openmc-plotter to establish the geometry as shown in the following figure, some areas are white, is there a problem with the geometry? But it can be calculated normally.

’‘’

pshriwise · November 28, 2022, 4:44pm

Hi @qianni!

Thanks for posting your question in the forum and thanks for using DAGMC with OpenMC – it’s a topic near and dear to me

For your first question above about post-processing, it’s a little tricky. We have methods for exporting data on a mesh in OpenMC, but your model is using a DAGMC geometry rather than a superimposed mesh. So we need to find a way to take the cell-based information and apply it to the tetrahedral mesh you’re representing in the file.

This is a little convoluted, but if you can convert your initial mesh format (.mphtxt it looks like) to VTK or an Exodus file, then you could use that as a mesh tally object in OpenMC to get a Python object with the vertices and connectivity post-simulation. At this point you can apply the cell-based data by using the openmc.UnstructuredMesh.write_data_to_vtk method to write a VTK file with the element-based data applied. So in this case you would have the following steps:

Generate the DAGMC geometry w/ vertices_to_h5m
Convert the .mphtxt file to Exodus or VTK and
Apply the tetrahedral mesh (the new Exodus/VTK file) as a mesh tally in OpenMC (only need to do this once to generate a statepoint file with the mesh information present)
Perform your depletion calculation with the tet-based DAGMC geometry
Load the unstructured mesh from the statepoint file resulting from step 3
Extract and apply cell-based information to the mesh loaded from step 5 using the write_data_to_vtk method.
Visualize with Paraview or Visit

Again, this assumes you can convert the original mesh into one of those formats.

Sorry for the complexity of this. FWIW, we’re working on an unstructured mesh tracking capability, in which case you’d simply be able to convert your original file to Exodus/VTK and run a problem directly on the tets.

Along the lines of @Shimwell’s comment on using one of the structured mesh type for tallying info, I’d be curious to know more about what information you’d like to gather as well. Are you looking to perform a depletion calculation for each tetrahedral element in the fuel region?

Interesting stuff. I hope this was helpful to you!

-Patrick

qianni · November 29, 2022, 7:38am

I think it is only a way to statistical the results whether the structured mesh or the unstructured mesh. The reason why I divide tetrahedron is that I think it is necessary to divide it into a serious of region to calculate the fuel depletion when taking a rod as the object. Therefore, tetrahedron is divided into regions, and each region counts the calculation results. Of course, I’m not sure the idea of divide many regions to depletion calculation is right

qianni · November 29, 2022, 7:39am

I hope you can give me your suggestions

pshriwise · December 5, 2022, 7:46pm

Hi @qianni,

For depletion calculations with cylindrical fuel rods, the rod is usually broken up into radial, angular, and axial regions (see pic below).
Screenshot 2022-12-05 at 1.44.25 PM

This can be done using CSG surfaces without a tet mesh fairly easily, and it will represent the fuel region without the approximation of linear surfaces near the boundary that you get with linear tetrahedra.

However, for fuel regions with irregular or complex shapes, I can see where a conformal tet mesh might be useful. Another benefit is its generality – if you can mesh the region and run it through your workflow, then you don’t have to worry about the fuel shape.

The main takeaway is that it depends on what you want to do, but, if you’re working on a problem with cylindrical fuel rods, you might want to consider doing the discretization manually.

I hope this gives some clarity and is helpful for whatever simulation you’re performing!

-Patrick

qianni · December 8, 2022, 1:25am

Thank you for your answer, which is very helpful to solve my problem .