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

Bug in Fixed Source Simulation (Calculation or Plotting) #177

Closed geogunow closed 8 years ago

geogunow commented 9 years ago

I have recently been working on some test cases which used a fixed cosine distribution for the neutron source. I have come across an error when an irregular lattice is used in conjunction with the fixed source solver and plotter. After the solver completes, I receive the following error messages when OpenMOC attempts to plot the flux distribution:

/usr/lib/python2.7/dist-packages/numpy/ma/core.py:3847: UserWarning: Warning: converting a masked element to nan.
  warnings.warn("Warning: converting a masked element to nan.")
/usr/lib/pymodules/python2.7/matplotlib/colorbar.py:581: RuntimeWarning: invalid value encountered in greater
  inrange = (ticks > -0.001) & (ticks < 1.001)
/usr/lib/pymodules/python2.7/matplotlib/colorbar.py:581: RuntimeWarning: invalid value encountered in less
  inrange = (ticks > -0.001) & (ticks < 1.001)
/usr/lib/pymodules/python2.7/matplotlib/colors.py:576: RuntimeWarning: invalid value encountered in less
  cbook._putmask(xa, xa < 0.0, -1)

Also, I commented the python code to print the fluxes in the plotting routine and it shows that the flux profile is entirely NaN's!

Important note: when I re-run the script, it runs just fine and produces the correct flux! It must be some kind of indexing error. My assumption was that this was due to loading the tracking file versus actually doing the segmentation but it turns out that is not the case. Even when I delete the tracking file it still runs correctly. The only way I know how to reproduce the error after it is initially run is to delete the tracks, log, and plots folder and re-run.

The input is given below:

import numpy
from openmoc import *
import openmoc.log as log
import openmoc.plotter as plotter
from openmoc.options import Options
from math import sin, cos, pi

###############################################################################
#######################   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, one directional'
    ' gradient...')

length = 100.0
num_cells_x = 10 
num_cells_y = 5

###############################################################################
#######################   Define Material Properties   ########################
###############################################################################

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...')

left = XPlane(x=-length/2, name='left')
right = XPlane(x=length/2, name='right')
top = YPlane(y=length/2, name='top')
bottom = YPlane(y=-length/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...')

lattice = Lattice(name='MxN lattice')
lattice.setWidth(width_x=length/num_cells_x, width_y=length/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,1)
#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)

print "begin setting source"
width_x = lattice.getWidthX()
width_y = lattice.getWidthY()
for i in range(num_cells_x):
  for j in range(num_cells_y):

    # calculate midpoint
    x_min = i * width_x - length/2
    x_max = x_min + width_x
    y_min = j * width_y - length/2
    y_max = y_min + width_y
    x_mid = (x_min + x_max) / 2
    y_mid = (y_min + y_max) / 2

    # extract cell from FSR ID
    coords = LocalCoords(x_mid, y_mid)
    coords.setUniverse(root_universe)
    geometry.findCellContainingCoords(coords)
    fsr_id = geometry.findFSRId(coords)

    # compute source based on cosine distribution
    x_comp = sin(x_max*pi/length) - sin(x_min*pi/length)
    y_comp = sin(y_max*pi/length) - sin(y_min*pi/length)
    source = x_comp * y_comp * length**2 / pi**2
    vol = (x_max - x_min) * (y_max - y_min)
    source /= vol

    # set source
    for g in range(basic_material.getNumEnergyGroups()):
      solver.setFixedSourceByFSR(fsr_id, g+1, source)

print "end setting source"

solver.computeFlux(max_iters)
solver.printTimerReport()

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

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

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

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

See #178 for the fix.

geogunow commented 9 years ago

178 fixed the issue in the script included above, but I found a new case that the updated code breaks on. When the number of cells in a given direction becomes exceedingly large, it seems I get a new issue. Code is given below. I am re-opening this issue.

import numpy
from openmoc import *
import openmoc.log as log
import openmoc.plotter as plotter
from openmoc.options import Options
from math import sin, cos, pi

###############################################################################
#######################   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, one directional'
    ' gradient...')

length = 100.0
num_cells_x = 10000
num_cells_y = 1

###############################################################################
#######################   Define Material Properties   ########################
###############################################################################

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...')

left = XPlane(x=-length/2, name='left')
right = XPlane(x=length/2, name='right')
top = YPlane(y=length/2, name='top')
bottom = YPlane(y=-length/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...')

lattice = Lattice(name='MxN lattice')
lattice.setWidth(width_x=length/num_cells_x, width_y=length/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,1)
#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)

print "begin setting source"
width_x = lattice.getWidthX()
width_y = lattice.getWidthY()
for i in range(num_cells_x):
  for j in range(num_cells_y):

    # calculate midpoint
    x_min = i * width_x - length/2
    x_max = x_min + width_x
    y_min = j * width_y - length/2
    y_max = y_min + width_y
    x_mid = (x_min + x_max) / 2
    y_mid = (y_min + y_max) / 2

    # extract cell from FSR ID
    coords = LocalCoords(x_mid, y_mid)
    coords.setUniverse(root_universe)
    geometry.findCellContainingCoords(coords)
    fsr_id = geometry.findFSRId(coords)

    # compute source based on cosine distribution
    x_comp = sin(x_max*pi/length) - sin(x_min*pi/length)
    y_comp = sin(y_max*pi/length) - sin(y_min*pi/length)
    source = x_comp * y_comp * length**2 / pi**2
    vol = (x_max - x_min) * (y_max - y_min)
    source /= vol

    # set source
    for g in range(basic_material.getNumEnergyGroups()):
      solver.setFixedSourceByFSR(fsr_id, g+1, source)

print "end setting source"

solver.computeFlux(max_iters)
solver.printTimerReport()

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

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

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

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

What issue do you come across? I've been unable to reproduce anything strange with the latest code on develop.

geogunow commented 9 years ago

So, I looked into this again and sometimes it runs fine. But other times I encounter the same error as before. My repository is up-to-date on the develop branch. And I just pulled & re-installed. The error persists. I get the following error:

/usr/lib/python2.7/dist-packages/numpy/ma/core.py:3847: UserWarning: Warning: converting a masked element to nan.
  warnings.warn("Warning: converting a masked element to nan.")
/usr/lib/pymodules/python2.7/matplotlib/colorbar.py:581: RuntimeWarning: invalid value encountered in greater
  inrange = (ticks > -0.001) & (ticks < 1.001)
/usr/lib/pymodules/python2.7/matplotlib/colorbar.py:581: RuntimeWarning: invalid value encountered in less
  inrange = (ticks > -0.001) & (ticks < 1.001)
/usr/lib/pymodules/python2.7/matplotlib/colors.py:576: RuntimeWarning: invalid value encountered in less
  cbook._putmask(xa, xa < 0.0, -1)

And it happens both when tracks are created or loaded from file.

wbinventor commented 9 years ago

I tried to reproduce this but am unable to do so. Of course that doesn't mean that it isn't an issue on some systems. But I think the onus may be on you @geogunow to dig this bug up. I imagine it is an initialization issue for one of the arrays in the Solver. On Jul 22, 2015 8:26 PM, "Geoffrey Gunow" notifications@github.com wrote:

So, I looked into this again and sometimes it runs fine. But other times I encounter the same error as before. My repository is up-to-date on the develop branch. And I just pulled & re-installed. The error persists.

— Reply to this email directly or view it on GitHub https://github.com/mit-crpg/OpenMOC/issues/177#issuecomment-123951313.

wbinventor commented 8 years ago

Based on the fixed source calculations I have been performing I have not observed anything like this issue, and believe it was resolved with changes made this past summer.