373 lines
14 KiB
Plaintext
373 lines
14 KiB
Plaintext
{
|
||
"cells": [
|
||
{
|
||
"attachments": {},
|
||
"cell_type": "markdown",
|
||
"id": "d63623a6",
|
||
"metadata": {},
|
||
"source": [
|
||
"# Source Convergence\n",
|
||
"\n",
|
||
"Monte Carlo \"run strategy\" generally consists of four choices:\n",
|
||
"\n",
|
||
"- Number of batches (`model.settings.batches`)\n",
|
||
"- Number of inactive batches (`model.settings.inactive`)\n",
|
||
"- Number of particles per batch (`model.settings.particles`)\n",
|
||
"- Number of generations per batch (`model.settings.generations_per_batch`), which defaults to 1\n",
|
||
"\n",
|
||
"In addition, for k-eigenvalue calculations you elect the initial source distribution. \n",
|
||
"\n",
|
||
"### What is a generation?\n",
|
||
"\n",
|
||
"During each generation in OpenMC:\n",
|
||
"\n",
|
||
"- `model.settings.particles` sites are sampled from a bank of fission sites\n",
|
||
"- Those particles are transported from birth until \"death\" (leakage or absorption)\n",
|
||
"- Any new neutrons which are produced during their lifetime (e.g., from fission) have their birth position, energy, and angle stored in a new fission bank.\n",
|
||
"\n",
|
||
"Therefore, each time a new generation begins, the sites used to sample the neutron positions are those which were \"banked\" from the previous generation.\n",
|
||
"\n",
|
||
"### What is a batch?\n",
|
||
"\n",
|
||
"A batch is a group of generations as one statistical realization. By default, `model.settings.generations_per_batch` is unity, in which case a batch is synonymous with a generation.\n",
|
||
"\n",
|
||
"### What is an inactive batch?\n",
|
||
"\n",
|
||
"An inactive batch involves the transport of neutrons, but not accumulation of any tallies.\n",
|
||
"\n",
|
||
"### What is an active batch?\n",
|
||
"\n",
|
||
"An active batch involves the transport of neutrons AND accumulation of tallies."
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"id": "09ac1456",
|
||
"metadata": {},
|
||
"source": [
|
||
"## Inactive batches and initial source\n",
|
||
"\n",
|
||
"For the very first generation run in OpenMC (the first generation of the first batch), we don't yet have a bank of sites to sample neutrons from. Therefore, the user has to provide a \"guess\" for the starting source, using `model.settings.source`. By default, this source distribution is taken to be a point source at the origin, isotropic in angle and with the Watt fission spectrum. In all k-eigenvalue cases, you don't know the source distribution, so whatever choice you make is a GUESS. However, we don't want to start accumulating tally statistics until our guessed source distribution undergoes enough generations so that it has statistically converged. Otherwise, even if you run a high number of particles/batch, your answers will be WRONG (but PRECISE!).\n",
|
||
"\n",
|
||
"The factors you choose are:\n",
|
||
"\n",
|
||
"- The starting source distribution. Ideally, you'd like this to be as close as possible to the true source distribution.\n",
|
||
"- The number of inactive batches. The higher this number is, the more likely you are to have evolved from your (WRONG) initial source to the true initial source."
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"id": "37c348fd",
|
||
"metadata": {},
|
||
"source": [
|
||
"To explore these effects, let's choose an initial source distribution which is clearly not the correct, converged source distribution. When choosing a source distribution, we will use the [`openmc.stats` module](https://docs.openmc.org/en/stable/pythonapi/stats.html), which contains probability distributions to sample from which we can use for energy, angle, and space."
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": null,
|
||
"id": "d8cc2d8b-dcc5-4edd-a028-3db368658768",
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"import openmc\n",
|
||
"import matplotlib.pyplot as plt\n",
|
||
"import pandas as pd\n",
|
||
"import numpy as np"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": null,
|
||
"id": "64693c53",
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"model = openmc.examples.pwr_pin_cell()\n",
|
||
"model.geometry.root_universe.plot()"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": null,
|
||
"id": "c9a90a0d",
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"model.geometry.bounding_box"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": null,
|
||
"id": "a73fe3c6",
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"bottom = openmc.ZPlane(z0=-150, boundary_type='vacuum')\n",
|
||
"top = openmc.ZPlane(z0=150, boundary_type='vacuum')\n",
|
||
"box = openmc.model.RectangularPrism(0.63*2, 0.63*2, boundary_type='reflective')\n",
|
||
"outer_cell = openmc.Cell(region=-box & +bottom & -top, fill=model.geometry.root_universe)\n",
|
||
"root_universe = openmc.Universe(cells=[outer_cell])\n",
|
||
"model.geometry.root_universe = root_universe"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"id": "95a8984d",
|
||
"metadata": {},
|
||
"source": [
|
||
"Let's make a source distribution which only exists over the lower half of the pincell -- as if the entire upper half had no fission at all. Clearly this is a very bad guess!"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": null,
|
||
"id": "2d9db573",
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": []
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"id": "8075e4b8",
|
||
"metadata": {},
|
||
"source": [
|
||
"In order to \"see\" the effect of having a wrong initial source distribution, we will need to add a tally to visualize a quantity affected by a wrong source distribution. Let's take a look at the heating distribution, with a mesh filter with 50 layers in the $z$ direction."
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": null,
|
||
"id": "18429b99",
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": []
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"id": "9d1b6e86",
|
||
"metadata": {},
|
||
"source": [
|
||
"Now, we will run this model with different choices for numbers of inactive batches. As long as we run enough inactive batches, we can evolve away from this very wrong source distribution to the true, cosine distribution we expect. We know that the source distribution should be symmetric about the midplane of the fuel, so we can monitor a quantity like the axial offset,\n",
|
||
"\n",
|
||
"$\\text{axial offset}\\equiv \\int_\\text{top half} \\dot{q}\\ dV-\\int_\\text{bottom half} \\dot{q}\\ dV$"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": null,
|
||
"id": "f1e03b6a",
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"model.settings.particles = 5000 \n",
|
||
"\n",
|
||
"fig, ax = plt.subplots()\n",
|
||
"axial_offset = []\n",
|
||
"inactive_batches = [1, 10, 50, 100, 200, 300]\n",
|
||
"\n",
|
||
"for i in inactive_batches:\n",
|
||
" \n",
|
||
" \n",
|
||
"plt.legend(ncol=3)\n",
|
||
"plt.grid()\n",
|
||
"plt.xlabel('Axial Position')\n",
|
||
"plt.ylabel('Kappa-Fission (unnormalized)')"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": null,
|
||
"id": "1d7d5bc0-fd6b-499f-b3b5-149cee89ed8b",
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"data = {'Inactive Batches' : inactive_batches, 'Axial Offset' : axial_offset}\n",
|
||
"df = pd.DataFrame(data=data)\n",
|
||
"df"
|
||
]
|
||
},
|
||
{
|
||
"attachments": {},
|
||
"cell_type": "markdown",
|
||
"id": "04a11a94",
|
||
"metadata": {},
|
||
"source": [
|
||
"We can see that we need more and more inactive batches if we want our fission source distribution to approach the true distribution. Based on monitoring the axial offset, we might conclude that somewhere around 100 inactive batches are required (other variations in axial offset arise from statistical noise).\n",
|
||
"The number of required inactive batches will depend on:\n",
|
||
"\n",
|
||
"- How big the problem is, and how far a neutron typically travels from birth to death. If the problem dimensions are large compared to the mean free path, it will take more generations for information to \"travel\" across the system.\n",
|
||
"- How realistic the initial source distribution is.\n",
|
||
"\n",
|
||
"If you increase the number of particles, the statistical noisiness will decrease, but it won't help us to get closer to the true distribution for the sampled source. As a rule of thumb on selecting the inactive batches:\n",
|
||
"\n",
|
||
"- Small problems (critical assemblies, pin cell): O(10) batches\n",
|
||
"- Medium problems (assemblies, fast reactor): O(10)–O(100) batches\n",
|
||
"- Large problems (PWR core, spent fuel pool): O(100)–O(1000) batches?\n",
|
||
"\n",
|
||
"The dominance ratio, or the ratio of the second-largest eigenvalue to the largest eigenvalue, dictates the convergence rate of power iteration (each generation in a Monte Carlo simulation is a power iteration).\n",
|
||
"\n",
|
||
"If we change our initial guess for the source distribution to something more realistic, we'll be able to use fewer inactive batches before we start collecting tally statistics."
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": null,
|
||
"id": "9aed66fb",
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": []
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": null,
|
||
"id": "873d004c",
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"fig, ax = plt.subplots()\n",
|
||
"axial_offset = []\n",
|
||
"inactive_batches = [1, 10, 50, 100, 200, 300]\n",
|
||
"\n",
|
||
"for i in inactive_batches:\n",
|
||
" model.settings.inactive = i\n",
|
||
" model.settings.batches = model.settings.inactive + 200\n",
|
||
" statepoint = model.run(output=False)\n",
|
||
" \n",
|
||
" with openmc.StatePoint(statepoint) as sp:\n",
|
||
" t = sp.get_tally(id=heat.id).get_values().squeeze()\n",
|
||
"\n",
|
||
" axial_offset.append(sum(t[25:]) - sum(t[:25]))\n",
|
||
" ax.plot(t, label='{} inactive'.format(i))\n",
|
||
" \n",
|
||
"plt.legend(ncol=3)\n",
|
||
"plt.grid()\n",
|
||
"plt.xlabel('Axial Position')\n",
|
||
"plt.ylabel('Kappa-Fission (unnormalized)')"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": null,
|
||
"id": "b63bd071-7719-47da-9cb1-85693182d3da",
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"data = {'Inactive Batches' : inactive_batches, 'Axial Offset' : axial_offset}\n",
|
||
"df = pd.DataFrame(data=data)\n",
|
||
"df"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"id": "b599b775",
|
||
"metadata": {},
|
||
"source": [
|
||
"## Shannon Entropy\n",
|
||
"\n",
|
||
"Of course, we know that using just 1 inactive batch is most likely inadequate - but our inspection of axial offset would fail to illuminate this. Instead, we should use a metric which considers the more detailed spatial variation (not just lumping the pincell into two halves).\n",
|
||
"\n",
|
||
"Shannon entropy is a concept from information theory that characterizes how much \"information\" a bit stream stores. In the context of eigenvalue calculations, it has been shown that when a source distibution is discretized over a mesh, the entropy of the source probability converges as the distribution itself reaches stationarity. Shannon entropy is defined as:\n",
|
||
"\n",
|
||
"$$ H = - \\sum_i p_i \\log_2 p_i $$\n",
|
||
"\n",
|
||
"where $p_i$ is the fraction of source particles in mesh cell $i$. In essence, we superimpose a mesh onto our geometry, and assess stationarity of the fission source by looking for stationarity in the Shannon entropy. Note that you will not want to run with Shannon entropy enabled all the time, because the simulation will be quite a bit slower.\n",
|
||
"\n",
|
||
"We'll need to create a `RegularMesh` object that will be used to count source particles over and calculate Shannon entropy. We can use the same mesh used for talling the fission source earlier (though this is not a requirement)."
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": null,
|
||
"id": "fd761c97",
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": []
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": null,
|
||
"id": "22e61eb2",
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"statepoint = model.run()"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": null,
|
||
"id": "fe53107b",
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"\n",
|
||
"\n",
|
||
"plt.plot(entropy)\n",
|
||
"plt.xlabel('Batch')\n",
|
||
"plt.ylabel('Shannon entropy')"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"id": "9ac5d111",
|
||
"metadata": {},
|
||
"source": [
|
||
"You should choose your number of inactive batches once the Shannon entropy plateaus to a constant value. It's important to use enough particles so that you can see the plateau. In the above, choosing 100 or higher looks to be a good choice, but we can re-run with higher particles to confirm. Based on this second run, 200 looks like a slightly better choice."
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": null,
|
||
"id": "bf37b32e",
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"\n",
|
||
" \n",
|
||
"plt.plot(entropy)\n",
|
||
"plt.xlabel('Batch')\n",
|
||
"plt.ylabel('Shannon entropy')"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"id": "2dd21030",
|
||
"metadata": {},
|
||
"source": [
|
||
"From the above, it looks like we can justify using about 200 inactive batches. It is common to use a moving window average in order to make this determination more quantitatively."
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": null,
|
||
"id": "f4dfea7b-4e86-4b9d-a1ca-e8e67849a833",
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": []
|
||
}
|
||
],
|
||
"metadata": {
|
||
"kernelspec": {
|
||
"display_name": "Python 3 (ipykernel)",
|
||
"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.12.4"
|
||
}
|
||
},
|
||
"nbformat": 4,
|
||
"nbformat_minor": 5
|
||
}
|