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.