From 3f5b0eea32e055a3f3ff518ee92474642ea4c987 Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Fri, 4 Jul 2025 14:58:38 +0200 Subject: [PATCH] Enhance reactor class with NaN checks and refactor methods for clarity --- src/reactor_class.js | 66 +++++++++++++++++++++++++------------------- 1 file changed, 37 insertions(+), 29 deletions(-) diff --git a/src/reactor_class.js b/src/reactor_class.js index 245f194..8473fe8 100644 --- a/src/reactor_class.js +++ b/src/reactor_class.js @@ -1,4 +1,4 @@ -const ASM3 = require('./asm3_class') +const ASM3 = require('./asm3_class.js') const { create, all } = require('mathjs') const config = { @@ -7,6 +7,12 @@ const config = { const math = create(all, config) +function assertNoNaN(arr, label="array") { + if (math.isNaN(arr)) { + throw new Error("NaN detected in ${label}!"); + } +} + class Reactor { /** * Reactor base class. @@ -57,11 +63,24 @@ class Reactor { * @param {number} T - Temperature in Celsius, default to 20 C. * @returns */ - calcOTR(S_O, T=20.0) { // caculate the OTR using basic correlation, default to temperature: 20 C + _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; return this.kla * (S_O_sat - S_O); } + /** + * Clip values in an array to zero. + * @param {Array} arr - Array of values to clip. + * @returns {Array} - New array with values clipped to zero. + */ + _arrayClip2Zero(arr) { + if (Array.isArray(arr)) { + return arr.map(clipToZero); + } else { + return arr < 0 ? 0 : arr; + } + } + /** * Update the reactor state based on the new time. * @param {number} newTime - New time to update reactor state to, in milliseconds since epoch. @@ -110,12 +129,11 @@ class Reactor_CSTR extends Reactor { const dC_in = math.multiply(math.divide([this.Fs], this.Vl), this.Cs_in)[0]; const dC_out = math.multiply(-1*math.sum(this.Fs)/this.Vl, this.state); const t_O = Array(13).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 + 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 const dC_total = math.multiply(math.add(dC_in, dC_out, r, t_O), time_step); - // clip value element-wise to each subarray to avoid negative concentrations - this.state = math.add(this.state, dC_total).map(val => val < 0 ? 0 : val); + this.state = this_.arrayClip2Zero(math.add(this.state, dC_total)); // clip value element-wise to avoid negative concentrations return this.state; } } @@ -141,8 +159,8 @@ 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(); + this.D_op = this._makeDoperator(true, true); + this.D2_op = this._makeD2operator(); } /** @@ -172,27 +190,20 @@ class Reactor_PFR extends Reactor { 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!"); - } + 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; }); } else { - transfer.forEach((x, i) => { x[0] = this.calcOTR(this.state[i][0]); }); + transfer.forEach((x, i) => { x[0] = this._calcOTR(this.state[i][0]); }); } 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!"); - } + + assertNoNaN(new_state, "new state"); // apply boundary conditions if (math.sum(this.Fs) > 0) { // Danckwerts BC @@ -210,11 +221,9 @@ class Reactor_PFR extends Reactor { // 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!"); - } + assertNoNaN(new_state, "new state post BC"); - this.state = new_state.map(row => row.map(val => val < 0 ? 0 : val)); // apply the new state + this.state = this._arrayClip2Zero(new_state); return new_state; } @@ -224,7 +233,7 @@ class Reactor_PFR extends Reactor { * @param {boolean} higher_order - Use higher order scheme if true, otherwise use first order scheme. * @returns {Array} - First derivative operator matrix. */ - makeDoperator(central=false, higher_order=false) { // create gradient operator + _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]); @@ -261,7 +270,7 @@ class Reactor_PFR extends Reactor { * Create central finite difference second derivative operator. * @returns {Array} - Second derivative operator matrix. */ - makeD2operator() { // create the central second derivative operator + _makeD2operator() { // create the central second derivative operator const I = math.diag(Array(this.n_x).fill(-2), 0); const A = math.resize(math.diag(Array(this.n_x).fill(1), 1), [this.n_x, this.n_x]); const B = math.resize(math.diag(Array(this.n_x).fill(1), -1), [this.n_x, this.n_x]); @@ -272,6 +281,7 @@ class Reactor_PFR extends Reactor { } } +module.exports = { Reactor_CSTR, Reactor_PFR }; // DEBUG // state: S_O, S_I, S_S, S_NH, S_N2, S_NO, S_HCO, X_I, X_S, X_H, X_STO, X_A, X_TS @@ -284,6 +294,4 @@ class Reactor_PFR extends Reactor { // while (N < 5000) { // console.log(Reactor.tick(0.001)); // N += 1; -// } - -module.exports = { Reactor_CSTR, Reactor_PFR }; \ No newline at end of file +// } \ No newline at end of file