From 4b49d1076368015a85b2eab2b38cee006aba4475 Mon Sep 17 00:00:00 2001
From: "p.vanderwilt" <p.van.der.wilt@brabantsedelta.nl>
Date: Fri, 4 Jul 2025 15:48:05 +0200
Subject: [PATCH] Fixed bug in NaN assertion

---
 src/reactor_class.js | 52 +++++++++++++++++++++++++-------------------
 1 file changed, 30 insertions(+), 22 deletions(-)

diff --git a/src/reactor_class.js b/src/reactor_class.js
index d483201..d41fd39 100644
--- a/src/reactor_class.js
+++ b/src/reactor_class.js
@@ -13,8 +13,14 @@ const math = create(all, config)
  * @param {string} label 
  */
 function assertNoNaN(arr, label = "array") {
-    if (math.isNaN(arr)) {
-        throw new Error("NaN detected in ${label}!");
+    if (Array.isArray(arr)) {
+        for (const el of arr) {
+            assertNoNaN(el, label);
+        }
+    } else {
+        if (Number.isNaN(arr)) {
+            throw new Error(`NaN detected in ${label}!`);
+        }
     }
 }
 
@@ -80,7 +86,7 @@ class Reactor {
      */
     _arrayClip2Zero(arr) {
         if (Array.isArray(arr)) {
-            return arr.map(clipToZero);
+            return arr.map(x => this._arrayClip2Zero(x));
         } else {
             return arr < 0 ? 0 : arr;
         }
@@ -166,6 +172,9 @@ class Reactor_PFR extends Reactor {
 
         this.D_op = this._makeDoperator(true, true);
         this.D2_op = this._makeD2operator();
+
+        assertNoNaN(this.D_op, "Derivative operator");
+        assertNoNaN(this.D2_op, "Second derivative operator");
     }
 
     /**
@@ -184,7 +193,7 @@ class Reactor_PFR extends Reactor {
         return { topic: "Fluent", payload: { inlet: 0, F: math.sum(this.Fs), C: this.state.at(-1) }, timestamp: this.currentTime };
     }
 
-    _applyBoundaryConditions(newState) {
+    _applyBoundaryConditions(state) {
         // apply boundary conditions
         if (math.sum(this.Fs) > 0) { // Danckwerts BC
             const BC_C_in = math.multiply(1 / math.sum(this.Fs), [this.Fs], this.Cs_in)[0];
@@ -192,14 +201,13 @@ class Reactor_PFR extends Reactor {
             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], stateNew)[0];
-            newState[0] = math.add(BC_C_in, BC_dispersion).map(val => val < 0 ? 0 : val);
+            const BC_dispersion = math.multiply((1 - (1 + 4 * this.volume / math.sum(this.Fs) / Pe) ^ 0.5) / Pe, [BC_gradient], state)[0];
+            state[0] = math.add(BC_C_in, BC_dispersion).map(val => val < 0 ? 0 : val);
         } else { // Neumann BC (no flux)
-            newState[0] = newState[1];
+            state[0] = state[1];
         }
         // Neumann BC (no flux)
-        newState[this.n_x - 1] = newState[this.n_x - 2]
-        return newState
+        state[this.n_x - 1] = state[this.n_x - 2]
     }
 
     /**
@@ -224,11 +232,11 @@ class Reactor_PFR extends Reactor {
         }
 
         const dC_total = math.multiply(math.add(dispersion, advection, reaction, transfer), time_step);
-        let stateNew = math.add(this.state, dC_total);
+        const stateNew = math.add(this.state, dC_total);
 
         assertNoNaN(stateNew, "new state");
 
-        stateNew = this._applyBoundaryConditions(stateNew);
+        this._applyBoundaryConditions(stateNew);
 
         assertNoNaN(stateNew, "new state post BC");
 
@@ -245,22 +253,22 @@ class Reactor_PFR extends Reactor {
     _makeDoperator(central = false, higher_order = false) { // create gradient operator
         if (higher_order) {
             if (central) {
-                const I = math.resize(math.diag(Array(this.n_x).fill(1 / 12), -2), [this.n_x, this.n_x]);
-                const A = math.resize(math.diag(Array(this.n_x).fill(-2 / 3), -1), [this.n_x, this.n_x]);
-                const B = math.resize(math.diag(Array(this.n_x).fill(2 / 3), 1), [this.n_x, this.n_x]);
-                const C = math.resize(math.diag(Array(this.n_x).fill(-1 / 12), 2), [this.n_x, this.n_x]);
+                const I = math.resize(math.diag(Array(this.n_x).fill(1/12), -2), [this.n_x, this.n_x]);
+                const A = math.resize(math.diag(Array(this.n_x).fill(-2/3), -1), [this.n_x, this.n_x]);
+                const B = math.resize(math.diag(Array(this.n_x).fill(2/3), 1), [this.n_x, this.n_x]);
+                const C = math.resize(math.diag(Array(this.n_x).fill(-1/12), 2), [this.n_x, this.n_x]);
                 const D = math.add(I, A, B, C);
                 const NearBoundary = Array(this.n_x).fill(0.0);
-                NearBoundary[0] = -1 / 4;
-                NearBoundary[1] = -5 / 6;
-                NearBoundary[2] = 3 / 2;
-                NearBoundary[3] = -1 / 2;
-                NearBoundary[4] = 1 / 12;
+                NearBoundary[0] = -1/4;
+                NearBoundary[1] = -5/6;
+                NearBoundary[2] = 3/2;
+                NearBoundary[3] = -1/2;
+                NearBoundary[4] = 1/12;
                 D[1] = NearBoundary;
                 NearBoundary.reverse();
-                D[this.n_x - 2] = math.multiply(-1, NearBoundary);
+                D[this.n_x-2] = math.multiply(-1, NearBoundary);
                 D[0] = Array(this.n_x).fill(0); // set by BCs elsewhere
-                D[this.n_x - 1] = Array(this.n_x).fill(0);
+                D[this.n_x-1] = Array(this.n_x).fill(0);
                 return D;
             } else {
                 throw new Error("Upwind higher order method not implemented! Use central scheme instead.");