Add example flow

Reviewed-on: #6

README.md

@@ -1,17 +1,20 @@
-# reactor
+# Reactor: Advanced Hydraulic Tank & Biological Process Simulator
 
-Reactor: Advanced Hydraulic Tank & Biological Process Simulator
+A comprehensive reactor class for wastewater treatment simulation featuring non-ideal plug flow hydraulics and ASM3 biological modeling.
 
-A comprehensive reactor class for wastewater treatment simulation featuring plug flow hydraulics, ASM1-ASM3 biological modeling, and multi-sectional concentration tracking. Implements hydraulic retention time calculations, dispersion modeling, and real-time biological reaction kinetics for accurate process simulation.
+## How to use this Node
+### Set Node Properties
+- Set reactor type: Continuously Stirred Tank Reactor (CSTR) or Plug Flow Reactor (PFR)
+- Configure reactor sizing: Set reactor volume [ $m^3$ ] and (for PFRs) set the reactor length [ $m$ ]
+- (For PFRs) set reactor spatial resolution: A value of 10 or 20 is good. A higher resolution means more accurate simulation, at higher computational expense. Note that connected reactors must have similar Δx values.
+- Set initial state of reactor: set the intial concentrations of all relevant reaction species.
+- (Optional) set $k_L a$ to calculate OTR internally, rather than providing it explicitly, using simple mass transfer model.
 
-Key Features:
-
-Plug Flow Hydraulics: Multi-section reactor with configurable sectioning factor and dispersion modeling
-ASM1 Integration: Complete biological nutrient removal modeling with 13 state variables (COD, nitrogen, phosphorus)
-Dynamic Volume Control: Automatic section management with overflow handling and retention time calculations
-Oxygen Transfer: Saturation-limited O2 transfer with Fick's law slowdown effects and solubility curves
-Real-time Kinetics: Continuous biological reaction rate calculations with configurable time acceleration
-Weighted Averaging: Volume-based concentration mixing for accurate mass balance calculations
-Child Registration: Integration with diffuser systems and upstream/downstream reactor networks
-Supports complex biological treatment train modeling with temperature compensation, sludge calculations, and comprehensive process monitoring for wastewater treatment plant optimization and regulatory compliance.
+### Accepted Node inputs
+- \{ topic: clock, payload: \<timestamp [ $ms$ ]\> \} - **required** clock signal to make reactor update state.
+- \{ topic: Fluent, payload: \{ F: \<flow rate [ $m^3 d^{-1}$ ]\>, C: \<array with concentrations\> \} \} - sets inflow composition and flow rate.
+- \{ topic: Dispersion, payload: \<dispersion coefficient in [ $m d^{-2}$ ]\> \} - sets PFR dispersion coefficient.
+- \{ topic: OTR, payload: \<oxygen transfer rate [ $ g d^{-1} m^{-3}$ ]\> \} - sets current oxygen transfer rate.
 
+## Troubleshooting
+Check for possible numerical warnings. These tell you which simulation parameters to change. If solutions appear to be oscillate, try reducing the time step. If solutions appear to be too dispersive, try increasing the reactor resolution.

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);
-};

example_flow/reactor_flows.json

