Refactor reactor class to improve NaN handling and add utility function for NaN assertions

src/reactor_class.js

@@ -1,5 +1,6 @@
 const ASM3 = require('./asm3_class.js');
 const { create, all } = require('mathjs');
+const { assertNoNaN } = require('./utils.js');
 
 const config = {
     matrix: 'Array' // use Array as the matrix type
@@ -7,26 +8,9 @@ const config = {
 
 const math = create(all, config);
 
-const OXYGEN_INDEX = 0;
+const S_O_INDEX = 0;
 const NUM_SPECIES = 13;
 
-/**
- * Assert that no NaN values are present in an array.
- * @param {Array} arr 
- * @param {string} label 
- */
-function assertNoNaN(arr, label = "array") {
-    if (Array.isArray(arr)) {
-        for (const el of arr) {
-            assertNoNaN(el, label);
-        }
-    } else {
-        if (Number.isNaN(arr)) {
-            throw new Error(`NaN detected in ${label}!`);
-        }
-    }
-}
-
 class Reactor {
     /**
      * Reactor base class.
@@ -35,16 +19,16 @@ class Reactor {
     constructor(config) {
         this.asm = new ASM3();
 
-        this.Vl = config.volume; // fluid volume reactor [m3]
+        this.volume = config.volume; // fluid volume reactor [m3]
 
-        this.Fs = Array(config.n_inlets).fill(0.0); // fluid debits per inlet [m3 d-1]
+        this.Fs = Array(config.n_inlets).fill(0); // fluid debits per inlet [m3 d-1]
-        this.Cs_in = Array.from(Array(config.n_inlets), () => new Array(NUM_SPECIES).fill(0.0)); // composition influents
+        this.Cs_in = Array.from(Array(config.n_inlets), () => new Array(NUM_SPECIES).fill(0)); // composition influents
         this.OTR = 0.0; // oxygen transfer rate [g O2 d-1]
 
         this.kla = config.kla; // if NaN, use externaly provided OTR [d-1]
 
         this.currentTime = Date.now(); // milliseconds since epoch [ms]
-        this.timeStep = 1 / (24 * 60 * 15); // time step [d]
+        this.timeStep = 1 / (24*60*15); // time step [d]
         this.speedUpFactor = 60; // speed up factor for simulation, 60 means 1 minute per simulated second
     }
 
@@ -78,7 +62,7 @@ class Reactor {
      * @returns 
      */
     _calcOTR(S_O, T = 20.0) { // caculate the OTR using basic correlation, default to temperature: 20 C
-        let S_O_sat = 14.652 - 4.1022e-1 * T + 7.9910e-3 * T * T + 7.7774e-5 * T * T * T;
+        let S_O_sat = 14.652 - 4.1022e-1 * T + 7.9910e-3 * T*T + 7.7774e-5 * T*T*T;
         return this.kla * (S_O_sat - S_O);
     }
 
@@ -102,7 +86,7 @@ class Reactor {
     updateState(newTime) { // expect update with timestamp
         const day2ms = 1000 * 60 * 60 * 24;
 
-        let n_iter = Math.floor(this.speedUpFactor * (newTime - this.currentTime) / (this.timeStep * day2ms));
+        let n_iter = Math.floor(this.speedUpFactor * (newTime-this.currentTime) / (this.timeStep*day2ms));
         if (n_iter) {
             let n = 0;
             while (n < n_iter) {
@@ -139,11 +123,11 @@ class Reactor_CSTR extends Reactor {
      * @returns {Array} - New reactor state.
      */
     tick(time_step) { // tick reactor state using forward Euler method
-        const inflow = math.multiply(math.divide([this.Fs], this.Vl), this.Cs_in)[0];
+        const inflow = math.multiply(math.divide([this.Fs], this.volume), this.Cs_in)[0];
-        const outflow = math.multiply(-1 * math.sum(this.Fs) / this.Vl, this.state);
+        const outflow = math.multiply(-1 * math.sum(this.Fs) / this.volume, this.state);
         const reaction = this.asm.compute_dC(this.state);
         const transfer = Array(NUM_SPECIES).fill(0.0);
-        transfer[OXYGEN_INDEX] = isNaN(this.kla) ? this.OTR : this._calcOTR(this.state[OXYGEN_INDEX]); // calculate OTR if kla is not NaN, otherwise use externaly calculated OTR
+        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");
@@ -166,7 +150,7 @@ class Reactor_PFR extends Reactor {
         this.n_x = config.resolution_L; // number of slices
 
         this.d_x = this.length / this.n_x;
-        this.A = this.Vl / this.length; // crosssectional area [m2]
+        this.A = this.volume / this.length; // crosssectional area [m2]
 
         this.state = Array.from(Array(this.n_x), () => config.initialState.slice())
 
@@ -198,21 +182,26 @@ class Reactor_PFR extends Reactor {
         return { topic: "Fluent", payload: { inlet: 0, F: math.sum(this.Fs), C: this.state.at(-1) }, timestamp: this.currentTime };
     }
 
+    /**
+     * Apply boundary conditions to the reactor state.
+     *  for inlet, apply generalised Danckwerts BC, if there is not flow, apply Neumann BC with no flux
+     *  for outlet, apply regular Danckwerts BC (Neumann BC with no flux)
+     * @param {Array} state - Current reactor state without enforced BCs.
+     */
     _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];
-            const BC_gradient = Array(this.n_x).fill(0.0);
+            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];
+            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)
             state[0] = state[1];
         }
         // Neumann BC (no flux)
-        state[this.n_x - 1] = state[this.n_x - 2]
+        state[this.n_x-1] = state[this.n_x-2]
     }
 
     /**
@@ -221,19 +210,19 @@ class Reactor_PFR extends Reactor {
      * @returns {Array} - New reactor state.
      */
     tick(time_step) {
-        const dispersion = math.multiply(this.D / (this.d_x * this.d_x), this.D2_op, this.state);
+        const dispersion = math.multiply(this.D / (this.d_x*this.d_x), this.D2_op, this.state);
-        const advection = math.multiply(-1 * math.sum(this.Fs) / (this.A * this.d_x), this.D_op, this.state);
+        const advection = math.multiply(-1 * math.sum(this.Fs) / (this.A*this.d_x), this.D_op, this.state);
         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.0));
+        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[OXYGEN_INDEX] = this.OTR; });
+            transfer.forEach((x) => { x[S_O_INDEX] = this.OTR; });
         } else {
-            transfer.forEach((x, i) => { x[OXYGEN_INDEX] = this._calcOTR(this.state[i][OXYGEN_INDEX]); });
+            transfer.forEach((x, i) => { x[S_O_INDEX] = this._calcOTR(this.state[i][S_O_INDEX]); });
         }
 
         const dC_total = math.multiply(math.add(dispersion, advection, reaction, transfer), time_step);
@@ -279,11 +268,11 @@ class Reactor_PFR extends Reactor {
                 throw new Error("Upwind higher order method not implemented! Use central scheme instead.");
             }
         } else {
-            const I = math.resize(math.diag(Array(this.n_x).fill(1 / (1 + central)), central), [this.n_x, this.n_x]);
+            const I = math.resize(math.diag(Array(this.n_x).fill(1 / (1+central)), central), [this.n_x, this.n_x]);
-            const A = math.resize(math.diag(Array(this.n_x).fill(-1 / (1 + central)), -1), [this.n_x, this.n_x]);
+            const A = math.resize(math.diag(Array(this.n_x).fill(-1 / (1+central)), -1), [this.n_x, this.n_x]);
             const D = math.add(I, A);
             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;
         }
     }

src/utils.js

@@ -0,0 +1,18 @@
+/**
+ * Assert that no NaN values are present in an array.
+ * @param {Array} arr 
+ * @param {string} label 
+ */
+function assertNoNaN(arr, label = "array") {
+    if (Array.isArray(arr)) {
+        for (const el of arr) {
+            assertNoNaN(el, label);
+        }
+    } else {
+        if (Number.isNaN(arr)) {
+            throw new Error(`NaN detected in ${label}!`);
+        }
+    }
+}
+
+module.exports = { assertNoNaN };