04b_Geometria.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.6.7"
    },
    "colab": {
      "name": "Geometria.ipynb",
      "version": "0.3.2",
      "provenance": []
    }
  },
  "cells": [
    {
      "cell_type": "code",
      "metadata": {
        "id": "NJvMsxAgkWuc",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#\n",
        "# Executing this cell you will install OpenMC and the nuclear\n",
        "# data libraries in this instance of the Google Colaboratory virtual machine.\n",
        "# The process takes about 2 minutes.\n",
        "#\n",
        "\n",
        "def install_openmc():\n",
        "  #\n",
        "  # Clone source code from Github, make and install\n",
        "  #\n",
        "  \n",
        "  import os\n",
        "  \n",
        "  if not os.path.isdir('/content'):\n",
        "    print(\"Esta función instala OpenMC en una instancia de Google Colaboratory.\")\n",
        "    print(\"Para instalar localmente siga las instrucciones de la documentacion:\")\n",
        "    print(\"http://docs.openmc.org/en/stable/quickinstall.html\")\n",
        "    return\n",
        "  \n",
        "  %cd -q /content\n",
        "  print(\"Obtaining OpenMC source code from Github...\")\n",
        "  !git  --no-pager clone https://github.com/mit-crpg/openmc.git &> /dev/null\n",
        "  %cd -q openmc\n",
        "  !git --no-pager checkout master &> /dev/null\n",
        "  !mkdir build\n",
        "  %cd -q build\n",
        "  print(\"Running cmake...\")\n",
        "  !cmake ../ -DPYTHON_EXECUTABLE=/usr/bin/python3 -DCMAKE_INSTALL_PREFIX=/usr/local &> /dev/null\n",
        "  print(\"Running make...\")\n",
        "  !make -j &> /dev/null\n",
        "  print(\"Running make install...\")\n",
        "  !make install &> /dev/null\n",
        "  import sys\n",
        "  sys.path.append('/usr/local/lib/python3.6/dist-packages/openmc-0.10.0-py3.6-linux-x86_64.egg')\n",
        "  sys.path.append('/usr/local/lib/python3.6/site-packages/')\n",
        "  %cd -q /content\n",
        "\n",
        "def install_data_from_onedrive():\n",
        "  #\n",
        "  # Download preprocessed HDF5 files from OneDrive (faster).\n",
        "  #\n",
        "  import os\n",
        "  \n",
        "  if not os.path.isdir('/content'):\n",
        "    print(\"Esta función instala los datos nucleares de OpenMC en una instancia de Google Colaboratory.\")\n",
        "    print(\"Para instalar localmente siga las instrucciones de la documentacion:\")\n",
        "    print(\"http://docs.openmc.org/en/stable/quickinstall.html\")\n",
        "    return\n",
        "  %cd -q /content\n",
        "  print(\"Obtaining HDF5 files from OneDrive...\")\n",
        "  !wget -O nndc_hdf5.tar.gz \"https://onedrive.live.com/download?cid=22422A8EEA2A85B9&resid=22422A8EEA2A85B9%21105&authkey=AHL6xwYFXDwEzkk\" &> /dev/null\n",
        "  print(\"Uncompressing...\")\n",
        "  !tar xzvf nndc_hdf5.tar.gz &> /dev/null\n",
        "    \n",
        "from time import time\n",
        "t1 = time()\n",
        "install_openmc()\n",
        "install_data_from_onedrive()\n",
        "t2 = time()\n",
        "print(\"Installed OpenMC in {:.2f} minutes\".format((t2-t1)/60.0))\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "uxDY3tJdkUkq",
        "colab_type": "text"
      },
      "source": [
        "# Ejemplos de definición de superficies y regiones en OpenMC"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "fGr8dpdskUkr",
        "colab_type": "text"
      },
      "source": [
        "Comenzamos importando el módulo `openmc`, contiene el API (applicationg programming interface) para generar el input de OpenMC."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "vsO9tCe1kUks",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import openmc"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "X1SmhgtLkUk0",
        "colab_type": "text"
      },
      "source": [
        "La definición de las distintas superficies se realiza con funciones del paquete `openmc`. Por ejemplo, para generar un plano perpendicular al eje $x$ hacemos:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "F30otsBHkUk3",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "plano = openmc.XPlane()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "XRtOK0uUkUk8",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "print(type(plano))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QWVXitDGkUlA",
        "colab_type": "text"
      },
      "source": [
        "El objeto de clase `XPlane` posee una serie de atributos y métodos que permiten modificarlo. Por ejemplo, para modificar el punto en que corta el eje $x$ hacemos:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "oNbwN6b-kUlB",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "plano.x0 = 10.0\n",
        "plano.name = \"Plano 1\""
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "fZjv-jh1kUlF",
        "colab_type": "text"
      },
      "source": [
        "Esto es equivalente a crear el objeto dándole los valores como parámetros a la función `XPlane()`:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "pVr38fr9kUlG",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "plano2 = openmc.XPlane(x0=10.0, name=\"Plano 2\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "hr-AFIEAkUlJ",
        "colab_type": "text"
      },
      "source": [
        "Una superficie perpendicular al eje $x$ está definida por la ecuación $x-x_0=0$, por lo que es necesario definir un único parámetro ($x_0$) para definir el plano. De forma similar para planos perpendiculares a los ejes $y$ y $z$:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "twBNVu1akUlK",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "plano3 = openmc.YPlane(y0=5)\n",
        "plano4 = openmc.ZPlane(z0=3)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "sQnZ19tkkUlN",
        "colab_type": "text"
      },
      "source": [
        "Para definir una esfera es necesario indicar la posición $x_0, y_0, z_0$ del centro, y su radio ya que está definida por $(x-x_0)^2+(y-y_0)^2+(z-z_0)^2-R^2=0$:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "5uz-TFfckUlO",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "esfera = openmc.Sphere(x0=10, R=100)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kVk1PQ23kUlR",
        "colab_type": "text"
      },
      "source": [
        "(los valores $y_0$ y $z_0$ que no fueron especificados toman su valor por defecto, 0)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "KNHyjaAEkUlS",
        "colab_type": "text"
      },
      "source": [
        "Se pueden definir cilindros paralelos a los ejes $x$, $y$, $z$ respectivamente utilizando las funciones `XCylinder()`, `YCylinder()` y `ZCylinder()`:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "1Vp9qRfxkUlU",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Un cilindro paralelo al eje x de 10 cm de radio, que pasa por y=2 cm, z=1 cm\n",
        "cil1 = openmc.XCylinder(y0=2, z0=1, R=10) \n",
        "\n",
        "# Un cilindro centrado en el eje y de radio 5 cm\n",
        "cil2 = openmc.YCylinder(R=5) \n",
        "\n",
        "# Un cilindro centrado en el eje z de radio 1 cm\n",
        "cil3 = openmc.ZCylinder() "
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6Wr8_7AMkUlW",
        "colab_type": "text"
      },
      "source": [
        "Además existen las funciones:\n",
        "* `Plane()`: define un plano con orientación arbitraria.\n",
        "* `XCone()`: define un cono con eje paralelo al eje $x$.\n",
        "* `YCone()`: define un cono con eje paralelo al eje $y$.\n",
        "* `ZCone()`: define un cono con eje paralelo al eje $z$.\n",
        "* `Quadric()`: Define una superficie cuadrática arbitraria que responde a la ecuación $Ax^2+By^2+Cz^2+Dxy+Eyz+Fxz+Gx+Hy+Jz+K=0$\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WTI90FjvkUlX",
        "colab_type": "text"
      },
      "source": [
        "Cada superficie que generamos tiene un identificador. Si no lo indicamos explícitamente al invocar la función que genera las superficies estos identificadores se asignan en forma secuencial:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "2ovYu0fQkUlX",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "print(plano.id)\n",
        "print(cil3.id)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Z3-4t3EYkUla",
        "colab_type": "text"
      },
      "source": [
        "Estos identificadores parecen en principio redundantes, ya que cada superficie tiene su objeto y potencialmente un nombre que la identifica. Pero, debido a que en otros códigos Monte Carlo suelen identificar las superficies con números enteros, es conveniente asignar los valores en la incialización de las superficies para facilitar la comparación:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "llCO2rI4kUlb",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "otro_plano = openmc.XPlane(x0=10, surface_id=999, name=\"Otro Plano\")\n",
        "print(otro_plano.id)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "AFdS9DtIkUlg",
        "colab_type": "text"
      },
      "source": [
        "## Regiones"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "GI3LyGvEkUlh",
        "colab_type": "text"
      },
      "source": [
        "Una vez que se creó una superficie, los operadores `+` y `-` permiten definir semiespacios:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "YH919jc8kUlh",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "adelante = +plano2\n",
        "atras = -plano2\n",
        "print(type(adelante))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dCDgGAlDkUll",
        "colab_type": "text"
      },
      "source": [
        "Estos semiespacios pueden combinarse para formar regiones utilizando los operadores `&` (intersección), `|` (unión) y `~` (complemento). Por ejemplo, si quisieramos definir la región dentro de una esfera y arriba de su plano medio podemos hacer:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "qPCMtyjMkUll",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Defino dos superficies:\n",
        "esfera = openmc.Sphere(R=10)\n",
        "plano_medio = openmc.ZPlane(z0=0)\n",
        "\n",
        "# Defino algunas semiespacios en base a estas superficies\n",
        "adentro_de_la_esfera = -esfera\n",
        "arriba_del_plano_medio = +plano_medio\n",
        "\n",
        "# Combino estas regiones para generar nuevas regiones\n",
        "hemisferio_norte = adentro_de_la_esfera & arriba_del_plano_medio\n",
        "print(type(hemisferio_norte))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZG-ppYE5kUlo",
        "colab_type": "text"
      },
      "source": [
        "La región dentro de la esfera pero abajo del plano medio es:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "60dPEwM5kUlp",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "hemisferio_sur = adentro_de_la_esfera & ~hemisferio_norte"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "e5mtuIjGkUlu",
        "colab_type": "text"
      },
      "source": [
        "A veces para hacer estas construcciones es útil recordar las reglas de De Morgan:\n",
        "* `~(A | B) = ~A & ~B`\n",
        "* `~(A & B) = ~A | ~B`"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Aif6hBq0kUlv",
        "colab_type": "text"
      },
      "source": [
        "Utilizando esto otra forma de definir esta región es:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "pNbBkQFekUly",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "hemisferio_sur = adentro_de_la_esfera & ~(adentro_de_la_esfera & arriba_del_plano_medio)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Jjbq2zXgkUl1",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "hemisferio_sur = adentro_de_la_esfera & (~adentro_de_la_esfera | ~arriba_del_plano_medio)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "EEcQgWoOkUl3",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "hemisferio_sur = -esfera & (~ -esfera | ~ +plano_medio)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "nFXQCrKqkUl6",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "hemisferio_sur = -esfera & ( +esfera | -plano_medio)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "08vwmZE7kUl9",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "hemisferio_sur = (-esfera &  +esfera) | (-esfera & -plano_medio)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "69nXRcrSkUl-",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "hemisferio_sur = (-esfera & -plano_medio)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "XzLZDdBPkUmB",
        "colab_type": "text"
      },
      "source": [
        "Es conveniente definir las regiones de la forma más sencilla posible para evitar errores de geometría y acelerar la ejecución del código."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BM3vOwC4kUmC",
        "colab_type": "text"
      },
      "source": [
        "## Celdas"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "wbaKnUmFkUmD",
        "colab_type": "text"
      },
      "source": [
        "Una vez que uno definió una region, es necesario definir una celda y asignarle materiales a esa region utilizando la función `.Cell()`:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "UqSSpIfvkUmD",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "umet = openmc.Material()\n",
        "umet.add_nuclide(\"U235\", 1.0, \"ao\")\n",
        "umet.set_density(\"g/cm3\", 19.0)\n",
        "cel_nucleo_del_reactor = openmc.Cell(fill=umet, region=hemisferio_sur)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Xh65ZoxTkUmG",
        "colab_type": "text"
      },
      "source": [
        "En este ejemplo el parámetro `fill=` utiliza un material, pero como veremos una celda también puede llenarse con otras celdas, reunidas en universos."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "K54aFxyKkUmG",
        "colab_type": "text"
      },
      "source": [
        "## Condiciones de contorno"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Gijw4TBIkUmH",
        "colab_type": "text"
      },
      "source": [
        "A diferencia de otros códigos Monte Carlo (como MCNP), en los que las condiciones de contorno se establecen en forma volumétrica, las condicionces de contorno en OpenMC se establecen sobre superficies con el parámetro `boundary_type=`. El valor por default para este parámetro es `transmission`, que permite que las partículas atraviesen la superficie sin ningún impedimento. Si esta parámetro se coloca en `vacuum`, OpenMC terminará la historia de las partículas que atraviesan la superficie, generando una condición de contorno de corriente entrante nula (esto es, vacío). Si se coloca el parámetro `reflective` la superficie reflejará las partículas que lleguen a ella, generando una condición de contorno reflejada."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "GSKkcsTpkUmI",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "esfera = openmc.Sphere(R=10, boundary_type='vacuum')"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "c4dLZNxnkUmL",
        "colab_type": "text"
      },
      "source": [
        "Como estas condiciones de contorno pueden aplicarse a cualquier superficie, incluyendo a superficies interiores, hay que tener cuidado de aplicar condiciones de contorno que tengan sentido físico."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "cr5Q2KsdkUmM",
        "colab_type": "text"
      },
      "source": [
        "## Universos"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Dskp3AQekUmN",
        "colab_type": "text"
      },
      "source": [
        "En OpenMC se denomina *universo* a una colección de celdas con un nombre asignado. Como mínimo debe haber un universo en el sistema, que se denomina `root universe`. Dentro de las celdas de este universo pueden haber otros universos, lo que permite generar estructuras repetidas (por ejemplo, conjuntos de placas combustibles en un reactor tipo MTR). En el siguiente ejemplo definiremos un cubo dentro de una esfera y en base a eso tres celdas: \n",
        "\n",
        "1) el interior del cubo, \n",
        "\n",
        "2) el exterior del cubo dentro de la esfera, y\n",
        "\n",
        "3) el exterior de la esfera.\n",
        "\n",
        "Al interior del cubo lo llenaremos de polietileno y a las esfera de agua. El exterior de la esfera será vacío."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "dcCwuUEekUmN",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "agua = openmc.Material()\n",
        "agua.add_nuclide(\"H1\", 2.0, 'ao')\n",
        "agua.add_nuclide(\"O16\", 1.0, 'ao')\n",
        "agua.set_density('g/cm3', 1.0)\n",
        "\n",
        "ple = openmc.Material()\n",
        "ple.add_nuclide(\"H1\", 2.0, 'ao')\n",
        "ple.add_nuclide(\"C0\", 1.0, 'ao')\n",
        "ple.set_density('g/cm3', 0.9)\n",
        "\n",
        "sup_plano1 = openmc.XPlane(x0=-7.5)\n",
        "sup_plano2 = openmc.XPlane(x0=+7.5)\n",
        "sup_plano3 = openmc.YPlane(y0=-7.5)\n",
        "sup_plano4 = openmc.YPlane(y0=+7.5)\n",
        "sup_plano5 = openmc.ZPlane(z0=-7.5)\n",
        "sup_plano6 = openmc.ZPlane(z0=+7.5)\n",
        "\n",
        "sup_esfera = openmc.Sphere(R=15)\n",
        "adentro_cubo = (+sup_plano1 & -sup_plano2) & \\\n",
        "               (+sup_plano3 & -sup_plano4) & \\\n",
        "               (+sup_plano5 & -sup_plano6) \n",
        "afuera_cubo = ~adentro_cubo\n",
        "adentro_esfera = -sup_esfera\n",
        "afuera_esfera = +sup_esfera\n",
        "\n",
        "cel_cubo = openmc.Cell(region=adentro_cubo, fill=ple)\n",
        "cel_esfera = openmc.Cell(region=adentro_esfera & afuera_cubo, fill=agua)\n",
        "cel_afuera = openmc.Cell(region=afuera_esfera & afuera_cubo, fill=None)\n",
        "\n",
        "universo = openmc.Universe(cells= [cel_cubo, cel_esfera, cel_afuera])"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "lRyV_qEnkUmR",
        "colab_type": "text"
      },
      "source": [
        "Los objetos de clase `Universe` pueden graficarse utilizando el método `.plot()`:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "MdRQ3Et8kUmS",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "universo.plot(width=(40,40))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Q8QYfsaMkUmW",
        "colab_type": "text"
      },
      "source": [
        "Aquí el parámetro `width=` toma una tupla que representa las dimensiones de la región a graficar en cm."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "alxugIEOk2Zb",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        ""
      ],
      "execution_count": 0,
      "outputs": []
    }
  ]
}