From f2d94b26c5141d976dac3e3cf22405fb8703cbd3 Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Tue, 24 Jun 2025 13:28:45 +0200 Subject: [PATCH 01/12] Add dispersion setting in advanced-reactor and initialize axial dispersion to zero in Reactor_PFR --- advanced-reactor.js | 3 +++ dependencies/reactor_class.js | 2 +- 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/advanced-reactor.js b/advanced-reactor.js index 8348137..e2eca51 100644 --- a/advanced-reactor.js +++ b/advanced-reactor.js @@ -78,6 +78,9 @@ module.exports = function(RED) { case "OTR": reactor.setOTR = msg; break; + case "Dispersion": + reactor.setDispersion = msg; + break; default: console.log("Unknown topic: " + msg.topic); } diff --git a/dependencies/reactor_class.js b/dependencies/reactor_class.js index 4027690..e030637 100644 --- a/dependencies/reactor_class.js +++ b/dependencies/reactor_class.js @@ -94,7 +94,7 @@ class Reactor_PFR { this.Fs = Array(n_inlets).fill(0.0); // fluid debits per inlet [m3 d-1] this.Cs_in = Array.from(Array(n_inlets), () => new Array(13).fill(0.0)); // composition influents this.OTR = 0.0; // oxygen transfer rate [g O2 d-1] - this.D = 0.1; // axial dispersion [m2 d-1] + this.D = 0.0; // axial dispersion [m2 d-1] this.kla = kla; // if NaN, use external OTR [d-1] From 2e76f733a8fed39a17f3f699a2c6b0506017a315 Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Tue, 24 Jun 2025 16:23:33 +0200 Subject: [PATCH 02/12] Working on fixing the Derivative operators and BCs --- dependencies/reactor_class.js | 40 ++++++++++++++++++++--------------- 1 file changed, 23 insertions(+), 17 deletions(-) diff --git a/dependencies/reactor_class.js b/dependencies/reactor_class.js index e030637..9e02fca 100644 --- a/dependencies/reactor_class.js +++ b/dependencies/reactor_class.js @@ -69,7 +69,8 @@ class Reactor_CSTR { const dC_total = math.multiply(math.add(dC_in, dC_out, r, t_O), time_step); - this.state = math.abs(math.add(this.state, dC_total)); // make sure that concentrations do not go negative + // clip value element-wise to each subarray to avoid negative concentrations + this.state = math.add(this.state, dC_total).map(val => val < 0 ? 0 : val); return this.state; } } @@ -99,8 +100,8 @@ class Reactor_PFR { this.kla = kla; // if NaN, use external OTR [d-1] this.currentTime = Date.now(); // milliseconds since epoch [ms] - this.timeStep = 1/(24*60*60); // time step [d] - this.speedUpFactor = 1; + this.timeStep = 1/(24*60*15); // time step [d] + this.speedUpFactor = 60; this.D_op = this.makeDoperator(); this.D2_op = this.makeD2operator(); @@ -147,7 +148,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 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)); reaction[0] = Array(13).fill(0.0); const transfer = Array.from(Array(this.n_x), () => new Array(13).fill(0.0)); @@ -158,36 +159,40 @@ class Reactor_PFR { transfer.forEach((x, i) => { x[0] = this.calcOTR(this.state[i][0]); }); } + transfer[0][0] = 0; + if (math.sum(this.Fs) > 0) { // Danckwerts BC - const BC_influx = math.multiply(math.divide([this.Fs], this.A), this.Cs_in)[0]; + 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.0); - BC_gradient[0] = 1; - BC_gradient[1] = -1; + 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)[0]; - this.state[0] = math.add(BC_influx, BC_dispersion); + this.state[0] = math.add(BC_C_in, BC_dispersion); + console.log(BC_dispersion); } const dC_total = math.multiply(math.add(dispersion, advection, reaction, transfer), time_step); - this.state = math.abs(math.add(this.state, dC_total)); // make sure that concentrations do not go negative + // clip value element-wise to each subarray to avoid negative concentrations + this.state = math.add(this.state, dC_total).map(row => row.map(val => val < 0 ? 0 : val)); return this.state; } makeDoperator() { // create the upwind scheme gradient operator - const I = math.identity(this.n_x); - const A = math.resize(math.diag(Array(this.n_x).fill(-1), 1), [this.n_x, this.n_x]); + const I = math.diag(Array(this.n_x).fill(-1), 0); + 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][1] = 1; + I[0][1] = -1; 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.identity(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.resize(math.diag(Array(this.n_x).fill(-1), -1), [this.n_x, this.n_x]); + 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]); I[0][0] = 0; - I[0][1] = 1; + I[0][1] = -1; return math.add(I, A, B); } } @@ -199,8 +204,9 @@ class Reactor_PFR { // 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 < 500) { +// while (N < 5000) { // console.log(Reactor.tick_fe(0.001)); // N += 1; // } From 9f1322978535e10c57eb83872d50cb416ce713bf Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Tue, 24 Jun 2025 16:38:07 +0200 Subject: [PATCH 03/12] Fix boundary conditions in gradient and second derivative operators for Reactor_PFR --- dependencies/reactor_class.js | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/dependencies/reactor_class.js b/dependencies/reactor_class.js index 9e02fca..adf744f 100644 --- a/dependencies/reactor_class.js +++ b/dependencies/reactor_class.js @@ -179,20 +179,20 @@ class Reactor_PFR { } makeDoperator() { // create the upwind scheme gradient operator - const I = math.diag(Array(this.n_x).fill(-1), 0); - const A = math.resize(math.diag(Array(this.n_x).fill(1), 1), [this.n_x, this.n_x]); + const I = math.diag(Array(this.n_x).fill(1), 0); + 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][1] = -1; - 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 + 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]); I[0][0] = 0; I[0][1] = -1; + I[0][0] = -1; // Dichelet boundary condition at outlet return math.add(I, A, B); } } From bb74fc86c212e8545a79e485146c1beaae0174d5 Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Fri, 27 Jun 2025 16:56:37 +0200 Subject: [PATCH 04/12] Refactor dispersion and boundary condition handling in Reactor_PFR --- dependencies/reactor_class.js | 49 +++++++++++++++++++---------------- 1 file changed, 27 insertions(+), 22 deletions(-) diff --git a/dependencies/reactor_class.js b/dependencies/reactor_class.js index adf744f..94844ee 100644 --- a/dependencies/reactor_class.js +++ b/dependencies/reactor_class.js @@ -103,7 +103,7 @@ class Reactor_PFR { this.timeStep = 1/(24*60*15); // time step [d] this.speedUpFactor = 60; - this.D_op = this.makeDoperator(); + this.D_op = this.makeDoperator(false); this.D2_op = this.makeD2operator(); } @@ -111,6 +111,8 @@ class Reactor_PFR { let index_in = input.payload.inlet; this.Fs[index_in] = input.payload.F; this.Cs_in[index_in] = input.payload.C; + // console.log("Pe " + this.d_x*math.sum(this.Fs)/(this.D*this.A)); + // console.log("Co " + math.sum(this.Fs)*this.timeStep/(this.A*this.d_x)); } set setOTR(input) { // setter for OTR (WIP) [g O2 d-1] @@ -150,7 +152,6 @@ class Reactor_PFR { 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)); - reaction[0] = Array(13).fill(0.0); const transfer = Array.from(Array(this.n_x), () => new Array(13).fill(0.0)); if (isNaN(this.kla)) { // calculate OTR if kla is not NaN, otherwise use externally calculated OTR @@ -159,41 +160,45 @@ class Reactor_PFR { transfer.forEach((x, i) => { x[0] = this.calcOTR(this.state[i][0]); }); } - transfer[0][0] = 0; + const dC_total = math.multiply(math.add(dispersion, advection, reaction, transfer), time_step); + const new_state = math.add(this.state, dC_total); + // apply boundary conditions 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_C_in = math.multiply(1/math.sum(this.Fs), [this.Fs], this.Cs_in)[0]; 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)[0]; - this.state[0] = math.add(BC_C_in, BC_dispersion); - console.log(BC_dispersion); + const BC_dispersion = math.multiply(this.D * this.A / (math.sum(this.Fs)*this.d_x), [BC_gradient], new_state)[0]; + new_state[0] = math.add(BC_C_in, BC_dispersion).map(val => val < 0 ? 0 : val); + console.log(new_state[0]) + } else { // Neumann BC (no flux) + new_state[0] = new_state[1]; } + // Neumann BC (no flux) + new_state[this.n_x-1] = new_state[this.n_x-2] - const dC_total = math.multiply(math.add(dispersion, advection, reaction, transfer), time_step); - - // clip value element-wise to each subarray to avoid negative concentrations - this.state = math.add(this.state, dC_total).map(row => row.map(val => val < 0 ? 0 : val)); - return this.state; + this.state = new_state.map(row => row.map(val => val < 0 ? 0 : val)); // apply the new state + return new_state; } - makeDoperator() { // create the upwind scheme gradient operator - const I = math.diag(Array(this.n_x).fill(1), 0); + makeDoperator(central=false) { // create the upwind scheme gradient operator + const I = math.resize(math.diag(Array(this.n_x).fill(1), central), [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][1] = -1; - return math.add(I, A); + 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; } 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]); - I[0][0] = 0; - I[0][1] = -1; - I[0][0] = -1; // Dichelet boundary condition at outlet - return math.add(I, A, B); + 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); + return D2; } } @@ -211,4 +216,4 @@ class Reactor_PFR { // N += 1; // } -module.exports = {Reactor_CSTR, Reactor_PFR}; \ No newline at end of file +module.exports = { Reactor_CSTR, Reactor_PFR }; \ No newline at end of file From 0cc653800324a0990508ebe38799c2f6a67ddf26 Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Fri, 27 Jun 2025 17:29:20 +0200 Subject: [PATCH 05/12] Handle division by zero in rate calculations for ASM3 --- dependencies/asm3_class.js | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/dependencies/asm3_class.js b/dependencies/asm3_class.js index 9b7762a..6f1e5fe 100644 --- a/dependencies/asm3_class.js +++ b/dependencies/asm3_class.js @@ -104,8 +104,8 @@ class ASM3 { // Heterotrophs rates[1] = k_STO * this._monod(S_O, K_O) * this._monod(S_S, K_S) * X_H; rates[2] = k_STO * nu_NO * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * this._monod(S_S, K_S) * X_H; - rates[3] = mu_H_max * this._monod(S_O, K_O) * this._monod(S_NH, K_NH) * this._monod(S_HCO, K_HCO) * this._monod(X_STO/X_H, K_STO) * X_H; - rates[4] = mu_H_max * nu_NO * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * this._monod(S_NH, K_NH) * this._monod(S_HCO, K_HCO) * this._monod(X_STO/X_H, K_STO) * X_H; + rates[3] = X_H == 0 ? 0 : mu_H_max * this._monod(S_O, K_O) * this._monod(S_NH, K_NH) * this._monod(S_HCO, K_HCO) * this._monod(X_STO/X_H, K_STO) * X_H; + rates[4] = X_H == 0 ? 0 : mu_H_max * nu_NO * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * this._monod(S_NH, K_NH) * this._monod(S_HCO, K_HCO) * this._monod(X_STO/X_H, K_STO) * X_H; rates[5] = b_H_O * this._monod(S_O, K_O) * X_H; rates[6] = b_H_NO * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * X_H; rates[7] = b_STO_O * this._monod(S_O, K_O) * X_H; From 8215c5ed9a6a736eed45a47dd8bbd6dd1772e888 Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Sat, 28 Jun 2025 19:19:38 +0200 Subject: [PATCH 06/12] Add checks for NaN values in Reactor_PFR calculations and update hydrolysis rate calculation to handle division by zero --- dependencies/asm3_class.js | 2 +- dependencies/reactor_class.js | 33 ++++++++++++++++++++++++++------- 2 files changed, 27 insertions(+), 8 deletions(-) diff --git a/dependencies/asm3_class.js b/dependencies/asm3_class.js index 6f1e5fe..f5ac1df 100644 --- a/dependencies/asm3_class.js +++ b/dependencies/asm3_class.js @@ -99,7 +99,7 @@ class ASM3 { const { k_H, K_X, k_STO, nu_NO, K_O, K_NO, K_S, K_STO, mu_H_max, K_NH, K_HCO, b_H_O, b_H_NO, b_STO_O, b_STO_NO, mu_A_max, K_A_NH, K_A_O, K_A_HCO, b_A_O, b_A_NO } = this.kin_params; // Hydrolysis - rates[0] = k_H * this._monod(X_S / X_H, K_X) * X_H; + rates[0] = X_H == 0 ? 0 : k_H * this._monod(X_S / X_H, K_X) * X_H; // Heterotrophs rates[1] = k_STO * this._monod(S_O, K_O) * this._monod(S_S, K_S) * X_H; diff --git a/dependencies/reactor_class.js b/dependencies/reactor_class.js index 94844ee..8212aaf 100644 --- a/dependencies/reactor_class.js +++ b/dependencies/reactor_class.js @@ -103,7 +103,7 @@ class Reactor_PFR { this.timeStep = 1/(24*60*15); // time step [d] this.speedUpFactor = 60; - this.D_op = this.makeDoperator(false); + this.D_op = this.makeDoperator(true); this.D2_op = this.makeD2operator(); } @@ -111,8 +111,10 @@ class Reactor_PFR { let index_in = input.payload.inlet; this.Fs[index_in] = input.payload.F; this.Cs_in[index_in] = input.payload.C; - // console.log("Pe " + this.d_x*math.sum(this.Fs)/(this.D*this.A)); - // console.log("Co " + math.sum(this.Fs)*this.timeStep/(this.A*this.d_x)); + console.log("Pe total " + this.length*math.sum(this.Fs)/(this.D*this.A)); + console.log("Pe local " + this.d_x*math.sum(this.Fs)/(this.D*this.A)); + console.log("Co ad " + math.sum(this.Fs)*this.timeStep/(this.A*this.d_x)); + console.log("Co D " + this.D*this.timeStep/(this.d_x*this.d_x)); } set setOTR(input) { // setter for OTR (WIP) [g O2 d-1] @@ -134,7 +136,6 @@ class Reactor_PFR { // expect update with timestamp updateState(newTime) { - const day2ms = 1000 * 60 * 60 * 24; let n_iter = Math.floor(this.speedUpFactor*(newTime - this.currentTime) / (this.timeStep * day2ms)); @@ -153,6 +154,16 @@ class Reactor_PFR { 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)); const transfer = Array.from(Array(this.n_x), () => new Array(13).fill(0.0)); + + if (dispersion.some(row => row.some(Number.isNaN))) { + throw new Error("NaN detected in dispersion!"); + } + if (advection.some(row => row.some(Number.isNaN))) { + throw new Error("NaN detected in advection!"); + } + if (reaction.some(row => row.some(Number.isNaN))) { + throw new Error("NaN detected in reaction!"); + } if (isNaN(this.kla)) { // calculate OTR if kla is not NaN, otherwise use externally calculated OTR transfer.forEach((x) => { x[0] = this.OTR; }); @@ -162,6 +173,9 @@ class Reactor_PFR { const dC_total = math.multiply(math.add(dispersion, advection, reaction, transfer), time_step); const new_state = math.add(this.state, dC_total); + if (new_state.some(row => row.some(Number.isNaN))) { + throw new Error("NaN detected in new_state after dC_total update!"); + } // apply boundary conditions if (math.sum(this.Fs) > 0) { // Danckwerts BC @@ -170,21 +184,26 @@ class Reactor_PFR { 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]; + 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); - console.log(new_state[0]) + } else { // Neumann BC (no flux) new_state[0] = new_state[1]; } // Neumann BC (no flux) new_state[this.n_x-1] = new_state[this.n_x-2] + if (new_state.some(row => row.some(Number.isNaN))) { + throw new Error("NaN detected in new_state after enforcing boundary conditions!"); + } + this.state = new_state.map(row => row.map(val => val < 0 ? 0 : val)); // apply the new state return new_state; } makeDoperator(central=false) { // create the upwind scheme gradient operator - const I = math.resize(math.diag(Array(this.n_x).fill(1), central), [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 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); From b2d32ba9f2011bbc878a0aa46db1925e5ed69ad4 Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Mon, 30 Jun 2025 12:50:02 +0200 Subject: [PATCH 07/12] Enhance makeDoperator to support higher-order central gradient schemes and improve boundary handling --- dependencies/reactor_class.js | 30 ++++++++++++++++++++++++++---- 1 file changed, 26 insertions(+), 4 deletions(-) diff --git a/dependencies/reactor_class.js b/dependencies/reactor_class.js index 8212aaf..e80f88b 100644 --- a/dependencies/reactor_class.js +++ b/dependencies/reactor_class.js @@ -201,10 +201,32 @@ class Reactor_PFR { return new_state; } - makeDoperator(central=false) { // create the upwind scheme gradient operator - 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); + 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); + 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) + } 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; From 3cc876533c64a2ac5b4dd0043af0aaf6b43e01f2 Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Mon, 30 Jun 2025 15:46:13 +0200 Subject: [PATCH 08/12] Changed the upper boundary to lower order scheme for now --- dependencies/reactor_class.js | 25 ++++++++++++++++--------- 1 file changed, 16 insertions(+), 9 deletions(-) diff --git a/dependencies/reactor_class.js b/dependencies/reactor_class.js index e80f88b..0d73705 100644 --- a/dependencies/reactor_class.js +++ b/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 From f4824b822cac7208cd87bd5fbb12a9c182f674a8 Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Tue, 1 Jul 2025 13:04:32 +0200 Subject: [PATCH 09/12] Improved wieghted finite differencing --- dependencies/reactor_class.js | 17 +++++++---------- 1 file changed, 7 insertions(+), 10 deletions(-) diff --git a/dependencies/reactor_class.js b/dependencies/reactor_class.js index 0d73705..2c8d3b8 100644 --- a/dependencies/reactor_class.js +++ b/dependencies/reactor_class.js @@ -211,18 +211,15 @@ class Reactor_PFR { 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); - 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; + 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] = NearBoundary; + 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; From e9847607e85d8f6776dae255337dfb04879f753d Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Tue, 1 Jul 2025 16:08:35 +0200 Subject: [PATCH 10/12] Use Generalized boundary condition by Nauman and Mallikarjun 1983 --- dependencies/reactor_class.js | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/dependencies/reactor_class.js b/dependencies/reactor_class.js index 2c8d3b8..6dd6168 100644 --- a/dependencies/reactor_class.js +++ b/dependencies/reactor_class.js @@ -2,7 +2,7 @@ const ASM3 = require('./asm3_class') const { create, all } = require('mathjs') const config = { - matrix: 'Array' // Choose 'Matrix' (default) or 'Array' + matrix: 'Array' // choose 'Matrix' (default) or 'Array' } const math = create(all, config) @@ -185,7 +185,8 @@ 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.alpha*this.D * this.A / (math.sum(this.Fs)*this.d_x), [BC_gradient], new_state)[0]; + let Pe = this.length*math.sum(this.Fs)/(this.D*this.A) + const BC_dispersion = math.multiply((1-(1+4*this.volume/math.sum(this.Fs)/Pe)^0.5)/Pe, 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); From dcc8562dbf5ac2800d6ce4bc10a2b6ff8f88cc1b Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Wed, 2 Jul 2025 10:37:02 +0200 Subject: [PATCH 11/12] Remove depreciated variable --- dependencies/reactor_class.js | 2 -- 1 file changed, 2 deletions(-) diff --git a/dependencies/reactor_class.js b/dependencies/reactor_class.js index 6dd6168..652ee26 100644 --- a/dependencies/reactor_class.js +++ b/dependencies/reactor_class.js @@ -105,8 +105,6 @@ class Reactor_PFR { 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) From f517b7764d083b49f8fac052a2f81e29da0325ab Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Thu, 3 Jul 2025 22:28:34 +0200 Subject: [PATCH 12/12] Remove mistake boundary condition --- dependencies/reactor_class.js | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dependencies/reactor_class.js b/dependencies/reactor_class.js index 652ee26..76a6340 100644 --- a/dependencies/reactor_class.js +++ b/dependencies/reactor_class.js @@ -184,7 +184,7 @@ class Reactor_PFR { BC_gradient[0] = -1; BC_gradient[1] = 1; let Pe = this.length*math.sum(this.Fs)/(this.D*this.A) - const BC_dispersion = math.multiply((1-(1+4*this.volume/math.sum(this.Fs)/Pe)^0.5)/Pe, this.D * this.A / (math.sum(this.Fs)*this.d_x), [BC_gradient], new_state)[0]; + const BC_dispersion = math.multiply((1-(1+4*this.volume/math.sum(this.Fs)/Pe)^0.5)/Pe, [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);