{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Firefly 003 \u2014 Google Colab edition\n",
    "\n",
    "Run the workshop in Google's cloud, with no local OpenMC installation.\n",
    "Sign in to Google and select **Runtime \u2192 Change runtime type \u2192 Python 3 / CPU**.\n",
    "No GPU or paid runtime is required; free runtime availability is controlled by Google.\n",
    "\n",
    "1. Run **A** once. Colab will disconnect briefly while CondaColab restarts Python.\n",
    "   This is expected. Wait for it to reconnect, then continue with **B**, not Run all.\n",
    "2. Run **B** to install OpenMC (allow several minutes).\n",
    "3. Run **C** and upload `openmc-reactor-circle.zip`, downloaded from the\n",
    "   Firefly 003 setup page. Leave the ZIP compressed.\n",
    "4. Run the original workshop cells below in order using their play buttons.\n",
    "\n",
    "On a phone, open Colab in Safari/Chrome; the wider desktop layout may be easier\n",
    "to use. Save a copy to Drive to keep notebook edits. Runtime files are temporary:\n",
    "after a runtime reset, repeat setup and upload the ZIP again. Download results\n",
    "with the final cell before disconnecting.\n",
    "\n",
    "This is an educational simulation, not a validated reactor engineering design.\n"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "# A \u2014 Install CondaColab. An automatic runtime restart is expected.\n",
    "import sys\n",
    "import subprocess\n",
    "subprocess.run([sys.executable, \"-m\", \"pip\", \"install\", \"-q\", \"condacolab==0.1.10\"], check=True)\n",
    "import condacolab\n",
    "condacolab.install()\n"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "# B \u2014 Run after the automatic restart.\n",
    "import condacolab\n",
    "condacolab.check()\n",
    "import os\n",
    "import sys\n",
    "import subprocess\n",
    "subprocess.run([\n",
    "    \"conda\", \"install\", \"-y\", \"--prefix\", sys.prefix,\n",
    "    \"--override-channels\", \"-c\", \"conda-forge\",\n",
    "    \"openmc=0.16\", \"matplotlib\"\n",
    "], check=True)\n",
    "os.environ[\"OMP_NUM_THREADS\"] = \"2\"\n",
    "import openmc\n",
    "print(\"OpenMC\", openmc.__version__, \"ready\")\n"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "# C \u2014 Upload the workshop ZIP; only its nuclear data will be extracted.\n",
    "from google.colab import files\n",
    "from pathlib import Path\n",
    "from hashlib import sha256\n",
    "from zipfile import ZipFile\n",
    "from io import BytesIO\n",
    "import os\n",
    "import xml.etree.ElementTree as ET\n",
    "\n",
    "work = Path(\"/content/firefly-003\")\n",
    "work.mkdir(exist_ok=True)\n",
    "os.chdir(work)\n",
    "uploaded = files.upload()\n",
    "if len(uploaded) != 1:\n",
    "    raise ValueError(\"Upload exactly one file: openmc-reactor-circle.zip\")\n",
    "payload = next(iter(uploaded.values()))\n",
    "if sha256(payload).hexdigest() != \"0fb4ef3efa0b9676d493326b78aacd8b47e6ff4e7d0b20633fe8d8c707931000\":\n",
    "    raise ValueError(\"Wrong ZIP. Download openmc-reactor-circle.zip from the Firefly 003 setup page.\")\n",
    "with ZipFile(BytesIO(payload)) as archive:\n",
    "    for member in archive.namelist():\n",
    "        prefix = \"openmc-reactor-circle/xs_data/\"\n",
    "        if member.startswith(prefix) and not member.endswith(\"/\"):\n",
    "            relative = Path(member.removeprefix(prefix))\n",
    "            if relative.is_absolute() or \"..\" in relative.parts:\n",
    "                raise ValueError(\"Unsafe archive path\")\n",
    "            target = work / \"xs_data\" / relative\n",
    "            target.parent.mkdir(parents=True, exist_ok=True)\n",
    "            target.write_bytes(archive.read(member))\n",
    "xs = work / \"xs_data/cross_sections.xml\"\n",
    "for library in ET.parse(xs).getroot().findall(\"library\"):\n",
    "    if not (xs.parent / library.attrib[\"path\"]).is_file():\n",
    "        raise FileNotFoundError(library.attrib[\"path\"])\n",
    "os.environ[\"OPENMC_CROSS_SECTIONS\"] = str(xs)\n",
    "del uploaded, payload\n",
    "print(\"Nuclear data ready. Continue with the workshop below.\")\n"
   ],
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "e1057945",
   "metadata": {},
   "source": [
    "# \ud83d\udd25 Firefly Reading Circle \u2014 Build a Nuclear Reactor in OpenMC\n",
    "\n",
    "**Goal (1 hour):** build, run, and visualize a simplified nuclear reactor\n",
    "model, understand k-effective, and experiment with what makes a reactor\n",
    "more or less reactive.\n",
    "\n",
    "Run each cell top-to-bottom with **Shift+Enter**. Every explanation lives\n",
    "in its own markdown (text) cell, right above the code it describes \u2014 code\n",
    "cells are kept short and mostly comment-free so the two stay visually\n",
    "distinct as you scroll.\n",
    "\n",
    "**Companion doc:** `READING_CIRCLE_GUIDE.md` in this folder has the full\n",
    "facilitator run-of-show, timing, and discussion prompts."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5213f574",
   "metadata": {},
   "source": [
    "---\n",
    "## 0. Setup check\n",
    "\n",
    "Run this first. It loads the libraries and points OpenMC at the nuclear\n",
    "data that ships with this folder (`xs_data/`), so there's nothing else to\n",
    "configure. It should print the OpenMC version and \"ready\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52aaf3f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import glob\n",
    "from pathlib import Path\n",
    "\n",
    "XS = Path(\"xs_data/cross_sections.xml\").resolve()\n",
    "os.environ[\"OPENMC_CROSS_SECTIONS\"] = str(XS)\n",
    "\n",
    "import openmc\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.image as mpimg\n",
    "\n",
    "openmc.config[\"cross_sections\"] = XS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d8f1812",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"OpenMC version:\", openmc.__version__)\n",
    "print(\"Nuclear data  :\", XS, \"\u2705 ready\" if XS.exists() else \"\u274c MISSING \u2014 open this notebook from the project folder\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85d498ed",
   "metadata": {},
   "source": [
    "---\n",
    "## 1. The big idea\n",
    "\n",
    "A nuclear reactor is just three ingredients in careful balance:\n",
    "\n",
    "- **Fuel** \u2014 releases neutrons when it fissions. We use UO2 enriched to\n",
    "  4% U-235 (typical PWR fuel).\n",
    "- **Moderator** \u2014 slows fast neutrons down so they're more likely to cause\n",
    "  another fission. We use ordinary (lightly borated) water.\n",
    "- **Control rods** \u2014 deliberately *absorb* neutrons to control the chain\n",
    "  reaction. We use boron carbide (B4C). Push them in to slow the reaction\n",
    "  down, pull them out to speed it up.\n",
    "\n",
    "The chain reaction is summarized by one number, **k-effective (keff)**:\n",
    "the average number of neutrons from one fission that go on to cause\n",
    "another fission.\n",
    "\n",
    "| keff | meaning |\n",
    "|---|---|\n",
    "| = 1.000 | **critical** \u2014 steady, self-sustaining chain reaction (what a running plant targets) |\n",
    "| > 1.000 | **supercritical** \u2014 neutron population (and power) growing |\n",
    "| < 1.000 | **subcritical** \u2014 neutron population (and power) shrinking |\n",
    "\n",
    "OpenMC is an open-source Monte Carlo particle transport code. It literally\n",
    "simulates many thousands of individual neutrons bouncing through the\n",
    "geometry we define, tracking every scatter/absorption/fission event, and\n",
    "from that statistics estimates keff and where the power is produced.\n",
    "\n",
    "No code in this section \u2014 just read, then continue to Section 2."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90e26fe2",
   "metadata": {},
   "source": [
    "---\n",
    "## 2. Materials\n",
    "\n",
    "Five materials make up this model. We'll define each one in its own cell.\n",
    "\n",
    "### 2a. Fuel \u2014 UO2, 4% enriched\n",
    "`enrichment=4.0` is the fissile U-235 weight percent \u2014 the main knob for\n",
    "\"experiment b\" later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d8af726",
   "metadata": {},
   "outputs": [],
   "source": [
    "fuel = openmc.Material(name=\"UO2 fuel (4.0% enriched)\")\n",
    "fuel.add_element('U', 1.0, enrichment=4.0)\n",
    "fuel.add_element('O', 2.0)\n",
    "fuel.set_density('g/cm3', 10.4)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f9e4a64",
   "metadata": {},
   "source": [
    "### 2b. Cladding \u2014 Zircaloy\n",
    "Simplified to pure zirconium for clarity in a teaching model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a687d6f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "clad = openmc.Material(name=\"Zircaloy cladding\")\n",
    "clad.add_element('Zr', 1.0)\n",
    "clad.set_density('g/cm3', 6.55)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ddee7b3",
   "metadata": {},
   "source": [
    "### 2c. Moderator/coolant \u2014 borated water\n",
    "~500 ppm natural boron and a hot PWR operating density (0.74 g/cm\u00b3) are\n",
    "both realistic plant values. Boron here is a deliberate neutron poison \u2014\n",
    "\"experiment c\" removes it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "efc2804f",
   "metadata": {},
   "outputs": [],
   "source": [
    "water = openmc.Material(name=\"Borated water moderator\")\n",
    "water.add_element('H', 2.0)\n",
    "water.add_element('O', 1.0)\n",
    "water.add_element('B', 500e-6)\n",
    "water.set_density('g/cm3', 0.74)\n",
    "water.add_s_alpha_beta('c_H_in_H2O')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "325b7b5d",
   "metadata": {},
   "source": [
    "### 2d. Control rod absorber \u2014 boron carbide (B4C)\n",
    "Deliberately a strong neutron absorber."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b3edbe8",
   "metadata": {},
   "outputs": [],
   "source": [
    "b4c = openmc.Material(name=\"B4C control rod absorber\")\n",
    "b4c.add_element('B', 4.0)\n",
    "b4c.add_element('C', 1.0)\n",
    "b4c.set_density('g/cm3', 2.52)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a6fd628",
   "metadata": {},
   "source": [
    "### 2e. Air \u2014 withdrawn control-rod gas gap\n",
    "Fills the thin gap inside the pin cladding."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc9623f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "air = openmc.Material(name=\"Air (void/instrument channel)\")\n",
    "air.add_element('N', 0.78)\n",
    "air.add_element('O', 0.21)\n",
    "air.add_element('Ar', 0.01)\n",
    "air.set_density('g/cm3', 0.001225)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "957937e7",
   "metadata": {},
   "source": [
    "### 2f. Collect and export\n",
    "OpenMC reads materials from an XML file it writes from this Python object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4487a61",
   "metadata": {},
   "outputs": [],
   "source": [
    "materials = openmc.Materials([fuel, clad, water, b4c, air])\n",
    "materials.export_to_xml()\n",
    "print(\"Materials defined:\", [m.name for m in materials])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85aaf5bf",
   "metadata": {},
   "source": [
    "---\n",
    "## 3. Geometry: pin cell \u2192 lattice \u2192 finite core\n",
    "\n",
    "We build a single **fuel pin cell** (concentric cylinders: fuel pellet \u2192\n",
    "gas gap \u2192 zircaloy clad \u2192 water), then repeat it in a **7\u00d77 lattice** to\n",
    "make a mini assembly. A cross of control-rod pins runs through the middle,\n",
    "with a couple of water-filled guide tubes. The whole thing sits inside a\n",
    "**water reflector** with a **vacuum boundary** \u2014 meaning neutrons that\n",
    "escape the edge are gone for good. This is a real *finite* reactor model,\n",
    "not an idealized infinite one.\n",
    "\n",
    "### 3a. Pin-cell dimensions and cylinders\n",
    "Standard PWR pin pitch and radii (cm)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "094e37a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "pin_pitch = 1.26\n",
    "fuel_or = 0.39\n",
    "clad_ir = 0.40\n",
    "clad_or = 0.46\n",
    "\n",
    "fuel_cyl = openmc.ZCylinder(r=fuel_or)\n",
    "clad_in_cyl = openmc.ZCylinder(r=clad_ir)\n",
    "clad_out_cyl = openmc.ZCylinder(r=clad_or)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e2f996e",
   "metadata": {},
   "source": [
    "### 3b. Three pin-cell \"universes\"\n",
    "A **fuel pin**, a **control-rod pin** (same clad/water, B4C instead of\n",
    "fuel), and a **water guide tube** (no fuel or absorber at all)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7f2da5b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_fuel_pin():\n",
    "    \"\"\"Standard UO2 fuel pin cell.\"\"\"\n",
    "    c_fuel = openmc.Cell(fill=fuel, region=-fuel_cyl)\n",
    "    c_gap = openmc.Cell(fill=air, region=+fuel_cyl & -clad_in_cyl)\n",
    "    c_clad = openmc.Cell(fill=clad, region=+clad_in_cyl & -clad_out_cyl)\n",
    "    c_water = openmc.Cell(fill=water, region=+clad_out_cyl)\n",
    "    return openmc.Universe(cells=[c_fuel, c_gap, c_clad, c_water])\n",
    "\n",
    "\n",
    "def make_control_pin():\n",
    "    \"\"\"Control rod pin: B4C absorber instead of fuel, same clad/water.\"\"\"\n",
    "    c_b4c = openmc.Cell(fill=b4c, region=-fuel_cyl)\n",
    "    c_gap = openmc.Cell(fill=air, region=+fuel_cyl & -clad_in_cyl)\n",
    "    c_clad = openmc.Cell(fill=clad, region=+clad_in_cyl & -clad_out_cyl)\n",
    "    c_water = openmc.Cell(fill=water, region=+clad_out_cyl)\n",
    "    return openmc.Universe(cells=[c_b4c, c_gap, c_clad, c_water])\n",
    "\n",
    "\n",
    "def make_water_pin():\n",
    "    \"\"\"Pure water pin (guide tube / instrument location).\"\"\"\n",
    "    c_water = openmc.Cell(fill=water, region=+clad_out_cyl)\n",
    "    c_clad = openmc.Cell(fill=clad, region=+clad_in_cyl & -clad_out_cyl)\n",
    "    c_inner = openmc.Cell(fill=water, region=-clad_in_cyl)\n",
    "    return openmc.Universe(cells=[c_inner, c_clad, c_water])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96ba5bb5",
   "metadata": {},
   "outputs": [],
   "source": [
    "fuel_u = make_fuel_pin()\n",
    "ctrl_u = make_control_pin()\n",
    "wtr_u = make_water_pin()\n",
    "print(\"Pin-cell universes built: fuel, control rod, water guide tube\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8962bc0e",
   "metadata": {},
   "source": [
    "### 3c. The core pattern\n",
    "\n",
    "`F` = fuel pin, `C` = control rod pin, `W` = water guide tube.\n",
    "\n",
    "**This is the main thing you'll edit in Section 8 to run experiments \u2014\n",
    "try changing a few `C`s to `F`s later and see what happens to keff!**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "244527e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 7\n",
    "F, C, W = fuel_u, ctrl_u, wtr_u\n",
    "pattern = [\n",
    "    [F, F, F, C, F, F, F],\n",
    "    [F, F, F, C, F, F, F],\n",
    "    [F, F, W, C, W, F, F],\n",
    "    [C, C, C, C, C, C, C],\n",
    "    [F, F, W, C, W, F, F],\n",
    "    [F, F, F, C, F, F, F],\n",
    "    [F, F, F, C, F, F, F],\n",
    "]\n",
    "\n",
    "lattice = openmc.RectLattice(name=\"7x7 simplified assembly\")\n",
    "lattice.lower_left = (-N * pin_pitch / 2, -N * pin_pitch / 2)\n",
    "lattice.pitch = (pin_pitch, pin_pitch)\n",
    "lattice.universes = pattern"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7729ae52",
   "metadata": {},
   "source": [
    "### 3d. Core boundary, reflector, and axial extent\n",
    "A 5 cm water reflector surrounds the lattice; outside that is **vacuum**\n",
    "(real leakage). Top/bottom are **reflective**, approximating an\n",
    "infinitely tall core so this stays a fast 2D-ish problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b39d3493",
   "metadata": {},
   "outputs": [],
   "source": [
    "core_half = N * pin_pitch / 2\n",
    "reflector_half = core_half + 5.0\n",
    "\n",
    "core_box = openmc.model.RectangularPrism(\n",
    "    width=2 * core_half, height=2 * core_half, boundary_type='transmission'\n",
    ")\n",
    "reflector_box = openmc.model.RectangularPrism(\n",
    "    width=2 * reflector_half, height=2 * reflector_half, boundary_type='vacuum'\n",
    ")\n",
    "z_min = openmc.ZPlane(z0=-10, boundary_type='reflective')\n",
    "z_max = openmc.ZPlane(z0=10, boundary_type='reflective')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "59aef303",
   "metadata": {},
   "source": [
    "### 3e. Assemble and export the geometry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8a0e6c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "core_cell = openmc.Cell(fill=lattice, region=-core_box & +z_min & -z_max)\n",
    "reflector_cell = openmc.Cell(\n",
    "    fill=water, region=+core_box & -reflector_box & +z_min & -z_max\n",
    ")\n",
    "\n",
    "root_universe = openmc.Universe(cells=[core_cell, reflector_cell])\n",
    "geometry = openmc.Geometry(root_universe)\n",
    "geometry.export_to_xml()\n",
    "\n",
    "print(f\"Core half-width:      {core_half:.2f} cm\")\n",
    "print(f\"Reflector half-width: {reflector_half:.2f} cm\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b4e6422",
   "metadata": {},
   "source": [
    "---\n",
    "## 4. Settings: a short criticality (k-eigenvalue) run\n",
    "\n",
    "We use a small particle count on purpose (2000 per batch) so this runs in\n",
    "a couple of seconds \u2014 perfect for a live demo where you'll rerun it many\n",
    "times. A production calculation would use millions of particles for a\n",
    "tighter uncertainty on keff."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9a5802d",
   "metadata": {},
   "outputs": [],
   "source": [
    "settings = openmc.Settings()\n",
    "settings.batches = 60\n",
    "settings.inactive = 15\n",
    "settings.particles = 2000\n",
    "settings.run_mode = 'eigenvalue'\n",
    "\n",
    "bounds = [-core_half, -core_half, -9, core_half, core_half, 9]\n",
    "uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)\n",
    "settings.source = openmc.IndependentSource(space=uniform_dist)\n",
    "settings.export_to_xml()\n",
    "print(\"Settings exported: 60 batches (15 inactive), 2000 particles/batch\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e5ea3de",
   "metadata": {},
   "source": [
    "---\n",
    "## 5. Tallies: where do neutrons go, and where does power come from?\n",
    "\n",
    "We lay a 100\u00d7100 mesh over the whole model and ask OpenMC to record\n",
    "**flux** (how many neutrons pass through each cell) and **fission**\n",
    "(how much power is produced in each cell)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d4573ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = openmc.RegularMesh()\n",
    "mesh.dimension = [100, 100]\n",
    "mesh.lower_left = [-reflector_half, -reflector_half]\n",
    "mesh.upper_right = [reflector_half, reflector_half]\n",
    "mesh_filter = openmc.MeshFilter(mesh)\n",
    "\n",
    "tally = openmc.Tally(name=\"flux_and_fission\")\n",
    "tally.filters = [mesh_filter]\n",
    "tally.scores = ['flux', 'fission']\n",
    "\n",
    "tallies = openmc.Tallies([tally])\n",
    "tallies.export_to_xml()\n",
    "print(\"Tallies exported\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3990fa96",
   "metadata": {},
   "source": [
    "---\n",
    "## 6. Define the geometry plots\n",
    "\n",
    "Two 2D slices: one colored by **material**, one by **cell** \u2014 useful for\n",
    "checking the geometry looks the way you expect before spending time on a\n",
    "full transport run."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4eff88df",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_xy = openmc.Plot()\n",
    "plot_xy.filename = 'plot_xy_materials'\n",
    "plot_xy.width = (2 * reflector_half + 2, 2 * reflector_half + 2)\n",
    "plot_xy.pixels = (1000, 1000)\n",
    "plot_xy.color_by = 'material'\n",
    "plot_xy.colors = {\n",
    "    fuel: (255, 60, 60),\n",
    "    clad: (150, 150, 150),\n",
    "    water: (120, 190, 255),\n",
    "    b4c: (40, 40, 40),\n",
    "    air: (255, 255, 255),\n",
    "}\n",
    "plot_xy.basis = 'xy'\n",
    "\n",
    "plot_xy_cell = openmc.Plot()\n",
    "plot_xy_cell.filename = 'plot_xy_cells'\n",
    "plot_xy_cell.width = plot_xy.width\n",
    "plot_xy_cell.pixels = (1000, 1000)\n",
    "plot_xy_cell.color_by = 'cell'\n",
    "plot_xy_cell.basis = 'xy'\n",
    "\n",
    "plots = openmc.Plots([plot_xy, plot_xy_cell])\n",
    "plots.export_to_xml()\n",
    "openmc.plot_geometry()\n",
    "print(\"Geometry plots rendered to plot_xy_materials.png / plot_xy_cells.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "53f7ec0f",
   "metadata": {},
   "source": [
    "---\n",
    "## 7. Run the simulation\n",
    "\n",
    "This launches the actual Monte Carlo transport calculation. Watch the live\n",
    "tally output scroll by \u2014 each line is one \"batch\" of 2000 neutron\n",
    "histories. The last several (\"active\") batches are averaged to estimate\n",
    "keff."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc79d536",
   "metadata": {},
   "outputs": [],
   "source": [
    "openmc.run()\n",
    "sp_path = sorted(glob.glob('statepoint.*.h5'))[-1]\n",
    "print(\"\\nStatepoint written to:\", sp_path)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "434d9900",
   "metadata": {},
   "source": [
    "---\n",
    "## 8. Look at the results\n",
    "\n",
    "### 8a. Geometry\n",
    "Red=fuel, grey=clad, blue=water, black=B4C control rods, white=air gap."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "641564f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 2, figsize=(13, 6.5))\n",
    "for ax, fname, title in [\n",
    "    (axes[0], 'plot_xy_materials.png', 'Materials\\n(red=fuel, grey=clad, blue=water, black=B4C)'),\n",
    "    (axes[1], 'plot_xy_cells.png', 'Cells (pin-cell structure)'),\n",
    "]:\n",
    "    img = mpimg.imread(fname)\n",
    "    ax.imshow(img)\n",
    "    ax.set_title(title, fontsize=10)\n",
    "    ax.axis('off')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b8b76d6",
   "metadata": {},
   "source": [
    "### 8b. Flux and fission (power) maps, plus keff\n",
    "Flux shows where neutrons travel; fission shows where power is actually\n",
    "produced \u2014 notice how the control-rod cross suppresses both locally."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2695eb2",
   "metadata": {},
   "outputs": [],
   "source": [
    "with openmc.StatePoint(sp_path) as sp:\n",
    "    tally_result = sp.get_tally(name='flux_and_fission')\n",
    "    flux = tally_result.get_slice(scores=['flux']).mean.reshape(100, 100)\n",
    "    fission = tally_result.get_slice(scores=['fission']).mean.reshape(100, 100)\n",
    "    keff = sp.keff\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, figsize=(13, 5.5))\n",
    "im0 = axes[0].imshow(flux, origin='lower', cmap='viridis')\n",
    "axes[0].set_title('Neutron flux distribution')\n",
    "axes[0].axis('off')\n",
    "plt.colorbar(im0, ax=axes[0], fraction=0.046)\n",
    "\n",
    "im1 = axes[1].imshow(fission, origin='lower', cmap='hot')\n",
    "axes[1].set_title('Fission (power) distribution')\n",
    "axes[1].axis('off')\n",
    "plt.colorbar(im1, ax=axes[1], fraction=0.046)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f7eba56",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"k-effective = {keff.nominal_value:.5f} +/- {keff.std_dev:.5f}\")\n",
    "if keff.nominal_value > 1.0:\n",
    "    print(\"  -> SUPERCRITICAL (neutron population growing)\")\n",
    "elif keff.nominal_value < 1.0:\n",
    "    print(\"  -> SUBCRITICAL (neutron population shrinking)\")\n",
    "else:\n",
    "    print(\"  -> CRITICAL (steady chain reaction)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d8ae617",
   "metadata": {},
   "source": [
    "---\n",
    "## 9. \ud83e\uddea Now experiment!\n",
    "\n",
    "Pick ONE change below, re-run **from the relevant section downward**\n",
    "(noted per experiment) \u2014 not from the very top, that just redefines\n",
    "things that don't need to change \u2014 and see what happens to keff and the\n",
    "power map.\n",
    "\n",
    "Discuss with your group before changing: do you think keff will go **up**\n",
    "or **down**, and why?\n",
    "\n",
    "| # | Change | Where | Re-run from |\n",
    "|---|---|---|---|\n",
    "| a | Pull a control rod: change one `C` to `F` in the pattern grid | \u00a73c | Section 3 |\n",
    "| b | Change enrichment: `enrichment=4.0` \u2192 `2.0` or `5.0` | \u00a72a | Section 2 |\n",
    "| c | Remove boron: delete `water.add_element('B', 500e-6)` | \u00a72c | Section 2 |\n",
    "| d | Bigger core: `N = 7` \u2192 `9` or `11` (extend the pattern grid too!) | \u00a73c | Section 3 |\n",
    "| e | Water density: `set_density('g/cm3', 1.0)` instead of `0.74` | \u00a72c | Section 2 |\n",
    "\n",
    "After editing the target cell, just re-run every cell from that section\n",
    "through Section 8 again (Cell menu \u2192 \"Run All Below\" works well here)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46cfb223",
   "metadata": {},
   "source": [
    "---\n",
    "## 10. Wrap-up discussion\n",
    "\n",
    "- Why do real reactors need **active, continuous control** \u2014 not just\n",
    "  \"reach critical and done\"?\n",
    "- What's the difference between what we modeled (a static k-eigenvalue\n",
    "  snapshot) and a real reactor **transient** \u2014 time-dependent behavior,\n",
    "  fuel temperature feedback, xenon poisoning building up over hours? OpenMC\n",
    "  can do depletion/burnup calculations too \u2014 that's a \"part 2\" topic.\n",
    "- This model is a **teaching toy**: the pin dimensions and general pattern\n",
    "  are realistic PWR-ish numbers, but it is **not** a licensed, validated\n",
    "  reactor design. Real reactor physics calculations require validated\n",
    "  nuclear data, full 3D geometry, depletion, thermal-hydraulic feedback,\n",
    "  and rigorous QA \u2014 well beyond what fits in an hour."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13273258",
   "metadata": {},
   "source": "---\n### Troubleshooting\n- `ModuleNotFoundError: No module named 'openmc'` \u2192 Jupyter wasn't started\n  by `setup.sh`. Close it and run `./setup.sh` again.\n- \"Nuclear data MISSING\" in Section 0 \u2192 the notebook was opened from\n  outside the project folder. Open it from the folder that contains `xs_data/`.\n- Geometry plot is all one color \u2192 check the `colors` dict in Section 6\n  matches the materials you actually defined."
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "# Optional: download your simulation outputs before Colab disconnects.\n",
    "from google.colab import files\n",
    "from zipfile import ZipFile, ZIP_DEFLATED\n",
    "from pathlib import Path\n",
    "with ZipFile(\"/content/firefly-results.zip\", \"w\", ZIP_DEFLATED) as archive:\n",
    "    for pattern in (\"*.xml\", \"*.png\", \"statepoint.*.h5\", \"summary.h5\"):\n",
    "        for path in Path(\"/content/firefly-003\").glob(pattern):\n",
    "            archive.write(path, arcname=path.name)\n",
    "files.download(\"/content/firefly-results.zip\")\n"
   ],
   "execution_count": null,
   "outputs": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python"
  },
  "colab": {
   "name": "Firefly_003_Colab.ipynb",
   "provenance": []
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
