From 8215c5ed9a6a736eed45a47dd8bbd6dd1772e888 Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Sat, 28 Jun 2025 19:19:38 +0200 Subject: [PATCH] Add checks for NaN values in Reactor_PFR calculations and update hydrolysis rate calculation to handle division by zero --- dependencies/asm3_class.js | 2 +- dependencies/reactor_class.js | 33 ++++++++++++++++++++++++++------- 2 files changed, 27 insertions(+), 8 deletions(-) diff --git a/dependencies/asm3_class.js b/dependencies/asm3_class.js index 6f1e5fe..f5ac1df 100644 --- a/dependencies/asm3_class.js +++ b/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; diff --git a/dependencies/reactor_class.js b/dependencies/reactor_class.js index 94844ee..8212aaf 100644 --- a/dependencies/reactor_class.js +++ b/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("Co " + math.sum(this.Fs)*this.timeStep/(this.A*this.d_x)); + console.log("Pe total " + this.length*math.sum(this.Fs)/(this.D*this.A)); + 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 A = math.resize(math.diag(Array(this.n_x).fill(-1), -1), [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 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);