Date: Mon, 16 Jun 2025 14:01:19 +0200
Subject: [PATCH 18/97] Add number of inlets input handling to advanced-reactor
node
---
advanced-reactor.html | 13 +++++++++++++
advanced-reactor.js | 13 ++++++-------
dependencies/reactor_class.js | 25 +++++++++++++------------
3 files changed, 32 insertions(+), 19 deletions(-)
diff --git a/advanced-reactor.html b/advanced-reactor.html
index a18e6a3..a1a7ea7 100644
--- a/advanced-reactor.html
+++ b/advanced-reactor.html
@@ -5,6 +5,7 @@
defaults: {
name: { value: "" },
volume: { value: 0., required: true},
+ n_inlets: { value: 1, required: true},
S_O_init: { value: 0., required: true },
S_I_init: { value: 30., required: true },
S_S_init: { value: 100., required: true },
@@ -31,6 +32,10 @@
type:"num",
types:["num"]
});
+ $("#node-input-n_inlets").typedInput({
+ type:"num",
+ types:["num"]
+ });
$(".concentrations").typedInput({
type:"num",
types:["num"]
@@ -41,6 +46,10 @@
if (isNaN(volume) || volume <= 0) {
RED.notify("Fluid volume not set correctly", {type: "error"});
}
+ let n_inlets = parseInt($("#node-input-n_inlets").typedInput("value"));
+ if (isNaN(n_inlets) || n_inlets < 1) {
+ RED.notify("Number of inlets not set correctly", {type: "error"});
+ }
}
});
@@ -55,6 +64,10 @@
+
+
+
+
Dissolved components
-
+
diff --git a/additional_nodes/recirculation-pump.js b/additional_nodes/recirculation-pump.js
index 489df01..cd9fc21 100644
--- a/additional_nodes/recirculation-pump.js
+++ b/additional_nodes/recirculation-pump.js
@@ -11,12 +11,12 @@ module.exports = function(RED) {
switch (msg.topic) {
case "Fluent":
// conserve volume flow debit
- let F1 = msg.payload.F;
- let F_diff = Math.max(F1 - F2, 0);
- let F2_corr = F1 < F2 ? F1 : F2;
+ let F_in = msg.payload.F;
+ let F1 = Math.max(F_in - F2, 0);
+ let F2_corr = F_in < F2 ? F_in : F2;
let msg_F1 = structuredClone(msg);
- msg_F1.payload.F = F_diff;
+ msg_F1.payload.F = F1;
let msg_F2 = {...msg};
msg_F2.payload.F = F2_corr;
diff --git a/additional_nodes/settling-basin.html b/additional_nodes/settling-basin.html
index da614bc..e8e8e8d 100644
--- a/additional_nodes/settling-basin.html
+++ b/additional_nodes/settling-basin.html
@@ -4,7 +4,7 @@
color: "#e4a363",
defaults: {
name: { value: "" },
- SVI: { value: 0.1, required: true },
+ TS_set: { value: 0.1, required: true },
inlet: { value: 1, required: true }
},
inputs: 1,
@@ -15,7 +15,7 @@
return this.name || "Settling basin";
},
oneditprepare: function() {
- $("#node-input-SVI").typedInput({
+ $("#node-input-TS_set").typedInput({
type:"num",
types:["num"]
});
@@ -25,9 +25,9 @@
});
},
oneditsave: function() {
- let SVI = parseFloat($("#node-input-SVI").typedInput("value"));
- if (isNaN(SVI) || SVI < 0) {
- RED.notify("SVI is not set correctly", {type: "error"});
+ let TS_set = parseFloat($("#node-input-TS_set").typedInput("value"));
+ if (isNaN(TS_set) || TS_set < 0) {
+ RED.notify("TS is not set correctly", {type: "error"});
}
let inlet = parseInt($("#node-input-n_inlets").typedInput("value"));
if (inlet < 1) {
@@ -43,8 +43,8 @@
-
-
+
+
Reactor properties
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/advanced-reactor.js b/advanced-reactor.js
index edfad14..68e9674 100644
--- a/advanced-reactor.js
+++ b/advanced-reactor.js
@@ -7,26 +7,60 @@ module.exports = function(RED) {
const Reactor = require('./dependencies/reactor_class');
- const reactor = new Reactor(
- 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)
- ]
- );
+ let new_reactor;
+
+ switch (config.reactor_type) {
+ case "CSTR":
+ new_reactor = new Reactor(
+ 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(
+ parseFloat(config.volume),
+ parseFloat(config.L),
+ 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;
diff --git a/dependencies/reactor_class.js b/dependencies/reactor_class.js
index 1b2e3d9..aa7f9dd 100644
--- a/dependencies/reactor_class.js
+++ b/dependencies/reactor_class.js
@@ -5,7 +5,6 @@ class Reactor_CSTR {
constructor(volume, n_inlets, kla, initial_state) {
this.state = initial_state;
- console.log(this.state);
this.asm = new ASM3();
this.Vl = volume; // fluid volume reactor [m3]
@@ -17,7 +16,7 @@ class Reactor_CSTR {
this.currentTime = Date.now(); // milliseconds since epoch [ms]
this.timeStep = 1/(24*60*15); // time step [d]
- this.speedUpFactor = 30;
+ this.speedUpFactor = 1;
}
set setInfluent(input) { // setter for C_in (WIP)
@@ -69,6 +68,113 @@ class Reactor_CSTR {
}
}
+class Reactor_PFR {
+
+ constructor(volume, length, resolution_L, n_inlets, kla, initial_state) {
+ this.asm = new ASM3();
+
+ 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;
+
+ this.A = volume / length; // crosssectional area [m2]
+
+ this.state = Array.from(Array(this.n_x), () => initial_state.slice())
+
+ 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 = 1;
+
+ this.D_op = makeDoperator();
+ this.D2_op = 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;
+ }
+
+ 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;
+ }
+
+ 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);
+ }
+
+ // 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
+ if (math.sum(this.Fs) > 0) {
+ this.state[0] = math.multiply(math.divide([this.Fs], this.A), this.Cs_in)[0] // Dichelet boundary condition
+ }
+
+ const dispersion = math.multiply(this.D / (this.d_x*this.d_x), this.D2_op, this.state);
+ const advection = math.multiply(math.sum(this.Fs)/(this.A*this.d_x), this.D_op, this.state);
+ const reaction = this.state.map(this.asm.compute_dC);
+ const transfer = Array.from(Array(this.n_x), () => new Array(13).fill(0.0))
+
+ 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]); });
+ }
+
+ const dC_total = math.multiply(math.add(dispersion, advection, reaction, transfer), time_step);
+
+ this.state = math.abs(math.add(this.state, dC_total)); // make sure that concentrations do not go negative
+ return this.state;
+ }
+
+ makeDoperator() { // create the upwind scheme gradient operator
+ const I = math.identity(this.n_x);
+ const A = math.diag(Array(this.n_x).fill(-1), 1).resize([this.n_x, this.n_x]);
+ I[this.n_x-1, this.n_x-1] = 0; // Neumann boundary condition at x=L
+ return math.add(I, A);
+ }
+
+ makeD2operator() { // create the upwind scheme second derivative operator
+ const I = math.diag(Array(this.n_x).fill(2), 0);
+ const A = math.diag(Array(this.n_x).fill(-1), 1).resize([this.n_x, this.n_x]);
+ const B = math.diag(Array(this.n_x).fill(-1), -1).resize([this.n_x, this.n_x]);
+ I[0, 0] = 1;
+ return math.add(I, A, B);
+ }
+}
+
+
// testing stuff
// 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];
From e6c1e21c166ac90751e8e320344bf741b9b4bfc8 Mon Sep 17 00:00:00 2001
From: "p.vanderwilt"
Date: Mon, 23 Jun 2025 17:46:55 +0200
Subject: [PATCH 29/97] Implement Danckwerts boundary condition in tick_fe
method for Reactor_PFR
---
dependencies/reactor_class.js | 21 +++++++++++++++------
1 file changed, 15 insertions(+), 6 deletions(-)
diff --git a/dependencies/reactor_class.js b/dependencies/reactor_class.js
index aa7f9dd..bb5a96c 100644
--- a/dependencies/reactor_class.js
+++ b/dependencies/reactor_class.js
@@ -137,14 +137,11 @@ class Reactor_PFR {
}
tick_fe(time_step) { // tick reactor state using forward Euler method
- if (math.sum(this.Fs) > 0) {
- this.state[0] = math.multiply(math.divide([this.Fs], this.A), this.Cs_in)[0] // Dichelet boundary condition
- }
-
const dispersion = math.multiply(this.D / (this.d_x*this.d_x), this.D2_op, this.state);
const advection = math.multiply(math.sum(this.Fs)/(this.A*this.d_x), this.D_op, this.state);
const reaction = this.state.map(this.asm.compute_dC);
- const transfer = Array.from(Array(this.n_x), () => new Array(13).fill(0.0))
+ reaction[0] = Array(13).fill(0.0);
+ const transfer = Array.from(Array(this.n_x), () => new Array(13).fill(0.0));
if (isNaN(this.kla)) { // calculate OTR if kla is not NaN, otherwise use externally calculated OTR
transfer.forEach((x) => { x[0] = this.OTR; });
@@ -152,6 +149,15 @@ class Reactor_PFR {
transfer.forEach((x, i) => { x[0] = this.calcOTR(this.state[i][0]); });
}
+ if (math.sum(this.Fs) > 0) { // Danckwerts BC
+ const BC_influx = math.multiply(math.divide([this.Fs], this.A), this.Cs_in)[0];
+ const BC_gradient = Array(this.n_x).fill(0.0);
+ BC_gradient[0] = 1;
+ BC_gradient[1] = -1;
+ const BC_dispersion = math.multiply(this.D * this.A / (math.sum(this.Fs)*this.d_x), [BC_gradient], this.state);
+ this.state[0] = math.add(BC_influx, BC_dispersion);
+ }
+
const dC_total = math.multiply(math.add(dispersion, advection, reaction, transfer), time_step);
this.state = math.abs(math.add(this.state, dC_total)); // make sure that concentrations do not go negative
@@ -161,6 +167,8 @@ class Reactor_PFR {
makeDoperator() { // create the upwind scheme gradient operator
const I = math.identity(this.n_x);
const A = math.diag(Array(this.n_x).fill(-1), 1).resize([this.n_x, this.n_x]);
+ I[0, 0] = 0;
+ I[0, 1] = 0;
I[this.n_x-1, this.n_x-1] = 0; // Neumann boundary condition at x=L
return math.add(I, A);
}
@@ -169,7 +177,8 @@ class Reactor_PFR {
const I = math.diag(Array(this.n_x).fill(2), 0);
const A = math.diag(Array(this.n_x).fill(-1), 1).resize([this.n_x, this.n_x]);
const B = math.diag(Array(this.n_x).fill(-1), -1).resize([this.n_x, this.n_x]);
- I[0, 0] = 1;
+ I[0, 0] = 0;
+ I[0, 1] = 0;
return math.add(I, A, B);
}
}
From e5c9010093e2cf8df67cead7d4ded041b1e758d1 Mon Sep 17 00:00:00 2001
From: "p.vanderwilt"
Date: Tue, 24 Jun 2025 11:20:28 +0200
Subject: [PATCH 30/97] Fixed various bugs
---
dependencies/asm3_class.js | 118 +++++++++++++++++-----------------
dependencies/reactor_class.js | 41 +++++++-----
2 files changed, 82 insertions(+), 77 deletions(-)
diff --git a/dependencies/asm3_class.js b/dependencies/asm3_class.js
index 9ad03d9..9b7762a 100644
--- a/dependencies/asm3_class.js
+++ b/dependencies/asm3_class.js
@@ -2,67 +2,65 @@ const math = require('mathjs')
class ASM3 {
- 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]
- // Heterotrophs
- k_STO: 5., // storage rate constant [g S_S g-1 X_H d-1]
- nu_NO: 0.6, // anoxic reduction factor [-]
- K_O: 0.2, // saturation constant S_0 [g O2 m-3]
- K_NO: 0.5, // saturation constant S_NO [g NO3-N m-3]
- K_S: 2., // saturation constant S_s [g COD m-3]
- K_STO: 1., // saturation constant X_STO [g X_STO g-1 X_H]
- mu_H_max: 2., // maximum specific growth rate [d-1]
- K_NH: 0.01, // saturation constant S_NH3 [g NH3-N m-3]
- K_HCO: 0.1, // saturation constant S_HCO [mole HCO3 m-3]
- b_H_O: 0.2, // aerobic respiration rate [d-1]
- b_H_NO: 0.1, // anoxic respiration rate [d-1]
- b_STO_O: 0.2, // aerobic respitation rate X_STO [d-1]
- b_STO_NO: 0.1, // anoxic respitation rate X_STO [d-1]
- // Autotrophs
- mu_A_max: 1.0, // maximum specific growth rate [d-1]
- K_A_NH: 1., // saturation constant S_NH3 [g NH3-N m-3]
- K_A_O: 0.5, // saturation constant S_0 [g O2 m-3]
- K_A_HCO: 0.5, // saturation constant S_HCO [mole HCO3 m-3]
- b_A_O: 0.15, // aerobic respiration rate [d-1]
- b_A_NO: 0.05 // anoxic respiration rate [d-1]
- }
-
- stoi_params = {
- // Stoichiometric and composition parameters
-
- 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
- Y_STO_O: 0.85, // aerobic yield X_STO per S_S [g X_STO g-1 S_S]
- Y_STO_NO: 0.80, // anoxic yield X_STO per S_S [g X_STO g-1 S_S]
- Y_H_O: 0.63, // aerobic yield X_H per X_STO [g X_H g-1 X_STO]
- Y_H_NO: 0.54, // anoxic yield X_H per X_STO [g X_H g-1 X_STO]
- Y_A: 0.24, // anoxic yield X_A per S_NO [g X_A g-1 NO3-N]
- // Composition (COD via DoR)
- i_CODN: -1.71, // COD content (DoR) [g COD g-1 N2-N]
- i_CODNO: -4.57, // COD content (DoR) [g COD g-1 NO3-N]
- // Composition (nitrogen)
- i_NSI: 0.01, // nitrogen content S_I [g N g-1 S_I]
- i_NSS: 0.03, // nitrogen content S_S [g N g-1 S_S]
- i_NXI: 0.02, // nitrogen content X_I [g N g-1 X_I]
- i_NXS: 0.04, // nitrogen content X_S [g N g-1 X_S]
- i_NBM: 0.07, // nitrogen content X_H / X_A [g N g-1 X_H / X_A]
- // Composition (TSS)
- i_TSXI: 0.75, // TSS content X_I [g TS g-1 X_I]
- i_TSXS: 0.75, // TSS content X_S [g TS g-1 X_S]
- i_TSBM: 0.90, // TSS content X_H / X_A [g TS g-1 X_H / X_A]
- i_TSSTO: 0.60, // TSS content X_STO (PHB based) [g TS g-1 X_STO]
- // Composition (charge)
- i_cNH: 1/14, // charge per S_NH [mole H+ g-1 NH3-N]
- i_cNO: -1/14 // charge per S_NO [mole H+ g-1 NO3-N]
- }
-
constructor() {
- this.stoi_matrix = this._initialise_stoi_matrix()
+ 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]
+ // Heterotrophs
+ k_STO: 5., // storage rate constant [g S_S g-1 X_H d-1]
+ nu_NO: 0.6, // anoxic reduction factor [-]
+ K_O: 0.2, // saturation constant S_0 [g O2 m-3]
+ K_NO: 0.5, // saturation constant S_NO [g NO3-N m-3]
+ K_S: 2., // saturation constant S_s [g COD m-3]
+ K_STO: 1., // saturation constant X_STO [g X_STO g-1 X_H]
+ mu_H_max: 2., // maximum specific growth rate [d-1]
+ K_NH: 0.01, // saturation constant S_NH3 [g NH3-N m-3]
+ K_HCO: 0.1, // saturation constant S_HCO [mole HCO3 m-3]
+ b_H_O: 0.2, // aerobic respiration rate [d-1]
+ b_H_NO: 0.1, // anoxic respiration rate [d-1]
+ b_STO_O: 0.2, // aerobic respitation rate X_STO [d-1]
+ b_STO_NO: 0.1, // anoxic respitation rate X_STO [d-1]
+ // Autotrophs
+ mu_A_max: 1.0, // maximum specific growth rate [d-1]
+ K_A_NH: 1., // saturation constant S_NH3 [g NH3-N m-3]
+ K_A_O: 0.5, // saturation constant S_0 [g O2 m-3]
+ K_A_HCO: 0.5, // saturation constant S_HCO [mole HCO3 m-3]
+ 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
+
+ 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
+ Y_STO_O: 0.85, // aerobic yield X_STO per S_S [g X_STO g-1 S_S]
+ Y_STO_NO: 0.80, // anoxic yield X_STO per S_S [g X_STO g-1 S_S]
+ Y_H_O: 0.63, // aerobic yield X_H per X_STO [g X_H g-1 X_STO]
+ Y_H_NO: 0.54, // anoxic yield X_H per X_STO [g X_H g-1 X_STO]
+ Y_A: 0.24, // anoxic yield X_A per S_NO [g X_A g-1 NO3-N]
+ // Composition (COD via DoR)
+ i_CODN: -1.71, // COD content (DoR) [g COD g-1 N2-N]
+ i_CODNO: -4.57, // COD content (DoR) [g COD g-1 NO3-N]
+ // Composition (nitrogen)
+ i_NSI: 0.01, // nitrogen content S_I [g N g-1 S_I]
+ i_NSS: 0.03, // nitrogen content S_S [g N g-1 S_S]
+ i_NXI: 0.02, // nitrogen content X_I [g N g-1 X_I]
+ i_NXS: 0.04, // nitrogen content X_S [g N g-1 X_S]
+ i_NBM: 0.07, // nitrogen content X_H / X_A [g N g-1 X_H / X_A]
+ // Composition (TSS)
+ i_TSXI: 0.75, // TSS content X_I [g TS g-1 X_I]
+ i_TSXS: 0.75, // TSS content X_S [g TS g-1 X_S]
+ i_TSBM: 0.90, // TSS content X_H / X_A [g TS g-1 X_H / X_A]
+ i_TSSTO: 0.60, // TSS content X_STO (PHB based) [g TS g-1 X_STO]
+ // Composition (charge)
+ i_cNH: 1/14, // charge per S_NH [mole H+ g-1 NH3-N]
+ i_cNO: -1/14 // charge per S_NO [mole H+ g-1 NO3-N]
+ };
+ this.stoi_matrix = this._initialise_stoi_matrix();
}
_initialise_stoi_matrix() { // initialise stoichiometric matrix
diff --git a/dependencies/reactor_class.js b/dependencies/reactor_class.js
index bb5a96c..af3289e 100644
--- a/dependencies/reactor_class.js
+++ b/dependencies/reactor_class.js
@@ -1,5 +1,11 @@
const ASM3 = require('./asm3_class')
-const math = require('mathjs')
+const { create, all } = require('mathjs')
+
+const config = {
+ matrix: 'Array' // Choose 'Matrix' (default) or 'Array'
+}
+
+const math = create(all, config)
class Reactor_CSTR {
@@ -93,8 +99,8 @@ class Reactor_PFR {
this.timeStep = 1/(24*60*15); // time step [d]
this.speedUpFactor = 1;
- this.D_op = makeDoperator();
- this.D2_op = makeD2operator();
+ this.D_op = this.makeDoperator();
+ this.D2_op = this.makeD2operator();
}
set setInfluent(input) { // setter for C_in (WIP)
@@ -139,7 +145,7 @@ class Reactor_PFR {
tick_fe(time_step) { // tick reactor state using forward Euler method
const dispersion = math.multiply(this.D / (this.d_x*this.d_x), this.D2_op, this.state);
const advection = math.multiply(math.sum(this.Fs)/(this.A*this.d_x), this.D_op, this.state);
- const reaction = this.state.map(this.asm.compute_dC);
+ const reaction = this.state.map((row) => this.asm.compute_dC(row));
reaction[0] = Array(13).fill(0.0);
const transfer = Array.from(Array(this.n_x), () => new Array(13).fill(0.0));
@@ -154,7 +160,7 @@ class Reactor_PFR {
const BC_gradient = Array(this.n_x).fill(0.0);
BC_gradient[0] = 1;
BC_gradient[1] = -1;
- const BC_dispersion = math.multiply(this.D * this.A / (math.sum(this.Fs)*this.d_x), [BC_gradient], this.state);
+ const BC_dispersion = math.multiply(this.D * this.A / (math.sum(this.Fs)*this.d_x), [BC_gradient], this.state)[0];
this.state[0] = math.add(BC_influx, BC_dispersion);
}
@@ -166,19 +172,19 @@ class Reactor_PFR {
makeDoperator() { // create the upwind scheme gradient operator
const I = math.identity(this.n_x);
- const A = math.diag(Array(this.n_x).fill(-1), 1).resize([this.n_x, this.n_x]);
- I[0, 0] = 0;
- I[0, 1] = 0;
- I[this.n_x-1, this.n_x-1] = 0; // Neumann boundary condition at x=L
+ const A = math.resize(math.diag(Array(this.n_x).fill(-1), 1), [this.n_x, this.n_x]);
+ I[0][0] = 0;
+ I[0][1] = 1;
+ I[this.n_x-1][this.n_x-1] = 0; // Neumann boundary condition at x=L
return math.add(I, A);
}
makeD2operator() { // create the upwind scheme second derivative operator
- const I = math.diag(Array(this.n_x).fill(2), 0);
- const A = math.diag(Array(this.n_x).fill(-1), 1).resize([this.n_x, this.n_x]);
- const B = math.diag(Array(this.n_x).fill(-1), -1).resize([this.n_x, this.n_x]);
- I[0, 0] = 0;
- I[0, 1] = 0;
+ const I = math.identity(this.n_x);
+ 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]);
+ I[0][0] = 0;
+ I[0][1] = 1;
return math.add(I, A, B);
}
}
@@ -187,9 +193,10 @@ class Reactor_PFR {
// testing stuff
// 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_CSTR(initial_state);
-// Reactor.C_in = [0.0, 30., 100., 16., 0., 0., 5., 25., 75., 30., 0., 0., 125.];
-// N = 0;
+// const Reactor = new Reactor_PFR(200, 10, 10, 1, 100, initial_state);
+// Reactor.Cs_in[0] = [0.0, 30., 100., 16., 0., 0., 5., 25., 75., 30., 0., 0., 125.];
+// Reactor.Fs[0] = 10;
+// let N = 0;
// while (N < 500) {
// console.log(Reactor.tick_fe(0.001));
// N += 1;
From 6b57a46aabd54511c7f3db31b9d480062bb18272 Mon Sep 17 00:00:00 2001
From: "p.vanderwilt"
Date: Tue, 24 Jun 2025 12:32:11 +0200
Subject: [PATCH 31/97] Add typed input fields for reactor length and
resolution in advanced-reactor, fixed NaN bug in reactor length
---
advanced-reactor.html | 8 ++++++++
advanced-reactor.js | 8 ++++----
dependencies/reactor_class.js | 15 +++++++++------
3 files changed, 21 insertions(+), 10 deletions(-)
diff --git a/advanced-reactor.html b/advanced-reactor.html
index 5b4d41b..512f3dd 100644
--- a/advanced-reactor.html
+++ b/advanced-reactor.html
@@ -40,6 +40,14 @@
type:"num",
types:["num"]
});
+ $("#node-input-length").typedInput({
+ type:"num",
+ types:["num"]
+ });
+ $("#node-input-resolution_L").typedInput({
+ type:"num",
+ types:["num"]
+ });
$("#node-input-kla").typedInput({
type:"num",
types:["num"]
diff --git a/advanced-reactor.js b/advanced-reactor.js
index 68e9674..8348137 100644
--- a/advanced-reactor.js
+++ b/advanced-reactor.js
@@ -5,13 +5,13 @@ module.exports = function(RED) {
let name = config.name;
- const Reactor = require('./dependencies/reactor_class');
+ const { Reactor_CSTR, Reactor_PFR } = require('./dependencies/reactor_class');
let new_reactor;
switch (config.reactor_type) {
case "CSTR":
- new_reactor = new Reactor(
+ new_reactor = new Reactor_CSTR(
parseFloat(config.volume),
parseInt(config.n_inlets),
parseFloat(config.kla),
@@ -33,9 +33,9 @@ module.exports = function(RED) {
);
break;
case "PFR":
- new_reactor = new Reactor(
+ new_reactor = new Reactor_PFR(
parseFloat(config.volume),
- parseFloat(config.L),
+ parseFloat(config.length),
parseInt(config.resolution_L),
parseInt(config.n_inlets),
parseFloat(config.kla),
diff --git a/dependencies/reactor_class.js b/dependencies/reactor_class.js
index af3289e..4027690 100644
--- a/dependencies/reactor_class.js
+++ b/dependencies/reactor_class.js
@@ -87,16 +87,19 @@ class Reactor_PFR {
this.A = volume / length; // crosssectional area [m2]
this.state = Array.from(Array(this.n_x), () => initial_state.slice())
+
+ // 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.D = 0.1; // 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.timeStep = 1/(24*60*60); // time step [d]
this.speedUpFactor = 1;
this.D_op = this.makeDoperator();
@@ -118,7 +121,7 @@ class Reactor_PFR {
}
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.at(-1)}, timestamp: this.currentTime};
}
calcOTR(S_O, T=20.0) { // caculate the OTR using basic correlation, default to temperature: 20 C
@@ -145,10 +148,10 @@ class Reactor_PFR {
tick_fe(time_step) { // tick reactor state using forward Euler method
const dispersion = math.multiply(this.D / (this.d_x*this.d_x), this.D2_op, this.state);
const advection = math.multiply(math.sum(this.Fs)/(this.A*this.d_x), this.D_op, this.state);
- const reaction = this.state.map((row) => this.asm.compute_dC(row));
+ const reaction = this.state.map((state_slice) => this.asm.compute_dC(state_slice));
reaction[0] = Array(13).fill(0.0);
const transfer = Array.from(Array(this.n_x), () => new Array(13).fill(0.0));
-
+
if (isNaN(this.kla)) { // calculate OTR if kla is not NaN, otherwise use externally calculated OTR
transfer.forEach((x) => { x[0] = this.OTR; });
} else {
@@ -202,4 +205,4 @@ class Reactor_PFR {
// N += 1;
// }
-module.exports = Reactor_CSTR;
\ No newline at end of file
+module.exports = {Reactor_CSTR, Reactor_PFR};
\ No newline at end of file
From f2d94b26c5141d976dac3e3cf22405fb8703cbd3 Mon Sep 17 00:00:00 2001
From: "p.vanderwilt"
Date: Tue, 24 Jun 2025 13:28:45 +0200
Subject: [PATCH 32/97] Add dispersion setting in advanced-reactor and
initialize axial dispersion to zero in Reactor_PFR
---
advanced-reactor.js | 3 +++
dependencies/reactor_class.js | 2 +-
2 files changed, 4 insertions(+), 1 deletion(-)
diff --git a/advanced-reactor.js b/advanced-reactor.js
index 8348137..e2eca51 100644
--- a/advanced-reactor.js
+++ b/advanced-reactor.js
@@ -78,6 +78,9 @@ module.exports = function(RED) {
case "OTR":
reactor.setOTR = msg;
break;
+ case "Dispersion":
+ reactor.setDispersion = msg;
+ break;
default:
console.log("Unknown topic: " + msg.topic);
}
diff --git a/dependencies/reactor_class.js b/dependencies/reactor_class.js
index 4027690..e030637 100644
--- a/dependencies/reactor_class.js
+++ b/dependencies/reactor_class.js
@@ -94,7 +94,7 @@ class Reactor_PFR {
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.1; // axial dispersion [m2 d-1]
+ this.D = 0.0; // axial dispersion [m2 d-1]
this.kla = kla; // if NaN, use external OTR [d-1]
From 2e76f733a8fed39a17f3f699a2c6b0506017a315 Mon Sep 17 00:00:00 2001
From: "p.vanderwilt"
Date: Tue, 24 Jun 2025 16:23:33 +0200
Subject: [PATCH 33/97] Working on fixing the Derivative operators and BCs
---
dependencies/reactor_class.js | 40 ++++++++++++++++++++---------------
1 file changed, 23 insertions(+), 17 deletions(-)
diff --git a/dependencies/reactor_class.js b/dependencies/reactor_class.js
index e030637..9e02fca 100644
--- a/dependencies/reactor_class.js
+++ b/dependencies/reactor_class.js
@@ -69,7 +69,8 @@ class Reactor_CSTR {
const dC_total = math.multiply(math.add(dC_in, dC_out, r, t_O), time_step);
- this.state = math.abs(math.add(this.state, dC_total)); // make sure that concentrations do not go negative
+ // 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);
return this.state;
}
}
@@ -99,8 +100,8 @@ class Reactor_PFR {
this.kla = kla; // if NaN, use external OTR [d-1]
this.currentTime = Date.now(); // milliseconds since epoch [ms]
- this.timeStep = 1/(24*60*60); // time step [d]
- this.speedUpFactor = 1;
+ this.timeStep = 1/(24*60*15); // time step [d]
+ this.speedUpFactor = 60;
this.D_op = this.makeDoperator();
this.D2_op = this.makeD2operator();
@@ -147,7 +148,7 @@ class Reactor_PFR {
tick_fe(time_step) { // tick reactor state using forward Euler method
const dispersion = math.multiply(this.D / (this.d_x*this.d_x), this.D2_op, this.state);
- const advection = math.multiply(math.sum(this.Fs)/(this.A*this.d_x), this.D_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));
reaction[0] = Array(13).fill(0.0);
const transfer = Array.from(Array(this.n_x), () => new Array(13).fill(0.0));
@@ -158,36 +159,40 @@ class Reactor_PFR {
transfer.forEach((x, i) => { x[0] = this.calcOTR(this.state[i][0]); });
}
+ transfer[0][0] = 0;
+
if (math.sum(this.Fs) > 0) { // Danckwerts BC
- const BC_influx = math.multiply(math.divide([this.Fs], this.A), this.Cs_in)[0];
+ 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;
+ BC_gradient[0] = -1;
+ BC_gradient[1] = 1;
const BC_dispersion = math.multiply(this.D * this.A / (math.sum(this.Fs)*this.d_x), [BC_gradient], this.state)[0];
- this.state[0] = math.add(BC_influx, BC_dispersion);
+ this.state[0] = math.add(BC_C_in, BC_dispersion);
+ console.log(BC_dispersion);
}
const dC_total = math.multiply(math.add(dispersion, advection, reaction, transfer), time_step);
- this.state = math.abs(math.add(this.state, dC_total)); // make sure that concentrations do not go negative
+ // clip value element-wise to each subarray to avoid negative concentrations
+ this.state = math.add(this.state, dC_total).map(row => row.map(val => val < 0 ? 0 : val));
return this.state;
}
makeDoperator() { // create the upwind scheme gradient operator
- const I = math.identity(this.n_x);
- const A = math.resize(math.diag(Array(this.n_x).fill(-1), 1), [this.n_x, this.n_x]);
+ const I = math.diag(Array(this.n_x).fill(-1), 0);
+ const A = math.resize(math.diag(Array(this.n_x).fill(1), 1), [this.n_x, this.n_x]);
I[0][0] = 0;
- I[0][1] = 1;
+ I[0][1] = -1;
I[this.n_x-1][this.n_x-1] = 0; // Neumann boundary condition at x=L
return math.add(I, A);
}
makeD2operator() { // create the upwind scheme second derivative operator
- const I = math.identity(this.n_x);
- 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]);
+ 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]);
I[0][0] = 0;
- I[0][1] = 1;
+ I[0][1] = -1;
return math.add(I, A, B);
}
}
@@ -199,8 +204,9 @@ class Reactor_PFR {
// const Reactor = new Reactor_PFR(200, 10, 10, 1, 100, initial_state);
// Reactor.Cs_in[0] = [0.0, 30., 100., 16., 0., 0., 5., 25., 75., 30., 0., 0., 125.];
// Reactor.Fs[0] = 10;
+// Reactor.D = 0.01;
// let N = 0;
-// while (N < 500) {
+// while (N < 5000) {
// console.log(Reactor.tick_fe(0.001));
// N += 1;
// }
From 9f1322978535e10c57eb83872d50cb416ce713bf Mon Sep 17 00:00:00 2001
From: "p.vanderwilt"
Date: Tue, 24 Jun 2025 16:38:07 +0200
Subject: [PATCH 34/97] Fix boundary conditions in gradient and second
derivative operators for Reactor_PFR
---
dependencies/reactor_class.js | 8 ++++----
1 file changed, 4 insertions(+), 4 deletions(-)
diff --git a/dependencies/reactor_class.js b/dependencies/reactor_class.js
index 9e02fca..adf744f 100644
--- a/dependencies/reactor_class.js
+++ b/dependencies/reactor_class.js
@@ -179,20 +179,20 @@ class Reactor_PFR {
}
makeDoperator() { // create the upwind scheme gradient operator
- const I = math.diag(Array(this.n_x).fill(-1), 0);
- const A = math.resize(math.diag(Array(this.n_x).fill(1), 1), [this.n_x, this.n_x]);
+ const I = math.diag(Array(this.n_x).fill(1), 0);
+ const A = math.resize(math.diag(Array(this.n_x).fill(-1), -1), [this.n_x, this.n_x]);
I[0][0] = 0;
I[0][1] = -1;
- I[this.n_x-1][this.n_x-1] = 0; // Neumann boundary condition at x=L
return math.add(I, A);
}
- makeD2operator() { // create the upwind scheme 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]);
I[0][0] = 0;
I[0][1] = -1;
+ I[0][0] = -1; // Dichelet boundary condition at outlet
return math.add(I, A, B);
}
}
From bb74fc86c212e8545a79e485146c1beaae0174d5 Mon Sep 17 00:00:00 2001
From: "p.vanderwilt"
Date: Fri, 27 Jun 2025 16:56:37 +0200
Subject: [PATCH 35/97] Refactor dispersion and boundary condition handling in
Reactor_PFR
---
dependencies/reactor_class.js | 49 +++++++++++++++++++----------------
1 file changed, 27 insertions(+), 22 deletions(-)
diff --git a/dependencies/reactor_class.js b/dependencies/reactor_class.js
index adf744f..94844ee 100644
--- a/dependencies/reactor_class.js
+++ b/dependencies/reactor_class.js
@@ -103,7 +103,7 @@ class Reactor_PFR {
this.timeStep = 1/(24*60*15); // time step [d]
this.speedUpFactor = 60;
- this.D_op = this.makeDoperator();
+ this.D_op = this.makeDoperator(false);
this.D2_op = this.makeD2operator();
}
@@ -111,6 +111,8 @@ class Reactor_PFR {
let index_in = input.payload.inlet;
this.Fs[index_in] = input.payload.F;
this.Cs_in[index_in] = input.payload.C;
+ // console.log("Pe " + this.d_x*math.sum(this.Fs)/(this.D*this.A));
+ // console.log("Co " + math.sum(this.Fs)*this.timeStep/(this.A*this.d_x));
}
set setOTR(input) { // setter for OTR (WIP) [g O2 d-1]
@@ -150,7 +152,6 @@ class Reactor_PFR {
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));
- reaction[0] = Array(13).fill(0.0);
const transfer = Array.from(Array(this.n_x), () => new Array(13).fill(0.0));
if (isNaN(this.kla)) { // calculate OTR if kla is not NaN, otherwise use externally calculated OTR
@@ -159,41 +160,45 @@ class Reactor_PFR {
transfer.forEach((x, i) => { x[0] = this.calcOTR(this.state[i][0]); });
}
- transfer[0][0] = 0;
+ const dC_total = math.multiply(math.add(dispersion, advection, reaction, transfer), time_step);
+ const new_state = math.add(this.state, dC_total);
+ // 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_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;
- const BC_dispersion = math.multiply(this.D * this.A / (math.sum(this.Fs)*this.d_x), [BC_gradient], this.state)[0];
- this.state[0] = math.add(BC_C_in, BC_dispersion);
- console.log(BC_dispersion);
+ const BC_dispersion = math.multiply(this.D * this.A / (math.sum(this.Fs)*this.d_x), [BC_gradient], new_state)[0];
+ new_state[0] = math.add(BC_C_in, BC_dispersion).map(val => val < 0 ? 0 : val);
+ console.log(new_state[0])
+ } else { // Neumann BC (no flux)
+ new_state[0] = new_state[1];
}
+ // Neumann BC (no flux)
+ new_state[this.n_x-1] = new_state[this.n_x-2]
- const dC_total = math.multiply(math.add(dispersion, advection, reaction, transfer), time_step);
-
- // clip value element-wise to each subarray to avoid negative concentrations
- this.state = math.add(this.state, dC_total).map(row => row.map(val => val < 0 ? 0 : val));
- return this.state;
+ this.state = new_state.map(row => row.map(val => val < 0 ? 0 : val)); // apply the new state
+ return new_state;
}
- makeDoperator() { // create the upwind scheme gradient operator
- const I = math.diag(Array(this.n_x).fill(1), 0);
+ makeDoperator(central=false) { // create the upwind scheme gradient operator
+ const I = math.resize(math.diag(Array(this.n_x).fill(1), central), [this.n_x, this.n_x]);
const A = math.resize(math.diag(Array(this.n_x).fill(-1), -1), [this.n_x, this.n_x]);
- I[0][0] = 0;
- I[0][1] = -1;
- return math.add(I, A);
+ 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);
+ return D;
}
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]);
- I[0][0] = 0;
- I[0][1] = -1;
- I[0][0] = -1; // Dichelet boundary condition at outlet
- return math.add(I, A, B);
+ 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);
+ return D2;
}
}
@@ -211,4 +216,4 @@ class Reactor_PFR {
// N += 1;
// }
-module.exports = {Reactor_CSTR, Reactor_PFR};
\ No newline at end of file
+module.exports = { Reactor_CSTR, Reactor_PFR };
\ No newline at end of file
From 0cc653800324a0990508ebe38799c2f6a67ddf26 Mon Sep 17 00:00:00 2001
From: "p.vanderwilt"
Date: Fri, 27 Jun 2025 17:29:20 +0200
Subject: [PATCH 36/97] Handle division by zero in rate calculations for ASM3
---
dependencies/asm3_class.js | 4 ++--
1 file changed, 2 insertions(+), 2 deletions(-)
diff --git a/dependencies/asm3_class.js b/dependencies/asm3_class.js
index 9b7762a..6f1e5fe 100644
--- a/dependencies/asm3_class.js
+++ b/dependencies/asm3_class.js
@@ -104,8 +104,8 @@ class ASM3 {
// Heterotrophs
rates[1] = k_STO * this._monod(S_O, K_O) * this._monod(S_S, K_S) * X_H;
rates[2] = k_STO * nu_NO * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * this._monod(S_S, K_S) * X_H;
- rates[3] = mu_H_max * this._monod(S_O, K_O) * this._monod(S_NH, K_NH) * this._monod(S_HCO, K_HCO) * this._monod(X_STO/X_H, K_STO) * X_H;
- rates[4] = mu_H_max * nu_NO * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * this._monod(S_NH, K_NH) * this._monod(S_HCO, K_HCO) * this._monod(X_STO/X_H, K_STO) * X_H;
+ rates[3] = X_H == 0 ? 0 : mu_H_max * this._monod(S_O, K_O) * this._monod(S_NH, K_NH) * this._monod(S_HCO, K_HCO) * this._monod(X_STO/X_H, K_STO) * X_H;
+ rates[4] = X_H == 0 ? 0 : mu_H_max * nu_NO * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * this._monod(S_NH, K_NH) * this._monod(S_HCO, K_HCO) * this._monod(X_STO/X_H, K_STO) * X_H;
rates[5] = b_H_O * this._monod(S_O, K_O) * X_H;
rates[6] = b_H_NO * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * X_H;
rates[7] = b_STO_O * this._monod(S_O, K_O) * X_H;
From 8215c5ed9a6a736eed45a47dd8bbd6dd1772e888 Mon Sep 17 00:00:00 2001
From: "p.vanderwilt"
Date: Sat, 28 Jun 2025 19:19:38 +0200
Subject: [PATCH 37/97] Add checks for NaN values in Reactor_PFR calculations
and update hydrolysis rate calculation to handle division by zero
---
dependencies/asm3_class.js | 2 +-
dependencies/reactor_class.js | 33 ++++++++++++++++++++++++++-------
2 files changed, 27 insertions(+), 8 deletions(-)
diff --git a/dependencies/asm3_class.js b/dependencies/asm3_class.js
index 6f1e5fe..f5ac1df 100644
--- a/dependencies/asm3_class.js
+++ b/dependencies/asm3_class.js
@@ -99,7 +99,7 @@ class ASM3 {
const { k_H, K_X, k_STO, nu_NO, K_O, K_NO, K_S, K_STO, mu_H_max, K_NH, K_HCO, b_H_O, b_H_NO, b_STO_O, b_STO_NO, mu_A_max, K_A_NH, K_A_O, K_A_HCO, b_A_O, b_A_NO } = this.kin_params;
// Hydrolysis
- rates[0] = k_H * this._monod(X_S / X_H, K_X) * X_H;
+ rates[0] = X_H == 0 ? 0 : k_H * this._monod(X_S / X_H, K_X) * X_H;
// Heterotrophs
rates[1] = k_STO * this._monod(S_O, K_O) * this._monod(S_S, K_S) * X_H;
diff --git a/dependencies/reactor_class.js b/dependencies/reactor_class.js
index 94844ee..8212aaf 100644
--- a/dependencies/reactor_class.js
+++ b/dependencies/reactor_class.js
@@ -103,7 +103,7 @@ class Reactor_PFR {
this.timeStep = 1/(24*60*15); // time step [d]
this.speedUpFactor = 60;
- this.D_op = this.makeDoperator(false);
+ this.D_op = this.makeDoperator(true);
this.D2_op = this.makeD2operator();
}
@@ -111,8 +111,10 @@ class Reactor_PFR {
let index_in = input.payload.inlet;
this.Fs[index_in] = input.payload.F;
this.Cs_in[index_in] = input.payload.C;
- // console.log("Pe " + this.d_x*math.sum(this.Fs)/(this.D*this.A));
- // console.log("Co " + math.sum(this.Fs)*this.timeStep/(this.A*this.d_x));
+ console.log("Pe total " + this.length*math.sum(this.Fs)/(this.D*this.A));
+ console.log("Pe local " + this.d_x*math.sum(this.Fs)/(this.D*this.A));
+ console.log("Co ad " + math.sum(this.Fs)*this.timeStep/(this.A*this.d_x));
+ console.log("Co D " + this.D*this.timeStep/(this.d_x*this.d_x));
}
set setOTR(input) { // setter for OTR (WIP) [g O2 d-1]
@@ -134,7 +136,6 @@ class Reactor_PFR {
// expect update with timestamp
updateState(newTime) {
-
const day2ms = 1000 * 60 * 60 * 24;
let n_iter = Math.floor(this.speedUpFactor*(newTime - this.currentTime) / (this.timeStep * day2ms));
@@ -153,6 +154,16 @@ class Reactor_PFR {
const advection = math.multiply(-1*math.sum(this.Fs)/(this.A*this.d_x), this.D_op, this.state);
const reaction = this.state.map((state_slice) => this.asm.compute_dC(state_slice));
const transfer = Array.from(Array(this.n_x), () => new Array(13).fill(0.0));
+
+ if (dispersion.some(row => row.some(Number.isNaN))) {
+ throw new Error("NaN detected in dispersion!");
+ }
+ if (advection.some(row => row.some(Number.isNaN))) {
+ throw new Error("NaN detected in advection!");
+ }
+ if (reaction.some(row => row.some(Number.isNaN))) {
+ throw new Error("NaN detected in reaction!");
+ }
if (isNaN(this.kla)) { // calculate OTR if kla is not NaN, otherwise use externally calculated OTR
transfer.forEach((x) => { x[0] = this.OTR; });
@@ -162,6 +173,9 @@ class Reactor_PFR {
const dC_total = math.multiply(math.add(dispersion, advection, reaction, transfer), time_step);
const new_state = math.add(this.state, dC_total);
+ if (new_state.some(row => row.some(Number.isNaN))) {
+ throw new Error("NaN detected in new_state after dC_total update!");
+ }
// apply boundary conditions
if (math.sum(this.Fs) > 0) { // Danckwerts BC
@@ -170,21 +184,26 @@ class Reactor_PFR {
BC_gradient[0] = -1;
BC_gradient[1] = 1;
const BC_dispersion = math.multiply(this.D * this.A / (math.sum(this.Fs)*this.d_x), [BC_gradient], new_state)[0];
+ console.log(math.add(BC_C_in, BC_dispersion));
new_state[0] = math.add(BC_C_in, BC_dispersion).map(val => val < 0 ? 0 : val);
- console.log(new_state[0])
+
} else { // Neumann BC (no flux)
new_state[0] = new_state[1];
}
// Neumann BC (no flux)
new_state[this.n_x-1] = new_state[this.n_x-2]
+ if (new_state.some(row => row.some(Number.isNaN))) {
+ throw new Error("NaN detected in new_state after enforcing boundary conditions!");
+ }
+
this.state = new_state.map(row => row.map(val => val < 0 ? 0 : val)); // apply the new state
return new_state;
}
makeDoperator(central=false) { // create the upwind scheme gradient operator
- const I = math.resize(math.diag(Array(this.n_x).fill(1), central), [this.n_x, this.n_x]);
- const A = math.resize(math.diag(Array(this.n_x).fill(-1), -1), [this.n_x, this.n_x]);
+ const I = math.resize(math.diag(Array(this.n_x).fill(1/(1+central)), central), [this.n_x, this.n_x]);
+ const A = math.resize(math.diag(Array(this.n_x).fill(-1/(1+central)), -1), [this.n_x, this.n_x]);
const D = math.add(I, A);
D[0] = Array(this.n_x).fill(0); // set by BCs elsewhere
D[this.n_x-1] = Array(this.n_x).fill(0);
From b2d32ba9f2011bbc878a0aa46db1925e5ed69ad4 Mon Sep 17 00:00:00 2001
From: "p.vanderwilt"
Date: Mon, 30 Jun 2025 12:50:02 +0200
Subject: [PATCH 38/97] Enhance makeDoperator to support higher-order central
gradient schemes and improve boundary handling
---
dependencies/reactor_class.js | 30 ++++++++++++++++++++++++++----
1 file changed, 26 insertions(+), 4 deletions(-)
diff --git a/dependencies/reactor_class.js b/dependencies/reactor_class.js
index 8212aaf..e80f88b 100644
--- a/dependencies/reactor_class.js
+++ b/dependencies/reactor_class.js
@@ -201,10 +201,32 @@ class Reactor_PFR {
return new_state;
}
- makeDoperator(central=false) { // create the upwind scheme gradient operator
- 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);
+ 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 D = math.add(I, A);
+ const NearBoundary = Array(this.n_x).fill(0.0);
+ NearBoundary[1] = -25/12;
+ NearBoundary[2] = 4;
+ NearBoundary[3] = -3;
+ NearBoundary[4] = 4/3;
+ NearBoundary[5] = -1/4;
+ D[1] = NearBoundary;
+ NearBoundary.reverse();
+ D[this.n_x-2] = math.multiply(-1, NearBoundary)
+ } 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 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);
return D;
From 3cc876533c64a2ac5b4dd0043af0aaf6b43e01f2 Mon Sep 17 00:00:00 2001
From: "p.vanderwilt"
Date: Mon, 30 Jun 2025 15:46:13 +0200
Subject: [PATCH 39/97] Changed the upper boundary to lower order scheme for
now
---
dependencies/reactor_class.js | 25 ++++++++++++++++---------
1 file changed, 16 insertions(+), 9 deletions(-)
diff --git a/dependencies/reactor_class.js b/dependencies/reactor_class.js
index e80f88b..0d73705 100644
--- a/dependencies/reactor_class.js
+++ b/dependencies/reactor_class.js
@@ -103,8 +103,10 @@ class Reactor_PFR {
this.timeStep = 1/(24*60*15); // time step [d]
this.speedUpFactor = 60;
- this.D_op = this.makeDoperator(true);
+ this.D_op = this.makeDoperator(true, true);
this.D2_op = this.makeD2operator();
+
+ this.alpha = 0.001; // boundary condition modifier
}
set setInfluent(input) { // setter for C_in (WIP)
@@ -183,7 +185,7 @@ class Reactor_PFR {
const BC_gradient = Array(this.n_x).fill(0.0);
BC_gradient[0] = -1;
BC_gradient[1] = 1;
- const BC_dispersion = math.multiply(this.D * this.A / (math.sum(this.Fs)*this.d_x), [BC_gradient], new_state)[0];
+ const BC_dispersion = math.multiply(this.alpha*this.D * this.A / (math.sum(this.Fs)*this.d_x), [BC_gradient], new_state)[0];
console.log(math.add(BC_C_in, BC_dispersion));
new_state[0] = math.add(BC_C_in, BC_dispersion).map(val => val < 0 ? 0 : val);
@@ -208,16 +210,22 @@ class Reactor_PFR {
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);
+ const D = math.add(I, A, B, C);
+ D[1][0] = -1;
+ D[1][1] = 0;
+ D[1][2] = 1;
+ D[1][3] = 0;
const NearBoundary = Array(this.n_x).fill(0.0);
NearBoundary[1] = -25/12;
NearBoundary[2] = 4;
NearBoundary[3] = -3;
NearBoundary[4] = 4/3;
NearBoundary[5] = -1/4;
- D[1] = NearBoundary;
NearBoundary.reverse();
- D[this.n_x-2] = math.multiply(-1, NearBoundary)
+ D[this.n_x-2] = NearBoundary;
+ D[0] = Array(this.n_x).fill(0); // set by BCs elsewhere
+ 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.");
}
@@ -225,11 +233,10 @@ class Reactor_PFR {
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);
+ return D;
}
-
- D[0] = Array(this.n_x).fill(0); // set by BCs elsewhere
- D[this.n_x-1] = Array(this.n_x).fill(0);
- return D;
}
makeD2operator() { // create the central second derivative operator
From f4824b822cac7208cd87bd5fbb12a9c182f674a8 Mon Sep 17 00:00:00 2001
From: "p.vanderwilt"
Date: Tue, 1 Jul 2025 13:04:32 +0200
Subject: [PATCH 40/97] Improved wieghted finite differencing
---
dependencies/reactor_class.js | 17 +++++++----------
1 file changed, 7 insertions(+), 10 deletions(-)
diff --git a/dependencies/reactor_class.js b/dependencies/reactor_class.js
index 0d73705..2c8d3b8 100644
--- a/dependencies/reactor_class.js
+++ b/dependencies/reactor_class.js
@@ -211,18 +211,15 @@ class Reactor_PFR {
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);
- D[1][0] = -1;
- D[1][1] = 0;
- D[1][2] = 1;
- D[1][3] = 0;
const NearBoundary = Array(this.n_x).fill(0.0);
- NearBoundary[1] = -25/12;
- NearBoundary[2] = 4;
- NearBoundary[3] = -3;
- NearBoundary[4] = 4/3;
- NearBoundary[5] = -1/4;
+ 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] = 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);
return D;
From e9847607e85d8f6776dae255337dfb04879f753d Mon Sep 17 00:00:00 2001
From: "p.vanderwilt"
Date: Tue, 1 Jul 2025 16:08:35 +0200
Subject: [PATCH 41/97] Use Generalized boundary condition by Nauman and
Mallikarjun 1983
---
dependencies/reactor_class.js | 5 +++--
1 file changed, 3 insertions(+), 2 deletions(-)
diff --git a/dependencies/reactor_class.js b/dependencies/reactor_class.js
index 2c8d3b8..6dd6168 100644
--- a/dependencies/reactor_class.js
+++ b/dependencies/reactor_class.js
@@ -2,7 +2,7 @@ const ASM3 = require('./asm3_class')
const { create, all } = require('mathjs')
const config = {
- matrix: 'Array' // Choose 'Matrix' (default) or 'Array'
+ matrix: 'Array' // choose 'Matrix' (default) or 'Array'
}
const math = create(all, config)
@@ -185,7 +185,8 @@ class Reactor_PFR {
const BC_gradient = Array(this.n_x).fill(0.0);
BC_gradient[0] = -1;
BC_gradient[1] = 1;
- const BC_dispersion = math.multiply(this.alpha*this.D * this.A / (math.sum(this.Fs)*this.d_x), [BC_gradient], new_state)[0];
+ 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, this.D * this.A / (math.sum(this.Fs)*this.d_x), [BC_gradient], new_state)[0];
console.log(math.add(BC_C_in, BC_dispersion));
new_state[0] = math.add(BC_C_in, BC_dispersion).map(val => val < 0 ? 0 : val);
From dcc8562dbf5ac2800d6ce4bc10a2b6ff8f88cc1b Mon Sep 17 00:00:00 2001
From: "p.vanderwilt"
Date: Wed, 2 Jul 2025 10:37:02 +0200
Subject: [PATCH 42/97] Remove depreciated variable
---
dependencies/reactor_class.js | 2 --
1 file changed, 2 deletions(-)
diff --git a/dependencies/reactor_class.js b/dependencies/reactor_class.js
index 6dd6168..652ee26 100644
--- a/dependencies/reactor_class.js
+++ b/dependencies/reactor_class.js
@@ -105,8 +105,6 @@ class Reactor_PFR {
this.D_op = this.makeDoperator(true, true);
this.D2_op = this.makeD2operator();
-
- this.alpha = 0.001; // boundary condition modifier
}
set setInfluent(input) { // setter for C_in (WIP)
From f517b7764d083b49f8fac052a2f81e29da0325ab Mon Sep 17 00:00:00 2001
From: "p.vanderwilt"
Date: Thu, 3 Jul 2025 22:28:34 +0200
Subject: [PATCH 43/97] Remove mistake boundary condition
---
dependencies/reactor_class.js | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/dependencies/reactor_class.js b/dependencies/reactor_class.js
index 652ee26..76a6340 100644
--- a/dependencies/reactor_class.js
+++ b/dependencies/reactor_class.js
@@ -184,7 +184,7 @@ class Reactor_PFR {
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, this.D * this.A / (math.sum(this.Fs)*this.d_x), [BC_gradient], new_state)[0];
+ 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);
From d0db1b416c36bf9b73238770dd269a6d15ac794a Mon Sep 17 00:00:00 2001
From: "p.vanderwilt"
Date: Fri, 4 Jul 2025 10:01:46 +0200
Subject: [PATCH 44/97] 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 45/97] 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 46/97] 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 47/97] 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 48/97] 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 49/97] 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 50/97] 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 51/97] 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 52/97] 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 53/97] 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 54/97] 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 55/97] 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 56/97] 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 57/97] 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 58/97] 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 59/97] 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 60/97] 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 61/97] 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.
From 01380c309f757622c5ad194f42ee938cedde2bc7 Mon Sep 17 00:00:00 2001
From: "p.vanderwilt"
Date: Mon, 7 Jul 2025 14:47:50 +0200
Subject: [PATCH 62/97] Add boundary condition input and update reactor
configuration handling
---
advanced-reactor.html | 17 +++++++++++++++++
src/nodeClass.js | 1 +
src/specificClass.js | 27 ++++++++++++++++++++++-----
3 files changed, 40 insertions(+), 5 deletions(-)
diff --git a/advanced-reactor.html b/advanced-reactor.html
index 512f3dd..6594947 100644
--- a/advanced-reactor.html
+++ b/advanced-reactor.html
@@ -8,6 +8,7 @@
volume: { value: 0., required: true },
length: { value: 0.},
resolution_L: { value: 0.},
+ boundary_condition: {value: "Generalised"},
n_inlets: { value: 1, required: true},
kla: { value: null },
S_O_init: { value: 0., required: true },
@@ -75,6 +76,18 @@
$(".PFR").show();
}
});
+ $("#node-input-boundary_condition").typedInput({
+ types: [
+ {
+ value: "Generalised",
+ options: [
+ { value: "Generalised", label: "Generalised"},
+ { value: "Diriclet", label: "Diriclet"},
+ { value: "Danckwerts", label: "Danckwerts"}
+ ]
+ }
+ ]
+ })
// Set initial visibility on dialog open
const initialType = $("#node-input-reactor_type").typedInput("value");
if (initialType === "CSTR") {
@@ -117,6 +130,10 @@
+
+
+
+
-