Merge pull request 'Switch temperature to Measurement node and finally fix boundary conditions properly' (#14) from experimental into main

Reviewed-on: p.vanderwilt/asm3#14

package-lock.json

@@ -9,6 +9,7 @@
       "version": "0.0.1",
       "license": "SEE LICENSE",
       "dependencies": {
+        "generalFunctions": "git+https://gitea.centraal.wbd-rd.nl/p.vanderwilt/generalFunctions.git",
         "mathjs": "^14.5.2"
       }
     },
@@ -59,6 +60,11 @@
         "url": "https://github.com/sponsors/rawify"
       }
     },
+    "node_modules/generalFunctions": {
+      "version": "1.0.0",
+      "resolved": "git+https://gitea.centraal.wbd-rd.nl/p.vanderwilt/generalFunctions.git#950ca2b6b4e91b37479aee90bff74b02c16f130e",
+      "license": "SEE LICENSE"
+    },
     "node_modules/javascript-natural-sort": {
       "version": "0.7.1",
       "resolved": "https://registry.npmjs.org/javascript-natural-sort/-/javascript-natural-sort-0.7.1.tgz",

package.json

@@ -27,6 +27,7 @@
     }
   },
   "dependencies": {
+    "generalFunctions": "git+https://gitea.centraal.wbd-rd.nl/p.vanderwilt/generalFunctions.git",
     "mathjs": "^14.5.2"
   }
 }

src/nodeClass.js

@@ -42,9 +42,18 @@ class nodeClass {
                 case "OTR":
                     this.reactor.setOTR = msg;
                     break;
+                case "Temperature":
+                    this.reactor.setTemperature = msg;
+                    break;
                 case "Dispersion":
                     this.reactor.setDispersion = msg;
                     break;
+                case 'registerChild':
+                    // Register this node as a child of the parent node
+                    const childId = msg.payload;
+                    const childObj = this.RED.nodes.getNode(childId); 
+                    this.reactor.childRegistrationUtils.registerChild(childObj.source, msg.positionVsParent);
+                    break;
                 default:
                     console.log("Unknown topic: " + msg.topic);
             }

src/specificClass.js

@@ -1,6 +1,7 @@
 const ASM3 = require('./reaction_modules/asm3_class.js');
 const { create, all } = require('mathjs');
 const { assertNoNaN } = require('./utils.js');
