diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.C
index 3662f1c354..18b0f90269 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.C
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.C
@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "GAMGSolver.H"
+#include "GAMGInterface.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -74,20 +75,212 @@ Foam::GAMGSolver::GAMGSolver
     interpolateCorrection_(false),
     scaleCorrection_(matrix.symmetric()),
     directSolveCoarsest_(false),
-    processorAgglomerate_(false),
     agglomeration_(GAMGAgglomeration::New(matrix_, controlDict_)),
 
     matrixLevels_(agglomeration_.size()),
     interfaceLevels_(agglomeration_.size()),
+    primitiveInterfaceLevels_(agglomeration_.size()),
     interfaceLevelsBouCoeffs_(agglomeration_.size()),
     interfaceLevelsIntCoeffs_(agglomeration_.size())
 {
     readControls();
 
-    forAll(agglomeration_, fineLevelIndex)
+    //if (agglomeration_.processorAgglomerate())
+    //{
+    //    procMatrixLevels_.setSize(agglomeration_.size());
+    //    procInterfaceLevels_.setSize(agglomeration_.size());
+    //    procPrimitiveInterfaces_.setSize(agglomeration_.size());
+    //    procInterfaceLevelsBouCoeffs_.setSize(agglomeration_.size());
+    //    procInterfaceLevelsIntCoeffs_.setSize(agglomeration_.size());
+    //}
+
+    if (agglomeration_.processorAgglomerate())
     {
-        agglomerateMatrix(fineLevelIndex);
+        forAll(agglomeration_, fineLevelIndex)
+        {
+            if (agglomeration_.hasMeshLevel(fineLevelIndex))
+            {
+                Pout<< "Level:" << fineLevelIndex
+                    << " agglomerating matrix." << endl;
+
+                if
+                (
+                    //   !agglomeration_.hasMeshLevel(fineLevelIndex+1)
+                    //||  (
+                    //        agglomeration_.nCells(fineLevelIndex)
+                    //     != agglomeration_.meshLevel(fineLevelIndex+1).
+                    //            lduAddr().size()
+                    //    )
+
+                    (fineLevelIndex+1) < agglomeration_.procBoundaryMap_.size()
+                 && agglomeration_.procBoundaryMap_.set(fineLevelIndex+1)
+                )
+                {
+                    Pout<< "Level:" << fineLevelIndex
+                        << " agglomerating onto dummy coarse mesh." << endl;
+
+                    // 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);
+
+                    Pout<< "Level:" << fineLevelIndex
+                        << " agglomerating onto procIDs:" << procIDs << endl;
+
+                    procAgglomerateMatrix
+                    (
+                        procAgglomMap,
+                        procIDs,
+                        fineLevelIndex
+                    );
+
+                    Pout<< "Level:" << fineLevelIndex
+                        << " DONE agglomerating onto procIDs:" << procIDs
+                        << endl;
+                }
+                else
+                {
+                    Pout<< "Level:" << fineLevelIndex
+                        << " agglomerating onto coarse mesh at level "
+                        << fineLevelIndex + 1 << endl;
+
+                    agglomerateMatrix
+                    (
+                        fineLevelIndex,
+                        agglomeration_.meshLevel(fineLevelIndex + 1),
+                        agglomeration_.interfaceLevel(fineLevelIndex + 1)
+                    );
+                    Pout<< "Level:" << fineLevelIndex
+                        << " DONE agglomerating onto coarse mesh at level "
+                        << fineLevelIndex + 1 << endl;
+
+                }
+            }
+            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())
     {
@@ -119,117 +312,118 @@ Foam::GAMGSolver::GAMGSolver
 
             UPstream::warnComm = oldWarn;
         }
-        else if (processorAgglomerate_)
+        else if (agglomeration_.processorAgglomerate())
         {
-            // Pick a level to processor agglomerate
-            label agglomLevel = matrixLevels_.size() - 1;//1;
-
-
-            // Get mesh and matrix at this level
-            const lduMatrix& levelMatrix = matrixLevels_[agglomLevel];
-            const lduMesh& levelMesh = levelMatrix.mesh();
-
-
-            label levelComm = levelMesh.comm();
-            label oldWarn = UPstream::warnComm;
-            UPstream::warnComm = levelComm;
-
-            Pout<< "Solve generic on mesh (level=" << agglomLevel
-                << ") using communicator " << levelComm << endl;
-
-            // 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 coarseI = procAgglomMap[procI];
-
-                    Map<label>::iterator fnd = agglomToMaster.find(coarseI);
-                    if (fnd == agglomToMaster.end())
-                    {
-                        agglomToMaster.insert(coarseI, procI);
-                    }
-                    else
-                    {
-                        fnd() = max(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(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;
-
-            // Allocate a communicator for the processor-agglomerated matrix
-            label procAgglomComm = UPstream::allocateCommunicator
-            (
-                levelComm,
-                masterProcs
-            );
-            Pout<< "** Allocated communicator " << procAgglomComm
-                << " for indices " << masterProcs
-                << " in processor list " << UPstream::procID(levelComm)
-                << endl;
-
-
-
-            // Gather matrix and mesh onto agglomProcIDs[0]
-            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-            procAgglomerateMatrix
-            (
-                // Agglomeration information
-                procAgglomMap,
-                agglomProcIDs,
-                procAgglomComm,
-
-                agglomLevel,  // level (coarse, not fine level!)
-
-                // Resulting matrix
-                allMatrixPtr_,
-                allInterfaceBouCoeffs_,
-                allInterfaceIntCoeffs_,
-                allPrimitiveInterfaces_,
-                allInterfaces_
-            );
-
-
-            UPstream::warnComm = oldWarn;
+            //// Pick a level to processor agglomerate
+            //label agglomLevel = matrixLevels_.size() - 1;//1;
+            //
+            //
+            //// Get mesh and matrix at this level
+            //const lduMatrix& levelMatrix = matrixLevels_[agglomLevel];
+            //const lduMesh& levelMesh = levelMatrix.mesh();
+            //
+            //
+            //label levelComm = levelMesh.comm();
+            //label oldWarn = UPstream::warnComm;
+            //UPstream::warnComm = levelComm;
+            //
+            //Pout<< "Solve generic on mesh (level=" << agglomLevel
+            //    << ") using communicator " << levelComm << endl;
+            //
+            //// 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 coarseI = procAgglomMap[procI];
+            //
+            //        Map<label>::iterator fnd = agglomToMaster.find(coarseI);
+            //        if (fnd == agglomToMaster.end())
+            //        {
+            //            agglomToMaster.insert(coarseI, procI);
+            //        }
+            //        else
+            //        {
+            //            fnd() = max(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(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;
+            //
+            //// Allocate a communicator for the processor-agglomerated matrix
+            //label procAgglomComm = UPstream::allocateCommunicator
+            //(
+            //    levelComm,
+            //    masterProcs
+            //);
+            //Pout<< "** Allocated communicator " << procAgglomComm
+            //    << " for indices " << masterProcs
+            //    << " in processor list " << UPstream::procID(levelComm)
+            //    << endl;
+            //
+            //
+            //
+            //// Gather matrix and mesh onto agglomProcIDs[0]
+            //// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+            //
+            //procAgglomerateMatrix
+            //(
+            //    // Agglomeration information
+            //    procAgglomMap,
+            //    agglomProcIDs,
+            //    procAgglomComm,
+            //
+            //    agglomLevel,  // level (coarse, not fine level!)
+            //
+            //    // Resulting matrix
+            //    allMatrixPtr_,
+            //    allInterfaceBouCoeffs_,
+            //    allInterfaceIntCoeffs_,
+            //    allPrimitiveInterfaces_,
+            //    allInterfaces_
+            //);
+            //
+            //
+            //UPstream::warnComm = oldWarn;
         }
     }
     else
@@ -258,19 +452,19 @@ 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);
-            }
-        }
-    }
+//    // 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_)
     {
@@ -305,7 +499,6 @@ void Foam::GAMGSolver::readControls()
     controlDict_.readIfPresent("interpolateCorrection", interpolateCorrection_);
     controlDict_.readIfPresent("scaleCorrection", scaleCorrection_);
     controlDict_.readIfPresent("directSolveCoarsest", directSolveCoarsest_);
-    controlDict_.readIfPresent("processorAgglomerate", processorAgglomerate_);
 
     Pout<< "GAMGSolver settings :"
         << " cacheAgglomeration:" << cacheAgglomeration_
@@ -319,7 +512,6 @@ void Foam::GAMGSolver::readControls()
         << " interpolateCorrection:" << interpolateCorrection_
         << " scaleCorrection:" << scaleCorrection_
         << " directSolveCoarsest:" << directSolveCoarsest_
-        << " processorAgglomerate:" << processorAgglomerate_
         << endl;
 }
 
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.H b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.H
index 955a798318..3497895241 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.H
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.H
@@ -122,6 +122,8 @@ class GAMGSolver
         //  Warning: Needs to be deleted explicitly.
         PtrList<lduInterfaceFieldPtrsList> interfaceLevels_;
 
