Changed the upper boundary to lower order scheme for now

dependencies/reactor_class.js

@@ -103,8 +103,10 @@ class Reactor_PFR {
         this.timeStep = 1/(24*60*15); // time step [d]
         this.speedUpFactor = 60;
 
-        this.D_op = this.makeDoperator(true);
+        this.D_op = this.makeDoperator(true, true);
         this.D2_op = this.makeD2operator();
+
+        this.alpha = 0.001; // boundary condition modifier
     }
 
     set setInfluent(input) { // setter for C_in (WIP)
@@ -183,7 +185,7 @@ class Reactor_PFR {
             const BC_gradient = Array(this.n_x).fill(0.0);
             BC_gradient[0] = -1;
             BC_gradient[1] = 1;
-            const BC_dispersion = math.multiply(this.D * this.A / (math.sum(this.Fs)*this.d_x), [BC_gradient], new_state)[0];
+            const BC_dispersion = math.multiply(this.alpha*this.D * this.A / (math.sum(this.Fs)*this.d_x), [BC_gradient], new_state)[0];
             console.log(math.add(BC_C_in, BC_dispersion));
             new_state[0] = math.add(BC_C_in, BC_dispersion).map(val => val < 0 ? 0 : val);
             
@@ -208,16 +210,22 @@ class Reactor_PFR {
                 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);
+                const D = math.add(I, A, B, C);
+                D[1][0] = -1;
+                D[1][1] = 0;
+                D[1][2] = 1;
+                D[1][3] = 0;
                 const NearBoundary = Array(this.n_x).fill(0.0);
                 NearBoundary[1] = -25/12;
                 NearBoundary[2] = 4;
                 NearBoundary[3] = -3;
                 NearBoundary[4] = 4/3;
                 NearBoundary[5] = -1/4;
-                D[1] = NearBoundary;
                 NearBoundary.reverse();
-                D[this.n_x-2] = math.multiply(-1, NearBoundary)
+                D[this.n_x-2] = 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.");
             }
@@ -225,11 +233,10 @@ class Reactor_PFR {
             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;
         }
-
-        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;
     }
 
     makeD2operator() { // create the central second derivative operator