Merge branch 'boundary-conditions'

package-lock.json

@@ -1,22 +1,22 @@
 {
-  "name": "asm3",
+  "name": "reactor",
   "version": "0.0.1",
   "lockfileVersion": 3,
   "requires": true,
   "packages": {
     "": {
-      "name": "asm3",
+      "name": "reactor",
       "version": "0.0.1",
       "license": "SEE LICENSE",
       "dependencies": {
-        "generalFunctions": "git+https://gitea.centraal.wbd-rd.nl/p.vanderwilt/generalFunctions.git#fix-missing-references",
+        "generalFunctions": "git+https://gitea.centraal.wbd-rd.nl/RnD/generalFunctions.git",
         "mathjs": "^14.5.2"
       }
     },
     "node_modules/@babel/runtime": {
-      "version": "7.28.3",
+      "version": "7.28.4",
-      "resolved": "https://registry.npmjs.org/@babel/runtime/-/runtime-7.28.3.tgz",
+      "resolved": "https://registry.npmjs.org/@babel/runtime/-/runtime-7.28.4.tgz",
-      "integrity": "sha512-9uIQ10o0WGdpP6GDhXcdOJPJuDgFtIDtN/9+ArJQ2NAfAmiuhTQdzkaTGR33v43GYS2UrSA0eX2pPPHoFVvpxA==",
+      "integrity": "sha512-Q/N6JNWvIvPnLDvjlE1OUBLPQHH6l3CltCEsHIujp45zQUSSh8K+gHnaEX45yAT1nyngnINhvWtzN+Nb9D8RAQ==",
       "license": "MIT",
       "engines": {
         "node": ">=6.9.0"
@@ -62,7 +62,7 @@
     },
     "node_modules/generalFunctions": {
       "version": "1.0.0",
-      "resolved": "git+https://gitea.centraal.wbd-rd.nl/p.vanderwilt/generalFunctions.git#302e12238745766a679ef11ca6ed5f4ea1548f87",
+      "resolved": "git+https://gitea.centraal.wbd-rd.nl/RnD/generalFunctions.git#efc97d6cd17399391b011298e47e8c1b1599592d",
       "license": "SEE LICENSE"
     },
     "node_modules/javascript-natural-sort": {
@@ -72,9 +72,9 @@
       "license": "MIT"
     },
     "node_modules/mathjs": {
-      "version": "14.7.0",
+      "version": "14.8.0",
-      "resolved": "https://registry.npmjs.org/mathjs/-/mathjs-14.7.0.tgz",
+      "resolved": "https://registry.npmjs.org/mathjs/-/mathjs-14.8.0.tgz",
-      "integrity": "sha512-RaMhb+9MSESjDZNox/FzzuFpIUI+oxGLyOy1t3BMoW53pGWnTzZtlucJ5cvbit0dIMYlCq00gNbW1giZX4/1Rg==",
+      "integrity": "sha512-DN4wmAjNzFVJ9vHqpAJ3vX0UF306u/1DgGKh7iVPuAFH19JDRd9NAaQS764MsKbSwDB6uBSkQEmgVmKdgYaCoQ==",
       "license": "Apache-2.0",
       "dependencies": {
         "@babel/runtime": "^7.26.10",

package.json

@@ -1,5 +1,5 @@
 {
-  "name": "asm3",
+  "name": "reactor",
   "version": "0.0.1",
   "description": "Implementation of the asm3 model for Node-Red",
   "repository": {
@@ -27,7 +27,7 @@
     }
   },
   "dependencies": {
-    "generalFunctions": "git+https://gitea.centraal.wbd-rd.nl/p.vanderwilt/generalFunctions.git#fix-missing-references",
+    "generalFunctions": "git+https://gitea.centraal.wbd-rd.nl/RnD/generalFunctions.git",
     "mathjs": "^14.5.2"
   }
 }

reactor.html

@@ -1,7 +1,7 @@
-<script src="/advancedReactor/menu.js"></script>
+<script src="/reactor/menu.js"></script>
 
 <script type="text/javascript">
-    RED.nodes.registerType("advancedReactor", {
+    RED.nodes.registerType("reactor", {
         category: "WWTP",
         color: "#c4cce0",
         defaults: {
@@ -41,13 +41,13 @@
         outputLabels: ["process", "dbase", "parent"],
         icon: "font-awesome/fa-recycle",
         label: function() {
-            return this.name || "advancedReactor";
+            return this.name || "Reactor";
         },
         oneditprepare: function() {
             // wait for the menu scripts to load
             const waitForMenuData = () => {
-                if (window.EVOLV?.nodes?.advancedReactor?.initEditor) {
+                if (window.EVOLV?.nodes?.reactor?.initEditor) {
-                    window.EVOLV.nodes.advancedReactor.initEditor(this);
+                    window.EVOLV.nodes.reactor.initEditor(this);
                 } else {
                     setTimeout(waitForMenuData, 50);
                 }
@@ -115,8 +115,8 @@
         },
         oneditsave: function() {
             // save logger fields
-            if (window.EVOLV?.nodes?.['advancedReactor']?.loggerMenu?.saveEditor) {
+            if (window.EVOLV?.nodes?.reactor?.loggerMenu?.saveEditor) {
-                window.EVOLV.nodes['advancedReactor'].loggerMenu.saveEditor(this);
+                window.EVOLV.nodes.reactor.loggerMenu.saveEditor(this);
             }
 
             // save position field
@@ -136,7 +136,7 @@
     });
 </script>
 
-<script type="text/html" data-template-name="advancedReactor">
+<script type="text/html" data-template-name="reactor">
     <div class="form-row">
         <label for="node-input-name"><i class="fa fa-tag"></i> Name</label>
         <input type="text" id="node-input-name" placeholder="Name">
@@ -243,6 +243,6 @@
 
 </script>
 
-<script type="text/html" data-help-name="advancedReactor">
+<script type="text/html" data-help-name="reactor">
     <p>New reactor node</p>
 </script>

reactor.js

@@ -1,4 +1,4 @@
-const nameOfNode = "advancedReactor"; // name of the node, should match file name and node type in Node-RED
+const nameOfNode = "reactor"; // name of the node, should match file name and node type in Node-RED
 const nodeClass = require('./src/nodeClass.js'); // node class
 const { MenuManager } = require('generalFunctions');
 

src/specificClass.js

@@ -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();
 
@@ -40,8 +45,8 @@ class Reactor {
     this.kla = config.kla; // if NaN, use externaly provided OTR [d-1]
 
     this.currentTime = Date.now(); // milliseconds since epoch [ms]
-    this.timeStep = 1 / (24*60*60) * this.config.timeStep; // time step [d]
+    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,19 +118,18 @@ class Reactor {
     }
   }
 
-  _connectMeasurement(measurement) {
+  _connectMeasurement(measurementChild) {
-    if (!measurement) {
+    if (!measurementChild) {
       this.logger.warn("Invalid measurement provided.");
       return;
     }
 
-    const position = measurement.config.functionality.positionVsParent;
+    const position = measurementChild.config.functionality.positionVsParent;
-    const measurementType = measurement.config.asset.type;
+    const measurementType = measurementChild.config.asset.type;
-    const key = `${measurementType}_${position}`;
     const eventName = `${measurementType}.measured.${position}`;
 
     // Register event listener for measurement updates
-    this.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
@@ -140,15 +144,24 @@ class Reactor {
   }
 
 
-  _connectReactor(reactor) {
+  _connectReactor(reactorChild) {
-    if (!reactor) {
+    if (!reactorChild) {
       this.logger.warn("Invalid reactor provided.");
       return;
     }
 
-    this.upstreamReactor = reactor;
+    if (reactorChild.functionality.positionVsParent != "upstream") {
+      this.logger.warn("Reactor children of reactors should always be upstream.");
+    }
 
-    reactor.emitter.on("stateChange", (data) => {
+    if (math.abs(reactorChild.d_x - this.d_x) / this.d_x < 0.025) {
+      this.logger.warn("Significant grid sizing discrepancies between adjacent reactors! Change resolutions to match reactors grid step, or implement boundary value interpolation.");
+    }
+
+    this.upstreamReactor = reactorChild;
+    reactorChild.downstreamReactor = this;
+
+    reactorChild.emitter.on("stateChange", (data) => {
       this.logger.debug(`State change of upstream reactor detected.`);
       this.updateState(data);
     });
@@ -241,11 +254,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();
@@ -257,15 +274,16 @@ class Reactor_PFR extends Reactor {
    * @param {object} input - Input object (msg) containing payload with dispersion value [m2 d-1].
    */
   set setDispersion(input) {
-    this.D = input.payload;
+    this.D = this._constrainDispersion(input.payload);
   }
 
   updateState(newTime) {
     super.updateState(newTime);
-    let Pe_local = this.d_x*math.sum(this.Fs)/(this.D*this.A)
+    // let Pe_local = this.d_x*math.sum(this.Fs)/(this.D*this.A)
+    this.D = this._constrainDispersion(this.D);
     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.`);
+    // (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) {
@@ -283,25 +301,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);
+    this._applyBoundaryConditions();
-    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 dispersion = math.multiply(this.D / (this.d_x*this.d_x), this.D2_op, this.extendedState);
-    const transfer = Array.from(Array(this.n_x), () => new Array(NUM_SPECIES).fill(0));
+    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++) {
+    for (let i = BC_PADDING+1; i < BC_PADDING+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);
+        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");
@@ -312,74 +331,70 @@ 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) {
     switch(measurementType) {
-      case "oxygen":
+      case "quantity (oxygen)":
-        grid_pos = Math.round(position * this.n_x);
+        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);
   }
 
   /**
    * 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) {
+  _applyBoundaryConditions() {
-    if (math.sum(this.Fs) > 0) { // Danckwerts BC
+    if (this.upstreamReactor) {
-      const BC_C_in = math.multiply(1 / math.sum(this.Fs), [this.Fs], this.Cs_in)[0];
+      for (let i = 0; i < BC_PADDING; i++) {
-      const BC_dispersion_term = (1-this.alpha)*this.D*this.A/(math.sum(this.Fs)*this.d_x);
+        this.extendedState[i] = this.upstreamReactor.state.at(i-BC_PADDING);
-      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];
+      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
+  _makeDoperator() { // create gradient operator
-    if (higher_order) {
+    const D_size = this.n_x+2*BC_PADDING;
-      if (central) {
+    const I = math.resize(math.diag(Array(D_size).fill(1/12), -2), [D_size, D_size]);
-        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(D_size).fill(-2/3), -1), [D_size, D_size]);
-        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(D_size).fill(2/3), 1), [D_size, D_size]);
-        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(D_size).fill(-1/12), 2), [D_size, D_size]);
-        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 D = math.add(I, A, B, C);
+    // set by BCs elsewhere
-        const NearBoundary = Array(this.n_x).fill(0.0);
+    D.forEach((row, i) => i < BC_PADDING || i >= this.n_x+BC_PADDING ? row.fill(0) : row);
-        NearBoundary[0] = -1/4;
+    return D;
-        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;
-    }
   }
 
   /**
@@ -387,14 +402,24 @@ 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 D_size = this.n_x+2*BC_PADDING;
-    const A = math.resize(math.diag(Array(this.n_x).fill(1), 1), [this.n_x, this.n_x]);
+    const I = math.diag(Array(D_size).fill(-2), 0);
-    const B = math.resize(math.diag(Array(this.n_x).fill(1), -1), [this.n_x, this.n_x]);
+    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
+    // set by BCs elsewhere
-    D2[this.n_x - 1] = Array(this.n_x).fill(0);
+    D2.forEach((row, i) => i < BC_PADDING || i >= this.n_x+BC_PADDING ? row.fill(0) : row);
     return D2;
   }
+
+  _constrainDispersion(D) {
+    const Dmin = math.sum(this.Fs) * this.d_x / (1.999 * this.A);
+    if (D < Dmin) {
+      this.logger.warn(`Local Péclet number too high! Constraining given dispersion (${D}) to minimal value (${Dmin}).`);
+      return Dmin;
+    }
+    return D;
+  }
 }
 
 module.exports = { Reactor_CSTR, Reactor_PFR };