ISAT::chemPointISAT: Use scalarMatrices::SVD
This commit is contained in:
Henry Weller
parent
5a26fb0338
commit
3420fee771
4 changed files with 132 additions and 469 deletions
|
|
@ -48,11 +48,11 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
|
|||
scalar scale = 0;
|
||||
scalar s = 0;
|
||||
scalar anorm = 0;
|
||||
label l=0;
|
||||
label l = 0;
|
||||
|
||||
for (label i=0; i<Un; i++)
|
||||
{
|
||||
l = i+2;
|
||||
l = i + 2;
|
||||
rv1[i] = scale*g;
|
||||
g = s = scale = 0;
|
||||
|
||||
|
|
@ -72,7 +72,7 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
|
|||
}
|
||||
|
||||
scalar f = U_(i, i);
|
||||
g = -sign(Foam::sqrt(s), f);
|
||||
g = -sign(sqrt(s), f);
|
||||
scalar h = f*g - s;
|
||||
U_(i, i) = f - g;
|
||||
|
||||
|
|
@ -118,7 +118,7 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
|
|||
}
|
||||
|
||||
scalar f = U_[i][l-1];
|
||||
g = -sign(Foam::sqrt(s),f);
|
||||
g = -sign(sqrt(s), f);
|
||||
scalar h = f*g - s;
|
||||
U_[i][l-1] = f - g;
|
||||
|
||||
|
|
@ -150,6 +150,8 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
|
|||
anorm = max(anorm, mag(S_[i]) + mag(rv1[i]));
|
||||
}
|
||||
|
||||
anorm *= SMALL;
|
||||
|
||||
for (label i=Un-1; i >= 0; i--)
|
||||
{
|
||||
if (i < Un-1)
|
||||
|
|
@ -178,7 +180,7 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
|
|||
|
||||
for (label j=l; j<Un;j++)
|
||||
{
|
||||
V_(i, j) = V_(j, i) = 0.0;
|
||||
V_(i, j) = V_(j, i) = 0;
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -194,7 +196,7 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
|
|||
|
||||
for (label j=l; j<Un; j++)
|
||||
{
|
||||
U_(i, j) = 0.0;
|
||||
U_(i, j) = 0;
|
||||
}
|
||||
|
||||
if (g != 0)
|
||||
|
|
@ -226,7 +228,7 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
|
|||
{
|
||||
for (label j=i; j<Um; j++)
|
||||
{
|
||||
U_(j, i) = 0.0;
|
||||
U_(j, i) = 0;
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -235,32 +237,40 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
|
|||
|
||||
for (label k=Un-1; k >= 0; k--)
|
||||
{
|
||||
for (label its = 0; its < 35; its++)
|
||||
for (label its = 0; its < 30; its++)
|
||||
{
|
||||
bool flag = true;
|
||||
|
||||
label mn;
|
||||
for (l = k; l >= 0; l--)
|
||||
{
|
||||
mn = l-1;
|
||||
if (mag(rv1[l]) + anorm == anorm)
|
||||
mn = l - 1;
|
||||
|
||||
if (l == 0 || mag(rv1[l]) <= anorm)
|
||||
{
|
||||
flag = false;
|
||||
break;
|
||||
}
|
||||
if (mag(S_[mn]) + anorm == anorm) break;
|
||||
|
||||
if (mag(S_[mn]) <= anorm)
|
||||
{
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
if (flag)
|
||||
{
|
||||
scalar c = 0.0;
|
||||
s = 1.0;
|
||||
scalar c = 0;
|
||||
s = 1;
|
||||
for (label i=l; i<k+1; i++)
|
||||
{
|
||||
scalar f = s*rv1[i];
|
||||
rv1[i] = c*rv1[i];
|
||||
|
||||
if (mag(f) + anorm == anorm) break;
|
||||
if (mag(f) <= anorm)
|
||||
{
|
||||
break;
|
||||
}
|
||||
|
||||
g = S_[i];
|
||||
scalar h = sqrtSumSqr(f, g);
|
||||
|
|
@ -283,18 +293,20 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
|
|||
|
||||
if (l == k)
|
||||
{
|
||||
if (z < 0.0)
|
||||
if (z < 0)
|
||||
{
|
||||
S_[k] = -z;
|
||||
|
||||
for (label j=0; j<Un; j++) V_(j, k) = -V_(j, k);
|
||||
for (label j=0; j<Un; j++)
|
||||
{
|
||||
V_(j, k) = -V_(j, k);
|
||||
}
|
||||
}
|
||||
break;
|
||||
}
|
||||
if (its == 34)
|
||||
if (its == 29)
|
||||
{
|
||||
WarningInFunction
|
||||
<< "No convergence in 35 SVD iterations"
|
||||
<< "No convergence in 30 SVD iterations"
|
||||
<< endl;
|
||||
}
|
||||
|
||||
|
|
@ -303,15 +315,15 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
|
|||
scalar y = S_[mn];
|
||||
g = rv1[mn];
|
||||
scalar h = rv1[k];
|
||||
scalar f = ((y - z)*(y + z) + (g - h)*(g + h))/(2.0*h*y);
|
||||
scalar f = ((y - z)*(y + z) + (g - h)*(g + h))/(2*h*y);
|
||||
g = sqrtSumSqr(f, scalar(1));
|
||||
f = ((x - z)*(x + z) + h*((y/(f + sign(g, f))) - h))/x;
|
||||
scalar c = 1.0;
|
||||
s = 1.0;
|
||||
scalar c = 1;
|
||||
s = 1;
|
||||
|
||||
for (label j=l; j <= mn; j++)
|
||||
{
|
||||
label i=j + 1;
|
||||
label i = j + 1;
|
||||
g = rv1[i];
|
||||
y = S_[i];
|
||||
h = s*g;
|
||||
|
|
@ -352,7 +364,7 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
|
|||
U_[jj][i] = z*c - y*s;
|
||||
}
|
||||
}
|
||||
rv1[l] = 0.0;
|
||||
rv1[l] = 0;
|
||||
rv1[k] = f;
|
||||
S_[k] = x;
|
||||
}
|
||||
|
|
@ -364,7 +376,6 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
|
|||
{
|
||||
if (S_[i] <= minS)
|
||||
{
|
||||
//Info<< "Removing " << S_[i] << " < " << minS << endl;
|
||||
S_[i] = 0;
|
||||
nZeros_++;
|
||||
}
|
||||
|
|
@ -372,28 +383,6 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
|
|||
|
||||
// Now multiply out to find the pseudo inverse of A, VSinvUt_
|
||||
multiply(VSinvUt_, V_, inv(S_), U_.T());
|
||||
|
||||
// test SVD
|
||||
/*
|
||||
scalarRectangularMatrix SVDA(A.m(), A.n());
|
||||
multiply(SVDA, U_, S_, transpose(V_));
|
||||
scalar maxDiff = 0;
|
||||
scalar diff = 0;
|
||||
for (label i=0; i<A.m(); i++)
|
||||
{
|
||||
for (label j=0; j<A.n(); j++)
|
||||
{
|
||||
diff = mag(A(i, j) - SVDA(i, j));
|
||||
if (diff > maxDiff) maxDiff = diff;
|
||||
}
|
||||
}
|
||||
Info<< "Maximum discrepancy between A and svd(A) = " << maxDiff << endl;
|
||||
|
||||
if (maxDiff > 4)
|
||||
{
|
||||
Info<< "singular values " << S_ << endl;
|
||||
}
|
||||
*/
|
||||
}
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -28,7 +28,8 @@ License
|
|||
template<class T>
|
||||
inline const T Foam::SVD::sign(const T& a, const T& b)
|
||||
{
|
||||
return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a);
|
||||
//return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a);
|
||||
return b >= 0 ? a : -a;
|
||||
}
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -24,7 +24,7 @@ License
|
|||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#include "chemPointISAT.H"
|
||||
#include <limits>
|
||||
#include "SVD.H"
|
||||
|
||||
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
|
||||
|
||||
|
|
@ -49,14 +49,14 @@ void Foam::chemPointISAT<CompType, ThermoType>::qrDecompose
|
|||
|
||||
for (label k=0; k<nCols-1; k++)
|
||||
{
|
||||
scale = 0.0;
|
||||
scale = 0;
|
||||
for (label i=k; i<nCols; i++)
|
||||
{
|
||||
scale=max(scale, fabs(R(i, k)));
|
||||
scale=max(scale, mag(R(i, k)));
|
||||
}
|
||||
if (scale == 0.0)
|
||||
if (scale == 0)
|
||||
{
|
||||
c[k] = d[k] = 0.0;
|
||||
c[k] = d[k] = 0;
|
||||
}
|
||||
else
|
||||
{
|
||||
|
|
@ -64,7 +64,7 @@ void Foam::chemPointISAT<CompType, ThermoType>::qrDecompose
|
|||
{
|
||||
R(i, k) /= scale;
|
||||
}
|
||||
sum = 0.0;
|
||||
sum = 0;
|
||||
for (label i=k; i<nCols; i++)
|
||||
{
|
||||
sum += sqr(R(i, k));
|
||||
|
|
@ -75,7 +75,7 @@ void Foam::chemPointISAT<CompType, ThermoType>::qrDecompose
|
|||
d[k] = -scale*sigma;
|
||||
for (label j=k+1; j<nCols; j++)
|
||||
{
|
||||
sum=0.0;
|
||||
sum=0;
|
||||
for ( label i=k; i<nCols; i++)
|
||||
{
|
||||
sum += R(i, k)*R(i, j);
|
||||
|
|
@ -88,14 +88,16 @@ void Foam::chemPointISAT<CompType, ThermoType>::qrDecompose
|
|||
}
|
||||
}
|
||||
}
|
||||
|
||||
d[nCols-1] = R(nCols-1, nCols-1);
|
||||
|
||||
// form R
|
||||
for (label i=0; i<nCols; i++)
|
||||
{
|
||||
R(i, i) = d[i];
|
||||
for ( label j=0; j<i; j++)
|
||||
{
|
||||
R(i, j)=0.0;
|
||||
R(i, j)=0;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -110,42 +112,47 @@ void Foam::chemPointISAT<CompType, ThermoType>::qrUpdate
|
|||
const Foam::scalarField &v
|
||||
)
|
||||
{
|
||||
label k, i;
|
||||
label k;
|
||||
|
||||
scalarField w(u);
|
||||
for (k=n-1;k>=0;k--)
|
||||
for (k=n-1; k>=0; k--)
|
||||
{
|
||||
if (w[k] != 0.0)
|
||||
if (w[k] != 0)
|
||||
{
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
if (k < 0)
|
||||
{
|
||||
k=0;
|
||||
k = 0;
|
||||
}
|
||||
for (i=k-1;i>=0;i--)
|
||||
|
||||
for (label i=k-1; i>=0; i--)
|
||||
{
|
||||
rotate(R, i, w[i],-w[i+1], n);
|
||||
if (w[i] == 0.0)
|
||||
if (w[i] == 0)
|
||||
{
|
||||
w[i] = fabs(w[i+1]);
|
||||
w[i] = mag(w[i+1]);
|
||||
}
|
||||
else if (fabs(w[i]) > fabs(w[i+1]))
|
||||
else if (mag(w[i]) > mag(w[i+1]))
|
||||
{
|
||||
w[i] = fabs(w[i])*sqrt(1.0+sqr(w[i+1]/w[i]));
|
||||
w[i] = mag(w[i])*sqrt(1 + sqr(w[i+1]/w[i]));
|
||||
}
|
||||
else
|
||||
{
|
||||
w[i] = fabs(w[i+1])*sqrt(1.0+sqr(w[i]/w[i+1]));
|
||||
w[i] = mag(w[i+1])*sqrt(1 + sqr(w[i]/w[i+1]));
|
||||
}
|
||||
}
|
||||
for (i=0;i<n;i++)
|
||||
|
||||
for (label i=0; i<n; i++)
|
||||
{
|
||||
R(0, i) += w[0]*v[i];
|
||||
}
|
||||
for (i=0;i<k;i++)
|
||||
|
||||
for (label i=0; i<k; i++)
|
||||
{
|
||||
rotate(R, i, R(i, i),-R(i+1, i), n);
|
||||
rotate(R, i, R(i, i), -R(i+1, i), n);
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -160,360 +167,36 @@ void Foam::chemPointISAT<CompType, ThermoType>::rotate
|
|||
label n
|
||||
)
|
||||
{
|
||||
label j;
|
||||
scalar c, fact, s, w, y;
|
||||
if (a == 0.0)
|
||||
if (a == 0)
|
||||
{
|
||||
c=0.0;
|
||||
s=(b >= 0.0 ? 1.0 : -1.0);
|
||||
c = 0;
|
||||
s = (b >= 0 ? 1.0 : -1.0);
|
||||
}
|
||||
else if (fabs(a) > fabs(b))
|
||||
else if (mag(a) > mag(b))
|
||||
{
|
||||
fact = b/a;
|
||||
c=sign(a)/sqrt(1.0+(fact*fact));
|
||||
s=fact*c;
|
||||
c = sign(a)/sqrt(1 + sqr(fact));
|
||||
s = fact*c;
|
||||
}
|
||||
else
|
||||
{
|
||||
fact=a/b;
|
||||
s=sign(b)/sqrt(1.0+(fact*fact));
|
||||
c=fact*s;
|
||||
fact = a/b;
|
||||
s = sign(b)/sqrt(1 + sqr(fact));
|
||||
c = fact*s;
|
||||
}
|
||||
for (j=i;j<n;j++)
|
||||
for (label j=i;j<n;j++)
|
||||
{
|
||||
y=R(i, j);
|
||||
w=R(i+1, j);
|
||||
R(i, j)=c*y-s*w;
|
||||
R(i+1, j)=s*y+c*w;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template<class CompType, class ThermoType>
|
||||
void Foam::chemPointISAT<CompType, ThermoType>::svd
|
||||
(
|
||||
scalarSquareMatrix& A,
|
||||
label m,
|
||||
label n,
|
||||
scalarDiagonalMatrix& d,
|
||||
scalarSquareMatrix& V
|
||||
)
|
||||
{
|
||||
// UPDATED VERSION NR3
|
||||
bool flag;
|
||||
label i, its, j, jj, k, l, nm;
|
||||
scalar anorm, c, f, g, h, s, scale, x, y, z;
|
||||
scalarField rv1(n);
|
||||
scalar eps = std::numeric_limits<scalar>::epsilon();
|
||||
g = scale = anorm = 0.0;
|
||||
|
||||
// Householder reduction to bidiagonal form
|
||||
for ( i = 0; i<n; i++)
|
||||
{
|
||||
l=i+2; // change from i+1 to i+2
|
||||
rv1[i] = scale*g;
|
||||
g=s=scale=0.0;
|
||||
if (i < m)
|
||||
{
|
||||
for (k=i;k<m;k++)
|
||||
{
|
||||
scale += fabs(A(k, i));
|
||||
}
|
||||
if (scale != 0.0)
|
||||
{
|
||||
for ( k=i;k<m;k++)
|
||||
{
|
||||
A(k, i) /= scale;
|
||||
s += A(k, i)*A(k, i);
|
||||
}
|
||||
f = A(i, i);
|
||||
g = -sign(f)*sqrt(s);
|
||||
h = f*g-s;
|
||||
A(i, i)=f-g;
|
||||
for (j=l-1;j<n;j++)
|
||||
{
|
||||
for (s=0.0,k=i;k<m;k++)
|
||||
{
|
||||
s += A(k, i)*A(k, j);
|
||||
}
|
||||
f = s/h;
|
||||
for (k=i; k<m;k++)
|
||||
{
|
||||
A(k, j) += f*A(k, i);
|
||||
}
|
||||
}
|
||||
for (k=i; k<m;k++)
|
||||
{
|
||||
A(k, i) *= scale;
|
||||
}
|
||||
}
|
||||
}
|
||||
d[i] = scale * g;
|
||||
g=s=scale=0.0;
|
||||
|
||||
if (i+1 <= m && i+1 != n)
|
||||
{
|
||||
for (k=l-1; k<n; k++)
|
||||
{
|
||||
scale += fabs(A(i, k));
|
||||
}
|
||||
if (scale != 0.0)
|
||||
{
|
||||
for (k=l-1; k<n; k++)
|
||||
{
|
||||
A(i, k) /= scale;
|
||||
s += A(i, k)*A(i, k);
|
||||
}
|
||||
f = A(i, l-1);
|
||||
g = -sign(f)*sqrt(s);
|
||||
h = f*g-s;
|
||||
A(i, l-1) = f-g;
|
||||
for (k=l-1; k<n; k++)
|
||||
{
|
||||
rv1[k] = A(i, k)/h;
|
||||
}
|
||||
for (j=l-1; j<m; j++)
|
||||
{
|
||||
for (s=0.0,k=l-1; k<n; k++)
|
||||
{
|
||||
s += A(j, k)*A(i, k);
|
||||
}
|
||||
for (k=l-1; k<n; k++)
|
||||
{
|
||||
A(j, k) += s*rv1[k];
|
||||
}
|
||||
}
|
||||
for (k=l-1; k<n; k++)
|
||||
{
|
||||
A(i, k) *= scale;
|
||||
}
|
||||
}
|
||||
}
|
||||
anorm = max(anorm, (fabs(d[i])+fabs(rv1[i])));
|
||||
}
|
||||
|
||||
// Accumulation of right-hand transformations
|
||||
for (i=n-1; i>=0; i--)
|
||||
{
|
||||
if (i < n-1)
|
||||
{
|
||||
if (g != 0.0)
|
||||
{
|
||||
for (j=l; j<n; j++)
|
||||
{
|
||||
V(j, i) = (A(i, j)/A(i, l))/g;
|
||||
}
|
||||
for (j=l; j<n; j++)
|
||||
{
|
||||
for (s=0.0,k=l; k<n; k++)
|
||||
{
|
||||
s += A(i, k)*V(k, j);
|
||||
}
|
||||
for (k=l; k<n; k++)
|
||||
{
|
||||
V(k, j) += s*V(k, i);
|
||||
}
|
||||
}
|
||||
}
|
||||
for (j=l; j<n; j++)
|
||||
{
|
||||
V(i, j)=V(j, i)=0.0;
|
||||
}
|
||||
}
|
||||
V(i, i) = 1.0;
|
||||
g = rv1[i];
|
||||
l = i;
|
||||
}
|
||||
// Accumulation of left-hand transformations
|
||||
for (i = min(m, n)-1; i>=0; i--)
|
||||
{
|
||||
l=i+1;
|
||||
g=d[i];
|
||||
for (j=l; j<n; j++)
|
||||
{
|
||||
A(i, j) = 0.0;
|
||||
}
|
||||
if (g != 0.0)
|
||||
{
|
||||
g = 1.0/g;
|
||||
for (j=l; j<n; j++)
|
||||
{
|
||||
for (s=0.0, k=l; k<m; k++)
|
||||
{
|
||||
s+= A(k, i)*A(k, j);
|
||||
}
|
||||
f = (s/A(i, i))*g;
|
||||
for (k=i; k<m; k++)
|
||||
{
|
||||
A(k, j) += f*A(k, i);
|
||||
}
|
||||
}
|
||||
for (j=i; j<m; j++)
|
||||
{
|
||||
A(j, i) *= g;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
for (j=i; j<m; j++)
|
||||
{
|
||||
A(j, i)=0.0;
|
||||
}
|
||||
}
|
||||
++A(i, i);
|
||||
}
|
||||
|
||||
// Diagonalization of the bidiagonal form :
|
||||
// Loop over singular values, and over allowed iteration
|
||||
for (k=n-1; k>=0; k--)
|
||||
{
|
||||
for (its=0; its<30; its++)
|
||||
{
|
||||
flag=true;
|
||||
// Test for splitting (rv1[1] always zero)
|
||||
for (l=k; l>=0; l--)
|
||||
{
|
||||
nm = l-1;
|
||||
if (l == 0 || fabs(rv1[l]) <= eps*anorm)
|
||||
{
|
||||
flag = false;
|
||||
break;
|
||||
}
|
||||
if (fabs(d[nm]) <= eps*anorm)
|
||||
{
|
||||
break;
|
||||
}
|
||||
}
|
||||
// Cancellation of rv1[l], if l>1
|
||||
if (flag)
|
||||
{
|
||||
c = 0.0;
|
||||
s = 1.0;
|
||||
for (i=l; i<k+1; i++)
|
||||
{
|
||||
f = s*rv1[i];
|
||||
rv1[i] = c*rv1[i];
|
||||
if (fabs(f) <= eps*anorm)
|
||||
{
|
||||
break;
|
||||
}
|
||||
g = d[i];
|
||||
h = pythag(f, g);
|
||||
d[i] = h;
|
||||
h = 1.0/h;
|
||||
c = g*h;
|
||||
s = -f*h;
|
||||
for (j=0; j<m; j++)
|
||||
{
|
||||
y = A(j, nm);
|
||||
z = A(j, i);
|
||||
A(j, nm) = y*c + z*s;
|
||||
A(j, i) = z*c - y*s;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
z = d[k];
|
||||
if (l == k) // Convergence
|
||||
{
|
||||
if (z < 0.0) // Singular value is made nonnegative
|
||||
{
|
||||
d[k] = -z;
|
||||
for (j=0; j<n; j++)
|
||||
{
|
||||
V(j, k) = -V(j, k);
|
||||
}
|
||||
}
|
||||
break;
|
||||
}
|
||||
if (its == 34)
|
||||
{
|
||||
WarningInFunction
|
||||
<< "no convergence in 35 SVD iterations"
|
||||
<< endl;
|
||||
}
|
||||
|
||||
x = d[l];
|
||||
nm = k-1;
|
||||
y = d[nm];
|
||||
g = rv1[nm];
|
||||
h = rv1[k];
|
||||
f = ((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
|
||||
g = pythag(f,1.0);
|
||||
f = ((x-z)*(x+z)+h*((y/(f+sign(f)*g))-h))/x;
|
||||
c=s=1.0;
|
||||
// Next QR transformation
|
||||
for (j=l; j<=nm; j++)
|
||||
{
|
||||
i = j+1;
|
||||
g = rv1[i];
|
||||
y = d[i];
|
||||
h = s*g;
|
||||
g = c*g;
|
||||
z = pythag(f, h);
|
||||
rv1[j] = z;
|
||||
c = f/z;
|
||||
s = h/z;
|
||||
f = x*c + g*s;
|
||||
g = g*c - x*s;
|
||||
h = y*s;
|
||||
y *= c;
|
||||
for (jj=0; jj<n; jj++)
|
||||
{
|
||||
x = V(jj, j);
|
||||
z = V(jj, i);
|
||||
V(jj, j) = x*c + z*s;
|
||||
V(jj, i) = z*c - x*s;
|
||||
}
|
||||
z = pythag(f, h);
|
||||
d[j] = z;
|
||||
if (z)
|
||||
{
|
||||
z = 1.0/z;
|
||||
c = f*z;
|
||||
s = h*z;
|
||||
}
|
||||
f = c*g + s*y;
|
||||
x = c*y - s*g;
|
||||
for (jj=0; jj<m; jj++)
|
||||
{
|
||||
y = A(jj, j);
|
||||
z = A(jj, i);
|
||||
A(jj, j) = y*c + z*s;
|
||||
A(jj, i) = z*c - y*s;
|
||||
}
|
||||
}
|
||||
rv1[l] = 0.0;
|
||||
rv1[k] = f;
|
||||
d[k] = x;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// pythag function used in svd
|
||||
// compute (a^2+b^2)^1/2 without descrutive underflow or overflow
|
||||
template<class CompType, class ThermoType>
|
||||
Foam::scalar
|
||||
Foam::chemPointISAT<CompType, ThermoType>::pythag(scalar a, scalar b)
|
||||
{
|
||||
scalar absa, absb;
|
||||
absa = fabs(a);
|
||||
absb = fabs(b);
|
||||
if (absa > absb)
|
||||
{
|
||||
return absa*sqrt(1.0+sqr(absb/absa));
|
||||
}
|
||||
else
|
||||
{
|
||||
return (absb == 0.0 ? 0.0 : absb*sqrt(1.0+sqr(absa/absb)));
|
||||
y = R(i, j);
|
||||
w = R(i+1, j);
|
||||
R(i, j) = c*y-s*w;
|
||||
R(i+1, j) = s*y+c*w;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
|
||||
|
||||
|
||||
template<class CompType, class ThermoType>
|
||||
Foam::chemPointISAT<CompType, ThermoType>::chemPointISAT
|
||||
(
|
||||
|
|
@ -565,37 +248,38 @@ Foam::chemPointISAT<CompType, ThermoType>::chemPointISAT
|
|||
label reduOrCompDim = completeSpaceSize;
|
||||
if (isMechRedActive)
|
||||
{
|
||||
reduOrCompDim = nActiveSpecies_+2;
|
||||
reduOrCompDim = nActiveSpecies_ + 2;
|
||||
}
|
||||
|
||||
// SVD decomposition A= U*D*V^T
|
||||
scalarSquareMatrix Atmp(A);// A computed in ISAT.C
|
||||
scalarSquareMatrix B(reduOrCompDim, Zero);
|
||||
DiagonalMatrix<scalar> diag(reduOrCompDim, Zero);
|
||||
svd(Atmp, reduOrCompDim, reduOrCompDim, diag, B);
|
||||
// SVD decomposition A = U*D*V^T
|
||||
SVD svdA(A);
|
||||
|
||||
// replace the value of vector diag by max(diag, 1/2), first ISAT paper,
|
||||
// Pope
|
||||
scalarDiagonalMatrix D(reduOrCompDim);
|
||||
const scalarDiagonalMatrix& S = svdA.S();
|
||||
|
||||
// Replace the value of vector D by max(D, 1/2), first ISAT paper
|
||||
for (label i=0; i<reduOrCompDim; i++)
|
||||
{
|
||||
diag[i] = max(diag[i], 0.5);
|
||||
D[i] = max(S[i], 0.5);
|
||||
}
|
||||
|
||||
// rebuild A with max length, tol and scale factor before QR decomposition
|
||||
scalarSquareMatrix Atilde(reduOrCompDim, reduOrCompDim);
|
||||
// Rebuild A with max length, tol and scale factor before QR decomposition
|
||||
scalarRectangularMatrix Atilde(reduOrCompDim);
|
||||
|
||||
// result stored in Atilde
|
||||
multiply(Atilde, Atmp, diag, B.T());
|
||||
// Result stored in Atilde
|
||||
multiply(Atilde, svdA.U(), D, svdA.V().T());
|
||||
|
||||
for (label i=0; i<reduOrCompDim; i++)// on species loop
|
||||
for (label i=0; i<reduOrCompDim; i++)
|
||||
{
|
||||
for (label j=0; j<reduOrCompDim; j++)// species, T and p loop
|
||||
for (label j=0; j<reduOrCompDim; j++)
|
||||
{
|
||||
label compi=i;
|
||||
label compi = i;
|
||||
|
||||
if (isMechRedActive)
|
||||
{
|
||||
compi = simplifiedToCompleteIndex(i);
|
||||
}
|
||||
|
||||
// SF*A/tolerance
|
||||
// (where SF is diagonal with inverse of scale factors)
|
||||
// SF*A is the same as dividing each line by the scale factor
|
||||
|
|
@ -639,6 +323,7 @@ Foam::chemPointISAT<CompType, ThermoType>::chemPointISAT
|
|||
tolerance_ = p.tolerance();
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
||||
|
||||
template<class CompType, class ThermoType>
|
||||
|
|
@ -647,27 +332,30 @@ bool Foam::chemPointISAT<CompType, ThermoType>::inEOA(const scalarField& phiq)
|
|||
scalarField dphi(phiq-phi());
|
||||
bool isMechRedActive = chemistry_.mechRed()->active();
|
||||
label dim = (isMechRedActive) ? nActiveSpecies_ : completeSpaceSize()-2;
|
||||
scalar epsTemp=0.0;
|
||||
List<scalar> propEps(completeSpaceSize(),0.0);
|
||||
scalar epsTemp=0;
|
||||
List<scalar> propEps(completeSpaceSize(),0);
|
||||
|
||||
for (label i=0; i<completeSpaceSize()-2; i++)
|
||||
{
|
||||
scalar temp(0.0);
|
||||
scalar temp = 0;
|
||||
|
||||
// When mechanism reduction is inactive OR on active species multiply L
|
||||
// by dphi to get the distance in the active species direction else (for
|
||||
// inactive species), just multiply the diagonal element and dphi
|
||||
if
|
||||
(
|
||||
!(isMechRedActive)
|
||||
||(isMechRedActive && completeToSimplifiedIndex_[i]!=-1)
|
||||
||(isMechRedActive && completeToSimplifiedIndex_[i] != -1)
|
||||
)
|
||||
{
|
||||
label si=(isMechRedActive) ? completeToSimplifiedIndex_[i] : i;
|
||||
|
||||
for (label j=si; j<dim; j++)// LT is upper triangular
|
||||
{
|
||||
label sj=(isMechRedActive) ? simplifiedToCompleteIndex_[j] : j;
|
||||
temp += LT_(si, j)*dphi[sj];
|
||||
}
|
||||
|
||||
temp += LT_(si, nActiveSpecies_)*dphi[completeSpaceSize()-2];
|
||||
temp += LT_(si, nActiveSpecies_+1)*dphi[completeSpaceSize()-1];
|
||||
}
|
||||
|
|
@ -675,12 +363,15 @@ bool Foam::chemPointISAT<CompType, ThermoType>::inEOA(const scalarField& phiq)
|
|||
{
|
||||
temp = dphi[i]/(tolerance_*scaleFactor_[i]);
|
||||
}
|
||||
|
||||
epsTemp += sqr(temp);
|
||||
|
||||
if (printProportion_)
|
||||
{
|
||||
propEps[i] = temp;
|
||||
}
|
||||
}
|
||||
|
||||
// Temperature
|
||||
epsTemp +=
|
||||
sqr
|
||||
|
|
@ -688,6 +379,7 @@ bool Foam::chemPointISAT<CompType, ThermoType>::inEOA(const scalarField& phiq)
|
|||
LT_(dim, dim)*dphi[completeSpaceSize()-2]
|
||||
+LT_(dim, dim+1)*dphi[completeSpaceSize()-1]
|
||||
);
|
||||
|
||||
// Pressure
|
||||
epsTemp += sqr(LT_(dim+1, dim+1)*dphi[completeSpaceSize()-1]);
|
||||
|
||||
|
|
@ -703,12 +395,12 @@ bool Foam::chemPointISAT<CompType, ThermoType>::inEOA(const scalarField& phiq)
|
|||
propEps[completeSpaceSize()-1] =
|
||||
sqr(LT_(dim+1, dim+1)*dphi[completeSpaceSize()-1]);
|
||||
}
|
||||
if (sqrt(epsTemp) > 1.0+tolerance_)
|
||||
if (sqrt(epsTemp) > 1 + tolerance_)
|
||||
{
|
||||
if (printProportion_)
|
||||
{
|
||||
scalar max=-1.0;
|
||||
label maxIndex=-1;
|
||||
scalar max = -1;
|
||||
label maxIndex = -1;
|
||||
for (label i=0; i<completeSpaceSize(); i++)
|
||||
{
|
||||
if(max < propEps[i])
|
||||
|
|
@ -736,7 +428,7 @@ bool Foam::chemPointISAT<CompType, ThermoType>::inEOA(const scalarField& phiq)
|
|||
Info<< "Direction maximum impact to error in ellipsoid: "
|
||||
<< propName << endl;
|
||||
Info<< "Proportion to the total error on the retrieve: "
|
||||
<< max / (epsTemp+SMALL) << endl;
|
||||
<< max/(epsTemp+SMALL) << endl;
|
||||
}
|
||||
return false;
|
||||
}
|
||||
|
|
@ -754,13 +446,13 @@ bool Foam::chemPointISAT<CompType, ThermoType>::checkSolution
|
|||
const scalarField& Rphiq
|
||||
)
|
||||
{
|
||||
scalar eps2 = 0.0;
|
||||
scalar eps2 = 0;
|
||||
scalarField dR(Rphiq - Rphi());
|
||||
scalarField dphi(phiq - phi());
|
||||
const scalarField& scaleFactorV(scaleFactor());
|
||||
const scalarSquareMatrix& Avar(A());
|
||||
bool isMechRedActive = chemistry_.mechRed()->active();
|
||||
scalar dRl = 0.0;
|
||||
scalar dRl = 0;
|
||||
label dim = completeSpaceSize()-2;
|
||||
if (isMechRedActive)
|
||||
{
|
||||
|
|
@ -771,12 +463,13 @@ bool Foam::chemPointISAT<CompType, ThermoType>::checkSolution
|
|||
// included
|
||||
for (label i=0; i<completeSpaceSize()-2; i++)
|
||||
{
|
||||
dRl = 0.0;
|
||||
dRl = 0;
|
||||
if (isMechRedActive)
|
||||
{
|
||||
label si = completeToSimplifiedIndex_[i];
|
||||
|
||||
// If this species is active
|
||||
if (si!=-1)
|
||||
if (si != -1)
|
||||
{
|
||||
for (label j=0; j<dim; j++)
|
||||
{
|
||||
|
|
@ -846,7 +539,7 @@ bool Foam::chemPointISAT<CompType, ThermoType>::grow(const scalarField& phiq)
|
|||
// corresponds to an inactive on the query side
|
||||
if
|
||||
(
|
||||
completeToSimplifiedIndex_[i]!=-1
|
||||
completeToSimplifiedIndex_[i] != -1
|
||||
&& chemistry_.completeToSimplifiedIndex()[i] == -1
|
||||
)
|
||||
{
|
||||
|
|
@ -860,7 +553,7 @@ bool Foam::chemPointISAT<CompType, ThermoType>::grow(const scalarField& phiq)
|
|||
(
|
||||
completeToSimplifiedIndex_[i] == -1
|
||||
&& chemistry_.completeToSimplifiedIndex()[i] == -1
|
||||
&& dphi[i] != 0.0
|
||||
&& dphi[i] != 0
|
||||
)
|
||||
{
|
||||
activeAdded++;
|
||||
|
|
@ -933,17 +626,17 @@ bool Foam::chemPointISAT<CompType, ThermoType>::grow(const scalarField& phiq)
|
|||
{
|
||||
LT_(i, i)=
|
||||
1.0
|
||||
/ (tolerance_*scaleFactor_[simplifiedToCompleteIndex_[i]]);
|
||||
A_(i, i)=1.0;
|
||||
/(tolerance_*scaleFactor_[simplifiedToCompleteIndex_[i]]);
|
||||
A_(i, i) = 1;
|
||||
}
|
||||
}
|
||||
|
||||
dim = nActiveSpecies_+2;
|
||||
dim = nActiveSpecies_ + 2;
|
||||
}
|
||||
|
||||
// beginning of grow algorithm
|
||||
scalarField phiTilde(dim, 0.0);
|
||||
scalar normPhiTilde = 0.0;
|
||||
scalarField phiTilde(dim, 0);
|
||||
scalar normPhiTilde = 0;
|
||||
// p' = L^T.(p-phi)
|
||||
|
||||
for (label i=0; i<dim; i++)
|
||||
|
|
@ -953,7 +646,7 @@ bool Foam::chemPointISAT<CompType, ThermoType>::grow(const scalarField& phiq)
|
|||
label sj = j;
|
||||
if (isMechRedActive)
|
||||
{
|
||||
sj=simplifiedToCompleteIndex_[j];
|
||||
sj = simplifiedToCompleteIndex_[j];
|
||||
}
|
||||
phiTilde[i] += LT_(i, j)*dphi[sj];
|
||||
}
|
||||
|
|
@ -961,12 +654,15 @@ bool Foam::chemPointISAT<CompType, ThermoType>::grow(const scalarField& phiq)
|
|||
phiTilde[i] += LT_(i, dim-1)*dphi[completeSpaceSize()-1];
|
||||
normPhiTilde += sqr(phiTilde[i]);
|
||||
}
|
||||
|
||||
scalar invSqrNormPhiTilde = 1.0/normPhiTilde;
|
||||
normPhiTilde = sqrt(normPhiTilde);
|
||||
|
||||
// gamma = (1/|p'| - 1)/|p'|^2
|
||||
scalar gamma = (1/normPhiTilde - 1)*invSqrNormPhiTilde;
|
||||
scalarField u(gamma*phiTilde);
|
||||
scalarField v(dim,0.0);
|
||||
scalarField v(dim, 0);
|
||||
|
||||
for ( label i=0; i<dim; i++)
|
||||
{
|
||||
for (register label j=0; j<=i;j++)
|
||||
|
|
|
|||
|
|
@ -208,8 +208,8 @@ class chemPointISAT
|
|||
(
|
||||
scalarSquareMatrix& R,
|
||||
const label n,
|
||||
const scalarField &u,
|
||||
const scalarField &v
|
||||
const scalarField& u,
|
||||
const scalarField& v
|
||||
);
|
||||
|
||||
void rotate
|
||||
|
|
@ -221,29 +221,6 @@ class chemPointISAT
|
|||
label n
|
||||
);
|
||||
|
||||
//- Singular Value Decomposition (SVD) for a square matrix
|
||||
// needed to compute the the length of the hyperellipsoid semi-axes
|
||||
// SVD decompose a matrix A into:
|
||||
// A = U * D * V^T ,
|
||||
// with the singular value in the diagonal matrix D and U and V
|
||||
// orthogonal A (scalarMatrix) the square matrix to apply the
|
||||
// decomposition on m (label) the number of line of the matrix A
|
||||
// n (label) the size of the matrix A
|
||||
// Output: U (scalarMatrix) replace A
|
||||
// V (scalarMatrix) not the transpose V^T
|
||||
// d (scalarField) the diagonal element of matrix D
|
||||
void svd
|
||||
(
|
||||
scalarSquareMatrix& A,
|
||||
label m,
|
||||
label n,
|
||||
scalarDiagonalMatrix& d,
|
||||
scalarSquareMatrix& V
|
||||
);
|
||||
|
||||
//- Function used in svd function
|
||||
scalar pythag(scalar a, scalar b);
|
||||
|
||||
|
||||
public:
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue