Merge branch 'master' of /home/dm4/OpenFOAM/OpenFOAM-dev

This commit is contained in:
mattijs
2011-10-19 18:27:18 +01:00
41 changed files with 501 additions and 379 deletions

View File

@ -66,6 +66,7 @@ Description
#include "Switch.H"
#include "bound.H"
#include "dynamicRefineFvMesh.H"
#include "pimpleControl.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -83,6 +84,8 @@ int main(int argc, char *argv[])
#include "compressibleCourantNo.H"
#include "setInitialDeltaT.H"
pimpleControl pimple(mesh);
scalar StCoNum = 0.0;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -94,7 +97,6 @@ int main(int argc, char *argv[])
while (runTime.run())
{
#include "readTimeControls.H"
#include "readPISOControls.H"
#include "compressibleCourantNo.H"
#include "setDeltaT.H"
@ -163,26 +165,35 @@ int main(int argc, char *argv[])
#include "rhoEqn.H"
#include "UEqn.H"
// --- PISO loop
for (int corr=1; corr<=nCorr; corr++)
// --- Pressure-velocity PIMPLE corrector loop
for (pimple.start(); pimple.loop(); pimple++)
{
#include "bEqn.H"
#include "ftEqn.H"
#include "huEqn.H"
#include "hEqn.H"
#include "UEqn.H"
if (!ign.ignited())
// --- PISO loop
for (int corr=1; corr<=pimple.nCorr(); corr++)
{
hu == h;
#include "bEqn.H"
#include "ftEqn.H"
#include "huEqn.H"
#include "hEqn.H"
if (!ign.ignited())
{
hu == h;
}
#include "pEqn.H"
}
#include "pEqn.H"
if (pimple.turbCorr())
{
turbulence->correct();
}
}
turbulence->correct();
runTime.write();
Info<< "\nExecutionTime = "

View File

@ -203,10 +203,6 @@ case ThirdParty:
switch ("$WM_COMPILER")
case Gcc:
case Gcc++0x:
set gcc_version=gcc-4.4.3
set gmp_version=gmp-5.0.1
set mpfr_version=mpfr-2.4.2
breaksw
case Gcc46:
case Gcc46++0x:
set gcc_version=gcc-4.6.1

View File

@ -222,12 +222,7 @@ fi
case "${foamCompiler}" in
OpenFOAM | ThirdParty)
case "$WM_COMPILER" in
Gcc | Gcc++0x)
gcc_version=gcc-4.4.3
gmp_version=gmp-5.0.1
mpfr_version=mpfr-2.4.2
;;
Gcc46 | Gcc46++0x)
Gcc | Gcc++0x | Gcc46 | Gcc46++0x)
gcc_version=gcc-4.6.1
gmp_version=gmp-5.0.2
mpfr_version=mpfr-3.0.1

View File