+        PtrList<PtrList<lduInterfaceField> > primitiveInterfaceLevels_;
+
         //- Hierarchy of interface boundary coefficients
         PtrList<FieldField<Field, scalar> > interfaceLevelsBouCoeffs_;
 
@@ -131,14 +133,14 @@ class GAMGSolver
         //- LU decompsed coarsest matrix
         autoPtr<LUscalarMatrix> coarsestLUMatrixPtr_;
 
-        // Processor matrix agglomeration
+        //// Processor matrix agglomeration
 
-            //- Combined coarsest level matrix for all processors
-            autoPtr<lduMatrix> allMatrixPtr_;
-            FieldField<Field, scalar> allInterfaceBouCoeffs_;
-            FieldField<Field, scalar> allInterfaceIntCoeffs_;
-            PtrList<lduInterfaceField> allPrimitiveInterfaces_;
-            lduInterfaceFieldPtrsList allInterfaces_;
+            ////- Combined coarsest level matrix for all processors
+            //PtrList<lduMatrix> procMatrixLevels_;
+            //PtrList<lduInterfaceFieldPtrsList> procInterfaceLevels_;
+            //PtrList<PtrList<lduInterfaceField> > procPrimitiveInterfaces_;
+            //PtrList<FieldField<Field, scalar> > procInterfaceLevelsBouCoeffs_;
+            //PtrList<FieldField<Field, scalar> > procInterfaceLevelsIntCoeffs_;
 
 
     // Private Member Functions