@@ -0,0 +1,838 @@
+[
+    {
+        "id": "a6b85e226d144df1",
+        "type": "tab",
+        "label": "Flow 3",
+        "disabled": false,
+        "info": "",
+        "env": []
+    },
+    {
+        "id": "a94c85d65c71b66a",
+        "type": "inject",
+        "z": "a6b85e226d144df1",
+        "name": "Influx composition 1",
+        "props": [
+            {
+                "p": "payload"
+            },
+            {
+                "p": "topic",
+                "vt": "str"
+            },
+            {
+                "p": "timestamp",
+                "v": "",
+                "vt": "date"
+            }
+        ],
+        "repeat": "864",
+        "crontab": "",
+        "once": true,
+        "onceDelay": "5",
+        "topic": "Fluent",
+        "payload": "{\"inlet\":0,\"F\":6600,\"C\":[0,30,100,16,0,0,5,25,75,30,0,0,125]}",
+        "payloadType": "json",
+        "x": 260,
+        "y": 240,
+        "wires": [
+            [
+                "b726c5be41c24dcb"
+            ]
+        ]
+    },
+    {
+        "id": "0ebfbbf57bba79f1",
+        "type": "inject",
+        "z": "a6b85e226d144df1",
+        "name": "",
+        "props": [
+            {
+                "p": "timestamp",
+                "v": "",
+                "vt": "date"
+            },
+            {
+                "p": "topic",
+                "vt": "str"
+            }
+        ],
+        "repeat": "1",
+        "crontab": "",
+        "once": true,
+        "onceDelay": 0.1,
+        "topic": "clock",
+        "x": 300,
+        "y": 280,
+        "wires": [
+            [
+                "b726c5be41c24dcb"
+            ]
+        ]
+    },
+    {
+        "id": "cd0e4a78d1a59a6e",
+        "type": "inject",
+        "z": "a6b85e226d144df1",
+        "name": "",
+        "props": [
+            {
+                "p": "payload"
+            },
+            {
+                "p": "topic",
+                "vt": "str"
+            }
+        ],
+        "repeat": "",
+        "crontab": "",
+        "once": true,
+        "onceDelay": 0.1,
+        "topic": "Dispersion",
+        "payload": "100",
+        "payloadType": "num",
+        "x": 260,
+        "y": 320,
+        "wires": [
+            [
+                "b726c5be41c24dcb",
+                "b0819ff9a4010227"
+            ]
+        ]
+    },
+    {
+        "id": "b726c5be41c24dcb",
+        "type": "reactor",
+        "z": "a6b85e226d144df1",
+        "name": "Anoxic 1",
+        "reactor_type": "PFR",
+        "volume": "730",
+        "length": "10",
+        "resolution_L": "20",
+        "kla": "",
+        "S_O_init": 0,
+        "S_I_init": 30,
+        "S_S_init": 100,
+        "S_NH_init": 16,
+        "S_N2_init": 0,
+        "S_NO_init": 0,
+        "S_HCO_init": 5,
+        "X_I_init": 25,
+        "X_S_init": 75,
+        "X_H_init": 30,
+        "X_STO_init": 0,
+        "X_A_init": "30",
+        "X_TS_init": "132",
+        "timeStep": "2",
+        "enableLog": true,
+        "logLevel": "info",
+        "positionVsParent": "upstream",
+        "x": 540,
+        "y": 240,
+        "wires": [
+            [
+                "057ab2dcd4739aef"
+            ],
+            [],
+            [
+                "b0819ff9a4010227",
+                "214aff7330e81d3f",
+                "2ec38c4ae9aa6e7e"
+            ]
+        ]
+    },
+    {
+        "id": "b0819ff9a4010227",
+        "type": "reactor",
+        "z": "a6b85e226d144df1",
+        "name": "Aerobic 1",
+        "reactor_type": "PFR",
+        "volume": "1460",
+        "length": "20",
+        "resolution_L": "40",
+        "kla": "2400",
+        "S_O_init": 0,
+        "S_I_init": 30,
+        "S_S_init": 100,
+        "S_NH_init": 16,
+        "S_N2_init": 0,
+        "S_NO_init": 0,
+        "S_HCO_init": 5,
+        "X_I_init": 25,
+        "X_S_init": 75,
+        "X_H_init": "500",
+        "X_STO_init": 0,
+        "X_A_init": "30",
+        "X_TS_init": "132",
+        "timeStep": "2",
+        "enableLog": true,
+        "logLevel": "info",
+        "positionVsParent": "upstream",
+        "x": 940,
+        "y": 240,
+        "wires": [
+            [
+                "e0049d66fefeb3e6"
+            ],
+            [],
+            [
+                "47120bc82aa7bf49"
+            ]
+        ]
+    },
+    {
+        "id": "214aff7330e81d3f",
+        "type": "rotatingMachine",
+        "z": "a6b85e226d144df1",
+        "speed": 1,
+        "startup": 0,
+        "warmup": 0,
+        "shutdown": 0,
+        "cooldown": 0,
+        "machineCurve": {},
+        "flowNumber": "1",
+        "uuid": "",
+        "supplier": "Hidrostal",
+        "category": "Pumps",
+        "assetType": "Centrifugal",
+        "model": "hidrostal-H05K-S03R",
+        "unit": "l/s",
+        "enableLog": true,
+        "logLevel": "info",
+        "positionVsParent": "downstream",
+        "positionIcon": "→",
+        "hasDistance": false,
+        "distance": "",
+        "x": 740,
+        "y": 340,
+        "wires": [
+            [],
+            [],
+            [
+                "b0819ff9a4010227"
+            ]
+        ]
+    },
+    {
+        "id": "0be74d5c77febec1",
+        "type": "measurement",
+        "z": "a6b85e226d144df1",
+        "name": "sensor",
+        "scaling": false,
+        "i_min": 0,
+        "i_max": 0,
+        "i_offset": 0,
+        "o_min": 900,
+        "o_max": 1200,
+        "simulator": true,
+        "smooth_method": "",
+        "count": 10,
+        "uuid": "",
+        "supplier": "Vega",
+        "category": "Sensor",
+        "assetType": "Flow",
+        "model": "VegaFlow 10",
+        "unit": "m³/h",
+        "enableLog": true,
+        "logLevel": "error",
+        "positionVsParent": "atEquipment",
+        "positionIcon": "⊥",
+        "hasDistance": false,
+        "distance": "",
+        "x": 600,
+        "y": 460,
+        "wires": [
+            [],
+            [],
+            [
+                "214aff7330e81d3f"
+            ]
+        ]
+    },
+    {
+        "id": "668f41a05698f21b",
+        "type": "inject",
+        "z": "a6b85e226d144df1",
+        "name": "",
+        "props": [
+            {
+                "p": "payload"
+            },
+            {
+                "p": "topic",
+                "vt": "str"
+            }
+        ],
+        "repeat": "",
+        "crontab": "",
+        "once": true,
+        "onceDelay": 0.1,
+        "topic": "execSequence",
+        "payload": "{\"source\":\"parent\",\"action\":\"execSequence\",\"parameter\":\"startup\"}",
+        "payloadType": "json",
+        "x": 330,
+        "y": 400,
+        "wires": [
+            [
+                "214aff7330e81d3f",
+                "2ec38c4ae9aa6e7e"
+            ]
+        ]
+    },
+    {
+        "id": "47120bc82aa7bf49",
+        "type": "settler",
+        "z": "a6b85e226d144df1",
+        "name": "",
+        "model": "mb-model",
+        "enableLog": true,
+        "logLevel": "info",
+        "positionVsParent": "atEquipment",
+        "x": 1100,
+        "y": 400,
+        "wires": [
+            [
+                "b1af49d1d0eb6783",
+                "6baffb7954cf0cce"
+            ],
+            [],
+            []
+        ]
+    },
+    {
+        "id": "2ec38c4ae9aa6e7e",
+        "type": "rotatingMachine",
+        "z": "a6b85e226d144df1",
+        "speed": 1,
+        "startup": 0,
+        "warmup": 0,
+        "shutdown": 0,
+        "cooldown": 0,
+        "machineCurve": {},
+        "flowNumber": "2",
+        "uuid": "",
+        "supplier": "Hidrostal",
+        "category": "Pumps",
+        "assetType": "Centrifugal",
+        "model": "hidrostal-H05K-S03R",
+        "unit": "l/s",
+        "enableLog": true,
+        "logLevel": "info",
+        "positionVsParent": "downstream",
+        "positionIcon": "→",
+        "hasDistance": false,
+        "distance": "",
+        "x": 860,
+        "y": 460,
+        "wires": [
+            [],
+            [],
+            [
+                "47120bc82aa7bf49"
+            ]
+        ]
+    },
+    {
+        "id": "21eff612371519c5",
+        "type": "measurement",
+        "z": "a6b85e226d144df1",
+        "name": "sensor",
+        "scaling": false,
+        "i_min": 0,
+        "i_max": 0,
+        "i_offset": 0,
+        "o_min": 500,
+        "o_max": 660,
+        "simulator": true,
+        "smooth_method": "",
+        "count": 10,
+        "uuid": "",
+        "supplier": "Vega",
+        "category": "Sensor",
+        "assetType": "Flow",
+        "model": "VegaFlow 10",
+        "unit": "m³/h",
+        "enableLog": true,
+        "logLevel": "error",
+        "positionVsParent": "atEquipment",
+        "positionIcon": "⊥",
+        "hasDistance": false,
+        "distance": "",
+        "x": 720,
+        "y": 580,
+        "wires": [
+            [],
+            [],
+            [
+                "2ec38c4ae9aa6e7e"
+            ]
+        ]
+    },
+    {
+        "id": "057ab2dcd4739aef",
+        "type": "function",
+        "z": "a6b85e226d144df1",
+        "name": "Data_converter",
+        "func": "if (msg.topic != \"Fluent\") {\n    return;\n}\n\nconst [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] = msg.payload.C;\n\nmsg = {payload: [\n    { \"Series\": \"S_O\",      \"Y\": S_O},\n    { \"Series\": \"S_I\",      \"Y\": S_I},\n    { \"Series\": \"S_S\",      \"Y\": S_S},\n    { \"Series\": \"S_NH\",     \"Y\": S_NH},\n    { \"Series\": \"S_N2\",     \"Y\": S_N2},\n    { \"Series\": \"S_NO\",     \"Y\": S_NO},\n    { \"Series\": \"S_HCO\",    \"Y\": S_HCO},\n    { \"Series\": \"X_I\",      \"Y\": X_I},\n    { \"Series\": \"X_S\",      \"Y\": X_S},\n    { \"Series\": \"X_H\",      \"Y\": X_H},\n    { \"Series\": \"X_STO\",    \"Y\": X_STO},\n    { \"Series\": \"X_A\",      \"Y\": X_A},\n    { \"Series\": \"X_TS\",     \"Y\": X_TS}\n    ]};\n\nreturn msg;",
+        "outputs": 1,
+        "timeout": 0,
+        "noerr": 0,
+        "initialize": "",
+        "finalize": "",
+        "libs": [],
+        "x": 720,
+        "y": 180,
+        "wires": [
+            [
+                "02a932317c76721e"
+            ]
+        ]
+    },
+    {
+        "id": "02a932317c76721e",
+        "type": "ui-chart",
+        "z": "a6b85e226d144df1",
+        "group": "ae38454098a37db0",
+        "name": "Anoxic reactor",
+        "label": "Anoxic reactor",
+        "order": 9007199254740991,
+        "chartType": "line",
+        "category": "Series",
+        "categoryType": "property",
+        "xAxisLabel": "",
+        "xAxisProperty": "",
+        "xAxisPropertyType": "timestamp",
+        "xAxisType": "time",
+        "xAxisFormat": "",
+        "xAxisFormatType": "auto",
+        "xmin": "",
+        "xmax": "",
+        "yAxisLabel": "",
+        "yAxisProperty": "Y",
+        "yAxisPropertyType": "property",
+        "ymin": "",
+        "ymax": "",
+        "bins": 10,
+        "action": "append",
+        "stackSeries": false,
+        "pointShape": "circle",
+        "pointRadius": 4,
+        "showLegend": true,
+        "removeOlder": "8",
+        "removeOlderUnit": "3600",
+        "removeOlderPoints": "2000",
+        "colors": [
+            "#0095ff",
+            "#ff0000",
+            "#ff7f0e",
+            "#2ca02c",
+            "#a347e1",
+            "#d62728",
+            "#ff9896",
+            "#9467bd",
+            "#c5b0d5"
+        ],
+        "textColor": [
+            "#666666"
+        ],
+        "textColorDefault": true,
+        "gridColor": [
+            "#e5e5e5"
+        ],
+        "gridColorDefault": true,
+        "width": 6,
+        "height": 8,
+        "className": "",
+        "interpolation": "linear",
+        "x": 920,
+        "y": 180,
+        "wires": [
+            []
+        ]
+    },
+    {
+        "id": "e0049d66fefeb3e6",
+        "type": "function",
+        "z": "a6b85e226d144df1",
+        "name": "Data_converter",
+        "func": "if (msg.topic != \"Fluent\") {\n    return;\n}\n\nconst [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] = msg.payload.C;\n\nmsg = {payload: [\n    { \"Series\": \"S_O\",      \"Y\": S_O},\n    { \"Series\": \"S_I\",      \"Y\": S_I},\n    { \"Series\": \"S_S\",      \"Y\": S_S},\n    { \"Series\": \"S_NH\",     \"Y\": S_NH},\n    { \"Series\": \"S_N2\",     \"Y\": S_N2},\n    { \"Series\": \"S_NO\",     \"Y\": S_NO},\n    { \"Series\": \"S_HCO\",    \"Y\": S_HCO},\n    { \"Series\": \"X_I\",      \"Y\": X_I},\n    { \"Series\": \"X_S\",      \"Y\": X_S},\n    { \"Series\": \"X_H\",      \"Y\": X_H},\n    { \"Series\": \"X_STO\",    \"Y\": X_STO},\n    { \"Series\": \"X_A\",      \"Y\": X_A},\n    { \"Series\": \"X_TS\",     \"Y\": X_TS}\n    ]};\n\nreturn msg;",
+        "outputs": 1,
+        "timeout": 0,
+        "noerr": 0,
+        "initialize": "",
+        "finalize": "",
+        "libs": [],
+        "x": 1140,
+        "y": 220,
+        "wires": [
+            [
+                "70eeb4c8caa2db77"
+            ]
+        ]
+    },
+    {
+        "id": "70eeb4c8caa2db77",
+        "type": "ui-chart",
+        "z": "a6b85e226d144df1",
+        "group": "ae38454098a37db0",
+        "name": "Aerobic reactor",
+        "label": "Aerobic reactor / recirculation",
+        "order": 9007199254740991,
+        "chartType": "line",
+        "category": "Series",
+        "categoryType": "property",
+        "xAxisLabel": "",
+        "xAxisProperty": "",
+        "xAxisPropertyType": "timestamp",
+        "xAxisType": "time",
+        "xAxisFormat": "",
+        "xAxisFormatType": "auto",
+        "xmin": "",
+        "xmax": "",
+        "yAxisLabel": "",
+        "yAxisProperty": "Y",
+        "yAxisPropertyType": "property",
+        "ymin": "",
+        "ymax": "",
+        "bins": 10,
+        "action": "append",
+        "stackSeries": false,
+        "pointShape": "circle",
+        "pointRadius": 4,
+        "showLegend": true,
+        "removeOlder": "8",
+        "removeOlderUnit": "3600",
+        "removeOlderPoints": "2000",
+        "colors": [
+            "#0095ff",
+            "#ff0000",
+            "#ff7f0e",
+            "#2ca02c",
+            "#a347e1",
+            "#d62728",
+            "#ff9896",
+            "#9467bd",
+            "#c5b0d5"
+        ],
+        "textColor": [
+            "#666666"
+        ],
+        "textColorDefault": true,
+        "gridColor": [
+            "#e5e5e5"
+        ],
+        "gridColorDefault": true,
+        "width": 6,
+        "height": 8,
+        "className": "",
+        "interpolation": "linear",
+        "x": 1340,
+        "y": 220,
+        "wires": [
+            []
+        ]
+    },
+    {
+        "id": "b1af49d1d0eb6783",
+        "type": "function",
+        "z": "a6b85e226d144df1",
+        "name": "Data_converter",
+        "func": "if (msg.topic != \"Fluent\" || msg.payload.inlet == 0) {\n    return;\n}\n\nconst [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] = msg.payload.C;\n\nmsg = {payload: [\n    { \"Series\": \"X_I\",      \"Y\": X_I},\n    { \"Series\": \"X_S\",      \"Y\": X_S},\n    { \"Series\": \"X_H\",      \"Y\": X_H},\n    { \"Series\": \"X_STO\",    \"Y\": X_STO},\n    { \"Series\": \"X_A\",      \"Y\": X_A},\n    { \"Series\": \"X_TS\",     \"Y\": X_TS}\n    ]};\n\nreturn msg;",
+        "outputs": 1,
+        "timeout": 0,
+        "noerr": 0,
+        "initialize": "",
+        "finalize": "",
+        "libs": [],
+        "x": 1320,
+        "y": 380,
+        "wires": [
+            [
+                "0a3ffa69046a4844"
+            ]
+        ]
+    },
+    {
+        "id": "0a3ffa69046a4844",
+        "type": "ui-chart",
+        "z": "a6b85e226d144df1",
+        "group": "de8b029d69f26c0e",
+        "name": "Sludge composition",
+        "label": "Sludge composition",
+        "order": 9007199254740991,
+        "chartType": "line",
+        "category": "Series",
+        "categoryType": "property",
+        "xAxisLabel": "",
+        "xAxisProperty": "",
+        "xAxisPropertyType": "timestamp",
+        "xAxisType": "time",
+        "xAxisFormat": "",
+        "xAxisFormatType": "auto",
+        "xmin": "",
+        "xmax": "",
+        "yAxisLabel": "",
+        "yAxisProperty": "Y",
+        "yAxisPropertyType": "property",
+        "ymin": "",
+        "ymax": "",
+        "bins": 10,
+        "action": "append",
+        "stackSeries": false,
+        "pointShape": "circle",
+        "pointRadius": 4,
+        "showLegend": true,
+        "removeOlder": "8",
+        "removeOlderUnit": "3600",
+        "removeOlderPoints": "2000",
+        "colors": [
+            "#0095ff",
+            "#ff0000",
+            "#ff7f0e",
+            "#2ca02c",
+            "#a347e1",
+            "#d62728",
+            "#ff9896",
+            "#9467bd",
+            "#c5b0d5"
+        ],
+        "textColor": [
+            "#666666"
+        ],
+        "textColorDefault": true,
+        "gridColor": [
+            "#e5e5e5"
+        ],
+        "gridColorDefault": true,
+        "width": 6,
+        "height": 8,
+        "className": "",
+        "interpolation": "linear",
+        "x": 1530,
+        "y": 380,
+        "wires": [
+            []
+        ]
+    },
+    {
+        "id": "272b4b0050479d13",
+        "type": "measurement",
+        "z": "a6b85e226d144df1",
+        "name": "sensor",
+        "scaling": false,
+        "i_min": 0,
+        "i_max": 0,
+        "i_offset": 0,
+        "o_min": 3200,
+        "o_max": 4000,
+        "simulator": true,
+        "smooth_method": "",
+        "count": 10,
+        "uuid": "",
+        "supplier": "Vega",
+        "category": "Sensor",
+        "assetType": "Quantity (TSS)",
+        "model": "VegaSolidsProbe",
+        "unit": "g/m³",
+        "enableLog": true,
+        "logLevel": "error",
+        "positionVsParent": "atEquipment",
+        "positionIcon": "⊥",
+        "hasDistance": false,
+        "distance": "",
+        "x": 930,
+        "y": 580,
+        "wires": [
+            [],
+            [],
+            [
+                "47120bc82aa7bf49"
+            ]
+        ]
+    },
+    {
+        "id": "6baffb7954cf0cce",
+        "type": "function",
+        "z": "a6b85e226d144df1",
+        "name": "Data_converter",
+        "func": "if (msg.topic != \"Fluent\" || msg.payload.inlet == 1 || msg.payload.inlet == 2) {\n    return;\n}\n\nconst [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] = msg.payload.C;\n\nmsg = {payload: [\n    { \"Series\": \"S_O\",      \"Y\": S_O},\n    { \"Series\": \"S_I\",      \"Y\": S_I},\n    { \"Series\": \"S_S\",      \"Y\": S_S},\n    { \"Series\": \"S_NH\",     \"Y\": S_NH},\n    { \"Series\": \"S_N2\",     \"Y\": S_N2},\n    { \"Series\": \"S_NO\",     \"Y\": S_NO},\n    { \"Series\": \"S_HCO\",    \"Y\": S_HCO}\n    ]};\n\nreturn msg;",
+        "outputs": 1,
+        "timeout": 0,
+        "noerr": 0,
+        "initialize": "",
+        "finalize": "",
+        "libs": [],
+        "x": 1320,
+        "y": 420,
+        "wires": [
+            [
+                "df288814a2a9a2b1"
+            ]
+        ]
+    },
+    {
+        "id": "df288814a2a9a2b1",
+        "type": "ui-chart",
+        "z": "a6b85e226d144df1",
+        "group": "de8b029d69f26c0e",
+        "name": "Effluent",
+        "label": "Effluent",
+        "order": 9007199254740991,
+        "chartType": "line",
+        "category": "Series",
+        "categoryType": "property",
+        "xAxisLabel": "",
+        "xAxisProperty": "",
+        "xAxisPropertyType": "timestamp",
+        "xAxisType": "time",
+        "xAxisFormat": "",
+        "xAxisFormatType": "auto",
+        "xmin": "",
+        "xmax": "",
+        "yAxisLabel": "",
+        "yAxisProperty": "Y",
+        "yAxisPropertyType": "property",
+        "ymin": "",
+        "ymax": "",
+        "bins": 10,
+        "action": "append",
+        "stackSeries": false,
+        "pointShape": "circle",
+        "pointRadius": 4,
+        "showLegend": true,
+        "removeOlder": "8",
+        "removeOlderUnit": "3600",
+        "removeOlderPoints": "2000",
+        "colors": [
+            "#0095ff",
+            "#ff0000",
+            "#ff7f0e",
+            "#2ca02c",
+            "#a347e1",
+            "#d62728",
+            "#ff9896",
+            "#9467bd",
+            "#c5b0d5"
+        ],
+        "textColor": [
+            "#666666"
+        ],
+        "textColorDefault": true,
+        "gridColor": [
+            "#e5e5e5"
+        ],
+        "gridColorDefault": true,
+        "width": 6,
+        "height": 8,
+        "className": "",
+        "interpolation": "linear",
+        "x": 1520,
+        "y": 420,
+        "wires": [
+            []
+        ]
+    },
+    {
+        "id": "ae38454098a37db0",
+        "type": "ui-group",
+        "name": "Group 3",
+        "page": "ca564642bfc5606c",
+        "width": 6,
+        "height": 1,
+        "order": -1,
+        "showTitle": true,
+        "className": "",
+        "visible": "true",
+        "disabled": "false",
+        "groupType": "default"
+    },
+    {
+        "id": "de8b029d69f26c0e",
+        "type": "ui-group",
+        "name": "Group 4",
+        "page": "ca564642bfc5606c",
+        "width": 6,
+        "height": 1,
+        "order": -1,
+        "showTitle": true,
+        "className": "",
+        "visible": "true",
+        "disabled": "false",
+        "groupType": "default"
+    },
+    {
+        "id": "ca564642bfc5606c",
+        "type": "ui-page",
+        "name": "PFR",
+        "ui": "90eb5f47d95b4087",
+        "path": "/page2",
+        "icon": "home",
+        "layout": "grid",
+        "theme": "2c8bcaa0046b4323",
+        "breakpoints": [
+            {
+                "name": "Default",
+                "px": "0",
+                "cols": "3"
+            },
+            {
+                "name": "Tablet",
+                "px": "576",
+                "cols": "6"
+            },
+            {
+                "name": "Small Desktop",
+                "px": "768",
+                "cols": "9"
+            },
+            {
+                "name": "Desktop",
+                "px": "1024",
+                "cols": "12"
+            }
+        ],
+        "order": -1,
+        "className": "",
+        "visible": "true",
+        "disabled": "false"
+    },
+    {
+        "id": "90eb5f47d95b4087",
+        "type": "ui-base",
+        "name": "Dashboard",
+        "path": "/dashboard",
+        "appIcon": "",
+        "includeClientData": true,
+        "acceptsClientConfig": [
+            "ui-notification",
+            "ui-control"
+        ],
+        "showPathInSidebar": false,
+        "headerContent": "page",
+        "navigationStyle": "default",
+        "titleBarStyle": "default",
+        "showReconnectNotification": true,
+        "notificationDisplayTime": 1,
+        "showDisconnectNotification": true,
+        "allowInstall": true
+    },
+    {
+        "id": "2c8bcaa0046b4323",
+        "type": "ui-theme",
+        "name": "Default",
+        "colors": {
+            "surface": "#ffffff",
+            "primary": "#0094ce",
+            "bgPage": "#eeeeee",
+            "groupBg": "#ffffff",
+            "groupOutline": "#cccccc"
+        },
+        "sizes": {
+            "density": "default",
+            "pagePadding": "12px",
+            "groupGap": "12px",
+            "groupBorderRadius": "4px",
+            "widgetGap": "12px"
+        }
+    }
+]

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>
@@ -154,20 +153,18 @@
         <label for="node-input-length"><i class="fa fa-tag"></i> Reactor length [m]</label>
         <input type="text" id="node-input-length" placeholder="m">
     </div>
