p.vanderwilt/reactor
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Add generalFunctions dependency and implement basic measurement child registration in nodeClass

Browse Source
This commit is contained in:
p.vanderwilt
2025-07-11 12:22:36 +02:00
parent 6fd86f71c8
commit 24de5a4c9f
4 changed files with 294 additions and 264 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
package-lock.json generated
View File

@@ -9,6 +9,7 @@
"version": "0.0.1", "version": "0.0.1",
"license": "SEE LICENSE", "license": "SEE LICENSE",
"dependencies": { "dependencies": {
"generalFunctions": "git+https://gitea.centraal.wbd-rd.nl/p.vanderwilt/generalFunctions.git",
"mathjs": "^14.5.2" "mathjs": "^14.5.2"
} }
}, },
@@ -59,6 +60,11 @@
"url": "https://github.com/sponsors/rawify" "url": "https://github.com/sponsors/rawify"
} }
}, },
"node_modules/generalFunctions": {
"version": "1.0.0",
"resolved": "git+https://gitea.centraal.wbd-rd.nl/p.vanderwilt/generalFunctions.git#950ca2b6b4e91b37479aee90bff74b02c16f130e",
"license": "SEE LICENSE"
},
"node_modules/javascript-natural-sort": { "node_modules/javascript-natural-sort": {
"version": "0.7.1", "version": "0.7.1",
"resolved": "https://registry.npmjs.org/javascript-natural-sort/-/javascript-natural-sort-0.7.1.tgz", "resolved": "https://registry.npmjs.org/javascript-natural-sort/-/javascript-natural-sort-0.7.1.tgz",

1
package.json
View File

@@ -27,6 +27,7 @@
} }
}, },
"dependencies": { "dependencies": {
"generalFunctions": "git+https://gitea.centraal.wbd-rd.nl/p.vanderwilt/generalFunctions.git",
"mathjs": "^14.5.2" "mathjs": "^14.5.2"
} }
} }

6
src/nodeClass.js
View File

@@ -48,6 +48,12 @@ class nodeClass {
case "Dispersion": case "Dispersion":
this.reactor.setDispersion = msg; this.reactor.setDispersion = msg;
break; break;
case 'registerChild':
// Register this node as a child of the parent node
const childId = msg.payload;
const childObj = this.RED.nodes.getNode(childId);
this.reactor.childRegistrationUtils.registerChild(childObj.source, msg.positionVsParent);
break;
default: default:
console.log("Unknown topic: " + msg.topic); console.log("Unknown topic: " + msg.topic);
} }

545
src/specificClass.js
View File

