Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update: include relaxFactors for gas and condensed species wihtin the EquilibriumSolver class and simplify procedure, solve minor error initialization procedure, and add miscellaneous changes #976

Merged
merged 6 commits into from
Aug 28, 2024
Merged

Update: include relaxFactors for gas and condensed species wihtin the EquilibriumSolver class and simplify procedure, solve minor error initialization procedure, and add miscellaneous changes #976

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
10f2510
Update: miscellaneous
AlbertoCuadra Aug 27, 2024
22b6192
Update: miscellaneous
AlbertoCuadra Aug 27, 2024
f16e93f
Solve: minor error initialization procedure
AlbertoCuadra Aug 27, 2024
bbaf5ef
Update: rename variables, generalize tolerances and solve error non-s…
AlbertoCuadra Aug 28, 2024
57712ba
Update: remove comments
AlbertoCuadra Aug 28, 2024
1861783
Update: include relaxFactors for gas and condensed species wihtin the…
AlbertoCuadra Aug 28, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion +combustiontoolbox/+core/@ChemicalSystem/ChemicalSystem.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@
% See also: :mat:func:`Database`, :mat:func:`NasaDatabase`, :mat:func:`findProducts`, :mat:func:`setListSpecies`

properties
species % Struct with class Species
species % Struct with Species objects
listSpecies % List of species
listElements % List of elements
stoichiometricMatrix % Stoichiometric matrix
Expand Down
45 changes: 39 additions & 6 deletions +combustiontoolbox/+equilibrium/@EquilibriumSolver/EquilibriumSolver.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,7 @@
tolTau = 1e-25 % Tolerance of the slack variables for condensed species
itMaxGibbs = 70 % Max number of iterations - Gibbs/Helmholtz minimization method
itMaxIons = 30 % Max number of iterations - charge balance (ions)
slackGuess = 1e-14 % Initial guess of the slack variables for condensed species
temperatureIons = 0 % Minimum temperature [K] to consider ionized species
tol0 = 1e-3 % Tolerance of the root finding algorithm
itMax = 30 % Max number of iterations - root finding method
Expand Down Expand Up @@ -74,6 +75,7 @@
addParameter(p, 'itMaxGibbs', obj.itMaxGibbs, @(x) isnumeric(x) && x > 0);
addParameter(p, 'itMaxIons', obj.itMaxIons, @(x) isnumeric(x) && x > 0);
addParameter(p, 'temperatureIons', obj.temperatureIons, @(x) isnumeric(x) && x >= 0);
addParameter(p, 'slackGuess', obj.slackGuess, @(x) isnumeric(x) && x > 0);
addParameter(p, 'tol0', obj.tol0, @(x) isnumeric(x) && x > 0);
addParameter(p, 'itMax', obj.itMax, @(x) isnumeric(x) && x > 0);
addParameter(p, 'rootMethod', obj.rootMethod, @(method) ismember(method, {@newton, @steff, @nsteff}));
Expand Down Expand Up @@ -103,6 +105,7 @@
obj.itMaxGibbs = p.Results.itMaxGibbs;
obj.itMaxIons = p.Results.itMaxIons;
obj.temperatureIons = p.Results.temperatureIons;
obj.slackGuess = p.Results.slackGuess;
obj.tol0 = p.Results.tol0;
obj.itMax = p.Results.itMax;
obj.rootMethod = p.Results.rootMethod;
Expand Down Expand Up @@ -431,15 +434,45 @@ function printTime(obj)
NS = length(index);
end

function delta = relaxFactor(NP, ni, eta, Delta_ln_NP, NG)
% Compute relaxation factor
FLAG = eta(1:NG) > 0;
FLAG_MINOR = ni(1:NG) / NP <= 1e-8 & FLAG;
delta1 = 2./max(5*abs(Delta_ln_NP), abs(eta(FLAG)));
delta2 = min(abs((-log(ni(FLAG_MINOR)/NP) - 9.2103404) ./ (eta(FLAG_MINOR) - Delta_ln_NP)));
function delta = relaxFactorGas(NP, nj_gas, Delta_ln_nj, Delta_ln_NP)
% Compute relaxation factor for gas species
FLAG = Delta_ln_nj > 0;
FLAG_MINOR = nj_gas / NP <= 1e-8 & FLAG;
delta1 = 2 / max(abs( [5 * Delta_ln_NP; Delta_ln_nj(FLAG)] ));
delta2 = min(abs( (-log(nj_gas(FLAG_MINOR) / NP) - 9.2103404) ./ (Delta_ln_nj(FLAG_MINOR) - Delta_ln_NP) ));
delta = min([1; delta1; delta2]);
end