+    <h2> Simulation parameters </h2>
+    <div class="form-row">
+        <label for="node-input-timeStep"><i class="fa fa-tag"></i> Time step [s]</label>
+        <input type="text" id="node-input-timeStep" placeholder="s">
+    </div>
     <div class="form-row PFR">
-        <label for="node-input-resolution_L"><i class="fa fa-tag"></i> Resolution</label>
+        <label for="node-input-resolution_L"><i class="fa fa-tag"></i> Spatial 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> Δx (length / resolution) [m]</label>
+        <span id="dx-output" style="display: inline-block; padding: 8px; font-weight: bold;">--</span>
     </div>
     <h3> Internal mass transfer calculation (optional) </h3>
     <div class="form-row">
@@ -228,11 +225,6 @@
         <label for="node-input-X_TS_init"><i class="fa fa-tag"></i> Initial total suspended solids [g TSS m-3]</label>
         <input type="text" id="node-input-X_TS_init" class="concentrations">
     </div>
-    <h2> Simulation parameters </h2>
-    <div class="form-row">
-        <label for="node-input-timeStep"><i class="fa fa-tag"></i> Time step [s]</label>
-        <input type="text" id="node-input-timeStep" placeholder="s">
-    </div>
 
     <!-- Logger fields injected here -->
 	<div id="logger-fields-placeholder"></div>
@@ -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 };