From c566766c4de982f6416ca87fd48d9855f7c19390 Mon Sep 17 00:00:00 2001
From: "p.vanderwilt" <p.van.der.wilt@brabantsedelta.nl>
Date: Tue, 8 Jul 2025 15:29:55 +0200
Subject: [PATCH] Refactor boundary condition application in Reactor_PFR class
 for improved clarity and efficiency

---
 src/specificClass.js | 39 ++++++++++++++++++---------------------
 1 file changed, 18 insertions(+), 21 deletions(-)

diff --git a/src/specificClass.js b/src/specificClass.js
index 18af512..1d715a7 100644
--- a/src/specificClass.js
+++ b/src/specificClass.js
@@ -196,27 +196,6 @@ class Reactor_PFR extends Reactor {
         return { topic: "Fluent", payload: { inlet: 0, F: math.sum(this.Fs), C: this.state.at(-1) }, timestamp: this.currentTime };
     }
 
-    /**
-     * 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_gradient = Array(this.n_x).fill(0);
-            BC_gradient[0] = -1;
-            BC_gradient[1] = 1;
-            const BC_dispersion = math.multiply((1 - this.alpha) * this.D*this.A / (math.sum(this.Fs) * this.d_x), [BC_gradient], state)[0];
-            state[0] = math.add(BC_C_in, BC_dispersion).map(val => val < 0 ? 0 : val);
-        } else { // Neumann BC (no flux)
-            state[0] = state[1];
-        }
-        // Neumann BC (no flux)
-        state[this.n_x-1] = state[this.n_x-2]
-    }
-
     /**
      * Tick the reactor state using explicit finite difference method.
      * @param {number} time_step - Time step for the simulation [d].
@@ -251,6 +230,24 @@ class Reactor_PFR extends Reactor {
         return stateNew;
     }
 
+    /**
+     * 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])));
+        } else {
+            state[0] = state[1];
+        }
+        // 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.