Godiva reactor problem

User Support
1

Hi,
I am trying to compute the Eigenvalue for a U235 sphere inside a given water sphere with the “search_for_keff” function.
My problem is that I got it working once with a bounding box out of 6 planes and “reflective” boundaries respectively.
Every time I use “vacuum” or “transmissive” boundaries I get a “f(a) and f(b) must have different signs” error. This problem persists also with the spherical surrounding water sphere(even with reflective boundaries).
Below is my source code for clarification. Thanks in advance for the interest or support.

import matplotlib.pyplot as plt
import numpy as np
import openmc
import openmc.model

def build_model(r_crit):
  fuel = openmc.Material(name='U235')
  fuel.set_density('g/cm3', 10.31341)
  fuel.add_nuclide('U-235', 1.)

  water = openmc.Material(name='water')
  water.set_density('g/cm3', 1.)
  water.add_element('H', 2.)
  water.add_element('O', 1.)

  materials = openmc.Materials([fuel, water])

  fuel_radius = openmc.Sphere(x0=0., y0=0., r=r_crit, boundary_type='reflective')
  moderator_radius = openmc.Sphere(x0=0., y0=0., r=15, boundary_type='reflective')

  fuel_cell = openmc.Cell(name = 'U235 fuel')
  fuel_cell.fill = fuel
  fuel_cell.region = -fuel_radius

  moderator_cell = openmc.Cell(name = 'Moderator')
  moderator_cell.fill = water
  moderator_cell.region = +fuel_radius & -moderator_radius

  root_universe = openmc.Universe(name = 'root universe')
  root_universe.add_cells([fuel_cell, moderator_cell])

  geometry openmc.Geometry(root_universe)

  settings = openmc.Settings()

  settings.batches = 300
  settings.inactive = 20
  settings.particles = 1000

  bounds = [-1., -1., -1., 1., 1., 1.]
  uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=true)
  settings.source = openmc.source.Source(space = uniform_dist)

  settings.output = {'tallies': False}
  
  model = openmc.Model(geometry, materials, settings)

  return model

r_crit, guesses, keffs = openmc.search_for_keff(build_model, bracket=[7., 11.], tol=1e-2,\
print_iterations = True)
print('U235-sphere critical radius: {:1.2f} cm'.format(r_crit))
2

@NotAnotherName The criticality search function relies on root-finding functionality from scipy.optimize. Many of the methods for finding the root of a scalar function require that the first two guesses (based on the bracket you pass) have opposite signs. For a criticality search where you want keff=1, this really means “opposite sides of 1”, so your the low value of your bracketing interval should produce a keff < 1 and the high value of your bracketing interval should produce a keff > 1. Alternatively, you can specify the initial_guess parameter instead of bracket — this will cause it to use the secant root-finding method which doesn’t have any requirement on the sign.

3

Thanks for the tip. The intervall was indeed too small!
The initial guess method worked in both cases though.