Immersed boundary method preparation

This commit is contained in:
Hrvoje Jasak 2012-04-25 10:22:35 +01:00
parent a3712d23f9
commit 913282318a
46 changed files with 231 additions and 92 deletions

View file

@ -38,7 +38,6 @@ Description
int main(int argc, char *argv[])
{
argList::validOptions.insert("writep", "");
# include "setRootCase.H"

View file

@ -299,7 +299,9 @@ void PDRkEpsilon::correct()
epsEqn().relax();
epsEqn().boundaryManipulate(epsilon_.boundaryField());
// No longer needed: matrix completes at the point of solution
// HJ, 17/Apr/2012
// epsEqn().completeAssembly();
solve(epsEqn);
bound(epsilon_, epsilon0_);

View file

@ -186,6 +186,7 @@ schemeSymbolList ({schemeSymbolListElement}+{space})
starStar ("**")
text ({space}({word}*{space})*)
anythingInBlock ([^)]*)
gridgenComment (({space}|{cspace})({word}*{space})*)
dateDDMMYYYY ({digit}{digit}"/"{digit}{digit}"/"{digit}{digit}{digit}{digit})
dateDDMonYYYY ((({digit}{digit}{space})|({digit}{space})){alpha}*{space}{digit}{digit}{digit}{digit})
@ -728,6 +729,9 @@ endOfSection {space}")"{space}
<ignoreBlock,ignoreEmbeddedBlock>{space}{text} {
}
<ignoreBlock,ignoreEmbeddedBlock>{gridgenComment} {
}
/* ------ Count newlines. ------ */

View file

@ -132,7 +132,7 @@ quote \"
dash "-"
dotColonDash [.:-]
schemeSpecialInitial [!$%&*/:<=>?~_^#.]
schemeSpecialInitial [!$%&*/\\:<=>?~_^#.@']
schemeSpecialSubsequent [.+-]
schemeSymbol (({some_space}|{alpha}|{quote}|{schemeSpecialInitial})({alpha}|{quote}|{digit}|{schemeSpecialInitial}|{schemeSpecialSubsequent})*)
@ -142,7 +142,7 @@ integer {decDigit}+
label [1-9]{decDigit}*
hexLabel {hexDigit}+
zeroLabel {digit}*
signedInteger [-+]?{integer}
word ({alpha}|{digit}|{dotColonDash})*
exponent_part [eE][-+]?{digit}+

View file

@ -72,6 +72,12 @@ public:
{}
// Destructor
virtual ~emptyPointPatch()
{}
// Member Functions
//- Accumulate the effect of constraint direction of this patch
@ -82,10 +88,6 @@ public:
) const;
// Destructor
virtual ~emptyPointPatch()
{}
};

View file

@ -191,7 +191,6 @@ public:
{
return polyPatch_.pointNormals();
}
};

View file

@ -258,26 +258,13 @@ bool Foam::primitiveMesh::checkClosedCells
// Calculate the aspect ration as the maximum of Cartesian component
// aspect ratio to the total area hydraulic area aspect ratio
scalar minCmpt = VGREAT;
scalar maxCmpt = -VGREAT;
for (direction dir = 0; dir < vector::nComponents; dir++)
{
if (meshD[dir] == 1)
{
minCmpt = min(minCmpt, sumMagClosed[cellI][dir]);
maxCmpt = max(maxCmpt, sumMagClosed[cellI][dir]);
}
}
scalar aspectRatio = maxCmpt/(minCmpt + VSMALL);
if (nDims == 3)
{
aspectRatio = max
scalar aspectRatio = max
(
aspectRatio,
1.0/6.0*cmptSum(sumMagClosed[cellI])/pow(vols[cellI], 2.0/3.0)
cmptMax(sumMagClosed[cellI])
/(cmptMin(sumMagClosed[cellI]) + VSMALL),
1.0/6.0*cmptSum(sumMagClosed[cellI])/
Foam::pow(Foam::max(vols[cellI], SMALL), 2.0/3.0)
);
}
maxAspectRatio = max(maxAspectRatio, aspectRatio);