function [N, psi_j, FLAG_UNSTABLE] = relaxFactorCondensed(NP, N, psi_j, Delta_nj, indexCondensed, NG, NS, SIZE, tau, RT)
% Compute and apply relaxation factor for condensed species

% Definitions
delta0 = 0.9999;

% Initialization
FLAG_UNSTABLE = false;

% Check if there are condensed species
if NS - NG == 0
return
end

% Compute and apply relaxation factor
FLAG_DELTA = N(indexCondensed) + Delta_nj < 0;
deltaCondensed = min([1; -delta0 * N(indexCondensed(FLAG_DELTA)) ./ Delta_nj(FLAG_DELTA)]);
N(indexCondensed) = N(indexCondensed) + deltaCondensed .* Delta_nj;

Delta_psi_j = (tau - psi_j(indexCondensed) .* Delta_nj) ./ N(indexCondensed) - psi_j(indexCondensed);
FLAG_DELTA = psi_j(indexCondensed) + Delta_psi_j < 0;
deltaCondensed = min([1; -delta0 * psi_j(indexCondensed(FLAG_DELTA)) ./ Delta_psi_j(FLAG_DELTA)]);
psi_j(indexCondensed) = psi_j(indexCondensed) + deltaCondensed .* Delta_psi_j;

% Check if there are unstable species
Omega_pi = exp(-psi_j(indexCondensed) / RT);
FLAG_UNSTABLE = (N(indexCondensed) / NP < exp(-SIZE)) | (abs(log10(Omega_pi)) > 1e-2);
N(indexCondensed(FLAG_UNSTABLE)) = 0;
end

function point = getPoint(x_vector, f_vector)
% Get point using the regula falsi method
%
Expand Down
2 changes: 1 addition & 1 deletion +combustiontoolbox/+equilibrium/@EquilibriumSolver/equilibriumCheckCondensed.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@
end

% Calculate dLdnj of the condensed species
dL_dnj(i) = (muRT(indexCondensed(i)) - dot(pi_i, A0(indexCondensed(i), :))) / W(i); % / W(i);
dL_dnj(i) = (muRT(indexCondensed(i)) - dot(pi_i, A0(indexCondensed(i), :))) / W(i);
end

% Get condensed species that may appear at chemical equilibrium
Expand Down
96 changes: 38 additions & 58 deletions +combustiontoolbox/+equilibrium/@EquilibriumSolver/equilibriumGibbs.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,13 +44,11 @@
R0 = combustiontoolbox.common.Constants.R0; % Universal gas constant [J/(K mol)]

% Definitions
% CHECK: ERROR IN N(:, 2) FOR CONDENSED SPECIES.
N = system.propertyVector; % Composition matrix [moles_i, phase_i]
N = system.propertyVector; % Composition vector [moles_i]
A0 = system.stoichiometricMatrix; % Stoichiometric matrix [a_ij]
RT = R0 * T; % [J/mol]
delta0 = 0.9999;
tau0RT = obj.tolTau;
opts.SYM = true; % Options linsolve method: real symmetric
tau0 = obj.tolTau; % Tolerance of the slack variables for condensed species
opts.SYM = true; % Options linsolve method: real symmetric

% Initialization
NatomE = mix.natomElementsReact;
Expand Down Expand Up @@ -108,7 +106,7 @@

% Initialization
psi_j = system.propertyVector;
tauRT = tau0RT .* min(NatomE);
tau = tau0 .* min(NatomE);

% Solve system
x = equilibriumLoop;
Expand All @@ -132,31 +130,23 @@
% NESTED FUNCTION
function x = equilibriumLoop
% Calculate composition at chemical equilibrium

% persistent totalIterations
%
% if isempty(totalIterations)
% totalIterations = 0;
% end

% Initialization
it = 0; counter_errors = 0;
itMax = obj.itMaxGibbs;
STOP = 1.0;
FLAG_UNSTABLE = false;
delta_j0 = ones(NS - NG, 1);

% Calculations
while STOP > obj.tolGibbs && it < itMax
it = it + 1;
% Chemical potentials
muRT(indexGas) = g0(indexGas) / RT + log(N(indexGas, 1) / NP) + log(p);
muRT(indexGas) = g0(indexGas) / RT + log(N(indexGas) / NP) + log(p);

