diff --git a/docs/source/pythonapi/index.rst b/docs/source/pythonapi/index.rst index 3c28a2185..8abb52528 100644 --- a/docs/source/pythonapi/index.rst +++ b/docs/source/pythonapi/index.rst @@ -22,7 +22,7 @@ there are many substantial benefits to using the Python API, including: - Ability to plot individual universes as geometry is being created - A :math:`k_\text{eff}` search function (:func:`openmc.search_for_keff`) - Random sphere packing for generating TRISO particle locations - (:func:`openmc.model.pack_trisos`) + (:func:`openmc.model.pack_spheres`) - Ability to create materials based on natural elements or uranium enrichment For those new to Python, there are many good tutorials available online. We diff --git a/docs/source/pythonapi/model.rst b/docs/source/pythonapi/model.rst index 9ed77f366..90e8b779d 100644 --- a/docs/source/pythonapi/model.rst +++ b/docs/source/pythonapi/model.rst @@ -37,7 +37,7 @@ Functions :template: myfunction.rst openmc.model.create_triso_lattice - openmc.model.pack_trisos + openmc.model.pack_spheres Model Container --------------- diff --git a/examples/jupyter/triso.ipynb b/examples/jupyter/triso.ipynb index bc7801e2d..36e0c1f14 100644 --- a/examples/jupyter/triso.ipynb +++ b/examples/jupyter/triso.ipynb @@ -101,7 +101,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Now that we have a universe that can be used for each TRISO particle, we need to randomly select locations. In this example, we will select locations at random within a 1 cm x 1 cm x 1 cm box centered at the origin with a packing fraction of 30%. Note that `pack_trisos` can handle up to the theoretical maximum of 60% (it will just be slow)." + "Next, we need a region to pack the TRISO particles in. We will use a 1 cm x 1 cm x 1 cm box centered at the origin." ] }, { @@ -111,16 +111,49 @@ "collapsed": false }, "outputs": [], + "source": [ + "min_x = openmc.XPlane(x0=-0.5, boundary_type='reflective')\n", + "max_x = openmc.XPlane(x0=0.5, boundary_type='reflective')\n", + "min_y = openmc.YPlane(y0=-0.5, boundary_type='reflective')\n", + "max_y = openmc.YPlane(y0=0.5, boundary_type='reflective')\n", + "min_z = openmc.ZPlane(z0=-0.5, boundary_type='reflective')\n", + "max_z = openmc.ZPlane(z0=0.5, boundary_type='reflective')\n", + "region = +min_x & -max_x & +min_y & -max_y & +min_z & -max_z" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we need to randomly select locations for the TRISO particles. In this example, we will select locations at random within the box with a packing fraction of 30%. Note that `pack_spheres` can handle up to the theoretical maximum of 60% (it will just be slow)." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], "source": [ "outer_radius = 425.*1e-4\n", - "\n", - "trisos = openmc.model.pack_trisos(\n", - " radius=outer_radius,\n", - " fill=triso_univ,\n", - " domain_shape='cube',\n", - " domain_length=1,\n", - " packing_fraction=0.3\n", - ")" + "centers = openmc.model.pack_spheres(radius=outer_radius, region=region, pf=0.3)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now that we have the locations of the TRISO particles determined and a universe that can be used for each particle, we can create the TRISO particles." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "trisos = [openmc.model.TRISO(outer_radius, triso_univ, c) for c in centers]" ] }, { @@ -132,7 +165,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 7, "metadata": { "collapsed": false }, @@ -142,10 +175,10 @@ "output_type": "stream", "text": [ "Cell\n", - "\tID =\t10005\n", + "\tID =\t6\n", "\tName =\t\n", - "\tFill =\t10000\n", - "\tRegion =\t-10004\n", + "\tFill =\t1\n", + "\tRegion =\t-11\n", "\tRotation =\tNone\n", "\tTranslation =\t[-0.33455672 0.31790187 0.24135378]\n", "\n" @@ -165,7 +198,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 8, "metadata": { "collapsed": false }, @@ -175,7 +208,7 @@ "output_type": "stream", "text": [ "[-0.45718713 -0.45730405 -0.45725048]\n", - "[ 0.45705454 0.45743843 0.45741142]\n" + "[0.45705454 0.45743843 0.45741142]\n" ] } ], @@ -194,7 +227,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 9, "metadata": { "collapsed": false }, @@ -205,7 +238,7 @@ "0.2996893513959326" ] }, - "execution_count": 7, + "execution_count": 9, "metadata": {}, "output_type": "execute_result" } @@ -218,42 +251,16 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Now that we have our TRISO particles created, we need to place them in a lattice to provide optimal tracking performance in OpenMC. We'll start by creating a box that the lattice will be placed within." + "Now that we have our TRISO particles created, we need to place them in a lattice to provide optimal tracking performance in OpenMC. We can use the box we created above to place the lattice in. Actually creating a lattice containing TRISO particles can be done with the `model.create_triso_lattice()` function. This function requires that we give it a list of TRISO particles, the lower-left coordinates of the lattice, the pitch of each lattice cell, the overall shape of the lattice (number of cells in each direction), and a background material." ] }, { "cell_type": "code", - "execution_count": 8, - "metadata": { - "collapsed": false, - "scrolled": true - }, - "outputs": [], - "source": [ - "min_x = openmc.XPlane(x0=-0.5, boundary_type='reflective')\n", - "max_x = openmc.XPlane(x0=0.5, boundary_type='reflective')\n", - "min_y = openmc.YPlane(y0=-0.5, boundary_type='reflective')\n", - "max_y = openmc.YPlane(y0=0.5, boundary_type='reflective')\n", - "min_z = openmc.ZPlane(z0=-0.5, boundary_type='reflective')\n", - "max_z = openmc.ZPlane(z0=0.5, boundary_type='reflective')\n", - "box = openmc.Cell(region=+min_x & -max_x & +min_y & -max_y & +min_z & -max_z)" - ] - }, - { - "cell_type": "markdown", + "execution_count": 10, "metadata": {}, - "source": [ - "Our last step is to actually create a lattice containing TRISO particles which can be done with `model.create_triso_lattice()` function. This function requires that we give it a list of TRISO particles, the lower-left coordinates of the lattice, the pitch of each lattice cell, the overall shape of the lattice (number of cells in each direction), and a background material." - ] - }, - { - "cell_type": "code", - "execution_count": 9, - "metadata": { - "collapsed": false - }, "outputs": [], "source": [ + "box = openmc.Cell(region=region)\n", "lower_left, upper_right = box.region.bounding_box\n", "shape = (3, 3, 3)\n", "pitch = (upper_right - lower_left)/shape\n", @@ -270,7 +277,7 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 11, "metadata": { "collapsed": true }, @@ -288,14 +295,14 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAGQBAMAAABykSv/AAAABGdBTUEAALGPC/xhBQAAACBjSFJN\nAAB6JgAAgIQAAPoAAACA6AAAdTAAAOpgAAA6mAAAF3CculE8AAAALVBMVEX///803tEsYIdnEy2T\nUVA4vPLpgJFyEhJNv8QsP88Otf3sYrFLC4epQXVfq1V9dXMyAAAAAWJLR0QAiAUdSAAAAAd0SU1F\nB+EEAxMJO4hQJO0AACJwSURBVHja7Z07dhpLEEC9BZ/DCmx2IDZguQmEhCJGOxDaAWgHsAAS5URS\nSqaUUCGpXkSsNTympz9V1VX9mRkk+53XmZFsz6W6/tU93779V9b37utH5/XzovP6H+R/kP9B/oMg\nI6XUfwFE6fX7rwdRqjqtGSH5+0BGJ475vCb5u0FOHHfzelWz3381iBrfzR+Wy+X8vlJ/N0g1X+o1\nv0ci+RSQ4cnI9AQyGt89LA0JEkl/IMPLGEecJB/ECWS5XNzPzgCiIg/b/OyyHxAnkOUS7a2eQLRt\nn/Ekw/qHs6hIskGanfW4Xq/13uodRI0rvdiv3Tqwy15A6p114livTnsLKklbEKUAyNDY9nuWxFD2\nA1KrSC2QWiSL+3wQBR8Yf/7Lg9QctfadLGK4gYYWRPUAoqqHRiANyCwXZIAeGH/uCYfju7mziJch\nSLPtZj2DrGptzwVRMRD3A28ST9/RmUHu2oOoFMjA+6haJB1BmFTjjCAKg9QCedQWkRFJGQiTaXya\nRJSzJLVIGJASZdf/43l0JA1yBy1JuLdi5nfYfOpARioukhG2WtnmV+VYrfGdNe3NPx6ARByiCV3y\nQcYZfoR55EGGHxkAkBUH0uwtdmcNKYhKgtx5zz5nPTvvMgasQL5Dz651fQ32Lfe0Mz6itJFmCciD\nJllpffzNgEjqkFwQZM2BRKLfACS1tUgYz4AMlKQPxSCX7PNyHN1A+MSq+dbaiISAzBmQH1I2MuRB\nZI56b82dQM4KwklEXIGyp/yIFomOIxakjAJAqp5A8jn03vqBQEYJjpNIqvnpP5zPsUAAyLjO4jqC\nrEpBfjS27Cd4ojhHXQ+qTO7zmwdRLUG+p/xIxvp5kV5WViMX8lwIIG3tbypE6QfEm4GRqQFc9A5S\nmXpA420LdL0EBLhK1rRF/EiuZxnYgsBqWWi0CkBGtRLHDIEMomO9HBKYWLVTkQwQZQK2JEgQohiV\nmmWQYG+b9eSkiJoDwjVFOBBSMBmMmwrP6e8mQYZOJIvMnUWLqEkQY6qyQOi3rEOBmiQpkiHwtlkC\nGdIwMhekKgYZ2JCmfrYUyA/tox60/LIEomiGcj4Qv+1PoUBKJKcnq7thJ2c7y9MQu08uW4KMspuh\nqMKTATIMHi2+sypSjigDYU2xANIIRHuGRXJvcbs+BdK0M1spO7/NeBDnrFfaEv1KgjQ5UhaG3lm1\nRQSFlSQINL/+r/9Ogoxd+LTK2FuZzw9B5uUg1iGOlP7b2jOkQAaowjOvzgNyXwLi9cI4Bm1PfydB\nfIpRK8kfAOKCxjqxihUfRBCdYvwBIDbkldq6PIgvTJ4RpEhHXGIl1YO+CKTYarnVCORxTdu6XwFS\n+5G7Mj/i5aJrpkxb9xNAhtRTtvDsAERo6wrK3idI6POD+nwBiNTWzTK/IkiTjGVwYK9fHv0CEKFB\nkuEQZT9iChYpkHEw91Ccj7iFWlawP90lRKlrYXVWn9CQcRgXF2eIPMgyCaLsxMdaGzkeZNBEoYlY\n0Wg2+a28nH0Uln1KQSrYqBVUZGCNz2UMRLEgZP2UOUgGJY4+ZCRWkq7XAsFpRb8gI1WFBYeRqwbm\nKDsZZpBUpKpMNB0FsWWlUpCR//fhp+MS81t/20YkcvFhYEaCEnurNUg9hTk//f+4DyJ2Q7+JT6kr\nPAu59jCIzjZ13lojWGNCIHw3VKxrNRFq/ZjCbwz8b5wDRMGqHxCJ1vbc6Pd7Y1vvTg9ZiSVTCzKP\ng7B+JAdE6HpKn4uTD85LfO8K0hTDC0FGuDLO/WCRB2Lq83K8mLm1uFgrB6RCvQq0txrdIRPMsVmU\neOMnU9lNOHJZCEK6R8huaZIFmSlvP2UKzG/0IX+kBcKBkH4eYrxjuqEdxmWbjDWj4Jsu3YUgtMOK\nDbNxTRe9gVR398moMWcFINrvgZ43ema2e9UBZGBjoctzgJDmPd53TDe0EAQZssbX9SCQOMg6AKlR\nfl90ATGTV04kTQf1C0DCHnURiAs3oHxSqe55QEYmg24FMmgCaxC2NB3tM4OsGJBwjqMExB7Tuied\n3i8ACXskBSCDsQ3tSaf3DCAxP3JheyQLOOtUAOKCnAURyTlAIp79okm6yPRZPsgAdHqrfkBc7MKA\niLGWEQidBywA8YE1yeMznviS+RA0UEMQF/2uafR7wfdI8kFARrMsA+HDX9jSZkBQPkJ0neuRFIOs\nw72V5BhzmZU1PHXMCUBGxg6Jc7EX7uhb8yjlICbrXwe99xQHG5INx8bw1CQ/PYYtyTnTQnJ2Jy3X\nI1GFIO5vr2j1McoxFMY7YJHk8qfjaOLa3xe+6bkIBIJKQtY0F4Ksi0EUn9rj0aifkMNkGk1daxkc\nkxaKdNkgA/c1rEtA3AgUORGHVeCndQ+VHwXgK40MyOxzQOC8FARBRum3eUA9m1wLoZk5YWq/zjQ/\nkvClBKTN1kI6ewn5QE9zZg3VHaqPmFlkytGArNuCtFR2uIUgiPtSfU9zhEIHZ8MCjo4S8c2fEvM7\nFI72KaYVOEKhQ/j8PelIM/p0+vsrySEOGQfebCFXMEQq8kh2Bwwd5gUgpeYX/z+/ApBhs6MvCQg6\no+g+R1G6iaWIx42BdPIjiaCRnwYcgson1HYAsjbf/4hYkyhIF8/uBk/5ML7xAdSBK1zCRSBrCgI2\nDFdu8KtbrAUSKwZk6Bw4mphtB7IKI3cM0in6dSJhU13V+DJ9XUlLEOQdEtreJR8xVetF2MWCDnyJ\nRVKiI2CSJwXSKUO0Veua4xcFERx4idUqAOmWs38HJxi+MyCc3yvxIyUgTRVl2bKKomuk3JkS7MDR\nEZICz16gIx3rWt9ttZd+6r545my+umNFxcRaBVara6WxQQk/s6qwXodn84X7IJjol/qRKEh4Yqyf\nC178F0xB/N7CKeI4yEdGJOG5uEiRtK/GiyDywTeQIVaJDJGE1xdl69wgP8BICfwY5ewGBEalpLb4\nB4DY63jSVZTsfORrQGxdq+LqWubw8A8LQjLErwCZy4dDhSur4GC5BrlA5f6vARlWkvltSLiRAXdf\nyqUtmZqZBltF+RIQG1KRINc9s3AFx7hC1XhU1/oiEOQWLn9kLZ/kGxBdyZpHz02eGQTW23JPucJl\nQEzk0YajNxD5uO4wY6LAgqRv4cgFGbQ51P6D1qTp9kkOefz0D5S68yEPxM4ylIMIx3XtIGDxmFM3\nkIEt6XCPG+C5D+pHMVeM0mvo9JAz4w3PDNKkKyxJuOcG7s/NV8/N0zl7irLGcLf1DAKmFH9xHMFN\nCr88CBuJDE1tZQEjxqHx9JfnA4G3FNJtFE5qEpAfzNDTEM5zw91GifsFGcCuXQAypsk6BWEiES49\nBN77XCB+9jmYIx9Y4xMDqTcN2vlc5YHdbX2D+O9vzoLgaWas7Nwa4sPkRiDMbusZJHJmxPXKRfPL\nLXKY/FLabT2DkOtN0iBuiSBuwscF90Kdq2eQOaxN/eoMMvR10ZUFIWeezgjyyN4q0BJkjg4BXloh\nPbra1eVZQHAHoCeQtQdRvJDOB8Ic2YufwcgAMRdvcUL6AhDpHqECEFRx93dxfdrW4jx7W5BH/FFv\nIE92bbCyXz3BZWKtJ3Y982tLdeT5WVEQ86svr50XBLnzfvh+hp91ozmuykGAQkxOIHcQZHk2EOgQ\nnxiSp0IQfCnJZ4E84UpI8LiSPESQZ3IGOwA509Z6QmGQ+NQFIMT7nUC2n6EjYG/NewLB8chzA7L8\nBBCQWBVwiCBbHCEyhmxyFpCnOtW1lZC2IFulJl4kd2CzTgJDptXmHCCbsauElOysJ4RRL7C3HqBA\nOP0/B8iTPWqvigTyBDmakMztLTBiivW/8bqe/rpXkI2NDYs4PIg+0Ff/AxP3wZ0tMDUfbfFcOpBi\nVxIEYhx42cbyIFtfpJu4T+wUV6D/Rm2erRR7BYk68CQIMBYTz4Y22zOaJJ3436i6iuRb6VPLIFs4\nq+weHKu//yWrNraVOOsokh5B0KjTxO83yGHFVvfdG4H4ybtuIukXBLs/3ksqM/RlBeJ91x8CAv34\nQgZ5dvk/EmO9Hf8YEDQsN4mQqLHdblt3NV7XvdUfiMI5vzPAWEXIZ1s/C3n/54A8gOkHA2KMliie\nLew9/ikgKIyqgLNz7jAB0klJzgFiSw2hYw/1xZvsjtqeB6KiDr95prFPYy2I9xmS7n8yiDk7flUI\nwrr6rwRx476FWwsVAXiRfKqyK39heR7Iqsliw0z3a0E2yp9Fu8oFqc3vdgw6b/ciyKf5Ed/1uxdq\nEtSPLBsQ9IGwt7YwRDk3CLyw/EoG2YJXLGidEPJzqu137rjG6+uufa6YBNmgC8ujIKjkSiqmUhip\n7IxmnZDsOmS9SRBUZeP3ln0mrNtb2lXglcSnx9c7O9PThiQNAuqeQZEeg2Bri6uKy3vJk7irlF79\nfYlnANmQrt+VDLKFL6qcoDqvvqPlWRKJydl3QDjnAMFFzgjIsz9nVcdWmSDPuhl2ihtegbrkgOyw\nMiVBcNmZVRKkuK4cBFsIHkSF+Umj4a+v8LxZDkbzt3JBgmYZ8ztb+3R+vz8LIKp5XVZIcno0kCwm\n99YuqIYlQWiPiRGaObBU7/exTz+4raVomcusl/rZYFcjBRLahRKQNQuycTGlK39OBJCtvtX7HtRU\nAQiqJiVAdn4TX7cECe2vv85+gtPx0PyCoTMWBLihhEBCdXIgQuqUBNk0tkbfNojtauAQ68BKJ1r0\nV0tBdkyo+c3tc756ndxa8P7HCQahIcoWHIOcdAHhfvUb3OezkCTowwYCgadXsDGiQSM6H9IBZIej\nOgjizydSkk3C/G7QQWb0PQdhvD9hS5lLQdCMxrUH0fv8oTkmxIDEPDu4tZaG6iaxegThsCS8MhB8\nfgiCwOyGPOmGtC+DH+PD/kRJSDgcNHlagpDqgAPZgMQ5iEHgMcfl/YzjBKcgsZL4XqgOh3NBFvN4\n/IvrNfceBB3YpyC4fRnqOj6XSgxwZZs6pjpEQbbG67wgFU45RHKG3IPgV3rQ7xwdc4yo0CpMzHUv\nVL9xTDvLMbULLkh7EZyDrOvA084siN8e+llZLbD/Pavra69CGASXTEMQF9BokNxWCQRZIxByp7zk\n8biZIehmmAoD6oWGW8s1GF+0Dmc2rwKQawNCX94eqIHv+l8VgqBRiABka45yGxAmEiyTiEIXSXDb\nx3X9n0pB9Lc+CczxvJkVMvebTcwIB0jfW4HQqz0oyEYH35zbzwHB5hg5RGVmVy3ILq/lviOXGluQ\nDQZchPsnMtkRV/YQxMZaBqRCIK8ufY8vej0oB7Ji6wsbccY0bn4pCBgWmHAgrzQR55fiHWIGiM5W\nnriFHOJCrJR4B2nL4RNma70GpREZBNT9i0CktSE3REVBnq3ZsFk9UfbctcM3TFz3AfIUCRoZkSib\nLUyMbYbmNx8EX17CgrDKHluRMJ7VEtjl3WKHmL9wlfnVgCTNb2JviYkVKxJlNpP1MT5EKQGBVeZr\nB5JyiCklkbKl+kEnIQmYH0BBY9neCmp56RAlubeE4oP/vumnmGtSDvIKj36+OpBE0JgUiVAO8l4i\ntpqfloIwhftkGJ8WCSzQYY7o0ANcpSCv7pys5YCJ1YpNrNIicRcrMZaWHXrYKq5kWkhi7ih5hSDR\nVDeHxLoDZGj9u1bJUzO6Uw4SRDO0+FA2TN6src853JcO37VKbUAFTXBbkBPKNfxTqhyUs5hNf0pr\nH3irvHXv2ewKgleqQNcSBL9taoI4jL2BnxqQXYcx5lTJtB3I1oZgKxq5mMR8ia31i9/218Ej5nXe\nU0XstiD4ijkoKs5/vpgHrm0GeehdXo6SbCu0A6GX/hFRWX1EIDubhxIO7lMZRGz0tAQhUQ8UyJIx\nAjWIO+pwjTgqezFILkiHxYOwNQl8p/cEgOz8peRQPzJLRDJIiXxCEHpVqX1iyQi8CFUtV31M994F\nkE2JwvAgMFfzIOA9xWBvvTQB7QOde3D3lqdb1gJIkeYXgJD3FCMQP//hQZhUsBwke2/lg0jW7EXv\nrKCEvcNXuv4lIGC01u8sdDL9+mtAOGWnL1xmQOAkCukLE5CsoZqOIJL5lfzLC5xE8SC0Ux9wKJxY\nCSBXvYCgSmpUIg7EfvU7OjuBOXCd+DzmV9KFqI4EEuF7bNbDjHGuKzvEbI6CoFE2v0hH7MPS+SIo\nkLEtCSRAShYHgiqpkxTgix87W2SCVBWOA3oA2TAndoTEaosuXvaAL2wwQkCg2dopd+pP9QWyYctw\nwlE+3AkHIFzNLQLiG6dVJsgmaYhVWHqAo7M401X8DM5L82hNT949b7C1AEjlizRZINbIATCWlBEJ\nWw7aopcJPAMQUDx8lUCgitCzAXGQjW1PXnmyECSsaqECHfyRUCZ6eeV6urvxg2R+YcacAbLxtuHK\ngVzh37Hvhg3UfcyWTHlJaRB7Qu7ag0gOcee2aB4InEWQJCLWqvkfbIG5eUYgYRV0h0OU9iAbL8DK\nloQ32SBSW2GM2m8AJKyXkKON161Bckqpke7BNrRmbsi5YkqmO1r3IcOErUGwMygHATykuZNZxHbD\ni813iRD9XkmCWIGslrFLeJIgKnA0DQeUlwgipbrQZmSAoKvHRZBxUFwPdT7s9yA8AYTtekLPvszy\n7Bvyag3JIrB+BFgpeKEAVhX31yQQNgFufuCvz0+C0ORBBIkd9CZXPAgqL4HwRaJXHwf4z2MgJM0T\ntT0C4o/r4l4pcZciyM6/mYEQOiueBqEvBJJ2IBc0eoGAQ4well4FI4Igf78DWbqb7cpIrMaZgwQb\nmQPNMEMQ8ObTOAhoNuyM678GhLOc8yO0piOBiPdr+avCtAUHe4u2eyIg3k3iljQ9fJkAyZpRiYDg\nQ4wTAOii8MlzXg8RnKLxgFlnrHK3lgyyJYcYAxCrORkgNhdeeEOF4pnuyh4FCW5qdCqCspIMEFeY\nF+abu5vfGAg9xEiNcj6ILxcJo/MJkLRDTIE8epB7AALC6kwQNKKfAbJxxQaFpoaqFiDSRZktQJr6\n0tpelZgGUc3FPrViZ08NZYHA47poa+WB+FP+QochuPGscg13E8av+eMYvYFQHQnSK6MiyP6lQDa+\ntfqUPzVUDgITPAJCHR1QEfcPLaoUCHgD9xU6OhK9S7MY5BkmeBjE3TJAQOjdoXEQcgnoZgy+uKuO\nIFDZg6LQi+dwR2cCXX/MB0Hnfq8aAdnC8lMLEMn82psoF0GsJQ0IBF9JAoTc+Lvx0fZVKxBQKsR9\nEjIoaEFATR6LJNikURB8mbRW93oU5D5ybVAKBB1i9G0Ed6khyUcU10ksByGXr1092RpJoi8qgjz7\nQ4wrkpCYvIhkiHKxoQwEF/b0bhJPjmSBcHewmp+gKtKL8xZS+adIR2g32OAl29QyyBYfYnz2l1fg\nMp0DQdcmtbVaG1r8TgFkgJCDbmDIH42mvvjHBXsb7K0iP7Kh3eDuIPR8mDt2MdFCcb9mQcRJhyLP\nvgHGcp0LUseXERB8B2v0zodXGBmucPW9LNba0E5XBkcd9MdAnsEhxsnWXznGXJXwGh3ZUAXRbwBy\nlSOQk0hiIO64rqpvGADZKgOyi4FY15CRj7SQyCYJAu502Y7ppeVxEPS4cHQrlSG20JEMEN9W2MJs\ntRhE5efsLaxWtFwKUSbPdhpiba8cK9tabBVlp9jpoBZ+JEcibsFslZy9TCt7M7RB61ogBcPvVsCS\nFYxUkbIDtZeqdc+M+WWcBa00kuovH2v5m0M2MIAP585VyvxCEP69CgAkMfo3ruCTN50J11jA0a+v\n1JpaL3ryjTlbi0muckGEah0Ekbq4ViO0+bu2W814qEZCOB+hb+kxNX3IwSQneSDSS1QAyA695yMM\n1aFOgLNPmvjbm197HEW/vTkfoH/sjrfNbt7Qes9aAQj42T8feh3xA3xE1hHd33P6AIK8obzm5m3v\n7vWqn3yvZbnUXroriL4XKgT5QJe83EZBYFytbgnIHg4W3rwBOwF+WMvyJgQ5aLl3B7Gp7n0cBL6V\nXv8qAnmDA+r1n1zA1/zswZU5AxA72NUN5OgGy+/jO+uDHD+nIHtfjbnRf3hwItij1IKCHFxvsouO\nnL5oXSHUzcXEzsJ3FlGQN32sxmjFHuYS0BJoBUIgB183mnYCOfp2dFQgx2qORroDkDc7NIX2Ui2C\nPTq4TkDgrKQEcsBpd8WCfLh2dFQgH2ike1kbBgry1nhvbMNqEKxdGOQAi9KSSA7Es09ZkA/zerw4\nxwcdvmdATgA31BhTkDkFQWfZYiAg1hJAPvQXmeAgIIsKgzhZNIrvM4DT93yj0Nt1MAjuP0h7i/Rq\nJJCsBW9FWFEQczu8fcY9Suhummzi0b4mCIIccEdI2luke/beBQTOAFAQp2YtQFCPLqokbrNOzwWy\nt4NVViT81loyW+tA7loQQaAB7wYS2VrYjRcpOzqnE1ESZM/fO4JIyo7dOAOyr0Q/QicLRE8C/g8n\nkEMdDHQEgeYX5fYOBDlE0bPTMaKpKBIfhQKOE0k5iOgQ0dd/g7Wdhig01qJjRCLIQaendUjmfkWb\n/GkxyJFcq+hCFO77JrsNZTLtQBqSkzymSCBKtQDBt0gBEEaVcSZlO7HahrUEMQ8+JX9OOnJeSer/\ndUXCeOy3rZKQxKoR0NztvRYg9WZCf2gmRlqAINtzC0CA17UgONV9AydC3gKQLGVnqGo33AbkiFT2\nw4PgG0uxszcB2F7ZQhIWSGB+i0DSsS7/YwVm10SQmQ9bxiY5MVsNcAkOMeJHWknkqJWIFYk1Rfqv\nx0GagNgKYI+5/C+Rke1CkFkM5Fj/nywJ2OkfaZDTw9/AP6BA3y6faayWJPVLraTVsrv7lkG0V958\niCDhs/pvn/nsYI9JaltYoiJJP0KeljCOQVIMQZYZINx6x3drxEEO+McJEF/dvOVIgDSRH3Ej01Y/\nbjJB0NGfhE7gMt4hGqIcQaWa33iODztEEErtFa8QHMgBjvZNoxxB8Us7ehEEOGFZjyAIql/XEcje\nTrLlgLzDQbIERzNsRH4ggcCwqBikCRr9Oc0ckIM/XRgRiFzGk0FAmTYLhGRRLswKCtY8yLvnjgkE\nvLMqG8RH3JkgOPmAxwlyQN5tQy8mEHgN+TQH5Ihus8wD2YM3utUNBTAqnRQJsEdRTZdtmwzyUAqC\nY3b8f2aBaEP6HlsRb9MnCIzZ9/i9GHkgqRXz/32CwJh9j4cbU3srE6SSI7IedeTN1uF1VovnpvoB\nicXI/YJoZ64fmo5A9AMSqX4JIDB7KgFxQPQ+815A6JGtaR7IQ5lnJyDkqvxMkKj9PcRKFCII6CS3\nAnnAIDdZIHFPcsDfTh7IEVzQ/tEChF6skgWiaPEtAnLasHkgJJ/9DJCDGxLvEeTDnC2+T0bxHUB8\nAt8IxB1nnqZBlrkgR/f1fJwLZG9uFjAgB7AJ+tMRW3xIdas/ZGVPgezdV2VAwB3fsvktBzGJeVIg\nrc2vOVWm85VmZ3lDORVBpHpkBOQjjyMJUv+fivuVO1/R1iDAdfEgyLMvYp69VX2eA8kIUWAOakF8\nKZYHyYy16hrpLEsIaRDygo9QRXD76h1mf2L9wTexV9H2tK0wJ2RCVYcHwS/4uOF/AeZIB0rGgchj\nBf8gDmPGb+MczfBPHGSPX/ARgthnWtsHzwCJDXp4kKOtx8TLJu7ahNsoyBs+xyeA+g5eHkhk9OYf\n8IQ5YeLRF5Zu4yDwZCWjIv66RZNaZIH4YahKqqIcaf9GEIgrLLlf4kH28Kwrx/mAutxToOwEBAYs\n8N6xdwEkK0s/giFSKzce5M38n0KBjs4dTLVjZ7bNgbRw3U1wUwEk6HEKAvFTtRZXANn78+CsCuGR\nlin2I/4pg2a0UVHCAUDAPY7i5O8RF3miIDaWmrGdBTpkpEHuQs9+oF2EQ9MyDJToH/SIj667LIH4\n0fellZsEohuGUoeEgNRmS4He5BQJhGmHvMsg+NJWuJ38H/xZAO1aEyBNvsH/iAE5KPueCaAiptUZ\nPngmyAzYW+DG6ctPEyCRxYC8+3AYqojuoheAkFfHAXkAElrk6Q1Ex00H804eqMkmjsgHUQTkFgrE\nhsRHehasPQhjtXSuS/pqh0wQ/4wSiJ28D0DW7rdagjyEwycHcxkdesIqvbUOYCRAADnaV14YEJLA\ndgAhr3OceouEnjpHRw5OkkkQk7kfq95AwljLkdA9M050f95R2C4ou62l9A8SRL/irkmcjgFmuw7/\nBD+iOoDsVWxcIMhHYiBxgVQwbD/CLFIEKdIRNBbEKAnJECP2KC4QlCLjEMWHvxSkwGr5po+0t5Lp\nhyWJqzp431UtAOWvTwVB41EESfmRvfEAMZG4KkrJWFO4s3DYzofxGKTEs++t4Y6IBNe1Wi68lW6F\nxAqb35JYi55c4ERyhyqN7dYhGN+FeY0ViFUSG2wVRL+uSCHuLVr7bQtCzS18ax8EASFKST6ivdR9\nDMRV4286gQQOkG/toKCxIEPcu9dmyyC0P9JqhbP6H8pF0V4gtK6Yl7OrG3nCl18dQILTE1JHBNVH\ns6oopy/6a0HwqCK/cupaJ5CbrK11NhCtEVGOrEpjDaKtVkLZzwiSsdwb7tzfkEHirr0/EKLsuSTp\narxSOQ6xJxAhtc0iyeiPGCcxy5uW7eRHhPJPzjqiP4nR79h4O3ntbaTfAsBExQe5IFe6vkUeMi4Q\nc2KiTYjiMscDn3/0CPKW5LDBlioPGkGZiM8/+gTRjjHOYS15cRjvTEmYj5wBJLHkpk2GQGzcMCWd\nrg4cbUH2Y7GNlhYIPB9aMCMXrrT5zRCIazLOS0XSHG9cNI0re353nhxaDNYRhWRtQVBGkKPerghx\nUOjosWvHlQoE5ygtQWAjPuO8Hmr4HOAUPhxYK95XMP1tCwKHkeU7BZwM4JCgO05g2qb2lTLFAsEj\nA+1A6JXySe0GRxnQgXO774pnaD78QY+uILCyFBeInUtsrAK46MjW9g7lhjfISNqB0CEyrNaEC1zC\nIICUH8xvBmkbI94BZC/OJ6pw+JdekNAPCHTIZwA5KKgNVkOQleoLBBxz6gLywM7wHvy8BNxZ+O6N\nQNlbgphbwOddQR4ZEHB//hSCoMI9Nb+tQI7wEFjvIAc4U+RVhFzrQhziF4OwOoKmvKYiiApvR+kd\n5KYLyAFfTiiB4KCxueGlvY6wIPG+IQK5w80vAwLvGJ/JIOial0P09HQMRLRa+1h6ixNG1iGSCzyn\nEghMrIznuW0B4l/DS0FUpAFKGNlYixwxE0FgqmvrEOUg3kDe8iCsSPaEEfWM7FE8Or8tgsAZFZW6\n9EHSdn8RIwXZm1kSTiRmhAFru3/no1cRb5IjIHCQwPSiykGs7w37I6AFxQgEV4Lx0VEOxEeSyCG6\nEMDEY4eghpu/t2qSignjbZmUBaFNXuYA8YFeci6EKGgd0NxPqUhgQygLJKxp7+Hx6ikHModxC/Xk\nPYCQBmkmSNhlgA7pPQrC3nPmQVRbkCO2d21B9uAKgjgITawYkDY6QgY2qdXilZ1ru4NLIVgQH20p\nnOrilXfjZ8ZCIGOhJcLOD9hN4Q7m0/OkRA/44cbasLfxIxGQiEPkO3Hu1y0Ib36ZC9vw3mrl2WMg\n1H3HrJb9C/q3ORBUWgknssnPugskcilrCCJFlOaR6AsF6PfOcrSNfqMgb1Lwm2gp2u8W38M/fc9b\n7fKROIi44i1FZ2ZNALbKKgjD9YkgsVzFfrXjshL9V4DEZzWd3mbfTvV1INEs2IH4N7WV7axPBIku\nr7bu3XllAvnTQFxZYV5lHLX4k0Hybqf6G0BcelHG8eeB5Iz1/x0g6SHyvwak1fof5D8I8i8rFUXT\nyFgUuwAAACV0RVh0ZGF0ZTpjcmVhdGUAMjAxNy0wNC0wM1QxNDowOTo1OS0wNTowMEJcIC0AAAAl\ndEVYdGRhdGU6bW9kaWZ5ADIwMTctMDQtMDNUMTQ6MDk6NTktMDU6MDAzAZiRAAAAAElFTkSuQmCC\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAGQBAMAAABykSv/AAAABGdBTUEAALGPC/xhBQAAACBjSFJNAAB6JgAAgIQAAPoAAACA6AAAdTAAAOpgAAA6mAAAF3CculE8AAAALVBMVEX///803tEsYIdnEy2TUVA4vPLpgJFyEhJNv8QsP88Otf3sYrFLC4epQXVfq1V9dXMyAAAAAWJLR0QAiAUdSAAAAAd0SU1FB+ILDgkoDSxqv5MAACJwSURBVHja7Z07dhpLEEC9BZ/DCmx2IDZguQmEhCJGOxDaAWgHsAAS5URSSqaUUCGpXkSsNTympz9V1VX9mRkk+53XmZFsz6W6/tU93779V9b37utH5/XzovP6H+R/kP9B/oMgI6XUfwFE6fX7rwdRqjqtGSH5+0BGJ475vCb5u0FOHHfzelWz3381iBrfzR+Wy+X8vlJ/N0g1X+o1v0ci+RSQ4cnI9AQyGt89LA0JEkl/IMPLGEecJB/ECWS5XNzPzgCiIg/b/OyyHxAnkOUS7a2eQLRtn/Ekw/qHs6hIskGanfW4Xq/13uodRI0rvdiv3Tqwy15A6p114livTnsLKklbEKUAyNDY9nuWxFD2A1KrSC2QWiSL+3wQBR8Yf/7Lg9QctfadLGK4gYYWRPUAoqqHRiANyCwXZIAeGH/uCYfju7mziJchSLPtZj2DrGptzwVRMRD3A28ST9/RmUHu2oOoFMjA+6haJB1BmFTjjCAKg9QCedQWkRFJGQiTaXyaRJSzJLVIGJASZdf/43l0JA1yBy1JuLdi5nfYfOpARioukhG2WtnmV+VYrfGdNe3NPx6ARByiCV3yQcYZfoR55EGGHxkAkBUH0uwtdmcNKYhKgtx5zz5nPTvvMgasQL5Dz651fQ32Lfe0Mz6itJFmCciDJllpffzNgEjqkFwQZM2BRKLfACS1tUgYz4AMlKQPxSCX7PNyHN1A+MSq+dbaiISAzBmQH1I2MuRBZI56b82dQM4KwklEXIGyp/yIFomOIxakjAJAqp5A8jn03vqBQEYJjpNIqvnpP5zPsUAAyLjO4jqCrEpBfjS27Cd4ojhHXQ+qTO7zmwdRLUG+p/xIxvp5kV5WViMX8lwIIG3tbypE6QfEm4GRqQFc9A5SmXpA420LdL0EBLhK1rRF/EiuZxnYgsBqWWi0CkBGtRLHDIEMomO9HBKYWLVTkQwQZQK2JEgQohiVmmWQYG+b9eSkiJoDwjVFOBBSMBmMmwrP6e8mQYZOJIvMnUWLqEkQY6qyQOi3rEOBmiQpkiHwtlkCGdIwMhekKgYZ2JCmfrYUyA/tox60/LIEomiGcj4Qv+1PoUBKJKcnq7thJ2c7y9MQu08uW4KMspuhqMKTATIMHi2+sypSjigDYU2xANIIRHuGRXJvcbs+BdK0M1spO7/NeBDnrFfaEv1KgjQ5UhaG3lm1RQSFlSQINL/+r/9Ogoxd+LTK2FuZzw9B5uUg1iGOlP7b2jOkQAaowjOvzgNyXwLi9cI4Bm1PfydBfIpRK8kfAOKCxjqxihUfRBCdYvwBIDbkldq6PIgvTJ4RpEhHXGIl1YO+CKTYarnVCORxTdu6XwFS+5G7Mj/i5aJrpkxb9xNAhtRTtvDsAERo6wrK3idI6POD+nwBiNTWzTK/IkiTjGVwYK9fHv0CEKFBkuEQZT9iChYpkHEw91Ccj7iFWlawP90lRKlrYXVWn9CQcRgXF2eIPMgyCaLsxMdaGzkeZNBEoYlY0Wg2+a28nH0Uln1KQSrYqBVUZGCNz2UMRLEgZP2UOUgGJY4+ZCRWkq7XAsFpRb8gI1WFBYeRqwbmKDsZZpBUpKpMNB0FsWWlUpCR//fhp+MS81t/20YkcvFhYEaCEnurNUg9hTk//f+4DyJ2Q7+JT6krPAu59jCIzjZ13lojWGNCIHw3VKxrNRFq/ZjCbwz8b5wDRMGqHxCJ1vbc6Pd7Y1vvTg9ZiSVTCzKPg7B+JAdE6HpKn4uTD85LfO8K0hTDC0FGuDLO/WCRB2Lq83K8mLm1uFgrB6RCvQq0txrdIRPMsVmUeOMnU9lNOHJZCEK6R8huaZIFmSlvP2UKzG/0IX+kBcKBkH4eYrxjuqEdxmWbjDWj4Jsu3YUgtMOKDbNxTRe9gVR398moMWcFINrvgZ43ema2e9UBZGBjoctzgJDmPd53TDe0EAQZssbX9SCQOMg6AKlRfl90ATGTV04kTQf1C0DCHnURiAs3oHxSqe55QEYmg24FMmgCaxC2NB3tM4OsGJBwjqMExB7Tuied3i8ACXskBSCDsQ3tSaf3DCAxP3JheyQLOOtUAOKCnAURyTlAIp79okm6yPRZPsgAdHqrfkBc7MKAiLGWEQidBywA8YE1yeMznviS+RA0UEMQF/2uafR7wfdI8kFARrMsA+HDX9jSZkBQPkJ0neuRFIOsw72V5BhzmZU1PHXMCUBGxg6Jc7EX7uhb8yjlICbrXwe99xQHG5INx8bw1CQ/PYYtyTnTQnJ2Jy3XI1GFIO5vr2j1McoxFMY7YJHk8qfjaOLa3xe+6bkIBIJKQtY0F4Ksi0EUn9rj0aifkMNkGk1daxkckxaKdNkgA/c1rEtA3AgUORGHVeCndQ+VHwXgK40MyOxzQOC8FARBRum3eUA9m1wLoZk5YWq/zjQ/kvClBKTN1kI6ewn5QE9zZg3VHaqPmFlkytGArNuCtFR2uIUgiPtSfU9zhEIHZ8MCjo4S8c2fEvM7FI72KaYVOEKhQ/j8PelIM/p0+vsrySEOGQfebCFXMEQq8kh2Bwwd5gUgpeYX/z+/ApBhs6MvCQg6o+g+R1G6iaWIx42BdPIjiaCRnwYcgson1HYAsjbf/4hYkyhIF8/uBk/5ML7xAdSBK1zCRSBrCgI2DFdu8KtbrAUSKwZk6Bw4mphtB7IKI3cM0in6dSJhU13V+DJ9XUlLEOQdEtreJR8xVetF2MWCDnyJRVKiI2CSJwXSKUO0Veua4xcFERx4idUqAOmWs38HJxi+MyCc3yvxIyUgTRVl2bKKomuk3JkS7MDREZICz16gIx3rWt9ttZd+6r545my+umNFxcRaBVara6WxQQk/s6qwXodn84X7IJjol/qRKEh4YqyfC178F0xB/N7CKeI4yEdGJOG5uEiRtK/GiyDywTeQIVaJDJGE1xdl69wgP8BICfwY5ewGBEalpLb4B4DY63jSVZTsfORrQGxdq+LqWubw8A8LQjLErwCZy4dDhSur4GC5BrlA5f6vARlWkvltSLiRAXdfyqUtmZqZBltF+RIQG1KRINc9s3AFx7hC1XhU1/oiEOQWLn9kLZ/kGxBdyZpHz02eGQTW23JPucJlQEzk0YajNxD5uO4wY6LAgqRv4cgFGbQ51P6D1qTp9kkOefz0D5S68yEPxM4ylIMIx3XtIGDxmFM3kIEt6XCPG+C5D+pHMVeM0mvo9JAz4w3PDNKkKyxJuOcG7s/NV8/N0zl7irLGcLf1DAKmFH9xHMFNCr88CBuJDE1tZQEjxqHx9JfnA4G3FNJtFE5qEpAfzNDTEM5zw91GifsFGcCuXQAypsk6BWEiES49BN77XCB+9jmYIx9Y4xMDqTcN2vlc5YHdbX2D+O9vzoLgaWas7Nwa4sPkRiDMbusZJHJmxPXKRfPLLXKY/FLabT2DkOtN0iBuiSBuwscF90Kdq2eQOaxN/eoMMvR10ZUFIWeezgjyyN4q0BJkjg4BXlohPbra1eVZQHAHoCeQtQdRvJDOB8Ic2YufwcgAMRdvcUL6AhDpHqECEFRx93dxfdrW4jx7W5BH/FFvIE92bbCyXz3BZWKtJ3Y982tLdeT5WVEQ86svr50XBLnzfvh+hp91ozmuykGAQkxOIHcQZHk2EOgQnxiSp0IQfCnJZ4E84UpI8LiSPESQZ3IGOwA509Z6QmGQ+NQFIMT7nUC2n6EjYG/NewLB8chzA7L8BBCQWBVwiCBbHCEyhmxyFpCnOtW1lZC2IFulJl4kd2CzTgJDptXmHCCbsauElOysJ4RRL7C3HqBAOP0/B8iTPWqvigTyBDmakMztLTBiivW/8bqe/rpXkI2NDYs4PIg+0Ff/AxP3wZ0tMDUfbfFcOpBiVxIEYhx42cbyIFtfpJu4T+wUV6D/Rm2erRR7BYk68CQIMBYTz4Y22zOaJJ3436i6iuRb6VPLIFs4q+weHKu//yWrNraVOOsokh5B0KjTxO83yGHFVvfdG4H4ybtuIukXBLs/3ksqM/RlBeJ91x8CAv34QgZ5dvk/EmO9Hf8YEDQsN4mQqLHdblt3NV7XvdUfiMI5vzPAWEXIZ1s/C3n/54A8gOkHA2KMliieLew9/ikgKIyqgLNz7jAB0klJzgFiSw2hYw/1xZvsjtqeB6KiDr95prFPYy2I9xmS7n8yiDk7flUIwrr6rwRx476FWwsVAXiRfKqyK39heR7Iqsliw0z3a0E2yp9Fu8oFqc3vdgw6b/ciyKf5Ed/1uxdqEtSPLBsQ9IGwt7YwRDk3CLyw/EoG2YJXLGidEPJzqu137rjG6+uufa6YBNmgC8ujIKjkSiqmUhip7IxmnZDsOmS9SRBUZeP3ln0mrNtb2lXglcSnx9c7O9PThiQNAuqeQZEeg2Bri6uKy3vJk7irlF79fYlnANmQrt+VDLKFL6qcoDqvvqPlWRKJydl3QDjnAMFFzgjIsz9nVcdWmSDPuhl2ihtegbrkgOywMiVBcNmZVRKkuK4cBFsIHkSF+Umj4a+v8LxZDkbzt3JBgmYZ8ztb+3R+vz8LIKp5XVZIcno0kCwm99YuqIYlQWiPiRGaObBU7/exTz+4raVomcusl/rZYFcjBRLahRKQNQuycTGlK39OBJCtvtX7HtRUAQiqJiVAdn4TX7cECe2vv85+gtPx0PyCoTMWBLihhEBCdXIgQuqUBNk0tkbfNojtauAQ68BKJ1r0V0tBdkyo+c3tc756ndxa8P7HCQahIcoWHIOcdAHhfvUb3OezkCTowwYCgadXsDGiQSM6H9IBZIejOgjizydSkk3C/G7QQWb0PQdhvD9hS5lLQdCMxrUH0fv8oTkmxIDEPDu4tZaG6iaxegThsCS8MhB8fgiCwOyGPOmGtC+DH+PD/kRJSDgcNHlagpDqgAPZgMQ5iEHgMcfl/YzjBKcgsZL4XqgOh3NBFvN4/IvrNfceBB3YpyC4fRnqOj6XSgxwZZs6pjpEQbbG67wgFU45RHKG3IPgV3rQ7xwdc4yo0CpMzHUvVL9xTDvLMbULLkh7EZyDrOvA084siN8e+llZLbD/Pavra69CGASXTEMQF9BokNxWCQRZIxByp7zk8biZIehmmAoD6oWGW8s1GF+0Dmc2rwKQawNCX94eqIHv+l8VgqBRiABka45yGxAmEiyTiEIXSXDbx3X9n0pB9Lc+CczxvJkVMvebTcwIB0jfW4HQqz0oyEYH35zbzwHB5hg5RGVmVy3ILq/lviOXGluQDQZchPsnMtkRV/YQxMZaBqRCIK8ufY8vej0oB7Ji6wsbccY0bn4pCBgWmHAgrzQR55fiHWIGiM5WnriFHOJCrJR4B2nL4RNma70GpREZBNT9i0CktSE3REVBnq3ZsFk9UfbctcM3TFz3AfIUCRoZkSibLUyMbYbmNx8EX17CgrDKHluRMJ7VEtjl3WKHmL9wlfnVgCTNb2JviYkVKxJlNpP1MT5EKQGBVeZrB5JyiCklkbKl+kEnIQmYH0BBY9neCmp56RAlubeE4oP/vumnmGtSDvIKj36+OpBE0JgUiVAO8l4itpqfloIwhftkGJ8WCSzQYY7o0ANcpSCv7pys5YCJ1YpNrNIicRcrMZaWHXrYKq5kWkhi7ih5hSDRVDeHxLoDZGj9u1bJUzO6Uw4SRDO0+FA2TN6src853JcO37VKbUAFTXBbkBPKNfxTqhyUs5hNf0prH3irvHXv2ewKgleqQNcSBL9taoI4jL2BnxqQXYcx5lTJtB3I1oZgKxq5mMR8ia31i9/218Ej5nXeU0XstiD4ijkoKs5/vpgHrm0GeehdXo6SbCu0A6GX/hFRWX1EIDubhxIO7lMZRGz0tAQhUQ8UyJIxAjWIO+pwjTgqezFILkiHxYOwNQl8p/cEgOz8peRQPzJLRDJIiXxCEHpVqX1iyQi8CFUtV31M994FkE2JwvAgMFfzIOA9xWBvvTQB7QOde3D3lqdb1gJIkeYXgJD3FCMQP//hQZhUsBwke2/lg0jW7EXvrKCEvcNXuv4lIGC01u8sdDL9+mtAOGWnL1xmQOAkCukLE5CsoZqOIJL5lfzLC5xE8SC0Ux9wKJxYCSBXvYCgSmpUIg7EfvU7OjuBOXCd+DzmV9KFqI4EEuF7bNbDjHGuKzvEbI6CoFE2v0hH7MPS+SIokLEtCSRAShYHgiqpkxTgix87W2SCVBWOA3oA2TAndoTEaosuXvaAL2wwQkCg2dopd+pP9QWyYctwwlE+3AkHIFzNLQLiG6dVJsgmaYhVWHqAo7M401X8DM5L82hNT949b7C1AEjlizRZINbIATCWlBEJWw7aopcJPAMQUDx8lUCgitCzAXGQjW1PXnmyECSsaqECHfyRUCZ6eeV6urvxg2R+YcacAbLxtuHKgVzh37Hvhg3UfcyWTHlJaRB7Qu7ag0gOcee2aB4InEWQJCLWqvkfbIG5eUYgYRV0h0OU9iAbL8DKloQ32SBSW2GM2m8AJKyXkKON161Bckqpke7BNrRmbsi5YkqmO1r3IcOErUGwMygHATykuZNZxHbDi813iRD9XkmCWIGslrFLeJIgKnA0DQeUlwgipbrQZmSAoKvHRZBxUFwPdT7s9yA8AYTtekLPvszy7Bvyag3JIrB+BFgpeKEAVhX31yQQNgFufuCvz0+C0ORBBIkd9CZXPAgqL4HwRaJXHwf4z2MgJM0TtT0C4o/r4l4pcZciyM6/mYEQOiueBqEvBJJ2IBc0eoGAQ4well4FI4Igf78DWbqb7cpIrMaZgwQbmQPNMEMQ8ObTOAhoNuyM678GhLOc8yO0piOBiPdr+avCtAUHe4u2eyIg3k3iljQ9fJkAyZpRiYDgQ4wTAOii8MlzXg8RnKLxgFlnrHK3lgyyJYcYAxCrORkgNhdeeEOF4pnuyh4FCW5qdCqCspIMEFeYF+abu5vfGAg9xEiNcj6ILxcJo/MJkLRDTIE8epB7AALC6kwQNKKfAbJxxQaFpoaqFiDSRZktQJr60tpelZgGUc3FPrViZ08NZYHA47poa+WB+FP+QochuPGscg13E8av+eMYvYFQHQnSK6MiyP6lQDa+tfqUPzVUDgITPAJCHR1QEfcPLaoUCHgD9xU6OhK9S7MY5BkmeBjE3TJAQOjdoXEQcgnoZgy+uKuOIFDZg6LQi+dwR2cCXX/MB0Hnfq8aAdnC8lMLEMn82psoF0GsJQ0IBF9JAoTc+Lvx0fZVKxBQKsR9EjIoaEFATR6LJNikURB8mbRW93oU5D5ybVAKBB1i9G0Ed6khyUcU10ksByGXr1092RpJoi8qgjz7Q4wrkpCYvIhkiHKxoQwEF/b0bhJPjmSBcHewmp+gKtKL8xZS+adIR2g32OAl29QyyBYfYnz2l1fgMp0DQdcmtbVaG1r8TgFkgJCDbmDIH42mvvjHBXsb7K0iP7Kh3eDuIPR8mDt2MdFCcb9mQcRJhyLPvgHGcp0LUseXERB8B2v0zodXGBmucPW9LNba0E5XBkcd9MdAnsEhxsnWXznGXJXwGh3ZUAXRbwBylSOQk0hiIO64rqpvGADZKgOyi4FY15CRj7SQyCYJAu502Y7ppeVxEPS4cHQrlSG20JEMEN9W2MJstRhE5efsLaxWtFwKUSbPdhpiba8cK9tabBVlp9jpoBZ+JEcibsFslZy9TCt7M7RB61ogBcPvVsCSFYxUkbIDtZeqdc+M+WWcBa00kuovH2v5m0M2MIAP585VyvxCEP69CgAkMfo3ruCTN50J11jA0a+v1JpaL3ryjTlbi0muckGEah0Ekbq4ViO0+bu2W814qEZCOB+hb+kxNX3IwSQneSDSS1QAyA695yMM1aFOgLNPmvjbm197HEW/vTkfoH/sjrfNbt7Qes9aAQj42T8feh3xA3xE1hHd33P6AIK8obzm5m3v7vWqn3yvZbnUXroriL4XKgT5QJe83EZBYFytbgnIHg4W3rwBOwF+WMvyJgQ5aLl3B7Gp7n0cBL6VXv8qAnmDA+r1n1zA1/zswZU5AxA72NUN5OgGy+/jO+uDHD+nIHtfjbnRf3hwItij1IKCHFxvsouOnL5oXSHUzcXEzsJ3FlGQN32sxmjFHuYS0BJoBUIgB183mnYCOfp2dFQgx2qORroDkDc7NIX2Ui2CPTq4TkDgrKQEcsBpd8WCfLh2dFQgH2ike1kbBgry1nhvbMNqEKxdGOQAi9KSSA7Es09ZkA/zerw4xwcdvmdATgA31BhTkDkFQWfZYiAg1hJAPvQXmeAgIIsKgzhZNIrvM4DT93yj0Nt1MAjuP0h7i/RqJJCsBW9FWFEQczu8fcY9Suhummzi0b4mCIIccEdI2luke/beBQTOAFAQp2YtQFCPLqokbrNOzwWyt4NVViT81loyW+tA7loQQaAB7wYS2VrYjRcpOzqnE1ESZM/fO4JIyo7dOAOyr0Q/QicLRE8C/g8nkEMdDHQEgeYX5fYOBDlE0bPTMaKpKBIfhQKOE0k5iOgQ0dd/g7Wdhig01qJjRCLIQaendUjmfkWb/GkxyJFcq+hCFO77JrsNZTLtQBqSkzymSCBKtQDBt0gBEEaVcSZlO7HahrUEMQ8+JX9OOnJeSer/dUXCeOy3rZKQxKoR0NztvRYg9WZCf2gmRlqAINtzC0CA17UgONV9AydC3gKQLGVnqGo33AbkiFT2w4PgG0uxszcB2F7ZQhIWSGB+i0DSsS7/YwVm10SQmQ9bxiY5MVsNcAkOMeJHWknkqJWIFYk1Rfqvx0GagNgKYI+5/C+Rke1CkFkM5Fj/nywJ2OkfaZDTw9/AP6BA3y6faayWJPVLraTVsrv7lkG0V958iCDhs/pvn/nsYI9JaltYoiJJP0KeljCOQVIMQZYZINx6x3drxEEO+McJEF/dvOVIgDSRH3Ej01Y/bjJB0NGfhE7gMt4hGqIcQaWa33iODztEEErtFa8QHMgBjvZNoxxB8Us7ehEEOGFZjyAIql/XEcjeTrLlgLzDQbIERzNsRH4ggcCwqBikCRr9Oc0ckIM/XRgRiFzGk0FAmTYLhGRRLswKCtY8yLvnjgkEvLMqG8RH3JkgOPmAxwlyQN5tQy8mEHgN+TQH5Ihus8wD2YM3utUNBTAqnRQJsEdRTZdtmwzyUAqCY3b8f2aBaEP6HlsRb9MnCIzZ9/i9GHkgqRXz/32CwJh9j4cbU3srE6SSI7IedeTN1uF1VovnpvoBicXI/YJoZ64fmo5A9AMSqX4JIDB7KgFxQPQ+815A6JGtaR7IQ5lnJyDkqvxMkKj9PcRKFCII6CS3AnnAIDdZIHFPcsDfTh7IEVzQ/tEChF6skgWiaPEtAnLasHkgJJ/9DJCDGxLvEeTDnC2+T0bxHUB8At8IxB1nnqZBlrkgR/f1fJwLZG9uFjAgB7AJ+tMRW3xIdas/ZGVPgezdV2VAwB3fsvktBzGJeVIgrc2vOVWm85VmZ3lDORVBpHpkBOQjjyMJUv+fivuVO1/R1iDAdfEgyLMvYp69VX2eA8kIUWAOakF8KZYHyYy16hrpLEsIaRDygo9QRXD76h1mf2L9wTexV9H2tK0wJ2RCVYcHwS/4uOF/AeZIB0rGgchjBf8gDmPGb+MczfBPHGSPX/ARgthnWtsHzwCJDXp4kKOtx8TLJu7ahNsoyBs+xyeA+g5eHkhk9OYf8IQ5YeLRF5Zu4yDwZCWjIv66RZNaZIH4YahKqqIcaf9GEIgrLLlf4kH28Kwrx/mAutxToOwEBAYs8N6xdwEkK0s/giFSKzce5M38n0KBjs4dTLVjZ7bNgbRw3U1wUwEk6HEKAvFTtRZXANn78+CsCuGRlin2I/4pg2a0UVHCAUDAPY7i5O8RF3miIDaWmrGdBTpkpEHuQs9+oF2EQ9MyDJToH/SIj667LIH40fellZsEohuGUoeEgNRmS4He5BQJhGmHvMsg+NJWuJ38H/xZAO1aEyBNvsH/iAE5KPueCaAiptUZPngmyAzYW+DG6ctPEyCRxYC8+3AYqojuoheAkFfHAXkAElrk6Q1Ex00H804eqMkmjsgHUQTkFgrEhsRHehasPQhjtXSuS/pqh0wQ/4wSiJ28D0DW7rdagjyEwycHcxkdesIqvbUOYCRAADnaV14YEJLAdgAhr3OceouEnjpHRw5OkkkQk7kfq95AwljLkdA9M050f95R2C4ou62l9A8SRL/irkmcjgFmuw7/BD+iOoDsVWxcIMhHYiBxgVQwbD/CLFIEKdIRNBbEKAnJECP2KC4QlCLjEMWHvxSkwGr5po+0t5LphyWJqzp431UtAOWvTwVB41EESfmRvfEAMZG4KkrJWFO4s3DYzofxGKTEs++t4Y6IBNe1Wi68lW6FxAqb35JYi55c4ERyhyqN7dYhGN+FeY0ViFUSG2wVRL+uSCHuLVr7bQtCzS18ax8EASFKST6ivdR9DMRV4286gQQOkG/toKCxIEPcu9dmyyC0P9JqhbP6H8pF0V4gtK6Yl7OrG3nCl18dQILTE1JHBNVHs6oopy/6a0HwqCK/cupaJ5CbrK11NhCtEVGOrEpjDaKtVkLZzwiSsdwb7tzfkEHirr0/EKLsuSTparxSOQ6xJxAhtc0iyeiPGCcxy5uW7eRHhPJPzjqiP4nR79h4O3ntbaTfAsBExQe5IFe6vkUeMi4Qc2KiTYjiMscDn3/0CPKW5LDBlioPGkGZiM8/+gTRjjHOYS15cRjvTEmYj5wBJLHkpk2GQGzcMCWdrg4cbUH2Y7GNlhYIPB9aMCMXrrT5zRCIazLOS0XSHG9cNI0re353nhxaDNYRhWRtQVBGkKPerghxUOjosWvHlQoE5ygtQWAjPuO8Hmr4HOAUPhxYK95XMP1tCwKHkeU7BZwM4JCgO05g2qb2lTLFAsEjA+1A6JXySe0GRxnQgXO774pnaD78QY+uILCyFBeInUtsrAK46MjW9g7lhjfISNqB0CEyrNaEC1zCIICUH8xvBmkbI94BZC/OJ6pw+JdekNAPCHTIZwA5KKgNVkOQleoLBBxz6gLywM7wHvy8BNxZ+O6NQNlbgphbwOddQR4ZEHB//hSCoMI9Nb+tQI7wEFjvIAc4U+RVhFzrQhziF4OwOoKmvKYiiApvR+kd5KYLyAFfTiiB4KCxueGlvY6wIPG+IQK5w80vAwLvGJ/JIOial0P09HQMRLRa+1h6ixNG1iGSCzynEghMrIznuW0B4l/DS0FUpAFKGNlYixwxE0FgqmvrEOUg3kDe8iCsSPaEEfWM7FE8Or8tgsAZFZW69EHSdn8RIwXZm1kSTiRmhAFru3/no1cRb5IjIHCQwPSiykGs7w37I6AFxQgEV4Lx0VEOxEeSyCG6EMDEY4eghpu/t2qSignjbZmUBaFNXuYA8YFeci6EKGgd0NxPqUhgQygLJKxp7+Hx6ikHModxC/XkPYCQBmkmSNhlgA7pPQrC3nPmQVRbkCO2d21B9uAKgjgITawYkDY6QgY2qdXilZ1ru4NLIVgQH20pnOrilXfjZ8ZCIGOhJcLOD9hN4Q7m0/OkRA/44cbasLfxIxGQiEPkO3Hu1y0Ib36ZC9vw3mrl2WMg1H3HrJb9C/q3ORBUWgknssnPugskcilrCCJFlOaR6AsF6PfOcrSNfqMgb1Lwm2gp2u8W38M/fc9b7fKROIi44i1FZ2ZNALbKKgjD9YkgsVzFfrXjshL9V4DEZzWd3mbfTvV1INEs2IH4N7WV7axPBIkur7bu3XllAvnTQFxZYV5lHLX4k0Hybqf6G0BcelHG8eeB5Iz1/x0g6SHyvwak1fof5D8I8i8rFUXTyFgUuwAAACV0RVh0ZGF0ZTpjcmVhdGUAMjAxOC0xMS0xNFQxNTo0MDoxMy0wNjowMByhwjAAAAAldEVYdGRhdGU6bW9kaWZ5ADIwMTgtMTEtMTRUMTU6NDA6MTMtMDY6MDBt/HqMAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] @@ -330,14 +337,14 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAGQBAMAAABykSv/AAAABGdBTUEAALGPC/xhBQAAACBjSFJN\nAAB6JgAAgIQAAPoAAACA6AAAdTAAAOpgAAA6mAAAF3CculE8AAAAFVBMVEX///+AgIA4vPKTUVDp\ngJFyEhJNv8S4mNHKAAAAAWJLR0QAiAUdSAAAAAd0SU1FB+EEAxMKAvx4/yYAAB1WSURBVHja7V3p\nWePKEiUFKwPLGWAywGQAZAD5h/CwpO6u5dTSLZmZud/TrxljSzq9VNd66unpv3Kd/iPX/4H8bdf/\ngfxt1z8JZJrn+b8AZF6u8z8PZJ6vP9ezQPLvAZl+cNxudyT/NpAfHC+3+3V9Pv/TQObLy+3t4+Pj\n9nplU/LvAbnePpbr9sqm5FeAIHE5eqvLy9vHhoRNyXFAprOH4ygkdUI+Pt5fnx8AZHZeFsr90efU\nCfn4YGvrICCLbH/GSKb7H5+PmZJ1ZX1+f38va+twIPPlulxw2PEBNgrkvrJ+cHx//awtuklGgbAB\nLrL9FSLZUB4D5L5F7hNyn5L31zwQa93zNX/Hcd99PxJRf30qQI5YW/P1bZ2QFQjZ7T6QydikE9vZ\nPwv3ViXiWX93XXbPBwP5uu/2LJDZA3Jud781ifhgIC/jQNDzGRByRolDagSIf+IcDWTmQO4T8rlI\nRDAlfUCCE+ehMzJXSXKfEgCkZ7MHSsCePRIDeaGSRK8tT/xO4tMpUAImLrXS4nfOSK3LSxHt8ubl\nJvaBKO8fArkkzhHw+ylxjkwEyBcCsq4tc2bP4r4BkJd2st/gyY6PjMm6aXuvZZC+4bptd352hEYn\nkLcFyZfcj0/mPbMXBfKNgDhvpx4aLS2hxgMg4R3yQPRNzFvvAoINq3kenRIB5AaAnKx3mzAQ92lF\nHXrnit3hQNCM2L9Vwx++xc+ULHrEu3CjECDXg4D0/Fg9MjTvp8v19vPAm9S0G5DLoBUXSi0fSeIT\n8fc7Em370BkZNEejc+Sgq8yVofLs3yOhinIcjoIE+QCOAHLd/AHraXt+DBByGEHR5pwj2ZOlOgS+\nPjqFVsc13TexN9I2kEXXyyChhtXDtshsK2wciBLg25bKeHHM09YB37uMQ5/SE/kq8/BcVg9PxrBr\nU/KeXFndTtRNVKWAiBFYPTwqDIGf0k7b1NtN3U7UKXSOGUCqSpN6t+WMerulUC/fjxb8cUDasr9l\nVssSDfs5bLOekmidBEAmMJ8YCPPwJIDES5h/+9rre2RAoCg2gKwTspwMqbO6a9Ub4czsSOFlhoHU\nw/orLYk6xNC8SsQuvzYRv/jnGMilqk9fubXVdW36ayeQIh/uMvKmTwYIZGIensPVJ0MRd6+2L7aD\nYZGn5PcGkGZi3DfJXwCkKo13w8pzPphAHmBiDAEpyoAV1sVAmmPygUA6g1ibemb5g/4QkG6pVa91\nQj6/ZVj3TwC5nyMvfedI++3iMwVh3V8Aonw+Ayc7AWKEdY3NfiQQZHsOB6zNsG5K/JpAMsc58Ab0\na78EiBEgSRyI9jmSclggj9l4UgcLWdH49B4V5e4LC636VcV7lrtkl8sGeQMNpbFkfHw7jpFNC41c\nnDPa2DmbfdJwe4FcaaB2th5zDay25VbzuIQCK9BMfUgYVtZev09IbFaMA9kGiq/KqXoDM5s9SmYo\n37km/CzVrTSCQ9/fjIZaQMqU2M6HkhIUrK1hIPcszNvP87lHw4yGWl6UzcPzbvse/Nym9j6DS8vw\nMZnRUNOvtWqoMimVPql94xFAZsPrt+z2rPZ7WmXry89L2m9ZgNx8IPAcyQAx/LDW52bmQ3hKpIFk\nneHi7oZnnPhnc0CKFHcelVpag5GX1SGFYhVl7wjXoZeLEgSKc5t9MIJvR4/uA/i2iKGDEpiJ+A3e\naUSvcuJ5eAB3pMuuFmvC4TuiHzoRVpzJtg/I9eX1oPRRcS3nHol5x+bMDiBQFzoMiBe8R+KjEwi7\nw3rWPWJCwnQKLT36gPDMq/7A02FAtG3WBUSqG+Ox+d1AJmlB9wDZFGuitjxoPuIEF53H0QOklGm9\njnjWjgWinZUdQLYjVWVKPeTyM3XWGMk7HdMOIFXJef+FKfEzdeZLjZGUMc0DmUik96AUIEdUuJk6\n7VXamHYAaYp1NsGhvRb6uiu9q/aLHFIoRpIHQiyazhwg/L5+SJvbI2KvoxhJN5Dv7rW1qEYaSRE8\nYr+VlGkvU2crffumaQAdQDar/zsbe684oEq22mXvAkkzXqpo0QtZlL7NnUDqr7+60rImI72DOknO\n7cuXYhfXoKd2SDGXUHmVTiDf3UBmbNrD1KiJOsFWvxbKT4JOujSQNgzfPUCqe0roA3ALFNALksn0\nZAogz78DxMiXmplQOm/fXXKTP7bZM+0dErVq6ksPkJGlZexZHtPc/vCjyzH/yJaLLHGsQL5HgQxu\ndrqEKJA6qC2mqVWHCR+Y+2akBX96xK8Vn5hBKDCtOuzbI5vT9XP1HiMgaPzWJaQdhvdR+RSrg6oO\ntw4gveI3es6EVrSoUayfMy1903jSJ+7OcySYeaw6Uc8n/RUBUopn8tJk58neamngylrPAHmAz4YL\nl7itGhCyYPzqjX26FjGsUEVuPcDZZI0B+Qr0633ab50SaOrO61km6Up6gLDTIdjte+yRzWv9DqNY\ntKSdhfw69gjJ5ImA7LIQqwWB4giW0t0htTqA7LPZT16E1jKDOs6RHiCrF+Vj0IvihNESgbLwZO/Y\nIzv9Wic7GWYbeFCbP7/AqQK6VofU2utpXO8AP70QnwcbTCOFAmi/8hxxgeiKsWMIXtoASyBtbXET\n8aLskU6DZ5c33ryc0GWzEK+BhTikXv8aECOlhNvss8Q2kMb+cCA1+Bt6Ufa4Mn8BSPFrXZFfixUP\nB1rpLwG52S5ng7JKJ5bvcvcfKLWg+F2RoNNHF+aUAEyU3vJIIEWlEkpufWeDguPC1QTm1/pDQICb\nKv6VNPIXT9ZtLHR/FHVbIpkzBW0suetAIHYQYOoJYI+TgTyBGw3cxizX7X6z0cA9BzKcy2CV65a1\n8vDgKQdSXTpwqM7eBxvFqKShW5KcwWn4YCCruQKR6DXH1z7Op6vylFmND8lgoUCcLEXMpHBm/9Wa\nyLT5VpjfZdpO+oORUCAOSyGsZkFlOqK6H+Rzu+SUhwBxonYoMovZcoCfkcuyPad3EkjLfVZ55Ej4\nANe7LAoDnge42o4G0sbvBoGIdRPKVFRMnqge2A3EqRmpsXL2YTSiqJg8U8+xE4hHbwKBhFfL8CGu\nkmOUMh/IzfBNDQIhtXPV45aqeToEyCdkFRgEclNFgKkqtL1AzLT6PUC+G5AZT9LjgIDBylXqOUC2\ndD40SX8ASKdagYAwj/uhXFy5pWXlXI0A+eQfPQCIs9mHqlnAHpl/B8hLO4dlUHVE9QYbgsR3dBHC\ngUDMA3HIBgakJL8CJPKE9KvdugabA3nQ0jpeDdKn369sdhqTOQiI0kceyIqRNKyGLq0hPpACR5q6\nbx8GIX/25alsa5HQGl6jguxBKsqpEuj1ceEIGKx2qaYRlDl+HE2UdActpfaj/rRJOCBKnKClmM4s\nn/C5od8NCTnohv2CS0EfdZAsAUQ6x8X6JbN0UKXWQS7T7ddyaao0bBAnPKh27iAn9nIRYXGu9xMJ\nCiyT9Ny+sdvNdWDXJJirLLXN+qWybUooca+b60AgONVJzvE6be8lTmhk3v1pIAkHybaRyrYxvPh/\nEgjTnh3VgNv/vTSovwLEbUNCkcyXstwaNd7etXUckBna/EgMks9ILuRORfVIIG8k+6EUIAQHUz9V\n8C8AYWrUqtiG/jCQJPQXAamenimkzyUie+duzwHJuFAuzYwtQNqZYRMp/SqQrXb87H5JA0nQOP8u\nkBwhk15aCSfAr272OUdYzjb73YrNxEJ+E0iJ+kVZVLM0/krA5cuJhfzmOdKifq/uKm7nyBafYh8Y\nrznBAtEHAUkSltMWC8s6SdnnhYZ4zabdYfWGQNKE5dLlKlkr8WzOJUfzvgP3GKghEOZl89cW39s5\n9mNqHveRtncDIX7Pd7/yiUtb7lU0HdaESinHlzgIpCPqJxpVUj+v04So2exZwq5RIHknZ6uzehbR\nKa+b0hIMu2eHkO2SASISRkIg3O3sbxKWJkVDCA0I0NqKe8KoNzNh9FUrZIJltQKCUYdBIDNql7X+\nmhqL4dqalDcsBBLHmErBkkglREvLI23simpoudAD5BsCofwG9CAAQFZWb4P80WGrAKBVrl8vEL1l\nGZ09Nce1+HWTzno6/YDtxJoNjQDZyife1RbVB+Kifr6Br/YCQarmU1kfhnIQLi2D//EEVBS3RrYH\nCPrq04bDUg6iOKxTvaL0er8+JA8EBtJXIHaTgCgOO1mtL09ajW8VtjD83QEE5Gg8rX9ZFi/iMpNp\nMGf1+Do4UoHZDKtPog67BfFpILx+iAJxrJsgfAmL/cEKWNVhtxisAwgiEXs6ccNZ3YKWOWo7r5WE\ngipIEgtd1OEskPebr/9yfw3hfGAF+yp12aVDcfkNqtKxmX8aCKWPz/shRA15A2K19Fjvz8ocnS30\nBbK87rHQt0KAtO0Z8qqsTDXth8CMgU+yYlhuZ39BzG4ROneZaiDs8Ep3jMNklE8nm1O+vKuM+ouF\n9/ltywKqJOoh4QHGbPAK9799OjnN27fftai/mvEACDMbFBBRyp3lpjZnZGZEEmj51Ki/umkEhNYC\nUO6MNVeo8ZttU5JiQjaBAGoPtT5eDMb+GAgbEX4gziJ3NRlyl6TGBYgA+G6y6cxQZX1LM05QXWsD\nIhS8ar77F6QHFUC+oH9BGcjtL674lUBIssAZAUny7s74QEwAsasDIdWVdQlqNrW0TrqWxgRC/P5d\nQOx3AwxR9htsYqNY9WKzpx+KGCZ2Ajk5SiN4g5lZC4pJo2v0RIFIZrP7g2yq8ejbLPY1jeR3L7cB\nEaSE+A1GxzSs4JQwqptR5zuK6SUOxACIYy3pl+Sp6dPYhEgvcwPiqijx6LxhnRWPN3vzaTAaAko/\nY6UxHB3DHZRr6DYCA1YGh2p8PDq43+S8qPDpnsH9DxXaDDGsvqBhlRidC3DoTu3MOOtf7M+OVRwl\noambQQKOA8JJJt/6mPRYFFbYm6IDhI/Va7WM5cB5rqGwO0TuoMHrZ9+9Yam8p8+md0UOutHRMrpN\nEXmDWUT3AdnXVxVc1W/3JTWXzTCH8TXLj55ahpETexSIpvPd3so8Pw2e4impj0VhhbFrtjp7Wu20\nTpYdmq47iAI9g0Asrcdx8KpShxXH1SQGMYAcepk+CcTp3V5Zh+XmpIvI6b63a34AVen6tqYQwF4t\ni4ovD2RULSVAkK3G+hQL/vW1ZoPrNJW3PA5ZWx0qHwTE6lMs8z/Kh3liBQfIjrVlATGlWXPMUmlm\nUbr+3UBQFY1Xmf5rQNBmNxsuQ+tVx4XVOz4aiCV+basacp66yYTiCH/MZpe8xKkZqUBqyYaXTCj9\nxI8Rv9ZecPeImhGvK5diebQPxB04TKXRFr9sj9Rv23kwa54TcQk8RkUR4fdz+9zSilHptgvkKlwe\nBwBBvgTDsOLEyxEvsADC2n7MteqP5aLsxAHccIaWa0fCgc/NAdICp1kC/NirCQ2ftuRlw2XLVdxi\n8kShEUuLALk2J00KiAqGKXJvwySjTyJ/Q80Eyn2q87B+ZO8R7S7xgejwJKrFQ0GOSQR1KhDbTaRi\nutPlzRK/2mJ2gUxKNmixDHvDnoB8dGfqRCrk6sf2gdiaK+SAgFwETaRs+arxH7S44T9gy5irKONA\npjaBTW1IA7HCCnaynhIbqLRxCEjGlepEDyYk7zSDNEfC7m3VBHUCSeW0xWEQ8X4dnier6Qfd7NcY\nSJmQrw8vPhgC0RHPFccUGwp23ZyWGS4Qj3q8vZUfLLfKdVMORLOSUXtePSBXS8Nj3/KD5YJQgHye\nSTsxa0s1fb4DxDYeBBBvyQuKB4Yjwb0BnUSnEwiIeUBAGzP0NQ9IG1LR7OKSizBanRm0FHeAyIZA\n5rOcxW4UMebLqeh5T4P0VXlKGFaXZCKBI36sIkbc+dRAUoINE6NultIyDWSIeqkUMS4SnPYacsI9\naJw25hGqwcq9GQAZzlHZXpgXMdY7jAQtZRWNWNIHLC1nMK0iRiduZU+M6kXOVvT+ze4CwdSljlVi\nT0hxzBv5zfvFrwfEKGK0DGF3TF6MQz4HJD4QIyCfDcgrARKr1RKIl6KPgDRddWZZQyPUElYR4wCQ\n1b/0bVM3SyDzSuxzPp32Zg0JINR7wJZWDkir8jciDIrxjBykl60b6cfHIO1zCoi874QUhbiIW3LQ\n0dDqzqwhB4g28DgO0NY2KuLmQHgHbtsJtRMIMPAaDtg5JOYO5fSG3JNM6pnHuD7MimXTKWR3c4+K\nuDnhpKj7LdljqQpzcJlFjFvoEyRCGgkCcRH3k/w6XUv7yAu4q5DHSYxEQaube6L2mXOZUjLpUhjx\n8jqcJSaKGMkW1f3Rto9x98ROIIB8bSetZSti/FLxDeSTsJ0NfUAgZ9me+C7iYN3+grxIjvuna4/g\naPCeMLUuYizjAo8LmzapS2rlqGX6gIhCN5Lkj0z9qlApPaTrHPGiwaOXqA+rZRfnE/CHO5kOXSc7\nFZbfWSCBE5dzsEaNhi6W971P18pSy7AXjSQBKWI8E8oxrCg4KRtzh/brxB6dCYmmpOgc8zm0Vj0i\nq3o0JOyRgRmZQiCE0yW0Vu0IFU/diizEgT2SANLCCsxa7QYy5232AamVqy3awDJrFQ2qy5GG1iXd\nofvOkcyM1Itaq7j20gMCWpKxM5X3VjBnFgxBvVcSSIJy7ML8NrI9r5TdXMvBula7z0TfU6sVHYpY\ngnIsSP2jp2mJTFQFmmu/zVNbGmjxajvdcii/shKUYz6zIy9irCbNZrlwe0QywGw+fYpj2DhJNFGZ\nWJ8PvYvosJLapwWxMKy4Fs15vWp525D9ngKS7wQ3Sf4eabNzXZXyeq1z+XEbNeBjyjFB8uIOF9Or\n7+8nvCgs+ZbLCZ9UKj5SkkCKqRsV+RNLZ/nqk/wzSVDnrBUBqVQYbs4AaYnl0bSr8vMndaPqNil3\nXadgktQTYj5EbNIHYp63q4fwPeRq1JxFT/JGq9tk8/0SW8JlrMn4jVJAzHC0/N71JlK6tTe+Jk21\ntVQ5Vi0gOlcSPZsBsdhmVfoZ/hojzbgLhid9pxpFZWtJ7S4+jrFvNUkmmmuPp5PvQcSqHjkmEEV2\nZNeyCSC+rtVGMsAhgNxH5Uneggxz4w74GeeZddfhD8rFH8JYTcdFWRG+JBDBDi8MutWa+KRtgsSE\nRBGhMHrWcdEcAAlEJlb1AEnF6Hg8c9eEeEBqYlUteoBLC71Gkmsh4KDrupylpcmn0pud1ek4i4bL\n851ArM3Oj3EAxGOWTGYWgGeMln7b4hfY9vJANBdGOo0IeBqnDhuTTa51IKpGWY6KIhdGOo1oMc1u\nItd+qAhK0yqWMnAw3mIl2Eyk+XwoZSxnq9U1EMkiVYGArcwtqRKJ1WpIR2LXJGZguI5+3o6uL6HG\n83O7bBJhWK0TpLte9GWoiSK8mjHSCUTKngqEnLrX6m1gFNa4IuS0K2fw7t0YohZRW7YAgYyl3Pmw\nOdZ1sq5VPpkCEuu6+M8zyV0zgTzT5zT7FXLrnE576KDiGbF8AdOW31mSJ3wg0kHHcJHBeRvNUGvJ\nsDaOyxUiESs9AMJdpsbo2GWS4RVKrbK6z+CnjPLG3SPoxuhtSpnk10fnXg/PEY+gh7MjUiAfCSAY\nncmtAV6N/zkA4no3UVgBpUzn5XsHnbV8bz/CIj3V+nb1p/xAJO/Sozyk019Baos7XPYhrC6ta62r\n3MhkM6dEVpyaOLqSjWy1KAGk5JyluHfLyM2ZFqb96V++mxYBAbkWXc1VU7h5z6okEEPjtoFw46O7\nJ2mCzLC/HVrA/Q+B8I5uc5rDnSCJ9PGeBiMVyFsvEK6z9z8zk5bScdrsAEJ1dtgXY+81cP4PAaE6\n+wQ7y+4FkiprlEC698iJ1qKaVct7rgEdeRBIY3WJUyBGgAzU1RDraYQobIJ85nuvkZKtjlosDOTW\nnxN4ChJTppEiutmPJMdA3jiQzDOjk0SMTlZstdMt+CoEYhOruDgu3tkuydBzQEzPzeOAGEni+4CU\n2uKEmjkORBjwER3qGKt7NDwHAJk4H0q4CIb2SEoZdYBcYiByqEyO73oNlv7OSUt1VPxuVWVVLMaC\nUlbMJoFk6fAiIF9GM8dSx32rjuHo6GInu9sYY8g/D0tcYxWF2qAFCHPFop/kdK27EHnuj5lAIKqx\nN5gQrl8TP4rd+JvxmJpDnqW1EEsOA0ENPtQXqI2UsByyaQU5VmvIG4+e+VGeCY258k417StjAqUS\nPao/xnebYCZ/8D0/SVJH8DJAUqk3OTUR91bAg+c5Hwjd4mZapIC0ZCjzJadYaKzvVxxLAU8jYk1k\n93kTUW7TlmOJ3HEPtJyVTpNIA6aa7ZmGg07nHRAWB/oC08xDuNdr4Gm0WQPxt0JWQEy13YDIlBZ2\njrTnq2B01DpsJjyOphbjND/FSExtTScZtRQQerIrj/60hgy9NfNWs2BtAimS+l4FuUleMTkRcAFk\nbWLMY5NtQjrCIbK/Fw8/tH/WtCpytNosHE66CwCydTHmCedhqDMA0vQx9jJm89PuCwBR6vAKJJMx\nSoGI1nF0SBoSo2ndfiDL5E5bTx66kxn9UuaajR54LPhr9gztv4DUQukEUxJIe0cLCMu8N1ttjgAB\n3N2TJqPjfcSMi6YEGEB4ywtswI4Bge0cZQpTbo+wpks+kG1usUthCIjWtSoSDjiRxMTUdmOz84P0\nSCDJLoiJEHcT28+zeY7MO4BM7kgqe8QD4k/IlartzIo0gXTtEf8VcrpdwgMiTGSjubgE0iG1Am9S\nOiYWVomTfleLptPoU2n/FRNIdI6EJ0BHJN+95NTiqRZac8fJHvYqUn6twUsS5eOpFt2pOnQtixKW\nDuVLT2aEPR5ic1t2DR3XDu1XOyn0nGXd5D4QKW5x1z5eHNhhjyyn1KvrWxptistvog5A7NXnGYN5\nC9Frmy1uvgcG8mLUzpaCBIk9LGezz2c7w/foC2mfOCLCy78yXpR7PccfBZLpYJnxa/0AOaeW1sOA\nJNrRZDyNdyArDcHrnwKS+eFFdrizgXRaqcNAxGbP/jL2xs9z5kA8CEhPn3r5mon4SK28eTQQ0Xyo\npzzg58fsf6b2mzFSh0pYyk/vv3Qccr2X08cqOO7mwP0ZDcFcwyrK/jgQSMIkCh3SHo4iSbD9cSSQ\nqNBxDxclESVZU3MHkOhl4qCNMwZVb9jN+7oXCOFr7bZIZGKuDkek7xSK38TL1CDjrXdK1vLGd071\nMUAbw6XmKJC+xGCqPQkC+BqO650QUdU0BoQG4hMHMvPJ8BLgcVNT1KIMAqHJyGFiMK/QreUE74UU\ndMjUlOJ/DEgfpTxXQVXB+ZipKQs9xoEoz5I5diUv8bnQrbwJgTv1C151jo0BsdsugMijIGEAQAau\nSVbqDHrjrfzEWdupkiDhGCDqQD4WCIrO81qhw4DIMqdBIG8whxfmS0jujYjRNAtkYwG/7QXyCYBw\n/vz2SKZNSfE7dOkisCOBwJwiVRgoDsQ/DATuEZjlpYHMgB3laCDJG2MgmJxQ83swpXGY4cUDEjvG\n2suB4BfnGH9u35VAOCfGmAe50f5pIK6qwMcNHoiYwFMD4Q7Z2D9qABFteFWxGL6nwAh1LVxiBqhj\nqNesp0ibvZAUkLDqDc0Hw8hiRiX4hfO3ARAaoezOhCpvJNtiNCAllwTdc0thILfRPR8N3mNI5tOG\nbItF9QMpZ6+OjzghKBXk5TEjBKRpkoBeqYmVSflws9dK93sFavxk10jpIC8o5rVIzl16sDB87E4J\nfd8UEO3TnnR5tVUHNHsn+TgQyZ6aA6KjDLrg3QICec7aQ4cbnAh5NwpEGTZ2ZZZmt5NAxtz+/CgV\nUgtvdhR2V+QIZkeDmZu64ta9+ZvWxYBcjLGB+QNyUZiNACYv/nUX7CPniAPEORDxm4iv240AnEU7\nfLJ7QOw0MSMSN7FvO40AvKSpQV3LA2LHRDhhmHU5DQUcFXdQ+3WBnKxb5kKKPg+/eY3ZIz4Q+x0z\nsqX2OvlKOYQPvrJAElvyIdwdBwNJ0iIPsKn8NpCUFUw6tf36yhqOIWK0tXfer0/IoUAIB9Sjk1ge\nDKSTlesvBrLDvPi7gBynOv1pID1cgn85kD91/R/I33b9D2BqKIQldvOpAAAAJXRFWHRkYXRlOmNy\nZWF0ZQAyMDE3LTA0LTAzVDE0OjEwOjAyLTA1OjAwNSCOuQAAACV0RVh0ZGF0ZTptb2RpZnkAMjAx\nNy0wNC0wM1QxNDoxMDowMi0wNTowMER9NgUAAAAASUVORK5CYII=\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAGQBAMAAABykSv/AAAABGdBTUEAALGPC/xhBQAAACBjSFJNAAB6JgAAgIQAAPoAAACA6AAAdTAAAOpgAAA6mAAAF3CculE8AAAAFVBMVEX///+AgIA4vPKTUVDpgJFyEhJNv8S4mNHKAAAAAWJLR0QAiAUdSAAAAAd0SU1FB+ILDgkoD8Jk3r8AAB1WSURBVHja7V3pWePKEiUFKwPLGWAywGQAZAD5h/CwpO6u5dTSLZmZud/TrxljSzq9VNd66unpv3Kd/iPX/4H8bdf/gfxt1z8JZJrn+b8AZF6u8z8PZJ6vP9ezQPLvAZl+cNxudyT/NpAfHC+3+3V9Pv/TQObLy+3t4+Pj9nplU/LvAbnePpbr9sqm5FeAIHE5eqvLy9vHhoRNyXFAprOH4ygkdUI+Pt5fnx8AZHZeFsr90efUCfn4YGvrICCLbH/GSKb7H5+PmZJ1ZX1+f38va+twIPPlulxw2PEBNgrkvrJ+cHx//awtuklGgbABLrL9FSLZUB4D5L5F7hNyn5L31zwQa93zNX/Hcd99PxJRf30qQI5YW/P1bZ2QFQjZ7T6QydikE9vZPwv3ViXiWX93XXbPBwP5uu/2LJDZA3Jud781ifhgIC/jQNDzGRByRolDagSIf+IcDWTmQO4T8rlIRDAlfUCCE+ehMzJXSXKfEgCkZ7MHSsCePRIDeaGSRK8tT/xO4tMpUAImLrXS4nfOSK3LSxHt8ublJvaBKO8fArkkzhHw+ylxjkwEyBcCsq4tc2bP4r4BkJd2st/gyY6PjMm6aXuvZZC+4bptd352hEYnkLcFyZfcj0/mPbMXBfKNgDhvpx4aLS2hxgMg4R3yQPRNzFvvAoINq3kenRIB5AaAnKx3mzAQ92lFHXrnit3hQNCM2L9Vwx++xc+ULHrEu3CjECDXg4D0/Fg9MjTvp8v19vPAm9S0G5DLoBUXSi0fSeIT8fc7Em370BkZNEejc+Sgq8yVofLs3yOhinIcjoIE+QCOAHLd/AHraXt+DBByGEHR5pwj2ZOlOgS+PjqFVsc13TexN9I2kEXXyyChhtXDtshsK2wciBLg25bKeHHM09YB37uMQ5/SE/kq8/BcVg9PxrBrU/KeXFndTtRNVKWAiBFYPTwqDIGf0k7b1NtN3U7UKXSOGUCqSpN6t+WMerulUC/fjxb8cUDasr9lVssSDfs5bLOekmidBEAmMJ8YCPPwJIDES5h/+9rre2RAoCg2gKwTspwMqbO6a9Ub4czsSOFlhoHUw/orLYk6xNC8SsQuvzYRv/jnGMilqk9fubXVdW36ayeQIh/uMvKmTwYIZGIensPVJ0MRd6+2L7aDYZGn5PcGkGZi3DfJXwCkKo13w8pzPphAHmBiDAEpyoAV1sVAmmPygUA6g1ibemb5g/4QkG6pVa91Qj6/ZVj3TwC5nyMvfedI++3iMwVh3V8Aonw+Ayc7AWKEdY3NfiQQZHsOB6zNsG5K/JpAMsc58Ab0a78EiBEgSRyI9jmSclggj9l4UgcLWdH49B4V5e4LC636VcV7lrtkl8sGeQMNpbFkfHw7jpFNC41cnDPa2DmbfdJwe4FcaaB2th5zDay25VbzuIQCK9BMfUgYVtZev09IbFaMA9kGiq/KqXoDM5s9SmYo37km/CzVrTSCQ9/fjIZaQMqU2M6HkhIUrK1hIPcszNvP87lHw4yGWl6UzcPzbvse/Nym9j6DS8vwMZnRUNOvtWqoMimVPql94xFAZsPrt+z2rPZ7WmXry89L2m9ZgNx8IPAcyQAx/LDW52bmQ3hKpIFkneHi7oZnnPhnc0CKFHcelVpag5GX1SGFYhVl7wjXoZeLEgSKc5t9MIJvR4/uA/i2iKGDEpiJ+A3eaUSvcuJ5eAB3pMuuFmvC4TuiHzoRVpzJtg/I9eX1oPRRcS3nHol5x+bMDiBQFzoMiBe8R+KjEwi7w3rWPWJCwnQKLT36gPDMq/7A02FAtG3WBUSqG+Ox+d1AJmlB9wDZFGuitjxoPuIEF53H0QOklGm9jnjWjgWinZUdQLYjVWVKPeTyM3XWGMk7HdMOIFXJef+FKfEzdeZLjZGUMc0DmUik96AUIEdUuJk67VXamHYAaYp1NsGhvRb6uiu9q/aLHFIoRpIHQiyazhwg/L5+SJvbI2KvoxhJN5Dv7rW1qEYaSRE8Yr+VlGkvU2crffumaQAdQDar/zsbe684oEq22mXvAkkzXqpo0QtZlL7NnUDqr7+60rImI72DOknO7cuXYhfXoKd2SDGXUHmVTiDf3UBmbNrD1KiJOsFWvxbKT4JOujSQNgzfPUCqe0roA3ALFNALksn0ZAogz78DxMiXmplQOm/fXXKTP7bZM+0dErVq6ksPkJGlZexZHtPc/vCjyzH/yJaLLHGsQL5HgQxudrqEKJA6qC2mqVWHCR+Y+2akBX96xK8Vn5hBKDCtOuzbI5vT9XP1HiMgaPzWJaQdhvdR+RSrg6oOtw4gveI3es6EVrSoUayfMy1903jSJ+7OcySYeaw6Uc8n/RUBUopn8tJk58neamngylrPAHmAz4YLl7itGhCyYPzqjX26FjGsUEVuPcDZZI0B+Qr0633ab50SaOrO61km6Up6gLDTIdjte+yRzWv9DqNYtKSdhfw69gjJ5ImA7LIQqwWB4giW0t0htTqA7LPZT16E1jKDOs6RHiCrF+Vj0IvihNESgbLwZO/YIzv9Wic7GWYbeFCbP7/AqQK6VofU2utpXO8AP70QnwcbTCOFAmi/8hxxgeiKsWMIXtoASyBtbXET8aLskU6DZ5c33ryc0GWzEK+BhTikXv8aECOlhNvss8Q2kMb+cCA1+Bt6Ufa4Mn8BSPFrXZFfixUPB1rpLwG52S5ng7JKJ5bvcvcfKLWg+F2RoNNHF+aUAEyU3vJIIEWlEkpufWeDguPC1QTm1/pDQICbKv6VNPIXT9ZtLHR/FHVbIpkzBW0suetAIHYQYOoJYI+TgTyBGw3cxizX7X6z0cA9BzKcy2CV65a18vDgKQdSXTpwqM7eBxvFqKShW5KcwWn4YCCruQKR6DXH1z7Op6vylFmND8lgoUCcLEXMpHBm/9WayLT5VpjfZdpO+oORUCAOSyGsZkFlOqK6H+Rzu+SUhwBxonYoMovZcoCfkcuyPad3EkjLfVZ55Ej4ANe7LAoDnge42o4G0sbvBoGIdRPKVFRMnqge2A3EqRmpsXL2YTSiqJg8U8+xE4hHbwKBhFfL8CGukmOUMh/IzfBNDQIhtXPV45aqeToEyCdkFRgEclNFgKkqtL1AzLT6PUC+G5AZT9LjgIDBylXqOUC2dD40SX8ASKdagYAwj/uhXFy5pWXlXI0A+eQfPQCIs9mHqlnAHpl/B8hLO4dlUHVE9QYbgsR3dBHCgUDMA3HIBgakJL8CJPKE9KvdugabA3nQ0jpeDdKn369sdhqTOQiI0kceyIqRNKyGLq0hPpACR5q6bx8GIX/25alsa5HQGl6jguxBKsqpEuj1ceEIGKx2qaYRlDl+HE2UdActpfaj/rRJOCBKnKClmM4sn/C5od8NCTnohv2CS0EfdZAsAUQ6x8X6JbN0UKXWQS7T7ddyaao0bBAnPKh27iAn9nIRYXGu9xMJCiyT9Ny+sdvNdWDXJJirLLXN+qWybUooca+b60AgONVJzvE6be8lTmhk3v1pIAkHybaRyrYxvPh/EgjTnh3VgNv/vTSovwLEbUNCkcyXstwaNd7etXUckBna/EgMks9ILuRORfVIIG8k+6EUIAQHUz9V8C8AYWrUqtiG/jCQJPQXAamenimkzyUie+duzwHJuFAuzYwtQNqZYRMp/SqQrXb87H5JA0nQOP8ukBwhk15aCSfAr272OUdYzjb73YrNxEJ+E0iJ+kVZVLM0/krA5cuJhfzmOdKifq/uKm7nyBafYh8YrznBAtEHAUkSltMWC8s6SdnnhYZ4zabdYfWGQNKE5dLlKlkr8WzOJUfzvgP3GKghEOZl89cW39s59mNqHveRtncDIX7Pd7/yiUtb7lU0HdaESinHlzgIpCPqJxpVUj+v04So2exZwq5RIHknZ6uzehbRKa+b0hIMu2eHkO2SASISRkIg3O3sbxKWJkVDCA0I0NqKe8KoNzNh9FUrZIJltQKCUYdBIDNql7X+mhqL4dqalDcsBBLHmErBkkglREvLI23simpoudAD5BsCofwG9CAAQFZWb4P80WGrAKBVrl8vEL1lGZ09Nce1+HWTzno6/YDtxJoNjQDZyife1RbVB+Kifr6Br/YCQarmU1kfhnIQLi2D//EEVBS3RrYHCPrq04bDUg6iOKxTvaL0er8+JA8EBtJXIHaTgCgOO1mtL09ajW8VtjD83QEE5Gg8rX9ZFi/iMpNpMGf1+Do4UoHZDKtPog67BfFpILx+iAJxrJsgfAmL/cEKWNVhtxisAwgiEXs6ccNZ3YKWOWo7r5WEgipIEgtd1OEskPebr/9yfw3hfGAF+yp12aVDcfkNqtKxmX8aCKWPz/shRA15A2K19Fjvz8ocnS30BbK87rHQt0KAtO0Z8qqsTDXth8CMgU+yYlhuZ39BzG4ROneZaiDs8Ep3jMNklE8nm1O+vKuM+ouF9/ltywKqJOoh4QHGbPAK9799OjnN27fftai/mvEACDMbFBBRyp3lpjZnZGZEEmj51Ki/umkEhNYCUO6MNVeo8ZttU5JiQjaBAGoPtT5eDMb+GAgbEX4gziJ3NRlyl6TGBYgA+G6y6cxQZX1LM05QXWsDIhS8ar77F6QHFUC+oH9BGcjtL674lUBIssAZAUny7s74QEwAsasDIdWVdQlqNrW0TrqWxgRC/P5dQOx3AwxR9htsYqNY9WKzpx+KGCZ2Ajk5SiN4g5lZC4pJo2v0RIFIZrP7g2yq8ejbLPY1jeR3L7cBEaSE+A1GxzSs4JQwqptR5zuK6SUOxACIYy3pl+Sp6dPYhEgvcwPiqijx6LxhnRWPN3vzaTAaAko/Y6UxHB3DHZRr6DYCA1YGh2p8PDq43+S8qPDpnsH9DxXaDDGsvqBhlRidC3DoTu3MOOtf7M+OVRwloambQQKOA8JJJt/6mPRYFFbYm6IDhI/Va7WM5cB5rqGwO0TuoMHrZ9+9Yam8p8+md0UOutHRMrpNEXmDWUT3AdnXVxVc1W/3JTWXzTCH8TXLj55ahpETexSIpvPd3so8Pw2e4impj0VhhbFrtjp7Wu20TpYdmq47iAI9g0Asrcdx8KpShxXH1SQGMYAcepk+CcTp3V5Zh+XmpIvI6b63a34AVen6tqYQwF4ti4ovD2RULSVAkK3G+hQL/vW1ZoPrNJW3PA5ZWx0qHwTE6lMs8z/Kh3liBQfIjrVlATGlWXPMUmlmUbr+3UBQFY1Xmf5rQNBmNxsuQ+tVx4XVOz4aiCV+basacp66yYTiCH/MZpe8xKkZqUBqyYaXTCj9xI8Rv9ZecPeImhGvK5diebQPxB04TKXRFr9sj9Rv23kwa54TcQk8RkUR4fdz+9zSilHptgvkKlweBwBBvgTDsOLEyxEvsADC2n7MteqP5aLsxAHccIaWa0fCgc/NAdICp1kC/NirCQ2ftuRlw2XLVdxi8kShEUuLALk2J00KiAqGKXJvwySjTyJ/Q80Eyn2q87B+ZO8R7S7xgejwJKrFQ0GOSQR1KhDbTaRiutPlzRK/2mJ2gUxKNmixDHvDnoB8dGfqRCrk6sf2gdiaK+SAgFwETaRs+arxH7S44T9gy5irKONApjaBTW1IA7HCCnaynhIbqLRxCEjGlepEDyYk7zSDNEfC7m3VBHUCSeW0xWEQ8X4dnier6Qfd7NcYSJmQrw8vPhgC0RHPFccUGwp23ZyWGS4Qj3q8vZUfLLfKdVMORLOSUXtePSBXS8Nj3/KD5YJQgHyeSTsxa0s1fb4DxDYeBBBvyQuKB4Yjwb0BnUSnEwiIeUBAGzP0NQ9IG1LR7OKSizBanRm0FHeAyIZA5rOcxW4UMebLqeh5T4P0VXlKGFaXZCKBI36sIkbc+dRAUoINE6NultIyDWSIeqkUMS4SnPYacsI9aJw25hGqwcq9GQAZzlHZXpgXMdY7jAQtZRWNWNIHLC1nMK0iRiduZU+M6kXOVvT+ze4CwdSljlViT0hxzBv5zfvFrwfEKGK0DGF3TF6MQz4HJD4QIyCfDcgrARKr1RKIl6KPgDRddWZZQyPUElYR4wCQ1b/0bVM3SyDzSuxzPp32Zg0JINR7wJZWDkir8jciDIrxjBykl60b6cfHIO1zCoi874QUhbiIW3LQ0dDqzqwhB4g28DgO0NY2KuLmQHgHbtsJtRMIMPAaDtg5JOYO5fSG3JNM6pnHuD7MimXTKWR3c4+KuDnhpKj7LdljqQpzcJlFjFvoEyRCGgkCcRH3k/w6XUv7yAu4q5DHSYxEQaube6L2mXOZUjLpUhjx8jqcJSaKGMkW1f3Rto9x98ROIIB8bSetZSti/FLxDeSTsJ0NfUAgZ9me+C7iYN3+grxIjvuna4/gaPCeMLUuYizjAo8LmzapS2rlqGX6gIhCN5Lkj0z9qlApPaTrHPGiwaOXqA+rZRfnE/CHO5kOXSc7FZbfWSCBE5dzsEaNhi6W971P18pSy7AXjSQBKWI8E8oxrCg4KRtzh/brxB6dCYmmpOgc8zm0Vj0iq3o0JOyRgRmZQiCE0yW0Vu0IFU/diizEgT2SANLCCsxa7QYy5232AamVqy3awDJrFQ2qy5GG1iXdofvOkcyM1Itaq7j20gMCWpKxM5X3VjBnFgxBvVcSSIJy7ML8NrI9r5TdXMvBula7z0TfU6sVHYpYgnIsSP2jp2mJTFQFmmu/zVNbGmjxajvdcii/shKUYz6zIy9irCbNZrlwe0QywGw+fYpj2DhJNFGZWJ8PvYvosJLapwWxMKy4Fs15vWp525D9ngKS7wQ3Sf4eabNzXZXyeq1z+XEbNeBjyjFB8uIOF9Or7+8nvCgs+ZbLCZ9UKj5SkkCKqRsV+RNLZ/nqk/wzSVDnrBUBqVQYbs4AaYnl0bSr8vMndaPqNil3XadgktQTYj5EbNIHYp63q4fwPeRq1JxFT/JGq9tk8/0SW8JlrMn4jVJAzHC0/N71JlK6tTe+Jk21tVQ5Vi0gOlcSPZsBsdhmVfoZ/hojzbgLhid9pxpFZWtJ7S4+jrFvNUkmmmuPp5PvQcSqHjkmEEV2ZNeyCSC+rtVGMsAhgNxH5Uneggxz4w74GeeZddfhD8rFH8JYTcdFWRG+JBDBDi8MutWa+KRtgsSERBGhMHrWcdEcAAlEJlb1AEnF6Hg8c9eEeEBqYlUteoBLC71Gkmsh4KDrupylpcmn0pud1ek4i4bL851ArM3Oj3EAxGOWTGYWgGeMln7b4hfY9vJANBdGOo0IeBqnDhuTTa51IKpGWY6KIhdGOo1oMc1uItd+qAhK0yqWMnAw3mIl2Eyk+XwoZSxnq9U1EMkiVYGArcwtqRKJ1WpIR2LXJGZguI5+3o6uL6HG83O7bBJhWK0TpLte9GWoiSK8mjHSCUTKngqEnLrX6m1gFNa4IuS0K2fw7t0YohZRW7YAgYyl3PmwOdZ1sq5VPpkCEuu6+M8zyV0zgTzT5zT7FXLrnE576KDiGbF8AdOW31mSJ3wg0kHHcJHBeRvNUGvJsDaOyxUiESs9AMJdpsbo2GWS4RVKrbK6z+CnjPLG3SPoxuhtSpnk10fnXg/PEY+gh7MjUiAfCSAYncmtAV6N/zkA4no3UVgBpUzn5XsHnbV8bz/CIj3V+nb1p/xAJO/Sozyk019Baos7XPYhrC6ta62r3MhkM6dEVpyaOLqSjWy1KAGk5JyluHfLyM2ZFqb96V++mxYBAbkWXc1VU7h5z6okEEPjtoFw46O7J2mCzLC/HVrA/Q+B8I5uc5rDnSCJ9PGeBiMVyFsvEK6z9z8zk5bScdrsAEJ1dtgXY+81cP4PAaE6+wQ7y+4FkiprlEC698iJ1qKaVct7rgEdeRBIY3WJUyBGgAzU1RDraYQobIJ85nuvkZKtjlosDOTWnxN4ChJTppEiutmPJMdA3jiQzDOjk0SMTlZstdMt+CoEYhOruDgu3tkuydBzQEzPzeOAGEni+4CU2uKEmjkORBjwER3qGKt7NDwHAJk4H0q4CIb2SEoZdYBcYiByqEyO73oNlv7OSUt1VPxuVWVVLMaCUlbMJoFk6fAiIF9GM8dSx32rjuHo6GInu9sYY8g/D0tcYxWF2qAFCHPFop/kdK27EHnuj5lAIKqxN5gQrl8TP4rd+JvxmJpDnqW1EEsOA0ENPtQXqI2UsByyaQU5VmvIG4+e+VGeCY258k417StjAqUSPao/xnebYCZ/8D0/SVJH8DJAUqk3OTUR91bAg+c5Hwjd4mZapIC0ZCjzJadYaKzvVxxLAU8jYk1k93kTUW7TlmOJ3HEPtJyVTpNIA6aa7ZmGg07nHRAWB/oC08xDuNdr4Gm0WQPxt0JWQEy13YDIlBZ2jrTnq2B01DpsJjyOphbjND/FSExtTScZtRQQerIrj/60hgy9NfNWs2BtAimS+l4FuUleMTkRcAFkbWLMY5NtQjrCIbK/Fw8/tH/WtCpytNosHE66CwCydTHmCedhqDMA0vQx9jJm89PuCwBR6vAKJJMxSoGI1nF0SBoSo2ndfiDL5E5bTx66kxn9UuaajR54LPhr9gztv4DUQukEUxJIe0cLCMu8N1ttjgAB3N2TJqPjfcSMi6YEGEB4ywtswI4Bge0cZQpTbo+wpks+kG1usUthCIjWtSoSDjiRxMTUdmOz84P0SCDJLoiJEHcT28+zeY7MO4BM7kgqe8QD4k/IlartzIo0gXTtEf8VcrpdwgMiTGSjubgE0iG1Am9SOiYWVomTfleLptPoU2n/FRNIdI6EJ0BHJN+95NTiqRZac8fJHvYqUn6twUsS5eOpFt2pOnQtixKWDuVLT2aEPR5ic1t2DR3XDu1XOyn0nGXd5D4QKW5x1z5eHNhhjyyn1KvrWxptistvog5A7NXnGYN5C9Frmy1uvgcG8mLUzpaCBIk9LGezz2c7w/foC2mfOCLCy78yXpR7PccfBZLpYJnxa/0AOaeW1sOAJNrRZDyNdyArDcHrnwKS+eFFdrizgXRaqcNAxGbP/jL2xs9z5kA8CEhPn3r5mon4SK28eTQQ0Xyopzzg58fsf6b2mzFSh0pYyk/vv3Qccr2X08cqOO7mwP0ZDcFcwyrK/jgQSMIkCh3SHo4iSbD9cSSQqNBxDxclESVZU3MHkOhl4qCNMwZVb9jN+7oXCOFr7bZIZGKuDkek7xSK38TL1CDjrXdK1vLGd071MUAbw6XmKJC+xGCqPQkC+BqO650QUdU0BoQG4hMHMvPJ8BLgcVNT1KIMAqHJyGFiMK/QreUE74UUdMjUlOJ/DEgfpTxXQVXB+ZipKQs9xoEoz5I5diUv8bnQrbwJgTv1C151jo0BsdsugMijIGEAQAauSVbqDHrjrfzEWdupkiDhGCDqQD4WCIrO81qhw4DIMqdBIG8whxfmS0jujYjRNAtkYwG/7QXyCYBw/vz2SKZNSfE7dOkisCOBwJwiVRgoDsQ/DATuEZjlpYHMgB3laCDJG2MgmJxQ83swpXGY4cUDEjvG2suB4BfnGH9u35VAOCfGmAe50f5pIK6qwMcNHoiYwFMD4Q7Z2D9qABFteFWxGL6nwAh1LVxiBqhjqNesp0ibvZAUkLDqDc0Hw8hiRiX4hfO3ARAaoezOhCpvJNtiNCAllwTdc0thILfRPR8N3mNI5tOGbItF9QMpZ6+OjzghKBXk5TEjBKRpkoBeqYmVSflws9dK93sFavxk10jpIC8o5rVIzl16sDB87E4Jfd8UEO3TnnR5tVUHNHsn+TgQyZ6aA6KjDLrg3QICec7aQ4cbnAh5NwpEGTZ2ZZZmt5NAxtz+/CgVUgtvdhR2V+QIZkeDmZu64ta9+ZvWxYBcjLGB+QNyUZiNACYv/nUX7CPniAPEORDxm4iv240AnEU7fLJ7QOw0MSMSN7FvO40AvKSpQV3LA2LHRDhhmHU5DQUcFXdQ+3WBnKxb5kKKPg+/eY3ZIz4Q+x0zsqX2OvlKOYQPvrJAElvyIdwdBwNJ0iIPsKn8NpCUFUw6tf36yhqOIWK0tXfer0/IoUAIB9Sjk1geDKSTlesvBrLDvPi7gBynOv1pID1cgn85kD91/R/I33b9D2BqKIQldvOpAAAAJXRFWHRkYXRlOmNyZWF0ZQAyMDE4LTExLTE0VDE1OjQwOjE1LTA2OjAwf3H3CgAAACV0RVh0ZGF0ZTptb2RpZnkAMjAxOC0xMS0xNFQxNTo0MDoxNS0wNjowMA4sT7YAAAAASUVORK5CYII=\n", "text/plain": [ "" ] @@ -356,7 +363,7 @@ "metadata": { "anaconda-cloud": {}, "kernelspec": { - "display_name": "Python [default]", + "display_name": "Python 3", "language": "python", "name": "python3" }, @@ -370,9 +377,9 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.5.2" + "version": "3.6.5" } }, "nbformat": 4, - "nbformat_minor": 0 + "nbformat_minor": 1 } diff --git a/openmc/model/triso.py b/openmc/model/triso.py index 1e4b4aa9b..e04d7f589 100644 --- a/openmc/model/triso.py +++ b/openmc/model/triso.py @@ -2,19 +2,23 @@ import copy import warnings import itertools import random -from collections import defaultdict +from abc import ABCMeta, abstractproperty, abstractmethod +from collections import Counter, defaultdict from collections.abc import Iterable -from numbers import Real -from random import uniform, gauss from heapq import heappush, heappop from math import pi, sin, cos, floor, log10, sqrt -from abc import ABCMeta, abstractproperty, abstractmethod +from numbers import Real +from random import uniform, gauss import numpy as np import scipy.spatial import openmc -import openmc.checkvalue as cv +from openmc.checkvalue import check_type + + +MAX_PF_RSP = 0.38 +MAX_PF_CRP = 0.64 class TRISO(openmc.Cell): @@ -55,7 +59,7 @@ class TRISO(openmc.Cell): @center.setter def center(self, center): - cv.check_type('TRISO center', center, Iterable, Real) + check_type('TRISO center', center, Iterable, Real) self._surface.x0 = center[0] self._surface.y0 = center[1] self._surface.z0 = center[2] @@ -91,43 +95,43 @@ class TRISO(openmc.Cell): k_min:k_max+1, j_min:j_max+1, i_min:i_max+1])) -class _Domain(metaclass=ABCMeta): - """Container in which to pack particles. +class _Container(metaclass=ABCMeta): + """Container in which to pack spheres. Parameters ---------- - particle_radius : float - Radius of particles to be packed in container. + sphere_radius : float + Radius of spheres to be packed in container. center : Iterable of float Cartesian coordinates of the center of the container. Default is - [0., 0., 0.] + (0., 0., 0.) Attributes ---------- - particle_radius : float - Radius of particles to be packed in container. + sphere_radius : float + Radius of spheres to be packed in container. center : list of float Cartesian coordinates of the center of the container. Default is - [0., 0., 0.] + (0., 0., 0.) cell_length : list of float Length in x-, y-, and z- directions of each cell in mesh overlaid on domain. limits : list of float - Constraint on where particle center can be placed. + Constraint on where sphere center can be placed. volume : float Volume of the container. """ - def __init__(self, particle_radius, center=[0., 0., 0.]): + def __init__(self, sphere_radius, center=(0., 0., 0.)): self._cell_length = None self._limits = None - self.particle_radius = particle_radius + self.sphere_radius = sphere_radius self.center = center @property - def particle_radius(self): - return self._particle_radius + def sphere_radius(self): + return self._sphere_radius @property def center(self): @@ -145,27 +149,24 @@ class _Domain(metaclass=ABCMeta): def volume(self): pass - @particle_radius.setter - def particle_radius(self, particle_radius): - self._particle_radius = float(particle_radius) + @sphere_radius.setter + def sphere_radius(self, sphere_radius): + self._sphere_radius = float(sphere_radius) self._limits = None self._cell_length = None @center.setter def center(self, center): - if np.asarray(center).size != 3: - raise ValueError('Unable to set domain center to {} since it must ' - 'be of length 3'.format(center)) - self._center = [float(x) for x in center] + self._center = center def mesh_cell(self, p): """Calculate the index of the cell in a mesh overlaid on the domain in - which the given particle center falls. + which the given sphere center falls. Parameters ---------- p : Iterable of float - Cartesian coordinates of particle center. + Cartesian coordinates of sphere center. Returns ------- @@ -177,12 +178,12 @@ class _Domain(metaclass=ABCMeta): def nearby_mesh_cells(self, p): """Calculates the indices of all cells in a mesh overlaid on the domain - within one diameter of the given particle. + within one diameter of the given sphere. Parameters ---------- p : Iterable of float - Cartesian coordinates of particle center. + Cartesian coordinates of sphere center. Returns ------- @@ -190,27 +191,41 @@ class _Domain(metaclass=ABCMeta): Indices of mesh cells. """ - d = 2*self.particle_radius + d = 2*self.sphere_radius r = [[a/self.cell_length[i] for a in [p[i]-d, p[i], p[i]+d]] for i in range(3)] return list(itertools.product(*({int(x) for x in y} for y in r))) @abstractmethod - def random_point(self): - """Generate Cartesian coordinates of center of a particle that is - contained entirely within the domain with uniform probability. + def from_region(self, region, sphere_radius): + """Create a container to pack spheres in based on a region. - Returns - ------- - list of float - Cartesian coordinates of particle center. + Parameters + ---------- + region : openmc.Region + Region to create container from. + sphere_radius : float + Outer radius of spheres. """ pass @abstractmethod - def repel_particles(self, p, q, d, d_new): - """Move particles p and q apart according to the following + def random_point(self): + """Generate Cartesian coordinates of center of a sphere that is + contained entirely within the domain with uniform probability. + + Returns + ------- + list of float + Cartesian coordinates of sphere center. + + """ + pass + + @abstractmethod + def repel_spheres(self, p, q, d, d_new): + """Move spheres p and q apart according to the following transformation (accounting for boundary conditions on domain): r_i^(n+1) = r_i^(n) + 1/2(d_out^(n+1) - d^(n)) @@ -219,78 +234,112 @@ class _Domain(metaclass=ABCMeta): Parameters ---------- p, q : numpy.ndarray - Cartesian coordinates of particle center. + Cartesian coordinates of sphere center. d : float - distance between centers of particles i and j. + distance between centers of spheres i and j. d_new : float - final distance between centers of particles i and j. + final distance between centers of spheres i and j. """ pass -class _CubicDomain(_Domain): - """Cubic container in which to pack particles. +class _RectangularPrism(_Container): + """Rectangular prism container in which to pack spheres. Parameters ---------- - length : float - Length of each side of the cubic container. - particle_radius : float - Radius of particles to be packed in container. + width : float + Prism length along the x-axis + depth : float + Prism length along the y-axis + height : float + Prism length along the z-axis + sphere_radius : float + Radius of spheres to be packed in container. center : Iterable of float Cartesian coordinates of the center of the container. Default is - [0., 0., 0.] + (0., 0., 0.) Attributes ---------- - length : float - Length of each side of the cubic container. - particle_radius : float - Radius of particles to be packed in container. + width : float + Prism length along the x-axis + depth : float + Prism length along the y-axis + height : float + Prism length along the z-axis + sphere_radius : float + Radius of spheres to be packed in container. center : list of float Cartesian coordinates of the center of the container. Default is - [0., 0., 0.] + (0., 0., 0.) cell_length : list of float Length in x-, y-, and z- directions of each cell in mesh overlaid on domain. limits : list of float - Maximum distance from center in x-, y-, or z-direction where particle + Minimum and maximum distance in x-, y-, and z-direction where sphere center can be placed. volume : float Volume of the container. """ - def __init__(self, length, particle_radius, center=[0., 0., 0.]): - super().__init__(particle_radius, center) - self.length = length + def __init__(self, width, depth, height, sphere_radius, center=(0., 0., 0.)): + super().__init__(sphere_radius, center) + self.width = width + self.depth = depth + self.height = height @property - def length(self): - return self._length + def width(self): + return self._width + + @property + def depth(self): + return self._depth + + @property + def height(self): + return self._height @property def limits(self): if self._limits is None: - self._limits = [self.length/2 - self.particle_radius] + c = self.center + r = self.sphere_radius + x, y, z = self.width/2, self.depth/2, self.height/2 + self._limits = [[c[0] - x + r, c[1] - y + r, c[2] - z + r], + [c[0] + x - r, c[1] + y - r, c[2] + z - r]] return self._limits @property def cell_length(self): if self._cell_length is None: - mesh_length = [self.length, self.length, self.length] - self._cell_length = [x/int(x/(4*self.particle_radius)) + mesh_length = [self.width, self.depth, self.height] + self._cell_length = [x/int(x/(4*self.sphere_radius)) for x in mesh_length] return self._cell_length @property def volume(self): - return self.length**3 + return self.width*self.depth*self.height - @length.setter - def length(self, length): - self._length = float(length) + @width.setter + def width(self, width): + self._width = float(width) + self._limits = None + self._cell_length = None + + @depth.setter + def depth(self, depth): + self._depth = float(depth) + self._limits = None + self._cell_length = None + + @height.setter + def height(self, height): + self._height = float(height) self._limits = None self._cell_length = None @@ -298,15 +347,54 @@ class _CubicDomain(_Domain): def limits(self, limits): self._limits = limits - def random_point(self): - x_max = self.limits[0] - return [uniform(-x_max, x_max), - uniform(-x_max, x_max), - uniform(-x_max, x_max)] + @classmethod + def from_region(self, region, sphere_radius): + check_type('region', region, openmc.Region) - def repel_particles(self, p, q, d, d_new): - # Moving each particle distance 's' away from the other along the line - # joining the particle centers will ensure their final distance is + # Assume the simplest case where the prism volume is the intersection + # of the half-spaces of six planes + if not isinstance(region, openmc.Intersection): + raise ValueError + + if any(not isinstance(node, openmc.Halfspace) for node in region): + raise ValueError + + if len(region) != 6: + raise ValueError + + # Sort half-spaces by surface type + px1, px2, py1, py2, pz1, pz2 = sorted(region, key=lambda x: x.surface.type) + + # Make sure the region consists of the correct surfaces + if (not isinstance(px1.surface, openmc.XPlane) or + not isinstance(px2.surface, openmc.XPlane) or + not isinstance(py1.surface, openmc.YPlane) or + not isinstance(py2.surface, openmc.YPlane) or + not isinstance(pz1.surface, openmc.ZPlane) or + not isinstance(pz2.surface, openmc.ZPlane)): + raise ValueError + + # Make sure the half-spaces are on the correct side of the surfaces + ll, ur = region.bounding_box + if any(x in ll or x in ur for x in (-np.inf, np.inf)): + raise ValueError + + # Calculate the parameters for the container + width, depth, height = ur - ll + center = ll + [width/2, depth/2, height/2] + + # The region is the volume of a rectangular prism, so create container + return _RectangularPrism(width, depth, height, sphere_radius, center) + + def random_point(self): + ll, ul = self.limits + return [uniform(ll[0], ul[0]), + uniform(ll[1], ul[1]), + uniform(ll[2], ul[2])] + + def repel_spheres(self, p, q, d, d_new): + # Moving each sphere distance 's' away from the other along the line + # joining the sphere centers will ensure their final distance is # equal to the outer diameter s = (d_new - d)/2 @@ -314,53 +402,64 @@ class _CubicDomain(_Domain): p += s*v q -= s*v - # Enforce the rigid boundary by moving each particle back along the + # Enforce the rigid boundary by moving each sphere back along the # surface normal until it is completely within the container if it # overlaps the surface - x_max = self.limits[0] - p[:] = np.clip(p, -x_max, x_max) - q[:] = np.clip(q, -x_max, x_max) + p[:] = np.clip(p, self.limits[0], self.limits[1]) + q[:] = np.clip(q, self.limits[0], self.limits[1]) -class _CylindricalDomain(_Domain): - """Cylindrical container in which to pack particles. +class _Cylinder(_Container): + """Cylindrical container in which to pack spheres. Parameters ---------- length : float - Length along z-axis of the cylindrical container. + Length of the cylindrical container. radius : float Radius of the cylindrical container. + axis : string + Axis along which the length of the cylinder is aligned. + sphere_radius : float + Radius of spheres to be packed in container. center : Iterable of float Cartesian coordinates of the center of the container. Default is - [0., 0., 0.] + (0., 0., 0.) Attributes ---------- length : float - Length along z-axis of the cylindrical container. + Length of the cylindrical container. radius : float Radius of the cylindrical container. - particle_radius : float - Radius of particles to be packed in container. + axis : string + Axis along which the length of the cylinder is aligned. + sphere_radius : float + Radius of spheres to be packed in container. center : list of float Cartesian coordinates of the center of the container. Default is - [0., 0., 0.] + (0., 0., 0.) + shift : list of int + Rolled indices of the x-, y-, and z- coordinates of a sphere so the + configuration is aligned with the correct axis. No shift corresponds to + a cylinder along the z-axis. cell_length : list of float Length in x-, y-, and z- directions of each cell in mesh overlaid on domain. limits : list of float - Maximum radial distance and maximum distance from center in z-direction - where particle center can be placed. + Maximum radial distance and minimum and maximum distance in the + direction parallel to the axis where sphere center can be placed. volume : float Volume of the container. """ - def __init__(self, length, radius, particle_radius, center=[0., 0., 0.]): - super().__init__(particle_radius, center) + def __init__(self, length, radius, axis, sphere_radius, center=(0., 0., 0.)): + super().__init__(sphere_radius, center) + self._shift = None self.length = length self.radius = radius + self.axis = axis @property def length(self): @@ -370,24 +469,44 @@ class _CylindricalDomain(_Domain): def radius(self): return self._radius + @property + def axis(self): + return self._axis + + @property + def shift(self): + if self._shift is None: + if self.axis == 'x': + self._shift = [1, 2, 0] + elif self.axis == 'y': + self._shift = [2, 0, 1] + else: + self._shift = [0, 1, 2] + return self._shift + @property def limits(self): if self._limits is None: - self._limits = [self.radius - self.particle_radius, - self.length/2 - self.particle_radius] + z0 = self.center[self.shift[2]] + z = self.length/2 + r = self.sphere_radius + self._limits = [[z0 - z + r], [z0 + z - r, self.radius - r]] return self._limits @property def cell_length(self): if self._cell_length is None: - mesh_length = [2*self.radius, 2*self.radius, self.length] - self._cell_length = [x/int(x/(4*self.particle_radius)) - for x in mesh_length] + h = 4*self.sphere_radius + i, j, k = self.shift + self._cell_length = [None]*3 + self._cell_length[i] = 2*self.radius/int(2*self.radius/h) + self._cell_length[j] = 2*self.radius/int(2*self.radius/h) + self._cell_length[k] = self.length/int(self.length/h) return self._cell_length @property def volume(self): - return self.length * pi * self.radius**2 + return self.length*pi*self.radius**2 @length.setter def length(self, length): @@ -401,20 +520,86 @@ class _CylindricalDomain(_Domain): self._limits = None self._cell_length = None + @axis.setter + def axis(self, axis): + self._axis = axis + self._shift = None + @limits.setter def limits(self, limits): self._limits = limits + @classmethod + def from_region(self, region, sphere_radius): + check_type('region', region, openmc.Region) + + # Assume the simplest case where the cylinder volume is the + # intersection of the half-spaces of a cylinder and two planes + if not isinstance(region, openmc.Intersection): + raise ValueError + + if any(not isinstance(node, openmc.Halfspace) for node in region): + raise ValueError + + if len(region) != 3: + raise ValueError + + # Identify the axis that the cylinder lies along + axis = region[0].surface.type[0] + + # Make sure the region is composed of a cylinder and two planes on the + # same axis + count = Counter(node.surface.type for node in region) + if count[axis + '-cylinder'] != 1 or count[axis + '-plane'] != 2: + raise ValueError + + # Sort the half-spaces by surface type + cyl, p1, p2 = sorted(region, key=lambda x: x.surface.type) + + # Calculate the parameters for a cylinder along the x-axis + if axis == 'x': + if p1.surface.x0 > p2.surface.x0: + p1, p2 = p2, p1 + length = p2.surface.x0 - p1.surface.x0 + center = (p1.surface.x0 + length/2, cyl.surface.y0, cyl.surface.z0) + + # Calculate the parameters for a cylinder along the y-axis + elif axis == 'y': + if p1.surface.y0 > p2.surface.y0: + p1, p2 = p2, p1 + length = p2.surface.y0 - p1.surface.y0 + center = (cyl.surface.x0, p1.surface.y0 + length/2, cyl.surface.z0) + + # Calculate the parameters for a cylinder along the z-axis + else: + if p1.surface.z0 > p2.surface.z0: + p1, p2 = p2, p1 + length = p2.surface.z0 - p1.surface.z0 + center = (cyl.surface.x0, cyl.surface.y0, p1.surface.z0 + length/2) + + # Make sure the half-spaces are on the correct side of the surfaces + if cyl.side != '-' or p1.side != '+' or p2.side != '-': + raise ValueError + + radius = cyl.surface.r + + # The region is the volume of a cylinder, so create container + return _Cylinder(length, radius, axis, sphere_radius, center) + def random_point(self): - r_max = self.limits[0] - z_max = self.limits[1] - r = sqrt(uniform(0, r_max**2)) + ll, ul = self.limits + r = sqrt(uniform(0, ul[1]**2)) t = uniform(0, 2*pi) - return [r*cos(t), r*sin(t), uniform(-z_max, z_max)] + i, j, k = self.shift + p = [None]*3 + p[i] = r*cos(t) + self.center[i] + p[j] = r*sin(t) + self.center[j] + p[k] = uniform(ll[0], ul[0]) + return p - def repel_particles(self, p, q, d, d_new): - # Moving each particle distance 's' away from the other along the line - # joining the particle centers will ensure their final distance is + def repel_spheres(self, p, q, d, d_new): + # Moving each sphere distance 's' away from the other along the line + # joining the sphere centers will ensure their final distance is # equal to the outer diameter s = (d_new - d)/2 @@ -422,78 +607,96 @@ class _CylindricalDomain(_Domain): p += s*v q -= s*v - # Enforce the rigid boundary by moving each particle back along the + # Enforce the rigid boundary by moving each sphere back along the # surface normal until it is completely within the container if it # overlaps the surface - r_max = self.limits[0] - z_max = self.limits[1] + ll, ul = self.limits + c = self.center + i, j, k = self.shift - r = sqrt(p[0]**2 + p[1]**2) - if r > r_max: - p[0:2] *= r_max/r - p[2] = np.clip(p[2], -z_max, z_max) + r = sqrt((p[i] - c[i])**2 + (p[j] - c[j])**2) + if r > ul[1]: + p[i] = (p[i] - c[i])*ul[1]/r + c[i] + p[j] = (p[j] - c[j])*ul[1]/r + c[j] + p[k] = np.clip(p[k], ll[0], ul[0]) - r = sqrt(q[0]**2 + q[1]**2) - if r > r_max: - q[0:2] *= r_max/r - q[2] = np.clip(q[2], -z_max, z_max) + r = sqrt((q[i] - c[i])**2 + (q[j] - c[j])**2) + if r > ul[1]: + q[i] = (q[i] - c[i])*ul[1]/r + c[i] + q[j] = (q[j] - c[j])*ul[1]/r + c[j] + q[k] = np.clip(q[k], ll[0], ul[0]) -class _SphericalDomain(_Domain): - """Spherical container in which to pack particles. +class _SphericalShell(_Container): + """Spherical shell container in which to pack spheres. Parameters ---------- radius : float - Radius of the spherical container. + Outer radius of the spherical shell container. + inner_radius : float + Inner radius of the spherical shell container. center : Iterable of float Cartesian coordinates of the center of the container. Default is - [0., 0., 0.] + (0., 0., 0.) Attributes ---------- radius : float - Radius of the spherical container. - particle_radius : float - Radius of particles to be packed in container. + Outer radius of the spherical shell container. + inner_radius : float + Inner radius of the spherical shell container. + sphere_radius : float + Radius of spheres to be packed in container. center : list of float Cartesian coordinates of the center of the container. Default is - [0., 0., 0.] + (0., 0., 0.) cell_length : list of float Length in x-, y-, and z- directions of each cell in mesh overlaid on domain. limits : list of float - Maximum radial distance where particle center can be placed. + Minimum and maximum radial distance where sphere center can be placed. volume : float Volume of the container. """ - def __init__(self, radius, particle_radius, center=[0., 0., 0.]): - super().__init__(particle_radius, center) + def __init__(self, radius, inner_radius, sphere_radius, + center=(0., 0., 0.)): + super().__init__(sphere_radius, center) self.radius = radius + self.inner_radius = inner_radius @property def radius(self): return self._radius + @property + def inner_radius(self): + return self._inner_radius + @property def limits(self): if self._limits is None: - self._limits = [self.radius - self.particle_radius] + r_max = self.radius - self.sphere_radius + if self.inner_radius == 0: + r_min = 0 + else: + r_min = self.inner_radius + self.sphere_radius + self._limits = [[r_min], [r_max]] return self._limits @property def cell_length(self): if self._cell_length is None: - mesh_length = [2*self.radius, 2*self.radius, 2*self.radius] - self._cell_length = [x/int(x/(4*self.particle_radius)) + mesh_length = 3*[2*self.radius] + self._cell_length = [x/int(x/(4*self.sphere_radius)) for x in mesh_length] return self._cell_length @property def volume(self): - return 4/3 * pi * self.radius**3 + return 4/3*pi*(self.radius**3 - self.inner_radius**3) @radius.setter def radius(self, radius): @@ -501,19 +704,71 @@ class _SphericalDomain(_Domain): self._limits = None self._cell_length = None + @inner_radius.setter + def inner_radius(self, inner_radius): + self._inner_radius = float(inner_radius) + self._limits = None + @limits.setter def limits(self, limits): self._limits = limits - def random_point(self): - r_max = self.limits[0] - x = (gauss(0, 1), gauss(0, 1), gauss(0, 1)) - r = (uniform(0, r_max**3)**(1/3) / sqrt(x[0]**2 + x[1]**2 + x[2]**2)) - return [r*s for s in x] + @classmethod + def from_region(self, region, sphere_radius): + check_type('region', region, openmc.Region) - def repel_particles(self, p, q, d, d_new): - # Moving each particle distance 's' away from the other along the line - # joining the particle centers will ensure their final distance is + # First check if the region is the volume inside a sphere. Assume the + # simplest case where the sphere volume is the negative half-space of a + # sphere. + if (isinstance(region, openmc.Halfspace) + and isinstance(region.surface, openmc.Sphere) + and region.side == '-'): + + # The region is the volume of a sphere, so create container + radius = region.surface.r + center = (region.surface.x0, region.surface.y0, region.surface.z0) + + return _SphericalShell(radius, 0., sphere_radius, center) + + # Next check for a spherical shell volume. Assume the simplest case + # where the spherical shell volume is the intersection of the + # half-spaces of two spheres. + if not isinstance(region, openmc.Intersection): + raise ValueError + + if any(not isinstance(node, openmc.Halfspace) for node in region): + raise ValueError + + if len(region) != 2: + raise ValueError + + if any(not isinstance(node.surface, openmc.Sphere) for node in region): + raise ValueError + + s1, s2 = sorted(region, key=lambda x: x.surface.r) + radius = s2.surface.r + inner_radius = s1.surface.r + center = (s1.surface.x0, s1.surface.y0, s1.surface.z0) + + if center != (s2.surface.x0, s2.surface.y0, s2.surface.z0): + raise ValueError + + if s1.side != '+' or s2.side != '-': + raise ValueError + + # The region is the volume of a spherical shell, so create container + return _SphericalShell(radius, inner_radius, sphere_radius, center) + + def random_point(self): + c = self.center + ll, ul = self.limits + x, y, z = (gauss(0, 1), gauss(0, 1), gauss(0, 1)) + r = (uniform(ll[0]**3, ul[0]**3)**(1/3)/sqrt(x**2 + y**2 + z**2)) + return [r*x + c[0], r*y + c[1], r*z + c[2]] + + def repel_spheres(self, p, q, d, d_new): + # Moving each sphere distance 's' away from the other along the line + # joining the sphere centers will ensure their final distance is # equal to the outer diameter s = (d_new - d)/2 @@ -521,18 +776,23 @@ class _SphericalDomain(_Domain): p += s*v q -= s*v - # Enforce the rigid boundary by moving each particle back along the + # Enforce the rigid boundary by moving each sphere back along the # surface normal until it is completely within the container if it # overlaps the surface - r_max = self.limits[0] + c = self.center + ll, ul = self.limits - r = sqrt(p[0]**2 + p[1]**2 + p[2]**2) - if r > r_max: - p *= r_max/r + r = sqrt((p[0] - c[0])**2 + (p[1] - c[1])**2 + (p[2] - c[2])**2) + if r > ul[0]: + p[:] = (p - c)*ul[0]/r + c + elif r < ll[0]: + p[:] = (p - c)*ll[0]/r + c - r = sqrt(q[0]**2 + q[1]**2 + q[2]**2) - if r > r_max: - q *= r_max/r + r = sqrt((q[0] - c[0])**2 + (q[1] - c[1])**2 + (q[2] - c[2])**2) + if r > ul[0]: + q[:] = (q - c)*ul[0]/r + c + elif r < ll[0]: + q[:] = (q - c)*ll[0]/r + c def create_triso_lattice(trisos, lower_left, pitch, shape, background): @@ -608,28 +868,28 @@ def create_triso_lattice(trisos, lower_left, pitch, shape, background): return lattice -def _random_sequential_pack(domain, n_particles): - """Random sequential packing of particles within a container. +def _random_sequential_pack(domain, num_spheres): + """Random sequential packing of spheres within a container. Parameters ---------- - domain : openmc.model._Domain - Container in which to pack particles. - n_particles : int - Number of particles to pack. + domain : openmc.model._Container + Container in which to pack spheres. + num_spheres : int + Number of spheres to pack. Returns ------ numpy.ndarray - Cartesian coordinates of centers of particles. + Cartesian coordinates of centers of spheres. """ - sqd = (2*domain.particle_radius)**2 - particles = [] + sqd = (2*domain.sphere_radius)**2 + spheres = [] mesh = defaultdict(list) - for i in range(n_particles): + for i in range(num_spheres): # Randomly sample new center coordinates while there are any overlaps while True: p = domain.random_point() @@ -639,23 +899,23 @@ def _random_sequential_pack(domain, n_particles): continue else: break - particles.append(p) + spheres.append(p) for idx in domain.nearby_mesh_cells(p): mesh[idx].append(p) - return np.array(particles) + return np.array(spheres) -def _close_random_pack(domain, particles, contraction_rate): - """Close random packing of particles using the Jodrey-Tory algorithm. +def _close_random_pack(domain, spheres, contraction_rate): + """Close random packing of spheres using the Jodrey-Tory algorithm. Parameters ---------- - domain : openmc.model._Domain - Container in which to pack particles. - particles : numpy.ndarray - Initial Cartesian coordinates of centers of particles. + domain : openmc.model._Container + Container in which to pack spheres. + spheres : numpy.ndarray + Initial Cartesian coordinates of centers of spheres. contraction_rate : float Contraction rate of outer diameter. @@ -667,9 +927,9 @@ def _close_random_pack(domain, particles, contraction_rate): Parameters ---------- d : float - distance between centers of particles i and j. + distance between centers of spheres i and j. i, j : int - Index of particles in particles array. + Index of spheres in spheres array. """ @@ -679,12 +939,12 @@ def _close_random_pack(domain, particles, contraction_rate): heappush(rods, rod) def remove_rod(i): - """Mark the rod containing particle i as removed. + """Mark the rod containing sphere i as removed. Parameters ---------- i : int - Index of particle in particles array. + Index of sphere in spheres array. """ @@ -700,9 +960,9 @@ def _close_random_pack(domain, particles, contraction_rate): Returns ------- d : float - distance between centers of particles i and j. + distance between centers of spheres i and j. i, j : int - Index of particles in particles array. + Index of spheres in spheres array. """ @@ -715,15 +975,15 @@ def _close_random_pack(domain, particles, contraction_rate): return None, None, None def create_rod_list(): - """Generate sorted list of rods (distances between particle centers). + """Generate sorted list of rods (distances between sphere centers). Rods are arranged in a heap where each element contains the rod length - and the particle indices. A rod between particles p and q is only + and the sphere indices. A rod between spheres p and q is only included if the distance between p and q could not be changed by the elimination of a greater overlap, i.e. q has no nearer neighbors than p. - A mapping of particle ids to rods is maintained in 'rods_map'. Each key - in the dict is the id of a particle that is in the rod list, and the + A mapping of sphere ids to rods is maintained in 'rods_map'. Each key + in the dict is the id of a sphere that is in the rod list, and the value is the id of its nearest neighbor and the rod that contains them. The dict is used to find rods in the priority queue and to mark removed rods so rods can be "removed" without breaking the heap structure @@ -732,25 +992,25 @@ def _close_random_pack(domain, particles, contraction_rate): """ # Create KD tree for quick nearest neighbor search - tree = scipy.spatial.cKDTree(particles) + tree = scipy.spatial.cKDTree(spheres) # Find distance to nearest neighbor and index of nearest neighbor for - # all particles - d, n = tree.query(particles, k=2) + # all spheres + d, n = tree.query(spheres, k=2) d = d[:, 1] n = n[:, 1] - # Array of particle indices, indices of nearest neighbors, and + # Array of sphere indices, indices of nearest neighbors, and # distances to nearest neighbors a = np.vstack((list(range(n.size)), n, d)).T # Sort along second column and swap first and second columns to create - # array of nearest neighbor indices, indices of particles they are + # array of nearest neighbor indices, indices of spheres they are # nearest neighbors of, and distances between them b = a[a[:, 1].argsort()] b[:, [0, 1]] = b[:, [1, 0]] - # Find the intersection between 'a' and 'b': a list of particles who + # Find the intersection between 'a' and 'b': a list of spheres who # are each other's nearest neighbors and the distance between them r = list({tuple(x) for x in a} & {tuple(x) for x in b}) @@ -766,30 +1026,30 @@ def _close_random_pack(domain, particles, contraction_rate): add_rod(d, i, j) def update_mesh(i): - """Update which mesh cells the particle is in based on new particle + """Update which mesh cells the sphere is in based on new sphere center coordinates. 'mesh'/'mesh_map' is a two way dictionary used to look up which - particles are located within one diameter of a given mesh cell and - which mesh cells a given particle center is within one diameter of. + spheres are located within one diameter of a given mesh cell and + which mesh cells a given sphere center is within one diameter of. This is used to speed up the nearest neighbor search. Parameters ---------- i : int - Index of particle in particles array. + Index of sphere in spheres array. """ - # Determine which mesh cells the particle is in and remove the - # particle id from those cells + # Determine which mesh cells the sphere is in and remove the + # sphere id from those cells for idx in mesh_map[i]: mesh[idx].remove(i) del mesh_map[i] - # Determine which mesh cells are within one diameter of particle's - # center and add this particle to the list of particles in those cells - for idx in domain.nearby_mesh_cells(particles[i]): + # Determine which mesh cells are within one diameter of sphere's + # center and add this sphere to the list of spheres in those cells + for idx in domain.nearby_mesh_cells(spheres[i]): mesh[idx].add(i) mesh_map[i].add(idx) @@ -798,7 +1058,7 @@ def _close_random_pack(domain, particles, contraction_rate): d_out^(i+1) = d_out^(i) - (1/2)^(j) * d_out0 * k / n, - where k is the contraction rate, n is the number of particles, and + where k is the contraction rate, n is the number of spheres, and j = floor(-log10(pf_out - pf_in)). @@ -809,26 +1069,26 @@ def _close_random_pack(domain, particles, contraction_rate): """ - inner_pf = 4/3*pi*(inner_diameter/2)**3*n_particles/domain.volume - outer_pf = 4/3*pi*(outer_diameter/2)**3*n_particles/domain.volume + inner_pf = 4/3*pi*(inner_diameter/2)**3*num_spheres/domain.volume + outer_pf = 4/3*pi*(outer_diameter/2)**3*num_spheres/domain.volume j = floor(-log10(outer_pf - inner_pf)) return (outer_diameter - 0.5**j * contraction_rate * - initial_outer_diameter / n_particles) + initial_outer_diameter / num_spheres) def nearest(i): - """Find index of nearest neighbor of particle i. + """Find index of nearest neighbor of sphere i. Parameters ---------- i : int - Index in particles array of particle for which to find nearest + Index in spheres array of sphere for which to find nearest neighbor. Returns ------- int - Index in particles array of nearest neighbor of i + Index in spheres array of nearest neighbor of i float distance between i and nearest neighbor. @@ -837,8 +1097,8 @@ def _close_random_pack(domain, particles, contraction_rate): # Need the second nearest neighbor of i since the nearest neighbor # will be itself. Using argpartition, the k-th nearest neighbor is # placed at index k. - idx = list(mesh[domain.mesh_cell(particles[i])]) - dists = scipy.spatial.distance.cdist([particles[i]], particles[idx])[0] + idx = list(mesh[domain.mesh_cell(spheres[i])]) + dists = scipy.spatial.distance.cdist([spheres[i]], spheres[idx])[0] if dists.size > 1: j = dists.argpartition(1)[1] return idx[j], dists[j] @@ -846,17 +1106,17 @@ def _close_random_pack(domain, particles, contraction_rate): return None, None def update_rod_list(i): - """Update the rod list with the new nearest neighbors of particle since + """Update the rod list with the new nearest neighbors of sphere since its overlap was eliminated. Parameters ---------- i : int - Index of particle in particles array. + Index of sphere in spheres array. """ - # If the nearest neighbor k of particle i has no nearer neighbors, + # If the nearest neighbor k of sphere i has no nearer neighbors, # remove the rod currently containing k from the rod list and add rod # k-i, keeping the rod list sorted k, d_ik = nearest(i) @@ -865,57 +1125,57 @@ def _close_random_pack(domain, particles, contraction_rate): remove_rod(k) add_rod(d_ik, i, k) - n_particles = len(particles) - diameter = 2*domain.particle_radius + num_spheres = len(spheres) + diameter = 2*domain.sphere_radius # Flag for marking rods that have been removed from priority queue removed = -1 # Outer diameter initially set to arbitrary value that yields pf of 1 - initial_outer_diameter = 2*(domain.volume/(n_particles*4/3*pi))**(1/3) + initial_outer_diameter = 2*(domain.volume/(num_spheres*4/3*pi))**(1/3) - # Inner and outer diameter of particles will change during packing + # Inner and outer diameter of spheres will change during packing outer_diameter = initial_outer_diameter inner_diameter = 0. - # List of rods arranged in a heap and mapping of particle ids to rods + # List of rods arranged in a heap and mapping of sphere ids to rods rods = [] rods_map = {} - # Initialize two-way dictionary that identifies which particles are near a - # given mesh cell and which mesh cells a particle is near + # Initialize two-way dictionary that identifies which spheres are near a + # given mesh cell and which mesh cells a sphere is near mesh = defaultdict(set) mesh_map = defaultdict(set) - for i in range(n_particles): - for idx in domain.nearby_mesh_cells(particles[i]): + for i in range(num_spheres): + for idx in domain.nearby_mesh_cells(spheres[i]): mesh[idx].add(i) mesh_map[i].add(idx) while True: - # Rebuild the sorted list of rods according to the current particle + # Rebuild the sorted list of rods according to the current sphere # configuration create_rod_list() # Set the inner diameter to the shortest center-to-center distance - # between any two particles + # between any two spheres if rods: inner_diameter = rods[0][0] - # Reached the desired particle radius + # Reached the desired sphere radius if inner_diameter >= diameter: break - # The algorithm converged before reaching the desired particle radius. + # The algorithm converged before reaching the desired sphere radius. # This can happen when the desired packing fraction is close to the # packing fraction limit. The packing fraction is a random variable - # that is determined by the particle locations and the contraction + # that is determined by the sphere locations and the contraction # rate. A higher packing fraction can be achieved with a smaller # contraction rate, though at the cost of a longer simulation time -- # the number of iterations needed to remove all overlaps is inversely # proportional to the contraction rate. if inner_diameter >= outer_diameter or not rods: warnings.warn('Close random pack converged before reaching true ' - 'particle radius; some particles may overlap. Try ' + 'sphere radius; some spheres may overlap. Try ' 'reducing contraction rate or packing fraction.') break @@ -924,7 +1184,7 @@ def _close_random_pack(domain, particles, contraction_rate): if not d: break outer_diameter = reduce_outer_diameter() - domain.repel_particles(particles[i], particles[j], d, outer_diameter) + domain.repel_spheres(spheres[i], spheres[j], d, outer_diameter) update_mesh(i) update_mesh(j) update_rod_list(i) @@ -936,54 +1196,46 @@ def _close_random_pack(domain, particles, contraction_rate): break -def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None, - domain_radius=None, domain_center=[0., 0., 0.], - n_particles=None, packing_fraction=None, - initial_packing_fraction=0.3, contraction_rate=1.e-3, seed=1): - """Generate a random, non-overlapping configuration of TRISO particles - within a container. +def pack_spheres(radius, region, pf=None, num_spheres=None, initial_pf=0.3, + contraction_rate=1.e-3, seed=1): + """Generate a random, non-overlapping configuration of spheres within a + container. Parameters ---------- radius : float - Outer radius of TRISO particles. - fill : openmc.Universe - Universe which contains all layers of the TRISO particle. - domain_shape : {'cube', 'cylinder', or 'sphere'} - Geometry of the container in which the TRISO particles are packed. - domain_length : float - Length of the container (if cube or cylinder). - domain_radius : float - Radius of the container (if cylinder or sphere). - domain_center : Iterable of float - Cartesian coordinates of the center of the container. - n_particles : int - Number of TRISO particles to pack in the domain. Exactly one of - 'n_particles' and 'packing_fraction' should be specified -- the other - will be calculated. - packing_fraction : float - Packing fraction of particles. Exactly one of 'n_particles' and - 'packing_fraction' should be specified -- the other will be calculated. - initial_packing_fraction : float, optional - Packing fraction used to initialize the configuration of particles in + Outer radius of spheres. + region : openmc.Region + Container in which the spheres are packed. Supported shapes are + rectangular prism, cylinder, sphere, and spherical shell. + pf : float + Packing fraction of the spheres. One of 'pf' and 'num_spheres' must + be specified; the other will be calculated. If both are specified, 'pf' + takes precedence over 'num_spheres'. + num_spheres : int + Number of spheres to pack in the domain. One of 'num_spheres' and 'pf' + must be specified; the other will be calculated. + initial_pf : float, optional + Packing fraction used to initialize the configuration of spheres in the domain. Default value is 0.3. It is not recommended to set the initial packing fraction much higher than 0.3 as the random sequential packing algorithm becomes prohibitively slow as it approaches its limit (~0.38). contraction_rate : float, optional - Contraction rate of outer diameter. This can affect the speed of the - close random packing algorithm. Default value is 1/400. + Contraction rate of the outer diameter. Higher packing fractions can be + reached using a smaller contraction rate, but the algorithm will take + longer to converge. seed : int, optional RNG seed. Returns - ------- - trisos : list of openmc.model.TRISO - List of TRISO particles in the domain. + ------ + numpy.ndarray + Cartesian coordinates of sphere centers. Notes ----- - The particle configuration is generated using a combination of random + The sphere configuration is generated using a combination of random sequential packing (RSP) and close random packing (CRP). RSP performs better than CRP for lower packing fractions (pf), but it becomes prohibitively slow as it approaches its packing limit (~0.38). CRP can @@ -991,23 +1243,23 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None, If the desired pf is below some threshold for which RSP will be faster than CRP ('initial_packing_fraction'), only RSP is used. If a higher pf is - required, particles with a radius smaller than the desired final radius + required, spheres with a radius smaller than the desired final radius (and therefore with a smaller pf) are initialized within the domain using - RSP. This initial configuration of particles is then used as a starting + RSP. This initial configuration of spheres is then used as a starting point for CRP using Jodrey and Tory's algorithm [1]_. - In RSP, particle centers are placed one by one at random, and placement - attempts for a particle are made until the particle is not overlapping any + In RSP, sphere centers are placed one by one at random, and placement + attempts for a sphere are made until the sphere is not overlapping any others. This implementation of the algorithm uses a mesh over the domain - to speed up the nearest neighbor search by only searching for a particle's + to speed up the nearest neighbor search by only searching for a sphere's neighbors within that mesh cell. - In CRP, each particle is assigned two diameters, an inner and an outer, + In CRP, each sphere is assigned two diameters, an inner and an outer, which approach each other during the simulation. The inner diameter, defined as the minimum center-to-center distance, is the true diameter of - the particles and defines the pf. At each iteration the worst overlap - between particles based on outer diameter is eliminated by moving the - particles apart along the line joining their centers. Iterations continue + the spheres and defines the pf. At each iteration the worst overlap + between spheres based on outer diameter is eliminated by moving the + spheres apart along the line joining their centers. Iterations continue until the two diameters converge or until the desired pf is reached. References @@ -1016,81 +1268,68 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None, packing of equal spheres", Phys. Rev. A 32 (1985) 2347-2351. """ - - # Check for valid container geometry and dimensions - if domain_shape not in ['cube', 'cylinder', 'sphere']: - raise ValueError('Unable to set domain_shape to "{}". Only "cube", ' - '"cylinder", and "sphere" are ' - 'supported."'.format(domain_shape)) - if not domain_length and domain_shape in ['cube', 'cylinder']: - raise ValueError('"domain_length" must be specified for {} domain ' - 'geometry '.format(domain_shape)) - if not domain_radius and domain_shape in ['cylinder', 'sphere']: - raise ValueError('"domain_radius" must be specified for {} domain ' - 'geometry '.format(domain_shape)) - - if domain_shape == 'cube': - domain = _CubicDomain(length=domain_length, particle_radius=radius, - center=domain_center) - elif domain_shape == 'cylinder': - domain = _CylindricalDomain(length=domain_length, radius=domain_radius, - particle_radius=radius, center=domain_center) - elif domain_shape == 'sphere': - domain = _SphericalDomain(radius=domain_radius, particle_radius=radius, - center=domain_center) - - # Calculate the packing fraction if the number of particles is specified; - # otherwise, calculate the number of particles from the packing fraction. - if ((n_particles is None and packing_fraction is None) or - (n_particles is not None and packing_fraction is not None)): - raise ValueError('Exactly one of "n_particles" and "packing_fraction" ' - 'must be specified.') - elif packing_fraction is None: - n_particles = int(n_particles) - packing_fraction = 4/3*pi*radius**3*n_particles / domain.volume - elif n_particles is None: - packing_fraction = float(packing_fraction) - n_particles = int(packing_fraction*domain.volume // (4/3*pi*radius**3)) - - # Check for valid packing fractions for each algorithm - if packing_fraction >= 0.64: - raise ValueError('Packing fraction of {} is greater than the ' - 'packing fraction limit for close random ' - 'packing (0.64)'.format(packing_fraction)) - if initial_packing_fraction >= 0.38: - raise ValueError('Initial packing fraction of {} is greater than the ' - 'packing fraction limit for random sequential' - 'packing (0.38)'.format(initial_packing_fraction)) - if initial_packing_fraction > packing_fraction: - initial_packing_fraction = packing_fraction - if packing_fraction > 0.3: - initial_packing_fraction = 0.3 - + # Seed RNG random.seed(seed) - # Calculate the particle radius used in the initial random sequential + # Create container with the correct shape based on the supplied region + domain = None + for cls in _Container.__subclasses__(): + try: + domain = cls.from_region(region, radius) + except ValueError: + pass + + if not domain: + raise ValueError('Could not map region {} to a container: supported ' + 'container shapes are rectangular prism, cylinder, ' + 'sphere, and spherical shell.'.format(region)) + + # Determine the packing fraction/number of spheres + volume = 4/3*pi*radius**3 + if pf is None and num_spheres is None: + raise ValueError('`pf` or `num_spheres` must be specified.') + elif pf is None: + num_spheres = int(num_spheres) + pf = volume*num_spheres/domain.volume + else: + pf = float(pf) + num_spheres = int(pf*domain.volume//volume) + + # Make sure initial packing fraction is less than packing fraction + if initial_pf > pf: + initial_pf = pf + + # Check packing fraction for close random packing + if pf > MAX_PF_CRP: + raise ValueError('Packing fraction {0} is greater than the limit for ' + 'close random packing, {1}'.format(pf, MAX_PF_CRP)) + + # Check packing fraction for random sequential packing + if initial_pf > MAX_PF_RSP: + raise ValueError('Initial packing fraction {0} is greater than the ' + 'limit for random sequential packing, ' + '{1}'.format(initial_pf, MAX_PF_RSP)) + + # Calculate the sphere radius used in the initial random sequential # packing from the initial packing fraction - initial_radius = (3/4 * initial_packing_fraction * domain.volume / - (pi * n_particles))**(1/3) - domain.particle_radius = initial_radius + initial_radius = (3/4*initial_pf*domain.volume/(pi*num_spheres))**(1/3) + domain.sphere_radius = initial_radius # Recalculate the limits for the initial random sequential packing using - # the desired final particle radius to ensure particles are fully contained + # the desired final sphere radius to ensure spheres are fully contained # within the domain during the close random pack - domain.limits = [x + initial_radius - radius for x in domain.limits] + domain.limits = [[x - initial_radius + radius for x in domain.limits[0]], + [x + initial_radius - radius for x in domain.limits[1]]] - # Generate non-overlapping particles for an initial inner radius using + # Generate non-overlapping spheres for an initial inner radius using # random sequential packing algorithm - particles = _random_sequential_pack(domain, n_particles) + spheres = _random_sequential_pack(domain, num_spheres) - # Use the particle configuration produced in random sequential packing as a - # starting point for close random pack with the desired final particle + # Use the sphere configuration produced in random sequential packing as a + # starting point for close random pack with the desired final sphere # radius - if initial_packing_fraction != packing_fraction: - domain.particle_radius = radius - _close_random_pack(domain, particles, contraction_rate) + if initial_pf != pf: + domain.sphere_radius = radius + _close_random_pack(domain, spheres, contraction_rate) - trisos = [] - for p in particles: - trisos.append(TRISO(radius, fill, [x + c for x, c in zip(p, domain.center)])) - return trisos + return spheres diff --git a/openmc/region.py b/openmc/region.py index f82db8553..ceffadbe3 100644 --- a/openmc/region.py +++ b/openmc/region.py @@ -475,7 +475,7 @@ class Complement(Region): >>> xr = openmc.XPlane(x0=10.0) >>> yl = openmc.YPlane(y0=-10.0) >>> yr = openmc.YPlane(y0=10.0) - >>> inside_box = +xl & -xr & +yl & -yl + >>> inside_box = +xl & -xr & +yl & -yr >>> outside_box = ~inside_box >>> type(outside_box) diff --git a/openmc/surface.py b/openmc/surface.py index 826ac5cc8..c861880c5 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -62,7 +62,7 @@ class Surface(IDManagerMixin): self.boundary_type = boundary_type # A dictionary of the quadratic surface coefficients - # Key - coefficeint name + # Key - coefficient name # Value - coefficient value self._coefficients = {} @@ -1661,7 +1661,7 @@ class Quadric(Surface): a, b, c, d, e, f, g, h, j, k : float, optional coefficients for the surface. All default to 0. name : str, optional - Name of the sphere. If not specified, the name will be the empty string. + Name of the surface. If not specified, the name will be the empty string. Attributes ---------- diff --git a/tests/regression_tests/triso/inputs_true.dat b/tests/regression_tests/triso/inputs_true.dat index a13187628..7e5007eca 100644 --- a/tests/regression_tests/triso/inputs_true.dat +++ b/tests/regression_tests/triso/inputs_true.dat @@ -5,7 +5,7 @@ - + @@ -220,12 +220,12 @@ - - - - - - + + + + + + diff --git a/tests/regression_tests/triso/test.py b/tests/regression_tests/triso/test.py index 174bba94f..9ede6dbb0 100644 --- a/tests/regression_tests/triso/test.py +++ b/tests/regression_tests/triso/test.py @@ -53,19 +53,21 @@ class TRISOTestHarness(PyAPITestHarness): c5 = openmc.Cell(fill=opyc, region=+spheres[3]) inner_univ = openmc.Universe(cells=[c1, c2, c3, c4, c5]) - outer_radius = 422.5*1e-4 - trisos = openmc.model.pack_trisos( - radius=outer_radius, fill=inner_univ, domain_shape='cube', - domain_length=1., domain_center=(0., 0., 0.), n_particles=100) - - # Define box to contain lattice + # Define box to contain lattice and to pack TRISO particles in min_x = openmc.XPlane(x0=-0.5, boundary_type='reflective') max_x = openmc.XPlane(x0=0.5, boundary_type='reflective') min_y = openmc.YPlane(y0=-0.5, boundary_type='reflective') max_y = openmc.YPlane(y0=0.5, boundary_type='reflective') min_z = openmc.ZPlane(z0=-0.5, boundary_type='reflective') max_z = openmc.ZPlane(z0=0.5, boundary_type='reflective') - box = openmc.Cell(region=+min_x & -max_x & +min_y & -max_y & +min_z & -max_z) + box_region = +min_x & -max_x & +min_y & -max_y & +min_z & -max_z + box = openmc.Cell(region=box_region) + + outer_radius = 422.5*1e-4 + centers = openmc.model.pack_spheres(radius=outer_radius, + region=box_region, num_spheres=100) + trisos = [openmc.model.TRISO(outer_radius, inner_univ, c) + for c in centers] # Create lattice ll, ur = box.region.bounding_box diff --git a/tests/unit_tests/test_model_triso.py b/tests/unit_tests/test_model_triso.py index 19c4f9081..1d2e0905c 100644 --- a/tests/unit_tests/test_model_triso.py +++ b/tests/unit_tests/test_model_triso.py @@ -10,21 +10,88 @@ import pytest import scipy.spatial +_RADIUS = 0.1 _PACKING_FRACTION = 0.35 -_RADIUS = 4.25e-2 -domain_params = [ - {'shape': 'cube', 'length': 0.75, 'radius': 0., 'volume': 0.75**3}, - {'shape': 'cylinder', 'length': 0.5, 'radius': 0.5, 'volume': 0.5*pi*0.5**2}, - {'shape': 'sphere', 'length': 0., 'radius': 0.5, 'volume': 4/3*pi*0.5**3} +_PARAMS = [ + {'shape': 'rectangular_prism', 'volume': 1**3}, + {'shape': 'x_cylinder', 'volume': 1*pi*1**2}, + {'shape': 'y_cylinder', 'volume': 1*pi*1**2}, + {'shape': 'z_cylinder', 'volume': 1*pi*1**2}, + {'shape': 'sphere', 'volume': 4/3*pi*1**3}, + {'shape': 'spherical_shell', 'volume': 4/3*pi*(1**3 - 0.5**3)} ] -@pytest.fixture(scope='module', params=domain_params, - ids=['cube', 'cylinder', 'sphere']) -def domain(request): +@pytest.fixture(scope='module', params=_PARAMS) +def container(request): return request.param +@pytest.fixture(scope='module') +def centers(request, container): + return request.getfixturevalue('centers_' + container['shape']) + + +@pytest.fixture(scope='module') +def centers_rectangular_prism(): + min_x = openmc.XPlane(x0=0) + max_x = openmc.XPlane(x0=1) + min_y = openmc.YPlane(y0=0) + max_y = openmc.YPlane(y0=1) + min_z = openmc.ZPlane(z0=0) + max_z = openmc.ZPlane(z0=1) + region = +min_x & -max_x & +min_y & -max_y & +min_z & -max_z + return openmc.model.pack_spheres(radius=_RADIUS, region=region, + pf=_PACKING_FRACTION, initial_pf=0.2) + + +@pytest.fixture(scope='module') +def centers_x_cylinder(): + cylinder = openmc.XCylinder(R=1, y0=1, z0=2) + min_x = openmc.XPlane(x0=0) + max_x = openmc.XPlane(x0=1) + region = +min_x & -max_x & -cylinder + return openmc.model.pack_spheres(radius=_RADIUS, region=region, + pf=_PACKING_FRACTION, initial_pf=0.2) + + +@pytest.fixture(scope='module') +def centers_y_cylinder(): + cylinder = openmc.YCylinder(R=1, x0=1, z0=2) + min_y = openmc.YPlane(y0=0) + max_y = openmc.YPlane(y0=1) + region = +min_y & -max_y & -cylinder + return openmc.model.pack_spheres(radius=_RADIUS, region=region, + pf=_PACKING_FRACTION, initial_pf=0.2) + + +@pytest.fixture(scope='module') +def centers_z_cylinder(): + cylinder = openmc.ZCylinder(R=1, x0=1, y0=2) + min_z = openmc.ZPlane(z0=0) + max_z = openmc.ZPlane(z0=1) + region = +min_z & -max_z & -cylinder + return openmc.model.pack_spheres(radius=_RADIUS, region=region, + pf=_PACKING_FRACTION, initial_pf=0.2) + + +@pytest.fixture(scope='module') +def centers_sphere(): + sphere = openmc.Sphere(R=1, x0=1, y0=2, z0=3) + region = -sphere + return openmc.model.pack_spheres(radius=_RADIUS, region=region, + pf=_PACKING_FRACTION, initial_pf=0.2) + + +@pytest.fixture(scope='module') +def centers_spherical_shell(): + sphere = openmc.Sphere(R=1, x0=1, y0=2, z0=3) + inner_sphere = openmc.Sphere(R=0.5, x0=1, y0=2, z0=3) + region = -sphere & +inner_sphere + return openmc.model.pack_spheres(radius=_RADIUS, region=region, + pf=_PACKING_FRACTION, initial_pf=0.2) + + @pytest.fixture(scope='module') def triso_universe(): sphere = openmc.Sphere(R=_RADIUS) @@ -33,76 +100,96 @@ def triso_universe(): return univ -@pytest.fixture(scope='module') -def trisos(domain, triso_universe): - trisos = openmc.model.pack_trisos( - radius=_RADIUS, - fill=triso_universe, - domain_shape=domain['shape'], - domain_length=domain['length'], - domain_radius=domain['radius'], - domain_center=(0., 0., 0.), - initial_packing_fraction=0.2, - packing_fraction=_PACKING_FRACTION - ) - return trisos - - -def test_overlap(trisos): - """Check that no TRISO particles overlap.""" - centers = [t.center for t in trisos] - +def test_overlap(centers): + """Check that none of the spheres in the packed configuration overlap.""" # Create KD tree for quick nearest neighbor search tree = scipy.spatial.cKDTree(centers) - # Find distance to nearest neighbor for all particles + # Find distance to nearest neighbor for all spheres d = tree.query(centers, k=2)[0] - # Get the smallest distance between any two particles + # Get the smallest distance between any two spheres d_min = min(d[:, 1]) assert d_min > 2*_RADIUS or d_min == pytest.approx(2*_RADIUS) -def test_contained(trisos, domain): - """Make sure all particles are entirely contained within the domain.""" - if domain['shape'] == 'cube': - x = max(np.hstack([abs(t.center) for t in trisos])) + _RADIUS - assert x < 0.5*domain['length'] or x == pytest.approx(0.5*domain['length']) - - elif domain['shape'] == 'cylinder': - r = max([norm(t.center[0:2]) for t in trisos]) + _RADIUS - z = max([abs(t.center[2]) for t in trisos]) + _RADIUS - assert r < domain['radius'] or r == pytest.approx(domain['radius']) - assert z < 0.5*domain['length'] or z == pytest.approx(0.5*domain['length']) - - elif domain['shape'] == 'sphere': - r = max([norm(t.center) for t in trisos]) + _RADIUS - assert r < domain['radius'] or r == pytest.approx(domain['radius']) +def test_contained_rectangular_prism(centers_rectangular_prism): + """Make sure all spheres are entirely contained within the domain.""" + d_max = np.amax(centers_rectangular_prism) + _RADIUS + d_min = np.amin(centers_rectangular_prism) - _RADIUS + assert d_max < 1 or d_max == pytest.approx(1) + assert d_min > 0 or d_min == pytest.approx(0) -def test_packing_fraction(trisos, domain): +def test_contained_x_cylinder(centers_x_cylinder): + """Make sure all spheres are entirely contained within the domain.""" + d = np.linalg.norm(centers_x_cylinder[:,[1,2]] - [1, 2], axis=1) + r_max = max(d) + _RADIUS + x_max = max(centers_x_cylinder[:,0]) + _RADIUS + x_min = min(centers_x_cylinder[:,0]) - _RADIUS + assert r_max < 1 or r_max == pytest.approx(1) + assert x_max < 1 or x_max == pytest.approx(1) + assert x_min > 0 or x_min == pytest.approx(0) + + +def test_contained_y_cylinder(centers_y_cylinder): + """Make sure all spheres are entirely contained within the domain.""" + d = np.linalg.norm(centers_y_cylinder[:,[0,2]] - [1, 2], axis=1) + r_max = max(d) + _RADIUS + y_max = max(centers_y_cylinder[:,1]) + _RADIUS + y_min = min(centers_y_cylinder[:,1]) - _RADIUS + assert r_max < 1 or r_max == pytest.approx(1) + assert y_max < 1 or y_max == pytest.approx(1) + assert y_min > 0 or y_min == pytest.approx(0) + + +def test_contained_z_cylinder(centers_z_cylinder): + """Make sure all spheres are entirely contained within the domain.""" + d = np.linalg.norm(centers_z_cylinder[:,[0,1]] - [1, 2], axis=1) + r_max = max(d) + _RADIUS + z_max = max(centers_z_cylinder[:,2]) + _RADIUS + z_min = min(centers_z_cylinder[:,2]) - _RADIUS + assert r_max < 1 or r_max == pytest.approx(1) + assert z_max < 1 or z_max == pytest.approx(1) + assert z_min > 0 or z_min == pytest.approx(0) + + +def test_contained_sphere(centers_sphere): + """Make sure all spheres are entirely contained within the domain.""" + d = np.linalg.norm(centers_sphere - [1, 2, 3], axis=1) + r_max = max(d) + _RADIUS + assert r_max < 1 or r_max == pytest.approx(1) + + +def test_contained_spherical_shell(centers_spherical_shell): + """Make sure all spheres are entirely contained within the domain.""" + d = np.linalg.norm(centers_spherical_shell - [1, 2, 3], axis=1) + r_max = max(d) + _RADIUS + r_min = min(d) - _RADIUS + assert r_max < 1 or r_max == pytest.approx(1) + assert r_min > 0.5 or r_min == pytest.approx(0.5) + + +def test_packing_fraction(container, centers): """Check that the actual PF is close to the requested PF.""" - pf = len(trisos)*4/3*pi*_RADIUS**3/domain['volume'] + pf = len(centers) * 4/3 * pi *_RADIUS**3 / container['volume'] assert pf == pytest.approx(_PACKING_FRACTION, rel=1e-2) -def test_n_particles(triso_universe): - """Check that the function returns the correct number of particles""" - trisos = openmc.model.pack_trisos( - radius=_RADIUS, fill=triso_universe, domain_shape='cube', - domain_length=1.0, n_particles=800 +def test_num_spheres(): + """Check that the function returns the correct number of spheres""" + centers = openmc.model.pack_spheres( + radius=_RADIUS, region=-openmc.Sphere(R=1), num_spheres=50 ) - assert len(trisos) == 800 + assert len(centers) == 50 -def test_triso_lattice(triso_universe): - trisos = openmc.model.pack_trisos( - radius=_RADIUS, fill=triso_universe, domain_shape='cube', - domain_length=1.0, domain_center=(0., 0., 0.), packing_fraction=0.2 - ) +def test_triso_lattice(triso_universe, centers_rectangular_prism): + trisos = [openmc.model.TRISO(_RADIUS, triso_universe, c) + for c in centers_rectangular_prism] - lower_left = np.array((-.5, -.5, -.5)) - upper_right = np.array((.5, .5, .5)) + lower_left = np.array((0, 0, 0)) + upper_right = np.array((1, 1, 1)) shape = (3, 3, 3) pitch = (upper_right - lower_left)/shape background = openmc.Material() @@ -112,50 +199,29 @@ def test_triso_lattice(triso_universe): ) -def test_domain_input(triso_universe): - # Invalid domain shape +def test_container_input(triso_universe): + # Invalid container shape with pytest.raises(ValueError): - trisos = openmc.model.pack_trisos( - radius=1, fill=triso_universe, n_particles=100, - domain_shape='circle' + centers = openmc.model.pack_spheres( + radius=_RADIUS, region=+openmc.Sphere(R=1), num_spheres=100 ) - # Don't specify domain length on a cube - with pytest.raises(ValueError): - trisos = openmc.model.pack_trisos( - radius=1, fill=triso_universe, n_particles=100, - domain_shape='cube' - ) - # Don't specify domain radius on a sphere - with pytest.raises(ValueError): - trisos = openmc.model.pack_trisos( - radius=1, fill=triso_universe, n_particles=100, - domain_shape='sphere' - ) - -def test_packing_fraction_input(triso_universe): - # Provide neither packing fraction nor number of particles + +def test_packing_fraction_input(): + # Provide neither packing fraction nor number of spheres with pytest.raises(ValueError): - trisos = openmc.model.pack_trisos( - radius=1, fill=triso_universe, domain_shape='cube', - domain_length=10 - ) - # Provide both packing fraction and number of particles - with pytest.raises(ValueError): - trisos = openmc.model.pack_trisos( - radius=1, fill=triso_universe, domain_shape='cube', - domain_length=10, n_particles=100, packing_fraction=0.2 + centers = openmc.model.pack_spheres( + radius=_RADIUS, region=-openmc.Sphere(R=1) ) + # Specify a packing fraction that is too high for CRP with pytest.raises(ValueError): - trisos = openmc.model.pack_trisos( - radius=1, fill=triso_universe, domain_shape='cube', - domain_length=10, packing_fraction=1 + centers = openmc.model.pack_spheres( + radius=_RADIUS, region=-openmc.Sphere(R=1), pf=1 ) + # Specify a packing fraction that is too high for RSP with pytest.raises(ValueError): - trisos = openmc.model.pack_trisos( - radius=1, fill=triso_universe, domain_shape='cube', - domain_length=10, packing_fraction=0.5, - initial_packing_fraction=0.4 + centers = openmc.model.pack_spheres( + radius=_RADIUS, region=-openmc.Sphere(R=1), pf=0.5, initial_pf=0.4 )