@@ -1,9 +1,10 @@
const ASM3 = require('./reaction_modules/asm3_class.js'); const ASM3 = require('./reaction_modules/asm3_class.js');
const { create, all } = require('mathjs'); const { create, all } = require('mathjs');
const { assertNoNaN } = require('./utils.js'); const { assertNoNaN } = require('./utils.js');
const { childRegistrationUtils, logger, MeasurementContainer } = require('generalFunctions');
const config = { const config = {
matrix: 'Array' // use Array as the matrix type matrix: 'Array' // use Array as the matrix type
}; };
const math = create(all, config); const math = create(all, config);
@@ -13,304 +14,320 @@ const NUM_SPECIES = 13;
const DEBUG = false; const DEBUG = false;
class Reactor { class Reactor {
/** /**
* Reactor base class. * Reactor base class.
* @param {object} config - Configuration object containing reactor parameters. * @param {object} config - Configuration object containing reactor parameters.
*/ */
constructor(config) { constructor(config) {
this.asm = new ASM3(); // EVOLV stuff
this.logger = new logger(); //TODO: attach config
this.measurements = new MeasurementContainer();
this.childRegistrationUtils = new childRegistrationUtils(this); // Child registration utility
this.volume = config.volume; // fluid volume reactor [m3] this.asm = new ASM3();
this.Fs = Array(config.n_inlets).fill(0); // fluid debits per inlet [m3 d-1] this.volume = config.volume; // fluid volume reactor [m3]
this.Cs_in = Array.from(Array(config.n_inlets), () => new Array(NUM_SPECIES).fill(0)); // composition influents
this.OTR = 0.0; // oxygen transfer rate [g O2 d-1]
this.temperature = 20; // temperature [C]
this.kla = config.kla; // if NaN, use externaly provided OTR [d-1] this.Fs = Array(config.n_inlets).fill(0); // fluid debits per inlet [m3 d-1]
this.Cs_in = Array.from(Array(config.n_inlets), () => new Array(NUM_SPECIES).fill(0)); // composition influents
this.OTR = 0.0; // oxygen transfer rate [g O2 d-1]
this.temperature = 20; // temperature [C]
this.currentTime = Date.now(); // milliseconds since epoch [ms] this.kla = config.kla; // if NaN, use externaly provided OTR [d-1]
this.timeStep = 1 / (24*60*15); // time step [d]
this.speedUpFactor = 60; // speed up factor for simulation, 60 means 1 minute per simulated second this.currentTime = Date.now(); // milliseconds since epoch [ms]
this.timeStep = 1 / (24*60*15); // time step [d]
this.speedUpFactor = 60; // speed up factor for simulation, 60 means 1 minute per simulated second
}
updateMeasurement(variant, subType, value, position) {
this.logger.debug(`---------------------- updating ${subType} ------------------ `);
switch (subType) {
case "temperature":
this.logger.debug(`no nothing`);
break;
default:
this.logger.error(`Type '${subType}' not recognized for measured update.`);
return;
} }
}
/** /**
* Setter for influent data. * Setter for influent data.
* @param {object} input - Input object (msg) containing payload with inlet index, flow rate, and concentrations. * @param {object} input - Input object (msg) containing payload with inlet index, flow rate, and concentrations.
*/ */
set setInfluent(input) { set setInfluent(input) {
let index_in = input.payload.inlet; let index_in = input.payload.inlet;
this.Fs[index_in] = input.payload.F; this.Fs[index_in] = input.payload.F;
this.Cs_in[index_in] = input.payload.C; this.Cs_in[index_in] = input.payload.C;
}
/**
* Setter for OTR (Oxygen Transfer Rate).
* @param {object} input - Input object (msg) containing payload with OTR value [g O2 d-1].
*/
set setOTR(input) {
this.OTR = input.payload;
}
/**
* Setter for temperature.
* @param {object} input - Input object (msg) containing payload with temperature value [C].
*/
set setTemperature(input) {
this.temperature = input.payload;
}
/**
* Calculate the oxygen transfer rate (OTR) based on the dissolved oxygen concentration and temperature.
* @param {number} S_O - Dissolved oxygen concentration [g O2 m-3].
* @param {number} T - Temperature in Celsius, default to 20 C.
* @returns {number} - Calculated OTR [g O2 d-1].
*/
_calcOTR(S_O, T = 20.0) { // caculate the OTR using basic correlation, default to temperature: 20 C
let S_O_sat = 14.652 - 4.1022e-1 * T + 7.9910e-3 * T*T + 7.7774e-5 * T*T*T;
return this.kla * (S_O_sat - S_O);
}
/**
* Clip values in an array to zero.
* @param {Array} arr - Array of values to clip.
* @returns {Array} - New array with values clipped to zero.
*/
_arrayClip2Zero(arr) {
if (Array.isArray(arr)) {
return arr.map(x => this._arrayClip2Zero(x));
} else {
return arr < 0 ? 0 : arr;
} }
}
/** /**
* Setter for OTR (Oxygen Transfer Rate). * Update the reactor state based on the new time.
* @param {object} input - Input object (msg) containing payload with OTR value [g O2 d-1]. * @param {number} newTime - New time to update reactor state to, in milliseconds since epoch.
*/ */
set setOTR(input) { updateState(newTime) { // expect update with timestamp
this.OTR = input.payload; const day2ms = 1000 * 60 * 60 * 24;
let n_iter = Math.floor(this.speedUpFactor * (newTime-this.currentTime) / (this.timeStep*day2ms));
if (n_iter) {
let n = 0;
while (n < n_iter) {
this.tick(this.timeStep);
n += 1;
} }
this.currentTime += n_iter * this.timeStep * day2ms / this.speedUpFactor;
/**
* Setter for temperature.
* @param {object} input - Input object (msg) containing payload with temperature value [C].
*/
set setTemperature(input) {
this.temperature = input.payload;
} }
}
/**
* Calculate the oxygen transfer rate (OTR) based on the dissolved oxygen concentration and temperature.
* @param {number} S_O - Dissolved oxygen concentration [g O2 m-3].
* @param {number} T - Temperature in Celsius, default to 20 C.
* @returns {number} - Calculated OTR [g O2 d-1].
*/
_calcOTR(S_O, T = 20.0) { // caculate the OTR using basic correlation, default to temperature: 20 C
let S_O_sat = 14.652 - 4.1022e-1 * T + 7.9910e-3 * T*T + 7.7774e-5 * T*T*T;
return this.kla * (S_O_sat - S_O);
}
/**
* Clip values in an array to zero.
* @param {Array} arr - Array of values to clip.
* @returns {Array} - New array with values clipped to zero.
*/
_arrayClip2Zero(arr) {
if (Array.isArray(arr)) {
return arr.map(x => this._arrayClip2Zero(x));
} else {
return arr < 0 ? 0 : arr;
}
}
/**
* Update the reactor state based on the new time.
* @param {number} newTime - New time to update reactor state to, in milliseconds since epoch.
*/
updateState(newTime) { // expect update with timestamp
const day2ms = 1000 * 60 * 60 * 24;
let n_iter = Math.floor(this.speedUpFactor * (newTime-this.currentTime) / (this.timeStep*day2ms));
if (n_iter) {
let n = 0;
while (n < n_iter) {
this.tick(this.timeStep);
n += 1;
}
this.currentTime += n_iter * this.timeStep * day2ms / this.speedUpFactor;
}
}
} }
class Reactor_CSTR extends Reactor { class Reactor_CSTR extends Reactor {
/** /**
* Reactor_CSTR class for Continuous Stirred Tank Reactor. * Reactor_CSTR class for Continuous Stirred Tank Reactor.
* @param {object} config - Configuration object containing reactor parameters. * @param {object} config - Configuration object containing reactor parameters.
*/ */
constructor(config) { constructor(config) {
super(config); super(config);
this.state = config.initialState; this.state = config.initialState;
} }
/** /**
* Getter for effluent data. * Getter for effluent data.
* @returns {object} Effluent data object (msg), defaults to inlet 0. * @returns {object} Effluent data object (msg), defaults to inlet 0.
*/ */
get getEffluent() { // getter for Effluent, defaults to inlet 0 get getEffluent() { // getter for Effluent, defaults to inlet 0
return { topic: "Fluent", payload: { inlet: 0, F: math.sum(this.Fs), C: this.state }, timestamp: this.currentTime }; return { topic: "Fluent", payload: { inlet: 0, F: math.sum(this.Fs), C: this.state }, timestamp: this.currentTime };
} }
/** /**
* Tick the reactor state using the forward Euler method. * Tick the reactor state using the forward Euler method.
* @param {number} time_step - Time step for the simulation [d]. * @param {number} time_step - Time step for the simulation [d].
* @returns {Array} - New reactor state. * @returns {Array} - New reactor state.
*/ */
tick(time_step) { // tick reactor state using forward Euler method tick(time_step) { // tick reactor state using forward Euler method
const inflow = math.multiply(math.divide([this.Fs], this.volume), this.Cs_in)[0]; const inflow = math.multiply(math.divide([this.Fs], this.volume), this.Cs_in)[0];
const outflow = math.multiply(-1 * math.sum(this.Fs) / this.volume, this.state); const outflow = math.multiply(-1 * math.sum(this.Fs) / this.volume, this.state);
const reaction = this.asm.compute_dC(this.state, this.temperature); const reaction = this.asm.compute_dC(this.state, this.temperature);
const transfer = Array(NUM_SPECIES).fill(0.0); const transfer = Array(NUM_SPECIES).fill(0.0);
transfer[S_O_INDEX] = isNaN(this.kla) ? this.OTR : this._calcOTR(this.state[S_O_INDEX], this.temperature); // calculate OTR if kla is not NaN, otherwise use externaly calculated OTR transfer[S_O_INDEX] = isNaN(this.kla) ? this.OTR : this._calcOTR(this.state[S_O_INDEX], this.temperature); // calculate OTR if kla is not NaN, otherwise use externaly calculated OTR
const dC_total = math.multiply(math.add(inflow, outflow, reaction, transfer), time_step) const dC_total = math.multiply(math.add(inflow, outflow, reaction, transfer), time_step)
this.state = this._arrayClip2Zero(math.add(this.state, dC_total)); // clip value element-wise to avoid negative concentrations this.state = this._arrayClip2Zero(math.add(this.state, dC_total)); // clip value element-wise to avoid negative concentrations
if(DEBUG){ if(DEBUG){
assertNoNaN(dC_total, "change in state"); assertNoNaN(dC_total, "change in state");
assertNoNaN(this.state, "new state"); assertNoNaN(this.state, "new state");
} }
return this.state; return this.state;
} }
} }
class Reactor_PFR extends Reactor { class Reactor_PFR extends Reactor {
/** /**
* Reactor_PFR class for Plug Flow Reactor. * Reactor_PFR class for Plug Flow Reactor.
* @param {object} config - Configuration object containing reactor parameters. * @param {object} config - Configuration object containing reactor parameters.
*/ */
constructor(config) { constructor(config) {
super(config); super(config);
this.length = config.length; // reactor length [m] this.length = config.length; // reactor length [m]
this.n_x = config.resolution_L; // number of slices this.n_x = config.resolution_L; // number of slices
this.d_x = this.length / this.n_x; this.d_x = this.length / this.n_x;
this.A = this.volume / this.length; // crosssectional area [m2] this.A = this.volume / this.length; // crosssectional area [m2]
this.alpha = config.alpha; this.alpha = config.alpha;
this.state = Array.from(Array(this.n_x), () => config.initialState.slice()) this.state = Array.from(Array(this.n_x), () => config.initialState.slice())
// console.log("Initial State: ") // console.log("Initial State: ")
// console.log(this.state) // console.log(this.state)
this.D = 0.0; // axial dispersion [m2 d-1] this.D = 0.0; // axial dispersion [m2 d-1]
this.D_op = this._makeDoperator(true, true); this.D_op = this._makeDoperator(true, true);
assertNoNaN(this.D_op, "Derivative operator"); assertNoNaN(this.D_op, "Derivative operator");
this.D2_op = this._makeD2operator(); this.D2_op = this._makeD2operator();
assertNoNaN(this.D2_op, "Second derivative operator"); assertNoNaN(this.D2_op, "Second derivative operator");
}
/**
* Setter for axial dispersion.
* @param {object} input - Input object (msg) containing payload with dispersion value [m2 d-1].
*/
set setDispersion(input) {
this.D = input.payload;
}
/**
* Getter for effluent data.
* @returns {object} Effluent data object (msg), defaults to inlet 0.
*/
get getEffluent() {
return { topic: "Fluent", payload: { inlet: 0, F: math.sum(this.Fs), C: this.state.at(-1) }, timestamp: this.currentTime };
}
updateState(newTime) {
super.updateState(newTime);
let Pe_local = this.d_x*math.sum(this.Fs)/(this.D*this.A)
let Co_D = this.D*this.timeStep/(this.d_x*this.d_x);
(Pe_local >= 2) && console.warn(`Local Péclet number (${Pe_local}) is too high! Increase reactor resolution.`);
(Co_D >= 0.5) && console.warn(`Courant number (${Co_D}) is too high! Reduce time step size.`);
if(DEBUG) {
console.log("Inlet state max " + math.max(this.state[0]))
console.log("Pe total " + this.length*math.sum(this.Fs)/(this.D*this.A));
console.log("Pe local " + Pe_local);
console.log("Co ad " + math.sum(this.Fs)*this.timeStep/(this.A*this.d_x));
console.log("Co D " + Co_D);
}
}
/**
* Tick the reactor state using explicit finite difference method.
* @param {number} time_step - Time step for the simulation [d].
* @returns {Array} - New reactor state.
*/
tick(time_step) {
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, this.temperature));
const transfer = Array.from(Array(this.n_x), () => new Array(NUM_SPECIES).fill(0));
if (isNaN(this.kla)) { // calculate OTR if kla is not NaN, otherwise use externally calculated OTR
transfer.forEach((x) => { x[S_O_INDEX] = this.OTR; });
} else {
transfer.forEach((x, i) => { x[S_O_INDEX] = this._calcOTR(this.state[i][S_O_INDEX], this.temperature); });
} }
/** const dC_total = math.multiply(math.add(dispersion, advection, reaction, transfer), time_step);
* Setter for axial dispersion.
* @param {object} input - Input object (msg) containing payload with dispersion value [m2 d-1]. const stateNew = math.add(this.state, dC_total);
*/ this._applyBoundaryConditions(stateNew);
set setDispersion(input) {
this.D = input.payload; if (DEBUG) {
assertNoNaN(dispersion, "dispersion");
assertNoNaN(advection, "advection");
assertNoNaN(reaction, "reaction");
assertNoNaN(dC_total, "change in state");
assertNoNaN(stateNew, "new state post BC");
} }
/** this.state = this._arrayClip2Zero(stateNew);
* Getter for effluent data. return stateNew;
* @returns {object} Effluent data object (msg), defaults to inlet 0. }
*/
get getEffluent() { /**
return { topic: "Fluent", payload: { inlet: 0, F: math.sum(this.Fs), C: this.state.at(-1) }, timestamp: this.currentTime }; * Apply boundary conditions to the reactor state.
* for inlet, apply generalised Danckwerts BC, if there is not flow, apply Neumann BC with no flux
* for outlet, apply regular Danckwerts BC (Neumann BC with no flux)
* @param {Array} state - Current reactor state without enforced BCs.
*/
_applyBoundaryConditions(state) {
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_dispersion_term = (1-this.alpha)*this.D*this.A/(math.sum(this.Fs)*this.d_x);
state[0] = math.multiply(1/(1+BC_dispersion_term), math.add(BC_C_in, math.multiply(BC_dispersion_term, state[1])));
} else {
state[0] = state[1];
} }
// Neumann BC (no flux)
state[this.n_x-1] = state[this.n_x-2];
}
updateState(newTime) { /**
super.updateState(newTime); * Create finite difference first derivative operator.
let Pe_local = this.d_x*math.sum(this.Fs)/(this.D*this.A) * @param {boolean} central - Use central difference scheme if true, otherwise use upwind scheme.
let Co_D = this.D*this.timeStep/(this.d_x*this.d_x); * @param {boolean} higher_order - Use higher order scheme if true, otherwise use first order scheme.
* @returns {Array} - First derivative operator matrix.
(Pe_local >= 2) && console.warn(`Local Péclet number (${Pe_local}) is too high! Increase reactor resolution.`); */
(Co_D >= 0.5) && console.warn(`Courant number (${Co_D}) is too high! Reduce time step size.`); _makeDoperator(central = false, higher_order = false) { // create gradient operator
if (higher_order) {
if(DEBUG) { if (central) {
console.log("Inlet state max " + math.max(this.state[0])) const I = math.resize(math.diag(Array(this.n_x).fill(1/12), -2), [this.n_x, this.n_x]);
console.log("Pe total " + this.length*math.sum(this.Fs)/(this.D*this.A)); const A = math.resize(math.diag(Array(this.n_x).fill(-2/3), -1), [this.n_x, this.n_x]);
console.log("Pe local " + Pe_local); const B = math.resize(math.diag(Array(this.n_x).fill(2/3), 1), [this.n_x, this.n_x]);
console.log("Co ad " + math.sum(this.Fs)*this.timeStep/(this.A*this.d_x)); const C = math.resize(math.diag(Array(this.n_x).fill(-1/12), 2), [this.n_x, this.n_x]);
console.log("Co D " + Co_D); 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;
D[1] = NearBoundary;
NearBoundary.reverse();
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;
} 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;
} }
}
/** /**
* Tick the reactor state using explicit finite difference method. * Create central finite difference second derivative operator.
* @param {number} time_step - Time step for the simulation [d]. * @returns {Array} - Second derivative operator matrix.
* @returns {Array} - New reactor state. */
*/ _makeD2operator() { // create the central second derivative operator
tick(time_step) { const I = math.diag(Array(this.n_x).fill(-2), 0);
const dispersion = math.multiply(this.D / (this.d_x*this.d_x), this.D2_op, this.state); const A = math.resize(math.diag(Array(this.n_x).fill(1), 1), [this.n_x, this.n_x]);
const advection = math.multiply(-1 * math.sum(this.Fs) / (this.A*this.d_x), this.D_op, this.state); const B = math.resize(math.diag(Array(this.n_x).fill(1), -1), [this.n_x, this.n_x]);
const reaction = this.state.map((state_slice) => this.asm.compute_dC(state_slice, this.temperature)); const D2 = math.add(I, A, B);
const transfer = Array.from(Array(this.n_x), () => new Array(NUM_SPECIES).fill(0)); D2[0] = Array(this.n_x).fill(0); // set by BCs elsewhere
D2[this.n_x - 1] = Array(this.n_x).fill(0);
if (isNaN(this.kla)) { // calculate OTR if kla is not NaN, otherwise use externally calculated OTR return D2;
transfer.forEach((x) => { x[S_O_INDEX] = this.OTR; }); }
} else {
transfer.forEach((x, i) => { x[S_O_INDEX] = this._calcOTR(this.state[i][S_O_INDEX], this.temperature); });
}
const dC_total = math.multiply(math.add(dispersion, advection, reaction, transfer), time_step);
const stateNew = math.add(this.state, dC_total);
this._applyBoundaryConditions(stateNew);
if (DEBUG) {
assertNoNaN(dispersion, "dispersion");
assertNoNaN(advection, "advection");
assertNoNaN(reaction, "reaction");
assertNoNaN(dC_total, "change in state");
assertNoNaN(stateNew, "new state post BC");
}
this.state = this._arrayClip2Zero(stateNew);
return stateNew;
}
/**
* Apply boundary conditions to the reactor state.
* for inlet, apply generalised Danckwerts BC, if there is not flow, apply Neumann BC with no flux
* for outlet, apply regular Danckwerts BC (Neumann BC with no flux)
* @param {Array} state - Current reactor state without enforced BCs.
*/
_applyBoundaryConditions(state) {
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_dispersion_term = (1-this.alpha)*this.D*this.A/(math.sum(this.Fs)*this.d_x);
state[0] = math.multiply(1/(1+BC_dispersion_term), math.add(BC_C_in, math.multiply(BC_dispersion_term, state[1])));
} else {
state[0] = state[1];
}
// Neumann BC (no flux)
state[this.n_x-1] = state[this.n_x-2]
}
/**
* Create finite difference first derivative operator.
* @param {boolean} central - Use central difference scheme if true, otherwise use upwind scheme.
* @param {boolean} higher_order - Use higher order scheme if true, otherwise use first order scheme.
* @returns {Array} - First derivative operator matrix.
*/
_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, 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;
D[1] = NearBoundary;
NearBoundary.reverse();
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;
} 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;
}
}
/**
* Create central finite difference second derivative operator.
* @returns {Array} - Second derivative operator matrix.
*/
_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]);
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;
}
} }
module.exports = { Reactor_CSTR, Reactor_PFR }; module.exports = { Reactor_CSTR, Reactor_PFR };