Commit dd4a5564 authored by mattijs's avatar mattijs
Browse files

BUG: lduPrimitiveMesh: incorrect upper-triangular ordering

parent 488cfd59
......@@ -29,6 +29,14 @@ License
#include "procLduInterface.H"
#include "cyclicLduInterface.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(LUscalarMatrix, 0);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::LUscalarMatrix::LUscalarMatrix(const scalarSquareMatrix& matrix)
......@@ -134,8 +142,39 @@ Foam::LUscalarMatrix::LUscalarMatrix
convert(ldum, interfaceCoeffs, interfaces);
}
if (Pstream::master(comm_))
if (debug && Pstream::master(comm_))
{
label nRows = n();
label nColumns = m();
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_);
}
......
......@@ -87,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
......
......@@ -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
......
......@@ -270,12 +270,19 @@ void Foam::GAMGAgglomeration::agglomerateLduAddressing
interfaceLevels_[fineLevelIndex];
// Create coarse-level interfaces
primitiveInterfaces_.set
(
fineLevelIndex + 1,
new PtrList<const lduInterface>(fineInterfaces.size())
);
PtrList<const lduInterface>& coarsePrimInterfaces =
primitiveInterfaces_[fineLevelIndex + 1];
interfaceLevels_.set
(
fineLevelIndex + 1,
new lduInterfacePtrsList(fineInterfaces.size())
);
lduInterfacePtrsList& coarseInterfaces =
interfaceLevels_[fineLevelIndex + 1];
......@@ -317,7 +324,7 @@ void Foam::GAMGAgglomeration::agglomerateLduAddressing
{
if (fineInterfaces.set(inti))
{
coarseInterfaces.set
coarsePrimInterfaces.set
(
inti,
GAMGInterface::New
......@@ -336,6 +343,7 @@ void Foam::GAMGAgglomeration::agglomerateLduAddressing
).ptr()
);
coarseInterfaces.set(inti, &coarsePrimInterfaces[inti]);
coarseInterfaceAddr[inti] = coarseInterfaces[inti].faceCells();
nPatchFaces[inti] = coarseInterfaceAddr[inti].size();
patchFineToCoarse[inti] = refCast<const GAMGInterface>
......@@ -408,13 +416,13 @@ void Foam::GAMGAgglomeration::gatherMeshes
// Combine all addressing
// ~~~~~~~~~~~~~~~~~~~~~~
Pout<< "Own Mesh " << myMesh.lduAddr().size() << endl;
forAll(otherMeshes, i)
{
Pout<< " otherMesh " << i << " "
<< otherMeshes[i].lduAddr().size()
<< endl;
}
//Pout<< "Own Mesh " << myMesh.lduAddr().size() << endl;
//forAll(otherMeshes, i)
//{
// Pout<< " otherMesh " << i << " "
// << otherMeshes[i].lduAddr().size()
// << endl;
//}
labelList procFaceOffsets;
......@@ -437,7 +445,7 @@ void Foam::GAMGAgglomeration::gatherMeshes
)
);
Pout<< "** Agglomerated Mesh " << allMesh().info();
//Pout<< "** Agglomerated Mesh " << allMesh().info();
}
UPstream::warnComm = oldWarn;
......@@ -459,7 +467,6 @@ void Foam::GAMGAgglomeration::procAgglomerateLduAddressing
<< " havefinemesh:" << meshLevels_.set(levelIndex-1) << endl;
const lduMesh& myMesh = meshLevels_[levelIndex-1];
//label meshComm = myMesh.comm();
Pout<< "my mesh :" << myMesh.info();
Pout<< "nCoarseCells :" << nCells_[levelIndex] << endl;
......@@ -510,8 +517,6 @@ void Foam::GAMGAgglomeration::procAgglomerateLduAddressing
// 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));
//procRestrictAddressing_.set(levelIndex, new labelField(0));
//procFaceRestrictAddressing_.set(levelIndex, new labelList(0));
procFaceMap_.set(levelIndex, new labelListList(0));
procBoundaryMap_.set(levelIndex, new labelListList(0));
procBoundaryFaceMap_.set(levelIndex, new labelListListList(0));
......@@ -575,46 +580,16 @@ void Foam::GAMGAgglomeration::procAgglomerateLduAddressing
}
else
{
const lduPrimitiveMesh& levelMesh = meshLevels_[levelIndex-1];
Pout<< "** DONE GAMGAgglomeration:procAgglomerateLduAddressing **"
<< endl;
Pout<< "my mesh:" << levelMesh.info();
Pout<< "nCoarseCells:" << nCells_[levelIndex] << endl;
Pout<< "restrictAddressing_:" << restrictAddressing_[levelIndex].size()
<< " max:" << max(restrictAddressing_[levelIndex]) << endl;
//Pout<< "faceRestrictAddressing_:"
// << faceRestrictAddressing_[levelIndex].size()
// << " max:" << max(faceRestrictAddressing_[levelIndex]) << endl;
//Pout<< "patchFaceRestrictAddressing_:" << endl;
//forAll(patchFaceRestrictAddressing_[levelIndex], patchI)
//{
// const labelList& map =
// patchFaceRestrictAddressing_[levelIndex][patchI];
//
// if (map.size())
// {
// Pout<< " patch:" << patchI
// << " size:" << map.size()
// << " max:" << max(map)
// << endl;
// }
//}
//Pout<< "interfaceLevels_:" << endl;
//forAll(interfaceLevels_[levelIndex], patchI)
//{
// if (interfaceLevels_[levelIndex].set(patchI))
// {
// Pout<< " patch:" << patchI
// << " interface:"
// << interfaceLevels_[levelIndex][patchI].type()
// << " size:"
// << interfaceLevels_[levelIndex][patchI].faceCells().size()
// << endl;
// }
//}
Pout<< "** DONE GAMGAgglomeration:procAgglomerateLduAddressing **"
<< endl;
//const lduPrimitiveMesh& levelMesh = meshLevels_[levelIndex-1];
//Pout<< "** DONE GAMGAgglomeration:procAgglomerateLduAddressing **"
// << endl;
//Pout<< "my mesh:" << levelMesh.info();
//Pout<< "nCoarseCells:" << nCells_[levelIndex] << endl;
//Pout<< "restrictAddressing_:"
// << restrictAddressing_[levelIndex].size()
// << " max:" << max(restrictAddressing_[levelIndex]) << endl;
//Pout<< "** DONE GAMGAgglomeration:procAgglomerateLduAddressing **"
// << endl;
}
UPstream::warnComm = oldWarn;
......@@ -628,14 +603,15 @@ void Foam::GAMGAgglomeration::procAgglomerateRestrictAddressing
const label levelIndex
)
{
Pout<< "** GAMGAgglomeration:procAgglomerateRestrictAddressing **" << endl;
Pout<< "level:" << levelIndex << endl;
Pout<< "procIDs:" << procIDs << endl;
Pout<< "myProcNo:" << UPstream::myProcNo(comm) << endl;
//Pout<< "** GAMGAgglomeration:procAgglomerateRestrictAddressing **"
// << endl;
//Pout<< "level:" << levelIndex << endl;
//Pout<< "procIDs:" << procIDs << endl;
//Pout<< "myProcNo:" << UPstream::myProcNo(comm) << endl;
//Pout<< "procCellOffsets_:" << procCellOffsets_[levelIndex] << endl;
//const lduMesh& levelMesh = meshLevel(levelIndex);
Pout<< "fine:" << restrictAddressing_[levelIndex].size() << endl;
//Pout<< "fine:" << restrictAddressing_[levelIndex].size() << endl;
// Collect number of cells
labelList nFineCells;
......@@ -648,7 +624,7 @@ void Foam::GAMGAgglomeration::procAgglomerateRestrictAddressing
);
//const labelList& offsets = procCellOffsets_[levelIndex];
//Pout<< "offsets:" << offsets << endl;
Pout<< "nFineCells:" << nFineCells << endl;
//Pout<< "nFineCells:" << nFineCells << endl;
labelList offsets(nFineCells.size()+1);
{
......@@ -658,10 +634,10 @@ void Foam::GAMGAgglomeration::procAgglomerateRestrictAddressing
offsets[i+1] = offsets[i] + nFineCells[i];
}
}
Pout<< "offsets:" << offsets << endl;
//Pout<< "offsets:" << offsets << endl;
// Combine and renumber nCoarseCells
Pout<< "Starting nCells_ ..." << endl;
//Pout<< "Starting nCells_ ..." << endl;
labelList nCoarseCells;
gatherList
(
......@@ -672,7 +648,7 @@ void Foam::GAMGAgglomeration::procAgglomerateRestrictAddressing
);
// (cell)restrictAddressing
Pout<< "Starting restrictAddressing_ ..." << endl;
//Pout<< "Starting restrictAddressing_ ..." << endl;
const globalIndex cellOffsetter(offsets);
labelList procRestrictAddressing;
......@@ -695,12 +671,13 @@ void Foam::GAMGAgglomeration::procAgglomerateRestrictAddressing
coarseCellOffsets[i+1] = coarseCellOffsets[i]+nCoarseCells[i];
}
}
label nOldCoarseCells = nCells_[levelIndex];
//label nOldCoarseCells = nCells_[levelIndex];
nCells_[levelIndex] = coarseCellOffsets.last();
Pout<< "Finished nCoarseCells_ ..."
<< " was:" << nOldCoarseCells
<< " now:" << nCells_[levelIndex] << endl;
//Pout<< "Finished nCoarseCells_ ..."
// << " was:" << nOldCoarseCells
// << " now:" << nCells_[levelIndex] << endl;
// Renumber consecutively
......@@ -718,17 +695,17 @@ void Foam::GAMGAgglomeration::procAgglomerateRestrictAddressing
}
}
Pout<< "coarseCellOffsets:" << coarseCellOffsets << endl;
//Pout<< "coarseCellOffsets:" << coarseCellOffsets << endl;
restrictAddressing_[levelIndex].transfer(procRestrictAddressing);
Pout<< "restrictAddressing_:"
<< restrictAddressing_[levelIndex].size()
<< " min:" << min(restrictAddressing_[levelIndex])
<< " max:" << max(restrictAddressing_[levelIndex])
<< endl;
Pout<< "Finished restrictAddressing_ ..." << endl;
//Pout<< "restrictAddressing_:"
// << restrictAddressing_[levelIndex].size()
// << " min:" << min(restrictAddressing_[levelIndex])
// << " max:" << max(restrictAddressing_[levelIndex])
// << endl;
//Pout<< "Finished restrictAddressing_ ..." << endl;
}
Pout<< "** DONE GAMGAgglomeration:procAgglomerateRestrictAddressing **"
<< endl;
//Pout<< "** DONE GAMGAgglomeration:procAgglomerateRestrictAddressing **"
// << endl;
}
......@@ -777,4 +754,55 @@ void Foam::GAMGAgglomeration::procAgglomerateRestrictAddressing
//}
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]);
}
// ************************************************************************* //
......@@ -52,6 +52,7 @@ void Foam::GAMGAgglomeration::compactLevels(const label nCreatedLevels)
nPatchFaces_.setSize(nCreatedLevels);
patchFaceRestrictAddressing_.setSize(nCreatedLevels);
meshLevels_.setSize(nCreatedLevels);
primitiveInterfaces_.setSize(nCreatedLevels + 1);
interfaceLevels_.setSize(nCreatedLevels + 1);
// Have procCommunicator_ always, even if not procAgglomerating
......@@ -62,9 +63,6 @@ void Foam::GAMGAgglomeration::compactLevels(const label nCreatedLevels)
<< " to " << nCreatedLevels << " levels" << endl;
procAgglomMap_.setSize(nCreatedLevels);
agglomProcIDs_.setSize(nCreatedLevels);
//procMeshLevels_.setSize(nCreatedLevels);
//procRestrictAddressing_.setSize(nCreatedLevels);
//procFaceRestrictAddressing_.setSize(nCreatedLevels);
procCellOffsets_.setSize(nCreatedLevels);
procFaceMap_.setSize(nCreatedLevels);
procBoundaryMap_.setSize(nCreatedLevels);
......@@ -165,60 +163,34 @@ void Foam::GAMGAgglomeration::compactLevels(const label nCreatedLevels)
// Processor restriction map: per processor the coarse processor
labelList procAgglomMap(UPstream::nProcs(levelComm));
// Master processor
labelList masterProcs;
// Local processors that agglomerate. agglomProcIDs[0] is in
// masterProc.
List<int> agglomProcIDs;
{
procAgglomMap[0] = 0;
procAgglomMap[1] = 0;
procAgglomMap[2] = 1;
//procAgglomMap[3] = 1;
// Determine the master processors
Map<label> agglomToMaster(procAgglomMap.size());
forAll(procAgglomMap, procI)
label half = (procAgglomMap.size()+1)/2;
for (label i = 0; i < half; i++)
{
label coarseI = procAgglomMap[procI];
Map<label>::iterator fnd = agglomToMaster.find(coarseI);
if (fnd == agglomToMaster.end())
{
agglomToMaster.insert(coarseI, procI);
}
else
{
fnd() = max(fnd(), procI);
}
procAgglomMap[i] = 0;
}
masterProcs.setSize(agglomToMaster.size());
forAllConstIter(Map<label>, agglomToMaster, iter)
for (label i = half; i < procAgglomMap.size(); i++)
{
masterProcs[iter.key()] = iter();
procAgglomMap[i] = 1;
}
}
// Master processor
labelList masterProcs;
// Local processors that agglomerate. agglomProcIDs[0] is in
// masterProc.
List<int> agglomProcIDs;
calculateRegionMaster
(
levelComm,
procAgglomMap,
masterProcs,
agglomProcIDs
);
// Collect all the processors in my agglomeration
label myProcID = Pstream::myProcNo(levelComm);
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]);
}
Pout<< "procAgglomMap:" << procAgglomMap << endl;
Pout<< "agglomProcIDs:" << agglomProcIDs << endl;
Pout<< "masterProcs:" << masterProcs << endl;
// Allocate a communicator for the processor-agglomerated matrix
......@@ -259,10 +231,6 @@ void Foam::GAMGAgglomeration::compactLevels(const label nCreatedLevels)
Pout<< "Starting procAgglomerateRestrictAddressing level:"
<< levelI << endl;
//procAgglomMap_.set(levelI, new labelList(procAgglomMap));
//agglomProcIDs_.set(levelI, new labelList(agglomProcIDs));
//procCommunicator_[levelI] = procAgglomComm;
procAgglomerateRestrictAddressing
(
levelComm,
......@@ -335,6 +303,9 @@ void Foam::GAMGAgglomeration::compactLevels(const label nCreatedLevels)
<< endl;
}
}
Pout<< fineMesh.info() << endl;
Pout<< endl;
}
else
......@@ -445,6 +416,7 @@ Foam::GAMGAgglomeration::GAMGAgglomeration
patchFaceRestrictAddressing_(maxLevels_),
meshLevels_(maxLevels_),
primitiveInterfaces_(maxLevels_ + 1),
interfaceLevels_(maxLevels_ + 1)
{
procCommunicator_.setSize(maxLevels_ + 1, -1);
......@@ -454,9 +426,6 @@ Foam::GAMGAgglomeration::GAMGAgglomeration
<< " levels" << endl;
procAgglomMap_.setSize(maxLevels_);
agglomProcIDs_.setSize(maxLevels_);
//procMeshLevels_.setSize(maxLevels_);
//procRestrictAddressing_.setSize(maxLevels_);
//procFaceRestrictAddressing_.setSize(maxLevels_);
procCellOffsets_.setSize(maxLevels_);
procFaceMap_.setSize(maxLevels_);
procBoundaryMap_.setSize(maxLevels_);
......@@ -573,42 +542,6 @@ const Foam::GAMGAgglomeration& Foam::GAMGAgglomeration::New
Foam::GAMGAgglomeration::~GAMGAgglomeration()
{
// // Temporary store the user-defined communicators so we can delete them
// labelHashSet communicators(meshLevels_.size());
// forAll(meshLevels_, leveli)
// {
// communicators.insert(meshLevels_[leveli].comm());
// }
//forAll(procMeshLevels_, leveli)
//{
// if (procMeshLevels_.set(leveli))
// {
// communicators.insert(procMeshLevels_[leveli].comm());
// }
//}
// Pout<< "~GAMGAgglomeration() : current communicators:" << communicators
// << endl;
// 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);
}
}
}
// forAllConstIter(labelHashSet, communicators, iter)
// {
// UPstream::freeCommunicator(iter.key());
// }
forAllReverse(procCommunicator_, i)
{