Fixed various bugs

dependencies/asm3_class.js

@@ -2,67 +2,65 @@ const math = require('mathjs')
 
 class ASM3 {
 
-  kin_params = {
-    // Kinetic parameters (20 C for now)
-
-    // 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]
-  }
-
-  stoi_params = {
-    // Stoichiometric and composition parameters
-
-    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]
-  }
-
   constructor() {
-    this.stoi_matrix = this._initialise_stoi_matrix()
+    this.kin_params = {
+      // Kinetic parameters (20 C for now)
+
+      // 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.stoi_params = {
+      // Stoichiometric and composition parameters
+
+      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_matrix = this._initialise_stoi_matrix();
   }
 
   _initialise_stoi_matrix() { // initialise stoichiometric matrix

dependencies/reactor_class.js

@@ -1,5 +1,11 @@
 const ASM3 = require('./asm3_class')
-const math = require('mathjs')
+const { create, all } = require('mathjs')
+
+const config = {
+  matrix: 'Array' // Choose 'Matrix' (default) or 'Array'
+}
+
+const math = create(all, config)
 
 class Reactor_CSTR {
     
@@ -93,8 +99,8 @@ class Reactor_PFR {
         this.timeStep = 1/(24*60*15); // time step [d]
         this.speedUpFactor = 1;
 
-        this.D_op = makeDoperator();
+        this.D_op = this.makeDoperator();
-        this.D2_op = makeD2operator();
+        this.D2_op = this.makeD2operator();
     }
 
     set setInfluent(input) { // setter for C_in (WIP)
@@ -139,7 +145,7 @@ class Reactor_PFR {
     tick_fe(time_step) { // tick reactor state using forward Euler method
         const dispersion = math.multiply(this.D / (this.d_x*this.d_x), this.D2_op, this.state);
         const advection = math.multiply(math.sum(this.Fs)/(this.A*this.d_x), this.D_op, this.state);
-        const reaction = this.state.map(this.asm.compute_dC);
+        const reaction = this.state.map((row) => this.asm.compute_dC(row));
         reaction[0] = Array(13).fill(0.0);
         const transfer = Array.from(Array(this.n_x), () => new Array(13).fill(0.0));
 
@@ -154,7 +160,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], this.state);
+            const BC_dispersion = math.multiply(this.D * this.A / (math.sum(this.Fs)*this.d_x), [BC_gradient], this.state)[0];
             this.state[0] = math.add(BC_influx, BC_dispersion);
         }
 
@@ -166,19 +172,19 @@ class Reactor_PFR {
 
     makeDoperator() { // create the upwind scheme gradient operator
         const I = math.identity(this.n_x);
-        const A = math.diag(Array(this.n_x).fill(-1), 1).resize([this.n_x, this.n_x]);
+        const A = math.resize(math.diag(Array(this.n_x).fill(-1), 1), [this.n_x, this.n_x]);
-        I[0, 0] = 0;
+        I[0][0] = 0;
-        I[0, 1] = 0;
+        I[0][1] = 1;
-        I[this.n_x-1, this.n_x-1] = 0; // Neumann boundary condition at x=L
+        I[this.n_x-1][this.n_x-1] = 0; // Neumann boundary condition at x=L
         return math.add(I, A);
     }
 
     makeD2operator() { // create the upwind scheme second derivative operator
-        const I = math.diag(Array(this.n_x).fill(2), 0);
+        const I = math.identity(this.n_x);
-        const A = math.diag(Array(this.n_x).fill(-1), 1).resize([this.n_x, this.n_x]);
+        const A = math.resize(math.diag(Array(this.n_x).fill(-1), 1), [this.n_x, this.n_x]);
-        const B = math.diag(Array(this.n_x).fill(-1), -1).resize([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]);
-        I[0, 0] = 0;
+        I[0][0] = 0;
-        I[0, 1] = 0;
+        I[0][1] = 1;
         return math.add(I, A, B);
     }
 }
@@ -187,9 +193,10 @@ class Reactor_PFR {
 // testing stuff
 // 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_CSTR(initial_state);
+// const Reactor = new Reactor_PFR(200, 10, 10, 1, 100, initial_state);
-// Reactor.C_in = [0.0, 30., 100., 16., 0., 0., 5., 25., 75., 30., 0., 0., 125.];
+// Reactor.Cs_in[0] = [0.0, 30., 100., 16., 0., 0., 5., 25., 75., 30., 0., 0., 125.];
-// N = 0;
+// Reactor.Fs[0] = 10;
+// let N = 0;
 // while (N < 500) {
 //     console.log(Reactor.tick_fe(0.001));
 //     N += 1;