Merge branch 'master' into feature/cvMesh

This commit is contained in:
laurence
2013-05-13 10:54:26 +01:00
39 changed files with 960 additions and 361 deletions

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -22,7 +22,7 @@ License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Application
uncoupledKinematicParcelFoam
icoUncoupledKinematicParcelFoam
Description
Transient solver for the passive transport of a single kinematic

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -26,7 +26,7 @@ Application
Description
Transient solver for the passive transport of a single kinematic
particle could.
particle cloud.
Uses a pre- calculated velocity field to evolve the cloud.

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -108,9 +108,8 @@ public:
// and a coefficient to multiply alphal for the vaporisation rate
virtual Pair<tmp<volScalarField> > mDotAlphal() const;
//- Return the mass condensation and vaporisation rates as an
// explicit term for the condensation rate and a coefficient to
// multiply (p - pSat) for the vaporisation rate
//- Return the mass condensation and vaporisation rates as coefficients
// to multiply (p - pSat)
virtual Pair<tmp<volScalarField> > mDotP() const;
//- Correct the Kunz phaseChange model

View File

@ -1,8 +1,8 @@
/*---------------------------------------------------------------------------*\
========Merkle= |
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -102,9 +102,8 @@ public:
// and a coefficient to multiply alphal for the vaporisation rate
virtual Pair<tmp<volScalarField> > mDotAlphal() const;
//- Return the mass condensation and vaporisation rates as an
// explicit term for the condensation rate and a coefficient to
// multiply (p - pSat) for the vaporisation rate
//- Return the mass condensation and vaporisation rates as coefficients
// to multiply (p - pSat)
virtual Pair<tmp<volScalarField> > mDotP() const;
//- Correct the Merkle phaseChange model

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -115,9 +115,8 @@ public:
// and a coefficient to multiply alphal for the vaporisation rate
virtual Pair<tmp<volScalarField> > mDotAlphal() const;
//- Return the mass condensation and vaporisation rates as an
// explicit term for the condensation rate and a coefficient to
// multiply (p - pSat) for the vaporisation rate
//- Return the mass condensation and vaporisation rates as coefficients
// to multiply (p - pSat)
virtual Pair<tmp<volScalarField> > mDotP() const;
//- Correct the SchnerrSauer phaseChange model

View File

@ -139,9 +139,8 @@ public:
// and a coefficient to multiply alphal for the vaporisation rate
virtual Pair<tmp<volScalarField> > mDotAlphal() const = 0;
//- Return the mass condensation and vaporisation rates as an
// explicit term for the condensation rate and a coefficient to
// multiply (p - pSat) for the vaporisation rate
//- Return the mass condensation and vaporisation rates as coefficients
// to multiply (p - pSat)
virtual Pair<tmp<volScalarField> > mDotP() const = 0;
//- Return the volumetric condensation and vaporisation rates as a
@ -149,9 +148,8 @@ public:
// and a coefficient to multiply alphal for the vaporisation rate
Pair<tmp<volScalarField> > vDotAlphal() const;
//- Return the volumetric condensation and vaporisation rates as an
// explicit term for the condensation rate and a coefficient to
// multiply (p - pSat) for the vaporisation rate
//- Return the volumetric condensation and vaporisation rates as
// coefficients to multiply (p - pSat)
Pair<tmp<volScalarField> > vDotP() const;
//- Correct the phaseChange model

View File

@ -111,11 +111,6 @@ Foam::multiphaseMixture::multiphaseMixture
{
calcAlphas();
alphas_.write();
forAllIter(PtrDictionary<phase>, phases_, iter)
{
alphaTable_.add(iter());
}
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -47,7 +47,6 @@ SourceFiles
#include "PtrDictionary.H"
#include "volFields.H"
#include "surfaceFields.H"
#include "multivariateSurfaceInterpolationScheme.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -156,9 +155,6 @@ private:
//- Conversion factor for degrees into radians
static const scalar convertToRad;
//- Phase-fraction field table for multivariate discretisation
multivariateSurfaceInterpolationScheme<scalar>::fieldTable alphaTable_;
// Private member functions

View File

@ -47,6 +47,10 @@ x 5;
varName x;
//Indirection for keys
key inlet_9;
boundaryField
{
Default_Boundary_Region
@ -63,6 +67,7 @@ boundaryField
inlet_6 { $.inactive } // Test scoping
inlet_7 { ${inactive}} // Test variable expansion
inlet_8 { $inactive }
${key} { $inactive }
#include "testDictInc"

View File

@ -402,9 +402,10 @@ addLayersControls
nLayerIter 50;
// Max number of iterations after which relaxed meshQuality controls
// get used. Up to nRelaxIter it uses the settings in
// get used. Up to nRelaxedIter it uses the settings in
// meshQualityControls,
// after nRelaxIter it uses the values in meshQualityControls::relaxed.
// after nRelaxedIter it uses the values in
// meshQualityControls::relaxed.
nRelaxedIter 20;
// Additional reporting: if there are just a few faces where there

View File

@ -488,7 +488,10 @@ void Foam::GAMGAgglomeration::procAgglomerateRestrictAddressing
comm,
procIDs,
restrictAddressing_[levelIndex],
procRestrictAddressing
procRestrictAddressing,
UPstream::msgType(),
Pstream::nonBlocking //Pstream::scheduled
);

View File

@ -126,8 +126,15 @@ void Foam::GAMGAgglomeration::restrictField
const List<int>& procIDs = agglomProcIDs(coarseLevelIndex);
const labelList& offsets = cellOffsets(coarseLevelIndex);
const globalIndex cellOffsetter(offsets);
cellOffsetter.gather(fineComm, procIDs, cf);
globalIndex::gather
(
offsets,
fineComm,
procIDs,
cf,
UPstream::msgType(),
Pstream::nonBlocking //Pstream::scheduled
);
}
}
@ -192,18 +199,18 @@ void Foam::GAMGAgglomeration::prolongField
const List<int>& procIDs = agglomProcIDs(coarseLevelIndex);
const labelList& offsets = cellOffsets(coarseLevelIndex);
const globalIndex cellOffsetter(offsets);
label localSize = nCells_[levelIndex];
Field<Type> allCf(localSize);
cellOffsetter.scatter
globalIndex::scatter
(
offsets,
coarseComm,
procIDs,
cf,
allCf
allCf,
UPstream::msgType(),
Pstream::nonBlocking //Pstream::scheduled
);
forAll(fineToCoarse, i)

View File

@ -25,6 +25,7 @@ License
#include "GAMGProcAgglomeration.H"
#include "GAMGAgglomeration.H"
#include "lduMesh.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -114,6 +115,152 @@ void Foam::GAMGProcAgglomeration::printStats
}
Foam::labelListList Foam::GAMGProcAgglomeration::globalCellCells
(
const lduMesh& mesh
)
{
const lduAddressing& addr = mesh.lduAddr();
lduInterfacePtrsList interfaces = mesh.interfaces();
const label myProcID = Pstream::myProcNo(mesh.comm());
globalIndex globalNumbering
(
addr.size(),
Pstream::msgType(),
mesh.comm(),
Pstream::parRun()
);
labelList globalIndices(addr.size());
forAll(globalIndices, cellI)
{
globalIndices[cellI] = globalNumbering.toGlobal(myProcID, cellI);
}
// Get the interface cells
PtrList<labelList> nbrGlobalCells(interfaces.size());
{
// Initialise transfer of restrict addressing on the interface
forAll(interfaces, inti)
{
if (interfaces.set(inti))
{
interfaces[inti].initInternalFieldTransfer
(
Pstream::nonBlocking,
globalIndices
);
}
}
if (Pstream::parRun())
{
Pstream::waitRequests();
}
forAll(interfaces, inti)
{
if (interfaces.set(inti))
{
nbrGlobalCells.set
(
inti,
new labelList
(
interfaces[inti].internalFieldTransfer
(
Pstream::nonBlocking,
globalIndices
)
)
);
}
}
}
// Scan the neighbour list to find out how many times the cell
// appears as a neighbour of the face. Done this way to avoid guessing
// and resizing list
labelList nNbrs(addr.size(), 1);
const labelUList& nbr = addr.upperAddr();
const labelUList& own = addr.lowerAddr();
{
forAll(nbr, faceI)
{
nNbrs[nbr[faceI]]++;
nNbrs[own[faceI]]++;
}
forAll(interfaces, inti)
{
if (interfaces.set(inti))
{
const labelUList& faceCells = interfaces[inti].faceCells();
forAll(faceCells, i)
{
nNbrs[faceCells[i]]++;
}
}
}
}
// Create cell-cells addressing
labelListList cellCells(addr.size());
forAll(cellCells, cellI)
{
cellCells[cellI].setSize(nNbrs[cellI], -1);
}
// Reset the list of number of neighbours to zero
nNbrs = 0;
// Scatter the neighbour faces
forAll(nbr, faceI)
{
label c0 = own[faceI];
label c1 = nbr[faceI];
cellCells[c0][nNbrs[c0]++] = globalIndices[c1];
cellCells[c1][nNbrs[c1]++] = globalIndices[c0];
}
forAll(interfaces, inti)
{
if (interfaces.set(inti))
{
const labelUList& faceCells = interfaces[inti].faceCells();
forAll(faceCells, i)
{
label c0 = faceCells[i];
cellCells[c0][nNbrs[c0]++] = nbrGlobalCells[inti][i];
}
}
}
forAll(cellCells, cellI)
{
Foam::stableSort(cellCells[cellI]);
}
// Replace the initial element (always -1) with the local cell
forAll(cellCells, cellI)
{
cellCells[cellI][0] = globalIndices[cellI];
}
return cellCells;
}
bool Foam::GAMGProcAgglomeration::agglomerate
(
const label fineLevelIndex,
@ -243,7 +390,4 @@ Foam::GAMGProcAgglomeration::~GAMGProcAgglomeration()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// ************************************************************************* //

View File

@ -44,6 +44,7 @@ namespace Foam
{
class GAMGAgglomeration;
class lduMesh;
/*---------------------------------------------------------------------------*\
Class GAMGProcAgglomeration Declaration
@ -74,6 +75,9 @@ protected:
const label procAgglomComm
);
//- Debug: calculate global cell-cells
static labelListList globalCellCells(const lduMesh&);
private:
// Private data
@ -139,7 +143,7 @@ public:
// Member Functions
//- Modify agglomeration. Return true if modified
//- Modify agglomeration. Return true if modified
virtual bool agglomerate() = 0;
};

View File

@ -109,7 +109,7 @@ bool Foam::manualGAMGProcAgglomeration::agglomerate()
// My processor id
const label myProcID = Pstream::myProcNo(levelMesh.comm());
const List<clusterAndMaster>& clusters =
const List<labelList>& clusters =
procAgglomMaps_[i].second();
// Coarse to fine master processor
@ -125,8 +125,8 @@ bool Foam::manualGAMGProcAgglomeration::agglomerate()
forAll(clusters, coarseI)
{
const labelList& cluster = clusters[coarseI].first();
coarseToMaster[coarseI] = clusters[coarseI].second();
const labelList& cluster = clusters[coarseI];
coarseToMaster[coarseI] = cluster[0];
forAll(cluster, i)
{

View File

@ -30,19 +30,23 @@ Description
In the GAMG control dictionary:
processorAgglomerator manual;
// List of level+procagglomeration where
// procagglomeration is a set of labelLists. Each labelList is
// a cluster of processor which gets combined onto the first element
// in the list.
processorAgglomeration
(
(
3 //at level 3
(
((0 1) 0) //coarse 0 from 0,1 (and moved onto 0)
((2 3) 3) //coarse 1 from 2,3 (and moved onto 3)
(0 1) //coarse 0 from 0,1 (and moved onto 0)
(3 2) //coarse 1 from 2,3 (and moved onto 3)
)
)
(
6 //at level6
(
((0 1) 0) //coarse 0 from 0,1 (and moved onto 0)
(0 1) //coarse 0 from 0,1 (and moved onto 0)
)
)
);
@ -76,10 +80,8 @@ class manualGAMGProcAgglomeration
{
// Private data
typedef Tuple2<labelList, label> clusterAndMaster;
//- Per level the agglomeration map
const List<Tuple2<label, List<clusterAndMaster> > > procAgglomMaps_;
const List<Tuple2<label, List<labelList> > > procAgglomMaps_;
//- Any allocated communicators
DynamicList<label> comms_;

View File

@ -114,13 +114,6 @@ Foam::solverPerformance Foam::GAMGSolver::solve
cmpt
);
//Pout<< "finestCorrection:" << finestCorrection << endl;
//Pout<< "finestResidual:" << finestResidual << endl;
//Pout<< "psi:" << psi << endl;
//Pout<< "Apsi:" << Apsi << endl;
// Calculate finest level residual field
matrix_.Amul(Apsi, psi, interfaceBouCoeffs_, interfaces_, cmpt);
finestResidual = source;
@ -205,7 +198,6 @@ void Foam::GAMGSolver::Vcycle
scratch1,
coarseCorrFields[leveli].size()
);
//scalarField ACf(coarseCorrFields[leveli].size(), VGREAT);
// Scale coarse-grid correction field
// but not on the coarsest level because it evaluates to 1
@ -218,8 +210,6 @@ void Foam::GAMGSolver::Vcycle
(
ACf.operator const scalarField&()
),
//ACf,
matrixLevels_[leveli],
interfaceLevelsBouCoeffs_[leveli],
interfaceLevels_[leveli],
@ -235,7 +225,6 @@ void Foam::GAMGSolver::Vcycle
(
ACf.operator const scalarField&()
),
//ACf,
coarseCorrFields[leveli],
interfaceLevelsBouCoeffs_[leveli],
interfaceLevels_[leveli],
@ -294,11 +283,6 @@ void Foam::GAMGSolver::Vcycle
scratch2,
coarseCorrFields[leveli].size()
);
//scalarField preSmoothedCoarseCorrField
//(
// coarseCorrFields[leveli].size(),
// VGREAT
//);
// Only store the preSmoothedCoarseCorrField if pre-smoothing is
// used
@ -328,12 +312,6 @@ void Foam::GAMGSolver::Vcycle
);
scalarField& ACfRef =
const_cast<scalarField&>(ACf.operator const scalarField&());
//scalarField ACfRef
//(
// coarseCorrFields[leveli].size(),
// VGREAT
//);
if (interpolateCorrection_)
{

View File

@ -29,7 +29,6 @@ Description
SourceFiles
processorGAMGInterface.C
processorGAMGInterfaceTemplates.C
\*---------------------------------------------------------------------------*/

