Merge remote-tracking branch 'origin/feature/procAgglom'

Conflicts:
	src/OpenFOAM/meshes/lduMesh/lduMesh.H
This commit is contained in:
mattijs
2013-05-02 14:44:59 +01:00
185 changed files with 10067 additions and 1238 deletions

10
README.txt Normal file
View File

@ -0,0 +1,10 @@
2013-01
- to UPstream added allocation of communicators
- added communicators to lduMesh,lduInterface and (indirectly)
to lduInterfaceField
- gamg agglomeration allocates new communicator for every level
- in all linear solvers/smoothers/preconditioners make they use
communicator
- added lots of warnings if using unexpected communicator (UPstream::warnComm)
- did LUScalarMatrix for 'directSolveCoarsest'

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

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

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

View File

@ -0,0 +1,3 @@
laplacianFoam.C
EXE = $(FOAM_USER_APPBIN)/laplacianFoam-communicators

View File

@ -0,0 +1,3 @@
EXE_INC = -I$(LIB_SRC)/finiteVolume/lnInclude
EXE_LIBS = -lfiniteVolume

View File

@ -0,0 +1,37 @@
Info<< "Reading field T\n" << endl;
volScalarField T
(
IOobject
(
"T",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading transportProperties\n" << endl;
IOdictionary transportProperties
(
IOobject
(
"transportProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
);
Info<< "Reading diffusivity DT\n" << endl;
dimensionedScalar DT
(
transportProperties.lookup("DT")
);

View File

@ -0,0 +1,817 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 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/>.
Application
laplacianFoam
Description
Solves a simple Laplace equation, e.g. for thermal diffusion in a solid.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "simpleControl.H"
#include "globalIndex.H"
#include "lduPrimitiveMesh.H"
#include "processorGAMGInterface.H"
#include "GAMGInterfaceField.H"
#include "processorLduInterfaceField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void checkUpperTriangular
(
const label size,
const labelUList& l,
const labelUList& u
)
{
forAll(l, faceI)
{
if (u[faceI] < l[faceI])
{
FatalErrorIn
(
"checkUpperTriangular"
"(const label, const labelUList&, const labelUList&)"
) << "Reversed face. Problem at face " << faceI
<< " l:" << l[faceI] << " u:" << u[faceI] << abort(FatalError);
}
if (l[faceI] < 0 || u[faceI] < 0 || u[faceI] >= size)
{
FatalErrorIn
(
"checkUpperTriangular"
"(const label, const labelUList&, const labelUList&)"
) << "Illegal cell label. Problem at face " << faceI
<< " l:" << l[faceI] << " u:" << u[faceI] << abort(FatalError);
}
}
for (label faceI=1; faceI < l.size(); faceI++)
{
if (l[faceI-1] > l[faceI])
{
FatalErrorIn
(
"checkUpperTriangular"
"(const label, const labelUList&, const labelUList&)"
) << "Lower not in incremental cell order."
<< " Problem at face " << faceI
<< " l:" << l[faceI] << " u:" << u[faceI]
<< " previous l:" << l[faceI-1] << abort(FatalError);
}
else if (l[faceI-1] == l[faceI])
{
// Same cell.
if (u[faceI-1] > u[faceI])
{
FatalErrorIn
(
"checkUpperTriangular"
"(const label, const labelUList&, const labelUList&)"
) << "Upper not in incremental cell order."
<< " Problem at face " << faceI
<< " l:" << l[faceI] << " u:" << u[faceI]
<< " previous u:" << u[faceI-1] << abort(FatalError);
}
}
}
}
void print(const string& msg, const lduMesh& mesh)
{
const lduAddressing& addressing = mesh.lduAddr();
const lduInterfacePtrsList interfaces = mesh.interfaces();
Pout<< "Mesh:" << msg.c_str() << nl
<< " cells:" << addressing.size() << nl
<< " faces:" << addressing.lowerAddr().size() << nl
<< " patches:" << interfaces.size() << nl;
const labelUList& l = addressing.lowerAddr();
const labelUList& startAddr = addressing.losortStartAddr();
const labelUList& addr = addressing.losortAddr();
forAll(addressing, cellI)
{
Pout<< " cell:" << cellI << nl;
label start = startAddr[cellI];
label end = startAddr[cellI+1];
for (label index = start; index < end; index++)
{
Pout<< " nbr:" << l[addr[index]] << nl;
}
}
Pout<< " Patches:" << nl;
forAll(interfaces, i)
{
if (interfaces.set(i))
{
if (isA<processorLduInterface>(interfaces[i]))
{
const processorLduInterface& pldui =
refCast<const processorLduInterface>(interfaces[i]);
Pout<< " " << i
<< " me:" << pldui.myProcNo()
<< " nbr:" << pldui.neighbProcNo()
<< " comm:" << pldui.comm()
<< " tag:" << pldui.tag()
<< nl;
}
{
Pout<< " " << i << " addressing:" << nl;
const labelUList& faceCells = interfaces[i].faceCells();
forAll(faceCells, i)
{
Pout<< "\t\t" << i << '\t' << faceCells[i] << nl;
}
}
}
}
}
template<class ProcPatch>
lduSchedule nonBlockingSchedule
(
const lduInterfacePtrsList& interfaces
)
{
lduSchedule schedule(2*interfaces.size());
label slotI = 0;
forAll(interfaces, i)
{
if (interfaces.set(i) && !isA<ProcPatch>(interfaces[i]))
{
schedule[slotI].patch = i;
schedule[slotI].init = true;
slotI++;
schedule[slotI].patch = i;
schedule[slotI].init = false;
slotI++;
}
}
forAll(interfaces, i)
{
if (interfaces.set(i) && isA<ProcPatch>(interfaces[i]))
{
schedule[slotI].patch = i;
schedule[slotI].init = true;
slotI++;
}
}
forAll(interfaces, i)
{
if (interfaces.set(i) && isA<ProcPatch>(interfaces[i]))
{
schedule[slotI].patch = i;
schedule[slotI].init = false;
slotI++;
}
}
return schedule;
}
void sendReceive
(
const label comm,
const label tag,
const globalIndex& offsets,
const scalarField& field,
scalarField& allField
)
{
label nProcs = Pstream::nProcs(comm);
if (Pstream::master(comm))
{
allField.setSize(offsets.size());
// Assign master slot
SubList<scalar>
(
allField,
offsets.localSize(0),
offsets.offset(0)
).assign(field);
// Assign slave slots
for (label procI = 1; procI < nProcs; procI++)
{
SubList<scalar> procSlot
(
allField,
offsets.localSize(procI),
offsets.offset(procI)
);
Pout<< "Receiving allField from " << procI
<< " at offset:" << offsets.offset(procI)
<< " size:" << offsets.size()
<< endl;
IPstream::read
(
Pstream::nonBlocking,
procI,
reinterpret_cast<char*>(procSlot.begin()),
procSlot.byteSize(),
tag,
comm
);
}
}
else
{
OPstream::write
(
Pstream::nonBlocking,
0, // master
reinterpret_cast<const char*>(field.begin()),
field.byteSize(),
tag,
comm
);
}
}
void sendReceive
(
const label comm,
const label tag,
const globalIndex& offsets,
const FieldField<Field, scalar>& field,
FieldField<Field, scalar>& allField
)
{
PstreamBuffers pBufs(Pstream::nonBlocking, Pstream::msgType(), comm);
if (!Pstream::master(comm))
{
UOPstream toMaster(Pstream::masterNo(), pBufs);
Pout<< "To 0 sending " << field.size()
<< " fields." << endl;
forAll(field, intI)
{
toMaster << field[intI];
}
}
pBufs.finishedSends();
if (Pstream::master(comm))
{
allField.setSize(offsets.size());
forAll(allField, i)
{
allField.set(i, new scalarField(0));
}
// Insert master values
forAll(field, intI)
{
allField[intI] = field[intI];
}
// Receive and insert slave values
label nProcs = Pstream::nProcs(comm);
for (label procI = 1; procI < nProcs; procI++)
{
UIPstream fromSlave(procI, pBufs);
label nSlaveInts = offsets.localSize(procI);
Pout<< "From " << procI << " receiving "
<< nSlaveInts << " fields." << endl;
for (label intI = 0; intI < nSlaveInts; intI++)
{
label slotI = offsets.toGlobal(procI, intI);
Pout<< " int:" << intI << " goes into slot " << slotI
<< endl;
fromSlave >> allField[slotI];
}
}
}
}
void collect
(
const label comm,
const globalIndex& cellOffsets,
const globalIndex& faceOffsets,
const scalarField& diagonal,
const scalarField& upper,
const scalarField& lower,
scalarField& allDiagonal,
scalarField& allUpper,
scalarField& allLower
)
{
label nOutstanding = Pstream::nRequests();
int allDiagonalTag = Pstream::allocateTag("allDiagonal:" __FILE__);
int allUpperTag = Pstream::allocateTag("allUpper:" __FILE__);
int allLowerTag = Pstream::allocateTag("allLower:" __FILE__);
sendReceive
(
comm,
allDiagonalTag,
cellOffsets,
diagonal,
allDiagonal
);
sendReceive
(
comm,
allUpperTag,
faceOffsets,
upper,
allUpper
);
sendReceive
(
comm,
allLowerTag,
faceOffsets,
lower,
allLower
);
Pstream::waitRequests(nOutstanding);
Pstream::freeTag("allDiagonal:" __FILE__, allDiagonalTag);
Pstream::freeTag("allUpper:" __FILE__, allUpperTag);
Pstream::freeTag("allLower:" __FILE__, allLowerTag);
}
void setCommunicator(fvMesh& mesh, const label newComm)
{
const polyBoundaryMesh& pbm = mesh.boundaryMesh();
// The current state is consistent with the mesh so check where the new
// communicator is and adjust accordingly.
forAll(pbm, patchI)
{
if (isA<processorPolyPatch>(pbm[patchI]))
{
processorPolyPatch& ppp = const_cast<processorPolyPatch&>
(
refCast
<
const processorPolyPatch
>(pbm[patchI])
);
label thisRank = UPstream::procNo
(
newComm,
ppp.comm(),
ppp.myProcNo()
);
label nbrRank = UPstream::procNo
(
newComm,
ppp.comm(),
ppp.neighbProcNo()
);
//ppp.comm() = newComm;
ppp.myProcNo() = thisRank;
ppp.neighbProcNo() = nbrRank;
}
}
mesh.polyMesh::comm() = newComm;
}
namespace Foam
{
typedef UPtrList<const GAMGInterfaceField> GAMGInterfaceFieldPtrsList;
}
// Gather matrices from processors procIDs[1..] on procIDs[0]
void gatherMatrices
(
const labelList& procIDs,
const PtrList<lduMesh>& procMeshes,
const lduMatrix& mat,
const FieldField<Field, scalar>& interfaceBouCoeffs,
const FieldField<Field, scalar>& interfaceIntCoeffs,
const lduInterfaceFieldPtrsList& interfaces,
PtrList<lduMatrix>& otherMats,
PtrList<FieldField<Field, scalar> >& otherBouCoeffs,
PtrList<FieldField<Field, scalar> >& otherIntCoeffs,
PtrList<GAMGInterfaceFieldPtrsList>& otherInterfaces
)
{
const label meshComm = mat.mesh().comm();
//lduInterfacePtrsList interfaces(mesh.interfaces());
if (Pstream::myProcNo(meshComm) == procIDs[0])
{
// Master.
otherMats.setSize(procIDs.size()-1);
otherBouCoeffs.setSize(procIDs.size()-1);
otherIntCoeffs.setSize(procIDs.size()-1);
otherInterfaces.setSize(procIDs.size()-1);
for (label i = 1; i < procIDs.size(); i++)
{
const lduMesh& procMesh = procMeshes[i];
const lduInterfacePtrsList procInterfaces = procMesh.interfaces();
IPstream fromSlave
(
Pstream::scheduled,
procIDs[i],
0, // bufSize
Pstream::msgType(),
meshComm
);
otherMats.set(i-1, new lduMatrix(procMesh, fromSlave));
// Receive number of/valid interfaces
boolList validTransforms(fromSlave);
List<int> validRanks(fromSlave);
// Size coefficients
otherBouCoeffs.set
(
i-1,
new FieldField<Field, scalar>(validTransforms.size())
);
otherIntCoeffs.set
(
i-1,
new FieldField<Field, scalar>(validTransforms.size())
);
otherInterfaces.set
(
i-1,
new GAMGInterfaceFieldPtrsList(validTransforms.size())
);
forAll(validTransforms, intI)
{
if (validTransforms[intI])
{
const processorGAMGInterface& interface =
refCast<const processorGAMGInterface>
(
procInterfaces[intI]
);
otherBouCoeffs[i-1].set(intI, new scalarField(fromSlave));
otherIntCoeffs[i-1].set(intI, new scalarField(fromSlave));
otherInterfaces[i-1].set
(
intI,
GAMGInterfaceField::New
(
interface, //procInterfaces[intI],
validTransforms[intI],
validRanks[intI]
).ptr()
);
}
}
}
}
else
{
// Send to master
OPstream toMaster
(
Pstream::scheduled,
procIDs[0],
0,
Pstream::msgType(),
meshComm
);
// Count valid interfaces
boolList validTransforms(interfaceBouCoeffs.size(), false);
List<int> validRanks(interfaceBouCoeffs.size(), -1);
forAll(interfaces, intI)
{
if (interfaces.set(intI))
{
const processorLduInterfaceField& interface =
refCast<const processorLduInterfaceField>
(
interfaces[intI]
);
validTransforms[intI] = interface.doTransform();
validRanks[intI] = interface.rank();
}
}
toMaster << mat << validTransforms << validRanks;
forAll(validTransforms, intI)
{
if (validTransforms[intI])
{
toMaster
<< interfaceBouCoeffs[intI]
<< interfaceIntCoeffs[intI];
}
}
}
}
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "createFields.H"
simpleControl simple(mesh);
//const polyBoundaryMesh& pbm = mesh.boundaryMesh();
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nCalculating temperature distribution\n" << endl;
// Get a subset of processors
labelList subProcs(3);
subProcs[0] = 0;
subProcs[1] = 1;
subProcs[2] = 2;
const UPstream::communicator newComm
(
UPstream::worldComm,
subProcs,
true
);
Pout<< "procIDs world :" << UPstream::procID(UPstream::worldComm) << endl;
Pout<< "procIDs newComm:" << UPstream::procID(newComm) << endl;
//// On ALL processors: get the interfaces:
//lduInterfacePtrsList interfaces(mesh.interfaces());
//PtrList<lduMesh> procMeshes;
//
//if (Pstream::myProcNo(newComm) != -1)
//{
// print("InitialMesh", mesh);
//
// labelList procIDs(3);
// procIDs[0] = 0;
// procIDs[1] = 1;
// procIDs[2] = 2;
//
////XXXXXX
// // Collect meshes from procs 0,1 (in newComm) on 1.
// lduPrimitiveMesh::gather(mesh, procIDs, procMeshes);
//
// if (Pstream::myProcNo(newComm) == procIDs[0])
// {
// // Print a bit
// forAll(procMeshes, i)
// {
// const lduMesh& pMesh = procMeshes[i];
// print("procMesh" + Foam::name(i), pMesh);
//
// const lduAddressing& addr = pMesh.lduAddr();
// checkUpperTriangular
// (
// addr.size(),
// addr.lowerAddr(),
// addr.upperAddr()
// );
// }
//
//
// // Combine
// labelList cellOffsets;
// labelListList faceMap;
// labelListList boundaryMap;
// labelListListList boundaryFaceMap;
// //autoPtr<lduPrimitiveMesh> allMeshPtr = combineMeshes
// //(
// // newComm,
// // procIDs,
// // procMeshes,
// //
// // cellOffsets, // per mesh the starting cell
// // faceMap, // per mesh the starting face
// // boundaryMap, // per mesh,per interface the starting face
// // boundaryFaceMap
// //);
// //const lduPrimitiveMesh& allMesh = allMeshPtr();
// const lduPrimitiveMesh allMesh
// (
// newComm,
// procIDs,
// procMeshes,
//
// cellOffsets,
// faceMap,
// boundaryMap,
// boundaryFaceMap
// );
//
//
// print("ALLMESH", allMesh);
//
// forAll(procMeshes, procMeshI)
// {
// const lduMesh& pMesh = procMeshes[procMeshI];
// //const lduAddressing& pAddressing = pMesh.lduAddr();
//
// Pout<< "procMesh" << procMeshI << endl
// << " cells start at:" << cellOffsets[procMeshI] << endl
// << " faces to to:" << faceMap[procMeshI] << endl;
//
// lduInterfacePtrsList interfaces = pMesh.interfaces();
// forAll(interfaces, intI)
// {
// Pout<< " patch:" << intI
// << " becomes patch:" << boundaryMap[procMeshI][intI]
// << endl;
//
// Pout<< " patch:" << intI
// << " faces become faces:"
// << boundaryFaceMap[procMeshI][intI]
// << endl;
// }
// }
// }
//
//
// // Construct test data
// // ~~~~~~~~~~~~~~~~~~~
//
// GAMGInterfaceFieldPtrsList interfaces(interfaces.size());
// FieldField<Field, scalar> interfaceBouCoeffs(interfaces.size());
// FieldField<Field, scalar> interfaceIntCoeffs(interfaces.size());
//
// forAll(interfaces, intI)
// {
// if (interfaces.set(intI))
// {
// label size = interfaces[intI].size();
//
// interfaces.set
// (
// intI,
// GAMGInterfaceField::New
// (
// mesh.interfaces()[intI],
// interfaces[intI]
// )
// );
// interfaceBouCoeffs.set(intI, new scalarField(size, 111));
// interfaceIntCoeffs.set(intI, new scalarField(size, 222));
// }
// }
//
//
// PtrList<lduMatrix> otherMats;
// PtrList<FieldField<Field, scalar> > otherBouCoeffs;
// PtrList<FieldField<Field, scalar> > otherIntCoeffs;
// PtrList<GAMGInterfaceFieldPtrsList> otherInterfaces;
// gatherMatrices
// (
// procIDs,
// procMeshes,
//
// mat,
// interfaceBouCoeffs,
// interfaceIntCoeffs,
// interfaces,
//
// otherMats,
// otherBouCoeffs,
// otherIntCoeffs,
// otherInterfaces
// );
////XXXXXX
//}
{
Pout<< "World:" << UPstream::worldComm
<< " procID:" << 2
<< " subComm:" << newComm
<< " rank1:" << UPstream::procNo(newComm, UPstream::worldComm, 1)
<< " rank2:" << UPstream::procNo(newComm, UPstream::worldComm, 2)
<< endl;
}
while (simple.loop())
{
Info<< "Time = " << runTime.timeName() << nl << endl;
while (simple.correctNonOrthogonal())
{
fvScalarMatrix Teqn
(
//fvm::ddt(T) - fvm::laplacian(DT, T)
fvm::laplacian(DT, T)
);
{
label oldWarn = UPstream::warnComm;
UPstream::warnComm = newComm;
label oldComm = mesh.comm();
setCommunicator(mesh, newComm);
Pout<< "** oldcomm:" << oldComm
<< " newComm:" << mesh.comm() << endl;
if (Pstream::myProcNo(mesh.comm()) != -1)
{
solve(Teqn);
}
setCommunicator(mesh, oldComm);
Pout<< "** reset mesh to:" << mesh.comm() << endl;
UPstream::warnComm = oldWarn;
}
}
#include "write.H"
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Pout<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

View File

@ -0,0 +1,46 @@
if (runTime.outputTime())
{
volVectorField gradT(fvc::grad(T));
volScalarField gradTx
(
IOobject
(
"gradTx",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
gradT.component(vector::X)
);
volScalarField gradTy
(
IOobject
(
"gradTy",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
gradT.component(vector::Y)
);
volScalarField gradTz
(
IOobject
(
"gradTz",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
gradT.component(vector::Z)
);
runTime.write();
}

View File

@ -0,0 +1,3 @@
Test-parallel-communicators.C
EXE = $(FOAM_USER_APPBIN)/Test-parallel-communicators

View File

@ -0,0 +1,203 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 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/>.
Application
Test-parallel-communicators
Description
Checks communication using user-defined communicators
\*---------------------------------------------------------------------------*/
#include "argList.H"
#include "Time.H"
#include "IPstream.H"
#include "OPstream.H"
#include "vector.H"
#include "IOstreams.H"
#include "PstreamReduceOps.H"
using namespace Foam;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
scalar sumReduce
(
const label comm,
const scalar localValue
)
{
scalar sum;
if (Pstream::parRun())
{
if (UPstream::master(comm))
{
// Add master value and all slaves
sum = localValue;
for
(
int slave=Pstream::firstSlave();
slave<=Pstream::lastSlave(comm);
slave++
)
{
scalar slaveValue;
UIPstream::read
(
Pstream::blocking,
slave,
reinterpret_cast<char*>(&slaveValue),
sizeof(scalar),
UPstream::msgType(), // tag
comm // communicator
);
sum += slaveValue;
}
// Send back to slaves
for
(
int slave=UPstream::firstSlave();
slave<=UPstream::lastSlave(comm);
slave++
)
{
UOPstream::write
(
UPstream::blocking,
slave,
reinterpret_cast<const char*>(&sum),
sizeof(scalar),
UPstream::msgType(), // tag
comm // communicator
);
}
}
else
{
{
UOPstream::write
(
UPstream::blocking,
UPstream::masterNo(),
reinterpret_cast<const char*>(&localValue),
sizeof(scalar),
UPstream::msgType(), // tag
comm // communicator
);
}
{
UIPstream::read
(
UPstream::blocking,
UPstream::masterNo(),
reinterpret_cast<char*>(&sum),
sizeof(scalar),
UPstream::msgType(), // tag
comm // communicator
);
}
}
}
return sum;
}
int main(int argc, char *argv[])
{
# include "setRootCase.H"
# include "createTime.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Allocate a communicator
label n = Pstream::nProcs(UPstream::worldComm);
DynamicList<label> bottom;
DynamicList<label> top;
for (label i = 0; i < n/2; i++)
{
bottom.append(i);
}
for (label i = n/2; i < n; i++)
{
top.append(i);
}
//Pout<< "bottom:" << bottom << endl;
Pout<< "top :" << top << endl;
scalar localValue = 111*UPstream::myProcNo(UPstream::worldComm);
Pout<< "localValue :" << localValue << endl;
label comm = Pstream::allocateCommunicator
(
UPstream::worldComm,
top
);
Pout<< "allocated comm :" << comm << endl;
Pout<< "comm myproc :" << Pstream::myProcNo(comm)
<< endl;
if (Pstream::myProcNo(comm) != -1)
{
//scalar sum = sumReduce(comm, localValue);
//scalar sum = localValue;
//reduce
//(
// UPstream::treeCommunication(comm),
// sum,
// sumOp<scalar>(),
// Pstream::msgType(),
// comm
//);
scalar sum = returnReduce
(
localValue,
sumOp<scalar>(),
Pstream::msgType(),
comm
);
Pout<< "sum :" << sum << endl;
}
Pstream::freeCommunicator(comm);
Pout<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

View File

@ -493,6 +493,7 @@ int main(int argc, char *argv[])
mesh.nFaces(), // start
patchI, // index
mesh.boundaryMesh(),// polyBoundaryMesh
Pstream::worldComm, // communicator
Pstream::myProcNo(),// myProcNo
nbrProcI // neighbProcNo
)

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
@ -448,6 +448,7 @@ bool Foam::domainDecomposition::writeDecomposition()
curStart,
nPatches,
procMesh.boundaryMesh(),
Pstream::worldComm,
procI,
curNeighbourProcessors[procPatchI]
);
@ -475,6 +476,7 @@ bool Foam::domainDecomposition::writeDecomposition()
curStart,
nPatches,
procMesh.boundaryMesh(),
Pstream::worldComm,
procI,
curNeighbourProcessors[procPatchI],
referPatch,

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

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

View File

@ -311,12 +311,34 @@ $(GAMGAgglomeration)/GAMGAgglomerateLduAddressing.C
pairGAMGAgglomeration = $(GAMGAgglomerations)/pairGAMGAgglomeration
$(pairGAMGAgglomeration)/pairGAMGAgglomeration.C
$(pairGAMGAgglomeration)/pairGAMGAgglomerate.C
$(pairGAMGAgglomeration)/pairGAMGAgglomerationCombineLevels.C
algebraicPairGAMGAgglomeration = $(GAMGAgglomerations)/algebraicPairGAMGAgglomeration
$(algebraicPairGAMGAgglomeration)/algebraicPairGAMGAgglomeration.C
dummyAgglomeration = $(GAMGAgglomerations)/dummyAgglomeration
$(dummyAgglomeration)/dummyAgglomeration.C
GAMGProcAgglomerations = $(GAMG)/GAMGProcAgglomerations
GAMGProcAgglomeration = $(GAMGProcAgglomerations)/GAMGProcAgglomeration
$(GAMGProcAgglomeration)/GAMGProcAgglomeration.C
masterCoarsestGAMGProcAgglomeration = $(GAMGProcAgglomerations)/masterCoarsestGAMGProcAgglomeration
$(masterCoarsestGAMGProcAgglomeration)/masterCoarsestGAMGProcAgglomeration.C
manualGAMGProcAgglomeration = $(GAMGProcAgglomerations)/manualGAMGProcAgglomeration
$(manualGAMGProcAgglomeration)/manualGAMGProcAgglomeration.C
eagerGAMGProcAgglomeration = $(GAMGProcAgglomerations)/eagerGAMGProcAgglomeration
$(eagerGAMGProcAgglomeration)/eagerGAMGProcAgglomeration.C
noneGAMGProcAgglomeration = $(GAMGProcAgglomerations)/noneGAMGProcAgglomeration
$(noneGAMGProcAgglomeration)/noneGAMGProcAgglomeration.C
/*
cellFaceRatioGAMGProcAgglomeration = $(GAMGProcAgglomerations)/cellFaceRatioGAMGProcAgglomeration
$(cellFaceRatioGAMGProcAgglomeration)/cellFaceRatioGAMGProcAgglomeration.C
*/
meshes/lduMesh/lduMesh.C
meshes/lduMesh/lduPrimitiveMesh.C
LduMatrix = matrices/LduMatrix
$(LduMatrix)/LduMatrix/lduMatrices.C

View File

@ -58,6 +58,7 @@ template<class T> class SubList;
// Forward declaration of friend functions and operators
template<class T> class UList;
template<class T> Ostream& operator<<(Ostream&, const UList<T>&);
template<class T> Istream& operator>>(Istream&, UList<T>&);
typedef UList<label> labelUList;
@ -349,6 +350,14 @@ public:
Ostream&,
const UList<T>&
);
//- Read UList contents from Istream. Requires size to have been set
// before.
friend Istream& operator>> <T>
(
Istream&,
UList<T>&
);
};
template<class T>

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,6 +26,7 @@ License
#include "UList.H"
#include "Ostream.H"
#include "token.H"
#include "SLList.H"
#include "contiguous.H"
// * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
@ -137,4 +138,155 @@ Foam::Ostream& Foam::operator<<(Foam::Ostream& os, const Foam::UList<T>& L)
}
template<class T>
Foam::Istream& Foam::operator>>(Istream& is, UList<T>& L)
{
is.fatalCheck("operator>>(Istream&, UList<T>&)");
token firstToken(is);
is.fatalCheck("operator>>(Istream&, UList<T>&) : reading first token");
if (firstToken.isCompound())
{
List<T> elems;
elems.transfer
(
dynamicCast<token::Compound<List<T> > >
(
firstToken.transferCompoundToken(is)
)
);
// Check list length
label s = elems.size();
if (s != L.size())
{
FatalIOErrorIn("operator>>(Istream&, UList<T>&)", is)
<< "incorrect length for UList. Read " << s
<< " expected " << L.size()
<< exit(FatalIOError);
}
for (register label i=0; i<s; i++)
{
L[i] = elems[i];
}
}
else if (firstToken.isLabel())
{
label s = firstToken.labelToken();
// Set list length to that read
if (s != L.size())
{
FatalIOErrorIn("operator>>(Istream&, UList<T>&)", is)
<< "incorrect length for UList. Read " << s
<< " expected " << L.size()
<< exit(FatalIOError);
}
// Read list contents depending on data format
if (is.format() == IOstream::ASCII || !contiguous<T>())
{
// Read beginning of contents
char delimiter = is.readBeginList("List");
if (s)
{
if (delimiter == token::BEGIN_LIST)
{
for (register label i=0; i<s; i++)
{
is >> L[i];
is.fatalCheck
(
"operator>>(Istream&, UList<T>&) : reading entry"
);
}
}
else
{
T element;
is >> element;
is.fatalCheck
(
"operator>>(Istream&, UList<T>&) : "
"reading the single entry"
);
for (register label i=0; i<s; i++)
{
L[i] = element;
}
}
}
// Read end of contents
is.readEndList("List");
}
else
{
if (s)
{
is.read(reinterpret_cast<char*>(L.data()), s*sizeof(T));
is.fatalCheck
(
"operator>>(Istream&, UList<T>&) : reading the binary block"
);
}
}
}
else if (firstToken.isPunctuation())
{
if (firstToken.pToken() != token::BEGIN_LIST)
{
FatalIOErrorIn("operator>>(Istream&, UList<T>&)", is)
<< "incorrect first token, expected '(', found "
<< firstToken.info()
<< exit(FatalIOError);
}
// Putback the opening bracket
is.putBack(firstToken);
// Now read as a singly-linked list
SLList<T> sll(is);
if (sll.size() != L.size())
{
FatalIOErrorIn("operator>>(Istream&, UList<T>&)", is)
<< "incorrect length for UList. Read " << sll.size()
<< " expected " << L.size()
<< exit(FatalIOError);
}
// Convert the singly-linked list to this list
label i = 0;
for
(
typename SLList<T>::const_iterator iter = sll.begin();
iter != sll.end();
++iter
)
{
L[i] = iter();
}
}
else
{
FatalIOErrorIn("operator>>(Istream&, UList<T>&)", is)
<< "incorrect first token, expected <int> or '(', found "
<< firstToken.info()
<< exit(FatalIOError);
}
return is;
}
// ************************************************************************* //

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
@ -78,9 +78,10 @@ void Foam::IOdictionary::readFile(const bool masterOnly)
(
comms,
const_cast<word&>(headerClassName()),
Pstream::msgType()
Pstream::msgType(),
Pstream::worldComm
);
Pstream::scatter(comms, note(), Pstream::msgType());
Pstream::scatter(comms, note(), Pstream::msgType(), Pstream::worldComm);
// Get my communication order
const Pstream::commsStruct& myComm = comms[Pstream::myProcNo()];

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
@ -33,6 +33,7 @@ Foam::IPstream::IPstream
const int fromProcNo,
const label bufSize,
const int tag,
const label comm,
streamFormat format,
versionNumber version
)
@ -45,6 +46,7 @@ Foam::IPstream::IPstream
buf_,
externalBufPosition_,
tag, // tag
comm,
false, // do not clear buf_ if at end
format,
version

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
@ -69,6 +69,7 @@ public:
const int fromProcNo,
const label bufSize = 0,
const int tag = UPstream::msgType(),
const label comm = UPstream::worldComm,
streamFormat format=BINARY,
versionNumber version=currentVersion
);

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
@ -33,12 +33,13 @@ Foam::OPstream::OPstream
const int toProcNo,
const label bufSize,
const int tag,
const label comm,
streamFormat format,
versionNumber version
)
:
Pstream(commsType, bufSize),
UOPstream(commsType, toProcNo, buf_, tag, true, format, version)
UOPstream(commsType, toProcNo, buf_, tag, comm, true, format, version)
{}

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
@ -66,6 +66,7 @@ public:
const int toProcNo,
const label bufSize = 0,
const int tag = UPstream::msgType(),
const label comm = UPstream::worldComm,
streamFormat format=BINARY,
versionNumber version=currentVersion
);

