Skip to content
GitLab
Commit eaaf9ad7 authored by Andreas Michael Bremen's avatar Andreas Michael Bremen
Browse files

Initial commit

parent 0420642a
Hide whitespace changes
Inline Side-by-side
ElectrolyteMedia/Media/Common/OneNonLinearEquation/package.mo 0 → 100644
View file @ eaaf9ad7
within ElectrolyteMedia.Media.Common;
package OneNonLinearEquation "Determine solution of a non-linear algebraic equation in one unknown without derivatives in a reliable and efficient way"
extends Modelica.Icons.Package;
replaceable partial function f_nonlinear
"Nonlinear algebraic equation in one unknown: y = f_nonlinear(x,p,X)"
extends Modelica.Icons.Function;
input Real x "Independent variable of function";
input Real p_=0.0 "Disregarded variables (here always used for pressure)";
input Real T_= 0 "Disregarded variables (here always used for temperature)";
input Real d_= 0 "Disregarded variables (here always used for density)";
input Real[:] X_=fill(0, 0)
"Disregarded variables (her always used for composition)";
output Real y "= f_nonlinear(x)";
// annotation(derivative(zeroDerivative=y)); // this must hold for all replaced functions
end f_nonlinear;
replaceable function solve
"Solve f_nonlinear(x_zero)=y_zero; f_nonlinear(x_min) - y_zero and f_nonlinear(x_max)-y_zero must have different sign"
import Modelica.Utilities.Streams.error;
extends Modelica.Icons.Function;
input Real y_zero
"Determine x_zero, such that f_nonlinear(x_zero) = y_zero";
input Real x_min "Minimum value of x";
input Real x_max "Maximum value of x";
input Real pressure=0.0
"Disregarded variables (here always used for pressure)";
input Real temperature= 0 "Disregarded variables (here always used for temperature)";
input Real density= 0 "Disregarded variables (here always used for density)";
input Real[:] X=fill(0, 0)
"Disregarded variables (here always used for composition)";
input Real x_tol=100*Modelica.Constants.eps
"Relative tolerance of the result";
output Real x_zero "f_nonlinear(x_zero) = y_zero";
protected
constant Real eps=Modelica.Constants.eps "Machine epsilon";
constant Real x_eps=1e-10
"Slight modification of x_min, x_max, since x_min, x_max are usually exactly at the borders T_min/h_min and then small numeric noise may make the interval invalid";
Real x_min2=x_min - x_eps;
Real x_max2=x_max + x_eps;
Real a=x_min2 "Current best minimum interval value";
Real b=x_max2 "Current best maximum interval value";
Real c "Intermediate point a <= c <= b";
Real d;
Real e "b - a";
Real m;
Real s;
Real p;
Real q;
Real r;
Real tol;
Real fa "= f_nonlinear(a) - y_zero";
Real fb "= f_nonlinear(b) - y_zero";
Real fc;
Boolean found=false;
algorithm
// Check that f(x_min) and f(x_max) have different sign
fa := f_nonlinear(
x_min2,
pressure,
temperature,
density,
X) - y_zero;
fb := f_nonlinear(
x_max2,
pressure,
temperature,
density,
X) - y_zero;
fc := fb;
if fa > 0.0 and fb > 0.0 or fa < 0.0 and fb < 0.0 then
error(
"The arguments x_min and x_max to OneNonLinearEquation.solve(..)\n"
+ "do not bracket the root of the single non-linear equation:\n" +
" x_min = " + String(x_min2) + "\n" + " x_max = " + String(x_max2)
+ "\n" + " y_zero = " + String(y_zero) + "\n" +
" fa = f(x_min) - y_zero = " + String(fa) + "\n" +
" fb = f(x_max) - y_zero = " + String(fb) + "\n" +
"fa and fb must have opposite sign which is not the case");
end if;
// Initialize variables
c := a;
fc := fa;
e := b - a;
d := e;
// Search loop
while not found loop
if abs(fc) < abs(fb) then
a := b;
b := c;
c := a;
fa := fb;
fb := fc;
fc := fa;
end if;
tol := 2*eps*abs(b) + x_tol;
m := (c - b)/2;
if abs(m) <= tol or fb == 0.0 then
// root found (interval is small enough)
found := true;
x_zero := b;
else
// Determine if a bisection is needed
if abs(e) < tol or abs(fa) <= abs(fb) then
e := m;
d := e;
else
s := fb/fa;
if a == c then
// linear interpolation
p := 2*m*s;
q := 1 - s;
else
// inverse quadratic interpolation
q := fa/fc;
r := fb/fc;
p := s*(2*m*q*(q - r) - (b - a)*(r - 1));
q := (q - 1)*(r - 1)*(s - 1);
end if;
if p > 0 then
q := -q;
else
p := -p;
end if;
s := e;
e := d;
if 2*p < 3*m*q - abs(tol*q) and p < abs(0.5*s*q) then
// interpolation successful
d := p/q;
else
// use bi-section
e := m;
d := e;
end if;
end if;
// Best guess value is defined as "a"
a := b;
fa := fb;
b := b + (if abs(d) > tol then d else if m > 0 then tol else -tol);
fb := f_nonlinear(
b,
pressure,
temperature,
density,
X) - y_zero;
if fb > 0 and fc > 0 or fb < 0 and fc < 0 then
// initialize variables
c := a;
fc := fa;
e := b - a;
d := e;
end if;
end if;
end while;
end solve;
annotation (Documentation(info="<html>
<p>
This function should currently only be used in Modelica.Media,
since it might be replaced in the future by another strategy,
where the tool is responsible for the solution of the non-linear
equation.
</p>
<p>
This library determines the solution of one non-linear algebraic equation \"y=f(x)\"
in one unknown \"x\" in a reliable way. As input, the desired value y of the
non-linear function has to be given, as well as an interval x_min, x_max that
contains the solution, i.e., \"f(x_min) - y\" and \"f(x_max) - y\" must
have a different sign. If possible, a smaller interval is computed by
inverse quadratic interpolation (interpolating with a quadratic polynomial
through the last 3 points and computing the zero). If this fails,
bisection is used, which always reduces the interval by a factor of 2.
The inverse quadratic interpolation method has superlinear convergence.
This is roughly the same convergence rate as a globally convergent Newton
method, but without the need to compute derivatives of the non-linear
function. The solver function is a direct mapping of the Algol 60 procedure
\"zero\" to Modelica, from:
</p>
<dl>
<dt> Brent R.P.:</dt>
<dd> <b>Algorithms for Minimization without derivatives</b>.
Prentice Hall, 1973, pp. 58-59.</dd>
</dl>
<p>
Due to current limitations of the
Modelica language (not possible to pass a function reference to a function),
the construction to use this solver on a user-defined function is a bit
complicated (this method is from Hans Olsson, Dassault Syst&egrave;mes AB). A user has to
provide a package in the following way:
</p>
<pre>
<b>package</b> MyNonLinearSolver
<b>extends</b> OneNonLinearEquation;
<b>redeclare record extends</b> Data
// Define data to be passed to user function
...
<b>end</b> Data;
<b>redeclare function extends</b> f_nonlinear
<b>algorithm</b>
// Compute the non-linear equation: y = f(x, Data)
<b>end</b> f_nonlinear;
// Dummy definition that has to be present for current Dymola
<b>redeclare function extends</b> solve
<b>end</b> solve;
<b>end</b> MyNonLinearSolver;
x_zero = MyNonLinearSolver.solve(y_zero, x_min, x_max, data=data);
</pre>
</html>"));
end OneNonLinearEquation;
ElectrolyteMedia/Media/Common/Reaction/calc_L.mo 0 → 100644
View file @ eaaf9ad7
within ElectrolyteMedia.Media.Common.Reaction;
function calc_L "calculates lower L of LU decomposition"
input Real[:,:] LU;
output Real[size(LU,1),size(LU,2)] L;
protected
Integer nF = size(LU,1);
Integer nR = size(LU,2);
algorithm
//compute lower L
for i in 1:nF loop
for j in 1:nR loop
if i == j then
L[i,j] :=1;
elseif j < i then
L[i,j] :=LU[i, j];
end if;
end for;
end for;
end calc_L;
ElectrolyteMedia/Media/Common/Reaction/calc_P.mo 0 → 100644
View file @ eaaf9ad7
within ElectrolyteMedia.Media.Common.Reaction;
function calc_P
"calculates permuation matrix of pivot vector from LAPACK"
input Integer[:] pivots_LAPACK_;
input Integer nF;
output Integer [nF,nF] P;
protected
Integer nR = size(pivots_LAPACK_,1);
Integer[nF] pivots;
Integer[nF] pivots_LAPACK;
Integer temp;
algorithm
// convert pivots_LAPACK_ to pivot vector
pivots :={i for i in 1:nF};
pivots_LAPACK[1:nR] :=pivots_LAPACK_;
pivots_LAPACK[nR+1:nF] :={i for i in nR + 1:nF};
for i in 1:nF-1 loop
temp :=pivots[pivots_LAPACK[i]];
pivots[pivots_LAPACK[i]] :=pivots[i];
pivots[i] :=temp;
end for;
// compute permutation matrix P
for i in 1:nF loop
P[pivots[i],i] :=1;
end for;
end calc_P;
ElectrolyteMedia/Media/Common/Reaction/calc_P_lambda_id.mo 0 → 100644
View file @ eaaf9ad7
within ElectrolyteMedia.Media.Common.Reaction;
function calc_P_lambda_id
"Calculates a stoichiometry matrix with pivoting so that lambdaout =(I|lambda_var) with permutation matrix Pout for reordering species elements of lambda_ordered = lambdaout*Pout"
input Real[:,:] lambda;
output Integer[size(lambda,1),size(lambda,1)] Pout;
protected
Real[size(lambda,1),size(lambda,2)] lambdaout;
algorithm
(lambdaout,Pout) :=calc_lambda_id(lambda);
end calc_P_lambda_id;
ElectrolyteMedia/Media/Common/Reaction/calc_P_nu_id.mo 0 → 100644
View file @ eaaf9ad7
within ElectrolyteMedia.Media.Common.Reaction;
function calc_P_nu_id
"Calculates a stoichiometry matrix with pivoting so that nuout =(I|nu_var) with permutation matrix Pout for reordering species elements of nu_ordered = nuout*Pout"
input Real[:,:] nu;
output Integer[size(nu,2),size(nu,2)] Pout;
protected
Real[size(nu,1),size(nu,2)] nuout;
algorithm
(nuout,Pout) :=calc_nu_id(nu);
end calc_P_nu_id;
ElectrolyteMedia/Media/Common/Reaction/calc_U.mo 0 → 100644
View file @ eaaf9ad7
within ElectrolyteMedia.Media.Common.Reaction;
function calc_U "calculates upper U of LU decomposition"
input Real[:,:] LU;
output Real[size(LU,2),size(LU,2)] U;
protected
Integer nF = size(LU,1);
Integer nR = size(LU,2);
algorithm
//compute upper U
for i in 1:nR loop
for j in 1:nR loop
if i < j+1 then
U[i,j] :=LU[i, j];
end if;
end for;
end for;
end calc_U;
ElectrolyteMedia/Media/Common/Reaction/calc_lambda.mo 0 → 100644
View file @ eaaf9ad7
within ElectrolyteMedia.Media.Common.Reaction;
function calc_lambda
"Calculates null vectors of stoichiometry matrix nu with positive elements only"
input Real[:,:] nu;
// input Real[nR,nF] nu;
input Integer nR;
input Integer nF;
output Real[nF,nF-nR] lambda;
protected
Integer[nF] sort;
Real[nF,nR+nF] total;
Real[1,nR+nF] total_ = -1*ones(1,nF+nR);
Real[1,nR+nF] total__;
Real[nF] iDummy;
Integer temp;
Real[nF-nR,nF] lambdaT;
Real[nF] iDummy1;
Integer[nF] sort1;
Integer[nF,nF] resort1;
Boolean[nF] resortindex = fill(false,nF);
Integer index;
algorithm
if nR > 0 then
// total = (transpose(nu)|identity)
total :=cat(2,transpose(nu),identity(nF));
// Gauss Jordan algorithm with total
for i in 1:nR loop
if not total[i,i] <> 0 then
(iDummy[i:nF],sort[i:nF]):=Modelica.Math.Vectors.sort(abs(total[i:nF, i]),false);
for j in i:nF loop
sort[j] :=sort[j] + (i-1);
end for;
total[i:nF,:] :=total[sort[i:nF], :];
end if;
total[i,:] :=total[i, :] / total[i, i];
for j in 1:nF loop
if i <> j then
total[j,:] :=total[j, :] - total[j, i]*total[i, :]/total[i, i];
end if;
end for;
end for;
// read transpose(lambda) = lambdaT from total
lambdaT :=total[nR + 1:nF, nR + 1:nR + nF];
// column changes (saved in resort1) and row operations to get rid of negativ elements in 1:nF-nR columns of lambdaT
for i in 1:nF-nR loop
if min(lambdaT) < 0 then
if not lambdaT[i,i] <> 0 then
(iDummy1[i:nF],sort1[i:nF]) :=Modelica.Math.Vectors.sort(abs(lambdaT[i,i:nF]), false);
for j in i:nF loop
sort1[j] :=sort1[j] + (i - 1);
end for;
lambdaT[:,i:nF] :=lambdaT[:, sort1[i:nF]];
resortindex[i] :=true;
(iDummy1[i:nF],resort1[i:nF,i]) :=Modelica.Math.Vectors.sort(sort1[i:nF], true);
for j in i:nF loop
resort1[j,i] :=resort1[j,i] + (i-1);
end for;
end if;
lambdaT[i,:] :=lambdaT[i, :] / lambdaT[i, i];
for j in 1:nF-nR loop
if i <> j then
lambdaT[j,:] :=lambdaT[j, :] - lambdaT[j, i]*lambdaT[i, :]/lambdaT[i, i];
end if;
end for;
else
break;
end if;
end for;
// get rid of negativ elements in nR+1:nF columns
for i in nR+1:nF loop
if min(lambdaT[:,i]) < 0 then
for j in 1:nF-nR loop
if lambdaT[j,i] < 0 then
index :=Modelica.Math.BooleanVectors.firstTrueIndex({lambdaT[k,i] > 0 for k in 1:nF-nR});
if index > 0 then
lambdaT[j,:] :=lambdaT[j, :] + lambdaT[index, :] * abs(lambdaT[j,i])/lambdaT[index,i];
index :=0;
end if;
end if;
end for;
end if;
end for;
// resort columns accordingly
for i in 2:nF loop
if resortindex[nF+1-i] then
lambdaT[:,nF+1-i:nF] :=lambdaT[:, resort1[nF + 1 - i:nF,nF+1-i]];
end if;
end for;
// transpose of lambdaT
lambda :=transpose(lambdaT);
//eliminate rounding errors
for i in 1:size(lambda,1) loop
for j in 1:size(lambda,2) loop
if lambda[i,j] < 1e-15 then
lambda[i,j] :=0;
end if;
end for;
end for;
else
lambda :=identity(nF);
end if;
// deprecated
// lambdaT :=total[nR + 1:nF, nR + 1:nR + nF];
// for i in 1:nF-nR loop
// if min(lambdaT) < 0 then
// if not lambdaT[i,i] <> 0 then
// (iDummy[i:nF-nR],sort[i:nF-nR]):=Modelica.Math.Vectors.sort(abs(lambdaT[i:nF-nR, i]),false);
// for j in i:nF-nR loop
// sort[j] :=sort[j] + (i-1);
// end for;
// lambdaT[i:nF-nR,:] :=lambdaT[sort[i:nF-nR], :];
// end if;
// lambdaT[i,:] :=lambdaT[i, :] / lambdaT[i, i];
// for j in 1:nF-nR loop
// if i <> j then
// lambdaT[j,:] :=lambdaT[j, :] - lambdaT[j, i]*lambdaT[i, :]/lambdaT[i, i];
// end if;
// end for;
// else
// break;
// end if;
// end for;
// for i in 1:nF-nR loop
// if min(lambdaT) < 0 then
// (iDummy[i:nF-nR],sort[i:nF-nR]):=Modelica.Math.Vectors.sort(lambdaT[i:nF-nR, i+nR],false);
// for j in i:nF-nR loop
// sort[j] :=sort[j] + (i-1);
// end for;
// lambdaT[i:nF-nR,:] :=lambdaT[sort[i:nF-nR], :];
// for j in 1:nF-nR loop
// if lambdaT[j,i+nR] < 0 then
// lambdaT[j,:] :=lambdaT[j, :] - lambdaT[j, i+nR]*lambdaT[i, :]/lambdaT[i, i+nR];
// end if;
// end for;
// else
// break;
// end if;
// end for;
end calc_lambda;
ElectrolyteMedia/Media/Common/Reaction/calc_lambda_id.mo 0 → 100644
View file @ eaaf9ad7
within ElectrolyteMedia.Media.Common.Reaction;
function calc_lambda_id
"calculates lambda with nX*nX identity matrix"
input Real[:,:] lambda;
output Real[size(lambda,1),size(lambda,2)] lambda_id;
output Integer[size(lambda,1),size(lambda,1)] P_id;
protected
Integer nF = size(lambda,1);
Integer nX = size(lambda,2);
Real[nF,nX] LU;
Integer[nX] pivots_LAPACK;
Real[nF,nX] L;
Integer[nF,nF] P;
Real[nX,nX] Lsquare;
Real[nX,nX] LsquareT;
Real[nX,nX] LU_;
Integer[nX] pivots_LAPACK_;
Real[nX,nX] U_;
algorithm
(LU,pivots_LAPACK) :=Modelica.Math.Matrices.LU(lambda);
L :=Media.Common.Reaction.calc_L(LU);
P :=Media.Common.Reaction.calc_P(pivots_LAPACK,nF);
Lsquare :=L[1:nX, :];
LsquareT :=transpose(Lsquare);
(LU_,pivots_LAPACK_) :=Modelica.Math.Matrices.LU(LsquareT);
U_ :=Media.Common.Reaction.calc_U(LU_);
lambda_id :=L*Modelica.Math.Matrices.inv(transpose(U_));
P_id :=P;
end calc_lambda_id;
ElectrolyteMedia/Media/Common/Reaction/calc_nu_id.mo 0 → 100644
View file @ eaaf9ad7
within ElectrolyteMedia.Media.Common.Reaction;
function calc_nu_id
"Calculates a stoichiometry matrix with pivoting so that nuout =(I|nu_var) with permutation matrix Pout for reordering species elements of nu_ordered = nuout*Pout"
input Real[:,:] nu;
output Real[size(nu,1),size(nu,2)] nuout;
output Integer[size(nu,2),size(nu,2)] Pout;
protected
parameter Integer nR = size(nu,1);
parameter Integer nF = size(nu,2);
Real[nF,nR] nuT = transpose(nu);
Real[nR,nF] nuvar;
Real[nF,nR] LU;
Integer[nR] pivots_LAPACK;
Integer info;
Integer[nF] pivots;
Real[nF,nR] L;
Real[nR,nR] U;
Integer[nF,nF] P;
Real[nR,nR] Lsquare;
Real[nR,nR] LsquareT;
Real[nR,nR] LU_;
Integer[nR] pivots_LAPACK_;
Real[nR,nR] U_;
Integer[nF,nF] P_;
Real[nF,nR] nuT_;
algorithm
(LU,pivots_LAPACK,info) :=Modelica.Math.Matrices.LU(nuT);
// compute lower L from LU
L :=calc_L(LU);
// compute permutation matrix P from pivots_LAPACK
P :=calc_P(pivots_LAPACK[1:nR],nF);
// compute identity matrix of square of L such that L_ = (I|L') with L = P_*L_*U_
Lsquare :=L[1:nR, :];
LsquareT :=transpose(Lsquare);
(LU_,pivots_LAPACK_) :=Modelica.Math.Matrices.LU(LsquareT);
U_ :=calc_U(LU_);
P_ :=identity(nF);
P_[1:nR,1:nR] :=calc_P(pivots_LAPACK_, nR);
nuT_ :=L*Modelica.Math.Matrices.inv(transpose(U_));
nuout :=transpose(nuT_);
Pout :=P*transpose(P_);
end calc_nu_id;
ElectrolyteMedia/Media/Common/Reaction/package.mo 0 → 100644
View file @ eaaf9ad7
within ElectrolyteMedia.Media.Common;
package Reaction "Contains functions regarding chemical reactions"
end Reaction;
ElectrolyteMedia/Media/Common/Reaction/package.order 0 → 100644
View file @ eaaf9ad7
calc_lambda
calc_nu_id
calc_P
calc_L
calc_U
calc_P_nu_id
calc_lambda_id
calc_P_lambda_id
ElectrolyteMedia/Media/Common/Types/GasModel.mo 0 → 100644
View file @ eaaf9ad7
within ElectrolyteMedia.Media.Common.Types;
type GasModel = enumeration(
Ideal "Ideal gas",
PengRobinson "Peng Robinson EOS") "Lists all types of gas models";
ElectrolyteMedia/Media/Common/Types/LiquidModel.mo 0 → 100644
View file @ eaaf9ad7
within ElectrolyteMedia.Media.Common.Types;
type LiquidModel = enumeration(
Ideal "Ideal liquid phase at infinite dilution",
Debye "Debye limiting law",
DebyeHueckel "Debye-Hueckel model",
Bromley "Bromley model",
Pitzer "Pitzer model") "Lists all types of liquid models";
ElectrolyteMedia/Media/Common/Types/package.mo 0 → 100644
View file @ eaaf9ad7
within ElectrolyteMedia.Media.Common;
package Types "Contains enumerations of all types of models"
end Types;