Enhance comments and documentation in Reactor classes for clarity and maintainability

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/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.");
-    }
-
-    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.logger.warn("Reactor children of other reactors should always be upstream!");
     }
 
     // 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
-    if (this.upstreamReactor) {
-      this.setInfluent = this.upstreamReactor.getEffluent[0]; // grab main effluent upstream reactor
+  updateState(newTime) {
+    if (!this.currentTime) { // initialise currentTime variable
+      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) {
-      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);
+
+    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); // 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) {
-      // Open boundary
+    if (this.upstreamReactor && this.upstreamReactor.config.reactor_type == "PFR") {
+      // 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) {
-      // Open boundary
+    if (this.downstreamReactor && this.downstreamReactor.config.reactor_type == "PFR") {
+      // 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) {