View File

@ -98,7 +98,8 @@ public:
const List<commsStruct>& comms,
T& Value,
const BinaryOp& bop,
const int tag
const int tag,
const label comm
);
//- Like above but switches between linear/tree communication
@ -107,7 +108,8 @@ public:
(
T& Value,
const BinaryOp& bop,
const int tag = Pstream::msgType()
const int tag = Pstream::msgType(),
const label comm = Pstream::worldComm
);
//- Scatter data. Distribute without modification. Reverse of gather
@ -116,13 +118,18 @@ public:
(
const List<commsStruct>& comms,
T& Value,
const int tag
const int tag,
const label comm
);
//- Like above but switches between linear/tree communication
template<class T>
static void scatter(T& Value, const int tag = Pstream::msgType());
template <class T>
static void scatter
(
T& Value,
const int tag = Pstream::msgType(),
const label comm = Pstream::worldComm
);
// Combine variants. Inplace combine values from processors.
// (Uses construct from Istream instead of <<)
@ -133,7 +140,8 @@ public:
const List<commsStruct>& comms,
T& Value,
const CombineOp& cop,
const int tag
const int tag,
const label comm
);
//- Like above but switches between linear/tree communication
@ -142,7 +150,8 @@ public:
(
T& Value,
const CombineOp& cop,
const int tag = Pstream::msgType()
const int tag = Pstream::msgType(),
const label comm = Pstream::worldComm
);
//- Scatter data. Reverse of combineGather
@ -151,7 +160,8 @@ public:
(
const List<commsStruct>& comms,
T& Value,
const int tag
const int tag,
const label comm
);
//- Like above but switches between linear/tree communication
@ -159,7 +169,8 @@ public:
static void combineScatter
(
T& Value,
const int tag = Pstream::msgType()
const int tag = Pstream::msgType(),
const label comm = Pstream::worldComm
);
// Combine variants working on whole List at a time.
@ -170,7 +181,8 @@ public:
const List<commsStruct>& comms,
List<T>& Value,
const CombineOp& cop,
const int tag
const int tag,
const label comm
);
//- Like above but switches between linear/tree communication
@ -179,7 +191,8 @@ public:
(
List<T>& Value,
const CombineOp& cop,
const int tag = Pstream::msgType()
const int tag = Pstream::msgType(),
const label comm = Pstream::worldComm
);
//- Scatter data. Reverse of combineGather
@ -188,7 +201,8 @@ public:
(
const List<commsStruct>& comms,
List<T>& Value,
const int tag
const int tag,
const label comm
);
//- Like above but switches between linear/tree communication
@ -196,7 +210,8 @@ public:
static void listCombineScatter
(
List<T>& Value,
const int tag = Pstream::msgType()
const int tag = Pstream::msgType(),
const label comm = Pstream::worldComm
);
// Combine variants working on whole map at a time. Container needs to
@ -208,7 +223,8 @@ public:
const List<commsStruct>& comms,
Container& Values,
const CombineOp& cop,
const int tag
const int tag,
const label comm
);
//- Like above but switches between linear/tree communication
@ -217,7 +233,8 @@ public:
(
Container& Values,
const CombineOp& cop,
const int tag = Pstream::msgType()
const int tag = Pstream::msgType(),
const label comm = UPstream::worldComm
);
//- Scatter data. Reverse of combineGather
@ -226,7 +243,8 @@ public:
(
const List<commsStruct>& comms,
Container& Values,
const int tag
const int tag,
const label comm
);
//- Like above but switches between linear/tree communication
@ -234,7 +252,8 @@ public:
static void mapCombineScatter
(
Container& Values,
const int tag = Pstream::msgType()
const int tag = Pstream::msgType(),
const label comm = UPstream::worldComm
);
@ -249,7 +268,8 @@ public:
(
const List<commsStruct>& comms,
List<T>& Values,
const int tag
const int tag,
const label comm
);
//- Like above but switches between linear/tree communication
@ -257,7 +277,8 @@ public:
static void gatherList
(
List<T>& Values,
const int tag = Pstream::msgType()
const int tag = Pstream::msgType(),
const label comm = UPstream::worldComm
);
//- Scatter data. Reverse of gatherList
@ -266,7 +287,8 @@ public:
(
const List<commsStruct>& comms,
List<T>& Values,
const int tag
const int tag,
const label comm
);
//- Like above but switches between linear/tree communication
@ -274,7 +296,8 @@ public:
static void scatterList
(
List<T>& Values,
const int tag = Pstream::msgType()
const int tag = Pstream::msgType(),
const label comm = UPstream::worldComm
);
@ -291,6 +314,7 @@ public:
List<Container >&,
labelListList& sizes,
const int tag = UPstream::msgType(),
const label comm = UPstream::worldComm,
const bool block = true
);

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
@ -40,17 +40,19 @@ Foam::PstreamBuffers::PstreamBuffers
(
const UPstream::commsTypes commsType,
const int tag,
const label comm,
IOstream::streamFormat format,
IOstream::versionNumber version
)
:
commsType_(commsType),
tag_(tag),
comm_(comm),
format_(format),
version_(version),
sendBuf_(UPstream::nProcs()),
recvBuf_(UPstream::nProcs()),
recvBufPos_(UPstream::nProcs(), 0),
sendBuf_(UPstream::nProcs(comm)),
recvBuf_(UPstream::nProcs(comm)),
recvBufPos_(UPstream::nProcs(comm), 0),
finishedSendsCalled_(false)
{}
@ -90,6 +92,7 @@ void Foam::PstreamBuffers::finishedSends(const bool block)
recvBuf_,
sizes,
tag_,
comm_,
block
);
}
@ -108,6 +111,7 @@ void Foam::PstreamBuffers::finishedSends(labelListList& sizes, const bool block)
recvBuf_,
sizes,
tag_,
comm_,
block
);
}
@ -123,9 +127,9 @@ void Foam::PstreamBuffers::finishedSends(labelListList& sizes, const bool block)
// Note: possible only if using different tag from write started
// by ~UOPstream. Needs some work.
//sizes.setSize(UPstream::nProcs());
//labelList& nsTransPs = sizes[UPstream::myProcNo()];
//nsTransPs.setSize(UPstream::nProcs());
//sizes.setSize(UPstream::nProcs(comm));
//labelList& nsTransPs = sizes[UPstream::myProcNo(comm)];
//nsTransPs.setSize(UPstream::nProcs(comm));
//
//forAll(sendBuf_, procI)
//{

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
@ -95,6 +95,8 @@ class PstreamBuffers
const int tag_;
const label comm_;
const IOstream::streamFormat format_;
const IOstream::versionNumber version_;
@ -127,6 +129,7 @@ public:
(
const UPstream::commsTypes commsType,
const int tag = UPstream::msgType(),
const label comm = UPstream::worldComm,
IOstream::streamFormat format=IOstream::BINARY,
IOstream::versionNumber version=IOstream::currentVersion
);

View File

@ -53,11 +53,12 @@ void combineReduce
const List<UPstream::commsStruct>& comms,
T& Value,
const CombineOp& cop,
const int tag
const int tag,
const label comm
)
{
Pstream::combineGather(comms, Value, cop, tag);
Pstream::combineScatter(comms, Value, tag);
Pstream::combineGather(comms, Value, cop, tag, comm);
Pstream::combineScatter(comms, Value, tag, comm);
}
@ -66,24 +67,45 @@ void combineReduce
(
T& Value,
const CombineOp& cop,
const int tag = Pstream::msgType()
const int tag = Pstream::msgType(),
const label comm = Pstream::worldComm
)
{
if (UPstream::nProcs() < UPstream::nProcsSimpleSum)
if (UPstream::nProcs(comm) < UPstream::nProcsSimpleSum)
{
Pstream::combineGather
(
UPstream::linearCommunication(),
UPstream::linearCommunication(comm),
Value,
cop,
tag
tag,
comm
);
Pstream::combineScatter
(
UPstream::linearCommunication(comm),
Value,
tag,
comm
);
Pstream::combineScatter(UPstream::linearCommunication(), Value, tag);
}
else
{
Pstream::combineGather(UPstream::treeCommunication(), Value, cop, tag);
Pstream::combineScatter(UPstream::treeCommunication(), Value, tag);
Pstream::combineGather
(
UPstream::treeCommunication(comm),
Value,
cop,
tag,
comm
);
Pstream::combineScatter
(
UPstream::treeCommunication(comm),
Value,
tag,
comm
);
}
}

View File

@ -44,11 +44,18 @@ void reduce
const List<UPstream::commsStruct>& comms,
T& Value,
const BinaryOp& bop,
const int tag
const int tag,
const label comm
)
{
Pstream::gather(comms, Value, bop, tag);
Pstream::scatter(comms, Value, tag);
if (UPstream::warnComm != -1 && comm != UPstream::warnComm)
{
Pout<< "** reducing:" << Value << " with comm:" << comm
<< endl;
error::printStack(Pout);
}
Pstream::gather(comms, Value, bop, tag, comm);
Pstream::scatter(comms, Value, tag, comm);
}
@ -58,16 +65,17 @@ void reduce
(
T& Value,
const BinaryOp& bop,
const int tag = Pstream::msgType()
const int tag = Pstream::msgType(),
const label comm = UPstream::worldComm
)
{
if (UPstream::nProcs() < UPstream::nProcsSimpleSum)
if (UPstream::nProcs(comm) < UPstream::nProcsSimpleSum)
{
reduce(UPstream::linearCommunication(), Value, bop, tag);
reduce(UPstream::linearCommunication(comm), Value, bop, tag, comm);
}
else
{
reduce(UPstream::treeCommunication(), Value, bop, tag);
reduce(UPstream::treeCommunication(comm), Value, bop, tag, comm);
}
}
@ -78,18 +86,33 @@ T returnReduce
(
const T& Value,
const BinaryOp& bop,
const int tag = Pstream::msgType()
const int tag = Pstream::msgType(),
const label comm = UPstream::worldComm
)
{
T WorkValue(Value);
if (UPstream::nProcs() < UPstream::nProcsSimpleSum)
if (UPstream::nProcs(comm) < UPstream::nProcsSimpleSum)
{
reduce(UPstream::linearCommunication(), WorkValue, bop, tag);
reduce
(
UPstream::linearCommunication(comm),
WorkValue,
bop,
tag,
comm
);
}
else
{
reduce(UPstream::treeCommunication(), WorkValue, bop, tag);
reduce
(
UPstream::treeCommunication(comm),
WorkValue,
bop,
tag,
comm
);
}
return WorkValue;
@ -102,11 +125,12 @@ void sumReduce
(
T& Value,
label& Count,
const int tag = Pstream::msgType()
const int tag = Pstream::msgType(),
const label comm = UPstream::worldComm
)
{
reduce(Value, sumOp<T>(), tag);
reduce(Count, sumOp<label>(), tag);
reduce(Value, sumOp<T>(), tag, comm);
reduce(Count, sumOp<label>(), tag, comm);
}
@ -117,10 +141,14 @@ void reduce
T& Value,
const BinaryOp& bop,
const int tag,
const label comm,
label& request
)
{
notImplemented("reduce(T&, const BinaryOp&, const int, label&");
notImplemented
(
"reduce(T&, const BinaryOp&, const int, const label, label&"
);
}
@ -129,28 +157,32 @@ void reduce
(
scalar& Value,
const sumOp<scalar>& bop,
const int tag = Pstream::msgType()
const int tag = Pstream::msgType(),
const label comm = UPstream::worldComm
);
void reduce
(
scalar& Value,
const minOp<scalar>& bop,
const int tag = Pstream::msgType()
const int tag = Pstream::msgType(),
const label comm = UPstream::worldComm
);
void reduce
(
vector2D& Value,
const sumOp<vector2D>& bop,
const int tag = Pstream::msgType()
const int tag = Pstream::msgType(),
const label comm = UPstream::worldComm
);
void sumReduce
(
scalar& Value,
label& Count,
const int tag = Pstream::msgType()
const int tag = Pstream::msgType(),
const label comm = UPstream::worldComm
);
void reduce
@ -158,6 +190,7 @@ void reduce
scalar& Value,
const sumOp<scalar>& bop,
const int tag,
const label comm,
label& request
);

View File

@ -337,7 +337,9 @@ Foam::Istream& Foam::UIPstream::rewind()
void Foam::UIPstream::print(Ostream& os) const
{
os << "Reading from processor " << fromProcNo_
<< " to processor " << myProcNo() << Foam::endl;
<< " using communicator " << comm_
<< " and tag " << tag_
<< Foam::endl;
}

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
@ -66,6 +66,8 @@ class UIPstream
const int tag_;
const label comm_;
const bool clearAtEnd_;
int messageSize_;
@ -97,6 +99,7 @@ public:
DynamicList<char>& externalBuf,
label& externalBufPosition,
const int tag = UPstream::msgType(),
const label comm = UPstream::worldComm,
const bool clearAtEnd = false, // destroy externalBuf if at end
streamFormat format=BINARY,
versionNumber version=currentVersion
@ -131,7 +134,8 @@ public:
const int fromProcNo,
char* buf,
const std::streamsize bufSize,
const int tag = UPstream::msgType()
const int tag = UPstream::msgType(),
const label communicator = 0
);
//- Return next token from stream

View File

@ -91,6 +91,7 @@ Foam::UOPstream::UOPstream
const int toProcNo,
DynamicList<char>& sendBuf,
const int tag,
const label comm,
const bool sendAtDestruct,
streamFormat format,
versionNumber version
@ -101,6 +102,7 @@ Foam::UOPstream::UOPstream
toProcNo_(toProcNo),
sendBuf_(sendBuf),
tag_(tag),
comm_(comm),
sendAtDestruct_(sendAtDestruct)
{
setOpened();
@ -115,6 +117,7 @@ Foam::UOPstream::UOPstream(const int toProcNo, PstreamBuffers& buffers)
toProcNo_(toProcNo),
sendBuf_(buffers.sendBuf_[toProcNo]),
tag_(buffers.tag_),
comm_(buffers.comm_),
sendAtDestruct_(buffers.commsType_ != UPstream::nonBlocking)
{
setOpened();
@ -136,7 +139,8 @@ Foam::UOPstream::~UOPstream()
toProcNo_,
sendBuf_.begin(),
sendBuf_.size(),
tag_
tag_,
comm_
)
)
{
@ -292,7 +296,8 @@ Foam::Ostream& Foam::UOPstream::write(const char* data, std::streamsize count)
void Foam::UOPstream::print(Ostream& os) const
{
os << "Writing from processor " << toProcNo_
<< " to processor " << myProcNo() << Foam::endl;
<< " to processor " << myProcNo() << " in communicator " << comm_
<< " and tag " << tag_ << Foam::endl;
}

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
@ -65,6 +65,8 @@ class UOPstream
const int tag_;
const label comm_;
const bool sendAtDestruct_;
@ -93,6 +95,7 @@ public:
const int toProcNo,
DynamicList<char>& sendBuf,
const int tag = UPstream::msgType(),
const label comm = UPstream::worldComm,
const bool sendAtDestruct = true,
streamFormat format=BINARY,
versionNumber version=currentVersion
@ -126,7 +129,8 @@ public:
const int toProcNo,
const char* buf,
const std::streamsize bufSize,
const int tag = UPstream::msgType()
const int tag = UPstream::msgType(),
const label communicator = 0
);
//- Write next token to stream

View File

@ -54,18 +54,33 @@ const Foam::NamedEnum<Foam::UPstream::commsTypes, 3>
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::UPstream::setParRun()
void Foam::UPstream::setParRun(const label nProcs)
{
parRun_ = true;
Pout.prefix() = '[' + name(myProcNo()) + "] ";
Perr.prefix() = '[' + name(myProcNo()) + "] ";
// Redo worldComm communicator (this has been created at static
// initialisation time)
freeCommunicator(UPstream::worldComm);
label comm = allocateCommunicator(-1, identity(nProcs), true);
if (comm != UPstream::worldComm)
{
FatalErrorIn("UPstream::setParRun(const label)")
<< "problem : comm:" << comm
<< " UPstream::worldComm:" << UPstream::worldComm
<< Foam::exit(FatalError);
}
Pout.prefix() = '[' + name(myProcNo(Pstream::worldComm)) + "] ";
Perr.prefix() = '[' + name(myProcNo(Pstream::worldComm)) + "] ";
}
void Foam::UPstream::calcLinearComm(const label nProcs)
Foam::List<Foam::UPstream::commsStruct> Foam::UPstream::calcLinearComm
(
const label nProcs
)
{
linearCommunication_.setSize(nProcs);
List<commsStruct> linearCommunication(nProcs);
// Master
labelList belowIDs(nProcs - 1);
@ -74,7 +89,7 @@ void Foam::UPstream::calcLinearComm(const label nProcs)
belowIDs[i] = i + 1;
}
linearCommunication_[0] = commsStruct
linearCommunication[0] = commsStruct
(
nProcs,
0,
@ -86,7 +101,7 @@ void Foam::UPstream::calcLinearComm(const label nProcs)
// Slaves. Have no below processors, only communicate up to master
for (label procID = 1; procID < nProcs; procID++)
{
linearCommunication_[procID] = commsStruct
linearCommunication[procID] = commsStruct
(
nProcs,
procID,
@ -95,6 +110,7 @@ void Foam::UPstream::calcLinearComm(const label nProcs)
labelList(0)
);
}
return linearCommunication;
}
@ -142,7 +158,10 @@ void Foam::UPstream::collectReceives
// 5 - 4
// 6 7 4
// 7 - 6
void Foam::UPstream::calcTreeComm(label nProcs)
Foam::List<Foam::UPstream::commsStruct> Foam::UPstream::calcTreeComm
(
label nProcs
)
{
label nLevels = 1;
while ((1 << nLevels) < nProcs)
@ -188,11 +207,11 @@ void Foam::UPstream::calcTreeComm(label nProcs)
}
treeCommunication_.setSize(nProcs);
List<commsStruct> treeCommunication(nProcs);
for (label procID = 0; procID < nProcs; procID++)
{
treeCommunication_[procID] = commsStruct
treeCommunication[procID] = commsStruct
(
nProcs,
procID,
@ -201,37 +220,198 @@ void Foam::UPstream::calcTreeComm(label nProcs)
allReceives[procID].shrink()
);
}
return treeCommunication;
}
// Callback from UPstream::init() : initialize linear and tree communication
// schedules now that nProcs is known.
void Foam::UPstream::initCommunicationSchedule()
Foam::label Foam::UPstream::allocateCommunicator
(
const label parentIndex,
const labelList& subRanks,
const bool doPstream
)
{
calcLinearComm(nProcs());
calcTreeComm(nProcs());
label index;
if (!freeComms_.empty())
{
index = freeComms_.pop();
}
else
{
// Extend storage
index = parentCommunicator_.size();
myProcNo_.append(-1);
procIDs_.append(List<int>(0));
parentCommunicator_.append(-1);
linearCommunication_.append(List<commsStruct>(0));
treeCommunication_.append(List<commsStruct>(0));
}
if (debug)
{
Pout<< "Communicators : Allocating communicator " << index << endl
<< " parent : " << parentIndex << endl
<< " procs : " << subRanks << endl
<< endl;
}
// Initialise; overwritten by allocatePstreamCommunicator
myProcNo_[index] = 0;
// Convert from label to int
procIDs_[index].setSize(subRanks.size());
forAll(procIDs_[index], i)
{
procIDs_[index][i] = subRanks[i];
// Enforce incremental order (so index is rank in next communicator)
if (i >= 1 && subRanks[i] <= subRanks[i-1])
{
FatalErrorIn
(
"UPstream::allocateCommunicator"
"(const label, const labelList&, const bool)"
) << "subranks not sorted : " << subRanks
<< " when allocating subcommunicator from parent "
<< parentIndex
<< Foam::abort(FatalError);
}
}
parentCommunicator_[index] = parentIndex;
linearCommunication_[index] = calcLinearComm(procIDs_[index].size());
treeCommunication_[index] = calcTreeComm(procIDs_[index].size());
if (doPstream && parRun())
{
allocatePstreamCommunicator(parentIndex, index);
}
return index;
}
void Foam::UPstream::freeCommunicator
(
const label communicator,
const bool doPstream
)
{
if (debug)
{
Pout<< "Communicators : Freeing communicator " << communicator << endl
<< " parent : " << parentCommunicator_[communicator] << endl
<< " myProcNo : " << myProcNo_[communicator] << endl
<< endl;
}
if (doPstream && parRun())
{
freePstreamCommunicator(communicator);
}
myProcNo_[communicator] = -1;
//procIDs_[communicator].clear();
parentCommunicator_[communicator] = -1;
linearCommunication_[communicator].clear();
treeCommunication_[communicator].clear();
freeComms_.push(communicator);
}
void Foam::UPstream::freeCommunicators(const bool doPstream)
{
forAll(myProcNo_, communicator)
{
if (myProcNo_[communicator] != -1)
{
freeCommunicator(communicator, doPstream);
}
}
}
int Foam::UPstream::baseProcNo(const label myComm, const int myProcID)
{
int procID = myProcID;
label comm = myComm;
while (parent(comm) != -1)
{
const List<int>& parentRanks = UPstream::procID(comm);
procID = parentRanks[procID];
comm = UPstream::parent(comm);
}
return procID;
}
Foam::label Foam::UPstream::procNo(const label myComm, const int baseProcID)
{
const List<int>& parentRanks = procID(myComm);
label parentComm = parent(myComm);
if (parentComm == -1)
{
return findIndex(parentRanks, baseProcID);
}
else
{
label parentRank = procNo(parentComm, baseProcID);
return findIndex(parentRanks, parentRank);
}
}
Foam::label Foam::UPstream::procNo
(
const label myComm,
const label currentComm,
const int currentProcID
)
{
label physProcID = UPstream::baseProcNo(currentComm, currentProcID);
return procNo(myComm, physProcID);
}
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
// Initialise my process number to 0 (the master)
int Foam::UPstream::myProcNo_(0);
// By default this is not a parallel run
bool Foam::UPstream::parRun_(false);
// Free communicators
Foam::LIFOStack<Foam::label> Foam::UPstream::freeComms_;
// My processor number
Foam::DynamicList<int> Foam::UPstream::myProcNo_(10);
// List of process IDs
Foam::List<int> Foam::UPstream::procIDs_(label(1), 0);
Foam::DynamicList<Foam::List<int> > Foam::UPstream::procIDs_(10);
// Parent communicator
Foam::DynamicList<Foam::label> Foam::UPstream::parentCommunicator_(10);
// Standard transfer message type
int Foam::UPstream::msgType_(1);
// Linear communication schedule
Foam::List<Foam::UPstream::commsStruct> Foam::UPstream::linearCommunication_(0);
Foam::DynamicList<Foam::List<Foam::UPstream::commsStruct> >
Foam::UPstream::linearCommunication_(10);
// Multi level communication schedule
Foam::List<Foam::UPstream::commsStruct> Foam::UPstream::treeCommunication_(0);
Foam::DynamicList<Foam::List<Foam::UPstream::commsStruct> >
Foam::UPstream::treeCommunication_(10);
// Allocate a serial communicator. This gets overwritten in parallel mode
// (by UPstream::setParRun())
Foam::UPstream::communicator serialComm(-1, Foam::labelList(1, 0), false);
// Should compact transfer be used in which floats replace doubles
// reducing the bandwidth requirement at the expense of some loss
@ -292,6 +472,14 @@ public:
addcommsTypeToOpt addcommsTypeToOpt_("commsType");
// Default communicator
Foam::label Foam::UPstream::worldComm(0);
// Warn for use of any communicator
Foam::label Foam::UPstream::warnComm(-1);
// Number of polling cycles in processor updates
int Foam::UPstream::nPollProcInterfaces
(

View File

@ -29,7 +29,6 @@ Description
SourceFiles
UPstream.C
UPstreamsPrint.C
UPstreamCommsStruct.C
gatherScatter.C
combineGatherScatter.C
@ -45,6 +44,8 @@ SourceFiles
#include "HashTable.H"
#include "string.H"
#include "NamedEnum.H"
#include "ListOps.H"
#include "LIFOStack.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -180,26 +181,30 @@ private:
// Private data
static int myProcNo_;
static bool parRun_;
static List<int> procIDs_;
static int msgType_;
static List<commsStruct> linearCommunication_;
static List<commsStruct> treeCommunication_;
// Communicator specific data
static LIFOStack<label> freeComms_;
static DynamicList<int> myProcNo_;
static DynamicList<List<int> > procIDs_;
static DynamicList<label> parentCommunicator_;
static DynamicList<List<commsStruct> > linearCommunication_;
static DynamicList<List<commsStruct> > treeCommunication_;
// Private Member Functions
//- Set data for parallel running
static void setParRun();
static void setParRun(const label nProcs);
//- Calculate linear communication schedule
static void calcLinearComm(const label nProcs);
static List<commsStruct> calcLinearComm(const label nProcs);
//- Calculate tree communication schedule
static void calcTreeComm(const label nProcs);
static List<commsStruct> calcTreeComm(const label nProcs);
//- Helper function for tree communication schedule determination
// Collects all processorIDs below a processor
@ -210,10 +215,18 @@ private:
DynamicList<label>& allReceives
);
//- Initialize all communication schedules. Callback from
// UPstream::init()
static void initCommunicationSchedule();
//- Allocate a communicator with index
static void allocatePstreamCommunicator
(
const label parentIndex,
const label index
);
//- Free a communicator
static void freePstreamCommunicator
(
const label index
);
protected:
@ -245,6 +258,13 @@ public:
//- Number of polling cycles in processor updates
static int nPollProcInterfaces;
//- Default communicator (all processors)
static label worldComm;
//- Debugging: warn for use of any communicator differing from warnComm
static label warnComm;
// Constructors
//- Construct given optional buffer size
@ -256,6 +276,73 @@ public:
// Member functions
//- Allocate a new communicator
static label allocateCommunicator
(
const label parent,
const labelList& subRanks,
const bool doPstream = true
);
//- Free a previously allocated communicator
static void freeCommunicator
(
const label communicator,
const bool doPstream = true
);
//- Free all communicators
static void freeCommunicators(const bool doPstream);
//- Helper class for allocating/freeing communicators
class communicator
{
label comm_;
//- Disallow copy and assignment
communicator(const communicator&);
void operator=(const communicator&);
public:
communicator
(
const label parent,
const labelList& subRanks,
const bool doPstream
)
:
comm_(allocateCommunicator(parent, subRanks, doPstream))
{}
~communicator()
{
freeCommunicator(comm_);
}
operator label() const
{
return comm_;
}
};
//- Return physical processor number (i.e. processor number in
// worldComm) given communicator and procssor
static int baseProcNo(const label myComm, const int procID);
//- Return processor number in communicator (given physical processor
// number) (= reverse of baseProcNo)
static label procNo(const label comm, const int baseProcID);
//- Return processor number in communicator (given processor number
// and communicator)
static label procNo
(
const label myComm,
const label currentComm,
const int currentProcID
);
//- Add the valid option this type of communications library
// adds/requires on the command line
static void addValidParOptions(HashTable<string>& validParOptions);
@ -281,6 +368,14 @@ public:
//- Non-blocking comms: has request i finished?
static bool finishedRequest(const label i);
static int allocateTag(const char*);
static int allocateTag(const word&);
static void freeTag(const char*, const int tag);
static void freeTag(const word&, const int tag);
//- Is this a parallel run?
static bool& parRun()
@ -289,15 +384,9 @@ public:
}
//- Number of processes in parallel run
static label nProcs()
static label nProcs(const label communicator = 0)
{
return procIDs_.size();
}
//- Am I the master process
static bool master()
{
return myProcNo_ == 0;
return procIDs_[communicator].size();
}
//- Process index of the master
@ -306,22 +395,27 @@ public:
return 0;
}
//- Number of this process (starting from masterNo() = 0)
static int myProcNo()
//- Am I the master process
static bool master(const label communicator = 0)
{
return myProcNo_;
return myProcNo_[communicator] == masterNo();
}
//- Process IDs
static const List<int>& procIDs()
//- Number of this process (starting from masterNo() = 0)
static int myProcNo(const label communicator = 0)
{
return procIDs_;
return myProcNo_[communicator];
}
static label parent(const label communicator)
{
return parentCommunicator_(communicator);
}
//- Process ID of given process index
static int procID(int procNo)
static List<int>& procID(int communicator)
{
return procIDs_[procNo];
return procIDs_[communicator];
}
//- Process index of first slave
@ -331,21 +425,27 @@ public:
}
//- Process index of last slave
static int lastSlave()
static int lastSlave(const label communicator = 0)
{
return nProcs() - 1;
return nProcs(communicator) - 1;
}
//- Communication schedule for linear all-to-master (proc 0)
static const List<commsStruct>& linearCommunication()
static const List<commsStruct>& linearCommunication
(
const label communicator = 0
)
{
return linearCommunication_;
return linearCommunication_[communicator];
}
//- Communication schedule for tree all-to-master (proc 0)
static const List<commsStruct>& treeCommunication()
static const List<commsStruct>& treeCommunication
(
const label communicator = 0
)
{
return treeCommunication_;
return treeCommunication_[communicator];
}
//- Message tag of standard messages

View File

@ -51,13 +51,14 @@ void Pstream::combineGather
const List<UPstream::commsStruct>& comms,
T& Value,
const CombineOp& cop,
const int tag
const int tag,
const label comm
)
{
if (UPstream::parRun())
{
// Get my communication order
const commsStruct& myComm = comms[UPstream::myProcNo()];
const commsStruct& myComm = comms[UPstream::myProcNo(comm)];
// Receive from my downstairs neighbours
forAll(myComm.below(), belowI)
@ -73,7 +74,8 @@ void Pstream::combineGather
belowID,
reinterpret_cast<char*>(&value),
sizeof(T),
tag
tag,
comm
);
if (debug & 2)
@ -86,7 +88,7 @@ void Pstream::combineGather
}
else
{
IPstream fromBelow(UPstream::scheduled, belowID, 0, tag);
IPstream fromBelow(UPstream::scheduled, belowID, 0, tag, comm);
T value(fromBelow);
if (debug & 2)
@ -116,12 +118,20 @@ void Pstream::combineGather
myComm.above(),
reinterpret_cast<const char*>(&Value),
sizeof(T),
tag
tag,
comm
);
}
else
{
OPstream toAbove(UPstream::scheduled, myComm.above(), 0, tag);
OPstream toAbove
(
UPstream::scheduled,
myComm.above(),
0,
tag,
comm
);
toAbove << Value;
}
}
@ -129,16 +139,36 @@ void Pstream::combineGather
}
template<class T, class CombineOp>
void Pstream::combineGather(T& Value, const CombineOp& cop, const int tag)
template <class T, class CombineOp>
void Pstream::combineGather
(
T& Value,
const CombineOp& cop,
const int tag,
const label comm
)
{
if (UPstream::nProcs() < UPstream::nProcsSimpleSum)
if (UPstream::nProcs(comm) < UPstream::nProcsSimpleSum)
{
combineGather(UPstream::linearCommunication(), Value, cop, tag);
combineGather
(
UPstream::linearCommunication(comm),
Value,
cop,
tag,
comm
);
}
else
{
combineGather(UPstream::treeCommunication(), Value, cop, tag);
combineGather
(
UPstream::treeCommunication(comm),
Value,
cop,
tag,
comm
);
}
}
@ -148,13 +178,14 @@ void Pstream::combineScatter
(
const List<UPstream::commsStruct>& comms,
T& Value,
const int tag
const int tag,
const label comm
)
{
if (UPstream::parRun())
{
// Get my communication order
const UPstream::commsStruct& myComm = comms[UPstream::myProcNo()];
const UPstream::commsStruct& myComm = comms[UPstream::myProcNo(comm)];
// Reveive from up
if (myComm.above() != -1)
@ -167,12 +198,20 @@ void Pstream::combineScatter
myComm.above(),
reinterpret_cast<char*>(&Value),
sizeof(T),
tag
tag,
comm
);
}
else
{
IPstream fromAbove(UPstream::scheduled, myComm.above(), 0, tag);
IPstream fromAbove
(
UPstream::scheduled,
myComm.above(),
0,
tag,
comm
);
Value = T(fromAbove);
}
@ -201,12 +240,13 @@ void Pstream::combineScatter
belowID,
reinterpret_cast<const char*>(&Value),
sizeof(T),
tag
tag,
comm
);
}
else
{
OPstream toBelow(UPstream::scheduled, belowID, 0, tag);
OPstream toBelow(UPstream::scheduled, belowID, 0, tag, comm);
toBelow << Value;
}
}
@ -214,16 +254,21 @@ void Pstream::combineScatter
}
template<class T>
void Pstream::combineScatter(T& Value, const int tag)
template <class T>
void Pstream::combineScatter
(
T& Value,
const int tag,
const label comm
)
{
if (UPstream::nProcs() < UPstream::nProcsSimpleSum)
if (UPstream::nProcs(comm) < UPstream::nProcsSimpleSum)
{
combineScatter(UPstream::linearCommunication(), Value, tag);
combineScatter(UPstream::linearCommunication(comm), Value, tag, comm);
}
else
{
combineScatter(UPstream::treeCommunication(), Value, tag);
combineScatter(UPstream::treeCommunication(comm), Value, tag, comm);
}
}
@ -238,13 +283,14 @@ void Pstream::listCombineGather
const List<UPstream::commsStruct>& comms,
List<T>& Values,
const CombineOp& cop,
const int tag
const int tag,
const label comm
)
{
if (UPstream::parRun())
{
// Get my communication order
const commsStruct& myComm = comms[UPstream::myProcNo()];
const commsStruct& myComm = comms[UPstream::myProcNo(comm)];
// Receive from my downstairs neighbours
forAll(myComm.below(), belowI)
@ -261,7 +307,8 @@ void Pstream::listCombineGather
belowID,
reinterpret_cast<char*>(receivedValues.begin()),
receivedValues.byteSize(),
tag
tag,
comm
);
if (debug & 2)
@ -277,7 +324,7 @@ void Pstream::listCombineGather
}
else
{
IPstream fromBelow(UPstream::scheduled, belowID, 0, tag);
IPstream fromBelow(UPstream::scheduled, belowID, 0, tag, comm);
List<T> receivedValues(fromBelow);
if (debug & 2)
@ -310,12 +357,20 @@ void Pstream::listCombineGather
myComm.above(),
reinterpret_cast<const char*>(Values.begin()),
Values.byteSize(),
tag
tag,
comm
);
}
else
{
OPstream toAbove(UPstream::scheduled, myComm.above(), 0, tag);
OPstream toAbove
(
UPstream::scheduled,
myComm.above(),
0,
tag,
comm
);
toAbove << Values;
}
}
@ -328,16 +383,31 @@ void Pstream::listCombineGather
(
List<T>& Values,
const CombineOp& cop,
const int tag
const int tag,
const label comm
)
{
if (UPstream::nProcs() < UPstream::nProcsSimpleSum)
if (UPstream::nProcs(comm) < UPstream::nProcsSimpleSum)
{
listCombineGather(UPstream::linearCommunication(), Values, cop, tag);
listCombineGather
(
UPstream::linearCommunication(comm),
Values,
cop,
tag,
comm
);
}
else
{
listCombineGather(UPstream::treeCommunication(), Values, cop, tag);
listCombineGather
(
UPstream::treeCommunication(comm),
Values,
cop,
tag,
comm
);
}
}
@ -347,13 +417,14 @@ void Pstream::listCombineScatter
(
const List<UPstream::commsStruct>& comms,
List<T>& Values,
const int tag
const int tag,
const label comm
)
{
if (UPstream::parRun())
{
// Get my communication order
const UPstream::commsStruct& myComm = comms[UPstream::myProcNo()];
const UPstream::commsStruct& myComm = comms[UPstream::myProcNo(comm)];
// Reveive from up
if (myComm.above() != -1)
@ -366,12 +437,20 @@ void Pstream::listCombineScatter
myComm.above(),
reinterpret_cast<char*>(Values.begin()),
Values.byteSize(),
tag
tag,
comm
);
}
else
{
IPstream fromAbove(UPstream::scheduled, myComm.above(), 0, tag);
IPstream fromAbove
(
UPstream::scheduled,
myComm.above(),
0,
tag,
comm
);
fromAbove >> Values;
}
@ -400,12 +479,13 @@ void Pstream::listCombineScatter
belowID,
reinterpret_cast<const char*>(Values.begin()),
Values.byteSize(),
tag
tag,
comm
);
}
else
{
OPstream toBelow(UPstream::scheduled, belowID, 0, tag);
OPstream toBelow(UPstream::scheduled, belowID, 0, tag, comm);
toBelow << Values;
}
}
@ -413,16 +493,33 @@ void Pstream::listCombineScatter
}
template<class T>
void Pstream::listCombineScatter(List<T>& Values, const int tag)
template <class T>
void Pstream::listCombineScatter
(
List<T>& Values,
const int tag,
const label comm
)
{
if (UPstream::nProcs() < UPstream::nProcsSimpleSum)
if (UPstream::nProcs(comm) < UPstream::nProcsSimpleSum)
{
listCombineScatter(UPstream::linearCommunication(), Values, tag);
listCombineScatter
(
UPstream::linearCommunication(comm),
Values,
tag,
comm
);
}
else
{
listCombineScatter(UPstream::treeCommunication(), Values, tag);
listCombineScatter
(
UPstream::treeCommunication(comm),
Values,
tag,
comm
);
}
}
@ -439,20 +536,21 @@ void Pstream::mapCombineGather
const List<UPstream::commsStruct>& comms,
Container& Values,
const CombineOp& cop,
const int tag
const int tag,
const label comm
)
{
if (UPstream::parRun())
{
// Get my communication order
const commsStruct& myComm = comms[UPstream::myProcNo()];
const commsStruct& myComm = comms[UPstream::myProcNo(comm)];
// Receive from my downstairs neighbours
forAll(myComm.below(), belowI)
{
label belowID = myComm.below()[belowI];
IPstream fromBelow(UPstream::scheduled, belowID, 0, tag);
IPstream fromBelow(UPstream::scheduled, belowID, 0, tag, comm);
Container receivedValues(fromBelow);
if (debug & 2)
@ -492,7 +590,7 @@ void Pstream::mapCombineGather
<< " data:" << Values << endl;
}
OPstream toAbove(UPstream::scheduled, myComm.above(), 0, tag);
OPstream toAbove(UPstream::scheduled, myComm.above(), 0, tag, comm);
toAbove << Values;
}
}
@ -504,16 +602,31 @@ void Pstream::mapCombineGather
(
Container& Values,
const CombineOp& cop,
const int tag
const int tag,
const label comm
)
{
if (UPstream::nProcs() < UPstream::nProcsSimpleSum)
if (UPstream::nProcs(comm) < UPstream::nProcsSimpleSum)
{
mapCombineGather(UPstream::linearCommunication(), Values, cop, tag);
mapCombineGather
(
UPstream::linearCommunication(comm),
Values,
cop,
tag,
comm
);
}
else
{
mapCombineGather(UPstream::treeCommunication(), Values, cop, tag);
mapCombineGather
(
UPstream::treeCommunication(comm),
Values,
cop,
tag,
comm
);
}
}
@ -523,18 +636,26 @@ void Pstream::mapCombineScatter
(
const List<UPstream::commsStruct>& comms,
Container& Values,
const int tag
const int tag,
const label comm
)
{
if (UPstream::parRun())
{
// Get my communication order
const UPstream::commsStruct& myComm = comms[UPstream::myProcNo()];
const UPstream::commsStruct& myComm = comms[UPstream::myProcNo(comm)];
// Reveive from up
if (myComm.above() != -1)
{
IPstream fromAbove(UPstream::scheduled, myComm.above(), 0, tag);
IPstream fromAbove
(
UPstream::scheduled,
myComm.above(),
0,
tag,
comm
);
fromAbove >> Values;
if (debug & 2)
@ -554,23 +675,40 @@ void Pstream::mapCombineScatter
Pout<< " sending to " << belowID << " data:" << Values << endl;
}
OPstream toBelow(UPstream::scheduled, belowID, 0, tag);
OPstream toBelow(UPstream::scheduled, belowID, 0, tag, comm);
toBelow << Values;
}
}
}
template<class Container>
void Pstream::mapCombineScatter(Container& Values, const int tag)
template <class Container>
void Pstream::mapCombineScatter
(
Container& Values,
const int tag,
const label comm
)
{
if (UPstream::nProcs() < UPstream::nProcsSimpleSum)
if (UPstream::nProcs(comm) < UPstream::nProcsSimpleSum)
{
mapCombineScatter(UPstream::linearCommunication(), Values, tag);
mapCombineScatter
(
UPstream::linearCommunication(comm),
Values,
tag,
comm
);
}
else
{
mapCombineScatter(UPstream::treeCommunication(), Values, tag);
mapCombineScatter
(
UPstream::treeCommunication(comm),
Values,
tag,
comm
);
}
}

View File

@ -46,6 +46,7 @@ void Pstream::exchange
List<Container>& recvBufs,
labelListList& sizes,
const int tag,
const label comm,
const bool block
)
{
@ -57,20 +58,20 @@ void Pstream::exchange
) << "Continuous data only." << Foam::abort(FatalError);
}
if (sendBufs.size() != UPstream::nProcs())
if (sendBufs.size() != UPstream::nProcs(comm))
{
FatalErrorIn
(
"Pstream::exchange(..)"
) << "Size of list:" << sendBufs.size()
<< " does not equal the number of processors:"
<< UPstream::nProcs()
<< UPstream::nProcs(comm)
<< Foam::abort(FatalError);
}
sizes.setSize(UPstream::nProcs());
labelList& nsTransPs = sizes[UPstream::myProcNo()];
nsTransPs.setSize(UPstream::nProcs());
sizes.setSize(UPstream::nProcs(comm));
labelList& nsTransPs = sizes[UPstream::myProcNo(comm)];
nsTransPs.setSize(UPstream::nProcs(comm));
forAll(sendBufs, procI)
{
@ -78,7 +79,7 @@ void Pstream::exchange
}
// Send sizes across. Note: blocks.
combineReduce(sizes, UPstream::listEq(), tag);
combineReduce(sizes, UPstream::listEq(), tag, comm);
if (Pstream::parRun())
{
@ -90,9 +91,9 @@ void Pstream::exchange
recvBufs.setSize(sendBufs.size());
forAll(sizes, procI)
{
label nRecv = sizes[procI][UPstream::myProcNo()];
label nRecv = sizes[procI][UPstream::myProcNo(comm)];
if (procI != Pstream::myProcNo() && nRecv > 0)
if (procI != Pstream::myProcNo(comm) && nRecv > 0)
{
recvBufs[procI].setSize(nRecv);
UIPstream::read
@ -101,7 +102,8 @@ void Pstream::exchange
procI,
reinterpret_cast<char*>(recvBufs[procI].begin()),
nRecv*sizeof(T),
tag
tag,
comm
);
}
}
@ -112,7 +114,7 @@ void Pstream::exchange
forAll(sendBufs, procI)
{
if (procI != Pstream::myProcNo() && sendBufs[procI].size() > 0)
if (procI != Pstream::myProcNo(comm) && sendBufs[procI].size() > 0)
{
if
(
@ -122,7 +124,8 @@ void Pstream::exchange
procI,
reinterpret_cast<const char*>(sendBufs[procI].begin()),
sendBufs[procI].size()*sizeof(T),
tag
tag,
comm
)
)
{
@ -146,7 +149,7 @@ void Pstream::exchange
}
// Do myself
recvBufs[Pstream::myProcNo()] = sendBufs[Pstream::myProcNo()];
recvBufs[Pstream::myProcNo(comm)] = sendBufs[Pstream::myProcNo(comm)];
}

View File

@ -48,13 +48,14 @@ void Pstream::gather
const List<UPstream::commsStruct>& comms,
T& Value,
const BinaryOp& bop,
const int tag
const int tag,
const label comm
)
{
if (UPstream::parRun())
{
// Get my communication order
const commsStruct& myComm = comms[UPstream::myProcNo()];
const commsStruct& myComm = comms[UPstream::myProcNo(comm)];
// Receive from my downstairs neighbours
forAll(myComm.below(), belowI)
@ -69,7 +70,8 @@ void Pstream::gather
myComm.below()[belowI],
reinterpret_cast<char*>(&value),
sizeof(T),
tag
tag,
comm
);
}
else
@ -79,7 +81,8 @@ void Pstream::gather
UPstream::scheduled,
myComm.below()[belowI],
0,
tag
tag,
comm
);
fromBelow >> value;
}
@ -98,12 +101,20 @@ void Pstream::gather
myComm.above(),
reinterpret_cast<const char*>(&Value),
sizeof(T),
tag
tag,
comm
);
}
else
{
OPstream toAbove(UPstream::scheduled, myComm.above(), 0, tag);
OPstream toAbove
(
UPstream::scheduled,
myComm.above(),
0,
tag,
comm
);
toAbove << Value;
}
}
@ -111,16 +122,22 @@ void Pstream::gather
}
template<class T, class BinaryOp>
void Pstream::gather(T& Value, const BinaryOp& bop, const int tag)
template <class T, class BinaryOp>
void Pstream::gather
(
T& Value,
const BinaryOp& bop,
const int tag,
const label comm
)
{
if (UPstream::nProcs() < UPstream::nProcsSimpleSum)
if (UPstream::nProcs(comm) < UPstream::nProcsSimpleSum)
{
gather(UPstream::linearCommunication(), Value, bop, tag);
gather(UPstream::linearCommunication(comm), Value, bop, tag, comm);
}
else
{
gather(UPstream::treeCommunication(), Value, bop, tag);
gather(UPstream::treeCommunication(comm), Value, bop, tag, comm);
}
}
@ -130,13 +147,14 @@ void Pstream::scatter
(
const List<UPstream::commsStruct>& comms,
T& Value,
const int tag
const int tag,
const label comm
)
{
if (UPstream::parRun())
{
// Get my communication order
const commsStruct& myComm = comms[UPstream::myProcNo()];
const commsStruct& myComm = comms[UPstream::myProcNo(comm)];
// Reveive from up
if (myComm.above() != -1)
@ -149,12 +167,20 @@ void Pstream::scatter
myComm.above(),
reinterpret_cast<char*>(&Value),
sizeof(T),
tag
tag,
comm
);
}
else
{
IPstream fromAbove(UPstream::scheduled, myComm.above(), 0, tag);
IPstream fromAbove
(
UPstream::scheduled,
myComm.above(),
0,
tag,
comm
);
fromAbove >> Value;
}
}
@ -170,7 +196,8 @@ void Pstream::scatter
myComm.below()[belowI],
reinterpret_cast<const char*>(&Value),
sizeof(T),
tag
tag,
comm
);
}
else
@ -180,7 +207,8 @@ void Pstream::scatter
UPstream::scheduled,
myComm.below()[belowI],
0,
tag
tag,
comm
);
toBelow << Value;
}
@ -189,16 +217,16 @@ void Pstream::scatter
}
template<class T>
void Pstream::scatter(T& Value, const int tag)
template <class T>
void Pstream::scatter(T& Value, const int tag, const label comm)
{
if (UPstream::nProcs() < UPstream::nProcsSimpleSum)
if (UPstream::nProcs(comm) < UPstream::nProcsSimpleSum)
{
scatter(UPstream::linearCommunication(), Value, tag);
scatter(UPstream::linearCommunication(comm), Value, tag, comm);
}
else
{
scatter(UPstream::treeCommunication(), Value, tag);
scatter(UPstream::treeCommunication(comm), Value, tag, comm);
}
}

View File

