Add checks for NaN values in Reactor_PFR calculations and update hydrolysis rate calculation to handle division by zero

dependencies/asm3_class.js

@@ -99,7 +99,7 @@ class ASM3 {
     const { k_H, K_X, k_STO, nu_NO, K_O, K_NO, K_S, K_STO, mu_H_max, K_NH, K_HCO, b_H_O, b_H_NO, b_STO_O, b_STO_NO, mu_A_max, K_A_NH, K_A_O, K_A_HCO, b_A_O, b_A_NO } = this.kin_params;
     
     // Hydrolysis
-    rates[0] = k_H * this._monod(X_S / X_H, K_X) * X_H;
+    rates[0] = X_H == 0 ? 0 : k_H * this._monod(X_S / X_H, K_X) * X_H;
 
     // Heterotrophs
     rates[1] = k_STO * this._monod(S_O, K_O) * this._monod(S_S, K_S) * X_H;

dependencies/reactor_class.js

@@ -103,7 +103,7 @@ class Reactor_PFR {
         this.timeStep = 1/(24*60*15); // time step [d]
         this.speedUpFactor = 60;
 
-        this.D_op = this.makeDoperator(false);
+        this.D_op = this.makeDoperator(true);
         this.D2_op = this.makeD2operator();
     }
 
@@ -111,8 +111,10 @@ class Reactor_PFR {
         let index_in = input.payload.inlet;
         this.Fs[index_in] = input.payload.F;
         this.Cs_in[index_in] = input.payload.C;
-        // console.log("Pe " + this.d_x*math.sum(this.Fs)/(this.D*this.A));
+        console.log("Pe total " + this.length*math.sum(this.Fs)/(this.D*this.A));
-        // console.log("Co " + math.sum(this.Fs)*this.timeStep/(this.A*this.d_x));
+        console.log("Pe local " + this.d_x*math.sum(this.Fs)/(this.D*this.A));
+        console.log("Co ad " + math.sum(this.Fs)*this.timeStep/(this.A*this.d_x));
+        console.log("Co D " + this.D*this.timeStep/(this.d_x*this.d_x));
     }
 
     set setOTR(input) { // setter for OTR (WIP) [g O2 d-1]
@@ -134,7 +136,6 @@ class Reactor_PFR {
 
     // expect update with timestamp
     updateState(newTime) {
-
         const day2ms = 1000 * 60 * 60 * 24; 
 
         let n_iter = Math.floor(this.speedUpFactor*(newTime - this.currentTime) / (this.timeStep * day2ms));
@@ -153,6 +154,16 @@ class Reactor_PFR {
         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));
+
+        if (dispersion.some(row => row.some(Number.isNaN))) {
+            throw new Error("NaN detected in dispersion!");
+        }
+        if (advection.some(row => row.some(Number.isNaN))) {
+            throw new Error("NaN detected in advection!");
+        }
+        if (reaction.some(row => row.some(Number.isNaN))) {
+            throw new Error("NaN detected in reaction!");
+        }
         
         if (isNaN(this.kla)) { // calculate OTR if kla is not NaN, otherwise use externally calculated OTR
             transfer.forEach((x) => { x[0] = this.OTR; });
@@ -162,6 +173,9 @@ class Reactor_PFR {
 
         const dC_total = math.multiply(math.add(dispersion, advection, reaction, transfer), time_step);
         const new_state = math.add(this.state, dC_total);
+        if (new_state.some(row => row.some(Number.isNaN))) {
+            throw new Error("NaN detected in new_state after dC_total update!");
+        }
 
         // apply boundary conditions
         if (math.sum(this.Fs) > 0) { // Danckwerts BC
@@ -170,21 +184,26 @@ class Reactor_PFR {
             BC_gradient[0] = -1;
             BC_gradient[1] = 1;
             const BC_dispersion = math.multiply(this.D * this.A / (math.sum(this.Fs)*this.d_x), [BC_gradient], new_state)[0];
+            console.log(math.add(BC_C_in, BC_dispersion));
             new_state[0] = math.add(BC_C_in, BC_dispersion).map(val => val < 0 ? 0 : val);
-            console.log(new_state[0])
+            
         } else { // Neumann BC (no flux)
             new_state[0] = new_state[1];
         }
         // Neumann BC (no flux)
         new_state[this.n_x-1] = new_state[this.n_x-2]
 
+        if (new_state.some(row => row.some(Number.isNaN))) {
+            throw new Error("NaN detected in new_state after enforcing boundary conditions!");
+        }
+
         this.state = new_state.map(row => row.map(val => val < 0 ? 0 : val)); // apply the new state
         return new_state;
     }
 
     makeDoperator(central=false) { // create the upwind scheme gradient operator
-        const I = math.resize(math.diag(Array(this.n_x).fill(1), 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), [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);