@@ -167,8 +169,25 @@ 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,
+            lduInterfaceFieldPtrsList& coarseInterfaces,
+            FieldField<Field, scalar>& coarseInterfaceBouCoeffs,
+            FieldField<Field, scalar>& coarseInterfaceIntCoeffs
+        ) const;
 
         //- Collect matrices from other processors
         void gatherMatrices
@@ -195,7 +214,6 @@ class GAMGSolver
             // Agglomeration information
             const labelList& procAgglomMap,
             const List<int>& agglomProcIDs,
-            label masterComm,
 
             const label levelI,
 
@@ -205,6 +223,14 @@ class GAMGSolver
             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
@@ -229,6 +255,7 @@ class GAMGSolver
             scalarField& field,
             scalarField& Acf,
             const lduMatrix& A,
+            const int comm,
             const FieldField<Field, scalar>& interfaceLevelBouCoeffs,
             const lduInterfaceFieldPtrsList& interfaceLevel,
             const scalarField& source,
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverAgglomerateMatrix.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverAgglomerateMatrix.C
index 6bd97630c9..e3754e47cc 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverAgglomerateMatrix.C
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverAgglomerateMatrix.C
@@ -30,27 +30,179 @@ License
 
 // * * * * * * * * * * * * * 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);
 
-    // Coarse matrix diagonal initialised by restricting the finer mesh diagonal
-    scalarField& coarseDiag = coarseMatrix.diag();
-    agglomeration_.restrictField(coarseDiag, fineMatrix.diag(), fineLevelIndex);
 