@ -26,7 +26,7 @@ Description
communication schedule (usually linear-to-master or tree-to-master).
The gathered data will be a list with element procID the data from processor
procID. Before calling every processor should insert its value into
Values[UPstream::myProcNo()].
Values[UPstream::myProcNo(comm)].
Note: after gather every processor only knows its own data and that of the
processors below it. Only the 'master' of the communication schedule holds
a fully filled List. Use scatter to distribute the data.
@ -49,12 +49,13 @@ void Pstream::gatherList
(
const List<UPstream::commsStruct>& comms,
List<T>& Values,
const int tag
const int tag,
const label comm
)
{
if (UPstream::parRun())
{
if (Values.size() != UPstream::nProcs())
if (Values.size() != UPstream::nProcs(comm))
{
FatalErrorIn
(
@ -62,12 +63,12 @@ void Pstream::gatherList
", List<T>)"
) << "Size of list:" << Values.size()
<< " does not equal the number of processors:"
<< UPstream::nProcs()
<< UPstream::nProcs(comm)
<< Foam::abort(FatalError);
}
// Get my communication order
const commsStruct& myComm = comms[UPstream::myProcNo()];
const commsStruct& myComm = comms[UPstream::myProcNo(comm)];
// Receive from my downstairs neighbours
forAll(myComm.below(), belowI)
@ -85,7 +86,8 @@ void Pstream::gatherList
belowID,
reinterpret_cast<char*>(receivedValues.begin()),
receivedValues.byteSize(),
tag
tag,
comm
);
Values[belowID] = receivedValues[0];
@ -97,7 +99,7 @@ void Pstream::gatherList
}
else
{
IPstream fromBelow(UPstream::scheduled, belowID, 0, tag);
IPstream fromBelow(UPstream::scheduled, belowID, 0, tag, comm);
fromBelow >> Values[belowID];
if (debug & 2)
@ -133,14 +135,14 @@ void Pstream::gatherList
if (debug & 2)
{
Pout<< " sending to " << myComm.above()
<< " data from me:" << UPstream::myProcNo()
<< " data:" << Values[UPstream::myProcNo()] << endl;
<< " data from me:" << UPstream::myProcNo(comm)
<< " data:" << Values[UPstream::myProcNo(comm)] << endl;
}
if (contiguous<T>())
{
List<T> sendingValues(belowLeaves.size() + 1);
sendingValues[0] = Values[UPstream::myProcNo()];
sendingValues[0] = Values[UPstream::myProcNo(comm)];
forAll(belowLeaves, leafI)
{
@ -153,13 +155,21 @@ void Pstream::gatherList
myComm.above(),
reinterpret_cast<const char*>(sendingValues.begin()),
sendingValues.byteSize(),
tag
tag,
comm
);
}
else
{
OPstream toAbove(UPstream::scheduled, myComm.above(), 0, tag);
toAbove << Values[UPstream::myProcNo()];
OPstream toAbove
(
UPstream::scheduled,
myComm.above(),
0,
tag,
comm
);
toAbove << Values[UPstream::myProcNo(comm)];
forAll(belowLeaves, leafI)
{
@ -179,16 +189,16 @@ void Pstream::gatherList
}
template<class T>
void Pstream::gatherList(List<T>& Values, const int tag)
template <class T>
void Pstream::gatherList(List<T>& Values, const int tag, const label comm)
{
if (UPstream::nProcs() < UPstream::nProcsSimpleSum)
if (UPstream::nProcs(comm) < UPstream::nProcsSimpleSum)
{
gatherList(UPstream::linearCommunication(), Values, tag);
gatherList(UPstream::linearCommunication(comm), Values, tag, comm);
}
else
{
gatherList(UPstream::treeCommunication(), Values, tag);
gatherList(UPstream::treeCommunication(comm), Values, tag, comm);
}
}
@ -198,12 +208,13 @@ void Pstream::scatterList
(
const List<UPstream::commsStruct>& comms,
List<T>& Values,
const int tag
const int tag,
const label comm
)
{
if (UPstream::parRun())
{
if (Values.size() != UPstream::nProcs())
if (Values.size() != UPstream::nProcs(comm))
{
FatalErrorIn
(
@ -211,12 +222,12 @@ void Pstream::scatterList
", List<T>)"
) << "Size of list:" << Values.size()
<< " does not equal the number of processors:"
<< UPstream::nProcs()
<< UPstream::nProcs(comm)
<< Foam::abort(FatalError);
}
// Get my communication order
const commsStruct& myComm = comms[UPstream::myProcNo()];
const commsStruct& myComm = comms[UPstream::myProcNo(comm)];
// Reveive from up
if (myComm.above() != -1)
@ -233,7 +244,8 @@ void Pstream::scatterList
myComm.above(),
reinterpret_cast<char*>(receivedValues.begin()),
receivedValues.byteSize(),
tag
tag,
comm
);
forAll(notBelowLeaves, leafI)
@ -243,7 +255,14 @@ void Pstream::scatterList
}
else
{
IPstream fromAbove(UPstream::scheduled, myComm.above(), 0, tag);
IPstream fromAbove
(
UPstream::scheduled,
myComm.above(),
0,
tag,
comm
);
forAll(notBelowLeaves, leafI)
{
@ -281,12 +300,13 @@ void Pstream::scatterList
belowID,
reinterpret_cast<const char*>(sendingValues.begin()),
sendingValues.byteSize(),
tag
tag,
comm
);
}
else
{
OPstream toBelow(UPstream::scheduled, belowID, 0, tag);
OPstream toBelow(UPstream::scheduled, belowID, 0, tag, comm);
// Send data destined for all other processors below belowID
forAll(notBelowLeaves, leafI)
@ -307,16 +327,16 @@ void Pstream::scatterList
}
template<class T>
void Pstream::scatterList(List<T>& Values, const int tag)
template <class T>
void Pstream::scatterList(List<T>& Values, const int tag, const label comm)
{
if (UPstream::nProcs() < UPstream::nProcsSimpleSum)
if (UPstream::nProcs(comm) < UPstream::nProcsSimpleSum)
{
scatterList(UPstream::linearCommunication(), Values, tag);
scatterList(UPstream::linearCommunication(comm), Values, tag, comm);
}
else
{
scatterList(UPstream::treeCommunication(), Values, tag);
scatterList(UPstream::treeCommunication(comm), Values, tag, comm);
}
}

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
@ -166,6 +166,60 @@ Foam::OSstream& Foam::messageStream::operator()
}
Foam::OSstream& Foam::messageStream::operator()(const label communicator)
{
if (UPstream::warnComm != -1 && communicator != UPstream::warnComm)
{
Pout<< "** messageStream with comm:" << communicator
<< endl;
error::printStack(Pout);
}
if (communicator == UPstream::worldComm)
{
return operator()();
}
else
{
bool master = UPstream::master(communicator);
if (level)
{
bool collect = (severity_ == INFO || severity_ == WARNING);
// Report the error
if (!master && collect)
{
return Snull;
}
else
{
if (title().size())
{
Pout<< title().c_str();
}
if (maxErrors_)
{
errorCount_++;
if (errorCount_ >= maxErrors_)
{
FatalErrorIn("messageStream::operator OSstream&()")
<< "Too many errors"
<< abort(FatalError);
}
}
return Pout;
}
}
}
return Snull;
}
Foam::messageStream::operator Foam::OSstream&()
{
if (level)

View File

@ -184,6 +184,11 @@ public:
const dictionary&
);
//- Convert to OSstream
// Use Info for default communicator, use Pout
// on master for non-default one.
OSstream& operator()(const label communicator);
//- Convert to OSstream for << operations
operator OSstream&();

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
@ -215,9 +215,10 @@ bool Foam::regIOobject::read()
(
comms,
const_cast<word&>(headerClassName()),
Pstream::msgType()
Pstream::msgType(),
Pstream::worldComm
);
Pstream::scatter(comms, note(), Pstream::msgType());
Pstream::scatter(comms, note(), Pstream::msgType(), Pstream::worldComm);
// Get my communication order

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

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

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

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

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

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
@ -477,10 +477,10 @@ TMP_UNARY_FUNCTION(Type, average)
#define G_UNARY_FUNCTION(ReturnType, gFunc, Func, rFunc) \
\
template<class Type> \
ReturnType gFunc(const UList<Type>& f) \
ReturnType gFunc(const UList<Type>& f, const int comm) \
{ \
ReturnType res = Func(f); \
reduce(res, rFunc##Op<Type>()); \
reduce(res, rFunc##Op<Type>(), Pstream::msgType(), comm); \
return res; \
} \
TMP_UNARY_FUNCTION(ReturnType, gFunc)
@ -495,27 +495,41 @@ G_UNARY_FUNCTION(Type, gSumCmptMag, sumCmptMag, sum)
#undef G_UNARY_FUNCTION
template<class Type>
scalar gSumProd(const UList<Type>& f1, const UList<Type>& f2)
scalar gSumProd
(
const UList<Type>& f1,
const UList<Type>& f2,
const int comm
)
{
scalar SumProd = sumProd(f1, f2);
reduce(SumProd, sumOp<scalar>());
reduce(SumProd, sumOp<scalar>(), Pstream::msgType(), comm);
return SumProd;
}
template<class Type>
Type gSumCmptProd(const UList<Type>& f1, const UList<Type>& f2)
Type gSumCmptProd
(
const UList<Type>& f1,
const UList<Type>& f2,
const int comm
)
{
Type SumProd = sumCmptProd(f1, f2);
reduce(SumProd, sumOp<Type>());
reduce(SumProd, sumOp<Type>(), Pstream::msgType(), comm);
return SumProd;
}
template<class Type>
Type gAverage(const UList<Type>& f)
Type gAverage
(
const UList<Type>& f,
const int comm
)
{
label n = f.size();
Type s = sum(f);
sumReduce(s, n);
sumReduce(s, n, Pstream::msgType(), comm);
if (n > 0)
{

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
@ -25,6 +25,7 @@ License
#define TEMPLATE template<class Type>
#include "FieldFunctionsM.H"
#include "UPstream.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -201,7 +202,7 @@ TMP_UNARY_FUNCTION(Type, average)
#define G_UNARY_FUNCTION(ReturnType, gFunc, Func, rFunc) \
\
template<class Type> \
ReturnType gFunc(const UList<Type>& f); \
ReturnType gFunc(const UList<Type>& f, const int comm = UPstream::worldComm); \
TMP_UNARY_FUNCTION(ReturnType, gFunc)
G_UNARY_FUNCTION(Type, gMax, max, max)
@ -214,13 +215,27 @@ G_UNARY_FUNCTION(Type, gSumCmptMag, sumCmptMag, sum)
#undef G_UNARY_FUNCTION
template<class Type>
scalar gSumProd(const UList<Type>& f1, const UList<Type>& f2);
scalar gSumProd
(
const UList<Type>& f1,
const UList<Type>& f2,
const int comm = UPstream::worldComm
);
template<class Type>
Type gSumCmptProd(const UList<Type>& f1, const UList<Type>& f2);
Type gSumCmptProd
(
const UList<Type>& f1,
const UList<Type>& f2,
const int comm = UPstream::worldComm
);
template<class Type>
Type gAverage(const UList<Type>& f);
Type gAverage
(
const UList<Type>& f,
const int comm = UPstream::worldComm
);
TMP_UNARY_FUNCTION(Type, gAverage)

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

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

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
@ -29,11 +29,20 @@ License
#include "procLduInterface.H"
#include "cyclicLduInterface.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(LUscalarMatrix, 0);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::LUscalarMatrix::LUscalarMatrix(const scalarSquareMatrix& matrix)
:
scalarSquareMatrix(matrix),
comm_(Pstream::worldComm),
pivotIndices_(n())
{
LUDecompose(*this, pivotIndices_);
@ -46,10 +55,12 @@ Foam::LUscalarMatrix::LUscalarMatrix
const FieldField<Field, scalar>& interfaceCoeffs,
const lduInterfaceFieldPtrsList& interfaces
)
:
comm_(ldum.mesh().comm())
{
if (Pstream::parRun())
{
PtrList<procLduMatrix> lduMatrices(Pstream::nProcs());
PtrList<procLduMatrix> lduMatrices(Pstream::nProcs(comm_));
label lduMatrixi = 0;
@ -64,25 +75,42 @@ Foam::LUscalarMatrix::LUscalarMatrix
)
);
if (Pstream::master())
if (Pstream::master(comm_))
{
for
(
int slave=Pstream::firstSlave();
slave<=Pstream::lastSlave();
slave<=Pstream::lastSlave(comm_);
slave++
)
{
lduMatrices.set
(
lduMatrixi++,
new procLduMatrix(IPstream(Pstream::scheduled, slave)())
new procLduMatrix
(
IPstream
(
Pstream::scheduled,
slave,
0, // bufSize
Pstream::msgType(),
comm_
)()
)
);
}
}
else
{
OPstream toMaster(Pstream::scheduled, Pstream::masterNo());
OPstream toMaster
(
Pstream::scheduled,
Pstream::masterNo(),
0, // bufSize
Pstream::msgType(),
comm_
);
procLduMatrix cldum
(
ldum,
@ -93,7 +121,7 @@ Foam::LUscalarMatrix::LUscalarMatrix
}
if (Pstream::master())
if (Pstream::master(comm_))
{
label nCells = 0;
forAll(lduMatrices, i)
@ -114,8 +142,44 @@ Foam::LUscalarMatrix::LUscalarMatrix
convert(ldum, interfaceCoeffs, interfaces);
}
if (Pstream::master())
if (Pstream::master(comm_))
{
label nRows = n();
label nColumns = m();
if (debug)
{
Pout<< "LUscalarMatrix : size:" << nRows << endl;
for (label rowI = 0; rowI < nRows; rowI++)
{
const scalar* row = operator[](rowI);
Pout<< "cell:" << rowI << " diagCoeff:" << row[rowI] << endl;
Pout<< " connects to upper cells :";
for (label columnI = rowI+1; columnI < nColumns; columnI++)
{
if (mag(row[columnI]) > SMALL)
{
Pout<< ' ' << columnI << " (coeff:" << row[columnI]
<< ")";
}
}
Pout<< endl;
Pout<< " connects to lower cells :";
for (label columnI = 0; columnI < rowI; columnI++)
{
if (mag(row[columnI]) > SMALL)
{
Pout<< ' ' << columnI << " (coeff:" << row[columnI]
<< ")";
}
}
Pout<< endl;
}
Pout<< endl;
}
pivotIndices_.setSize(n());
LUDecompose(*this, pivotIndices_);
}

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
@ -58,6 +58,9 @@ class LUscalarMatrix
{
// Private data
//- Communicator to use
const label comm_;
//- Processor matrix offsets
labelList procOffsets_;
@ -84,6 +87,9 @@ class LUscalarMatrix
public:
// Declare name of the class and its debug switch
ClassName("LUscalarMatrix");
// Constructors
//- Construct from scalarSquareMatrix and perform LU decomposition

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
@ -34,7 +34,7 @@ void Foam::LUscalarMatrix::solve(Field<Type>& sourceSol) const
{
Field<Type> completeSourceSol(n());
if (Pstream::master())
if (Pstream::master(comm_))
{
typename Field<Type>::subField
(
@ -45,7 +45,7 @@ void Foam::LUscalarMatrix::solve(Field<Type>& sourceSol) const
for
(
int slave=Pstream::firstSlave();
slave<=Pstream::lastSlave();
slave<=Pstream::lastSlave(comm_);
slave++
)
{
@ -57,7 +57,9 @@ void Foam::LUscalarMatrix::solve(Field<Type>& sourceSol) const
(
&(completeSourceSol[procOffsets_[slave]])
),
(procOffsets_[slave + 1] - procOffsets_[slave])*sizeof(Type)
(procOffsets_[slave+1]-procOffsets_[slave])*sizeof(Type),
Pstream::msgType(),
comm_
);
}
}
@ -68,11 +70,13 @@ void Foam::LUscalarMatrix::solve(Field<Type>& sourceSol) const
Pstream::scheduled,
Pstream::masterNo(),
reinterpret_cast<const char*>(sourceSol.begin()),
sourceSol.byteSize()
sourceSol.byteSize(),
Pstream::msgType(),
comm_
);
}
if (Pstream::master())
if (Pstream::master(comm_))
{
LUBacksubstitute(*this, pivotIndices_, completeSourceSol);
@ -85,7 +89,7 @@ void Foam::LUscalarMatrix::solve(Field<Type>& sourceSol) const
for
(
int slave=Pstream::firstSlave();
slave<=Pstream::lastSlave();
slave<=Pstream::lastSlave(comm_);
slave++
)
{
@ -97,7 +101,9 @@ void Foam::LUscalarMatrix::solve(Field<Type>& sourceSol) const
(
&(completeSourceSol[procOffsets_[slave]])
),
(procOffsets_[slave + 1] - procOffsets_[slave])*sizeof(Type)
(procOffsets_[slave + 1]-procOffsets_[slave])*sizeof(Type),
Pstream::msgType(),
comm_
);
}
}
@ -108,7 +114,9 @@ void Foam::LUscalarMatrix::solve(Field<Type>& sourceSol) const
Pstream::scheduled,
Pstream::masterNo(),
reinterpret_cast<char*>(sourceSol.begin()),
sourceSol.byteSize()
sourceSol.byteSize(),
Pstream::msgType(),
comm_
);
}
}

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
@ -40,7 +40,8 @@ Foam::procLduInterface::procLduInterface
coeffs_(coeffs),
myProcNo_(-1),
neighbProcNo_(-1),
tag_(-1)
tag_(-1),
comm_(-1)
{
if (isA<processorLduInterface>(interface.interface()))
{
@ -50,6 +51,7 @@ Foam::procLduInterface::procLduInterface
myProcNo_ = pldui.myProcNo();
neighbProcNo_ = pldui.neighbProcNo();
tag_ = pldui.tag();
comm_ = pldui.comm();
}
else if (isA<cyclicLduInterface>(interface.interface()))
{
@ -73,7 +75,8 @@ Foam::procLduInterface::procLduInterface(Istream& is)
coeffs_(is),
myProcNo_(readLabel(is)),
neighbProcNo_(readLabel(is)),
tag_(readLabel(is))
tag_(readLabel(is)),
comm_(readLabel(is))
{}
@ -85,7 +88,8 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const procLduInterface& cldui)
<< cldui.coeffs_
<< cldui.myProcNo_
<< cldui.neighbProcNo_
<< cldui.tag_;
<< cldui.tag_
<< cldui.comm_;
return os;
}

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
@ -25,7 +25,7 @@ Class
Foam::procLduInterface
Description
Foam::procLduInterface
IO interface for processorLduInterface
SourceFiles
procLduInterface.C
@ -58,6 +58,7 @@ class procLduInterface
label myProcNo_;
label neighbProcNo_;
label tag_;
label comm_;
// Private Member Functions

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
@ -25,7 +25,7 @@ Class
Foam::procLduMatrix
Description
Foam::procLduMatrix
I/O for lduMatrix and interface values.
SourceFiles
procLduMatrix.C

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
@ -24,6 +24,7 @@ License
\*---------------------------------------------------------------------------*/
#include "SolverPerformance.H"
#include "IOstreams.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -64,7 +65,8 @@ bool Foam::SolverPerformance<Type>::checkConvergence
{
if (debug >= 2)
{
Info<< solverName_
//Info<< solverName_
Pout<< solverName_
<< ": Iteration " << noIterations_
<< " residual = " << finalResidual_
<< endl;

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
@ -39,7 +39,7 @@ Description
order but only for groups of edges belonging to each point. An example
is given below:
\verbatim
owner eighbour
owner neighbour
0 1
0 20
1 2

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
@ -84,10 +84,13 @@ public:
// Access
//- Return processor number
//- Return communicator used for parallel communication
virtual int comm() const = 0;
//- Return processor number (rank in communicator)
virtual int myProcNo() const = 0;
//- Return neigbour processor number
//- Return neigbour processor number (rank in communicator)
virtual int neighbProcNo() const = 0;
//- Return face transformation tensor

View File

@ -46,7 +46,8 @@ void Foam::processorLduInterface::send
neighbProcNo(),
reinterpret_cast<const char*>(f.begin()),
nBytes,
tag()
tag(),
comm()
);
}
else if (commsType == Pstream::nonBlocking)
@ -59,7 +60,8 @@ void Foam::processorLduInterface::send
neighbProcNo(),
receiveBuf_.begin(),
nBytes,
tag()
tag(),
comm()
);
resizeBuf(sendBuf_, nBytes);
@ -71,7 +73,8 @@ void Foam::processorLduInterface::send
neighbProcNo(),
sendBuf_.begin(),
nBytes,
tag()
tag(),
comm()
);
}
else
@ -98,7 +101,8 @@ void Foam::processorLduInterface::receive
neighbProcNo(),
reinterpret_cast<char*>(f.begin()),
f.byteSize(),
tag()
tag(),
comm()
);
}
else if (commsType == Pstream::nonBlocking)
@ -162,7 +166,8 @@ void Foam::processorLduInterface::compressedSend
neighbProcNo(),
sendBuf_.begin(),
nBytes,
tag()
tag(),
comm()
);
}
else if (commsType == Pstream::nonBlocking)
@ -175,7 +180,8 @@ void Foam::processorLduInterface::compressedSend
neighbProcNo(),
receiveBuf_.begin(),
nBytes,
tag()
tag(),
comm()
);
OPstream::write
@ -184,7 +190,8 @@ void Foam::processorLduInterface::compressedSend
neighbProcNo(),
sendBuf_.begin(),
nBytes,
tag()
tag(),
comm()
);
}
else
@ -225,7 +232,8 @@ void Foam::processorLduInterface::compressedReceive
neighbProcNo(),
receiveBuf_.begin(),
nBytes,
tag()
tag(),
comm()
);
}
else if (commsType != Pstream::nonBlocking)

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
@ -71,6 +71,9 @@ public:
// Access
//- Return communicator used for comms
virtual int comm() const = 0;
//- Return processor number
virtual int myProcNo() const = 0;

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
@ -25,6 +25,7 @@ License
#include "lduMatrix.H"
#include "IOstreams.H"
#include "Switch.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -116,17 +117,30 @@ Foam::lduMatrix::lduMatrix(lduMatrix& A, bool reUse)
}
Foam::lduMatrix::lduMatrix
(
const lduMesh& mesh,
Istream& is
)
Foam::lduMatrix::lduMatrix(const lduMesh& mesh, Istream& is)
:
lduMesh_(mesh),
lowerPtr_(new scalarField(is)),
diagPtr_(new scalarField(is)),
upperPtr_(new scalarField(is))
{}
lowerPtr_(NULL),
diagPtr_(NULL),
upperPtr_(NULL)
{
Switch hasLow(is);
Switch hasDiag(is);
Switch hasUp(is);
if (hasLow)
{
lowerPtr_ = new scalarField(is);
}
if (hasDiag)
{
diagPtr_ = new scalarField(is);
}
if (hasUp)
{
upperPtr_ = new scalarField(is);
}
}
Foam::lduMatrix::~lduMatrix()
@ -195,6 +209,53 @@ Foam::scalarField& Foam::lduMatrix::upper()
}
Foam::scalarField& Foam::lduMatrix::lower(const label nCoeffs)
{
if (!lowerPtr_)
{
if (upperPtr_)
{
lowerPtr_ = new scalarField(*upperPtr_);
}
else
{
lowerPtr_ = new scalarField(nCoeffs, 0.0);
}
}
return *lowerPtr_;
}
Foam::scalarField& Foam::lduMatrix::diag(const label size)
{
if (!diagPtr_)
{
diagPtr_ = new scalarField(size, 0.0);
}
return *diagPtr_;
}
Foam::scalarField& Foam::lduMatrix::upper(const label nCoeffs)
{
if (!upperPtr_)
{
if (lowerPtr_)
{
upperPtr_ = new scalarField(*lowerPtr_);
}
else
{
upperPtr_ = new scalarField(nCoeffs, 0.0);
}
}
return *upperPtr_;
}
const Foam::scalarField& Foam::lduMatrix::lower() const
{
if (!lowerPtr_ && !upperPtr_)
@ -252,25 +313,26 @@ const Foam::scalarField& Foam::lduMatrix::upper() const
Foam::Ostream& Foam::operator<<(Ostream& os, const lduMatrix& ldum)
{
if (ldum.lowerPtr_)
Switch hasLow = ldum.hasLower();
Switch hasDiag = ldum.hasDiag();
Switch hasUp = ldum.hasUpper();
os << hasLow << token::SPACE << hasDiag << token::SPACE
<< hasUp << token::SPACE;
if (hasLow)
{
os << "Lower triangle = "
<< *ldum.lowerPtr_
<< endl << endl;
os << ldum.lower();
}
if (ldum.diagPtr_)
if (hasDiag)
{
os << "diagonal = "
<< *ldum.diagPtr_
<< endl << endl;
os << ldum.diag();
}
if (ldum.upperPtr_)
if (hasUp)
{
os << "Upper triangle = "
<< *ldum.upperPtr_
<< endl << endl;
os << ldum.upper();
}
os.check("Ostream& operator<<(Ostream&, const lduMatrix&");
@ -279,4 +341,64 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const lduMatrix& ldum)
}
Foam::Ostream& Foam::operator<<(Ostream& os, const InfoProxy<lduMatrix>& ip)
{
const lduMatrix& ldum = ip.t_;
Switch hasLow = ldum.hasLower();
Switch hasDiag = ldum.hasDiag();
Switch hasUp = ldum.hasUpper();
os << "Lower:" << hasLow
<< " Diag:" << hasDiag
<< " Upper:" << hasUp << endl;
if (hasLow)
{
os << "lower:" << ldum.lower().size() << endl;
}
if (hasDiag)
{
os << "diag :" << ldum.diag().size() << endl;
}
if (hasUp)
{
os << "upper:" << ldum.upper().size() << endl;
}
//if (hasLow)
//{
// os << "lower contents:" << endl;
// forAll(ldum.lower(), i)
// {
// os << "i:" << i << "\t" << ldum.lower()[i] << endl;
// }
// os << endl;
//}
//if (hasDiag)
//{
// os << "diag contents:" << endl;
// forAll(ldum.diag(), i)
// {
// os << "i:" << i << "\t" << ldum.diag()[i] << endl;
// }
// os << endl;
//}
//if (hasUp)
//{
// os << "upper contents:" << endl;
// forAll(ldum.upper(), i)
// {
// os << "i:" << i << "\t" << ldum.upper()[i] << endl;
// }
// os << endl;
//}
os.check("Ostream& operator<<(Ostream&, const lduMatrix&");
return os;
}
// ************************************************************************* //

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
@ -58,6 +58,7 @@ SourceFiles
#include "autoPtr.H"
#include "runTimeSelectionTables.H"
#include "solverPerformance.H"
#include "InfoProxy.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -562,6 +563,14 @@ public:
scalarField& diag();
scalarField& upper();
// Size with externally provided sizes (for constructing with 'fake'
// mesh in GAMG)
scalarField& lower(const label size);
scalarField& diag(const label nCoeffs);
scalarField& upper(const label nCoeffs);
const scalarField& lower() const;
const scalarField& diag() const;
const scalarField& upper() const;
@ -691,6 +700,16 @@ public:
tmp<Field<Type> > faceH(const tmp<Field<Type> >&) const;
// Info
//- Return info proxy.
// Used to print matrix information to a stream
InfoProxy<lduMatrix> info() const
{
return *this;
}
// Member operators
void operator=(const lduMatrix&);
@ -707,6 +726,7 @@ public:
// Ostream operator
friend Ostream& operator<<(Ostream&, const lduMatrix&);
friend Ostream& operator<<(Ostream&, const InfoProxy<lduMatrix>&);
};

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
@ -186,10 +186,15 @@ Foam::scalar Foam::lduMatrix::solver::normFactor
{
// --- Calculate A dot reference value of psi
matrix_.sumA(tmpField, interfaceBouCoeffs_, interfaces_);
tmpField *= gAverage(psi);
tmpField *= gAverage(psi, matrix_.lduMesh_.comm());
return
gSum(mag(Apsi - tmpField) + mag(source - tmpField))
gSum
(
(mag(Apsi - tmpField) + mag(source - tmpField))(),
matrix_.lduMesh_.comm()
)
+ solverPerformance::small_;
// At convergence this simpler method is equivalent to the above

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
@ -99,8 +99,20 @@ void Foam::GAMGPreconditioner::precondition
// Create the smoothers for all levels
PtrList<lduMatrix::smoother> smoothers;
// Scratch fields if processor-agglomerated coarse level meshes
// are bigger than original. Usually not needed
scalarField ApsiScratch;
scalarField finestCorrectionScratch;
// Initialise the above data structures
initVcycle(coarseCorrFields, coarseSources, smoothers);
initVcycle
(
coarseCorrFields,
coarseSources,
smoothers,
ApsiScratch,
finestCorrectionScratch
);
for (label cycle=0; cycle<nVcycles_; cycle++)
{
@ -112,6 +124,14 @@ void Foam::GAMGPreconditioner::precondition
AwA,
finestCorrection,
finestResidual,
(ApsiScratch.size() ? ApsiScratch : AwA),
(
finestCorrectionScratch.size()
? finestCorrectionScratch
: finestCorrection
),
coarseCorrFields,
coarseSources,
cmpt

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
@ -25,6 +25,7 @@ License
#include "GAMGAgglomeration.H"
#include "GAMGInterface.H"
#include "processorGAMGInterface.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
@ -85,7 +86,8 @@ void Foam::GAMGAgglomeration::agglomerateLduAddressing
labelList initCoarseNeighb(nFineFaces);
// Counter for coarse faces
label nCoarseFaces = 0;
label& nCoarseFaces = nFaces_[fineLevelIndex];
nCoarseFaces = 0;
// Loop through all fine faces
forAll(upperAddr, fineFacei)
@ -197,6 +199,58 @@ void Foam::GAMGAgglomeration::agglomerateLduAddressing
}
}
// Create face-flip status
faceFlipMap_.set(fineLevelIndex, new boolList(nFineFaces, false));
boolList& faceFlipMap = faceFlipMap_[fineLevelIndex];
label nFlipped = 0;
label nDissapear = 0;
forAll(faceRestrictAddr, fineFacei)
{
label coarseFacei = faceRestrictAddr[fineFacei];
if (coarseFacei >= 0)
{
// Maps to coarse face
label cOwn = coarseOwner[coarseFacei];
label cNei = coarseNeighbour[coarseFacei];
label rmUpperAddr = restrictMap[upperAddr[fineFacei]];
label rmLowerAddr = restrictMap[lowerAddr[fineFacei]];
if (cOwn == rmUpperAddr && cNei == rmLowerAddr)
{
faceFlipMap[fineFacei] = true;
nFlipped++;
}
else if (cOwn == rmLowerAddr && cNei == rmUpperAddr)
{
//faceFlipMap[fineFacei] = false;
}
else
{
FatalErrorIn("GAMGAgglomeration::agglomerateLduAddressing(..)")
<< "problem."
<< " fineFacei:" << fineFacei
<< " rmUpperAddr:" << rmUpperAddr
<< " rmLowerAddr:" << rmLowerAddr
<< " coarseFacei:" << coarseFacei
<< " cOwn:" << cOwn
<< " cNei:" << cNei
<< exit(FatalError);
}
}
else
{
nDissapear++;
}
}
// Clear the temporary storage for the coarse cell data
cCellnFaces.setSize(0);
cCellFaces.setSize(0);
@ -207,20 +261,19 @@ void Foam::GAMGAgglomeration::agglomerateLduAddressing
// Create coarse-level interfaces
// Get reference to fine-level interfaces
const lduInterfacePtrsList& fineInterfaces =
interfaceLevels_[fineLevelIndex];
const lduInterfacePtrsList& fineInterfaces = interfaceLevel(fineLevelIndex);
// Create coarse-level interfaces
interfaceLevels_.set
nPatchFaces_.set(fineLevelIndex, new labelList(fineInterfaces.size(), 0));
labelList& nPatchFaces = nPatchFaces_[fineLevelIndex];
patchFaceRestrictAddressing_.set
(
fineLevelIndex + 1,
new lduInterfacePtrsList(fineInterfaces.size())
fineLevelIndex,
new labelListList(fineInterfaces.size())
);
labelListList& patchFineToCoarse =
patchFaceRestrictAddressing_[fineLevelIndex];
lduInterfacePtrsList& coarseInterfaces =
interfaceLevels_[fineLevelIndex + 1];
labelListList coarseInterfaceAddr(fineInterfaces.size());
// Initialise transfer of restrict addressing on the interface
forAll(fineInterfaces, inti)
@ -240,7 +293,23 @@ void Foam::GAMGAgglomeration::agglomerateLduAddressing
Pstream::waitRequests();
}
// Add the coarse level
meshLevels_.set
(
fineLevelIndex,
new lduPrimitiveMesh
(
nCoarseCells,
coarseOwner,
coarseNeighbour,
fineMesh.comm(),
true
)
);
lduInterfacePtrsList coarseInterfaces(fineInterfaces.size());
forAll(fineInterfaces, inti)
{
if (fineInterfaces.set(inti))
@ -251,7 +320,7 @@ void Foam::GAMGAgglomeration::agglomerateLduAddressing
GAMGInterface::New
(
inti,
coarseInterfaces,
meshLevels_[fineLevelIndex].rawInterfaces(),
fineInterfaces[inti],
fineInterfaces[inti].interfaceInternalField(restrictMap),
fineInterfaces[inti].internalFieldTransfer
@ -259,30 +328,394 @@ void Foam::GAMGAgglomeration::agglomerateLduAddressing
Pstream::nonBlocking,
restrictMap
),
fineLevelIndex
fineLevelIndex,
fineMesh.comm()
).ptr()
);
coarseInterfaceAddr[inti] = coarseInterfaces[inti].faceCells();
nPatchFaces[inti] = coarseInterfaces[inti].faceCells().size();
patchFineToCoarse[inti] = refCast<const GAMGInterface>
(
coarseInterfaces[inti]
).faceRestrictAddressing();
}
}
// Set the coarse ldu addressing onto the list
meshLevels_.set
meshLevels_[fineLevelIndex].addInterfaces
(
fineLevelIndex,
new lduPrimitiveMesh
coarseInterfaces,
lduPrimitiveMesh::nonBlockingSchedule<processorGAMGInterface>
(
nCoarseCells,
coarseOwner,
coarseNeighbour,
coarseInterfaceAddr,
coarseInterfaces,
fineMeshAddr.patchSchedule(),
true
coarseInterfaces
)
);
if (debug)
{
Pout<< "GAMGAgglomeration :"
<< " agglomerated level " << fineLevelIndex
<< " from nCells:" << fineMeshAddr.size()
<< " nFaces:" << upperAddr.size()
<< " to nCells:" << nCoarseCells
<< " nFaces:" << nCoarseFaces
<< endl;
}
}
void Foam::GAMGAgglomeration::procAgglomerateLduAddressing
(
const label meshComm,
const labelList& procAgglomMap,
const labelList& procIDs,
const label allMeshComm,
const label levelIndex
)
{
const lduMesh& myMesh = meshLevels_[levelIndex-1];
label oldWarn = UPstream::warnComm;
UPstream::warnComm = meshComm;
procAgglomMap_.set(levelIndex, new labelList(procAgglomMap));
agglomProcIDs_.set(levelIndex, new labelList(procIDs));
procCommunicator_[levelIndex] = allMeshComm;
// These could only be set on the master procs but it is
// quite convenient to also have them on the slaves
procCellOffsets_.set(levelIndex, new labelList(0));
procFaceMap_.set(levelIndex, new labelListList(0));
procBoundaryMap_.set(levelIndex, new labelListList(0));
procBoundaryFaceMap_.set(levelIndex, new labelListListList(0));
// Collect meshes
PtrList<lduPrimitiveMesh> otherMeshes;
lduPrimitiveMesh::gather(meshComm, myMesh, procIDs, otherMeshes);
if (Pstream::myProcNo(meshComm) == procIDs[0])
{
// Combine all addressing
labelList procFaceOffsets;
meshLevels_.set
(
levelIndex-1,
new lduPrimitiveMesh
(
allMeshComm,
procAgglomMap,
procIDs,
myMesh,
otherMeshes,
procCellOffsets_[levelIndex],
procFaceOffsets,
procFaceMap_[levelIndex],
procBoundaryMap_[levelIndex],
procBoundaryFaceMap_[levelIndex]
)
);
}
// Combine restrict addressing
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
procAgglomerateRestrictAddressing
(
meshComm,
procIDs,
levelIndex
);
if (Pstream::myProcNo(meshComm) != procIDs[0])
{
clearLevel(levelIndex);
}
UPstream::warnComm = oldWarn;
}
void Foam::GAMGAgglomeration::procAgglomerateRestrictAddressing
(
const label comm,
const labelList& procIDs,
const label levelIndex
)
{
// Collect number of cells
labelList nFineCells;
gatherList
(
comm,
procIDs,
restrictAddressing_[levelIndex].size(),
nFineCells
);
labelList offsets(nFineCells.size()+1);
{
offsets[0] = 0;
forAll(nFineCells, i)
{
offsets[i+1] = offsets[i] + nFineCells[i];
}
}
// Combine and renumber nCoarseCells
labelList nCoarseCells;
gatherList
(
comm,
procIDs,
nCells_[levelIndex],
nCoarseCells
);
// (cell)restrictAddressing
const globalIndex cellOffsetter(offsets);
labelList procRestrictAddressing;
cellOffsetter.gather
(
comm,
procIDs,
restrictAddressing_[levelIndex],
procRestrictAddressing
);
if (Pstream::myProcNo(comm) == procIDs[0])
{
labelList coarseCellOffsets(procIDs.size()+1);
{
coarseCellOffsets[0] = 0;
forAll(procIDs, i)
{
coarseCellOffsets[i+1] = coarseCellOffsets[i]+nCoarseCells[i];
}
}
nCells_[levelIndex] = coarseCellOffsets.last();
// Renumber consecutively
for (label procI = 1; procI < procIDs.size(); procI++)
{
SubList<label> procSlot
(
procRestrictAddressing,
offsets[procI+1]-offsets[procI],
offsets[procI]
);
forAll(procSlot, i)
{
procSlot[i] += coarseCellOffsets[procI];
}
}
restrictAddressing_[levelIndex].transfer(procRestrictAddressing);
}
}
void Foam::GAMGAgglomeration::combineLevels(const label curLevel)
{
label prevLevel = curLevel - 1;
// Set the previous level nCells to the current
nCells_[prevLevel] = nCells_[curLevel];
nFaces_[prevLevel] = nFaces_[curLevel];
// Map the restrictAddressing from the coarser level into the previous
// finer level
const labelList& curResAddr = restrictAddressing_[curLevel];
labelList& prevResAddr = restrictAddressing_[prevLevel];
const labelList& curFaceResAddr = faceRestrictAddressing_[curLevel];
labelList& prevFaceResAddr = faceRestrictAddressing_[prevLevel];
const boolList& curFaceFlipMap = faceFlipMap_[curLevel];
boolList& prevFaceFlipMap = faceFlipMap_[prevLevel];
forAll(prevFaceResAddr, i)
{
if (prevFaceResAddr[i] >= 0)
{
label fineFaceI = prevFaceResAddr[i];
prevFaceResAddr[i] = curFaceResAddr[fineFaceI];
prevFaceFlipMap[i] = curFaceFlipMap[fineFaceI];
}
else
{
label fineFaceI = -prevFaceResAddr[i] - 1;
prevFaceResAddr[i] = -curResAddr[fineFaceI] - 1;
prevFaceFlipMap[i] = curFaceFlipMap[fineFaceI];
}
}
// Delete the restrictAddressing for the coarser level
faceRestrictAddressing_.set(curLevel, NULL);
faceFlipMap_.set(curLevel, NULL);
forAll(prevResAddr, i)
{
prevResAddr[i] = curResAddr[prevResAddr[i]];
}
const labelListList& curPatchFaceResAddr =
patchFaceRestrictAddressing_[curLevel];
labelListList& prevPatchFaceResAddr =
patchFaceRestrictAddressing_[prevLevel];
forAll(prevPatchFaceResAddr, inti)
{
const labelList& curResAddr = curPatchFaceResAddr[inti];
labelList& prevResAddr = prevPatchFaceResAddr[inti];
forAll(prevResAddr, i)
{
label fineFaceI = prevResAddr[i];
prevResAddr[i] = curResAddr[fineFaceI];
}
}
// Delete the restrictAddressing for the coarser level
restrictAddressing_.set(curLevel, NULL);
// Patch faces
nPatchFaces_[prevLevel] = nPatchFaces_[curLevel];
// Adapt the restrict addressing for the patches
const lduInterfacePtrsList& curInterLevel =
meshLevels_[curLevel].rawInterfaces();
const lduInterfacePtrsList& prevInterLevel =
meshLevels_[prevLevel].rawInterfaces();
forAll(prevInterLevel, inti)
{
if (prevInterLevel.set(inti))
{
GAMGInterface& prevInt = refCast<GAMGInterface>
(
const_cast<lduInterface&>
(
prevInterLevel[inti]
)
);
const GAMGInterface& curInt = refCast<const GAMGInterface>
(
curInterLevel[inti]
);
prevInt.combine(curInt);
}
}
// Delete the matrix addressing and coefficients from the previous level
// and replace with the corresponding entry from the coarser level
meshLevels_.set(prevLevel, meshLevels_.set(curLevel, NULL));
}
//void Foam::GAMGAgglomeration::gatherList
//(
// const label comm,
// const labelList& procIDs,
//
// const label myVal,
// labelList& vals,
// const int tag
//)
//{
// vals.setSize(procIDs.size());
//
// if (Pstream::myProcNo(comm) == procIDs[0])
// {
// vals[0] = myVal;
//
// for (label i = 1; i < procIDs.size(); i++)
// {
// label& slaveVal = vals[i];
// IPstream::read
// (
// Pstream::scheduled,
// procIDs[i],
// reinterpret_cast<char*>(&slaveVal),
// sizeof(slaveVal),
// tag,
// comm
// );
// }
// }
// else
// {
// OPstream::write
// (
// Pstream::scheduled,
// procIDs[0],
// reinterpret_cast<const char*>(&myVal),
// sizeof(myVal),
// tag,
// comm
// );
// }
//}
void Foam::GAMGAgglomeration::calculateRegionMaster
(
const label comm,
const labelList& procAgglomMap,
labelList& masterProcs,
List<int>& agglomProcIDs
)
{
// Determine the master processors
Map<label> agglomToMaster(procAgglomMap.size());
forAll(procAgglomMap, procI)
{
label coarseI = procAgglomMap[procI];
Map<label>::iterator fnd = agglomToMaster.find(coarseI);
if (fnd == agglomToMaster.end())
{
agglomToMaster.insert(coarseI, procI);
}
else
{
fnd() = min(fnd(), procI);
}
}
masterProcs.setSize(agglomToMaster.size());
forAllConstIter(Map<label>, agglomToMaster, iter)
{
masterProcs[iter.key()] = iter();
}
// Collect all the processors in my agglomeration
label myProcID = Pstream::myProcNo(comm);
label myAgglom = procAgglomMap[myProcID];
// Get all processors agglomerating to the same coarse
// processor
agglomProcIDs = findIndices(procAgglomMap, myAgglom);
// Make sure the master is the first element.
label index = findIndex
(
agglomProcIDs,
agglomToMaster[myAgglom]
);
Swap(agglomProcIDs[0], agglomProcIDs[index]);
}

