Refactor measurement position assignment and update grid position calculation in Reactor classes to align with new generalFunctions

Functioning code, requires improved sequencing
Fix boundary conditions for advection
2025-10-10 11:27:55 +02:00 · 2025-10-02 17:39:31 +02:00 · 2025-10-02 13:48:47 +02:00 · 2025-10-01 16:50:48 +02:00 · 2025-09-30 15:36:25 +02:00 · 2025-09-29 16:58:46 +02:00
6 changed files with 107 additions and 98 deletions
--- a/additional_nodes/recirculation-pump.js
+++ b/additional_nodes/recirculation-pump.js
@@ -3,12 +3,13 @@ module.exports = function(RED) {
        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": {
+                case "Fluent":
                    // conserve volume flow debit
                    let F_in = msg.payload.F;
                    let F1 = Math.max(F_in - F2, 0);
@@ -23,7 +24,6 @@ module.exports = function(RED) {

                    send([msg_F1, msg_F2]);
                    break;
-                }
                case "clock":
                    break;
                default:
--- a/additional_nodes/settling-basin.js
+++ b/additional_nodes/settling-basin.js
@@ -3,12 +3,13 @@ module.exports = function(RED) {
        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": {
+                case "Fluent":
                    // conserve volume flow debit
                    let F_in = msg.payload.F;
                    let C_in = msg.payload.C;
@@ -40,7 +41,6 @@ module.exports = function(RED) {

                    send([msg_F1, msg_F2]);
                    break;
-                }
                case "clock":
                    break;
                default:
--- a/src/nodeClass.js
+++ b/src/nodeClass.js
@@ -1,5 +1,4 @@
 const { Reactor_CSTR, Reactor_PFR } = require('./specificClass.js');
-const { configManager } = require('generalFunctions');


 class nodeClass {
@@ -49,13 +48,12 @@ class nodeClass {
                case "Dispersion":
                    this.source.setDispersion = msg;
                    break;
-                case 'registerChild': {
+                case 'registerChild':
                    // Register this node as a parent of the child node
                    const childId = msg.payload;
-                    const childObj = this.RED.nodes.getNode(childId);
+                    const childObj = this.RED.nodes.getNode(childId); 
                    this.source.childRegistrationUtils.registerChild(childObj.source, msg.positionVsParent);
                    break;
-                }
                default:
                    console.log("Unknown topic: " + msg.topic);
            }
@@ -71,10 +69,20 @@ class nodeClass {
     * @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, {
+        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),
@@ -98,7 +106,7 @@ class nodeClass {
                parseFloat(uiConfig.X_TS_init)
            ],
            timeStep: parseFloat(uiConfig.timeStep)
-        });
+        }
    }

    /**
@@ -129,7 +137,7 @@ class nodeClass {
                new_reactor = new Reactor_PFR(this.config);
                break;
            default:
-                console.warn("Unknown reactor type: " + this.config.reactor_type);
+                console.warn("Unknown reactor type: " + uiConfig.reactor_type);
        }

        this.source = new_reactor; // protect from reassignment
--- a/src/reaction_modules/asm3_class
+++ b/src/reaction_modules/asm3_class
@@ -19,7 +19,7 @@ class ASM3 {
      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.0, // saturation constant S_s        [g COD 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]
@@ -171,7 +171,7 @@ class ASM3 {
  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 [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;

--- a/src/reaction_modules/asm3_class.js
+++ b/src/reaction_modules/asm3_class.js
@@ -171,7 +171,7 @@ class ASM3 {
  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 [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;

--- a/src/specificClass.js
+++ b/src/specificClass.js
@@ -1,7 +1,7 @@
 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 { childRegistrationUtils, logger, MeasurementContainer } = require('generalFunctions');
 const EventEmitter = require('events');

 const mathConfig = {
@@ -12,6 +12,7 @@ const math = create(all, mathConfig);

 const S_O_INDEX = 0;
 const NUM_SPECIES = 13;
+const BC_PADDING = 2;
 const DEBUG = false;

 class Reactor {
@@ -27,6 +28,10 @@ class Reactor {
    this.measurements = new MeasurementContainer();
    this.upstreamReactor = null;
    this.childRegistrationUtils = new childRegistrationUtils(this); // Child registration utility
+    this.parent = []; // Gets assigned via child registration
+
+    this.upstreamReactor = null;
+    this.downstreamReactor = null;

    this.asm = new ASM3();

@@ -41,7 +46,7 @@ class Reactor {

    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 = 60; // speed up factor for simulation, 60 means 1 minute per simulated second
+    this.speedUpFactor = 100; // speed up factor for simulation, 60 means 1 minute per simulated second
  }

  /**
@@ -113,23 +118,18 @@ class Reactor {
    }
  }

-  _connectMeasurement(measurement) {
-    if (!measurement) {
+  _connectMeasurement(measurementChild) {
+    if (!measurementChild) {
      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 position = measurementChild.config.functionality.positionVsParent;
+    const measurementType = measurementChild.config.asset.type;
    const eventName = `${measurementType}.measured.${position}`;

    // Register event listener for measurement updates
-    measurement.measurements.emitter.on(eventName, (eventData) => {
+    measurementChild.measurements.emitter.on(eventName, (eventData) => {
      this.logger.debug(`${position} ${measurementType} from ${eventData.childName}: ${eventData.value} ${eventData.unit}`);

      // Store directly in parent's measurement container
@@ -144,26 +144,27 @@ class Reactor {
  }


-  _connectReactor(reactor) {
-    if (!reactor) {
+  _connectReactor(reactorChild) {
+    if (!reactorChild) {
      this.logger.warn("Invalid reactor provided.");
      return;
    }

-    this.upstreamReactor = reactor;
+    this.upstreamReactor = reactorChild;
+    reactorChild.downstreamReactor = this;

-    reactor.emitter.on("stateChange", (data) => {
+    reactorChild.emitter.on("stateChange", (data) => {
      this.logger.debug(`State change of upstream reactor detected.`);
      this.updateState(data);
    });
  }

  
-  _updateMeasurement(measurementType, value, position, _context) {
+  _updateMeasurement(measurementType, value, position, context) {
    this.logger.debug(`---------------------- updating ${measurementType} ------------------ `);
    switch (measurementType) {
      case "temperature":
-        if (position == POSITIONS.AT_EQUIPMENT) {
+        if (position == "atEquipment") {
          this.temperature = value;
        }
        break;
@@ -245,11 +246,15 @@ class Reactor_PFR extends Reactor {

    this.alpha = config.alpha;

-    this.state = Array.from(Array(this.n_x), () => config.initialState.slice())
+    this.state = Array.from(Array(this.n_x), () => config.initialState.slice());
+    this.extendedState = Array.from(Array(this.n_x + 2*BC_PADDING), () => new Array(NUM_SPECIES).fill(0));
+
+    // initialise extended state
+    this.state.forEach((row, i) => this.extendedState[i+BC_PADDING] = row);

    this.D = 0.0; // axial dispersion [m2 d-1]

-    this.D_op = this._makeDoperator(true, true);
+    this.D_op = this._makeDoperator();
    assertNoNaN(this.D_op, "Derivative operator");

    this.D2_op = this._makeD2operator();
@@ -287,25 +292,26 @@ class Reactor_PFR extends Reactor {
   * @returns {Array} - New reactor state.
   */
  tick(time_step) {
-    const dispersion = math.multiply(this.D / (this.d_x*this.d_x), this.D2_op, this.state);
-    const advection = math.multiply(-1 * math.sum(this.Fs) / (this.A*this.d_x), this.D_op, this.state);
-    const 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));
+    this._applyBoundaryConditions();
+
+    const dispersion = math.multiply(this.D / (this.d_x*this.d_x), this.D2_op, this.extendedState);
+    const advection = math.multiply(-1 * math.sum(this.Fs) / (this.A*this.d_x), this.D_op, this.extendedState);
+    const reaction = this.extendedState.map((state_slice) => this.asm.compute_dC(state_slice, this.temperature));
+    const transfer = Array.from(Array(this.n_x+2*BC_PADDING), () => 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++) {
+      for (let i = BC_PADDING+1; i < BC_PADDING+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);
+    for (let i = BC_PADDING+1; i < BC_PADDING+this.n_x - 1; i++) {
+        transfer[i][S_O_INDEX] = this._calcOTR(this.extendedState[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 dC_total = math.multiply(math.add(dispersion, advection, reaction, transfer).slice(BC_PADDING, this.n_x+BC_PADDING), time_step);

    const stateNew = math.add(this.state, dC_total);
-    this._applyBoundaryConditions(stateNew);

    if (DEBUG) {
      assertNoNaN(dispersion, "dispersion");
@@ -316,18 +322,18 @@ class Reactor_PFR extends Reactor {
    }

    this.state = this._arrayClip2Zero(stateNew);
+    this.state.forEach((row, i) => this.extendedState[i+BC_PADDING] = row);
    return stateNew;
  }

-  _updateMeasurement(measurementType, value, position, _context) {
+  _updateMeasurement(measurementType, value, position, context) {
    switch(measurementType) {
-      case "quantity (oxygen)": {
-        let grid_pos = Math.round(position / this.config.length * this.n_x);
+      case "quantity (oxygen)":
+        let grid_pos = Math.round(context.distance / this.config.length * this.n_x);
        this.state[grid_pos][S_O_INDEX] = value; // naive approach for reconciling measurements and simulation
        break;
-      }
      default:
-        super._updateMeasurement(measurementType, value, position, _context);
+        super._updateMeasurement(measurementType, value, position, context);
    }
  }

@@ -335,57 +341,51 @@ class Reactor_PFR extends Reactor {
   * 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])));
+  _applyBoundaryConditions() {
+    if (this.upstreamReactor) {
+      for (let i = 0; i < BC_PADDING; i++) {
+        this.extendedState[i] = this.upstreamReactor.state.at(i-BC_PADDING);
+      }
    } else {
-      state[0] = state[1];
+      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);
+        this.extendedState[BC_PADDING] = math.multiply(1/(1+BC_dispersion_term), math.add(BC_C_in, math.multiply(BC_dispersion_term, this.extendedState[BC_PADDING+1])));
+        this.extendedState[BC_PADDING-1] = math.add(math.multiply(2, this.extendedState[BC_PADDING]), math.multiply(-2, this.extendedState[BC_PADDING+2]), this.extendedState[BC_PADDING+3]);
+      } else {
+        for (let i = 0; i < BC_PADDING; i++) {
+          this.extendedState[i] = this.extendedState[BC_PADDING];
+        }
+      }
+    }
+
+    if (this.downstreamReactor) {
+      for (let i = 0; i < BC_PADDING; i++) {
+        this.extendedState[this.n_x+BC_PADDING+i] = this.downstreamReactor.state[i];
+      }
+    } else {
+      // Neumann BC (no flux)
+      for (let i = 0; i < BC_PADDING; i++) {
+        this.extendedState[BC_PADDING+this.n_x+i] = this.extendedState.at(-1-BC_PADDING);
+      }
    }
-    // 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;
-    }
+  _makeDoperator() { // create gradient operator
+    const D_size = this.n_x+2*BC_PADDING;
+    const I = math.resize(math.diag(Array(D_size).fill(1/12), -2), [D_size, D_size]);
+    const A = math.resize(math.diag(Array(D_size).fill(-2/3), -1), [D_size, D_size]);
+    const B = math.resize(math.diag(Array(D_size).fill(2/3), 1), [D_size, D_size]);
+    const C = math.resize(math.diag(Array(D_size).fill(-1/12), 2), [D_size, D_size]);
+    const D = math.add(I, A, B, C);
+    // set by BCs elsewhere
+    D.forEach((row, i) => i < BC_PADDING || i >= this.n_x+BC_PADDING ? row.fill(0) : row);
+    return D;
  }

  /**
@@ -393,12 +393,13 @@ class Reactor_PFR extends Reactor {
   * @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 D_size = this.n_x+2*BC_PADDING;
+    const I = math.diag(Array(D_size).fill(-2), 0);
+    const A = math.resize(math.diag(Array(D_size).fill(1), 1), [D_size, D_size]);
+    const B = math.resize(math.diag(Array(D_size).fill(1), -1), [D_size, D_size]);
    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);
+    // set by BCs elsewhere
+    D2.forEach((row, i) => i < BC_PADDING || i >= this.n_x+BC_PADDING ? row.fill(0) : row);
    return D2;
  }
 }
Author	SHA1	Message	Date
p.vanderwilt	2a520be33b	Refactor measurement position assignment and update grid position calculation in Reactor classes to align with new generalFunctions	2025-10-10 11:27:55 +02:00
p.vanderwilt	baecf2f599	Functioning code, requires improved sequencing	2025-10-02 17:39:31 +02:00
p.vanderwilt	cd3a19e66f	Fix boundary conditions for advection	2025-10-02 13:48:47 +02:00
p.vanderwilt	3aea0e55c4	Rewrite for improved boundary condition	2025-10-01 16:50:48 +02:00
p.vanderwilt	d9511dc3c7	Implement simple BCs	2025-09-30 15:36:25 +02:00
p.vanderwilt	993482f8c0	Deal with mulitple parents and set downstreamReactor for improved boundary conditions	2025-09-29 16:58:46 +02:00
p.vanderwilt	5f4ebdc2af	Fix reference error and improve child variable naming	2025-09-29 15:45:07 +02:00
p.vanderwilt	6c79d0ef9b	Use improved boundary conditions for upstream and downstream reactors	2025-09-29 15:34:54 +02:00