Refactor reactor class to improve NaN handling and removed magic numbers

src/reactor_class.js

@@ -1,11 +1,14 @@
-const ASM3 = require('./asm3_class.js')
+const ASM3 = require('./asm3_class.js');
-const { create, all } = require('mathjs')
+const { create, all } = require('mathjs');
 
 const config = {
     matrix: 'Array' // use Array as the matrix type
-}
+};
 
-const math = create(all, config)
+const math = create(all, config);
+
+const OXYGEN_INDEX = 0;
+const NUM_SPECIES = 13;
 
 /**
  * Assert that no NaN values are present in an array.
@@ -35,7 +38,7 @@ class Reactor {
         this.Vl = config.volume; // fluid volume reactor [m3]
 
         this.Fs = Array(config.n_inlets).fill(0.0); // fluid debits per inlet [m3 d-1]
-        this.Cs_in = Array.from(Array(config.n_inlets), () => new Array(13).fill(0.0)); // composition influents
+        this.Cs_in = Array.from(Array(config.n_inlets), () => new Array(NUM_SPECIES).fill(0.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]
@@ -136,15 +139,17 @@ class Reactor_CSTR extends Reactor {
      * @returns {Array} - New reactor state.
      */
     tick(time_step) { // tick reactor state using forward Euler method
-        const r = this.asm.compute_dC(this.state);
+        const inflow = math.multiply(math.divide([this.Fs], this.Vl), this.Cs_in)[0];
-        const dC_in = math.multiply(math.divide([this.Fs], this.Vl), this.Cs_in)[0];
+        const outflow = math.multiply(-1 * math.sum(this.Fs) / this.Vl, this.state);
-        const dC_out = math.multiply(-1 * math.sum(this.Fs) / this.Vl, this.state);
+        const reaction = this.asm.compute_dC(this.state);
-        const t_O = Array(13).fill(0.0);
+        const transfer = Array(NUM_SPECIES).fill(0.0);
-        t_O[0] = isNaN(this.kla) ? this.OTR : this._calcOTR(this.state[0]); // calculate OTR if kla is not NaN, otherwise use externaly calculated OTR
+        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
 
-        const dC_total = math.multiply(math.add(dC_in, dC_out, r, t_O), time_step);
+        const dC_total = math.multiply(math.add(inflow, outflow, reaction, transfer), time_step);
+        assertNoNaN(dC_total, "change in state");
 
         this.state = this._arrayClip2Zero(math.add(this.state, dC_total)); // clip value element-wise to avoid negative concentrations
+        assertNoNaN(this.state, "new state");
         return this.state;
     }
 }
@@ -171,9 +176,9 @@ class Reactor_PFR extends Reactor {
         this.D = 0.0; // axial dispersion [m2 d-1]
 
         this.D_op = this._makeDoperator(true, true);
-        this.D2_op = this._makeD2operator();
-
         assertNoNaN(this.D_op, "Derivative operator");
+        
+        this.D2_op = this._makeD2operator();
         assertNoNaN(this.D2_op, "Second derivative operator");
     }
 
@@ -219,25 +224,25 @@ class Reactor_PFR extends Reactor {
         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 reaction = this.state.map((state_slice) => this.asm.compute_dC(state_slice));
-        const transfer = Array.from(Array(this.n_x), () => new Array(13).fill(0.0));
+        const transfer = Array.from(Array(this.n_x), () => new Array(NUM_SPECIES).fill(0.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[0] = this.OTR; });
+            transfer.forEach((x) => { x[OXYGEN_INDEX] = this.OTR; });
         } else {
-            transfer.forEach((x, i) => { x[0] = this._calcOTR(this.state[i][0]); });
+            transfer.forEach((x, i) => { x[OXYGEN_INDEX] = this._calcOTR(this.state[i][OXYGEN_INDEX]); });
         }
 
         const dC_total = math.multiply(math.add(dispersion, advection, reaction, transfer), time_step);
-        const stateNew = math.add(this.state, dC_total);
+        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");
 
         this.state = this._arrayClip2Zero(stateNew);