View File

@ -360,6 +360,7 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
+ procMesh.lduAddr().size();
}
// Faces initially get added in order, sorted later
labelList internalFaceOffsets(nMeshes+1);
internalFaceOffsets[0] = 0;
@ -818,9 +819,37 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
{
const labelPairList& elems = iter();
// Sort processors in increasing order
labelList order(identity(elems.size()));
Foam::sort(order, procLess(elems));
// Sort processors in increasing order so both sides walk through in
// same order.
labelPairList procPairs(elems.size());
forAll(elems, i)
{
const labelPair& elem = elems[i];
label procMeshI = elem[0];
label interfaceI = elem[1];
const lduInterfacePtrsList interfaces = mesh
(
myMesh,
otherMeshes,
procMeshI
).interfaces();
const processorLduInterface& pldui =
refCast<const processorLduInterface>
(
interfaces[interfaceI]
);
label myProcNo = pldui.myProcNo();
label nbrProcNo = pldui.neighbProcNo();
procPairs[i] = labelPair
(
min(myProcNo, nbrProcNo),
max(myProcNo, nbrProcNo)
);
}
labelList order;
sortedOrder(procPairs, order);
// Count
label n = 0;
@ -837,11 +866,6 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
procMeshI
).interfaces();
//Pout<< " adding interface:" << " procMeshI:" << procMeshI
// << " interface:" << interfaceI
// << " size:" << interfaces[interfaceI].faceCells().size()
// << endl;
n += interfaces[interfaceI].faceCells().size();
}
@ -868,6 +892,7 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
const labelUList& l = interfaces[interfaceI].faceCells();
bfMap.setSize(l.size());
forAll(l, faceI)
{
allFaceCells[n] = cellOffsets[procMeshI]+l[faceI];
@ -947,7 +972,10 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
patchSchedule_ = nonBlockingSchedule<processorGAMGInterface>(interfaces_);
checkUpperTriangular(cellOffsets.last(), lowerAddr_, upperAddr_);
if (debug)
{
checkUpperTriangular(cellOffsets.last(), lowerAddr_, upperAddr_);
}
}

