search

mit-crpg / OpenMOC

A method of characteristics code for nuclear reactor physics calculations.
https://mit-crpg.github.io/OpenMOC/
Other
148 stars 83 forks source link

Problem with Setting Boundary Conditions #155

Closed geogunow closed 9 years ago

geogunow commented 9 years ago

I have been writing a test case today to showcase an example where the user might want cells of unequal x and y lengths. That part seems to be working fine - OpenMOC seems to handle uneven mesh properly now. However, it does seem to have a problem with applying the correct boundary conditions. When I run my case with one vacuum boundary (the left to be specific), I would expect to see a pure one-directional gradient as for the flux. However, I instead see the flux has the profile shown below:

fsr-flux-group-1

I suspect this should be an easy to fix bug in the setting of boundary conditions. Here is the source code for this case:

import numpy
from openmoc import *
import openmoc.log as log
import openmoc.plotter as plotter
from openmoc.options import Options

###############################################################################
#######################   Main Simulation Parameters   ########################
###############################################################################

options = Options()

num_threads = options.getNumThreads()
track_spacing = options.getTrackSpacing()
num_azim = options.getNumAzimAngles()
tolerance = options.getTolerance()
max_iters = options.getMaxIterations()

log.set_log_level('NORMAL')

log.py_printf('TITLE', 'Simulating a one group homogeneous infinite medium...')
log.py_printf('HEADER', 'The reference keff = 1.43...')

###############################################################################
#######################   Main Simulation Parameters   ########################
###############################################################################

log.py_printf('NORMAL', 'Creating materials...')

basic_material = Material(name='1-group infinite medium')
basic_material.setNumEnergyGroups(1)
basic_material.setSigmaA(numpy.array([0.069389522]))
basic_material.setSigmaF(numpy.array([0.0414198575]))
basic_material.setNuSigmaF(numpy.array([0.0994076580]))
basic_material.setSigmaS(numpy.array([0.383259177]))
basic_material.setChi(numpy.array([1.0]))
basic_material.setSigmaT(numpy.array([0.452648699]))

###############################################################################
###########################   Creating Surfaces   #############################
###############################################################################

log.py_printf('NORMAL', 'Creating surfaces...')

L = 100.0
left = XPlane(x=-L/2, name='left')
right = XPlane(x=L/2, name='right')
top = YPlane(y=L/2, name='top')
bottom = YPlane(y=-L/2, name='bottom')

left.setBoundaryType(VACUUM)
right.setBoundaryType(REFLECTIVE)
top.setBoundaryType(REFLECTIVE)
bottom.setBoundaryType(REFLECTIVE)

###############################################################################
#############################   Creating Cells   ##############################
###############################################################################

log.py_printf('NORMAL', 'Creating cells...')

fill = Cell(name='fill')
fill.setFill(basic_material)

root_cell = Cell(name='root cell')
root_cell.addSurface(halfspace=+1, surface=left)
root_cell.addSurface(halfspace=-1, surface=right)
root_cell.addSurface(halfspace=+1, surface=bottom)
root_cell.addSurface(halfspace=-1, surface=top)

###############################################################################
###########################    Creating Universes   ###########################
###############################################################################

log.py_printf('NORMAL', 'Creating universes...')

fill_universe = Universe(name='homogeneous fill cell')
fill_universe.addCell(fill)

root_universe = Universe(name='root universe')
root_universe.addCell(root_cell)

###############################################################################
############################    Creating Lattices   ###########################
###############################################################################

log.py_printf('NORMAL', 'Creating 100 x 1 lattice...')

num_cells_x = 10
num_cells_y = 1
lattice = Lattice(name='MxN lattice')
lattice.setWidth(width_x=L/num_cells_x, width_y=L/num_cells_y)
lattice.setUniverses([[fill_universe] * num_cells_x]*num_cells_y)
root_cell.setFill(lattice)

###############################################################################
##########################   Creating the Geometry   ##########################
###############################################################################

log.py_printf('NORMAL', 'Creating geometry...')