% Construction of matrix J
J = update_matrix_J(A0_T, J22, N, NP, indexGas, indexCondensed, NS - NG, psi_j);

% Construction of vector b
b = update_vector_b(A0, N, NP, NatomE, ind_E, index, indexGas, indexCondensed, indexIons, muRT, tauRT);
b = update_vector_b(A0, N, NP, NatomE, ind_E, index, indexGas, indexCondensed, indexIons, muRT, tau);

% Solve the linear system J*x = b
[x, ~] = linsolve(J, b, opts);
Expand All @@ -171,10 +161,14 @@

NG = length(indexGas);
NS = length(index);

if NS - NG > 0
J22 = zeros(NS - NG + 1);
end

% Reset removed species to 1e-6 to try the avoid singular matrix
N( N(index, 1) < obj.tolMoles, 1) = 1e-6;
psi_j(indexCondensed) = tauRT ./ N(indexCondensed, 1);
% Reset removed species to tolMolesGuess to try the avoid singular matrix
N( N(index) < obj.tolMoles ) = obj.tolMolesGuess;
psi_j(indexCondensed) = obj.slackGuess;

if counter_errors > 2
x = NaN;
Expand All @@ -193,50 +187,34 @@
% Compute correction moles of gases
Delta_ln_nj = update_Delta_ln_nj(A0, pi_i, Delta_ln_NP, muRT, indexGas);

% Calculate correction factor
delta = obj.relaxFactor(NP, N(index, 1), [Delta_ln_nj; Delta_nj], Delta_ln_NP, NG);
% Calculate correction factor for gases
deltaGas = obj.relaxFactorGas(NP, N(indexGas), Delta_ln_nj, Delta_ln_NP);

% Apply correction gaseous species and total moles in the mixture
N(indexGas, 1) = N(indexGas, 1) .* exp(delta * Delta_ln_nj);
NP = NP * exp(delta * Delta_ln_NP);

% Apply correction condensed species
if NS - NG > 0
delta_j = delta_j0;
FLAG_DELTA = N(indexCondensed, 1) + Delta_nj < 0;
delta_j(FLAG_DELTA) = -delta0 * N(indexCondensed(FLAG_DELTA), 1) ./ Delta_nj(FLAG_DELTA);
N(indexCondensed, 1) = N(indexCondensed, 1) + min(delta_j) .* Delta_nj;

delta_j = delta_j0;
Delta_psi_j = (tauRT - psi_j(indexCondensed) .* Delta_nj) ./ N(indexCondensed, 1) - psi_j(indexCondensed);
FLAG_DELTA = psi_j(indexCondensed) + Delta_psi_j < 0;
delta_j(FLAG_DELTA) = -delta0 * psi_j(indexCondensed(FLAG_DELTA)) ./ Delta_psi_j(FLAG_DELTA);
psi_j(indexCondensed) = psi_j(indexCondensed) + min(delta_j) .* Delta_psi_j;

Omega_pi = exp(-psi_j(indexCondensed));
FLAG_UNSTABLE = (N(indexCondensed, 1) / NP < exp(-SIZE)) | (abs(log10(Omega_pi)) > 1e-2);
N(indexCondensed(FLAG_UNSTABLE), 1) = 0;
end
N(indexGas) = N(indexGas) .* exp(deltaGas * Delta_ln_nj);
NP = NP * exp(deltaGas * Delta_ln_NP);

% Calculate and apply correction condensed species
[N, psi_j, FLAG_UNSTABLE] = obj.relaxFactorCondensed(NP, N, psi_j, Delta_nj, indexCondensed, NG, NS, SIZE, tau, RT);

% Compute STOP criteria
STOP = compute_STOP(NP, Delta_ln_NP, N(index, 1), [Delta_ln_nj; Delta_nj], NG, A0(index, :), NatomE, max_NatomE, obj.tolE);
STOP = compute_STOP(NP, Delta_ln_NP, N(index), [Delta_ln_nj; Delta_nj], NG, A0(index, :), NatomE, max_NatomE, obj.tolE);

% Update temp values in order to remove species with moles < tolerance
[index, indexCondensed, indexGas, indexIons, NG, NS, N] = obj.updateTemp(N, index, indexCondensed, indexGas, indexIons, NP, NG, NS, SIZE);

