point kinetics

docs/source/examples/point_kinetics.ipynb

@@ -0,0 +1,152 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Reactor Point Kinetics\n",
+    "\n",
+    "Simulating the transient neutron population in a nuclear reactor using the point kinetics equations (PKE) with six delayed neutron precursor groups."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The point kinetics model couples the neutron density $n$ with $G$ delayed neutron precursor groups $C_i$:\n",
+    "\n",
+    "$$\\frac{dn}{dt} = \\frac{\\rho - \\beta}{\\Lambda} \\, n + \\sum_{i=1}^{G} \\lambda_i \\, C_i + S$$\n",
+    "\n",
+    "$$\\frac{dC_i}{dt} = \\frac{\\beta_i}{\\Lambda} \\, n - \\lambda_i \\, C_i \\qquad i = 1, \\ldots, G$$\n",
+    "\n",
+    "where $\\beta = \\sum_i \\beta_i$ is the total delayed neutron fraction, $\\Lambda$ is the prompt neutron generation time, and $\\rho$ is the reactivity."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": "import matplotlib.pyplot as plt\n\nplt.style.use('../pathsim_docs.mplstyle')\n\nfrom pathsim import Simulation, Connection\nfrom pathsim.blocks import Source, Scope\nfrom pathsim.solvers import GEAR52A\n\nfrom pathsim_chem.neutronics import PointKinetics"
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": "The point kinetics system is very stiff — the prompt neutron generation time $\\Lambda \\sim 10^{-5}\\,\\text{s}$ creates eigenvalues on the order of $10^5$. The variable-order BDF solver GEAR52A is ideal here: it requires only one implicit solve per step and adapts both step size and order to the smooth exponential dynamics. The default fixed-point tolerance (`1e-9`) is unnecessarily tight for this problem; relaxing it to `1e-6` gives a ~70x speedup with negligible loss in accuracy.\n\n## 1. Delayed Supercritical Step\n\nInsert a step reactivity of $\\rho = 0.003$ (about $0.46\\beta$). Since $\\rho < \\beta$, the reactor is delayed supercritical — the power rises on a slow time scale governed by the delayed neutrons."
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": "reactor = PointKinetics(n0=1.0)\n\nrho_step = 0.003\nsrc_rho = Source(lambda t: rho_step if t > 0 else 0.0)\nsrc_s = Source(lambda t: 0.0)\nsco = Scope(labels=[\"n\"])\n\nsim = Simulation(\n    [src_rho, src_s, reactor, sco],\n    [\n        Connection(src_rho, reactor),       # rho\n        Connection(src_s, reactor[1]),       # S\n        Connection(reactor, sco),            # n\n    ],\n    dt=0.01,\n    Solver=GEAR52A,\n    tolerance_fpi=1e-6,\n    log=True,\n)\n\nsim.run(100)"
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sco.plot(lw=1.5)\n",
+    "plt.yscale('log')\n",
+    "plt.title('Delayed Supercritical (ρ = 0.003)')\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The neutron density rises exponentially on a time scale of seconds — much slower than the prompt neutron lifetime ($\\Lambda \\sim 10^{-5}$ s) because the delayed neutrons control the dynamics when $\\rho < \\beta$."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## 2. Prompt Supercritical\n",
+    "\n",
+    "Insert $\\rho = 0.008 > \\beta \\approx 0.0065$. Now the reactor is prompt supercritical — the power rises on the prompt neutron time scale, producing a rapid excursion."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": "reactor = PointKinetics(n0=1.0)\n\nrho_prompt = 0.008\nsrc_rho = Source(lambda t: rho_prompt if t > 0 else 0.0)\nsrc_s = Source(lambda t: 0.0)\nsco = Scope(labels=[\"n\"])\n\nsim = Simulation(\n    [src_rho, src_s, reactor, sco],\n    [\n        Connection(src_rho, reactor),\n        Connection(src_s, reactor[1]),\n        Connection(reactor, sco),\n    ],\n    dt=0.001,\n    Solver=GEAR52A,\n    tolerance_fpi=1e-6,\n    log=True,\n)\n\nsim.run(0.5)"
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sco.plot(lw=1.5)\n",
+    "plt.yscale('log')\n",
+    "plt.title('Prompt Supercritical (ρ = 0.008 > β)')\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The power rises orders of magnitude within milliseconds. This is why prompt criticality must be avoided in reactor design — the delayed neutrons are the key safety mechanism that keeps power transients manageable."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## 3. Subcritical with External Source\n",
+    "\n",
+    "A subcritical assembly ($\\rho = -0.05$) with a constant external neutron source. The system reaches an equilibrium where the source multiplication produces a steady neutron population:\n",
+    "\n",
+    "$$n_{ss} = \\frac{-S \\, \\Lambda}{\\rho}$$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": "rho_sub = -0.05\ns_ext = 1e5\nLambda = 1e-5\n\nreactor = PointKinetics(n0=0.001, Lambda=Lambda)\n\nsrc_rho = Source(lambda t: rho_sub)\nsrc_s = Source(lambda t: s_ext)\nsco = Scope(labels=[\"n\"])\n\nsim = Simulation(\n    [src_rho, src_s, reactor, sco],\n    [\n        Connection(src_rho, reactor),\n        Connection(src_s, reactor[1]),\n        Connection(reactor, sco),\n    ],\n    dt=0.01,\n    Solver=GEAR52A,\n    tolerance_fpi=1e-6,\n    log=True,\n)\n\nsim.run(50)"
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "n_ss = -s_ext * Lambda / rho_sub\n",
+    "\n",
+    "sco.plot(lw=1.5)\n",
+    "plt.axhline(n_ss, color='r', ls='--', lw=1, label=f'$n_{{ss}}$ = {n_ss:.4f}')\n",
+    "plt.legend()\n",
+    "plt.title('Subcritical with Source (ρ = −0.05)')\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The neutron density converges to the expected source-multiplied equilibrium value, confirming the subcritical multiplication physics."
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "name": "python",
+   "version": "3.11.0"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}

src/pathsim_chem/__init__.py

@@ -16,3 +16,4 @@
 from .tritium import *
 from .thermodynamics import *
 from .process import *
+from .neutronics import *

src/pathsim_chem/neutronics/__init__.py

@@ -0,0 +1 @@
+from .point_kinetics import *

src/pathsim_chem/neutronics/point_kinetics.py

@@ -0,0 +1,146 @@
+#########################################################################################
+##
+##                         Reactor Point Kinetics Block
+##
+#########################################################################################
+
+# IMPORTS ===============================================================================
+
+import numpy as np
+
+from pathsim.blocks.dynsys import DynamicalSystem
+
+# CONSTANTS =============================================================================
+
+# Keepin U-235 thermal 6-group delayed neutron data
+BETA_U235 = [0.000215, 0.001424, 0.001274, 0.002568, 0.000748, 0.000273]
+LAMBDA_U235 = [0.0124, 0.0305, 0.111, 0.301, 1.14, 3.01]
+
+# BLOCKS ================================================================================
+
+class PointKinetics(DynamicalSystem):
+    """Reactor point kinetics equations with delayed neutron precursors.
+
+    Models the time-dependent neutron population in a nuclear reactor using
+    the point kinetics equations (PKE). The neutron density responds to
+    reactivity changes through prompt fission and delayed neutron emission
+    from G precursor groups.
+
+    Mathematical Formulation
+    -------------------------
+    The state vector is :math:`[n, C_1, C_2, \\ldots, C_G]` where :math:`n`
+    is the neutron density (or power) and :math:`C_i` are the delayed neutron
+    precursor concentrations.
+
+    .. math::
+
+        \\frac{dn}{dt} = \\frac{\\rho - \\beta}{\\Lambda} \\, n
+            + \\sum_{i=1}^{G} \\lambda_i \\, C_i + S
+
+    .. math::
+
+        \\frac{dC_i}{dt} = \\frac{\\beta_i}{\\Lambda} \\, n
+            - \\lambda_i \\, C_i \\qquad i = 1, \\ldots, G
+
+    where :math:`\\beta = \\sum_i \\beta_i` is the total delayed neutron
+    fraction and :math:`\\Lambda` is the prompt neutron generation time.
+
+    Parameters
+    ----------
+    n0 : float
+        Initial neutron density [-]. Default 1.0 (normalized).
+    Lambda : float
+        Prompt neutron generation time [s].
+    beta : array_like
+        Delayed neutron fractions per group [-].
+    lam : array_like
+        Precursor decay constants per group [1/s].
+    """
+
+    input_port_labels = {
+        "rho": 0,
+        "S":   1,
+    }
+
+    output_port_labels = {
+        "n": 0,
+    }
+
+    def __init__(self, n0=1.0, Lambda=1e-5,
+                 beta=None, lam=None):
+
+        # defaults
+        if beta is None:
+            beta = BETA_U235
+        if lam is None:
+            lam = LAMBDA_U235
+
+        beta = np.asarray(beta, dtype=float)
+        lam = np.asarray(lam, dtype=float)
+
+        # input validation
+        if Lambda <= 0:
+            raise ValueError(f"'Lambda' must be positive but is {Lambda}")
+        if len(beta) != len(lam):
+            raise ValueError(
+                f"'beta' and 'lam' must have the same length "
+                f"but got {len(beta)} and {len(lam)}"
+            )
+        if len(beta) == 0:
+            raise ValueError("'beta' and 'lam' must have at least one group")
+
+        # store parameters
+        G = len(beta)
+        beta_total = float(np.sum(beta))
+
+        self.n0 = n0
+        self.Lambda = Lambda
+        self.beta = beta
+        self.lam = lam
+        self.G = G
+        self.beta_total = beta_total
+
+        # initial conditions: steady state at rho = 0
+        x0 = np.empty(G + 1)
+        x0[0] = n0
+        for i in range(G):
+            x0[1 + i] = beta[i] / (Lambda * lam[i]) * n0
+
+        # rhs of point kinetics ode system
+        def _fn_d(x, u, t):
+            n = x[0]
+            C = x[1:]
+
+            rho = u[0]
+            S = u[1] if len(u) > 1 else 0.0
+
+            dn = (rho - beta_total) / Lambda * n + np.dot(lam, C) + S
+
+            dC = beta / Lambda * n - lam * C
+
+            return np.concatenate(([dn], dC))
+
+        # jacobian of rhs wrt state [n, C_1, ..., C_G]
+        def _jc_d(x, u, t):
+            rho = u[0]
+
+            J = np.zeros((G + 1, G + 1))
+
+            J[0, 0] = (rho - beta_total) / Lambda
+            for i in range(G):
+                J[0, 1 + i] = lam[i]
+                J[1 + i, 0] = beta[i] / Lambda
+                J[1 + i, 1 + i] = -lam[i]
+
+            return J
+
+        # output function: neutron density only
+        def _fn_a(x, u, t):
+            return x[:1].copy()
+
+        super().__init__(
+            func_dyn=_fn_d,
+            jac_dyn=_jc_d,
+            func_alg=_fn_a,
+            initial_value=x0,
+        )