View file

@ -45,7 +45,7 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
Class emptyFvPatch Declaration
Class emptyFvPatchField Declaration
\*---------------------------------------------------------------------------*/
template<class Type>

View file

@ -166,6 +166,18 @@ void fvsPatchField<Type>::check(const fvsPatchField<Type>& ptf) const
}
template<class Type>
void fvsPatchField<Type>::check(const fvPatchField<Type>& ptf) const
{
if (&patch_ != &(ptf.patch()))
{
FatalErrorIn("PatchField<Type>::check(const fvsPatchField<Type>&)")
<< "different patches for fvsPatchField<Type>s"
<< abort(FatalError);
}
}
// Map from self
template<class Type>
void fvsPatchField<Type>::autoMap
@ -280,6 +292,17 @@ void fvsPatchField<Type>::operator/=
}
template<class Type>
void fvsPatchField<Type>::operator=
(
const fvPatchField<Type>& ptf
)
{
check(ptf);
Field<Type>::operator=(ptf);
}
template<class Type>
void fvsPatchField<Type>::operator+=
(

View file

@ -60,6 +60,8 @@ class dictionary;
class fvPatchFieldMapper;
class surfaceMesh;
template<class Type>
class fvPatchField;
// Forward declaration of friend functions and operators
@ -313,6 +315,23 @@ public:
);
// Evaluation functions
//- Initialise the evaluation of the patch field
virtual void initEvaluate
(
const Pstream::commsTypes commsType=Pstream::blocking
)
{}
//- Evaluate the patch field, sets Updated to false
virtual void evaluate
(
const Pstream::commsTypes commsType=Pstream::blocking
)
{}
//- Write
virtual void write(Ostream&) const;
@ -322,6 +341,9 @@ public:
//- Check fvsPatchField<Type> against given fvsPatchField<Type>
void check(const fvsPatchField<Type>&) const;
//- Check fvsPatchField<Type> against given fvPatchField<Type>
void check(const fvPatchField<Type>&) const;
// Member operators
@ -333,6 +355,8 @@ public:
virtual void operator*=(const fvsPatchField<scalar>&);
virtual void operator/=(const fvsPatchField<scalar>&);
virtual void operator=(const fvPatchField<Type>&);
virtual void operator+=(const Field<Type>&);
virtual void operator-=(const Field<Type>&);

View file