View File

@ -28,6 +28,8 @@ License
#include "lduMatrix.H"
#include "Time.H"
#include "dlLibraryTable.H"
#include "GAMGInterface.H"
#include "GAMGProcAgglomeration.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -36,6 +38,7 @@ namespace Foam
defineTypeNameAndDebug(GAMGAgglomeration, 0);
defineRunTimeSelectionTable(GAMGAgglomeration, lduMesh);
defineRunTimeSelectionTable(GAMGAgglomeration, lduMatrix);
defineRunTimeSelectionTable(GAMGAgglomeration, geometry);
}
@ -45,8 +48,44 @@ void Foam::GAMGAgglomeration::compactLevels(const label nCreatedLevels)
{
nCells_.setSize(nCreatedLevels);
restrictAddressing_.setSize(nCreatedLevels);
nFaces_.setSize(nCreatedLevels);
faceRestrictAddressing_.setSize(nCreatedLevels);
faceFlipMap_.setSize(nCreatedLevels);
nPatchFaces_.setSize(nCreatedLevels);
patchFaceRestrictAddressing_.setSize(nCreatedLevels);
meshLevels_.setSize(nCreatedLevels);
interfaceLevels_.setSize(nCreatedLevels + 1);
// Have procCommunicator_ always, even if not procAgglomerating
procCommunicator_.setSize(nCreatedLevels + 1);
if (processorAgglomerate())
{
procAgglomMap_.setSize(nCreatedLevels);
agglomProcIDs_.setSize(nCreatedLevels);
procCellOffsets_.setSize(nCreatedLevels);
procFaceMap_.setSize(nCreatedLevels);
procBoundaryMap_.setSize(nCreatedLevels);
procBoundaryFaceMap_.setSize(nCreatedLevels);
procAgglomeratorPtr_().agglomerate();
}
if (debug)
{
for (label levelI = 0; levelI <= size(); levelI++)
{
if (hasMeshLevel(levelI))
{
const lduMesh& fineMesh = meshLevel(levelI);
Pout<< "Level " << levelI << " fine mesh:"<< nl;
Pout<< fineMesh.info() << endl;
}
else
{
Pout<< "Level " << levelI << " has no fine mesh:" << nl
<< endl;
}
}
}
}
@ -57,7 +96,7 @@ bool Foam::GAMGAgglomeration::continueAgglomerating
{
// Check the need for further agglomeration on all processors
bool contAgg = nCoarseCells >= nCellsInCoarsestLevel_;
reduce(contAgg, andOp<bool>());
mesh().reduce(contAgg, andOp<bool>());
return contAgg;
}
@ -78,14 +117,42 @@ Foam::GAMGAgglomeration::GAMGAgglomeration
(
readLabel(controlDict.lookup("nCellsInCoarsestLevel"))
),
procAgglomeratorPtr_
(
(
(UPstream::nProcs(mesh.comm()) > 1)
&& controlDict.found("processorAgglomerator")
)
? GAMGProcAgglomeration::New
(
controlDict.lookup("processorAgglomerator"),
*this,
controlDict
)
: autoPtr<GAMGProcAgglomeration>(NULL)
),
nCells_(maxLevels_),
restrictAddressing_(maxLevels_),
nFaces_(maxLevels_),
faceRestrictAddressing_(maxLevels_),
faceFlipMap_(maxLevels_),
nPatchFaces_(maxLevels_),
patchFaceRestrictAddressing_(maxLevels_),
meshLevels_(maxLevels_),
interfaceLevels_(maxLevels_ + 1)
{}
meshLevels_(maxLevels_)
{
procCommunicator_.setSize(maxLevels_ + 1, -1);
if (processorAgglomerate())
{
procAgglomMap_.setSize(maxLevels_);
agglomProcIDs_.setSize(maxLevels_);
procCellOffsets_.setSize(maxLevels_);
procFaceMap_.setSize(maxLevels_);
procBoundaryMap_.setSize(maxLevels_);
procBoundaryFaceMap_.setSize(maxLevels_);
}
}
const Foam::GAMGAgglomeration& Foam::GAMGAgglomeration::New
@ -192,25 +259,56 @@ const Foam::GAMGAgglomeration& Foam::GAMGAgglomeration::New
}
Foam::autoPtr<Foam::GAMGAgglomeration> Foam::GAMGAgglomeration::New
(
const lduMesh& mesh,
const scalarField& cellVolumes,
const vectorField& faceAreas,
const dictionary& controlDict
)
{
const word agglomeratorType(controlDict.lookup("agglomerator"));
const_cast<Time&>(mesh.thisDb().time()).libs().open
(
controlDict,
"geometricGAMGAgglomerationLibs",
geometryConstructorTablePtr_
);
geometryConstructorTable::iterator cstrIter =
geometryConstructorTablePtr_->find(agglomeratorType);
if (cstrIter == geometryConstructorTablePtr_->end())
{
FatalErrorIn
(
"GAMGAgglomeration::New"
"(const lduMesh& mesh, const dictionary& controlDict)"
) << "Unknown GAMGAgglomeration type "
<< agglomeratorType << ".\n"
<< "Valid geometric GAMGAgglomeration types are :"
<< geometryConstructorTablePtr_->sortedToc()
<< exit(FatalError);
}
return autoPtr<GAMGAgglomeration>
(
cstrIter()
(
mesh,
cellVolumes,
faceAreas,
controlDict
)
);
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::GAMGAgglomeration::~GAMGAgglomeration()
{
// Clear the interface storage by hand.
// It is a list of ptrs not a PtrList for consistency of the interface
for (label leveli=1; leveli<interfaceLevels_.size(); leveli++)
{
lduInterfacePtrsList& curLevel = interfaceLevels_[leveli];
forAll(curLevel, i)
{
if (curLevel.set(i))
{
delete curLevel(i);
}
}
}
}
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
@ -231,12 +329,118 @@ const Foam::lduMesh& Foam::GAMGAgglomeration::meshLevel
}
bool Foam::GAMGAgglomeration::hasMeshLevel(const label i) const
{
if (i == 0)
{
return true;
}
else
{
return meshLevels_.set(i - 1);
}
}
const Foam::lduInterfacePtrsList& Foam::GAMGAgglomeration::interfaceLevel
(
const label i
) const
{
return interfaceLevels_[i];
if (i == 0)
{
return meshInterfaces_;
}
else
{
return meshLevels_[i - 1].rawInterfaces();
}
}
void Foam::GAMGAgglomeration::clearLevel(const label i)
{
if (hasMeshLevel(i))
{
meshLevels_.set(i - 1, NULL);
if (i < nCells_.size())
{
nCells_[i] = -555;
restrictAddressing_.set(i, NULL);
nFaces_[i] = -666;
faceRestrictAddressing_.set(i, NULL);
faceFlipMap_.set(i, NULL);
nPatchFaces_.set(i, NULL);
patchFaceRestrictAddressing_.set(i, NULL);
}
}
}
const Foam::labelList& Foam::GAMGAgglomeration::procAgglomMap
(
const label leveli
) const
{
return procAgglomMap_[leveli];
}
const Foam::labelList& Foam::GAMGAgglomeration::agglomProcIDs
(
const label leveli
) const
{
return agglomProcIDs_[leveli];
}
bool Foam::GAMGAgglomeration::hasProcMesh(const label leveli) const
{
return procCommunicator_[leveli] != -1;
}
Foam::label Foam::GAMGAgglomeration::procCommunicator(const label leveli) const
{
return procCommunicator_[leveli];
}
const Foam::labelList& Foam::GAMGAgglomeration::cellOffsets
(
const label leveli
) const
{
return procCellOffsets_[leveli];
}
const Foam::labelListList& Foam::GAMGAgglomeration::faceMap
(
const label leveli
) const
{
return procFaceMap_[leveli];
}
const Foam::labelListList& Foam::GAMGAgglomeration::boundaryMap
(
const label leveli
) const
{
return procBoundaryMap_[leveli];
}
const Foam::labelListListList& Foam::GAMGAgglomeration::boundaryFaceMap
(
const label leveli
) const
{
return procBoundaryFaceMap_[leveli];
}

View File