View File

@ -150,6 +150,20 @@ public:
// Other
//- Collect data in processor order on master (== procIDs[0]).
// Needs offsets only on master.
template<class Type>
static void gather
(
const labelUList& offsets,
const label comm,
const labelList& procIDs,
const UList<Type>& fld,
List<Type>& allFld,
const int tag = UPstream::msgType(),
const Pstream::commsTypes commsType=Pstream::nonBlocking
);
//- Collect data in processor order on master (== procIDs[0]).
// Needs offsets only on master.
template<class Type>
@ -159,8 +173,25 @@ public:
const labelList& procIDs,
const UList<Type>& fld,
List<Type>& allFld,
const int tag = UPstream::msgType()
) const;
const int tag = UPstream::msgType(),
const Pstream::commsTypes commsType=Pstream::nonBlocking
) const
{
gather(offsets_, comm, procIDs, fld, allFld, tag, commsType);
}
//- Inplace collect data in processor order on master
// (== procIDs[0]). Needs offsets only on master.
template<class Type>
static void gather
(
const labelUList& offsets,
const label comm,
const labelList& procIDs,
List<Type>& fld,
const int tag = UPstream::msgType(),
const Pstream::commsTypes commsType=Pstream::nonBlocking
);
//- Inplace collect data in processor order on master
// (== procIDs[0]). Needs offsets only on master.
@ -170,8 +201,25 @@ public:
const label comm,
const labelList& procIDs,
List<Type>& fld,
const int tag = UPstream::msgType()
) const;
const int tag = UPstream::msgType(),
const Pstream::commsTypes commsType=Pstream::nonBlocking
) const
{
gather(offsets_, comm, procIDs, fld, tag, commsType);
}
//- Distribute data in processor order. Requires fld to be sized!
template<class Type>
static void scatter
(
const labelUList& offsets,
const label comm,
const labelList& procIDs,
const UList<Type>& allFld,
UList<Type>& fld,
const int tag = UPstream::msgType(),
const Pstream::commsTypes commsType=Pstream::nonBlocking
);
//- Distribute data in processor order. Requires fld to be sized!
template<class Type>
@ -181,8 +229,12 @@ public:
const labelList& procIDs,
const UList<Type>& allFld,
UList<Type>& fld,
const int tag = UPstream::msgType()
) const;
const int tag = UPstream::msgType(),
const Pstream::commsTypes commsType=Pstream::nonBlocking
) const
{
scatter(offsets_, comm, procIDs, allFld, fld, tag, commsType);
}
// IOstream Operators

View File

