Merge pull request 'Final bug fixes and documentation' (#6) from dev-Pieter into main

Reviewed-on: #6

reactor.html

@@ -102,6 +102,19 @@
             } else {
                 $(".PFR").show();
             }
+
+            const updateDx = () => {
+                const length = parseFloat($("#node-input-length").val()) || 0;
+                const resolution = parseFloat($("#node-input-resolution_L").val()) || 1;
+                const dx = resolution > 0 ? (length / resolution).toFixed(6) : "N/A";
+                $("#dx-output").text(dx + " m");
+            };
+
+            // Set up event listeners for real-time updates
+            $("#node-input-length, #node-input-resolution_L").on("change keyup", updateDx);
+            
+            // Initial calculation
+            updateDx();
         },
         oneditsave: function() {
             // save logger fields
@@ -144,6 +157,10 @@
         <label for="node-input-resolution_L"><i class="fa fa-tag"></i> Resolution</label>
         <input type="text" id="node-input-resolution_L" placeholder="#">
     </div>
+    <div class="form-row PFR">
+        <label for="node-input-dx"><i class="fa fa-tag"></i> dx (length / resolution) [m]</label>
+        <span id="dx-output" style="display: inline-block; padding: 8px; font-weight: bold; color: #333;">--</span>
+    </div>
     <h3> Internal mass transfer calculation (optional) </h3>
     <div class="form-row">
         <label for="node-input-kla"><i class="fa fa-tag"></i> kLa [d-1]</label>

src/nodeClass.js

@@ -34,7 +34,6 @@ class nodeClass {
             switch (msg.topic) {
                 case "clock":
                     this.source.updateState(msg.timestamp);
-                    send([msg, null, null]);
                     break;
                 case "Fluent":
                     this.source.setInfluent = msg;
@@ -42,9 +41,6 @@ class nodeClass {
                 case "OTR":
                     this.source.setOTR = msg;
                     break;
-                case "Temperature":
-                    this.source.setTemperature = msg;
-                    break;
                 case "Dispersion":
                     this.source.setDispersion = msg;
                     break;

src/reaction_modules/asm3_class Koch.js

@@ -7,17 +7,7 @@ const ASM_CONSTANTS = {
   NUM_SPECIES: 13
 };
 
-/**
+const KINETIC_CONSTANTS = {
- * 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] 
@@ -42,13 +32,9 @@ class ASM3 {
   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]
-    };
+};
 
-    /**
+const STOICHIOMETRIC_CONSTANTS = {
-     * 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]
@@ -75,7 +61,25 @@ class ASM3 {
   // 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]
-    };
+};
+
+/**
+ * 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 = KINETIC_CONSTANTS;
+
+    /**
+     * Stoichiometric and composition parameters for ASM3. Using Koch et al. 2000 parameters.
+     * @property {Object} stoi_params - Stoichiometric parameters
+     */
+    this.stoi_params = STOICHIOMETRIC_CONSTANTS;
 
     /**
      * Temperature theta parameters for ASM3. Using Koch et al. 2000 parameters.
@@ -215,4 +219,4 @@ class ASM3 {
   }
 }
 
-module.exports = { ASM3, ASM_CONSTANTS };
+module.exports = { ASM3, ASM_CONSTANTS, KINETIC_CONSTANTS, STOICHIOMETRIC_CONSTANTS };

src/reaction_modules/asm3_class.js

@@ -7,17 +7,7 @@ const ASM_CONSTANTS = {
   NUM_SPECIES: 13
 };
 
-/**
+const KINETIC_CONSTANTS = {
- * 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] 
@@ -42,13 +32,9 @@ class ASM3 {
   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]
-    };
+};
 
-    /**
+const STOICHIOMETRIC_CONSTANTS = {
-     * 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]
@@ -75,7 +61,25 @@ class ASM3 {
   // 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]
-    };
+};
+
+/**
+ * ASM3 class for the Activated Sludge Model No. 3 (ASM3).
+ */
+class ASM3 {
+
+  constructor() {
+    /**
+     * Kinetic parameters for ASM3 at 20 C.
+     * @property {Object} kin_params - Kinetic parameters
+     */
+    this.kin_params = KINETIC_CONSTANTS;
+
+    /**
+     * Stoichiometric and composition parameters for ASM3.
+     * @property {Object} stoi_params - Stoichiometric parameters
+     */
+    this.stoi_params = STOICHIOMETRIC_CONSTANTS;
 
     /**
      * Temperature theta parameters for ASM3.
@@ -215,4 +219,4 @@ class ASM3 {
   }
 }
 
-module.exports = { ASM3, ASM_CONSTANTS };
+module.exports = { ASM3, ASM_CONSTANTS, KINETIC_CONSTANTS, STOICHIOMETRIC_CONSTANTS };

src/specificClass.js

@@ -10,9 +10,9 @@ const mathConfig = {
 
 const math = create(all, mathConfig);
 
-const BC_PADDING = 2;
+const BC_PADDING = 2; // Boundary Condition padding for open boundaries in extendedState variable
 const DEBUG = false;
-const DAY2MS = 1000 * 60 * 60 * 24;
+const DAY2MS = 1000 * 60 * 60 * 24; // convert between days and milliseconds
 
 class Reactor {
   /**
@@ -25,13 +25,14 @@ class Reactor {
     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.childRegistrationUtils = new childRegistrationUtils(this); // Child registration utility
+    this.childRegistrationUtils = new childRegistrationUtils(this); // child registration utility
 
+    // placeholder variables for children and parents
     this.upstreamReactor = null;
     this.downstreamReactor = null;
     this.returnPump = null;
 
-    this.asm = new ASM3();
+    this.asm = new ASM3(); // Reaction model
 
     this.volume = config.volume; // fluid volume reactor [m3]
 
@@ -42,9 +43,9 @@ class Reactor {
 
     this.kla = config.kla; // if NaN, use externaly provided OTR [d-1]
 
-    this.currentTime = Date.now(); // milliseconds since epoch [ms]
+    this.currentTime = null; // milliseconds since epoch [ms]
     this.timeStep = 1 / (24*60*60) * this.config.timeStep; // time step in seconds, converted to days.
-    this.speedUpFactor = 100; // speed up factor for simulation, 60 means 1 minute per simulated second
+    this.speedUpFactor = 1; // speed up factor for simulation, 60 means 1 minute per simulated second
   }
 
   /**
@@ -113,6 +114,11 @@ class Reactor {
     }
   }
 
+  /**
+   * Register child function required for child registration.
+   * @param {object} child 
+   * @param {string} softwareType 
+   */
   registerChild(child, softwareType) {
     if(!child) {
       this.logger.error(`Invalid ${softwareType} child provided.`);
@@ -161,18 +167,14 @@ class Reactor {
 
   _connectReactor(reactorChild) {
     if (reactorChild.config.functionality.positionVsParent != "upstream") {
-      this.logger.warn("Reactor children of reactors should always be upstream.");
+      this.logger.warn("Reactor children of other reactors should always be upstream!");
-    }
-
-    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.");
     }
 
     // set upstream and downstream reactor variable in current and child nodes respectively for easy access
     this.upstreamReactor = reactorChild;
     reactorChild.downstreamReactor = this;
 
-    reactorChild.emitter.on("stateChange", (eventData) => {
+    reactorChild.emitter.on("stateChange", (eventData) => { // Triggers state update in downstream reactor.
       this.logger.debug(`State change of upstream reactor detected.`);
       this.updateState(eventData);
     });
@@ -203,20 +205,32 @@ class Reactor {
    * Update the reactor state based on the new time.
    * @param {number} newTime - New time to update reactor state to, in milliseconds since epoch.
    */
-  updateState(newTime) { // expect update with timestamp
+  updateState(newTime) {
-    if (this.upstreamReactor) {
+    if (!this.currentTime) { // initialise currentTime variable
-      this.setInfluent = this.upstreamReactor.getEffluent[0]; // grab main effluent upstream reactor
+      this.currentTime = newTime;
+      return;
+    }
+
+    if (this.upstreamReactor) { // grab main effluent upstream reactor
+      this.setInfluent = this.upstreamReactor.getEffluent[0];
     }
 
     const n_iter = Math.floor(this.speedUpFactor * (newTime-this.currentTime) / (this.timeStep*DAY2MS));
-    if (n_iter) {
+
+    if (n_iter == 0) { // no update required, change in currentTime smaller than time step
+      return;
+    }
+
     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);
+    this.emitter.emit("stateChange", this.currentTime); // update downstream reactors
+
+    if (this.returnPump) { // update recirculation pump state
+      this.returnPump.updateSourceSink();
     }
   }
 }
@@ -231,6 +245,23 @@ class Reactor_CSTR extends Reactor {
       this.state = config.initialState;
   }
 
+  _updateMeasurement(measurementType, value, position, context) {
+
+    switch(measurementType) {
+      case "quantity (oxygen)":
+        this.state[ASM_CONSTANTS.S_O_INDEX] = value;
+        break;
+      case "quantity (ammonium)":
+        this.state[ASM_CONSTANTS.S_NH_INDEX] = value;
+        break;
+      case "quantity (nox)":
+        this.state[ASM_CONSTANTS.S_NO_INDEX] = value;
+        break;
+      default:
+        super._updateMeasurement(measurementType, value, position, context);
+    }
+  }
+
   /**
    * Tick the reactor state using the forward Euler method.
    * @param {number} time_step - Time step for the simulation [d].
@@ -241,7 +272,7 @@ class Reactor_CSTR extends Reactor {
     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(ASM_CONSTANTS.NUM_SPECIES).fill(0.0);
-    transfer[ASM_CONSTANTS.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
+    transfer[ASM_CONSTANTS.S_O_INDEX] = isNaN(this.kla) ? this.OTR : this._calcOTR(this.state[ASM_CONSTANTS.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
@@ -290,13 +321,23 @@ class Reactor_PFR extends Reactor {
     this.D = this._constrainDispersion(input.payload);
   }
 
+  _connectReactor(reactorChild) {
+    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.");
+    }
+    super._connectReactor(reactorChild);
+  }
+
+  /**
+   * Update the reactor state based on the new time. Performs checks specific to PFR.
+   * @param {number} newTime - New time to update reactor state to, in milliseconds since epoch.
+   */
   updateState(newTime) {
     super.updateState(newTime);
-    // let Pe_local = this.d_x*math.sum(this.Fs)/(this.D*this.A)
+    
-    this.D = this._constrainDispersion(this.D);
+    this.D = this._constrainDispersion(this.D); // constrains D to minimum dispersion, so that local Péclet number is always above 2
     const 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) {
@@ -375,8 +416,8 @@ class Reactor_PFR extends Reactor {
    */
   _applyBoundaryConditions() {
     // Upstream BC
-    if (this.upstreamReactor) {
+    if (this.upstreamReactor && this.upstreamReactor.config.reactor_type == "PFR") {
-      // Open boundary
+      // Open boundary, if upstream reactor is PFR
       this.extendedState.splice(0, BC_PADDING, ...this.upstreamReactor.state.slice(-BC_PADDING)); 
     } else {
       if (math.sum(this.Fs) > 0) {
@@ -384,7 +425,7 @@ class Reactor_PFR extends Reactor {
         const BC_C_in = math.multiply(1 / math.sum(this.Fs), [this.Fs], this.Cs_in)[0];
         const BC_dispersion_term = 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])));
-        // Numerical boundary condition
+        // Numerical boundary condition (first-order accurate)
         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 {
         // Neumann BC (no flux)
@@ -393,8 +434,8 @@ class Reactor_PFR extends Reactor {
     }
 
     // Downstream BC
-    if (this.downstreamReactor) {
+    if (this.downstreamReactor && this.downstreamReactor.config.reactor_type == "PFR") {
-      // Open boundary
+      // Open boundary, if downstream reactor is PFR
       this.extendedState.splice(this.n_x+BC_PADDING, BC_PADDING, ...this.downstreamReactor.state.slice(0, BC_PADDING));
     } else {
       // Neumann BC (no flux)
@@ -404,7 +445,6 @@ class Reactor_PFR extends Reactor {
 
   /**
    * Create finite difference first derivative operator.
-   * @returns {Array} - First derivative operator matrix.
    */
   _makeDoperator() { // create gradient operator
     const D_size = this.n_x+2*BC_PADDING;
@@ -420,7 +460,6 @@ class Reactor_PFR extends Reactor {
 
   /**
    * Create central finite difference second derivative operator.
-   * @returns {Array} - Second derivative operator matrix.
    */
   _makeD2operator() { // create the central second derivative operator
     const D_size = this.n_x+2*BC_PADDING;
@@ -433,6 +472,9 @@ class Reactor_PFR extends Reactor {
     return D2;
   }
 
+  /**
+   * Constrains dispersion so that local Péclet number stays below 2. Needed for stable central differencing method.
+   */
   _constrainDispersion(D) {
     const Dmin = math.sum(this.Fs) * this.d_x / (1.999 * this.A);
     if (D < Dmin) {