@ -164,10 +164,13 @@ Foam::label Foam::objectRegistry::getEvent() const
if (event_ == labelMax)
{
WarningIn("objectRegistry::getEvent() const")
<< "Event counter has overflowed. "
<< "Resetting counter on all dependent objects." << nl
<< "This might cause extra evaluations." << endl;
if (objectRegistry::debug)
{
WarningIn("objectRegistry::getEvent() const")
<< "Event counter has overflowed. "
<< "Resetting counter on all dependent objects." << nl
<< "This might cause extra evaluations." << endl;
}
// Reset event counter
curEvent = 1;

View File

@ -334,6 +334,7 @@ $(snGradSchemes)/snGradScheme/snGradSchemes.C
$(snGradSchemes)/correctedSnGrad/correctedSnGrads.C
$(snGradSchemes)/limitedSnGrad/limitedSnGrads.C
$(snGradSchemes)/uncorrectedSnGrad/uncorrectedSnGrads.C
$(snGradSchemes)/orthogonalSnGrad/orthogonalSnGrads.C
/*
$(snGradSchemes)/quadraticFitSnGrad/quadraticFitSnGradData.C
$(snGradSchemes)/quadraticFitSnGrad/quadraticFitSnGrads.C

View File

@ -110,38 +110,34 @@ void Foam::extendedLeastSquaresVectors::makeLeastSquaresVectors() const
const labelUList& owner = mesh_.owner();
const labelUList& neighbour = mesh_.neighbour();
// Build the d-vectors
surfaceVectorField d
(
mesh_.Sf()/(mesh_.magSf()*mesh_.deltaCoeffs())
);
if (!mesh_.orthogonal())
{
d -= mesh_.correctionVectors()/mesh_.deltaCoeffs();
}
const volVectorField& C = mesh.C();
// Set up temporary storage for the dd tensor (before inversion)
symmTensorField dd(mesh_.nCells(), symmTensor::zero);
forAll(owner, faceI)
forAll(owner, facei)
{
const symmTensor wdd(1.0/magSqr(d[faceI])*sqr(d[faceI]));
label own = owner[facei];
label nei = neighbour[facei];
dd[owner[faceI]] += wdd;
dd[neighbour[faceI]] += wdd;
vector d = C[nei] - C[own];
const symmTensor wdd(1.0/magSqr(d[facei])*sqr(d[facei]));
dd[own] += wdd;
dd[nei] += wdd;
}
// Visit the boundaries. Coupled boundaries are taken into account
// in the construction of d vectors.
forAll(d.boundaryField(), patchI)
forAll(lsP.boundaryField(), patchi)
{
const fvsPatchVectorField& pd = d.boundaryField()[patchI];
const fvPatch& p = pd.patch();
const fvPatch& p = lsP.boundaryField()[patchi].patch();
const labelUList& faceCells = p.faceCells();
// Build the d-vectors
vectorField pd = p.delta();
forAll(pd, patchFaceI)
{
dd[faceCells[patchFaceI]] +=
@ -232,24 +228,30 @@ void Foam::extendedLeastSquaresVectors::makeLeastSquaresVectors() const
// Revisit all faces and calculate the lsP and lsN vectors
forAll(owner, faceI)
forAll(owner, facei)
{
lsP[faceI] =
(1.0/magSqr(d[faceI]))*(invDd[owner[faceI]] & d[faceI]);
label own = owner[facei];
label nei = neighbour[facei];
lsN[faceI] =
((-1.0)/magSqr(d[faceI]))*(invDd[neighbour[faceI]] & d[faceI]);
vector d = C[nei] - C[own];
lsP[facei] =
(1.0/magSqr(d[facei]))*(invDd[owner[facei]] & d);
lsN[facei] =
((-1.0)/magSqr(d[facei]))*(invDd[neighbour[facei]] & d);
}
forAll(lsP.boundaryField(), patchI)
{
const fvsPatchVectorField& pd = d.boundaryField()[patchI];
fvsPatchVectorField& patchLsP = lsP.boundaryField()[patchI];
const fvPatch& p = patchLsP.patch();
const labelUList& faceCells = p.faceCells();
// Build the d-vectors
vectorField pd = p.delta();
forAll(p, patchFaceI)
{
patchLsP[patchFaceI] =

View File

@ -121,24 +121,12 @@ Foam::fv::fourthGrad<Type>::calcGrad
const scalarField& lambdap = lambda.boundaryField()[patchi];
const fvPatch& p = fGrad.boundaryField()[patchi].patch();
const labelUList& faceCells = p.faceCells();
// Build the d-vectors
vectorField pd
(
mesh.Sf().boundaryField()[patchi]
/ (
mesh.magSf().boundaryField()[patchi]
* mesh.deltaCoeffs().boundaryField()[patchi]
)
);
if (!mesh.orthogonal())
{
pd -= mesh.correctionVectors().boundaryField()[patchi]
/mesh.deltaCoeffs().boundaryField()[patchi];
}
const labelUList& faceCells =
fGrad.boundaryField()[patchi].patch().faceCells();
vectorField pd = p.delta();
const Field<GradType> neighbourSecondfGrad
(

View File

@ -130,20 +130,7 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const
const labelUList& faceCells = p.patch().faceCells();
// Build the d-vectors
vectorField pd
(
mesh.Sf().boundaryField()[patchi]
/ (
mesh.magSf().boundaryField()[patchi]
* mesh.deltaCoeffs().boundaryField()[patchi]
)
);
if (!mesh.orthogonal())
{
pd -= mesh.correctionVectors().boundaryField()[patchi]
/mesh.deltaCoeffs().boundaryField()[patchi];
}
vectorField pd = p.delta();
if (p.coupled())
{
@ -196,21 +183,7 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const
const labelUList& faceCells = p.faceCells();
// Build the d-vectors
vectorField pd
(
mesh.Sf().boundaryField()[patchi]
/(
mesh.magSf().boundaryField()[patchi]
*mesh.deltaCoeffs().boundaryField()[patchi]
)
);
if (!mesh.orthogonal())
{
pd -= mesh.correctionVectors().boundaryField()[patchi]
/mesh.deltaCoeffs().boundaryField()[patchi];
}
vectorField pd = p.delta();
if (p.coupled())
{

View File

@ -211,7 +211,7 @@ public:
// Add the patch constructor functions to the hash tables
#define makeFvLaplacianTypeScheme(SS, Type, GType) \
#define makeFvLaplacianTypeScheme(SS, GType, Type) \
\
typedef SS<Type, GType> SS##Type##GType; \
defineNamedTemplateTypeNameAndDebug(SS##Type##GType, 0); \
@ -224,13 +224,20 @@ public:
#define makeFvLaplacianScheme(SS) \
\
makeFvLaplacianTypeScheme(SS, scalar, scalar) \
makeFvLaplacianTypeScheme(SS, scalar, symmTensor) \
makeFvLaplacianTypeScheme(SS, scalar, tensor) \
makeFvLaplacianTypeScheme(SS, vector, scalar) \
makeFvLaplacianTypeScheme(SS, sphericalTensor, scalar) \
makeFvLaplacianTypeScheme(SS, symmTensor, scalar) \
makeFvLaplacianTypeScheme(SS, tensor, scalar) \
makeFvLaplacianTypeScheme(SS, scalar, vector) \
makeFvLaplacianTypeScheme(SS, symmTensor, vector) \
makeFvLaplacianTypeScheme(SS, tensor, vector) \
makeFvLaplacianTypeScheme(SS, scalar, sphericalTensor) \
makeFvLaplacianTypeScheme(SS, symmTensor, sphericalTensor) \
makeFvLaplacianTypeScheme(SS, tensor, sphericalTensor) \
makeFvLaplacianTypeScheme(SS, scalar, symmTensor) \
makeFvLaplacianTypeScheme(SS, symmTensor, symmTensor) \
makeFvLaplacianTypeScheme(SS, tensor, symmTensor) \
makeFvLaplacianTypeScheme(SS, scalar, tensor) \
makeFvLaplacianTypeScheme(SS, symmTensor, tensor) \
makeFvLaplacianTypeScheme(SS, tensor, tensor)
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View File

@ -50,7 +50,7 @@ Foam::fv::correctedSnGrad<Type>::fullGradCorrection
// construct GeometricField<Type, fvsPatchField, surfaceMesh>
tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > tssf =
mesh.correctionVectors()
mesh.nonOrthCorrectionVectors()
& linear<typename outerProduct<vector, Type>::type>(mesh).interpolate
(
gradScheme<Type>::New
@ -88,7 +88,7 @@ Foam::fv::correctedSnGrad<Type>::correction
IOobject::NO_WRITE
),
mesh,
vf.dimensions()*mesh.deltaCoeffs().dimensions()
vf.dimensions()*mesh.nonOrthDeltaCoeffs().dimensions()
)
);
GeometricField<Type, fvsPatchField, surfaceMesh>& ssf = tssf();
@ -98,7 +98,7 @@ Foam::fv::correctedSnGrad<Type>::correction
ssf.replace
(
cmpt,
mesh.correctionVectors()
mesh.nonOrthCorrectionVectors()
& linear
<
typename

View File

@ -96,13 +96,13 @@ public:
const GeometricField<Type, fvPatchField, volMesh>&
) const
{
return this->mesh().deltaCoeffs();
return this->mesh().nonOrthDeltaCoeffs();
}
//- Return true if this scheme uses an explicit correction
virtual bool corrected() const
{
return !this->mesh().orthogonal();
return true;
}
//- Return the explicit correction to the correctedSnGrad

View File

@ -120,13 +120,13 @@ public:
const GeometricField<Type, fvPatchField, volMesh>&
) const
{
return this->mesh().deltaCoeffs();
return this->mesh().nonOrthDeltaCoeffs();
}
//- Return true if this scheme uses an explicit correction
virtual bool corrected() const
{
return !this->mesh().orthogonal();
return true;
}
//- Return the explicit correction to the limitedSnGrad

View File

@ -0,0 +1,76 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ 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 3 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, see <http://www.gnu.org/licenses/>.
Description
Simple central-difference snGrad scheme without non-orthogonal correction.
\*---------------------------------------------------------------------------*/
#include "orthogonalSnGrad.H"
#include "volFields.H"
#include "surfaceFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace fv
{
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class Type>
orthogonalSnGrad<Type>::~orthogonalSnGrad()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
orthogonalSnGrad<Type>::correction
(
const GeometricField<Type, fvPatchField, volMesh>&
) const
{
notImplemented
(
"orthogonalSnGrad<Type>::correction"
"(const GeometricField<Type, fvPatchField, volMesh>&)"
);
return tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >(NULL);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace fv
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View File

@ -0,0 +1,133 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ 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 3 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, see <http://www.gnu.org/licenses/>.
Class
Foam::fv::orthogonalSnGrad
Description
Simple central-difference snGrad scheme without non-orthogonal correction.
SourceFiles
orthogonalSnGrad.C
\*---------------------------------------------------------------------------*/
#ifndef orthogonalSnGrad_H
#define orthogonalSnGrad_H
#include "snGradScheme.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace fv
{
/*---------------------------------------------------------------------------*\
Class orthogonalSnGrad Declaration
\*---------------------------------------------------------------------------*/
template<class Type>
class orthogonalSnGrad
:
public snGradScheme<Type>
{
// Private Member Functions
//- Disallow default bitwise assignment
void operator=(const orthogonalSnGrad&);
public:
//- Runtime type information
TypeName("uncorrected");
// Constructors
//- Construct from mesh
orthogonalSnGrad(const fvMesh& mesh)
:
snGradScheme<Type>(mesh)
{}
//- Construct from mesh and data stream
orthogonalSnGrad(const fvMesh& mesh, Istream&)
:
snGradScheme<Type>(mesh)
{}
//- Destructor
virtual ~orthogonalSnGrad();
// Member Functions
//- Return the interpolation weighting factors for the given field
virtual tmp<surfaceScalarField> deltaCoeffs
(
const GeometricField<Type, fvPatchField, volMesh>&
) const
{
return this->mesh().deltaCoeffs();
}
//- Return true if this scheme uses an explicit correction
virtual bool corrected() const
{
return false;
}
//- Return the explicit correction to the orthogonalSnGrad
// for the given field
virtual tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
correction(const GeometricField<Type, fvPatchField, volMesh>&) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace fv
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "orthogonalSnGrad.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,42 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ 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 3 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, see <http://www.gnu.org/licenses/>.
Description
Simple central-difference snGrad scheme without non-orthogonal correction.
\*---------------------------------------------------------------------------*/
#include "orthogonalSnGrad.H"
#include "fvMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace fv
{
makeSnGradScheme(orthogonalSnGrad)
}
}
// ************************************************************************* //

View File

@ -109,7 +109,7 @@ public:
const GeometricField<Type, fvPatchField, volMesh>&
) const
{
return this->mesh().deltaCoeffs();
return this->mesh().nonOrthDeltaCoeffs();
}
//- Return true if this scheme uses an explicit correction

View File

@ -247,7 +247,7 @@ Foam::label Foam::quadraticFitSnGradData::calcFit
scalarList singVals(minSize_);
label nSVDzeros = 0;
const scalar deltaCoeff = mesh().deltaCoeffs()[faci];
const scalar deltaCoeff = deltaCoeffs()[faci];
bool goodFit = false;
for (int iIt = 0; iIt < 10 && !goodFit; iIt++)

View File

@ -96,7 +96,7 @@ public:
const GeometricField<Type, fvPatchField, volMesh>&
) const
{
return this->mesh().deltaCoeffs();
return this->mesh().nonOrthDeltaCoeffs();
}
//- Return true if this scheme uses an explicit correction

View File

@ -66,9 +66,6 @@ protected:
//- Make patch weighting factors
virtual void makeWeights(scalarField&) const = 0;
//- Make patch face - neighbour cell distances
virtual void makeDeltaCoeffs(scalarField&) const = 0;
public:

View File

@ -57,25 +57,6 @@ void Foam::cyclicFvPatch::makeWeights(scalarField& w) const
}
// Make patch face - neighbour cell distances
void Foam::cyclicFvPatch::makeDeltaCoeffs(scalarField& dc) const
{
//const cyclicPolyPatch& nbrPatch = cyclicPolyPatch_.neighbPatch();
const cyclicFvPatch& nbrPatch = neighbFvPatch();
const scalarField deltas(nf() & fvPatch::delta());
const scalarField nbrDeltas(nbrPatch.nf() & nbrPatch.fvPatch::delta());
forAll(deltas, facei)
{
scalar di = deltas[facei];
scalar dni = nbrDeltas[facei];
dc[facei] = 1.0/(di + dni);
}
}
// Return delta (P to N) vectors across coupled patch
Foam::tmp<Foam::vectorField> Foam::cyclicFvPatch::delta() const
{

View File

@ -66,9 +66,6 @@ protected:
//- Make patch weighting factors
void makeWeights(scalarField&) const;
//- Make patch face - neighbour cell distances
void makeDeltaCoeffs(scalarField&) const;
public:

View File

@ -60,27 +60,6 @@ void Foam::cyclicAMIFvPatch::makeWeights(scalarField& w) const
}
void Foam::cyclicAMIFvPatch::makeDeltaCoeffs(scalarField& dc) const
{
const cyclicAMIFvPatch& nbrPatch = neighbFvPatch();
const scalarField deltas(nf() & fvPatch::delta());
const scalarField nbrDeltas
(
interpolate(nbrPatch.nf() & nbrPatch.fvPatch::delta())
);
forAll(deltas, faceI)
{
scalar di = deltas[faceI];
scalar dni = nbrDeltas[faceI];
dc[faceI] = 1.0/(di + dni);
}
}
Foam::tmp<Foam::vectorField> Foam::cyclicAMIFvPatch::delta() const
{
const vectorField patchD(fvPatch::delta());

View File

@ -66,9 +66,6 @@ protected:
//- Make patch weighting factors
void makeWeights(scalarField&) const;
//- Make patch face - neighbour cell distances
void makeDeltaCoeffs(scalarField&) const;
public:

View File

@ -65,19 +65,6 @@ void processorFvPatch::makeWeights(scalarField& w) const
}
void processorFvPatch::makeDeltaCoeffs(scalarField& dc) const
{
if (Pstream::parRun())
{
dc = (1.0 - weights())/(nf() & fvPatch::delta());
}
else
{
dc = 1.0/(nf() & fvPatch::delta());
}
}
tmp<vectorField> processorFvPatch::delta() const
{
if (Pstream::parRun())

View File

@ -65,9 +65,6 @@ protected:
//- Make patch weighting factors
void makeWeights(scalarField&) const;
//- Make patch face - neighbour cell distances
void makeDeltaCoeffs(scalarField&) const;
public:

View File

@ -150,12 +150,6 @@ void Foam::fvPatch::makeWeights(scalarField& w) const
}
void Foam::fvPatch::makeDeltaCoeffs(scalarField& dc) const
{
dc = 1.0/(nf() & delta());
}
void Foam::fvPatch::initMovePoints()
{}

View File

@ -86,9 +86,6 @@ protected:
//- Make patch weighting factors
virtual void makeWeights(scalarField&) const;
//- Make patch face - neighbour cell distances
virtual void makeDeltaCoeffs(scalarField&) const;
//- Initialise the patches for moving points
virtual void initMovePoints();

View File

@ -122,20 +122,7 @@ Foam::LimitedScheme<Type, Limiter, LimitFunc>::limiter
);
// Build the d-vectors
vectorField pd
(
mesh.Sf().boundaryField()[patchi]
/ (
mesh.magSf().boundaryField()[patchi]
* mesh.deltaCoeffs().boundaryField()[patchi]
)
);
if (!mesh.orthogonal())
{
pd -= mesh.correctionVectors().boundaryField()[patchi]
/mesh.deltaCoeffs().boundaryField()[patchi];
}
vectorField pd = CDweights.boundaryField()[patchi].patch().delta();
forAll(pLim, face)
{

View File

@ -111,20 +111,7 @@ Foam::linearUpwind<Type>::correction
);
// Build the d-vectors
vectorField pd
(
mesh.Sf().boundaryField()[patchi]
/ (
mesh.magSf().boundaryField()[patchi]
* mesh.deltaCoeffs().boundaryField()[patchi]
)
);
if (!mesh.orthogonal())
{
pd -= mesh.correctionVectors().boundaryField()[patchi]
/mesh.deltaCoeffs().boundaryField()[patchi];
}
vectorField pd = Cf.boundaryField()[patchi].patch().delta();
forAll(pOwner, facei)
{

View File

@ -82,21 +82,17 @@ void Foam::skewCorrectionVectors::makeSkewCorrectionVectors() const
const surfaceVectorField& Sf = mesh_.Sf();
const labelUList& owner = mesh_.owner();
const labelUList& neighbour = mesh_.neighbour();
// Build the d-vectors
surfaceVectorField d(Sf/(mesh_.magSf()*mesh_.deltaCoeffs()));
if (!mesh_.orthogonal())
forAll(owner, facei)
{
d -= mesh_.correctionVectors()/mesh_.deltaCoeffs();
}
label own = owner[facei];
label nei = neighbour[facei];
forAll(owner, faceI)
{
vector Cpf = Cf[faceI] - C[owner[faceI]];
vector d = C[nei] - C[own];
vector Cpf = Cf[facei] - C[own];
SkewCorrVecs[faceI] =
Cpf - ((Sf[faceI] & Cpf)/(Sf[faceI] & d[faceI]))*d[faceI];
SkewCorrVecs[facei] = Cpf - ((Sf[facei] & Cpf)/(Sf[facei] & d))*d;
}
@ -115,7 +111,7 @@ void Foam::skewCorrectionVectors::makeSkewCorrectionVectors() const
const labelUList& faceCells = p.faceCells();
const vectorField& patchFaceCentres = Cf.boundaryField()[patchI];
const vectorField& patchSf = Sf.boundaryField()[patchI];
const vectorField& patchD = d.boundaryField()[patchI];
const vectorField patchD = p.delta();
forAll(p, patchFaceI)
{
@ -136,7 +132,7 @@ void Foam::skewCorrectionVectors::makeSkewCorrectionVectors() const
if (Sf.internalField().size())
{
skewCoeff = max(mag(SkewCorrVecs)/mag(d)).value();
skewCoeff = max(mag(SkewCorrVecs)*mesh_.deltaCoeffs()).value();
}
if (debug)
@ -182,7 +178,7 @@ const Foam::surfaceVectorField& Foam::skewCorrectionVectors::operator()() const
if (!skew())
{
FatalErrorIn("skewCorrectionVectors::operator()()")
<< "Cannot return correctionVectors; mesh is not skewed"
<< "Cannot return skewCorrectionVectors; mesh is not skewed"
<< abort(FatalError);
}

View File

@ -31,43 +31,38 @@ Description
#include "surfaceFields.H"
#include "demandDrivenData.H"
#include "coupledFvPatch.H"
#include "unitConversion.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(surfaceInterpolation, 0);
defineTypeNameAndDebug(Foam::surfaceInterpolation, 0);
// * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * //
void surfaceInterpolation::clearOut()
void Foam::surfaceInterpolation::clearOut()
{
deleteDemandDrivenData(weightingFactors_);
deleteDemandDrivenData(differenceFactors_);
deleteDemandDrivenData(correctionVectors_);
deleteDemandDrivenData(weights_);
deleteDemandDrivenData(deltaCoeffs_);
deleteDemandDrivenData(nonOrthDeltaCoeffs_);
deleteDemandDrivenData(nonOrthCorrectionVectors_);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
surfaceInterpolation::surfaceInterpolation(const fvMesh& fvm)
Foam::surfaceInterpolation::surfaceInterpolation(const fvMesh& fvm)
:
mesh_(fvm),
weightingFactors_(NULL),
differenceFactors_(NULL),
orthogonal_(false),
correctionVectors_(NULL)
weights_(NULL),
deltaCoeffs_(NULL),
nonOrthDeltaCoeffs_(NULL),
nonOrthCorrectionVectors_(NULL)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
surfaceInterpolation::~surfaceInterpolation()
Foam::surfaceInterpolation::~surfaceInterpolation()
{
clearOut();
}
@ -75,66 +70,67 @@ surfaceInterpolation::~surfaceInterpolation()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
const surfaceScalarField& surfaceInterpolation::weights() const
const Foam::surfaceScalarField&
Foam::surfaceInterpolation::weights() const
{
if (!weightingFactors_)
if (!weights_)
{
makeWeights();
}
return (*weightingFactors_);
return (*weights_);
}
const surfaceScalarField& surfaceInterpolation::deltaCoeffs() const
const Foam::surfaceScalarField&
Foam::surfaceInterpolation::deltaCoeffs() const
{
if (!differenceFactors_)
if (!deltaCoeffs_)
{
makeDeltaCoeffs();
}
return (*differenceFactors_);
return (*deltaCoeffs_);
}
bool surfaceInterpolation::orthogonal() const
const Foam::surfaceScalarField&
Foam::surfaceInterpolation::nonOrthDeltaCoeffs() const
{
if (orthogonal_ == false && !correctionVectors_)
if (!nonOrthDeltaCoeffs_)
{
makeCorrectionVectors();
makeNonOrthDeltaCoeffs();
}
return orthogonal_;
return (*nonOrthDeltaCoeffs_);
}
const surfaceVectorField& surfaceInterpolation::correctionVectors() const
const Foam::surfaceVectorField&
Foam::surfaceInterpolation::nonOrthCorrectionVectors() const
{
if (orthogonal())
if (!nonOrthCorrectionVectors_)
{
FatalErrorIn("surfaceInterpolation::correctionVectors()")
<< "cannot return correctionVectors; mesh is orthogonal"
<< abort(FatalError);
makeNonOrthCorrectionVectors();
}
return (*correctionVectors_);
return (*nonOrthCorrectionVectors_);
}
// Do what is neccessary if the mesh has moved
bool surfaceInterpolation::movePoints()
bool Foam::surfaceInterpolation::movePoints()
{
deleteDemandDrivenData(weightingFactors_);
deleteDemandDrivenData(differenceFactors_);
orthogonal_ = false;
deleteDemandDrivenData(correctionVectors_);
deleteDemandDrivenData(weights_);
deleteDemandDrivenData(deltaCoeffs_);
deleteDemandDrivenData(nonOrthDeltaCoeffs_);
deleteDemandDrivenData(nonOrthCorrectionVectors_);
return true;
}
void surfaceInterpolation::makeWeights() const
void Foam::surfaceInterpolation::makeWeights() const
{
if (debug)
{
@ -143,20 +139,18 @@ void surfaceInterpolation::makeWeights() const
<< endl;
}
weightingFactors_ = new surfaceScalarField
weights_ = new surfaceScalarField
(
IOobject
(
"weightingFactors",
"weights",
mesh_.pointsInstance(),
mesh_
),
mesh_,
dimless
);
surfaceScalarField& weightingFactors = *weightingFactors_;
surfaceScalarField& weights = *weights_;
// Set local references to mesh data
// (note that we should not use fvMesh sliced fields at this point yet
@ -170,7 +164,7 @@ void surfaceInterpolation::makeWeights() const
const vectorField& Sf = mesh_.faceAreas();
// ... and reference to the internal field of the weighting factors
scalarField& w = weightingFactors.internalField();
scalarField& w = weights.internalField();
forAll(owner, facei)
{
@ -188,11 +182,10 @@ void surfaceInterpolation::makeWeights() const
{
mesh_.boundary()[patchi].makeWeights
(
weightingFactors.boundaryField()[patchi]
weights.boundaryField()[patchi]
);
}
if (debug)
{
Info<< "surfaceInterpolation::makeWeights() : "
@ -202,7 +195,7 @@ void surfaceInterpolation::makeWeights() const
}
void surfaceInterpolation::makeDeltaCoeffs() const
void Foam::surfaceInterpolation::makeDeltaCoeffs() const
{
if (debug)
{
@ -215,18 +208,63 @@ void surfaceInterpolation::makeDeltaCoeffs() const
// needed to make sure deltaCoeffs are calculated for parallel runs.
weights();
differenceFactors_ = new surfaceScalarField
deltaCoeffs_ = new surfaceScalarField
(
IOobject
(
"differenceFactors_",
"deltaCoeffs",
mesh_.pointsInstance(),
mesh_
),
mesh_,
dimless/dimLength
);
surfaceScalarField& DeltaCoeffs = *differenceFactors_;
surfaceScalarField& DeltaCoeffs = *deltaCoeffs_;
// Set local references to mesh data
const volVectorField& C = mesh_.C();
const labelUList& owner = mesh_.owner();
const labelUList& neighbour = mesh_.neighbour();
forAll(owner, facei)
{
DeltaCoeffs[facei] = 1.0/mag(C[neighbour[facei]] - C[owner[facei]]);
}
forAll(DeltaCoeffs.boundaryField(), patchi)
{
DeltaCoeffs.boundaryField()[patchi] =
1.0/mag(mesh_.boundary()[patchi].delta());
}
}
void Foam::surfaceInterpolation::makeNonOrthDeltaCoeffs() const
{
if (debug)
{
Info<< "surfaceInterpolation::makeNonOrthDeltaCoeffs() : "
<< "Constructing differencing factors array for face gradient"
<< endl;
}
// Force the construction of the weighting factors
// needed to make sure deltaCoeffs are calculated for parallel runs.
weights();
nonOrthDeltaCoeffs_ = new surfaceScalarField
(
IOobject
(
"nonOrthDeltaCoeffs",
mesh_.pointsInstance(),
mesh_
),
mesh_,
dimless/dimLength
);
surfaceScalarField& nonOrthDeltaCoeffs = *nonOrthDeltaCoeffs_;
// Set local references to mesh data
@ -242,49 +280,49 @@ void surfaceInterpolation::makeDeltaCoeffs() const
vector unitArea = Sf[facei]/magSf[facei];
// Standard cell-centre distance form
//DeltaCoeffs[facei] = (unitArea & delta)/magSqr(delta);
//NonOrthDeltaCoeffs[facei] = (unitArea & delta)/magSqr(delta);
// Slightly under-relaxed form
//DeltaCoeffs[facei] = 1.0/mag(delta);
//NonOrthDeltaCoeffs[facei] = 1.0/mag(delta);
// More under-relaxed form
//DeltaCoeffs[facei] = 1.0/(mag(unitArea & delta) + VSMALL);
//NonOrthDeltaCoeffs[facei] = 1.0/(mag(unitArea & delta) + VSMALL);
// Stabilised form for bad meshes
DeltaCoeffs[facei] = 1.0/max(unitArea & delta, 0.05*mag(delta));
nonOrthDeltaCoeffs[facei] = 1.0/max(unitArea & delta, 0.05*mag(delta));
}
forAll(DeltaCoeffs.boundaryField(), patchi)
forAll(nonOrthDeltaCoeffs.boundaryField(), patchi)
{
mesh_.boundary()[patchi].makeDeltaCoeffs
(
DeltaCoeffs.boundaryField()[patchi]
);
vectorField delta = mesh_.boundary()[patchi].delta();
nonOrthDeltaCoeffs.boundaryField()[patchi] =
1.0/max(mesh_.boundary()[patchi].nf() & delta, 0.05*mag(delta));
}
}
void surfaceInterpolation::makeCorrectionVectors() const
void Foam::surfaceInterpolation::makeNonOrthCorrectionVectors() const
{
if (debug)
{
Info<< "surfaceInterpolation::makeCorrectionVectors() : "
Info<< "surfaceInterpolation::makeNonOrthCorrectionVectors() : "
<< "Constructing non-orthogonal correction vectors"
<< endl;
}
correctionVectors_ = new surfaceVectorField
nonOrthCorrectionVectors_ = new surfaceVectorField
(
IOobject
(
"correctionVectors",
"nonOrthCorrectionVectors",
mesh_.pointsInstance(),
mesh_
),
mesh_,
dimless
);
surfaceVectorField& corrVecs = *correctionVectors_;
surfaceVectorField& corrVecs = *nonOrthCorrectionVectors_;
// Set local references to mesh data
const volVectorField& C = mesh_.C();
@ -292,14 +330,14 @@ void surfaceInterpolation::makeCorrectionVectors() const
const labelUList& neighbour = mesh_.neighbour();
const surfaceVectorField& Sf = mesh_.Sf();
const surfaceScalarField& magSf = mesh_.magSf();
const surfaceScalarField& DeltaCoeffs = deltaCoeffs();
const surfaceScalarField& NonOrthDeltaCoeffs = nonOrthDeltaCoeffs();
forAll(owner, facei)
{
vector unitArea = Sf[facei]/magSf[facei];
vector delta = C[neighbour[facei]] - C[owner[facei]];
corrVecs[facei] = unitArea - delta*DeltaCoeffs[facei];
corrVecs[facei] = unitArea - delta*NonOrthDeltaCoeffs[facei];
}
// Boundary correction vectors set to zero for boundary patches
@ -308,18 +346,18 @@ void surfaceInterpolation::makeCorrectionVectors() const
forAll(corrVecs.boundaryField(), patchi)
{
fvsPatchVectorField& patchcorrVecs = corrVecs.boundaryField()[patchi];
fvsPatchVectorField& patchCorrVecs = corrVecs.boundaryField()[patchi];
if (!patchcorrVecs.coupled())
if (!patchCorrVecs.coupled())
{
patchcorrVecs = vector::zero;
patchCorrVecs = vector::zero;
}
else
{
const fvsPatchScalarField& patchDeltaCoeffs
= DeltaCoeffs.boundaryField()[patchi];
const fvsPatchScalarField& patchNonOrthDeltaCoeffs
= NonOrthDeltaCoeffs.boundaryField()[patchi];
const fvPatch& p = patchcorrVecs.patch();
const fvPatch& p = patchCorrVecs.patch();
const vectorField patchDeltas(mesh_.boundary()[patchi].delta());
@ -331,60 +369,19 @@ void surfaceInterpolation::makeCorrectionVectors() const
const vector& delta = patchDeltas[patchFacei];
patchcorrVecs[patchFacei] =
unitArea - delta*patchDeltaCoeffs[patchFacei];
patchCorrVecs[patchFacei] =
unitArea - delta*patchNonOrthDeltaCoeffs[patchFacei];
}
}
}
scalar NonOrthogCoeff = 0.0;
// Calculate the non-orthogonality for meshes with 1 face or more
if (returnReduce(magSf.size(), sumOp<label>()) > 0)
{
NonOrthogCoeff = radToDeg
(
asin
(
min
(
(sum(magSf*mag(corrVecs))/sum(magSf)).value(),
1.0
)
)
);
}
if (debug)
{
Info<< "surfaceInterpolation::makeCorrectionVectors() : "
<< "non-orthogonality coefficient = " << NonOrthogCoeff << " deg."
<< endl;
}
//NonOrthogCoeff = 0.0;
if (NonOrthogCoeff < 0.1)
{
orthogonal_ = true;
deleteDemandDrivenData(correctionVectors_);
}
else
{
orthogonal_ = false;
}
if (debug)
{
Info<< "surfaceInterpolation::makeCorrectionVectors() : "
Info<< "surfaceInterpolation::makeNonOrthCorrectionVectors() : "
<< "Finished constructing non-orthogonal correction vectors"
<< endl;
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View File

@ -59,17 +59,17 @@ class surfaceInterpolation
// Demand-driven data
//- Central-differencing weighting factors
mutable surfaceScalarField* weightingFactors_;
//- Linear difference weighting factors
mutable surfaceScalarField* weights_;
//- Face-gradient difference factors
mutable surfaceScalarField* differenceFactors_;
//- Cell-centre difference coefficients
mutable surfaceScalarField* deltaCoeffs_;
//- Is mesh orthogonal
mutable bool orthogonal_;
//- Non-orthogonal cell-centre difference coefficients
mutable surfaceScalarField* nonOrthDeltaCoeffs_;
//- Non-orthogonality correction vectors
mutable surfaceVectorField* correctionVectors_;
mutable surfaceVectorField* nonOrthCorrectionVectors_;
// Private Member Functions
@ -80,8 +80,11 @@ class surfaceInterpolation
//- Construct face-gradient difference factors
void makeDeltaCoeffs() const;
//- Construct face-gradient difference factors
void makeNonOrthDeltaCoeffs() const;
//- Construct non-orthogonality correction vectors
void makeCorrectionVectors() const;
void makeNonOrthCorrectionVectors() const;
protected:
@ -112,17 +115,18 @@ public:
// Member functions
//- Return reference to weighting factors array
//- Return reference to linear difference weighting factors
const surfaceScalarField& weights() const;
//- Return reference to difference factors array
//- Return reference to cell-centre difference coefficients
const surfaceScalarField& deltaCoeffs() const;
//- Return whether mesh is orthogonal or not
bool orthogonal() const;
//- Return reference to non-orthogonal cell-centre difference
// coefficients
const surfaceScalarField& nonOrthDeltaCoeffs() const;
//- Return reference to non-orthogonality correction vectors array
const surfaceVectorField& correctionVectors() const;
//- Return reference to non-orthogonality correction vectors
const surfaceVectorField& nonOrthCorrectionVectors() const;
//- Do what is neccessary if the mesh has moved
bool movePoints();

View File

@ -86,7 +86,7 @@ atmBoundaryLayerInletEpsilonFvPatchScalarField
z_(dict.lookup("z")),
z0_(readScalar(dict.lookup("z0"))),
kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
zGround_(readScalar(dict.lookup("zGround")))
zGround_("zGround", dict, p.size())
{
if (mag(z_) < SMALL)
{

View File

@ -90,7 +90,7 @@ class atmBoundaryLayerInletEpsilonFvPatchScalarField
const scalar kappa_;
//- Minimum corrdinate value in z direction
const scalar zGround_;
const scalarField zGround_;
public:

View File

@ -115,7 +115,7 @@ atmBoundaryLayerInletVelocityFvPatchVectorField
n_ /= mag(n_);
z_ /= mag(z_);
Ustar_ = kappa_*Uref_/(log((Href_ + z0_)/min(z0_ , 0.001)));
Ustar_ = kappa_*Uref_/(log((Href_ + z0_)/max(z0_ , 0.001)));
evaluate();
}
@ -150,9 +150,11 @@ void atmBoundaryLayerInletVelocityFvPatchVectorField::updateCoeffs()
forAll(coord, i)
{
if ((coord[i] - zGround_) < Href_)
if ((coord[i] - zGround_[i]) < Href_)
{
Un[i] = (Ustar_/kappa_)*log((coord[i] - zGround_ + z0_)/z0_);
Un[i] =
(Ustar_/kappa_)
* log((coord[i] - zGround_[i] + z0_)/max(z0_, 0.001));
}
else
{
@ -181,8 +183,7 @@ void atmBoundaryLayerInletVelocityFvPatchVectorField::write(Ostream& os) const
<< Uref_ << token::END_STATEMENT << nl;
os.writeKeyword("Href")
<< Href_ << token::END_STATEMENT << nl;
os.writeKeyword("zGround")
<< zGround_ << token::END_STATEMENT << nl;
zGround_.writeEntry("zGround", os) ;
writeEntry("value", os);
}

View File

@ -113,7 +113,7 @@ class atmBoundaryLayerInletVelocityFvPatchVectorField
const scalar Href_;
//- Minimum corrdinate value in z direction
const scalar zGround_;
const scalarField zGround_;
public:

View File

@ -190,7 +190,7 @@ relaxationFactors
h 0.95;
k 0.9;
epsilon 0.9;
}w
}
}
relaxationFactors0

View File

@ -13,5 +13,5 @@ z0 0.1;
turbulentKE 1.3;
windDirection (1 0 0);
zDirection (0 0 1);
zGround 935.0;
zGround uniform 935.0;
// ************************************************************************* //

View File

@ -94,6 +94,7 @@ PIMPLE
nCorrectors 2;
nNonOrthogonalCorrectors 0;
nAlphaCorr 2;
correctAlpha no;
}
relaxationFactors

View File

@ -1,3 +1,2 @@
c++DBUG =
#c++OPT = -xSSE3 -O3 -no-prec-div
c++OPT = -xSSE3 -O1 -no-prec-div
c++OPT = -xSSE3 -O2 -no-prec-div

View File

@ -1,5 +1,2 @@
c++DBUG =
#c++OPT = -O3 -xP -no-prec-div
c++OPT = -ansi-alias -O3 -ftz -fno-alias \
-fargument-noalias-global \
-unroll0
c++OPT = -xSSE3 -O2 -no-prec-div