@ -30,7 +30,6 @@ Description
SourceFiles
GAMGAgglomeration.C
GAMGAgglomerationTemplates.C
GAMGAgglomerate.C
GAMGAgglomerateLduAddressing.C
\*---------------------------------------------------------------------------*/
@ -44,6 +43,8 @@ SourceFiles
#include "primitiveFields.H"
#include "runTimeSelectionTables.H"
#include "boolList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
@ -51,6 +52,8 @@ namespace Foam
class lduMesh;
class lduMatrix;
class mapDistribute;
class GAMGProcAgglomeration;
/*---------------------------------------------------------------------------*\
Class GAMGAgglomeration Declaration
@ -65,10 +68,15 @@ protected:
// Protected data
//- Max number of levels
label maxLevels_;
const label maxLevels_;
//- Number of cells in coarsest level
label nCellsInCoarsestLevel_;
const label nCellsInCoarsestLevel_;
autoPtr<GAMGProcAgglomeration> procAgglomeratorPtr_;
//- Cached mesh interfaces
lduInterfacePtrsList meshInterfaces_;
//- The number of cells in each level
labelList nCells_;
@ -77,6 +85,10 @@ protected:
// Maps from the finer to the coarser level.
PtrList<labelField> restrictAddressing_;
//- The number of (coarse) faces in each level.
// max(faceRestrictAddressing)+1.
labelList nFaces_;
//- Face restriction addressing array.
// Maps from the finer to the coarser level.
// Positive indices map the finer faces which form part of the boundary
@ -85,22 +97,107 @@ protected:
// coarser cells to minus the corresponding coarser cell index minus 1.
PtrList<labelList> faceRestrictAddressing_;
//- Face flip: for faces mapped to internal faces stores whether
// the face is reversed or not. This is used to avoid having
// to access the coarse mesh at all when mapping
PtrList<boolList> faceFlipMap_;
//- The number of (coarse) patch faces in each level.
// max(patchFaceRestrictAddressing_)+1.
PtrList<labelList> nPatchFaces_;
//- Patch-local face restriction addressing array.
// Maps from the finer to the coarser level. Always positive.
// Extracted from GAMGInterfaces after agglomeration.
PtrList<labelListList> patchFaceRestrictAddressing_;
//- Hierarchy of mesh addressing
PtrList<lduPrimitiveMesh> meshLevels_;
//- Hierarchy interfaces.
// Warning: Needs to be deleted explicitly.
PtrList<lduInterfacePtrsList> interfaceLevels_;
// Processor agglomeration
//- Per level, per processor the processor it agglomerates into
mutable PtrList<labelList> procAgglomMap_;
//- Per level the set of processors to agglomerate. Element 0 is
// the 'master' of the cluster.
mutable PtrList<labelList> agglomProcIDs_;
//- Communicator for given level
mutable labelList procCommunicator_;
//- Mapping from processor to procMeshLevel cells
mutable PtrList<labelList> procCellOffsets_;
//- Mapping from processor to procMeshLevel face
mutable PtrList<labelListList> procFaceMap_;
//- Mapping from processor to procMeshLevel boundary
mutable PtrList<labelListList> procBoundaryMap_;
//- Mapping from processor to procMeshLevel boundary face
mutable PtrList<labelListListList> procBoundaryFaceMap_;
// Protected Member Functions
//- Assemble coarse mesh addressing
void agglomerateLduAddressing(const label fineLevelIndex);
//- Combine a level with the previous one
void combineLevels(const label curLevel);
//- Shrink the number of levels to that specified
void compactLevels(const label nCreatedLevels);
//- Check the need for further agglomeration
bool continueAgglomerating(const label nCoarseCells) const;
//- Gather value from all procIDs onto procIDs[0]
template<class Type>
static void gatherList
(
const label comm,
const labelList& procIDs,
const Type& myVal,
List<Type>& allVals,
const int tag = Pstream::msgType()
);
void clearLevel(const label leveli);
// Processor agglomeration
//- Collect and combine processor meshes into allMesh.
// - allMeshComm : communicator for combined mesh.
// - procAgglomMap : per processor the new agglomerated processor
// (rank in allMeshComm!). Global information.
// - procIDs : local information: same for all in
// agglomerated processor.
void procAgglomerateLduAddressing
(
const label comm,
const labelList& procAgglomMap,
const labelList& procIDs,
const label allMeshComm,
const label levelIndex
);
//- Collect and combine basic restriction addressing:
// nCells_
// restrictAddressing_
void procAgglomerateRestrictAddressing
(
const label comm,
const labelList& procIDs,
const label levelIndex
);
private:
// Private Member Functions
@ -113,6 +210,9 @@ protected:
public:
//- Declare friendship with GAMGProcAgglomeration
friend class GAMGProcAgglomeration;
//- Runtime type information
TypeName("GAMGAgglomeration");
@ -152,6 +252,27 @@ public:
)
);
//- Runtime selection table for matrix or mixed geometric/matrix
// agglomerators
declareRunTimeSelectionTable
(
autoPtr,
GAMGAgglomeration,
geometry,
(
const lduMesh& mesh,
const scalarField& cellVolumes,
const vectorField& faceAreas,
const dictionary& controlDict
),
(
mesh,
cellVolumes,
faceAreas,
controlDict
)
);
// Constructors
@ -179,6 +300,15 @@ public:
const dictionary& controlDict
);
//- Return the selected geometric agglomerator
static autoPtr<GAMGAgglomeration> New
(
const lduMesh& mesh,
const scalarField& cellVolumes,
const vectorField& faceAreas,
const dictionary& controlDict
);
//- Destructor
~GAMGAgglomeration();
@ -196,6 +326,9 @@ public:
//- Return LDU mesh of given level
const lduMesh& meshLevel(const label leveli) const;
//- Do we have mesh for given level?
bool hasMeshLevel(const label leveli) const;
//- Return LDU interface addressing of given level
const lduInterfacePtrsList& interfaceLevel
(
@ -214,6 +347,37 @@ public:
return faceRestrictAddressing_[leveli];
}
const labelListList& patchFaceRestrictAddressing(const label leveli)
const
{
return patchFaceRestrictAddressing_[leveli];
}
//- Return face flip map of given level
const boolList& faceFlipMap(const label leveli) const
{
return faceFlipMap_[leveli];
}
//- Return number of coarse cells (before processor agglomeration)
label nCells(const label leveli) const
{
return nCells_[leveli];
}
//- Return number of coarse faces (before processor agglomeration)
label nFaces(const label leveli) const
{
return nFaces_[leveli];
}
//- Return number of coarse patch faces (before processor
// agglomeration)
const labelList& nPatchFaces(const label leveli) const
{
return nPatchFaces_[leveli];
}
// Restriction and prolongation
@ -223,7 +387,8 @@ public:
(
Field<Type>& cf,
const Field<Type>& ff,
const label fineLevelIndex
const label fineLevelIndex,
const bool procAgglom
) const;
//- Restrict (integrate by summation) face field
@ -235,14 +400,80 @@ public:
const label fineLevelIndex
) const;
//- Restrict (integrate by summation) cell field
template<class Type>
void restrictField
(
Field<Type>& cf,
const Field<Type>& ff,
const labelList& fineToCoarse
) const;
//- Prolong (interpolate by injection) cell field
template<class Type>
void prolongField
(
Field<Type>& ff,
const Field<Type>& cf,
const label coarseLevelIndex
const label coarseLevelIndex,
const bool procAgglom
) const;
// Procesor agglomeration. Note that the mesh and agglomeration is
// stored per fineLevel (even though it is the coarse level mesh that
// has been agglomerated). This is just for convenience and consistency
// with GAMGSolver notation.
//- Given fine to coarse processor map determine:
// - for each coarse processor a master (minimum of the fine
// processors)
// - for each coarse processor the set of fine processors
// (element 0 is the master processor)
static void calculateRegionMaster
(
const label comm,
const labelList& procAgglomMap,
labelList& masterProcs,
List<int>& agglomProcIDs
);
//- Whether to agglomerate across processors
bool processorAgglomerate() const
{
return procAgglomeratorPtr_.valid();
}
//- Mapping from processor to agglomerated processor (global, all
// processors have the same information). Note that level is
// the fine, not the coarse, level index. This is to be
// consistent with the way the restriction is stored
const labelList& procAgglomMap(const label fineLeveli) const;
//- Set of processors to agglomerate. Element 0 is the
// master processor. (local, same only on those processors that
// agglomerate)
const labelList& agglomProcIDs(const label fineLeveli) const;
//- Check that level has combined mesh
bool hasProcMesh(const label fineLeveli) const;
//- Communicator for current level or -1
label procCommunicator(const label fineLeveli) const;
//- Mapping from processor to procMesh cells
const labelList& cellOffsets(const label fineLeveli) const;
//- Mapping from processor to procMesh face
const labelListList& faceMap(const label fineLeveli) const;
//- Mapping from processor to procMesh boundary
const labelListList& boundaryMap(const label fineLeveli) const;
//- Mapping from processor to procMesh boundary face
const labelListListList& boundaryFaceMap(const label fineLeveli)
const;
};

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
@ -24,20 +24,85 @@ License
\*---------------------------------------------------------------------------*/
#include "GAMGAgglomeration.H"
#include "mapDistribute.H"
#include "globalIndex.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class Type>
void Foam::GAMGAgglomeration::gatherList
(
const label comm,
const labelList& procIDs,
const Type& myVal,
List<Type>& allVals,
const int tag
)
{
if (Pstream::myProcNo(comm) == procIDs[0])
{
allVals.setSize(procIDs.size());
allVals[0] = myVal;
for (label i = 1; i < procIDs.size(); i++)
{
IPstream fromSlave
(
Pstream::scheduled,
procIDs[i],
0,
tag,
comm
);
fromSlave >> allVals[i];
}
}
else
{
OPstream toMaster
(
Pstream::scheduled,
procIDs[0],
0,
tag,
comm
);
toMaster << myVal;
}
}
template<class Type>
void Foam::GAMGAgglomeration::restrictField
(
Field<Type>& cf,
const Field<Type>& ff,
const label fineLevelIndex
const labelList& fineToCoarse
) const
{
cf = pTraits<Type>::zero;
forAll(ff, i)
{
cf[fineToCoarse[i]] += ff[i];
}
}
template<class Type>
void Foam::GAMGAgglomeration::restrictField
(
Field<Type>& cf,
const Field<Type>& ff,
const label fineLevelIndex,
const bool procAgglom
) const
{
const labelList& fineToCoarse = restrictAddressing_[fineLevelIndex];
if (ff.size() != fineToCoarse.size())
if (!procAgglom && ff.size() != fineToCoarse.size())
{
FatalErrorIn
(
@ -50,11 +115,19 @@ void Foam::GAMGAgglomeration::restrictField
<< abort(FatalError);
}
cf = pTraits<Type>::zero;
restrictField(cf, ff, fineToCoarse);
forAll(ff, i)
label coarseLevelIndex = fineLevelIndex+1;
if (procAgglom && hasProcMesh(coarseLevelIndex))
{
cf[fineToCoarse[i]] += ff[i];
label fineComm = UPstream::parent(procCommunicator_[coarseLevelIndex]);
const List<int>& procIDs = agglomProcIDs(coarseLevelIndex);
const labelList& offsets = cellOffsets(coarseLevelIndex);
const globalIndex cellOffsetter(offsets);
cellOffsetter.gather(fineComm, procIDs, cf);
}
}
@ -69,6 +142,19 @@ void Foam::GAMGAgglomeration::restrictFaceField
{
const labelList& fineToCoarse = faceRestrictAddressing_[fineLevelIndex];
if (ff.size() != fineToCoarse.size())
{
FatalErrorIn
(
"void GAMGAgglomeration::restrictFaceField"
"(Field<Type>& cf, const Field<Type>& ff, "
"const label fineLevelIndex) const"
) << "field does not correspond to level " << fineLevelIndex
<< " sizes: field = " << ff.size()
<< " level = " << fineToCoarse.size()
<< abort(FatalError);
}
cf = pTraits<Type>::zero;
forAll(fineToCoarse, ffacei)
@ -88,14 +174,49 @@ void Foam::GAMGAgglomeration::prolongField
(
Field<Type>& ff,
const Field<Type>& cf,
const label coarseLevelIndex
const label levelIndex,
const bool procAgglom
) const
{
const labelList& fineToCoarse = restrictAddressing_[coarseLevelIndex];
const labelList& fineToCoarse = restrictAddressing_[levelIndex];
forAll(fineToCoarse, i)
label coarseLevelIndex = levelIndex+1;
if (procAgglom && hasProcMesh(coarseLevelIndex))
{
ff[i] = cf[fineToCoarse[i]];
label coarseComm = UPstream::parent
(
procCommunicator_[coarseLevelIndex]
);
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
(
coarseComm,
procIDs,
cf,
allCf
);
forAll(fineToCoarse, i)
{
ff[i] = allCf[fineToCoarse[i]];
}
}
else
{
forAll(fineToCoarse, i)
{
ff[i] = cf[fineToCoarse[i]];
}
}
}

View File

@ -0,0 +1,82 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 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/>.
\*---------------------------------------------------------------------------*/
#include "dummyAgglomeration.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(dummyAgglomeration, 0);
addToRunTimeSelectionTable
(
GAMGAgglomeration,
dummyAgglomeration,
lduMesh
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::dummyAgglomeration::dummyAgglomeration
(
const lduMesh& mesh,
const dictionary& controlDict
)
:
GAMGAgglomeration(mesh, controlDict),
nLevels_(readLabel(controlDict.lookup("nLevels")))
{
// Get the finest-level interfaces from the mesh
meshInterfaces_ = mesh.interfaces();
const label nCoarseCells = mesh.lduAddr().size();
for
(
label nCreatedLevels = 0;
nCreatedLevels < nLevels_;
nCreatedLevels++
)
{
nCells_[nCreatedLevels] = nCoarseCells;
restrictAddressing_.set
(
nCreatedLevels,
new labelField(identity(nCoarseCells))
);
agglomerateLduAddressing(nCreatedLevels);
}
// Shrink the storage of the levels to those created
compactLevels(nLevels_);
}
// ************************************************************************* //

View File

@ -0,0 +1,93 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 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::dummyAgglomeration
Description
Agglomerate without combining cells. Used for testing.
SourceFiles
dummyAgglomeration.C
\*---------------------------------------------------------------------------*/
#ifndef dummyAgglomeration_H
#define dummyAgglomeration_H
#include "GAMGAgglomeration.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class dummyAgglomeration Declaration
\*---------------------------------------------------------------------------*/
class dummyAgglomeration
:
public GAMGAgglomeration
{
// Private data
//- Preset number of levels
label nLevels_;
// Private Member Functions
//- Disallow default bitwise copy construct
dummyAgglomeration(const dummyAgglomeration&);
//- Disallow default bitwise assignment
void operator=(const dummyAgglomeration&);
public:
//- Runtime type information
TypeName("dummy");
// Constructors
//- Construct given mesh and controls
dummyAgglomeration
(
const lduMesh& mesh,
const dictionary& controlDict
);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

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
@ -205,11 +205,7 @@ void Foam::pairGAMGAgglomeration::agglomerate
)
{
// Get the finest-level interfaces from the mesh
interfaceLevels_.set
(
0,
new lduInterfacePtrsList(mesh.interfaces())
);
meshInterfaces_ = mesh.interfaces();
// Start geometric agglomeration from the given faceWeights
scalarField* faceWeightsPtr = const_cast<scalarField*>(&faceWeights);

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
@ -30,7 +30,6 @@ Description
SourceFiles
pairGAMGAgglomeration.C
pairGAMGAgglomerate.C
pairGAMGAgglomerationCombineLevels.C
\*---------------------------------------------------------------------------*/
@ -77,8 +76,6 @@ protected:
const scalarField& faceWeights
);
void combineLevels(const label curLevel);
//- Disallow default bitwise copy construct
pairGAMGAgglomeration(const pairGAMGAgglomeration&);

View File

@ -1,98 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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/>.
\*---------------------------------------------------------------------------*/
#include "pairGAMGAgglomeration.H"
#include "lduInterfacePtrsList.H"
#include "GAMGInterface.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::pairGAMGAgglomeration::combineLevels(const label curLevel)
{
label prevLevel = curLevel - 1;
// Set the previous level nCells to the current
nCells_[prevLevel] = nCells_[curLevel];
// Map the restrictAddressing from the coarser level into the previous
// finer level
const labelList& curResAddr = restrictAddressing_[curLevel];
labelList& prevResAddr = restrictAddressing_[prevLevel];
const labelList& curFaceResAddr = faceRestrictAddressing_[curLevel];
labelList& prevFaceResAddr = faceRestrictAddressing_[prevLevel];
forAll(prevFaceResAddr, i)
{
if (prevFaceResAddr[i] >= 0)
{
prevFaceResAddr[i] = curFaceResAddr[prevFaceResAddr[i]];
}
else
{
prevFaceResAddr[i] = -curResAddr[-prevFaceResAddr[i] - 1] - 1;
}
}
// Delete the restrictAddressing for the coarser level
faceRestrictAddressing_.set(curLevel, NULL);
forAll(prevResAddr, i)
{
prevResAddr[i] = curResAddr[prevResAddr[i]];
}
// Delete the restrictAddressing for the coarser level
restrictAddressing_.set(curLevel, NULL);
// Delete the matrix addressing and coefficients from the previous level
// and replace with the corresponding entried from the coarser level
meshLevels_.set(prevLevel, meshLevels_.set(curLevel, NULL));
// Same for the lduInterfaceFields taking care to delete the sub-entries
// held on List<T*>
const lduInterfacePtrsList& curInterLevel = interfaceLevels_[curLevel+1];
lduInterfacePtrsList& prevInterLevel = interfaceLevels_[prevLevel+1];
forAll(prevInterLevel, inti)
{
if (prevInterLevel.set(inti))
{
refCast<GAMGInterface>(const_cast<lduInterface&>
(
prevInterLevel[inti]
)).combine(refCast<const GAMGInterface>(curInterLevel[inti]));
delete curInterLevel(inti);
}
}
interfaceLevels_.set(curLevel+1, NULL);
}
// ************************************************************************* //

View File

@ -0,0 +1,249 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 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/>.
\*---------------------------------------------------------------------------*/
#include "GAMGProcAgglomeration.H"
#include "GAMGAgglomeration.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(GAMGProcAgglomeration, 0);
defineRunTimeSelectionTable(GAMGProcAgglomeration, GAMGAgglomeration);
}
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
void Foam::GAMGProcAgglomeration::printStats
(
Ostream& os,
GAMGAgglomeration& agglom
) const
{
for (label levelI = 0; levelI <= agglom.size(); levelI++)
{
if (agglom.hasMeshLevel(levelI))
{
os << agglom.meshLevel(levelI).info() << endl;
}
else
{
os << "Level " << levelI << " has no fine mesh:" << endl;
}
if
(
levelI < agglom.restrictAddressing_.size()
&& agglom.restrictAddressing_.set(levelI)
)
{
const labelList& cellRestrict =
agglom.restrictAddressing(levelI);
const labelList& faceRestrict =
agglom.faceRestrictAddressing(levelI);
os << "Level " << levelI << " agglomeration:" << nl
<< " nCoarseCells:" << agglom.nCells(levelI) << nl
<< " nCoarseFaces:" << agglom.nFaces(levelI) << nl
<< " cellRestriction:"
<< " size:" << cellRestrict.size()
<< " max:" << max(cellRestrict)
<< nl
<< " faceRestriction:"
<< " size:" << faceRestrict.size()
<< " max:" << max(faceRestrict)
<< nl;
const labelListList& patchFaceRestrict =
agglom.patchFaceRestrictAddressing(levelI);
forAll(patchFaceRestrict, i)
{
if (patchFaceRestrict[i].size())
{
const labelList& faceRestrict =
patchFaceRestrict[i];
os << " " << i
<< " size:" << faceRestrict.size()
<< " max:" << max(faceRestrict)
<< nl;
}
}
}
if
(
levelI < agglom.procCellOffsets_.size()
&& agglom.procCellOffsets_.set(levelI)
)
{
os << " procCellOffsets:" << agglom.procCellOffsets_[levelI]
<< nl
<< " procAgglomMap:" << agglom.procAgglomMap_[levelI]
<< nl
<< " procIDs:" << agglom.agglomProcIDs_[levelI]
<< nl
<< " comm:" << agglom.procCommunicator_[levelI]
<< endl;
}
os << endl;
}
os << endl;
}
bool Foam::GAMGProcAgglomeration::agglomerate
(
const label fineLevelIndex,
const labelList& procAgglomMap,
const labelList& masterProcs,
const List<int>& agglomProcIDs,
const label procAgglomComm
)
{
const lduMesh& levelMesh = agglom_.meshLevels_[fineLevelIndex];
label levelComm = levelMesh.comm();
if (Pstream::myProcNo(levelComm) != -1)
{
// Collect meshes and restrictAddressing onto master
// Overwrites the fine mesh (meshLevels_[index-1]) and addressing
// from fine mesh to coarse mesh (restrictAddressing_[index]).
agglom_.procAgglomerateLduAddressing
(
levelComm,
procAgglomMap,
agglomProcIDs,
procAgglomComm,
fineLevelIndex //fine level index
);
// Combine restrict addressing only onto master
for
(
label levelI = fineLevelIndex+1;
levelI < agglom_.meshLevels_.size();
levelI++
)
{
agglom_.procAgglomerateRestrictAddressing
(
levelComm,
agglomProcIDs,
levelI
);
}
if (Pstream::myProcNo(levelComm) == agglomProcIDs[0])
{
// On master. Recreate coarse meshes from restrict addressing
for
(
label levelI = fineLevelIndex;
levelI < agglom_.meshLevels_.size();
levelI++
)
{
agglom_.agglomerateLduAddressing(levelI);
}
}
else
{
// Agglomerated away. Clear mesh storage.
for
(
label levelI = fineLevelIndex+1;
levelI <= agglom_.size();
levelI++
)
{
agglom_.clearLevel(levelI);
}
}
}
// Should check!
return true;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::GAMGProcAgglomeration::GAMGProcAgglomeration
(
GAMGAgglomeration& agglom,
const dictionary& controlDict
)
:
agglom_(agglom)
{}
Foam::autoPtr<Foam::GAMGProcAgglomeration> Foam::GAMGProcAgglomeration::New
(
const word& type,
GAMGAgglomeration& agglom,
const dictionary& controlDict
)
{
if (debug)
{
Info<< "GAMGProcAgglomeration::New(const word&, GAMGAgglomeration&"
", const dictionary&) : "
"constructing GAMGProcAgglomeration"
<< endl;
}
GAMGAgglomerationConstructorTable::iterator cstrIter =
GAMGAgglomerationConstructorTablePtr_->find(type);
if (cstrIter == GAMGAgglomerationConstructorTablePtr_->end())
{
FatalErrorIn
(
"GAMGProcAgglomeration::New(const word&, GAMGAgglomeration&"
", const dictionary&) "
) << "Unknown GAMGProcAgglomeration type "
<< type << " for GAMGAgglomeration " << agglom.type() << nl << nl
<< "Valid GAMGProcAgglomeration types are :" << endl
<< GAMGAgglomerationConstructorTablePtr_->sortedToc()
<< exit(FatalError);
}
return autoPtr<GAMGProcAgglomeration>(cstrIter()(agglom, controlDict));
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::GAMGProcAgglomeration::~GAMGProcAgglomeration()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// ************************************************************************* //

View File

@ -0,0 +1,156 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 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::GAMGProcAgglomeration
Description
Processor agglomeration of GAMGAgglomerations.
SourceFiles
GAMGProcAgglomeration.C
\*---------------------------------------------------------------------------*/
#ifndef GAMGProcAgglomeration_H
#define GAMGProcAgglomeration_H
#include "runTimeSelectionTables.H"
#include "labelList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
class GAMGAgglomeration;
/*---------------------------------------------------------------------------*\
Class GAMGProcAgglomeration Declaration
\*---------------------------------------------------------------------------*/
class GAMGProcAgglomeration
{
protected:
// Protected data
//- Reference to agglomeration
GAMGAgglomeration& agglom_;
// Protected Member Functions
//- Debug: write agglomeration info
void printStats(Ostream& os, GAMGAgglomeration& agglom) const;
//- Agglomerate a level. Return true if anything has changed
bool agglomerate
(
const label fineLevelIndex,
const labelList& procAgglomMap,
const labelList& masterProcs,
const List<int>& agglomProcIDs,
const label procAgglomComm
);
private:
// Private data
// Private Member Functions
//- Disallow default bitwise copy construct
GAMGProcAgglomeration(const GAMGProcAgglomeration&);
//- Disallow default bitwise assignment
void operator=(const GAMGProcAgglomeration&);
public:
//- Runtime type information
TypeName("GAMGProcAgglomeration");
// Declare run-time constructor selection tables
//- Runtime selection table for pure geometric agglomerators
declareRunTimeSelectionTable
(
autoPtr,
GAMGProcAgglomeration,
GAMGAgglomeration,
(
GAMGAgglomeration& agglom,
const dictionary& controlDict
),
(
agglom,
controlDict
)
);
// Constructors
//- Construct given agglomerator and controls
GAMGProcAgglomeration
(
GAMGAgglomeration& agglom,
const dictionary& controlDict
);
// Selectors
//- Return the selected agglomerator
static autoPtr<GAMGProcAgglomeration> New
(
const word& type,
GAMGAgglomeration& agglom,
const dictionary& controlDict
);
//- Destructor
virtual ~GAMGProcAgglomeration();
// Member Functions
//- Modify agglomeration. Return true if modified
virtual bool agglomerate() = 0;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,168 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 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/>.
\*---------------------------------------------------------------------------*/
#include "eagerGAMGProcAgglomeration.H"
#include "addToRunTimeSelectionTable.H"
#include "GAMGAgglomeration.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(eagerGAMGProcAgglomeration, 0);
addToRunTimeSelectionTable
(
GAMGProcAgglomeration,
eagerGAMGProcAgglomeration,
GAMGAgglomeration
);
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::eagerGAMGProcAgglomeration::eagerGAMGProcAgglomeration
(
GAMGAgglomeration& agglom,
const dictionary& controlDict
)
:
GAMGProcAgglomeration(agglom, controlDict),
mergeLevels_(controlDict.lookupOrDefault<label>("mergeLevels", 1))
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::eagerGAMGProcAgglomeration::
~eagerGAMGProcAgglomeration()
{
forAllReverse(comms_, i)
{
if (comms_[i] != -1)
{
UPstream::freeCommunicator(comms_[i]);
}
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool Foam::eagerGAMGProcAgglomeration::agglomerate()
{
if (debug)
{
Pout<< nl << "Starting mesh overview" << endl;
printStats(Pout, agglom_);
}
if (agglom_.size() >= 1)
{
// Agglomerate one but last level (since also agglomerating
// restrictAddressing)
for
(
label fineLevelIndex = 2;
fineLevelIndex < agglom_.size();
fineLevelIndex++
)
{
if (agglom_.hasMeshLevel(fineLevelIndex))
{
// Get the fine mesh
const lduMesh& levelMesh = agglom_.meshLevel(fineLevelIndex);
label levelComm = levelMesh.comm();
label nProcs = UPstream::nProcs(levelComm);
if (nProcs > 1)
{
// Processor restriction map: per processor the coarse
// processor
labelList procAgglomMap(nProcs);
forAll(procAgglomMap, procI)
{
procAgglomMap[procI] = procI/(1<<mergeLevels_);
}
// Master processor
labelList masterProcs;
// Local processors that agglomerate. agglomProcIDs[0]
// is in masterProc.
List<int> agglomProcIDs;
GAMGAgglomeration::calculateRegionMaster
(
levelComm,
procAgglomMap,
masterProcs,
agglomProcIDs
);
// Allocate a communicator for the processor-agglomerated
// matrix
comms_.append
(
UPstream::allocateCommunicator
(
levelComm,
masterProcs
)
);
// Use procesor agglomeration maps to do the actual
// collecting.
if (Pstream::myProcNo(levelComm) != -1)
{
GAMGProcAgglomeration::agglomerate
(
fineLevelIndex,
procAgglomMap,
masterProcs,
agglomProcIDs,
comms_.last()
);
}
}
}
}
}
// Print a bit
if (debug)
{
Pout<< nl << "Agglomerated mesh overview" << endl;
printStats(Pout, agglom_);
}
return true;
}
// ************************************************************************* //

View File

@ -0,0 +1,113 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 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::eagerGAMGProcAgglomeration
Description
'Eager' processor agglomeration of GAMGAgglomerations: at every
level agglomerates 'mergeLevels' number of processors onto the
minimum processor number.
SourceFiles
eagerGAMGProcAgglomeration.C
\*---------------------------------------------------------------------------*/
#ifndef eagerGAMGProcAgglomeration_H
#define eagerGAMGProcAgglomeration_H
#include "GAMGProcAgglomeration.H"
#include "DynamicList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
class GAMGAgglomeration;
/*---------------------------------------------------------------------------*\
Class eagerGAMGProcAgglomeration Declaration
\*---------------------------------------------------------------------------*/
class eagerGAMGProcAgglomeration
:
public GAMGProcAgglomeration
{
// Private data
//- Agglpmeration level
const label mergeLevels_;
DynamicList<label> comms_;
// Private Member Functions
//- Disallow default bitwise copy construct
eagerGAMGProcAgglomeration
(
const eagerGAMGProcAgglomeration&
);
//- Disallow default bitwise assignment
void operator=(const eagerGAMGProcAgglomeration&);
public:
//- Runtime type information
TypeName("eager");
// Constructors
//- Construct given agglomerator and controls
eagerGAMGProcAgglomeration
(
GAMGAgglomeration& agglom,
const dictionary& controlDict
);
//- Destructor
virtual ~eagerGAMGProcAgglomeration();
// Member Functions
//- Modify agglomeration. Return true if modified
virtual bool agglomerate();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,219 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 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/>.
\*---------------------------------------------------------------------------*/
#include "manualGAMGProcAgglomeration.H"
#include "addToRunTimeSelectionTable.H"
#include "GAMGAgglomeration.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(manualGAMGProcAgglomeration, 0);
addToRunTimeSelectionTable
(
GAMGProcAgglomeration,
manualGAMGProcAgglomeration,
GAMGAgglomeration
);
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::manualGAMGProcAgglomeration::manualGAMGProcAgglomeration
(
GAMGAgglomeration& agglom,
const dictionary& controlDict
)
:
GAMGProcAgglomeration(agglom, controlDict),
procAgglomMaps_(controlDict.lookup("processorAgglomeration"))
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::manualGAMGProcAgglomeration::
~manualGAMGProcAgglomeration()
{
forAllReverse(comms_, i)
{
if (comms_[i] != -1)
{
UPstream::freeCommunicator(comms_[i]);
}
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool Foam::manualGAMGProcAgglomeration::agglomerate()
{
if (debug)
{
Pout<< nl << "Starting mesh overview" << endl;
printStats(Pout, agglom_);
}
if (agglom_.size() >= 1)
{
forAll(procAgglomMaps_, i)
{
const label fineLevelIndex = procAgglomMaps_[i].first();
if (fineLevelIndex >= agglom_.size())
{
WarningIn("manualGAMGProcAgglomeration::agglomerate()")
<< "Ignoring specification for level " << fineLevelIndex
<< " since outside agglomeration." << endl;
continue;
}
if (agglom_.hasMeshLevel(fineLevelIndex))
{
// Get the fine mesh
const lduMesh& levelMesh = agglom_.meshLevel(fineLevelIndex);
label nProcs = UPstream::nProcs(levelMesh.comm());
if (nProcs > 1)
{
// My processor id
const label myProcID = Pstream::myProcNo(levelMesh.comm());
const List<clusterAndMaster>& clusters =
procAgglomMaps_[i].second();
// Coarse to fine master processor
labelList coarseToMaster(clusters.size());
// Fine to coarse map
labelList procAgglomMap(nProcs, -1);
// Cluster for my processor (with master index first)
labelList agglomProcIDs;
forAll(clusters, coarseI)
{
const labelList& cluster = clusters[coarseI].first();
coarseToMaster[coarseI] = clusters[coarseI].second();
forAll(cluster, i)
{
procAgglomMap[cluster[i]] = coarseI;
}
label masterIndex = findIndex
(
cluster,
coarseToMaster[coarseI]
);
if (masterIndex == -1)
{
FatalErrorIn
(
"manualGAMGProcAgglomeration::agglomerate()"
) << "At level " << fineLevelIndex
<< " the master processor "
<< coarseToMaster[coarseI]
<< " is not in the cluster "
<< cluster
<< exit(FatalError);
}
if (findIndex(cluster, myProcID) != -1)
{
// This is my cluster. Make sure master index is
// first
agglomProcIDs = cluster;
Swap(agglomProcIDs[0], agglomProcIDs[masterIndex]);
}
}
// Check that we've done all processors
if (findIndex(procAgglomMap, -1) != -1)
{
FatalErrorIn
(
"manualGAMGProcAgglomeration::agglomerate()"
) << "At level " << fineLevelIndex
<< " processor "
<< findIndex(procAgglomMap, -1)
<< " is not in any cluster"
<< exit(FatalError);
}
// Allocate a communicator for the processor-agglomerated
// matrix
comms_.append
(
UPstream::allocateCommunicator
(
levelMesh.comm(),
coarseToMaster
)
);
// Use procesor agglomeration maps to do the actual
// collecting
if (Pstream::myProcNo(levelMesh.comm()) != -1)
{
GAMGProcAgglomeration::agglomerate
(
fineLevelIndex,
procAgglomMap,
coarseToMaster,
agglomProcIDs,
comms_.last()
);
}
}
}
}
// Print a bit
if (debug)
{
Pout<< nl << "Agglomerated mesh overview" << endl;
printStats(Pout, agglom_);
}
}
return true;
}
// ************************************************************************* //

View File

@ -0,0 +1,136 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 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::manualGAMGProcAgglomeration
Description
Manual processor agglomeration of GAMGAgglomerations.
In the GAMG control dictionary:
processorAgglomerator manual;
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)
)
)
(
6 //at level6
(
((0 1) 0) //coarse 0 from 0,1 (and moved onto 0)
)
)
);
SourceFiles
manualGAMGProcAgglomeration.C
\*---------------------------------------------------------------------------*/
#ifndef manualGAMGProcAgglomeration_H
#define manualGAMGProcAgglomeration_H
#include "GAMGProcAgglomeration.H"
#include "DynamicList.H"
#include "Tuple2.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
class GAMGAgglomeration;
/*---------------------------------------------------------------------------*\
Class manualGAMGProcAgglomeration Declaration
\*---------------------------------------------------------------------------*/
class manualGAMGProcAgglomeration
:
public GAMGProcAgglomeration
{
// Private data
typedef Tuple2<labelList, label> clusterAndMaster;
//- Per level the agglomeration map
const List<Tuple2<label, List<clusterAndMaster> > > procAgglomMaps_;
//- Any allocated communicators
DynamicList<label> comms_;
// Private Member Functions
//- Disallow default bitwise copy construct
manualGAMGProcAgglomeration
(
const manualGAMGProcAgglomeration&
);
//- Disallow default bitwise assignment
void operator=(const manualGAMGProcAgglomeration&);
public:
//- Runtime type information
TypeName("manual");
// Constructors
//- Construct given agglomerator and controls
manualGAMGProcAgglomeration
(
GAMGAgglomeration& agglom,
const dictionary& controlDict
);
//- Destructor
virtual ~manualGAMGProcAgglomeration();
// Member Functions
//- Modify agglomeration. Return true if modified
virtual bool agglomerate();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,153 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 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/>.
\*---------------------------------------------------------------------------*/
#include "masterCoarsestGAMGProcAgglomeration.H"
#include "addToRunTimeSelectionTable.H"
#include "GAMGAgglomeration.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(masterCoarsestGAMGProcAgglomeration, 0);
addToRunTimeSelectionTable
(
GAMGProcAgglomeration,
masterCoarsestGAMGProcAgglomeration,
GAMGAgglomeration
);
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::masterCoarsestGAMGProcAgglomeration::masterCoarsestGAMGProcAgglomeration
(
GAMGAgglomeration& agglom,
const dictionary& controlDict
)
:
GAMGProcAgglomeration(agglom, controlDict)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::masterCoarsestGAMGProcAgglomeration::
~masterCoarsestGAMGProcAgglomeration()
{
forAllReverse(comms_, i)
{
if (comms_[i] != -1)
{
UPstream::freeCommunicator(comms_[i]);
}
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool Foam::masterCoarsestGAMGProcAgglomeration::agglomerate()
{
if (debug)
{
Pout<< nl << "Starting mesh overview" << endl;
printStats(Pout, agglom_);
}
if (agglom_.size() >= 1)
{
// Agglomerate one but last level (since also agglomerating
// restrictAddressing)
label fineLevelIndex = agglom_.size()-1;
if (agglom_.hasMeshLevel(fineLevelIndex))
{
// Get the fine mesh
const lduMesh& levelMesh = agglom_.meshLevel(fineLevelIndex);
label levelComm = levelMesh.comm();
label nProcs = UPstream::nProcs(levelComm);
if (nProcs > 1)
{
// Processor restriction map: per processor the coarse processor
labelList procAgglomMap(nProcs, 0);
// Master processor
labelList masterProcs;
// Local processors that agglomerate. agglomProcIDs[0] is in
// masterProc.
List<int> agglomProcIDs;
GAMGAgglomeration::calculateRegionMaster
(
levelComm,
procAgglomMap,
masterProcs,
agglomProcIDs
);
// Allocate a communicator for the processor-agglomerated matrix
comms_.append
(
UPstream::allocateCommunicator
(
levelComm,
masterProcs
)
);
// Use procesor agglomeration maps to do the actual collecting.
if (Pstream::myProcNo(levelComm) != -1)
{
GAMGProcAgglomeration::agglomerate
(
fineLevelIndex,
procAgglomMap,
masterProcs,
agglomProcIDs,
comms_.last()
);
}
}
}
}
// Print a bit
if (debug)
{
Pout<< nl << "Agglomerated mesh overview" << endl;
printStats(Pout, agglom_);
}
return true;
}
// ************************************************************************* //

View File

@ -0,0 +1,108 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 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::masterCoarsestGAMGProcAgglomeration
Description
Processor agglomeration of GAMGAgglomerations.
SourceFiles
masterCoarsestGAMGProcAgglomeration.C
\*---------------------------------------------------------------------------*/
#ifndef masterCoarsestGAMGProcAgglomeration_H
#define masterCoarsestGAMGProcAgglomeration_H
#include "GAMGProcAgglomeration.H"
#include "DynamicList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
class GAMGAgglomeration;
/*---------------------------------------------------------------------------*\
Class masterCoarsestGAMGProcAgglomeration Declaration
\*---------------------------------------------------------------------------*/
class masterCoarsestGAMGProcAgglomeration
:
public GAMGProcAgglomeration
{
// Private data
DynamicList<label> comms_;
// Private Member Functions
//- Disallow default bitwise copy construct
masterCoarsestGAMGProcAgglomeration
(
const masterCoarsestGAMGProcAgglomeration&
);
//- Disallow default bitwise assignment
void operator=(const masterCoarsestGAMGProcAgglomeration&);
public:
//- Runtime type information
TypeName("masterCoarsest");
// Constructors
//- Construct given agglomerator and controls
masterCoarsestGAMGProcAgglomeration
(
GAMGAgglomeration& agglom,
const dictionary& controlDict
);
//- Destructor
virtual ~masterCoarsestGAMGProcAgglomeration();
// Member Functions
//- Modify agglomeration. Return true if modified
virtual bool agglomerate();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,70 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 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/>.
\*---------------------------------------------------------------------------*/
#include "noneGAMGProcAgglomeration.H"
#include "addToRunTimeSelectionTable.H"
#include "GAMGAgglomeration.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(noneGAMGProcAgglomeration, 0);
addNamedToRunTimeSelectionTable
(
GAMGProcAgglomeration,
noneGAMGProcAgglomeration,
GAMGAgglomeration,
none
);
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::noneGAMGProcAgglomeration::noneGAMGProcAgglomeration
(
GAMGAgglomeration& agglom,
const dictionary& controlDict
)
:
GAMGProcAgglomeration(agglom, controlDict)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::noneGAMGProcAgglomeration::~noneGAMGProcAgglomeration()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// ************************************************************************* //

View File

@ -0,0 +1,107 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 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::noneGAMGProcAgglomeration
Description
Processor agglomeration of GAMGAgglomerations.
SourceFiles
noneGAMGProcAgglomeration.C
\*---------------------------------------------------------------------------*/
#ifndef noneGAMGProcAgglomeration_H
#define noneGAMGProcAgglomeration_H
#include "GAMGProcAgglomeration.H"
#include "DynamicList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
class GAMGAgglomeration;
/*---------------------------------------------------------------------------*\
Class noneGAMGProcAgglomeration Declaration
\*---------------------------------------------------------------------------*/
class noneGAMGProcAgglomeration
:
public GAMGProcAgglomeration
{
// Private Member Functions
//- Disallow default bitwise copy construct
noneGAMGProcAgglomeration
(
const noneGAMGProcAgglomeration&
);
//- Disallow default bitwise assignment
void operator=(const noneGAMGProcAgglomeration&);
public:
//- Runtime type information
TypeName("noneGAMGProcAgglomeration");
// Constructors
//- Construct given agglomerator and controls
noneGAMGProcAgglomeration
(
GAMGAgglomeration& agglom,
const dictionary& controlDict
);
//- Destructor
virtual ~noneGAMGProcAgglomeration();
// Member Functions
//- Modify agglomeration. Return true if modified
virtual bool agglomerate()
{
return false;
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -24,6 +24,7 @@ License
\*---------------------------------------------------------------------------*/
#include "GAMGSolver.H"
#include "GAMGInterface.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -77,32 +78,200 @@ Foam::GAMGSolver::GAMGSolver
agglomeration_(GAMGAgglomeration::New(matrix_, controlDict_)),
matrixLevels_(agglomeration_.size()),
primitiveInterfaceLevels_(agglomeration_.size()),
interfaceLevels_(agglomeration_.size()),
interfaceLevelsBouCoeffs_(agglomeration_.size()),
interfaceLevelsIntCoeffs_(agglomeration_.size())
{
readControls();
forAll(agglomeration_, fineLevelIndex)
if (agglomeration_.processorAgglomerate())
{
agglomerateMatrix(fineLevelIndex);
forAll(agglomeration_, fineLevelIndex)
{
if (agglomeration_.hasMeshLevel(fineLevelIndex))
{
if
(
(fineLevelIndex+1) < agglomeration_.size()
&& agglomeration_.hasProcMesh(fineLevelIndex+1)
)
{
// Construct matrix without referencing the coarse mesh so
// construct a dummy mesh instead. This will get overwritten
// by the call to procAgglomerateMatrix so is only to get
// it through agglomerateMatrix
const lduInterfacePtrsList& fineMeshInterfaces =
agglomeration_.interfaceLevel(fineLevelIndex);
PtrList<GAMGInterface> dummyPrimMeshInterfaces
(
fineMeshInterfaces.size()
);
lduInterfacePtrsList dummyMeshInterfaces
(
dummyPrimMeshInterfaces.size()
);
forAll(fineMeshInterfaces, intI)
{
if (fineMeshInterfaces.set(intI))
{
OStringStream os;
refCast<const GAMGInterface>
(
fineMeshInterfaces[intI]
).write(os);
IStringStream is(os.str());
dummyPrimMeshInterfaces.set
(
intI,
GAMGInterface::New
(
fineMeshInterfaces[intI].type(),
intI,
dummyMeshInterfaces,
is
)
);
}
}
forAll(dummyPrimMeshInterfaces, intI)
{
if (dummyPrimMeshInterfaces.set(intI))
{
dummyMeshInterfaces.set
(
intI,
&dummyPrimMeshInterfaces[intI]
);
}
}
// So:
// - pass in incorrect mesh (= fine mesh instead of coarse)
// - pass in dummy interfaces
agglomerateMatrix
(
fineLevelIndex,
agglomeration_.meshLevel(fineLevelIndex),
dummyMeshInterfaces
);
const labelList& procAgglomMap =
agglomeration_.procAgglomMap(fineLevelIndex+1);
const List<int>& procIDs =
agglomeration_.agglomProcIDs(fineLevelIndex+1);
procAgglomerateMatrix
(
procAgglomMap,
procIDs,
fineLevelIndex
);
}
else
{
agglomerateMatrix
(
fineLevelIndex,
agglomeration_.meshLevel(fineLevelIndex + 1),
agglomeration_.interfaceLevel(fineLevelIndex + 1)
);
}
}
else
{
// No mesh. Not involved in calculation anymore
}
}
}
else
{
forAll(agglomeration_, fineLevelIndex)
{
// Agglomerate on to coarse level mesh
agglomerateMatrix
(
fineLevelIndex,
agglomeration_.meshLevel(fineLevelIndex + 1),
agglomeration_.interfaceLevel(fineLevelIndex + 1)
);
}
}
if (debug)
{
for
(
label fineLevelIndex = 0;
fineLevelIndex <= matrixLevels_.size();
fineLevelIndex++
)
{
if (fineLevelIndex == 0 || matrixLevels_.set(fineLevelIndex-1))
{
const lduMatrix& matrix = matrixLevel(fineLevelIndex);
const lduInterfaceFieldPtrsList& interfaces =
interfaceLevel(fineLevelIndex);
Pout<< "level:" << fineLevelIndex << nl
<< " nCells:" << matrix.diag().size() << nl
<< " nFaces:" << matrix.lower().size() << nl
<< " nInterfaces:" << interfaces.size()
<< endl;
forAll(interfaces, i)
{
if (interfaces.set(i))
{
Pout<< " " << i
<< "\ttype:" << interfaces[i].type()
<< endl;
}
}
}
else
{
Pout<< "level:" << fineLevelIndex << " : no matrix" << endl;
}
}
Pout<< endl;
}
if (matrixLevels_.size())
{
const label coarsestLevel = matrixLevels_.size() - 1;
if (directSolveCoarsest_)
{
coarsestLUMatrixPtr_.set
(
new LUscalarMatrix
const label coarsestLevel = matrixLevels_.size() - 1;
if (matrixLevels_.set(coarsestLevel))
{
const lduMesh& coarsestMesh =
matrixLevels_[coarsestLevel].mesh();
label coarseComm = coarsestMesh.comm();
label oldWarn = UPstream::warnComm;
UPstream::warnComm = coarseComm;
coarsestLUMatrixPtr_.set
(
matrixLevels_[coarsestLevel],
interfaceLevelsBouCoeffs_[coarsestLevel],
interfaceLevels_[coarsestLevel]
)
);
new LUscalarMatrix
(
matrixLevels_[coarsestLevel],
interfaceLevelsBouCoeffs_[coarsestLevel],
interfaceLevels_[coarsestLevel]
)
);
UPstream::warnComm = oldWarn;
}
}
}
else
@ -131,20 +300,6 @@ Foam::GAMGSolver::GAMGSolver
Foam::GAMGSolver::~GAMGSolver()
{
// Clear the the lists of pointers to the interfaces
forAll(interfaceLevels_, leveli)
{
lduInterfaceFieldPtrsList& curLevel = interfaceLevels_[leveli];
forAll(curLevel, i)
{
if (curLevel.set(i))
{
delete curLevel(i);
}
}
}
if (!cacheAgglomeration_)
{
delete &agglomeration_;
@ -178,6 +333,23 @@ void Foam::GAMGSolver::readControls()
controlDict_.readIfPresent("interpolateCorrection", interpolateCorrection_);
controlDict_.readIfPresent("scaleCorrection", scaleCorrection_);
controlDict_.readIfPresent("directSolveCoarsest", directSolveCoarsest_);
if (debug)
{
Pout<< "GAMGSolver settings :"
<< " cacheAgglomeration:" << cacheAgglomeration_
<< " nPreSweeps:" << nPreSweeps_
<< " preSweepsLevelMultiplier:" << preSweepsLevelMultiplier_
<< " maxPreSweeps:" << maxPreSweeps_
<< " nPostSweeps:" << nPostSweeps_
<< " postSweepsLevelMultiplier:" << postSweepsLevelMultiplier_
<< " maxPostSweeps:" << maxPostSweeps_
<< " nFinestSweeps:" << nFinestSweeps_
<< " interpolateCorrection:" << interpolateCorrection_
<< " scaleCorrection:" << scaleCorrection_
<< " directSolveCoarsest:" << directSolveCoarsest_
<< endl;
}
}

View File

@ -43,9 +43,7 @@ Description
SourceFiles
GAMGSolver.C
GAMGSolverCalcAgglomeration.C
GAMGSolverMakeCoarseMatrix.C
GAMGSolverOperations.C
GAMGSolverAgglomerateMatrix.C
GAMGSolverInterpolate.C
GAMGSolverScale.C
GAMGSolverSolve.C
@ -118,7 +116,9 @@ class GAMGSolver
PtrList<lduMatrix> matrixLevels_;
//- Hierarchy of interfaces.
// Warning: Needs to be deleted explicitly.
PtrList<PtrList<lduInterfaceField> > primitiveInterfaceLevels_;
//- Hierarchy of interfaces in lduInterfaceFieldPtrs form
PtrList<lduInterfaceFieldPtrsList> interfaceLevels_;
//- Hierarchy of interface boundary coefficients
@ -157,8 +157,70 @@ class GAMGSolver
const label i
) const;
//- Agglomerate coarse matrix
void agglomerateMatrix(const label fineLevelIndex);
//- Agglomerate coarse matrix. Supply mesh to use - so we can
// construct temporary matrix on the fine mesh (instead of the coarse
// mesh)
void agglomerateMatrix
(
const label fineLevelIndex,
const lduMesh& coarseMesh,
const lduInterfacePtrsList& coarseMeshInterfaces
);
//- Agglomerate coarse interface coefficients
void agglomerateInterfaceCoefficients
(
const label fineLevelIndex,
const lduInterfacePtrsList& coarseMeshInterfaces,
PtrList<lduInterfaceField>& coarsePrimInterfaces,
lduInterfaceFieldPtrsList& coarseInterfaces,
FieldField<Field, scalar>& coarseInterfaceBouCoeffs,
FieldField<Field, scalar>& coarseInterfaceIntCoeffs
) const;
//- Collect matrices from other processors
void gatherMatrices
(
const labelList& procIDs,
const lduMesh& dummyMesh,
const label meshComm,
const lduMatrix& mat,
const FieldField<Field, scalar>& interfaceBouCoeffs,
const FieldField<Field, scalar>& interfaceIntCoeffs,
const lduInterfaceFieldPtrsList& interfaces,
PtrList<lduMatrix>& otherMats,
PtrList<FieldField<Field, scalar> >& otherBouCoeffs,
PtrList<FieldField<Field, scalar> >& otherIntCoeffs,
List<boolList>& otherTransforms,
List<List<int> >& otherRanks
) const;
//- Agglomerate processor matrices
void procAgglomerateMatrix
(
// Agglomeration information
const labelList& procAgglomMap,
const List<int>& agglomProcIDs,
const label levelI,
// Resulting matrix
autoPtr<lduMatrix>& allMatrixPtr,
FieldField<Field, scalar>& allInterfaceBouCoeffs,
FieldField<Field, scalar>& allInterfaceIntCoeffs,
PtrList<lduInterfaceField>& allPrimitiveInterfaces,
lduInterfaceFieldPtrsList& allInterfaces
) const;
//- Agglomerate processor matrices
void procAgglomerateMatrix
(
const labelList& procAgglomMap,
const List<int>& agglomProcIDs,
const label levelI
);
//- Interpolate the correction after injected prolongation
void interpolate
@ -193,7 +255,9 @@ class GAMGSolver
(
PtrList<scalarField>& coarseCorrFields,
PtrList<scalarField>& coarseSources,
PtrList<lduMatrix::smoother>& smoothers
PtrList<lduMatrix::smoother>& smoothers,
scalarField& scratch1,
scalarField& scratch2
) const;
@ -206,6 +270,10 @@ class GAMGSolver
scalarField& Apsi,
scalarField& finestCorrection,
scalarField& finestResidual,
scalarField& scratch1,
scalarField& scratch2,
PtrList<scalarField>& coarseCorrFields,
PtrList<scalarField>& coarseSources,
const direction cmpt=0

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
@ -25,30 +25,185 @@ License
#include "GAMGSolver.H"
#include "GAMGInterfaceField.H"
#include "processorLduInterfaceField.H"
#include "processorGAMGInterfaceField.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::GAMGSolver::agglomerateMatrix(const label fineLevelIndex)
void Foam::GAMGSolver::agglomerateMatrix
(
const label fineLevelIndex,
const lduMesh& coarseMesh,
const lduInterfacePtrsList& coarseMeshInterfaces
)
{
// Get fine matrix
const lduMatrix& fineMatrix = matrixLevel(fineLevelIndex);
// Set the coarse level matrix
matrixLevels_.set
(
fineLevelIndex,
new lduMatrix(agglomeration_.meshLevel(fineLevelIndex + 1))
);
lduMatrix& coarseMatrix = matrixLevels_[fineLevelIndex];
if (UPstream::myProcNo(fineMatrix.mesh().comm()) != -1)
{
const label nCoarseFaces = agglomeration_.nFaces(fineLevelIndex);
const label nCoarseCells = agglomeration_.nCells(fineLevelIndex);
// Get face restriction map for current level
const labelList& faceRestrictAddr =
agglomeration_.faceRestrictAddressing(fineLevelIndex);
// Set the coarse level matrix
matrixLevels_.set
(
fineLevelIndex,
new lduMatrix(coarseMesh)
);
lduMatrix& coarseMatrix = matrixLevels_[fineLevelIndex];
// Coarse matrix diagonal initialised by restricting the finer mesh diagonal
scalarField& coarseDiag = coarseMatrix.diag();
agglomeration_.restrictField(coarseDiag, fineMatrix.diag(), fineLevelIndex);
// Coarse matrix diagonal initialised by restricting the finer mesh
// diagonal. Note that we size with the cached coarse nCells and not
// the actual coarseMesh size since this might be dummy when processor
// agglomerating.
scalarField& coarseDiag = coarseMatrix.diag(nCoarseCells);
agglomeration_.restrictField
(
coarseDiag,
fineMatrix.diag(),
fineLevelIndex,
false // no processor agglomeration
);
// Get reference to fine-level interfaces
const lduInterfaceFieldPtrsList& fineInterfaces =
interfaceLevel(fineLevelIndex);
// Create coarse-level interfaces
primitiveInterfaceLevels_.set
(
fineLevelIndex,
new PtrList<lduInterfaceField>(fineInterfaces.size())
);
PtrList<lduInterfaceField>& coarsePrimInterfaces =
primitiveInterfaceLevels_[fineLevelIndex];
interfaceLevels_.set
(
fineLevelIndex,
new lduInterfaceFieldPtrsList(fineInterfaces.size())
);
lduInterfaceFieldPtrsList& coarseInterfaces =
interfaceLevels_[fineLevelIndex];
// Set coarse-level boundary coefficients
interfaceLevelsBouCoeffs_.set
(
fineLevelIndex,
new FieldField<Field, scalar>(fineInterfaces.size())
);
FieldField<Field, scalar>& coarseInterfaceBouCoeffs =
interfaceLevelsBouCoeffs_[fineLevelIndex];
// Set coarse-level internal coefficients
interfaceLevelsIntCoeffs_.set
(
fineLevelIndex,
new FieldField<Field, scalar>(fineInterfaces.size())
);
FieldField<Field, scalar>& coarseInterfaceIntCoeffs =
interfaceLevelsIntCoeffs_[fineLevelIndex];
// Add the coarse level
agglomerateInterfaceCoefficients
(
fineLevelIndex,
coarseMeshInterfaces,
coarsePrimInterfaces,
coarseInterfaces,
coarseInterfaceBouCoeffs,
coarseInterfaceIntCoeffs
);
// Get face restriction map for current level
const labelList& faceRestrictAddr =
agglomeration_.faceRestrictAddressing(fineLevelIndex);
const boolList& faceFlipMap =
agglomeration_.faceFlipMap(fineLevelIndex);
// Check if matrix is asymetric and if so agglomerate both upper
// and lower coefficients ...
if (fineMatrix.hasLower())
{
// Get off-diagonal matrix coefficients
const scalarField& fineUpper = fineMatrix.upper();
const scalarField& fineLower = fineMatrix.lower();
// Coarse matrix upper coefficients. Note passed in size
scalarField& coarseUpper = coarseMatrix.upper(nCoarseFaces);
scalarField& coarseLower = coarseMatrix.lower(nCoarseFaces);
forAll(faceRestrictAddr, fineFacei)
{
label cFace = faceRestrictAddr[fineFacei];
if (cFace >= 0)
{
// Check the orientation of the fine-face relative to the
// coarse face it is being agglomerated into
if (!faceFlipMap[fineFacei])
{
coarseUpper[cFace] += fineUpper[fineFacei];
coarseLower[cFace] += fineLower[fineFacei];
}
else
{
coarseUpper[cFace] += fineLower[fineFacei];
coarseLower[cFace] += fineUpper[fineFacei];
}
}
else
{
// Add the fine face coefficients into the diagonal.
coarseDiag[-1 - cFace] +=
fineUpper[fineFacei] + fineLower[fineFacei];
}
}
}
else // ... Otherwise it is symmetric so agglomerate just the upper
{
// Get off-diagonal matrix coefficients
const scalarField& fineUpper = fineMatrix.upper();
// Coarse matrix upper coefficients
scalarField& coarseUpper = coarseMatrix.upper(nCoarseFaces);
forAll(faceRestrictAddr, fineFacei)
{
label cFace = faceRestrictAddr[fineFacei];
if (cFace >= 0)
{
coarseUpper[cFace] += fineUpper[fineFacei];
}
else
{
// Add the fine face coefficient into the diagonal.
coarseDiag[-1 - cFace] += 2*fineUpper[fineFacei];
}
}
}
}
}
// Agglomerate only the interface coefficients.
void Foam::GAMGSolver::agglomerateInterfaceCoefficients
(
const label fineLevelIndex,
const lduInterfacePtrsList& coarseMeshInterfaces,
PtrList<lduInterfaceField>& coarsePrimInterfaces,
lduInterfaceFieldPtrsList& coarseInterfaces,
FieldField<Field, scalar>& coarseInterfaceBouCoeffs,
FieldField<Field, scalar>& coarseInterfaceIntCoeffs
) const
{
// Get reference to fine-level interfaces
const lduInterfaceFieldPtrsList& fineInterfaces =
interfaceLevel(fineLevelIndex);
@ -61,34 +216,12 @@ void Foam::GAMGSolver::agglomerateMatrix(const label fineLevelIndex)
const FieldField<Field, scalar>& fineInterfaceIntCoeffs =
interfaceIntCoeffsLevel(fineLevelIndex);
const labelListList& patchFineToCoarse =
agglomeration_.patchFaceRestrictAddressing(fineLevelIndex);
// Create coarse-level interfaces
interfaceLevels_.set
(
fineLevelIndex,
new lduInterfaceFieldPtrsList(fineInterfaces.size())
);
const labelList& nPatchFaces =
agglomeration_.nPatchFaces(fineLevelIndex);
lduInterfaceFieldPtrsList& coarseInterfaces =
interfaceLevels_[fineLevelIndex];
// Set coarse-level boundary coefficients
interfaceLevelsBouCoeffs_.set
(
fineLevelIndex,
new FieldField<Field, scalar>(fineInterfaces.size())
);
FieldField<Field, scalar>& coarseInterfaceBouCoeffs =
interfaceLevelsBouCoeffs_[fineLevelIndex];
// Set coarse-level internal coefficients
interfaceLevelsIntCoeffs_.set
(
fineLevelIndex,
new FieldField<Field, scalar>(fineInterfaces.size())
);
FieldField<Field, scalar>& coarseInterfaceIntCoeffs =
interfaceLevelsIntCoeffs_[fineLevelIndex];
// Add the coarse level
forAll(fineInterfaces, inti)
@ -98,10 +231,10 @@ void Foam::GAMGSolver::agglomerateMatrix(const label fineLevelIndex)
const GAMGInterface& coarseInterface =
refCast<const GAMGInterface>
(
agglomeration_.interfaceLevel(fineLevelIndex + 1)[inti]
coarseMeshInterfaces[inti]
);
coarseInterfaces.set
coarsePrimInterfaces.set
(
inti,
GAMGInterfaceField::New
@ -110,101 +243,553 @@ void Foam::GAMGSolver::agglomerateMatrix(const label fineLevelIndex)
fineInterfaces[inti]
).ptr()
);
coarseInterfaces.set
(
inti,
&coarsePrimInterfaces[inti]
);
const labelList& faceRestrictAddressing = patchFineToCoarse[inti];
coarseInterfaceBouCoeffs.set
(
inti,
coarseInterface.agglomerateCoeffs(fineInterfaceBouCoeffs[inti])
new scalarField(nPatchFaces[inti], 0.0)
);
agglomeration_.restrictField
(
coarseInterfaceBouCoeffs[inti],
fineInterfaceBouCoeffs[inti],
faceRestrictAddressing
);
coarseInterfaceIntCoeffs.set
(
inti,
coarseInterface.agglomerateCoeffs(fineInterfaceIntCoeffs[inti])
new scalarField(nPatchFaces[inti], 0.0)
);
agglomeration_.restrictField
(
coarseInterfaceIntCoeffs[inti],
fineInterfaceIntCoeffs[inti],
faceRestrictAddressing
);
}
}
}
// Check if matrix is assymetric and if so agglomerate both upper and lower
// coefficients ...
if (fineMatrix.hasLower())
// Gather matrices.
// Note: matrices get constructed with dummy mesh
void Foam::GAMGSolver::gatherMatrices
(
const labelList& procIDs,
const lduMesh& dummyMesh,
const label meshComm,
const lduMatrix& mat,
const FieldField<Field, scalar>& interfaceBouCoeffs,
const FieldField<Field, scalar>& interfaceIntCoeffs,
const lduInterfaceFieldPtrsList& interfaces,
PtrList<lduMatrix>& otherMats,
PtrList<FieldField<Field, scalar> >& otherBouCoeffs,
PtrList<FieldField<Field, scalar> >& otherIntCoeffs,
List<boolList>& otherTransforms,
List<List<int> >& otherRanks
) const
{
if (debug)
{
// Get off-diagonal matrix coefficients
const scalarField& fineUpper = fineMatrix.upper();
const scalarField& fineLower = fineMatrix.lower();
Pout<< "GAMGSolver::gatherMatrices :"
<< " collecting matrices from procs:" << procIDs
<< " using comm:" << meshComm << endl;
}
// Coarse matrix upper coefficients
scalarField& coarseUpper = coarseMatrix.upper();
scalarField& coarseLower = coarseMatrix.lower();
if (Pstream::myProcNo(meshComm) == procIDs[0])
{
// Master.
otherMats.setSize(procIDs.size()-1);
otherBouCoeffs.setSize(procIDs.size()-1);
otherIntCoeffs.setSize(procIDs.size()-1);
otherTransforms.setSize(procIDs.size()-1);
otherRanks.setSize(procIDs.size()-1);
const labelList& restrictAddr =
agglomeration_.restrictAddressing(fineLevelIndex);
const labelUList& l = fineMatrix.lduAddr().lowerAddr();
const labelUList& cl = coarseMatrix.lduAddr().lowerAddr();
const labelUList& cu = coarseMatrix.lduAddr().upperAddr();
forAll(faceRestrictAddr, fineFacei)
for (label procI = 1; procI < procIDs.size(); procI++)
{
label cFace = faceRestrictAddr[fineFacei];
label otherI = procI-1;
if (cFace >= 0)
IPstream fromSlave
(
Pstream::scheduled,
procIDs[procI],
0, // bufSize
Pstream::msgType(),
meshComm
);
otherMats.set(otherI, new lduMatrix(dummyMesh, fromSlave));
// Receive number of/valid interfaces
boolList& procTransforms = otherTransforms[otherI];
List<int>& procRanks = otherRanks[otherI];
fromSlave >> procTransforms;
fromSlave >> procRanks;
// Size coefficients
otherBouCoeffs.set
(
otherI,
new FieldField<Field, scalar>(procRanks.size())
);
otherIntCoeffs.set
(
otherI,
new FieldField<Field, scalar>(procRanks.size())
);
forAll(procRanks, intI)
{
// Check the orientation of the fine-face relative to the
// coarse face it is being agglomerated into
if (cl[cFace] == restrictAddr[l[fineFacei]])
if (procRanks[intI] != -1)
{
coarseUpper[cFace] += fineUpper[fineFacei];
coarseLower[cFace] += fineLower[fineFacei];
}
else if (cu[cFace] == restrictAddr[l[fineFacei]])
{
coarseUpper[cFace] += fineLower[fineFacei];
coarseLower[cFace] += fineUpper[fineFacei];
}
else
{
FatalErrorIn
otherBouCoeffs[otherI].set
(
"GAMGSolver::agglomerateMatrix(const label)"
) << "Inconsistent addressing between "
"fine and coarse grids"
<< exit(FatalError);
intI,
new scalarField(fromSlave)
);
otherIntCoeffs[otherI].set
(
intI,
new scalarField(fromSlave)
);
}
}
else
{
// Add the fine face coefficients into the diagonal.
coarseDiag[-1 - cFace] +=
fineUpper[fineFacei] + fineLower[fineFacei];
}
}
}
else // ... Otherwise it is symmetric so agglomerate just the upper
else
{
// Get off-diagonal matrix coefficients
const scalarField& fineUpper = fineMatrix.upper();
// Send to master
// Coarse matrix upper coefficients
scalarField& coarseUpper = coarseMatrix.upper();
forAll(faceRestrictAddr, fineFacei)
// Count valid interfaces
boolList procTransforms(interfaceBouCoeffs.size(), false);
List<int> procRanks(interfaceBouCoeffs.size(), -1);
forAll(interfaces, intI)
{
label cFace = faceRestrictAddr[fineFacei];
if (interfaces.set(intI))
{
const processorLduInterfaceField& interface =
refCast<const processorLduInterfaceField>
(
interfaces[intI]
);
if (cFace >= 0)
{
coarseUpper[cFace] += fineUpper[fineFacei];
procTransforms[intI] = interface.doTransform();
procRanks[intI] = interface.rank();
}
else
}
OPstream toMaster
(
Pstream::scheduled,
procIDs[0],
0,
Pstream::msgType(),
meshComm
);
toMaster << mat << procTransforms << procRanks;
forAll(procRanks, intI)
{
if (procRanks[intI] != -1)
{
// Add the fine face coefficient into the diagonal.
coarseDiag[-1 - cFace] += 2*fineUpper[fineFacei];
toMaster
<< interfaceBouCoeffs[intI]
<< interfaceIntCoeffs[intI];
}
}
}
}
void Foam::GAMGSolver::procAgglomerateMatrix
(
// Agglomeration information
const labelList& procAgglomMap,
const List<int>& agglomProcIDs,
const label levelI,
// Resulting matrix
autoPtr<lduMatrix>& allMatrixPtr,
FieldField<Field, scalar>& allInterfaceBouCoeffs,
FieldField<Field, scalar>& allInterfaceIntCoeffs,
PtrList<lduInterfaceField>& allPrimitiveInterfaces,
lduInterfaceFieldPtrsList& allInterfaces
) const
{
const lduMatrix& coarsestMatrix = matrixLevels_[levelI];
const lduInterfaceFieldPtrsList& coarsestInterfaces =
interfaceLevels_[levelI];
const FieldField<Field, scalar>& coarsestBouCoeffs =
interfaceLevelsBouCoeffs_[levelI];
const FieldField<Field, scalar>& coarsestIntCoeffs =
interfaceLevelsIntCoeffs_[levelI];
const lduMesh& coarsestMesh = coarsestMatrix.mesh();
label coarseComm = coarsestMesh.comm();
label oldWarn = UPstream::warnComm;
UPstream::warnComm = coarseComm;
// Gather all matrix coefficients onto agglomProcIDs[0]
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
PtrList<lduMatrix> otherMats;
PtrList<FieldField<Field, scalar> > otherBouCoeffs;
PtrList<FieldField<Field, scalar> > otherIntCoeffs;
List<boolList> otherTransforms;
List<List<int> > otherRanks;
gatherMatrices
(
agglomProcIDs,
coarsestMesh,
coarseComm,
coarsestMatrix,
coarsestBouCoeffs,
coarsestIntCoeffs,
coarsestInterfaces,
otherMats,
otherBouCoeffs,
otherIntCoeffs,
otherTransforms,
otherRanks
);
if (Pstream::myProcNo(coarseComm) == agglomProcIDs[0])
{
// Agglomerate all matrix
// ~~~~~~~~~~~~~~~~~~~~~~
//Pout<< "Own matrix:" << coarsestMatrix.info() << endl;
//
//forAll(otherMats, i)
//{
// Pout<< "** otherMats " << i << " "
// << otherMats[i].info()
// << endl;
//}
//Pout<< endl;
const lduMesh& allMesh = agglomeration_.meshLevel(levelI+1);
const labelList& cellOffsets = agglomeration_.cellOffsets(levelI+1);
const labelListList& faceMap = agglomeration_.faceMap(levelI+1);
const labelListList& boundaryMap = agglomeration_.boundaryMap(levelI+1);
const labelListListList& boundaryFaceMap =
agglomeration_.boundaryFaceMap(levelI+1);
allMatrixPtr.reset(new lduMatrix(allMesh));
lduMatrix& allMatrix = allMatrixPtr();
if (coarsestMatrix.hasDiag())
{
scalarField& allDiag = allMatrix.diag();
SubList<scalar>
(
allDiag,
coarsestMatrix.diag().size()
).assign
(
coarsestMatrix.diag()
);
forAll(otherMats, i)
{
SubList<scalar>
(
allDiag,
otherMats[i].diag().size(),
cellOffsets[i+1]
).assign
(
otherMats[i].diag()
);
}
}
if (coarsestMatrix.hasLower())
{
scalarField& allLower = allMatrix.lower();
UIndirectList<scalar>
(
allLower,
faceMap[0]
) = coarsestMatrix.lower();
forAll(otherMats, i)
{
UIndirectList<scalar>
(
allLower,
faceMap[i+1]
) = otherMats[i].lower();
}
}
if (coarsestMatrix.hasUpper())
{
scalarField& allUpper = allMatrix.upper();
UIndirectList<scalar>
(
allUpper,
faceMap[0]
) = coarsestMatrix.upper();
forAll(otherMats, i)
{
UIndirectList<scalar>
(
allUpper,
faceMap[i+1]
) = otherMats[i].upper();
}
}
// Agglomerate interface fields and coefficients
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
lduInterfacePtrsList allMeshInterfaces = allMesh.interfaces();
allInterfaceBouCoeffs.setSize(allMeshInterfaces.size());
allInterfaceIntCoeffs.setSize(allMeshInterfaces.size());
allPrimitiveInterfaces.setSize(allMeshInterfaces.size());
allInterfaces.setSize(allMeshInterfaces.size());
forAll(allMeshInterfaces, intI)
{
const lduInterface& patch = allMeshInterfaces[intI];
label size = patch.faceCells().size();
allInterfaceBouCoeffs.set(intI, new scalarField(size));
allInterfaceIntCoeffs.set(intI, new scalarField(size));
}
labelList nBounFaces(allMeshInterfaces.size());
forAll(boundaryMap, procI)
{
const FieldField<Field, scalar>& procBouCoeffs
(
(procI == 0)
? coarsestBouCoeffs
: otherBouCoeffs[procI-1]
);
const FieldField<Field, scalar>& procIntCoeffs
(
(procI == 0)
? coarsestIntCoeffs
: otherIntCoeffs[procI-1]
);
const labelList& bMap = boundaryMap[procI];
forAll(bMap, procIntI)
{
label allIntI = bMap[procIntI];
if (allIntI != -1)
{
// So this boundary has been preserved. Copy
// data across.
if (!allInterfaces.set(allIntI))
{
// Construct lduInterfaceField
bool doTransform = false;
int rank = -1;
if (procI == 0)
{
const processorGAMGInterfaceField& procInt =
refCast
<
const processorGAMGInterfaceField
>
(
coarsestInterfaces[procIntI]
);
doTransform = procInt.doTransform();
rank = procInt.rank();
}
else
{
doTransform =
otherTransforms[procI-1][procIntI];
rank = otherRanks[procI-1][procIntI];
}
allPrimitiveInterfaces.set
(
allIntI,
GAMGInterfaceField::New
(
refCast<const GAMGInterface>
(
allMeshInterfaces[allIntI]
),
doTransform,
rank
).ptr()
);
allInterfaces.set
(
allIntI,
&allPrimitiveInterfaces[allIntI]
);
}
// Map data from processor to complete mesh
scalarField& allBou = allInterfaceBouCoeffs[allIntI];
scalarField& allInt = allInterfaceIntCoeffs[allIntI];
const labelList& map = boundaryFaceMap[procI][procIntI];
const scalarField& procBou = procBouCoeffs[procIntI];
const scalarField& procInt = procIntCoeffs[procIntI];
forAll(map, i)
{
label allFaceI = map[i];
if (allFaceI < 0)
{
FatalErrorIn("GAMGSolver::GAMGSolver()")
<< "problem." << abort(FatalError);
}
allBou[allFaceI] = procBou[i];
allInt[allFaceI] = procInt[i];
}
}
else if (procBouCoeffs.set(procIntI))
{
// Boundary has become internal face
const labelList& map = boundaryFaceMap[procI][procIntI];
const scalarField& procBou = procBouCoeffs[procIntI];
const scalarField& procInt = procIntCoeffs[procIntI];
forAll(map, i)
{
if (map[i] >= 0)
{
label allFaceI = map[i];
if (coarsestMatrix.hasUpper())
{
allMatrix.upper()[allFaceI] = -procBou[i];
}
if (coarsestMatrix.hasLower())
{
allMatrix.lower()[allFaceI] = -procInt[i];
}
}
else
{
label allFaceI = -map[i]-1;
if (coarsestMatrix.hasUpper())
{
allMatrix.upper()[allFaceI] = -procInt[i];
}
if (coarsestMatrix.hasLower())
{
allMatrix.lower()[allFaceI] = -procBou[i];
}
}
}
}
}
}
//Pout<< "** Assembled allMatrix:" << allMatrix.info() << endl;
//
//forAll(allInterfaces, intI)
//{
// if (allInterfaces.set(intI))
// {
// Pout<< " patch:" << intI
// << " type:" << allInterfaces[intI].type()
// << " size:"
// << allInterfaces[intI].interface().
// faceCells().size()
// << endl;
//
// //const scalarField& bouCoeffs = allInterfaceBouCoeffs[intI];
// //const scalarField& intCoeffs = allInterfaceIntCoeffs[intI];
// //forAll(bouCoeffs, faceI)
// //{
// // Pout<< " " << faceI
// // << "\tbou:" << bouCoeffs[faceI]
// // << "\tint:" << intCoeffs[faceI]
// // << endl;
// //}
// }
//}
}
UPstream::warnComm = oldWarn;
}
void Foam::GAMGSolver::procAgglomerateMatrix
(
const labelList& procAgglomMap,
const List<int>& agglomProcIDs,
const label levelI
)
{
autoPtr<lduMatrix> allMatrixPtr;
autoPtr<FieldField<Field, scalar> > allInterfaceBouCoeffs
(
new FieldField<Field, scalar>(0)
);
autoPtr<FieldField<Field, scalar> > allInterfaceIntCoeffs
(
new FieldField<Field, scalar>(0)
);
autoPtr<PtrList<lduInterfaceField> > allPrimitiveInterfaces
(
new PtrList<lduInterfaceField>(0)
);
autoPtr<lduInterfaceFieldPtrsList> allInterfaces
(
new lduInterfaceFieldPtrsList(0)
);
procAgglomerateMatrix
(
// Agglomeration information
procAgglomMap,
agglomProcIDs,
levelI,
// Resulting matrix
allMatrixPtr,
allInterfaceBouCoeffs(),
allInterfaceIntCoeffs(),
allPrimitiveInterfaces(),
allInterfaces()
);
matrixLevels_.set(levelI, allMatrixPtr);
interfaceLevelsBouCoeffs_.set(levelI, allInterfaceBouCoeffs);
interfaceLevelsIntCoeffs_.set(levelI, allInterfaceIntCoeffs);
primitiveInterfaceLevels_.set(levelI, allPrimitiveInterfaces);
interfaceLevels_.set(levelI, allInterfaces);
}
// ************************************************************************* //

View File

@ -58,7 +58,8 @@ void Foam::GAMGSolver::scale
}
vector2D scalingVector(scalingFactorNum, scalingFactorDenom);
reduce(scalingVector, sumOp<vector2D>());
A.mesh().reduce(scalingVector, sumOp<vector2D>());
scalar sf = scalingVector.x()/stabilise(scalingVector.y(), VSMALL);
if (debug >= 2)

View File

@ -28,7 +28,6 @@ License
#include "BICCG.H"
#include "SubField.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::solverPerformance Foam::GAMGSolver::solve
@ -61,7 +60,11 @@ Foam::solverPerformance Foam::GAMGSolver::solve
scalarField finestResidual(source - Apsi);
// Calculate normalised residual for convergence test
solverPerf.initialResidual() = gSumMag(finestResidual)/normFactor;
solverPerf.initialResidual() = gSumMag
(
finestResidual,
matrix().mesh().comm()
)/normFactor;
solverPerf.finalResidual() = solverPerf.initialResidual();
@ -77,8 +80,20 @@ Foam::solverPerformance Foam::GAMGSolver::solve
// Create the smoothers for all levels
PtrList<lduMatrix::smoother> smoothers;
// Scratch fields if processor-agglomerated coarse level meshes
// are bigger than original. Usually not needed
scalarField scratch1;
scalarField scratch2;
// Initialise the above data structures
initVcycle(coarseCorrFields, coarseSources, smoothers);
initVcycle
(
coarseCorrFields,
coarseSources,
smoothers,
scratch1,
scratch2
);
do
{
@ -90,21 +105,36 @@ Foam::solverPerformance Foam::GAMGSolver::solve
Apsi,
finestCorrection,
finestResidual,
(scratch1.size() ? scratch1 : Apsi),
(scratch2.size() ? scratch2 : finestCorrection),
coarseCorrFields,
coarseSources,
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;
finestResidual -= Apsi;
solverPerf.finalResidual() = gSumMag(finestResidual)/normFactor;
solverPerf.finalResidual() = gSumMag
(
finestResidual,
matrix().mesh().comm()
)/normFactor;
if (debug >= 2)
{
solverPerf.print(Info);
solverPerf.print(Info(matrix().mesh().comm()));
}
} while
(
@ -125,6 +155,10 @@ void Foam::GAMGSolver::Vcycle
scalarField& Apsi,
scalarField& finestCorrection,
scalarField& finestResidual,
scalarField& scratch1,
scalarField& scratch2,
PtrList<scalarField>& coarseCorrFields,
PtrList<scalarField>& coarseSources,
const direction cmpt
@ -134,8 +168,8 @@ void Foam::GAMGSolver::Vcycle
const label coarsestLevel = matrixLevels_.size() - 1;
// Restrict finest grid residual for the next level up
agglomeration_.restrictField(coarseSources[0], finestResidual, 0);
// Restrict finest grid residual for the next level up.
agglomeration_.restrictField(coarseSources[0], finestResidual, 0, true);
if (debug >= 2 && nPreSweeps_)
{
@ -146,66 +180,80 @@ void Foam::GAMGSolver::Vcycle
// Residual restriction (going to coarser levels)
for (label leveli = 0; leveli < coarsestLevel; leveli++)
{
// If the optional pre-smoothing sweeps are selected
// smooth the coarse-grid field for the restriced source
if (nPreSweeps_)
if (coarseSources.set(leveli + 1))
{
coarseCorrFields[leveli] = 0.0;
smoothers[leveli + 1].smooth
(
coarseCorrFields[leveli],
coarseSources[leveli],
cmpt,
min
(
nPreSweeps_ + preSweepsLevelMultiplier_*leveli,
maxPreSweeps_
)
);
scalarField::subField ACf
(
Apsi,
coarseCorrFields[leveli].size()
);
// Scale coarse-grid correction field
// but not on the coarsest level because it evaluates to 1
if (scaleCorrection_ && leveli < coarsestLevel - 1)
// If the optional pre-smoothing sweeps are selected
// smooth the coarse-grid field for the restriced source
if (nPreSweeps_)
{
scale
coarseCorrFields[leveli] = 0.0;
smoothers[leveli + 1].smooth
(
coarseCorrFields[leveli],
const_cast<scalarField&>(ACf.operator const scalarField&()),
matrixLevels_[leveli],
coarseSources[leveli],
cmpt,
min
(
nPreSweeps_ + preSweepsLevelMultiplier_*leveli,
maxPreSweeps_
)
);
scalarField::subField ACf
(
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
if (scaleCorrection_ && leveli < coarsestLevel - 1)
{
scale
(
coarseCorrFields[leveli],
const_cast<scalarField&>
(
ACf.operator const scalarField&()
),
//ACf,
matrixLevels_[leveli],
interfaceLevelsBouCoeffs_[leveli],
interfaceLevels_[leveli],
coarseSources[leveli],
cmpt
);
}
// Correct the residual with the new solution
matrixLevels_[leveli].Amul
(
const_cast<scalarField&>
(
ACf.operator const scalarField&()
),
//ACf,
coarseCorrFields[leveli],
interfaceLevelsBouCoeffs_[leveli],
interfaceLevels_[leveli],
coarseSources[leveli],
cmpt
);
coarseSources[leveli] -= ACf;
}
// Correct the residual with the new solution
matrixLevels_[leveli].Amul
// Residual is equal to source
agglomeration_.restrictField
(
const_cast<scalarField&>(ACf.operator const scalarField&()),
coarseCorrFields[leveli],
interfaceLevelsBouCoeffs_[leveli],
interfaceLevels_[leveli],
cmpt
coarseSources[leveli + 1],
coarseSources[leveli],
leveli + 1,
true
);
coarseSources[leveli] -= ACf;
}
// Residual is equal to source
agglomeration_.restrictField
(
coarseSources[leveli + 1],
coarseSources[leveli],
leveli + 1
);
}
if (debug >= 2 && nPreSweeps_)
@ -215,12 +263,14 @@ void Foam::GAMGSolver::Vcycle
// Solve Coarsest level with either an iterative or direct solver
solveCoarsestLevel
(
coarseCorrFields[coarsestLevel],
coarseSources[coarsestLevel]
);
if (coarseCorrFields.set(coarsestLevel))
{
solveCoarsestLevel
(
coarseCorrFields[coarsestLevel],
coarseSources[coarsestLevel]
);
}
if (debug >= 2)
{
@ -229,87 +279,111 @@ void Foam::GAMGSolver::Vcycle
// Smoothing and prolongation of the coarse correction fields
// (going to finer levels)
scalarField dummyField(0);
for (label leveli = coarsestLevel - 1; leveli >= 0; leveli--)
{
// Create a field for the pre-smoothed correction field
// as a sub-field of the finestCorrection which is not
// currently being used
scalarField::subField preSmoothedCoarseCorrField
(
finestCorrection,
coarseCorrFields[leveli].size()
);
// Only store the preSmoothedCoarseCorrField if pre-smoothing is used
if (nPreSweeps_)
if (coarseCorrFields.set(leveli))
{
preSmoothedCoarseCorrField.assign(coarseCorrFields[leveli]);
}
// Create a field for the pre-smoothed correction field
// as a sub-field of the finestCorrection which is not
// currently being used
scalarField::subField preSmoothedCoarseCorrField
(
scratch2,
coarseCorrFields[leveli].size()
);
//scalarField preSmoothedCoarseCorrField
//(
// coarseCorrFields[leveli].size(),
// VGREAT
//);
agglomeration_.prolongField
(
coarseCorrFields[leveli],
coarseCorrFields[leveli + 1],
leveli + 1
);
// Only store the preSmoothedCoarseCorrField if pre-smoothing is
// used
if (nPreSweeps_)
{
preSmoothedCoarseCorrField.assign(coarseCorrFields[leveli]);
}
// Create A.psi for this coarse level as a sub-field of Apsi
scalarField::subField ACf
(
Apsi,
coarseCorrFields[leveli].size()
);
scalarField& ACfRef =
const_cast<scalarField&>(ACf.operator const scalarField&());
if (interpolateCorrection_)
{
interpolate
agglomeration_.prolongField
(
coarseCorrFields[leveli],
ACfRef,
matrixLevels_[leveli],
interfaceLevelsBouCoeffs_[leveli],
interfaceLevels_[leveli],
coarseSources[leveli],
cmpt
(
coarseCorrFields.set(leveli + 1)
? coarseCorrFields[leveli + 1]
: dummyField // dummy value
),
leveli + 1,
true
);
}
// Scale coarse-grid correction field
// but not on the coarsest level because it evaluates to 1
if (scaleCorrection_ && leveli < coarsestLevel - 1)
{
scale
// Create A.psi for this coarse level as a sub-field of Apsi
scalarField::subField ACf
(
scratch1,
coarseCorrFields[leveli].size()
);
scalarField& ACfRef =
const_cast<scalarField&>(ACf.operator const scalarField&());
//scalarField ACfRef
//(
// coarseCorrFields[leveli].size(),
// VGREAT
//);
if (interpolateCorrection_)
{
interpolate
(
coarseCorrFields[leveli],
ACfRef,
matrixLevels_[leveli],
interfaceLevelsBouCoeffs_[leveli],
interfaceLevels_[leveli],
coarseSources[leveli],
cmpt
);
}
// Scale coarse-grid correction field
// but not on the coarsest level because it evaluates to 1
if (scaleCorrection_ && leveli < coarsestLevel - 1)
{
scale
(
coarseCorrFields[leveli],
ACfRef,
matrixLevels_[leveli],
interfaceLevelsBouCoeffs_[leveli],
interfaceLevels_[leveli],
coarseSources[leveli],
cmpt
);
}
// Only add the preSmoothedCoarseCorrField if pre-smoothing is
// used
if (nPreSweeps_)
{
coarseCorrFields[leveli] += preSmoothedCoarseCorrField;
}
smoothers[leveli + 1].smooth
(
coarseCorrFields[leveli],
ACfRef,
matrixLevels_[leveli],
interfaceLevelsBouCoeffs_[leveli],
interfaceLevels_[leveli],
coarseSources[leveli],
cmpt
cmpt,
min
(
nPostSweeps_ + postSweepsLevelMultiplier_*leveli,
maxPostSweeps_
)
);
}
// Only add the preSmoothedCoarseCorrField if pre-smoothing is used
if (nPreSweeps_)
{
coarseCorrFields[leveli] += preSmoothedCoarseCorrField;
}
smoothers[leveli + 1].smooth
(
coarseCorrFields[leveli],
coarseSources[leveli],
cmpt,
min
(
nPostSweeps_ + postSweepsLevelMultiplier_*leveli,
maxPostSweeps_
)
);
}
// Prolong the finest level correction
@ -317,7 +391,8 @@ void Foam::GAMGSolver::Vcycle
(
finestCorrection,
coarseCorrFields[0],
0
0,
true
);
if (interpolateCorrection_)
@ -368,9 +443,13 @@ void Foam::GAMGSolver::initVcycle
(
PtrList<scalarField>& coarseCorrFields,
PtrList<scalarField>& coarseSources,
PtrList<lduMatrix::smoother>& smoothers
PtrList<lduMatrix::smoother>& smoothers,
scalarField& scratch1,
scalarField& scratch2
) const
{
label maxSize = matrix_.diag().size();
coarseCorrFields.setSize(matrixLevels_.size());
coarseSources.setSize(matrixLevels_.size());
smoothers.setSize(matrixLevels_.size() + 1);
@ -392,37 +471,44 @@ void Foam::GAMGSolver::initVcycle
forAll(matrixLevels_, leveli)
{
coarseCorrFields.set
(
leveli,
new scalarField
(
agglomeration_.meshLevel(leveli + 1).lduAddr().size()
)
);
if (agglomeration_.nCells(leveli) >= 0)
{
label nCoarseCells = agglomeration_.nCells(leveli);
coarseSources.set
(
leveli,
new scalarField
(
agglomeration_.meshLevel(leveli + 1).lduAddr().size()
)
);
coarseSources.set(leveli, new scalarField(nCoarseCells));
}
smoothers.set
(
leveli + 1,
lduMatrix::smoother::New
if (matrixLevels_.set(leveli))
{
const lduMatrix& mat = matrixLevels_[leveli];
label nCoarseCells = mat.diag().size();
maxSize = max(maxSize, nCoarseCells);
coarseCorrFields.set(leveli, new scalarField(nCoarseCells));
smoothers.set
(
fieldName_,
matrixLevels_[leveli],
interfaceLevelsBouCoeffs_[leveli],
interfaceLevelsIntCoeffs_[leveli],
interfaceLevels_[leveli],
controlDict_
)
);
leveli + 1,
lduMatrix::smoother::New
(
fieldName_,
matrixLevels_[leveli],
interfaceLevelsBouCoeffs_[leveli],
interfaceLevelsIntCoeffs_[leveli],
interfaceLevels_[leveli],
controlDict_
)
);
}
}
if (maxSize > matrix_.diag().size())
{
// Allocate some scratch storage
scratch1.setSize(maxSize);
scratch2.setSize(maxSize);
}
}
@ -433,14 +519,124 @@ void Foam::GAMGSolver::solveCoarsestLevel
const scalarField& coarsestSource
) const
{
const label coarsestLevel = matrixLevels_.size() - 1;
label coarseComm = matrixLevels_[coarsestLevel].mesh().comm();
label oldWarn = UPstream::warnComm;
UPstream::warnComm = coarseComm;
if (directSolveCoarsest_)
{
coarsestCorrField = coarsestSource;
coarsestLUMatrixPtr_->solve(coarsestCorrField);
}
//else if
//(
// agglomeration_.processorAgglomerate()
// && procMatrixLevels_.set(coarsestLevel)
//)
//{
// //const labelList& agglomProcIDs = agglomeration_.agglomProcIDs
// //(
// // coarsestLevel
// //);
// //
// //scalarField allSource;
// //
// //globalIndex cellOffsets;
// //if (Pstream::myProcNo(coarseComm) == agglomProcIDs[0])
// //{
// // cellOffsets.offsets() =
// // agglomeration_.cellOffsets(coarsestLevel);
// //}
// //
// //cellOffsets.gather
// //(
// // coarseComm,
// // agglomProcIDs,
// // coarsestSource,
// // allSource
// //);
// //
// //scalarField allCorrField;
// //solverPerformance coarseSolverPerf;
//
// label solveComm = agglomeration_.procCommunicator(coarsestLevel);
// label oldWarn = UPstream::warnComm;
// UPstream::warnComm = solveComm;
//
//
// coarsestCorrField = 0;
// solverPerformance coarseSolverPerf;
//
// if (Pstream::myProcNo(solveComm) != -1)
// {
// const lduMatrix& allMatrix = procMatrixLevels_[coarsestLevel];
//
// {
// Pout<< "** Master:Solving on comm:" << solveComm
// << " with procs:" << UPstream::procID(solveComm) << endl;
//
// if (allMatrix.asymmetric())
// {
// coarseSolverPerf = BICCG
// (
// "coarsestLevelCorr",
// allMatrix,
// procInterfaceLevelsBouCoeffs_[coarsestLevel],
// procInterfaceLevelsIntCoeffs_[coarsestLevel],
// procInterfaceLevels_[coarsestLevel],
// tolerance_,
// relTol_
// ).solve
// (
// coarsestCorrField,
// coarsestSource
// );
// }
// else
// {
// coarseSolverPerf = ICCG
// (
// "coarsestLevelCorr",
// allMatrix,
// procInterfaceLevelsBouCoeffs_[coarsestLevel],
// procInterfaceLevelsIntCoeffs_[coarsestLevel],
// procInterfaceLevels_[coarsestLevel],
// tolerance_,
// relTol_
// ).solve
// (
// coarsestCorrField,
// coarsestSource
// );
// }
// }
// }
//
// UPstream::warnComm = oldWarn;
// Pout<< "done master solve." << endl;
//
// //// Scatter to all processors
// //coarsestCorrField.setSize(coarsestSource.size());
// //cellOffsets.scatter
// //(
// // coarseComm,
// // agglomProcIDs,
// // allCorrField,
// // coarsestCorrField
// //);
//
// if (debug >= 2)
// {
// coarseSolverPerf.print(Info(coarseComm));
// }
//
// Pout<< "procAgglom: coarsestSource :" << coarsestSource << endl;
// Pout<< "procAgglom: coarsestCorrField:" << coarsestCorrField << endl;
//}
else
{
const label coarsestLevel = matrixLevels_.size() - 1;
coarsestCorrField = 0;
solverPerformance coarseSolverPerf;
@ -481,9 +677,11 @@ void Foam::GAMGSolver::solveCoarsestLevel
if (debug >= 2)
{
coarseSolverPerf.print(Info);
coarseSolverPerf.print(Info(coarseComm));
}
}
UPstream::warnComm = oldWarn;
}

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
@ -32,6 +32,7 @@ namespace Foam
{
defineTypeNameAndDebug(GAMGInterfaceField, 0);
defineRunTimeSelectionTable(GAMGInterfaceField, lduInterface);
defineRunTimeSelectionTable(GAMGInterfaceField, lduInterfaceField);
}
// ************************************************************************* //

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
@ -29,7 +29,7 @@ Description
SourceFiles
GAMGInterfaceField.C
newAmgInterfaceField.C
GAMGInterfaceFieldNew.C
\*---------------------------------------------------------------------------*/
@ -80,7 +80,7 @@ public:
(
autoPtr,
GAMGInterfaceField,
lduInterface,
lduInterfaceField,
(
const GAMGInterface& GAMGCp,
const lduInterfaceField& fineInterface
@ -88,6 +88,19 @@ public:
(GAMGCp, fineInterface)
);
declareRunTimeSelectionTable
(
autoPtr,
GAMGInterfaceField,
lduInterface,
(
const GAMGInterface& GAMGCp,
const bool doTransform,
const int rank
),
(GAMGCp, doTransform, rank)
);
// Selectors
@ -99,6 +112,15 @@ public:
const lduInterfaceField& fineInterface
);
//- Return a pointer to a new interface created on freestore given
// the fine interface
static autoPtr<GAMGInterfaceField> New
(
const GAMGInterface& GAMGCp,
const bool doTransform,
const int rank
);
// Constructors
@ -112,6 +134,30 @@ public:
lduInterfaceField(GAMGCp),
interface_(GAMGCp)
{}
//- Construct from GAMG interface and fine level interface field
GAMGInterfaceField
(
const GAMGInterface& GAMGCp,
const bool doTransform,
const int rank
)
:
lduInterfaceField(GAMGCp),
interface_(GAMGCp)
{}
// Member Functions
// Access
//- Return interface
const GAMGInterface& interface() const
{
return interface_;
}
};

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
@ -35,10 +35,10 @@ Foam::autoPtr<Foam::GAMGInterfaceField> Foam::GAMGInterfaceField::New
{
const word coupleType(fineInterface.interfaceFieldType());
lduInterfaceConstructorTable::iterator cstrIter =
lduInterfaceConstructorTablePtr_->find(coupleType);
lduInterfaceFieldConstructorTable::iterator cstrIter =
lduInterfaceFieldConstructorTablePtr_->find(coupleType);
if (cstrIter == lduInterfaceConstructorTablePtr_->end())
if (cstrIter == lduInterfaceFieldConstructorTablePtr_->end())
{
FatalErrorIn
(
@ -48,7 +48,7 @@ Foam::autoPtr<Foam::GAMGInterfaceField> Foam::GAMGInterfaceField::New
) << "Unknown GAMGInterfaceField type "
<< coupleType << nl
<< "Valid GAMGInterfaceField types are :"
<< lduInterfaceConstructorTablePtr_->sortedToc()
<< lduInterfaceFieldConstructorTablePtr_->sortedToc()
<< exit(FatalError);
}
@ -56,4 +56,33 @@ Foam::autoPtr<Foam::GAMGInterfaceField> Foam::GAMGInterfaceField::New
}
Foam::autoPtr<Foam::GAMGInterfaceField> Foam::GAMGInterfaceField::New
(
const GAMGInterface& GAMGCp,
const bool doTransform,
const int rank
)
{
const word coupleType(GAMGCp.type());
lduInterfaceConstructorTable::iterator cstrIter =
lduInterfaceConstructorTablePtr_->find(coupleType);
if (cstrIter == lduInterfaceConstructorTablePtr_->end())
{
FatalErrorIn
(
"GAMGInterfaceField::New"
"(const word&, const GAMGInterface&, const bool, const int)"
) << "Unknown GAMGInterfaceField type "
<< coupleType << nl
<< "Valid GAMGInterfaceField types are :"
<< lduInterfaceConstructorTablePtr_->sortedToc()
<< exit(FatalError);
}
return autoPtr<GAMGInterfaceField>(cstrIter()(GAMGCp, doTransform, rank));
}
// ************************************************************************* //

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
@ -38,6 +38,12 @@ namespace Foam
cyclicGAMGInterfaceField,
lduInterface
);
addToRunTimeSelectionTable
(
GAMGInterfaceField,
cyclicGAMGInterfaceField,
lduInterfaceField
);
}
@ -62,6 +68,20 @@ Foam::cyclicGAMGInterfaceField::cyclicGAMGInterfaceField
}
Foam::cyclicGAMGInterfaceField::cyclicGAMGInterfaceField
(
const GAMGInterface& GAMGCp,
const bool doTransform,
const int rank
)
:
GAMGInterfaceField(GAMGCp, doTransform, rank),
cyclicInterface_(refCast<const cyclicGAMGInterface>(GAMGCp)),
doTransform_(doTransform),
rank_(rank)
{}
// * * * * * * * * * * * * * * * * Desstructor * * * * * * * * * * * * * * * //
Foam::cyclicGAMGInterfaceField::~cyclicGAMGInterfaceField()

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
@ -89,6 +89,14 @@ public:
const lduInterfaceField& fineInterfaceField
);
//- Construct from GAMG interface and fine level interface field
cyclicGAMGInterfaceField
(
const GAMGInterface& GAMGCp,
const bool doTransform,
const int rank
);
//- Destructor
virtual ~cyclicGAMGInterfaceField();

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
@ -38,6 +38,12 @@ namespace Foam
processorCyclicGAMGInterfaceField,
lduInterface
);
addToRunTimeSelectionTable
(
GAMGInterfaceField,
processorCyclicGAMGInterfaceField,
lduInterfaceField
);
}
@ -53,6 +59,17 @@ Foam::processorCyclicGAMGInterfaceField::processorCyclicGAMGInterfaceField
{}
Foam::processorCyclicGAMGInterfaceField::processorCyclicGAMGInterfaceField
(
const GAMGInterface& GAMGCp,
const bool doTransform,
const int rank
)
:
processorGAMGInterfaceField(GAMGCp, doTransform, rank)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::processorCyclicGAMGInterfaceField::~processorCyclicGAMGInterfaceField()

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
@ -77,6 +77,14 @@ public:
const lduInterfaceField& fineInterface
);
//- Construct from GAMG interface and fine level interface field
processorCyclicGAMGInterfaceField
(
const GAMGInterface& GAMGCp,
const bool doTransform,
const int rank
);
// Destructor

View File

@ -38,6 +38,12 @@ namespace Foam
processorGAMGInterfaceField,
lduInterface
);
addToRunTimeSelectionTable
(
GAMGInterfaceField,
processorGAMGInterfaceField,
lduInterfaceField
);
}
@ -62,6 +68,20 @@ Foam::processorGAMGInterfaceField::processorGAMGInterfaceField
}
Foam::processorGAMGInterfaceField::processorGAMGInterfaceField
(
const GAMGInterface& GAMGCp,
const bool doTransform,
const int rank
)
:
GAMGInterfaceField(GAMGCp, doTransform, rank),
procInterface_(refCast<const processorGAMGInterface>(GAMGCp)),
doTransform_(doTransform),
rank_(rank)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::processorGAMGInterfaceField::~processorGAMGInterfaceField()
@ -79,6 +99,9 @@ void Foam::processorGAMGInterfaceField::initInterfaceMatrixUpdate
const Pstream::commsTypes commsType
) const
{
label oldWarn = UPstream::warnComm;
UPstream::warnComm = comm();
procInterface_.interfaceInternalField(psiInternal, scalarSendBuf_);
if (commsType == Pstream::nonBlocking && !Pstream::floatTransfer)
@ -92,7 +115,8 @@ void Foam::processorGAMGInterfaceField::initInterfaceMatrixUpdate
procInterface_.neighbProcNo(),
reinterpret_cast<char*>(scalarReceiveBuf_.begin()),
scalarReceiveBuf_.byteSize(),
procInterface_.tag()
procInterface_.tag(),
comm()
);
outstandingSendRequest_ = UPstream::nRequests();
@ -102,7 +126,8 @@ void Foam::processorGAMGInterfaceField::initInterfaceMatrixUpdate
procInterface_.neighbProcNo(),
reinterpret_cast<const char*>(scalarSendBuf_.begin()),
scalarSendBuf_.byteSize(),
procInterface_.tag()
procInterface_.tag(),
comm()
);
}
else
@ -111,6 +136,8 @@ void Foam::processorGAMGInterfaceField::initInterfaceMatrixUpdate
}
const_cast<processorGAMGInterfaceField&>(*this).updatedMatrix() = false;
UPstream::warnComm = oldWarn;
}
@ -128,6 +155,9 @@ void Foam::processorGAMGInterfaceField::updateInterfaceMatrix
return;
}
label oldWarn = UPstream::warnComm;
UPstream::warnComm = comm();
const labelUList& faceCells = procInterface_.faceCells();
if (commsType == Pstream::nonBlocking && !Pstream::floatTransfer)
@ -171,6 +201,8 @@ void Foam::processorGAMGInterfaceField::updateInterfaceMatrix
}
const_cast<processorGAMGInterfaceField&>(*this).updatedMatrix() = true;
UPstream::warnComm = oldWarn;
}

