zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

ENH: nonBlockingGaussSmoother: block as late as possible.

Browse Source 
renumberMesh: option for sorting proc cells last to benefit from this.
This commit is contained in:
mattijs
2011-12-20 09:09:51 +00:00
parent 97e4abdc87
commit f80b3a793a
5 changed files with 524 additions and 52 deletions
Download Patch File Download Diff File Expand all files Collapse all files

172
applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
View File

@ -464,23 +464,6 @@ int main(int argc, char *argv[])
"Renumber mesh to minimise bandwidth"
);
argList::addOption
(
"blockSize",
"block size",
"order cells into blocks (using decomposition) before ordering"
);
argList::addBoolOption
(
"orderPoints",
"order points into internal and boundary points"
);
argList::addBoolOption
(
"writeMaps",
"write cellMap, faceMap, pointMap in polyMesh/"
);
# include "addRegionOption.H"
# include "addOverwriteOption.H"
# include "addTimeOptions.H"
@ -507,41 +490,6 @@ int main(int argc, char *argv[])
const bool readDict = args.optionFound("dict");
label blockSize = 0;
args.optionReadIfPresent("blockSize", blockSize, 0);
if (blockSize > 0)
{
Info<< "Ordering cells into regions of size " << blockSize
<< " (using decomposition);"
<< " ordering faces into region-internal and region-external." << nl
<< endl;
if (blockSize < 0 || blockSize >= mesh.nCells())
{
FatalErrorIn(args.executable())
<< "Block size " << blockSize << " should be positive integer"
<< " and less than the number of cells in the mesh."
<< exit(FatalError);
}
}
const bool orderPoints = args.optionFound("orderPoints");
if (orderPoints)
{
Info<< "Ordering points into internal and boundary points." << nl
<< endl;
}
const bool writeMaps = args.optionFound("writeMaps");
if (writeMaps)
{
Info<< "Writing renumber maps (new to old) to polyMesh." << nl
<< endl;
}
const bool overwrite = args.optionFound("overwrite");
label band = getBand(mesh.faceOwner(), mesh.faceNeighbour());
@ -551,6 +499,11 @@ int main(int argc, char *argv[])
<< returnReduce(band, maxOp<label>()) << nl << endl;
bool sortCoupledFaceCells = false;
bool writeMaps = false;
bool orderPoints = false;
label blockSize = 0;
// Construct renumberMethod
autoPtr<renumberMethod> renumberPtr;
@ -571,6 +524,51 @@ int main(int argc, char *argv[])
);
renumberPtr = renumberMethod::New(renumberDict);
sortCoupledFaceCells = renumberDict.lookupOrDefault
(
"sortCoupledFaceCells",
false
);
if (sortCoupledFaceCells)
{
Info<< "Sorting cells on coupled boundaries to be last." << nl
<< endl;
}
blockSize = renumberDict.lookupOrDefault("blockSize", 0);
if (blockSize > 0)
{
Info<< "Ordering cells into regions of size " << blockSize
<< " (using decomposition);"
<< " ordering faces into region-internal and region-external."
<< nl << endl;
if (blockSize < 0 || blockSize >= mesh.nCells())
{
FatalErrorIn(args.executable())
<< "Block size " << blockSize
<< " should be positive integer"
<< " and less than the number of cells in the mesh."
<< exit(FatalError);
}
}
orderPoints = renumberDict.lookupOrDefault("orderPoints", false);
if (orderPoints)
{
Info<< "Ordering points into internal and boundary points." << nl
<< endl;
}
writeMaps = readLabel(renumberDict.lookup("writeMaps"));
if (writeMaps)
{
Info<< "Writing renumber maps (new to old) to polyMesh." << nl
<< endl;
}
}
else
{
@ -765,6 +763,76 @@ int main(int argc, char *argv[])
cellOrder = invert(mesh.nCells(), reverseCellOrder);
if (sortCoupledFaceCells)
{
// Change order so all coupled patch faceCells are at the end.
const polyBoundaryMesh& pbm = mesh.boundaryMesh();
// Collect all boundary cells on coupled patches
label nBndCells = 0;
forAll(pbm, patchI)
{
if (pbm[patchI].coupled())
{
nBndCells += pbm[patchI].size();
}
}
labelList bndCellMap(nBndCells);
labelList bndCells(bndCellMap);
nBndCells = 0;
forAll(pbm, patchI)
{
if (pbm[patchI].coupled())
{
const labelUList& faceCells = pbm[patchI].faceCells();
forAll(faceCells, i)
{
label cellI = faceCells[i];
bndCells[nBndCells] = cellI;
bndCellMap[nBndCells++] = reverseCellOrder[cellI];
}
}
}
bndCells.setSize(nBndCells);
bndCellMap.setSize(nBndCells);
// Sort
labelList order;
sortedOrder(bndCellMap, order);
// Redo newReverseCellOrder
labelList newReverseCellOrder(mesh.nCells(), -1);
label sortedI = mesh.nCells();
forAllReverse(order, i)
{
label origCellI = bndCells[order[i]];
newReverseCellOrder[origCellI] = --sortedI;
}
Info<< "Ordered all " << nBndCells << " cells with a coupled face"
<< " to the end of the cell list, starting at " << sortedI
<< endl;
// Compact
sortedI = 0;
forAll(cellOrder, newCellI)
{
label origCellI = cellOrder[newCellI];
if (newReverseCellOrder[origCellI] == -1)
{
newReverseCellOrder[origCellI] = sortedI++;
}
}
// Update sorted back to original (unsorted) map
cellOrder = invert(mesh.nCells(), newReverseCellOrder);
}
// Determine new to old face order with new cell numbering
faceOrder = getFaceOrder
(

17
applications/utilities/mesh/manipulation/renumberMesh/renumberMeshDict
View File

@ -15,6 +15,23 @@ FoamFile
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Write maps from renumbered back to original mesh
writeMaps true;
// Optional entry: sort cells on coupled boundaries to last for use with
// e.g. nonBlockingGaussSeidel.
sortCoupledFaceCells false;
// Optional entry: renumber on a block-by-block basis. This can be used on
// large cases to keep the blocks fitting in cache with all the the cache
// missed bunched at the end.
//blockSize 0;
// Optional entry: sort points into internal and boundary points
//orderPoints false;
method CuthillMcKee;
//method manual;
//method random;

1
src/OpenFOAM/Make/files
View File

@ -250,6 +250,7 @@ $(lduMatrix)/solvers/ICCG/ICCG.C
$(lduMatrix)/solvers/BICCG/BICCG.C
$(lduMatrix)/smoothers/GaussSeidel/GaussSeidelSmoother.C
$(lduMatrix)/smoothers/nonBlockingGaussSeidel/nonBlockingGaussSeidelSmoother.C
$(lduMatrix)/smoothers/DIC/DICSmoother.C
$(lduMatrix)/smoothers/DICGaussSeidel/DICGaussSeidelSmoother.C
$(lduMatrix)/smoothers/DILU/DILUSmoother.C

266
src/OpenFOAM/matrices/lduMatrix/smoothers/nonBlockingGaussSeidel/nonBlockingGaussSeidelSmoother.C Normal file
View File

@ -0,0 +1,266 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 "nonBlockingGaussSeidelSmoother.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(nonBlockingGaussSeidelSmoother, 0);
lduMatrix::smoother::
addsymMatrixConstructorToTable<nonBlockingGaussSeidelSmoother>
addnonBlockingGaussSeidelSmootherSymMatrixConstructorToTable_;
lduMatrix::smoother::
addasymMatrixConstructorToTable<nonBlockingGaussSeidelSmoother>
addnonBlockingGaussSeidelSmootherAsymMatrixConstructorToTable_;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::nonBlockingGaussSeidelSmoother::nonBlockingGaussSeidelSmoother
(
const word& fieldName,
const lduMatrix& matrix,
const FieldField<Field, scalar>& interfaceBouCoeffs,
const FieldField<Field, scalar>& interfaceIntCoeffs,
const lduInterfaceFieldPtrsList& interfaces
)
:
lduMatrix::smoother
(
fieldName,
matrix,
interfaceBouCoeffs,
interfaceIntCoeffs,
interfaces
)
{
// Check that all interface addressing is sorted to be after the
// non-interface addressing.
register const label nCells = matrix.diag().size();
blockStart_ = nCells;
labelList startCellI(interfaceBouCoeffs.size(), -1);
forAll(interfaces, patchi)
{
if (interfaces.set(patchi))
{
const labelUList& faceCells = matrix_.lduAddr().patchAddr(patchi);
blockStart_ = min(blockStart_, min(faceCells));
}
}
if (debug)
{
Pout<< "nonBlockingGaussSeidelSmoother :"
<< " Starting block on cell " << blockStart_
<< " out of " << nCells << endl;
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::nonBlockingGaussSeidelSmoother::smooth
(
const word& fieldName_,
scalarField& psi,
const lduMatrix& matrix_,
const label blockStart,
const scalarField& source,
const FieldField<Field, scalar>& interfaceBouCoeffs_,
const lduInterfaceFieldPtrsList& interfaces_,
const direction cmpt,
const label nSweeps
)
{
register scalar* __restrict__ psiPtr = psi.begin();
register const label nCells = psi.size();
scalarField bPrime(nCells);
register scalar* __restrict__ bPrimePtr = bPrime.begin();
register const scalar* const __restrict__ diagPtr = matrix_.diag().begin();
register const scalar* const __restrict__ upperPtr =
matrix_.upper().begin();
register const scalar* const __restrict__ lowerPtr =
matrix_.lower().begin();
register const label* const __restrict__ uPtr =
matrix_.lduAddr().upperAddr().begin();
register const label* const __restrict__ ownStartPtr =
matrix_.lduAddr().ownerStartAddr().begin();
// Parallel boundary initialisation. The parallel boundary is treated
// as an effective jacobi interface in the boundary.
// Note: there is a change of sign in the coupled
// interface update. The reason for this is that the
// internal coefficients are all located at the l.h.s. of
// the matrix whereas the "implicit" coefficients on the
// coupled boundaries are all created as if the
// coefficient contribution is of a source-kind (i.e. they
// have a sign as if they are on the r.h.s. of the matrix.
// To compensate for this, it is necessary to turn the
// sign of the contribution.
FieldField<Field, scalar>& mBouCoeffs =
const_cast<FieldField<Field, scalar>&>
(
interfaceBouCoeffs_
);
forAll(mBouCoeffs, patchi)
{
if (interfaces_.set(patchi))
{
mBouCoeffs[patchi].negate();
}
}
for (label sweep=0; sweep<nSweeps; sweep++)
{
bPrime = source;
matrix_.initMatrixInterfaces
(
mBouCoeffs,
interfaces_,
psi,
bPrime,
cmpt
);
register scalar curPsi;
register label fStart;
register label fEnd = ownStartPtr[0];
for (register label cellI=0; cellI<blockStart; cellI++)
{
// Start and end of this row
fStart = fEnd;
fEnd = ownStartPtr[cellI + 1];
// Get the accumulated neighbour side
curPsi = bPrimePtr[cellI];
// Accumulate the owner product side
for (register label curFace=fStart; curFace<fEnd; curFace++)
{
curPsi -= upperPtr[curFace]*psiPtr[uPtr[curFace]];
}
// Finish current psi
curPsi /= diagPtr[cellI];
// Distribute the neighbour side using current psi
for (register label curFace=fStart; curFace<fEnd; curFace++)
{
bPrimePtr[uPtr[curFace]] -= lowerPtr[curFace]*curPsi;
}
psiPtr[cellI] = curPsi;
}
matrix_.updateMatrixInterfaces
(
mBouCoeffs,
interfaces_,
psi,
bPrime,
cmpt
);
// Update rest of the cells
for (label cellI=blockStart; cellI < nCells; cellI++)
{
// Start and end of this row
fStart = fEnd;
fEnd = ownStartPtr[cellI + 1];
// Get the accumulated neighbour side
curPsi = bPrimePtr[cellI];
// Accumulate the owner product side
for (register label curFace=fStart; curFace<fEnd; curFace++)
{
curPsi -= upperPtr[curFace]*psiPtr[uPtr[curFace]];
}
// Finish current psi
curPsi /= diagPtr[cellI];
// Distribute the neighbour side using current psi
for (register label curFace=fStart; curFace<fEnd; curFace++)
{
bPrimePtr[uPtr[curFace]] -= lowerPtr[curFace]*curPsi;
}
psiPtr[cellI] = curPsi;
}
}
// Restore interfaceBouCoeffs_
forAll(mBouCoeffs, patchi)
{
if (interfaces_.set(patchi))
{
mBouCoeffs[patchi].negate();
}
}
}
void Foam::nonBlockingGaussSeidelSmoother::smooth
(
scalarField& psi,
const scalarField& source,
const direction cmpt,
const label nSweeps
) const
{
smooth
(
fieldName_,
psi,
matrix_,
blockStart_,
source,
interfaceBouCoeffs_,
interfaces_,
cmpt,
nSweeps
);
}
// ************************************************************************* //

120
src/OpenFOAM/matrices/lduMatrix/smoothers/nonBlockingGaussSeidel/nonBlockingGaussSeidelSmoother.H Normal file
View File

@ -0,0 +1,120 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::nonBlockingGaussSeidelSmoother
Description
Variant of gaussSeidelSmoother that expects processor boundary
cells to be sorted last and so can block later. Only when the
cells are actually visited does it need the results to be present.
It is expected that there is little benefit to be gained from doing
this on a patch by patch basis since the number of processor interfaces
is quite small and the overhead of checking whether a processor interface
is finished might be quite high (call into mpi). Also this would
require a dynamic memory allocation to store the state of the outstanding
requests.
SourceFiles
nonBlockingGaussSeidelSmoother.C
\*---------------------------------------------------------------------------*/
#ifndef nonBlockingGaussSeidelSmoother_H
#define nonBlockingGaussSeidelSmoother_H
#include "lduMatrix.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class nonBlockingGaussSeidelSmoother Declaration
\*---------------------------------------------------------------------------*/
class nonBlockingGaussSeidelSmoother
:
public lduMatrix::smoother
{
// Private data
//- Starting cell when to block
label blockStart_;
public:
//- Runtime type information
TypeName("nonBlockingGaussSeidel");
// Constructors
//- Construct from components
nonBlockingGaussSeidelSmoother
(
const word& fieldName,
const lduMatrix& matrix,
const FieldField<Field, scalar>& interfaceBouCoeffs,
const FieldField<Field, scalar>& interfaceIntCoeffs,
const lduInterfaceFieldPtrsList& interfaces
);
// Member Functions
//- Smooth for the given number of sweeps
static void smooth
(
const word& fieldName,
scalarField& psi,
const lduMatrix& matrix,
const label blockStart,
const scalarField& source,
const FieldField<Field, scalar>& interfaceBouCoeffs,
const lduInterfaceFieldPtrsList& interfaces,
const direction cmpt,
const label nSweeps
);
//- Smooth the solution for a given number of sweeps
virtual void smooth
(
scalarField& psi,
const scalarField& Source,
const direction cmpt,
const label nSweeps
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //