From b2d32ba9f2011bbc878a0aa46db1925e5ed69ad4 Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Mon, 30 Jun 2025 12:50:02 +0200 Subject: [PATCH] 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;