const Reaction = require('./ReactionKinetics.js');
const Species = require('./Species.js');
const clamp = require('../util/clamp.js');
const SpeciesNode = require('./SpeciesNode.js');
const Constants = require('./Constants.js');
const RungeKutta4 = require('./RungeKutta4.js');
const logger = require('../util/logger.js');
function computeRelChangeInConc(speciesNodes, oldConcentrations) {
// Compute the relative change in concentration, and return it.
// Candidate for simplification
var relChangeInConc = 0.0;
// const change = speciesNodes.map( (x, k) => {
// newConc = current.getConcentration();
// if (x.infinite || (newConc <= Constants.MINIMUM_CONCENTRATION)) {
// return 0;
// } else {
// return Math.abs(oldConcentrations[k] - newConc) / newConc;
// }
// }).reduce((x,y)=>x+y, 0)
for (k = 0; k < speciesNodes.length; k++) {
current = speciesNodes[k];
if (!current.infinite) {
if (current.getConcentration() > Constants.MINIMUM_CONCENTRATION) {
relChange = Math.abs(oldConcentrations[k] - current.getConcentration()) /
current.getConcentration();
relChangeInConc += relChange;
}
}
}
return relChangeInConc;
}
/**
* This class models an individual reaction, including references to the
* quantities of products and reactants, as well as the model for the general
* reaction.
*
* Variables
* archetype {Reaction} Reference to our defining Reaction.
*
* speciesNodes {Array} List of SpeciesNodes that correspond one-to-one
* with the species in our archetype.
*
* finiteReactants {Boolean} Set to true if the Reactants of our reaction have finite limits.
*
* finiteProducts {Boolean} Set to true if the Products of our reaction have finite limits.
*
* leChatlierDenominator {Number} If one of the sides of our reaction is infinite, the maximum distance the
* reaction can travel will be defined by LeChatlier's principle.
*
* While that number is not constant, a large portion of it is, and can be
* cached here.
*
* leChatlierRootPower {Number} Similarly, the total sum of the coefficients is cached for computing
* leChatelier's constant, and is used as the root power.
*
* unitActivityReactants {Boolean} Set to true if there are reactants that have unit activity.
*
* unitActivityProducts {Boolean} Set to true if there are products that have unit activity.
*
* xLow {Number} The maximum amount that the reaction can progress in the foward
* direction.
*
* xHigh {Number} The maximum amount that the reaction can progress in the backward
* direction.
*
* xPreviousIteration {Number} Amount reaction was pushed during most recent call to iteration()
*
* relaxation {Number} Used to speed up the iteration's use of setRangeOfX. The relaxation is
* used to compute a range based on the previous iteration's push.
*
* xConcentrations {Array} Cached concentrations of our species, scaled by their coefficients.
*
* maxChangeAllowed {Number} Control value to avoid jumps in temperature caused by reaction advance
*
* isShiftLimited {Boolean} Control whether the shifting has been limited by maximum allowed changed
*
* MAX_TEMP_CHANGE {Number} Specifies the maximal temperature change that can be caused by one iteration of
* this reaction node. That avoids formation of bucles in sharp changes of K
* with temperature.
*
* solventTemperature {Number} Cached variables to increase equilibrium calculations
*
* currentK {Number} Cached variables to increase equilibrium calculations
*
* coef {Array} Cached variables to increase equilibrium calculations
*
* polarity {Array} Cached variables to increase equilibrium calculations
*
* isInfinite {Array} Cached variables to increase equilibrium calculations
*
* hasUnitActivity {Array} Cached variables to increase equilibrium calculations
*/
/**
* Creates a new ReactionNode from the template and SpeciesNodes. It is
* required that the SpeciesNodes within the Vector of nodes provided
* correspond exactly to the Vector of species within the Reaction.
*
* @param {Reaction} reactionTemplate
* @param {Array} speciesNodes
*
* @constructor
*/
class ReactionNodeKinetics {
constructor(reactionTemplate, speciesNodes) {
if (reactionTemplate === null || speciesNodes === null) {
throw "NullPointerException : arguments is null in ReactionNode constructor";
}
if (reactionTemplate.getSpeciesCount() !== speciesNodes.length) {
throw "IllegalArgumentException : length of arguments does not match in ReactionNode constructor";
}
/** @type {Reaction} */
this.archetype = reactionTemplate;
this.speciesNodes = speciesNodes;
this.xMoles = [];
this.coef = [];
this.polarity = [];
this.isInfinite = [];
this.hasUnitActivity = [];
this.finiteReactants = false;
this.finiteProducts = false;
this.leChatlierDenominator = 1.0;
this.leChatlierRootPower = 0.0;
this.unitActivityReactants = false;
this.unitActivityProducts = false;
this.xLow = 1e-100;
this.xHigh = 1e100;
this.xPreviousIteration = 1e-2;
this.relaxation = 10;
this.maxChangeAllowed = 1;
this.isShiftLimited = false;
this.MAX_TEMP_CHANGE = 2;
this.solventTemperature = this.getSolvent().getTemperature();
this.kinetic_reaction = false;
this.currentK = this.archetype.getK(this.solventTemperature);
var i = void 0;
var current = void 0;
var xMol = void 0;
i = 0;
while (i < this.speciesNodes.length) {
current = this.speciesNodes[i];
this.coef[i] = Math.abs(this.archetype.getCoefficientAt(i));
this.polarity[i] = this.coef[i] / this.archetype.getCoefficientAt(i);
this.isInfinite[i] = current.infinite;
this.hasUnitActivity[i] = current.unitActivity;
if (i !== this.archetype.getSpeciesIndex(current.archetype)) {
throw "Exception : Wrong order in ReactionNode.constructor";
}
if (!this.isInfinite[i]) {
if (this.polarity[i] < 0) {
this.finiteReactants = true;
} else {
this.finiteProducts = true;
}
this.leChatlierRootPower += this.coef[i];
this.leChatlierDenominator *= Math.pow(this.coef[i], this.coef[i]);
}
if (this.hasUnitActivity[i] && !this.isInfinite[i]) {
if (this.polarity[i] < 0) {
this.unitActivityReactants = true;
} else {
this.unitActivityProducts = true;
}
}
if (this.isInfinite[i]) {
xMol = Infinity/this.coef[i];
} else {
xMol = current.moles / this.coef[i];
}
this.setXMolesAt(i, xMol);
i++;
}
this.maxChangeAllowed = Math.max(Math.abs(this.MAX_TEMP_CHANGE / -this.archetype.getHo() / this.getSolvent().liquidVolume / 1000 * (this.getSolvent().getSolutionWeight() * this.getSolvent().getSpecificHeat()) * Constants.CAL_TO_J), 0.01);
if (reactionTemplate.kinetics) {
this.kinetic_reaction = true;
const concentrations = this.speciesNodes.map(x => x.getConcentration());
let f_kinetic;
if (reactionTemplate.useEquilibrium) {
f_kinetic = (t, y) => {
const k = this.archetype.getRateConstant(this.solventTemperature);
const K = this.archetype.getK(this.solventTemperature);
const k_rev = k/K;
let rate_fwd = k;
this.archetype.orders.forEach( (order, index) =>
{
if (order !== 0) {
rate_fwd *= Math.pow(y[index], order);
}
});
let rate_rev = k_rev;
this.archetype.coefficients.forEach( (coeff, index) => {
const order = coeff + this.archetype.orders[index];
if (order !== 0) {
rate_rev *= Math.pow(y[index], order);
}
});
const rate = rate_fwd - rate_rev;
return this.archetype.coefficients.map(x=> rate*x);
}
} else {
f_kinetic = (t, y) => {
const k = this.archetype.getRateConstant(this.solventTemperature);
let rate = k;
this.archetype.orders.forEach( (order, index) =>
{
if (order !== 0) {
rate *= Math.pow(y[index], order);
}
});
return this.archetype.coefficients.map(x=> rate*x);
}
}
this.rk = new RungeKutta4(f_kinetic, 0.0, concentrations, 0.1);
}
}
getMolesShifted(dX) {
const liquidVolume = this.getSolvent().liquidVolume;
return Object.fromEntries(this.speciesNodes.map((x, i) => [x.archetype.id, x.moles + this.coef[i] * this.polarity[i] * dX ]));
}
/**
* Evaluates the equilibrium when the reaction advances by the argument
* double and returns a measure of the difference between it's current value
* and its equilibrium value.
*
* @param dX {Number}
* @returns {Number}
*/
f(dX) {
//
const logK = Math.log(this.currentK);
let logQ = 0;
let xNewConc,
raised,
i;
let newMolesObj;
const liquidUnitActivity = this.getSolvent().equilibriumSettings.liquidUnitActivity
if (!liquidUnitActivity) {
newMolesObj = this.getMolesShifted(dX);
}
for (i = 0; i < this.speciesNodes.length; i++) {
const specie = this.speciesNodes[i];
xNewConc = 0;
if (!this.isInfinite[i]) {
xNewConc = specie.getConcentration() + this.polarity[i] * dX * this.coef[i] * specie.dConc_dmoles();
if (xNewConc < 1e-40) { // Or some other very small positive constant like Constants.MINIMUM_CONCENTRATION
// console.warn("Clamping concentration for " + specie.archetype.simpleName + " from " + xNewConc + " to 1e-30 in f(dX=" + dX + ")");
xNewConc = 1e-40;
}
}
if (!this.hasUnitActivity[i] && specie.archetype.state !== 'l') {
raised = this.coef[i]*Math.log(xNewConc);
if (this.polarity[i] < 0) {
logQ -= raised;
} else {
logQ += raised;
}
} else if (!liquidUnitActivity && specie.archetype.state == 'l') {
raised = Math.log(specie.getMoleFraction(newMolesObj))* this.coef[i];
if (this.polarity[i] < 0) {
logQ -= raised;
} else {
logQ += raised;
}
}
}
return logK - logQ;
};
/**
* Advances the reaction by the argument double.
*
* @param dX
*/
shift(dX) {
// Local Variables
var current,
xConc,
i;
for (i = 0; i < this.speciesNodes.length; i++) {
current = this.speciesNodes[i];
if (!this.isInfinite[i]) {
current.moles += this.polarity[i] * dX * this.coef[i];
}
}
};
/**
* Computes the maximum amount forwards and backwards that the reaction can
* progress, and stores these values in the variables xLow & xHigh. When
* both sides are finite, these values are computed via limiting reagents.
* However, if either side is infinite, then we are bound by LeChatlier's
* principle.
*
* This method will update the cached XConcentrations.
*/
setRangeOfX() {
this.xHigh = Number.POSITIVE_INFINITY;
this.xLow = Number.NEGATIVE_INFINITY;
// Local Variables
var current,
i,
k,
kcd,
root,
leChatlier;
for (i = 0; i < this.speciesNodes.length; i++) {
current = this.speciesNodes[i];
if (!this.isInfinite[i]) {
// Update cached Concentrations.
this.setXMolesAt(i, current.moles / this.coef[i]);
// The range is computed by limiting reagents.
if (this.polarity[i] < 0) {
this.xHigh = Math.min(this.xMoles[i], this.xHigh);
}
else {
this.xLow = Math.max(-this.xMoles[i], this.xLow);
}
}
}
// If one of the sides is infinite, then it's bound will still be set
// to either Double.MAX_VALUE or Double.MIN_VALUE. However, we can
// set a bound on this value by using Le Chatlier's principle. For
// efficiency, this value was partially computed at initialization if
// either of these cases is applicable.
if (!this.finiteReactants || !this.finiteProducts) {
k = this.currentK;
kcd = k / this.leChatlierDenominator;
root = 1 / this.leChatlierRootPower;
leChatlier = Math.pow(kcd, root);
if (this.finiteReactants) {
this.xLow = -1 * leChatlier;
}
else {
this.xHigh = leChatlier;
}
}
};
/**
* Brings the solution to equilibrium.
*
* @param {boolean} constantTemperature Keep temperature constant during the iteration
* @param {number|null} [dT=null] The difference in time (in seconds) since previous equilibrium-needed for reactions using kinetics
* @returns {number}
*/
iteration(constantTemperature, dT=null) {
this.solventTemperature = this.getSolvent().getTemperature();
this.currentK = this.archetype.getK(this.solventTemperature);
this.isShiftLimited = false;
this.depletionState = 0;
// Local Variables
var oldConcentrations,
current,
amtLastPushed,
xLowTry,
xHighTry,
fLow,
fHigh,
zeroConcSolid,
i,
k,
change,
// Store the old concentrations.
oldConcentrations = [];
for (i = 0; i < this.speciesNodes.length; i++) {
current = this.speciesNodes[i];
if (!current.infinite) {
oldConcentrations[i] = current.getConcentration();
}
}
// NEW: Check if the reaction is already at equilibrium
const f0 = this.f(0);
// Define a constant for equilibrium tolerance
const EQUILIBRIUM_TOLERANCE_F = 1e-10
const STOICHIOMETRIC_LIMIT_TOLERANCE_DX = 1e-23; // For comparing xHigh/xLow to 0 effectively
if (Math.abs(f0) < EQUILIBRIUM_TOLERANCE_F) {
// The reaction is already at equilibrium, no need to shift
this.xPreviousIteration = 0;
return 0; // Return zero change
}
// Here try to introduce kinetics code...
if (dT !== null) {
this.setRangeOfX();
this.rk.y = oldConcentrations.slice(0); // Copy
const newC = this.rk.deltaT(dT); // Copy
const volume = this.getSolvent().liquidVolume;
let extent = newC.map((x, i) => ((x - oldConcentrations[i]) * volume) / this.archetype.coefficients[i])
.reduce((x, y) => x + y, 0) / newC.length;
if (extent > this.xHigh) {
extent = this.xHigh;
} else if (extent < this.xLow) {
extent = this.xLow;
}
this.shift(extent);
this.xPreviousIteration = extent;
return 0;
}
// Attempt to bracket based off of stoichiometry.
this.setRangeOfX();
let doBisection = false;
// Sanity check for invalid stoichiometric range
if (this.xLow > this.xHigh + STOICHIOMETRIC_LIMIT_TOLERANCE_DX) {
this.xPreviousIteration = 0;
return 0;
} else {
let bisectionLowX, bisectionHighX, fAtBracketLow, fAtBracketHigh; // For bisection setup
// Reminder: f(dX) > 0 means K > Q (FORWARD shift for dX > 0)
// f(dX) < 0 means K < Q (REVERSE shift for dX < 0)
if (f0 > 0) { // Reaction wants to shift FORWARD (K > Q_initial)
const fAtXHigh = this.f(this.xHigh);
// If K >= Q even at the max forward extent (xHigh), then reaction goes to xHigh.
// This covers the case where xHigh is ~0 (then fAtXHigh ~ f0, which is >0).
if (fAtXHigh >= -EQUILIBRIUM_TOLERANCE_F) {
this.xPreviousIteration = this.xHigh;
this.depletionState = 1; // Out of reactants
} else { // K < Q at xHigh, so equilibrium was crossed. Root is between 0 and xHigh.
bisectionLowX = 0;
fAtBracketLow = f0;
bisectionHighX = this.xHigh;
fAtBracketHigh = fAtXHigh;
doBisection = true;
}
} else { // Reaction wants to shift REVERSE (f0 < 0, K < Q_initial)
const fAtXLow = this.f(this.xLow);
// If K <= Q even at the max reverse extent (xLow), then reaction goes to xLow.
// This covers the case where xLow is ~0 (then fAtXLow ~ f0, which is <0).
if (fAtXLow <= EQUILIBRIUM_TOLERANCE_F) {
this.xPreviousIteration = this.xLow;
this.depletionState = -1; // Out of products
} else { // K > Q at xLow, so equilibrium was crossed. Root is between xLow and 0.
bisectionLowX = this.xLow;
fAtBracketLow = fAtXLow;
bisectionHighX = 0;
fAtBracketHigh = f0;
doBisection = true;
}
}
if (doBisection) {
// --- Bisection Method ---
let currentBisectionXLow = bisectionLowX;
let currentBisectionXHigh = bisectionHighX;
let fAtCurrentXLow = fAtBracketLow; // This is f(currentBisectionXLow)
var nIterations = (currentBisectionXHigh - currentBisectionXLow <= 0) ? 10 :
Math.ceil(Math.log((currentBisectionXHigh - currentBisectionXLow) / Constants.REACTION_BRACKET_ACCURACY) / Math.log(2));
nIterations = Math.max(nIterations, 10);
var xmid = (currentBisectionXLow + currentBisectionXHigh) / 2.0;
var fmid = NaN; // Ensure fmid is defined
// Check if bracket ends are already the root (unlikely if f0 !=0 and previous checks passed, but for safety)
if (Math.abs(fAtCurrentXLow) < EQUILIBRIUM_TOLERANCE_F) {
this.xPreviousIteration = currentBisectionXLow;
} else if (Math.abs(fAtBracketHigh) < EQUILIBRIUM_TOLERANCE_F) { // fAtBracketHigh is f(original bisectionHighX)
this.xPreviousIteration = currentBisectionXHigh;
} else {
for (let j = 0; j < nIterations; j++) {
if ((currentBisectionXHigh - currentBisectionXLow) <= Constants.REACTION_BRACKET_ACCURACY) {
break;
}
xmid = (currentBisectionXHigh + currentBisectionXLow) / 2.0;
if (xmid === currentBisectionXLow || xmid === currentBisectionXHigh) break; // Precision limit
fmid = this.f(xmid);
if (Math.abs(fmid) < EQUILIBRIUM_TOLERANCE_F) {
currentBisectionXLow = xmid; // Converged
currentBisectionXHigh = xmid;
break;
}
if (fAtCurrentXLow * fmid < 0.0) { // Root is between currentBisectionXLow and xmid
currentBisectionXHigh = xmid;
} else { // Root is between xmid and currentBisectionXHigh
currentBisectionXLow = xmid;
fAtCurrentXLow = fmid;
}
}
this.xPreviousIteration = (currentBisectionXHigh + currentBisectionXLow) / 2.0;
}
}
} // End of main equilibrium logic (after xLow > xHigh check)
// Limit the expected change in temperature to allow for
// re-equilibration, if temperature is not constant
if (!constantTemperature
&& Math.abs(this.xPreviousIteration) > Math.abs(this.maxChangeAllowed)) {
this.isShiftLimited = true;
this.xPreviousIteration = (this.xPreviousIteration < 0) ? -this.maxChangeAllowed
: this.maxChangeAllowed;
}
this.shift(this.xPreviousIteration);
return computeRelChangeInConc(this.speciesNodes, oldConcentrations);
};
/**
* Converts the ReactionNode into a String.
*/
toString() {
var returnValue = "Reaction:",
i;
for (i = 0; i < this.speciesNodes.length; i++) {
returnValue += this.speciesNodes[i].toString() + "\n\t" + this.archetype.getCoefficientAt(i) + "\t";
}
return returnValue;
};
/**
* Returns a reference to the solution that this reactionNode is in.
*
* @returns {Solution}
*/
getSolvent() {
return this.speciesNodes[0].solvent;
};
/**
* Returns the species node at the specified index. This species node will
* correspond to the Reaction archetype's species at the same index.
*
* @param index
* @returns {SpeciesNode}
*/
getSpeciesNodeAt(index) {
return this.speciesNodes[index];
}
/**
* Sets the concentration (relative to the specie's coefficient) for the
* species at the specified index.
*
* @param index
* @param value
*/
setXMolesAt(index, value) {
this.xMoles[index] = value;
// FIX - RD - This is probably not needed.
this.maxChangeAllowed = Math.abs(this.MAX_TEMP_CHANGE
/ -this.archetype.getHo()
/ this.getSolvent().liquidVolume
/ 1000
* (this.getSolvent().getSolutionWeight() * this.getSolvent()
.getSpecificHeat()) * Constants.CAL_TO_J);
}
/**
* Determines whether a ReactionNode is equal by comparing the template.
*
* @param other
*
* @returns {Boolean}
*/
equals(other) {
if (other instanceof ReactionNodeKinetics) {
return this.archetype.equals(other.archetype);
}
else if (other instanceof Reaction) {
return this.archetype.equals(other);
} else {
return false;
}
}
}
module.exports = ReactionNodeKinetics;