From fe3add40071e9ecb1c9f02cedd9061a4e7139c12 Mon Sep 17 00:00:00 2001
From: "p.vanderwilt" <p.van.der.wilt@brabantsedelta.nl>
Date: Fri, 4 Jul 2025 17:42:31 +0200
Subject: [PATCH] Add debug assertions for state changes in Reactor classes

---
 src/reactor_class.js | 31 +++++++++++++++++--------------
 1 file changed, 17 insertions(+), 14 deletions(-)

diff --git a/src/reactor_class.js b/src/reactor_class.js
index dd5daf7..d12e158 100644
--- a/src/reactor_class.js
+++ b/src/reactor_class.js
@@ -10,6 +10,7 @@ const math = create(all, config);
 
 const S_O_INDEX = 0;
 const NUM_SPECIES = 13;
+const DEBUG = false;
 
 class Reactor {
     /**
@@ -129,11 +130,12 @@ class Reactor_CSTR extends Reactor {
         const transfer = Array(NUM_SPECIES).fill(0.0);
         transfer[S_O_INDEX] = isNaN(this.kla) ? this.OTR : this._calcOTR(this.state[S_O_INDEX]); // calculate OTR if kla is not NaN, otherwise use externaly calculated OTR
 
-        const dC_total = math.multiply(math.add(inflow, outflow, reaction, transfer), time_step);
-        assertNoNaN(dC_total, "change in state");
-
+        const dC_total = math.multiply(math.add(inflow, outflow, reaction, transfer), time_step)
         this.state = this._arrayClip2Zero(math.add(this.state, dC_total)); // clip value element-wise to avoid negative concentrations
-        assertNoNaN(this.state, "new state");
+        if(DEBUG){
+          assertNoNaN(dC_total, "change in state");
+          assertNoNaN(this.state, "new state");
+        }
         return this.state;
     }
 }
@@ -194,8 +196,9 @@ class Reactor_PFR extends Reactor {
             const BC_gradient = Array(this.n_x).fill(0);
             BC_gradient[0] = -1;
             BC_gradient[1] = 1;
-            let Pe = this.length * math.sum(this.Fs) / (this.D * this.A)
-            const BC_dispersion = math.multiply((1 - (1 + 4*this.volume/math.sum(this.Fs)/Pe)^0.5) / Pe, [BC_gradient], state)[0];
+            let Pe = this.length * math.sum(this.Fs) / (this.D * this.A);
+            let residence_time = this.volume/math.sum(this.Fs);
+            const BC_dispersion = math.multiply((1 - (1 + 4*residence_time/Pe)^0.5) / (Pe*this.d_x), [BC_gradient], state)[0];
             state[0] = math.add(BC_C_in, BC_dispersion).map(val => val < 0 ? 0 : val);
         } else { // Neumann BC (no flux)
             state[0] = state[1];
@@ -215,10 +218,6 @@ class Reactor_PFR extends Reactor {
         const reaction = this.state.map((state_slice) => this.asm.compute_dC(state_slice));
         const transfer = Array.from(Array(this.n_x), () => new Array(NUM_SPECIES).fill(0));
 
-        assertNoNaN(dispersion, "dispersion");
-        assertNoNaN(advection, "advection");
-        assertNoNaN(reaction, "reaction");
-
         if (isNaN(this.kla)) { // calculate OTR if kla is not NaN, otherwise use externally calculated OTR
             transfer.forEach((x) => { x[S_O_INDEX] = this.OTR; });
         } else {
@@ -226,13 +225,17 @@ class Reactor_PFR extends Reactor {
         }
 
         const dC_total = math.multiply(math.add(dispersion, advection, reaction, transfer), time_step);
-        assertNoNaN(dC_total, "change in state");
 
         const stateNew = math.add(this.state, dC_total);
-        assertNoNaN(stateNew, "new state");
-
         this._applyBoundaryConditions(stateNew);
-        assertNoNaN(stateNew, "new state post BC");
+
+        if (DEBUG) {
+          assertNoNaN(dispersion, "dispersion");
+          assertNoNaN(advection, "advection");
+          assertNoNaN(reaction, "reaction");
+          assertNoNaN(dC_total, "change in state");
+          assertNoNaN(stateNew, "new state post BC");
+        }
 
         this.state = this._arrayClip2Zero(stateNew);
         return stateNew;