% Update psi_j vector
if sum(FLAG_UNSTABLE)
J22 = zeros(NS - NG + 1);
FLAG_UNSTABLE(:) = false;
delta_j0 = ones(NS - NG, 1);
end

% Debug
% aux_delta(it) = delta;
% aux_delta(it) = min(deltaGas, deltaCondensed);
% aux_STOP(it) = STOP;
end

% totalIterations = totalIterations + it;

% Debug
% debug_plot_error(it, aux_STOP, aux_delta);

% Check convergence of charge balance (ionized species)
[N, STOP_ions, FLAG_ION] = equilibriumCheckIons(obj, N, A0, ind_E, indexGas, indexIons);
Expand All @@ -263,9 +241,6 @@
% Remove element E from matrix
indexElements(ind_E) = [];
NE = NE - 1;

% Debug
% debug_plot_error(it, aux_STOP, aux_delta);
end

function x = equilibriumLoopCondensed(x)
Expand Down Expand Up @@ -340,16 +315,16 @@
N_backup = N;

% Check if there are non initialized condensed species
N(indexCondensed_add(N(indexCondensed_add, 1) == 0), 1) = 1e-5;
N(indexCondensed_add(N(indexCondensed_add) == 0)) = obj.tolMolesGuess;

% Initialize Lagrange multiplier vector psi
psi_j(indexCondensed_add) = 1e-15 ./ N(indexCondensed_add, 1);
psi_j(indexCondensed_add) = obj.slackGuess;

% Compute chemical equilibrium considering condensed species
x0 = equilibriumLoop;

% Debug
% aux2 = N(indexCondensed_add, 1);
% aux2 = N(indexCondensed_add);
% fprintf('\n n0 n\n');
% for k = 1:NC_add
% fprintf('%10s %1.3e %1.3e\n', system.listSpecies{indexCondensed_add(k)}, aux1(k), aux2(k));
Expand Down Expand Up @@ -399,6 +374,11 @@
function J11 = update_matrix_J11(A0_T, N, indexGas)
% Compute submatrix J11
J11 = A0_T(:, indexGas) * (A0_T(:, indexGas) .* N(indexGas)')';

% J11 is expected to be symmetric. However, due to precision errors,
% slight asymmetries may occur. To enforce symmetry, we explicitly
% symmetrize the matrix by averaging it with its transpose
J11 = (J11 + J11') / 2;
end

function J12 = update_matrix_J12(A0_T, N, indexGas, indexCondensed)
Expand All @@ -410,7 +390,7 @@

function J22 = update_matrix_J22(J22, N, NP, indexGas)
% Compute submatrix J22
J22(end, end) = sum(N(indexGas, 1)) - NP;
J22(end, end) = sum(N(indexGas)) - NP;
end

function J = update_matrix_J(A0_T, J22, N, NP, indexGas, indexCondensed, NC, psi_j)
Expand All @@ -422,16 +402,16 @@
J = [J11, J12; J12', J22];
end

function b = update_vector_b(A0, N, NP, NatomE, ind_E, index, indexGas, indexCondensed, indexIons, muRT, tauRT)
function b = update_vector_b(A0, N, NP, NatomE, ind_E, index, indexGas, indexCondensed, indexIons, muRT, tau)
% Compute vector b
bi = N(index)' * A0(index, :);

if any(indexIons)
bi(ind_E) = NatomE(ind_E);
end

b1 = (NatomE - bi + sum(A0(indexGas, :) .* N(indexGas, 1) .* muRT(indexGas)))';
b2 = muRT(indexCondensed) - tauRT ./ N(indexCondensed);
b1 = (NatomE - bi + sum(A0(indexGas, :) .* N(indexGas) .* muRT(indexGas)))';
b2 = muRT(indexCondensed) - tau ./ N(indexCondensed);
b3 = NP + sum(N(indexGas) .* muRT(indexGas) - N(indexGas));

b = [b1; b2; b3];
Expand Down
2 changes: 1 addition & 1 deletion +combustiontoolbox/+equilibrium/@EquilibriumSolver/equilibriumGuess.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,7 @@
N = getSimplex(N, A0, muRT, b0, index, indexIons, NG);
catch
% Get molar composition using a uniform distribution
N(index) = NP/NG;
N(indexGas) = NP/NG;
end

% Recompute mol gaseous species
Expand Down
Loading
Loading