From d0db1b416c36bf9b73238770dd269a6d15ac794a Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Fri, 4 Jul 2025 10:01:46 +0200 Subject: [PATCH 01/18] Remove debug messages --- dependencies/reactor_class.js | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/dependencies/reactor_class.js b/dependencies/reactor_class.js index 76a6340..c1e8989 100644 --- a/dependencies/reactor_class.js +++ b/dependencies/reactor_class.js @@ -111,10 +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 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)); + // 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] @@ -185,7 +185,6 @@ class Reactor_PFR { 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]; - 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); } else { // Neumann BC (no flux) From 25cd728b68e04a3ceca70d0cf4288bffd8ca3d00 Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Fri, 4 Jul 2025 10:44:54 +0200 Subject: [PATCH 02/18] Refactor reactor node registration --- advanced-reactor.js | 105 ++---------------------- {dependencies => src}/asm3_class.js | 0 src/nodeClass.js | 109 +++++++++++++++++++++++++ {dependencies => src}/reactor_class.js | 0 4 files changed, 118 insertions(+), 96 deletions(-) rename {dependencies => src}/asm3_class.js (100%) create mode 100644 src/nodeClass.js rename {dependencies => src}/reactor_class.js (100%) diff --git a/advanced-reactor.js b/advanced-reactor.js index e2eca51..ee954cd 100644 --- a/advanced-reactor.js +++ b/advanced-reactor.js @@ -1,99 +1,12 @@ +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) { - function reactor(config) { + // Register the node type + RED.nodes.registerType(nameOfNode, function(config) { + // Initialize the Node-RED node first RED.nodes.createNode(this, config); - var node = this; - - let name = config.name; - - const { Reactor_CSTR, Reactor_PFR } = require('./dependencies/reactor_class'); - - let new_reactor; - - switch (config.reactor_type) { - case "CSTR": - new_reactor = new Reactor_CSTR( - parseFloat(config.volume), - parseInt(config.n_inlets), - parseFloat(config.kla), - [ - parseFloat(config.S_O_init), - parseFloat(config.S_I_init), - parseFloat(config.S_S_init), - parseFloat(config.S_NH_init), - parseFloat(config.S_N2_init), - parseFloat(config.S_NO_init), - parseFloat(config.S_HCO_init), - parseFloat(config.X_I_init), - parseFloat(config.X_S_init), - parseFloat(config.X_H_init), - parseFloat(config.X_STO_init), - parseFloat(config.X_A_init), - parseFloat(config.X_TS_init) - ] - ); - break; - case "PFR": - new_reactor = new Reactor_PFR( - parseFloat(config.volume), - parseFloat(config.length), - parseInt(config.resolution_L), - parseInt(config.n_inlets), - parseFloat(config.kla), - [ - parseFloat(config.S_O_init), - parseFloat(config.S_I_init), - parseFloat(config.S_S_init), - parseFloat(config.S_NH_init), - parseFloat(config.S_N2_init), - parseFloat(config.S_NO_init), - parseFloat(config.S_HCO_init), - parseFloat(config.X_I_init), - parseFloat(config.X_S_init), - parseFloat(config.X_H_init), - parseFloat(config.X_STO_init), - parseFloat(config.X_A_init), - parseFloat(config.X_TS_init) - ] - ); - break; - default: - console.warn("Unknown reactor type: " + config.reactor_type); - } - - const reactor = new_reactor; // protect from reassignment - - node.on('input', function(msg, send, done) { - let toggleUpdate = false; - - switch (msg.topic) { - case "clock": - toggleUpdate = true; - break; - case "Fluent": - reactor.setInfluent = msg; - if (msg.payload.inlet == 0) { - toggleUpdate = true; - } - break; - case "OTR": - reactor.setOTR = msg; - break; - case "Dispersion": - reactor.setDispersion = msg; - break; - default: - console.log("Unknown topic: " + msg.topic); - } - - if (toggleUpdate) { - reactor.updateState(msg.timestamp); - send(reactor.getEffluent); - } - - if (done) { - done(); - } - }); - } - RED.nodes.registerType("advanced-reactor", reactor); + // Then create your custom class and attach it + this.nodeClass = new nodeClass(config, RED, this, nameOfNode); + }); }; diff --git a/dependencies/asm3_class.js b/src/asm3_class.js similarity index 100% rename from dependencies/asm3_class.js rename to src/asm3_class.js diff --git a/src/nodeClass.js b/src/nodeClass.js new file mode 100644 index 0000000..9281a45 --- /dev/null +++ b/src/nodeClass.js @@ -0,0 +1,109 @@ +const { Reactor_CSTR, Reactor_PFR } = require('./reactor_class.js'); + + +class nodeClass { + /** + * Create a ReactorNode. + * @param {object} uiConfig - Node-RED node configuration. + * @param {object} RED - Node-RED runtime API. + * @param {object} nodeInstance - The Node-RED node instance. + * @param {string} nameOfNode - The name of the node, used for + */ + constructor(uiConfig, RED, nodeInstance, nameOfNode) { + // Preserve RED reference for HTTP endpoints if needed + this.node = nodeInstance; + this.RED = RED; + this.name = nameOfNode; + + let new_reactor; + + switch (uiConfig.reactor_type) { + case "CSTR": + new_reactor = new Reactor_CSTR( + parseFloat(uiConfig.volume), + parseInt(uiConfig.n_inlets), + parseFloat(uiConfig.kla), + [ + parseFloat(uiConfig.S_O_init), + parseFloat(uiConfig.S_I_init), + parseFloat(uiConfig.S_S_init), + parseFloat(uiConfig.S_NH_init), + parseFloat(uiConfig.S_N2_init), + parseFloat(uiConfig.S_NO_init), + parseFloat(uiConfig.S_HCO_init), + parseFloat(uiConfig.X_I_init), + parseFloat(uiConfig.X_S_init), + parseFloat(uiConfig.X_H_init), + parseFloat(uiConfig.X_STO_init), + parseFloat(uiConfig.X_A_init), + parseFloat(uiConfig.X_TS_init) + ] + ); + break; + case "PFR": + new_reactor = new Reactor_PFR( + parseFloat(uiConfig.volume), + parseFloat(uiConfig.length), + parseInt(uiConfig.resolution_L), + parseInt(uiConfig.n_inlets), + parseFloat(uiConfig.kla), + [ + parseFloat(uiConfig.S_O_init), + parseFloat(uiConfig.S_I_init), + parseFloat(uiConfig.S_S_init), + parseFloat(uiConfig.S_NH_init), + parseFloat(uiConfig.S_N2_init), + parseFloat(uiConfig.S_NO_init), + parseFloat(uiConfig.S_HCO_init), + parseFloat(uiConfig.X_I_init), + parseFloat(uiConfig.X_S_init), + parseFloat(uiConfig.X_H_init), + parseFloat(uiConfig.X_STO_init), + parseFloat(uiConfig.X_A_init), + parseFloat(uiConfig.X_TS_init) + ] + ); + break; + default: + console.warn("Unknown reactor type: " + uiConfig.reactor_type); + } + + const reactor = new_reactor; // protect from reassignment + + this.node.on('input', function(msg, send, done) { + let toggleUpdate = false; + + switch (msg.topic) { + case "clock": + toggleUpdate = true; + break; + case "Fluent": + reactor.setInfluent = msg; + if (msg.payload.inlet == 0) { + toggleUpdate = true; + } + break; + case "OTR": + reactor.setOTR = msg; + break; + case "Dispersion": + reactor.setDispersion = msg; + break; + default: + console.log("Unknown topic: " + msg.topic); + } + + if (toggleUpdate) { + reactor.updateState(msg.timestamp); + send(reactor.getEffluent); + } + + if (done) { + done(); + } + }); + } + +} + +module.exports = nodeClass; \ No newline at end of file diff --git a/dependencies/reactor_class.js b/src/reactor_class.js similarity index 100% rename from dependencies/reactor_class.js rename to src/reactor_class.js From 1cda956d832cecc2a079cdc14941601c0f15e170 Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Fri, 4 Jul 2025 11:42:34 +0200 Subject: [PATCH 03/18] Refactor Reactor class structure and include inheritance for CSTR and PFR --- src/reactor_class.js | 66 +++++++++++++++++--------------------------- 1 file changed, 25 insertions(+), 41 deletions(-) diff --git a/src/reactor_class.js b/src/reactor_class.js index c1e8989..a81ff7f 100644 --- a/src/reactor_class.js +++ b/src/reactor_class.js @@ -7,43 +7,56 @@ const config = { const math = create(all, config) -class Reactor_CSTR { - - constructor(volume, n_inlets, kla, initial_state) { - this.state = initial_state; +class Reactor { + + constructor(volume, n_inlets, kla){ this.asm = new ASM3(); this.Vl = volume; // fluid volume reactor [m3] + this.Fs = Array(n_inlets).fill(0.0); // fluid debits per inlet [m3 d-1] this.Cs_in = Array.from(Array(n_inlets), () => new Array(13).fill(0.0)); // composition influents this.OTR = 0.0; // oxygen transfer rate [g O2 d-1] - + this.kla = kla; // if NaN, use external OTR [d-1] this.currentTime = Date.now(); // milliseconds since epoch [ms] this.timeStep = 1/(24*60*15); // time step [d] - this.speedUpFactor = 1; + this.speedUpFactor = 60; } set setInfluent(input) { // setter for C_in (WIP) let index_in = input.payload.inlet; this.Fs[index_in] = input.payload.F; this.Cs_in[index_in] = input.payload.C; + // 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] this.OTR = input.payload; } - 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}; - } - 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); } +} + +class Reactor_CSTR extends Reactor { + + constructor(volume, n_inlets, kla, initial_state) { + super(volume, n_inlets, kla); + this.state = initial_state; + } + + 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}; + } + // expect update with timestamp updateState(newTime) { @@ -75,12 +88,11 @@ class Reactor_CSTR { } } -class Reactor_PFR { +class Reactor_PFR extends Reactor { constructor(volume, length, resolution_L, n_inlets, kla, initial_state) { - this.asm = new ASM3(); + super(volume, n_inlets, kla); - this.Vl = volume; // fluid volume reactor [m3] this.length = length; // reactor length [m] this.n_x = resolution_L; // number of slices this.d_x = length / resolution_L; @@ -92,35 +104,12 @@ class Reactor_PFR { // console.log("Initial State: ") // console.log(this.state) - this.Fs = Array(n_inlets).fill(0.0); // fluid debits per inlet [m3 d-1] - this.Cs_in = Array.from(Array(n_inlets), () => new Array(13).fill(0.0)); // composition influents - this.OTR = 0.0; // oxygen transfer rate [g O2 d-1] this.D = 0.0; // axial dispersion [m2 d-1] - - this.kla = kla; // if NaN, use external OTR [d-1] - - this.currentTime = Date.now(); // milliseconds since epoch [ms] - this.timeStep = 1/(24*60*15); // time step [d] - this.speedUpFactor = 60; this.D_op = this.makeDoperator(true, true); this.D2_op = this.makeD2operator(); } - set setInfluent(input) { // setter for C_in (WIP) - let index_in = input.payload.inlet; - this.Fs[index_in] = input.payload.F; - this.Cs_in[index_in] = input.payload.C; - // 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] - this.OTR = input.payload; - } - set setDispersion(input) { // setter for Axial dispersion [m2 d-1] this.D = input.payload; } @@ -129,11 +118,6 @@ class Reactor_PFR { return {topic: "Fluent", payload: {inlet: 0, F: math.sum(this.Fs), C:this.state.at(-1)}, timestamp: this.currentTime}; } - 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); - } - // expect update with timestamp updateState(newTime) { const day2ms = 1000 * 60 * 60 * 24; From fee6881f1b9847ca68602022871ccfcf2141ea39 Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Fri, 4 Jul 2025 12:06:58 +0200 Subject: [PATCH 04/18] Refactor nodeClass for to mostly allign with the standard EVOLV structure --- src/nodeClass.js | 78 ++++++++++++++++++++++++++---------------------- 1 file changed, 43 insertions(+), 35 deletions(-) diff --git a/src/nodeClass.js b/src/nodeClass.js index 9281a45..ad11f69 100644 --- a/src/nodeClass.js +++ b/src/nodeClass.js @@ -15,6 +15,48 @@ class nodeClass { this.RED = RED; this.name = nameOfNode; + this._setupClass(uiConfig, this.node); + + this._attachInputHandler(); + + } + + _attachInputHandler() { // Handle input messages + this.node.on('input', function(msg, send, done) { + let toggleUpdate = false; + + switch (msg.topic) { + case "clock": + toggleUpdate = true; + break; + case "Fluent": + this.reactor.setInfluent = msg; + if (msg.payload.inlet == 0) { + toggleUpdate = true; + } + break; + case "OTR": + this.reactor.setOTR = msg; + break; + case "Dispersion": + this.reactor.setDispersion = msg; + break; + default: + console.log("Unknown topic: " + msg.topic); + } + + if (toggleUpdate) { + this.reactor.updateState(msg.timestamp); + send(this.reactor.getEffluent); + } + + if (done) { + done(); + } + }); + } + + _setupClass(uiConfig, node) { let new_reactor; switch (uiConfig.reactor_type) { @@ -68,42 +110,8 @@ class nodeClass { console.warn("Unknown reactor type: " + uiConfig.reactor_type); } - const reactor = new_reactor; // protect from reassignment - - this.node.on('input', function(msg, send, done) { - let toggleUpdate = false; - - switch (msg.topic) { - case "clock": - toggleUpdate = true; - break; - case "Fluent": - reactor.setInfluent = msg; - if (msg.payload.inlet == 0) { - toggleUpdate = true; - } - break; - case "OTR": - reactor.setOTR = msg; - break; - case "Dispersion": - reactor.setDispersion = msg; - break; - default: - console.log("Unknown topic: " + msg.topic); - } - - if (toggleUpdate) { - reactor.updateState(msg.timestamp); - send(reactor.getEffluent); - } - - if (done) { - done(); - } - }); + node.reactor = new_reactor; // protect from reassignment } - } module.exports = nodeClass; \ No newline at end of file From c23818c10805bf60cdd7386c68e2ed7bc0e6c0ac Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Fri, 4 Jul 2025 12:16:08 +0200 Subject: [PATCH 05/18] Remove unnecessary node parameter _setupClass --- src/nodeClass.js | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/nodeClass.js b/src/nodeClass.js index ad11f69..4cbb4bf 100644 --- a/src/nodeClass.js +++ b/src/nodeClass.js @@ -15,7 +15,7 @@ class nodeClass { this.RED = RED; this.name = nameOfNode; - this._setupClass(uiConfig, this.node); + this._setupClass(uiConfig); this._attachInputHandler(); @@ -56,7 +56,7 @@ class nodeClass { }); } - _setupClass(uiConfig, node) { + _setupClass(uiConfig) { let new_reactor; switch (uiConfig.reactor_type) { @@ -110,7 +110,7 @@ class nodeClass { console.warn("Unknown reactor type: " + uiConfig.reactor_type); } - node.reactor = new_reactor; // protect from reassignment + this.reactor = new_reactor; // protect from reassignment } } From 530dac5c771b7dd6110c64965f3cc96cf52cd1aa Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Fri, 4 Jul 2025 12:51:37 +0200 Subject: [PATCH 06/18] Refactor nodeClass to streamline configuration loading and reactor setup --- src/nodeClass.js | 94 ++++++++++++++++++++++++------------------------ 1 file changed, 47 insertions(+), 47 deletions(-) diff --git a/src/nodeClass.js b/src/nodeClass.js index 4cbb4bf..917f1c0 100644 --- a/src/nodeClass.js +++ b/src/nodeClass.js @@ -3,11 +3,11 @@ const { Reactor_CSTR, Reactor_PFR } = require('./reactor_class.js'); class nodeClass { /** - * Create a ReactorNode. - * @param {object} uiConfig - Node-RED node configuration. - * @param {object} RED - Node-RED runtime API. - * @param {object} nodeInstance - The Node-RED node instance. - * @param {string} nameOfNode - The name of the node, used for + * Create ReactorNode. + * @param {object} uiConfig - Node-RED node configuration + * @param {object} RED - Node-RED runtime API + * @param {object} nodeInstance - Node-RED node instance + * @param {string} nameOfNode - Name of the node */ constructor(uiConfig, RED, nodeInstance, nameOfNode) { // Preserve RED reference for HTTP endpoints if needed @@ -15,14 +15,16 @@ class nodeClass { this.RED = RED; this.name = nameOfNode; - this._setupClass(uiConfig); + this._loadConfig(uiConfig) + + this._setupClass(); this._attachInputHandler(); } _attachInputHandler() { // Handle input messages - this.node.on('input', function(msg, send, done) { + this.node.on('input', (msg, send, done) => { let toggleUpdate = false; switch (msg.topic) { @@ -56,54 +58,52 @@ class nodeClass { }); } - _setupClass(uiConfig) { + _loadConfig(uiConfig) { + this.config = { + reactor_type: uiConfig.reactor_type, + volume: parseFloat(uiConfig.volume), + length: parseFloat(uiConfig.length), + resolution_L: parseInt(uiConfig.resolution_L), + n_inlets: parseInt(uiConfig.n_inlets), + kla: parseFloat(uiConfig.kla), + initialState: [ + parseFloat(uiConfig.S_O_init), + parseFloat(uiConfig.S_I_init), + parseFloat(uiConfig.S_S_init), + parseFloat(uiConfig.S_NH_init), + parseFloat(uiConfig.S_N2_init), + parseFloat(uiConfig.S_NO_init), + parseFloat(uiConfig.S_HCO_init), + parseFloat(uiConfig.X_I_init), + parseFloat(uiConfig.X_S_init), + parseFloat(uiConfig.X_H_init), + parseFloat(uiConfig.X_STO_init), + parseFloat(uiConfig.X_A_init), + parseFloat(uiConfig.X_TS_init) + ] + } + } + + _setupClass() { let new_reactor; - switch (uiConfig.reactor_type) { + switch (this.config.reactor_type) { case "CSTR": new_reactor = new Reactor_CSTR( - parseFloat(uiConfig.volume), - parseInt(uiConfig.n_inlets), - parseFloat(uiConfig.kla), - [ - parseFloat(uiConfig.S_O_init), - parseFloat(uiConfig.S_I_init), - parseFloat(uiConfig.S_S_init), - parseFloat(uiConfig.S_NH_init), - parseFloat(uiConfig.S_N2_init), - parseFloat(uiConfig.S_NO_init), - parseFloat(uiConfig.S_HCO_init), - parseFloat(uiConfig.X_I_init), - parseFloat(uiConfig.X_S_init), - parseFloat(uiConfig.X_H_init), - parseFloat(uiConfig.X_STO_init), - parseFloat(uiConfig.X_A_init), - parseFloat(uiConfig.X_TS_init) - ] + this.config.volume, + this.config.n_inlets, + this.config.kla, + this.config.initialState ); break; case "PFR": new_reactor = new Reactor_PFR( - parseFloat(uiConfig.volume), - parseFloat(uiConfig.length), - parseInt(uiConfig.resolution_L), - parseInt(uiConfig.n_inlets), - parseFloat(uiConfig.kla), - [ - parseFloat(uiConfig.S_O_init), - parseFloat(uiConfig.S_I_init), - parseFloat(uiConfig.S_S_init), - parseFloat(uiConfig.S_NH_init), - parseFloat(uiConfig.S_N2_init), - parseFloat(uiConfig.S_NO_init), - parseFloat(uiConfig.S_HCO_init), - parseFloat(uiConfig.X_I_init), - parseFloat(uiConfig.X_S_init), - parseFloat(uiConfig.X_H_init), - parseFloat(uiConfig.X_STO_init), - parseFloat(uiConfig.X_A_init), - parseFloat(uiConfig.X_TS_init) - ] + this.config.volume, + this.config.length, + this.config.resolution_L, + this.config.n_inlets, + this.config.kla, + this.config.initialState ); break; default: From 348307d99919688cebae08e457d6a099e65ba539 Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Fri, 4 Jul 2025 13:09:20 +0200 Subject: [PATCH 07/18] Add documentation --- src/nodeClass.js | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/src/nodeClass.js b/src/nodeClass.js index 917f1c0..62a6226 100644 --- a/src/nodeClass.js +++ b/src/nodeClass.js @@ -3,7 +3,7 @@ const { Reactor_CSTR, Reactor_PFR } = require('./reactor_class.js'); class nodeClass { /** - * Create ReactorNode. + * Construct ReactorNode. * @param {object} uiConfig - Node-RED node configuration * @param {object} RED - Node-RED runtime API * @param {object} nodeInstance - Node-RED node instance @@ -23,7 +23,10 @@ class nodeClass { } - _attachInputHandler() { // Handle input messages + /** + * Handle node-red input messages + */ + _attachInputHandler() { this.node.on('input', (msg, send, done) => { let toggleUpdate = false; @@ -58,6 +61,10 @@ class nodeClass { }); } + /** + * Parse node configuration + * @param {object} uiConfig Config set in UI in node-red + */ _loadConfig(uiConfig) { this.config = { reactor_type: uiConfig.reactor_type, @@ -84,6 +91,9 @@ class nodeClass { } } + /** + * Setup reactor class based on config + */ _setupClass() { let new_reactor; From c239b71ad87b9ee9243a858ab4483e0a803cb19d Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Fri, 4 Jul 2025 13:16:49 +0200 Subject: [PATCH 08/18] Refactor reactor constructors to accept a config object for improved clarity and maintainability --- src/nodeClass.js | 16 ++-------------- src/reactor_class.js | 30 +++++++++++++++--------------- 2 files changed, 17 insertions(+), 29 deletions(-) diff --git a/src/nodeClass.js b/src/nodeClass.js index 62a6226..15a959b 100644 --- a/src/nodeClass.js +++ b/src/nodeClass.js @@ -99,22 +99,10 @@ class nodeClass { switch (this.config.reactor_type) { case "CSTR": - new_reactor = new Reactor_CSTR( - this.config.volume, - this.config.n_inlets, - this.config.kla, - this.config.initialState - ); + new_reactor = new Reactor_CSTR(this.config); break; case "PFR": - new_reactor = new Reactor_PFR( - this.config.volume, - this.config.length, - this.config.resolution_L, - this.config.n_inlets, - this.config.kla, - this.config.initialState - ); + new_reactor = new Reactor_PFR(this.config); break; default: console.warn("Unknown reactor type: " + uiConfig.reactor_type); diff --git a/src/reactor_class.js b/src/reactor_class.js index a81ff7f..220fbf9 100644 --- a/src/reactor_class.js +++ b/src/reactor_class.js @@ -9,16 +9,16 @@ const math = create(all, config) class Reactor { - constructor(volume, n_inlets, kla){ + constructor(config){ this.asm = new ASM3(); - this.Vl = volume; // fluid volume reactor [m3] + this.Vl = config.volume; // fluid volume reactor [m3] - this.Fs = Array(n_inlets).fill(0.0); // fluid debits per inlet [m3 d-1] - this.Cs_in = Array.from(Array(n_inlets), () => new Array(13).fill(0.0)); // composition influents + 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(13).fill(0.0)); // composition influents this.OTR = 0.0; // oxygen transfer rate [g O2 d-1] - this.kla = kla; // if NaN, use external OTR [d-1] + this.kla = config.kla; // if NaN, use external OTR [d-1] this.currentTime = Date.now(); // milliseconds since epoch [ms] this.timeStep = 1/(24*60*15); // time step [d] @@ -48,9 +48,9 @@ class Reactor { class Reactor_CSTR extends Reactor { - constructor(volume, n_inlets, kla, initial_state) { - super(volume, n_inlets, kla); - this.state = initial_state; + constructor(config) { + super(config); + this.state = config.initialState; } get getEffluent() { // getter for Effluent, defaults to inlet 0 @@ -90,16 +90,16 @@ class Reactor_CSTR extends Reactor { class Reactor_PFR extends Reactor { - constructor(volume, length, resolution_L, n_inlets, kla, initial_state) { - super(volume, n_inlets, kla); + constructor(config) { + super(config); - this.length = length; // reactor length [m] - this.n_x = resolution_L; // number of slices - this.d_x = length / resolution_L; + this.length = config.length; // reactor length [m] + this.n_x = config.resolution_L; // number of slices - this.A = volume / length; // crosssectional area [m2] + this.d_x = this.length / this.n_x; + this.A = this.Vl / this.length; // crosssectional area [m2] - this.state = Array.from(Array(this.n_x), () => initial_state.slice()) + this.state = Array.from(Array(this.n_x), () => config.initialState.slice()) // console.log("Initial State: ") // console.log(this.state) From 09e7072d16e64b38e2fd6ff11cd77f961705289b Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Fri, 4 Jul 2025 13:52:28 +0200 Subject: [PATCH 09/18] Refactor documentation in nodeClass and reactor_class for clarity and consistency --- src/nodeClass.js | 4 +- src/reactor_class.js | 133 ++++++++++++++++++++++++++++--------------- 2 files changed, 90 insertions(+), 47 deletions(-) diff --git a/src/nodeClass.js b/src/nodeClass.js index 15a959b..b9485be 100644 --- a/src/nodeClass.js +++ b/src/nodeClass.js @@ -3,9 +3,9 @@ const { Reactor_CSTR, Reactor_PFR } = require('./reactor_class.js'); class nodeClass { /** - * Construct ReactorNode. + * Node-RED node class for advanced-reactor. * @param {object} uiConfig - Node-RED node configuration - * @param {object} RED - Node-RED runtime API + * @param {object} RED - Node-RED runtime API * @param {object} nodeInstance - Node-RED node instance * @param {string} nameOfNode - Name of the node */ diff --git a/src/reactor_class.js b/src/reactor_class.js index 220fbf9..245f194 100644 --- a/src/reactor_class.js +++ b/src/reactor_class.js @@ -8,7 +8,10 @@ const config = { const math = create(all, config) class Reactor { - + /** + * Reactor base class. + * @param {object} config - Configuration object containing reactor parameters. + */ constructor(config){ this.asm = new ASM3(); @@ -18,62 +21,91 @@ class Reactor { this.Cs_in = Array.from(Array(config.n_inlets), () => new Array(13).fill(0.0)); // composition influents this.OTR = 0.0; // oxygen transfer rate [g O2 d-1] - this.kla = config.kla; // if NaN, use external OTR [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.speedUpFactor = 60; + this.speedUpFactor = 60; // speed up factor for simulation, 60 means 1 minute per simulated second } - set setInfluent(input) { // setter for C_in (WIP) + /** + * Setter for influent data. + * @param {object} input - Input object (msg) containing payload with inlet index, flow rate, and concentrations. + */ + set setInfluent(input) { let index_in = input.payload.inlet; this.Fs[index_in] = input.payload.F; this.Cs_in[index_in] = input.payload.C; + // DEBUG // 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] + /** + * Setter for OTR (Oxygen Transfer Rate). + * @param {object} input - Input object (msg) containing payload with OTR value [g O2 d-1]. + */ + set setOTR(input) { this.OTR = input.payload; } + /** + * + * @param {number} S_O - Dissolved oxygen concentration [g O2 m-3]. + * @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; return this.kla * (S_O_sat - S_O); } -} - -class Reactor_CSTR extends Reactor { - - constructor(config) { - super(config); - this.state = config.initialState; - } - - 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}; - } - - // expect update with timestamp - updateState(newTime) { - + /** + * Update the reactor state based on the new time. + * @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; let n_iter = Math.floor(this.speedUpFactor*(newTime - this.currentTime) / (this.timeStep * day2ms)); if (n_iter) { let n = 0; while (n < n_iter) { - this.tick_fe(this.timeStep); + this.tick(this.timeStep); n += 1; } this.currentTime += n_iter * this.timeStep * day2ms / this.speedUpFactor; } } - tick_fe(time_step) { // tick reactor state using forward Euler method +} + +class Reactor_CSTR extends Reactor { + /** + * Reactor_CSTR class for Continuous Stirred Tank Reactor. + * @param {object} config - Configuration object containing reactor parameters. + */ + constructor(config) { + super(config); + this.state = config.initialState; + } + + /** + * Getter for effluent data. + * @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}; + } + + /** + * Tick the reactor state using the forward Euler method. + * @param {number} time_step - Time step for the simulation [d]. + * @returns {Array} - New reactor state. + */ + 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); @@ -89,7 +121,10 @@ class Reactor_CSTR extends Reactor { } class Reactor_PFR extends Reactor { - + /** + * Reactor_PFR class for Plug Flow Reactor. + * @param {object} config - Configuration object containing reactor parameters. + */ constructor(config) { super(config); @@ -109,31 +144,29 @@ class Reactor_PFR extends Reactor { this.D_op = this.makeDoperator(true, true); this.D2_op = this.makeD2operator(); } - - set setDispersion(input) { // setter for Axial dispersion [m2 d-1] + + /** + * Setter for axial dispersion. + * @param {object} input - Input object (msg) containing payload with dispersion value [m2 d-1]. + */ + set setDispersion(input) { this.D = input.payload; } - get getEffluent() { // getter for Effluent, defaults to inlet 0 + /** + * Getter for effluent data. + * @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}; } - // 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)); - if (n_iter) { - let n = 0; - while (n < n_iter) { - this.tick_fe(this.timeStep); - n += 1; - } - this.currentTime += n_iter * this.timeStep * day2ms / this.speedUpFactor; - } - } - - tick_fe(time_step) { // tick reactor state using forward Euler method + /** + * Tick the reactor state using explicit finite difference method. + * @param {number} time_step - Time step for the simulation [d]. + * @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 reaction = this.state.map((state_slice) => this.asm.compute_dC(state_slice)); @@ -185,6 +218,12 @@ class Reactor_PFR extends Reactor { return new_state; } + /** + * Create finite difference first derivative operator. + * @param {boolean} central - Use central difference scheme if true, otherwise use upwind scheme. + * @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 if (higher_order) { if (central) { @@ -218,6 +257,10 @@ 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 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]); @@ -230,7 +273,7 @@ class Reactor_PFR extends Reactor { } -// testing stuff +// 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 // let initial_state = [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]; // const Reactor = new Reactor_PFR(200, 10, 10, 1, 100, initial_state); @@ -239,7 +282,7 @@ class Reactor_PFR extends Reactor { // Reactor.D = 0.01; // let N = 0; // while (N < 5000) { -// console.log(Reactor.tick_fe(0.001)); +// console.log(Reactor.tick(0.001)); // N += 1; // } From 3f5b0eea32e055a3f3ff518ee92474642ea4c987 Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Fri, 4 Jul 2025 14:58:38 +0200 Subject: [PATCH 10/18] 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 From c6b0cab06785efb050e4a978728d374321c17c0e Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Fri, 4 Jul 2025 15:14:03 +0200 Subject: [PATCH 11/18] 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; } } From 4b49d1076368015a85b2eab2b38cee006aba4475 Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Fri, 4 Jul 2025 15:48:05 +0200 Subject: [PATCH 12/18] Fixed bug in NaN assertion --- src/reactor_class.js | 52 +++++++++++++++++++++++++------------------- 1 file changed, 30 insertions(+), 22 deletions(-) diff --git a/src/reactor_class.js b/src/reactor_class.js index d483201..d41fd39 100644 --- a/src/reactor_class.js +++ b/src/reactor_class.js @@ -13,8 +13,14 @@ const math = create(all, config) * @param {string} label */ function assertNoNaN(arr, label = "array") { - if (math.isNaN(arr)) { - throw new Error("NaN detected in ${label}!"); + if (Array.isArray(arr)) { + for (const el of arr) { + assertNoNaN(el, label); + } + } else { + if (Number.isNaN(arr)) { + throw new Error(`NaN detected in ${label}!`); + } } } @@ -80,7 +86,7 @@ class Reactor { */ _arrayClip2Zero(arr) { if (Array.isArray(arr)) { - return arr.map(clipToZero); + return arr.map(x => this._arrayClip2Zero(x)); } else { return arr < 0 ? 0 : arr; } @@ -166,6 +172,9 @@ class Reactor_PFR extends Reactor { this.D_op = this._makeDoperator(true, true); this.D2_op = this._makeD2operator(); + + assertNoNaN(this.D_op, "Derivative operator"); + assertNoNaN(this.D2_op, "Second derivative operator"); } /** @@ -184,7 +193,7 @@ class Reactor_PFR extends Reactor { return { topic: "Fluent", payload: { inlet: 0, F: math.sum(this.Fs), C: this.state.at(-1) }, timestamp: this.currentTime }; } - _applyBoundaryConditions(newState) { + _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]; @@ -192,14 +201,13 @@ class Reactor_PFR extends Reactor { 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); + 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) - newState[0] = newState[1]; + state[0] = state[1]; } // Neumann BC (no flux) - newState[this.n_x - 1] = newState[this.n_x - 2] - return newState + state[this.n_x - 1] = state[this.n_x - 2] } /** @@ -224,11 +232,11 @@ class Reactor_PFR extends Reactor { } const dC_total = math.multiply(math.add(dispersion, advection, reaction, transfer), time_step); - let stateNew = math.add(this.state, dC_total); + const stateNew = math.add(this.state, dC_total); assertNoNaN(stateNew, "new state"); - stateNew = this._applyBoundaryConditions(stateNew); + this._applyBoundaryConditions(stateNew); assertNoNaN(stateNew, "new state post BC"); @@ -245,22 +253,22 @@ class Reactor_PFR extends Reactor { _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."); From 6755f2bd28b98dfd26f0107650bf950d0a9d94f5 Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Fri, 4 Jul 2025 16:03:42 +0200 Subject: [PATCH 13/18] Refactor reactor class to improve NaN handling and removed magic numbers --- src/reactor_class.js | 41 +++++++++++++++++++++++------------------ 1 file changed, 23 insertions(+), 18 deletions(-) diff --git a/src/reactor_class.js b/src/reactor_class.js index d41fd39..2062d8f 100644 --- a/src/reactor_class.js +++ b/src/reactor_class.js @@ -1,11 +1,14 @@ -const ASM3 = require('./asm3_class.js') -const { create, all } = require('mathjs') +const ASM3 = require('./asm3_class.js'); +const { create, all } = require('mathjs'); const config = { matrix: 'Array' // use Array as the matrix type -} +}; -const math = create(all, config) +const math = create(all, config); + +const OXYGEN_INDEX = 0; +const NUM_SPECIES = 13; /** * Assert that no NaN values are present in an array. @@ -35,7 +38,7 @@ class Reactor { this.Vl = 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(13).fill(0.0)); // composition influents + this.Cs_in = Array.from(Array(config.n_inlets), () => new Array(NUM_SPECIES).fill(0.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] @@ -136,15 +139,17 @@ class Reactor_CSTR extends Reactor { * @returns {Array} - New reactor state. */ 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 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 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 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 - const dC_total = math.multiply(math.add(dC_in, dC_out, r, t_O), time_step); + const dC_total = math.multiply(math.add(inflow, outflow, reaction, transfer), time_step); + assertNoNaN(dC_total, "change in state"); this.state = this._arrayClip2Zero(math.add(this.state, dC_total)); // clip value element-wise to avoid negative concentrations + assertNoNaN(this.state, "new state"); return this.state; } } @@ -171,9 +176,9 @@ 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(); - assertNoNaN(this.D_op, "Derivative operator"); + + this.D2_op = this._makeD2operator(); assertNoNaN(this.D2_op, "Second derivative operator"); } @@ -219,25 +224,25 @@ class Reactor_PFR extends Reactor { 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)); + const transfer = Array.from(Array(this.n_x), () => new Array(NUM_SPECIES).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; }); + transfer.forEach((x) => { x[OXYGEN_INDEX] = this.OTR; }); } else { - transfer.forEach((x, i) => { x[0] = this._calcOTR(this.state[i][0]); }); + transfer.forEach((x, i) => { x[OXYGEN_INDEX] = this._calcOTR(this.state[i][OXYGEN_INDEX]); }); } const dC_total = math.multiply(math.add(dispersion, advection, reaction, transfer), time_step); - const stateNew = math.add(this.state, dC_total); + assertNoNaN(dC_total, "change in state"); + const stateNew = math.add(this.state, dC_total); assertNoNaN(stateNew, "new state"); this._applyBoundaryConditions(stateNew); - assertNoNaN(stateNew, "new state post BC"); this.state = this._arrayClip2Zero(stateNew); From a2cfb20e2ca405671ca86b99c6b3a92b93ee6038 Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Fri, 4 Jul 2025 16:28:35 +0200 Subject: [PATCH 14/18] 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 From fe3add40071e9ecb1c9f02cedd9061a4e7139c12 Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Fri, 4 Jul 2025 17:42:31 +0200 Subject: [PATCH 15/18] Add debug assertions for state changes in Reactor classes --- src/reactor_class.js | 31 +++++++++++++++++-------------- 1 file changed, 17 insertions(+), 14 deletions(-) diff --git a/src/reactor_class.js b/src/reactor_class.js index dd5daf7..d12e158 100644 --- a/src/reactor_class.js +++ b/src/reactor_class.js @@ -10,6 +10,7 @@ const math = create(all, config); const S_O_INDEX = 0; const NUM_SPECIES = 13; +const DEBUG = false; class Reactor { /** @@ -129,11 +130,12 @@ class Reactor_CSTR extends Reactor { const transfer = Array(NUM_SPECIES).fill(0.0); 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"); - + const dC_total = math.multiply(math.add(inflow, outflow, reaction, transfer), time_step) this.state = this._arrayClip2Zero(math.add(this.state, dC_total)); // clip value element-wise to avoid negative concentrations - assertNoNaN(this.state, "new state"); + if(DEBUG){ + assertNoNaN(dC_total, "change in state"); + assertNoNaN(this.state, "new state"); + } return this.state; } } @@ -194,8 +196,9 @@ class Reactor_PFR extends Reactor { 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]; + let Pe = this.length * math.sum(this.Fs) / (this.D * this.A); + let residence_time = this.volume/math.sum(this.Fs); + const BC_dispersion = math.multiply((1 - (1 + 4*residence_time/Pe)^0.5) / (Pe*this.d_x), [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]; @@ -215,10 +218,6 @@ 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(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[S_O_INDEX] = this.OTR; }); } else { @@ -226,13 +225,17 @@ class Reactor_PFR extends Reactor { } const dC_total = math.multiply(math.add(dispersion, advection, reaction, transfer), time_step); - assertNoNaN(dC_total, "change in state"); const stateNew = math.add(this.state, dC_total); - assertNoNaN(stateNew, "new state"); - this._applyBoundaryConditions(stateNew); - assertNoNaN(stateNew, "new state post BC"); + + if (DEBUG) { + assertNoNaN(dispersion, "dispersion"); + assertNoNaN(advection, "advection"); + assertNoNaN(reaction, "reaction"); + assertNoNaN(dC_total, "change in state"); + assertNoNaN(stateNew, "new state post BC"); + } this.state = this._arrayClip2Zero(stateNew); return stateNew; From b4ddb6b8dfbb431bf23f9e07840f39018a400144 Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Mon, 7 Jul 2025 11:08:11 +0200 Subject: [PATCH 16/18] Enhance documentation for ASM3 class and its parameters --- src/asm3_class.js | 45 ++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 40 insertions(+), 5 deletions(-) diff --git a/src/asm3_class.js b/src/asm3_class.js index f5ac1df..ab272dc 100644 --- a/src/asm3_class.js +++ b/src/asm3_class.js @@ -1,11 +1,16 @@ const math = require('mathjs') +/** + * ASM3 class for the Activated Sludge Model No. 3 (ASM3). + */ class ASM3 { constructor() { + /** + * Kinetic parameters for ASM3 at 20 C. + * @property {Object} kin_params - Kinetic parameters + */ this.kin_params = { - // Kinetic parameters (20 C for now) - // Hydrolysis k_H: 3., // hydrolysis rate constant [g X_S g-1 X_H d-1] K_X: 1., // hydrolysis saturation constant [g X_S g-1 X_H] @@ -31,9 +36,13 @@ class ASM3 { b_A_O: 0.15, // aerobic respiration rate [d-1] b_A_NO: 0.05 // anoxic respiration rate [d-1] }; - this.stoi_params = { - // Stoichiometric and composition parameters + /** + * Stoichiometric and composition parameters for ASM3. + * @property {Object} stoi_params - Stoichiometric parameters + */ + this.stoi_params = { + // Fractions f_SI: 0., // fraction S_I from hydrolysis [g S_I g-1 X_S] f_XI: 0.2, // fraction X_I from decomp X_H [g X_I g-1 X_H] // Yields @@ -63,6 +72,10 @@ class ASM3 { this.stoi_matrix = this._initialise_stoi_matrix(); } + /** + * Initialises the stoichiometric matrix for ASM3. + * @returns {Array} - The stoichiometric matrix for ASM3. (2D array) + */ _initialise_stoi_matrix() { // initialise stoichiometric matrix const { f_SI, f_XI, Y_STO_O, Y_STO_NO, Y_H_O, Y_H_NO, Y_A, i_CODN, i_CODNO, i_NSI, i_NSS, i_NXI, i_NXS, i_NBM, i_TSXI, i_TSXS, i_TSBM, i_TSSTO, i_cNH, i_cNO } = this.stoi_params; @@ -84,15 +97,32 @@ class ASM3 { return stoi_matrix[0].map((col, i) => stoi_matrix.map(row => row[i])); // transpose matrix } + /** + * Computes the Monod equation rate value for a given concentration and half-saturation constant. + * @param {number} c - Concentration of reaction species. + * @param {number} K - Half-saturation constant for the reaction species. + * @returns {number} - Monod equation rate value for the given concentration and half-saturation constant. + */ _monod(c, K){ return c / (K + c); } + /** + * Computes the inverse Monod equation rate value for a given concentration and half-saturation constant. Used for inhibition. + * @param {number} c - Concentration of reaction species. + * @param {number} K - Half-saturation constant for the reaction species. + * @returns {number} - Inverse Monod equation rate value for the given concentration and half-saturation constant. + */ _inv_monod(c, K){ return K / (K + c); } - compute_rates(state) { // computes reaction rates. state is optional + /** + * Computes the reaction rates for each process reaction based on the current state. + * @param {Array} state - State vector containing concentrations of reaction species. + * @returns {Array} - Reaction rates for each process reaction. + */ + compute_rates(state) { // 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 const rates = Array(12); const [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] = state; @@ -119,6 +149,11 @@ class ASM3 { return rates; } + /** + * Computes the change in concentrations of reaction species based on the current state. + * @param {Array} state - State vector containing concentrations of reaction species. + * @returns {Array} - Change in reaction species concentrations. + */ compute_dC(state) { // compute changes in concentrations // 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 return math.multiply(this.stoi_matrix, this.compute_rates(state)); From deb5269d1a59d7eeb14e19ad7aed9c6514d08537 Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Mon, 7 Jul 2025 11:59:11 +0200 Subject: [PATCH 17/18] Change file structure to align with project --- src/nodeClass.js | 2 +- src/{ => reaction_modules}/asm3_class.js | 0 src/{reactor_class.js => specificClass.js} | 2 +- 3 files changed, 2 insertions(+), 2 deletions(-) rename src/{ => reaction_modules}/asm3_class.js (100%) rename src/{reactor_class.js => specificClass.js} (99%) diff --git a/src/nodeClass.js b/src/nodeClass.js index a873284..7919df2 100644 --- a/src/nodeClass.js +++ b/src/nodeClass.js @@ -1,4 +1,4 @@ -const { Reactor_CSTR, Reactor_PFR } = require('./reactor_class.js'); +const { Reactor_CSTR, Reactor_PFR } = require('./specificClass.js'); class nodeClass { diff --git a/src/asm3_class.js b/src/reaction_modules/asm3_class.js similarity index 100% rename from src/asm3_class.js rename to src/reaction_modules/asm3_class.js diff --git a/src/reactor_class.js b/src/specificClass.js similarity index 99% rename from src/reactor_class.js rename to src/specificClass.js index d12e158..ce5c786 100644 --- a/src/reactor_class.js +++ b/src/specificClass.js @@ -1,4 +1,4 @@ -const ASM3 = require('./asm3_class.js'); +const ASM3 = require('./reaction_modules/asm3_class.js'); const { create, all } = require('mathjs'); const { assertNoNaN } = require('./utils.js'); From 302780726a8b6d850f7dd4b64440c83d63dcbd41 Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Mon, 7 Jul 2025 12:24:15 +0200 Subject: [PATCH 18/18] Refactor Reactor class to remove debug logs and enhance setter for influent data with conditional logging --- src/specificClass.js | 23 ++++++++++++++++------- 1 file changed, 16 insertions(+), 7 deletions(-) diff --git a/src/specificClass.js b/src/specificClass.js index ce5c786..384aaf0 100644 --- a/src/specificClass.js +++ b/src/specificClass.js @@ -41,11 +41,6 @@ class Reactor { let index_in = input.payload.inlet; this.Fs[index_in] = input.payload.F; this.Cs_in[index_in] = input.payload.C; - // DEBUG - // 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)); } /** @@ -57,10 +52,10 @@ class Reactor { } /** - * + * Calculate the oxygen transfer rate (OTR) based on the dissolved oxygen concentration and temperature. * @param {number} S_O - Dissolved oxygen concentration [g O2 m-3]. * @param {number} T - Temperature in Celsius, default to 20 C. - * @returns + * @returns {number} - Calculated OTR [g O2 d-1]. */ _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; @@ -176,6 +171,20 @@ class Reactor_PFR extends Reactor { this.D = input.payload; } + /** + * Setter for influent data. + * @param {object} input - Input object (msg) containing payload with inlet index, flow rate, and concentrations. + */ + set setInfluent(input) { + super.setInfluent = input; + if(DEBUG) { + 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)); + } + } + /** * Getter for effluent data. * @returns {object} Effluent data object (msg), defaults to inlet 0.