geometry = Geometry()
geometry.setRootUniverse(root_universe)
geometry.initializeFlatSourceRegions()

###############################################################################
########################   Creating the TrackGenerator   ######################
###############################################################################

log.py_printf('NORMAL', 'Initializing the track generator...')

track_generator = TrackGenerator(geometry, num_azim, track_spacing)
track_generator.setNumThreads(num_threads)
track_generator.generateTracks()

###############################################################################
###########################   Running a Simulation   ##########################
###############################################################################

solver = CPUSolver(track_generator)
solver.setNumThreads(num_threads)
solver.setConvergenceThreshold(tolerance)
solver.computeEigenvalue(max_iters)
solver.printTimerReport()

###############################################################################
############################    Generating Plots   ############################
###############################################################################

log.py_printf('NORMAL', 'Plotting data...')

plotter.plot_materials(geometry, gridsize=500)
plotter.plot_cells(geometry, gridsize=500)
plotter.plot_flat_source_regions(geometry, gridsize=500)
plotter.plot_spatial_fluxes(solver, energy_groups=[1])

log.py_printf('TITLE', 'Finished')
wbinventor commented 9 years ago

@geogunow - the problem isn't with the boundary condition treatment in OpenMOC. Rather, the issue is that you were using a very coarse spatial discretization of the geometry. As a result, the flux gradient did not show up in your plot. Simply changing num_cells_x = 1000 from 10 in your script above will produce the following flux gradient:

fsr-flux-group-1

You'll note that this is ever so slightly asymmetric and tilted to the right. If you set L = 10.0 instead of 100 this effect of the vacuum BC will be even more noticeable as expected:

fsr-flux-group-1

geogunow commented 9 years ago

Ok, but neither of those plots are right either. It should be a pure gradient. As in the maximum should occur on the right hand side. With 1000 mesh cells I would hope that OpenMOC would be able to crush this very simple problem.

nathangibson14 commented 9 years ago

If you set all of the boundary conditions to reflective, you don't get a flat flux. So something is screwy in the file. I'm looking at it now.

nathangibson14 commented 9 years ago

Alright, I take that back. The matplotlib plots are just a little screwy looking. You do get a flat flux with all reflective boundary conditions.

After playing with @geogunow's file a bit and playing with @wbinventor's simple-lattice sample input, I'm convinced that setting one side's boundary condition affects both that side and the opposite side. So our boundary conditions aren't working as expected, and I'm reopening this issue.

geogunow commented 9 years ago

Update: when I turn on CMFD, somehow this problem goes away (or at least mostly). I'm not sure why. But here's the updated case I'm running (now 100 x 100) without CMFD:

import numpy
from openmoc import *
import openmoc.log as log
import openmoc.plotter as plotter
from openmoc.options import Options

###############################################################################
#######################   Main Simulation Parameters   ########################
###############################################################################

options = Options()

num_threads = options.getNumThreads()
track_spacing = 0.1
num_azim = 16
tolerance = options.getTolerance()
max_iters = options.getMaxIterations()

log.set_log_level('NORMAL')

log.py_printf('TITLE', 'Simulating a one group homogeneous infinite medium...')
log.py_printf('HEADER', 'The reference keff = 1.387...')

###############################################################################
#######################   Main Simulation Parameters   ########################
###############################################################################

log.py_printf('NORMAL', 'Creating materials...')

basic_material = Material(name='1-group infinite medium')
basic_material.setNumEnergyGroups(1)
basic_material.setSigmaA(numpy.array([0.069389522]))
basic_material.setSigmaF(numpy.array([0.0414198575]))
basic_material.setNuSigmaF(numpy.array([0.0994076580]))
basic_material.setSigmaS(numpy.array([0.383259177]))
basic_material.setChi(numpy.array([1.0]))
basic_material.setSigmaT(numpy.array([0.452648699]))

###############################################################################
###########################   Creating Surfaces   #############################
###############################################################################

log.py_printf('NORMAL', 'Creating surfaces...')