@ -193,6 +193,7 @@ Foam::fvMatrix<Type>::fvMatrix
source_(psi.size(), pTraits<Type>::zero),
internalCoeffs_(psi.mesh().boundary().size()),
boundaryCoeffs_(psi.mesh().boundary().size()),
assemblyCompleted_(false),
faceFluxCorrectionPtr_(NULL)
{
if (debug)
@ -241,6 +242,7 @@ Foam::fvMatrix<Type>::fvMatrix(const fvMatrix<Type>& fvm)
source_(fvm.source_),
internalCoeffs_(fvm.internalCoeffs_),
boundaryCoeffs_(fvm.boundaryCoeffs_),
assemblyCompleted_(fvm.assemblyCompleted_),
faceFluxCorrectionPtr_(NULL)
{
if (debug)
@ -288,6 +290,7 @@ Foam::fvMatrix<Type>::fvMatrix(const tmp<fvMatrix<Type> >& tfvm)
const_cast<fvMatrix<Type>&>(tfvm()).boundaryCoeffs_,
tfvm.isTmp()
),
assemblyCompleted_(tfvm().assemblyCompleted()),
faceFluxCorrectionPtr_(NULL)
{
if (debug)
@ -332,6 +335,7 @@ Foam::fvMatrix<Type>::fvMatrix
source_(is),
internalCoeffs_(psi.mesh().boundary().size()),
boundaryCoeffs_(psi.mesh().boundary().size()),
assemblyCompleted_(false),
faceFluxCorrectionPtr_(NULL)
{
if (debug)
@ -614,12 +618,25 @@ void Foam::fvMatrix<Type>::relax()
template<class Type>
void Foam::fvMatrix<Type>::boundaryManipulate
(
typename GeometricField<Type, fvPatchField, volMesh>::
GeometricBoundaryField& bFields
)
void Foam::fvMatrix<Type>::completeAssembly()
{
if (assemblyCompleted_)
{
return;
}
if (debug)
{
InfoIn("void Foam::fvMatrix<Type>::completeAssembly()")
<< "Completing matrix for equation " << this->psi().name()
<< endl;
}
assemblyCompleted_ = true;
typename GeometricField<Type, fvPatchField, volMesh>::
GeometricBoundaryField& bFields = psi_.boundaryField();
forAll(bFields, patchI)
{
bFields[patchI].manipulateMatrix(*this);

View file

@ -132,6 +132,8 @@ class fvMatrix
// for boundary cells
FieldField<Field, Type> boundaryCoeffs_;
//- Has assembly been completed?
bool assemblyCompleted_;
//- Face flux field for non-orthogonal correction
mutable GeometricField<Type, fvsPatchField, surfaceMesh>
@ -290,6 +292,12 @@ public:
return boundaryCoeffs_;
}
//- Return true if matrix assemble has been completed
bool assemblyCompleted() const
{
return assemblyCompleted_;
}
//- Declare return type of the faceFluxCorrectionPtr() function
typedef GeometricField<Type, fvsPatchField, surfaceMesh>
@ -315,7 +323,7 @@ public:
void addBoundarySource
(
Field<Type>& source,
const bool couples=true
const bool couples = true
) const;
@ -357,12 +365,9 @@ public:
// alpha is read from controlDict
void relax();
//- Manipulate based on a boundary field
void boundaryManipulate
(
typename GeometricField<Type, fvPatchField, volMesh>::
GeometricBoundaryField& values
);
//- Complete matrix assembly for solution:
// Manipulate based on a boundary field
void completeAssembly();
//- Construct and return the solver
// Use the given solver controls

View file

@ -63,6 +63,9 @@ Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::solve
<< endl;
}
// Complete matrix assembly. HJ, 17/Apr/2012
this->completeAssembly();
lduSolverPerformance solverPerfVec
(
"fvMatrix<Type>::solve",

View file

@ -104,6 +104,9 @@ Foam::fvMatrix<Foam::scalar>::fvSolver::solve
const dictionary& solverControls
)
{
// Complete matrix assembly. HJ, 17/Apr/2012
fvMat_.completeAssembly();
scalarField saveDiag = fvMat_.diag();
fvMat_.addBoundaryDiag(fvMat_.diag(), 0);
@ -138,6 +141,9 @@ Foam::lduMatrix::solverPerformance Foam::fvMatrix<Foam::scalar>::solve
<< endl;
}
// Complete matrix assembly. HJ, 17/Apr/2012
completeAssembly();
scalarField saveDiag = diag();
addBoundaryDiag(diag(), 0);

View file

@ -55,7 +55,7 @@ class emptyFvPatch
{
// Private data
//- face-cell addressing
//- Dummy face-cell addressing
const labelList::subList faceCells_;
@ -80,7 +80,7 @@ public:
return 0;
}
//- Return faceCells
//- Return faceCells of zero size
virtual const unallocLabelList& faceCells() const;
};

View file

@ -56,6 +56,7 @@ class processorFvPatch
{
// Private Data
//- Reference to processor patch
const processorPolyPatch& procPolyPatch_;

View file

@ -87,6 +87,7 @@ private:
//- Is surface closed
mutable label surfaceClosed_;
// Private Member Functions
////- Helper: find instance of files without header

View file

@ -99,7 +99,7 @@ public:
// Destructor
~twoPhaseMixture()
virtual ~twoPhaseMixture()
{}

View file

@ -91,10 +91,10 @@ public:
// Member Functions
//- Return the laminar viscosity
const volScalarField& nu() const;
virtual const volScalarField& nu() const;
//- Correct the laminar viscosity
void correct();
virtual void correct();
//- Read transportProperties dictionary
virtual bool read();

View file

@ -93,26 +93,26 @@ public:
// Destructor
~BirdCarreau()
virtual ~BirdCarreau()
{}
// Member Functions
//- Return the laminar viscosity
const volScalarField& nu() const
virtual const volScalarField& nu() const
{
return nu_;
}
//- Correct the laminar viscosity
void correct()
virtual void correct()
{
nu_ = calcNu();
}
//- Read transportProperties dictionary
bool read(const dictionary& viscosityProperties);
virtual bool read(const dictionary& viscosityProperties);
};

View file

@ -92,26 +92,26 @@ public:
// Destructor
~CrossPowerLaw()
virtual ~CrossPowerLaw()
{}
// Member Functions
//- Return the laminar viscosity
const volScalarField& nu() const
virtual const volScalarField& nu() const
{
return nu_;
}
//- Correct the laminar viscosity
void correct()
virtual void correct()
{
nu_ = calcNu();
}
//- Read transportProperties dictionary
bool read(const dictionary& viscosityProperties);
virtual bool read(const dictionary& viscosityProperties);
};

View file

@ -82,24 +82,24 @@ public:
// Destructor
~Newtonian()
virtual ~Newtonian()
{}
// Member Functions
//- Return the laminar viscosity
const volScalarField& nu() const
virtual const volScalarField& nu() const
{
return nu_;
}
//- Correct the laminar viscosity (not appropriate, viscosity constant)
void correct()
virtual void correct()
{}
//- Read transportProperties dictionary
bool read(const dictionary& viscosityProperties);
virtual bool read(const dictionary& viscosityProperties);
};