@ -30,34 +30,75 @@ License
template<class Type>
void Foam::globalIndex::gather
(
const labelUList& off,
const label comm,
const labelList& procIDs,
const UList<Type>& fld,
List<Type>& allFld,
const int tag
) const
const int tag,
const Pstream::commsTypes commsType
)
{
if (Pstream::myProcNo(comm) == procIDs[0])
{
allFld.setSize(size());
allFld.setSize(off.last());
// Assign my local data
SubList<Type>(allFld, fld.size(), 0).assign(fld);
for (label i = 1; i < procIDs.size(); i++)
if (commsType == Pstream::scheduled || commsType == Pstream::blocking)
{
SubList<Type> procSlot
(
allFld,
offsets_[i+1]-offsets_[i],
offsets_[i]
);
if (contiguous<Type>())
for (label i = 1; i < procIDs.size(); i++)
{
SubList<Type> procSlot(allFld, off[i+1]-off[i], off[i]);
if (contiguous<Type>())
{
IPstream::read
(
commsType,
procIDs[i],
reinterpret_cast<char*>(procSlot.begin()),
procSlot.byteSize(),
tag,
comm
);
}
else
{
IPstream fromSlave
(
commsType,
procIDs[i],
0,
tag,
comm
);
fromSlave >> procSlot;
}
}
}
else
{
// nonBlocking
if (!contiguous<Type>())
{
FatalErrorIn("globalIndex::gather(..)")
<< "nonBlocking not supported for non-contiguous data"
<< exit(FatalError);
}
label startOfRequests = Pstream::nRequests();
// Set up reads
for (label i = 1; i < procIDs.size(); i++)
{
SubList<Type> procSlot(allFld, off[i+1]-off[i], off[i]);
IPstream::read
(
Pstream::scheduled,
commsType,
procIDs[i],
reinterpret_cast<char*>(procSlot.begin()),
procSlot.byteSize(),
@ -65,45 +106,66 @@ void Foam::globalIndex::gather
comm
);
}
else
{
IPstream fromSlave
(
Pstream::scheduled,
procIDs[i],
0,
tag,
comm
);
fromSlave >> procSlot;
}
// Wait for all to finish
Pstream::waitRequests(startOfRequests);
}
}
else
{
if (contiguous<Type>())
if (commsType == Pstream::scheduled || commsType == Pstream::blocking)
{
if (contiguous<Type>())
{
OPstream::write
(
commsType,
procIDs[0],
reinterpret_cast<const char*>(fld.begin()),
fld.byteSize(),
tag,
comm
);
}
else
{
OPstream toMaster
(
commsType,
procIDs[0],
0,
tag,
comm
);
toMaster << fld;
}
}
else
{
// nonBlocking
if (!contiguous<Type>())
{
FatalErrorIn("globalIndex::gather(..)")
<< "nonBlocking not supported for non-contiguous data"
<< exit(FatalError);
}
label startOfRequests = Pstream::nRequests();
// Set up write
OPstream::write
(
Pstream::scheduled,
commsType,
procIDs[0],
reinterpret_cast<const char*>(fld.begin()),
fld.byteSize(),
tag,
comm
);
}
else
{
OPstream toMaster
(
Pstream::scheduled,
procIDs[0],
0,
tag,
comm
);
toMaster << fld;
// Wait for all to finish
Pstream::waitRequests(startOfRequests);
}
}
}
@ -112,121 +174,104 @@ void Foam::globalIndex::gather
template<class Type>
void Foam::globalIndex::gather
(
const labelUList& off,
const label comm,
const labelList& procIDs,
List<Type>& fld,
const int tag
) const
const int tag,
const Pstream::commsTypes commsType
)
{
List<Type> allFld;
gather(off, comm, procIDs, fld, allFld, tag, commsType);
if (Pstream::myProcNo(comm) == procIDs[0])
{
List<Type> allFld(size());
// Assign my local data
SubList<Type>(allFld, fld.size(), 0).assign(fld);
for (label i = 1; i < procIDs.size(); i++)
{
SubList<Type> procSlot
(
allFld,
offsets_[i+1]-offsets_[i],
offsets_[i]
);
if (contiguous<Type>())
{
IPstream::read
(
Pstream::scheduled,
procIDs[i],
reinterpret_cast<char*>(procSlot.begin()),
procSlot.byteSize(),
tag,
comm
);
}
else
{
IPstream fromSlave
(
Pstream::scheduled,
procIDs[i],
0,
tag,
comm
);
fromSlave >> procSlot;
}
}
fld.transfer(allFld);
}
else
{
if (contiguous<Type>())
{
OPstream::write
(
Pstream::scheduled,
procIDs[0],
reinterpret_cast<const char*>(fld.begin()),
fld.byteSize(),
tag,
comm
);
}
else
{
OPstream toMaster
(
Pstream::scheduled,
procIDs[0],
0,
tag,
comm
);
toMaster << fld;
}
}
}
template<class Type>
void Foam::globalIndex::scatter
(
const labelUList& off,
const label comm,
const labelList& procIDs,
const UList<Type>& allFld,
UList<Type>& fld,
const int tag
) const
const int tag,
const Pstream::commsTypes commsType
)
{
if (Pstream::myProcNo(comm) == procIDs[0])
{
fld.assign
(
SubList<Type>
(
allFld,
offsets_[1]-offsets_[0]
)
);
fld.assign(SubList<Type>(allFld, off[1]-off[0]));
for (label i = 1; i < procIDs.size(); i++)
if (commsType == Pstream::scheduled || commsType == Pstream::blocking)
{
const SubList<Type> procSlot
(
allFld,
offsets_[i+1]-offsets_[i],
offsets_[i]
);
if (contiguous<Type>())
for (label i = 1; i < procIDs.size(); i++)
{
const SubList<Type> procSlot
(
allFld,
off[i+1]-off[i],
off[i]
);
if (contiguous<Type>())
{
OPstream::write
(
commsType,
procIDs[i],
reinterpret_cast<const char*>(procSlot.begin()),
procSlot.byteSize(),
tag,
comm
);
}
else
{
OPstream toSlave
(
commsType,
procIDs[i],
0,
tag,
comm
);
toSlave << procSlot;
}
}
}
else
{
// nonBlocking
if (!contiguous<Type>())
{
FatalErrorIn("globalIndex::scatter(..)")
<< "nonBlocking not supported for non-contiguous data"
<< exit(FatalError);
}
label startOfRequests = Pstream::nRequests();
// Set up writes
for (label i = 1; i < procIDs.size(); i++)
{
const SubList<Type> procSlot
(
allFld,
off[i+1]-off[i],
off[i]
);
OPstream::write
(
Pstream::scheduled,
commsType,
procIDs[i],
reinterpret_cast<const char*>(procSlot.begin()),
procSlot.byteSize(),
@ -234,45 +279,66 @@ void Foam::globalIndex::scatter
comm
);
}
else
{
OPstream toSlave
(
Pstream::scheduled,
procIDs[i],
0,
tag,
comm
);
toSlave << procSlot;
}
// Wait for all to finish
Pstream::waitRequests(startOfRequests);
}
}
else
{
if (contiguous<Type>())
if (commsType == Pstream::scheduled || commsType == Pstream::blocking)
{
if (contiguous<Type>())
{
IPstream::read
(
commsType,
procIDs[0],
reinterpret_cast<char*>(fld.begin()),
fld.byteSize(),
tag,
comm
);
}
else
{
IPstream fromMaster
(
commsType,
procIDs[0],
0,
tag,
comm
);
fromMaster >> fld;
}
}
else
{
// nonBlocking
if (!contiguous<Type>())
{
FatalErrorIn("globalIndex::scatter(..)")
<< "nonBlocking not supported for non-contiguous data"
<< exit(FatalError);
}
label startOfRequests = Pstream::nRequests();
// Set up read
IPstream::read
(
Pstream::scheduled,
commsType,
procIDs[0],
reinterpret_cast<char*>(fld.begin()),
fld.byteSize(),
tag,
comm
);
}
else
{
IPstream fromMaster
(
Pstream::scheduled,
procIDs[0],
0,
tag,
comm
);
fromMaster >> fld;
// Wait for all to finish
Pstream::waitRequests(startOfRequests);
}
}
}

View File