L = 100.0
left = XPlane(x=-L/2, name='left')
right = XPlane(x=L/2, name='right')
top = YPlane(y=L/2, name='top')
bottom = YPlane(y=-L/2, name='bottom')

left.setBoundaryType(VACUUM)
right.setBoundaryType(REFLECTIVE)
top.setBoundaryType(REFLECTIVE)
bottom.setBoundaryType(REFLECTIVE)

###############################################################################
#############################   Creating Cells   ##############################
###############################################################################

log.py_printf('NORMAL', 'Creating cells...')

fill = Cell(name='fill')
fill.setFill(basic_material)

root_cell = Cell(name='root cell')
root_cell.addSurface(halfspace=+1, surface=left)
root_cell.addSurface(halfspace=-1, surface=right)
root_cell.addSurface(halfspace=+1, surface=bottom)
root_cell.addSurface(halfspace=-1, surface=top)

###############################################################################
###########################    Creating Universes   ###########################
###############################################################################

log.py_printf('NORMAL', 'Creating universes...')

fill_universe = Universe(name='homogeneous fill cell')
fill_universe.addCell(fill)

root_universe = Universe(name='root universe')
root_universe.addCell(root_cell)

###############################################################################
############################    Creating Lattices   ###########################
###############################################################################

log.py_printf('NORMAL', 'Creating 100 x 1 lattice...')

num_cells_x = 100
num_cells_y = 100
lattice = Lattice(name='MxN lattice')
lattice.setWidth(width_x=L/num_cells_x, width_y=L/num_cells_y)
lattice.setUniverses([[fill_universe] * num_cells_x]*num_cells_y)
root_cell.setFill(lattice)

###############################################################################
##########################   Creating the Geometry   ##########################
###############################################################################

log.py_printf('NORMAL', 'Creating geometry...')

geometry = Geometry()
geometry.setRootUniverse(root_universe)
geometry.initializeFlatSourceRegions()

###############################################################################
########################   Creating the TrackGenerator   ######################
###############################################################################

log.py_printf('NORMAL', 'Initializing the track generator...')

track_generator = TrackGenerator(geometry, num_azim, track_spacing)
track_generator.setNumThreads(num_threads)
track_generator.generateTracks()

###############################################################################
###########################   Running a Simulation   ##########################
###############################################################################

solver = CPUSolver(track_generator)
solver.setNumThreads(num_threads)
solver.setConvergenceThreshold(tolerance)
solver.computeEigenvalue(max_iters)
solver.printTimerReport()

###############################################################################
############################    Generating Plots   ############################
###############################################################################

log.py_printf('NORMAL', 'Plotting data...')

plotter.plot_materials(geometry, gridsize=10)
plotter.plot_cells(geometry, gridsize=10)
plotter.plot_flat_source_regions(geometry, gridsize=10)
plotter.plot_spatial_fluxes(solver, energy_groups=[1])

log.py_printf('TITLE', 'Finished')

which produces the following flux plot:

fsr-flux-group-1-no-cmfd

which is obviously incorrect (should be a pure gradient). Now I turn on CMFD using this code:

import numpy
from openmoc import *
import openmoc.log as log
import openmoc.plotter as plotter
from openmoc.options import Options

###############################################################################
#######################   Main Simulation Parameters   ########################
###############################################################################

options = Options()

num_threads = options.getNumThreads()
track_spacing = 0.1
num_azim = 16
tolerance = options.getTolerance()
max_iters = options.getMaxIterations()

log.set_log_level('NORMAL')

log.py_printf('TITLE', 'Simulating a one group homogeneous infinite medium...')
log.py_printf('HEADER', 'The reference keff = 1.387...')

###############################################################################
#######################   Main Simulation Parameters   ########################
###############################################################################

log.py_printf('NORMAL', 'Creating materials...')