+const { childRegistrationUtils, logger, MeasurementContainer } = require('generalFunctions');
 
 const config = {
   matrix: 'Array' // use Array as the matrix type
@@ -10,7 +11,7 @@ const math = create(all, config);
 
 const S_O_INDEX = 0;
 const NUM_SPECIES = 13;
-const DEBUG = true;
+const DEBUG = false;
 
 class Reactor {
   /**
@@ -18,6 +19,11 @@ class Reactor {
    * @param {object} config - Configuration object containing reactor parameters.
    */
   constructor(config) {
+    // EVOLV stuff
+    this.logger = new logger(); //TODO: attach config
+    this.measurements = new MeasurementContainer();
+    this.childRegistrationUtils = new childRegistrationUtils(this); // Child registration utility
+
     this.asm = new ASM3();
 
     this.volume = config.volume; // fluid volume reactor [m3]
@@ -25,12 +31,25 @@ class Reactor {
     this.Fs = Array(config.n_inlets).fill(0); // fluid debits per inlet [m3 d-1]
     this.Cs_in = Array.from(Array(config.n_inlets), () => new Array(NUM_SPECIES).fill(0)); // composition influents
     this.OTR = 0.0; // oxygen transfer rate [g O2 d-1]
+    this.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*15); // time step [d]
-        this.speedUpFactor = 1; // speed up factor for simulation, 60 means 1 minute per simulated second
+    this.speedUpFactor = 60; // speed up factor for simulation, 60 means 1 minute per simulated second
+  }
+
+  updateMeasurement(variant, subType, value, position) {
+    this.logger.debug(`---------------------- updating ${subType} ------------------ `);
+    switch (subType) {
+      case "temperature":
+        this.temperature = value;
+        break;
+      default:
+        this.logger.error(`Type '${subType}' not recognized for measured update.`);
+        return;
+    }
   }
 
   /**
@@ -92,7 +111,6 @@ class Reactor {
       this.currentTime += n_iter * this.timeStep * day2ms / this.speedUpFactor;
     }
   }
-
 }
 
 class Reactor_CSTR extends Reactor {
@@ -121,9 +139,9 @@ class Reactor_CSTR extends Reactor {
   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);
+      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]); // calculate OTR if kla is not NaN, otherwise use externaly calculated OTR
+      transfer[S_O_INDEX] = isNaN(this.kla) ? this.OTR : this._calcOTR(this.state[S_O_INDEX], 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._arrayClip2Zero(math.add(this.state, dC_total)); // clip value element-wise to avoid negative concentrations
@@ -173,21 +191,6 @@ class Reactor_PFR extends Reactor {
     this.D = input.payload;
   }
 
-    /**
-     * Setter for influent data.
-     * @param {object} input - Input object (msg) containing payload with inlet index, flow rate, and concentrations.
-    */
-    set setInfluent(input) {
-        super.setInfluent = input;
-        if(DEBUG) {
-            console.log("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 " + this.d_x*math.sum(this.Fs)/(this.D*this.A));
-            console.log("Co ad " + math.sum(this.Fs)*this.timeStep/(this.A*this.d_x));
-            console.log("Co D " + this.D*this.timeStep/(this.d_x*this.d_x));
-        }
-    }
-
   /**
    * Getter for effluent data.
    * @returns {object} Effluent data object (msg), defaults to inlet 0.
@@ -196,25 +199,21 @@ class Reactor_PFR extends Reactor {
     return { topic: "Fluent", payload: { inlet: 0, F: math.sum(this.Fs), C: this.state.at(-1) }, timestamp: this.currentTime };
   }
 
-    /**
+  updateState(newTime) {
-     * Apply boundary conditions to the reactor state.
+    super.updateState(newTime);
-     *  for inlet, apply generalised Danckwerts BC, if there is not flow, apply Neumann BC with no flux
+    let Pe_local = this.d_x*math.sum(this.Fs)/(this.D*this.A)
-     *  for outlet, apply regular Danckwerts BC (Neumann BC with no flux)
+    let Co_D = this.D*this.timeStep/(this.d_x*this.d_x);
-     * @param {Array} state - Current reactor state without enforced BCs.
+
-     */
+    (Pe_local >= 2) && console.warn(`Local Péclet number (${Pe_local}) is too high! Increase reactor resolution.`);
-    _applyBoundaryConditions(state) {
+    (Co_D >= 0.5) && console.warn(`Courant number (${Co_D}) is too high! Reduce time step size.`);
-        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];
+    if(DEBUG) {
-            const BC_gradient = Array(this.n_x).fill(0);
+      console.log("Inlet state max " + math.max(this.state[0]))
-            BC_gradient[0] = -1;
+      console.log("Pe total " + this.length*math.sum(this.Fs)/(this.D*this.A));
-            BC_gradient[1] = 1;
+      console.log("Pe local " + Pe_local);
-            const BC_dispersion = math.multiply((1 - this.alpha) * this.D*this.A / (math.sum(this.Fs) * this.d_x), [BC_gradient], state)[0];
+      console.log("Co ad " + math.sum(this.Fs)*this.timeStep/(this.A*this.d_x));
-            state[0] = math.add(BC_C_in, BC_dispersion).map(val => val < 0 ? 0 : val);
+      console.log("Co D " + Co_D);
-        } else { // Neumann BC (no flux)
-            state[0] = state[1];
     }
-        // Neumann BC (no flux)
-        state[this.n_x-1] = state[this.n_x-2]
   }
 
   /**
@@ -225,13 +224,13 @@ class Reactor_PFR extends Reactor {
   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));
+    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
       transfer.forEach((x) => { x[S_O_INDEX] = this.OTR; });
     } else {
-            transfer.forEach((x, i) => { x[S_O_INDEX] = this._calcOTR(this.state[i][S_O_INDEX]); });
+      transfer.forEach((x, i) => { x[S_O_INDEX] = this._calcOTR(this.state[i][S_O_INDEX], this.temperature); });
     }
 
     const dC_total = math.multiply(math.add(dispersion, advection, reaction, transfer), time_step);
@@ -251,6 +250,24 @@ class Reactor_PFR extends Reactor {
     return stateNew;
   }
 
+  /**
+   * 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.