+        Pout<< "agglomerateMatrix : I have fine mesh:"
+            << fineMatrix.mesh().lduAddr().size()
+            << " at fine level:" << fineLevelIndex
+            << " constricting to coarse cells:"
+            << nCoarseCells
+            << " coarse nFaces:" << nCoarseFaces << endl;
+
+        // 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. 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
+        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,
+            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];
+                }
+            }
+        }
+    }
+}
+
+
+//XXXXX
+// Agglomerate only the interface coefficients.
+void Foam::GAMGSolver::agglomerateInterfaceCoefficients
+(
+    const label fineLevelIndex,
+    const lduInterfacePtrsList& coarseMeshInterfaces,
+    lduInterfaceFieldPtrsList& coarseInterfaces,
+    FieldField<Field, scalar>& coarseInterfaceBouCoeffs,
+    FieldField<Field, scalar>& coarseInterfaceIntCoeffs
+) const
+{
     // Get reference to fine-level interfaces
     const lduInterfaceFieldPtrsList& fineInterfaces =
         interfaceLevel(fineLevelIndex);
@@ -63,34 +215,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)
@@ -100,7 +230,7 @@ void Foam::GAMGSolver::agglomerateMatrix(const label fineLevelIndex)
             const GAMGInterface& coarseInterface =
                 refCast<const GAMGInterface>
                 (
-                    agglomeration_.interfaceLevel(fineLevelIndex + 1)[inti]
+                    coarseMeshInterfaces[inti]
                 );
 
             coarseInterfaces.set
@@ -113,97 +243,31 @@ void Foam::GAMGSolver::agglomerateMatrix(const label fineLevelIndex)
                 ).ptr()
             );
 
+            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 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
-        scalarField& coarseUpper = coarseMatrix.upper();
-        scalarField& coarseLower = coarseMatrix.lower();
-
-        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)
-        {
-            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 (cl[cFace] == restrictAddr[l[fineFacei]])
-                {
-                    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
-                    (
-                        "GAMGSolver::agglomerateMatrix(const label)"
-                    )   << "Inconsistent addressing between "
-                           "fine and coarse grids"
-                        << exit(FatalError);
-                }
-            }
-            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();
-
-        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];
-            }
         }
     }
 }
@@ -355,7 +419,6 @@ void Foam::GAMGSolver::procAgglomerateMatrix
     // Agglomeration information
     const labelList& procAgglomMap,
     const List<int>& agglomProcIDs,
-    label masterComm,
 
     const label levelI,
 
@@ -365,7 +428,7 @@ void Foam::GAMGSolver::procAgglomerateMatrix
     FieldField<Field, scalar>& allInterfaceIntCoeffs,
     PtrList<lduInterfaceField>& allPrimitiveInterfaces,
     lduInterfaceFieldPtrsList& allInterfaces
