diff --git a/additional_nodes/recirculation-pump.js b/additional_nodes/recirculation-pump.js
index 5e6d50f..94c7018 100644
--- a/additional_nodes/recirculation-pump.js
+++ b/additional_nodes/recirculation-pump.js
@@ -1,40 +1,40 @@
-module.exports = function(RED) {
- function recirculation(config) {
- RED.nodes.createNode(this, config);
- var node = this;
-
- let name = config.name;
- let F2 = parseFloat(config.F2);
- const inlet_F2 = parseInt(config.inlet);
-
- node.on('input', function(msg, send, done) {
- switch (msg.topic) {
- case "Fluent":
- // conserve volume flow debit
- 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 = F1;
-
- let msg_F2 = {...msg};
- msg_F2.payload.F = F2_corr;
- msg_F2.payload.inlet = inlet_F2;
-
- send([msg_F1, msg_F2]);
- break;
- case "clock":
- break;
- default:
- console.log("Unknown topic: " + msg.topic);
- }
-
- if (done) {
- done();
- }
- });
-
- }
- RED.nodes.registerType("recirculation-pump", recirculation);
-};
+module.exports = function(RED) {
+ function recirculation(config) {
+ RED.nodes.createNode(this, config);
+ var node = this;
+
+ let F2 = parseFloat(config.F2);
+ const inlet_F2 = parseInt(config.inlet);
+
+ node.on('input', function(msg, send, done) {
+ switch (msg.topic) {
+ case "Fluent": {
+ // conserve volume flow debit
+ 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 = F1;
+
+ let msg_F2 = {...msg};
+ msg_F2.payload.F = F2_corr;
+ msg_F2.payload.inlet = inlet_F2;
+
+ send([msg_F1, msg_F2]);
+ break;
+ }
+ case "clock":
+ break;
+ default:
+ console.log("Unknown topic: " + msg.topic);
+ }
+
+ if (done) {
+ done();
+ }
+ });
+
+ }
+ RED.nodes.registerType("recirculation-pump", recirculation);
+};
diff --git a/additional_nodes/settling-basin.js b/additional_nodes/settling-basin.js
index 66bc57c..80affcf 100644
--- a/additional_nodes/settling-basin.js
+++ b/additional_nodes/settling-basin.js
@@ -1,57 +1,57 @@
-module.exports = function(RED) {
- function settler(config) {
- RED.nodes.createNode(this, config);
- var node = this;
-
- let name = config.name;
- let TS_set = parseFloat(config.TS_set);
- const inlet_sludge = parseInt(config.inlet);
-
- node.on('input', function(msg, send, done) {
- switch (msg.topic) {
- case "Fluent":
- // conserve volume flow debit
- let F_in = msg.payload.F;
- let C_in = msg.payload.C;
- let F2 = (F_in * C_in[12]) / TS_set;
-
- 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 = F1;
- msg_F1.payload.C[7] = 0;
- msg_F1.payload.C[8] = 0;
- msg_F1.payload.C[9] = 0;
- msg_F1.payload.C[10] = 0;
- msg_F1.payload.C[11] = 0;
- msg_F1.payload.C[12] = 0;
-
- let msg_F2 = {...msg};
- msg_F2.payload.F = F2_corr;
- if (F2_corr > 0) {
- msg_F2.payload.C[7] = F_in * C_in[7] / F2;
- msg_F2.payload.C[8] = F_in * C_in[8] / F2;
- msg_F2.payload.C[9] = F_in * C_in[9] / F2;
- msg_F2.payload.C[10] = F_in * C_in[10] / F2;
- msg_F2.payload.C[11] = F_in * C_in[11] / F2;
- msg_F2.payload.C[12] = F_in * C_in[12] / F2;
- }
- msg_F2.payload.inlet = inlet_sludge;
-
- send([msg_F1, msg_F2]);
- break;
- case "clock":
- break;
- default:
- console.log("Unknown topic: " + msg.topic);
- }
-
- if (done) {
- done();
- }
- });
-
- }
- RED.nodes.registerType("settling-basin", settler);
-};
+module.exports = function(RED) {
+ function settler(config) {
+ RED.nodes.createNode(this, config);
+ var node = this;
+
+ let TS_set = parseFloat(config.TS_set);
+ const inlet_sludge = parseInt(config.inlet);
+
+ node.on('input', function(msg, send, done) {
+ switch (msg.topic) {
+ case "Fluent": {
+ // conserve volume flow debit
+ let F_in = msg.payload.F;
+ let C_in = msg.payload.C;
+ let F2 = (F_in * C_in[12]) / TS_set;
+
+ 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 = F1;
+ msg_F1.payload.C[7] = 0;
+ msg_F1.payload.C[8] = 0;
+ msg_F1.payload.C[9] = 0;
+ msg_F1.payload.C[10] = 0;
+ msg_F1.payload.C[11] = 0;
+ msg_F1.payload.C[12] = 0;
+
+ let msg_F2 = {...msg};
+ msg_F2.payload.F = F2_corr;
+ if (F2_corr > 0) {
+ msg_F2.payload.C[7] = F_in * C_in[7] / F2;
+ msg_F2.payload.C[8] = F_in * C_in[8] / F2;
+ msg_F2.payload.C[9] = F_in * C_in[9] / F2;
+ msg_F2.payload.C[10] = F_in * C_in[10] / F2;
+ msg_F2.payload.C[11] = F_in * C_in[11] / F2;
+ msg_F2.payload.C[12] = F_in * C_in[12] / F2;
+ }
+ msg_F2.payload.inlet = inlet_sludge;
+
+ send([msg_F1, msg_F2]);
+ break;
+ }
+ case "clock":
+ break;
+ default:
+ console.log("Unknown topic: " + msg.topic);
+ }
+
+ if (done) {
+ done();
+ }
+ });
+
+ }
+ RED.nodes.registerType("settling-basin", settler);
+};
diff --git a/reactor.html b/reactor.html
index 82734f3..3ccff6b 100644
--- a/reactor.html
+++ b/reactor.html
@@ -1,267 +1,286 @@
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
diff --git a/src/nodeClass.js b/src/nodeClass.js
index a6a7701..4573cf3 100644
--- a/src/nodeClass.js
+++ b/src/nodeClass.js
@@ -1,5 +1,5 @@
const { Reactor_CSTR, Reactor_PFR } = require('./specificClass.js');
-const { outputUtils } = require('generalFunctions');
+const { outputUtils, configManager } = require('generalFunctions');
const REACTOR_SPECIES = [
'S_O',
@@ -16,23 +16,23 @@ const REACTOR_SPECIES = [
'X_A',
'X_TS'
];
-
-
-class nodeClass {
- /**
- * Node-RED node class for advanced-reactor.
- * @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
- this.node = nodeInstance;
- this.RED = RED;
- this.name = nameOfNode;
- this.source = null;
-
+
+
+class nodeClass {
+ /**
+ * Node-RED node class for advanced-reactor.
+ * @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
+ this.node = nodeInstance;
+ this.RED = RED;
+ this.name = nameOfNode;
+ this.source = null;
+
this._loadConfig(uiConfig)
this._setupClass();
this._output = new outputUtils();
@@ -40,143 +40,133 @@ class nodeClass {
this._attachInputHandler();
this._registerChild();
this._startTickLoop();
- this._attachCloseHandler();
- }
-
- /**
- * Handle node-red input messages
- */
- _attachInputHandler() {
- this.node.on('input', (msg, send, done) => {
- try {
- switch (msg.topic) {
- case "clock":
- this.source.updateState(msg.timestamp);
- send([msg, null, null]);
- break;
- case "Fluent":
- this.source.setInfluent = msg;
- break;
- case "OTR":
- this.source.setOTR = msg;
- break;
- case "Temperature":
- this.source.setTemperature = msg;
- break;
- case "Dispersion":
- this.source.setDispersion = msg;
- break;
- case 'registerChild': {
- const childId = msg.payload;
- const childObj = this.RED.nodes.getNode(childId);
- if (!childObj || !childObj.source) {
- this.source?.logger?.warn(`registerChild skipped: missing child/source for id=${childId}`);
- break;
- }
- this.source.childRegistrationUtils.registerChild(childObj.source, msg.positionVsParent);
- break;
- }
- default:
- this.source?.logger?.warn(`Unknown topic: ${msg.topic}`);
- }
- } catch (error) {
- this.source?.logger?.error(`Input handler failure: ${error.message}`);
- }
-
- if (typeof done === 'function') {
- done();
- }
- });
- }
-
- /**
- * Parse node configuration
- * @param {object} uiConfig Config set in UI in node-red
- */
- _loadConfig(uiConfig) {
- this.config = {
- general: {
- name: uiConfig.name || this.name,
- id: this.node.id,
- unit: null,
- logging: {
- enabled: uiConfig.enableLog,
- logLevel: uiConfig.logLevel
- }
- },
- functionality: {
- positionVsParent: uiConfig.positionVsParent || 'atEquipment', // Default to 'atEquipment' if not specified
- softwareType: "reactor" // should be set in config manager
- },
- reactor_type: uiConfig.reactor_type,
- volume: parseFloat(uiConfig.volume),
- length: parseFloat(uiConfig.length),
- resolution_L: parseInt(uiConfig.resolution_L),
- alpha: parseFloat(uiConfig.alpha),
- 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)
- ],
- timeStep: parseFloat(uiConfig.timeStep),
- speedUpFactor: Number(uiConfig.speedUpFactor) || 1
- }
- }
-
- /**
- * Register this node as a child upstream and downstream.
- * Delayed to avoid Node-RED startup race conditions.
- */
- _registerChild() {
- setTimeout(() => {
- this.node.send([
- null,
- null,
- { topic: 'registerChild', payload: this.node.id, positionVsParent: this.config?.functionality?.positionVsParent || 'atEquipment' }
- ]);
- }, 100);
- }
-
- /**
- * Setup reactor class based on config
- */
- _setupClass() {
- let new_reactor;
-
- switch (this.config.reactor_type) {
- case "CSTR":
- new_reactor = new Reactor_CSTR(this.config);
- break;
- case "PFR":
- new_reactor = new Reactor_PFR(this.config);
- break;
- default:
- this.node.warn("Unknown reactor type: " + this.config.reactor_type + ". Falling back to CSTR.");
- new_reactor = new Reactor_CSTR(this.config);
- }
-
- this.source = new_reactor; // protect from reassignment
- this.node.source = this.source;
- }
-
- _startTickLoop() {
- setTimeout(() => {
- this._tickInterval = setInterval(() => this._tick(), 1000);
- }, 1000);
- }
-
+ this._attachCloseHandler();
+ }
+
+ /**
+ * Handle node-red input messages
+ */
+ _attachInputHandler() {
+ this.node.on('input', (msg, send, done) => {
+ try {
+ switch (msg.topic) {
+ case "clock":
+ this.source.updateState(msg.timestamp);
+ send([msg, null, null]);
+ break;
+ case "Fluent":
+ this.source.setInfluent = msg;
+ break;
+ case "OTR":
+ this.source.setOTR = msg;
+ break;
+ case "Temperature":
+ this.source.setTemperature = msg;
+ break;
+ case "Dispersion":
+ this.source.setDispersion = msg;
+ break;
+ case 'registerChild': {
+ const childId = msg.payload;
+ const childObj = this.RED.nodes.getNode(childId);
+ if (!childObj || !childObj.source) {
+ this.source?.logger?.warn(`registerChild skipped: missing child/source for id=${childId}`);
+ break;
+ }
+ this.source.childRegistrationUtils.registerChild(childObj.source, msg.positionVsParent);
+ break;
+ }
+ default:
+ this.source?.logger?.warn(`Unknown topic: ${msg.topic}`);
+ }
+ } catch (error) {
+ this.source?.logger?.error(`Input handler failure: ${error.message}`);
+ }
+
+ if (typeof done === 'function') {
+ done();
+ }
+ });
+ }
+
+ /**
+ * Parse node configuration using ConfigManager
+ * @param {object} uiConfig Config set in UI in node-red
+ */
+ _loadConfig(uiConfig) {
+ const cfgMgr = new configManager();
+
+ // Build config: base sections + reactor-specific domain config
+ this.config = cfgMgr.buildConfig('reactor', uiConfig, this.node.id, {
+ reactor_type: uiConfig.reactor_type,
+ volume: parseFloat(uiConfig.volume),
+ length: parseFloat(uiConfig.length),
+ resolution_L: parseInt(uiConfig.resolution_L),
+ alpha: parseFloat(uiConfig.alpha),
+ 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)
+ ],
+ timeStep: parseFloat(uiConfig.timeStep),
+ speedUpFactor: Number(uiConfig.speedUpFactor) || 1
+ });
+ }
+
+ /**
+ * Register this node as a child upstream and downstream.
+ * Delayed to avoid Node-RED startup race conditions.
+ */
+ _registerChild() {
+ setTimeout(() => {
+ this.node.send([
+ null,
+ null,
+ { topic: 'registerChild', payload: this.node.id, positionVsParent: this.config?.functionality?.positionVsParent || 'atEquipment' }
+ ]);
+ }, 100);
+ }
+
+ /**
+ * Setup reactor class based on config
+ */
+ _setupClass() {
+ let new_reactor;
+
+ switch (this.config.reactor_type) {
+ case "CSTR":
+ new_reactor = new Reactor_CSTR(this.config);
+ break;
+ case "PFR":
+ new_reactor = new Reactor_PFR(this.config);
+ break;
+ default:
+ this.node.warn("Unknown reactor type: " + this.config.reactor_type + ". Falling back to CSTR.");
+ new_reactor = new Reactor_CSTR(this.config);
+ }
+
+ this.source = new_reactor; // protect from reassignment
+ this.node.source = this.source;
+ }
+
+ _startTickLoop() {
+ setTimeout(() => {
+ this._tickInterval = setInterval(() => this._tick(), 1000);
+ }, 1000);
+ }
+
_tick(){
const gridProfile = this.source.getGridProfile;
if (gridProfile) {
@@ -209,10 +199,10 @@ class nodeClass {
_attachCloseHandler() {
this.node.on('close', (done) => {
- clearInterval(this._tickInterval);
- if (typeof done === 'function') done();
- });
- }
-}
-
-module.exports = nodeClass;
+ clearInterval(this._tickInterval);
+ if (typeof done === 'function') done();
+ });
+ }
+}
+
+module.exports = nodeClass;
diff --git a/src/reaction_modules/asm3_class Koch.js b/src/reaction_modules/asm3_class Koch.js
index 11cedb2..6e4e788 100644
--- a/src/reaction_modules/asm3_class Koch.js
+++ b/src/reaction_modules/asm3_class Koch.js
@@ -1,211 +1,211 @@
-const math = require('mathjs')
-
-/**
- * ASM3 class for the Activated Sludge Model No. 3 (ASM3). Using Koch et al. 2000 parameters.
- */
-class ASM3 {
-
- constructor() {
- /**
- * Kinetic parameters for ASM3 at 20 C. Using Koch et al. 2000 parameters.
- * @property {Object} kin_params - Kinetic parameters
- */
- this.kin_params = {
- // Hydrolysis
- k_H: 9., // 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: 12., // storage rate constant [g S_S g-1 X_H d-1]
- nu_NO: 0.5, // 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: 10., // saturation constant S_s [g COD m-3]
- K_STO: 0.1, // saturation constant X_STO [g X_STO g-1 X_H]
- mu_H_max: 3., // 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.3, // aerobic respiration rate [d-1]
- b_H_NO: 0.15, // anoxic respiration rate [d-1]
- b_STO_O: 0.3, // aerobic respitation rate X_STO [d-1]
- b_STO_NO: 0.15, // anoxic respitation rate X_STO [d-1]
- // Autotrophs
- mu_A_max: 1.3, // maximum specific growth rate [d-1]
- K_A_NH: 1.4, // 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.20, // aerobic respiration rate [d-1]
- b_A_NO: 0.10 // anoxic respiration rate [d-1]
- };
-
- /**
- * Stoichiometric and composition parameters for ASM3. Using Koch et al. 2000 parameters.
- * @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
- Y_STO_O: 0.80, // aerobic yield X_STO per S_S [g X_STO g-1 S_S]
- Y_STO_NO: 0.70, // anoxic yield X_STO per S_S [g X_STO g-1 S_S]
- Y_H_O: 0.80, // aerobic yield X_H per X_STO [g X_H g-1 X_STO]
- Y_H_NO: 0.65, // 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.04, // nitrogen content X_I [g N g-1 X_I]
- i_NXS: 0.03, // 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]
- };
-
- /**
- * Temperature theta parameters for ASM3. Using Koch et al. 2000 parameters.
- * These parameters are used to adjust reaction rates based on temperature.
- * @property {Object} temp_params - Temperature theta parameters
- */
- this.temp_params = {
- // Hydrolysis
- theta_H: 0.04,
- // Heterotrophs
- theta_STO: 0.07,
- theta_mu_H: 0.07,
- theta_b_H_O: 0.07,
- theta_b_H_NO: 0.07,
- theta_b_STO_O: this._compute_theta(0.1, 0.3, 10, 20),
- theta_b_STO_NO: this._compute_theta(0.05, 0.15, 10, 20),
- // Autotrophs
- theta_mu_A: 0.105,
- theta_b_A_O: 0.105,
- theta_b_A_NO: 0.105
- };
-
- 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;
-
- const stoi_matrix = Array(12);
- // 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
- stoi_matrix[0] = [0., f_SI, 1.-f_SI, i_NXS-(1.-f_SI)*i_NSS-f_SI*i_NSI, 0., 0., (i_NXS-(1.-f_SI)*i_NSS-f_SI*i_NSI)*i_cNH, 0., -1., 0., 0., 0., -i_TSXS];
- stoi_matrix[1] = [-(1.-Y_STO_O), 0, -1., i_NSS, 0., 0., i_NSS*i_cNH, 0., 0., 0., Y_STO_O, 0., Y_STO_O*i_TSSTO];
- stoi_matrix[2] = [0., 0., -1., i_NSS, -(1.-Y_STO_NO)/(i_CODNO-i_CODN), (1.-Y_STO_NO)/(i_CODNO-i_CODN), i_NSS*i_cNH + (1.-Y_STO_NO)/(i_CODNO-i_CODN)*i_cNO, 0., 0., 0., Y_STO_NO, 0., Y_STO_NO*i_TSSTO];
- stoi_matrix[3] = [-(1.-Y_H_O)/Y_H_O, 0., 0., -i_NBM, 0., 0., -i_NBM*i_cNH, 0., 0., 1., -1./Y_H_O, 0., i_TSBM-i_TSSTO/Y_H_O];
- stoi_matrix[4] = [0., 0., 0., -i_NBM, -(1.-Y_H_NO)/(Y_H_NO*(i_CODNO-i_CODN)), (1.-Y_H_NO)/(Y_H_NO*(i_CODNO-i_CODN)), -i_NBM*i_cNH+(1.-Y_H_NO)/(Y_H_NO*(i_CODNO-i_CODN))*i_cNO, 0., 0., 1., -1./Y_H_NO, 0., i_TSBM-i_TSSTO/Y_H_NO];
- stoi_matrix[5] = [f_XI-1., 0., 0., i_NBM-f_XI*i_NXI, 0., 0., (i_NBM-f_XI*i_NXI)*i_cNH, f_XI, 0., -1., 0., 0., f_XI*i_TSXI-i_TSBM];
- stoi_matrix[6] = [0., 0., 0., i_NBM-f_XI*i_NXI, -(1.-f_XI)/(i_CODNO-i_CODN), (1.-f_XI)/(i_CODNO-i_CODN), (i_NBM-f_XI*i_NXI)*i_cNH+(1-f_XI)/(i_CODNO-i_CODN)*i_cNO, f_XI, 0., -1., 0., 0., f_XI*i_TSXI-i_TSBM];
- stoi_matrix[7] = [-1., 0., 0., 0., 0., 0., 0., 0., 0., 0., -1., 0., -i_TSSTO];
- stoi_matrix[8] = [0., 0., 0., 0., -1./(i_CODNO-i_CODN), 1./(i_CODNO-i_CODN), i_cNO/(i_CODNO-i_CODN), 0., 0., 0., -1., 0., -i_TSSTO];
- stoi_matrix[9] = [1.+i_CODNO/Y_A, 0., 0., -1./Y_A-i_NBM, 0., 1./Y_A, (-1./Y_A-i_NBM)*i_cNH+i_cNO/Y_A, 0., 0., 0., 0., 1., i_TSBM];
- stoi_matrix[10] = [f_XI-1., 0., 0., i_NBM-f_XI*i_NXI, 0., 0., (i_NBM-f_XI*i_NXI)*i_cNH, f_XI, 0., 0., 0., -1., f_XI*i_TSXI-i_TSBM];
- stoi_matrix[11] = [0., 0., 0., i_NBM-f_XI*i_NXI, -(1.-f_XI)/(i_CODNO-i_CODN), (1.-f_XI)/(i_CODNO-i_CODN), (i_NBM-f_XI*i_NXI)*i_cNH+(1-f_XI)/(i_CODNO-i_CODN)*i_cNO, 0., 0., 0., 0., -1., f_XI*i_TSXI-i_TSBM];
-
- 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);
- }
-
- /**
- * Adjust the rate parameter for temperature T using simplied Arrhenius equation based on rate constant at 20 degrees Celsius and theta parameter.
- * @param {number} k - Rate constant at 20 degrees Celcius.
- * @param {number} theta - Theta parameter.
- * @param {number} T - Temperature in Celcius.
- * @returns {number} - Adjusted rate parameter at temperature T based on the Arrhenius equation.
- */
- _arrhenius(k, theta, T) {
- return k * Math.exp(theta*(T-20));
- }
-
- /**
- * Computes the temperature theta parameter based on two rate constants and their corresponding temperatures.
- * @param {number} k1 - Rate constant at temperature T1.
- * @param {number} k2 - Rate constant at temperature T2.
- * @param {number} T1 - Temperature T1 in Celcius.
- * @param {number} T2 - Temperature T2 in Celcius.
- * @returns {number} - Theta parameter.
- */
- _compute_theta(k1, k2, T1, T2) {
- return Math.log(k1/k2)/(T1-T2);
- }
-
- /**
- * Computes the reaction rates for each process reaction based on the current state and temperature.
- * @param {Array} state - State vector containing concentrations of reaction species.
- * @param {number} [T=20] - Temperature in degrees Celsius (default is 20).
- * @returns {Array} - Reaction rates for each process reaction.
- */
- compute_rates(state, T = 20) {
- // 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;
- 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;
- const { theta_H, theta_STO, theta_mu_H, theta_b_H_O, theta_b_H_NO, theta_b_STO_O, theta_b_STO_NO, theta_mu_A, theta_b_A_O, theta_b_A_NO } = this.temp_params;
-
- // Hydrolysis
- rates[0] = X_H == 0 ? 0 : this._arrhenius(k_H, theta_H, T) * this._monod(X_S / X_H, K_X) * X_H;
-
- // Heterotrophs
- rates[1] = this._arrhenius(k_STO, theta_STO, T) * this._monod(S_O, K_O) * this._monod(S_S, K_S) * X_H;
- rates[2] = this._arrhenius(k_STO, theta_STO, T) * 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] = X_H == 0 ? 0 : this._arrhenius(mu_H_max, theta_mu_H, T) * 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 : this._arrhenius(mu_H_max, theta_mu_H, T) * 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] = this._arrhenius(b_H_O, theta_b_H_O, T) * this._monod(S_O, K_O) * X_H;
- rates[6] = this._arrhenius(b_H_NO, theta_b_H_NO, T) * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * X_H;
- rates[7] = this._arrhenius(b_STO_O, theta_b_STO_O, T) * this._monod(S_O, K_O) * X_H;
- rates[8] = this._arrhenius(b_STO_NO, theta_b_STO_NO, T) * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * X_STO;
-
- // Autotrophs
- rates[9] = this._arrhenius(mu_A_max, theta_mu_A, T) * this._monod(S_O, K_A_O) * this._monod(S_NH, K_A_NH) * this._monod(S_HCO, K_A_HCO) * X_A;
- rates[10] = this._arrhenius(b_A_O, theta_b_A_O, T) * this._monod(S_O, K_O) * X_A;
- rates[11] = this._arrhenius(b_A_NO, theta_b_A_NO, T) * this._inv_monod(S_O, K_A_O) * this._monod(S_NO, K_NO) * X_A;
-
- return rates;
- }
-
- /**
- * Computes the change in concentrations of reaction species based on the current state and temperature.
- * @param {Array} state - State vector containing concentrations of reaction species.
- * @param {number} [T=20] - Temperature in degrees Celsius (default is 20).
- * @returns {Array} - Change in reaction species concentrations.
- */
- compute_dC(state, T = 20) { // 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, T));
- }
-}
-
-module.exports = ASM3;
\ No newline at end of file
+const math = require('mathjs')
+
+/**
+ * ASM3 class for the Activated Sludge Model No. 3 (ASM3). Using Koch et al. 2000 parameters.
+ */
+class ASM3 {
+
+ constructor() {
+ /**
+ * Kinetic parameters for ASM3 at 20 C. Using Koch et al. 2000 parameters.
+ * @property {Object} kin_params - Kinetic parameters
+ */
+ this.kin_params = {
+ // Hydrolysis
+ k_H: 9., // 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: 12., // storage rate constant [g S_S g-1 X_H d-1]
+ nu_NO: 0.5, // 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: 10., // saturation constant S_s [g COD m-3]
+ K_STO: 0.1, // saturation constant X_STO [g X_STO g-1 X_H]
+ mu_H_max: 3., // 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.3, // aerobic respiration rate [d-1]
+ b_H_NO: 0.15, // anoxic respiration rate [d-1]
+ b_STO_O: 0.3, // aerobic respitation rate X_STO [d-1]
+ b_STO_NO: 0.15, // anoxic respitation rate X_STO [d-1]
+ // Autotrophs
+ mu_A_max: 1.3, // maximum specific growth rate [d-1]
+ K_A_NH: 1.4, // 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.20, // aerobic respiration rate [d-1]
+ b_A_NO: 0.10 // anoxic respiration rate [d-1]
+ };
+
+ /**
+ * Stoichiometric and composition parameters for ASM3. Using Koch et al. 2000 parameters.
+ * @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
+ Y_STO_O: 0.80, // aerobic yield X_STO per S_S [g X_STO g-1 S_S]
+ Y_STO_NO: 0.70, // anoxic yield X_STO per S_S [g X_STO g-1 S_S]
+ Y_H_O: 0.80, // aerobic yield X_H per X_STO [g X_H g-1 X_STO]
+ Y_H_NO: 0.65, // 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.04, // nitrogen content X_I [g N g-1 X_I]
+ i_NXS: 0.03, // 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]
+ };
+
+ /**
+ * Temperature theta parameters for ASM3. Using Koch et al. 2000 parameters.
+ * These parameters are used to adjust reaction rates based on temperature.
+ * @property {Object} temp_params - Temperature theta parameters
+ */
+ this.temp_params = {
+ // Hydrolysis
+ theta_H: 0.04,
+ // Heterotrophs
+ theta_STO: 0.07,
+ theta_mu_H: 0.07,
+ theta_b_H_O: 0.07,
+ theta_b_H_NO: 0.07,
+ theta_b_STO_O: this._compute_theta(0.1, 0.3, 10, 20),
+ theta_b_STO_NO: this._compute_theta(0.05, 0.15, 10, 20),
+ // Autotrophs
+ theta_mu_A: 0.105,
+ theta_b_A_O: 0.105,
+ theta_b_A_NO: 0.105
+ };
+
+ 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;
+
+ const stoi_matrix = Array(12);
+ // 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
+ stoi_matrix[0] = [0., f_SI, 1.-f_SI, i_NXS-(1.-f_SI)*i_NSS-f_SI*i_NSI, 0., 0., (i_NXS-(1.-f_SI)*i_NSS-f_SI*i_NSI)*i_cNH, 0., -1., 0., 0., 0., -i_TSXS];
+ stoi_matrix[1] = [-(1.-Y_STO_O), 0, -1., i_NSS, 0., 0., i_NSS*i_cNH, 0., 0., 0., Y_STO_O, 0., Y_STO_O*i_TSSTO];
+ stoi_matrix[2] = [0., 0., -1., i_NSS, -(1.-Y_STO_NO)/(i_CODNO-i_CODN), (1.-Y_STO_NO)/(i_CODNO-i_CODN), i_NSS*i_cNH + (1.-Y_STO_NO)/(i_CODNO-i_CODN)*i_cNO, 0., 0., 0., Y_STO_NO, 0., Y_STO_NO*i_TSSTO];
+ stoi_matrix[3] = [-(1.-Y_H_O)/Y_H_O, 0., 0., -i_NBM, 0., 0., -i_NBM*i_cNH, 0., 0., 1., -1./Y_H_O, 0., i_TSBM-i_TSSTO/Y_H_O];
+ stoi_matrix[4] = [0., 0., 0., -i_NBM, -(1.-Y_H_NO)/(Y_H_NO*(i_CODNO-i_CODN)), (1.-Y_H_NO)/(Y_H_NO*(i_CODNO-i_CODN)), -i_NBM*i_cNH+(1.-Y_H_NO)/(Y_H_NO*(i_CODNO-i_CODN))*i_cNO, 0., 0., 1., -1./Y_H_NO, 0., i_TSBM-i_TSSTO/Y_H_NO];
+ stoi_matrix[5] = [f_XI-1., 0., 0., i_NBM-f_XI*i_NXI, 0., 0., (i_NBM-f_XI*i_NXI)*i_cNH, f_XI, 0., -1., 0., 0., f_XI*i_TSXI-i_TSBM];
+ stoi_matrix[6] = [0., 0., 0., i_NBM-f_XI*i_NXI, -(1.-f_XI)/(i_CODNO-i_CODN), (1.-f_XI)/(i_CODNO-i_CODN), (i_NBM-f_XI*i_NXI)*i_cNH+(1-f_XI)/(i_CODNO-i_CODN)*i_cNO, f_XI, 0., -1., 0., 0., f_XI*i_TSXI-i_TSBM];
+ stoi_matrix[7] = [-1., 0., 0., 0., 0., 0., 0., 0., 0., 0., -1., 0., -i_TSSTO];
+ stoi_matrix[8] = [0., 0., 0., 0., -1./(i_CODNO-i_CODN), 1./(i_CODNO-i_CODN), i_cNO/(i_CODNO-i_CODN), 0., 0., 0., -1., 0., -i_TSSTO];
+ stoi_matrix[9] = [1.+i_CODNO/Y_A, 0., 0., -1./Y_A-i_NBM, 0., 1./Y_A, (-1./Y_A-i_NBM)*i_cNH+i_cNO/Y_A, 0., 0., 0., 0., 1., i_TSBM];
+ stoi_matrix[10] = [f_XI-1., 0., 0., i_NBM-f_XI*i_NXI, 0., 0., (i_NBM-f_XI*i_NXI)*i_cNH, f_XI, 0., 0., 0., -1., f_XI*i_TSXI-i_TSBM];
+ stoi_matrix[11] = [0., 0., 0., i_NBM-f_XI*i_NXI, -(1.-f_XI)/(i_CODNO-i_CODN), (1.-f_XI)/(i_CODNO-i_CODN), (i_NBM-f_XI*i_NXI)*i_cNH+(1-f_XI)/(i_CODNO-i_CODN)*i_cNO, 0., 0., 0., 0., -1., f_XI*i_TSXI-i_TSBM];
+
+ 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);
+ }
+
+ /**
+ * Adjust the rate parameter for temperature T using simplied Arrhenius equation based on rate constant at 20 degrees Celsius and theta parameter.
+ * @param {number} k - Rate constant at 20 degrees Celcius.
+ * @param {number} theta - Theta parameter.
+ * @param {number} T - Temperature in Celcius.
+ * @returns {number} - Adjusted rate parameter at temperature T based on the Arrhenius equation.
+ */
+ _arrhenius(k, theta, T) {
+ return k * Math.exp(theta*(T-20));
+ }
+
+ /**
+ * Computes the temperature theta parameter based on two rate constants and their corresponding temperatures.
+ * @param {number} k1 - Rate constant at temperature T1.
+ * @param {number} k2 - Rate constant at temperature T2.
+ * @param {number} T1 - Temperature T1 in Celcius.
+ * @param {number} T2 - Temperature T2 in Celcius.
+ * @returns {number} - Theta parameter.
+ */
+ _compute_theta(k1, k2, T1, T2) {
+ return Math.log(k1/k2)/(T1-T2);
+ }
+
+ /**
+ * Computes the reaction rates for each process reaction based on the current state and temperature.
+ * @param {Array} state - State vector containing concentrations of reaction species.
+ * @param {number} [T=20] - Temperature in degrees Celsius (default is 20).
+ * @returns {Array} - Reaction rates for each process reaction.
+ */
+ compute_rates(state, T = 20) {
+ // 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_S, S_NH, , S_NO, S_HCO, , X_S, X_H, X_STO, X_A] = state;
+ 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;
+ const { theta_H, theta_STO, theta_mu_H, theta_b_H_O, theta_b_H_NO, theta_b_STO_O, theta_b_STO_NO, theta_mu_A, theta_b_A_O, theta_b_A_NO } = this.temp_params;
+
+ // Hydrolysis
+ rates[0] = X_H == 0 ? 0 : this._arrhenius(k_H, theta_H, T) * this._monod(X_S / X_H, K_X) * X_H;
+
+ // Heterotrophs
+ rates[1] = this._arrhenius(k_STO, theta_STO, T) * this._monod(S_O, K_O) * this._monod(S_S, K_S) * X_H;
+ rates[2] = this._arrhenius(k_STO, theta_STO, T) * 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] = X_H == 0 ? 0 : this._arrhenius(mu_H_max, theta_mu_H, T) * 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 : this._arrhenius(mu_H_max, theta_mu_H, T) * 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] = this._arrhenius(b_H_O, theta_b_H_O, T) * this._monod(S_O, K_O) * X_H;
+ rates[6] = this._arrhenius(b_H_NO, theta_b_H_NO, T) * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * X_H;
+ rates[7] = this._arrhenius(b_STO_O, theta_b_STO_O, T) * this._monod(S_O, K_O) * X_H;
+ rates[8] = this._arrhenius(b_STO_NO, theta_b_STO_NO, T) * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * X_STO;
+
+ // Autotrophs
+ rates[9] = this._arrhenius(mu_A_max, theta_mu_A, T) * this._monod(S_O, K_A_O) * this._monod(S_NH, K_A_NH) * this._monod(S_HCO, K_A_HCO) * X_A;
+ rates[10] = this._arrhenius(b_A_O, theta_b_A_O, T) * this._monod(S_O, K_O) * X_A;
+ rates[11] = this._arrhenius(b_A_NO, theta_b_A_NO, T) * this._inv_monod(S_O, K_A_O) * this._monod(S_NO, K_NO) * X_A;
+
+ return rates;
+ }
+
+ /**
+ * Computes the change in concentrations of reaction species based on the current state and temperature.
+ * @param {Array} state - State vector containing concentrations of reaction species.
+ * @param {number} [T=20] - Temperature in degrees Celsius (default is 20).
+ * @returns {Array} - Change in reaction species concentrations.
+ */
+ compute_dC(state, T = 20) { // 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, T));
+ }
+}
+
+module.exports = ASM3;
diff --git a/src/reaction_modules/asm3_class.js b/src/reaction_modules/asm3_class.js
index c10c00c..b2bc3bf 100644
--- a/src/reaction_modules/asm3_class.js
+++ b/src/reaction_modules/asm3_class.js
@@ -1,211 +1,211 @@
-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 = {
- // 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]
- };
-
- /**
- * 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
- 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]
- };
-
- /**
- * Temperature theta parameters for ASM3.
- * These parameters are used to adjust reaction rates based on temperature.
- * @property {Object} temp_params - Temperature theta parameters
- */
- this.temp_params = {
- // Hydrolysis
- theta_H: this._compute_theta(2, 3, 10, 20),
- // Heterotrophs
- theta_STO: this._compute_theta(2.5, 5, 10, 20),
- theta_mu_H: this._compute_theta(1, 2, 10, 20),
- theta_b_H_O: this._compute_theta(0.1, 0.2, 10, 20),
- theta_b_H_NO: this._compute_theta(0.05, 0.1, 10, 20),
- theta_b_STO_O: this._compute_theta(0.1, 0.2, 10, 20),
- theta_b_STO_NO: this._compute_theta(0.05, 0.1, 10, 20),
- // Autotrophs
- theta_mu_A: this._compute_theta(0.35, 1, 10, 20),
- theta_b_A_O: this._compute_theta(0.05, 0.15, 10, 20),
- theta_b_A_NO: this._compute_theta(0.02, 0.05, 10, 20)
- };
-
- 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;
-
- const stoi_matrix = Array(12);
- // 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
- stoi_matrix[0] = [0., f_SI, 1.-f_SI, i_NXS-(1.-f_SI)*i_NSS-f_SI*i_NSI, 0., 0., (i_NXS-(1.-f_SI)*i_NSS-f_SI*i_NSI)*i_cNH, 0., -1., 0., 0., 0., -i_TSXS];
- stoi_matrix[1] = [-(1.-Y_STO_O), 0, -1., i_NSS, 0., 0., i_NSS*i_cNH, 0., 0., 0., Y_STO_O, 0., Y_STO_O*i_TSSTO];
- stoi_matrix[2] = [0., 0., -1., i_NSS, -(1.-Y_STO_NO)/(i_CODNO-i_CODN), (1.-Y_STO_NO)/(i_CODNO-i_CODN), i_NSS*i_cNH + (1.-Y_STO_NO)/(i_CODNO-i_CODN)*i_cNO, 0., 0., 0., Y_STO_NO, 0., Y_STO_NO*i_TSSTO];
- stoi_matrix[3] = [-(1.-Y_H_O)/Y_H_O, 0., 0., -i_NBM, 0., 0., -i_NBM*i_cNH, 0., 0., 1., -1./Y_H_O, 0., i_TSBM-i_TSSTO/Y_H_O];
- stoi_matrix[4] = [0., 0., 0., -i_NBM, -(1.-Y_H_NO)/(Y_H_NO*(i_CODNO-i_CODN)), (1.-Y_H_NO)/(Y_H_NO*(i_CODNO-i_CODN)), -i_NBM*i_cNH+(1.-Y_H_NO)/(Y_H_NO*(i_CODNO-i_CODN))*i_cNO, 0., 0., 1., -1./Y_H_NO, 0., i_TSBM-i_TSSTO/Y_H_NO];
- stoi_matrix[5] = [f_XI-1., 0., 0., i_NBM-f_XI*i_NXI, 0., 0., (i_NBM-f_XI*i_NXI)*i_cNH, f_XI, 0., -1., 0., 0., f_XI*i_TSXI-i_TSBM];
- stoi_matrix[6] = [0., 0., 0., i_NBM-f_XI*i_NXI, -(1.-f_XI)/(i_CODNO-i_CODN), (1.-f_XI)/(i_CODNO-i_CODN), (i_NBM-f_XI*i_NXI)*i_cNH+(1-f_XI)/(i_CODNO-i_CODN)*i_cNO, f_XI, 0., -1., 0., 0., f_XI*i_TSXI-i_TSBM];
- stoi_matrix[7] = [-1., 0., 0., 0., 0., 0., 0., 0., 0., 0., -1., 0., -i_TSSTO];
- stoi_matrix[8] = [0., 0., 0., 0., -1./(i_CODNO-i_CODN), 1./(i_CODNO-i_CODN), i_cNO/(i_CODNO-i_CODN), 0., 0., 0., -1., 0., -i_TSSTO];
- stoi_matrix[9] = [1.+i_CODNO/Y_A, 0., 0., -1./Y_A-i_NBM, 0., 1./Y_A, (-1./Y_A-i_NBM)*i_cNH+i_cNO/Y_A, 0., 0., 0., 0., 1., i_TSBM];
- stoi_matrix[10] = [f_XI-1., 0., 0., i_NBM-f_XI*i_NXI, 0., 0., (i_NBM-f_XI*i_NXI)*i_cNH, f_XI, 0., 0., 0., -1., f_XI*i_TSXI-i_TSBM];
- stoi_matrix[11] = [0., 0., 0., i_NBM-f_XI*i_NXI, -(1.-f_XI)/(i_CODNO-i_CODN), (1.-f_XI)/(i_CODNO-i_CODN), (i_NBM-f_XI*i_NXI)*i_cNH+(1-f_XI)/(i_CODNO-i_CODN)*i_cNO, 0., 0., 0., 0., -1., f_XI*i_TSXI-i_TSBM];
-
- 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);
- }
-
- /**
- * Adjust the rate parameter for temperature T using simplied Arrhenius equation based on rate constant at 20 degrees Celsius and theta parameter.
- * @param {number} k - Rate constant at 20 degrees Celcius.
- * @param {number} theta - Theta parameter.
- * @param {number} T - Temperature in Celcius.
- * @returns {number} - Adjusted rate parameter at temperature T based on the Arrhenius equation.
- */
- _arrhenius(k, theta, T) {
- return k * Math.exp(theta*(T-20));
- }
-
- /**
- * Computes the temperature theta parameter based on two rate constants and their corresponding temperatures.
- * @param {number} k1 - Rate constant at temperature T1.
- * @param {number} k2 - Rate constant at temperature T2.
- * @param {number} T1 - Temperature T1 in Celcius.
- * @param {number} T2 - Temperature T2 in Celcius.
- * @returns {number} - Theta parameter.
- */
- _compute_theta(k1, k2, T1, T2) {
- return Math.log(k1/k2)/(T1-T2);
- }
-
- /**
- * Computes the reaction rates for each process reaction based on the current state and temperature.
- * @param {Array} state - State vector containing concentrations of reaction species.
- * @param {number} [T=20] - Temperature in degrees Celsius (default is 20).
- * @returns {Array} - Reaction rates for each process reaction.
- */
- compute_rates(state, T = 20) {
- // 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;
- 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;
- const { theta_H, theta_STO, theta_mu_H, theta_b_H_O, theta_b_H_NO, theta_b_STO_O, theta_b_STO_NO, theta_mu_A, theta_b_A_O, theta_b_A_NO } = this.temp_params;
-
- // Hydrolysis
- rates[0] = X_H == 0 ? 0 : this._arrhenius(k_H, theta_H, T) * this._monod(X_S / X_H, K_X) * X_H;
-
- // Heterotrophs
- rates[1] = this._arrhenius(k_STO, theta_STO, T) * this._monod(S_O, K_O) * this._monod(S_S, K_S) * X_H;
- rates[2] = this._arrhenius(k_STO, theta_STO, T) * 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] = X_H == 0 ? 0 : this._arrhenius(mu_H_max, theta_mu_H, T) * 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 : this._arrhenius(mu_H_max, theta_mu_H, T) * 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] = this._arrhenius(b_H_O, theta_b_H_O, T) * this._monod(S_O, K_O) * X_H;
- rates[6] = this._arrhenius(b_H_NO, theta_b_H_NO, T) * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * X_H;
- rates[7] = this._arrhenius(b_STO_O, theta_b_STO_O, T) * this._monod(S_O, K_O) * X_H;
- rates[8] = this._arrhenius(b_STO_NO, theta_b_STO_NO, T) * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * X_STO;
-
- // Autotrophs
- rates[9] = this._arrhenius(mu_A_max, theta_mu_A, T) * this._monod(S_O, K_A_O) * this._monod(S_NH, K_A_NH) * this._monod(S_HCO, K_A_HCO) * X_A;
- rates[10] = this._arrhenius(b_A_O, theta_b_A_O, T) * this._monod(S_O, K_O) * X_A;
- rates[11] = this._arrhenius(b_A_NO, theta_b_A_NO, T) * this._inv_monod(S_O, K_A_O) * this._monod(S_NO, K_NO) * X_A;
-
- return rates;
- }
-
- /**
- * Computes the change in concentrations of reaction species based on the current state and temperature.
- * @param {Array} state - State vector containing concentrations of reaction species.
- * @param {number} [T=20] - Temperature in degrees Celsius (default is 20).
- * @returns {Array} - Change in reaction species concentrations.
- */
- compute_dC(state, T = 20) { // 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, T));
- }
-}
-
-module.exports = ASM3;
\ No newline at end of file
+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 = {
+ // 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]
+ };
+
+ /**
+ * 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
+ 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]
+ };
+
+ /**
+ * Temperature theta parameters for ASM3.
+ * These parameters are used to adjust reaction rates based on temperature.
+ * @property {Object} temp_params - Temperature theta parameters
+ */
+ this.temp_params = {
+ // Hydrolysis
+ theta_H: this._compute_theta(2, 3, 10, 20),
+ // Heterotrophs
+ theta_STO: this._compute_theta(2.5, 5, 10, 20),
+ theta_mu_H: this._compute_theta(1, 2, 10, 20),
+ theta_b_H_O: this._compute_theta(0.1, 0.2, 10, 20),
+ theta_b_H_NO: this._compute_theta(0.05, 0.1, 10, 20),
+ theta_b_STO_O: this._compute_theta(0.1, 0.2, 10, 20),
+ theta_b_STO_NO: this._compute_theta(0.05, 0.1, 10, 20),
+ // Autotrophs
+ theta_mu_A: this._compute_theta(0.35, 1, 10, 20),
+ theta_b_A_O: this._compute_theta(0.05, 0.15, 10, 20),
+ theta_b_A_NO: this._compute_theta(0.02, 0.05, 10, 20)
+ };
+
+ 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;
+
+ const stoi_matrix = Array(12);
+ // 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
+ stoi_matrix[0] = [0., f_SI, 1.-f_SI, i_NXS-(1.-f_SI)*i_NSS-f_SI*i_NSI, 0., 0., (i_NXS-(1.-f_SI)*i_NSS-f_SI*i_NSI)*i_cNH, 0., -1., 0., 0., 0., -i_TSXS];
+ stoi_matrix[1] = [-(1.-Y_STO_O), 0, -1., i_NSS, 0., 0., i_NSS*i_cNH, 0., 0., 0., Y_STO_O, 0., Y_STO_O*i_TSSTO];
+ stoi_matrix[2] = [0., 0., -1., i_NSS, -(1.-Y_STO_NO)/(i_CODNO-i_CODN), (1.-Y_STO_NO)/(i_CODNO-i_CODN), i_NSS*i_cNH + (1.-Y_STO_NO)/(i_CODNO-i_CODN)*i_cNO, 0., 0., 0., Y_STO_NO, 0., Y_STO_NO*i_TSSTO];
+ stoi_matrix[3] = [-(1.-Y_H_O)/Y_H_O, 0., 0., -i_NBM, 0., 0., -i_NBM*i_cNH, 0., 0., 1., -1./Y_H_O, 0., i_TSBM-i_TSSTO/Y_H_O];
+ stoi_matrix[4] = [0., 0., 0., -i_NBM, -(1.-Y_H_NO)/(Y_H_NO*(i_CODNO-i_CODN)), (1.-Y_H_NO)/(Y_H_NO*(i_CODNO-i_CODN)), -i_NBM*i_cNH+(1.-Y_H_NO)/(Y_H_NO*(i_CODNO-i_CODN))*i_cNO, 0., 0., 1., -1./Y_H_NO, 0., i_TSBM-i_TSSTO/Y_H_NO];
+ stoi_matrix[5] = [f_XI-1., 0., 0., i_NBM-f_XI*i_NXI, 0., 0., (i_NBM-f_XI*i_NXI)*i_cNH, f_XI, 0., -1., 0., 0., f_XI*i_TSXI-i_TSBM];
+ stoi_matrix[6] = [0., 0., 0., i_NBM-f_XI*i_NXI, -(1.-f_XI)/(i_CODNO-i_CODN), (1.-f_XI)/(i_CODNO-i_CODN), (i_NBM-f_XI*i_NXI)*i_cNH+(1-f_XI)/(i_CODNO-i_CODN)*i_cNO, f_XI, 0., -1., 0., 0., f_XI*i_TSXI-i_TSBM];
+ stoi_matrix[7] = [-1., 0., 0., 0., 0., 0., 0., 0., 0., 0., -1., 0., -i_TSSTO];
+ stoi_matrix[8] = [0., 0., 0., 0., -1./(i_CODNO-i_CODN), 1./(i_CODNO-i_CODN), i_cNO/(i_CODNO-i_CODN), 0., 0., 0., -1., 0., -i_TSSTO];
+ stoi_matrix[9] = [1.+i_CODNO/Y_A, 0., 0., -1./Y_A-i_NBM, 0., 1./Y_A, (-1./Y_A-i_NBM)*i_cNH+i_cNO/Y_A, 0., 0., 0., 0., 1., i_TSBM];
+ stoi_matrix[10] = [f_XI-1., 0., 0., i_NBM-f_XI*i_NXI, 0., 0., (i_NBM-f_XI*i_NXI)*i_cNH, f_XI, 0., 0., 0., -1., f_XI*i_TSXI-i_TSBM];
+ stoi_matrix[11] = [0., 0., 0., i_NBM-f_XI*i_NXI, -(1.-f_XI)/(i_CODNO-i_CODN), (1.-f_XI)/(i_CODNO-i_CODN), (i_NBM-f_XI*i_NXI)*i_cNH+(1-f_XI)/(i_CODNO-i_CODN)*i_cNO, 0., 0., 0., 0., -1., f_XI*i_TSXI-i_TSBM];
+
+ 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);
+ }
+
+ /**
+ * Adjust the rate parameter for temperature T using simplied Arrhenius equation based on rate constant at 20 degrees Celsius and theta parameter.
+ * @param {number} k - Rate constant at 20 degrees Celcius.
+ * @param {number} theta - Theta parameter.
+ * @param {number} T - Temperature in Celcius.
+ * @returns {number} - Adjusted rate parameter at temperature T based on the Arrhenius equation.
+ */
+ _arrhenius(k, theta, T) {
+ return k * Math.exp(theta*(T-20));
+ }
+
+ /**
+ * Computes the temperature theta parameter based on two rate constants and their corresponding temperatures.
+ * @param {number} k1 - Rate constant at temperature T1.
+ * @param {number} k2 - Rate constant at temperature T2.
+ * @param {number} T1 - Temperature T1 in Celcius.
+ * @param {number} T2 - Temperature T2 in Celcius.
+ * @returns {number} - Theta parameter.
+ */
+ _compute_theta(k1, k2, T1, T2) {
+ return Math.log(k1/k2)/(T1-T2);
+ }
+
+ /**
+ * Computes the reaction rates for each process reaction based on the current state and temperature.
+ * @param {Array} state - State vector containing concentrations of reaction species.
+ * @param {number} [T=20] - Temperature in degrees Celsius (default is 20).
+ * @returns {Array} - Reaction rates for each process reaction.
+ */
+ compute_rates(state, T = 20) {
+ // 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_S, S_NH, , S_NO, S_HCO, , X_S, X_H, X_STO, X_A] = state;
+ 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;
+ const { theta_H, theta_STO, theta_mu_H, theta_b_H_O, theta_b_H_NO, theta_b_STO_O, theta_b_STO_NO, theta_mu_A, theta_b_A_O, theta_b_A_NO } = this.temp_params;
+
+ // Hydrolysis
+ rates[0] = X_H == 0 ? 0 : this._arrhenius(k_H, theta_H, T) * this._monod(X_S / X_H, K_X) * X_H;
+
+ // Heterotrophs
+ rates[1] = this._arrhenius(k_STO, theta_STO, T) * this._monod(S_O, K_O) * this._monod(S_S, K_S) * X_H;
+ rates[2] = this._arrhenius(k_STO, theta_STO, T) * 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] = X_H == 0 ? 0 : this._arrhenius(mu_H_max, theta_mu_H, T) * 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 : this._arrhenius(mu_H_max, theta_mu_H, T) * 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] = this._arrhenius(b_H_O, theta_b_H_O, T) * this._monod(S_O, K_O) * X_H;
+ rates[6] = this._arrhenius(b_H_NO, theta_b_H_NO, T) * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * X_H;
+ rates[7] = this._arrhenius(b_STO_O, theta_b_STO_O, T) * this._monod(S_O, K_O) * X_H;
+ rates[8] = this._arrhenius(b_STO_NO, theta_b_STO_NO, T) * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * X_STO;
+
+ // Autotrophs
+ rates[9] = this._arrhenius(mu_A_max, theta_mu_A, T) * this._monod(S_O, K_A_O) * this._monod(S_NH, K_A_NH) * this._monod(S_HCO, K_A_HCO) * X_A;
+ rates[10] = this._arrhenius(b_A_O, theta_b_A_O, T) * this._monod(S_O, K_O) * X_A;
+ rates[11] = this._arrhenius(b_A_NO, theta_b_A_NO, T) * this._inv_monod(S_O, K_A_O) * this._monod(S_NO, K_NO) * X_A;
+
+ return rates;
+ }
+
+ /**
+ * Computes the change in concentrations of reaction species based on the current state and temperature.
+ * @param {Array} state - State vector containing concentrations of reaction species.
+ * @param {number} [T=20] - Temperature in degrees Celsius (default is 20).
+ * @returns {Array} - Change in reaction species concentrations.
+ */
+ compute_dC(state, T = 20) { // 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, T));
+ }
+}
+
+module.exports = ASM3;
diff --git a/src/specificClass.js b/src/specificClass.js
index fb4eab6..17725e7 100644
--- a/src/specificClass.js
+++ b/src/specificClass.js
@@ -1,104 +1,104 @@
-const ASM3 = require('./reaction_modules/asm3_class.js');
-const { create, all, isArray } = require('mathjs');
-const { assertNoNaN } = require('./utils.js');
-const { childRegistrationUtils, logger, MeasurementContainer } = require('generalFunctions');
-const EventEmitter = require('events');
-
-const mathConfig = {
- matrix: 'Array' // use Array as the matrix type
-};
-
-const math = create(all, mathConfig);
-
-const S_O_INDEX = 0;
-const NUM_SPECIES = 13;
-const DEBUG = false;
-
-class Reactor {
- /**
- * Reactor base class.
- * @param {object} config - Configuration object containing reactor parameters.
- */
- constructor(config) {
- this.config = config;
- // EVOLV stuff
- this.logger = new logger(this.config.general.logging.enabled, this.config.general.logging.logLevel, config.general.name);
- this.emitter = new EventEmitter();
- this.measurements = new MeasurementContainer();
- this.upstreamReactor = null;
- this.childRegistrationUtils = new childRegistrationUtils(this); // Child registration utility
-
- this.asm = new ASM3();
-
- this.volume = config.volume; // fluid volume reactor [m3]
-
- 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 m-3]
- this.temperature = 20; // temperature [C]
-
- 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*60) * this.config.timeStep; // time step in seconds, converted to days.
- this.speedUpFactor = config.speedUpFactor ?? 1; // speed up factor for simulation
- }
-
- /**
- * 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;
- }
-
- /**
- * Setter for OTR (Oxygen Transfer Rate).
- * @param {object} input - Input object (msg) containing payload with OTR value [g O2 d-1 m-3].
- */
- set setOTR(input) {
- this.OTR = input.payload;
- }
-
- /**
- * Setter for reactor temperature [C].
- * Accepts either a direct numeric payload or { value } object payload.
- * @param {object} input - Input object (msg)
- */
- set setTemperature(input) {
- const payload = input?.payload;
- const rawValue = (payload && typeof payload === 'object' && payload.value !== undefined)
- ? payload.value
- : payload;
- const parsedValue = Number(rawValue);
- if (!Number.isFinite(parsedValue)) {
- this.logger.warn(`Invalid temperature input: ${rawValue}`);
- return;
- }
- this.temperature = parsedValue;
- }
-
- /**
- * Getter for effluent data.
- * @returns {object} Effluent data object (msg), defaults to inlet 0.
- */
- get getEffluent() { // getter for Effluent, defaults to inlet 0
- if (isArray(this.state.at(-1))) {
- 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 }, timestamp: this.currentTime };
- }
-
- get getGridProfile() { return null; }
-
- /**
- * 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 {number} - Calculated OTR [g O2 d-1 m-3].
- */
+const ASM3 = require('./reaction_modules/asm3_class.js');
+const { create, all, isArray } = require('mathjs');
+const { assertNoNaN } = require('./utils.js');
+const { childRegistrationUtils, logger, MeasurementContainer, POSITIONS } = require('generalFunctions');
+const EventEmitter = require('events');
+
+const mathConfig = {
+ matrix: 'Array' // use Array as the matrix type
+};
+
+const math = create(all, mathConfig);
+
+const S_O_INDEX = 0;
+const NUM_SPECIES = 13;
+const DEBUG = false;
+
+class Reactor {
+ /**
+ * Reactor base class.
+ * @param {object} config - Configuration object containing reactor parameters.
+ */
+ constructor(config) {
+ this.config = config;
+ // EVOLV stuff
+ this.logger = new logger(this.config.general.logging.enabled, this.config.general.logging.logLevel, config.general.name);
+ this.emitter = new EventEmitter();
+ this.measurements = new MeasurementContainer();
+ this.upstreamReactor = null;
+ this.childRegistrationUtils = new childRegistrationUtils(this); // Child registration utility
+
+ this.asm = new ASM3();
+
+ this.volume = config.volume; // fluid volume reactor [m3]
+
+ 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 m-3]
+ this.temperature = 20; // temperature [C]
+
+ 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*60) * this.config.timeStep; // time step in seconds, converted to days.
+ this.speedUpFactor = config.speedUpFactor ?? 1; // speed up factor for simulation
+ }
+
+ /**
+ * 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;
+ }
+
+ /**
+ * Setter for OTR (Oxygen Transfer Rate).
+ * @param {object} input - Input object (msg) containing payload with OTR value [g O2 d-1 m-3].
+ */
+ set setOTR(input) {
+ this.OTR = input.payload;
+ }
+
+ /**
+ * Setter for reactor temperature [C].
+ * Accepts either a direct numeric payload or { value } object payload.
+ * @param {object} input - Input object (msg)
+ */
+ set setTemperature(input) {
+ const payload = input?.payload;
+ const rawValue = (payload && typeof payload === 'object' && payload.value !== undefined)
+ ? payload.value
+ : payload;
+ const parsedValue = Number(rawValue);
+ if (!Number.isFinite(parsedValue)) {
+ this.logger.warn(`Invalid temperature input: ${rawValue}`);
+ return;
+ }
+ this.temperature = parsedValue;
+ }
+
+ /**
+ * Getter for effluent data.
+ * @returns {object} Effluent data object (msg), defaults to inlet 0.
+ */
+ get getEffluent() { // getter for Effluent, defaults to inlet 0
+ if (isArray(this.state.at(-1))) {
+ 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 }, timestamp: this.currentTime };
+ }
+
+ get getGridProfile() { return null; }
+
+ /**
+ * 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 {number} - Calculated OTR [g O2 d-1 m-3].
+ */
_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);
@@ -126,357 +126,356 @@ class Reactor {
}
return capRow(state);
}
-
- /**
- * 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(x => this._arrayClip2Zero(x));
- } else {
- return arr < 0 ? 0 : arr;
- }
- }
-
- registerChild(child, softwareType) {
- switch (softwareType) {
- case "measurement":
- this.logger.debug(`Registering measurement child.`);
- this._connectMeasurement(child);
- break;
- case "reactor":
- this.logger.debug(`Registering reactor child.`);
- this._connectReactor(child);
- break;
-
- default:
- this.logger.error(`Unrecognized softwareType: ${softwareType}`);
- }
- }
-
- _connectMeasurement(measurement) {
- if (!measurement) {
- this.logger.warn("Invalid measurement provided.");
- return;
- }
-
- let position;
- if (measurement.config.functionality.distance !== 'undefined') {
- position = measurement.config.functionality.distance;
- } else {
- position = measurement.config.functionality.positionVsParent;
- }
- const measurementType = measurement.config.asset.type;
- const key = `${measurementType}_${position}`;
- const eventName = `${measurementType}.measured.${position}`;
-
- // Register event listener for measurement updates
- measurement.measurements.emitter.on(eventName, (eventData) => {
- this.logger.debug(`${position} ${measurementType} from ${eventData.childName}: ${eventData.value} ${eventData.unit}`);
-
- // Store directly in parent's measurement container
- this.measurements
- .type(measurementType)
- .variant("measured")
- .position(position)
- .value(eventData.value, eventData.timestamp, eventData.unit);
-
- this._updateMeasurement(measurementType, eventData.value, position, eventData);
- });
- }
-
-
- _connectReactor(reactor) {
- if (!reactor) {
- this.logger.warn("Invalid reactor provided.");
- return;
- }
-
- this.upstreamReactor = reactor;
-
- reactor.emitter.on("stateChange", (data) => {
- this.logger.debug(`State change of upstream reactor detected.`);
- this.updateState(data);
- });
- }
-
-
- _updateMeasurement(measurementType, value, position, context) {
- this.logger.debug(`---------------------- updating ${measurementType} ------------------ `);
- switch (measurementType) {
- case "temperature":
- if (position == "atEquipment") {
- this.temperature = value;
- }
- break;
- default:
- this.logger.error(`Type '${measurementType}' not recognized for measured update.`);
- return;
- }
- }
-
- /**
- * 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 = Date.now()) { // expect update with timestamp
- const day2ms = 1000 * 60 * 60 * 24;
-
- if (this.upstreamReactor) {
- this.setInfluent = this.upstreamReactor.getEffluent;
- }
-
- 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(this.timeStep);
- n += 1;
- }
- this.currentTime += n_iter * this.timeStep * day2ms / this.speedUpFactor;
- this.emitter.emit("stateChange", this.currentTime);
- }
- }
-}
-
-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;
- }
-
- /**
- * 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 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, this.temperature);
- const transfer = Array(NUM_SPECIES).fill(0.0);
- transfer[S_O_INDEX] = isNaN(this.kla) ? this.OTR : this._calcOTR(this.state[S_O_INDEX], this.temperature); // 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)
+
+ /**
+ * 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(x => this._arrayClip2Zero(x));
+ } else {
+ return arr < 0 ? 0 : arr;
+ }
+ }
+
+ registerChild(child, softwareType) {
+ switch (softwareType) {
+ case "measurement":
+ this.logger.debug(`Registering measurement child.`);
+ this._connectMeasurement(child);
+ break;
+ case "reactor":
+ this.logger.debug(`Registering reactor child.`);
+ this._connectReactor(child);
+ break;
+
+ default:
+ this.logger.error(`Unrecognized softwareType: ${softwareType}`);
+ }
+ }
+
+ _connectMeasurement(measurement) {
+ if (!measurement) {
+ this.logger.warn("Invalid measurement provided.");
+ return;
+ }
+
+ let position;
+ if (measurement.config.functionality.distance !== 'undefined') {
+ position = measurement.config.functionality.distance;
+ } else {
+ position = measurement.config.functionality.positionVsParent;
+ }
+ const measurementType = measurement.config.asset.type;
+ const eventName = `${measurementType}.measured.${position}`;
+
+ // Register event listener for measurement updates
+ measurement.measurements.emitter.on(eventName, (eventData) => {
+ this.logger.debug(`${position} ${measurementType} from ${eventData.childName}: ${eventData.value} ${eventData.unit}`);
+
+ // Store directly in parent's measurement container
+ this.measurements
+ .type(measurementType)
+ .variant("measured")
+ .position(position)
+ .value(eventData.value, eventData.timestamp, eventData.unit);
+
+ this._updateMeasurement(measurementType, eventData.value, position, eventData);
+ });
+ }
+
+
+ _connectReactor(reactor) {
+ if (!reactor) {
+ this.logger.warn("Invalid reactor provided.");
+ return;
+ }
+
+ this.upstreamReactor = reactor;
+
+ reactor.emitter.on("stateChange", (data) => {
+ this.logger.debug(`State change of upstream reactor detected.`);
+ this.updateState(data);
+ });
+ }
+
+
+ _updateMeasurement(measurementType, value, position, _context) {
+ this.logger.debug(`---------------------- updating ${measurementType} ------------------ `);
+ switch (measurementType) {
+ case "temperature":
+ if (position == POSITIONS.AT_EQUIPMENT) {
+ this.temperature = value;
+ }
+ break;
+ default:
+ this.logger.error(`Type '${measurementType}' not recognized for measured update.`);
+ return;
+ }
+ }
+
+ /**
+ * 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 = Date.now()) { // expect update with timestamp
+ const day2ms = 1000 * 60 * 60 * 24;
+
+ if (this.upstreamReactor) {
+ this.setInfluent = this.upstreamReactor.getEffluent;
+ }
+
+ 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(this.timeStep);
+ n += 1;
+ }
+ this.currentTime += n_iter * this.timeStep * day2ms / this.speedUpFactor;
+ this.emitter.emit("stateChange", this.currentTime);
+ }
+ }
+}
+
+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;
+ }
+
+ /**
+ * 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 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, this.temperature);
+ const transfer = Array(NUM_SPECIES).fill(0.0);
+ transfer[S_O_INDEX] = isNaN(this.kla) ? this.OTR : this._calcOTR(this.state[S_O_INDEX], this.temperature); // 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)
this.state = this._capDissolvedOxygen(this._arrayClip2Zero(math.add(this.state, dC_total))); // clip concentrations and enforce physical DO saturation
- if(DEBUG){
- assertNoNaN(dC_total, "change in state");
- assertNoNaN(this.state, "new state");
- }
- return this.state;
- }
-}
-
-class Reactor_PFR extends Reactor {
- /**
- * Reactor_PFR class for Plug Flow Reactor.
- * @param {object} config - Configuration object containing reactor parameters.
- */
- constructor(config) {
- super(config);
-
- this.length = config.length; // reactor length [m]
- this.n_x = config.resolution_L; // number of slices
-
- this.d_x = this.length / this.n_x;
- this.A = this.volume / this.length; // crosssectional area [m2]
-
- this.alpha = config.alpha;
-
- this.state = Array.from(Array(this.n_x), () => config.initialState.slice())
-
- this.D = 0.0; // axial dispersion [m2 d-1]
-
- 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");
- }
-
- get getGridProfile() {
- return {
- grid: this.state.map(row => row.slice()),
- n_x: this.n_x,
- d_x: this.d_x,
- length: this.length,
- species: ['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'],
- timestamp: this.currentTime
- };
- }
-
- /**
- * 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;
- }
-
- updateState(newTime) {
- super.updateState(newTime);
- let Pe_local = this.d_x*math.sum(this.Fs)/(this.D*this.A)
- let Co_D = this.D*this.timeStep/(this.d_x*this.d_x);
-
- (Pe_local >= 2) && this.logger.warn(`Local Péclet number (${Pe_local}) is too high! Increase reactor resolution.`);
- (Co_D >= 0.5) && this.logger.warn(`Courant number (${Co_D}) is too high! Reduce time step size.`);
-
- if(DEBUG) {
- console.log("Inlet state max " + math.max(this.state[0]))
- console.log("Pe total " + this.length*math.sum(this.Fs)/(this.D*this.A));
- console.log("Pe local " + Pe_local);
- console.log("Co ad " + math.sum(this.Fs)*this.timeStep/(this.A*this.d_x));
- console.log("Co D " + Co_D);
- }
- }
-
- /**
- * 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, this.temperature));
- const transfer = Array.from(Array(this.n_x), () => new Array(NUM_SPECIES).fill(0));
-
- if (isNaN(this.kla)) { // calculate OTR if kla is not NaN, otherwise use externally calculated OTR
- for (let i = 1; i < this.n_x - 1; i++) {
- transfer[i][S_O_INDEX] = this.OTR * this.n_x/(this.n_x-2);
- }
- } else {
- for (let i = 1; i < this.n_x - 1; i++) {
- transfer[i][S_O_INDEX] = this._calcOTR(this.state[i][S_O_INDEX], this.temperature) * this.n_x/(this.n_x-2);
- }
- }
-
- const dC_total = math.multiply(math.add(dispersion, advection, reaction, transfer), time_step);
-
- const stateNew = math.add(this.state, dC_total);
- this._applyBoundaryConditions(stateNew);
-
- if (DEBUG) {
- assertNoNaN(dispersion, "dispersion");
- assertNoNaN(advection, "advection");
- assertNoNaN(reaction, "reaction");
- assertNoNaN(dC_total, "change in state");
- assertNoNaN(stateNew, "new state post BC");
- }
-
+ if(DEBUG){
+ assertNoNaN(dC_total, "change in state");
+ assertNoNaN(this.state, "new state");
+ }
+ return this.state;
+ }
+}
+
+class Reactor_PFR extends Reactor {
+ /**
+ * Reactor_PFR class for Plug Flow Reactor.
+ * @param {object} config - Configuration object containing reactor parameters.
+ */
+ constructor(config) {
+ super(config);
+
+ this.length = config.length; // reactor length [m]
+ this.n_x = config.resolution_L; // number of slices
+
+ this.d_x = this.length / this.n_x;
+ this.A = this.volume / this.length; // crosssectional area [m2]
+
+ this.alpha = config.alpha;
+
+ this.state = Array.from(Array(this.n_x), () => config.initialState.slice())
+
+ this.D = 0.0; // axial dispersion [m2 d-1]
+
+ 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");
+ }
+
+ get getGridProfile() {
+ return {
+ grid: this.state.map(row => row.slice()),
+ n_x: this.n_x,
+ d_x: this.d_x,
+ length: this.length,
+ species: ['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'],
+ timestamp: this.currentTime
+ };
+ }
+
+ /**
+ * 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;
+ }
+
+ updateState(newTime) {
+ super.updateState(newTime);
+ let Pe_local = this.d_x*math.sum(this.Fs)/(this.D*this.A)
+ let Co_D = this.D*this.timeStep/(this.d_x*this.d_x);
+
+ (Pe_local >= 2) && this.logger.warn(`Local Peclet number (${Pe_local}) is too high! Increase reactor resolution.`);
+ (Co_D >= 0.5) && this.logger.warn(`Courant number (${Co_D}) is too high! Reduce time step size.`);
+
+ if(DEBUG) {
+ console.log("Inlet state max " + math.max(this.state[0]))
+ console.log("Pe total " + this.length*math.sum(this.Fs)/(this.D*this.A));
+ console.log("Pe local " + Pe_local);
+ console.log("Co ad " + math.sum(this.Fs)*this.timeStep/(this.A*this.d_x));
+ console.log("Co D " + Co_D);
+ }
+ }
+
+ /**
+ * 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, this.temperature));
+ const transfer = Array.from(Array(this.n_x), () => new Array(NUM_SPECIES).fill(0));
+
+ if (isNaN(this.kla)) { // calculate OTR if kla is not NaN, otherwise use externally calculated OTR
+ for (let i = 1; i < this.n_x - 1; i++) {
+ transfer[i][S_O_INDEX] = this.OTR * this.n_x/(this.n_x-2);
+ }
+ } else {
+ for (let i = 1; i < this.n_x - 1; i++) {
+ transfer[i][S_O_INDEX] = this._calcOTR(this.state[i][S_O_INDEX], this.temperature) * this.n_x/(this.n_x-2);
+ }
+ }
+
+ const dC_total = math.multiply(math.add(dispersion, advection, reaction, transfer), time_step);
+
+ const stateNew = math.add(this.state, dC_total);
+ this._applyBoundaryConditions(stateNew);
+
+ 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._capDissolvedOxygen(this._arrayClip2Zero(stateNew));
return stateNew;
}
-
- _updateMeasurement(measurementType, value, position, context) {
- switch(measurementType) {
- case "quantity (oxygen)":
- if (!Number.isFinite(position) || !Number.isFinite(value) || this.config.length <= 0) {
- this.logger.warn(`Ignoring oxygen measurement update with invalid data (position=${position}, value=${value}).`);
- break;
- }
- {
- // Clamp sensor-derived position to valid PFR grid bounds.
- const rawIndex = Math.round(position / this.config.length * this.n_x);
- const grid_pos = Math.max(0, Math.min(this.n_x - 1, rawIndex));
- this.state[grid_pos][S_O_INDEX] = value; // reconcile measured oxygen concentration into nearest grid cell
- }
- break;
- default:
- super._updateMeasurement(measurementType, value, position, context);
- }
- }
-
- /**
- * 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) {
- 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_dispersion_term = (1-this.alpha)*this.D*this.A/(math.sum(this.Fs)*this.d_x);
- state[0] = math.multiply(1/(1+BC_dispersion_term), math.add(BC_C_in, math.multiply(BC_dispersion_term, state[1])));
- } else {
- state[0] = state[1];
- }
- // Neumann BC (no flux)
- state[this.n_x-1] = state[this.n_x-2];
- }
-
- /**
- * 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) {
- 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;
- D[1] = NearBoundary;
- NearBoundary.reverse();
- 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;
- } 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;
- }
- }
-
- /**
- * 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]);
- 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);
- return D2;
- }
-}
-
-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
-// 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);
-// 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 < 5000) {
-// console.log(Reactor.tick(0.001));
-// N += 1;
-// }
+
+ _updateMeasurement(measurementType, value, position, context) {
+ switch(measurementType) {
+ case "quantity (oxygen)":
+ if (!Number.isFinite(position) || !Number.isFinite(value) || this.config.length <= 0) {
+ this.logger.warn(`Ignoring oxygen measurement update with invalid data (position=${position}, value=${value}).`);
+ break;
+ }
+ {
+ // Clamp sensor-derived position to valid PFR grid bounds.
+ const rawIndex = Math.round(position / this.config.length * this.n_x);
+ const grid_pos = Math.max(0, Math.min(this.n_x - 1, rawIndex));
+ this.state[grid_pos][S_O_INDEX] = value; // reconcile measured oxygen concentration into nearest grid cell
+ }
+ break;
+ default:
+ super._updateMeasurement(measurementType, value, position, context);
+ }
+ }
+
+ /**
+ * 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) {
+ 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_dispersion_term = (1-this.alpha)*this.D*this.A/(math.sum(this.Fs)*this.d_x);
+ state[0] = math.multiply(1/(1+BC_dispersion_term), math.add(BC_C_in, math.multiply(BC_dispersion_term, state[1])));
+ } else {
+ state[0] = state[1];
+ }
+ // Neumann BC (no flux)
+ state[this.n_x-1] = state[this.n_x-2];
+ }
+
+ /**
+ * 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) {
+ 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;
+ D[1] = NearBoundary;
+ NearBoundary.reverse();
+ 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;
+ } 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;
+ }
+ }
+
+ /**
+ * 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]);
+ 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);
+ return D2;
+ }
+}
+
+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
+// 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);
+// 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 < 5000) {
+// console.log(Reactor.tick(0.001));
+// N += 1;
+// }
diff --git a/test/specificClass.test.js b/test/specificClass.test.js
new file mode 100644
index 0000000..93c4a4c
--- /dev/null
+++ b/test/specificClass.test.js
@@ -0,0 +1,346 @@
+/**
+ * Tests for reactor specificClass (domain logic).
+ *
+ * Two reactor classes are exported: Reactor_CSTR and Reactor_PFR.
+ * Both extend a base Reactor class.
+ *
+ * Key methods tested:
+ * - _calcOTR: oxygen transfer rate calculation
+ * - _arrayClip2Zero: clip negative values to zero
+ * - setInfluent / getEffluent: influent/effluent data flow
+ * - setOTR: external OTR override
+ * - tick (CSTR): forward Euler state update
+ * - tick (PFR): finite difference state update
+ * - registerChild: dispatches to measurement / reactor handlers
+ */
+
+const { Reactor_CSTR, Reactor_PFR } = require('../src/specificClass');
+
+// --------------- helpers ---------------
+
+const NUM_SPECIES = 13;
+
+function makeCSTRConfig(overrides = {}) {
+ return {
+ general: {
+ name: 'TestCSTR',
+ id: 'cstr-test-1',
+ logging: { enabled: false, logLevel: 'error' },
+ },
+ functionality: {
+ softwareType: 'reactor',
+ positionVsParent: 'atEquipment',
+ },
+ volume: 1000,
+ n_inlets: 1,
+ kla: 240,
+ timeStep: 1, // 1 second
+ initialState: new Array(NUM_SPECIES).fill(1.0),
+ ...overrides,
+ };
+}
+
+function makePFRConfig(overrides = {}) {
+ return {
+ general: {
+ name: 'TestPFR',
+ id: 'pfr-test-1',
+ logging: { enabled: false, logLevel: 'error' },
+ },
+ functionality: {
+ softwareType: 'reactor',
+ positionVsParent: 'atEquipment',
+ },
+ volume: 200,
+ length: 10,
+ resolution_L: 10,
+ n_inlets: 1,
+ kla: 240,
+ alpha: 0.5,
+ timeStep: 1,
+ initialState: new Array(NUM_SPECIES).fill(0.1),
+ ...overrides,
+ };
+}
+
+// --------------- CSTR tests ---------------
+
+describe('Reactor_CSTR', () => {
+
+ describe('constructor / initialization', () => {
+ it('should create an instance and set state from initialState', () => {
+ const r = new Reactor_CSTR(makeCSTRConfig());
+ expect(r).toBeDefined();
+ expect(r.state).toEqual(new Array(NUM_SPECIES).fill(1.0));
+ });
+
+ it('should initialize Fs and Cs_in arrays based on n_inlets', () => {
+ const r = new Reactor_CSTR(makeCSTRConfig({ n_inlets: 3 }));
+ expect(r.Fs).toHaveLength(3);
+ expect(r.Cs_in).toHaveLength(3);
+ expect(r.Fs.every(v => v === 0)).toBe(true);
+ });
+
+ it('should store volume from config', () => {
+ const r = new Reactor_CSTR(makeCSTRConfig({ volume: 500 }));
+ expect(r.volume).toBe(500);
+ });
+
+ it('should initialize temperature to 20', () => {
+ const r = new Reactor_CSTR(makeCSTRConfig());
+ expect(r.temperature).toBe(20);
+ });
+ });
+
+ describe('_calcOTR()', () => {
+ let r;
+ beforeAll(() => { r = new Reactor_CSTR(makeCSTRConfig({ kla: 240 })); });
+
+ it('should return a positive value when S_O < saturation', () => {
+ const otr = r._calcOTR(0, 20);
+ expect(otr).toBeGreaterThan(0);
+ });
+
+ it('should return approximately zero when S_O equals saturation', () => {
+ // S_O_sat at T=20: 14.652 - 4.1022e-1*20 + 7.9910e-3*400 + 7.7774e-5*8000
+ const T = 20;
+ const S_O_sat = 14.652 - 4.1022e-1 * T + 7.9910e-3 * T * T + 7.7774e-5 * T * T * T;
+ const otr = r._calcOTR(S_O_sat, T);
+ expect(otr).toBeCloseTo(0, 5);
+ });
+
+ it('should return a negative value when S_O > saturation (supersaturated)', () => {
+ const otr = r._calcOTR(100, 20);
+ expect(otr).toBeLessThan(0);
+ });
+
+ it('should use T=20 as default temperature', () => {
+ const otr1 = r._calcOTR(0);
+ const otr2 = r._calcOTR(0, 20);
+ expect(otr1).toBe(otr2);
+ });
+ });
+
+ describe('_arrayClip2Zero()', () => {
+ let r;
+ beforeAll(() => { r = new Reactor_CSTR(makeCSTRConfig()); });
+
+ it('should clip negative values to zero', () => {
+ expect(r._arrayClip2Zero([-5, 3, -1, 0, 7])).toEqual([0, 3, 0, 0, 7]);
+ });
+
+ it('should leave all-positive arrays unchanged', () => {
+ expect(r._arrayClip2Zero([1, 2, 3])).toEqual([1, 2, 3]);
+ });
+
+ it('should handle nested arrays (2D)', () => {
+ const result = r._arrayClip2Zero([[-1, 2], [3, -4]]);
+ expect(result).toEqual([[0, 2], [3, 0]]);
+ });
+
+ it('should handle a single scalar', () => {
+ expect(r._arrayClip2Zero(-5)).toBe(0);
+ expect(r._arrayClip2Zero(5)).toBe(5);
+ });
+ });
+
+ describe('setInfluent / getEffluent', () => {
+ it('should store influent data via setter', () => {
+ const r = new Reactor_CSTR(makeCSTRConfig({ n_inlets: 2 }));
+ const input = {
+ payload: {
+ inlet: 0,
+ F: 100,
+ C: new Array(NUM_SPECIES).fill(5),
+ },
+ };
+ r.setInfluent = input;
+ expect(r.Fs[0]).toBe(100);
+ expect(r.Cs_in[0]).toEqual(new Array(NUM_SPECIES).fill(5));
+ });
+
+ it('should return effluent with the sum of Fs and the current state', () => {
+ const r = new Reactor_CSTR(makeCSTRConfig());
+ r.Fs[0] = 50;
+ const eff = r.getEffluent;
+ expect(eff.topic).toBe('Fluent');
+ expect(eff.payload.F).toBe(50);
+ expect(eff.payload.C).toEqual(r.state);
+ });
+ });
+
+ describe('setOTR', () => {
+ it('should set the OTR value', () => {
+ const r = new Reactor_CSTR(makeCSTRConfig({ kla: NaN }));
+ r.setOTR = { payload: 42 };
+ expect(r.OTR).toBe(42);
+ });
+ });
+
+ describe('tick()', () => {
+ it('should return a new state array of correct length', () => {
+ const r = new Reactor_CSTR(makeCSTRConfig());
+ const result = r.tick(0.001);
+ expect(result).toHaveLength(NUM_SPECIES);
+ });
+
+ it('should not produce NaN values', () => {
+ const r = new Reactor_CSTR(makeCSTRConfig());
+ r.Fs[0] = 10;
+ r.Cs_in[0] = new Array(NUM_SPECIES).fill(5);
+ const result = r.tick(0.001);
+ result.forEach(v => expect(Number.isNaN(v)).toBe(false));
+ });
+
+ it('should not produce negative concentrations', () => {
+ const r = new Reactor_CSTR(makeCSTRConfig());
+ // Run multiple ticks
+ for (let i = 0; i < 100; i++) {
+ r.tick(0.001);
+ }
+ r.state.forEach(v => expect(v).toBeGreaterThanOrEqual(0));
+ });
+
+ it('should reach steady state with zero flow (concentrations change only via reaction)', () => {
+ const r = new Reactor_CSTR(makeCSTRConfig());
+ // No inflow
+ const initial = [...r.state];
+ r.tick(0.0001);
+ // State should have changed due to reaction/OTR
+ const changed = r.state.some((v, i) => v !== initial[i]);
+ expect(changed).toBe(true);
+ });
+ });
+
+ describe('registerChild()', () => {
+ it('should not throw for "measurement" software type', () => {
+ const r = new Reactor_CSTR(makeCSTRConfig());
+ // Passing null child will trigger warn but not crash
+ expect(() => r.registerChild(null, 'measurement')).not.toThrow();
+ });
+
+ it('should not throw for "reactor" software type', () => {
+ const r = new Reactor_CSTR(makeCSTRConfig());
+ expect(() => r.registerChild(null, 'reactor')).not.toThrow();
+ });
+
+ it('should not throw for unknown software type', () => {
+ const r = new Reactor_CSTR(makeCSTRConfig());
+ expect(() => r.registerChild(null, 'unknown')).not.toThrow();
+ });
+ });
+});
+
+// --------------- PFR tests ---------------
+
+describe('Reactor_PFR', () => {
+
+ describe('constructor / initialization', () => {
+ it('should create an instance with 2D state grid', () => {
+ const r = new Reactor_PFR(makePFRConfig());
+ expect(r).toBeDefined();
+ expect(r.state).toHaveLength(10); // resolution_L = 10
+ expect(r.state[0]).toHaveLength(NUM_SPECIES);
+ });
+
+ it('should compute d_x = length / n_x', () => {
+ const r = new Reactor_PFR(makePFRConfig({ length: 10, resolution_L: 5 }));
+ expect(r.d_x).toBe(2);
+ });
+
+ it('should compute cross-sectional area A = volume / length', () => {
+ const r = new Reactor_PFR(makePFRConfig({ volume: 200, length: 10 }));
+ expect(r.A).toBe(20);
+ });
+
+ it('should initialize D (dispersion) to 0', () => {
+ const r = new Reactor_PFR(makePFRConfig());
+ expect(r.D).toBe(0);
+ });
+
+ it('should create derivative operators of correct size', () => {
+ const r = new Reactor_PFR(makePFRConfig({ resolution_L: 8 }));
+ expect(r.D_op).toHaveLength(8);
+ expect(r.D_op[0]).toHaveLength(8);
+ expect(r.D2_op).toHaveLength(8);
+ expect(r.D2_op[0]).toHaveLength(8);
+ });
+ });
+
+ describe('setDispersion', () => {
+ it('should set the axial dispersion value', () => {
+ const r = new Reactor_PFR(makePFRConfig());
+ r.setDispersion = { payload: 0.5 };
+ expect(r.D).toBe(0.5);
+ });
+ });
+
+ describe('tick()', () => {
+ it('should return a 2D state grid of correct dimensions', () => {
+ const r = new Reactor_PFR(makePFRConfig());
+ r.D = 0.01;
+ const result = r.tick(0.0001);
+ expect(result).toHaveLength(10);
+ expect(result[0]).toHaveLength(NUM_SPECIES);
+ });
+
+ it('should not produce NaN values with small time step and dispersion', () => {
+ const r = new Reactor_PFR(makePFRConfig());
+ r.D = 0.01;
+ r.Fs[0] = 10;
+ r.Cs_in[0] = new Array(NUM_SPECIES).fill(5);
+ const result = r.tick(0.0001);
+ result.forEach(row => {
+ row.forEach(v => expect(Number.isNaN(v)).toBe(false));
+ });
+ });
+
+ it('should not produce negative concentrations', () => {
+ const r = new Reactor_PFR(makePFRConfig());
+ r.D = 0.01;
+ for (let i = 0; i < 10; i++) {
+ r.tick(0.0001);
+ }
+ r.state.forEach(row => {
+ row.forEach(v => expect(v).toBeGreaterThanOrEqual(0));
+ });
+ });
+ });
+
+ describe('_applyBoundaryConditions()', () => {
+ it('should apply Neumann BC at outlet (last = second to last)', () => {
+ const r = new Reactor_PFR(makePFRConfig({ resolution_L: 5 }));
+ const state = Array.from({ length: 5 }, () => new Array(NUM_SPECIES).fill(1));
+ state[3] = new Array(NUM_SPECIES).fill(7);
+ r._applyBoundaryConditions(state);
+ // outlet BC: state[4] = state[3]
+ expect(state[4]).toEqual(new Array(NUM_SPECIES).fill(7));
+ });
+
+ it('should apply Neumann BC at inlet when no flow', () => {
+ const r = new Reactor_PFR(makePFRConfig({ resolution_L: 5 }));
+ r.Fs[0] = 0;
+ const state = Array.from({ length: 5 }, () => new Array(NUM_SPECIES).fill(1));
+ state[1] = new Array(NUM_SPECIES).fill(3);
+ r._applyBoundaryConditions(state);
+ // No flow: state[0] = state[1]
+ expect(state[0]).toEqual(new Array(NUM_SPECIES).fill(3));
+ });
+ });
+
+ describe('_arrayClip2Zero() (inherited)', () => {
+ it('should clip 2D arrays correctly', () => {
+ const r = new Reactor_PFR(makePFRConfig());
+ const result = r._arrayClip2Zero([[-1, 2], [3, -4]]);
+ expect(result).toEqual([[0, 2], [3, 0]]);
+ });
+ });
+
+ describe('_calcOTR() (inherited)', () => {
+ it('should work the same as in CSTR', () => {
+ const r = new Reactor_PFR(makePFRConfig({ kla: 240 }));
+ const otr = r._calcOTR(0, 20);
+ expect(otr).toBeGreaterThan(0);
+ });
+ });
+});