Enhance comments and documentation in Reactor classes for clarity and maintainability

Reviewed-on: #4

additional_nodes/recirculation-pump.html

@@ -1,57 +0,0 @@
-<script type="text/javascript">
-    RED.nodes.registerType("recirculation-pump", {
-        category: "WWTP",
-        color: "#e4a363",
-        defaults: {
-            name: { value: "" },
-            F2: { value: 0, required: true },
-            inlet: { value: 1, required: true }
-        },
-        inputs: 1,
-        outputs: 2,
-        outputLabels: ["Main effluent", "Recirculation effluent"],
-        icon: "font-awesome/fa-random",
-        label: function() {
-            return this.name || "Recirculation pump";
-        },
-        oneditprepare: function() {
-            $("#node-input-F2").typedInput({
-                type:"num",
-                types:["num"]
-            });
-            $("#node-input-inlet").typedInput({
-                type:"num",
-                types:["num"]
-            });
-        },
-        oneditsave: function() {
-            let debit = parseFloat($("#node-input-F2").typedInput("value"));
-            if (isNaN(debit) || debit < 0) {
-                RED.notify("Debit is not set correctly", {type: "error"});
-            }
-            let inlet = parseInt($("#node-input-n_inlets").typedInput("value"));
-            if (inlet < 1) {
-                RED.notify("Number of inlets not set correctly", {type: "error"});
-            }
-        }
-    });
-</script>
-
-<script type="text/html" data-template-name="recirculation-pump">
-    <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">
-    </div>
-    <div class="form-row">
-        <label for="node-input-F2"><i class="fa fa-tag"></i> Recirculation debit [m3 d-1]</label>
-        <input type="text" id="node-input-F2" placeholder="m3 s-1">
-    </div>
-    <div class="form-row">
-        <label for="node-input-inlet"><i class="fa fa-tag"></i> Assigned inlet recirculation</label>
-        <input type="text" id="node-input-inlet" placeholder="#">
-    </div>
-</script>
-
-<script type="text/html" data-help-name="recirculation-pump">
-    <p>Recirculation-pump for splitting streams</p>
-</script>

additional_nodes/recirculation-pump.js

@@ -1,40 +0,0 @@
-module.exports = function(RED) {
-    function recirculation(config) {
-        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":
-                    // conserve volume flow debit
-                    let F_in = msg.payload.F;
-                    let F1 = Math.max(F_in - F2, 0);
-                    let F2_corr = F_in < F2 ? F_in : F2;
-
-                    let msg_F1 = structuredClone(msg);
-                    msg_F1.payload.F = F1;
-
-                    let msg_F2 = {...msg};
-                    msg_F2.payload.F = F2_corr;
-                    msg_F2.payload.inlet = inlet_F2;
-
-                    send([msg_F1, msg_F2]);
-                    break;
-                case "clock":
-                    break;
-                default:
-                    console.log("Unknown topic: " + msg.topic);
-            }
-
-            if (done) {
-                done();
-            }
-        });
-
-    }
-    RED.nodes.registerType("recirculation-pump", recirculation);
-};

additional_nodes/settling-basin.html

@@ -1,57 +0,0 @@
-<script type="text/javascript">
-    RED.nodes.registerType("settling-basin", {
-        category: "WWTP",
-        color: "#e4a363",
-        defaults: {
-            name: { value: "" },
-            TS_set: { value: 0.1, required: true },
-            inlet: { value: 1, required: true }
-        },
-        inputs: 1,
-        outputs: 2,
-        outputLabels: ["Main effluent", "Sludge effluent"],
-        icon: "font-awesome/fa-random",
-        label: function() {
-            return this.name || "Settling basin";
-        },
-        oneditprepare: function() {
-            $("#node-input-TS_set").typedInput({
-                type:"num",
-                types:["num"]
-            });
-            $("#node-input-inlet").typedInput({
-                type:"num",
-                types:["num"]
-            });
-        },
-        oneditsave: function() {
-            let TS_set = parseFloat($("#node-input-TS_set").typedInput("value"));
-            if (isNaN(TS_set) || TS_set < 0) {
-                RED.notify("TS is not set correctly", {type: "error"});
-            }
-            let inlet = parseInt($("#node-input-n_inlets").typedInput("value"));
-            if (inlet < 1) {
-                RED.notify("Number of inlets not set correctly", {type: "error"});
-            }
-        }
-    });
-</script>
-
-<script type="text/html" data-template-name="settling-basin">
-    <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">
-    </div>
-    <div class="form-row">
-        <label for="node-input-TS_set"><i class="fa fa-tag"></i> Total Solids set point [g m-3]</label>
-        <input type="text" id="node-input-TS_set" placeholder="">
-    </div>
-    <div class="form-row">
-        <label for="node-input-inlet"><i class="fa fa-tag"></i> Assigned inlet return line</label>
-        <input type="text" id="node-input-inlet" placeholder="#">
-    </div>
-</script>
-
-<script type="text/html" data-help-name="settling-basin">
-    <p>Settling tank</p>
-</script>

additional_nodes/settling-basin.js

