From 5213c098247f5c4d4b1f2904c6ed615ecb372e3a Mon Sep 17 00:00:00 2001 From: Milan Rother Date: Mon, 2 Mar 2026 15:46:24 +0100 Subject: [PATCH 1/4] Add neutronics module with PointKinetics block --- docs/source/examples/point_kinetics.ipynb | 228 ++++++++++++++++++ src/pathsim_chem/__init__.py | 1 + src/pathsim_chem/neutronics/__init__.py | 1 + src/pathsim_chem/neutronics/point_kinetics.py | 146 +++++++++++ tests/neutronics/__init__.py | 0 tests/neutronics/test_point_kinetics.py | 146 +++++++++++ 6 files changed, 522 insertions(+) create mode 100644 docs/source/examples/point_kinetics.ipynb create mode 100644 src/pathsim_chem/neutronics/__init__.py create mode 100644 src/pathsim_chem/neutronics/point_kinetics.py create mode 100644 tests/neutronics/__init__.py create mode 100644 tests/neutronics/test_point_kinetics.py diff --git a/docs/source/examples/point_kinetics.ipynb b/docs/source/examples/point_kinetics.ipynb new file mode 100644 index 0000000..f115068 --- /dev/null +++ b/docs/source/examples/point_kinetics.ipynb @@ -0,0 +1,228 @@ +{ + "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", + "\n", + "plt.style.use('../pathsim_docs.mplstyle')\n", + "\n", + "from pathsim import Simulation, Connection\n", + "from pathsim.blocks import Source, Scope\n", + "\n", + "from pathsim_chem.neutronics import PointKinetics" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 1. Delayed Supercritical Step\n", + "\n", + "Insert 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", + "\n", + "rho_step = 0.003\n", + "src_rho = Source(lambda t: rho_step if t > 0 else 0.0)\n", + "src_s = Source(lambda t: 0.0)\n", + "sco = Scope(labels=[\"n\"])\n", + "\n", + "sim = 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", + " log=True,\n", + ")\n", + "\n", + "sim.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", + "\n", + "rho_prompt = 0.008\n", + "src_rho = Source(lambda t: rho_prompt if t > 0 else 0.0)\n", + "src_s = Source(lambda t: 0.0)\n", + "sco = Scope(labels=[\"n\"])\n", + "\n", + "sim = 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", + " log=True,\n", + ")\n", + "\n", + "sim.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\n", + "s_ext = 1e5\n", + "Lambda = 1e-5\n", + "\n", + "reactor = PointKinetics(n0=0.001, Lambda=Lambda)\n", + "\n", + "src_rho = Source(lambda t: rho_sub)\n", + "src_s = Source(lambda t: s_ext)\n", + "sco = Scope(labels=[\"n\"])\n", + "\n", + "sim = 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", + " log=True,\n", + ")\n", + "\n", + "sim.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 +} diff --git a/src/pathsim_chem/__init__.py b/src/pathsim_chem/__init__.py index 39dc189..51a8900 100644 --- a/src/pathsim_chem/__init__.py +++ b/src/pathsim_chem/__init__.py @@ -16,3 +16,4 @@ from .tritium import * from .thermodynamics import * from .process import * +from .neutronics import * diff --git a/src/pathsim_chem/neutronics/__init__.py b/src/pathsim_chem/neutronics/__init__.py new file mode 100644 index 0000000..2b92bfd --- /dev/null +++ b/src/pathsim_chem/neutronics/__init__.py @@ -0,0 +1 @@ +from .point_kinetics import * diff --git a/src/pathsim_chem/neutronics/point_kinetics.py b/src/pathsim_chem/neutronics/point_kinetics.py new file mode 100644 index 0000000..c32a264 --- /dev/null +++ b/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, + ) diff --git a/tests/neutronics/__init__.py b/tests/neutronics/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/tests/neutronics/test_point_kinetics.py b/tests/neutronics/test_point_kinetics.py new file mode 100644 index 0000000..bf6492c --- /dev/null +++ b/tests/neutronics/test_point_kinetics.py @@ -0,0 +1,146 @@ +######################################################################################## +## +## TESTS FOR +## 'neutronics.point_kinetics.py' +## +######################################################################################## + +# IMPORTS ============================================================================== + +import unittest +import numpy as np + +from pathsim_chem.neutronics import PointKinetics +from pathsim_chem.neutronics.point_kinetics import BETA_U235, LAMBDA_U235 + +from pathsim.solvers import EUF + + +# TESTS ================================================================================ + +class TestPointKinetics(unittest.TestCase): + """Test the PointKinetics (reactor point kinetics equations) block.""" + + def test_init_default(self): + """Test default initialization with Keepin U-235 data.""" + pk = PointKinetics() + self.assertEqual(pk.n0, 1.0) + self.assertEqual(pk.Lambda, 1e-5) + self.assertEqual(pk.G, 6) + np.testing.assert_array_equal(pk.beta, BETA_U235) + np.testing.assert_array_equal(pk.lam, LAMBDA_U235) + self.assertAlmostEqual(pk.beta_total, sum(BETA_U235)) + + def test_init_custom(self): + """Test custom initialization with user-specified groups.""" + beta = [0.003, 0.004] + lam = [0.1, 1.0] + pk = PointKinetics(n0=2.0, Lambda=1e-4, beta=beta, lam=lam) + + self.assertEqual(pk.n0, 2.0) + self.assertEqual(pk.Lambda, 1e-4) + self.assertEqual(pk.G, 2) + np.testing.assert_array_almost_equal(pk.beta, beta) + np.testing.assert_array_almost_equal(pk.lam, lam) + + pk.set_solver(EUF, parent=None) + x0 = pk.engine.initial_value + self.assertEqual(len(x0), 3) + self.assertAlmostEqual(x0[0], 2.0) + # C_i(0) = beta_i / (Lambda * lam_i) * n0 + self.assertAlmostEqual(x0[1], 0.003 / (1e-4 * 0.1) * 2.0) + self.assertAlmostEqual(x0[2], 0.004 / (1e-4 * 1.0) * 2.0) + + def test_init_validation(self): + """Test input validation.""" + with self.assertRaises(ValueError): + PointKinetics(Lambda=0) + with self.assertRaises(ValueError): + PointKinetics(Lambda=-1e-5) + with self.assertRaises(ValueError): + PointKinetics(beta=[0.001, 0.002], lam=[0.1]) + with self.assertRaises(ValueError): + PointKinetics(beta=[], lam=[]) + + def test_port_labels(self): + """Test port label definitions.""" + self.assertEqual(PointKinetics.input_port_labels["rho"], 0) + self.assertEqual(PointKinetics.input_port_labels["S"], 1) + self.assertEqual(PointKinetics.output_port_labels["n"], 0) + + def test_state_size(self): + """State vector has G+1 entries (neutron density + G precursor groups).""" + for G in [1, 2, 6]: + beta = [0.001] * G + lam = [0.1] * G + pk = PointKinetics(beta=beta, lam=lam) + pk.set_solver(EUF, parent=None) + self.assertEqual(len(pk.engine.initial_value), G + 1) + + def test_steady_state(self): + """At rho=0 and S=0, the initial conditions are steady state (dn/dt ~ 0).""" + pk = PointKinetics() + pk.set_solver(EUF, parent=None) + + x0 = pk.engine.initial_value.copy() + u = np.array([0.0, 0.0]) # rho=0, S=0 + dx = pk.op_dyn(x0, u, 0) + + np.testing.assert_allclose(dx, 0.0, atol=1e-10) + + def test_precursor_equilibrium(self): + """Verify initial precursor concentrations satisfy equilibrium.""" + n0 = 5.0 + Lambda = 2e-5 + beta = [0.0003, 0.001, 0.002] + lam = [0.05, 0.5, 2.0] + + pk = PointKinetics(n0=n0, Lambda=Lambda, beta=beta, lam=lam) + pk.set_solver(EUF, parent=None) + + x0 = pk.engine.initial_value + for i in range(3): + expected = beta[i] / (Lambda * lam[i]) * n0 + self.assertAlmostEqual(x0[1 + i], expected, places=8) + + def test_subcritical_source(self): + """With rho<0 and external source, system converges to steady state. + + At equilibrium: n_ss ≈ S * Lambda / (beta - rho) (for |rho| >> 0) + """ + beta = [0.003, 0.004] + lam = [0.1, 1.0] + Lambda = 1e-4 + rho = -0.05 + S = 1e6 + beta_total = sum(beta) + + # Start from source-free equilibrium, then apply source + negative rho + pk = PointKinetics(n0=0.001, Lambda=Lambda, beta=beta, lam=lam) + pk.set_solver(EUF, parent=None) + + # Expected steady state: dn/dt=0, dC_i/dt=0 + # From dC_i/dt=0: C_i = beta_i/(Lambda*lam_i) * n_ss + # Substituting into dn/dt=0: + # 0 = (rho-beta)/Lambda * n_ss + sum(lam_i * beta_i/(Lambda*lam_i) * n_ss) + S + # 0 = (rho-beta)/Lambda * n_ss + beta/Lambda * n_ss + S + # 0 = rho/Lambda * n_ss + S + # n_ss = -S * Lambda / rho + n_ss = -S * Lambda / rho + self.assertGreater(n_ss, 0) # sanity check: should be positive for rho < 0 + + # Verify by plugging n_ss into the equations directly + x_ss = np.empty(3) + x_ss[0] = n_ss + for i in range(2): + x_ss[1 + i] = beta[i] / (Lambda * lam[i]) * n_ss + + u = np.array([rho, S]) + dx = pk.op_dyn(x_ss, u, 0) + np.testing.assert_allclose(dx, 0.0, atol=1e-6) + + +# RUN TESTS LOCALLY ==================================================================== + +if __name__ == '__main__': + unittest.main(verbosity=2) From 52f0578ab3fd0e0d88bf5c96670d81c8cea38a39 Mon Sep 17 00:00:00 2001 From: Milan Rother Date: Mon, 2 Mar 2026 15:53:50 +0100 Subject: [PATCH 2/4] Use implicit solver (ESDIRK54) in point kinetics notebook for stiff PKE system --- docs/source/examples/point_kinetics.ipynb | 88 ++--------------------- 1 file changed, 6 insertions(+), 82 deletions(-) diff --git a/docs/source/examples/point_kinetics.ipynb b/docs/source/examples/point_kinetics.ipynb index f115068..baa4bd2 100644 --- a/docs/source/examples/point_kinetics.ipynb +++ b/docs/source/examples/point_kinetics.ipynb @@ -27,52 +27,19 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": [ - "import matplotlib.pyplot as plt\n", - "\n", - "plt.style.use('../pathsim_docs.mplstyle')\n", - "\n", - "from pathsim import Simulation, Connection\n", - "from pathsim.blocks import Source, Scope\n", - "\n", - "from pathsim_chem.neutronics import PointKinetics" - ] + "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 ESDIRK54\n\nfrom pathsim_chem.neutronics import PointKinetics" }, { "cell_type": "markdown", "metadata": {}, - "source": [ - "## 1. Delayed Supercritical Step\n", - "\n", - "Insert 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." - ] + "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$. An implicit solver (here ESDIRK54) is essential for stable integration at practical time steps.\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", - "\n", - "rho_step = 0.003\n", - "src_rho = Source(lambda t: rho_step if t > 0 else 0.0)\n", - "src_s = Source(lambda t: 0.0)\n", - "sco = Scope(labels=[\"n\"])\n", - "\n", - "sim = 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", - " log=True,\n", - ")\n", - "\n", - "sim.run(100)" - ] + "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=ESDIRK54,\n log=True,\n)\n\nsim.run(100)" }, { "cell_type": "code", @@ -107,27 +74,7 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": [ - "reactor = PointKinetics(n0=1.0)\n", - "\n", - "rho_prompt = 0.008\n", - "src_rho = Source(lambda t: rho_prompt if t > 0 else 0.0)\n", - "src_s = Source(lambda t: 0.0)\n", - "sco = Scope(labels=[\"n\"])\n", - "\n", - "sim = 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", - " log=True,\n", - ")\n", - "\n", - "sim.run(0.5)" - ] + "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=ESDIRK54,\n log=True,\n)\n\nsim.run(0.5)" }, { "cell_type": "code", @@ -164,30 +111,7 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": [ - "rho_sub = -0.05\n", - "s_ext = 1e5\n", - "Lambda = 1e-5\n", - "\n", - "reactor = PointKinetics(n0=0.001, Lambda=Lambda)\n", - "\n", - "src_rho = Source(lambda t: rho_sub)\n", - "src_s = Source(lambda t: s_ext)\n", - "sco = Scope(labels=[\"n\"])\n", - "\n", - "sim = 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", - " log=True,\n", - ")\n", - "\n", - "sim.run(50)" - ] + "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=ESDIRK54,\n log=True,\n)\n\nsim.run(50)" }, { "cell_type": "code", @@ -225,4 +149,4 @@ }, "nbformat": 4, "nbformat_minor": 4 -} +} \ No newline at end of file From d90bc85edd36a1729f01403621f54dbca4b25a28 Mon Sep 17 00:00:00 2001 From: Milan Rother Date: Mon, 2 Mar 2026 16:08:36 +0100 Subject: [PATCH 3/4] Switch point kinetics notebook to GEAR52A solver for stiff PKE system --- docs/source/examples/point_kinetics.ipynb | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/docs/source/examples/point_kinetics.ipynb b/docs/source/examples/point_kinetics.ipynb index baa4bd2..6570394 100644 --- a/docs/source/examples/point_kinetics.ipynb +++ b/docs/source/examples/point_kinetics.ipynb @@ -27,19 +27,19 @@ "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 ESDIRK54\n\nfrom pathsim_chem.neutronics import PointKinetics" + "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$. An implicit solver (here ESDIRK54) is essential for stable integration at practical time steps.\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." + "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.\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=ESDIRK54,\n log=True,\n)\n\nsim.run(100)" + "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 log=True,\n)\n\nsim.run(100)" }, { "cell_type": "code", @@ -74,7 +74,7 @@ "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=ESDIRK54,\n log=True,\n)\n\nsim.run(0.5)" + "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 log=True,\n)\n\nsim.run(0.5)" }, { "cell_type": "code", @@ -111,7 +111,7 @@ "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=ESDIRK54,\n log=True,\n)\n\nsim.run(50)" + "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 log=True,\n)\n\nsim.run(50)" }, { "cell_type": "code", From bc5b399456d02cce29f500e1528b58e47174a1c4 Mon Sep 17 00:00:00 2001 From: Milan Rother Date: Mon, 2 Mar 2026 16:13:50 +0100 Subject: [PATCH 4/4] Relax fixed-point tolerance to 1e-6 in point kinetics notebook for 70x speedup --- docs/source/examples/point_kinetics.ipynb | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/source/examples/point_kinetics.ipynb b/docs/source/examples/point_kinetics.ipynb index 6570394..d2a2b4d 100644 --- a/docs/source/examples/point_kinetics.ipynb +++ b/docs/source/examples/point_kinetics.ipynb @@ -32,14 +32,14 @@ { "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.\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." + "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 log=True,\n)\n\nsim.run(100)" + "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", @@ -74,7 +74,7 @@ "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 log=True,\n)\n\nsim.run(0.5)" + "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", @@ -111,7 +111,7 @@ "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 log=True,\n)\n\nsim.run(50)" + "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",