diff --git a/CMakeLists.txt b/CMakeLists.txt index b761929523..c003f39ac1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -165,7 +165,6 @@ endif() add_subdirectory(vendor/xtl) set(xtl_DIR ${CMAKE_CURRENT_BINARY_DIR}/vendor/xtl) add_subdirectory(vendor/xtensor) -target_link_libraries(xtensor INTERFACE xtl) #=============================================================================== # GSL header-only library diff --git a/docs/source/usersguide/materials.rst b/docs/source/usersguide/materials.rst index 8472768489..85584a6c46 100644 --- a/docs/source/usersguide/materials.rst +++ b/docs/source/usersguide/materials.rst @@ -40,6 +40,12 @@ of an element, you specify the element itself. For example, mat.add_element('C', 1.0) +This method can also accept case-insensitive element names such as + +:: + + mat.add_element('aluminium', 1.0) + Internally, OpenMC stores data on the atomic masses and natural abundances of all known isotopes and then uses this data to determine what isotopes should be added to the material. When the material is later exported to XML for use by the @@ -55,6 +61,20 @@ following would add 3.2% enriched uranium to a material:: In addition to U235 and U238, concentrations of U234 and U236 will be present and are determined through a correlation based on measured data. +It is also possible to perform enrichment of any element that is composed +of two naturally-occurring isotopes (e.g., Li or B) in terms of atomic percent. +To invoke this, provide the additional argument `enrichment_target` to +:meth:`Material.add_element`. For example the following would enrich B10 +to 30ao%:: + + mat.add_element('B', 1.0, enrichment=30.0, enrichment_target='B10') + +In order to enrich an isotope in terms of mass percent (wo%), provide the extra +argument `enrichment_type`. For example the following would enrich Li6 to 15wo%:: + + mat.add_element('Li', 1.0, enrichment=15.0, enrichment_target='Li6', + enrichment_type='wo') + Often, cross section libraries don't actually have all naturally-occurring isotopes for a given element. For example, in ENDF/B-VII.1, cross section evaluations are given for O16 and O17 but not for O18. If OpenMC is aware of @@ -135,6 +155,33 @@ attribute, e.g., :attr:`Material.temperature` or :attr:`Cell.temperature` attributes, respectively. +----------------- +Material Mixtures +----------------- + +In OpenMC it is possible to mix any number of materials to create a new material +with the correct nuclide composition and density. The +:meth:`Material.mix_materials` method takes a list of materials and +a list of their mixing fractions. Mixing fractions can be provided as atomic +fractions, weight fractions, or volume fractions. The fraction type +can be specified by passing 'ao', 'wo', or 'vo' as the third argument, respectively. +For example, assuming the required materials have already been defined, a MOX +material with 3% plutonium oxide by weight could be created using the following: + +:: + + mox = openmc.Material.mix_materials([uo2, puo2], [0.97, 0.03], 'wo') + +It should be noted that, if mixing fractions are specifed as atomic or weight +fractions, the supplied fractions should sum to one. If the fractions are specified +as volume fractions, and the sum of the fractions is less than one, then the remaining +fraction is set as void material. + +.. warning:: Materials with :math:`S(\alpha,\beta)` thermal scattering data + cannot be used in :meth:`Material.mix_materials`. However, thermal + scattering data can be added to a material created by + :meth:`Material.mix_materials`. + -------------------- Material Collections -------------------- diff --git a/examples/jupyter/pincell.ipynb b/examples/jupyter/pincell.ipynb index e50a4a322c..88fe0ca359 100644 --- a/examples/jupyter/pincell.ipynb +++ b/examples/jupyter/pincell.ipynb @@ -173,7 +173,7 @@ "name": "stderr", "output_type": "stream", "text": [ - "/home/romano/openmc/openmc/mixin.py:71: IDWarning: Another Material instance already exists with id=2.\n", + "/home/sam/openmc/openmc/openmc/mixin.py:71: IDWarning: Another Material instance already exists with id=2.\n", " warn(msg, IDWarning)\n" ] } @@ -435,6 +435,39 @@ "uo2_three.set_density('g/cc', 10.0)" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Mixtures\n", + "\n", + "In OpenMC it is also possible to define materials by mixing existing materials. For example, if we wanted to create MOX fuel out of a mixture of UO2 (97 wt%) and PuO2 (3 wt%) we could do the following:" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [], + "source": [ + "# Create PuO2 material\n", + "puo2 = openmc.Material()\n", + "puo2.add_nuclide('Pu239', 0.94)\n", + "puo2.add_nuclide('Pu240', 0.06)\n", + "puo2.add_nuclide('O16', 2.0)\n", + "puo2.set_density('g/cm3', 11.5)\n", + "\n", + "# Create the mixture\n", + "mox = openmc.Material.mix_materials([uo2, puo2], [0.97, 0.03], 'wo')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The 'wo' argument in the `mix_materials()` method specifies that the fractions are weight fractions. Materials can also be mixed by atomic and volume fractions with 'ao' and 'vo', respectively. For 'ao' and 'wo' the fractions must sum to one. For 'vo', if fractions do not sum to one, the remaining fraction is set as void." + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -456,7 +489,7 @@ }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 16, "metadata": {}, "outputs": [], "source": [ @@ -474,7 +507,7 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 17, "metadata": {}, "outputs": [], "source": [ @@ -491,7 +524,7 @@ }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 18, "metadata": {}, "outputs": [ { @@ -517,7 +550,7 @@ }, { "cell_type": "code", - "execution_count": 18, + "execution_count": 19, "metadata": {}, "outputs": [], "source": [ @@ -534,7 +567,7 @@ }, { "cell_type": "code", - "execution_count": 19, + "execution_count": 20, "metadata": {}, "outputs": [ { @@ -543,7 +576,7 @@ "(array([-1., -1., 0.]), array([1., 1., 1.]))" ] }, - "execution_count": 19, + "execution_count": 20, "metadata": {}, "output_type": "execute_result" } @@ -561,7 +594,7 @@ }, { "cell_type": "code", - "execution_count": 20, + "execution_count": 21, "metadata": {}, "outputs": [], "source": [ @@ -581,7 +614,7 @@ }, { "cell_type": "code", - "execution_count": 21, + "execution_count": 22, "metadata": {}, "outputs": [], "source": [ @@ -604,7 +637,7 @@ }, { "cell_type": "code", - "execution_count": 22, + "execution_count": 23, "metadata": {}, "outputs": [], "source": [ @@ -624,22 +657,22 @@ }, { "cell_type": "code", - "execution_count": 23, + "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "" + "" ] }, - "execution_count": 23, + "execution_count": 24, "metadata": {}, "output_type": "execute_result" }, { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAARUAAAD8CAYAAABZ0jAcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAE2FJREFUeJzt3XusZWV9xvHvUyhDY6PMMARHMFyUOGBoQI9oS2KtoqJpZqBiHU0rGAy1LW1T4wVCUhuUFvAPrKkVJmDFSxiUhnSMGMq1/qGDnEmRGYbgDAOtM0UYGaBpQOjgr3+s9+CaM3ufs89e7163/XySnbP3uuzzrrPe9Zx3rX35KSIwM8vl15pugJn1i0PFzLJyqJhZVg4VM8vKoWJmWTlUzCyrLKEi6SuSnpC0dch8SfqipB2S7pf0htK8cyVtT7dzc7THzJqTa6TyVeDMBea/Bzgh3S4AvgwgaQXwGeDNwGnAZyQtz9QmM2tAllCJiO8DexdYZC3wtShsAg6TtAp4N3BbROyNiKeA21g4nMys5Q6u6fccBfy09HhXmjZs+gEkXUAxyuFlL3vZG1evXj2ZltrInt71wkSe97CjD5nI89roNm/e/POIOGKcdesKlcoiYj2wHmBmZiZmZ2cbblH/3fypR5tuwkBnX3ls003oPUn/Oe66dYXKbuDVpcdHp2m7gbfNm353TW2ykrYGyCCD2uqgaY+6QmUjcKGkDRQXZZ+JiMck3Qr8Xeni7LuAi2tq01TrUoiMYv72OGSakyVUJN1AMeJYKWkXxSs6vw4QEVcDtwDvBXYAzwIfSfP2SvoscG96qksjYqELvjaGvgXIKDyaaY66+NUHvqayuGkMkqVwwCxM0uaImBlnXb+j1syy6syrPzYaj1BGM/d38oglP4dKDzhIxlf+2zlg8nCodJSDJD8HTB4OlQ5xkNTHATM+h0oHOEya5esvS+NQaTGHSbs4XEbjl5RbyoHSXt43C/NIpUXcWbvDo5bhHCot4DDpLl/QPZBPfxrmQOkP78uCRyoNcQfsJ58WeaRiZpl5pFIzj1CmwzSPWDxSqZEDZfpM4z53qNRkGjuXFaZt3/v0Z8KmrUPZYNN0OuRQmRCHiQ0yDeGSq+zpmZIeSmVNLxow/ypJ96XbTyQ9XZr3YmnexhztaZoDxRbT5z5SeaQi6SDgS8A7KYqB3StpY0Rsm1smIv66tPxfAKeWnuK5iDilajvaos+dxfK6+VOP9nLEkuP05zRgR0TsBEhlONYC24Ys/0GKb9vvFYeJjaOPp0M5Tn+WUrr0GOA44M7S5EMlzUraJOmsDO2pnQPFqupTH6r7JeV1wE0R8WJp2jGpFMCHgC9Ies2gFSVdkMJnds+ePXW0dSR96gzWrL70pRyhMqyk6SDrgBvKEyJid/q5k6Lk6akHrlbUUo6ImYiYOeKIsepGZ9eXTmDt0Yc+lSNU7gVOkHScpEMoguOAV3EkrQaWAz8sTVsuaVm6vxI4neHXYsysAyqHSkTsAy4EbgUeBL4VEQ9IulTSmtKi64ANsX9JxBOBWUk/Bu4CLi+/atRmffiPYu3U9b7lsqdj6PpOt25o8hUhlz2tkQPF6tLVvua36Y+oqzvYuq2L72PxSGUEDhRrWpf6oENlEV3amdZvXemLDpUFdGUn2vToQp90qJhZVg6VIbrwH8GmU9v7pkNlgLbvNLM291GHyjxt3llmZW3tqw4VM8vKoVLS1uQ3G6aNfdahYmZZOVSSNia+2Sja1ncdKrRvp5gtVZv68NSHSpt2hlkVbenLUx8qZpbX1H71QVtS3SynNnxVgkcqZpbVVIaKRynWd0328bpqKZ8naU+pZvJHS/POlbQ93c7N0Z6FOFBsWjTV12uppZzcGBEXzlt3BUUJ1BkggM1p3aeqtsvMmpFjpPJSLeWIeAGYq6U8incDt0XE3hQktwFnZmjTQB6l2LRpos/XWUv5fZLul3STpLmKhkupw9zKsqdmtr+6LtR+Bzg2In6LYjRy/VKfoI1lT83sQLXUUo6IJyPi+fTwWuCNo66bi099bFrV3fdrqaUsaVXp4RqK8qhQlEp9V6qpvBx4V5qWlQPFpl2dx0DlV38iYp+kuVrKBwFfmaulDMxGxEbgL1Nd5X3AXuC8tO5eSZ+lCCaASyNib9U2mVlzel9L2aMUs18Z9e37rqVsZq3R61DxKMVsf3UcE70OFTOrn0PFzLLqbaj41MdssEkfG70NFTNrhkPFzLLqZaj41MdsYZM8RnoZKmbWHIeKmWXVu1DxqY/ZaCZ1rPQuVMysWb2p++MRitnSTaJOkEcqZpaVQ8XMsupFqPjUx6yanMdQL0LFzNrDoWJmWdVV9vTjkraluj93SDqmNO/FUjnUjfPXNbNuqavs6X8AMxHxrKQ/Ba4EPpDmPRcRp1Rth5m1Qy1lTyPiroh4Nj3cRFHfJwtfpDXLI9exVGfZ0znnA98rPT40lTPdJOmsYSu57KlZN9T6jlpJfwTMAL9bmnxMROyWdDxwp6QtEfHw/HUjYj2wHooSHbU02MyWrJaypwCSzgAuAdaUSqASEbvTz53A3cCpGdpkZg2pq+zpqcA1FIHyRGn6cknL0v2VwOlA+QLvgnw9xSyvHMdUXWVPPw/8JvBtSQD/FRFrgBOBayT9kiLgLp/3qpGZdUyWayoRcQtwy7xpf1O6f8aQ9X4AnJyjDWbWDn5HrZll5VAxs6w6Gyq+SGs2GVWPrc6Gipm1k0PFzLJyqJhZVg4VM8vKoWJmWTlUzCwrh4qZZdXJUHl61wtNN8Gs115z5MlvHHfdToaKmbWXQ8XMsnKomFlWDhUzy8qhYmZZOVTMLCuHipllVVfZ02WSbkzz75F0bGnexWn6Q5LenaM9ZtacyqFSKnv6HuAk4IOSTpq32PnAUxHxWuAq4Iq07kkU377/euBM4J/S85lZR9VS9jQ9vj7dvwl4h4qv1V8LbIiI5yPiEWBHej4z66i6yp6+tExE7AOeAQ4fcV1g/7Kn//PckxmabWaT0JkLtRGxPiJmImLm5b9xeNPNMbMh6ip7+tIykg4GXgE8OeK6ZtYhtZQ9TY/PTffPAe6MiEjT16VXh44DTgB+lKFNZtaQusqeXgd8XdIOYC9F8JCW+xZF/eR9wJ9HxItV22Rmzamr7OkvgPcPWfcy4LIc7TCz5nXmQq2ZdYNDxcyycqiYWVYOFTPLyqFiZll1MlQOO/qQpptg1msPP75l87jrdjJUzKy9HCpmlpVDxcyycqiYWVYOFTPLyqFiZlk5VMwsq86GytlXHtt0E8x6qeqx1dlQMbN2cqiYWVYOFTPLyqFiZllVChVJKyTdJml7+rl8wDKnSPqhpAck3S/pA6V5X5X0iKT70u2Upfx+X6w1yyvHMVV1pHIRcEdEnADckR7P9yzw4YiYK236BUmHleZ/MiJOSbf7KrbHzBpWNVTK5UyvB86av0BE/CQitqf7/w08ARxR8feaWUtVDZUjI+KxdP9nwJELLSzpNOAQ4OHS5MvSadFVkpYtsO5LZU/37NlTsdlmNimLhoqk2yVtHXDbrwh7Kg4WCzzPKuDrwEci4pdp8sXAauBNwArg08PWL5c9PeKIXw10fF3FLI9cx9KidX8i4oxh8yQ9LmlVRDyWQuOJIcu9HPgucElEbCo999wo53lJ/wx8YkmtN7PWqXr6Uy5nei7wr/MXSKVQbwa+FhE3zZu3Kv0UxfWYrRXbY2YNqxoqlwPvlLQdOCM9RtKMpGvTMn8IvBU4b8BLx9+UtAXYAqwEPlexPWbWsEplTyPiSeAdA6bPAh9N978BfGPI+m+v8vvNrH168Y5aX6w1qybnMdSLUDGz9nComFlWla6ptMnc8O3mTz3aaDvMumQSlw48UjGzrHoXKr5oazaaSR0rvQsVM2uWQ8XMsuplqPgUyGxhkzxGehkqZtYch4qZZdXbUPEpkNlgkz42ehsqZtYMh4qZZdXrUPEpkNn+6jgmeh0qZla/3oeKRytmhbqOhd6HipnVa+JlT9NyL5a+n3Zjafpxku6RtEPSjelLsrPzaMWmXZ3HQB1lTwGeK5U2XVOafgVwVUS8FngKOL9ie4ZysNi0qrvvT7zs6TCpLMfbgbmyHUta38zaqa6yp4emkqWbJM0Fx+HA0xGxLz3eBRw17BflKHvq0YpNmyb6/KJfJynpduCVA2ZdUn4QESFpWNnTYyJit6TjgTtTrZ9nltLQiFgPrAeYmZkZWl7VzJpVS9nTiNidfu6UdDdwKvAvwGGSDk6jlaOB3WNsg5m1SB1lT5dLWpburwROB7algu53AecstH5uPgWyadFUX6+j7OmJwKykH1OEyOURsS3N+zTwcUk7KK6xXFexPSNxsFjfNdnH6yh7+gPg5CHr7wROq9IGM2uX3tT9WSrXCbI+asMo3G/TN7Ospj5U2pDsZjm0pS9PfahAe3aG2bja1IcdKkmbdorZUrSt7zpUzCwrh0pJ2xLfbDFt7LMOFTPLyqEyTxuT32yQtvZVh8oAbd1ZZnPa3EcdKkO0eafZdGt733SomFlWDpUFtP0/gk2fLvRJh8oiurATbTp0pS86VEbQlZ1p/dWlPji1X32wVP6qBGtCl8JkjkcqS9TFnWzd1NW+5lAZQ1d3tnVHl/vYxMueSvq9UsnT+yT9Yq72j6SvSnqkNO+UKu2pU5d3urVb1/vWxMueRsRdcyVPKSoSPgv8W2mRT5ZKot5XsT1m1rC6y56eA3wvIp6t+Htboev/Uax9+tCn6ip7OmcdcMO8aZdJul/SVXP1gbqkD53A2qEvfWnRUJF0u6StA25ry8ul4mBDy5GmCoYnA7eWJl8MrAbeBKygqAM0bP3KtZQnpS+dwZrTpz6kIgvGXFl6CHhbqezp3RHxuiHL/hXw+oi4YMj8twGfiIjfX+z3zszMxOzs7NjtniS/j8WWoq1hImlzRMyMs+7Ey56WfJB5pz4piJAkiusxWyu2p3Ft7STWPn3tK3WUPUXSscCrgX+ft/43JW0BtgArgc9VbE8r9LWzWD597iOVTn+a0ubTn/l8OmRlXQmTKqc//uzPhPkzQwbdCZMc/Db9mkxTp7L9Tdu+d6jUaNo6l03nPvfpT818OjQdpjFM5nikYmZZeaTSEI9Y+mmaRyhzPFJpmDthf3hfFjxSaQGPWrrLQXIgh0qLlDuoA6bdHCbD+fSnpdxp28v7ZmEeqbSYT4vaxWEyGodKBzhcmuUwWRqHSof4mkt9HCTjc6h0lAMmPwdJHg6VHnDAjM9Bkp9DpWd8/WU0DpPJ8UvKZpaVRyo95VOiA3l0Ug+HyhQYdDD1PWgcIM2pFCqS3g/8LXAicFpEDPziWElnAv8AHARcGxFzX5B9HLABOBzYDPxxRLxQpU02mvkHXddDxiHSHlVHKluBPwCuGbaApIOALwHvBHYB90raGBHbgCuAqyJig6SrgfOBL1dsk42hS6MZB0i7VQqViHgQoCjbM9RpwI6I2JmW3QCslfQgRcH2D6XlrqcY9ThUWmKxg3dSoePQ6LY6rqkcBfy09HgX8GaKU56nI2JfafpRw55E0gXAXHXD5yV1vvDYACuBnzfdiAkZfds+P9mGZNbXfTaw0ugoFg0VSbcDrxww65KIWKgiYVYRsR5Yn9o0O25Nkjbr63ZBf7etz9s17rqLhkpEnDHukye7KaoTzjk6TXsSOEzSwWm0MjfdzDqsjje/3QucIOk4SYcA64CNUZRGvAs4Jy23WC1mM+uASqEi6WxJu4DfBr4r6dY0/VWSbgFIo5ALgVuBB4FvRcQD6Sk+DXxc0g6KayzXjfir11dpd4v1dbugv9vm7Zqnk7WUzay9/NkfM8vKoWJmWXUiVCS9X9IDkn4paejLd5LOlPSQpB2SLqqzjeOQtELSbZK2p5/Lhyz3oqT70m1j3e0c1WJ/f0nLJN2Y5t8j6dj6WzmeEbbtPEl7Svvpo020cykkfUXSE8Pe86XCF9M23y/pDSM9cUS0/kbx2aLXAXcDM0OWOQh4GDgeOAT4MXBS021fZLuuBC5K9y8Crhiy3P823dYRtmXRvz/wZ8DV6f464Mam251x284D/rHpti5xu94KvAHYOmT+e4HvAQLeAtwzyvN2YqQSEQ9GxEOLLPbSxwGi+FDiBmDt5FtXyVqKjyeQfp7VYFuqGuXvX97em4B3aJHPeLREF/vWoiLi+8DeBRZZC3wtCpso3le2arHn7USojGjQxwGGvu2/JY6MiMfS/Z8BRw5Z7lBJs5I2SWpr8Izy939pmSjeavAMxVsJ2m7UvvW+dJpwk6RXD5jfNWMdU635PpW2fBwgt4W2q/wgIkLSsNf3j4mI3ZKOB+6UtCUiHs7dVqvkO8ANEfG8pD+hGJG9veE2NaI1oRKT+zhAoxbaLkmPS1oVEY+lYeUTQ55jd/q5U9LdwKkU5/htMsrff26ZXZIOBl5B8XGNtlt02yKivB3XUlwv67qxjqk+nf4M/DhAw21azEaKjyfAkI8pSFouaVm6vxI4HdhWWwtHN8rfv7y95wB3Rroi2HKLbtu8aw1rKN493nUbgQ+nV4HeAjxTOl0frukr0CNepT6b4nzueeBx4NY0/VXALfOuVv+E4r/4JU23e4TtOhy4A9gO3A6sSNNnKL4hD+B3gC0UrzhsAc5vut0LbM8Bf3/gUmBNun8o8G1gB/Aj4Pim25xx2/4eeCDtp7uA1U23eYRtugF4DPi/dHydD3wM+FiaL4ovWHs49b2Br7zOv/lt+maWVZ9Of8ysBRwqZpaVQ8XMsnKomFlWDhUzy8qhYmZZOVTMLKv/BxrmRHimLrvyAAAAAElFTkSuQmCC\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAARUAAAD8CAYAAABZ0jAcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAATXElEQVR4nO3dfbBcdX3H8fenUEJrp5IQBiMqkIoGqh3AK9oyY1tFRP9IsKJGpxUcHGpb2pk6OsIwUzuoLdg/sE6tmsEHfBiC0nEaRxzksf6hQW4skgCFhDjVpCiRIJ1OMJrw7R/nt/HkZnfv3t3fnqf9vGZ27u552Ps79/zO5/7O2YevIgIzs1x+re4GmFm3OFTMLCuHipll5VAxs6wcKmaWlUPFzLLKEiqSPiPpcUnbBsyXpI9J2iHpfklnl+ZdLGl7ul2coz1mVp9cI5XPARcMmf964LR0uwz4BICkFcAHgFcA5wAfkLQ8U5vMrAZZQiUivgXsHbLIOuDzUdgMHCdpFfA64LaI2BsRTwK3MTyczKzhjq7o95wE/Kj0eFeaNmj6ESRdRjHK4VnPetbL1qxZM52W2sh2bntwKs+7+iVnTOV5bXRbtmz5aUScMM66VYXKxCJiA7ABYG5uLubn52tuUfetf9HZQ+e/8AVTCvb/fWbo7I2PfG86v9cOkfTf465bVajsBp5fevy8NG038EcLpt9dUZusZLEAaZJ+bXXQNEdVobIJuFzSRoqLsk9FxGOSbgX+oXRx9nzgyoraNNPaFCKjWLg9Dpn6ZAkVSTdSjDhWStpF8YrOrwNExCeBW4A3ADuAfcA707y9kj4I3Jue6uqIGHbB18bQtQAZhUcz9VEbv/rA11QWN4tBshQOmOEkbYmIuXHW9TtqzSyr1rz6Y6PxCGU0vb+TRyz5OVQ6wEEyvvLfzgGTh0OlpRwk+Tlg8nCotIiDpDoOmPE5VFrAYVIvX39ZGodKgzlMmsXhMhq/pNxQDpTm8r4ZziOVBnFnbQ+PWgZzqDSAw6S9fEH3SD79qZkDpTu8LwseqdTEHbCbfFrkkYqZZeaRSsU8QpkNszxi8UilQg6U2TOL+9yhUpFZ7FxWmLV979OfKZu1DmX9zdLpkENlShwm1s8shEuusqcXSHo4lTW9os/86yTdl26PSPpZad7B0rxNOdpTNweKLabLfWTikYqko4CPA6+lKAZ2r6RNEXGo0lRE/G1p+b8Gzio9xdMRceak7WiKLncWy2v9i87u5Iglx+nPOcCOiNgJkMpwrAMGla97G8W37XeKw8TG0cXToRynP0spXXoycCpwZ2nysZLmJW2WdGGG9lTOgWKT6lIfqvol5fXAzRFxsDTt5FQK4O3ARyX9Tr8VJV2Wwmd+z549VbR1JF3qDFavrvSlHKEyqKRpP+uBG8sTImJ3+rmTouTpWUeuVtRSjoi5iJg74YSx6kZn15VOYM3RhT6VI1TuBU6TdKqkYyiC44hXcSStAZYD3ylNWy5pWbq/EjiXwddizKwFJg6ViDgAXA7cCjwEfDkiHpB0taS1pUXXAxvj8JKIpwPzkr4P3AVcU37VqMm68B/FmqntfctlT8fQ9p1u7VDnK0Iue1ohB4pVpa19zW/TH1Fbd7C1Wxvfx+KRyggcKFa3NvVBh8oi2rQzrdva0hcdKkO0ZSfa7GhDn3SomFlWDpUB2vAfwWZT0/umQ6WPpu80syb3UYfKAk3eWWZlTe2rDhUzy8qhUtLU5DcbpIl91qFiZlk5VJImJr7ZKJrWdx0qNG+nmC1Vk/rwzIdKk3aG2SSa0pdnPlTMLK+Z/eqDpqS6WU5N+KoEj1TMLKuZDBWPUqzr6uzjVdVSvkTSnlLN5HeV5l0saXu6XZyjPcM4UGxW1NXXK6mlnNwUEZcvWHcFRQnUOSCALWndJydtl5nVI8dI5VAt5Yj4BdCrpTyK1wG3RcTeFCS3ARdkaFNfHqXYrKmjz1dZS/lNku6XdLOkXkXDpdRhbmTZUzM7XFUXar8GnBIRv0cxGrlhqU/QxLKnZnakSmopR8QTEbE/PbweeNmo6+biUx+bVVX3/UpqKUtaVXq4lqI8KhSlUs9PNZWXA+enaVk5UGzWVXkMTPzqT0QckNSrpXwU8JleLWVgPiI2AX+T6iofAPYCl6R190r6IEUwAVwdEXsnbZOZ1afztZQ9SjH7lVHfvu9aymbWGJ0OFY9SzA5XxTHR6VAxs+o5VMwsq86Gik99zPqb9rHR2VAxs3o4VMwsq06Gik99zIab5jHSyVAxs/o4VMwsq86Fik99zEYzrWOlc6FiZvXqTN0fj1DMlm4adYI8UjGzrBwqZpZVJ0LFpz5mk8l5DHUiVMysORwqZpZVVWVP3yPpwVT35w5JJ5fmHSyVQ920cF0za5eqyp7+JzAXEfsk/QXwEeCtad7TEXHmpO0ws2aopOxpRNwVEfvSw80U9X2y8EVaszxyHUtVlj3tuRT4Runxsamc6WZJFw5ayWVPzdqh0nfUSvpTYA74w9LkkyNit6TVwJ2StkbEowvXjYgNwAYoSnRU0mAzW7JKyp4CSDoPuApYWyqBSkTsTj93AncDZ2Vok5nVpKqyp2cBn6IIlMdL05dLWpburwTOBcoXeIfy9RSzvHIcU1WVPf0n4LeAr0gC+GFErAVOBz4l6RmKgLtmwatGZtYyWa6pRMQtwC0Lpv1d6f55A9b7NvDSHG0ws2bwO2rNLCuHipll1dpQ8UVas+mY9NhqbaiYWTM5VMwsK4eKmWXlUDGzrBwqZpaVQ8XMsnKomFlWrQyVndv88SCzaVq+7DdeNu66rQwVM2suh4qZZeVQMbOsHCpmlpVDxcyycqiYWVYOFTPLqqqyp8sk3ZTm3yPplNK8K9P0hyW9Lkd7zKw+E4dKqezp64EzgLdJOmPBYpcCT0bEC4HrgGvTumdQfPv+7wIXAP+ans/MWqqSsqfp8Q3p/s3Aa1R8rf46YGNE7I+IHwA70vOZWUtVVfb00DIRcQB4Cjh+xHWBw8ue7j94IEOzzWwaWnOhNiI2RMRcRMwtO6rSaq1mtgRVlT09tIyko4FnA0+MuK6ZtUglZU/T44vT/YuAOyMi0vT16dWhU4HTgO9maJOZ1aSqsqefBr4gaQewlyJ4SMt9maJ+8gHgryLi4KRtMrP6qBgwtMuKY38zzn/BmrqbYdZZ3/zhf7H35/s0zrqtuVBrZu3gUDGzrBwqZpaVQ8XMsnKomFlWrQyV1S9Z+HlFM8vpyf1Pbxl33VaGipk1l0PFzLJyqJhZVg4VM8vKoWJmWTlUzCwrh4qZZdXaUNn4yPfqboJZJ016bLU2VMysmRwqZpaVQ8XMsnKomFlWE4WKpBWSbpO0Pf1c3meZMyV9R9IDku6X9NbSvM9J+oGk+9LtzKX8fl+sNcsrxzE16UjlCuCOiDgNuCM9Xmgf8I6I6JU2/aik40rz3xcRZ6bbfRO2x8xqNmmolMuZ3gBcuHCBiHgkIran+/8DPA6cMOHvNbOGmjRUToyIx9L9HwMnDltY0jnAMcCjpckfTqdF10laNmTdQ2VP9+zZM2GzzWxaFg0VSbdL2tbndlgR9lQcbGC9D0mrgC8A74yIZ9LkK4E1wMuBFcD7B61fLnt6wgm/Guj4uopZHrmOpUWLiUXEeYPmSfqJpFUR8VgKjccHLPfbwNeBqyJic+m5e6Oc/ZI+C7x3Sa03s8aZ9PSnXM70YuDfFy6QSqF+Ffh8RNy8YN6q9FMU12O2TdgeM6vZpKFyDfBaSduB89JjJM1Juj4t8xbgVcAlfV46/pKkrcBWYCXwoQnbY2Y1m6iWckQ8Abymz/R54F3p/heBLw5Y/9WT/H4za55OvKPWF2vNJpPzGOpEqJhZczhUzCyria6pNElv+Lb+RWfX3BKz9pjGpQOPVMwsq86Fii/amo1mWsdK50LFzOrlUDGzrDoZKj4FMhtumsdIJ0PFzOrjUDGzrDobKj4FMutv2sdGZ0PFzOrhUDGzrDodKj4FMjtcFcdEp0PFzKrX+VDxaMWsUNWx0PlQMbNqTb3saVruYOn7aTeVpp8q6R5JOyTdlL4kOzuPVmzWVXkMVFH2FODpUmnTtaXp1wLXRcQLgSeBSydsz0AOFptVVff9qZc9HSSV5Xg10CvbsaT1zayZqip7emwqWbpZUi84jgd+FhEH0uNdwEmDflGOsqcerdisqaPPL/p1kpJuB57TZ9ZV5QcREZIGlT09OSJ2S1oN3Jlq/Ty1lIZGxAZgA8Dc3NzA8qpmVq9Kyp5GxO70c6eku4GzgH8DjpN0dBqtPA/YPcY2mFmDVFH2dLmkZen+SuBc4MFU0P0u4KJh6+fmUyCbFXX19SrKnp4OzEv6PkWIXBMRD6Z57wfeI2kHxTWWT0/YnpE4WKzr6uzjVZQ9/Tbw0gHr7wTOmaQNZtYsnan7s1SuE2Rd1IRRuN+mb2ZZzXyoNCHZzXJoSl+e+VCB5uwMs3E1qQ87VJIm7RSzpWha33WomFlWDpWSpiW+2WKa2GcdKmaWlUNlgSYmv1k/Te2rDpU+mrqzzHqa3EcdKgM0eafZbGt633SomFlWDpUhmv4fwWZPG/qkQ2URbdiJNhva0hcdKiNoy8607mpTH5zZrz5YKn9VgtWhTWHS45HKErVxJ1s7tbWvOVTG0Nadbe3R5j429bKnkv64VPL0Pkk/79X+kfQ5ST8ozTtzkvZUqc073Zqt7X1r6mVPI+KuXslTioqE+4BvlhZ5X6kk6n0TtsfMalZ12dOLgG9ExL4Jf28jtP0/ijVPF/pUVWVPe9YDNy6Y9mFJ90u6rlcfqE260AmsGbrSlxYNFUm3S9rW57auvFwqDjawHGmqYPhS4NbS5CuBNcDLgRUUdYAGrT9xLeVp6UpnsPp0qQ9VUvY0eQvw1Yj4Zem5e6Oc/ZI+C7x3SDsaXUvZ72OxcXQpTHqmXva05G0sOPVJQYQkUVyP2TZhe2rXxU5i09HVvlJF2VMknQI8H/iPBet/SdJWYCuwEvjQhO1phK52Fsuny31ExaWQdpmbm4v5+fm6mzESnw5ZWVvCRNKWiJgbZ11/9mfKfK3FoD1hkoPfpl+RWepUdrhZ2/cOlQrNWuey2dznPv2pmE+HZsMshkmPRypmlpVHKjXxiKWbZnmE0uORSs3cCbvD+7LgkUoDeNTSXg6SIzlUGqTcQR0wzeYwGcynPw3lTttc3jfDeaTSYD4tahaHyWgcKi3gcKmXw2RpHCot4msu1XGQjM+h0lIOmPwcJHk4VDrAATM+B0l+DpWO8fWX0ThMpscvKZtZVh6pdJRPiY7k0Uk1HCozoN/B1PWgcYDUZ6JQkfRm4O+B04FzIqLvF8dKugD4Z+Ao4PqI6H1B9qnARuB4YAvwZxHxi0naZKNZeNC1PWQcIs0x6UhlG/AnwKcGLSDpKODjwGuBXcC9kjZFxIPAtcB1EbFR0ieBS4FPTNgmG0ObRjMOkGabKFQi4iGAomzPQOcAOyJiZ1p2I7BO0kMUBdvfnpa7gWLU41BpiMUO3mmFjkOj3aq4pnIS8KPS413AKyhOeX4WEQdK008a9CSSLgMuSw/3S2p94bE+VgI/rbsRUzLytt00/J9U03R1n7143BUXDRVJtwPP6TPrqogYVpEwq3LZU0nz49YkabKubhd0d9u6vF3jrjtRLeUR7aaoTtjzvDTtCeA4SUen0Upvupm1WBVvfrsXOE3SqZKOAdYDm6IojXgXcFFabrFazGbWAhOFiqQ3StoF/D7wdUm3punPlXQLQBqFXA7cCjwEfDkiHkhP8X7gPZJ2UFxj+fSIv3rDJO1usK5uF3R327xdC7SylrKZNZc/+2NmWTlUzCyrVoSKpDdLekDSM5IGvnwn6QJJD0vaIemKKts4DkkrJN0maXv6uXzAcgcl3Zdum6pu56gW+/tLWibppjT/HkmnVN/K8YywbZdI2lPaT++qo51LIekzkh4f9J4vFT6Wtvl+SaO92zEiGn+j+GzRi4G7gbkByxwFPAqsBo4Bvg+cUXfbF9mujwBXpPtXANcOWO7/6m7rCNuy6N8f+Evgk+n+euCmutudcdsuAf6l7rYucbteBZwNbBsw/w3ANwABrwTuGeV5WzFSiYiHIuLhRRY79HGAKD6UuBFYN/3WTWQdxccTSD8vrLEtkxrl71/e3puB12iRz3g0RBv71qIi4lvA3iGLrAM+H4XNFO8rW7XY87YiVEbU7+MAA9/23xAnRsRj6f6PgRMHLHespHlJmyU1NXhG+fsfWiaKtxo8RfFWgqYbtW+9KZ0m3Czp+X3mt81Yx1Rjvk+lKR8HyG3YdpUfRERIGvT6/skRsVvSauBOSVsj4tHcbbWJfA24MSL2S/pzihHZq2tuUy0aEyoxvY8D1GrYdkn6iaRVEfFYGlY+PuA5dqefOyXdDZxFcY7fJKP8/XvL7JJ0NPBsio9rNN2i2xYR5e24nuJ6WduNdUx16fSn78cBam7TYjZRfDwBBnxMQdJyScvS/ZXAucCDlbVwdKP8/cvbexFwZ6Qrgg236LYtuNawluLd4223CXhHehXolcBTpdP1weq+Aj3iVeo3UpzP7Qd+Atyapj8XuGXB1epHKP6LX1V3u0fYruOBO4DtwO3AijR9juIb8gD+ANhK8YrDVuDSuts9ZHuO+PsDVwNr0/1jga8AO4DvAqvrbnPGbftH4IG0n+4C1tTd5hG26UbgMeCX6fi6FHg38O40XxRfsPZo6nt9X3ldePPb9M0sqy6d/phZAzhUzCwrh4qZZeVQMbOsHCpmlpVDxcyycqiYWVb/D5MtR6Xki1ByAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] @@ -663,22 +696,22 @@ }, { "cell_type": "code", - "execution_count": 24, + "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "" + "" ] }, - "execution_count": 24, + "execution_count": 25, "metadata": {}, "output_type": "execute_result" }, { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAARUAAAD8CAYAAABZ0jAcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAEW9JREFUeJzt3X2MXNV9xvHvU6iNRNVgY4uaFxsTUI2rtAa2kBYpbYEEyB82aWhiR1VMZOSmLapUlAgQEkhO0kL6h6O0aRMLCBBFBuI2qqOAXINNkdKYsFYNNhCwMW1j18EOBqTK1InNr3/cs/Syntmd3Tl7X2afjzSamfsye+7MmWfPvXdmfooIzMxy+aW6G2Bmg8WhYmZZOVTMLCuHipll5VAxs6wcKmaWVZZQkXSfpIOSdnWZL0lflbRH0nOSLi7NWylpd7qszNEeM6tPrpHK/cA1Y8y/FrggXVYD/wAgaTZwJ3AZcClwp6RZmdpkZjXIEioR8RRweIxFlgEPRmEbcJqkecDVwOaIOBwRbwCbGTuczKzhTq7o75wF/KR0f1+a1m36CSStphjlcOqpp16yaNGiqWmp9ezwgR9PyePOnufXtm7bt2//WUTMncy6VYVK3yJiHbAOYGhoKIaHh2tu0eBbv+bycZb4rUraMdqKO35Qy9+dTiT952TXrSpU9gPnlO6fnabtB35/1PQnK2qTlYwfIM3Rqa0OmuaoKlQ2AjdJeojioOxbEXFA0ibgr0oHZz8C3FZRm6a1NoVIL0Zvj0OmPllCRdJ6ihHHHEn7KM7o/DJARHwdeBT4KLAHOAJ8Js07LOkLwDPpodZExFgHfG0SBi1AeuHRTH3Uxp8+8DGV8U3HIJkIB8zYJG2PiKHJrOtP1JpZVq05+2O98QilNyPPk0cs+TlUBoCDZPLKz50DJg+HSks5SPJzwOThUGkRB0l1HDCT51BpAYdJvXz8ZWIcKg3mMGkWh0tvfEq5oRwozeXXZmweqTSIO2t7eNTSnUOlARwm7eUDuify7k/NHCiDw69lwSOVmrgDDibvFnmkYmaZeaRSMY9QpofpPGLxSKVCDpTpZzq+5g6VikzHzmWF6fbae/dnik23DmWdTafdIYfKFHGYWCfTIVxylT29RtJLqazprR3mr5W0I11elvRmad7x0ryNOdpTNweKjWeQ+0jfoSLpJOBrFKVNFwMrJC0uLxMRfxkRSyJiCfC3wD+VZr89Mi8ilvbbnroNcmexvAa1r+TY/bkU2BMRewFSGY5lwAtdll9B8Wv7A2VQO4hNrUHcHcqx+zOR0qULgIXAltLkUyQNS9om6boM7amcA8X6NUh9qOpTysuBDRFxvDRtQSoF8CngK5Le32lFSatT+AwfOnSoirb2ZJA6g9VrUPpSjlDpVtK0k+XA+vKEiNifrvdSlDy9qNOKEbEuIoYiYmju3EnVjc5uUDqBNccg9KkcofIMcIGkhZJmUATHCWdxJC0CZgE/LE2bJWlmuj0HuJzux2LMrAX6DpWIOAbcBGwCXgQeiYjnJa2RVD6bsxx4KN5bEvFCYFjSs8BW4K6IaEWoDMJ/FGumtvctlz2dhLa/6NYOdZ4RctnTCjlQrCpt7Wv+mH6P2voCW7u18XMsHqn0wIFidWtTH3SojKNNL6YNtrb0RYfKGNryItr00YY+6VAxs6wcKl204T+CTU9N75sOlQ6a/qKZNbmPOlRGafKLZVbW1L7qUDGzrBwqJU1NfrNumthnHSpmlpVDJWli4pv1oml916FC814Us4lqUh+e9qHSpBfDrB9N6cvTPlTMLK9p+9MHTUl1s5ya8FMJHqmYWVbTMlQ8SrFBV2cfr6qW8g2SDpVqJt9YmrdS0u50WZmjPWNxoNh0UVdf7/uYSqmW8ocpqhM+I2ljh1/Ffzgibhq17myKEqhDQADb07pv9NsuM6tHjpHKu7WUI+LnwEgt5V5cDWyOiMMpSDYD12RoU0cepdh0U0efr7KW8sclPSdpg6SRioYTqcPcyLKnZvZeVR2o/R5wbkT8JsVo5IGJPkATy56a2YkqqaUcEa9HxNF09x7gkl7XzcW7PjZdVd33K6mlLGle6e5SivKoUJRK/UiqqTwL+EialpUDxaa7Kt8DfZ/9iYhjkkZqKZ8E3DdSSxkYjoiNwF+kusrHgMPADWndw5K+QBFMAGsi4nC/bTKz+gx8LWWPUsz+X68f33ctZTNrjIEOFY9SzN6rivfEQIeKmVXPoWJmWQ1sqHjXx6yzqX5vDGyomFk9HCpmltVAhop3fczGNpXvkYEMFTOrj0PFzLIauFDxro9Zb6bqvTJwoWJm9RqYuj8eoZhN3FTUCfJIxcyycqiYWVYDESre9THrT8730ECEipk1h0PFzLKqquzpzZJeSHV/npC0oDTveKkc6sbR65pZu1RV9vTfgaGIOCLpT4EvA59M896OiCX9tsPMmqGSsqcRsTUijqS72yjq+2Thg7RmeeR6L1VZ9nTEKuCx0v1TUjnTbZKu67aSy56atUOln6iV9MfAEPB7pckLImK/pPOALZJ2RsQro9eNiHXAOihKdFTSYDObsErKngJIugq4HVhaKoFKROxP13uBJ4GLMrTJzGpSVdnTi4BvUATKwdL0WZJmpttzgMuB8gHeMfl4illeOd5TVZU9/RvgV4DvSAL4r4hYClwIfEPSOxQBd9eos0Zm1jJZjqlExKPAo6Om3VG6fVWX9f4N+ECONphZM/gTtWaWlUPFzLJqbaj4IK3Z1Oj3vdXaUDGzZnKomFlWDhUzy8qhYmZZOVTMLCuHipll5VAxs6xaGSqHD/y47iaYDbSFZ556yWTXbWWomFlzOVTMLCuHipll5VAxs6wcKmaWlUPFzLJyqJhZVlWVPZ0p6eE0/2lJ55bm3ZamvyTp6hztMbP69B0qpbKn1wKLgRWSFo9abBXwRkScD6wF7k7rLqb49f3fAK4B/j49npm1VCVlT9P9B9LtDcCVKn5WfxnwUEQcjYhXgT3p8cyspXL8mn6nsqeXdVsmlfR4Czg9Td82at2OJVMlrQZWA8yfP58Vd/wgQ9PNrJNP3antk123NQdqI2JdRAxFxNDcuXPrbo6ZdVFV2dN3l5F0MvA+4PUe1zWzFqmk7Gm6vzLdvh7YEhGRpi9PZ4cWAhcAP8rQJjOrSVVlT+8FviVpD3CYInhIyz1CUT/5GPDnEXG83zaZWX1UDBjaZWhoKIaHh+tuhtnAkrQ9IoYms25rDtSaWTs4VMwsK4eKmWXlUDGzrBwqZpaVQ8XMsnKomFlWDhUzy8qhYmZZOVTMLCuHipll5VAxs6wcKmaWlUPFzLJyqJhZVg4VM8vKoWJmWTlUzCyrvkJF0mxJmyXtTtezOiyzRNIPJT0v6TlJnyzNu1/Sq5J2pMuSftpjZvXrd6RyK/BERFwAPJHuj3YE+HREjJQ2/Yqk00rzPx8RS9JlR5/tMbOa9Rsq5XKmDwDXjV4gIl6OiN3p9n8DBwFXAzMbUP2GyhkRcSDd/ilwxlgLS7oUmAG8Upr8pbRbtFbSzDHWXS1pWNLwoUOH+my2mU2VcUNF0uOSdnW4vKcIeyoO1rXeh6R5wLeAz0TEO2nybcAi4LeB2cAt3dZ32VOzdhi3mFhEXNVtnqTXJM2LiAMpNA52We5Xge8Dt0fEuwXZS6Oco5K+CXxuQq03s8bpd/enXM50JfDPoxdIpVC/CzwYERtGzZuXrkVxPGZXn+0xs5r1Gyp3AR+WtBu4Kt1H0pCke9IynwA+BNzQ4dTxtyXtBHYCc4Av9tkeM6uZy56a2Qlc9tTMGsOhYmZZOVTMLCuHipll5VAxs6wcKmaWlUPFzLJyqJhZVg4VM8vKoWJmWTlUzCwrh4qZZeVQMbOsHCpmlpVDxcyycqiYWVYOFTPLyqFiZllNednTtNzx0u/TbixNXyjpaUl7JD2cfiTbzFqsirKnAG+XSpsuLU2/G1gbEecDbwCr+myPmdVsysuedpPKclwBjJTtmND6ZtZMVZU9PSWVLN0maSQ4TgfejIhj6f4+4Kxuf8hlT83aYdwKhZIeB36tw6zby3ciIiR1q/exICL2SzoP2JJq/bw1kYZGxDpgHRQlOiayrplVp5KypxGxP13vlfQkcBHwj8Bpkk5Oo5Wzgf2T2AYza5Aqyp7OkjQz3Z4DXA68kAq6bwWuH2t9M2uXKsqeXggMS3qWIkTuiogX0rxbgJsl7aE4xnJvn+0xs5q57KmZncBlT82sMRwqZpaVQ8XMsnKomFlWDhUzy8qhYmZZOVTMLCuHipll5VAxs6wcKmaWlUPFzLJyqJhZVg4VM8vKoWJmWTlUzCwrh4qZZeVQMbOsHCpmltWUlz2V9Aelkqc7JP3vSO0fSfdLerU0b0k/7TGz+k152dOI2DpS8pSiIuER4F9Ki3y+VBJ1R5/tMbOaVV329HrgsYg40uffNbOGqqrs6YjlwPpR074k6TlJa0fqA5lZe1VV9pRUwfADwKbS5NsowmgGRUnTW4A1XdZfDawGmD9//njNNrOaVFL2NPkE8N2I+EXpsUdGOUclfRP43BjtcC1lsxaY8rKnJSsYteuTgghJojges6vP9phZzaooe4qkc4FzgH8dtf63Je0EdgJzgC/22R4zq9m4uz9jiYjXgSs7TB8Gbizd/w/grA7LXdHP3zez5vEnas0sK4eKmWXlUDGzrBwqZpaVQ8XMsnKomFlWDhUzy8qhYmZZOVTMLCuHipll5VAxs6wcKmaWlUPFzLJyqJhZVg4VM8vKoWJmWTlUzCwrh4qZZeVQMbOs+q2l/EeSnpf0jqShMZa7RtJLkvZIurU0faGkp9P0hyXN6Kc9Zla/fkcqu4A/BJ7qtoCkk4CvAdcCi4EVkhan2XcDayPifOANYFWf7TGzmvUVKhHxYkS8NM5ilwJ7ImJvRPwceAhYlmr9XAFsSMv1UovZzBqurxIdPToL+Enp/j7gMuB04M2IOFaafkIZjxHlsqcUFQ0HsfDYHOBndTdiigzqtg3qdv36ZFfsq5ZyRIxVkTCrctlTScMR0fUYTlsN6nbB4G7bIG/XZNftq5Zyj/ZTVCcccXaa9jpwmqST02hlZLqZtVgVp5SfAS5IZ3pmAMuBjRERwFbg+rTceLWYzawF+j2l/DFJ+4DfAb4vaVOafqakRwHSKOQmYBPwIvBIRDyfHuIW4GZJeyiOsdzb459e10+7G2xQtwsGd9u8XaOoGDCYmeXhT9SaWVYOFTPLqhWh0u/XAZpK0mxJmyXtTtezuix3XNKOdNlYdTt7Nd7zL2lm+jrGnvT1jHOrb+Xk9LBtN0g6VHqdbqyjnRMh6T5JB7t95kuFr6Ztfk7SxT09cEQ0/gJcSPFhnCeBoS7LnAS8ApwHzACeBRbX3fZxtuvLwK3p9q3A3V2W+5+629rDtoz7/AN/Bnw93V4OPFx3uzNu2w3A39Xd1glu14eAi4FdXeZ/FHgMEPBB4OleHrcVI5Xo4+sAU9+6viyj+HoCtP9rCr08/+Xt3QBcmb6u0XRt7FvjioingMNjLLIMeDAK2yg+VzZvvMdtRaj0qNPXAbp+7L8hzoiIA+n2T4Ezuix3iqRhSdskNTV4enn+310mio8avEXxUYKm67VvfTztJmyQdE6H+W0zqfdUFd/96UlTvg6Q21jbVb4TESGp2/n9BRGxX9J5wBZJOyPildxttb58D1gfEUcl/QnFiOyKmttUi8aESkzd1wFqNdZ2SXpN0ryIOJCGlQe7PMb+dL1X0pPARRT7+E3Sy/M/ssw+SScD76P4ukbTjbttEVHejnsojpe13aTeU4O0+9Px6wA1t2k8Gym+ngBdvqYgaZakmen2HOBy4IXKWti7Xp7/8vZeD2yJdESw4cbdtlHHGpZSfHq87TYCn05ngT4IvFXaXe+u7iPQPR6l/hjF/txR4DVgU5p+JvDoqKPVL1P8F7+97nb3sF2nA08Au4HHgdlp+hBwT7r9u8BOijMOO4FVdbd7jO054fkH1gBL0+1TgO8Ae4AfAefV3eaM2/bXwPPpddoKLKq7zT1s03rgAPCL9P5aBXwW+GyaL4ofWHsl9b2OZ15HX/wxfTPLapB2f8ysARwqZpaVQ8XMsnKomFlWDhUzy8qhYmZZOVTMLKv/A7NFXJz1zhCeAAAAAElFTkSuQmCC\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAARUAAAD8CAYAAABZ0jAcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAARZ0lEQVR4nO3de4xc5X3G8e9TXBuJqsHGDnUAGxxQjau0NmwhLVLaAuGSP2zS0MRUVWxk5KYtrVREBAipkUjSQprKUdS0iQUEiCIucRXVUUCuwab80ZiwVg02RuDFNIldBzs2IFWmTmx+/eO8Sw/rnd3ZnXfPZfb5SKOZc5t9z8yZZ99zmfkpIjAzy+WX6m6AmfUXh4qZZeVQMbOsHCpmlpVDxcyycqiYWVZZQkXS/ZIOStrVYbokfVXSkKQXJF1UmrZK0p50W5WjPWZWn1w9lQeAa8aYfi1wQbqtBf4ZQNIc4HPApcAlwOckzc7UJjOrQZZQiYhngCNjzLICeCgK24DTJc0HrgY2R8SRiHgD2MzY4WRmDTejor9zFvCT0vC+NK7T+JNIWkvRy+G00067ePHixVPTUuva4ef3TcnznvFbZ0/J81r3tm/f/rOImDeZZasKlZ5FxHpgPcDAwEAMDg7W3KL+98D7bx17hqnaUf3x2JNXH/zyFP1hGybpR5NdtqpQ2Q+cUxo+O43bD/z+iPFPV9QmKxk3QBpktLY6aJqjqlDZCNws6RGKg7JvRcQBSZuAvy0dnL0KuKOiNk1rbQqRboxcH4dMfbKEiqSHKXoccyXtozij88sAEfF14HHgY8AQcBS4MU07IunzwHPpqe6KiLEO+Nok9FuAdMO9mfqojT994GMq45uOQTIRDpixSdoeEQOTWdZX1JpZVq05+2PdcQ+lO8Ovk3ss+TlU+oCDZPLKr50DJg+HSks5SPJzwOThUGkRB0l1HDCT51BpAYdJvXz8ZWIcKg3mMGkWh0t3fEq5oRwozeX3ZmzuqTSIN9b2cK+lM4dKAzhM2ssHdE/m3Z+aOVD6h9/LgnsqNfEG2J+8W+Seipll5p5KxdxDmR6mc4/FPZUKOVCmn+n4njtUKjIdNy4rTLf33rs/U2y6bVA2uum0O+RQmSIOExvNdAiXXGVPr5H0ciprevso09dJ2pFur0h6szTtRGnaxhztqZsDxcbTz9tIzz0VSacAXwM+SlEM7DlJGyNi9/A8EfHXpfn/ElhWeoq3I2Jpr+1oin7eWCyvB95/a1/2WHLs/lwCDEXEXoBUhmMFsLvD/DdQ/Np+X3GY2GT04+5Qjt2fiZQuXQicB2wpjT5V0qCkbZKuy9CeyjlQrFf9tA1VfUp5JbAhIk6Uxi1MpQD+GPiKpA+OtqCktSl8Bg8dOlRFW7vSTxuD1atftqUcodKppOloVgIPl0dExP50v5ei5OmykxcrailHxEBEDMybN6m60dn1y0ZgzdEP21SOUHkOuEDSeZJmUgTHSWdxJC2mKOn9g9K42ZJmpcdzgcvofCzGzFqg51CJiOPAzcAm4CXgsYh4UdJdkpaXZl0JPBLvLYl4ITAo6XlgK3B3+axRk/XDfxRrprZvWy57Ogltf9OtHeo8I+SypxVyoFhV2rqt+TL9LrX1DbZ2a+N1LO6pdMGBYnVr0zboUBlHm95M629t2RYdKmNoy5to00cbtkmHipll5VDpoA3/EWx6avq26VAZRdPfNLMmb6MOlRGa/GaZlTV1W3WomFlWDpWSpia/WSdN3GYdKmaWlUMlaWLim3WjaduuQ4XmvSlmE9WkbXjah0qT3gyzXjRlW572oWJmeU3bnz5oSqqb5dSEn0pwT8XMspqWoeJeivW7Orfxqmopr5Z0qFQz+abStFWS9qTbqhztGYsDxaaLurb1SmopJ49GxM0jlp1DUQJ1AAhge1r2jV7bZWb1yNFTebeWckT8HBiupdyNq4HNEXEkBclm4JoMbRqVeyk23dSxzVdZS/kTkl6QtEHScEXDidRhbmTZUzN7r6oO1H4PODcifpOiN/LgRJ+giWVPzexkldRSjojDEXEsDd4LXNztsrl418emq6q3/UpqKUuaXxpcTlEeFYpSqVelmsqzgavSuKwcKDbdVfkZ6PnsT0QclzRcS/kU4P7hWsrAYERsBP4q1VU+DhwBVqdlj0j6PEUwAdwVEUd6bZOZ1afvaym7l2L2/7q9fN+1lM2sMfo6VNxLMXuvKj4TfR0qZlY9h4qZZdW3oeJdH7PRTfVno29Dxczq4VAxs6z6MlS862M2tqn8jPRlqJhZfRwqZpZV34WKd33MujNVn5W+CxUzq1ff1P1xD8Vs4qaiTpB7KmaWlUPFzLLqi1Dxro9Zb3J+hvoiVMysORwqZpZVVWVPb5G0O9X9eUrSwtK0E6VyqBtHLmtm7VJV2dP/BAYi4qikPwO+BHwqTXs7Ipb22g4za4ZKyp5GxNaIOJoGt1HU98nCB2nN8sj1Waqy7OmwNcATpeFTUznTbZKu67SQy56atUOlV9RK+hNgAPi90uiFEbFf0iJgi6SdEfHqyGUjYj2wHooSHZU02MwmrJKypwCSrgTuBJaXSqASEfvT/V7gaWBZhjaZWU2qKnu6DPgGRaAcLI2fLWlWejwXuAwoH+Adk4+nmOWV4zNVVdnTvwd+BfiOJIAfR8Ry4ELgG5LeoQi4u0ecNTKzlslyTCUiHgceHzHub0qPr+yw3H8AH8rRBjNrBl9Ra2ZZOVTMLKvWhooP0ppNjV4/W60NFTNrJoeKmWXlUDGzrBwqZpaVQ8XMsnKomFlWDhUzy6qVoXL4+X11N8Gsr50748yLJ7tsK0PFzJrLoWJmWTlUzCwrh4qZZeVQMbOsHCpmlpVDxcyyqqrs6SxJj6bpz0o6tzTtjjT+ZUlX52iPmdWn51AplT29FlgC3CBpyYjZ1gBvRMT5wDrgnrTsEopf3/8N4Brgn9LzmVlLVVL2NA0/mB5vAK5Q8bP6K4BHIuJYRLwGDKXnM7OWyvFr+qOVPb200zyppMdbwBlp/LYRy45aMlXSWmAtwIIFC1j9oy9naLqZjeZG/cP2yS7bmgO1EbE+IgYiYmDevHl1N8fMOqiq7Om780iaAbwPONzlsmbWIpWUPU3Dq9Lj64EtERFp/Mp0dug84ALghxnaZGY1qars6X3AtyQNAUcogoc032MU9ZOPA38RESd6bZOZ1UdFh6FdBgYGYnBwsO5mmPUtSdsjYmAyy7bmQK2ZtYNDxcyycqiYWVYOFTPLyqFiZlk5VMwsK4eKmWXlUDGzrBwqZpaVQ8XMsnKomFlWDhUzy8qhYmZZOVTMLCuHipll5VAxs6wcKmaWlUPFzLLqKVQkzZG0WdKedD97lHmWSvqBpBclvSDpU6VpD0h6TdKOdFvaS3vMrH699lRuB56KiAuAp9LwSEeBT0fEcGnTr0g6vTT9sxGxNN129NgeM6tZr6FSLmf6IHDdyBki4pWI2JMe/zdwEHA1MLM+1WuonBkRB9LjnwJnjjWzpEuAmcCrpdFfTLtF6yTNGmPZtZIGJQ0eOnSox2ab2VQZN1QkPSlp1yi39xRhT8XBOtb7kDQf+BZwY0S8k0bfASwGfhuYA9zWaXmXPTVrh3GLiUXElZ2mSXpd0vyIOJBC42CH+X4V+D5wZ0S8W5C91Ms5JumbwK0Tar2ZNU6vuz/lcqargH8dOUMqhfpd4KGI2DBi2vx0L4rjMbt6bI+Z1azXULkb+KikPcCVaRhJA5LuTfN8EvgIsHqUU8fflrQT2AnMBb7QY3vMrGYue2pmJ3HZUzNrDIeKmWXlUDGzrBwqZpaVQ8XMsnKomFlWDhUzy8qhYmZZOVTMLCuHipll5VAxs6wcKmaWlUPFzLJyqJhZVg4VM8vKoWJmWTlUzCwrh4qZZTXlZU/TfCdKv0+7sTT+PEnPShqS9Gj6kWwza7Eqyp4CvF0qbbq8NP4eYF1EnA+8AazpsT1mVrMpL3vaSSrLcTkwXLZjQsubWTNVVfb01FSydJuk4eA4A3gzIo6n4X3AWZ3+kMuemrXDuBUKJT0J/Nook+4sD0RESOpU72NhROyXtAjYkmr9vDWRhkbEemA9FCU6JrKsmVWnkrKnEbE/3e+V9DSwDPgX4HRJM1Jv5Wxg/yTWwcwapIqyp7MlzUqP5wKXAbtTQfetwPVjLW9m7VJF2dMLgUFJz1OEyN0RsTtNuw24RdIQxTGW+3psj5nVzGVPzewkLntqZo3hUDGzrBwqZpaVQ8XMsnKomFlWDhUzy8qhYmZZOVTMLCuHipll5VAxs6wcKmaWlUPFzLJyqJhZVg4VM8vKoWJmWTlUzCwrh4qZZeVQMbOsprzsqaQ/KJU83SHpf4dr/0h6QNJrpWlLe2mPmdVvysueRsTW4ZKnFBUJjwL/Vprls6WSqDt6bI+Z1azqsqfXA09ExNEe/66ZNVRVZU+HrQQeHjHui5JekLRuuD6QmbVXVWVPSRUMPwRsKo2+gyKMZlKUNL0NuKvD8muBtQALFiwYr9lmVpNKyp4mnwS+GxG/KD33cC/nmKRvAreO0Q7XUjZrgSkve1pyAyN2fVIQIUkUx2N29dgeM6tZFWVPkXQucA7w7yOW/7akncBOYC7whR7bY2Y1G3f3ZywRcRi4YpTxg8BNpeH/As4aZb7Le/n7ZtY8vqLWzLJyqJhZVg4VM8vKoWJmWTlUzCwrh4qZZeVQMbOsHCpmlpVDxcyycqiYWVYOFTPLyqFiZlk5VMwsK4eKmWXlUDGzrBwqZpaVQ8XMsnKomFlWDhUzy6rXWsp/JOlFSe9IGhhjvmskvSxpSNLtpfHnSXo2jX9U0sxe2mNm9eu1p7IL+EPgmU4zSDoF+BpwLbAEuEHSkjT5HmBdRJwPvAGs6bE9ZlaznkIlIl6KiJfHme0SYCgi9kbEz4FHgBWp1s/lwIY0Xze1mM2s4Xoq0dGls4CflIb3AZcCZwBvRsTx0viTyngMK5c9paho2I+Fx+YCP6u7EVOkX9etX9fr1ye7YE+1lCNirIqEWZXLnkoajIiOx3Daql/XC/p33fp5vSa7bE+1lLu0n6I64bCz07jDwOmSZqTeyvB4M2uxKk4pPwdckM70zARWAhsjIoCtwPVpvvFqMZtZC/R6SvnjkvYBvwN8X9KmNP4Dkh4HSL2Qm4FNwEvAYxHxYnqK24BbJA1RHGO5r8s/vb6XdjdYv64X9O+6eb1GUNFhMDPLw1fUmllWDhUzy6oVodLr1wGaStIcSZsl7Un3szvMd0LSjnTbWHU7uzXe6y9pVvo6xlD6esa51bdycrpYt9WSDpXep5vqaOdESLpf0sFO13yp8NW0zi9IuqirJ46Ixt+ACykuxnkaGOgwzynAq8AiYCbwPLCk7raPs15fAm5Pj28H7ukw3//U3dYu1mXc1x/4c+Dr6fFK4NG6251x3VYD/1h3Wye4Xh8BLgJ2dZj+MeAJQMCHgWe7ed5W9FSih68DTH3rerKC4usJ0P6vKXTz+pfXdwNwRfq6RtO1cdsaV0Q8AxwZY5YVwENR2EZxXdn88Z63FaHSpdG+DtDxsv+GODMiDqTHPwXO7DDfqZIGJW2T1NTg6eb1f3eeKC41eIviUoKm63bb+kTaTdgg6ZxRprfNpD5TVXz3pytN+TpAbmOtV3kgIkJSp/P7CyNiv6RFwBZJOyPi1dxttZ58D3g4Io5J+lOKHtnlNbepFo0JlZi6rwPUaqz1kvS6pPkRcSB1Kw92eI796X6vpKeBZRT7+E3Szes/PM8+STOA91F8XaPpxl23iCivx70Ux8vablKfqX7a/Rn16wA1t2k8Gym+ngAdvqYgabakWenxXOAyYHdlLexeN69/eX2vB7ZEOiLYcOOu24hjDcsprh5vu43Ap9NZoA8Db5V21zur+wh0l0epP06xP3cMeB3YlMZ/AHh8xNHqVyj+i99Zd7u7WK8zgKeAPcCTwJw0fgC4Nz3+XWAnxRmHncCauts9xvqc9PoDdwHL0+NTge8AQ8APgUV1tznjuv0d8GJ6n7YCi+tucxfr9DBwAPhF+nytAT4DfCZNF8UPrL2atr1Rz7yOvPkyfTPLqp92f8ysARwqZpaVQ8XMsnKomFlWDhUzy8qhYmZZOVTMLKv/A+YXQ8OS8hxTAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] @@ -702,22 +735,22 @@ }, { "cell_type": "code", - "execution_count": 25, + "execution_count": 26, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "" + "" ] }, - "execution_count": 25, + "execution_count": 26, "metadata": {}, "output_type": "execute_result" }, { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAARUAAAD8CAYAAABZ0jAcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAESBJREFUeJzt3V2MXOV9x/Hvr1AbiarBxhZ1ABssUI2rVgZvgRYpbYHwkgtDGpoYqcJERm7a0EpFiQBxgUQSFdILV1HTBgsIIYqAxFVURwG5BkO5iQlr1cHGCPxC0+A62MGAVJk6wfx7cZ6lh/XO7uyeZ87L7O8jjXbmvMw8Z+ac3z7nnDnzV0RgZpbLrzXdADMbLg4VM8vKoWJmWTlUzCwrh4qZZeVQMbOssoSKpIckHZK0q8d4SfqapL2SXpR0UWncGkl70m1NjvaYWXNy9VQeBq6ZZPy1wPnptg74ZwBJ84G7gUuAi4G7Jc3L1CYza0CWUImI54Ajk0xyHfBIFLYBp0laBFwNbImIIxHxFrCFycPJzFru5Jpe50zgZ6XHr6dhvYafQNI6il4Op5566sply5YNpqXWv+0Det6VA3pe69v27dt/ERELZzJvXaFSWURsADYAjIyMxOjoaMMtmgXU0OtOFVa+smTgJP10pvPWdfbnAHB26fFZaViv4VY3TXBrqy61dRaqK1Q2ATels0CXAu9ExEFgM3CVpHnpAO1VaZgN2rBtlMO2PB2WZfdH0qPAHwMLJL1OcUbn1wEi4hvAE8AngL3AUeCzadwRSV8CXkhPdU9ETHbA12ZiNm5kEy2zd5tqkSVUIuLGKcYH8Pke4x4CHsrRDiuZjUEylfJ74oAZGH+j1syy6szZH+uTeyj9GXuf3GPJzqEyDBwkM+ddouwcKl3lIMnPAZOFQ6VLHCT1ccDMmEOlCxwmzfLxl2lxqLSZw6RdHC598SnltnKgtJc/m0m5p9ImXlm7w72WnhwqbeAw6S4f0D2Bd3+a5kAZHv4sAfdUmuMVcDh5t8g9FTPLyz2VurmHMjvM4h6Leyp1cqDMPrPwM3eo1GUWrlyWzLLP3rs/gzbLVijrYRbtDjlUBsVhYhOZBeGSq+zpNZJeSWVN75hg/HpJO9LtVUlvl8YdL43blKM9jXOg2FSGeB2p3FORdBLwdeDjFMXAXpC0KSJ2j00TEX9bmv6vgQtLT/FuRKyo2o7WGOKVxTITQ9ljydFTuRjYGxH7I+KXwGMUZU57uRF4NMPrtotLQ9hMDOF6kyNUplO6dAlwLrC1NPgUSaOStkm6PkN76jdkK4U1YIjWoboP1K4GNkbE8dKwJRFxQNJSYKuknRGxb/yM5VrKixcvrqe1/RiilcEaNiS7Qzl6KtMpXbqacbs+EXEg/d0PPMuHj7eUp9sQESMRMbJw4YzqRufnQLHchmCdyhEqLwDnSzpX0hyK4DjhLI6kZcA84EelYfMkzU33FwCXAbvHz2tm3VF59yci3pN0K0UN5JOAhyLiJUn3AKMRMRYwq4HHUrXCMRcA90t6nyLg7i2fNWq1IfiPYi3V8d0gfXgb74aRkZEYHR1trgEOFKtDg5umpO0RMTKTeX3tz3Q5UKwuHV3X/DX9fnX0A7aO6+DX+t1T6YcDxZrWoXXQoTKVDn2YNuQ6si46VCbTkQ/RZpEOrJMOFTPLyqHSSwf+I9gs1fJ106EykZZ/aGZtXkcdKuO1+MMy+5CWrqsOFTPLyqFS1tLkN+upheusQ8XMsnKojGlh4pv1pWXrrkMFWvehmE1bi9Zhh0qLPgyzSlqyLjtUzCyr2fvTBy1JdbOsWvBTCe6pmFlWszNU3EuxYdfgOl5XLeWbJR0u1Uy+pTRujaQ96bYmR3smb+zAX8GsHRpa12uppZw8HhG3jpt3PnA3MEKxF7g9zftW1XaZWTOaqKVcdjWwJSKOpCDZAlyToU0Tcy/FZpsG1vk6ayl/StKLkjZKGqtoOJ06zOtSzeXRw4cPZ2i2mQ1CXQdqfwCcExG/R9Eb+dZ0n6CVZU/N7AS11FKOiDcj4lh6+ACwst95s/Guj81WNa/7tdRSlrSo9HAV8HK6vxm4KtVUngdclYbl5UCx2a7GbaCuWsp/I2kV8B5wBLg5zXtE0pcoggngnog4UrVNZtac4a+l7F6K2f/rc3N3LWUza43hDhX3Usw+rIZtYrhDxcxq51Axs6yGN1S862M2sQFvG8MbKmbWCIeKmWU1nKHiXR+zyQ1wGxnOUDGzxjhUzCyr4QsV7/qY9WdA28rwhYqZNWp46v64h2I2fQOoE+Seipll5VAxs6yGI1S862NWTcZtaDhCxcxaw6FiZlnVVfb0Nkm7U92fpyUtKY07XiqHumn8vGbWLXWVPf0PYCQijkr6S+CrwGfSuHcjYkXVdphZO9RS9jQinomIo+nhNor6Pnn4IK1ZHpm2pTrLno5ZCzxZenxKKme6TdL1vWZy2VOzbqj1G7WS/hwYAf6oNHhJRByQtBTYKmlnROwbP29EbAA2QFGio5YGm9m01VL2FEDSlcBdwKpSCVQi4kD6ux94FrgwQ5vMrCF1lT29ELifIlAOlYbPkzQ33V8AXAaUD/BOzsdTzPLKsE3VVfb074HfAL4nCeC/ImIVcAFwv6T3KQLu3nFnjcysY7pd9tQ9FbP8wmVPzaxFHCpmllV3Q8W7PmaDUXHb6m6omFkrOVTMLCuHipll5VAxs6wcKmaWlUPFzLJyqJhZVt0Mle1NN8BsuK1k5cqZztvNUDGz1nKomFlWDhUzy8qhYmZZOVTMLCuHipll5VAxs6zqKns6V9Ljafzzks4pjbszDX9F0tU52mNmzakcKqWyp9cCy4EbJS0fN9la4K2IOA9YD9yX5l1O8ev7vwNcA/xTej4z66hayp6mx99K9zcCV6j4Wf3rgMci4lhEvAbsTc9nZh2Vo0LhRGVPL+k1TSrp8Q5wehq+bdy8E5ZMlbQOWAewePFi+GmGlpvZhLZr+4wvhunMgdqI2BARIxExsnDhwqabY2Y91FX29INpJJ0MfAR4s895zaxDail7mh6vSfdvALZGUcVsE7A6nR06Fzgf+HGGNplZQ+oqe/og8G1Je4EjFMFDmu67FPWT3wM+HxHHq7bJzJrT7bKnZjYQLntqZq3hUDGzrBwqZpaVQ8XMsnKomFlWDhUzy8qhYmZZOVTMLCuHipll5VAxs6wcKmaWlUPFzLJyqJhZVg4VM8vKoWJmWTlUzCwrh4qZZeVQMbOsKoWKpPmStkjak/7Om2CaFZJ+JOklSS9K+kxp3MOSXpO0I91WVGmPmTWvak/lDuDpiDgfeDo9Hu8ocFNEjJU2/QdJp5XGfzEiVqTbjortMbOGVQ2VcjnTbwHXj58gIl6NiD3p/n8DhwBXAzMbUlVD5YyIOJju/xw4Y7KJJV0MzAH2lQZ/Je0WrZc0d5J510kalTR6+PDhis02s0GZMlQkPSVp1wS3DxVhT8XBetb7kLQI+Dbw2Yh4Pw2+E1gG/D4wH7i91/wue2rWDVMWE4uIK3uNk/SGpEURcTCFxqEe0/0m8EPgroj4oCB7qZdzTNI3gS9Mq/Vm1jpVd3/K5UzXAP86foJUCvX7wCMRsXHcuEXpryiOx+yq2B4za1jVULkX+LikPcCV6TGSRiQ9kKb5NPAx4OYJTh1/R9JOYCewAPhyxfaYWcNc9tTMTuCyp2bWGg4VM8vKoWJmWTlUzCwrh4qZZeVQMbOsHCpmlpVDxcyycqiYWVYOFTPLyqFiZlk5VMwsK4eKmWXlUDGzrBwqZpaVQ8XMsnKomFlWDhUzy2rgZU/TdMdLv0+7qTT8XEnPS9or6fH0I9lm1mF1lD0FeLdU2nRVafh9wPqIOA94C1hbsT1m1rCBlz3tJZXluBwYK9sxrfnNrJ3qKnt6SipZuk3SWHCcDrwdEe+lx68DZ/Z6IZc9NeuGKSsUSnoK+K0JRt1VfhARIalXvY8lEXFA0lJga6r18850GhoRG4ANUJTomM68ZlafWsqeRsSB9He/pGeBC4F/AU6TdHLqrZwFHJjBMphZi9RR9nSepLnp/gLgMmB3Kuj+DHDDZPObWbfUUfb0AmBU0k8oQuTeiNidxt0O3CZpL8UxlgcrtsfMGuayp2Z2Apc9NbPWcKiYWVYOFTPLyqFiZlk5VMwsK4eKmWXlUDGzrBwqZpaVQ8XMsnKomFlWDhUzy8qhYmZZOVTMLCuHipll5VAxs6wcKmaWlUPFzLJyqJhZVgMveyrpT0olT3dI+t+x2j+SHpb0WmnciirtMbPmDbzsaUQ8M1bylKIi4VHg30qTfLFUEnVHxfaYWcPqLnt6A/BkRByt+Lpm1lJ1lT0dsxp4dNywr0h6UdL6sfpAZtZddZU9JVUw/F1gc2nwnRRhNIeipOntwD095l8HrANYvHjxVM02s4bUUvY0+TTw/Yj4Vem5x3o5xyR9E/jCJO1wLWWzDhh42dOSGxm365OCCEmiOB6zq2J7zKxhdZQ9RdI5wNnAv4+b/zuSdgI7gQXAlyu2x8waNuXuz2Qi4k3gigmGjwK3lB7/J3DmBNNdXuX1zax9/I1aM8vKoWJmWTlUzCwrh4qZZeVQMbOsHCpmlpVDxcyycqiYWVYOFTPLyqFiZlk5VMwsK4eKmWXlUDGzrBwqZpaVQ8XMsnKomFlWDhUzy8qhYmZZOVTMLKuqtZT/TNJLkt6XNDLJdNdIekXSXkl3lIafK+n5NPxxSXOqtMfMmle1p7IL+FPguV4TSDoJ+DpwLbAcuFHS8jT6PmB9RJwHvAWsrdgeM2tYpVCJiJcj4pUpJrsY2BsR+yPil8BjwHWp1s/lwMY0XT+1mM2s5SqV6OjTmcDPSo9fBy4BTgfejoj3SsNPKOMxplz2lKKi4TAWHlsA/KLpRgzIsC7bsC7Xb890xkq1lCNisoqEWZXLnkoajYiex3C6aliXC4Z32YZ5uWY6b6Vayn06QFGdcMxZadibwGmSTk69lbHhZtZhdZxSfgE4P53pmQOsBjZFRADPADek6aaqxWxmHVD1lPInJb0O/AHwQ0mb0/CPSnoCIPVCbgU2Ay8D342Il9JT3A7cJmkvxTGWB/t86Q1V2t1iw7pcMLzL5uUaR0WHwcwsD3+j1syycqiYWVadCJWqlwO0laT5krZI2pP+zusx3XFJO9JtU93t7NdU77+kuelyjL3p8oxz6m/lzPSxbDdLOlz6nG5pop3TIekhSYd6fedLha+lZX5R0kV9PXFEtP4GXEDxZZxngZEe05wE7AOWAnOAnwDLm277FMv1VeCOdP8O4L4e0/1P023tY1mmfP+BvwK+ke6vBh5vut0Zl+1m4B+bbus0l+tjwEXArh7jPwE8CQi4FHi+n+ftRE8lKlwOMPjWVXIdxeUJ0P3LFPp5/8vLuxG4Il2u0XZdXLemFBHPAUcmmeQ64JEobKP4XtmiqZ63E6HSp4kuB+j5tf+WOCMiDqb7PwfO6DHdKZJGJW2T1Nbg6ef9/2CaKL5q8A7FVwnart9161NpN2GjpLMnGN81M9qm6rj2py9tuRwgt8mWq/wgIkJSr/P7SyLigKSlwFZJOyNiX+62WiU/AB6NiGOS/oKiR3Z5w21qRGtCJQZ3OUCjJlsuSW9IWhQRB1O38lCP5ziQ/u6X9CxwIcU+fpv08/6PTfO6pJOBj1BcrtF2Uy5bRJSX4wGK42VdN6Ntaph2fya8HKDhNk1lE8XlCdDjMgVJ8yTNTfcXAJcBu2trYf/6ef/Ly3sDsDXSEcGWm3LZxh1rWEXx7fGu2wTclM4CXQq8U9pd763pI9B9HqX+JMX+3DHgDWBzGv5R4IlxR6tfpfgvflfT7e5juU4Hngb2AE8B89PwEeCBdP8PgZ0UZxx2Amubbvcky3PC+w/cA6xK908BvgfsBX4MLG26zRmX7e+Al9Ln9AywrOk297FMjwIHgV+l7Wst8Dngc2m8KH5gbV9a9yY88zr+5q/pm1lWw7T7Y2Yt4FAxs6wcKmaWlUPFzLJyqJhZVg4VM8vKoWJmWf0f1/avloPSXUoAAAAASUVORK5CYII=\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAARUAAAD8CAYAAABZ0jAcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAARIElEQVR4nO3dXYxc5X3H8e+vUBuJqsHGFnUAGyxQjatWBm+BFiltgfCSC0MamhipwkRGbtrQSkWJAHGBRBIV0gtXUdMGCwghioDEVVRHAbkGQ7mJCWvVwcYI/ELT4DrYwYBUmTrB/HtxnqWH9c7u7J5nzsvs7yONdua8zDxn5pzfPuecOfNXRGBmlsuvNd0AMxsuDhUzy8qhYmZZOVTMLCuHipll5VAxs6yyhIqkhyQdkrSrx3hJ+pqkvZJelHRRadwaSXvSbU2O9phZc3L1VB4Grplk/LXA+em2DvhnAEnzgbuBS4CLgbslzcvUJjNrQJZQiYjngCOTTHId8EgUtgGnSVoEXA1siYgjEfEWsIXJw8nMWu7kml7nTOBnpcevp2G9hp9A0jqKXg6nnnrqymXLlg2mpda/7QN63pUDel7r2/bt238REQtnMm9doVJZRGwANgCMjIzE6Ohowy2aBdTQ604VVr6yZOAk/XSm89Z19ucAcHbp8VlpWK/hVjdNcGurLrV1FqorVDYBN6WzQJcC70TEQWAzcJWkeekA7VVpmA3asG2Uw7Y8HZZl90fSo8AfAwskvU5xRufXASLiG8ATwCeAvcBR4LNp3BFJXwJeSE91T0RMdsDXZmI2bmQTLbN3m2qRJVQi4sYpxgfw+R7jHgIeytEOK5mNQTKV8nvigBkYf6PWzLLqzNkf65N7KP0Ze5/cY8nOoTIMHCQz512i7BwqXeUgyc8Bk4VDpUscJPVxwMyYQ6ULHCbN8vGXaXGotJnDpF0cLn3xKeW2cqC0lz+bSbmn0iZeWbvDvZaeHCpt4DDpLh/QPYF3f5rmQBke/iwB91Sa4xVwOHm3yD0VM8vLPZW6uYcyO8ziHot7KnVyoMw+s/Azd6jUZRauXJbMss/euz+DNstWKOthFu0OOVQGxWFiE5kF4ZKr7Ok1kl5JZU3vmGD8ekk70u1VSW+Xxh0vjduUoz2Nc6DYVIZ4HancU5F0EvB14OMUxcBekLQpInaPTRMRf1ua/q+BC0tP8W5ErKjajtYY4pXFMhND2WPJ0VO5GNgbEfsj4pfAYxRlTnu5EXg0w+u2i0tD2EwM4XqTI1SmU7p0CXAusLU0+BRJo5K2Sbo+Q3vqN2QrhTVgiNahug/UrgY2RsTx0rAlEXFA0lJgq6SdEbFv/IzlWsqLFy+up7X9GKKVwRo2JLtDOXoq0ylduppxuz4RcSD93Q88y4ePt5Sn2xARIxExsnDhjOpG5+dAsdyGYJ3KESovAOdLOlfSHIrgOOEsjqRlwDzgR6Vh8yTNTfcXAJcBu8fPa2bdUXn3JyLek3QrRQ3kk4CHIuIlSfcAoxExFjCrgcdStcIxFwD3S3qfIuDuLZ81arUh+I9iLdXx3SB9eBvvhpGRkRgdHW2uAQ4Uq0ODm6ak7RExMpN5fe3PdDlQrC4dXdf8Nf1+dfQDto7r4Nf63VPphwPFmtahddChMpUOfZg25DqyLjpUJtORD9FmkQ6skw4VM8vKodJLB/4j2CzV8nXToTKRln9oZm1eRx0q47X4wzL7kJauqw4VM8vKoVLW0uQ366mF66xDxcyycqiMaWHim/WlZeuuQwVa96GYTVuL1mGHSos+DLNKWrIuO1TMLKvZ+9MHLUl1s6xa8FMJ7qmYWVazM1TcS7Fh1+A6Xlct5ZslHS7VTL6lNG6NpD3ptiZHeyZv7MBfwawdGlrXa6mlnDweEbeOm3c+cDcwQrEXuD3N+1bVdplZM5qopVx2NbAlIo6kINkCXJOhTRNzL8VmmwbW+TprKX9K0ouSNkoaq2g4nTrM61LN5dHDhw9naLaZDUJdB2p/AJwTEb9H0Rv51nSfoJVlT83sBLXUUo6INyPiWHr4ALCy33mz8a6PzVY1r/u11FKWtKj0cBXwcrq/Gbgq1VSeB1yVhuXlQLHZrsZtoK5ayn8jaRXwHnAEuDnNe0TSlyiCCeCeiDhStU1m1pzhr6XsXorZ/+tzc3ctZTNrjeEOFfdSzD6shm1iuEPFzGrnUDGzrIY3VLzrYzaxAW8bwxsqZtYIh4qZZTWcoeJdH7PJDXAbGc5QMbPGOFTMLKvhCxXv+pj1Z0DbyvCFipk1anjq/riHYjZ9A6gT5J6KmWXlUDGzrIYjVLzrY1ZNxm1oOELFzFrDoWJmWdVV9vQ2SbtT3Z+nJS0pjTteKoe6afy8ZtYtdZU9/Q9gJCKOSvpL4KvAZ9K4dyNiRdV2mFk71FL2NCKeiYij6eE2ivo+efggrVkembalOsuejlkLPFl6fEoqZ7pN0vW9ZnLZU7NuqPUbtZL+HBgB/qg0eElEHJC0FNgqaWdE7Bs/b0RsADZAUaKjlgab2bTVUvYUQNKVwF3AqlIJVCLiQPq7H3gWuDBDm8ysIXWVPb0QuJ8iUA6Vhs+TNDfdXwBcBpQP8E7Ox1PM8sqwTdVV9vTvgd8AvicJ4L8iYhVwAXC/pPcpAu7ecWeNzKxjul321D0Vs/zCZU/NrEUcKmaWVXdDxbs+ZoNRcdvqbqiYWSs5VMwsK4eKmWXlUDGzrBwqZpaVQ8XMsnKomFlW3QyV7U03wGy4rWTlypnO281QMbPWcqiYWVYOFTPLyqFiZlk5VMwsK4eKmWXlUDGzrOoqezpX0uNp/POSzimNuzMNf0XS1TnaY2bNqRwqpbKn1wLLgRslLR832VrgrYg4D1gP3JfmXU7x6/u/A1wD/FN6PjPrqFrKnqbH30r3NwJXqPhZ/euAxyLiWES8BuxNz2dmHZWjQuFEZU8v6TVNKunxDnB6Gr5t3LwTlkyVtA5YB7B48WL4aYaWm9mEtmv7jC+G6cyB2ojYEBEjETGycOHCpptjZj3UVfb0g2kknQx8BHizz3nNrENqKXuaHq9J928AtkZRxWwTsDqdHToXOB/4cYY2mVlD6ip7+iDwbUl7gSMUwUOa7rsU9ZPfAz4fEcertsnMmtPtsqdmNhAue2pmreFQMbOsHCpmlpVDxcyycqiYWVYOFTPLyqFiZlk5VMwsK4eKmWXlUDGzrBwqZpaVQ8XMsnKomFlWDhUzy8qhYmZZOVTMLCuHipll5VAxs6wqhYqk+ZK2SNqT/s6bYJoVkn4k6SVJL0r6TGncw5Jek7Qj3VZUaY+ZNa9qT+UO4OmIOB94Oj0e7yhwU0SMlTb9B0mnlcZ/MSJWpNuOiu0xs4ZVDZVyOdNvAdePnyAiXo2IPen+fwOHAFcDMxtSVUPljIg4mO7/HDhjsoklXQzMAfaVBn8l7RatlzR3knnXSRqVNHr48OGKzTazQZkyVCQ9JWnXBLcPFWFPxcF61vuQtAj4NvDZiHg/Db4TWAb8PjAfuL3X/C57atYNUxYTi4gre42T9IakRRFxMIXGoR7T/SbwQ+CuiPigIHupl3NM0jeBL0yr9WbWOlV3f8rlTNcA/zp+glQK9fvAIxGxcdy4RemvKI7H7KrYHjNrWNVQuRf4uKQ9wJXpMZJGJD2Qpvk08DHg5glOHX9H0k5gJ7AA+HLF9phZw1z21MxO4LKnZtYaDhUzy8qhYmZZOVTMLCuHipll5VAxs6wcKmaWlUPFzLJyqJhZVg4VM8vKoWJmWTlUzCwrh4qZZeVQMbOsHCpmlpVDxcyycqiYWVYOFTPLauBlT9N0x0u/T7upNPxcSc9L2ivp8fQj2WbWYXWUPQV4t1TadFVp+H3A+og4D3gLWFuxPWbWsIGXPe0lleW4HBgr2zGt+c2sneoqe3pKKlm6TdJYcJwOvB0R76XHrwNn9nohlz0164YpKxRKegr4rQlG3VV+EBEhqVe9jyURcUDSUmBrqvXzznQaGhEbgA1QlOiYzrxmVp9ayp5GxIH0d7+kZ4ELgX8BTpN0cuqtnAUcmMEymFmL1FH2dJ6kuen+AuAyYHcq6P4McMNk85tZt9RR9vQCYFTSTyhC5N6I2J3G3Q7cJmkvxTGWByu2x8wa5rKnZnYClz01s9ZwqJhZVg4VM8vKoWJmWTlUzCwrh4qZZeVQMbOsHCpmlpVDxcyycqiYWVYOFTPLyqFiZlk5VMwsK4eKmWXlUDGzrBwqZpaVQ8XMsnKomFlWAy97KulPSiVPd0j637HaP5IelvRaadyKKu0xs+YNvOxpRDwzVvKUoiLhUeDfSpN8sVQSdUfF9phZw+oue3oD8GREHK34umbWUnWVPR2zGnh03LCvSHpR0vqx+kBm1l11lT0lVTD8XWBzafCdFGE0h6Kk6e3APT3mXwesA1i8ePFUzTazhtRS9jT5NPD9iPhV6bnHejnHJH0T+MIk7XAtZbMOGHjZ05IbGbfrk4IISaI4HrOrYnvMrGF1lD1F0jnA2cC/j5v/O5J2AjuBBcCXK7bHzBo25e7PZCLiTeCKCYaPAreUHv8ncOYE011e5fXNrH38jVozy8qhYmZZOVTMLCuHipll5VAxs6wcKmaWlUPFzLJyqJhZVg4VM8vKoWJmWTlUzCwrh4qZZeVQMbOsHCpmlpVDxcyycqiYWVYOFTPLyqFiZlk5VMwsq6q1lP9M0kuS3pc0Msl010h6RdJeSXeUhp8r6fk0/HFJc6q0x8yaV7Wnsgv4U+C5XhNIOgn4OnAtsBy4UdLyNPo+YH1EnAe8Bayt2B4za1ilUImIlyPilSkmuxjYGxH7I+KXwGPAdanWz+XAxjRdP7WYzazlKpXo6NOZwM9Kj18HLgFOB96OiPdKw08o4zGmXPaUoqLhMBYeWwD8oulGDMiwLtuwLtdvz3TGSrWUI2KyioRZlcueShqNiJ7HcLpqWJcLhnfZhnm5ZjpvpVrKfTpAUZ1wzFlp2JvAaZJOTr2VseFm1mF1nFJ+ATg/nemZA6wGNkVEAM8AN6TppqrFbGYdUPWU8iclvQ78AfBDSZvT8I9KegIg9UJuBTYDLwPfjYiX0lPcDtwmaS/FMZYH+3zpDVXa3WLDulwwvMvm5RpHRYfBzCwPf6PWzLJyqJhZVp0IlaqXA7SVpPmStkjak/7O6zHdcUk70m1T3e3s11Tvv6S56XKMvenyjHPqb+XM9LFsN0s6XPqcbmmindMh6SFJh3p950uFr6VlflHSRX09cUS0/gZcQPFlnGeBkR7TnATsA5YCc4CfAMubbvsUy/VV4I50/w7gvh7T/U/Tbe1jWaZ8/4G/Ar6R7q8GHm+63RmX7WbgH5tu6zSX62PARcCuHuM/ATwJCLgUeL6f5+1ETyUqXA4w+NZVch3F5QnQ/csU+nn/y8u7EbgiXa7Rdl1ct6YUEc8BRyaZ5DrgkShso/he2aKpnrcTodKniS4H6Pm1/5Y4IyIOpvs/B87oMd0pkkYlbZPU1uDp5/3/YJoovmrwDsVXCdqu33XrU2k3YaOksycY3zUz2qbquPanL225HCC3yZar/CAiQlKv8/tLIuKApKXAVkk7I2Jf7rZaJT8AHo2IY5L+gqJHdnnDbWpEa0IlBnc5QKMmWy5Jb0haFBEHU7fyUI/nOJD+7pf0LHAhxT5+m/Tz/o9N87qkk4GPUFyu0XZTLltElJfjAYrjZV03o21qmHZ/JrwcoOE2TWUTxeUJ0OMyBUnzJM1N9xcAlwG7a2th//p5/8vLewOwNdIRwZabctnGHWtYRfHt8a7bBNyUzgJdCrxT2l3vrekj0H0epf4kxf7cMeANYHMa/lHgiXFHq1+l+C9+V9Pt7mO5TgeeBvYATwHz0/AR4IF0/w+BnRRnHHYCa5tu9yTLc8L7D9wDrEr3TwG+B+wFfgwsbbrNGZft74CX0uf0DLCs6Tb3sUyPAgeBX6Xtay3wOeBzabwofmBtX1r3JjzzOv7mr+mbWVbDtPtjZi3gUDGzrBwqZpaVQ8XMsnKomFlWDhUzy8qhYmZZ/R/X9q+Wg9JdSgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] @@ -750,7 +783,7 @@ }, { "cell_type": "code", - "execution_count": 26, + "execution_count": 27, "metadata": {}, "outputs": [], "source": [ @@ -768,7 +801,7 @@ }, { "cell_type": "code", - "execution_count": 27, + "execution_count": 28, "metadata": {}, "outputs": [], "source": [ @@ -786,16 +819,16 @@ }, { "cell_type": "code", - "execution_count": 28, + "execution_count": 29, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ - "/home/romano/openmc/openmc/mixin.py:71: IDWarning: Another Cell instance already exists with id=1.\n", + "/home/sam/openmc/openmc/openmc/mixin.py:71: IDWarning: Another Cell instance already exists with id=1.\n", " warn(msg, IDWarning)\n", - "/home/romano/openmc/openmc/mixin.py:71: IDWarning: Another Cell instance already exists with id=2.\n", + "/home/sam/openmc/openmc/openmc/mixin.py:71: IDWarning: Another Cell instance already exists with id=2.\n", " warn(msg, IDWarning)\n" ] } @@ -822,7 +855,7 @@ }, { "cell_type": "code", - "execution_count": 29, + "execution_count": 30, "metadata": {}, "outputs": [], "source": [ @@ -842,7 +875,7 @@ }, { "cell_type": "code", - "execution_count": 30, + "execution_count": 31, "metadata": {}, "outputs": [], "source": [ @@ -862,7 +895,7 @@ }, { "cell_type": "code", - "execution_count": 31, + "execution_count": 32, "metadata": {}, "outputs": [ { @@ -871,7 +904,7 @@ "openmc.region.Intersection" ] }, - "execution_count": 31, + "execution_count": 32, "metadata": {}, "output_type": "execute_result" } @@ -891,7 +924,7 @@ }, { "cell_type": "code", - "execution_count": 32, + "execution_count": 33, "metadata": {}, "outputs": [], "source": [ @@ -907,7 +940,7 @@ }, { "cell_type": "code", - "execution_count": 33, + "execution_count": 34, "metadata": {}, "outputs": [ { @@ -954,7 +987,7 @@ }, { "cell_type": "code", - "execution_count": 34, + "execution_count": 35, "metadata": {}, "outputs": [], "source": [ @@ -971,7 +1004,7 @@ }, { "cell_type": "code", - "execution_count": 35, + "execution_count": 36, "metadata": {}, "outputs": [], "source": [ @@ -984,7 +1017,7 @@ }, { "cell_type": "code", - "execution_count": 36, + "execution_count": 37, "metadata": {}, "outputs": [ { @@ -1026,7 +1059,7 @@ }, { "cell_type": "code", - "execution_count": 37, + "execution_count": 38, "metadata": {}, "outputs": [], "source": [ @@ -1045,7 +1078,7 @@ }, { "cell_type": "code", - "execution_count": 38, + "execution_count": 39, "metadata": {}, "outputs": [], "source": [ @@ -1062,7 +1095,7 @@ }, { "cell_type": "code", - "execution_count": 39, + "execution_count": 40, "metadata": {}, "outputs": [ { @@ -1100,7 +1133,7 @@ }, { "cell_type": "code", - "execution_count": 40, + "execution_count": 41, "metadata": { "scrolled": true }, @@ -1134,30 +1167,34 @@ " %%%%%%%%%%%\n", "\n", " | The OpenMC Monte Carlo Code\n", - " Copyright | 2011-2019 MIT and OpenMC contributors\n", + " Copyright | 2011-2020 MIT and OpenMC contributors\n", " License | http://openmc.readthedocs.io/en/latest/license.html\n", - " Version | 0.11.0-dev\n", - " Git SHA1 | 61c911cffdae2406f9f4bc667a9a6954748bb70c\n", - " Date/Time | 2019-07-19 06:20:10\n", - " OpenMP Threads | 4\n", + " Version | 0.12.0-dev\n", + " Git SHA1 | dd74b2f43f1d2060486de2823ee30058b8971dbd\n", + " Date/Time | 2020-03-03 13:58:53\n", + " MPI Processes | 1\n", + " OpenMP Threads | 8\n", "\n", " Reading settings XML file...\n", " Reading cross sections XML file...\n", " Reading materials XML file...\n", " Reading geometry XML file...\n", - " Reading U235 from /opt/data/hdf5/nndc_hdf5_v15/U235.h5\n", - " Reading U238 from /opt/data/hdf5/nndc_hdf5_v15/U238.h5\n", - " Reading O16 from /opt/data/hdf5/nndc_hdf5_v15/O16.h5\n", - " Reading Zr90 from /opt/data/hdf5/nndc_hdf5_v15/Zr90.h5\n", - " Reading Zr91 from /opt/data/hdf5/nndc_hdf5_v15/Zr91.h5\n", - " Reading Zr92 from /opt/data/hdf5/nndc_hdf5_v15/Zr92.h5\n", - " Reading Zr94 from /opt/data/hdf5/nndc_hdf5_v15/Zr94.h5\n", - " Reading Zr96 from /opt/data/hdf5/nndc_hdf5_v15/Zr96.h5\n", - " Reading H1 from /opt/data/hdf5/nndc_hdf5_v15/H1.h5\n", - " Reading O17 from /opt/data/hdf5/nndc_hdf5_v15/O17.h5\n", - " Reading c_H_in_H2O from /opt/data/hdf5/nndc_hdf5_v15/c_H_in_H2O.h5\n", + " Reading U235 from /home/sam/openmc/libs/nndc_hdf5/U235.h5\n", + " Reading U238 from /home/sam/openmc/libs/nndc_hdf5/U238.h5\n", + " Reading O16 from /home/sam/openmc/libs/nndc_hdf5/O16.h5\n", + " Reading Zr90 from /home/sam/openmc/libs/nndc_hdf5/Zr90.h5\n", + " Reading Zr91 from /home/sam/openmc/libs/nndc_hdf5/Zr91.h5\n", + " Reading Zr92 from /home/sam/openmc/libs/nndc_hdf5/Zr92.h5\n", + " Reading Zr94 from /home/sam/openmc/libs/nndc_hdf5/Zr94.h5\n", + " Reading Zr96 from /home/sam/openmc/libs/nndc_hdf5/Zr96.h5\n", + " Reading H1 from /home/sam/openmc/libs/nndc_hdf5/H1.h5\n", + " Reading O17 from /home/sam/openmc/libs/nndc_hdf5/O17.h5\n", + " Reading c_H_in_H2O from /home/sam/openmc/libs/nndc_hdf5/c_H_in_H2O.h5\n", " Maximum neutron transport energy: 20000000.000000 eV for U235\n", + " Minimum neutron data temperature: 294.000000 K\n", + " Maximum neutron data temperature: 294.000000 K\n", " Reading tallies XML file...\n", + " Preparing distributed cell instances...\n", " Writing summary.h5 file...\n", " Initializing source particles...\n", "\n", @@ -1269,20 +1306,20 @@ "\n", " =======================> TIMING STATISTICS <=======================\n", "\n", - " Total time for initialization = 7.5853e-01 seconds\n", - " Reading cross sections = 7.3383e-01 seconds\n", - " Total time in simulation = 6.5719e+00 seconds\n", - " Time in transport only = 5.9772e+00 seconds\n", - " Time in inactive batches = 4.8850e-01 seconds\n", - " Time in active batches = 6.0834e+00 seconds\n", - " Time synchronizing fission bank = 5.9939e-03 seconds\n", - " Sampling source sites = 5.1295e-03 seconds\n", - " SEND/RECV source sites = 7.3640e-04 seconds\n", - " Time accumulating tallies = 9.9301e-05 seconds\n", - " Total time for finalization = 1.1585e-04 seconds\n", - " Total time elapsed = 7.3346e+00 seconds\n", - " Calculation Rate (inactive) = 20471.0 particles/second\n", - " Calculation Rate (active) = 14794.4 particles/second\n", + " Total time for initialization = 9.7663e-01 seconds\n", + " Reading cross sections = 8.7037e-01 seconds\n", + " Total time in simulation = 1.8046e+00 seconds\n", + " Time in transport only = 1.7523e+00 seconds\n", + " Time in inactive batches = 2.0849e-01 seconds\n", + " Time in active batches = 1.5961e+00 seconds\n", + " Time synchronizing fission bank = 6.7422e-03 seconds\n", + " Sampling source sites = 4.9715e-03 seconds\n", + " SEND/RECV source sites = 1.1323e-03 seconds\n", + " Time accumulating tallies = 5.3810e-04 seconds\n", + " Total time for finalization = 5.1810e-04 seconds\n", + " Total time elapsed = 2.7828e+00 seconds\n", + " Calculation Rate (inactive) = 47962.8 particles/second\n", + " Calculation Rate (active) = 56386.7 particles/second\n", "\n", " ============================> RESULTS <============================\n", "\n", @@ -1308,7 +1345,7 @@ }, { "cell_type": "code", - "execution_count": 41, + "execution_count": 42, "metadata": {}, "outputs": [ { @@ -1341,7 +1378,7 @@ }, { "cell_type": "code", - "execution_count": 42, + "execution_count": 43, "metadata": {}, "outputs": [], "source": [ @@ -1362,7 +1399,7 @@ }, { "cell_type": "code", - "execution_count": 43, + "execution_count": 44, "metadata": {}, "outputs": [ { @@ -1397,7 +1434,7 @@ }, { "cell_type": "code", - "execution_count": 44, + "execution_count": 45, "metadata": {}, "outputs": [ { @@ -1429,18 +1466,20 @@ " %%%%%%%%%%%\n", "\n", " | The OpenMC Monte Carlo Code\n", - " Copyright | 2011-2019 MIT and OpenMC contributors\n", + " Copyright | 2011-2020 MIT and OpenMC contributors\n", " License | http://openmc.readthedocs.io/en/latest/license.html\n", - " Version | 0.11.0-dev\n", - " Git SHA1 | 61c911cffdae2406f9f4bc667a9a6954748bb70c\n", - " Date/Time | 2019-07-19 06:20:18\n", - " OpenMP Threads | 4\n", + " Version | 0.12.0-dev\n", + " Git SHA1 | dd74b2f43f1d2060486de2823ee30058b8971dbd\n", + " Date/Time | 2020-03-03 13:58:56\n", + " MPI Processes | 1\n", + " OpenMP Threads | 8\n", "\n", " Reading settings XML file...\n", " Reading cross sections XML file...\n", " Reading materials XML file...\n", " Reading geometry XML file...\n", " Reading tallies XML file...\n", + " Preparing distributed cell instances...\n", " Reading plot XML file...\n", "\n", " =======================> PLOTTING SUMMARY <========================\n", @@ -1449,11 +1488,11 @@ "Plot file: pinplot.ppm\n", "Universe depth: -1\n", "Plot Type: Slice\n", - "Origin: 0 0 0\n", + "Origin: 0.0 0.0 0.0\n", "Width: 1.26 1.26\n", "Coloring: Materials\n", "Basis: XY\n", - "Pixels: 200 200 \n", + "Pixels: 200 200\n", "\n", " Processing plot 1: pinplot.ppm...\n" ] @@ -1472,7 +1511,7 @@ }, { "cell_type": "code", - "execution_count": 45, + "execution_count": 46, "metadata": {}, "outputs": [], "source": [ @@ -1488,19 +1527,19 @@ }, { "cell_type": "code", - "execution_count": 46, + "execution_count": 47, "metadata": { "scrolled": false }, "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAMgAAADIAgMAAADQNkYNAAAABGdBTUEAALGPC/xhBQAAACBjSFJNAAB6JgAAgIQAAPoAAACA6AAAdTAAAOpgAAA6mAAAF3CculE8AAAADFBMVEUAAP9yEhL//////wAZPRNOAAAAAWJLR0QCZgt8ZAAAAAd0SU1FB+MHEwEUEnBdK8cAAAI7SURBVGje7ZmxkcIwEEUhcAnuhxIIEIGvAwioxiU4QASU4GougA4gsM6SfNxh766sfx7N3Ix+zJv9X7KxtLtaZWVlZWWl1Fr12sQQpXLaRRMRTKFe2kbk+Na8POpNkbbmWlMqtkwxRsJllIotMxQ56F7neWX8LmrdGtNpPWc//Z4cLsbpfp6xN85XNRA904Sd2V/stXlJ16EFcL4O7Q/SnUPOrK/9xfzSvQ44K0dFhjLCmq0nRYYyvDPrq2rfka4RnVlfoyJ9GdGZ9XUdI89aWGYbpTITNUKYgvLlnXFhSsqXd7bjEcKXc8YhvYEjhdzY/Db9lUKebP4+/b6lkK7m8nNRhDAlE8WFoREuig/DpKej+DAbGmGiuDAUUsgItWR8ejZ/yaZ3+SmET+/zL4JIC8Ys2VpI7/JPkSKEbCnkyiNPCilDyI5A+AVzSzZFVAhRBCKssVvlBRB5W8iNwZCThDymiLyT5F5iSCsh3RSRN5/cfgwxoiikkpFmgqgwopZAPmTkcxnkKCO3BZDQs088/RhykpFHRjLyP5AE70uadz/Z/1iK/2TgYwF9klJ8K4GPeJrTBXDsSXMeS3O2RA69wNEaOMAD1wTgMgJceYCLVZobH3IVBS68wLUauLwDLQKgEYG0O4CmCtC6ARpEQBsKaXYBLTWgcQe0B4EmJNLqBBqqQNsWaA4jLWig0Q2005GmPTAaAAYQyJgDGKYgIxtgMISMn4AhFzJKAwZ2yFgQGT4iI86srKysrD/rC4LWcCSWwIp+AAAAJXRFWHRkYXRlOmNyZWF0ZQAyMDE5LTA3LTE5VDA2OjIwOjE4LTA1OjAwIrpoEwAAACV0RVh0ZGF0ZTptb2RpZnkAMjAxOS0wNy0xOVQwNjoyMDoxOC0wNTowMFPn0K8AAAAASUVORK5CYII=\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAMgAAADIAgMAAADQNkYNAAAABGdBTUEAALGPC/xhBQAAACBjSFJNAAB6JgAAgIQAAPoAAACA6AAAdTAAAOpgAAA6mAAAF3CculE8AAAADFBMVEUAAP9yEhL//////wAZPRNOAAAAAWJLR0QCZgt8ZAAAAAd0SU1FB+QDAw06OBFn074AAAI7SURBVGje7ZmxkcIwEEUhcAnuhxIIEIGvAwioxiU4QASU4GougA4gsM6SfNxh766sfx7N3Ix+zJv9X7KxtLtaZWVlZWWl1Fr12sQQpXLaRRMRTKFe2kbk+Na8POpNkbbmWlMqtkwxRsJllIotMxQ56F7neWX8LmrdGtNpPWc//Z4cLsbpfp6xN85XNRA904Sd2V/stXlJ16EFcL4O7Q/SnUPOrK/9xfzSvQ44K0dFhjLCmq0nRYYyvDPrq2rfka4RnVlfoyJ9GdGZ9XUdI89aWGYbpTITNUKYgvLlnXFhSsqXd7bjEcKXc8YhvYEjhdzY/Db9lUKebP4+/b6lkK7m8nNRhDAlE8WFoREuig/DpKej+DAbGmGiuDAUUsgItWR8ejZ/yaZ3+SmET+/zL4JIC8Ys2VpI7/JPkSKEbCnkyiNPCilDyI5A+AVzSzZFVAhRBCKssVvlBRB5W8iNwZCThDymiLyT5F5iSCsh3RSRN5/cfgwxoiikkpFmgqgwopZAPmTkcxnkKCO3BZDQs088/RhykpFHRjLyP5AE70uadz/Z/1iK/2TgYwF9klJ8K4GPeJrTBXDsSXMeS3O2RA69wNEaOMAD1wTgMgJceYCLVZobH3IVBS68wLUauLwDLQKgEYG0O4CmCtC6ARpEQBsKaXYBLTWgcQe0B4EmJNLqBBqqQNsWaA4jLWig0Q2005GmPTAaAAYQyJgDGKYgIxtgMISMn4AhFzJKAwZ2yFgQGT4iI86srKysrD/rC4LWcCSWwIp+AAAAJXRFWHRkYXRlOmNyZWF0ZQAyMDIwLTAzLTAzVDEzOjU4OjU2KzAwOjAw73cO+QAAACV0RVh0ZGF0ZTptb2RpZnkAMjAyMC0wMy0wM1QxMzo1ODo1NiswMDowMJ4qtkUAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, - "execution_count": 46, + "execution_count": 47, "metadata": {}, "output_type": "execute_result" } @@ -1519,17 +1558,17 @@ }, { "cell_type": "code", - "execution_count": 47, + "execution_count": 48, "metadata": {}, "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAMgAAADIAgMAAADQNkYNAAAABGdBTUEAALGPC/xhBQAAACBjSFJNAAB6JgAAgIQAAPoAAACA6AAAdTAAAOpgAAA6mAAAF3CculE8AAAADFBMVEUAAP9yEhL//////wAZPRNOAAAAAWJLR0QCZgt8ZAAAAAd0SU1FB+MHEwEUEnBdK8cAAAI7SURBVGje7ZmxkcIwEEUhcAnuhxIIEIGvAwioxiU4QASU4GougA4gsM6SfNxh766sfx7N3Ix+zJv9X7KxtLtaZWVlZWWl1Fr12sQQpXLaRRMRTKFe2kbk+Na8POpNkbbmWlMqtkwxRsJllIotMxQ56F7neWX8LmrdGtNpPWc//Z4cLsbpfp6xN85XNRA904Sd2V/stXlJ16EFcL4O7Q/SnUPOrK/9xfzSvQ44K0dFhjLCmq0nRYYyvDPrq2rfka4RnVlfoyJ9GdGZ9XUdI89aWGYbpTITNUKYgvLlnXFhSsqXd7bjEcKXc8YhvYEjhdzY/Db9lUKebP4+/b6lkK7m8nNRhDAlE8WFoREuig/DpKej+DAbGmGiuDAUUsgItWR8ejZ/yaZ3+SmET+/zL4JIC8Ys2VpI7/JPkSKEbCnkyiNPCilDyI5A+AVzSzZFVAhRBCKssVvlBRB5W8iNwZCThDymiLyT5F5iSCsh3RSRN5/cfgwxoiikkpFmgqgwopZAPmTkcxnkKCO3BZDQs088/RhykpFHRjLyP5AE70uadz/Z/1iK/2TgYwF9klJ8K4GPeJrTBXDsSXMeS3O2RA69wNEaOMAD1wTgMgJceYCLVZobH3IVBS68wLUauLwDLQKgEYG0O4CmCtC6ARpEQBsKaXYBLTWgcQe0B4EmJNLqBBqqQNsWaA4jLWig0Q2005GmPTAaAAYQyJgDGKYgIxtgMISMn4AhFzJKAwZ2yFgQGT4iI86srKysrD/rC4LWcCSWwIp+AAAAJXRFWHRkYXRlOmNyZWF0ZQAyMDE5LTA3LTE5VDA2OjIwOjE4LTA1OjAwIrpoEwAAACV0RVh0ZGF0ZTptb2RpZnkAMjAxOS0wNy0xOVQwNjoyMDoxOC0wNTowMFPn0K8AAAAASUVORK5CYII=\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAMgAAADIAgMAAADQNkYNAAAABGdBTUEAALGPC/xhBQAAACBjSFJNAAB6JgAAgIQAAPoAAACA6AAAdTAAAOpgAAA6mAAAF3CculE8AAAADFBMVEUAAP9yEhL//////wAZPRNOAAAAAWJLR0QCZgt8ZAAAAAd0SU1FB+QDAw06OBFn074AAAI7SURBVGje7ZmxkcIwEEUhcAnuhxIIEIGvAwioxiU4QASU4GougA4gsM6SfNxh766sfx7N3Ix+zJv9X7KxtLtaZWVlZWWl1Fr12sQQpXLaRRMRTKFe2kbk+Na8POpNkbbmWlMqtkwxRsJllIotMxQ56F7neWX8LmrdGtNpPWc//Z4cLsbpfp6xN85XNRA904Sd2V/stXlJ16EFcL4O7Q/SnUPOrK/9xfzSvQ44K0dFhjLCmq0nRYYyvDPrq2rfka4RnVlfoyJ9GdGZ9XUdI89aWGYbpTITNUKYgvLlnXFhSsqXd7bjEcKXc8YhvYEjhdzY/Db9lUKebP4+/b6lkK7m8nNRhDAlE8WFoREuig/DpKej+DAbGmGiuDAUUsgItWR8ejZ/yaZ3+SmET+/zL4JIC8Ys2VpI7/JPkSKEbCnkyiNPCilDyI5A+AVzSzZFVAhRBCKssVvlBRB5W8iNwZCThDymiLyT5F5iSCsh3RSRN5/cfgwxoiikkpFmgqgwopZAPmTkcxnkKCO3BZDQs088/RhykpFHRjLyP5AE70uadz/Z/1iK/2TgYwF9klJ8K4GPeJrTBXDsSXMeS3O2RA69wNEaOMAD1wTgMgJceYCLVZobH3IVBS68wLUauLwDLQKgEYG0O4CmCtC6ARpEQBsKaXYBLTWgcQe0B4EmJNLqBBqqQNsWaA4jLWig0Q2005GmPTAaAAYQyJgDGKYgIxtgMISMn4AhFzJKAwZ2yFgQGT4iI86srKysrD/rC4LWcCSWwIp+AAAAJXRFWHRkYXRlOmNyZWF0ZQAyMDIwLTAzLTAzVDEzOjU4OjU2KzAwOjAw73cO+QAAACV0RVh0ZGF0ZTptb2RpZnkAMjAyMC0wMy0wM1QxMzo1ODo1NiswMDowMJ4qtkUAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, - "execution_count": 47, + "execution_count": 48, "metadata": {}, "output_type": "execute_result" } @@ -1556,7 +1595,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.7.0" + "version": "3.6.9" } }, "nbformat": 4, diff --git a/include/openmc/geometry_aux.h b/include/openmc/geometry_aux.h index 800baef265..5ff05e7a7d 100644 --- a/include/openmc/geometry_aux.h +++ b/include/openmc/geometry_aux.h @@ -86,11 +86,15 @@ void count_cell_instances(int32_t univ_indx); //! Recursively search through universes and count universe instances. //! \param search_univ The index of the universe to begin searching from. //! \param target_univ_id The ID of the universe to be counted. +//! \param univ_count_memo Memoized counts that make this function faster for +//! large systems. The first call to this function for each target_univ_id +//! should start with an empty memo. //! \return The number of instances of target_univ_id in the geometry tree under //! search_univ. //============================================================================== -int count_universe_instances(int32_t search_univ, int32_t target_univ_id); +int count_universe_instances(int32_t search_univ, int32_t target_univ_id, + std::unordered_map& univ_count_memo); //============================================================================== //! Build a character array representing the path to a distribcell instance. diff --git a/include/openmc/lattice.h b/include/openmc/lattice.h index cd73b96408..95d30aed5a 100644 --- a/include/openmc/lattice.h +++ b/include/openmc/lattice.h @@ -77,7 +77,8 @@ public: {offsets_.resize(n_maps * universes_.size(), C_NONE);} //! Populate the distribcell offset tables. - int32_t fill_offset_table(int32_t offset, int32_t target_univ_id, int map); + int32_t fill_offset_table(int32_t offset, int32_t target_univ_id, int map, + std::unordered_map& univ_count_memo); //! \brief Check lattice indices. //! \param i_xyz[3] The indices for a lattice tile. diff --git a/include/openmc/mgxs.h b/include/openmc/mgxs.h index 8d8ca4d8e6..f9c4722e94 100644 --- a/include/openmc/mgxs.h +++ b/include/openmc/mgxs.h @@ -40,8 +40,8 @@ class Mgxs { xt::xtensor kTs; // temperature in eV (k * T) AngleDistributionType scatter_format; // flag for if this is legendre, histogram, or tabular - int num_delayed_groups; // number of delayed neutron groups int num_groups; // number of energy groups + int num_delayed_groups; // number of delayed neutron groups std::vector xs; // Cross section data // MGXS Incoming Flux Angular grid information bool is_isotropic; // used to skip search for angle indices if isotropic @@ -185,6 +185,9 @@ class Mgxs { //! @param u Incoming particle direction. void set_angle_index(Direction u); + + //! \brief Provide const access to list of XsData held by this + const std::vector& get_xsdata() const { return xs; } }; } // namespace openmc diff --git a/include/openmc/scattdata.h b/include/openmc/scattdata.h index a7ab243fa0..75fd7a2fdb 100644 --- a/include/openmc/scattdata.h +++ b/include/openmc/scattdata.h @@ -43,8 +43,8 @@ class ScattData { double_2dvec energy; // Normalized p0 matrix for sampling Eout double_2dvec mult; // nu-scatter multiplication (nu-scatt/scatt) double_3dvec dist; // Angular distribution - xt::xtensor gmin; // minimum outgoing group - xt::xtensor gmax; // maximum outgoing group + xt::xtensor gmin; // minimum outgoing group + xt::xtensor gmax; // maximum outgoing group xt::xtensor scattxs; // Isotropic Sigma_{s,g_{in}} //! \brief Calculates the value of normalized f(mu). diff --git a/openmc/data/data.py b/openmc/data/data.py index e813209f0a..5a0b16c7e3 100644 --- a/openmc/data/data.py +++ b/openmc/data/data.py @@ -10,6 +10,7 @@ from warnings import warn # pp. 293-306 (2013). The "representative isotopic abundance" values from # column 9 are used except where an interval is given, in which case the # "best measurement" is used. +# Note that the abundances are given as atomic fractions! NATURAL_ABUNDANCE = { 'H1': 0.99984426, 'H2': 0.00015574, 'He3': 0.000002, 'He4': 0.999998, 'Li6': 0.07589, 'Li7': 0.92411, @@ -110,6 +111,49 @@ NATURAL_ABUNDANCE = { 'U238': 0.992742 } +# Dictionary to give element symbols from IUPAC names +# (and some common mispellings) +ELEMENT_SYMBOL = {'neutron': 'n', 'hydrogen': 'H', 'helium': 'He', + 'lithium': 'Li', 'beryllium': 'Be', 'boron': 'B', + 'carbon': 'C', 'nitrogen': 'N', 'oxygen': 'O', 'fluorine': 'F', + 'neon': 'Ne', 'sodium': 'Na', 'magnesium': 'Mg', + 'aluminium': 'Al', 'aluminum': 'Al', 'silicon': 'Si', + 'phosphorus': 'P', 'sulfur': 'S', 'sulphur': 'S', + 'chlorine': 'Cl', 'argon': 'Ar', 'potassium': 'K', + 'calcium': 'Ca', 'scandium': 'Sc', 'titanium': 'Ti', + 'vanadium': 'V', 'chromium': 'Cr', 'manganese': 'Mn', + 'iron': 'Fe', 'cobalt': 'Co', 'nickel': 'Ni', 'copper': 'Cu', + 'zinc': 'Zn', 'gallium': 'Ga', 'germanium': 'Ge', + 'arsenic': 'As', 'selenium': 'Se', 'bromine': 'Br', + 'krypton': 'Kr', 'rubidium': 'Rb', 'strontium': 'Sr', + 'yttrium': 'Y', 'zirconium': 'Zr', 'niobium': 'Nb', + 'molybdenum': 'Mo', 'technetium': 'Tc', 'ruthenium': 'Ru', + 'rhodium': 'Rh', 'palladium': 'Pd', 'silver': 'Ag', + 'cadmium': 'Cd', 'indium': 'In', 'tin': 'Sn', 'antimony': 'Sb', + 'tellurium': 'Te', 'iodine': 'I', 'xenon': 'Xe', + 'caesium': 'Cs', 'cesium': 'Cs', 'barium': 'Ba', + 'lanthanum': 'La', 'cerium': 'Ce', 'praseodymium': 'Pr', + 'neodymium': 'Nd', 'promethium': 'Pm', 'samarium': 'Sm', + 'europium': 'Eu', 'gadolinium': 'Gd', 'terbium': 'Tb', + 'dysprosium': 'Dy', 'holmium': 'Ho', 'erbium': 'Er', + 'thulium': 'Tm', 'ytterbium': 'Yb', 'lutetium': 'Lu', + 'hafnium': 'Hf', 'tantalum': 'Ta', 'tungsten': 'W', + 'wolfram': 'W', 'rhenium': 'Re', 'osmium': 'Os', + 'iridium': 'Ir', 'platinum': 'Pt', 'gold': 'Au', + 'mercury': 'Hg', 'thallium': 'Tl', 'lead': 'Pb', + 'bismuth': 'Bi', 'polonium': 'Po', 'astatine': 'At', + 'radon': 'Rn', 'francium': 'Fr', 'radium': 'Ra', + 'actinium': 'Ac', 'thorium': 'Th', 'protactinium': 'Pa', + 'uranium': 'U', 'neptunium': 'Np', 'plutonium': 'Pu', + 'americium': 'Am', 'curium': 'Cm', 'berkelium': 'Bk', + 'californium': 'Cf', 'einsteinium': 'Es', 'fermium': 'Fm', + 'mendelevium': 'Md', 'nobelium': 'No', 'lawrencium': 'Lr', + 'rutherfordium': 'Rf', 'dubnium': 'Db', 'seaborgium': 'Sg', + 'bohrium': 'Bh', 'hassium': 'Hs', 'meitnerium': 'Mt', + 'darmstadtium': 'Ds', 'roentgenium': 'Rg', 'copernicium': 'Cn', + 'nihonium': 'Nh', 'flerovium': 'Fl', 'moscovium': 'Mc', + 'livermorium': 'Lv', 'tennessine': 'Ts', 'oganesson': 'Og'} + ATOMIC_SYMBOL = {0: 'n', 1: 'H', 2: 'He', 3: 'Li', 4: 'Be', 5: 'B', 6: 'C', 7: 'N', 8: 'O', 9: 'F', 10: 'Ne', 11: 'Na', 12: 'Mg', 13: 'Al', 14: 'Si', 15: 'P', 16: 'S', 17: 'Cl', 18: 'Ar', 19: 'K', diff --git a/openmc/element.py b/openmc/element.py index e364603348..3df2ff1f24 100644 --- a/openmc/element.py +++ b/openmc/element.py @@ -4,6 +4,7 @@ import os from xml.etree import ElementTree as ET import openmc.checkvalue as cv +from numbers import Real from openmc.data import NATURAL_ABUNDANCE, atomic_mass @@ -35,6 +36,7 @@ class Element(str): return self def expand(self, percent, percent_type, enrichment=None, + enrichment_target=None, enrichment_type=None, cross_sections=None): """Expand natural element into its naturally-occurring isotopes. @@ -52,9 +54,15 @@ class Element(str): percent_type : {'ao', 'wo'} 'ao' for atom percent and 'wo' for weight percent enrichment : float, optional - Enrichment for U235 in weight percent. For example, input 4.95 for - 4.95 weight percent enriched U. Default is None - (natural composition). + Enrichment of an enrichment_taget nuclide in percent (ao or wo). + If enrichment_taget is not supplied then it is enrichment for U235 + in weight percent. For example, input 4.95 for 4.95 weight percent + enriched U. Default is None (natural composition). + enrichment_target: str, optional + Single nuclide name to enrich from a natural composition (e.g., 'O16') + enrichment_type: {'ao', 'wo'}, optional + 'ao' for enrichment as atom percent and 'wo' for weight percent. + Default is: 'ao' for two-isotope enrichment; 'wo' for U enrichment cross_sections : str, optional Location of cross_sections.xml file. Default is None. @@ -65,16 +73,47 @@ class Element(str): is a tuple consisting of a nuclide string, the atom/weight percent, and the string 'ao' or 'wo'. + Raises + ------ + ValueError + No data is available for any of natural isotopes of the element + ValueError + If only some natural isotopes are available in the cross-section data + library and the element is not O, W, or Ta + ValueError + If a non-naturally-occurring isotope is requested + ValueError + If enrichment is requested of an element with more than two + naturally-occurring isotopes. + ValueError + If enrichment procedure for Uranium is used when element is not + Uranium. + ValueError + Uranium enrichment is requested with enrichment_type=='ao' + Notes ----- When the `enrichment` argument is specified, a correlation from `ORNL/CSD/TM-244 `_ is used to calculate the weight fractions of U234, U235, U236, and U238. Namely, the weight fraction of U234 and U236 are taken to be 0.89% and 0.46%, - respectively, of the U235 weight fraction. The remainder of the isotopic - weight is assigned to U238. + respectively, of the U235 weight fraction. The remainder of the + isotopic weight is assigned to U238. + + When the `enrichment` argument is specified with `enrichment_target`, a + general enrichment procedure is used for elements composed of exactly + two naturally-occurring isotopes. `enrichment` is interpreted as atom + percent by default but can be controlled by the `enrichment_type` + argument. """ + # Check input + if enrichment_type is not None: + cv.check_value('enrichment_type', enrichment_type, {'ao', 'wo'}) + + if enrichment is not None: + cv.check_less_than('enrichment', enrichment, 100.0, equality=True) + cv.check_greater_than('enrichment', enrichment, 0., equality=True) # Get the nuclides present in nature natural_nuclides = set() @@ -110,8 +149,8 @@ class Element(str): mutual_nuclides = sorted(list(mutual_nuclides)) absent_nuclides = sorted(list(absent_nuclides)) - # If all natural nuclides are present in the library, expand element - # using all natural nuclides + # If all natural nuclides are present in the library, + # expand element using all natural nuclides if len(absent_nuclides) == 0: for nuclide in mutual_nuclides: abundances[nuclide] = NATURAL_ABUNDANCE[nuclide] @@ -164,7 +203,20 @@ class Element(str): abundances[nuclide] = NATURAL_ABUNDANCE[nuclide] # Modify mole fractions if enrichment provided - if enrichment is not None: + # Old treatment for Uranium + if enrichment is not None and enrichment_target is None: + + # Check that the element is Uranium + if self.name != 'U': + msg = ('Enrichment procedure for Uranium was requested, ' + 'but the isotope is {} not U'.format(self)) + raise ValueError(msg) + + # Check that enrichment_type is not 'ao' + if enrichment_type == 'ao': + msg = ('Enrichment procedure for Uranium requires that ' + 'enrichment value is provided as wo%.') + raise ValueError(msg) # Calculate the mass fractions of isotopes abundances['U234'] = 0.0089 * enrichment @@ -181,6 +233,73 @@ class Element(str): for nuclide in abundances.keys(): abundances[nuclide] /= sum_abundances + # Modify mole fractions if enrichment provided + # New treatment for arbitrary element + elif enrichment is not None and enrichment_target is not None: + + # Provide more informative error message for U235 + if enrichment_target == 'U235': + msg = ("There is a special procedure for enrichment of U235 " + "in U. To invoke it, the arguments 'enrichment_target'" + "and 'enrichment_type' should be omitted. Provide " + "a value only for 'enrichment' in weight percent.") + raise ValueError(msg) + + # Check if it is two-isotope mixture + if len(abundances) != 2: + msg = ('Element {} does not consist of two naturally-occurring ' + 'isotopes. Please enter isotopic abundances manually.' + .format(self)) + raise ValueError(msg) + + # Check if the target nuclide is present in the mixture + if enrichment_target not in abundances: + msg = ('The target nuclide {} is not one of the naturally-occurring ' + 'isotopes ({})'.format(enrichment_target, list(abundances))) + raise ValueError(msg) + + # If weight percent enrichment is requested convert to mass fractions + if enrichment_type == 'wo': + # Convert the atomic abundances to weight fractions + # Compute the element atomic mass + element_am = sum(atomic_mass(nuc)*abundances[nuc] for nuc in abundances) + + # Convert Molar Fractions to mass fractions + for nuclide in abundances: + abundances[nuclide] *= atomic_mass(nuclide) / element_am + + # Normalize to one + sum_abundances = sum(abundances.values()) + for nuclide in abundances: + abundances[nuclide] /= sum_abundances + + # Enrich the mixture + # The procedure is more generic that it needs to be. It allows + # to enrich mixtures of more then 2 isotopes, keeping the ratios + # of non-enriched nuclides the same as in natural composition + + # Get fraction of non-enriched isotopes in nat. composition + non_enriched = 1.0 - abundances[enrichment_target] + tail_fraction = 1.0 - enrichment / 100.0 + + # Enrich all nuclides + # Do bogus operation for enrichment target but overwrite immediatly + # to avoid if statement in the loop + for nuclide, fraction in abundances.items(): + abundances[nuclide] = tail_fraction * fraction / non_enriched + abundances[enrichment_target] = enrichment / 100.0 + + # Convert back to atomic fractions if requested + if enrichment_type == 'wo': + # Convert the mass fractions to mole fractions + for nuclide in abundances: + abundances[nuclide] /= atomic_mass(nuclide) + + # Normalize the mole fractions to one + sum_abundances = sum(abundances.values()) + for nuclide in abundances: + abundances[nuclide] /= sum_abundances + # Compute the ratio of the nuclide atomic masses to the element # atomic mass if percent_type == 'wo': diff --git a/openmc/material.py b/openmc/material.py index 996d6a3327..3952fe4e99 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -499,34 +499,62 @@ class Material(IDManagerMixin): if macroscopic == self._macroscopic: self._macroscopic = None - def add_element(self, element, percent, percent_type='ao', enrichment=None): + def add_element(self, element, percent, percent_type='ao', enrichment=None, + enrichment_target=None, enrichment_type=None): """Add a natural element to the material Parameters ---------- element : str - Element to add, e.g., 'Zr' + Element to add, e.g., 'Zr' or 'Zirconium' percent : float Atom or weight percent percent_type : {'ao', 'wo'}, optional 'ao' for atom percent and 'wo' for weight percent. Defaults to atom percent. enrichment : float, optional - Enrichment for U235 in weight percent. For example, input 4.95 for - 4.95 weight percent enriched U. Default is None - (natural composition). + Enrichment of an enrichment_taget nuclide in percent (ao or wo). + If enrichment_taget is not supplied then it is enrichment for U235 + in weight percent. For example, input 4.95 for 4.95 weight percent + enriched U. + Default is None (natural composition). + enrichment_target: str, optional + Single nuclide name to enrich from a natural composition (e.g., 'O16') + enrichment_type: {'ao', 'wo'}, optional + 'ao' for enrichment as atom percent and 'wo' for weight percent. + Default is: 'ao' for two-isotope enrichment; 'wo' for U enrichment + + Notes + ----- + General enrichment procedure is allowed only for elements composed of + two isotopes. If `enrichment_target` is given without `enrichment` + natural composition is added to the material. """ + cv.check_type('nuclide', element, str) cv.check_type('percent', percent, Real) cv.check_value('percent type', percent_type, {'ao', 'wo'}) + # Make sure element name is just that + if not element.isalpha(): + raise ValueError("Element name should be given by the " + "element's symbol or name, e.g., 'Zr', 'zirconium'") + + # Allow for element identifier to be given as a symbol or name + if len(element) > 2: + el = element.lower() + element = openmc.data.ELEMENT_SYMBOL.get(el) + if element is None: + msg = 'Element name "{}" not recognised'.format(el) + raise ValueError(msg) + if self._macroscopic is not None: msg = 'Unable to add an Element to Material ID="{}" as a ' \ 'macroscopic data-set has already been added'.format(self._id) raise ValueError(msg) - if enrichment is not None: + if enrichment is not None and enrichment_target is None: if not isinstance(enrichment, Real): msg = 'Unable to add an Element to Material ID="{}" with a ' \ 'non-floating point enrichment value "{}"'\ @@ -551,14 +579,13 @@ class Material(IDManagerMixin): format(enrichment, self._id) warnings.warn(msg) - # Make sure element name is just that - if not element.isalpha(): - raise ValueError("Element name should be given by the " - "element's symbol, e.g., 'Zr'") - # Add naturally-occuring isotopes element = openmc.Element(element) - for nuclide in element.expand(percent, percent_type, enrichment): + for nuclide in element.expand(percent, + percent_type, + enrichment, + enrichment_target, + enrichment_type): self.add_nuclide(*nuclide) def add_s_alpha_beta(self, name, fraction=1.0): @@ -927,7 +954,7 @@ class Material(IDManagerMixin): Fractions of each material to be combined percent_type : {'ao', 'wo', 'vo'} Type of percentage, must be one of 'ao', 'wo', or 'vo', to signify atom - percent (molar percent), weight percent, or volume percent, + percent (molar percent), weight percent, or volume percent, optional. Defaults to 'ao' name : str The name for the new material, optional. Defaults to concatenated @@ -996,7 +1023,7 @@ class Material(IDManagerMixin): zip(materials, fracs)]) new_mat = openmc.Material(name=name) - # Compute atom fractions of nuclides and add them to the new material + # Compute atom fractions of nuclides and add them to the new material tot_nuclides_per_cc = np.sum([dens for dens in nuclides_per_cc.values()]) for nuc, atom_dens in nuclides_per_cc.items(): new_mat.add_nuclide(nuc, atom_dens/tot_nuclides_per_cc, 'ao') diff --git a/openmc/surface.py b/openmc/surface.py index 021f134b20..76be3b75dd 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -27,6 +27,35 @@ will not accept positional parameters for superclass arguments.\ """ +class SurfaceCoefficient: + """Descriptor class for surface coefficients. + + Parameters + ----------- + value : float or str + Value of the coefficient (float) or the name of the coefficient that + it is equivalent to (str). + + """ + def __init__(self, value): + self.value = value + + def __get__(self, instance, owner=None): + if instance is None: + return self + else: + if isinstance(self.value, str): + return instance._coefficients[self.value] + else: + return self.value + + def __set__(self, instance, value): + if isinstance(self.value, Real): + raise AttributeError('This coefficient is read-only') + check_type('{} coefficient'.format(self.value), value, Real) + instance._coefficients[self.value] = value + + def _future_kwargs_warning_helper(cls, *args, **kwargs): # Warn if Surface parameters are passed by position, not by keyword argsdict = dict(zip(('boundary_type', 'name', 'surface_id'), args)) @@ -484,7 +513,7 @@ class PlaneMixin(metaclass=ABCMeta): if np.any(np.isclose(np.abs(nhat), 1., rtol=0., atol=self._atol)): sign = nhat.sum() a, b, c, d = self._get_base_coeffs() - vals = [d/val if not np.isclose(val, 0., rtol=0., atol=self._atol) + vals = [d/val if not np.isclose(val, 0., rtol=0., atol=self._atol) else np.nan for val in (a, b, c)] if side == '-': if sign > 0: @@ -674,41 +703,10 @@ class Plane(PlaneMixin, Surface): return True return NotImplemented - @property - def a(self): - return self.coefficients['a'] - - @property - def b(self): - return self.coefficients['b'] - - @property - def c(self): - return self.coefficients['c'] - - @property - def d(self): - return self.coefficients['d'] - - @a.setter - def a(self, a): - check_type('A coefficient', a, Real) - self._coefficients['a'] = a - - @b.setter - def b(self, b): - check_type('B coefficient', b, Real) - self._coefficients['b'] = b - - @c.setter - def c(self, c): - check_type('C coefficient', c, Real) - self._coefficients['c'] = c - - @d.setter - def d(self, d): - check_type('D coefficient', d, Real) - self._coefficients['d'] = d + a = SurfaceCoefficient('a') + b = SurfaceCoefficient('b') + c = SurfaceCoefficient('c') + d = SurfaceCoefficient('d') @classmethod def from_points(cls, p1, p2, p3, **kwargs): @@ -792,30 +790,11 @@ class XPlane(PlaneMixin, Surface): super().__init__(**kwargs) self.x0 = x0 - @property - def x0(self): - return self.coefficients['x0'] - - @property - def a(self): - return 1. - - @property - def b(self): - return 0. - - @property - def c(self): - return 0. - - @property - def d(self): - return self.x0 - - @x0.setter - def x0(self, x0): - check_type('x0 coefficient', x0, Real) - self._coefficients['x0'] = x0 + x0 = SurfaceCoefficient('x0') + a = SurfaceCoefficient(1.) + b = SurfaceCoefficient(0.) + c = SurfaceCoefficient(0.) + d = x0 def evaluate(self, point): return point[0] - self.x0 @@ -870,30 +849,11 @@ class YPlane(PlaneMixin, Surface): super().__init__(**kwargs) self.y0 = y0 - @property - def y0(self): - return self.coefficients['y0'] - - @property - def a(self): - return 0. - - @property - def b(self): - return 1. - - @property - def c(self): - return 0. - - @property - def d(self): - return self.y0 - - @y0.setter - def y0(self, y0): - check_type('y0 coefficient', y0, Real) - self._coefficients['y0'] = y0 + y0 = SurfaceCoefficient('y0') + a = SurfaceCoefficient(0.) + b = SurfaceCoefficient(1.) + c = SurfaceCoefficient(0.) + d = y0 def evaluate(self, point): return point[1] - self.y0 @@ -948,30 +908,11 @@ class ZPlane(PlaneMixin, Surface): super().__init__(**kwargs) self.z0 = z0 - @property - def z0(self): - return self.coefficients['z0'] - - @property - def a(self): - return 0. - - @property - def b(self): - return 0. - - @property - def c(self): - return 1. - - @property - def d(self): - return self.z0 - - @z0.setter - def z0(self, z0): - check_type('z0 coefficient', z0, Real) - self._coefficients['z0'] = z0 + z0 = SurfaceCoefficient('z0') + a = SurfaceCoefficient(0.) + b = SurfaceCoefficient(0.) + c = SurfaceCoefficient(1.) + d = z0 def evaluate(self, point): return point[2] - self.z0 @@ -1095,7 +1036,7 @@ class QuadricMixin(metaclass=ABCMeta): return surf def rotate(self, rotation, pivot=(0., 0., 0.), order='xyz', inplace=False): - # Get pivot and rotation matrix + # Get pivot and rotation matrix pivot = np.asarray(pivot) rotation = np.asarray(rotation, dtype=float) @@ -1225,68 +1166,13 @@ class Cylinder(QuadricMixin, Surface): return True return NotImplemented - @property - def x0(self): - return self.coefficients['x0'] - - @property - def y0(self): - return self.coefficients['y0'] - - @property - def z0(self): - return self.coefficients['z0'] - - @property - def r(self): - return self.coefficients['r'] - - @property - def dx(self): - return self.coefficients['dx'] - - @property - def dy(self): - return self.coefficients['dy'] - - @property - def dz(self): - return self.coefficients['dz'] - - @x0.setter - def x0(self, x0): - check_type('x0 coefficient', x0, Real) - self._coefficients['x0'] = x0 - - @y0.setter - def y0(self, y0): - check_type('y0 coefficient', y0, Real) - self._coefficients['y0'] = y0 - - @z0.setter - def z0(self, z0): - check_type('z0 coefficient', z0, Real) - self._coefficients['z0'] = z0 - - @r.setter - def r(self, r): - check_type('r coefficient', r, Real) - self._coefficients['r'] = r - - @dx.setter - def dx(self, dx): - check_type('dx coefficient', dx, Real) - self._coefficients['dx'] = dx - - @dy.setter - def dy(self, dy): - check_type('dy coefficient', dy, Real) - self._coefficients['dy'] = dy - - @dz.setter - def dz(self, dz): - check_type('dz coefficient', dz, Real) - self._coefficients['dz'] = dz + x0 = SurfaceCoefficient('x0') + y0 = SurfaceCoefficient('y0') + z0 = SurfaceCoefficient('z0') + r = SurfaceCoefficient('r') + dx = SurfaceCoefficient('dx') + dy = SurfaceCoefficient('dy') + dz = SurfaceCoefficient('dz') def bounding_box(self, side): if side == '-': @@ -1305,7 +1191,7 @@ class Cylinder(QuadricMixin, Surface): # Get x, y, z coordinates of two points x1, y1, z1 = self._origin x2, y2, z2 = self._origin + self._axis - r = self.r + r = self.r # Define intermediate terms dx = x2 - x1 @@ -1375,7 +1261,7 @@ class Cylinder(QuadricMixin, Surface): with catch_warnings(): simplefilter('ignore', IDWarning) kwargs = {'boundary_type': self.boundary_type, 'name': self.name, - 'surface_id': self.id} + 'surface_id': self.id} quad_rep = Quadric(*self._get_base_coeffs(), **kwargs) return quad_rep.to_xml_element() @@ -1440,48 +1326,13 @@ class XCylinder(QuadricMixin, Surface): for key, val in zip(self._coeff_keys, (y0, z0, r)): setattr(self, key, val) - @property - def y0(self): - return self.coefficients['y0'] - - @property - def z0(self): - return self.coefficients['z0'] - - @property - def r(self): - return self.coefficients['r'] - - @property - def x0(self): - return 0. - - @property - def dx(self): - return 1. - - @property - def dy(self): - return 0. - - @property - def dz(self): - return 0. - - @y0.setter - def y0(self, y0): - check_type('y0 coefficient', y0, Real) - self._coefficients['y0'] = y0 - - @z0.setter - def z0(self, z0): - check_type('z0 coefficient', z0, Real) - self._coefficients['z0'] = z0 - - @r.setter - def r(self, r): - check_type('r coefficient', r, Real) - self._coefficients['r'] = r + x0 = SurfaceCoefficient(0.) + y0 = SurfaceCoefficient('y0') + z0 = SurfaceCoefficient('z0') + r = SurfaceCoefficient('r') + dx = SurfaceCoefficient(1.) + dy = SurfaceCoefficient(0.) + dz = SurfaceCoefficient(0.) def _get_base_coeffs(self): y0, z0, r = self.y0, self.z0, self.r @@ -1566,48 +1417,13 @@ class YCylinder(QuadricMixin, Surface): for key, val in zip(self._coeff_keys, (x0, z0, r)): setattr(self, key, val) - @property - def x0(self): - return self.coefficients['x0'] - - @property - def z0(self): - return self.coefficients['z0'] - - @property - def r(self): - return self.coefficients['r'] - - @property - def y0(self): - return 0. - - @property - def dx(self): - return 0. - - @property - def dy(self): - return 1. - - @property - def dz(self): - return 0. - - @x0.setter - def x0(self, x0): - check_type('x0 coefficient', x0, Real) - self._coefficients['x0'] = x0 - - @z0.setter - def z0(self, z0): - check_type('z0 coefficient', z0, Real) - self._coefficients['z0'] = z0 - - @r.setter - def r(self, r): - check_type('r coefficient', r, Real) - self._coefficients['r'] = r + x0 = SurfaceCoefficient('x0') + y0 = SurfaceCoefficient(0.) + z0 = SurfaceCoefficient('z0') + r = SurfaceCoefficient('r') + dx = SurfaceCoefficient(0.) + dy = SurfaceCoefficient(1.) + dz = SurfaceCoefficient(0.) def _get_base_coeffs(self): x0, z0, r = self.x0, self.z0, self.r @@ -1692,48 +1508,13 @@ class ZCylinder(QuadricMixin, Surface): for key, val in zip(self._coeff_keys, (x0, y0, r)): setattr(self, key, val) - @property - def x0(self): - return self.coefficients['x0'] - - @property - def y0(self): - return self.coefficients['y0'] - - @property - def r(self): - return self.coefficients['r'] - - @property - def z0(self): - return 0. - - @property - def dx(self): - return 0. - - @property - def dy(self): - return 0. - - @property - def dz(self): - return 1. - - @x0.setter - def x0(self, x0): - check_type('x0 coefficient', x0, Real) - self._coefficients['x0'] = x0 - - @y0.setter - def y0(self, y0): - check_type('y0 coefficient', y0, Real) - self._coefficients['y0'] = y0 - - @r.setter - def r(self, r): - check_type('r coefficient', r, Real) - self._coefficients['r'] = r + x0 = SurfaceCoefficient('x0') + y0 = SurfaceCoefficient('y0') + z0 = SurfaceCoefficient(0.) + r = SurfaceCoefficient('r') + dx = SurfaceCoefficient(0.) + dy = SurfaceCoefficient(0.) + dz = SurfaceCoefficient(1.) def _get_base_coeffs(self): x0, y0, r = self.x0, self.y0, self.r @@ -1820,41 +1601,10 @@ class Sphere(QuadricMixin, Surface): for key, val in zip(self._coeff_keys, (x0, y0, z0, r)): setattr(self, key, val) - @property - def x0(self): - return self.coefficients['x0'] - - @property - def y0(self): - return self.coefficients['y0'] - - @property - def z0(self): - return self.coefficients['z0'] - - @property - def r(self): - return self.coefficients['r'] - - @x0.setter - def x0(self, x0): - check_type('x0 coefficient', x0, Real) - self._coefficients['x0'] = x0 - - @y0.setter - def y0(self, y0): - check_type('y0 coefficient', y0, Real) - self._coefficients['y0'] = y0 - - @z0.setter - def z0(self, z0): - check_type('z0 coefficient', z0, Real) - self._coefficients['z0'] = z0 - - @r.setter - def r(self, r): - check_type('r coefficient', r, Real) - self._coefficients['r'] = r + x0 = SurfaceCoefficient('x0') + y0 = SurfaceCoefficient('y0') + z0 = SurfaceCoefficient('z0') + r = SurfaceCoefficient('r') def _get_base_coeffs(self): x0, y0, z0, r = self.x0, self.y0, self.z0, self.r @@ -1966,68 +1716,13 @@ class Cone(QuadricMixin, Surface): return True return NotImplemented - @property - def x0(self): - return self.coefficients['x0'] - - @property - def y0(self): - return self.coefficients['y0'] - - @property - def z0(self): - return self.coefficients['z0'] - - @property - def r2(self): - return self.coefficients['r2'] - - @property - def dx(self): - return self.coefficients['dx'] - - @property - def dy(self): - return self.coefficients['dy'] - - @property - def dz(self): - return self.coefficients['dz'] - - @x0.setter - def x0(self, x0): - check_type('x0 coefficient', x0, Real) - self._coefficients['x0'] = x0 - - @y0.setter - def y0(self, y0): - check_type('y0 coefficient', y0, Real) - self._coefficients['y0'] = y0 - - @z0.setter - def z0(self, z0): - check_type('z0 coefficient', z0, Real) - self._coefficients['z0'] = z0 - - @r2.setter - def r2(self, r2): - check_type('r^2 coefficient', r2, Real) - self._coefficients['r2'] = r2 - - @dx.setter - def dx(self, dx): - check_type('dx coefficient', dx, Real) - self._coefficients['dx'] = dx - - @dy.setter - def dy(self, dy): - check_type('dy coefficient', dy, Real) - self._coefficients['dy'] = dy - - @dz.setter - def dz(self, dz): - check_type('dz coefficient', dz, Real) - self._coefficients['dz'] = dz + x0 = SurfaceCoefficient('x0') + y0 = SurfaceCoefficient('y0') + z0 = SurfaceCoefficient('z0') + r2 = SurfaceCoefficient('r2') + dx = SurfaceCoefficient('dx') + dy = SurfaceCoefficient('dy') + dz = SurfaceCoefficient('dz') def _get_base_coeffs(self): # The equation for a general cone with vertex at point p = (x0, y0, z0) @@ -2074,7 +1769,7 @@ class Cone(QuadricMixin, Surface): with catch_warnings(): simplefilter('ignore', IDWarning) kwargs = {'boundary_type': self.boundary_type, 'name': self.name, - 'surface_id': self.id} + 'surface_id': self.id} quad_rep = Quadric(*self._get_base_coeffs(), **kwargs) return quad_rep.to_xml_element() @@ -2142,53 +1837,13 @@ class XCone(QuadricMixin, Surface): for key, val in zip(self._coeff_keys, (x0, y0, z0, r2)): setattr(self, key, val) - @property - def x0(self): - return self.coefficients['x0'] - - @property - def y0(self): - return self.coefficients['y0'] - - @property - def z0(self): - return self.coefficients['z0'] - - @property - def r2(self): - return self.coefficients['r2'] - - @property - def dx(self): - return 1. - - @property - def dy(self): - return 0. - - @property - def dz(self): - return 0. - - @x0.setter - def x0(self, x0): - check_type('x0 coefficient', x0, Real) - self._coefficients['x0'] = x0 - - @y0.setter - def y0(self, y0): - check_type('y0 coefficient', y0, Real) - self._coefficients['y0'] = y0 - - @z0.setter - def z0(self, z0): - check_type('z0 coefficient', z0, Real) - self._coefficients['z0'] = z0 - - @r2.setter - def r2(self, r2): - check_type('r^2 coefficient', r2, Real) - self._coefficients['r2'] = r2 + x0 = SurfaceCoefficient('x0') + y0 = SurfaceCoefficient('y0') + z0 = SurfaceCoefficient('z0') + r2 = SurfaceCoefficient('r2') + dx = SurfaceCoefficient(1.) + dy = SurfaceCoefficient(0.) + dz = SurfaceCoefficient(0.) def _get_base_coeffs(self): x0, y0, z0, r2 = self.x0, self.y0, self.z0, self.r2 @@ -2271,53 +1926,13 @@ class YCone(QuadricMixin, Surface): for key, val in zip(self._coeff_keys, (x0, y0, z0, r2)): setattr(self, key, val) - @property - def x0(self): - return self.coefficients['x0'] - - @property - def y0(self): - return self.coefficients['y0'] - - @property - def z0(self): - return self.coefficients['z0'] - - @property - def r2(self): - return self.coefficients['r2'] - - @property - def dx(self): - return 0. - - @property - def dy(self): - return 1. - - @property - def dz(self): - return 0. - - @x0.setter - def x0(self, x0): - check_type('x0 coefficient', x0, Real) - self._coefficients['x0'] = x0 - - @y0.setter - def y0(self, y0): - check_type('y0 coefficient', y0, Real) - self._coefficients['y0'] = y0 - - @z0.setter - def z0(self, z0): - check_type('z0 coefficient', z0, Real) - self._coefficients['z0'] = z0 - - @r2.setter - def r2(self, r2): - check_type('r^2 coefficient', r2, Real) - self._coefficients['r2'] = r2 + x0 = SurfaceCoefficient('x0') + y0 = SurfaceCoefficient('y0') + z0 = SurfaceCoefficient('z0') + r2 = SurfaceCoefficient('r2') + dx = SurfaceCoefficient(0.) + dy = SurfaceCoefficient(1.) + dz = SurfaceCoefficient(0.) def _get_base_coeffs(self): x0, y0, z0, r2 = self.x0, self.y0, self.z0, self.r2 @@ -2400,53 +2015,13 @@ class ZCone(QuadricMixin, Surface): for key, val in zip(self._coeff_keys, (x0, y0, z0, r2)): setattr(self, key, val) - @property - def x0(self): - return self.coefficients['x0'] - - @property - def y0(self): - return self.coefficients['y0'] - - @property - def z0(self): - return self.coefficients['z0'] - - @property - def r2(self): - return self.coefficients['r2'] - - @property - def dx(self): - return 0. - - @property - def dy(self): - return 0. - - @property - def dz(self): - return 1. - - @x0.setter - def x0(self, x0): - check_type('x0 coefficient', x0, Real) - self._coefficients['x0'] = x0 - - @y0.setter - def y0(self, y0): - check_type('y0 coefficient', y0, Real) - self._coefficients['y0'] = y0 - - @z0.setter - def z0(self, z0): - check_type('z0 coefficient', z0, Real) - self._coefficients['z0'] = z0 - - @r2.setter - def r2(self, r2): - check_type('r^2 coefficient', r2, Real) - self._coefficients['r2'] = r2 + x0 = SurfaceCoefficient('x0') + y0 = SurfaceCoefficient('y0') + z0 = SurfaceCoefficient('z0') + r2 = SurfaceCoefficient('r2') + dx = SurfaceCoefficient(0.) + dy = SurfaceCoefficient(0.) + dz = SurfaceCoefficient(1.) def _get_base_coeffs(self): x0, y0, z0, r2 = self.x0, self.y0, self.z0, self.r2 @@ -2513,95 +2088,16 @@ class Quadric(QuadricMixin, Surface): for key, val in zip(self._coeff_keys, (a, b, c, d, e, f, g, h, j, k)): setattr(self, key, val) - @property - def a(self): - return self.coefficients['a'] - - @property - def b(self): - return self.coefficients['b'] - - @property - def c(self): - return self.coefficients['c'] - - @property - def d(self): - return self.coefficients['d'] - - @property - def e(self): - return self.coefficients['e'] - - @property - def f(self): - return self.coefficients['f'] - - @property - def g(self): - return self.coefficients['g'] - - @property - def h(self): - return self.coefficients['h'] - - @property - def j(self): - return self.coefficients['j'] - - @property - def k(self): - return self.coefficients['k'] - - @a.setter - def a(self, a): - check_type('a coefficient', a, Real) - self._coefficients['a'] = a - - @b.setter - def b(self, b): - check_type('b coefficient', b, Real) - self._coefficients['b'] = b - - @c.setter - def c(self, c): - check_type('c coefficient', c, Real) - self._coefficients['c'] = c - - @d.setter - def d(self, d): - check_type('d coefficient', d, Real) - self._coefficients['d'] = d - - @e.setter - def e(self, e): - check_type('e coefficient', e, Real) - self._coefficients['e'] = e - - @f.setter - def f(self, f): - check_type('f coefficient', f, Real) - self._coefficients['f'] = f - - @g.setter - def g(self, g): - check_type('g coefficient', g, Real) - self._coefficients['g'] = g - - @h.setter - def h(self, h): - check_type('h coefficient', h, Real) - self._coefficients['h'] = h - - @j.setter - def j(self, j): - check_type('j coefficient', j, Real) - self._coefficients['j'] = j - - @k.setter - def k(self, k): - check_type('k coefficient', k, Real) - self._coefficients['k'] = k + a = SurfaceCoefficient('a') + b = SurfaceCoefficient('b') + c = SurfaceCoefficient('c') + d = SurfaceCoefficient('d') + e = SurfaceCoefficient('e') + f = SurfaceCoefficient('f') + g = SurfaceCoefficient('g') + h = SurfaceCoefficient('h') + j = SurfaceCoefficient('j') + k = SurfaceCoefficient('k') def _get_base_coeffs(self): return tuple(getattr(self, c) for c in self._coeff_keys) @@ -2734,7 +2230,7 @@ class Halfspace(Region): Parameters ---------- redundant_surfaces : dict - Dictionary mapping redundant surface IDs to surface IDs for the + Dictionary mapping redundant surface IDs to surface IDs for the :class:`openmc.Surface` instances that should replace them. """ diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index ae6d894aa8..c58af934b4 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -385,8 +385,10 @@ prepare_distribcell() } // Fill the cell and lattice offset tables. + #pragma omp parallel for for (int map = 0; map < target_univ_ids.size(); map++) { auto target_univ_id = target_univ_ids[map]; + std::unordered_map univ_count_memo; for (const auto& univ : model::universes) { int32_t offset = 0; for (int32_t cell_indx : univ->cells_) { @@ -395,11 +397,13 @@ prepare_distribcell() if (c.type_ == Fill::UNIVERSE) { c.offset_[map] = offset; int32_t search_univ = c.fill_; - offset += count_universe_instances(search_univ, target_univ_id); + offset += count_universe_instances(search_univ, target_univ_id, + univ_count_memo); } else if (c.type_ == Fill::LATTICE) { Lattice& lat = *model::lattices[c.fill_]; - offset = lat.fill_offset_table(offset, target_univ_id, map); + offset = lat.fill_offset_table(offset, target_univ_id, map, + univ_count_memo); } } } @@ -411,7 +415,6 @@ prepare_distribcell() void count_cell_instances(int32_t univ_indx) { - const auto univ_counts = model::universe_cell_counts.find(univ_indx); if (univ_counts != model::universe_cell_counts.end()) { for (const auto& it : univ_counts->second) { @@ -442,30 +445,42 @@ count_cell_instances(int32_t univ_indx) //============================================================================== int -count_universe_instances(int32_t search_univ, int32_t target_univ_id) +count_universe_instances(int32_t search_univ, int32_t target_univ_id, + std::unordered_map& univ_count_memo) { - // If this is the target, it can't contain itself. + // If this is the target, it can't contain itself. if (model::universes[search_univ]->id_ == target_univ_id) { return 1; } + // If we have already counted the number of instances, reuse that value. + auto search = univ_count_memo.find(search_univ); + if (search != univ_count_memo.end()) { + return search->second; + } + int count {0}; for (int32_t cell_indx : model::universes[search_univ]->cells_) { Cell& c = *model::cells[cell_indx]; if (c.type_ == Fill::UNIVERSE) { int32_t next_univ = c.fill_; - count += count_universe_instances(next_univ, target_univ_id); + count += count_universe_instances(next_univ, target_univ_id, + univ_count_memo); } else if (c.type_ == Fill::LATTICE) { Lattice& lat = *model::lattices[c.fill_]; for (auto it = lat.begin(); it != lat.end(); ++it) { int32_t next_univ = *it; - count += count_universe_instances(next_univ, target_univ_id); + count += count_universe_instances(next_univ, target_univ_id, + univ_count_memo); } } } + // Remember the number of instances in this universe. + univ_count_memo[search_univ] = count; + return count; } diff --git a/src/lattice.cpp b/src/lattice.cpp index 8f88d8c591..eba8807c1e 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -93,11 +93,12 @@ Lattice::adjust_indices() //============================================================================== int32_t -Lattice::fill_offset_table(int32_t offset, int32_t target_univ_id, int map) +Lattice::fill_offset_table(int32_t offset, int32_t target_univ_id, int map, + std::unordered_map& univ_count_memo) { for (LatticeIter it = begin(); it != end(); ++it) { offsets_[map * universes_.size() + it.indx_] = offset; - offset += count_universe_instances(*it, target_univ_id); + offset += count_universe_instances(*it, target_univ_id, univ_count_memo); } return offset; } diff --git a/src/main.cpp b/src/main.cpp index dea347ce4b..ccd2005bf5 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -43,6 +43,8 @@ int main(int argc, char* argv[]) { case RunMode::VOLUME: err = openmc_calculate_volumes(); break; + default: + break; } if (err) fatal_error(openmc_err_msg); diff --git a/src/particle.cpp b/src/particle.cpp index ea0d03cd01..5f9713941e 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -675,6 +675,8 @@ Particle::write_restart() const case RunMode::PARTICLE: write_dataset(file_id, "run_mode", "particle restart"); break; + default: + break; } write_dataset(file_id, "id", id_); write_dataset(file_id, "type", static_cast(type_)); diff --git a/src/scattdata.cpp b/src/scattdata.cpp index 11e5c4c4c8..4fd7a7b7c4 100644 --- a/src/scattdata.cpp +++ b/src/scattdata.cpp @@ -227,6 +227,8 @@ ScattData::get_xs(MgxsType xstype, int gin, const int* gout, const double* mu) fatal_error("Invalid call to get_xs"); } break; + default: + break; } return val; } diff --git a/src/state_point.cpp b/src/state_point.cpp index 503b1354e6..d077590365 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -87,6 +87,8 @@ openmc_statepoint_write(const char* filename, bool* write_source) case RunMode::EIGENVALUE: write_dataset(file_id, "run_mode", "eigenvalue"); break; + default: + break; } write_attribute(file_id, "photon_transport", settings::photon_transport); write_dataset(file_id, "n_particles", settings::n_particles); diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 81f6c504eb..e0cd5d095c 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -27,7 +27,7 @@ namespace openmc { //============================================================================== FilterBinIter::FilterBinIter(const Tally& tally, Particle* p) - : tally_{tally}, filter_matches_{p->filter_matches_} + : filter_matches_{p->filter_matches_}, tally_{tally} { // Find all valid bins in each relevant filter if they have not already been // found for this event. @@ -57,7 +57,7 @@ FilterBinIter::FilterBinIter(const Tally& tally, Particle* p) FilterBinIter::FilterBinIter(const Tally& tally, bool end, std::vector* particle_filter_matches) - : tally_{tally}, filter_matches_{*particle_filter_matches} + : filter_matches_{*particle_filter_matches}, tally_{tally} { // Handle the special case for an iterator that points to the end. if (end) { diff --git a/src/volume_calc.cpp b/src/volume_calc.cpp index 8cae7743d3..97594c3b5d 100644 --- a/src/volume_calc.cpp +++ b/src/volume_calc.cpp @@ -295,6 +295,8 @@ std::vector VolumeCalculation::execute() const case TriggerMetric::variance: val = result.volume[1] * result.volume[1]; break; + default: + break; } // update max if entry is valid if (val > 0.0) { trigger_val = std::max(trigger_val, val); } @@ -374,6 +376,8 @@ void VolumeCalculation::to_hdf5(const std::string& filename, case TriggerMetric::relative_error: trigger_str = "rel_err"; break; + default: + break; } write_attribute(file_id, "trigger_type", trigger_str); } else { diff --git a/tests/unit_tests/test_element.py b/tests/unit_tests/test_element.py new file mode 100644 index 0000000000..bacb988b9a --- /dev/null +++ b/tests/unit_tests/test_element.py @@ -0,0 +1,81 @@ +import openmc +from pytest import approx, raises + +from openmc.data import NATURAL_ABUNDANCE, atomic_mass + + +def test_expand_no_enrichment(): + """ Expand Li in natural compositions""" + lithium = openmc.Element('Li') + + # Verify the expansion into ATOMIC fraction against natural composition + for isotope in lithium.expand(100.0, 'ao'): + assert isotope[1] == approx(NATURAL_ABUNDANCE[isotope[0]] * 100.0) + + # Verify the expansion into WEIGHT fraction against natural composition + natural = {'Li6': NATURAL_ABUNDANCE['Li6'] * atomic_mass('Li6'), + 'Li7': NATURAL_ABUNDANCE['Li7'] * atomic_mass('Li7')} + li_am = sum(natural.values()) + for key in natural: + natural[key] /= li_am + + for isotope in lithium.expand(100.0, 'wo'): + assert isotope[1] == approx(natural[isotope[0]] * 100.0) + + +def test_expand_enrichment(): + """ Expand and verify enrichment of Li """ + lithium = openmc.Element('Li') + + # Verify the enrichment by atoms + ref = {'Li6': 75.0, 'Li7': 25.0} + for isotope in lithium.expand(100.0, 'ao', 25.0, 'Li7', 'ao'): + assert isotope[1] == approx(ref[isotope[0]]) + + # Verify the enrichment by weight + for isotope in lithium.expand(100.0, 'wo', 25.0, 'Li7', 'wo'): + assert isotope[1] == approx(ref[isotope[0]]) + + +def test_expand_exceptions(): + """ Test that correct exceptions are raised for invalid input """ + + # 1 Isotope Element + with raises(ValueError): + element = openmc.Element('Be') + element.expand(70.0, 'ao', 4.0, 'Be9') + + # 3 Isotope Element + with raises(ValueError): + element = openmc.Element('Cr') + element.expand(70.0, 'ao', 4.0, 'Cr52') + + # Non-present Enrichment Target + with raises(ValueError): + element = openmc.Element('H') + element.expand(70.0, 'ao', 4.0, 'H4') + + # Enrichment Procedure for Uranium if not Uranium + with raises(ValueError): + element = openmc.Element('Li') + element.expand(70.0, 'ao', 4.0) + + # Missing Enrichment Target + with raises(ValueError): + element = openmc.Element('Li') + element.expand(70.0, 'ao', 4.0, enrichment_type='ao') + + # Invalid Enrichment Type Entry + with raises(ValueError): + element = openmc.Element('Li') + element.expand(70.0, 'ao', 4.0, 'Li7', 'Grand Moff Tarkin') + + # Trying to enrich Uranium + with raises(ValueError): + element = openmc.Element('U') + element.expand(70.0, 'ao', 4.0, 'U235', 'wo') + + # Trying to enrich Uranium with wrong enrichment_target + with raises(ValueError): + element = openmc.Element('U') + element.expand(70.0, 'ao', 4.0, enrichment_type='ao') diff --git a/tests/unit_tests/test_material.py b/tests/unit_tests/test_material.py index cbff16e403..6f51135ab2 100644 --- a/tests/unit_tests/test_material.py +++ b/tests/unit_tests/test_material.py @@ -29,11 +29,34 @@ def test_elements(): m = openmc.Material() m.add_element('Zr', 1.0) m.add_element('U', 1.0, enrichment=4.5) + m.add_element('Li', 1.0, enrichment=60.0, enrichment_target='Li7') + m.add_element('H', 1.0, enrichment=50.0, enrichment_target='H2', + enrichment_type='wo') with pytest.raises(ValueError): m.add_element('U', 1.0, enrichment=100.0) with pytest.raises(ValueError): m.add_element('Pu', 1.0, enrichment=3.0) + with pytest.raises(ValueError): + m.add_element('U', 1.0, enrichment=70.0, enrichment_target='U235') + with pytest.raises(ValueError): + m.add_element('He', 1.0, enrichment=17.0, enrichment_target='He6') +def test_elements_by_name(): + """Test adding elements by name""" + m = openmc.Material() + m.add_element('woLfrAm', 1.0) + with pytest.raises(ValueError): + m.add_element('uranum', 1.0) + m.add_element('uRaNiUm', 1.0) + m.add_element('Aluminium', 1.0) + a = openmc.Material() + b = openmc.Material() + c = openmc.Material() + a.add_element('sulfur', 1.0) + b.add_element('SulPhUR', 1.0) + c.add_element('S', 1.0) + assert a._nuclides == b._nuclides + assert b._nuclides == c._nuclides def test_density(): m = openmc.Material() diff --git a/vendor/xtensor b/vendor/xtensor index ef091807f7..31acec1e90 160000 --- a/vendor/xtensor +++ b/vendor/xtensor @@ -1 +1 @@ -Subproject commit ef091807f7ed0e5ba7e251a6c46f4af7bba79e2e +Subproject commit 31acec1e90bbea6d4bc17af0710a123bd5da6689 diff --git a/vendor/xtl b/vendor/xtl index f5d13e6c4f..0024346605 160000 --- a/vendor/xtl +++ b/vendor/xtl @@ -1 +1 @@ -Subproject commit f5d13e6c4f856becc178939365fcdcf9a657ffb5 +Subproject commit 0024346605bd92bcc4009caad7f4be88687e063a