diff --git a/docs/source/examples/index.rst b/docs/source/examples/index.rst index c1117b4630..720909284c 100644 --- a/docs/source/examples/index.rst +++ b/docs/source/examples/index.rst @@ -23,4 +23,5 @@ features via the :ref:`pythonapi`. mg-mode-part-iii mdgxs-part-i mdgxs-part-ii + search nuclear-data diff --git a/docs/source/examples/search.ipynb b/docs/source/examples/search.ipynb new file mode 100644 index 0000000000..69359470ce --- /dev/null +++ b/docs/source/examples/search.ipynb @@ -0,0 +1,238 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This Notebook illustrates the usage of the OpenMC Python API's generic eigenvalue search capability. In this Notebook, we will do a critical boron concentration search of a typical PWR pin cell.\n", + "\n", + "To use the search functionality, we must create a function which creates our model according to the input parameter we wish to search for (in this case, the boron concentration). \n", + "\n", + "This notebook will first create that function, and then, run the search." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# Initialize third-party libraries and the OpenMC Python API\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "\n", + "import openmc\n", + "import openmc.model\n", + "\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Create Parametrized Model\n", + "\n", + "To perform the search we will use the `openmc.search_for_keff` function. This function requires a different function be defined which creates an parametrized model to analyze. This model is required to be stored in an `openmc.model.Model` object. The first parameter of this function will be modified during the search process for our critical eigenvalue.\n", + "\n", + "Our model will be a pin-cell from the [Multi-Group Mode Part II](./mg-mode-part-ii.rst) assembly, except this time the entire model building process will be contained within a function, and the Boron concentration will be parametrized." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "# Create the model. `ppm_Boron` will be the parametric variable.\n", + "\n", + "def build_model(ppm_Boron):\n", + " # Create the pin materials\n", + " fuel = openmc.Material(name='1.6% Fuel')\n", + " fuel.set_density('g/cm3', 10.31341)\n", + " fuel.add_element('U', 1., enrichment=1.6)\n", + " fuel.add_element('O', 2.)\n", + "\n", + " zircaloy = openmc.Material(name='Zircaloy')\n", + " zircaloy.set_density('g/cm3', 6.55)\n", + " zircaloy.add_element('Zr', 1.)\n", + "\n", + " water = openmc.Material(name='Borated Water')\n", + " water.set_density('g/cm3', 0.741)\n", + " water.add_element('H', 2.)\n", + " water.add_element('O', 1.)\n", + "\n", + " # Include the amount of boron in the water based on the ppm,\n", + " # neglecting the other constituents of boric acid\n", + " water.add_element('B', ppm_Boron * 1E-6)\n", + " \n", + " # Instantiate a Materials object\n", + " materials = openmc.Materials((fuel, zircaloy, water))\n", + " \n", + " # Create cylinders for the fuel and clad\n", + " fuel_outer_radius = openmc.ZCylinder(R=0.39218)\n", + " clad_outer_radius = openmc.ZCylinder(R=0.45720)\n", + "\n", + " # Create boundary planes to surround the geometry\n", + " min_x = openmc.XPlane(x0=-0.63, boundary_type='reflective')\n", + " max_x = openmc.XPlane(x0=+0.63, boundary_type='reflective')\n", + " min_y = openmc.YPlane(y0=-0.63, boundary_type='reflective')\n", + " max_y = openmc.YPlane(y0=+0.63, boundary_type='reflective')\n", + "\n", + " # Create fuel Cell\n", + " fuel_cell = openmc.Cell(name='1.6% Fuel')\n", + " fuel_cell.fill = fuel\n", + " fuel_cell.region = -fuel_outer_radius\n", + "\n", + " # Create a clad Cell\n", + " clad_cell = openmc.Cell(name='1.6% Clad')\n", + " clad_cell.fill = zircaloy\n", + " clad_cell.region = +fuel_outer_radius & -clad_outer_radius\n", + "\n", + " # Create a moderator Cell\n", + " moderator_cell = openmc.Cell(name='1.6% Moderator')\n", + " moderator_cell.fill = water\n", + " moderator_cell.region = +clad_outer_radius & (+min_x & -max_x & +min_y & -max_y)\n", + "\n", + " # Create root Universe\n", + " root_universe = openmc.Universe(name='root universe', universe_id=0)\n", + " root_universe.add_cells([fuel_cell, clad_cell, moderator_cell])\n", + "\n", + " # Create Geometry and set root universe\n", + " geometry = openmc.Geometry(root_universe)\n", + " \n", + " # Finish with the settings file\n", + " settings = openmc.Settings()\n", + " settings.batches = 300\n", + " settings.inactive = 20\n", + " settings.particles = 1000\n", + " settings.run_mode = 'eigenvalue'\n", + "\n", + " # Create an initial uniform spatial source distribution over fissionable zones\n", + " bounds = [-0.63, -0.63, -10, 0.63, 0.63, 10.]\n", + " uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)\n", + " settings.source = openmc.source.Source(space=uniform_dist)\n", + "\n", + " # We dont need a tallies file so dont waste the disk input/output time\n", + " settings.output = {'tallies': False}\n", + " \n", + " model = openmc.model.Model(geometry, materials, settings)\n", + " \n", + " return model" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Search for the Critical Boron Concentration\n", + "\n", + "To perform the search we imply call the `openmc.search_for_keff` function and pass in the relvant arguments. For our purposes we will be passing in the model building function (`build_model` defined above), a bracketed range for the expected critical Boron concentration (1,000 to 2,500 ppm), the tolerance, and the method we wish to use. \n", + "\n", + "Instead of the bracketed range we could have used a single initial guess, but have elected not to in this example. Finally, due to the high noise inherent in using as few histories as are used in this example, our tolerance on the final keff value will be rather large (1.e-2) and a bisection method will be used for the search." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Iteration: 1; Guess of 1.00e+03 produced a keff of 1.08721 +/- 0.00158\n", + "Iteration: 2; Guess of 2.50e+03 produced a keff of 0.95263 +/- 0.00147\n", + "Iteration: 3; Guess of 1.75e+03 produced a keff of 1.01466 +/- 0.00163\n", + "Iteration: 4; Guess of 2.12e+03 produced a keff of 0.98475 +/- 0.00167\n", + "Iteration: 5; Guess of 1.94e+03 produced a keff of 0.99954 +/- 0.00154\n", + "Iteration: 6; Guess of 1.84e+03 produced a keff of 1.00428 +/- 0.00162\n", + "Iteration: 7; Guess of 1.89e+03 produced a keff of 1.00633 +/- 0.00166\n", + "Iteration: 8; Guess of 1.91e+03 produced a keff of 1.00388 +/- 0.00166\n", + "Iteration: 9; Guess of 1.93e+03 produced a keff of 0.99813 +/- 0.00142\n", + "Critical Boron Concentration: 1926 ppm\n" + ] + } + ], + "source": [ + "# Perform the search\n", + "crit_ppm, guesses, keffs = openmc.search_for_keff(build_model, bracket=[1000., 2500.],\n", + " tol=1.E-2, bracketed_method='bisect',\n", + " print_iterations=True)\n", + "\n", + "print('Critical Boron Concentration: {:4.0f} ppm'.format(crit_ppm))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Finally, the `openmc.search_for_keff` function also provided us with `List`s of the guesses and corresponding keff values generated during the search process with OpenMC. Let's use that information to make a quick plot of the value of keff versus the boron concentration." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfsAAAEyCAYAAAD9bHmuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+8VXWd7/HXW8A8qXhUyARUrNTCVLTjr24qWVfQmUnS\nysy5/qiJadK6/ZAZGbvZ2JST2M1xdHJwImVqNDP0WmloJtJUmgdREB0UfxQHSFCDIskAP/eP9d24\n2Ox9zuZw1tlnr/N+Ph77wVrf7/rx+e512J+9vuu711JEYGZmZuW1Q7MDMDMzs2I52ZuZmZWck72Z\nmVnJOdmbmZmVnJO9mZlZyTnZm5mZlZyTvQ0Iks6SdFez4+iOpLmS/qrZcZg1QtKdks5pdhw2MDjZ\nW7+R9Kyk9ZLW5V5XA0TEtyPipGbHaD1LX3r+mI7fWknzJB3S7LgqJB0o6buSnk/xLZT0GUlDmh1b\nLX3xJVLSFyR9K18WESdHxA3bF52VhZO99be/iIhdcq8Lmh1QWUga2o+7uyAidgH2BOYC/9GbjfR1\nzJLeCDwALAMOiYjdgPcDHcCufbmv/tLPx9VKysneBgRJ50r6r9z8SZKWpDOzf5V0X/7sR9KHJT0u\n6beS5kjaL1cXkj4m6clUf40yr5G0RtJbc8uOTL0Nr5O0u6QfSFqd1vuBpDF14t3iTErS2LTfoWl+\nN0nfkLRS0nJJ/1jrzFLSqLT/PXJlh6ez0mENtvV8SU8CT6Z2fk3SqtxZ7VvTslucQebf8+7W605E\nbARuAsbltvsaSVdKWpFeV0p6TaqbIKlL0t9J+g3wzVT+UUlLJb0o6XZJo3o6nnVC+gfg5xHxmYhY\nmWJcEhEfiog1aXvvkbQ4/S3MlfSW3L6elXRhav9aSd+RtFOu/lRJD0v6naSnJE3q6XhX3mdJV6T4\nn5F0cqr7EnAccLVyPV3VxzWV/bOkZWnf8yUdl8onAX8PnJG28Uj18Za0g6TPSfpVOsazJO2W6ip/\nu+dI+nX627u4p2NvrcXJ3gYcSSOAW4BpZGeOS4C35+onk324nQaMBH4K3Fi1mT8HjgQOAz4ATIyI\nl4HZwJm55T4A3BcRq8j+P3wT2A/YF1gPXN3LZtwAbATeBBwOnARs1VUbESuAXwCn54o/BNwSERsa\nbOtk4GiyhHsScDxwINAOnAG80EC8vVpP0o7AWcD9ueKLgWOA8WTv/1HA53L1rwf2IHufp0g6EbiM\n7FjsDfyK7AtE3lbHs05I7yb726kX74Fk79+nyN7PO4Dvp3ZUfACYBOwPHAqcm9Y9CpgFTCV7j44H\nnk3r9HS8jyb7Ox4BXA58Q5Ii4mKyY3pBjZ6u/HEFeJDsPd0D+E/gu5J2iogfAV8GvpO2cViNpp+b\nXu8E3gDswtZ/2+8ADgLeBXw+/yXISiAi/PKrX15kH4zrgDW510dT3bnAf6Xps4Ff5NYTWbfsX6X5\nO4GP5Op3AF4C9kvzAbwjV38zcFGafjfwdK7uZ8DZdeIdD/w2Nz83F8MXgG/l6sam/Q4F9gJeBtpy\n9WcC99bZz18BP6lq6/Hb0NYTc/UnAk+QJdsdqvazOf4a73nd9WrEOzfFsAb4E7AWeFeu/inglNz8\nRODZND0hrbNTrv4bwOW5+V2ADcDYno5njdg2AJO6if3/ADdXvZ/LgQm5v9G/zNVfDlybpv8N+FqN\nbXZ7vNP7vDRX99rUptfXOi61jmudtvwWOKzW32ONv9d7gI/n6g5K79VQXv3bHZOr/yXwwe39P+/X\nwHn5zN762+SIaM+9rquxzCiyhAdAZJ8+Xbn6/YB/Tt2wa4AXyZLk6Nwyv8lNv0SWQAB+ArRJOjp1\nh48HbgWQ9FpJ/5a6On8HzAPate0Du/YDhgErczH+G/C6OsvfAhybuq6PJ/vg/ek2tDX/Xv2E7Izt\nGuA5STMkDe8p4F6s98mIaAd2IjvrvkXSoaluFNnZecWvUlnF6oj4Y25+i+UjYh1Zr0Ijx7PaC2S9\nA/VU7+sVsvevkX3tQ/ZFplojx3vzNiPipTRZrw0Vy/Izkj6r7HLO2rSP3ch6ChpR65hUvphuFSPd\nv8fWgpzsbSBaCWy+Vp6uz+avnS8D/rrqS0NbRPy8pw2nD/ebyc68PgT8ICJ+n6o/S3bGc3REDCdL\nvJAl12p/IDtDq3h9VXwvAyNy8Q2PiIPrxLQGuIus+/hDwI3pC06jbY2q7V0VEW8DDibrlp/aQMzd\nrVdXRLwSET8FlpJ1XQOsIEuAFfumsprxVi8vaWeyyzfLe9p/DT9my0si1ar3JbIk3si+lgFvrFPe\n8PGuod6jRzeXp+vzf0f2N7J7+qK1llf/Nnt6fGmtY7IReK7BGK3FOdnbQPRD4BBJk5UNeDufLRPT\ntcA0SQfD5sFR79+G7f8n2TXps9J0xa5k1+nXKBswd0k323gYOF7Svmmg07RKRWQDw+4CvippeBoc\n9UZJJ/QQ09lkiSof0za1VdKRqddiGFly/yOwKRfzaakH403ARxpcr1uSjiW7rrw4Fd0IfE7Z4McR\nwOeBb9VbP7X3PEnjlQ3k+zLwQEQ828j+q1wCvF3SdEmvT/G9SdK3JLWTfdH7M0nvSm39LFmi7vGL\nItnlhvPSujtIGi3pzb083nnPkV1H786uZMl5NTBU0ueBfM/Lc8BYSfU+028EPi1pf0m78Oo1/o0N\nxmgtzsne+tv3teXv7G+tXiAinif7udTlZN2y44BOsg9lIuJW4CvATam7/VHg5EYDiIgHyBLaKLJr\n4hVXAm3A82QDzn7UzTbuBr4DLATmAz+oWuRsYEfgMbJrq7fQfffy7cABwHMR8UhuP9va1uHAdWmf\nvyJ7/65IdV8ju17+HNmAsm83uF4tldHj68h+dve5iKi8l/9IdrwWAouAh1JZTRFxD9m19O+R9eq8\nEfhgN/uuKyKeAo4luw69WNLatN1O4PcRsQT4S+BfyI7zX5D9HPRPDWz7l8B5ZO/jWuA+Xj1b3tbj\nnffPwPuUjdS/qs4yc8j+Vp8gOz5/ZMtu/u+mf1+Q9FCN9WeSHad5wDNp/U80GJ+VgF7tLTQbmNLZ\nShdwVkTc2+x4zMxajc/sbUCSNFFSe+rW/Xuya5P397CamZnV4GRvA9WxZCOfK12tkyNifXNDMjNr\nTe7GNzMzKzmf2ZuZmZWck72ZmVnJleZpSiNGjIixY8c2OwwzM7N+M3/+/OcjYmRPy5Um2Y8dO5bO\nzs5mh2FmZtZvJP2q56XcjW9mZlZ6TvZmZmYl52RvZmZWck72ZmZmJedkb2ZmVnKFJXtJMyWtkvRo\nnfo3S/qFpJclXVhVN0nSEklLJV1UVIxmZmaDQZFn9tcDk7qpfxH4JFWP0ZQ0BLiG7DGe44AzJY0r\nKEYzM7PSKyzZR8Q8soRer35VRDwIbKiqOgpYGhFPp2dM3wScWlScZmZmZTcQr9mPBpbl5rtSmZmZ\nmfXCQEz2qlFW89F8kqZI6pTUuXr16oLDMjMza00DMdl3Afvk5scAK2otGBEzIqIjIjpGjuzx1sBm\nZmaD0kBM9g8CB0jaX9KOwAeB25sck5mZWcsq7EE4km4EJgAjJHUBlwDDACLiWkmvBzqB4cArkj4F\njIuI30m6AJgDDAFmRsTiouI0MzMru8KSfUSc2UP9b8i66GvV3QHcUURcZmZmg81A7MY3MzOzPuRk\nb2ZmVnJO9mZmZiXnZG9mZlZyTvZmZmYl52RvZmZWck72ZmZmJedkb2ZmVnJO9mZmZiVX2B30Wtlt\nC5Yzfc4SVqxZz6j2NqZOPIjJh/spu2Zm1pqc7KvctmA502YvYv2GTQAsX7OeabMXATjhm5lZS3I3\nfpXpc5ZsTvQV6zdsYvqcJU2KyMzMbPs42VdZsWb9NpWbmZkNdE72VUa1t21TuZmZ2UDnZF9l6sSD\naBs2ZIuytmFDmDrxoCZFZGZmtn08QK9KZRCeR+ObmVlZONnXMPnw0U7uZmZWGu7GNzMzKzknezMz\ns5IrLNlLmilplaRH69RL0lWSlkpaKOmIXN3lkhZLejwto6LiNDMzK7siz+yvByZ1U38ycEB6TQG+\nDiDp7cD/AA4F3gocCZxQYJxmZmalVliyj4h5wIvdLHIqMCsy9wPtkvYGAtgJ2BF4DTAMeK6oOM3M\nzMqumdfsRwPLcvNdwOiI+AVwL7AyveZExONNiM/MzKwUmpnsa12HD0lvAt4CjCH7QnCipONrbkCa\nIqlTUufq1asLDNXMzKx1NTPZdwH75ObHACuA9wL3R8S6iFgH3AkcU2sDETEjIjoiomPkyJGFB2xm\nZtaKmpnsbwfOTqPyjwHWRsRK4NfACZKGShpGNjjP3fhmZma9VNgd9CTdCEwARkjqAi4hG2xHRFwL\n3AGcAiwFXgLOS6veApwILCIbrPejiPh+UXGamZmVXWHJPiLO7KE+gPNrlG8C/rqouMzMzAYb30HP\nzMys5JzszczMSs7J3szMrOSc7M3MzErOyd7MzKzknOzNzMxKzsnezMys5JzszczMSs7J3szMrOSc\n7M3MzErOyd7MzKzknOzNzMxKzsnezMys5JzszczMSs7J3szMrOSc7M3MzErOyd7MzKzknOzNzMxK\nrrBkL2mmpFWSHq1TL0lXSVoqaaGkI3J1+0q6S9Ljkh6TNLaoOM3MzMquyDP764FJ3dSfDByQXlOA\nr+fqZgHTI+ItwFHAqoJiNDMzK72hRW04Iub1cEZ+KjArIgK4X1K7pL2B3YGhEXF32s66omI0MzMb\nDJp5zX40sCw335XKDgTWSJotaYGk6ZKGNCVCMzOzEmhmsleNsiDrbTgOuBA4EngDcG7NDUhTJHVK\n6ly9enVRcZqZmbW0Zib7LmCf3PwYYEUqXxART0fERuA24Iga6xMRMyKiIyI6Ro4cWXjAZmZmraiZ\nyf524Ow0Kv8YYG1ErAQeBHaXVMneJwKPNStIMzOzVlfYAD1JNwITgBGSuoBLgGEAEXEtcAdwCrAU\neAk4L9VtknQhcI8kAfOB64qK08zMrOyKHI1/Zg/1AZxfp+5u4NAi4jIzMxtsfAc9MzOzknOyNzMz\nKzknezMzs5JzsjczMys5J3szM7OSc7I3MzMrOSd7MzOzknOyNzMzKzknezMzs5JzsjczMys5J3sz\nM7OSc7I3MzMrOSd7MzOzknOyNzMzKzknezMzs5JzsjczMys5J3szM7OSc7I3MzMrOSd7MzOzkiss\n2UuaKWmVpEfr1EvSVZKWSloo6Yiq+uGSlku6uqgYzczMBoMiz+yvByZ1U38ycEB6TQG+XlX/ReC+\nQiIzMzMbRApL9hExD3ixm0VOBWZF5n6gXdLeAJLeBuwF3FVUfGZmZoNFM6/ZjwaW5ea7gNGSdgC+\nCkxtSlRmZmYl08xkrxplAXwcuCMiltWo33ID0hRJnZI6V69e3ecBmpmZlcHQJu67C9gnNz8GWAEc\nCxwn6ePALsCOktZFxEXVG4iIGcAMgI6Ojig+ZDMzs9bTzGR/O3CBpJuAo4G1EbESOKuygKRzgY5a\nid7MzMwa01Cyl7QX8GVgVEScLGkccGxEfKObdW4EJgAjJHUBlwDDACLiWuAO4BRgKfAScN52tMPM\nzMzqUETPvd+S7gS+CVwcEYdJGgosiIhDig6wUR0dHdHZ2dnsMMzMzPqNpPkR0dHTco0O0BsRETcD\nrwBExEZg03bEZ2ZmZv2k0WT/B0l7ko2WR9IxwNrCojIzM7M+0+gAvc+QDah7o6SfASOB9xUWlZmZ\nmfWZhpJ9RDwk6QTgILLfxy+JiA2FRmZmZmZ9otHR+GdXFR0hiYiYVUBMZmZm1oca7cY/Mje9E/Au\n4CHAyd7MzGyAa7Qb/xP5eUm7Af9RSERmZmbWp3p7b/yXyB5Na2ZmZgNco9fsv0/62R3ZF4RxwM1F\nBWVmZmZ9p9Fr9lfkpjcCv4qIrgLiMTMzsz7W6DX7+4oOxMzMzIrRbbKX9Hte7b7fogqIiBheSFRm\nZmbWZ7pN9hGxa38FYmZmZsXYpufZS3od2e/sAYiIX/d5RGZmZtanGvrpnaT3SHoSeAa4D3gWuLPA\nuMzMzKyPNPo7+y8CxwBPRMT+ZHfQ+1lhUZmZmVmfaTTZb4iIF4AdJO0QEfcC4wuMy8zMzPpIo9fs\n10jaBZgHfFvSKrLf25uZmdkA1+iZ/alkt8j9NPAj4CngL4oKyszMzPpOo8l+CjAqIjZGxA0RcVXq\n1q9L0kxJqyQ9Wqdekq6StFTSQklHpPLxkn4haXEqP2PbmmRmZmZ5jSb74cAcST+VdL6kvRpY53pg\nUjf1J5M9TOcAsi8TX0/lLwFnR8TBaf0rJbU3GKeZmZlVaSjZR8Q/pOR7PjAKuE/Sj3tYZx7wYjeL\nnArMisz9QLukvSPiiYh4Mm1jBbAKGNlInGZmZra1bX3E7SrgN8ALwOu2c9+jgWW5+a5Utpmko4Ad\nycYImJmZWS80elOdv5E0F7gHGAF8NCIO3c59q0bZ5vvwS9ob+A/gvIh4pU5cUyR1SupcvXr1doZj\nZmZWTo3+9G4/4FMR8XAf7rsL2Cc3PwZYASBpOPBD4HOpi7+miJgBzADo6Oio9cAeMzOzQa/RR9xe\nJGmIpFH5dbbz3vi3AxdIugk4GlgbESsl7QjcSnY9/7vbsX0zMzOjwWQv6QLgC8BzQKVLPYC6XfmS\nbgQmACMkdQGXAMMAIuJa4A7gFGAp2Qj889KqHwCOB/aUdG4qO7ePexXMzMwGDUX03PstaSlwdE+/\nrW+mjo6O6OzsbHYYZmZm/UbS/Ijo6Gm5Rq/ZLwPWbl9IZlY2ty1YzvQ5S1ixZj2j2tuYOvEgJh8+\nuucVzaxfNZrsnwbmSvoh8HKlMCL+byFRmdmAd9uC5UybvYj1GzYBsHzNeqbNXgTghG82wDT6O/tf\nA3eT/eZ919zLzAap6XOWbE70Fes3bGL6nCVNisjM6ml0NP4/AEjaOSL+UGxIZtYKVqxZv03lZtY8\njd5U51hJjwGPp/nDJP1roZGZ2YA2qr1tm8rNrHka7ca/EphIdptcIuIRsp/HmdkgNXXiQbQNG7JF\nWduwIUydeFCTIjKzehodoEdELJO2uMPtpnrLmln5VQbheTS+2cDX8E/vJL0diHSHu0+SuvTNbPCa\nfPhoJ3ezFtBoN/7HyB5vO5rsnvbj07yZmZkNcI2Oxn8eOKvgWMzMzKwAjd4b/6oaxWuBzoj4f30b\nkpmZmfWlRq/Z7wS8Gag8he50YDHwEUnvjIhPFRGcmVlv+Da+ZltqNNm/CTgxIjYCSPo6cBfwP4FF\nBcVmZrbNfBtfs601OkBvNLBzbn5nYFREbCJ3r3wzs2bzbXzNttbomf3lwMOS5gIiu6HOlyXtDPy4\noNjMzLa5S9638TXbWqOj8b8h6Q7gKLJk//cRsSJVTy0qODMb3HrTJT+qvY3lNRL7qPY2X8u3Qavb\nbnxJb07/HgHsTfZc+18Dr09lZmaF6U2XfL3b+L7zzSOZNnsRy9esJ3j1i8NtC5YXEbrZgNLTmf1n\ngY8CX61RF8CJfR6RmVnSmy75erfx7e6Lg8/urey6TfYR8dH07zv7Jxwzs1d11yXfnVq38f30dx6u\nuayv5dtg0FM3/t/mpt9fVfflHtadKWmVpEfr1EvSVZKWSlqYvywg6RxJT6bXOY01xczKpi+frOdH\n8tpg1tNP7z6Ym55WVTeph3Wv72GZk4ED0msK8HUASXsAlwBHkw0IvETS7j3sy8xKaPLho7nstEMY\n3d6GgNHtbVx22iG96nb3I3ltMOvpmr3qTNea30JEzJM0tptFTgVmRUQA90tql7Q3MAG4OyJeBJB0\nN9mXhht7iNXMSqivnqzX3SN5PUrfyq6nZB91pmvNb6vRZKP7K7pSWb1yM7PtUuuLg++4Z4NBT8n+\nMEm/IzuLb0vTpPmdtnPftXoGopvyrTcgTSG7BMC+++67neGY2WBSOZuvNQDQo/StbLq9Zh8RQyJi\neETsGhFD03Rlfth27rsL2Cc3PwZY0U15rfhmRERHRHSMHDlyO8Mxs8GicjZfK9FXeJS+lUmj98Yv\nwu3A2WlU/jHA2ohYCcwBTpK0exqYd1IqMzPrE7V+c1/No/StTBq9N/42k3Qj2WC7EZK6yEbYDwOI\niGuBO4BTgKXAS8B5qe5FSV8EHkyburQyWM/MrC/0dNbuUfpWNoUl+4g4s4f6AM6vUzcTmFlEXGZm\n9W7WA9nP+zwa38qmmd34ZmZNUe8391eeMZ6fXXSiE72VTmFn9mZmA1V3v7k3KyMnezMblPrqZj1m\nrcDd+GZmZiXnZG9mZlZyTvZmZmYl52RvZmZWck72ZmZmJedkb2ZmVnJO9mZmZiXnZG9mZlZyTvZm\nZmYl52RvZmZWck72ZmZmJed745uZlchtC5b7AT+2FSd7M7OSuG3BcqbNXsT6DZsAWL5mPdNmLwJw\nwh/k3I1vZlYS0+cs2ZzoK9Zv2MT0OUuaFJENFE72ZmYlsWLN+m0qt8HDyd7MrCRGtbdtU7kNHoUm\ne0mTJC2RtFTSRTXq95N0j6SFkuZKGpOru1zSYkmPS7pKkoqM1cys1U2deBBtw4ZsUdY2bAhTJx7U\npIhsoCgs2UsaAlwDnAyMA86UNK5qsSuAWRFxKHApcFla9+3A/wAOBd4KHAmcUFSsZmZlMPnw0Vx2\n2iGMbm9DwOj2Ni477RAPzrNCR+MfBSyNiKcBJN0EnAo8lltmHPDpNH0vcFuaDmAnYEdAwDDguQJj\nNTMrhcmHj3Zyt60U2Y0/GliWm+9KZXmPAKen6fcCu0raMyJ+QZb8V6bXnIh4vMBYzczMSqvIZF/r\nGntUzV8InCBpAVk3/XJgo6Q3AW8BxpB9QThR0vFb7UCaIqlTUufq1av7NnozM7OSKDLZdwH75ObH\nACvyC0TEiog4LSIOBy5OZWvJzvLvj4h1EbEOuBM4pnoHETEjIjoiomPkyJFFtcPMzKylFZnsHwQO\nkLS/pB2BDwK35xeQNEJSJYZpwMw0/WuyM/6hkoaRnfW7G9/MzKwXCkv2EbERuACYQ5aob46IxZIu\nlfSetNgEYImkJ4C9gC+l8luAp4BFZNf1H4mI7xcVq5mZWZkpovoyemvq6OiIzs7OZodhZmbWbyTN\nj4iOnpbzHfTMzMxKzsnezMys5JzszczMSs7J3szMrOSc7M3MzErOyd7MzKzknOzNzMxKzsnezMys\n5JzszczMSs7J3szMrOSc7M3MzErOyd7MzKzknOzNzMxKzsnezMys5JzszczMSs7J3szMrOSc7M3M\nzErOyd7MzKzknOzNzMxKrtBkL2mSpCWSlkq6qEb9fpLukbRQ0lxJY3J1+0q6S9Ljkh6TNLbIWM3M\nzMqqsGQvaQhwDXAyMA44U9K4qsWuAGZFxKHApcBlubpZwPSIeAtwFLCqqFjNzMzKrMgz+6OApRHx\ndET8CbgJOLVqmXHAPWn63kp9+lIwNCLuBoiIdRHxUoGxmpmZlVaRyX40sCw335XK8h4BTk/T7wV2\nlbQncCCwRtJsSQskTU89BWZmZraNikz2qlEWVfMXAidIWgCcACwHNgJDgeNS/ZHAG4Bzt9qBNEVS\np6TO1atX92HoZmZm5VFksu8C9snNjwFW5BeIiBURcVpEHA5cnMrWpnUXpEsAG4HbgCOqdxARMyKi\nIyI6Ro4cWVQ7zMzMWlqRyf5B4ABJ+0vaEfggcHt+AUkjJFVimAbMzK27u6RKBj8ReKzAWM3MzEqr\nsGSfzsgvAOYAjwM3R8RiSZdKek9abAKwRNITwF7Al9K6m8i68O+RtIjsksB1RcVqZmZWZoqovoze\nmjo6OqKzs7PZYZiZmfUbSfMjoqOn5XwHPTMzs5JzsjczMys5J3szM7OSG9rsAMzMzMrutgXLmT5n\nCSvWrGdUextTJx7E5MOr7zNXHCd7MzOzAt22YDnTZi9i/YZNACxfs55psxcB9FvCdze+mZlZgabP\nWbI50Ves37CJ6XOW9FsMTvZmZmYFWrFm/TaVF8HJ3szMrECj2tu2qbwITvZmZmYFmjrxINqGbfng\n1rZhQ5g68aB+i8ED9MzMzApUGYTn0fhmZmYlNvnw0f2a3Ku5G9/MzKzknOzNzMxKzsnezMys5Jzs\nzczMSs7J3szMrOSc7M3MzErOyd7MzKzknOzNzMxKrtBkL2mSpCWSlkq6qEb9fpLukbRQ0lxJY6rq\nh0taLunqIuM0MzMrs8KSvaQhwDXAycA44ExJ46oWuwKYFRGHApcCl1XVfxG4r6gYzczMBoMiz+yP\nApZGxNMR8SfgJuDUqmXGAfek6Xvz9ZLeBuwF3FVgjGZmZqVXZLIfDSzLzXelsrxHgNPT9HuBXSXt\nKWkH4KvA1ALjMzMzGxSKTPaqURZV8xcCJ0haAJwALAc2Ah8H7oiIZXRD0hRJnZI6V69e3Rcxm5mZ\nlU6RT73rAvbJzY8BVuQXiIgVwGkAknYBTo+ItZKOBY6T9HFgF2BHSesi4qKq9WcAMwA6Ojqqv0iY\nmZkZxSb7B4EDJO1Pdsb+QeBD+QUkjQBejIhXgGnATICIOCu3zLlAR3WiNzMzs8YU1o0fERuBC4A5\nwOPAzRGxWNKlkt6TFpsALJH0BNlgvC8VFY+ZmdlgpYhy9H53dHREZ2dns8MwMzPrN5LmR0RHT8v5\nDnpmZmYl52RvZmZWck72ZmZmJedkb2ZmVnJO9mZmZiXnZG9mZlZypfnpnaTVwK/6eLMjgOf7eJsD\ngdvVWtyu1uJ2tZ5Wbtt+ETGyp4VKk+yLIKmzkd8vthq3q7W4Xa3F7Wo9ZW5bhbvxzczMSs7J3szM\nrOSc7Ls3o9kBFMTtai1uV2txu1pPmdsG+Jq9mZlZ6fnM3szMrOQGVbKXNFPSKkmP5sr2kHS3pCfT\nv7unckm6StJSSQslHZFb55y0/JOSzmlGW/LqtGu6pP9Osd8qqT1XNy21a4mkibnySalsqaSL+rsd\ntdRqW67uQkkhaUSab+ljlso/kY7BYkmX58pb4pjV+VscL+l+SQ9L6pR0VCpvieMlaR9J90p6PB2X\n/53Ky/DZUa9tLf35Ua9dufqW/ezotYgYNC/geOAI4NFc2eXARWn6IuArafoU4E5AwDHAA6l8D+Dp\n9O/uaXqpn5qZAAAItElEQVT3Adiuk4ChaforuXaNAx4BXgPsDzwFDEmvp4A3ADumZcYNxGOWyvcB\n5pDdW2FESY7ZO4EfA69J869rtWNWp113ASfnjtHcVjpewN7AEWl6V+CJdEzK8NlRr20t/flRr11p\nvqU/O3r7GlRn9hExD3ixqvhU4IY0fQMwOVc+KzL3A+2S9gYmAndHxIsR8VvgbmBS8dHXV6tdEXFX\nRGxMs/cDY9L0qcBNEfFyRDwDLAWOSq+lEfF0RPwJuCkt21R1jhnA14C/BfKDTlr6mAF/A/xTRLyc\nllmVylvmmNVpVwDD0/RuwIo03RLHKyJWRsRDafr3wOPAaMrx2VGzba3++dHNMYMW/+zorUGV7OvY\nKyJWQvYHArwulY8GluWW60pl9coHsg+TfWuFErRL0nuA5RHxSFVVq7ftQOA4SQ9Iuk/Skam81dv1\nKWC6pGXAFcC0VN5y7ZI0FjgceICSfXZUtS2vpT8/8u0q8WdHj4Y2O4ABTDXKopvyAUnSxcBG4NuV\nohqLBbW/+A24dkl6LXAxWTfjVtU1ylrpmA0l6yo8BjgSuFnSG2jxY0bWY/HpiPiepA8A3wDeTYsd\nL0m7AN8DPhURv5NqhZktWqNswLYLtm5brrylPz/y7SJrR1k/O3rkM3t4LnXXkP6tdJ12kV3bqRhD\n1v1Yr3zASYNJ/hw4K9IFKFq/XW8ku1b4iKRnyeJ8SNLraf22dQGzU1fiL4FXyO7Z3ertOgeYnaa/\nS9blCy3ULknDyJLGtyOi0pZSfHbUaVvLf37UaFeZPzt61uxBA/39Asay5eCh6Ww5yObyNP1nbDlg\n45fx6oCNZ8jOwHZP03sMwHZNAh4DRlYtdzBbDrB5mmxwzdA0vT+vDrA5uNntqtW2qrpneXWQTasf\ns48Bl6bpA8m6D9Vqx6xGux4HJqTpdwHzW+l4pfhmAVdWlbf8Z0c3bWvpz4967apapmU/O3r1njQ7\ngH7+A7gRWAlsIPvG9hFgT+Ae4Mn07x65P5ZryEaYLgI6ctv5MNnAlKXAeQO0XUvJksXD6XVtbvmL\nU7uWkEZJp/JTyEatPgVc3Ox21WtbVX3+P2yrH7MdgW8BjwIPASe22jGr0653APNTAngAeFsrHa8U\nfwALc/+fTinJZ0e9trX050e9dlUt05KfHb19+Q56ZmZmJedr9mZmZiXnZG9mZlZyTvZmZmYl52Rv\nZmZWck72ZmZmJedkb9ZLkjalJ7k9IukhSW9vQgxnS3o0PdnrMUkX9ncMVfGMl3RKL9YbK+lDufkO\nSVf1UUyV4zSqL7bXzX6+LelFSe8rcj9mveFkb9Z76yNifEQcRna/98saXVHSkO3duaSTyW4DelJE\nHEz2tLm127vd7TSe7PfWW5HU3e25xwKbk31EdEbEJ/sopspxKvTOZxFxFnB7kfsw6y0ne7O+MRz4\nLWx+Nvb0dMa9SNIZqXxCesb2f5LduANJn0nLPSrpU6lsbHoO93XpjP0uSW019jkNuLCSxCLijxFx\nXdpG5RnyleeRV561PlfSVyT9UtITko5L5UMkXZHiXSjpE6n8bemhPPMlzcndHnar7UjaEbgUOCOd\nSZ8h6QuSZki6C5iV2vbT1BOS7w35J7KHAD0s6dPpvfpB2tcekm5Lcd0v6dBU/gVJM1MsT0tq6MuB\npHWSvpr2f4+kkbk2XSnp5+l4HJXbzw3pODwr6TRJl6f36kfptqxmA1uz7+rjl1+t+gI2kd2Z67/J\nzqgrd4Y7nexRmEOAvYBfkz1fewLwB2D/tNzbyJL+zsAuwGKyp3ONJXtox/i03M3AX9bY/4vAbnVi\nWwickKYvJd02FJgLfDVNnwL8OE3/Ddl9xCvPMN8DGAb8nHTLVOAMYGYP2zkXuDoXxxfI7p7XluZf\nC+yUpg8AOtP0BOAHufU2zwP/AlySpk8EHs5t++dkt24dAbwADKvxXqyrmg+y+70DfL4Sb2rTdWn6\neNItf9N+/iu9H4cBL5HuHAfcCkzObft64H3N/tv0y6/ql596Z9Z76yNiPICkY8nOXN9KdqvOGyNi\nE9nDUu4je4rd78juuf1MWv8dwK0R8Ye0jdnAcWRdwc9ExMNpuflkXwAaImk3oD0i7ktFN5A9gKai\n8rCT/HbfTXZL1I0AEfFiastbgbuVPeFtCNmtcLvbTi23R8T6ND0MuFrSeLIvSwc20KR3kH2BIiJ+\nImnP1EaAH0bEy8DLklaRfbnq6mF7rwDfSdPfyrUDstv9EhHzJA2X1J7K74yIDZIWkb0PP0rli9iG\nY2PWLE72Zn0gIn4haQQwktqPxaz4Q266u+Vezk1vAmp14y8m6x34SaNxVm17E69+BoitH90pYHFE\nHLsN26kl3+ZPA8+RnSHvAPyxgXi7e8xo9fvUm8+0qDO91X4i4hVJGyKiUv5KL/dp1q98zd6sD0h6\nM9kZ3wvAPLLr1kPS9eDjgV/WWG0eMFnSayXtDLwX+Ok27PYy4HJlj+hE0mskfTIi1gK/rVyPB/4X\ncF+9jSR3AR+rDKKTtAfZg05Gpl4LJA2TdHAP2/k9sGs39bsBKyPilRRXZaBid+vNA85KMUwAno/c\nM9d7YQegMmL+Q2Rd9BWV8RXvANam99Ks5fkbqVnvtUmqdLULOCciNkm6FTiW7ClvAfxtRPwmfSHY\nLCIeknQ9r34R+PeIWCBpbCM7j4g7JO0F/FhZP3sAM1P1OcC1kl5L9ujR83rY3L+TdakvlLSB7Nr1\n1cp+RnZV6jYfClxJ1qNQz73ARel9qfXrhH8Fvifp/WnZyln/QmCjpEfIrnsvyK3zBeCbkhaSXS8/\np4e29OQPwMGS5pONtTgjV/dbST8nG3D54e3cj9mA4afemVmpSVoXEbvUm8+VzyX7dUPnduzrerKB\nhbf0dhtmRXA3vpmV3e/UTzfVAU6gsXEIZv3KZ/ZmZmYl5zN7MzOzknOyNzMzKzknezMzs5Jzsjcz\nMys5J3szM7OSc7I3MzMruf8PsJgBT1rI6FgAAAAASUVORK5CYII=\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plt.figure(figsize=(8, 4.5))\n", + "plt.title('Eigenvalue versus Boron Concentration')\n", + "# Create a scatter plot using the mean value of keff\n", + "plt.scatter(guesses, [keffs[i][0] for i in range(len(keffs))])\n", + "plt.xlabel('Boron Concentration [ppm]')\n", + "plt.ylabel('Eigenvalue')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": true + }, + "source": [ + "We see a nearly linear reactivity coefficient for the boron concentration, exactly as one would expect for a pure 1/v absorber at small concentrations." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.6.0" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/docs/source/examples/search.rst b/docs/source/examples/search.rst new file mode 100644 index 0000000000..9bb75b5830 --- /dev/null +++ b/docs/source/examples/search.rst @@ -0,0 +1,13 @@ +.. _notebook_search: + +================== +Criticality Search +================== + +.. only:: html + + .. notebook:: search.ipynb + +.. only:: latex + + IPython notebooks must be viewed in the online HTML documentation. diff --git a/docs/source/pythonapi/index.rst b/docs/source/pythonapi/index.rst index c23bf19211..161b8564ab 100644 --- a/docs/source/pythonapi/index.rst +++ b/docs/source/pythonapi/index.rst @@ -173,6 +173,7 @@ Running OpenMC openmc.calculate_volumes openmc.plot_geometry openmc.plot_inline + openmc.search_for_keff Post-processing --------------- @@ -337,6 +338,19 @@ Functions openmc.model.create_triso_lattice openmc.model.pack_trisos +Model Container +--------------- + +Classes ++++++++ + +.. autosummary:: + :toctree: generated + :nosignatures: + :template: myclass.rst + + openmc.model.Model + -------------------------------------------- :mod:`openmc.data` -- Nuclear Data Interface -------------------------------------------- diff --git a/openmc/__init__.py b/openmc/__init__.py index 53fbec3f22..d8e7e24822 100644 --- a/openmc/__init__.py +++ b/openmc/__init__.py @@ -26,3 +26,4 @@ from openmc.summary import * from openmc.particle_restart import * from openmc.mixin import * from openmc.plotter import * +from openmc.search import * \ No newline at end of file diff --git a/openmc/model/__init__.py b/openmc/model/__init__.py index ffb1f42820..557effcefa 100644 --- a/openmc/model/__init__.py +++ b/openmc/model/__init__.py @@ -1 +1,2 @@ from .triso import * +from .model import * diff --git a/openmc/model/model.py b/openmc/model/model.py new file mode 100644 index 0000000000..7716cae97b --- /dev/null +++ b/openmc/model/model.py @@ -0,0 +1,164 @@ +import openmc +from openmc.checkvalue import check_type + + +class Model(object): + """OpenMC model container for the openmc.Geometry, openmc.Materials, + openmc.Settings, openmc.Tallies, openmc.CMFD objects, and openmc.Plot + objects + + Parameters + ---------- + geometry : openmc.Geometry + Geometry information + materials : openmc.Materials + Materials information + settings : openmc.Settings + Settings information + tallies : openmc.Tallies + Tallies information, optional + cmfd : openmc.CMFD + CMFD information, optional + plots : openmc.Plots + Plot information, optional + + Attributes + ---------- + geometry : openmc.Geometry + Geometry information + materials : openmc.Materials + Materials information + settings : openmc.Settings + Settings information + tallies : openmc.Tallies + Tallies information + cmfd : openmc.CMFD + CMFD information + plots : openmc.Plots + Plot information + + """ + + def __init__(self, geometry, materials, settings, tallies=None, cmfd=None, + plots=None): + self.geometry = geometry + self.materials = materials + self.settings = settings + if tallies: + self.tallies = tallies + else: + self._tallies = openmc.Tallies() + if cmfd: + self.cmfd = cmfd + else: + self._cmfd = None + if plots: + self.plots = plots + else: + self.plots = openmc.Plots() + + self.sp = None + + @property + def geometry(self): + return self._geometry + + @property + def materials(self): + return self._materials + + @property + def settings(self): + return self._settings + + @property + def tallies(self): + return self._tallies + + @property + def cmfd(self): + return self._cmfd + + @property + def plots(self): + return self._plots + + @geometry.setter + def geometry(self, geometry): + check_type('geometry', geometry, openmc.Geometry) + self._geometry = geometry + + @materials.setter + def materials(self, materials): + check_type('materials', materials, openmc.Materials) + self._materials = materials + + @settings.setter + def settings(self, settings): + check_type('settings', settings, openmc.Settings) + self._settings = settings + + @tallies.setter + def tallies(self, tallies): + check_type('tallies', tallies, openmc.Tallies) + self._tallies = tallies + + @cmfd.setter + def cmfd(self, cmfd): + check_type('cmfd', cmfd, openmc.CMFD) + self._cmfd = cmfd + + @plots.setter + def plots(self, plots): + check_type('plots', plots, openmc.Plots) + self._plots = plots + + def export_to_xml(self): + """Export model settings to XML files. + """ + + self.geometry.export_to_xml() + self.materials.export_to_xml() + self.settings.export_to_xml() + self.tallies.export_to_xml() + if self.cmfd is not None: + self.cmfd.export_to_xml() + self.plots.export_to_xml() + + def run(self, **kwargs): + """Creates the XML files, runs OpenMC, and loads the statepoint. + + Parameters + ---------- + **kwargs + All keyword arguments are passed to openmc.run + + Returns + ------- + 2-tuple of float + k_combined from the statepoint + + """ + + self.export_to_xml() + + return_code = openmc.run(**kwargs) + + assert (return_code == 0), "OpenMC did not execute successfully" + + statepoint_batches = self.settings.batches + if self.settings.statepoint is not None: + if 'batches' in self.settings.statepoint: + statepoint_batches = self.settings.statepoint['batches'][-1] + self.sp = \ + openmc.StatePoint('statepoint.{}.h5'.format(statepoint_batches)) + + return self.sp.k_combined + + def close(self): + """Close the statepoint and summary files + """ + + if self.sp is not None: + self.sp._f.close() + self.sp.summary._f.close() diff --git a/openmc/search.py b/openmc/search.py new file mode 100644 index 0000000000..dd2c4345af --- /dev/null +++ b/openmc/search.py @@ -0,0 +1,201 @@ +from collections import Callable +from numbers import Real + +import openmc +import openmc.model +import openmc.checkvalue as cv + + +_SCALAR_BRACKETED_METHODS = ['brentq', 'brenth', 'ridder', 'bisect'] + + +def _search_keff(guess, target, model_builder, model_args, print_iterations, + print_output, guesses, results): + """Function which will actually create our model, run the calculation, and + obtain the result. This function will be passed to the root finding + algorithm + + Parameters + ---------- + guess : Real + Current guess for the parameter to be searched in `model_builder`. + target_keff : Real + Value to search for + model_builder : collections.Callable + Callable function which builds a model according to a passed + parameter. This function must return an openmc.model.Model object. + model_args : dict + Keyword-based arguments to pass to the `model_builder` method. + print_iterations : bool + Whether or not to print the guess and the resultant keff during the + iteration process. + print_output : bool + Whether or not to print the OpenMC output during the iterations. + guesses : Iterable of Real + Running list of guesses thus far, to be updated during the execution of + this function. + results : Iterable of Real + Running list of results thus far, to be updated during the execution of + this function. + + Returns + ------- + float + Value of the model for the current guess compared to the target value. + + """ + + # Build the model + model = model_builder(guess, **model_args) + + # Run the model and obtain keff + keff = model.run(output=print_output) + + # Close the model to ensure HDF5 will allow access during the next + # OpenMC execution + model.close() + + # Record the history + guesses.append(guess) + results.append(keff) + + if print_iterations: + text = 'Iteration: {}; Guess of {:.2e} produced a keff of ' + \ + '{:1.5f} +/- {:1.5f}' + print(text.format(len(guesses), guess, keff[0], keff[1])) + + return (keff[0] - target) + + +def search_for_keff(model_builder, initial_guess=None, target=1.0, + bracket=None, model_args=None, tol=None, + bracketed_method='bisect', print_iterations=False, + print_output=False, **kwargs): + """Function to perform a keff search by modifying a model parametrized by a + single independent variable. + + Parameters + ---------- + model_builder : collections.Callable + Callable function which builds a model according to a passed + parameter. This function must return an openmc.model.Model object. + initial_guess : Real, optional + Initial guess for the parameter to be searched in + `model_builder`. One of `guess` or `bracket` must be provided. + target : Real, optional + keff value to search for, defaults to 1.0. + bracket : None or Iterable of Real, optional + Bracketing interval to search for the solution; if not provided, + a generic non-bracketing method is used. If provided, the brackets + are used. Defaults to no brackets provided. One of `guess` or `bracket` + must be provided. If both are provided, the bracket will be + preferentially used. + model_args : dict, optional + Keyword-based arguments to pass to the `model_builder` method. Defaults + to no arguments. + tol : float + Tolerance to pass to the search method + bracketed_method : {'brentq', 'brenth', 'ridder', 'bisect'}, optional + Solution method to use; only applies if + `bracket` is set, otherwise the Secant method is used. + Defaults to 'bisect'. + print_iterations : bool + Whether or not to print the guess and the result during the iteration + process. Defaults to False. + print_output : bool + Whether or not to print the OpenMC output during the iterations. + Defaults to False. + **kwargs + All remaining keyword arguments are passed to the root-finding + method. + + Returns + ------- + zero_value : float + Estimated value of the variable parameter where keff is the + targeted value + guesses : List of Real + List of guesses attempted by the search + results : List of 2-tuple of Real + List of keffs and uncertainties corresponding to the guess attempted by + the search + + """ + + if initial_guess is not None: + cv.check_type('initial_guess', initial_guess, Real) + if bracket is not None: + cv.check_iterable_type('bracket', bracket, Real) + cv.check_length('bracket', bracket, 2) + cv.check_less_than('bracket values', bracket[0], bracket[1]) + if model_args is None: + model_args = {} + else: + cv.check_type('model_args', model_args, dict) + cv.check_type('target', target, Real) + cv.check_type('tol', tol, Real) + cv.check_value('bracketed_method', bracketed_method, + _SCALAR_BRACKETED_METHODS) + cv.check_type('print_iterations', print_iterations, bool) + cv.check_type('print_output', print_output, bool) + cv.check_type('model_builder', model_builder, Callable) + + # Run the model builder function once to make sure it provides the correct + # output type + if bracket is not None: + model = model_builder(bracket[0], **model_args) + elif initial_guess is not None: + model = model_builder(initial_guess, **model_args) + cv.check_type('model_builder return', model, openmc.model.Model) + + import scipy.optimize as sopt + + # Set the iteration data storage variables + guesses = [] + results = [] + + # Set the searching function (for easy replacement should a later + # generic function be added. + search_function = _search_keff + + if bracket is not None: + # Generate our arguments + args = {'f': search_function, 'a': bracket[0], 'b': bracket[1]} + if tol is not None: + args['rtol'] = tol + + # Set the root finding method + if bracketed_method == 'brentq': + root_finder = sopt.brentq + elif bracketed_method == 'brenth': + root_finder = sopt.brenth + elif bracketed_method == 'ridder': + root_finder = sopt.ridder + elif bracketed_method == 'bisect': + root_finder = sopt.bisect + + elif initial_guess is not None: + + # Generate our arguments + args = {'func': search_function, 'x0': initial_guess} + if tol is not None: + args['tol'] = tol + + # Set the root finding method + root_finder = sopt.newton + + else: + raise ValueError("Either the 'bracket' or 'initial_guess' parameters " + "must be set") + + # Add information to be passed to the searching function + args['args'] = (target, model_builder, model_args, print_iterations, + print_output, guesses, results) + + # Create a new dictionary with the arguments from args and kwargs + args.update(kwargs) + + # Perform the search + zero_value = root_finder(**args) + + return zero_value, guesses, results