basic_material = Material(name='1-group infinite medium')
basic_material.setNumEnergyGroups(1)
basic_material.setSigmaA(numpy.array([0.069389522]))
basic_material.setSigmaF(numpy.array([0.0414198575]))
basic_material.setNuSigmaF(numpy.array([0.0994076580]))
basic_material.setSigmaS(numpy.array([0.383259177]))
basic_material.setChi(numpy.array([1.0]))
basic_material.setSigmaT(numpy.array([0.452648699]))

###############################################################################
###########################   Creating Surfaces   #############################
###############################################################################

log.py_printf('NORMAL', 'Creating surfaces...')

L = 100.0
left = XPlane(x=-L/2, name='left')
right = XPlane(x=L/2, name='right')
top = YPlane(y=L/2, name='top')
bottom = YPlane(y=-L/2, name='bottom')

left.setBoundaryType(VACUUM)
right.setBoundaryType(REFLECTIVE)
top.setBoundaryType(REFLECTIVE)
bottom.setBoundaryType(REFLECTIVE)

###############################################################################
#############################   Creating Cells   ##############################
###############################################################################

log.py_printf('NORMAL', 'Creating cells...')

fill = Cell(name='fill')
fill.setFill(basic_material)

root_cell = Cell(name='root cell')
root_cell.addSurface(halfspace=+1, surface=left)
root_cell.addSurface(halfspace=-1, surface=right)
root_cell.addSurface(halfspace=+1, surface=bottom)
root_cell.addSurface(halfspace=-1, surface=top)

###############################################################################
###########################    Creating Universes   ###########################
###############################################################################

log.py_printf('NORMAL', 'Creating universes...')

fill_universe = Universe(name='homogeneous fill cell')
fill_universe.addCell(fill)

root_universe = Universe(name='root universe')
root_universe.addCell(root_cell)

###############################################################################
############################    Creating Lattices   ###########################
###############################################################################

log.py_printf('NORMAL', 'Creating 100 x 1 lattice...')

num_cells_x = 100
num_cells_y = 100
lattice = Lattice(name='MxN lattice')
lattice.setWidth(width_x=L/num_cells_x, width_y=L/num_cells_y)
lattice.setUniverses([[fill_universe] * num_cells_x]*num_cells_y)
root_cell.setFill(lattice)

###############################################################################
###########################   Creating CMFD Mesh    ###########################
###############################################################################

log.py_printf('NORMAL', 'Creating Cmfd mesh...')

cmfd = Cmfd()
cmfd.setMOCRelaxationFactor(0.6)
cmfd.setSORRelaxationFactor(1.5)
cmfd.setLatticeStructure(51,51)
cmfd.setGroupStructure([1])

###############################################################################
##########################   Creating the Geometry   ##########################
###############################################################################

log.py_printf('NORMAL', 'Creating geometry...')

geometry = Geometry()
geometry.setRootUniverse(root_universe)
geometry.setCmfd(cmfd)
geometry.initializeFlatSourceRegions()

###############################################################################
########################   Creating the TrackGenerator   ######################
###############################################################################

log.py_printf('NORMAL', 'Initializing the track generator...')

track_generator = TrackGenerator(geometry, num_azim, track_spacing)
track_generator.setNumThreads(num_threads)
track_generator.generateTracks()

###############################################################################
###########################   Running a Simulation   ##########################
###############################################################################

solver = CPUSolver(track_generator)
solver.setNumThreads(num_threads)
solver.setConvergenceThreshold(tolerance)
solver.computeEigenvalue(max_iters)
solver.printTimerReport()

###############################################################################
############################    Generating Plots   ############################
###############################################################################

log.py_printf('NORMAL', 'Plotting data...')

plotter.plot_materials(geometry, gridsize=10)
plotter.plot_cells(geometry, gridsize=10)
plotter.plot_flat_source_regions(geometry, gridsize=10)
plotter.plot_spatial_fluxes(solver, energy_groups=[1])

log.py_printf('TITLE', 'Finished')

and I get:

fsr-flux-group-1

which looks a lot better! Still not sure why this is the case. Though there does seems to still be something a bit weird going on near the right side.