@@ -1,57 +0,0 @@
-module.exports = function(RED) {
-    function settler(config) {
-        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":
-                    // conserve volume flow debit
-                    let F_in = msg.payload.F;
-                    let C_in = msg.payload.C;
-                    let F2 = (F_in * C_in[12]) / TS_set;
-
-                    let F1 = Math.max(F_in - F2, 0);
-                    let F2_corr = F_in < F2 ? F_in : F2;
-
-                    let msg_F1 = structuredClone(msg);
-                    msg_F1.payload.F = F1;
-                    msg_F1.payload.C[7] = 0;
-                    msg_F1.payload.C[8] = 0;
-                    msg_F1.payload.C[9] = 0;
-                    msg_F1.payload.C[10] = 0;
-                    msg_F1.payload.C[11] = 0;
-                    msg_F1.payload.C[12] = 0;
-
-                    let msg_F2 = {...msg};
-                    msg_F2.payload.F = F2_corr;
-                    if (F2_corr > 0) {
-                        msg_F2.payload.C[7] = F_in * C_in[7] / F2;
-                        msg_F2.payload.C[8] = F_in * C_in[8] / F2;
-                        msg_F2.payload.C[9] = F_in * C_in[9] / F2;
-                        msg_F2.payload.C[10] = F_in * C_in[10] / F2;
-                        msg_F2.payload.C[11] = F_in * C_in[11] / F2;
-                        msg_F2.payload.C[12] = F_in * C_in[12] / F2;
-                    }
-                    msg_F2.payload.inlet = inlet_sludge;
-
-                    send([msg_F1, msg_F2]);
-                    break;
-                case "clock":
-                    break;
-                default:
-                    console.log("Unknown topic: " + msg.topic);
-            }
-
-            if (done) {
-                done();
-            }
-        });
-
-    }
-    RED.nodes.registerType("settling-basin", settler);
-};

package-lock.json

@@ -1,119 +0,0 @@
-{
-  "name": "reactor",
-  "version": "0.0.1",
-  "lockfileVersion": 3,
-  "requires": true,
-  "packages": {
-    "": {
-      "name": "reactor",
-      "version": "0.0.1",
-      "license": "SEE LICENSE",
-      "dependencies": {
-        "generalFunctions": "git+https://gitea.centraal.wbd-rd.nl/RnD/generalFunctions.git",
-        "mathjs": "^14.5.2"
-      }
-    },
-    "node_modules/@babel/runtime": {
-      "version": "7.28.4",
-      "resolved": "https://registry.npmjs.org/@babel/runtime/-/runtime-7.28.4.tgz",
-      "integrity": "sha512-Q/N6JNWvIvPnLDvjlE1OUBLPQHH6l3CltCEsHIujp45zQUSSh8K+gHnaEX45yAT1nyngnINhvWtzN+Nb9D8RAQ==",
-      "license": "MIT",
-      "engines": {
-        "node": ">=6.9.0"
-      }
-    },
-    "node_modules/complex.js": {
-      "version": "2.4.2",
-      "resolved": "https://registry.npmjs.org/complex.js/-/complex.js-2.4.2.tgz",
-      "integrity": "sha512-qtx7HRhPGSCBtGiST4/WGHuW+zeaND/6Ld+db6PbrulIB1i2Ev/2UPiqcmpQNPSyfBKraC0EOvOKCB5dGZKt3g==",
-      "license": "MIT",
-      "engines": {
-        "node": "*"
-      },
-      "funding": {
-        "type": "github",
-        "url": "https://github.com/sponsors/rawify"
-      }
-    },
-    "node_modules/decimal.js": {
-      "version": "10.6.0",
-      "resolved": "https://registry.npmjs.org/decimal.js/-/decimal.js-10.6.0.tgz",
-      "integrity": "sha512-YpgQiITW3JXGntzdUmyUR1V812Hn8T1YVXhCu+wO3OpS4eU9l4YdD3qjyiKdV6mvV29zapkMeD390UVEf2lkUg==",
-      "license": "MIT"
-    },
-    "node_modules/escape-latex": {
-      "version": "1.2.0",
-      "resolved": "https://registry.npmjs.org/escape-latex/-/escape-latex-1.2.0.tgz",
-      "integrity": "sha512-nV5aVWW1K0wEiUIEdZ4erkGGH8mDxGyxSeqPzRNtWP7ataw+/olFObw7hujFWlVjNsaDFw5VZ5NzVSIqRgfTiw==",
-      "license": "MIT"
-    },
-    "node_modules/fraction.js": {
-      "version": "5.3.4",
-      "resolved": "https://registry.npmjs.org/fraction.js/-/fraction.js-5.3.4.tgz",
-      "integrity": "sha512-1X1NTtiJphryn/uLQz3whtY6jK3fTqoE3ohKs0tT+Ujr1W59oopxmoEh7Lu5p6vBaPbgoM0bzveAW4Qi5RyWDQ==",
-      "license": "MIT",
-      "engines": {
-        "node": "*"
-      },
-      "funding": {
-        "type": "github",
-        "url": "https://github.com/sponsors/rawify"
-      }
-    },
-    "node_modules/generalFunctions": {
-      "version": "1.0.0",
-      "resolved": "git+https://gitea.centraal.wbd-rd.nl/RnD/generalFunctions.git#efc97d6cd17399391b011298e47e8c1b1599592d",
-      "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",
-      "integrity": "sha512-nO6jcEfZWQXDhOiBtG2KvKyEptz7RVbpGP4vTD2hLBdmNQSsCiicO2Ioinv6UI4y9ukqnBpy+XZ9H6uLNgJTlw==",
-      "license": "MIT"
-    },
-    "node_modules/mathjs": {
-      "version": "14.8.0",
-      "resolved": "https://registry.npmjs.org/mathjs/-/mathjs-14.8.0.tgz",
-      "integrity": "sha512-DN4wmAjNzFVJ9vHqpAJ3vX0UF306u/1DgGKh7iVPuAFH19JDRd9NAaQS764MsKbSwDB6uBSkQEmgVmKdgYaCoQ==",
-      "license": "Apache-2.0",
-      "dependencies": {
-        "@babel/runtime": "^7.26.10",
-        "complex.js": "^2.2.5",
-        "decimal.js": "^10.4.3",
-        "escape-latex": "^1.2.0",
-        "fraction.js": "^5.2.1",
-        "javascript-natural-sort": "^0.7.1",
-        "seedrandom": "^3.0.5",
-        "tiny-emitter": "^2.1.0",
-        "typed-function": "^4.2.1"
-      },
-      "bin": {
-        "mathjs": "bin/cli.js"
-      },
-      "engines": {
-        "node": ">= 18"
-      }
-    },
-    "node_modules/seedrandom": {
-      "version": "3.0.5",
-      "resolved": "https://registry.npmjs.org/seedrandom/-/seedrandom-3.0.5.tgz",
-      "integrity": "sha512-8OwmbklUNzwezjGInmZ+2clQmExQPvomqjL7LFqOYqtmuxRgQYqOD3mHaU+MvZn5FLUeVxVfQjwLZW/n/JFuqg==",
-      "license": "MIT"
-    },
-    "node_modules/tiny-emitter": {
-      "version": "2.1.0",
-      "resolved": "https://registry.npmjs.org/tiny-emitter/-/tiny-emitter-2.1.0.tgz",
-      "integrity": "sha512-NB6Dk1A9xgQPMoGqC5CVXn123gWyte215ONT5Pp5a0yt4nlEoO1ZWeCwpncaekPHXO60i47ihFnZPiRPjRMq4Q==",
-      "license": "MIT"
-    },
-    "node_modules/typed-function": {
-      "version": "4.2.1",
-      "resolved": "https://registry.npmjs.org/typed-function/-/typed-function-4.2.1.tgz",
-      "integrity": "sha512-EGjWssW7Tsk4DGfE+5yluuljS1OGYWiI1J6e8puZz9nTMM51Oug8CD5Zo4gWMsOhq5BI+1bF+rWTm4Vbj3ivRA==",
-      "license": "MIT",
-      "engines": {
-        "node": ">= 18"
-      }
-    }
-  }
-}