@ -52,8 +52,8 @@ Foam::porosityModels::DarcyForchheimer::DarcyForchheimer
)
:
porosityModel(name, modelType, mesh, dict, cellZoneName),
D_(cellZoneIds_.size()),
F_(cellZoneIds_.size()),
D_(cellZoneIDs_.size()),
F_(cellZoneIDs_.size()),
rhoName_(coeffs_.lookupOrDefault<word>("rho", "rho")),
muName_(coeffs_.lookupOrDefault<word>("mu", "thermo:mu")),
nuName_(coeffs_.lookupOrDefault<word>("nu", "nu"))
@ -67,7 +67,7 @@ Foam::porosityModels::DarcyForchheimer::DarcyForchheimer
if (coordSys_.R().uniform())
{
forAll (cellZoneIds_, zoneI)
forAll (cellZoneIDs_, zoneI)
{
D_[zoneI].setSize(1, tensor::zero);
F_[zoneI].setSize(1, tensor::zero);
@ -89,9 +89,9 @@ Foam::porosityModels::DarcyForchheimer::DarcyForchheimer
}
else
{
forAll(cellZoneIds_, zoneI)
forAll(cellZoneIDs_, zoneI)
{
const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]];
const labelList& cells = mesh_.cellZones()[cellZoneIDs_[zoneI]];
D_[zoneI].setSize(cells.size(), tensor::zero);
F_[zoneI].setSize(cells.size(), tensor::zero);
@ -122,6 +122,24 @@ Foam::porosityModels::DarcyForchheimer::~DarcyForchheimer()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::porosityModels::DarcyForchheimer::calcForce
(
const volVectorField& U,
const volScalarField& rho,
const volScalarField& mu,
vectorField& force
) const
{
scalarField Udiag(U.size());
vectorField Usource(U.size());
const scalarField& V = mesh_.V();
apply(Udiag, Usource, V, rho, mu, U);
force = Udiag*U - Usource;
}
void Foam::porosityModels::DarcyForchheimer::correct
(
fvVectorMatrix& UEqn
@ -194,10 +212,12 @@ void Foam::porosityModels::DarcyForchheimer::correct
}
void Foam::porosityModels::DarcyForchheimer::writeData(Ostream& os) const
bool Foam::porosityModels::DarcyForchheimer::writeData(Ostream& os) const
{
os << indent << name_ << endl;
dict_.write(os);
return true;
}

View File

@ -143,6 +143,15 @@ public:
// Member Functions
//- Calculate the porosity force
virtual void calcForce
(
const volVectorField& U,
const volScalarField& rho,
const volScalarField& mu,
vectorField& force
) const;
//- Add resistance
virtual void correct(fvVectorMatrix& UEqn) const;
@ -165,7 +174,7 @@ public:
// I-O
//- Write
void writeData(Ostream& os) const;
bool writeData(Ostream& os) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View File

@ -36,12 +36,12 @@ void Foam::porosityModels::DarcyForchheimer::apply
const vectorField& U
) const
{
forAll(cellZoneIds_, zoneI)
forAll(cellZoneIDs_, zoneI)
{
const tensorField& dZones = D_[zoneI];
const tensorField& fZones = F_[zoneI];
const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]];
const labelList& cells = mesh_.cellZones()[cellZoneIDs_[zoneI]];
forAll(cells, i)
{
@ -68,12 +68,12 @@ void Foam::porosityModels::DarcyForchheimer::apply
const vectorField& U
) const
{
forAll(cellZoneIds_, zoneI)
forAll(cellZoneIDs_, zoneI)
{
const tensorField& dZones = D_[zoneI];
const tensorField& fZones = F_[zoneI];
const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]];
const labelList& cells = mesh_.cellZones()[cellZoneIDs_[zoneI]];
forAll(cells, i)
{

View File

@ -50,12 +50,12 @@ void Foam::porosityModels::fixedCoeff::apply
const scalar rho
) const
{
forAll(cellZoneIds_, zoneI)
forAll(cellZoneIDs_, zoneI)
{
const tensorField& alphaZones = alpha_[zoneI];
const tensorField& betaZones = beta_[zoneI];
const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]];
const labelList& cells = mesh_.cellZones()[cellZoneIDs_[zoneI]];
forAll(cells, i)
{
@ -79,12 +79,12 @@ void Foam::porosityModels::fixedCoeff::apply
) const
{
forAll(cellZoneIds_, zoneI)
forAll(cellZoneIDs_, zoneI)
{
const tensorField& alphaZones = alpha_[zoneI];
const tensorField& betaZones = beta_[zoneI];
const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]];
const labelList& cells = mesh_.cellZones()[cellZoneIDs_[zoneI]];
forAll(cells, i)
{
@ -111,8 +111,8 @@ Foam::porosityModels::fixedCoeff::fixedCoeff
)
:
porosityModel(name, modelType, mesh, dict, cellZoneName),
alpha_(cellZoneIds_.size()),
beta_(cellZoneIds_.size())
alpha_(cellZoneIDs_.size()),
beta_(cellZoneIDs_.size())
{
dimensionedVector alpha(coeffs_.lookup("alpha"));
dimensionedVector beta(coeffs_.lookup("beta"));
@ -122,7 +122,7 @@ Foam::porosityModels::fixedCoeff::fixedCoeff
if (coordSys_.R().uniform())
{
forAll (cellZoneIds_, zoneI)
forAll (cellZoneIDs_, zoneI)
{
alpha_[zoneI].setSize(1, tensor::zero);
beta_[zoneI].setSize(1, tensor::zero);
@ -140,9 +140,9 @@ Foam::porosityModels::fixedCoeff::fixedCoeff
}
else
{
forAll(cellZoneIds_, zoneI)
forAll(cellZoneIDs_, zoneI)
{
const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]];
const labelList& cells = mesh_.cellZones()[cellZoneIDs_[zoneI]];
alpha_[zoneI].setSize(cells.size(), tensor::zero);
beta_[zoneI].setSize(cells.size(), tensor::zero);
@ -175,6 +175,25 @@ Foam::porosityModels::fixedCoeff::~fixedCoeff()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::porosityModels::fixedCoeff::calcForce
(
const volVectorField& U,
const volScalarField& rho,
const volScalarField& mu,
vectorField& force
) const
{
scalarField Udiag(U.size());
vectorField Usource(U.size());
const scalarField& V = mesh_.V();
scalar rhoRef = readScalar(coeffs_.lookup("rhoRef"));
apply(Udiag, Usource, V, U, rhoRef);
force = Udiag*U - Usource;
}
void Foam::porosityModels::fixedCoeff::correct
(
fvVectorMatrix& UEqn
@ -235,10 +254,12 @@ void Foam::porosityModels::fixedCoeff::correct
}
void Foam::porosityModels::fixedCoeff::writeData(Ostream& os) const
bool Foam::porosityModels::fixedCoeff::writeData(Ostream& os) const
{
os << indent << name_ << endl;
dict_.write(os);
return true;
}

View File

@ -119,6 +119,15 @@ public:
// Member Functions
//- Calculate the porosity force
virtual void calcForce
(
const volVectorField& U,
const volScalarField& rho,
const volScalarField& mu,
vectorField& force
) const;
//- Add resistance
virtual void correct(fvVectorMatrix& UEqn) const;
@ -141,7 +150,7 @@ public:
// I-O
//- Write
void writeData(Ostream& os) const;
bool writeData(Ostream& os) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View File

@ -88,13 +88,14 @@ Foam::porosityModel::porosityModel
const word& cellZoneName
)
:
MeshObject<fvMesh, Foam::UpdateableMeshObject, porosityModel>(mesh),
name_(name),
mesh_(mesh),
dict_(dict),
coeffs_(dict.subDict(modelType + "Coeffs")),
active_(true),
zoneName_(cellZoneName),
cellZoneIds_(),
cellZoneIDs_(),
coordSys_(coordinateSystem::New(mesh, coeffs_))
{
if (zoneName_ == word::null)
@ -103,11 +104,11 @@ Foam::porosityModel::porosityModel
dict_.lookup("cellZone") >> zoneName_;
}
cellZoneIds_ = mesh_.cellZones().findIndices(zoneName_);
cellZoneIDs_ = mesh_.cellZones().findIndices(zoneName_);
Info<< " creating porous zone: " << zoneName_ << endl;
bool foundZone = !cellZoneIds_.empty();
bool foundZone = !cellZoneIDs_.empty();
reduce(foundZone, orOp<bool>());
if (!foundZone && Pstream::master())
@ -136,12 +137,30 @@ Foam::porosityModel::~porosityModel()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::tmp<Foam::vectorField> Foam::porosityModel::porosityModel::force
(
const volVectorField& U,
const volScalarField& rho,
const volScalarField& mu
) const
{
tmp<vectorField> tforce(new vectorField(U.size(), vector::zero));
if (!cellZoneIDs_.empty())
{
this->calcForce(U, rho, mu, tforce());
}
return tforce;
}
void Foam::porosityModel::porosityModel::addResistance
(
fvVectorMatrix& UEqn
) const
{
if (cellZoneIds_.empty())
if (cellZoneIDs_.empty())
{
return;
}
@ -157,7 +176,7 @@ void Foam::porosityModel::porosityModel::addResistance
const volScalarField& mu
) const
{
if (cellZoneIds_.empty())
if (cellZoneIDs_.empty())
{
return;
}
@ -173,7 +192,7 @@ void Foam::porosityModel::porosityModel::addResistance
bool correctAUprocBC
) const
{
if (cellZoneIds_.empty())
if (cellZoneIDs_.empty())
{
return;
}
@ -190,12 +209,30 @@ void Foam::porosityModel::porosityModel::addResistance
}
bool Foam::porosityModel::porosityModel::movePoints()
{
// no updates necessary; all member data independent of mesh
return true;
}
void Foam::porosityModel::porosityModel::updateMesh(const mapPolyMesh& mpm)
{
// no updates necessary; all member data independent of mesh
}
bool Foam::porosityModel::writeData(Ostream& os) const
{
return true;
}
bool Foam::porosityModel::read(const dictionary& dict)
{
active_ = readBool(dict.lookup("active"));
coeffs_ = dict.subDict(type() + "Coeffs");
dict.lookup("cellZone") >> zoneName_;
cellZoneIds_ = mesh_.cellZones().findIndices(zoneName_);
cellZoneIDs_ = mesh_.cellZones().findIndices(zoneName_);
return true;
}

View File

@ -36,6 +36,7 @@ SourceFiles
#ifndef porosityModel_H
#define porosityModel_H
#include "MeshObject.H"
#include "fvMesh.H"
#include "dictionary.H"
#include "fvMatricesFwd.H"
@ -54,6 +55,8 @@ namespace Foam
\*---------------------------------------------------------------------------*/
class porosityModel
:
public MeshObject<fvMesh, UpdateableMeshObject, porosityModel>
{
private:
@ -88,8 +91,8 @@ protected:
//- Name(s) of cell-zone
keyType zoneName_;
//- Cell zone Ids
labelList cellZoneIds_;
//- Cell zone IDs
labelList cellZoneIDs_;
//- Local co-ordinate system
coordinateSystem coordSys_;
@ -100,6 +103,15 @@ protected:
//- Adjust negative resistance values to be multiplier of max value
void adjustNegativeResistance(dimensionedVector& resist);
//- Calculate the porosity force
virtual void calcForce
(
const volVectorField& U,
const volScalarField& rho,
const volScalarField& mu,
vectorField& force
) const = 0;
virtual void correct(fvVectorMatrix& UEqn) const = 0;
virtual void correct
@ -207,6 +219,17 @@ public:
//- Return const access to the porosity active flag
inline bool active() const;
//- Return const access to the cell zone IDs
inline const labelList& cellZoneIDs() const;
//- Return the force over the cell zone(s)
virtual tmp<vectorField> force
(
const volVectorField& U,
const volScalarField& rho,
const volScalarField& mu
) const;
//- Add resistance
virtual void addResistance(fvVectorMatrix& UEqn) const;
@ -227,10 +250,19 @@ public:
) const;
// Topology change
//- Move points
virtual bool movePoints();
//- Update on meshUpdate
virtual void updateMesh(const mapPolyMesh& mpm);
// I-O
//- Write
virtual void writeData(Ostream& os) const = 0;
virtual bool writeData(Ostream& os) const;
//- Read porosity dictionary
virtual bool read(const dictionary& dict);

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2012-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -35,4 +35,10 @@ inline bool Foam::porosityModel::active() const
}
inline const Foam::labelList& Foam::porosityModel::cellZoneIDs() const
{
return cellZoneIDs_;
}
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2012-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -66,6 +66,23 @@ Foam::porosityModels::powerLaw::~powerLaw()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::porosityModels::powerLaw::calcForce
(
const volVectorField& U,
const volScalarField& rho,
const volScalarField& mu,
vectorField& force
) const
{
scalarField Udiag(U.size());
const scalarField& V = mesh_.V();
apply(Udiag, V, rho, U);
force = Udiag*U;
}
void Foam::porosityModels::powerLaw::correct
(
fvVectorMatrix& UEqn
@ -126,10 +143,12 @@ void Foam::porosityModels::powerLaw::correct
}
void Foam::porosityModels::powerLaw::writeData(Ostream& os) const
bool Foam::porosityModels::powerLaw::writeData(Ostream& os) const
{
os << indent << name_ << endl;
dict_.write(os);
return true;
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2012-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -127,6 +127,15 @@ public:
// Member Functions
//- Calculate the porosity force
virtual void calcForce
(
const volVectorField& U,
const volScalarField& rho,
const volScalarField& mu,
vectorField& force
) const;
//- Add resistance
virtual void correct(fvVectorMatrix& UEqn) const;
@ -149,7 +158,7 @@ public:
// I-O
//- Write
void writeData(Ostream& os) const;
bool writeData(Ostream& os) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2012-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -37,9 +37,9 @@ void Foam::porosityModels::powerLaw::apply
const scalar C0 = C0_;
const scalar C1m1b2 = (C1_ - 1.0)/2.0;
forAll(cellZoneIds_, zoneI)
forAll(cellZoneIDs_, zoneI)
{
const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]];
const labelList& cells = mesh_.cellZones()[cellZoneIDs_[zoneI]];
forAll(cells, i)
{
@ -63,9 +63,9 @@ void Foam::porosityModels::powerLaw::apply
const scalar C0 = C0_;
const scalar C1m1b2 = (C1_ - 1.0)/2.0;
forAll(cellZoneIds_, zoneI)
forAll(cellZoneIDs_, zoneI)
{
const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]];
const labelList& cells = mesh_.cellZones()[cellZoneIDs_[zoneI]];
forAll(cells, i)
{

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -43,7 +43,7 @@ tmp<GeometricField<Type, fvPatchField, volMesh> >
Su
(
const GeometricField<Type, fvPatchField, volMesh>& su,
GeometricField<Type, fvPatchField, volMesh>& vf
const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
return su;
@ -54,7 +54,7 @@ tmp<GeometricField<Type, fvPatchField, volMesh> >
Su
(
const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu,
GeometricField<Type, fvPatchField, volMesh>& vf
const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
return tsu;
@ -66,7 +66,7 @@ tmp<GeometricField<Type, fvPatchField, volMesh> >
Sp
(
const volScalarField& sp,
GeometricField<Type, fvPatchField, volMesh>& vf
const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
return sp*vf;
@ -77,7 +77,7 @@ tmp<GeometricField<Type, fvPatchField, volMesh> >
Sp
(
const tmp<volScalarField>& tsp,
GeometricField<Type, fvPatchField, volMesh>& vf
const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
return tsp*vf;
@ -89,7 +89,7 @@ tmp<GeometricField<Type, fvPatchField, volMesh> >
Sp
(
const dimensionedScalar& sp,
GeometricField<Type, fvPatchField, volMesh>& vf
const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
return sp*vf;
@ -101,7 +101,7 @@ tmp<GeometricField<Type, fvPatchField, volMesh> >
SuSp
(
const volScalarField& sp,
GeometricField<Type, fvPatchField, volMesh>& vf
const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
return sp*vf;
@ -112,7 +112,7 @@ tmp<GeometricField<Type, fvPatchField, volMesh> >
SuSp
(
const tmp<volScalarField>& tsp,
GeometricField<Type, fvPatchField, volMesh>& vf
const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
return tsp*vf;

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -55,14 +55,14 @@ namespace fvc
tmp<GeometricField<Type, fvPatchField, volMesh> > Su
(
const GeometricField<Type, fvPatchField, volMesh>&,
GeometricField<Type, fvPatchField, volMesh>&
const GeometricField<Type, fvPatchField, volMesh>&
);
template<class Type>
tmp<GeometricField<Type, fvPatchField, volMesh> > Su
(
const tmp<GeometricField<Type, fvPatchField, volMesh> >&,
GeometricField<Type, fvPatchField, volMesh>&
const GeometricField<Type, fvPatchField, volMesh>&
);
@ -72,14 +72,14 @@ namespace fvc
tmp<GeometricField<Type, fvPatchField, volMesh> > Sp
(
const volScalarField&,
GeometricField<Type, fvPatchField, volMesh>&
const GeometricField<Type, fvPatchField, volMesh>&
);
template<class Type>
tmp<GeometricField<Type, fvPatchField, volMesh> > Sp
(
const tmp<volScalarField>&,
GeometricField<Type, fvPatchField, volMesh>&
const GeometricField<Type, fvPatchField, volMesh>&
);
@ -87,7 +87,7 @@ namespace fvc
tmp<GeometricField<Type, fvPatchField, volMesh> > Sp
(
const dimensionedScalar&,
GeometricField<Type, fvPatchField, volMesh>&
const GeometricField<Type, fvPatchField, volMesh>&
);
@ -97,14 +97,14 @@ namespace fvc
tmp<GeometricField<Type, fvPatchField, volMesh> > SuSp
(
const volScalarField&,
GeometricField<Type, fvPatchField, volMesh>&
const GeometricField<Type, fvPatchField, volMesh>&
);
template<class Type>
tmp<GeometricField<Type, fvPatchField, volMesh> > SuSp
(
const tmp<volScalarField>&,
GeometricField<Type, fvPatchField, volMesh>&
const GeometricField<Type, fvPatchField, volMesh>&
);
}

View File

@ -120,8 +120,8 @@ void Foam::forceCoeffs::write()
scalar pDyn = 0.5*rhoRef_*magUInf_*magUInf_;
Field<vector> totForce(force_[0] + force_[1]);
Field<vector> totMoment(moment_[0] + moment_[1]);
Field<vector> totForce(force_[0] + force_[1] + force_[2]);
Field<vector> totMoment(moment_[0] + moment_[1] + moment_[2]);
List<Field<scalar> > coeffs(3);
coeffs[0].setSize(nBin_);

View File

@ -37,6 +37,7 @@ License
#include "compressible/RAS/RASModel/RASModel.H"
#include "compressible/LES/LESModel/LESModel.H"
#include "porosityModel.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -52,14 +53,15 @@ void Foam::forces::writeFileHeader(const label i)
{
file()
<< "# Time" << tab
<< "forces(pressure, viscous) moment(pressure, viscous)";
<< "forces(pressure, viscous, porous) "
<< "moment(pressure, viscous, porous)";
if (localSystem_)
{
file()
<< tab
<< "local forces(pressure, viscous) "
<< "local moment(pressure, viscous)";
<< "local forces(pressure, viscous, porous) "
<< "local moment(pressure, viscous, porous)";
}
file()<< endl;
@ -132,7 +134,7 @@ Foam::tmp<Foam::volSymmTensorField> Foam::forces::devRhoReff() const
else
{
FatalErrorIn("forces::devRhoReff()")
<< "No valid model for viscous stress calculation."
<< "No valid model for viscous stress calculation"
<< exit(FatalError);
return volSymmTensorField::null();
@ -140,6 +142,46 @@ Foam::tmp<Foam::volSymmTensorField> Foam::forces::devRhoReff() const
}
Foam::tmp<Foam::volScalarField> Foam::forces::mu() const
{
if (obr_.foundObject<fluidThermo>("thermophysicalProperties"))
{
const fluidThermo& thermo =
obr_.lookupObject<fluidThermo>("thermophysicalProperties");
return thermo.mu();
}
else if
(
obr_.foundObject<singlePhaseTransportModel>("transportProperties")
)
{
const singlePhaseTransportModel& laminarT =
obr_.lookupObject<singlePhaseTransportModel>
("transportProperties");
return rho()*laminarT.nu();
}
else if (obr_.foundObject<dictionary>("transportProperties"))
{
const dictionary& transportProperties =
obr_.lookupObject<dictionary>("transportProperties");
dimensionedScalar nu(transportProperties.lookup("nu"));
return rho()*nu;
}
else
{
FatalErrorIn("forces::mu()")
<< "No valid model for dynamic viscosity calculation"
<< exit(FatalError);
return volScalarField::null();
}
}
Foam::tmp<Foam::volScalarField> Foam::forces::rho() const
{
if (rhoName_ == "rhoInf")
@ -190,32 +232,35 @@ Foam::scalar Foam::forces::rho(const volScalarField& p) const
void Foam::forces::applyBins
(
const label patchI,
const vectorField fN,
const vectorField Md,
const vectorField fT
const vectorField& Md,
const vectorField& fN,
const vectorField& fT,
const vectorField& fP,
const vectorField& d
)
{
if (nBin_ == 1)
{
force_[0][0] += sum(fN);
force_[1][0] += sum(fT);
moment_[0][0] += sum(Md ^ fN);
moment_[1][0] += sum(Md ^ fT);
force_[2][0] += sum(fP);
moment_[0][0] += sum(Md^fN);
moment_[1][0] += sum(Md^fT);
moment_[2][0] += sum(Md^fP);
}
else
{
const fvMesh& mesh = refCast<const fvMesh>(obr_);
const vectorField& Cf = mesh.C().boundaryField()[patchI];
scalarField d((Cf & binDir_) - binMin_);
scalarField dd((d & binDir_) - binMin_);
forAll(d, i)
forAll(dd, i)
{
label binI = floor(d[i]/binDx_);
label binI = floor(dd[i]/binDx_);
force_[0][binI] += fN[i];
force_[1][binI] += fT[i];
force_[2][binI] += fP[i];
moment_[0][binI] += Md[i]^fN[i];
moment_[1][binI] += Md[i]^fT[i];
moment_[2][binI] += Md[i]^fP[i];
}
}
}
@ -239,7 +284,7 @@ void Foam::forces::writeBins() const
Info<< " Writing bins to " << forcesDir << endl;
}
wordList fieldNames(IStringStream("(pressure viscous)")());
wordList fieldNames(IStringStream("(pressure viscous porous)")());
List<Field<vector> > f(force_);
List<Field<vector> > m(moment_);
@ -250,8 +295,11 @@ void Foam::forces::writeBins() const
{
f[0][i] += f[0][i-1];
f[1][i] += f[1][i-1];
f[2][i] += f[2][i-1];
m[0][i] += m[0][i-1];
m[1][i] += m[1][i-1];
m[2][i] += m[2][i-1];
}
}
@ -264,12 +312,14 @@ void Foam::forces::writeBins() const
if (localSystem_)
{
List<Field<vector> > localForce(2);
List<Field<vector> > localMoment(2);
List<Field<vector> > localForce(3);
List<Field<vector> > localMoment(3);
localForce[0] = coordSys_.localVector(force_[0]);
localForce[1] = coordSys_.localVector(force_[1]);
localForce[2] = coordSys_.localVector(force_[2]);
localMoment[0] = coordSys_.localVector(moment_[0]);
localMoment[1] = coordSys_.localVector(moment_[1]);
localMoment[2] = coordSys_.localVector(moment_[2]);
if (binCumulative_)
{
@ -277,8 +327,10 @@ void Foam::forces::writeBins() const
{
localForce[0][i] += localForce[0][i-1];
localForce[1][i] += localForce[1][i-1];
localForce[2][i] += localForce[2][i-1];
localMoment[0][i] += localMoment[0][i-1];
localMoment[1][i] += localMoment[1][i-1];
localMoment[2][i] += localMoment[2][i-1];
}
}
@ -306,8 +358,8 @@ Foam::forces::forces
obr_(obr),
active_(true),
log_(false),
force_(2),
moment_(2),
force_(3),
moment_(3),
patchSet_(),
pName_(word::null),
UName_(word::null),
@ -318,6 +370,7 @@ Foam::forces::forces
pRef_(0),
coordSys_(),
localSystem_(false),
porosity_(false),
nBin_(1),
binDir_(vector::zero),
binDx_(0.0),
@ -365,8 +418,8 @@ Foam::forces::forces
obr_(obr),
active_(true),
log_(false),
force_(2),
moment_(2),
force_(3),
moment_(3),
patchSet_(patchSet),
pName_(pName),
UName_(UName),
@ -377,6 +430,7 @@ Foam::forces::forces
pRef_(pRef),
coordSys_(coordSys),
localSystem_(false),
porosity_(false),
nBin_(1),
binDir_(vector::zero),
binDx_(0.0),
@ -475,6 +529,18 @@ void Foam::forces::read(const dictionary& dict)
localSystem_ = true;
}
if (dict.readIfPresent("porosity", porosity_))
{
if (porosity_)
{
Info<< "Including porosity effects" << endl;
}
else
{
Info<< "Not including porosity effects" << endl;
}
}
if (dict.found("binData"))
{
const dictionary& binDict(dict.subDict("binData"));
@ -491,10 +557,11 @@ void Foam::forces::read(const dictionary& dict)
else if ((nBin_ == 0) || (nBin_ == 1))
{
nBin_ = 1;
force_[0].setSize(1);
force_[1].setSize(1);
moment_[0].setSize(1);
moment_[1].setSize(1);
forAll(force_, i)
{
force_[i].setSize(1);
moment_[i].setSize(1);
}
}
if (nBin_ > 1)
@ -546,8 +613,10 @@ void Foam::forces::read(const dictionary& dict)
// allocate storage for forces and moments
force_[0].setSize(1);
force_[1].setSize(1);
force_[2].setSize(1);
moment_[0].setSize(1);
moment_[1].setSize(1);
moment_[2].setSize(1);
}
}
}
@ -581,28 +650,50 @@ void Foam::forces::write()
if (log_)
{
Info<< type() << " output:" << nl
<< " forces(pressure,viscous)"
<< "(" << sum(force_[0]) << "," << sum(force_[1]) << ")" << nl
<< " moment(pressure,viscous)"
<< "(" << sum(moment_[0]) << "," << sum(moment_[1]) << ")"
<< " forces(pressure,viscous,porous) = ("
<< sum(force_[0]) << ","
<< sum(force_[1]) << ","
<< sum(force_[2]) << ")" << nl
<< " moment(pressure,viscous,porous) = ("
<< sum(moment_[0]) << ","
<< sum(moment_[1]) << ","
<< sum(moment_[2]) << ")"
<< nl;
}
file() << obr_.time().value() << tab
<< "(" << sum(force_[0]) << "," << sum(force_[1]) << ") "
<< "(" << sum(moment_[0]) << "," << sum(moment_[1]) << ")"
<< "("
<< sum(force_[0]) << ","
<< sum(force_[1]) << ","
<< sum(force_[2])
<< ") "
<< "("
<< sum(moment_[0]) << ","
<< sum(moment_[1]) << ","
<< sum(moment_[2])
<< ")"
<< endl;
if (localSystem_)
{
vectorField localForceP(coordSys_.localVector(force_[0]));
vectorField localForceV(coordSys_.localVector(force_[1]));
vectorField localMomentP(coordSys_.localVector(moment_[0]));
vectorField localMomentV(coordSys_.localVector(moment_[1]));
vectorField localForceN(coordSys_.localVector(force_[0]));
vectorField localForceT(coordSys_.localVector(force_[1]));
vectorField localForceP(coordSys_.localVector(force_[2]));
vectorField localMomentN(coordSys_.localVector(moment_[0]));
vectorField localMomentT(coordSys_.localVector(moment_[1]));
vectorField localMomentP(coordSys_.localVector(moment_[2]));
file() << obr_.time().value() << tab
<< "(" << sum(localForceP) << "," << sum(localForceV) << ") "
<< "(" << sum(localMomentP) << "," << sum(localMomentV) << ")"
<< "("
<< sum(localForceN) << ","
<< sum(localForceT) << ","
<< sum(localForceP)
<< ") "
<< "("
<< sum(localMomentN) << ","
<< sum(localMomentT) << ","
<< sum(localMomentP)
<< ")"
<< endl;
}
@ -620,9 +711,11 @@ void Foam::forces::calcForcesMoment()
{
force_[0] = vector::zero;
force_[1] = vector::zero;
force_[2] = vector::zero;
moment_[0] = vector::zero;
moment_[1] = vector::zero;
moment_[2] = vector::zero;
if (directForceDensity_)
{
@ -656,7 +749,10 @@ void Foam::forces::calcForcesMoment()
// Tangential force (total force minus normal fN)
vectorField fT(sA*fD.boundaryField()[patchI] - fN);
applyBins(patchI, fN, Md, fT);
//- Porous force
vectorField fP(Md.size(), vector::zero);
applyBins(Md, fN, fT, fP, mesh.C().boundaryField()[patchI]);
}
}
else
@ -685,14 +781,51 @@ void Foam::forces::calcForcesMoment()
mesh.C().boundaryField()[patchI] - coordSys_.origin()
);
vectorField pf
vectorField fN
(
rho(p)*Sfb[patchI]*(p.boundaryField()[patchI] - pRef)
);
vectorField vf(Sfb[patchI] & devRhoReffb[patchI]);
vectorField fT(Sfb[patchI] & devRhoReffb[patchI]);
applyBins(patchI, pf, Md, vf);
vectorField fP(Md.size(), vector::zero);
applyBins(Md, fN, fT, fP, mesh.C().boundaryField()[patchI]);
}
}
if (porosity_)
{
const volVectorField& U = obr_.lookupObject<volVectorField>(UName_);
const volScalarField rho(this->rho());
const volScalarField mu(this->mu());
const fvMesh& mesh = U.mesh();
const HashTable<const porosityModel*> models =
obr_.lookupClass<porosityModel>();
forAllConstIter(HashTable<const porosityModel*>, models, iter)
{
const porosityModel& pm = *iter();
vectorField fPTot(pm.force(U, rho, mu));
const labelList& cellZoneIDs = pm.cellZoneIDs();
forAll(cellZoneIDs, i)
{
label zoneI = cellZoneIDs[i];
const cellZone& cZone = mesh.cellZones()[zoneI];
const vectorField d(mesh.C(), cZone);
const vectorField fP(fPTot, cZone);
const vectorField Md(d - coordSys_.origin());
const vectorField fDummy(Md.size(), vector::zero);
applyBins(Md, fDummy, fDummy, fP, d);
}
}
}
@ -703,13 +836,13 @@ void Foam::forces::calcForcesMoment()
Foam::vector Foam::forces::forceEff() const
{
return sum(force_[0]) + sum(force_[1]);
return sum(force_[0]) + sum(force_[1]) + sum(force_[2]);
}
Foam::vector Foam::forces::momentEff() const
{
return sum(moment_[0]) + sum(moment_[1]);
return sum(moment_[0]) + sum(moment_[1]) + sum(moment_[2]);
}

View File

@ -149,10 +149,10 @@ protected:
//- Switch to send output to Info as well as to file
Switch log_;
//- Pressure and viscous force per bin
//- Pressure, viscous and porous force per bin
List<Field<vector> > force_;
//- Pressure and viscous pressure per bin
//- Pressure, viscous and porous moment per bin
List<Field<vector> > moment_;
@ -188,6 +188,9 @@ protected:
//- Flag to indicate whether we are using a local co-ordinate sys
bool localSystem_;
//- Flag to include porosity effects
bool porosity_;
// Bin information
@ -221,6 +224,9 @@ protected:
//- Return the effective viscous stress (laminar + turbulent).
tmp<volSymmTensorField> devRhoReff() const;
//- Dynamic viscosity field
tmp<volScalarField> mu() const;
//- Return rho if rhoName is specified otherwise rhoRef
tmp<volScalarField> rho() const;
@ -231,10 +237,11 @@ protected:
//- Accumulate bin data
void applyBins
(
const label patchI,
const vectorField fN,
const vectorField Md,
const vectorField fT
const vectorField& Md,
const vectorField& fN,
const vectorField& fT,
const vectorField& fP,
const vectorField& d
);
//- Helper function to write bin data

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2012-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -241,6 +241,20 @@ Foam::Ostream& Foam::OBJstream::write(const linePointRef& ln)
}
Foam::Ostream& Foam::OBJstream::write
(
const linePointRef& ln,
const vector& n0,
const vector& n1
)
{
write(ln.start(), n0);
write(ln.end(), n1);
write("l ") << nVertices_-1 << ' ' << nVertices_ << nl;
return *this;
}
Foam::Ostream& Foam::OBJstream::write
(
const triPointRef& f,

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2012-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -92,7 +92,7 @@ public:
// Access
//- Return the name of the stream
//- Return the number of vertices written
label nVertices() const
{
return nVertices_;
@ -135,6 +135,14 @@ public:
//- Write line
Ostream& write(const linePointRef&);
//- Write line with points and vector normals ('vn')
Ostream& write
(
const linePointRef&,
const vector& n0,
const vector& n1
);
//- Write triangle as points with lines or filled polygon
Ostream& write(const triPointRef&, const bool lines = true);