From b722041fff3129f208eccde20cac205556be1002 Mon Sep 17 00:00:00 2001 From: henry Date: Sun, 28 Sep 2008 22:46:33 +0100 Subject: [PATCH 1/7] Corrected reading of "LESProperties". --- .../execFlowFunctionObjects/execFlowFunctionObjects.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/applications/utilities/postProcessing/miscellaneous/execFlowFunctionObjects/execFlowFunctionObjects.C b/applications/utilities/postProcessing/miscellaneous/execFlowFunctionObjects/execFlowFunctionObjects.C index b848e330ea..3f7a99040b 100644 --- a/applications/utilities/postProcessing/miscellaneous/execFlowFunctionObjects/execFlowFunctionObjects.C +++ b/applications/utilities/postProcessing/miscellaneous/execFlowFunctionObjects/execFlowFunctionObjects.C @@ -139,7 +139,7 @@ void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh) IOobject LESPropertiesHeader ( - "RASProperties", + "LESProperties", runTime.constant(), mesh, IOobject::MUST_READ, From 318b71e206e9e635aa649938e9d1c5dfbc0b7f6b Mon Sep 17 00:00:00 2001 From: henry Date: Sun, 28 Sep 2008 22:50:57 +0100 Subject: [PATCH 2/7] Made Matrix generic base-class for SquareMatrix and RectangularMatrix by introducing the type it will become as a template argument. Brought everything into line with this change. --- src/ODE/ODE/ODE.H | 4 +- src/ODE/ODESolvers/KRR4/KRR4.C | 19 +-- src/ODE/ODESolvers/KRR4/KRR4.H | 6 +- src/ODE/ODESolvers/SIBS/SIBS.H | 10 +- src/ODE/ODESolvers/SIBS/SIMPR.C | 13 +- src/ODE/ODESolvers/SIBS/polyExtrapolate.C | 2 +- src/OpenFOAM/Make/files | 6 +- .../IOstreams/Pstreams/PstreamCommsStruct.C | 2 - .../matrices/DiagonalMatrix/DiagonalMatrix.C | 3 +- .../matrices/DiagonalMatrix/DiagonalMatrix.H | 5 +- .../matrices/LUscalarMatrix/LUscalarMatrix.C | 8 +- .../matrices/LUscalarMatrix/LUscalarMatrix.H | 10 +- src/OpenFOAM/matrices/Matrix/Matrix.C | 106 ++++++------ src/OpenFOAM/matrices/Matrix/Matrix.H | 104 ++++++------ src/OpenFOAM/matrices/Matrix/MatrixI.H | 46 +++--- src/OpenFOAM/matrices/Matrix/MatrixIO.C | 24 +-- .../RectangularMatrix/RectangularMatrix.H | 92 +++++++++++ .../RectangularMatrix/RectangularMatrixI.H | 70 ++++++++ .../matrices/SquareMatrix/SquareMatrix.H | 96 +++++++++++ .../matrices/SquareMatrix/SquareMatrixI.H | 75 +++++++++ .../SVD/SVD.C | 13 +- .../SVD/SVD.H | 22 +-- .../SVD/SVDI.H | 8 +- .../scalarMatrices.C} | 99 +++++------ .../matrices/scalarMatrices/scalarMatrices.H | 137 ++++++++++++++++ .../scalarMatricesTemplates.C} | 25 +-- .../matrices/scalarMatrix/scalarMatrix.H | 155 ------------------ .../matrices/simpleMatrix/simpleMatrix.C | 28 ++-- .../matrices/simpleMatrix/simpleMatrix.H | 72 +++----- .../schemes/quadraticFit/quadraticFitData.C | 14 +- .../chemistryModel/chemistryModel.C | 38 ++--- .../chemistryModel/chemistryModel.H | 3 +- 32 files changed, 791 insertions(+), 524 deletions(-) create mode 100644 src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H create mode 100644 src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H create mode 100644 src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H create mode 100644 src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H rename src/OpenFOAM/matrices/{scalarMatrix => scalarMatrices}/SVD/SVD.C (96%) rename src/OpenFOAM/matrices/{scalarMatrix => scalarMatrices}/SVD/SVD.H (86%) rename src/OpenFOAM/matrices/{scalarMatrix => scalarMatrices}/SVD/SVDI.H (88%) rename src/OpenFOAM/matrices/{scalarMatrix/scalarMatrix.C => scalarMatrices/scalarMatrices.C} (76%) create mode 100644 src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.H rename src/OpenFOAM/matrices/{scalarMatrix/scalarMatrixTemplates.C => scalarMatrices/scalarMatricesTemplates.C} (90%) delete mode 100644 src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.H diff --git a/src/ODE/ODE/ODE.H b/src/ODE/ODE/ODE.H index 6b3c69011c..fe57cb91f4 100644 --- a/src/ODE/ODE/ODE.H +++ b/src/ODE/ODE/ODE.H @@ -34,7 +34,7 @@ Description #define ODE_H #include "scalarField.H" -#include "Matrix.H" +#include "scalarMatrices.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -80,7 +80,7 @@ public: const scalar x, const scalarField& y, scalarField& dfdx, - Matrix& dfdy + scalarSquareMatrix& dfdy ) const = 0; }; diff --git a/src/ODE/ODESolvers/KRR4/KRR4.C b/src/ODE/ODESolvers/KRR4/KRR4.C index fa00919a42..ec3db04574 100644 --- a/src/ODE/ODESolvers/KRR4/KRR4.C +++ b/src/ODE/ODESolvers/KRR4/KRR4.C @@ -25,7 +25,6 @@ License \*---------------------------------------------------------------------------*/ #include "KRR4.H" -#include "simpleMatrix.H" #include "addToRunTimeSelectionTable.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -36,11 +35,11 @@ namespace Foam { addToRunTimeSelectionTable(ODESolver, KRR4, ODE); -const scalar +const scalar KRR4::safety = 0.9, KRR4::grow = 1.5, KRR4::pgrow = -0.25, KRR4::shrink = 0.5, KRR4::pshrink = (-1.0/3.0), KRR4::errcon = 0.1296; -const scalar +const scalar KRR4::gamma = 1.0/2.0, KRR4::a21 = 2.0, KRR4::a31 = 48.0/25.0, KRR4::a32 = 6.0/25.0, KRR4::c21 = -8.0, KRR4::c31 = 372.0/25.0, KRR4::c32 = 12.0/5.0, @@ -81,8 +80,8 @@ void Foam::KRR4::solve const ODE& ode, scalar& x, scalarField& y, - scalarField& dydx, - const scalar eps, + scalarField& dydx, + const scalar eps, const scalarField& yScale, const scalar hTry, scalar& hDid, @@ -109,14 +108,14 @@ void Foam::KRR4::solve a_[i][i] += 1.0/(gamma*h); } - simpleMatrix::LUDecompose(a_, pivotIndices_); + LUDecompose(a_, pivotIndices_); for (register label i=0; i::LUBacksubstitute(a_, pivotIndices_, g1_); + LUBacksubstitute(a_, pivotIndices_, g1_); for (register label i=0; i::LUBacksubstitute(a_, pivotIndices_, g2_); + LUBacksubstitute(a_, pivotIndices_, g2_); for (register label i=0; i::LUBacksubstitute(a_, pivotIndices_, g3_); + LUBacksubstitute(a_, pivotIndices_, g3_); for (register label i=0; i::LUBacksubstitute(a_, pivotIndices_, g4_); + LUBacksubstitute(a_, pivotIndices_, g4_); for (register label i=0; i dfdy_; - mutable Matrix a_; + mutable scalarSquareMatrix dfdy_; + mutable scalarSquareMatrix a_; mutable labelList pivotIndices_; static const int maxtry = 40; static const scalar safety, grow, pgrow, shrink, pshrink, errcon; - static const scalar + static const scalar gamma, a21, a31, a32, c21, c31, c32, c41, c42, c43, diff --git a/src/ODE/ODESolvers/SIBS/SIBS.H b/src/ODE/ODESolvers/SIBS/SIBS.H index 75383bc38a..164b92d299 100644 --- a/src/ODE/ODESolvers/SIBS/SIBS.H +++ b/src/ODE/ODESolvers/SIBS/SIBS.H @@ -60,8 +60,8 @@ class SIBS static const scalar safe1, safe2, redMax, redMin, scaleMX; mutable scalarField a_; - mutable Matrix alpha_; - mutable Matrix d_p_; + mutable scalarSquareMatrix alpha_; + mutable scalarSquareMatrix d_p_; mutable scalarField x_p_; mutable scalarField err_; @@ -69,7 +69,7 @@ class SIBS mutable scalarField ySeq_; mutable scalarField yErr_; mutable scalarField dfdx_; - mutable Matrix dfdy_; + mutable scalarSquareMatrix dfdy_; mutable label first_, kMax_, kOpt_; mutable scalar epsOld_, xNew_; @@ -82,7 +82,7 @@ void SIMPR const scalarField& y, const scalarField& dydx, const scalarField& dfdx, - const Matrix& dfdy, + const scalarSquareMatrix& dfdy, const scalar deltaX, const label nSteps, scalarField& yEnd @@ -97,7 +97,7 @@ void polyExtrapolate scalarField& yz, scalarField& dy, scalarField& x_p, - Matrix& d_p + scalarSquareMatrix& d_p ) const; diff --git a/src/ODE/ODESolvers/SIBS/SIMPR.C b/src/ODE/ODESolvers/SIBS/SIMPR.C index f540f26d11..671f1acec9 100644 --- a/src/ODE/ODESolvers/SIBS/SIMPR.C +++ b/src/ODE/ODESolvers/SIBS/SIMPR.C @@ -25,7 +25,6 @@ License \*---------------------------------------------------------------------------*/ #include "SIBS.H" -#include "simpleMatrix.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -36,7 +35,7 @@ void Foam::SIBS::SIMPR const scalarField& y, const scalarField& dydx, const scalarField& dfdx, - const Matrix& dfdy, + const scalarSquareMatrix& dfdy, const scalar deltaX, const label nSteps, scalarField& yEnd @@ -44,7 +43,7 @@ void Foam::SIBS::SIMPR { scalar h = deltaX/nSteps; - Matrix a(n_, n_); + scalarSquareMatrix a(n_); for (register label i=0; i::LUDecompose(a, pivotIndices); + LUDecompose(a, pivotIndices); for (register label i=0; i::LUBacksubstitute(a, pivotIndices, yEnd); + LUBacksubstitute(a, pivotIndices, yEnd); scalarField del(yEnd); scalarField ytemp(n_); @@ -83,7 +82,7 @@ void Foam::SIBS::SIMPR yEnd[i] = h*yEnd[i] - del[i]; } - simpleMatrix::LUBacksubstitute(a, pivotIndices, yEnd); + LUBacksubstitute(a, pivotIndices, yEnd); for (register label i=0; i::LUBacksubstitute(a, pivotIndices, yEnd); + LUBacksubstitute(a, pivotIndices, yEnd); for (register label i=0; i& d + scalarSquareMatrix& d ) const { label n = yz.size(); diff --git a/src/OpenFOAM/Make/files b/src/OpenFOAM/Make/files index 57dc7f8dc7..cf4e331e25 100644 --- a/src/OpenFOAM/Make/files +++ b/src/OpenFOAM/Make/files @@ -177,9 +177,9 @@ dimensionedTypes/dimensionedTensor/dimensionedTensor.C matrices/solution/solution.C -scalarMatrix = matrices/scalarMatrix -$(scalarMatrix)/scalarMatrix.C -$(scalarMatrix)/SVD/SVD.C +scalarMatrices = matrices/scalarMatrices +$(scalarMatrices)/scalarMatrices.C +$(scalarMatrices)/SVD/SVD.C LUscalarMatrix = matrices/LUscalarMatrix $(LUscalarMatrix)/LUscalarMatrix.C diff --git a/src/OpenFOAM/db/IOstreams/Pstreams/PstreamCommsStruct.C b/src/OpenFOAM/db/IOstreams/Pstreams/PstreamCommsStruct.C index 6b078aa25d..7a1a680643 100644 --- a/src/OpenFOAM/db/IOstreams/Pstreams/PstreamCommsStruct.C +++ b/src/OpenFOAM/db/IOstreams/Pstreams/PstreamCommsStruct.C @@ -22,8 +22,6 @@ License along with OpenFOAM; if not, write to the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA -Description - \*---------------------------------------------------------------------------*/ #include "Pstream.H" diff --git a/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C b/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C index c023e59c5a..a31058eace 100644 --- a/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C +++ b/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C @@ -29,7 +29,8 @@ License // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // template -Foam::DiagonalMatrix::DiagonalMatrix(const Matrix& a) +template +Foam::DiagonalMatrix::DiagonalMatrix(const Matrix& a) : List(min(a.n(), a.m())) { diff --git a/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.H b/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.H index 36138fc25d..ae8ceffc1c 100644 --- a/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.H +++ b/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.H @@ -46,7 +46,7 @@ namespace Foam // * * * * * * * * * * * * Class Forward declaration * * * * * * * * * * * // -template class Matrix; +template class Matrix; /*---------------------------------------------------------------------------*\ Class DiagonalMatrix Declaration @@ -62,7 +62,8 @@ public: // Constructors //- Construct from diagonal component of a Matrix - DiagonalMatrix(const Matrix&); + template + DiagonalMatrix(const Matrix&); //- Construct empty from size DiagonalMatrix(const label size); diff --git a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C index 5de8693f2f..f425accf72 100644 --- a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C +++ b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C @@ -31,9 +31,9 @@ License // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -Foam::LUscalarMatrix::LUscalarMatrix(const Matrix& matrix) +Foam::LUscalarMatrix::LUscalarMatrix(const scalarSquareMatrix& matrix) : - scalarMatrix(matrix), + scalarSquareMatrix(matrix), pivotIndices_(n()) { LUDecompose(*this, pivotIndices_); @@ -101,7 +101,7 @@ Foam::LUscalarMatrix::LUscalarMatrix nCells += lduMatrices[i].size(); } - Matrix m(nCells, nCells, 0.0); + scalarSquareMatrix m(nCells, 0.0); transfer(m); convert(lduMatrices); } @@ -109,7 +109,7 @@ Foam::LUscalarMatrix::LUscalarMatrix else { label nCells = ldum.lduAddr().size(); - Matrix m(nCells, nCells, 0.0); + scalarSquareMatrix m(nCells, 0.0); transfer(m); convert(ldum, interfaceCoeffs, interfaces); } diff --git a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H index d5e79faf19..2255cee236 100644 --- a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H +++ b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H @@ -36,7 +36,7 @@ SourceFiles #ifndef LUscalarMatrix_H #define LUscalarMatrix_H -#include "scalarMatrix.H" +#include "scalarMatrices.H" #include "labelList.H" #include "FieldField.H" #include "lduInterfaceFieldPtrsList.H" @@ -55,7 +55,7 @@ class procLduMatrix; class LUscalarMatrix : - public scalarMatrix + public scalarSquareMatrix { // Private data @@ -78,7 +78,7 @@ class LUscalarMatrix void convert(const PtrList& lduMatrices); - //- Print the ratio of the mag-sum of the off-diagonal coefficients + //- Print the ratio of the mag-sum of the off-diagonal coefficients // to the mag-diagonal void printDiagonalDominance() const; @@ -87,8 +87,8 @@ public: // Constructors - //- Construct from Matrix and perform LU decomposition - LUscalarMatrix(const Matrix&); + //- Construct from scalarSquareMatrix and perform LU decomposition + LUscalarMatrix(const scalarSquareMatrix&); //- Construct from lduMatrix and perform LU decomposition LUscalarMatrix diff --git a/src/OpenFOAM/matrices/Matrix/Matrix.C b/src/OpenFOAM/matrices/Matrix/Matrix.C index 4fbef71c2b..577cbcabcb 100644 --- a/src/OpenFOAM/matrices/Matrix/Matrix.C +++ b/src/OpenFOAM/matrices/Matrix/Matrix.C @@ -28,8 +28,8 @@ License // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // -template -void Foam::Matrix::allocate() +template +void Foam::Matrix::allocate() { if (n_ && m_) { @@ -46,8 +46,8 @@ void Foam::Matrix::allocate() // * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * * // -template -Foam::Matrix::~Matrix() +template +Foam::Matrix::~Matrix() { if (v_) { @@ -59,16 +59,16 @@ Foam::Matrix::~Matrix() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -template -const Foam::Matrix& Foam::Matrix::null() +template +const Foam::Matrix& Foam::Matrix::null() { - Matrix* nullPtr = reinterpret_cast*>(NULL); + Matrix* nullPtr = reinterpret_cast*>(NULL); return *nullPtr; } -template -Foam::Matrix::Matrix(const label n, const label m) +template +Foam::Matrix::Matrix(const label n, const label m) : n_(n), m_(m), @@ -76,7 +76,7 @@ Foam::Matrix::Matrix(const label n, const label m) { if (n_ < 0 || m_ < 0) { - FatalErrorIn("Matrix::Matrix(const label n, const label m)") + FatalErrorIn("Matrix::Matrix(const label n, const label m)") << "bad n, m " << n_ << ", " << m_ << abort(FatalError); } @@ -85,8 +85,8 @@ Foam::Matrix::Matrix(const label n, const label m) } -template -Foam::Matrix::Matrix(const label n, const label m, const Type& a) +template +Foam::Matrix::Matrix(const label n, const label m, const Type& a) : n_(n), m_(m), @@ -96,7 +96,7 @@ Foam::Matrix::Matrix(const label n, const label m, const Type& a) { FatalErrorIn ( - "Matrix::Matrix(const label n, const label m, const T&)" + "Matrix::Matrix(const label n, const label m, const T&)" ) << "bad n, m " << n_ << ", " << m_ << abort(FatalError); } @@ -117,8 +117,8 @@ Foam::Matrix::Matrix(const label n, const label m, const Type& a) } -template -Foam::Matrix::Matrix(const Matrix& a) +template +Foam::Matrix::Matrix(const Matrix& a) : n_(a.n_), m_(a.m_), @@ -139,8 +139,8 @@ Foam::Matrix::Matrix(const Matrix& a) } -template -void Foam::Matrix::clear() +template +void Foam::Matrix::clear() { if (v_) { @@ -153,8 +153,8 @@ void Foam::Matrix::clear() } -template -void Foam::Matrix::transfer(Matrix& a) +template +void Foam::Matrix::transfer(Matrix& a) { clear(); @@ -169,13 +169,11 @@ void Foam::Matrix::transfer(Matrix& a) } -// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // - -template -Foam::Matrix Foam::Matrix::T() const +template +Form Foam::Matrix::T() const { - const Matrix& A = *this; - Matrix At(m(), n()); + const Matrix& A = *this; + Form At(m(), n()); for (register label i=0; i Foam::Matrix::T() const } -template -void Foam::Matrix::operator=(const Type& t) +// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // + +template +void Foam::Matrix::operator=(const Type& t) { if (v_) { @@ -206,12 +206,12 @@ void Foam::Matrix::operator=(const Type& t) // Assignment operator. Takes linear time. -template -void Foam::Matrix::operator=(const Matrix& a) +template +void Foam::Matrix::operator=(const Matrix& a) { if (this == &a) { - FatalErrorIn("Matrix::operator=(const Matrix&)") + FatalErrorIn("Matrix::operator=(const Matrix&)") << "attempted assignment to self" << abort(FatalError); } @@ -240,8 +240,8 @@ void Foam::Matrix::operator=(const Matrix& a) // * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * // -template -const Type& Foam::max(const Matrix& a) +template +const Type& Foam::max(const Matrix& a) { label nm = a.n_*a.m_; @@ -262,7 +262,7 @@ const Type& Foam::max(const Matrix& a) } else { - FatalErrorIn("max(const Matrix&)") + FatalErrorIn("max(const Matrix&)") << "matrix is empty" << abort(FatalError); @@ -272,8 +272,8 @@ const Type& Foam::max(const Matrix& a) } -template -const Type& Foam::min(const Matrix& a) +template +const Type& Foam::min(const Matrix& a) { label nm = a.n_*a.m_; @@ -294,7 +294,7 @@ const Type& Foam::min(const Matrix& a) } else { - FatalErrorIn("min(const Matrix&)") + FatalErrorIn("min(const Matrix&)") << "matrix is empty" << abort(FatalError); @@ -306,10 +306,10 @@ const Type& Foam::min(const Matrix& a) // * * * * * * * * * * * * * * * Global Operators * * * * * * * * * * * * * // -template -Foam::Matrix Foam::operator-(const Matrix& a) +template +Form Foam::operator-(const Matrix& a) { - Matrix na(a.n_, a.m_); + Form na(a.n_, a.m_); if (a.n_ && a.m_) { @@ -327,14 +327,14 @@ Foam::Matrix Foam::operator-(const Matrix& a) } -template -Foam::Matrix Foam::operator+(const Matrix& a, const Matrix& b) +template +Form Foam::operator+(const Matrix& a, const Matrix& b) { if (a.n_ != b.n_) { FatalErrorIn ( - "Matrix::operator+(const Matrix&, const Matrix&)" + "Matrix::operator+(const Matrix&, const Matrix&)" ) << "attempted add matrices with different number of rows: " << a.n_ << ", " << b.n_ << abort(FatalError); @@ -344,13 +344,13 @@ Foam::Matrix Foam::operator+(const Matrix& a, const Matrix& b) { FatalErrorIn ( - "Matrix::operator+(const Matrix&, const Matrix&)" + "Matrix::operator+(const Matrix&, const Matrix&)" ) << "attempted add matrices with different number of columns: " << a.m_ << ", " << b.m_ << abort(FatalError); } - Matrix ab(a.n_, a.m_); + Form ab(a.n_, a.m_); Type* abv = ab.v_[0]; const Type* av = a.v_[0]; @@ -366,14 +366,14 @@ Foam::Matrix Foam::operator+(const Matrix& a, const Matrix& b) } -template -Foam::Matrix Foam::operator-(const Matrix& a, const Matrix& b) +template +Form Foam::operator-(const Matrix& a, const Matrix& b) { if (a.n_ != b.n_) { FatalErrorIn ( - "Matrix::operator-(const Matrix&, const Matrix&)" + "Matrix::operator-(const Matrix&, const Matrix&)" ) << "attempted add matrices with different number of rows: " << a.n_ << ", " << b.n_ << abort(FatalError); @@ -383,13 +383,13 @@ Foam::Matrix Foam::operator-(const Matrix& a, const Matrix& b) { FatalErrorIn ( - "Matrix::operator-(const Matrix&, const Matrix&)" + "Matrix::operator-(const Matrix&, const Matrix&)" ) << "attempted add matrices with different number of columns: " << a.m_ << ", " << b.m_ << abort(FatalError); } - Matrix ab(a.n_, a.m_); + Form ab(a.n_, a.m_); Type* abv = ab.v_[0]; const Type* av = a.v_[0]; @@ -405,17 +405,17 @@ Foam::Matrix Foam::operator-(const Matrix& a, const Matrix& b) } -template -Foam::Matrix Foam::operator*(const scalar s, const Matrix& a) +template +Form Foam::operator*(const scalar s, const Matrix& a) { - Matrix sa(a.n_, a.m_); + Form sa(a.n_, a.m_); if (a.n_ && a.m_) { Type* sav = sa.v_[0]; const Type* av = a.v_[0]; - label nm = a.n_*a.m_;; + label nm = a.n_*a.m_; for (register label i=0; i class Matrix; +template class Matrix; -template const Type& max(const Matrix&); -template const Type& min(const Matrix&); - -template Matrix operator-(const Matrix&); -template Matrix operator+ +template Istream& operator>> ( - const Matrix&, - const Matrix& -); -template Matrix operator- -( - const Matrix&, - const Matrix& -); -template Matrix operator* -( - const scalar, - const Matrix& + Istream&, + Matrix& ); -template Istream& operator>>(Istream&, Matrix&); -template Ostream& operator<<(Ostream&, const Matrix&); +template Ostream& operator<< +( + Ostream&, + const Matrix& +); /*---------------------------------------------------------------------------*\ Class Matrix Declaration \*---------------------------------------------------------------------------*/ -template +template class Matrix { // Private data @@ -112,13 +101,13 @@ public: Matrix(const label n, const label m, const Type&); //- Copy constructor. - Matrix(const Matrix&); + Matrix(const Matrix&); //- Construct from Istream. Matrix(Istream&); //- Clone - inline autoPtr > clone() const; + inline autoPtr > clone() const; // Destructor @@ -129,7 +118,7 @@ public: // Member functions //- Return a null Matrix - static const Matrix& null(); + static const Matrix& null(); // Access @@ -160,7 +149,11 @@ public: //- Transfer the contents of the argument Matrix into this Matrix // and annull the argument Matrix. - void transfer(Matrix&); + void transfer(Matrix&); + + + //- Return the transpose of the matrix + Form T() const; // Member operators @@ -171,52 +164,49 @@ public: //- Return subscript-checked element of constant Matrix. inline const Type* operator[](const label) const; - //- Return the transpose of the matrix - Matrix T() const; - //- Assignment operator. Takes linear time. - void operator=(const Matrix&); + void operator=(const Matrix&); //- Assignment of all entries to the given value void operator=(const Type&); - // Friend functions - - friend const Type& max(const Matrix&); - friend const Type& min(const Matrix&); - - - // Friend operators - - friend Matrix operator-(const Matrix&); - friend Matrix operator+ - ( - const Matrix&, - const Matrix& - ); - friend Matrix operator- - ( - const Matrix&, - const Matrix& - ); - friend Matrix operator* - ( - const scalar, - const Matrix& - ); - - // IOstream operators //- Read Matrix from Istream, discarding contents of existing Matrix. - friend Istream& operator>> (Istream&, Matrix&); + friend Istream& operator>> (Istream&, Matrix&); // Write Matrix to Ostream. - friend Ostream& operator<< (Ostream&, const Matrix&); + friend Ostream& operator<< (Ostream&, const Matrix&); }; +// Global functions and operators + +template const Type& max(const Matrix&); +template const Type& min(const Matrix&); + +template Form operator-(const Matrix&); + +template Form operator+ +( + const Matrix&, + const Matrix& +); + +template Form operator- +( + const Matrix&, + const Matrix& +); + +template Form operator* +( + const scalar, + const Matrix& +); + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace Foam diff --git a/src/OpenFOAM/matrices/Matrix/MatrixI.H b/src/OpenFOAM/matrices/Matrix/MatrixI.H index 8323174d4f..8d758cf5b1 100644 --- a/src/OpenFOAM/matrices/Matrix/MatrixI.H +++ b/src/OpenFOAM/matrices/Matrix/MatrixI.H @@ -26,8 +26,8 @@ License // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -template -inline Foam::Matrix::Matrix() +template +inline Foam::Matrix::Matrix() : n_(0), m_(0), @@ -35,67 +35,67 @@ inline Foam::Matrix::Matrix() {} -template -inline Foam::autoPtr > Foam::Matrix::clone() const +template +inline Foam::autoPtr > Foam::Matrix::clone() const { - return autoPtr >(new Matrix(*this)); + return autoPtr >(new Matrix(*this)); } // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // //- Return the number of rows -template -inline Foam::label Foam::Matrix::n() const +template +inline Foam::label Foam::Matrix::n() const { return n_; } -template -inline Foam::label Foam::Matrix::m() const +template +inline Foam::label Foam::Matrix::m() const { return m_; } -template -inline Foam::label Foam::Matrix::size() const +template +inline Foam::label Foam::Matrix::size() const { return n_*m_; } -template -inline void Foam::Matrix::checki(const label i) const +template +inline void Foam::Matrix::checki(const label i) const { if (!n_) { - FatalErrorIn("Matrix::checki(const label)") + FatalErrorIn("Matrix::checki(const label)") << "attempt to access element from zero sized row" << abort(FatalError); } else if (i<0 || i>=n_) { - FatalErrorIn("Matrix::checki(const label)") + FatalErrorIn("Matrix::checki(const label)") << "index " << i << " out of range 0 ... " << n_-1 << abort(FatalError); } } -template -inline void Foam::Matrix::checkj(const label j) const +template +inline void Foam::Matrix::checkj(const label j) const { if (!m_) { - FatalErrorIn("Matrix::checkj(const label)") + FatalErrorIn("Matrix::checkj(const label)") << "attempt to access element from zero sized column" << abort(FatalError); } else if (j<0 || j>=m_) { - FatalErrorIn("Matrix::checkj(const label)") + FatalErrorIn("Matrix::checkj(const label)") << "index " << j << " out of range 0 ... " << m_-1 << abort(FatalError); } @@ -104,8 +104,8 @@ inline void Foam::Matrix::checkj(const label j) const // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // -template -inline Type* Foam::Matrix::operator[](const label i) +template +inline Type* Foam::Matrix::operator[](const label i) { # ifdef FULLDEBUG checki(i); @@ -114,8 +114,8 @@ inline Type* Foam::Matrix::operator[](const label i) } -template -inline const Type* Foam::Matrix::operator[](const label i) const +template +inline const Type* Foam::Matrix::operator[](const label i) const { # ifdef FULLDEBUG checki(i); diff --git a/src/OpenFOAM/matrices/Matrix/MatrixIO.C b/src/OpenFOAM/matrices/Matrix/MatrixIO.C index 9b5a922f0d..15f0d5345e 100644 --- a/src/OpenFOAM/matrices/Matrix/MatrixIO.C +++ b/src/OpenFOAM/matrices/Matrix/MatrixIO.C @@ -32,8 +32,8 @@ License // * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * // -template -Foam::Matrix::Matrix(Istream& is) +template +Foam::Matrix::Matrix(Istream& is) : n_(0), m_(0), @@ -43,17 +43,17 @@ Foam::Matrix::Matrix(Istream& is) } -template -Foam::Istream& Foam::operator>>(Istream& is, Matrix& M) +template +Foam::Istream& Foam::operator>>(Istream& is, Matrix& M) { // Anull matrix M.clear(); - is.fatalCheck("operator>>(Istream&, Matrix&)"); + is.fatalCheck("operator>>(Istream&, Matrix&)"); token firstToken(is); - is.fatalCheck("operator>>(Istream&, Matrix&) : reading first token"); + is.fatalCheck("operator>>(Istream&, Matrix&) : reading first token"); if (firstToken.isLabel()) { @@ -88,7 +88,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix& M) is.fatalCheck ( - "operator>>(Istream&, Matrix&) : " + "operator>>(Istream&, Matrix&) : " "reading entry" ); } @@ -103,7 +103,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix& M) is.fatalCheck ( - "operator>>(Istream&, Matrix&) : " + "operator>>(Istream&, Matrix&) : " "reading the single entry" ); @@ -128,7 +128,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix& M) is.fatalCheck ( - "operator>>(Istream&, Matrix&) : " + "operator>>(Istream&, Matrix&) : " "reading the binary block" ); } @@ -136,7 +136,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix& M) } else { - FatalIOErrorIn("operator>>(Istream&, Matrix&)", is) + FatalIOErrorIn("operator>>(Istream&, Matrix&)", is) << "incorrect first token, expected , found " << firstToken.info() << exit(FatalIOError); @@ -146,8 +146,8 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix& M) } -template -Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix& M) +template +Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix& M) { label nm = M.n_*M.m_; diff --git a/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H b/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H new file mode 100644 index 0000000000..e5e6b16777 --- /dev/null +++ b/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H @@ -0,0 +1,92 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software; you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation; either version 2 of the License, or (at your + option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Class + Foam::RectangularMatrix + +Description + A templated 2D rectangular matrix of objects of \, where the n x n matrix + dimension is known and used for subscript bounds checking, etc. + +SourceFiles + RectangularMatrixI.H + RectangularMatrix.C + +\*---------------------------------------------------------------------------*/ + +#ifndef RectangularMatrix_H +#define RectangularMatrix_H + +#include "Matrix.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class Matrix Declaration +\*---------------------------------------------------------------------------*/ + +template +class RectangularMatrix +: + public Matrix, Type> +{ + +public: + + // Constructors + + //- Null constructor. + inline RectangularMatrix(); + + //- Construct given number of rows and columns, + inline RectangularMatrix(const label m, const label n); + + //- Construct with given number of rows and columns + // and value for all elements. + inline RectangularMatrix(const label m, const label n, const Type&); + + //- Construct from Istream. + inline RectangularMatrix(Istream&); + + //- Clone + inline autoPtr > clone() const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +# include "RectangularMatrixI.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H b/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H new file mode 100644 index 0000000000..1cacd2f11a --- /dev/null +++ b/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H @@ -0,0 +1,70 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software; you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation; either version 2 of the License, or (at your + option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +inline Foam::RectangularMatrix::RectangularMatrix() +: + Matrix, Type>() +{} + +template +inline Foam::RectangularMatrix::RectangularMatrix +( + const label m, + const label n +) +: + Matrix, Type>(m, n) +{} + +template +inline Foam::RectangularMatrix::RectangularMatrix +( + const label m, + const label n, + const Type& t +) +: + Matrix, Type>(m, n, t) +{} + +template +inline Foam::RectangularMatrix::RectangularMatrix(Istream& is) +: + Matrix, Type>(is) +{} + +template +inline Foam::autoPtr > +Foam::RectangularMatrix::clone() const +{ + return autoPtr >(new RectangularMatrix(*this)); +} + + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H new file mode 100644 index 0000000000..6239781d85 --- /dev/null +++ b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H @@ -0,0 +1,96 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software; you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation; either version 2 of the License, or (at your + option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Class + Foam::SquareMatrix + +Description + A templated 2D square matrix of objects of \, where the n x n matrix + dimension is known and used for subscript bounds checking, etc. + +SourceFiles + SquareMatrixI.H + SquareMatrix.C + +\*---------------------------------------------------------------------------*/ + +#ifndef SquareMatrix_H +#define SquareMatrix_H + +#include "Matrix.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class Matrix Declaration +\*---------------------------------------------------------------------------*/ + +template +class SquareMatrix +: + public Matrix, Type> +{ + +public: + + // Constructors + + //- Null constructor. + inline SquareMatrix(); + + //- Construct given number of rows/columns. + inline SquareMatrix(const label n); + + //- Construct given number of rows and columns, + // It checks that m == n. + inline SquareMatrix(const label m, const label n); + + //- Construct with given number of rows/columns + // and value for all elements. + inline SquareMatrix(const label n, const Type&); + + //- Construct from Istream. + inline SquareMatrix(Istream&); + + //- Clone + inline autoPtr > clone() const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +# include "SquareMatrixI.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H new file mode 100644 index 0000000000..a7371a6a2e --- /dev/null +++ b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H @@ -0,0 +1,75 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software; you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation; either version 2 of the License, or (at your + option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +inline Foam::SquareMatrix::SquareMatrix() +: + Matrix, Type>() +{} + +template +inline Foam::SquareMatrix::SquareMatrix(const label n) +: + Matrix, Type>(n, n) +{} + +template +inline Foam::SquareMatrix::SquareMatrix(const label m, const label n) +: + Matrix, Type>(m, n) +{ + if (m != n) + { + FatalErrorIn + ( + "SquareMatrix::SquareMatrix(const label m, const label n)" + ) << "m != n for constructing a square matrix" << exit(FatalError); + } +} + +template +inline Foam::SquareMatrix::SquareMatrix(const label n, const Type& t) +: + Matrix, Type>(n, t) +{} + +template +inline Foam::SquareMatrix::SquareMatrix(Istream& is) +: + Matrix, Type>(is) +{} + +template +inline Foam::autoPtr > +Foam::SquareMatrix::clone() const +{ + return autoPtr >(new SquareMatrix(*this)); +} + + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/scalarMatrix/SVD/SVD.C b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.C similarity index 96% rename from src/OpenFOAM/matrices/scalarMatrix/SVD/SVD.C rename to src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.C index a1355c393d..13ef29ed9b 100644 --- a/src/OpenFOAM/matrices/scalarMatrix/SVD/SVD.C +++ b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.C @@ -26,12 +26,12 @@ License #include "SVD.H" #include "scalarList.H" -#include "scalarMatrix.H" +#include "scalarMatrices.H" #include "ListOps.H" // * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * // -Foam::SVD::SVD(const Matrix& A, const scalar minCondition) +Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition) : U_(A), V_(A.m(), A.m()), @@ -294,8 +294,11 @@ Foam::SVD::SVD(const Matrix& A, const scalar minCondition) } if (its == 34) { - WarningIn("SVD::SVD(Matrix& A)") - << "no convergence in 35 SVD iterations" + WarningIn + ( + "SVD::SVD" + "(scalarRectangularMatrix& A, const scalar minCondition)" + ) << "no convergence in 35 SVD iterations" << endl; } @@ -375,7 +378,7 @@ Foam::SVD::SVD(const Matrix& A, const scalar minCondition) multiply(VSinvUt_, V_, inv(S_), U_.T()); // test SVD - /*Matrix SVDA(A.n(), A.m()); + /*scalarRectangularMatrix SVDA(A.n(), A.m()); multiply(SVDA, U_, S_, transpose(V_)); scalar maxDiff = 0; scalar diff = 0; diff --git a/src/OpenFOAM/matrices/scalarMatrix/SVD/SVD.H b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.H similarity index 86% rename from src/OpenFOAM/matrices/scalarMatrix/SVD/SVD.H rename to src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.H index b597027732..48e581198a 100644 --- a/src/OpenFOAM/matrices/scalarMatrix/SVD/SVD.H +++ b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.H @@ -23,12 +23,13 @@ License Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA Class - SVD + Foam::SVD Description Singular value decomposition of a rectangular matrix. SourceFiles + SVDI.H SVD.C \*---------------------------------------------------------------------------*/ @@ -36,8 +37,7 @@ SourceFiles #ifndef SVD_H #define SVD_H -#include "DiagonalMatrix.H" -#include "Matrix.H" +#include "scalarMatrices.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -56,16 +56,16 @@ class SVD // Private data //- Rectangular matrix with the same dimensions as the input - Matrix U_; + scalarRectangularMatrix U_; //- square matrix V - Matrix V_; + scalarRectangularMatrix V_; //- The singular values DiagonalMatrix S_; //- The matrix product V S^(-1) U^T - Matrix VSinvUt_; + scalarRectangularMatrix VSinvUt_; //- The number of zero singular values label nZeros_; @@ -88,22 +88,22 @@ public: // Constructors //- Construct from a rectangular Matrix - SVD(const Matrix& A, const scalar minCondition = 0); + SVD(const scalarRectangularMatrix& A, const scalar minCondition = 0); // Access functions //- Return U - inline const Matrix& U() const; + inline const scalarRectangularMatrix& U() const; //- Return the square matrix V - inline const Matrix& V() const; + inline const scalarRectangularMatrix& V() const; //- Return the singular values - inline const DiagonalMatrix& S() const; + inline const scalarDiagonalMatrix& S() const; //- Return VSinvUt (the pseudo inverse) - inline const Matrix& VSinvUt() const; + inline const scalarRectangularMatrix& VSinvUt() const; //- Return the number of zero singular values inline label nZeros() const; diff --git a/src/OpenFOAM/matrices/scalarMatrix/SVD/SVDI.H b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVDI.H similarity index 88% rename from src/OpenFOAM/matrices/scalarMatrix/SVD/SVDI.H rename to src/OpenFOAM/matrices/scalarMatrices/SVD/SVDI.H index 337647edfb..298aac1aeb 100644 --- a/src/OpenFOAM/matrices/scalarMatrix/SVD/SVDI.H +++ b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVDI.H @@ -35,22 +35,22 @@ inline const T Foam::SVD::sign(const T& a, const T& b) // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -inline const Foam::Matrix& Foam::SVD::U() const +inline const Foam::scalarRectangularMatrix& Foam::SVD::U() const { return U_; } -inline const Foam::Matrix& Foam::SVD::V() const +inline const Foam::scalarRectangularMatrix& Foam::SVD::V() const { return V_; } -inline const Foam::DiagonalMatrix& Foam::SVD::S() const +inline const Foam::scalarDiagonalMatrix& Foam::SVD::S() const { return S_; } -inline const Foam::Matrix& Foam::SVD::VSinvUt() const +inline const Foam::scalarRectangularMatrix& Foam::SVD::VSinvUt() const { return VSinvUt_; } diff --git a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.C b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C similarity index 76% rename from src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.C rename to src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C index b88e167b6f..8ddd01a82c 100644 --- a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.C +++ b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C @@ -24,39 +24,14 @@ License \*---------------------------------------------------------------------------*/ -#include "scalarMatrix.H" +#include "scalarMatrices.H" #include "SVD.H" - -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -Foam::scalarMatrix::scalarMatrix() -{} - - -Foam::scalarMatrix::scalarMatrix(const label mSize) -: - Matrix(mSize, mSize, 0.0) -{} - - -Foam::scalarMatrix::scalarMatrix(const Matrix& matrix) -: - Matrix(matrix) -{} - - -Foam::scalarMatrix::scalarMatrix(Istream& is) -: - Matrix(is) -{} - - // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -void Foam::scalarMatrix::LUDecompose +void Foam::LUDecompose ( - Matrix& matrix, + scalarSquareMatrix& matrix, labelList& pivotIndices ) { @@ -81,8 +56,8 @@ void Foam::scalarMatrix::LUDecompose { FatalErrorIn ( - "scalarMatrix::LUdecompose" - "(Matrix& matrix, labelList& rowIndices)" + "LUdecompose" + "(scalarSquareMatrix& matrix, labelList& rowIndices)" ) << "Singular matrix" << exit(FatalError); } @@ -164,9 +139,9 @@ void Foam::scalarMatrix::LUDecompose void Foam::multiply ( - Matrix& ans, // value changed in return - const Matrix& A, - const Matrix& B + scalarRectangularMatrix& ans, // value changed in return + const scalarRectangularMatrix& A, + const scalarRectangularMatrix& B ) { if (A.m() != B.n()) @@ -174,15 +149,15 @@ void Foam::multiply FatalErrorIn ( "multiply(" - "Matrix& answer " - "const Matrix& A, " - "const Matrix& B)" + "scalarRectangularMatrix& answer " + "const scalarRectangularMatrix& A, " + "const scalarRectangularMatrix& B)" ) << "A and B must have identical inner dimensions but A.m = " << A.m() << " and B.n = " << B.n() << abort(FatalError); } - ans = Matrix(A.n(), B.m(), scalar(0)); + ans = scalarRectangularMatrix(A.n(), B.m(), scalar(0)); for(register label i = 0; i < A.n(); i++) { @@ -199,10 +174,10 @@ void Foam::multiply void Foam::multiply ( - Matrix& ans, // value changed in return - const Matrix& A, - const Matrix& B, - const Matrix& C + scalarRectangularMatrix& ans, // value changed in return + const scalarRectangularMatrix& A, + const scalarRectangularMatrix& B, + const scalarRectangularMatrix& C ) { if (A.m() != B.n()) @@ -210,10 +185,10 @@ void Foam::multiply FatalErrorIn ( "multiply(" - "const Matrix& A, " - "const Matrix& B, " - "const Matrix& C, " - "Matrix& answer)" + "const scalarRectangularMatrix& A, " + "const scalarRectangularMatrix& B, " + "const scalarRectangularMatrix& C, " + "scalarRectangularMatrix& answer)" ) << "A and B must have identical inner dimensions but A.m = " << A.m() << " and B.n = " << B.n() << abort(FatalError); @@ -224,16 +199,16 @@ void Foam::multiply FatalErrorIn ( "multiply(" - "const Matrix& A, " - "const Matrix& B, " - "const Matrix& C, " - "Matrix& answer)" + "const scalarRectangularMatrix& A, " + "const scalarRectangularMatrix& B, " + "const scalarRectangularMatrix& C, " + "scalarRectangularMatrix& answer)" ) << "B and C must have identical inner dimensions but B.m = " << B.m() << " and C.n = " << C.n() << abort(FatalError); } - ans = Matrix(A.n(), C.m(), scalar(0)); + ans = scalarRectangularMatrix(A.n(), C.m(), scalar(0)); for(register label i = 0; i < A.n(); i++) { @@ -255,10 +230,10 @@ void Foam::multiply void Foam::multiply ( - Matrix& ans, // value changed in return - const Matrix& A, + scalarRectangularMatrix& ans, // value changed in return + const scalarRectangularMatrix& A, const DiagonalMatrix& B, - const Matrix& C + const scalarRectangularMatrix& C ) { if (A.m() != B.size()) @@ -266,10 +241,10 @@ void Foam::multiply FatalErrorIn ( "multiply(" - "const Matrix& A, " + "const scalarRectangularMatrix& A, " "const DiagonalMatrix& B, " - "const Matrix& C, " - "Matrix& answer)" + "const scalarRectangularMatrix& C, " + "scalarRectangularMatrix& answer)" ) << "A and B must have identical inner dimensions but A.m = " << A.m() << " and B.n = " << B.size() << abort(FatalError); @@ -280,16 +255,16 @@ void Foam::multiply FatalErrorIn ( "multiply(" - "const Matrix& A, " + "const scalarRectangularMatrix& A, " "const DiagonalMatrix& B, " - "const Matrix& C, " - "Matrix& answer)" + "const scalarRectangularMatrix& C, " + "scalarRectangularMatrix& answer)" ) << "B and C must have identical inner dimensions but B.m = " << B.size() << " and C.n = " << C.n() << abort(FatalError); } - ans = Matrix(A.n(), C.m(), scalar(0)); + ans = scalarRectangularMatrix(A.n(), C.m(), scalar(0)); for(register label i = 0; i < A.n(); i++) { @@ -304,9 +279,9 @@ void Foam::multiply } -Foam::Matrix Foam::SVDinv +Foam::RectangularMatrix Foam::SVDinv ( - const Matrix& A, + const scalarRectangularMatrix& A, scalar minCondition ) { diff --git a/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.H b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.H new file mode 100644 index 0000000000..17d8727516 --- /dev/null +++ b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.H @@ -0,0 +1,137 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software; you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation; either version 2 of the License, or (at your + option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Class + scalarMatrices + +Description + Scalar matrices + +SourceFiles + scalarMatrices.C + scalarMatricesTemplates.C + +\*---------------------------------------------------------------------------*/ + +#ifndef scalarMatrices_H +#define scalarMatrices_H + +#include "RectangularMatrix.H" +#include "SquareMatrix.H" +#include "DiagonalMatrix.H" +#include "scalarField.H" +#include "labelList.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +typedef RectangularMatrix scalarRectangularMatrix; +typedef SquareMatrix scalarSquareMatrix; +typedef DiagonalMatrix scalarDiagonalMatrix; + +//- Solve the matrix using Gaussian elimination with pivoting, +// returning the solution in the source +template +void solve(scalarSquareMatrix& matrix, Field& source); + +//- Solve the matrix using Gaussian elimination with pivoting +// and return the solution +template +void solve +( + Field& psi, + const scalarSquareMatrix& matrix, + const Field& source +); + +//- LU decompose the matrix with pivoting +void LUDecompose +( + scalarSquareMatrix& matrix, + labelList& pivotIndices +); + +//- LU back-substitution with given source, returning the solution +// in the source +template +void LUBacksubstitute +( + const scalarSquareMatrix& luMmatrix, + const labelList& pivotIndices, + Field& source +); + +//- Solve the matrix using LU decomposition with pivoting +// returning the LU form of the matrix and the solution in the source +template +void LUsolve(scalarSquareMatrix& matrix, Field& source); + +void multiply +( + scalarRectangularMatrix& answer, // value changed in return + const scalarRectangularMatrix& A, + const scalarRectangularMatrix& B +); + +void multiply +( + scalarRectangularMatrix& answer, // value changed in return + const scalarRectangularMatrix& A, + const scalarRectangularMatrix& B, + const scalarRectangularMatrix& C +); + +void multiply +( + scalarRectangularMatrix& answer, // value changed in return + const scalarRectangularMatrix& A, + const DiagonalMatrix& B, + const scalarRectangularMatrix& C +); + +//- Return the inverse of matrix A using SVD +scalarRectangularMatrix SVDinv +( + const scalarRectangularMatrix& A, + scalar minCondition = 0 +); + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "scalarMatricesTemplates.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrixTemplates.C b/src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C similarity index 90% rename from src/OpenFOAM/matrices/scalarMatrix/scalarMatrixTemplates.C rename to src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C index fb48a9b26d..131856b3fc 100644 --- a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrixTemplates.C +++ b/src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C @@ -24,15 +24,15 @@ License \*---------------------------------------------------------------------------*/ -#include "scalarMatrix.H" +#include "scalarMatrices.H" #include "Swap.H" // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // template -void Foam::scalarMatrix::solve +void Foam::solve ( - Matrix& tmpMatrix, + scalarSquareMatrix& tmpMatrix, Field& sourceSol ) { @@ -68,7 +68,7 @@ void Foam::scalarMatrix::solve // Check that the system of equations isn't singular if (mag(tmpMatrix[i][i]) < 1e-20) { - FatalErrorIn("scalarMatrix::solve()") + FatalErrorIn("solve(scalarSquareMatrix&, Field& sourceSol)") << "Singular Matrix" << exit(FatalError); } @@ -102,18 +102,23 @@ void Foam::scalarMatrix::solve template -void Foam::scalarMatrix::solve(Field& psi, const Field& source) const +void Foam::solve +( + Field& psi, + const scalarSquareMatrix& matrix, + const Field& source +) { - Matrix tmpMatrix = *this; + scalarSquareMatrix tmpMatrix = matrix; psi = source; solve(tmpMatrix, psi); } template -void Foam::scalarMatrix::LUBacksubstitute +void Foam::LUBacksubstitute ( - const Matrix& luMatrix, + const scalarSquareMatrix& luMatrix, const labelList& pivotIndices, Field& sourceSol ) @@ -160,9 +165,9 @@ void Foam::scalarMatrix::LUBacksubstitute template -void Foam::scalarMatrix::LUsolve +void Foam::LUsolve ( - Matrix& matrix, + scalarSquareMatrix& matrix, Field& sourceSol ) { diff --git a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.H b/src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.H deleted file mode 100644 index 2d22e71d4e..0000000000 --- a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.H +++ /dev/null @@ -1,155 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM is free software; you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation; either version 2 of the License, or (at your - option) any later version. - - OpenFOAM is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License - for more details. - - You should have received a copy of the GNU General Public License - along with OpenFOAM; if not, write to the Free Software Foundation, - Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -Class - Foam::scalarMatrix - -Description - Foam::scalarMatrix - -SourceFiles - scalarMatrix.C - -\*---------------------------------------------------------------------------*/ - -#ifndef scalarMatrix_H -#define scalarMatrix_H - -#include "Matrix.H" -#include "DiagonalMatrix.H" -#include "scalarField.H" -#include "labelList.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -/*---------------------------------------------------------------------------*\ - Class scalarMatrix Declaration -\*---------------------------------------------------------------------------*/ - -class scalarMatrix -: - public Matrix -{ - -public: - - // Constructors - - //- Construct null - scalarMatrix(); - - //- Construct given size - scalarMatrix(const label); - - //- Construct from Matrix - scalarMatrix(const Matrix&); - - //- Construct from Istream - scalarMatrix(Istream&); - - - // Member Functions - - //- Solve the matrix using Gaussian elimination with pivoting, - // returning the solution in the source - template - static void solve(Matrix& matrix, Field& source); - - //- Solve the matrix using Gaussian elimination with pivoting - // and return the solution - template - void solve(Field& psi, const Field& source) const; - - - //- LU decompose the matrix with pivoting - static void LUDecompose - ( - Matrix& matrix, - labelList& pivotIndices - ); - - //- LU back-substitution with given source, returning the solution - // in the source - template - static void LUBacksubstitute - ( - const Matrix& luMmatrix, - const labelList& pivotIndices, - Field& source - ); - - //- Solve the matrix using LU decomposition with pivoting - // returning the LU form of the matrix and the solution in the source - template - static void LUsolve(Matrix& matrix, Field& source); -}; - - -// Global functions - -void multiply -( - Matrix& answer, // value changed in return - const Matrix& A, - const Matrix& B -); - -void multiply -( - Matrix& answer, // value changed in return - const Matrix& A, - const Matrix& B, - const Matrix& C -); - -void multiply -( - Matrix& answer, // value changed in return - const Matrix& A, - const DiagonalMatrix& B, - const Matrix& C -); - -//- Return the inverse of matrix A using SVD -Matrix SVDinv(const Matrix& A, scalar minCondition = 0); - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#ifdef NoRepository -# include "scalarMatrixTemplates.C" -#endif - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#endif - -// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C b/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C index 1265b87ee5..f0f80750a5 100644 --- a/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C +++ b/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C @@ -31,7 +31,7 @@ License template Foam::simpleMatrix::simpleMatrix(const label mSize) : - scalarMatrix(mSize), + scalarSquareMatrix(mSize), source_(mSize, pTraits::zero) {} @@ -39,11 +39,11 @@ Foam::simpleMatrix::simpleMatrix(const label mSize) template Foam::simpleMatrix::simpleMatrix ( - const scalarMatrix& matrix, + const scalarSquareMatrix& matrix, const Field& source ) : - scalarMatrix(matrix), + scalarSquareMatrix(matrix), source_(source) {} @@ -51,7 +51,7 @@ Foam::simpleMatrix::simpleMatrix template Foam::simpleMatrix::simpleMatrix(Istream& is) : - scalarMatrix(is), + scalarSquareMatrix(is), source_(is) {} @@ -61,10 +61,10 @@ Foam::simpleMatrix::simpleMatrix(Istream& is) template Foam::Field Foam::simpleMatrix::solve() const { - scalarMatrix tmpMatrix = *this; + scalarSquareMatrix tmpMatrix = *this; Field sourceSol = source_; - scalarMatrix::solve(tmpMatrix, sourceSol); + Foam::solve(tmpMatrix, sourceSol); return sourceSol; } @@ -73,10 +73,10 @@ Foam::Field Foam::simpleMatrix::solve() const template Foam::Field Foam::simpleMatrix::LUsolve() const { - scalarMatrix luMatrix = *this; + scalarSquareMatrix luMatrix = *this; Field sourceSol = source_; - scalarMatrix::LUsolve(luMatrix, sourceSol); + Foam::LUsolve(luMatrix, sourceSol); return sourceSol; } @@ -108,7 +108,7 @@ void Foam::simpleMatrix::operator=(const simpleMatrix& m) << abort(FatalError); } - scalarMatrix::operator=(m); + scalarSquareMatrix::operator=(m); source_ = m.source_; } @@ -124,8 +124,8 @@ Foam::simpleMatrix Foam::operator+ { return simpleMatrix ( - static_cast(m1) - + static_cast(m2), + static_cast(m1) + + static_cast(m2), m1.source_ + m2.source_ ); } @@ -140,8 +140,8 @@ Foam::simpleMatrix Foam::operator- { return simpleMatrix ( - static_cast(m1) - - static_cast(m2), + static_cast(m1) + - static_cast(m2), m1.source_ - m2.source_ ); } @@ -159,7 +159,7 @@ Foam::simpleMatrix Foam::operator*(const scalar s, const simpleMatrix Foam::Ostream& Foam::operator<<(Ostream& os, const simpleMatrix& m) { - os << static_cast(m) << nl << m.source_; + os << static_cast(m) << nl << m.source_; return os; } diff --git a/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.H b/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.H index d7b1a027be..a972c787c8 100644 --- a/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.H +++ b/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.H @@ -36,7 +36,7 @@ SourceFiles #ifndef simpleMatrix_H #define simpleMatrix_H -#include "scalarMatrix.H" +#include "scalarMatrices.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -48,27 +48,6 @@ namespace Foam template class simpleMatrix; -template -simpleMatrix operator+ -( - const simpleMatrix&, - const simpleMatrix& -); - -template -simpleMatrix operator- -( - const simpleMatrix&, - const simpleMatrix& -); - -template -simpleMatrix operator* -( - const scalar, - const simpleMatrix& -); - template Ostream& operator<< ( @@ -84,7 +63,7 @@ Ostream& operator<< template class simpleMatrix : - public scalarMatrix + public scalarSquareMatrix { // Private data @@ -99,7 +78,7 @@ public: simpleMatrix(const label); //- Construct from components - simpleMatrix(const scalarMatrix&, const Field&); + simpleMatrix(const scalarSquareMatrix&, const Field&); //- Construct from Istream simpleMatrix(Istream&); @@ -137,27 +116,6 @@ public: void operator=(const simpleMatrix&); - // Friend Operators - - friend simpleMatrix operator+ - ( - const simpleMatrix&, - const simpleMatrix& - ); - - friend simpleMatrix operator- - ( - const simpleMatrix&, - const simpleMatrix& - ); - - friend simpleMatrix operator* - ( - const scalar, - const simpleMatrix& - ); - - // Ostream Operator friend Ostream& operator<< @@ -168,6 +126,30 @@ public: }; +// Global operators + +template +simpleMatrix operator+ +( + const simpleMatrix&, + const simpleMatrix& +); + +template +simpleMatrix operator- +( + const simpleMatrix&, + const simpleMatrix& +); + +template +simpleMatrix operator* +( + const scalar, + const simpleMatrix& +); + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace Foam diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C index fd0e55fcdd..0caebcd421 100644 --- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C @@ -49,11 +49,11 @@ Foam::quadraticFitData::quadraticFitData : MeshObject(mesh), centralWeight_(cWeight), -#ifdef SPHERICAL_GEOMETRY +# ifdef SPHERICAL_GEOMETRY dim_(2), -#else +# else dim_(mesh.nGeometricD()), -#endif +# endif minSize_ ( dim_ == 1 ? 3 : @@ -204,7 +204,7 @@ Foam::label Foam::quadraticFitData::calcFit scalar scale = 0; // calculate the matrix of the polynomial components - Matrix B(C.size(), minSize_, scalar(0)); + scalarRectangularMatrix B(C.size(), minSize_, scalar(0)); for(label ip = 0; ip < C.size(); ip++) { @@ -212,11 +212,11 @@ Foam::label Foam::quadraticFitData::calcFit scalar px = (p - p0)&idir; scalar py = (p - p0)&jdir; -#ifndef SPHERICAL_GEOMETRY +# ifndef SPHERICAL_GEOMETRY scalar pz = (p - p0)&kdir; -#else +# else scalar pz = mag(p) - mag(p0); -#endif +# endif if (ip == 0) { diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel.C index 3dd49e126e..f2da51c746 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel.C +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel.C @@ -114,7 +114,7 @@ Foam::scalarField Foam::chemistryModel::omega forAll(reactions_, i) { const reaction& R = reactions_[i]; - + scalar omegai = omega ( R, c, T, p, pf, cf, lRef, pr, cr, rRef @@ -164,13 +164,13 @@ Foam::scalar Foam::chemistryModel::omega pf = 1.0; pr = 1.0; - + label Nl = R.lhs().size(); label Nr = R.rhs().size(); - + label slRef = 0; lRef = R.lhs()[slRef].index; - + pf = kf; for(label s=1; s& dfdc + scalarSquareMatrix& dfdc ) const { scalar T = c[Ns_]; scalar p = c[Ns_ + 1]; - + scalarField c2(Ns(), 0.0); for(label i=0; i Foam::chemistryModel::tc() const scalar pi = thermo_.p()[celli]; scalarField c(Ns_); scalar cSum = 0.0; - + for(label i=0; i& dfdc + scalarSquareMatrix& dfdc ) const; //- Calculates the reaction rates From 5a9cfc954b226cb5316e8b4af4b5ea77d1e66183 Mon Sep 17 00:00:00 2001 From: andy Date: Mon, 29 Sep 2008 10:14:41 +0100 Subject: [PATCH 3/7] == assignments for nut/mut --- src/turbulenceModels/RAS/compressible/LRR/LRR.C | 6 +++--- .../RAS/compressible/LaunderGibsonRSTM/LaunderGibsonRSTM.C | 6 +++--- .../RAS/compressible/RNGkEpsilon/RNGkEpsilon.C | 6 +++--- src/turbulenceModels/RAS/compressible/kEpsilon/kEpsilon.C | 6 +++--- src/turbulenceModels/RAS/compressible/kOmegaSST/kOmegaSST.C | 6 +++--- .../RAS/compressible/realizableKE/realizableKE.C | 6 +++--- src/turbulenceModels/RAS/incompressible/LRR/LRR.C | 4 ++-- .../incompressible/LaunderGibsonRSTM/LaunderGibsonRSTM.C | 4 ++-- .../RAS/incompressible/LienCubicKE/LienCubicKE.C | 4 ++-- .../RAS/incompressible/RNGkEpsilon/RNGkEpsilon.C | 4 ++-- src/turbulenceModels/RAS/incompressible/kEpsilon/kEpsilon.C | 4 ++-- src/turbulenceModels/RAS/incompressible/kOmega/kOmega.C | 4 ++-- .../RAS/incompressible/kOmegaSST/kOmegaSST.C | 4 ++-- .../RAS/incompressible/realizableKE/realizableKE.C | 4 ++-- 14 files changed, 34 insertions(+), 34 deletions(-) diff --git a/src/turbulenceModels/RAS/compressible/LRR/LRR.C b/src/turbulenceModels/RAS/compressible/LRR/LRR.C index 176e2a63c8..31ae4ceba9 100644 --- a/src/turbulenceModels/RAS/compressible/LRR/LRR.C +++ b/src/turbulenceModels/RAS/compressible/LRR/LRR.C @@ -206,7 +206,7 @@ LRR::LRR autoCreateMut("mut", mesh_) ) { - mut_ = Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_); + mut_ == Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_); mut_.correctBoundaryConditions(); if (couplingFactor_.value() < 0.0 || couplingFactor_.value() > 1.0) @@ -310,7 +310,7 @@ void LRR::correct() if (!turbulence_) { // Re-calculate viscosity - mut_ = rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_); + mut_ == rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_); mut_.correctBoundaryConditions(); return; } @@ -399,7 +399,7 @@ void LRR::correct() // Re-calculate viscosity - mut_ = rho_*Cmu_*sqr(k_)/epsilon_; + mut_ == rho_*Cmu_*sqr(k_)/epsilon_; mut_.correctBoundaryConditions(); diff --git a/src/turbulenceModels/RAS/compressible/LaunderGibsonRSTM/LaunderGibsonRSTM.C b/src/turbulenceModels/RAS/compressible/LaunderGibsonRSTM/LaunderGibsonRSTM.C index 0736481d09..4943d73613 100644 --- a/src/turbulenceModels/RAS/compressible/LaunderGibsonRSTM/LaunderGibsonRSTM.C +++ b/src/turbulenceModels/RAS/compressible/LaunderGibsonRSTM/LaunderGibsonRSTM.C @@ -228,7 +228,7 @@ LaunderGibsonRSTM::LaunderGibsonRSTM autoCreateMut("mut", mesh_) ) { - mut_ = Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_); + mut_ == Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_); mut_.correctBoundaryConditions(); if (couplingFactor_.value() < 0.0 || couplingFactor_.value() > 1.0) @@ -335,7 +335,7 @@ void LaunderGibsonRSTM::correct() if (!turbulence_) { // Re-calculate viscosity - mut_ = rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_); + mut_ == rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_); mut_.correctBoundaryConditions(); return; } @@ -438,7 +438,7 @@ void LaunderGibsonRSTM::correct() // Re-calculate turbulent viscosity - mut_ = Cmu_*rho_*sqr(k_)/epsilon_; + mut_ == Cmu_*rho_*sqr(k_)/epsilon_; mut_.correctBoundaryConditions(); diff --git a/src/turbulenceModels/RAS/compressible/RNGkEpsilon/RNGkEpsilon.C b/src/turbulenceModels/RAS/compressible/RNGkEpsilon/RNGkEpsilon.C index 6d54e97f70..8851355ee1 100644 --- a/src/turbulenceModels/RAS/compressible/RNGkEpsilon/RNGkEpsilon.C +++ b/src/turbulenceModels/RAS/compressible/RNGkEpsilon/RNGkEpsilon.C @@ -176,7 +176,7 @@ RNGkEpsilon::RNGkEpsilon autoCreateMut("mut", mesh_) ) { - mut_ = Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_); + mut_ == Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_); mut_.correctBoundaryConditions(); printCoeffs(); @@ -263,7 +263,7 @@ void RNGkEpsilon::correct() if (!turbulence_) { // Re-calculate viscosity - mut_ = rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_); + mut_ == rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_); mut_.correctBoundaryConditions(); return; } @@ -330,7 +330,7 @@ void RNGkEpsilon::correct() // Re-calculate viscosity - mut_ = rho_*Cmu_*sqr(k_)/epsilon_; + mut_ == rho_*Cmu_*sqr(k_)/epsilon_; mut_.correctBoundaryConditions(); } diff --git a/src/turbulenceModels/RAS/compressible/kEpsilon/kEpsilon.C b/src/turbulenceModels/RAS/compressible/kEpsilon/kEpsilon.C index 80a79ab88f..6c10f02466 100644 --- a/src/turbulenceModels/RAS/compressible/kEpsilon/kEpsilon.C +++ b/src/turbulenceModels/RAS/compressible/kEpsilon/kEpsilon.C @@ -157,7 +157,7 @@ kEpsilon::kEpsilon autoCreateMut("mut", mesh_) ) { - mut_ = Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_); + mut_ == Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_); mut_.correctBoundaryConditions(); printCoeffs(); @@ -243,7 +243,7 @@ void kEpsilon::correct() if (!turbulence_) { // Re-calculate viscosity - mut_ = rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_); + mut_ == rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_); mut_.correctBoundaryConditions(); return; } @@ -303,7 +303,7 @@ void kEpsilon::correct() // Re-calculate viscosity - mut_ = rho_*Cmu_*sqr(k_)/epsilon_; + mut_ == rho_*Cmu_*sqr(k_)/epsilon_; mut_.correctBoundaryConditions(); } diff --git a/src/turbulenceModels/RAS/compressible/kOmegaSST/kOmegaSST.C b/src/turbulenceModels/RAS/compressible/kOmegaSST/kOmegaSST.C index 854dcb1bbf..ec44d17c93 100644 --- a/src/turbulenceModels/RAS/compressible/kOmegaSST/kOmegaSST.C +++ b/src/turbulenceModels/RAS/compressible/kOmegaSST/kOmegaSST.C @@ -260,7 +260,7 @@ kOmegaSST::kOmegaSST autoCreateMut("mut", mesh_) ) { - mut_ = a1_*rho_*k_/max(a1_*omega_, F2()*sqrt(magSqr(symm(fvc::grad(U_))))); + mut_ == a1_*rho_*k_/max(a1_*omega_, F2()*sqrt(magSqr(symm(fvc::grad(U_))))); mut_.correctBoundaryConditions(); printCoeffs(); @@ -351,7 +351,7 @@ void kOmegaSST::correct() if (!turbulence_) { // Re-calculate viscosity - mut_ = + mut_ == a1_*rho_*k_ /max(a1_*omega_, F2()*sqrt(magSqr(symm(fvc::grad(U_))))); mut_.correctBoundaryConditions(); @@ -430,7 +430,7 @@ void kOmegaSST::correct() // Re-calculate viscosity - mut_ = a1_*rho_*k_/max(a1_*omega_, F2()*sqrt(S2)); + mut_ == a1_*rho_*k_/max(a1_*omega_, F2()*sqrt(S2)); mut_.correctBoundaryConditions(); } diff --git a/src/turbulenceModels/RAS/compressible/realizableKE/realizableKE.C b/src/turbulenceModels/RAS/compressible/realizableKE/realizableKE.C index fae47c1d40..3806e47ab2 100644 --- a/src/turbulenceModels/RAS/compressible/realizableKE/realizableKE.C +++ b/src/turbulenceModels/RAS/compressible/realizableKE/realizableKE.C @@ -192,7 +192,7 @@ realizableKE::realizableKE bound(k_, k0_); bound(epsilon_, epsilon0_); - mut_ = rCmu(fvc::grad(U_))*rho_*sqr(k_)/(epsilon_ + epsilonSmall_); + mut_ == rCmu(fvc::grad(U_))*rho_*sqr(k_)/(epsilon_ + epsilonSmall_); mut_.correctBoundaryConditions(); printCoeffs(); @@ -276,7 +276,7 @@ void realizableKE::correct() if (!turbulence_) { // Re-calculate viscosity - mut_ = rCmu(fvc::grad(U_))*rho_*sqr(k_)/epsilon_; + mut_ == rCmu(fvc::grad(U_))*rho_*sqr(k_)/epsilon_; mut_.correctBoundaryConditions(); return; } @@ -342,7 +342,7 @@ void realizableKE::correct() bound(k_, k0_); // Re-calculate viscosity - mut_ = rCmu(gradU, S2, magS)*rho_*sqr(k_)/epsilon_; + mut_ == rCmu(gradU, S2, magS)*rho_*sqr(k_)/epsilon_; mut_.correctBoundaryConditions(); } diff --git a/src/turbulenceModels/RAS/incompressible/LRR/LRR.C b/src/turbulenceModels/RAS/incompressible/LRR/LRR.C index 7eebd91b5a..7bf38a9982 100644 --- a/src/turbulenceModels/RAS/incompressible/LRR/LRR.C +++ b/src/turbulenceModels/RAS/incompressible/LRR/LRR.C @@ -187,7 +187,7 @@ LRR::LRR autoCreateNut("nut", mesh_) ) { - nut_ = Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_); + nut_ == Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_); nut_.correctBoundaryConditions(); if (couplingFactor_.value() < 0.0 || couplingFactor_.value() > 1.0) @@ -381,7 +381,7 @@ void LRR::correct() // Re-calculate viscosity - nut_ = Cmu_*sqr(k_)/epsilon_; + nut_ == Cmu_*sqr(k_)/epsilon_; nut_.correctBoundaryConditions(); diff --git a/src/turbulenceModels/RAS/incompressible/LaunderGibsonRSTM/LaunderGibsonRSTM.C b/src/turbulenceModels/RAS/incompressible/LaunderGibsonRSTM/LaunderGibsonRSTM.C index 4e33412681..d67c91a227 100644 --- a/src/turbulenceModels/RAS/incompressible/LaunderGibsonRSTM/LaunderGibsonRSTM.C +++ b/src/turbulenceModels/RAS/incompressible/LaunderGibsonRSTM/LaunderGibsonRSTM.C @@ -216,7 +216,7 @@ LaunderGibsonRSTM::LaunderGibsonRSTM autoCreateNut("nut", mesh_) ) { - nut_ = Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_); + nut_ == Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_); nut_.correctBoundaryConditions(); if (couplingFactor_.value() < 0.0 || couplingFactor_.value() > 1.0) @@ -422,7 +422,7 @@ void LaunderGibsonRSTM::correct() // Re-calculate turbulent viscosity - nut_ = Cmu_*sqr(k_)/epsilon_; + nut_ == Cmu_*sqr(k_)/epsilon_; nut_.correctBoundaryConditions(); diff --git a/src/turbulenceModels/RAS/incompressible/LienCubicKE/LienCubicKE.C b/src/turbulenceModels/RAS/incompressible/LienCubicKE/LienCubicKE.C index 5000d6a135..b347693adc 100644 --- a/src/turbulenceModels/RAS/incompressible/LienCubicKE/LienCubicKE.C +++ b/src/turbulenceModels/RAS/incompressible/LienCubicKE/LienCubicKE.C @@ -228,7 +228,7 @@ LienCubicKE::LienCubicKE ) ) { - nut_ = Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_) + C5viscosity_; + nut_ == Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_) + C5viscosity_; nut_.correctBoundaryConditions(); printCoeffs(); @@ -385,7 +385,7 @@ void LienCubicKE::correct() - 2.0*pow(Cmu_, 3.0)*pow(k_, 4.0)/pow(epsilon_, 3.0) *(magSqr(gradU_ + gradU_.T()) - magSqr(gradU_ - gradU_.T())); - nut_ = Cmu_*sqr(k_)/epsilon_ + C5viscosity_; + nut_ == Cmu_*sqr(k_)/epsilon_ + C5viscosity_; nut_.correctBoundaryConditions(); nonlinearStress_ = symm diff --git a/src/turbulenceModels/RAS/incompressible/RNGkEpsilon/RNGkEpsilon.C b/src/turbulenceModels/RAS/incompressible/RNGkEpsilon/RNGkEpsilon.C index eebcbbdb5e..3c98bfec9e 100644 --- a/src/turbulenceModels/RAS/incompressible/RNGkEpsilon/RNGkEpsilon.C +++ b/src/turbulenceModels/RAS/incompressible/RNGkEpsilon/RNGkEpsilon.C @@ -156,7 +156,7 @@ RNGkEpsilon::RNGkEpsilon autoCreateNut("nut", mesh_) ) { - nut_ = Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_); + nut_ == Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_); nut_.correctBoundaryConditions(); printCoeffs(); @@ -297,7 +297,7 @@ void RNGkEpsilon::correct() // Re-calculate viscosity - nut_ = Cmu_*sqr(k_)/epsilon_; + nut_ == Cmu_*sqr(k_)/epsilon_; nut_.correctBoundaryConditions(); } diff --git a/src/turbulenceModels/RAS/incompressible/kEpsilon/kEpsilon.C b/src/turbulenceModels/RAS/incompressible/kEpsilon/kEpsilon.C index 7e9f971122..c387604301 100644 --- a/src/turbulenceModels/RAS/incompressible/kEpsilon/kEpsilon.C +++ b/src/turbulenceModels/RAS/incompressible/kEpsilon/kEpsilon.C @@ -130,7 +130,7 @@ kEpsilon::kEpsilon autoCreateNut("nut", mesh_) ) { - nut_ = Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_); + nut_ == Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_); nut_.correctBoundaryConditions(); printCoeffs(); @@ -267,7 +267,7 @@ void kEpsilon::correct() // Re-calculate viscosity - nut_ = Cmu_*sqr(k_)/epsilon_; + nut_ == Cmu_*sqr(k_)/epsilon_; nut_.correctBoundaryConditions(); } diff --git a/src/turbulenceModels/RAS/incompressible/kOmega/kOmega.C b/src/turbulenceModels/RAS/incompressible/kOmega/kOmega.C index 47a68d2e18..ff2dceca5a 100644 --- a/src/turbulenceModels/RAS/incompressible/kOmega/kOmega.C +++ b/src/turbulenceModels/RAS/incompressible/kOmega/kOmega.C @@ -141,7 +141,7 @@ kOmega::kOmega autoCreateNut("nut", mesh_) ) { - nut_ = k_/(omega_ + omegaSmall_); + nut_ == k_/(omega_ + omegaSmall_); nut_.correctBoundaryConditions(); printCoeffs(); @@ -273,7 +273,7 @@ void kOmega::correct() // Re-calculate viscosity - nut_ = k_/omega_; + nut_ == k_/omega_; nut_.correctBoundaryConditions(); } diff --git a/src/turbulenceModels/RAS/incompressible/kOmegaSST/kOmegaSST.C b/src/turbulenceModels/RAS/incompressible/kOmegaSST/kOmegaSST.C index 2c66a36f30..4760c17152 100644 --- a/src/turbulenceModels/RAS/incompressible/kOmegaSST/kOmegaSST.C +++ b/src/turbulenceModels/RAS/incompressible/kOmegaSST/kOmegaSST.C @@ -250,7 +250,7 @@ kOmegaSST::kOmegaSST autoCreateNut("nut", mesh_) ) { - nut_ = + nut_ == a1_*k_ /max ( @@ -411,7 +411,7 @@ void kOmegaSST::correct() // Re-calculate viscosity - nut_ = a1_*k_/max(a1_*omega_, F2()*sqrt(S2)); + nut_ == a1_*k_/max(a1_*omega_, F2()*sqrt(S2)); nut_.correctBoundaryConditions(); } diff --git a/src/turbulenceModels/RAS/incompressible/realizableKE/realizableKE.C b/src/turbulenceModels/RAS/incompressible/realizableKE/realizableKE.C index 294bf9c915..ffeba065d7 100644 --- a/src/turbulenceModels/RAS/incompressible/realizableKE/realizableKE.C +++ b/src/turbulenceModels/RAS/incompressible/realizableKE/realizableKE.C @@ -182,7 +182,7 @@ realizableKE::realizableKE bound(k_, k0_); bound(epsilon_, epsilon0_); - nut_ = rCmu(fvc::grad(U_))*sqr(k_)/(epsilon_ + epsilonSmall_); + nut_ == rCmu(fvc::grad(U_))*sqr(k_)/(epsilon_ + epsilonSmall_); nut_.correctBoundaryConditions(); printCoeffs(); @@ -326,7 +326,7 @@ void realizableKE::correct() // Re-calculate viscosity - nut_ = rCmu(gradU, S2, magS)*sqr(k_)/epsilon_; + nut_ == rCmu(gradU, S2, magS)*sqr(k_)/epsilon_; nut_.correctBoundaryConditions(); } From 071b8ffe89451c5c3f85d86ad00de22b475a3ea3 Mon Sep 17 00:00:00 2001 From: andy Date: Tue, 30 Sep 2008 14:26:46 +0100 Subject: [PATCH 4/7] adding run-time selectable thermal wall function capability to compressible RAS models --- .../RAS/compressible/LRR/LRR.C | 30 +++- .../RAS/compressible/LRR/LRR.H | 3 +- .../LaunderGibsonRSTM/LaunderGibsonRSTM.C | 29 +++- .../LaunderGibsonRSTM/LaunderGibsonRSTM.H | 3 +- .../RAS/compressible/Make/files | 3 + .../RAS/compressible/RASModel/RASModel.C | 29 +++- .../RAS/compressible/RASModel/RASModel.H | 7 + .../compressible/RNGkEpsilon/RNGkEpsilon.C | 24 +++ .../compressible/RNGkEpsilon/RNGkEpsilon.H | 3 +- .../backwardsCompatibilityWallFunctions.C | 71 ++++++++ .../backwardsCompatibilityWallFunctions.H | 7 + .../alphatWallFunctionFvPatchScalarField.C | 132 +++++++++++++++ .../alphatWallFunctionFvPatchScalarField.H | 155 ++++++++++++++++++ .../RAS/compressible/kEpsilon/kEpsilon.C | 24 +++ .../RAS/compressible/kEpsilon/kEpsilon.H | 4 +- .../RAS/compressible/kOmegaSST/kOmegaSST.C | 23 +++ .../RAS/compressible/kOmegaSST/kOmegaSST.H | 3 +- .../compressible/realizableKE/realizableKE.C | 24 +++ .../compressible/realizableKE/realizableKE.H | 3 +- .../RAS/incompressible/RASModel/RASModel.C | 15 +- 20 files changed, 563 insertions(+), 29 deletions(-) create mode 100644 src/turbulenceModels/RAS/compressible/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatWallFunction/alphatWallFunctionFvPatchScalarField.C create mode 100644 src/turbulenceModels/RAS/compressible/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatWallFunction/alphatWallFunctionFvPatchScalarField.H diff --git a/src/turbulenceModels/RAS/compressible/LRR/LRR.C b/src/turbulenceModels/RAS/compressible/LRR/LRR.C index 31ae4ceba9..a22d495249 100644 --- a/src/turbulenceModels/RAS/compressible/LRR/LRR.C +++ b/src/turbulenceModels/RAS/compressible/LRR/LRR.C @@ -204,11 +204,20 @@ LRR::LRR IOobject::AUTO_WRITE ), autoCreateMut("mut", mesh_) + ), + alphat_ + ( + IOobject + ( + "alphat", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + autoCreateAlphat("alphat", mesh_) ) { - mut_ == Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_); - mut_.correctBoundaryConditions(); - if (couplingFactor_.value() < 0.0 || couplingFactor_.value() > 1.0) { FatalErrorIn @@ -221,6 +230,12 @@ LRR::LRR << exit(FatalError); } + mut_ == Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_); + mut_.correctBoundaryConditions(); + + alphat_ == mut_/Prt_; + alphat_.correctBoundaryConditions(); + printCoeffs(); } @@ -312,6 +327,11 @@ void LRR::correct() // Re-calculate viscosity mut_ == rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_); mut_.correctBoundaryConditions(); + + // Re-calculate thermal diffusivity + alphat_ = mut_/Prt_; + alphat_.correctBoundaryConditions(); + return; } @@ -402,6 +422,10 @@ void LRR::correct() mut_ == rho_*Cmu_*sqr(k_)/epsilon_; mut_.correctBoundaryConditions(); + // Re-calculate thermal diffusivity + alphat_ = mut_/Prt_; + alphat_.correctBoundaryConditions(); + // Correct wall shear stresses diff --git a/src/turbulenceModels/RAS/compressible/LRR/LRR.H b/src/turbulenceModels/RAS/compressible/LRR/LRR.H index 553b11be96..0dc2ed2ad8 100644 --- a/src/turbulenceModels/RAS/compressible/LRR/LRR.H +++ b/src/turbulenceModels/RAS/compressible/LRR/LRR.H @@ -97,6 +97,7 @@ class LRR volScalarField k_; volScalarField epsilon_; volScalarField mut_; + volScalarField alphat_; public: @@ -152,7 +153,7 @@ public: { return tmp ( - new volScalarField("alphaEff", alphah_*mut_ + alpha()) + new volScalarField("alphaEff", alphah_*alphat_ + alpha()) ); } diff --git a/src/turbulenceModels/RAS/compressible/LaunderGibsonRSTM/LaunderGibsonRSTM.C b/src/turbulenceModels/RAS/compressible/LaunderGibsonRSTM/LaunderGibsonRSTM.C index 4943d73613..0181605fb1 100644 --- a/src/turbulenceModels/RAS/compressible/LaunderGibsonRSTM/LaunderGibsonRSTM.C +++ b/src/turbulenceModels/RAS/compressible/LaunderGibsonRSTM/LaunderGibsonRSTM.C @@ -226,11 +226,20 @@ LaunderGibsonRSTM::LaunderGibsonRSTM IOobject::AUTO_WRITE ), autoCreateMut("mut", mesh_) + ), + alphat_ + ( + IOobject + ( + "alphat", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + autoCreateAlphat("alphat", mesh_) ) { - mut_ == Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_); - mut_.correctBoundaryConditions(); - if (couplingFactor_.value() < 0.0 || couplingFactor_.value() > 1.0) { FatalErrorIn @@ -243,6 +252,12 @@ LaunderGibsonRSTM::LaunderGibsonRSTM << exit(FatalError); } + mut_ == Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_); + mut_.correctBoundaryConditions(); + + alphat_ == mut_/Prt_; + alphat_.correctBoundaryConditions(); + printCoeffs(); } @@ -337,6 +352,11 @@ void LaunderGibsonRSTM::correct() // Re-calculate viscosity mut_ == rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_); mut_.correctBoundaryConditions(); + + // Re-calculate thermal diffusivity + alphat_ = mut_/Prt_; + alphat_.correctBoundaryConditions(); + return; } @@ -441,6 +461,9 @@ void LaunderGibsonRSTM::correct() mut_ == Cmu_*rho_*sqr(k_)/epsilon_; mut_.correctBoundaryConditions(); + // Re-calculate thermal diffusivity + alphat_ = mut_/Prt_; + alphat_.correctBoundaryConditions(); // Correct wall shear stresses diff --git a/src/turbulenceModels/RAS/compressible/LaunderGibsonRSTM/LaunderGibsonRSTM.H b/src/turbulenceModels/RAS/compressible/LaunderGibsonRSTM/LaunderGibsonRSTM.H index 2be69e96b5..4fa0e9315e 100644 --- a/src/turbulenceModels/RAS/compressible/LaunderGibsonRSTM/LaunderGibsonRSTM.H +++ b/src/turbulenceModels/RAS/compressible/LaunderGibsonRSTM/LaunderGibsonRSTM.H @@ -104,6 +104,7 @@ class LaunderGibsonRSTM volScalarField k_; volScalarField epsilon_; volScalarField mut_; + volScalarField alphat_; public: @@ -161,7 +162,7 @@ public: { return tmp ( - new volScalarField("alphaEff", alphah_*mut_ + alpha()) + new volScalarField("alphaEff", alphah_*alphat_ + alpha()) ); } diff --git a/src/turbulenceModels/RAS/compressible/Make/files b/src/turbulenceModels/RAS/compressible/Make/files index ff6c31794b..0598d06a00 100644 --- a/src/turbulenceModels/RAS/compressible/Make/files +++ b/src/turbulenceModels/RAS/compressible/Make/files @@ -14,6 +14,9 @@ kOmegaSST/kOmegaSST.C /* Wall functions */ wallFunctions = derivedFvPatchFields/wallFunctions +alphatWallFunctions = $(wallFunctions)/alphatWallFunctions +$(alphatWallFunctions)/alphatWallFunction/alphatWallFunctionFvPatchScalarField.C + mutWallFunctions = $(wallFunctions)/mutWallFunctions $(mutWallFunctions)/mutWallFunction/mutWallFunctionFvPatchScalarField.C $(mutWallFunctions)/mutRoughWallFunction/mutRoughWallFunctionFvPatchScalarField.C diff --git a/src/turbulenceModels/RAS/compressible/RASModel/RASModel.C b/src/turbulenceModels/RAS/compressible/RASModel/RASModel.C index a38da49b0e..c1bbdc1684 100644 --- a/src/turbulenceModels/RAS/compressible/RASModel/RASModel.C +++ b/src/turbulenceModels/RAS/compressible/RASModel/RASModel.C @@ -47,7 +47,8 @@ void RASModel::printCoeffs() { if (printCoeffs_) { - Info<< type() << "Coeffs" << coeffDict_ << endl; + Info<< type() << "Coeffs" << coeffDict_ << nl + << "wallFunctionCoeffs" << wallFunctionDict_ << endl; } } @@ -115,6 +116,15 @@ RASModel::RASModel 0.09 ) ), + Prt_ + ( + dimensioned::lookupOrAddToDict + ( + "Prt", + wallFunctionDict_, + 0.85 + ) + ), yPlusLam_(yPlusLam(kappa_.value(), E_.value())), @@ -148,11 +158,9 @@ tmp RASModel::yPlus(const label patchNo) const tmp tYp(new scalarField(curPatch.size())); scalarField& Yp = tYp(); - if (typeid(curPatch) == typeid(wallFvPatch)) + if (isType(curPatch)) { - scalar Cmu(readScalar(coeffDict_.lookup("Cmu"))); - - Yp = pow(Cmu, 0.25) + Yp = pow(Cmu_.value(), 0.25) *y_[patchNo] *sqrt(k()().boundaryField()[patchNo].patchInternalField()) /( @@ -165,8 +173,8 @@ tmp RASModel::yPlus(const label patchNo) const WarningIn ( "tmp RASModel::yPlus(const label patchNo) const" - ) << "Patch " << patchNo << " is not a wall. Returning blank field" - << endl; + ) << "Patch " << patchNo << " is not a wall. Returning null field" + << nl << endl; Yp.setSize(0); } @@ -191,8 +199,11 @@ bool RASModel::read() lookup("turbulence") >> turbulence_; coeffDict_ = subDict(type() + "Coeffs"); - kappa_.readIfPresent(subDict("wallFunctionCoeffs")); - E_.readIfPresent(subDict("wallFunctionCoeffs")); + wallFunctionDict_ = subDict("wallFunctionCoeffs"); + kappa_.readIfPresent(wallFunctionDict_); + E_.readIfPresent(wallFunctionDict_); + Cmu_.readIfPresent(wallFunctionDict_); + Prt_.readIfPresent(wallFunctionDict_); yPlusLam_ = yPlusLam(kappa_.value(), E_.value()); diff --git a/src/turbulenceModels/RAS/compressible/RASModel/RASModel.H b/src/turbulenceModels/RAS/compressible/RASModel/RASModel.H index 51b4e706da..1be29e6e79 100644 --- a/src/turbulenceModels/RAS/compressible/RASModel/RASModel.H +++ b/src/turbulenceModels/RAS/compressible/RASModel/RASModel.H @@ -95,6 +95,7 @@ protected: dimensionedScalar kappa_; dimensionedScalar E_; dimensionedScalar Cmu_; + dimensionedScalar Prt_; scalar yPlusLam_; @@ -244,6 +245,12 @@ public: return Cmu_; } + //- Return turbulent Prandtl number for use in wall-functions + dimensionedScalar Prt() const + { + return Prt_; + } + //- Return the near wall distances const nearWallDist& y() const { diff --git a/src/turbulenceModels/RAS/compressible/RNGkEpsilon/RNGkEpsilon.C b/src/turbulenceModels/RAS/compressible/RNGkEpsilon/RNGkEpsilon.C index 8851355ee1..8f6ba82321 100644 --- a/src/turbulenceModels/RAS/compressible/RNGkEpsilon/RNGkEpsilon.C +++ b/src/turbulenceModels/RAS/compressible/RNGkEpsilon/RNGkEpsilon.C @@ -174,11 +174,26 @@ RNGkEpsilon::RNGkEpsilon IOobject::AUTO_WRITE ), autoCreateMut("mut", mesh_) + ), + alphat_ + ( + IOobject + ( + "alphat", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + autoCreateAlphat("alphat", mesh_) ) { mut_ == Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_); mut_.correctBoundaryConditions(); + alphat_ == mut_/Prt_; + alphat_.correctBoundaryConditions(); + printCoeffs(); } @@ -265,6 +280,11 @@ void RNGkEpsilon::correct() // Re-calculate viscosity mut_ == rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_); mut_.correctBoundaryConditions(); + + // Re-calculate thermal diffusivity + alphat_ = mut_/Prt_; + alphat_.correctBoundaryConditions(); + return; } @@ -332,6 +352,10 @@ void RNGkEpsilon::correct() // Re-calculate viscosity mut_ == rho_*Cmu_*sqr(k_)/epsilon_; mut_.correctBoundaryConditions(); + + // Re-calculate thermal diffusivity + alphat_ = mut_/Prt_; + alphat_.correctBoundaryConditions(); } diff --git a/src/turbulenceModels/RAS/compressible/RNGkEpsilon/RNGkEpsilon.H b/src/turbulenceModels/RAS/compressible/RNGkEpsilon/RNGkEpsilon.H index 5063615ccb..98445ced1d 100644 --- a/src/turbulenceModels/RAS/compressible/RNGkEpsilon/RNGkEpsilon.H +++ b/src/turbulenceModels/RAS/compressible/RNGkEpsilon/RNGkEpsilon.H @@ -87,6 +87,7 @@ class RNGkEpsilon volScalarField k_; volScalarField epsilon_; volScalarField mut_; + volScalarField alphat_; public: @@ -142,7 +143,7 @@ public: { return tmp ( - new volScalarField("alphaEff", alphah_*mut_ + alpha()) + new volScalarField("alphaEff", alphah_*alphat_ + alpha()) ); } diff --git a/src/turbulenceModels/RAS/compressible/backwardsCompatibilityWallFunctions/backwardsCompatibilityWallFunctions.C b/src/turbulenceModels/RAS/compressible/backwardsCompatibilityWallFunctions/backwardsCompatibilityWallFunctions.C index a7c93d4ad6..19df95a3d5 100644 --- a/src/turbulenceModels/RAS/compressible/backwardsCompatibilityWallFunctions/backwardsCompatibilityWallFunctions.C +++ b/src/turbulenceModels/RAS/compressible/backwardsCompatibilityWallFunctions/backwardsCompatibilityWallFunctions.C @@ -27,6 +27,7 @@ License #include "backwardsCompatibilityWallFunctions.H" #include "calculatedFvPatchField.H" +#include "alphatWallFunctionFvPatchScalarField.H" #include "mutWallFunctionFvPatchScalarField.H" #include "epsilonWallFunctionFvPatchScalarField.H" #include "kQRWallFunctionFvPatchField.H" @@ -41,6 +42,76 @@ namespace compressible // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +tmp autoCreateAlphat +( + const word& fieldName, + const fvMesh& mesh +) +{ + IOobject alphatHeader + ( + fieldName, + mesh.time().timeName(), + mesh, + IOobject::MUST_READ, + IOobject::NO_WRITE, + false + ); + + if (alphatHeader.headerOk()) + { + return tmp(new volScalarField(alphatHeader, mesh)); + } + else + { + Info<< "--> Upgrading " << fieldName << " to employ run-time " + << "selectable wall functions" << endl; + + const fvBoundaryMesh& bm = mesh.boundary(); + + wordList alphatBoundaryTypes(bm.size()); + + forAll(bm, patchI) + { + if (isType(bm[patchI])) + { + alphatBoundaryTypes[patchI] = + RASModels::alphatWallFunctionFvPatchScalarField::typeName; + } + else + { + alphatBoundaryTypes[patchI] = + calculatedFvPatchField::typeName; + } + } + + tmp alphat + ( + new volScalarField + ( + IOobject + ( + fieldName, + mesh.time().timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + mesh, + dimensionedScalar("zero", dimDensity*dimArea/dimTime, 0.0), + alphatBoundaryTypes + ) + ); + + Info<< " Writing updated " << fieldName << endl; + alphat().write(); + + return alphat; + } +} + + tmp autoCreateMut ( const word& fieldName, diff --git a/src/turbulenceModels/RAS/compressible/backwardsCompatibilityWallFunctions/backwardsCompatibilityWallFunctions.H b/src/turbulenceModels/RAS/compressible/backwardsCompatibilityWallFunctions/backwardsCompatibilityWallFunctions.H index a97b3090cd..bcc812a6b4 100644 --- a/src/turbulenceModels/RAS/compressible/backwardsCompatibilityWallFunctions/backwardsCompatibilityWallFunctions.H +++ b/src/turbulenceModels/RAS/compressible/backwardsCompatibilityWallFunctions/backwardsCompatibilityWallFunctions.H @@ -53,6 +53,13 @@ namespace compressible const fvMesh& mesh ); + //- alphat + tmp autoCreateAlphat + ( + const word& fieldName, + const fvMesh& mesh + ); + //- epsilon tmp autoCreateEpsilon ( diff --git a/src/turbulenceModels/RAS/compressible/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatWallFunction/alphatWallFunctionFvPatchScalarField.C b/src/turbulenceModels/RAS/compressible/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatWallFunction/alphatWallFunctionFvPatchScalarField.C new file mode 100644 index 0000000000..a34631cc3c --- /dev/null +++ b/src/turbulenceModels/RAS/compressible/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatWallFunction/alphatWallFunctionFvPatchScalarField.C @@ -0,0 +1,132 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software; you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation; either version 2 of the License, or (at your + option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "alphatWallFunctionFvPatchScalarField.H" +#include "RASModel.H" +#include "fvPatchFieldMapper.H" +#include "volFields.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace compressible +{ +namespace RASModels +{ + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +alphatWallFunctionFvPatchScalarField:: +alphatWallFunctionFvPatchScalarField +( + const fvPatch& p, + const DimensionedField& iF +) +: + fixedValueFvPatchScalarField(p, iF) +{} + + +alphatWallFunctionFvPatchScalarField:: +alphatWallFunctionFvPatchScalarField +( + const alphatWallFunctionFvPatchScalarField& ptf, + const fvPatch& p, + const DimensionedField& iF, + const fvPatchFieldMapper& mapper +) +: + fixedValueFvPatchScalarField(ptf, p, iF, mapper) +{} + + +alphatWallFunctionFvPatchScalarField:: +alphatWallFunctionFvPatchScalarField +( + const fvPatch& p, + const DimensionedField& iF, + const dictionary& dict +) +: + fixedValueFvPatchScalarField(p, iF, dict) +{} + + +alphatWallFunctionFvPatchScalarField:: +alphatWallFunctionFvPatchScalarField +( + const alphatWallFunctionFvPatchScalarField& awfpsf +) +: + fixedValueFvPatchScalarField(awfpsf) +{} + + +alphatWallFunctionFvPatchScalarField:: +alphatWallFunctionFvPatchScalarField +( + const alphatWallFunctionFvPatchScalarField& awfpsf, + const DimensionedField& iF +) +: + fixedValueFvPatchScalarField(awfpsf, iF) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void alphatWallFunctionFvPatchScalarField::updateCoeffs() +{ + const RASModel& ras = db().lookupObject("RASProperties"); + const scalar Prt = ras.Prt().value(); + + const scalarField& mutw = + patch().lookupPatchField("mut"); + + operator==(mutw/Prt); +} + + +void alphatWallFunctionFvPatchScalarField::write(Ostream& os) const +{ + fvPatchField::write(os); + writeEntry("value", os); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +makePatchTypeField(fvPatchScalarField, alphatWallFunctionFvPatchScalarField); + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace RASModels +} // End namespace compressible +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/turbulenceModels/RAS/compressible/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatWallFunction/alphatWallFunctionFvPatchScalarField.H b/src/turbulenceModels/RAS/compressible/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatWallFunction/alphatWallFunctionFvPatchScalarField.H new file mode 100644 index 0000000000..d94775a342 --- /dev/null +++ b/src/turbulenceModels/RAS/compressible/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatWallFunction/alphatWallFunctionFvPatchScalarField.H @@ -0,0 +1,155 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software; you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation; either version 2 of the License, or (at your + option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Class + Foam::compressible::RASModels::alphatWallFunctionFvPatchScalarField + +Description + Boundary condition for turbulent thermal diffusivity when using wall + functions + - replicates OpenFOAM v1.5 (and earlier) behaviour + +SourceFiles + alphatWallFunctionFvPatchScalarField.C + +\*---------------------------------------------------------------------------*/ + +#ifndef alphatWallFunctionFvPatchScalarField_H +#define alphatWallFunctionFvPatchScalarField_H + +#include "fixedValueFvPatchFields.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace compressible +{ +namespace RASModels +{ + +/*---------------------------------------------------------------------------*\ + Class alphatWallFunctionFvPatchScalarField Declaration +\*---------------------------------------------------------------------------*/ + +class alphatWallFunctionFvPatchScalarField +: + public fixedValueFvPatchScalarField +{ + +public: + + //- Runtime type information + TypeName("alphatWallFunction"); + + + // Constructors + + //- Construct from patch and internal field + alphatWallFunctionFvPatchScalarField + ( + const fvPatch&, + const DimensionedField& + ); + + //- Construct from patch, internal field and dictionary + alphatWallFunctionFvPatchScalarField + ( + const fvPatch&, + const DimensionedField&, + const dictionary& + ); + + //- Construct by mapping given + // alphatWallFunctionFvPatchScalarField + // onto a new patch + alphatWallFunctionFvPatchScalarField + ( + const alphatWallFunctionFvPatchScalarField&, + const fvPatch&, + const DimensionedField&, + const fvPatchFieldMapper& + ); + + //- Construct as copy + alphatWallFunctionFvPatchScalarField + ( + const alphatWallFunctionFvPatchScalarField& + ); + + //- Construct and return a clone + virtual tmp clone() const + { + return tmp + ( + new alphatWallFunctionFvPatchScalarField(*this) + ); + } + + //- Construct as copy setting internal field reference + alphatWallFunctionFvPatchScalarField + ( + const alphatWallFunctionFvPatchScalarField&, + const DimensionedField& + ); + + //- Construct and return a clone setting internal field reference + virtual tmp clone + ( + const DimensionedField& iF + ) const + { + return tmp + ( + new alphatWallFunctionFvPatchScalarField(*this, iF) + ); + } + + + // Member functions + + // Evaluation functions + + //- Update the coefficients associated with the patch field + virtual void updateCoeffs(); + + + // I-O + + //- Write + void write(Ostream&) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace RASModels +} // End namespace compressible +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/turbulenceModels/RAS/compressible/kEpsilon/kEpsilon.C b/src/turbulenceModels/RAS/compressible/kEpsilon/kEpsilon.C index 6c10f02466..cb0fecc03d 100644 --- a/src/turbulenceModels/RAS/compressible/kEpsilon/kEpsilon.C +++ b/src/turbulenceModels/RAS/compressible/kEpsilon/kEpsilon.C @@ -155,11 +155,26 @@ kEpsilon::kEpsilon IOobject::AUTO_WRITE ), autoCreateMut("mut", mesh_) + ), + alphat_ + ( + IOobject + ( + "alphat", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + autoCreateAlphat("alphat", mesh_) ) { mut_ == Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_); mut_.correctBoundaryConditions(); + alphat_ == mut_/Prt_; + alphat_.correctBoundaryConditions(); + printCoeffs(); } @@ -245,6 +260,11 @@ void kEpsilon::correct() // Re-calculate viscosity mut_ == rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_); mut_.correctBoundaryConditions(); + + // Re-calculate thermal diffusivity + alphat_ = mut_/Prt_; + alphat_.correctBoundaryConditions(); + return; } @@ -305,6 +325,10 @@ void kEpsilon::correct() // Re-calculate viscosity mut_ == rho_*Cmu_*sqr(k_)/epsilon_; mut_.correctBoundaryConditions(); + + // Re-calculate thermal diffusivity + alphat_ = mut_/Prt_; + alphat_.correctBoundaryConditions(); } diff --git a/src/turbulenceModels/RAS/compressible/kEpsilon/kEpsilon.H b/src/turbulenceModels/RAS/compressible/kEpsilon/kEpsilon.H index d42bc73a66..76d952d40d 100644 --- a/src/turbulenceModels/RAS/compressible/kEpsilon/kEpsilon.H +++ b/src/turbulenceModels/RAS/compressible/kEpsilon/kEpsilon.H @@ -74,7 +74,6 @@ class kEpsilon // Model coefficients -// dimensionedScalar Cmu; dimensionedScalar Cmu_; dimensionedScalar C1_; dimensionedScalar C2_; @@ -88,6 +87,7 @@ class kEpsilon volScalarField k_; volScalarField epsilon_; volScalarField mut_; + volScalarField alphat_; public: @@ -144,7 +144,7 @@ public: { return tmp ( - new volScalarField("alphaEff", alphah_*mut_ + alpha()) + new volScalarField("alphaEff", alphah_*alphat_ + alpha()) ); } diff --git a/src/turbulenceModels/RAS/compressible/kOmegaSST/kOmegaSST.C b/src/turbulenceModels/RAS/compressible/kOmegaSST/kOmegaSST.C index ec44d17c93..1181169f87 100644 --- a/src/turbulenceModels/RAS/compressible/kOmegaSST/kOmegaSST.C +++ b/src/turbulenceModels/RAS/compressible/kOmegaSST/kOmegaSST.C @@ -258,11 +258,26 @@ kOmegaSST::kOmegaSST IOobject::AUTO_WRITE ), autoCreateMut("mut", mesh_) + ), + alphat_ + ( + IOobject + ( + "alphat", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + autoCreateAlphat("alphat", mesh_) ) { mut_ == a1_*rho_*k_/max(a1_*omega_, F2()*sqrt(magSqr(symm(fvc::grad(U_))))); mut_.correctBoundaryConditions(); + alphat_ == mut_/Prt_; + alphat_.correctBoundaryConditions(); + printCoeffs(); } @@ -356,6 +371,10 @@ void kOmegaSST::correct() /max(a1_*omega_, F2()*sqrt(magSqr(symm(fvc::grad(U_))))); mut_.correctBoundaryConditions(); + // Re-calculate thermal diffusivity + alphat_ = mut_/Prt_; + alphat_.correctBoundaryConditions(); + return; } @@ -432,6 +451,10 @@ void kOmegaSST::correct() // Re-calculate viscosity mut_ == a1_*rho_*k_/max(a1_*omega_, F2()*sqrt(S2)); mut_.correctBoundaryConditions(); + + // Re-calculate thermal diffusivity + alphat_ = mut_/Prt_; + alphat_.correctBoundaryConditions(); } diff --git a/src/turbulenceModels/RAS/compressible/kOmegaSST/kOmegaSST.H b/src/turbulenceModels/RAS/compressible/kOmegaSST/kOmegaSST.H index 4836d02611..3f9c1486d3 100644 --- a/src/turbulenceModels/RAS/compressible/kOmegaSST/kOmegaSST.H +++ b/src/turbulenceModels/RAS/compressible/kOmegaSST/kOmegaSST.H @@ -134,6 +134,7 @@ class kOmegaSST volScalarField k_; volScalarField omega_; volScalarField mut_; + volScalarField alphat_; // Private member functions @@ -238,7 +239,7 @@ public: { return tmp ( - new volScalarField("alphaEff", alphah_*mut_ + alpha()) + new volScalarField("alphaEff", alphah_*alphat_ + alpha()) ); } diff --git a/src/turbulenceModels/RAS/compressible/realizableKE/realizableKE.C b/src/turbulenceModels/RAS/compressible/realizableKE/realizableKE.C index 3806e47ab2..1e762433f3 100644 --- a/src/turbulenceModels/RAS/compressible/realizableKE/realizableKE.C +++ b/src/turbulenceModels/RAS/compressible/realizableKE/realizableKE.C @@ -187,6 +187,18 @@ realizableKE::realizableKE IOobject::AUTO_WRITE ), autoCreateMut("mut", mesh_) + ), + alphat_ + ( + IOobject + ( + "alphat", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + autoCreateAlphat("alphat", mesh_) ) { bound(k_, k0_); @@ -195,6 +207,9 @@ realizableKE::realizableKE mut_ == rCmu(fvc::grad(U_))*rho_*sqr(k_)/(epsilon_ + epsilonSmall_); mut_.correctBoundaryConditions(); + alphat_ == mut_/Prt_; + alphat_.correctBoundaryConditions(); + printCoeffs(); } @@ -278,6 +293,11 @@ void realizableKE::correct() // Re-calculate viscosity mut_ == rCmu(fvc::grad(U_))*rho_*sqr(k_)/epsilon_; mut_.correctBoundaryConditions(); + + // Re-calculate thermal diffusivity + alphat_ = mut_/Prt_; + alphat_.correctBoundaryConditions(); + return; } @@ -344,6 +364,10 @@ void realizableKE::correct() // Re-calculate viscosity mut_ == rCmu(gradU, S2, magS)*rho_*sqr(k_)/epsilon_; mut_.correctBoundaryConditions(); + + // Re-calculate thermal diffusivity + alphat_ = mut_/Prt_; + alphat_.correctBoundaryConditions(); } diff --git a/src/turbulenceModels/RAS/compressible/realizableKE/realizableKE.H b/src/turbulenceModels/RAS/compressible/realizableKE/realizableKE.H index 0eb770f083..f997efdd5f 100644 --- a/src/turbulenceModels/RAS/compressible/realizableKE/realizableKE.H +++ b/src/turbulenceModels/RAS/compressible/realizableKE/realizableKE.H @@ -91,6 +91,7 @@ class realizableKE volScalarField k_; volScalarField epsilon_; volScalarField mut_; + volScalarField alphat_; tmp rCmu ( @@ -157,7 +158,7 @@ public: { return tmp ( - new volScalarField("alphaEff", alphah_*mut_ + alpha()) + new volScalarField("alphaEff", alphah_*alphat_ + alpha()) ); } diff --git a/src/turbulenceModels/RAS/incompressible/RASModel/RASModel.C b/src/turbulenceModels/RAS/incompressible/RASModel/RASModel.C index 628da387da..46236828f2 100644 --- a/src/turbulenceModels/RAS/incompressible/RASModel/RASModel.C +++ b/src/turbulenceModels/RAS/incompressible/RASModel/RASModel.C @@ -45,7 +45,8 @@ void RASModel::printCoeffs() { if (printCoeffs_) { - Info<< type() << "Coeffs" << coeffDict_ << endl;; + Info<< type() << "Coeffs" << coeffDict_ << nl + << "wallFunctionCoeffs" << wallFunctionDict_ << endl; } } @@ -148,9 +149,10 @@ tmp RASModel::yPlus(const label patchNo) const tmp tYp(new scalarField(curPatch.size())); scalarField& Yp = tYp(); - if (typeid(curPatch) == typeid(wallFvPatch)) + if (isType(curPatch)) { - Yp = pow(Cmu_.value(), 0.25)*y_[patchNo] + Yp = pow(Cmu_.value(), 0.25) + *y_[patchNo] *sqrt(k()().boundaryField()[patchNo].patchInternalField()) /nu().boundaryField()[patchNo]; } @@ -158,9 +160,8 @@ tmp RASModel::yPlus(const label patchNo) const { WarningIn ( - "tmp RASModel::yPlus(const label patchNo)" - ) << "const : " << nl - << "Patch " << patchNo << " is not a wall. Returning zero field" + "tmp RASModel::yPlus(const label patchNo) const" + ) << "Patch " << patchNo << " is not a wall. Returning null field" << nl << endl; Yp.setSize(0); @@ -185,8 +186,8 @@ bool RASModel::read() { lookup("turbulence") >> turbulence_; coeffDict_ = subDict(type() + "Coeffs"); - wallFunctionDict_ = subDict("wallFunctionCoeffs"); + wallFunctionDict_ = subDict("wallFunctionCoeffs"); kappa_.readIfPresent(wallFunctionDict_); E_.readIfPresent(wallFunctionDict_); Cmu_.readIfPresent(wallFunctionDict_); From cc0dbb05499e55521d6b291a7d14c2fb55e61f0e Mon Sep 17 00:00:00 2001 From: andy Date: Tue, 30 Sep 2008 15:57:55 +0100 Subject: [PATCH 5/7] adding -fields option --- .../reconstructPar/fvFieldReconstructor.H | 10 ++-- .../fvFieldReconstructorReconstructFields.C | 56 ++++++++++--------- .../reconstructPar/reconstructPar.C | 44 +++++++++++++-- 3 files changed, 74 insertions(+), 36 deletions(-) diff --git a/applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructor.H b/applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructor.H index 201ce66811..71ac334133 100644 --- a/applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructor.H +++ b/applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructor.H @@ -151,18 +151,20 @@ public: const IOobject& fieldIoObject ); - //- Reconstruct and write all volume fields + //- Reconstruct and write all/selected volume fields template void reconstructFvVolumeFields ( - const IOobjectList& objects + const IOobjectList& objects, + const HashSet& selectedFields ); - //- Reconstruct and write all volume fields + //- Reconstruct and write all/selected volume fields template void reconstructFvSurfaceFields ( - const IOobjectList& objects + const IOobjectList& objects, + const HashSet& selectedFields ); }; diff --git a/applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructorReconstructFields.C b/applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructorReconstructFields.C index 06f4f6a634..2ad0735ccf 100644 --- a/applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructorReconstructFields.C +++ b/applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructorReconstructFields.C @@ -131,7 +131,7 @@ Foam::fvFieldReconstructor::reconstructFvVolumeField forAll(cp, faceI) { // Subtract one to take into account offsets for - // face direction. + // face direction. reverseAddressing[faceI] = cp[faceI] - 1 - curPatchStart; } @@ -151,7 +151,7 @@ Foam::fvFieldReconstructor::reconstructFvVolumeField forAll(cp, faceI) { // Subtract one to take into account offsets for - // face direction. + // face direction. label curF = cp[faceI] - 1; // Is the face on the boundary? @@ -282,7 +282,7 @@ Foam::fvFieldReconstructor::reconstructFvSurfaceField // It is necessary to create a copy of the addressing array to // take care of the face direction offset trick. - // + // { labelList curAddr(faceProcAddressing_[procI]); @@ -342,7 +342,7 @@ Foam::fvFieldReconstructor::reconstructFvSurfaceField forAll(cp, faceI) { // Subtract one to take into account offsets for - // face direction. + // face direction. reverseAddressing[faceI] = cp[faceI] - 1 - curPatchStart; } @@ -452,11 +452,12 @@ Foam::fvFieldReconstructor::reconstructFvSurfaceField } -// Reconstruct and write all volume fields +// Reconstruct and write all/selected volume fields template void Foam::fvFieldReconstructor::reconstructFvVolumeFields ( - const IOobjectList& objects + const IOobjectList& objects, + const HashSet& selectedFields ) { const word& fieldClassName = @@ -468,27 +469,29 @@ void Foam::fvFieldReconstructor::reconstructFvVolumeFields { Info<< " Reconstructing " << fieldClassName << "s\n" << endl; - for - ( - IOobjectList::const_iterator fieldIter = fields.begin(); - fieldIter != fields.end(); - ++fieldIter - ) + forAllConstIter(IOobjectList, fields, fieldIter) { - Info<< " " << fieldIter()->name() << endl; + if + ( + !selectedFields.size() + || selectedFields.found(fieldIter()->name()) + ) + { + Info<< " " << fieldIter()->name() << endl; - reconstructFvVolumeField(*fieldIter())().write(); + reconstructFvVolumeField(*fieldIter())().write(); + } } - Info<< endl; } } -// Reconstruct and write all surface fields +// Reconstruct and write all/selected surface fields template void Foam::fvFieldReconstructor::reconstructFvSurfaceFields ( - const IOobjectList& objects + const IOobjectList& objects, + const HashSet& selectedFields ) { const word& fieldClassName = @@ -500,18 +503,19 @@ void Foam::fvFieldReconstructor::reconstructFvSurfaceFields { Info<< " Reconstructing " << fieldClassName << "s\n" << endl; - for - ( - IOobjectList::const_iterator fieldIter = fields.begin(); - fieldIter != fields.end(); - ++fieldIter - ) + forAllConstIter(IOobjectList, fields, fieldIter) { - Info<< " " << fieldIter()->name() << endl; + if + ( + !selectedFields.size() + || selectedFields.found(fieldIter()->name()) + ) + { + Info<< " " << fieldIter()->name() << endl; - reconstructFvSurfaceField(*fieldIter())().write(); + reconstructFvSurfaceField(*fieldIter())().write(); + } } - Info<< endl; } } diff --git a/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C b/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C index db14efe3e5..50818f4db9 100644 --- a/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C +++ b/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C @@ -48,9 +48,17 @@ int main(int argc, char *argv[]) argList::noParallel(); timeSelector::addOptions(); # include "addRegionOption.H" + argList::validOptions.insert("fields", "\"(list of fields)\""); + # include "setRootCase.H" # include "createTime.H" + HashSet selectedFields; + if (args.options().found("fields")) + { + IStringStream(args.options()["fields"])() >> selectedFields; + } + // determine the processor count directly label nProcs = 0; while (dir(args.path()/(word("processor") + name(nProcs)))) @@ -184,13 +192,37 @@ int main(int argc, char *argv[]) procMeshes.boundaryProcAddressing() ); - fvReconstructor.reconstructFvVolumeFields(objects); - fvReconstructor.reconstructFvVolumeFields(objects); - fvReconstructor.reconstructFvVolumeFields(objects); - fvReconstructor.reconstructFvVolumeFields(objects); - fvReconstructor.reconstructFvVolumeFields(objects); + fvReconstructor.reconstructFvVolumeFields + ( + objects, + selectedFields + ); + fvReconstructor.reconstructFvVolumeFields + ( + objects, + selectedFields + ); + fvReconstructor.reconstructFvVolumeFields + ( + objects, + selectedFields + ); + fvReconstructor.reconstructFvVolumeFields + ( + objects, + selectedFields + ); + fvReconstructor.reconstructFvVolumeFields + ( + objects, + selectedFields + ); - fvReconstructor.reconstructFvSurfaceFields(objects); + fvReconstructor.reconstructFvSurfaceFields + ( + objects, + selectedFields + ); } else { From d1f908df2ca3660a2cc35c9d189f8b7471eac381 Mon Sep 17 00:00:00 2001 From: andy Date: Tue, 30 Sep 2008 16:55:44 +0100 Subject: [PATCH 6/7] correcting y[] in updateCoeffs() --- .../omegaWallFunction/omegaWallFunctionFvPatchScalarField.C | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/turbulenceModels/RAS/incompressible/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C b/src/turbulenceModels/RAS/incompressible/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C index 2fed92c105..84793a2c42 100644 --- a/src/turbulenceModels/RAS/incompressible/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C +++ b/src/turbulenceModels/RAS/incompressible/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C @@ -160,7 +160,7 @@ void omegaWallFunctionFvPatchScalarField::updateCoeffs() scalar yPlus = Cmu25*y[faceI]*sqrt(k[faceCellI])/nuw[faceI]; - omega[faceCellI] = sqrt(k[faceCellI])/(Cmu25*kappa*y[faceCellI]); + omega[faceCellI] = sqrt(k[faceCellI])/(Cmu25*kappa*y[faceI]); if (yPlus > yPlusLam) { @@ -168,7 +168,7 @@ void omegaWallFunctionFvPatchScalarField::updateCoeffs() (nutw[faceI] + nuw[faceI]) *magGradUw[faceI] *Cmu25*sqrt(k[faceCellI]) - /(kappa*y[faceCellI]); + /(kappa*y[faceI]); } else { From 3ac92339f1f8a869ea8e1e8fe7a8917b805494af Mon Sep 17 00:00:00 2001 From: henry Date: Tue, 30 Sep 2008 17:35:58 +0100 Subject: [PATCH 7/7] Changed the SquareMatrix construction from initial size and value and corrected usages of it. --- .../matrices/LUscalarMatrix/LUscalarMatrix.C | 4 ++-- .../matrices/SquareMatrix/SquareMatrix.H | 5 +++-- .../matrices/SquareMatrix/SquareMatrixI.H | 20 ++++++++++++++++--- 3 files changed, 22 insertions(+), 7 deletions(-) diff --git a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C index f425accf72..3e8fc3fe0f 100644 --- a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C +++ b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C @@ -101,7 +101,7 @@ Foam::LUscalarMatrix::LUscalarMatrix nCells += lduMatrices[i].size(); } - scalarSquareMatrix m(nCells, 0.0); + scalarSquareMatrix m(nCells, nCells, 0.0); transfer(m); convert(lduMatrices); } @@ -109,7 +109,7 @@ Foam::LUscalarMatrix::LUscalarMatrix else { label nCells = ldum.lduAddr().size(); - scalarSquareMatrix m(nCells, 0.0); + scalarSquareMatrix m(nCells, nCells, 0.0); transfer(m); convert(ldum, interfaceCoeffs, interfaces); } diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H index 6239781d85..27160dc239 100644 --- a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H +++ b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H @@ -69,9 +69,10 @@ public: // It checks that m == n. inline SquareMatrix(const label m, const label n); - //- Construct with given number of rows/columns + //- Construct with given number of rows and rows // and value for all elements. - inline SquareMatrix(const label n, const Type&); + // It checks that m == n. + inline SquareMatrix(const label m, const label n, const Type&); //- Construct from Istream. inline SquareMatrix(Istream&); diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H index a7371a6a2e..d3ee1cd6bf 100644 --- a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H +++ b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H @@ -53,10 +53,24 @@ inline Foam::SquareMatrix::SquareMatrix(const label m, const label n) } template -inline Foam::SquareMatrix::SquareMatrix(const label n, const Type& t) +inline Foam::SquareMatrix::SquareMatrix +( + const label m, + const label n, + const Type& t +) : - Matrix, Type>(n, t) -{} + Matrix, Type>(m, n, t) +{ + if (m != n) + { + FatalErrorIn + ( + "SquareMatrix::SquareMatrix" + "(const label m, const label n, const Type&)" + ) << "m != n for constructing a square matrix" << exit(FatalError); + } +} template inline Foam::SquareMatrix::SquareMatrix(Istream& is)