RnD/reactor
3
0
Fork 0
Code Issues 1 Pull Requests Actions Packages Projects Releases Wiki Activity

Fixed bug in NaN assertion

Browse Source
This commit is contained in:
p.vanderwilt
2025-07-04 15:48:05 +02:00
parent c6b0cab067
commit 4b49d10763
1 changed files with 30 additions and 22 deletions
Download Patch File Download Diff File Expand all files Collapse all files

52
src/reactor_class.js
View File

@@ -13,8 +13,14 @@ const math = create(all, config)
* @param {string} label
*/
function assertNoNaN(arr, label = "array") {
if (math.isNaN(arr)) {
throw new Error("NaN detected in ${label}!");
if (Array.isArray(arr)) {
for (const el of arr) {
assertNoNaN(el, label);
}
} else {
if (Number.isNaN(arr)) {
throw new Error(`NaN detected in ${label}!`);
}
}
}
@@ -80,7 +86,7 @@ class Reactor {
*/
_arrayClip2Zero(arr) {
if (Array.isArray(arr)) {
return arr.map(clipToZero);
return arr.map(x => this._arrayClip2Zero(x));
} else {
return arr < 0 ? 0 : arr;
}
@@ -166,6 +172,9 @@ class Reactor_PFR extends Reactor {
this.D_op = this._makeDoperator(true, true);
this.D2_op = this._makeD2operator();
assertNoNaN(this.D_op, "Derivative operator");
assertNoNaN(this.D2_op, "Second derivative operator");
}
/**
@@ -184,7 +193,7 @@ class Reactor_PFR extends Reactor {
return { topic: "Fluent", payload: { inlet: 0, F: math.sum(this.Fs), C: this.state.at(-1) }, timestamp: this.currentTime };
}
_applyBoundaryConditions(newState) {
_applyBoundaryConditions(state) {
// 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];
@@ -192,14 +201,13 @@ class Reactor_PFR extends Reactor {
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, [BC_gradient], stateNew)[0];
newState[0] = math.add(BC_C_in, BC_dispersion).map(val => val < 0 ? 0 : val);
const BC_dispersion = math.multiply((1 - (1 + 4 * this.volume / math.sum(this.Fs) / Pe) ^ 0.5) / Pe, [BC_gradient], state)[0];
state[0] = math.add(BC_C_in, BC_dispersion).map(val => val < 0 ? 0 : val);
} else { // Neumann BC (no flux)
newState[0] = newState[1];
state[0] = state[1];
}
// Neumann BC (no flux)
newState[this.n_x - 1] = newState[this.n_x - 2]
return newState
state[this.n_x - 1] = state[this.n_x - 2]
}
/**
@@ -224,11 +232,11 @@ class Reactor_PFR extends Reactor {
}
const dC_total = math.multiply(math.add(dispersion, advection, reaction, transfer), time_step);
let stateNew = math.add(this.state, dC_total);
const stateNew = math.add(this.state, dC_total);
assertNoNaN(stateNew, "new state");
stateNew = this._applyBoundaryConditions(stateNew);
this._applyBoundaryConditions(stateNew);
assertNoNaN(stateNew, "new state post BC");
@@ -245,22 +253,22 @@ class Reactor_PFR extends Reactor {
_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 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, B, C);
const NearBoundary = Array(this.n_x).fill(0.0);
NearBoundary[0] = -1 / 4;
NearBoundary[1] = -5 / 6;
NearBoundary[2] = 3 / 2;
NearBoundary[3] = -1 / 2;
NearBoundary[4] = 1 / 12;
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] = math.multiply(-1, 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);
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.");