View File

@ -105,6 +105,14 @@ public:
const lduInterfaceField& fineInterface
);
//- Construct from GAMG interface and fine level interface field
processorGAMGInterfaceField
(
const GAMGInterface& GAMGCp,
const bool doTransform,
const int rank
);
//- Destructor
virtual ~processorGAMGInterfaceField();
@ -146,6 +154,12 @@ public:
//- Processor interface functions
//- Return communicator used for comms
virtual int comm() const
{
return procInterface_.comm();
}
//- Return processor number
virtual int myProcNo() const
{

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
@ -31,9 +31,26 @@ namespace Foam
{
defineTypeNameAndDebug(GAMGInterface, 0);
defineRunTimeSelectionTable(GAMGInterface, lduInterface);
defineRunTimeSelectionTable(GAMGInterface, Istream);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::GAMGInterface::GAMGInterface
(
const label index,
const lduInterfacePtrsList& coarseInterfaces,
Istream& is
)
:
index_(index),
coarseInterfaces_(coarseInterfaces),
faceCells_(is),
faceRestrictAddressing_(is)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::GAMGInterface::combine(const GAMGInterface& coarseGi)
@ -66,6 +83,27 @@ Foam::tmp<Foam::scalarField> Foam::GAMGInterface::agglomerateCoeffs
tmp<scalarField> tcoarseCoeffs(new scalarField(size(), 0.0));
scalarField& coarseCoeffs = tcoarseCoeffs();
if (fineCoeffs.size() != faceRestrictAddressing_.size())
{
FatalErrorIn
(
"GAMGInterface::agglomerateCoeffs(const scalarField&) const"
) << "Size of coefficients " << fineCoeffs.size()
<< " does not correspond to the size of the restriction "
<< faceRestrictAddressing_.size()
<< abort(FatalError);
}
if (debug && max(faceRestrictAddressing_) > size())
{
FatalErrorIn
(
"GAMGInterface::agglomerateCoeffs(const scalarField&) const"
) << "Face restrict addressing addresses outside of coarse interface"
<< " size. Max addressing:" << max(faceRestrictAddressing_)
<< " coarse size:" << size()
<< abort(FatalError);
}
forAll(faceRestrictAddressing_, ffi)
{
coarseCoeffs[faceRestrictAddressing_[ffi]] += fineCoeffs[ffi];
@ -75,4 +113,10 @@ Foam::tmp<Foam::scalarField> Foam::GAMGInterface::agglomerateCoeffs
}
void Foam::GAMGInterface::write(Ostream& os) const
{
os << faceCells_ << token::SPACE << faceRestrictAddressing_;
}
// ************************************************************************* //

Some files were not shown because too many files have changed in this diff Show More