package.json

@@ -11,6 +11,7 @@
     "activated sludge",
     "wastewater",
     "biological model",
+    "EVOLV",
     "node-red"
   ],
   "license": "SEE LICENSE",
@@ -21,9 +22,7 @@
   },
   "node-red": {
     "nodes": {
-      "reactor": "reactor.js",
-      "recirculation-pump": "additional_nodes/recirculation-pump.js",
-      "settling-basin": "additional_nodes/settling-basin.js"
+      "reactor": "reactor.js"
     }
   },
   "dependencies": {

reactor.html

@@ -2,7 +2,7 @@
 
 <script type="text/javascript">
     RED.nodes.registerType("reactor", {
-        category: "WWTP",
+        category: "EVOLV",
         color: "#c4cce0",
         defaults: {
             name: { value: "" },
@@ -10,8 +10,6 @@
             volume: { value: 0., required: true },
             length: { value: 0.},
             resolution_L: { value: 0.},
-            alpha: {value: 0},
-            n_inlets: { value: 1, required: true},
             kla: { value: null },
             
             S_O_init: { value: 0., required: true },
@@ -33,7 +31,7 @@
             enableLog: { value: false },
             logLevel: { value: "error" },
 
-            positionVsParent: { value: "" },
+            positionVsParent: { value: "" }
         },
         inputs: 1,
         outputs: 3,
@@ -58,10 +56,6 @@
                 type:"num",
                 types:["num"]
             });
-            $("#node-input-n_inlets").typedInput({
-                type:"num",
-                types:["num"]
-            });
             $("#node-input-length").typedInput({
                 type:"num",
                 types:["num"]
@@ -97,10 +91,6 @@
                     $(".PFR").show();
                 }
             });
-            $("#node-input-alpha").typedInput({
-                type:"num",
-                types:["num"]
-            })
             $("#node-input-timeStep").typedInput({
                 type:"num",
                 types:["num"]
@@ -112,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
@@ -128,10 +131,6 @@
             if (isNaN(volume) || volume <= 0) {
                 RED.notify("Fluid volume not set correctly", {type: "error"});
             }
-            let n_inlets = parseInt($("#node-input-n_inlets").typedInput("value"));
-            if (isNaN(n_inlets) || n_inlets < 1) {
-                RED.notify("Number of inlets not set correctly", {type: "error"});
-            }
         }
     });
 </script>
