From c6b0cab06785efb050e4a978728d374321c17c0e Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Fri, 4 Jul 2025 15:14:03 +0200 Subject: [PATCH] Refactor advanced-reactor and nodeClass for improved readability and consistency --- advanced-reactor.js | 4 +- src/nodeClass.js | 1 - src/reactor_class.js | 115 +++++++++++++++++++++++-------------------- 3 files changed, 64 insertions(+), 56 deletions(-) diff --git a/advanced-reactor.js b/advanced-reactor.js index ee954cd..9ca6959 100644 --- a/advanced-reactor.js +++ b/advanced-reactor.js @@ -1,9 +1,9 @@ const nameOfNode = "advanced-reactor"; // name of the node, should match file name and node type in Node-RED const nodeClass = require('./src/nodeClass.js'); // node class -module.exports = function(RED) { +module.exports = function (RED) { // Register the node type - RED.nodes.registerType(nameOfNode, function(config) { + RED.nodes.registerType(nameOfNode, function (config) { // Initialize the Node-RED node first RED.nodes.createNode(this, config); // Then create your custom class and attach it diff --git a/src/nodeClass.js b/src/nodeClass.js index b9485be..a873284 100644 --- a/src/nodeClass.js +++ b/src/nodeClass.js @@ -20,7 +20,6 @@ class nodeClass { this._setupClass(); this._attachInputHandler(); - } /** diff --git a/src/reactor_class.js b/src/reactor_class.js index 8473fe8..d483201 100644 --- a/src/reactor_class.js +++ b/src/reactor_class.js @@ -2,12 +2,17 @@ const ASM3 = require('./asm3_class.js') const { create, all } = require('mathjs') const config = { - matrix: 'Array' // choose 'Matrix' (default) or 'Array' + matrix: 'Array' // use Array as the matrix type } const math = create(all, config) -function assertNoNaN(arr, label="array") { +/** + * Assert that no NaN values are present in an array. + * @param {Array} arr + * @param {string} label + */ +function assertNoNaN(arr, label = "array") { if (math.isNaN(arr)) { throw new Error("NaN detected in ${label}!"); } @@ -18,7 +23,7 @@ class Reactor { * Reactor base class. * @param {object} config - Configuration object containing reactor parameters. */ - constructor(config){ + constructor(config) { this.asm = new ASM3(); this.Vl = config.volume; // fluid volume reactor [m3] @@ -30,7 +35,7 @@ class Reactor { 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 } @@ -63,8 +68,8 @@ 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 - let S_O_sat = 14.652 - 4.1022e-1*T + 7.9910e-3*T*T + 7.7774e-5*T*T*T; + _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); } @@ -86,9 +91,9 @@ class Reactor { * @param {number} newTime - New time to update reactor state to, in milliseconds since epoch. */ updateState(newTime) { // expect update with timestamp - const day2ms = 1000 * 60 * 60 * 24; + 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) { @@ -116,7 +121,7 @@ class Reactor_CSTR extends Reactor { * @returns {object} Effluent data object (msg), defaults to inlet 0. */ get getEffluent() { // getter for Effluent, defaults to inlet 0 - return {topic: "Fluent", payload: {inlet: 0, F: math.sum(this.Fs), C:this.state}, timestamp: this.currentTime}; + return { topic: "Fluent", payload: { inlet: 0, F: math.sum(this.Fs), C: this.state }, timestamp: this.currentTime }; } /** @@ -127,13 +132,13 @@ class Reactor_CSTR extends Reactor { tick(time_step) { // tick reactor state using forward Euler method const r = this.asm.compute_dC(this.state); 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 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 const dC_total = math.multiply(math.add(dC_in, dC_out, r, t_O), time_step); - this.state = this_.arrayClip2Zero(math.add(this.state, dC_total)); // clip value element-wise to avoid negative concentrations + this.state = this._arrayClip2Zero(math.add(this.state, dC_total)); // clip value element-wise to avoid negative concentrations return this.state; } } @@ -153,7 +158,7 @@ class Reactor_PFR extends Reactor { this.A = this.Vl / this.length; // crosssectional area [m2] this.state = Array.from(Array(this.n_x), () => config.initialState.slice()) - + // console.log("Initial State: ") // console.log(this.state) @@ -162,7 +167,7 @@ class Reactor_PFR extends Reactor { this.D_op = this._makeDoperator(true, true); this.D2_op = this._makeD2operator(); } - + /** * Setter for axial dispersion. * @param {object} input - Input object (msg) containing payload with dispersion value [m2 d-1]. @@ -176,7 +181,25 @@ class Reactor_PFR extends Reactor { * @returns {object} Effluent data object (msg), defaults to inlet 0. */ get getEffluent() { - return {topic: "Fluent", payload: {inlet: 0, F: math.sum(this.Fs), C:this.state.at(-1)}, timestamp: this.currentTime}; + return { topic: "Fluent", payload: { inlet: 0, F: math.sum(this.Fs), C: this.state.at(-1) }, timestamp: this.currentTime }; + } + + _applyBoundaryConditions(newState) { + // 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); + 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], stateNew)[0]; + newState[0] = math.add(BC_C_in, BC_dispersion).map(val => val < 0 ? 0 : val); + } else { // Neumann BC (no flux) + newState[0] = newState[1]; + } + // Neumann BC (no flux) + newState[this.n_x - 1] = newState[this.n_x - 2] + return newState } /** @@ -185,15 +208,15 @@ 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(13).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; }); } else { @@ -201,30 +224,16 @@ class Reactor_PFR extends Reactor { } const dC_total = math.multiply(math.add(dispersion, advection, reaction, transfer), time_step); - const new_state = math.add(this.state, dC_total); + let stateNew = math.add(this.state, dC_total); - assertNoNaN(new_state, "new state"); + assertNoNaN(stateNew, "new 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); - 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], new_state)[0]; - new_state[0] = math.add(BC_C_in, BC_dispersion).map(val => val < 0 ? 0 : val); - - } 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] + stateNew = this._applyBoundaryConditions(stateNew); - assertNoNaN(new_state, "new state post BC"); + assertNoNaN(stateNew, "new state post BC"); - this.state = this._arrayClip2Zero(new_state); - return new_state; + this.state = this._arrayClip2Zero(stateNew); + return stateNew; } /** @@ -233,35 +242,35 @@ 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]); - const A = math.resize(math.diag(Array(this.n_x).fill(-2/3), -1), [this.n_x, this.n_x]); - const B = math.resize(math.diag(Array(this.n_x).fill(2/3), 1), [this.n_x, this.n_x]); - const C = math.resize(math.diag(Array(this.n_x).fill(-1/12), 2), [this.n_x, this.n_x]); + const I = math.resize(math.diag(Array(this.n_x).fill(1 / 12), -2), [this.n_x, this.n_x]); + const A = math.resize(math.diag(Array(this.n_x).fill(-2 / 3), -1), [this.n_x, this.n_x]); + const B = math.resize(math.diag(Array(this.n_x).fill(2 / 3), 1), [this.n_x, this.n_x]); + const C = math.resize(math.diag(Array(this.n_x).fill(-1 / 12), 2), [this.n_x, this.n_x]); const D = math.add(I, A, B, C); const NearBoundary = Array(this.n_x).fill(0.0); - NearBoundary[0] = -1/4; - NearBoundary[1] = -5/6; - NearBoundary[2] = 3/2; - NearBoundary[3] = -1/2; - NearBoundary[4] = 1/12; + NearBoundary[0] = -1 / 4; + NearBoundary[1] = -5 / 6; + NearBoundary[2] = 3 / 2; + NearBoundary[3] = -1 / 2; + NearBoundary[4] = 1 / 12; D[1] = NearBoundary; NearBoundary.reverse(); - D[this.n_x-2] = math.multiply(-1, NearBoundary); + D[this.n_x - 2] = math.multiply(-1, NearBoundary); 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; } else { 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; } } @@ -276,7 +285,7 @@ class Reactor_PFR extends Reactor { const B = math.resize(math.diag(Array(this.n_x).fill(1), -1), [this.n_x, this.n_x]); const D2 = math.add(I, A, B); D2[0] = Array(this.n_x).fill(0); // set by BCs elsewhere - D2[this.n_x-1] = Array(this.n_x).fill(0); + D2[this.n_x - 1] = Array(this.n_x).fill(0); return D2; } }