From f80b3a793abec398078b7500873f816f484e8edf Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 20 Dec 2011 09:09:51 +0000
Subject: [PATCH] ENH: nonBlockingGaussSmoother: block as late as possible.

renumberMesh: option for sorting proc cells last to benefit from this.
---
 .../manipulation/renumberMesh/renumberMesh.C  | 172 +++++++----
 .../renumberMesh/renumberMeshDict             |  17 ++
 src/OpenFOAM/Make/files                       |   1 +
 .../nonBlockingGaussSeidelSmoother.C          | 266 ++++++++++++++++++
 .../nonBlockingGaussSeidelSmoother.H          | 120 ++++++++
 5 files changed, 524 insertions(+), 52 deletions(-)
 create mode 100644 src/OpenFOAM/matrices/lduMatrix/smoothers/nonBlockingGaussSeidel/nonBlockingGaussSeidelSmoother.C
 create mode 100644 src/OpenFOAM/matrices/lduMatrix/smoothers/nonBlockingGaussSeidel/nonBlockingGaussSeidelSmoother.H

diff --git a/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C b/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
index a8c18c0712..7e090fa30e 100644
--- a/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
+++ b/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
@@ -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
         (
diff --git a/applications/utilities/mesh/manipulation/renumberMesh/renumberMeshDict b/applications/utilities/mesh/manipulation/renumberMesh/renumberMeshDict
index 7ff33fc530..a852fad7db 100644
--- a/applications/utilities/mesh/manipulation/renumberMesh/renumberMeshDict
+++ b/applications/utilities/mesh/manipulation/renumberMesh/renumberMeshDict
@@ -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;
diff --git a/src/OpenFOAM/Make/files b/src/OpenFOAM/Make/files
index cbd9888573..4bab8ad32a 100644
--- a/src/OpenFOAM/Make/files
+++ b/src/OpenFOAM/Make/files
@@ -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
diff --git a/src/OpenFOAM/matrices/lduMatrix/smoothers/nonBlockingGaussSeidel/nonBlockingGaussSeidelSmoother.C b/src/OpenFOAM/matrices/lduMatrix/smoothers/nonBlockingGaussSeidel/nonBlockingGaussSeidelSmoother.C
new file mode 100644
index 0000000000..77c1ba5430
--- /dev/null
+++ b/src/OpenFOAM/matrices/lduMatrix/smoothers/nonBlockingGaussSeidel/nonBlockingGaussSeidelSmoother.C
@@ -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
+    );
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/lduMatrix/smoothers/nonBlockingGaussSeidel/nonBlockingGaussSeidelSmoother.H b/src/OpenFOAM/matrices/lduMatrix/smoothers/nonBlockingGaussSeidel/nonBlockingGaussSeidelSmoother.H
new file mode 100644
index 0000000000..0562d875cf
--- /dev/null
+++ b/src/OpenFOAM/matrices/lduMatrix/smoothers/nonBlockingGaussSeidel/nonBlockingGaussSeidelSmoother.H
@@ -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
+
+// ************************************************************************* //