@@ -158,16 +157,9 @@
         <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="PFR">
-        <p> Inlet boundary condition parameter &alpha; (&alpha; = 0: Danckwerts BC / &alpha; = 1: Dirichlet BC) </p>
-        <div class="form-row">
-            <label for="node-input-alpha"><i class="fa fa-tag"></i>Adjustable parameter BC</label>
-            <input type="text" id="node-input-alpha">
-        </div>
-    </div>
-    <div class="form-row">
-        <label for="node-input-n_inlets"><i class="fa fa-tag"></i> Number of inlets</label>
-        <input type="text" id="node-input-n_inlets" placeholder="#">
+    <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">
@@ -240,7 +232,6 @@
     <!-- Position fields will be injected here -->
     <div id="position-fields-placeholder"></div>
 
-
 </script>
 
 <script type="text/html" data-help-name="reactor">

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;
@@ -87,8 +83,6 @@ class nodeClass {
             volume: parseFloat(uiConfig.volume),
             length: parseFloat(uiConfig.length),
             resolution_L: parseInt(uiConfig.resolution_L),
-            alpha: parseFloat(uiConfig.alpha),
-            n_inlets: parseInt(uiConfig.n_inlets),
             kla: parseFloat(uiConfig.kla),
             initialState: [
                 parseFloat(uiConfig.S_O_init),

src/reaction_modules/asm3_class Koch.js

@@ -1,4 +1,67 @@
-const math = require('mathjs')
+const math = require('mathjs');
+
+const ASM_CONSTANTS = {
+  S_O_INDEX: 0,
+  S_NH_INDEX: 3,
+  S_NO_INDEX: 5,
+  NUM_SPECIES: 13
+};
+
+const KINETIC_CONSTANTS = {
+  // 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] 
+  // Heterotrophs
+  k_STO:    12.,  // storage rate constant          [g S_S g-1 X_H d-1]
+  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.,  // 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]
+  K_HCO:    0.1,  // saturation constant S_HCO      [mole HCO3 m-3]
+  b_H_O:    0.3,  // aerobic respiration rate       [d-1]
+  b_H_NO:   0.15, // anoxic respiration rate        [d-1]
+  b_STO_O:  0.3,  // aerobic respitation rate X_STO [d-1]
+  b_STO_NO: 0.15, // anoxic respitation rate X_STO  [d-1]
+  // Autotrophs
+  mu_A_max: 1.3,  // maximum specific growth rate   [d-1]
+  K_A_NH:   1.4,  // saturation constant S_NH3      [g NH3-N m-3] 
+  K_A_O:    0.5,  // saturation constant S_0        [g O2 m-3]
+  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 = {
+  // 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]
+  // Yields
+  Y_STO_O:  0.80, // aerobic yield X_STO per S_S    [g X_STO g-1 S_S]
+  Y_STO_NO: 0.70, // anoxic yield X_STO per S_S     [g X_STO g-1 S_S]
+  Y_H_O:    0.80, // aerobic yield X_H per X_STO    [g X_H g-1 X_STO]
+  Y_H_NO:   0.65, // anoxic yield X_H per X_STO     [g X_H g-1 X_STO]
+  Y_A:      0.24, // anoxic yield X_A per S_NO      [g X_A g-1 NO3-N]
+  // Composition (COD via DoR)
+  i_CODN:  -1.71, // COD content (DoR)              [g COD g-1 N2-N]
+  i_CODNO: -4.57, // COD content (DoR)              [g COD g-1 NO3-N]
+  // Composition (nitrogen)
+  i_NSI:    0.01, // nitrogen content S_I           [g N g-1 S_I]
+  i_NSS:    0.03, // nitrogen content S_S           [g N g-1 S_S]
+  i_NXI:    0.04, // nitrogen content X_I           [g N g-1 X_I]
+  i_NXS:    0.03, // nitrogen content X_S           [g N g-1 X_S]
+  i_NBM:    0.07, // nitrogen content X_H / X_A     [g N g-1 X_H / X_A]
+  // Composition (TSS)
+  i_TSXI:   0.75, // TSS content X_I                [g TS g-1 X_I]
+  i_TSXS:   0.75, // TSS content X_S                [g TS g-1 X_S]
+  i_TSBM:   0.90, // TSS content X_H / X_A          [g TS g-1 X_H / X_A]
+  i_TSSTO:  0.60, // TSS content X_STO (PHB based)  [g TS g-1 X_STO]
+  // 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.
@@ -10,65 +73,13 @@ class ASM3 {
      * 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] 
-      // Heterotrophs
-      k_STO:    12.,  // storage rate constant          [g S_S g-1 X_H d-1]
-      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.,  // 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]
-      K_HCO:    0.1,  // saturation constant S_HCO      [mole HCO3 m-3]
-      b_H_O:    0.3,  // aerobic respiration rate       [d-1]
-      b_H_NO:   0.15, // anoxic respiration rate        [d-1]
-      b_STO_O:  0.3,  // aerobic respitation rate X_STO [d-1]
-      b_STO_NO: 0.15, // anoxic respitation rate X_STO  [d-1]
-      // Autotrophs
-      mu_A_max: 1.3,  // maximum specific growth rate   [d-1]
-      K_A_NH:   1.4,  // saturation constant S_NH3      [g NH3-N m-3] 
-      K_A_O:    0.5,  // saturation constant S_0        [g O2 m-3]
-      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]
-    };
+    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 = {
-      // 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]
-      // Yields
-      Y_STO_O:  0.80, // aerobic yield X_STO per S_S    [g X_STO g-1 S_S]
-      Y_STO_NO: 0.70, // anoxic yield X_STO per S_S     [g X_STO g-1 S_S]
-      Y_H_O:    0.80, // aerobic yield X_H per X_STO    [g X_H g-1 X_STO]
-      Y_H_NO:   0.65, // anoxic yield X_H per X_STO     [g X_H g-1 X_STO]
-      Y_A:      0.24, // anoxic yield X_A per S_NO      [g X_A g-1 NO3-N]
-      // Composition (COD via DoR)
-      i_CODN:  -1.71, // COD content (DoR)              [g COD g-1 N2-N]
-      i_CODNO: -4.57, // COD content (DoR)              [g COD g-1 NO3-N]
-      // Composition (nitrogen)
-      i_NSI:    0.01, // nitrogen content S_I           [g N g-1 S_I]
-      i_NSS:    0.03, // nitrogen content S_S           [g N g-1 S_S]
-      i_NXI:    0.04, // nitrogen content X_I           [g N g-1 X_I]
-      i_NXS:    0.03, // nitrogen content X_S           [g N g-1 X_S]
-      i_NBM:    0.07, // nitrogen content X_H / X_A     [g N g-1 X_H / X_A]
-      // Composition (TSS)
-      i_TSXI:   0.75, // TSS content X_I                [g TS g-1 X_I]
-      i_TSXS:   0.75, // TSS content X_S                [g TS g-1 X_S]
-      i_TSBM:   0.90, // TSS content X_H / X_A          [g TS g-1 X_H / X_A]
-      i_TSSTO:  0.60, // TSS content X_STO (PHB based)  [g TS g-1 X_STO]
-      // 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]
-    };
+    this.stoi_params = STOICHIOMETRIC_CONSTANTS;
 
     /**
      * Temperature theta parameters for ASM3. Using Koch et al. 2000 parameters.
@@ -208,4 +219,4 @@ class ASM3 {
   }
 }
 
-module.exports = ASM3;
+module.exports = { ASM3, ASM_CONSTANTS, KINETIC_CONSTANTS, STOICHIOMETRIC_CONSTANTS };

src/reaction_modules/asm3_class.js

@@ -1,4 +1,67 @@
-const math = require('mathjs')
+const math = require('mathjs');
+
+const ASM_CONSTANTS = {
+  S_O_INDEX: 0,
+  S_NH_INDEX: 3,
+  S_NO_INDEX: 5,
+  NUM_SPECIES: 13
+};
+
+const KINETIC_CONSTANTS = {
+  // 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] 
+  // Heterotrophs
+  k_STO:    5.,   // storage rate constant          [g S_S g-1 X_H d-1]
+  nu_NO:    0.6,  // 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:      2.,   // saturation constant S_s        [g COD m-3]
+  K_STO:    1.,   // saturation constant X_STO      [g X_STO g-1 X_H]
+  mu_H_max: 2.,   // maximum specific growth rate   [d-1]
+  K_NH:     0.01, // saturation constant S_NH3      [g NH3-N m-3]
+  K_HCO:    0.1,  // saturation constant S_HCO      [mole HCO3 m-3]
+  b_H_O:    0.2,  // aerobic respiration rate       [d-1]
+  b_H_NO:   0.1,  // anoxic respiration rate        [d-1]
+  b_STO_O:  0.2,  // aerobic respitation rate X_STO [d-1]
+  b_STO_NO: 0.1,  // anoxic respitation rate X_STO  [d-1]
+  // Autotrophs
+  mu_A_max: 1.0,  // maximum specific growth rate   [d-1]
+  K_A_NH:   1.,   // saturation constant S_NH3      [g NH3-N m-3] 
+  K_A_O:    0.5,  // saturation constant S_0        [g O2 m-3]
+  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 = {
+  // 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]
+  // Yields
+  Y_STO_O:  0.85, // aerobic yield X_STO per S_S    [g X_STO g-1 S_S]
+  Y_STO_NO: 0.80, // anoxic yield X_STO per S_S     [g X_STO g-1 S_S]
+  Y_H_O:    0.63, // aerobic yield X_H per X_STO    [g X_H g-1 X_STO]
+  Y_H_NO:   0.54, // anoxic yield X_H per X_STO     [g X_H g-1 X_STO]
+  Y_A:      0.24, // anoxic yield X_A per S_NO      [g X_A g-1 NO3-N]
+  // Composition (COD via DoR)
+  i_CODN:  -1.71, // COD content (DoR)              [g COD g-1 N2-N]
+  i_CODNO: -4.57, // COD content (DoR)              [g COD g-1 NO3-N]
+  // Composition (nitrogen)
+  i_NSI:    0.01, // nitrogen content S_I           [g N g-1 S_I]
+  i_NSS:    0.03, // nitrogen content S_S           [g N g-1 S_S]
+  i_NXI:    0.02, // nitrogen content X_I           [g N g-1 X_I]
+  i_NXS:    0.04, // nitrogen content X_S           [g N g-1 X_S]
+  i_NBM:    0.07, // nitrogen content X_H / X_A     [g N g-1 X_H / X_A]
+  // Composition (TSS)
+  i_TSXI:   0.75, // TSS content X_I                [g TS g-1 X_I]
+  i_TSXS:   0.75, // TSS content X_S                [g TS g-1 X_S]
+  i_TSBM:   0.90, // TSS content X_H / X_A          [g TS g-1 X_H / X_A]
+  i_TSSTO:  0.60, // TSS content X_STO (PHB based)  [g TS g-1 X_STO]
+  // 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).
@@ -10,65 +73,13 @@ class ASM3 {
      * 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] 
-      // Heterotrophs
-      k_STO:    5.,   // storage rate constant          [g S_S g-1 X_H d-1]
-      nu_NO:    0.6,  // 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:      2.,   // saturation constant S_s        [g COD m-3]
-      K_STO:    1.,   // saturation constant X_STO      [g X_STO g-1 X_H]
-      mu_H_max: 2.,   // maximum specific growth rate   [d-1]
-      K_NH:     0.01, // saturation constant S_NH3      [g NH3-N m-3]
-      K_HCO:    0.1,  // saturation constant S_HCO      [mole HCO3 m-3]
-      b_H_O:    0.2,  // aerobic respiration rate       [d-1]
-      b_H_NO:   0.1,  // anoxic respiration rate        [d-1]
-      b_STO_O:  0.2,  // aerobic respitation rate X_STO [d-1]
-      b_STO_NO: 0.1,  // anoxic respitation rate X_STO  [d-1]
-      // Autotrophs
-      mu_A_max: 1.0,  // maximum specific growth rate   [d-1]
-      K_A_NH:   1.,   // saturation constant S_NH3      [g NH3-N m-3] 
-      K_A_O:    0.5,  // saturation constant S_0        [g O2 m-3]
-      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]
-    };
+    this.kin_params = KINETIC_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]
-      // Yields
-      Y_STO_O:  0.85, // aerobic yield X_STO per S_S    [g X_STO g-1 S_S]
-      Y_STO_NO: 0.80, // anoxic yield X_STO per S_S     [g X_STO g-1 S_S]
-      Y_H_O:    0.63, // aerobic yield X_H per X_STO    [g X_H g-1 X_STO]
-      Y_H_NO:   0.54, // anoxic yield X_H per X_STO     [g X_H g-1 X_STO]
-      Y_A:      0.24, // anoxic yield X_A per S_NO      [g X_A g-1 NO3-N]
-      // Composition (COD via DoR)
-      i_CODN:  -1.71, // COD content (DoR)              [g COD g-1 N2-N]
-      i_CODNO: -4.57, // COD content (DoR)              [g COD g-1 NO3-N]
-      // Composition (nitrogen)
-      i_NSI:    0.01, // nitrogen content S_I           [g N g-1 S_I]
-      i_NSS:    0.03, // nitrogen content S_S           [g N g-1 S_S]
-      i_NXI:    0.02, // nitrogen content X_I           [g N g-1 X_I]
-      i_NXS:    0.04, // nitrogen content X_S           [g N g-1 X_S]
-      i_NBM:    0.07, // nitrogen content X_H / X_A     [g N g-1 X_H / X_A]
-      // Composition (TSS)
-      i_TSXI:   0.75, // TSS content X_I                [g TS g-1 X_I]
-      i_TSXS:   0.75, // TSS content X_S                [g TS g-1 X_S]
-      i_TSBM:   0.90, // TSS content X_H / X_A          [g TS g-1 X_H / X_A]
-      i_TSSTO:  0.60, // TSS content X_STO (PHB based)  [g TS g-1 X_STO]
-      // 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]
-    };
+    this.stoi_params = STOICHIOMETRIC_CONSTANTS;
 
     /**
      * Temperature theta parameters for ASM3.
@@ -208,4 +219,4 @@ class ASM3 {
   }
 }
 
-module.exports = ASM3;
+module.exports = { ASM3, ASM_CONSTANTS, KINETIC_CONSTANTS, STOICHIOMETRIC_CONSTANTS };

src/specificClass.js

@@ -1,4 +1,4 @@
-const ASM3 = require('./reaction_modules/asm3_class.js');
+const { ASM3, ASM_CONSTANTS } = require('./reaction_modules/asm3_class.js');
 const { create, all, isArray } = require('mathjs');
 const { assertNoNaN } = require('./utils.js');
 const { childRegistrationUtils, logger, MeasurementContainer } = require('generalFunctions');
@@ -10,9 +10,9 @@ const mathConfig = {
 
 const math = create(all, mathConfig);
 
-const S_O_INDEX = 0;
-const NUM_SPECIES = 13;
+const BC_PADDING = 2; // Boundary Condition padding for open boundaries in extendedState variable
 const DEBUG = false;
+const DAY2MS = 1000 * 60 * 60 * 24; // convert between days and milliseconds
 
 class Reactor {
   /**
@@ -25,23 +25,27 @@ 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.upstreamReactor = null;
-    this.childRegistrationUtils = new childRegistrationUtils(this); // Child registration utility
+    this.childRegistrationUtils = new childRegistrationUtils(this); // child registration utility
 
-    this.asm = new ASM3();
+    // placeholder variables for children and parents
+    this.upstreamReactor = null;
+    this.downstreamReactor = null;
+    this.returnPump = null;
+
+    this.asm = new ASM3(); // Reaction model
 
     this.volume = config.volume; // fluid volume reactor [m3]
 
-    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.Fs = [0]; // fluid debits per inlet [m3 d-1]
+    this.Cs_in = [Array(ASM_CONSTANTS.NUM_SPECIES).fill(0)]; // composition influents
     this.OTR = 0.0; // oxygen transfer rate [g O2 d-1 m-3]
     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.currentTime = null; // 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 = 1; // speed up factor for simulation, 60 means 1 minute per simulated second
   }
 
   /**
@@ -49,9 +53,15 @@ class Reactor {
    * @param {object} input - Input object (msg) containing payload with inlet index, flow rate, and concentrations.
   */
   set setInfluent(input) {
-    let index_in = input.payload.inlet;
-    this.Fs[index_in] = input.payload.F;
-    this.Cs_in[index_in] = input.payload.C;
+    const i_in = input.payload.inlet;
+    if (this.Fs.length <= i_in) {
+      this.logger.debug(`Adding new inlet index ${i_in}.`);
+      this.Fs.push(0);
+      this.Cs_in.push(Array(ASM_CONSTANTS.NUM_SPECIES).fill(0));
+      this.setInfluent = input;
+    }
+    this.Fs[i_in] = input.payload.F;
+    this.Cs_in[i_in] = input.payload.C;
   }
 
   /**
@@ -64,13 +74,20 @@ class Reactor {
 
   /**
    * Getter for effluent data.
-   * @returns {object} Effluent data object (msg), defaults to inlet 0.
+   * @returns {object} Effluent data object (msg).
    */
-  get getEffluent() { // getter for Effluent, defaults to inlet 0
-    if (isArray(this.state.at(-1))) {
-      return { topic: "Fluent", payload: { inlet: 0, F: math.sum(this.Fs), C: this.state.at(-1) }, timestamp: this.currentTime };
+  get getEffluent() {
+    const Cs = isArray(this.state.at(-1)) ? this.state.at(-1) : this.state;
+    const effluent = [{ topic: "Fluent", payload: { inlet: 0, F: math.sum(this.Fs), C: Cs }, timestamp: this.currentTime }];
+    if (this.returnPump) {
+      const recirculationFlow = this.returnPump.measurements.type("flow").variant("measured").position("atEquipment").getCurrentValue();
+      // constrain flow to prevent negatives
+      const F_main = Math.max(effluent[0].payload.F - recirculationFlow, 0);
+      const F_sidestream = effluent[0].payload.F < recirculationFlow ? effluent[0].payload.F : recirculationFlow;
+      effluent[0].payload.F = F_main;
+      effluent.push({ topic: "Fluent", payload: { inlet: 1, F: F_sidestream, C: Cs }, timestamp: this.currentTime });
     }
-    return { topic: "Fluent", payload: { inlet: 0, F: math.sum(this.Fs), C: this.state }, timestamp: this.currentTime };
+    return effluent;
   }
 
   /**
@@ -80,7 +97,7 @@ class Reactor {
    * @returns {number} - Calculated OTR [g O2 d-1 m-3].
    */
   _calcOTR(S_O, T = 20.0) { // caculate the OTR using basic correlation, default to temperature: 20 C
-    let S_O_sat = 14.652 - 4.1022e-1 * T + 7.9910e-3 * T*T + 7.7774e-5 * T*T*T;
+    const S_O_sat = 14.652 - 4.1022e-1 * T + 7.9910e-3 * T*T + 7.7774e-5 * T*T*T;
     return this.kla * (S_O_sat - S_O);
   }
 
@@ -97,40 +114,43 @@ 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.`);
+      return;
+    }
+
     switch (softwareType) {
       case "measurement":
-        this.logger.debug(`Registering measurement child.`);
+        this.logger.debug(`Registering measurement child...`);
         this._connectMeasurement(child);
         break;
       case "reactor":
-        this.logger.debug(`Registering reactor child.`);
+        this.logger.debug(`Registering reactor child...`);
         this._connectReactor(child);
         break;
+      case "machine":
+        this.logger.debug(`Registering machine child...`);
+        this._connectMachine(child);
+        break;
 
       default:
         this.logger.error(`Unrecognized softwareType: ${softwareType}`);
     }
   }
 
-  _connectMeasurement(measurement) {
-    if (!measurement) {
-      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 key = `${measurementType}_${position}`;
+  _connectMeasurement(measurementChild) {
+    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
@@ -145,20 +165,27 @@ class Reactor {
   }
 
 
-  _connectReactor(reactor) {
-    if (!reactor) {
-      this.logger.warn("Invalid reactor provided.");
-      return;
+  _connectReactor(reactorChild) {
+    if (reactorChild.config.functionality.positionVsParent != "upstream") {
+      this.logger.warn("Reactor children of other reactors should always be upstream!");
     }
 
-    this.upstreamReactor = reactor;
+    // set upstream and downstream reactor variable in current and child nodes respectively for easy access
+    this.upstreamReactor = reactorChild;
+    reactorChild.downstreamReactor = this;
 
-    reactor.emitter.on("stateChange", (data) => {
+    reactorChild.emitter.on("stateChange", (eventData) => { // Triggers state update in downstream reactor.
       this.logger.debug(`State change of upstream reactor detected.`);
-      this.updateState(data);
+      this.updateState(eventData);
     });
   }
 
+  _connectMachine(machineChild) {
+    if (machineChild.config.functionality.positionVsParent == "downstream") {
+      machineChild.upstreamSource = this;
+      this.returnPump = machineChild;
+    }
+  }
   
   _updateMeasurement(measurementType, value, position, context) {
     this.logger.debug(`---------------------- updating ${measurementType} ------------------ `);
@@ -178,22 +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 = Date.now()) { // expect update with timestamp
-    const day2ms = 1000 * 60 * 60 * 24;
-
-    if (this.upstreamReactor) {
-      this.setInfluent = this.upstreamReactor.getEffluent;
+  updateState(newTime) {
+    if (!this.currentTime) { // initialise currentTime variable
+      this.currentTime = newTime;
+      return;
     }
 
-    let 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 (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 == 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();
     }
   }
 }
@@ -208,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].
@@ -217,8 +271,8 @@ class Reactor_CSTR extends Reactor {
     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, 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], this.temperature); // calculate OTR if kla is not NaN, otherwise use externaly calculated OTR
+    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[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
@@ -244,13 +298,15 @@ class Reactor_PFR extends Reactor {
     this.d_x = this.length / this.n_x;
     this.A = this.volume / this.length; // crosssectional area [m2]
 
-    this.alpha = config.alpha;
+    this.state = Array.from(Array(this.n_x), () => config.initialState.slice());
+    this.extendedState = Array.from(Array(this.n_x + 2*BC_PADDING), () => new Array(ASM_CONSTANTS.NUM_SPECIES).fill(0));
 
-    this.state = Array.from(Array(this.n_x), () => config.initialState.slice())
+    // 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();
@@ -262,15 +318,26 @@ 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);
   }
 
+  _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)
-    let Co_D = this.D*this.timeStep/(this.d_x*this.d_x);
+    
+    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) {
@@ -288,25 +355,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(ASM_CONSTANTS.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++) {
-        transfer[i][S_O_INDEX] = this.OTR * this.n_x/(this.n_x-2);
+      for (let i = BC_PADDING+1; i < BC_PADDING+this.n_x - 1; i++) {
+        transfer[i][ASM_CONSTANTS.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][ASM_CONSTANTS.S_O_INDEX] = this._calcOTR(this.extendedState[i][ASM_CONSTANTS.S_O_INDEX], this.temperature);
       }
     }
 
-    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");
@@ -317,14 +385,24 @@ 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) {
+    const grid_pos = Math.round(context.distance / this.config.length * this.n_x);
+
+    // naive approach for reconciling measurements and simulation
+    // could benefit from Kalman filter?
     switch(measurementType) {
       case "quantity (oxygen)":
-        let grid_pos = Math.round(position / this.config.length * this.n_x);
-        this.state[grid_pos][S_O_INDEX] = value; // naive approach for reconciling measurements and simulation
+        this.state[grid_pos][ASM_CONSTANTS.S_O_INDEX] = value;
+        break;
+      case "quantity (ammonium)":
+        this.state[grid_pos][ASM_CONSTANTS.S_NH_INDEX] = value;
+        break;
+      case "quantity (nox)":
+        this.state[grid_pos][ASM_CONSTANTS.S_NO_INDEX] = value;
         break;
       default:
         super._updateMeasurement(measurementType, value, position, context);
@@ -335,85 +413,76 @@ 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() {
+    // Upstream BC
+    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 {
-      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 = 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 (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)
+        this.extendedState.fill(this.extendedState[BC_PADDING], 0, BC_PADDING);
+      }
+    }
+
+    // Downstream BC
+    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)
+      this.extendedState.fill(this.extendedState.at(-1-BC_PADDING), BC_PADDING+this.n_x);
     }
-    // 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;
   }
 
   /**
    * Create central finite difference second derivative operator.
-   * @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;
   }
+
+  /**
+   * 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) {
+      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 };
-
-// DEBUG
-// 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
-// let initial_state = [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1];
-// const Reactor = new Reactor_PFR(200, 10, 10, 1, 100, initial_state);
-// Reactor.Cs_in[0] = [0.0, 30., 100., 16., 0., 0., 5., 25., 75., 30., 0., 0., 125.];
-// Reactor.Fs[0] = 10;
-// Reactor.D = 0.01;
-// let N = 0;
-// while (N < 5000) {
-//     console.log(Reactor.tick(0.001));
-//     N += 1;
-// }
+module.exports = { Reactor_CSTR, Reactor_PFR };