From 24de5a4c9f2129e1443a5f90ad4ca438c98ae79e Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Fri, 11 Jul 2025 12:22:36 +0200 Subject: [PATCH] Add generalFunctions dependency and implement basic measurement child registration in nodeClass --- package-lock.json | 6 + package.json | 1 + src/nodeClass.js | 6 + src/specificClass.js | 545 ++++++++++++++++++++++--------------------- 4 files changed, 294 insertions(+), 264 deletions(-) diff --git a/package-lock.json b/package-lock.json index e2e89fc..19b3188 100644 --- a/package-lock.json +++ b/package-lock.json @@ -9,6 +9,7 @@ "version": "0.0.1", "license": "SEE LICENSE", "dependencies": { + "generalFunctions": "git+https://gitea.centraal.wbd-rd.nl/p.vanderwilt/generalFunctions.git", "mathjs": "^14.5.2" } }, @@ -59,6 +60,11 @@ "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": { "version": "0.7.1", "resolved": "https://registry.npmjs.org/javascript-natural-sort/-/javascript-natural-sort-0.7.1.tgz", diff --git a/package.json b/package.json index e60fa97..b6d04c5 100644 --- a/package.json +++ b/package.json @@ -27,6 +27,7 @@ } }, "dependencies": { + "generalFunctions": "git+https://gitea.centraal.wbd-rd.nl/p.vanderwilt/generalFunctions.git", "mathjs": "^14.5.2" } } diff --git a/src/nodeClass.js b/src/nodeClass.js index 2e00ced..f4bbe25 100644 --- a/src/nodeClass.js +++ b/src/nodeClass.js @@ -48,6 +48,12 @@ class nodeClass { case "Dispersion": this.reactor.setDispersion = msg; 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: console.log("Unknown topic: " + msg.topic); } diff --git a/src/specificClass.js b/src/specificClass.js index 62236d1..509f460 100644 --- a/src/specificClass.js +++ b/src/specificClass.js @@ -1,9 +1,10 @@ const ASM3 = require('./reaction_modules/asm3_class.js'); const { create, all } = require('mathjs'); const { assertNoNaN } = require('./utils.js'); +const { childRegistrationUtils, logger, MeasurementContainer } = require('generalFunctions'); const config = { - matrix: 'Array' // use Array as the matrix type + matrix: 'Array' // use Array as the matrix type }; const math = create(all, config); @@ -13,304 +14,320 @@ const NUM_SPECIES = 13; const DEBUG = false; class Reactor { - /** - * Reactor base class. - * @param {object} config - Configuration object containing reactor parameters. - */ - constructor(config) { - this.asm = new ASM3(); + /** + * Reactor base class. + * @param {object} config - Configuration object containing reactor parameters. + */ + constructor(config) { + // 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.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.volume = config.volume; // fluid volume reactor [m3] - 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.timeStep = 1 / (24*60*15); // time step [d] - this.speedUpFactor = 60; // speed up factor for simulation, 60 means 1 minute per simulated second + this.kla = config.kla; // if NaN, use externaly provided OTR [d-1] + + 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. - * @param {object} input - Input object (msg) containing payload with inlet index, flow rate, and concentrations. - */ - set setInfluent(input) { - let index_in = input.payload.inlet; - this.Fs[index_in] = input.payload.F; - this.Cs_in[index_in] = input.payload.C; + /** + * Setter for influent data. + * @param {object} input - Input object (msg) containing payload with inlet index, flow rate, and concentrations. + */ + set setInfluent(input) { + let index_in = input.payload.inlet; + this.Fs[index_in] = input.payload.F; + 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). - * @param {object} input - Input object (msg) containing payload with OTR value [g O2 d-1]. - */ - set setOTR(input) { - this.OTR = input.payload; + /** + * 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; } - - /** - * Setter for temperature. - * @param {object} input - Input object (msg) containing payload with temperature value [C]. - */ - set setTemperature(input) { - this.temperature = input.payload; + this.currentTime += n_iter * this.timeStep * day2ms / this.speedUpFactor; } - - /** - * 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 { - /** - * Reactor_CSTR class for Continuous Stirred Tank Reactor. - * @param {object} config - Configuration object containing reactor parameters. - */ - constructor(config) { - super(config); - this.state = config.initialState; - } + /** + * Reactor_CSTR class for Continuous Stirred Tank Reactor. + * @param {object} config - Configuration object containing reactor parameters. + */ + constructor(config) { + super(config); + this.state = config.initialState; + } - /** - * Getter for effluent data. - * @returns {object} Effluent data object (msg), 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 }; - } + /** + * Getter for effluent data. + * @returns {object} Effluent data object (msg), 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 }; + } - /** - * Tick the reactor state using the forward Euler method. - * @param {number} time_step - Time step for the simulation [d]. - * @returns {Array} - New reactor state. - */ - 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 outflow = math.multiply(-1 * math.sum(this.Fs) / this.volume, this.state); - const reaction = this.asm.compute_dC(this.state, this.temperature); - 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 + /** + * Tick the reactor state using the forward Euler method. + * @param {number} time_step - Time step for the simulation [d]. + * @returns {Array} - New reactor state. + */ + 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 outflow = math.multiply(-1 * math.sum(this.Fs) / this.volume, this.state); + const reaction = this.asm.compute_dC(this.state, this.temperature); + 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 - 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 - if(DEBUG){ - assertNoNaN(dC_total, "change in state"); - assertNoNaN(this.state, "new state"); - } - return this.state; - } + 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 + if(DEBUG){ + assertNoNaN(dC_total, "change in state"); + assertNoNaN(this.state, "new state"); + } + return this.state; + } } class Reactor_PFR extends Reactor { - /** - * Reactor_PFR class for Plug Flow Reactor. - * @param {object} config - Configuration object containing reactor parameters. - */ - constructor(config) { - super(config); + /** + * Reactor_PFR class for Plug Flow Reactor. + * @param {object} config - Configuration object containing reactor parameters. + */ + constructor(config) { + super(config); - this.length = config.length; // reactor length [m] - this.n_x = config.resolution_L; // number of slices + this.length = config.length; // reactor length [m] + this.n_x = config.resolution_L; // number of slices - this.d_x = this.length / this.n_x; - this.A = this.volume / this.length; // crosssectional area [m2] + this.d_x = this.length / this.n_x; + 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(this.state) + // console.log("Initial 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); - assertNoNaN(this.D_op, "Derivative operator"); + this.D_op = this._makeDoperator(true, true); + assertNoNaN(this.D_op, "Derivative operator"); - this.D2_op = this._makeD2operator(); - assertNoNaN(this.D2_op, "Second derivative operator"); + this.D2_op = this._makeD2operator(); + 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); }); } - /** - * 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; + 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"); } - /** - * 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 }; + 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]; + } - 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); - } + /** + * 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; } + } - /** - * 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); - - 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; - } + /** + * 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 };