-)
+) const
 {
     const lduMatrix& coarsestMatrix = matrixLevels_[levelI];
     const lduInterfaceFieldPtrsList& coarsestInterfaces =
@@ -387,15 +450,6 @@ void Foam::GAMGSolver::procAgglomerateMatrix
     Pout<< "procAgglomerateMatrix :" << " level:" << levelI << endl;
 
 
-    agglomeration_.procAgglomerateLduAddressing
-    (
-        procAgglomMap,
-        agglomProcIDs,
-        masterComm,
-        levelI     // coarsest level
-    );
-
-
     // Gather all matrix coefficients onto agglomProcIDs[0]
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -439,12 +493,15 @@ void Foam::GAMGSolver::procAgglomerateMatrix
         Pout<< endl;
 
 
-        const lduMesh& allMesh = agglomeration_.procMeshLevel(levelI);
-        const labelList& cellOffsets = agglomeration_.cellOffsets(levelI);
-        const labelListList& faceMap = agglomeration_.faceMap(levelI);
-        const labelListList& boundaryMap = agglomeration_.boundaryMap(levelI);
+        //const lduMesh& allMesh = agglomeration_.procMeshLevel(levelI);
+        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);
+            agglomeration_.boundaryFaceMap(levelI+1);
+
+Pout<< "Agglomerating onto mesh:" << allMesh.info() << endl;
 
 
         allMatrixPtr.reset(new lduMatrix(allMesh));
@@ -740,4 +797,54 @@ Pout<< "    from proc:" << procI << " interface:" << procIntI
 //XXXX
 
 
+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);
+}
+
+
 // ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverScale.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverScale.C
index cd787047bb..d7f73c1266 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverScale.C
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverScale.C
@@ -33,6 +33,7 @@ void Foam::GAMGSolver::scale
     scalarField& field,
     scalarField& Acf,
     const lduMatrix& A,
+    const int comm,
     const FieldField<Field, scalar>& interfaceLevelBouCoeffs,
     const lduInterfaceFieldPtrsList& interfaceLevel,
     const scalarField& source,
@@ -58,11 +59,13 @@ void Foam::GAMGSolver::scale
     }
 
     vector2D scalingVector(scalingFactorNum, scalingFactorDenom);
-    matrix().mesh().reduce
-    (
-        scalingVector,
-        sumOp<vector2D>()
-    );
+    //A.mesh().reduce
+    //(
+    //    scalingVector,
+    //    sumOp<vector2D>()
+    //);
+    Foam::reduce(scalingVector, sumOp<vector2D>(), Pstream::msgType(), comm);
+
     scalar sf = scalingVector.x()/stabilise(scalingVector.y(), VSMALL);
 
     if (debug >= 2)
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C
index f62a4ee769..61022c353d 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C
@@ -142,8 +142,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_)
     {
@@ -154,66 +154,81 @@ 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_)
+        Pout<< "Restriction for level:" << leveli << endl;
+
+        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
+                (
+                    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)
+                {
+                    int comm = matrixLevels_[leveli].mesh().comm();
+
+                    scale
+                    (
+                        coarseCorrFields[leveli],
+                        const_cast<scalarField&>
+                        (
+                            ACf.operator const scalarField&()
+                        ),
+                        matrixLevels_[leveli],
+                        comm,
+                        interfaceLevelsBouCoeffs_[leveli],
+                        interfaceLevels_[leveli],
+                        coarseSources[leveli],
+                        cmpt
+                    );
+                }
+
+                // Correct the residual with the new solution
+                matrixLevels_[leveli].Amul
+                (
+                    const_cast<scalarField&>
+                    (
+                        ACf.operator const scalarField&()
+                    ),
+                    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_)
@@ -223,12 +238,15 @@ void Foam::GAMGSolver::Vcycle
 
 
     // Solve Coarsest level with either an iterative or direct solver
-    solveCoarsestLevel
-    (
-        coarseCorrFields[coarsestLevel],
-        coarseSources[coarsestLevel]
-    );
-
+    if (coarseCorrFields.set(coarsestLevel))
+    {
+        Pout<< "Coarsest solve for level:" << coarsestLevel << endl;
+        solveCoarsestLevel
+        (
+            coarseCorrFields[coarsestLevel],
+            coarseSources[coarsestLevel]
+        );
+    }
 
     if (debug >= 2)
     {
@@ -237,99 +255,141 @@ 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()
-        );
+        Pout<< "Smoothing and prolongation for level:" << leveli << endl;
 