View file

@ -89,26 +89,26 @@ public:
// Destructor
~freeSurface()
virtual ~freeSurface()
{}
// Member Functions
//- Return the laminar viscosity
const volScalarField& nu() const
virtual const volScalarField& nu() const
{
return nu_;
}
//- Correct the laminar viscosity
void correct()
virtual void correct()
{
nu_ = gamma_*nu1_ + (scalar(1) - gamma_)*nu2_;
}
//- Read transportProperties dictionary
bool read(const dictionary& viscosityProperties);
virtual bool read(const dictionary& viscosityProperties);
};

View file

@ -93,26 +93,26 @@ public:
// Destructor
~powerLaw()
virtual ~powerLaw()
{}
// Member Functions
//- Return the laminar viscosity
const volScalarField& nu() const
virtual const volScalarField& nu() const
{
return nu_;
}
//- Correct the laminar viscosity
void correct()
virtual void correct()
{
nu_ = calcNu();
}
//- Read transportProperties dictionary
bool read(const dictionary& viscosityProperties);
virtual bool read(const dictionary& viscosityProperties);
};

View file

@ -65,6 +65,7 @@ class geometricSurfacePatch
//- Index of patch in boundary
label index_;
public:
//- Runtime type information

View file

@ -44,11 +44,16 @@ void triSurface::writeSTLASCII(Ostream& os) const
surfacePatchList myPatches(calcPatches(faceMap));
label faceIndex = 0;
forAll(myPatches, patchI)
forAll (myPatches, patchI)
{
// Print all faces belonging to this region
const surfacePatch& patch = myPatches[patchI];
if (patch.size() == 0)
{
continue;
}
os << "solid " << patch.name() << endl;
for

View file

@ -201,9 +201,10 @@ class triSurface
const pointField&
);
//- read non-comment line
//- Read non-comment line
static string getLineNoComment(IFstream&);
protected:
// Protected Member Functions
@ -220,6 +221,7 @@ protected:
return static_cast<List<Face>&>(*this);
}
public:
// Public typedefs

View file

