From a2cfb20e2ca405671ca86b99c6b3a92b93ee6038 Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Fri, 4 Jul 2025 16:28:35 +0200 Subject: [PATCH] Refactor reactor class to improve NaN handling and add utility function for NaN assertions --- src/reactor_class.js | 71 +++++++++++++++++++------------------------- src/utils.js | 18 +++++++++++ 2 files changed, 48 insertions(+), 41 deletions(-) create mode 100644 src/utils.js diff --git a/src/reactor_class.js b/src/reactor_class.js index 2062d8f..dd5daf7 100644 --- a/src/reactor_class.js +++ b/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.Cs_in = Array.from(Array(config.n_inlets), () => new Array(NUM_SPECIES).fill(0.0)); // composition influents + 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)); // 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 outflow = math.multiply(-1 * math.sum(this.Fs) / this.Vl, this.state); + const inflow = math.multiply(math.divide([this.Fs], this.volume), this.Cs_in)[0]; + 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()) @@ -177,7 +161,7 @@ class Reactor_PFR extends Reactor { this.D_op = this._makeDoperator(true, true); assertNoNaN(this.D_op, "Derivative operator"); - + this.D2_op = this._makeD2operator(); assertNoNaN(this.D2_op, "Second derivative operator"); } @@ -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 advection = math.multiply(-1 * math.sum(this.Fs) / (this.A * this.d_x), this.D_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 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 A = math.resize(math.diag(Array(this.n_x).fill(-1 / (1 + central)), -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); + D[this.n_x-1] = Array(this.n_x).fill(0); return D; } } diff --git a/src/utils.js b/src/utils.js new file mode 100644 index 0000000..2ca1896 --- /dev/null +++ b/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 }; \ No newline at end of file