-        // Only store the preSmoothedCoarseCorrField if pre-smoothing is used
-        if (nPreSweeps_)
+        if (coarseCorrFields.set(leveli))
         {
-            preSmoothedCoarseCorrField.assign(coarseCorrFields[leveli]);
-        }
+            Pout<< "Prolonging from " << leveli + 1 << " up to "
+                << leveli << endl;
+            //// 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()
+            //);
+            scalarField preSmoothedCoarseCorrField;
 
-        agglomeration_.prolongField
-        (
-            coarseCorrFields[leveli],
-            coarseCorrFields[leveli + 1],
-            leveli + 1
-        );
+            // Only store the preSmoothedCoarseCorrField if pre-smoothing is
+            // used
+            if (nPreSweeps_)
+            {
+                //preSmoothedCoarseCorrField.assign(coarseCorrFields[leveli]);
+                preSmoothedCoarseCorrField = 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
+            Pout<< "Doing stuff at level " << leveli << endl;
+
+            //// 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&());
+            scalarField ACfRef(coarseCorrFields[leveli].size());
+
+
+            if (interpolateCorrection_)
+            {
+                Pout<< "doing interpolate." << endl;
+                interpolate
+                (
+                    coarseCorrFields[leveli],
+                    ACfRef,
+                    matrixLevels_[leveli],
+                    interfaceLevelsBouCoeffs_[leveli],
+                    interfaceLevels_[leveli],
+                    coarseSources[leveli],
+                    cmpt
+                );
+                Pout<< "done interpolate." << endl;
+            }
+
+            // Scale coarse-grid correction field
+            // but not on the coarsest level because it evaluates to 1
+            if (scaleCorrection_ && leveli < coarsestLevel - 1)
+            {
+                //int comm =
+                //(
+                //    matrixLevels_.set(leveli+1)
+                //  ? matrixLevels_[leveli+1].mesh().comm()
+                //  : matrixLevels_[leveli].mesh().comm()
+                //);
+                int comm = matrixLevels_[leveli].mesh().comm();
+
+
+                Pout<< "doing scale with comm:" << comm << endl;
+                scale
+                (
+                    coarseCorrFields[leveli],
+                    ACfRef,
+                    matrixLevels_[leveli],
+                    comm,
+                    interfaceLevelsBouCoeffs_[leveli],
+                    interfaceLevels_[leveli],
+                    coarseSources[leveli],
+                    cmpt
+                );
+                Pout<< "done scale with comm:" << comm << endl;
+            }
+
+            // Only add the preSmoothedCoarseCorrField if pre-smoothing is
+            // used
+            if (nPreSweeps_)
+            {
+                coarseCorrFields[leveli] += preSmoothedCoarseCorrField;
+            }
+
+            Pout<< "doing smooth." << endl;
+            smoothers[leveli + 1].smooth
             (
                 coarseCorrFields[leveli],
-                ACfRef,
-                matrixLevels_[leveli],
-                interfaceLevelsBouCoeffs_[leveli],
-                interfaceLevels_[leveli],
                 coarseSources[leveli],
-                cmpt
+                cmpt,
+                min
+                (
+                    nPostSweeps_ + postSweepsLevelMultiplier_*leveli,
+                    maxPostSweeps_
+                )
             );
+            Pout<< "done smooth." << endl;
         }
-
-        // 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
+    Pout<< "Doing Prolong to finest level" << endl;
     agglomeration_.prolongField
     (
         finestCorrection,
         coarseCorrFields[0],
-        0
+        0,
+        true            //false               // no proc agglomeration for now
     );