@ -363,7 +363,9 @@ void LRR::correct()
epsEqn().relax();
epsEqn().boundaryManipulate(epsilon_.boundaryField());
// No longer needed: matrix completes at the point of solution
// HJ, 17/Apr/2012
// epsEqn().completeAssembly();
solve(epsEqn);
bound(epsilon_, epsilon0_);
@ -382,8 +384,13 @@ void LRR::correct()
forAll(curPatch, facei)
{
label faceCelli = curPatch.faceCells()[facei];
P[faceCelli]
*= min(G[faceCelli]/(0.5*mag(tr(P[faceCelli])) + SMALL), 100.0);
P[faceCelli] *=
min
(
G[faceCelli]/(0.5*mag(tr(P[faceCelli])) + SMALL),
100.0
);
}
}
}

View file

@ -401,7 +401,9 @@ void LaunderGibsonRSTM::correct()
epsEqn().relax();
epsEqn().boundaryManipulate(epsilon_.boundaryField());
// No longer needed: matrix completes at the point of solution
// HJ, 17/Apr/2012
// epsEqn().completeAssembly();
solve(epsEqn);
bound(epsilon_, epsilon0_);

View file

@ -332,7 +332,9 @@ void RNGkEpsilon::correct()
epsEqn().relax();
epsEqn().boundaryManipulate(epsilon_.boundaryField());
// No longer needed: matrix completes at the point of solution
// HJ, 17/Apr/2012
// epsEqn().completeAssembly();
solve(epsEqn);
bound(epsilon_, epsilon0_);

View file

@ -177,6 +177,11 @@ void epsilonWallFunctionFvPatchScalarField::updateCoeffs()
// HJ, 20/Mar/2011
if (!db().foundObject<volScalarField>(GName_))
{
InfoIn("void epsilonWallFunctionFvPatchScalarField::updateCoeffs()")
<< "Cannot access " << GName_ << " field. for patch "
<< patch().name() << ". Evaluating as zeroGradient"
<< endl;
zeroGradientFvPatchScalarField::updateCoeffs();
return;

View file

@ -182,6 +182,11 @@ void omegaWallFunctionFvPatchScalarField::updateCoeffs()
// HJ, 20/Mar/2011
if (!db().foundObject<volScalarField>(GName_))
{
InfoIn("void omegaWallFunctionFvPatchScalarField::updateCoeffs()")
<< "Cannot access " << GName_ << " field. for patch "
<< patch().name() << ". Evaluating as zeroGradient"
<< endl;
zeroGradientFvPatchScalarField::evaluate();
return;

View file

@ -305,7 +305,9 @@ void kEpsilon::correct()
epsEqn().relax();
epsEqn().boundaryManipulate(epsilon_.boundaryField());
// No longer needed: matrix completes at the point of solution
// HJ, 17/Apr/2012
// epsEqn().completeAssembly();
solve(epsEqn);
bound(epsilon_, epsilon0_);

View file

@ -429,7 +429,9 @@ void kOmegaSST::correct()
omegaEqn().relax();
omegaEqn().boundaryManipulate(omega_.boundaryField());
// No longer needed: matrix completes at the point of solution
// HJ, 17/Apr/2012
// omegaEqn().completeAssembly();
solve(omegaEqn);
bound(omega_, omega0_);

View file

@ -346,7 +346,9 @@ void realizableKE::correct()
epsEqn().relax();
epsEqn().boundaryManipulate(epsilon_.boundaryField());
// No longer needed: matrix completes at the point of solution
// HJ, 17/Apr/2012
// epsEqn().completeAssembly();
solve(epsEqn);
bound(epsilon_, epsilon0_);

View file

@ -324,7 +324,9 @@ void LRR::correct()
epsEqn().relax();
epsEqn().boundaryManipulate(epsilon_.boundaryField());
// No longer needed: matrix completes at the point of solution
// HJ, 17/Apr/2012
// epsEqn().completeAssembly();
solve(epsEqn);
bound(epsilon_, epsilon0_);
@ -343,8 +345,13 @@ void LRR::correct()
forAll(curPatch, facei)
{
label faceCelli = curPatch.faceCells()[facei];
P[faceCelli]
*= min(G[faceCelli]/(0.5*mag(tr(P[faceCelli])) + SMALL), 1.0);
P[faceCelli] *=
min
(
G[faceCelli]/(0.5*mag(tr(P[faceCelli])) + SMALL),
1.0
);
}
}
}

View file

@ -365,7 +365,9 @@ void LaunderGibsonRSTM::correct()
epsEqn().relax();
epsEqn().boundaryManipulate(epsilon_.boundaryField());
// No longer needed: matrix completes at the point of solution
// HJ, 17/Apr/2012
// epsEqn().completeAssembly();
solve(epsEqn);
bound(epsilon_, epsilon0_);

View file

@ -356,7 +356,9 @@ void LienCubicKE::correct()
epsEqn().relax();
epsEqn().boundaryManipulate(epsilon_.boundaryField());
// No longer needed: matrix completes at the point of solution
// HJ, 17/Apr/2012
// epsEqn().completeAssembly();
solve(epsEqn);
bound(epsilon_, epsilon0_);

View file

@ -280,7 +280,9 @@ void RNGkEpsilon::correct()
epsEqn().relax();
epsEqn().boundaryManipulate(epsilon_.boundaryField());
// No longer needed: matrix completes at the point of solution
// HJ, 17/Apr/2012
// epsEqn().completeAssembly();
solve(epsEqn);
bound(epsilon_, epsilon0_);

View file

@ -172,6 +172,11 @@ void epsilonWallFunctionFvPatchScalarField::updateCoeffs()
// HJ, 20/Mar/2011
if (!db().foundObject<volScalarField>(GName_))
{
InfoIn("void epsilonWallFunctionFvPatchScalarField::updateCoeffs()")
<< "Cannot access " << GName_ << " field. for patch "
<< patch().name() << ". Evaluating as zeroGradient"
<< endl;
zeroGradientFvPatchScalarField::evaluate();
return;

View file

@ -177,6 +177,11 @@ void omegaWallFunctionFvPatchScalarField::updateCoeffs()
// HJ, 20/Mar/2011
if (!db().foundObject<volScalarField>(GName_))
{
InfoIn("void omegaWallFunctionFvPatchScalarField::updateCoeffs()")
<< "Cannot access " << GName_ << " field. for patch "
<< patch().name() << ". Evaluating as zeroGradient"
<< endl;
zeroGradientFvPatchScalarField::evaluate();
return;

View file

@ -242,7 +242,9 @@ void kEpsilon::correct()
epsEqn().relax();
epsEqn().boundaryManipulate(epsilon_.boundaryField());
// No longer needed: matrix completes at the point of solution
// HJ, 17/Apr/2012
// epsEqn().completeAssembly();
solve(epsEqn);
bound(epsilon_, epsilon0_);

View file

@ -251,7 +251,9 @@ void kOmega::correct()
omegaEqn().relax();
omegaEqn().boundaryManipulate(omega_.boundaryField());
// No longer needed: matrix completes at the point of solution
// HJ, 17/Apr/2012
// omegaEqn().completeAssembly();
solve(omegaEqn);
bound(omega_, omega0_);

View file

@ -375,7 +375,9 @@ void kOmegaSST::correct()
omegaEqn().relax();
omegaEqn().boundaryManipulate(omega_.boundaryField());
// No longer needed: matrix completes at the point of solution
// HJ, 17/Apr/2012
// omegaEqn().completeAssembly();
solve(omegaEqn);
bound(omega_, omega0_);

View file

@ -308,7 +308,9 @@ void realizableKE::correct()
epsEqn().relax();
epsEqn().boundaryManipulate(epsilon_.boundaryField());
// No longer needed: matrix completes at the point of solution
// HJ, 17/Apr/2012
// epsEqn().completeAssembly();
solve(epsEqn);
bound(epsilon_, epsilon0_);