+    Pout<< "Done Prolong to finest level" << endl;
 
     if (interpolateCorrection_)
     {
+        Pout<< "doing interpolate on finest level" << endl;
         interpolate
         (
             finestCorrection,
@@ -340,21 +400,26 @@ void Foam::GAMGSolver::Vcycle
             finestResidual,
             cmpt
         );
+        Pout<< "done interpolate on finest level" << endl;
     }
 
     if (scaleCorrection_)
     {
         // Scale the finest level correction
+        int comm = matrix_.mesh().comm();
+        Pout<< "doing scale on finest level with comm:" << comm << endl;
         scale
         (
             finestCorrection,
             Apsi,
             matrix_,
+            comm,
             interfaceBouCoeffs_,
             interfaces_,
             finestResidual,
             cmpt
         );
+        Pout<< "done scale on finest level with comm:" << comm << endl;
     }
 
     forAll(psi, i)
@@ -362,6 +427,7 @@ void Foam::GAMGSolver::Vcycle
         psi[i] += finestCorrection[i];
     }
 
+    Pout<< "Doing smooth on finest level" << endl;
     smoothers[0].smooth
     (
         psi,
@@ -369,6 +435,7 @@ void Foam::GAMGSolver::Vcycle
         cmpt,
         nFinestSweeps_
     );
+    Pout<< "Done smooth on finest level" << endl;
 }
 
 
@@ -400,37 +467,47 @@ 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()
-            )
-        );
+            Pout<< "initVCucle level:" << leveli << " nCoarseCells:"
+                << nCoarseCells << endl;
 
-        smoothers.set
-        (
-            leveli + 1,
-            lduMatrix::smoother::New
+            coarseSources.set(leveli, new scalarField(nCoarseCells));
+
+            //if (!matrixLevels_.set(leveli))
+            //{
+            //    coarseCorrFields.set(leveli, new scalarField(nCoarseCells));
+            //}
+        }
+
+        if (matrixLevels_.set(leveli))
+        {
+            const lduMatrix& mat = matrixLevels_[leveli];
+
+            label nCoarseCells = mat.diag().size();
+            Pout<< "initVCucle level:" << leveli << " matrix size:"
+                << nCoarseCells << endl;
+
+            coarseCorrFields.set(leveli, new scalarField(nCoarseCells));
+            //coarseCorrFields.set(leveli, new scalarField(nCoarseCells));
+            //coarseSources.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_
+                )
+            );
+        }
     }
 }
 
@@ -456,106 +533,111 @@ void Foam::GAMGSolver::solveCoarsestLevel
         coarsestCorrField = coarsestSource;
         coarsestLUMatrixPtr_->solve(coarsestCorrField);
     }
-    else if (processorAgglomerate_)
-    {
-        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;
-
-        if (Pstream::myProcNo(solveComm) != -1)
-        {
-            const lduMatrix& allMatrix = allMatrixPtr_();
-
-            allCorrField.setSize(allSource.size(), 0);
-
-            {
-                Pout<< "** Master:Solving on comm:" << solveComm
-                    << " with procs:" << UPstream::procID(solveComm) << endl;
-
-                if (allMatrix.asymmetric())
-                {
-                    coarseSolverPerf = BICCG
-                    (
-                        "coarsestLevelCorr",
-                        allMatrix,
-                        allInterfaceBouCoeffs_,
-                        allInterfaceIntCoeffs_,
-                        allInterfaces_,
-                        tolerance_,
-                        relTol_
-                    ).solve
-                    (
-                        allCorrField,
-                        allSource
-                    );
-                }
-                else
-                {
-                    coarseSolverPerf = ICCG
-                    (
-                        "coarsestLevelCorr",
-                        allMatrix,
-                        allInterfaceBouCoeffs_,
-                        allInterfaceIntCoeffs_,
-                        allInterfaces_,
-                        tolerance_,
-                        relTol_
-                    ).solve
-                    (
-                        allCorrField,
-                        allSource
-                    );
-                }
-            }
-        }
-
-        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 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
     {
         coarsestCorrField = 0;