Merge branch 'PstreamUpdate' into development

ThirdParty/AllMake.stage2

@@ -40,8 +40,8 @@
 # Requirements:
 #     1: Your foam-extend environment must be properly initialized
 #     2: AllMake.stage1 if you are overriding your system compiler
-#     3: The file etc/prefs.sh should be used for setting the variables enabling
-#        the compilation of the various packages
+#     3: The file etc/prefs.sh should be used for setting the variables
+#        enabling the compilation of the various packages
 #
 # Author:
 #     Martin Beaudoin, Hydro-Quebec, (2015)
@@ -104,6 +104,13 @@ then
 	( rpm_make -p openmpi-1.8.8 -s openmpi-1.8.8.spec  -u http://www.open-mpi.org/software/ompi/v1.8/downloads/openmpi-1.8.8.tar.gz \
           -f --define '_configureAdditionalArgs "$WM_THIRD_PARTY_USE_OPENMPI_188_ConfigureAdditionalArgs"')
     }
+    # mvipich2-2.2
+    #
+    [ ! -z "$WM_THIRD_PARTY_USE_MVAPICH2_22" ] && {
+        echo "Building mvapich2 2.2"
+        ( rpm_make -p mvapich2-2.2 -s mvapich2-2.2.spec -u http://mvapich.cse.ohio-state.edu/download/mvapich/mv2/mvapich2-2.2.tar.gz \
+          -f --define '_configureAdditionalArgs "$WM_THIRD_PARTY_USE_MVAPICH2_22_ConfigureAdditionalArgs"')
+    }
 else
     echo "Using system installed OpenMPI"
     echo ""

ThirdParty/rpmBuild/SPECS/mvapich2-2.2.spec

@@ -0,0 +1,268 @@
+#------------------------------------------------------------------------------
+# =========                 |
+# \\      /  F ield         | foam-extend: Open Source CFD
+#  \\    /   O peration     |
+#   \\  /    A nd           | For copyright notice see file Copyright
+#    \\/     M anipulation  |
+#------------------------------------------------------------------------------
+# License
+#     This file is part of foam-extend.
+#
+#     foam-extend 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.
+#
+#     foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+#
+# Script
+#     RPM spec file for mvapich2-2.2
+#
+# Description
+#     RPM spec file for creating a relocatable RPM
+#
+# Author:
+#     Martin Beaudoin, Hydro-Quebec, (2015)
+#
+#------------------------------------------------------------------------------
+
+# We grab the value of WM_THIRD_PARTY and WM_OPTIONS from the environment variable
+%{expand:%%define _WM_THIRD_PARTY_DIR %(echo $WM_THIRD_PARTY_DIR)}
+%{expand:%%define _WM_OPTIONS         %(echo $WM_OPTIONS)}
+
+# Disable the generation of debuginfo packages
+%define debug_package %{nil}
+
+# The topdir needs to point to the $WM_THIRD_PARTY/rpmbuild directory
+%define _topdir	 	%{_WM_THIRD_PARTY_DIR}/rpmBuild
+%define _tmppath	%{_topdir}/tmp
+
+# Will install the package directly $WM_THIRD_PARTY_DIR
+#   Some comments about package relocation:
+#   By using this prefix for the Prefix:  parameter in this file, you will make this
+#   package relocatable.
+#
+#   This is fine, as long as your software is itself relocatable.
+#
+#   Simply take note that libraries built with libtool are not relocatable because the
+#   prefix we specify will be hard-coded in the library .la files.
+#   Ref: http://sourceware.org/autobook/autobook/autobook_80.html
+#
+#   In that case, if you ever change the value of the $WM_THIRD_PARTY_DIR, you will
+#   not be able to reutilize this RPM, even though it is relocatable. You will need to
+#   regenerate the RPM.
+#
+%define _prefix         %{_WM_THIRD_PARTY_DIR}
+
+%define name		mvapich2
+%define release		%{_WM_OPTIONS}
+%define version 	2.2
+
+%define buildroot       %{_topdir}/BUILD/%{name}-%{version}-root
+
+BuildRoot:	        %{buildroot}
+Summary: 		mvapich2
+License: 		BSD
+Name: 			%{name}
+Version: 		%{version}
+Release: 		%{release}
+URL:            http://mvapich.cse.ohio-state.edu/download/mvapich/mv2/
+Source: 		%url/%{name}-%{version}.tar.gz
+Prefix: 		%{_prefix}
+Group: 			Development/Tools
+
+
+%define _installPrefix  %{_prefix}/packages/%{name}-%{version}/platforms/%{_WM_OPTIONS}
+
+#--------------------------------------------------------------------------
+#
+# Here, we define default compiling options for mvapich2
+#
+# One can override the option on the commande line : --define='MACRO EXPR'
+#
+%{!?_configureAdditionalArgs: %define _configureAdditionalArgs Undefined}
+
+
+%description
+%{summary}
+
+%prep
+%setup -q
+
+
+%build
+    # export WM settings in a form that GNU configure recognizes
+    [ -n "$WM_CC" ]         &&  export CC="$WM_CC"
+    [ -n "$WM_CXX" ]        &&  export CXX="$WM_CXX"
+    [ -n "$WM_CFLAGS" ]     &&  export CFLAGS="$WM_CFLAGS"
+    [ -n "$WM_CXXFLAGS" ]   &&  export CXXFLAGS="$WM_CXXFLAGS"
+    [ -n "$WM_LDFLAGS" ]    &&  export LDFLAGS="$WM_LDFLAGS"
+
+    set +x
+    echo ""
+    echo "Compilation options:"
+    echo "     _configureAdditionalArgs : %{_configureAdditionalArgs}"
+    echo ""
+    set -x
+
+    unset mpiWith
+    # Enable GridEngine if it appears to be in use
+    # If you don't want any integration with SGE, simply unset the SGE
+    # environment variable
+    if [ -n "$SGE_ROOT" ]
+    then
+        mpiWith="$mpiWith --with-sge"
+    else
+        mpiWith="$mpiWith --without-sge"
+	mpiWith="$mpiWith --enable-mca-no-build=ras-gridengine,pls-gridengine"
+    fi
+
+    # Infiniband support
+    # Adjust according to your local paths
+    # if [ -d /usr/local/ofed -a -d /usr/local/ofed/lib64 ]
+    # then
+    #     mpiWith="$mpiWith --with-openib=/usr/local/ofed"
+    #     mpiWith="$mpiWith --with-openib-libdir=/usr/local/ofed/lib64"
+    # fi
+
+    # NOTE. Option --disable-mcast disables Infiniband when libs are not present
+
+    ./configure \
+        --prefix=%{_installPrefix}  \
+        --exec_prefix=%{_installPrefix}  \
+        --enable-mpirun-prefix-by-default \
+        --enable-orterun-prefix-by-default \
+        --enable-shared \
+        --enable-mpi-cxx \
+        --disable-static \
+        --disable-mpi-f77 \
+        --disable-mpi-f90 \
+        --without-slurm \
+        --enable-mpi-profile $mpiWith \
+        --disable-vt \
+        echo %{?_configureAdditionalArgs}`
+
+    [ -z "$WM_NCOMPPROCS" ] && WM_NCOMPPROCS=1
+    make -j $WM_NCOMPPROCS
+
+%install
+    make install DESTDIR=$RPM_BUILD_ROOT
+
+    # Creation of foam-extend specific .csh and .sh files"
+
+    echo ""
+    echo "Generating foam-extend specific .csh and .sh files for the package %{name}-%{version}"
+    echo ""
+    #
+    # Generate package specific .sh file for foam-extend
+    #
+mkdir -p $RPM_BUILD_ROOT/%{_installPrefix}/etc
+cat << DOT_SH_EOF > $RPM_BUILD_ROOT/%{_installPrefix}/etc/%{name}-%{version}.sh
+# Load %{name}-%{version} libraries and binaries
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+export MVAPICH2_DIR=\$WM_THIRD_PARTY_DIR/packages/%{name}-%{version}/platforms/\$WM_OPTIONS
+export MVAPICH2_BIN_DIR=\$MVAPICH2_DIR/bin
+export MVAPICH2_LIB_DIR=\$MVAPICH2_DIR/lib
+
+# Enable access to the package runtime applications and libraries
+[ -d \$MVAPICH2_BIN_DIR ] && _foamAddPath \$MVAPICH2_BIN_DIR
+[ -d \$MVAPICH2_LIB_DIR ] && _foamAddLib  \$MVAPICH2_LIB_DIR
+
+export MPI_HOME=\$MVAPICH2_DIR
+export MPI_ARCH_PATH=\$MPI_HOME
+export OPAL_PREFIX=\$MPI_ARCH_PATH
+
+# We initialize the rest of the environment using mpicc --showme:
+# This does not work with MVAPICH2
+#export MVAPICH2_INCLUDE_DIR="\`mpicc --showme:incdirs\`"
+#export MVAPICH2_COMPILE_FLAGS="\`mpicc --showme:compile\`"
+#export MVAPICH2_LINK_FLAGS="\`mpicc --showme:link\`"
+
+# Set the foam-extend compilation flags
+export PINC=\$MVAPICH2_COMPILE_FLAGS
+export PLIBS=\$MVAPICH2_LINK_FLAGS
+
+
+if [ "\$FOAM_VERBOSE" -a "\$PS1" ]
+then
+    echo "  Environment variables defined for OpenMPI:"
+    echo "    MVAPICH2_BIN_DIR       : \$MVAPICH2_BIN_DIR"
+    echo "    MVAPICH2_LIB_DIR       : \$MVAPICH2_LIB_DIR"
+    echo "    MVAPICH2_INCLUDE_DIR   : \$MVAPICH2_INCLUDE_DIR"
+    echo "    MVAPICH2_COMPILE_FLAGS : \$MVAPICH2_COMPILE_FLAGS"
+    echo "    MVAPICH2_LINK_FLAGS    : \$MVAPICH2_LINK_FLAGS"
+    echo ""
+    echo "    MPI_HOME              : \$MPI_HOME"
+    echo "    MPI_ARCH_PATH         : \$MPI_ARCH_PATH"
+    echo "    OPAL_PREFIX           : \$OPAL_PREFIX"
+    echo "    PINC                  : \$PINC"
+    echo "    PLIBS                 : \$PLIBS"
+fi
+DOT_SH_EOF
+
+    #
+    # Generate package specific .csh file for foam-extend
+    #
+cat << DOT_CSH_EOF > $RPM_BUILD_ROOT/%{_installPrefix}/etc/%{name}-%{version}.csh
+# Load %{name}-%{version} libraries and binaries
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+setenv MVAPICH2_DIR \$WM_THIRD_PARTY_DIR/packages/%{name}-%{version}/platforms/\$WM_OPTIONS
+setenv MVAPICH2_BIN_DIR \$MVAPICH2_DIR/bin
+setenv MVAPICH2_LIB_DIR \$MVAPICH2_DIR/lib
+
+# Enable access to the package runtime applications and libraries
+if ( -e \$MVAPICH2_BIN_DIR ) then
+    _foamAddPath \$MVAPICH2_BIN_DIR
+endif
+
+if ( -e \$MVAPICH2_LIB_DIR ) then
+    _foamAddLib \$MVAPICH2_LIB_DIR
+endif
+
+setenv MPI_HOME \$MVAPICH2_DIR
+setenv MPI_ARCH_PATH \$MPI_HOME
+setenv OPAL_PREFIX \$MPI_ARCH_PATH
+
+# We initialize the rest of the environment using mpicc --showme:
+setenv MVAPICH2_INCLUDE_DIR "\`mpicc --showme:incdirs\`"
+setenv MVAPICH2_COMPILE_FLAGS "\`mpicc --showme:compile\`"
+setenv MVAPICH2_LINK_FLAGS "\`mpicc --showme:link\`"
+
+# Set the foam-extend compilation flags
+setenv PINC "\$MVAPICH2_COMPILE_FLAGS"
+setenv PLIBS "\$MVAPICH2_LINK_FLAGS"
+
+
+if (\$?FOAM_VERBOSE && \$?prompt) then
+    echo "  Environment variables defined for OpenMPI:"
+    echo "    MVAPICH2_BIN_DIR       : \$MVAPICH2_BIN_DIR"
+    echo "    MVAPICH2_LIB_DIR       : \$MVAPICH2_LIB_DIR"
+    echo "    MVAPICH2_INCLUDE_DIR   : \$MVAPICH2_INCLUDE_DIR"
+    echo "    MVAPICH2_COMPILE_FLAGS : \$MVAPICH2_COMPILE_FLAGS"
+    echo "    MVAPICH2_LINK_FLAGS    : \$MVAPICH2_LINK_FLAGS"
+    echo ""
+    echo "    MPI_HOME              : \$MPI_HOME"
+    echo "    MPI_ARCH_PATH         : \$MPI_ARCH_PATH"
+    echo "    OPAL_PREFIX           : \$OPAL_PREFIX"
+    echo "    PINC                  : \$PINC"
+    echo "    PLIBS                 : \$PLIBS"
+endif
+DOT_CSH_EOF
+
+    #finally, generate a .tgz file for systems where using rpm for installing packages
+    # as a non-root user might be a problem.
+    (mkdir -p  %{_topdir}/TGZS/%{_target_cpu}; cd $RPM_BUILD_ROOT/%{_prefix}; tar -zcvf %{_topdir}/TGZS/%{_target_cpu}/%{name}-%{version}.tgz  packages/%{name}-%{version})
+
+
+%clean
+rm -rf %{buildroot}
+
+%files
+%defattr(-,root,root)
+%{_installPrefix}

ThirdParty/tools/makeThirdPartyFunctionsForRPM

@@ -82,6 +82,7 @@ rpm_make()
     echo "Package name      : $_PACKAGE"
     echo "Package URL       : $_PACKAGE_URL"
     echo "RPM spec file name: $_SPECFILE"
+    echo "RPM file name     : $_RPMFILENAME"
     echo "Additional flags  : $_ADDITIONALFLAGS"
 
     # Shift options

applications/solvers/coupled/MRFPorousFoam/UEqn.H

@@ -5,6 +5,9 @@
         fvc::div(phi)
     );
 
+    // Update boundary velocity for consistency with the flux
+    mrfZones.correctBoundaryVelocity(U);
+
     // Momentum equation
     fvVectorMatrix UEqn
     (
@@ -67,11 +70,20 @@
 
         const scalarField& V = mesh.V().field();
 
-        // Note: this insertion should only happen in porous cell zones
-        // A rewrite of the porousZones class interface is needed.
+        // Note: insertion should only happen in porous cell zones
         // HJ, 14/Mar/2016
-        forAll (TUIn, cellI)
+
+        register label cellI;
+
+        forAll (pZones, pZoneI)
         {
+            const labelList& curZoneCells = pZones[pZoneI].zone();
+
+            // Loop over all cells in the zone
+            forAll (curZoneCells, zcI)
+            {
+                cellI = curZoneCells[zcI];
+
                 const scalar& cellV = V[cellI];
 
                 const tensor& cellTU = TUIn[cellI];
@@ -91,4 +103,5 @@
                 cellDD(2, 2) += cellV*cellTU.zz();
             }
         }
+    }
 }

applications/solvers/incompressible/MRFSimpleFoam/MRFSimpleFoam.C

@@ -61,53 +61,13 @@ int main(int argc, char *argv[])
 
         // Pressure-velocity SIMPLE corrector
         {
-            // Momentum predictor
-
-            tmp<fvVectorMatrix> UEqn
-            (
-                fvm::div(phi, U)
-              + turbulence->divDevReff()
-            );
-            mrfZones.addCoriolis(UEqn());
-
-            UEqn().relax();
-
-            solve(UEqn() == -fvc::grad(p));
-
-            volScalarField rAU = 1.0/UEqn().A();
-            U = rAU*UEqn().H();
-            UEqn.clear();
-
-            phi = fvc::interpolate(U, "interpolate(HbyA)") & mesh.Sf();
-            mrfZones.relativeFlux(phi);
-            adjustPhi(phi, U, p);
-
-            // Non-orthogonal pressure corrector loop
-            while (simple.correctNonOrthogonal())
-            {
-                fvScalarMatrix pEqn
-                (
-                    fvm::laplacian(rAU, p) == fvc::div(phi)
-                );
-
-                pEqn.setReference(pRefCell, pRefValue);
-                pEqn.solve();
-
-                if (simple.finalNonOrthogonalIter())
-                {
-                    phi -= pEqn.flux();
-                }
+#           include "UEqn.H"
+#           include "pEqn.H"
         }
 
-#           include "continuityErrs.H"
-
-            // Explicitly relax pressure for momentum corrector
-            p.relax();
-
-            // Momentum corrector
-            U -= rAU*fvc::grad(p);
-            U.correctBoundaryConditions();
-        }
+        // Calculate relative velocity
+        Urel == U;
+        mrfZones.relativeVelocity(Urel);
 
         turbulence->correct();
 

applications/solvers/incompressible/MRFSimpleFoam/UEqn.H

@@ -0,0 +1,19 @@
+    // Solve the momentum equation
+
+    tmp<fvVectorMatrix> HUEqn
+    (
+        fvm::div(phi, U)
+      + turbulence->divDevReff()
+    );
+    mrfZones.addCoriolis(HUEqn());
+
+    // Get under-relaxation factor
+    const scalar UUrf = mesh.solutionDict().equationRelaxationFactor(U.name());
+
+    // Momentum solution
+    solve
+    (
+        relax(HUEqn(), UUrf)
+     ==
+        -fvc::grad(p)
+    );

applications/solvers/incompressible/MRFSimpleFoam/createFields.H

@@ -44,3 +44,19 @@
 
     MRFZones mrfZones(mesh);
     mrfZones.correctBoundaryVelocity(U);
+
+    // Create Urel as a permanent field to make it available for on-the-fly
+    // post-processing operations
+    volVectorField Urel
+    (
+        IOobject
+        (
+            "Urel",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        U
+    );
+    mrfZones.relativeVelocity(Urel);

applications/solvers/incompressible/MRFSimpleFoam/pEqn.H

@@ -0,0 +1,50 @@
+{
+    p.boundaryField().updateCoeffs();
+
+    // Prepare clean 1/Ap without contribution from under-relaxation
+    // HJ, 26/Oct/2015
+    volScalarField rAU
+    (
+        "(1|A(U))",
+        1/HUEqn().A()
+    );
+
+    // Store velocity under-relaxation point before using U for
+    // the flux precursor
+    U.storePrevIter();
+
+    U = rAU*HUEqn().H();
+    HUEqn.clear();
+
+     phi = fvc::interpolate(U) & mesh.Sf();
+     mrfZones.relativeFlux(phi);
+     adjustPhi(phi, U, p);
+
+     // Non-orthogonal pressure corrector loop
+     while (simple.correctNonOrthogonal())
+     {
+         fvScalarMatrix pEqn
+         (
+             fvm::laplacian(rAU, p) == fvc::div(phi)
+         );
+
+         pEqn.setReference(pRefCell, pRefValue);
+         pEqn.solve();
+
+         if (simple.finalNonOrthogonalIter())
+         {
+             phi -= pEqn.flux();
+         }
+     }
+
+#   include "continuityErrs.H"
+
+    // Explicitly relax pressure for momentum corrector
+    p.relax();
+
+    // Momentum corrector
+    // Note: since under-relaxation does not change aU, H/a in U can be
+    // re-used.  HJ, 22/Jan/2016
+    U = UUrf*(U - rAU*fvc::grad(p)) + (1 - UUrf)*U.prevIter();
+    U.correctBoundaryConditions();
+}

applications/solvers/incompressible/simpleFoam/pEqn.H

@@ -1,3 +1,4 @@
+{
     p.boundaryField().updateCoeffs();
 
     // Prepare clean 1/Ap without contribution from under-relaxation
@@ -14,6 +15,7 @@
 
     U = rUA*HUEqn().H();
     HUEqn.clear();
+
     phi = fvc::interpolate(U) & mesh.Sf();
     adjustPhi(phi, U, p);
 
@@ -22,7 +24,7 @@
     {
         fvScalarMatrix pEqn
         (
-            fvm::laplacian(rUA, p) == fvc::div(phi)
+            fvm::laplacian(rAU, p) == fvc::div(phi)
         );
 
         pEqn.setReference(pRefCell, pRefValue);
@@ -43,6 +45,6 @@
     // Momentum corrector
     // Note: since under-relaxation does not change aU, H/a in U can be
     // re-used.  HJ, 22/Jan/2016
-    U = UUrf*(U - rUA*fvc::grad(p)) + (1 - UUrf)*U.prevIter();
+    U = UUrf*(U - rAU*fvc::grad(p)) + (1 - UUrf)*U.prevIter();
     U.correctBoundaryConditions();
-
+}

bin/plotForces.py

@@ -53,6 +53,19 @@ for i in range(1,len(t)):
         my[i] += mvy[i]
         mz[i] += mvz[i]
 
+# write clean data file
+outForces=open('forces.dat','w')
+
+for data in zip(t,fx,fy,fz):
+        outForces.write(' '.join([str(d) for d in data])+'\n')
+outForces.close()
+
+outMoments=open('moments.dat','w')
+
+for data in zip(t,mx,my,mz):
+        outMoments.write(' '.join([str(d) for d in data])+'\n')
+outMoments.close()
+
 # plot forces
 import pylab
 pylab.xlabel('iteration')

bin/plotResidual.py

@@ -62,6 +62,7 @@ for line in lines:
                 tepsilon.append(iepsilon)
                 epsilon.append(float(matchepsilon.group(2)))
 
+# write clean data file
 outfile=open('residual.dat','w')
 
 if iomega > 0:

etc/settings.csh

@@ -300,12 +300,17 @@ case SYSTEMOPENMPI:
     breaksw
 
 case MVAPICH2:
-    set mpi_version=mvapich2-2.2
-    if ($?WM_THIRD_PARTY_USE_MVAPICH2_22 != 0 && -d $WM_THIRD_PARTY_DIR/packages/mvapich2-2.2/platforms/$WM_OPTIONS ) then
-        if ($?FOAM_VERBOSE && $?prompt) then
-            echo "Using mvapich2-2.2 from the ThirdParty package: $WM_THIRD_PARTY_DIR/packages/$mpi_version"
+    set mpi_version=mvapich2
+
+    if ($?MVAPICH2_BIN_DIR != 0) then
+        if (-d "${MVAPICH2_BIN_DIR}" ) then
+        _foamAddPath $MVAPICH2_BIN_DIR
+        endif
+    else
+        set mpicc_cmd=`which mpicc`
+        setenv MVAPICH2_BIN_DIR `dirname $mpicc_cmd`
+        unset mpicc_cmd
     endif
-        _foamSource $WM_THIRD_PARTY_DIR/packages/$mpi_version/platforms/$WM_OPTIONS/etc/$mpi_version.csh
 
     setenv MPI_HOME `dirname $MVAPICH2_BIN_DIR`
     setenv MPI_ARCH_PATH $MPI_HOME
@@ -571,8 +576,8 @@ endif
 
 # zoltan
 # ~~~~~
-if ( $?ZOLTAN_SYSTEM == 0 && $?WM_THIRD_PARTY_USE_ZOLTAN_36 != 0 && -e "$WM_THIRD_PARTY_DIR"/packages/zoltan-3.6/platforms/$WM_OPTIONS ) then
-    _foamSource $WM_THIRD_PARTY_DIR/packages/zoltan-3.6/platforms/$WM_OPTIONS/etc/zoltan-3.6.csh
+if ( $?ZOLTAN_SYSTEM == 0 && $?WM_THIRD_PARTY_USE_ZOLTAN_35 != 0 && -e "$WM_THIRD_PARTY_DIR"/packages/zoltan-3.5/platforms/$WM_OPTIONS ) then
+    _foamSource $WM_THIRD_PARTY_DIR/packages/zoltan-3.5/platforms/$WM_OPTIONS/etc/zoltan-3.5.csh
 endif
 
 # Python

etc/settings.sh

@@ -375,22 +375,16 @@ SYSTEMOPENMPI)
     unset mpi_version
     ;;
 
+
 MVAPICH2)
-    mpi_version=mvapich2
-
-    if [ -n "${MVAPICH2_BIN_DIR}" ] && [ -d "${MVAPICH2_BIN_DIR}" ]
+    mpi_version=mvapich2-2.2
+    if [ ! -z $WM_THIRD_PARTY_USE_MVAPICH2_20 ] && [ -e $WM_THIRD_PARTY_DIR/packages/mvapich2-2.2/platforms/$WM_OPTIONS ]
         then
-        _foamAddPath $MVAPICH2_BIN_DIR
-    else
-        MVAPICH2_BIN_DIR=$(dirname `which mpicc`)
+        if [ "$FOAM_VERBOSE" -a "$PS1" ]
+        then
+            echo "Using mvapich2 from the ThirdParty package: $WM_THIRD_PARTY_DIR/packages/$mpi_version"
         fi
-
-    if which mpicc >/dev/null
-    then
-        mpicc -v 2>/dev/null | grep -q "mpicc for MVAPICH2" ||
-            echo "Warning: `which mpicc` does not appear to be for MVAPICH2"
-    else
-        echo "Warning: mpicc not available"
+        _foamSource $WM_THIRD_PARTY_DIR/packages/$mpi_version/platforms/$WM_OPTIONS/etc/$mpi_version.sh
     fi
 
     export MPI_HOME=`dirname $MVAPICH2_BIN_DIR`

src/finiteArea/Make/files

@@ -48,6 +48,8 @@ derivedFaPatchFields = $(faPatchFields)/derived
 $(derivedFaPatchFields)/fixedValueOutflow/fixedValueOutflowFaPatchFields.C
 $(derivedFaPatchFields)/inletOutlet/inletOutletFaPatchFields.C
 $(derivedFaPatchFields)/slip/slipFaPatchFields.C
+$(derivedFaPatchFields)/edgeNormalFixedValue/edgeNormalFixedValueFaPatchVectorField.C
+$(derivedFaPatchFields)/timeVaryingUniformFixedValue/timeVaryingUniformFixedValueFaPatchFields.C
 
 faePatchFields = fields/faePatchFields
 $(faePatchFields)/faePatchField/faePatchFields.C
@@ -92,6 +94,7 @@ $(ddtSchemes)/backwardFaDdtScheme/backwardFaDdtSchemes.C
 $(ddtSchemes)/boundedBackwardFaDdtScheme/boundedBackwardFaDdtScheme.C
 
 divSchemes = finiteArea/divSchemes
+finiteArea/fam/vectorFamDiv.C
 $(divSchemes)/faDivScheme/faDivSchemes.C
 $(divSchemes)/gaussFaDivScheme/gaussFaDivSchemes.C
 

src/finiteArea/faMesh/faMesh.C

@@ -1047,6 +1047,18 @@ void Foam::faMesh::addFaPatches(const List<faPatch*>& p)
 }
 
 
+Foam::label Foam::faMesh::comm() const
+{
+    return comm_;
+}
+
+
+Foam::label& Foam::faMesh::comm()
+{
+    return comm_;
+}
+
+
 const Foam::objectRegistry& Foam::faMesh::thisDb() const
 {
     return mesh().thisDb();

src/finiteArea/faMesh/faMesh.H

@@ -110,6 +110,12 @@ class faMesh
             mutable label nFaces_;
 
 
+        // Communication support
+
+            //- Communicator used for parallel communication
+            label comm_;
+
+
     // Demand-driven data
 
         //- Primitive patch
@@ -287,8 +293,7 @@ public:
         );
 
 
-    // Destructor
-
+    //- Destructor
     virtual ~faMesh();
 
 
@@ -369,6 +374,15 @@ public:
                 }
 
 
+            // Communication support
+
+                //- Return communicator used for parallel communication
+                label comm() const;
+
+                //- Return communicator used for parallel communication
+                label& comm();
+
+
         // Access
 
             //- Return reference to the mesh database

src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.C

@@ -55,6 +55,18 @@ processorFaPatch::~processorFaPatch()
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
+Foam::label Foam::processorFaPatch::comm() const
+{
+    return boundaryMesh().mesh().comm();
+}
+
+
+int Foam::processorFaPatch::tag() const
+{
+    return Pstream::msgType();
+}
+
+
 void processorFaPatch::makeNonGlobalPatchPoints() const
 {
     // If it is not runing parallel or there are no global points

src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.H

@@ -37,7 +37,6 @@ SourceFiles
 
 #include "coupledFaPatch.H"
 #include "processorLduInterface.H"
-// #include "processorPolyPatch.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -55,7 +54,10 @@ class processorFaPatch
 {
     // Private data
 
+        //- My processro number
         int myProcNo_;
+
+        //- Neighbour processor number
         int neighbProcNo_;
 
         //- Processor-neighbbour patch edge centres
@@ -75,6 +77,7 @@ class processorFaPatch
         //  non-global, i.e. present in this processor patch
         mutable labelList* nonGlobalPatchPointsPtr_;
 
+
 protected:
 
     // Protected Member functions
@@ -88,9 +91,8 @@ protected:
         //- Find non-globa patch points
         void makeNonGlobalPatchPoints() const;
 
-protected:
 
-    // Protected Member functions
+        // Geometry functions
 
             //- Initialise the calculation of the patch geometry
             void initGeometry();
@@ -110,6 +112,7 @@ protected:
             //- Update of the patch topology
             virtual void updateMesh();
 
+
 public:
 
     //- Runtime type information
@@ -160,8 +163,8 @@ public:
             nonGlobalPatchPointsPtr_(NULL)
         {}
 
-    // Destructor
 
+    //- Destructor
     virtual ~processorFaPatch();
 
 
@@ -192,6 +195,16 @@ public:
             }
         }
 
+
+        // Communications support
+
+            //- Return communicator used for communication
+            virtual label comm() const;
+
+            //- Return message tag to use for communication
+            virtual int tag() const;
+
+
         //- Return face transformation tensor
         virtual const tensorField& forwardT() const
         {

src/finiteArea/fields/faPatchFields/basic/transform/transformFaPatchFields.C

@@ -28,7 +28,6 @@ Description
 #include "faPatchFields.H"
 #include "transformFaPatchFields.H"
 #include "addToRunTimeSelectionTable.H"
-// #include "volFields.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -37,10 +36,6 @@ namespace Foam
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-// defineNamedTemplateTypeNameAndDebug(transformFaPatchScalarField, 0);
-// defineNamedTemplateTypeNameAndDebug(transformFaPatchVectorField, 0);
-// defineNamedTemplateTypeNameAndDebug(transformFaPatchTensorField, 0);
-
 makeFaPatchFieldsTypeName(transform);
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/finiteArea/fields/faPatchFields/constraint/processor/processorFaPatchField.C

@@ -176,7 +176,7 @@ void processorFaPatchField<Type>::initEvaluate
 {
     if (Pstream::parRun())
     {
-        procPatch_.compressedSend(commsType, this->patchInternalField()());
+        procPatch_.send(commsType, this->patchInternalField()());
     }
 }
 
@@ -189,7 +189,7 @@ void processorFaPatchField<Type>::evaluate
 {
     if (Pstream::parRun())
     {
-        procPatch_.compressedReceive<Type>(commsType, *this);
+        procPatch_.receive<Type>(commsType, *this);
 
         if (doTransform())
         {
@@ -218,7 +218,7 @@ void processorFaPatchField<Type>::initInterfaceMatrixUpdate
     const bool switchToLhs
 ) const
 {
-    procPatch_.compressedSend
+    procPatch_.send
     (
         commsType,
         this->patch().patchInternalField(psiInternal)()
@@ -240,7 +240,7 @@ void processorFaPatchField<Type>::updateInterfaceMatrix
 {
     scalarField pnf
     (
-        procPatch_.compressedReceive<scalar>(commsType, this->size())()
+        procPatch_.receive<scalar>(commsType, this->size())()
     );
 
     // Transform according to the transformation tensor

src/finiteArea/fields/faPatchFields/constraint/processor/processorFaPatchScalarField.C

@@ -52,7 +52,7 @@ void processorFaPatchField<scalar>::initInterfaceMatrixUpdate
     const bool switchToLhs
 ) const
 {
-    procPatch_.compressedSend
+    procPatch_.send
     (
         commsType,
         patch().patchInternalField(psiInternal)()
@@ -74,7 +74,7 @@ void processorFaPatchField<scalar>::updateInterfaceMatrix
 {
     scalarField pnf
     (
-        procPatch_.compressedReceive<scalar>(commsType, this->size())()
+        procPatch_.receive<scalar>(commsType, this->size())()
     );
 
     const unallocLabelList& edgeFaces = patch().edgeFaces();

src/finiteArea/fields/faPatchFields/constraint/wedge/wedgeFaPatchFields.C

@@ -28,7 +28,6 @@ Description
 #include "faPatchFields.H"
 #include "wedgeFaPatchFields.H"
 #include "addToRunTimeSelectionTable.H"
-// #include "volFields.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -37,10 +36,6 @@ namespace Foam
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-// defineNamedTemplateTypeNameAndDebug(wedgeFaPatchScalarField, 0);
-// defineNamedTemplateTypeNameAndDebug(wedgeFaPatchVectorField, 0);
-// defineNamedTemplateTypeNameAndDebug(wedgeFaPatchTensorField, 0);
-
 makeFaPatchFields(wedge);
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/finiteArea/fields/faPatchFields/derived/edgeNormalFixedValue/edgeNormalFixedValueFaPatchVectorField.C

@@ -0,0 +1,155 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "edgeNormalFixedValueFaPatchVectorField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "areaFields.H"
+#include "faPatchFieldMapper.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::edgeNormalFixedValueFaPatchVectorField::
+edgeNormalFixedValueFaPatchVectorField
+(
+    const faPatch& p,
+    const DimensionedField<vector, areaMesh>& iF
+)
+:
+    fixedValueFaPatchVectorField(p, iF),
+    refValue_(p.size(), 0)
+{}
+
+
+Foam::edgeNormalFixedValueFaPatchVectorField::
+edgeNormalFixedValueFaPatchVectorField
+(
+    const edgeNormalFixedValueFaPatchVectorField& ptf,
+    const faPatch& p,
+    const DimensionedField<vector, areaMesh>& iF,
+    const faPatchFieldMapper& mapper
+)
+:
+    fixedValueFaPatchVectorField(ptf, p, iF, mapper),
+    refValue_(ptf.refValue_, mapper)
+{}
+
+
+Foam::edgeNormalFixedValueFaPatchVectorField::
+edgeNormalFixedValueFaPatchVectorField
+(
+    const faPatch& p,
+    const DimensionedField<vector, areaMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFaPatchVectorField(p, iF, dict),
+    refValue_("refValue", dict, p.size())
+{}
+
+
+Foam::edgeNormalFixedValueFaPatchVectorField::
+edgeNormalFixedValueFaPatchVectorField
+(
+    const edgeNormalFixedValueFaPatchVectorField& pivpvf
+)
+:
+    fixedValueFaPatchVectorField(pivpvf),
+    refValue_(pivpvf.refValue_)
+{}
+
+
+Foam::edgeNormalFixedValueFaPatchVectorField::
+edgeNormalFixedValueFaPatchVectorField
+(
+    const edgeNormalFixedValueFaPatchVectorField& pivpvf,
+    const DimensionedField<vector, areaMesh>& iF
+)
+:
+    fixedValueFaPatchVectorField(pivpvf, iF),
+    refValue_(pivpvf.refValue_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::edgeNormalFixedValueFaPatchVectorField::autoMap
+(
+    const faPatchFieldMapper& m
+)
+{
+    fixedValueFaPatchVectorField::autoMap(m);
+    refValue_.autoMap(m);
+}
+
+
+void Foam::edgeNormalFixedValueFaPatchVectorField::rmap
+(
+    const faPatchVectorField& ptf,
+    const labelList& addr
+)
+{
+    fixedValueFaPatchVectorField::rmap(ptf, addr);
+
+    const edgeNormalFixedValueFaPatchVectorField& tiptf =
+        refCast<const edgeNormalFixedValueFaPatchVectorField>(ptf);
+
+    refValue_.rmap(tiptf.refValue_, addr);
+}
+
+
+void Foam::edgeNormalFixedValueFaPatchVectorField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    // Bug fix: update for moving mesh.  HJ, 15/Oct/2010
+    operator==(refValue_*patch().edgeNormals());
+}
+
+
+void Foam::edgeNormalFixedValueFaPatchVectorField::write(Ostream& os) const
+{
+    fixedValueFaPatchVectorField::write(os);
+    refValue_.writeEntry("refValue", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+makeFaPatchTypeField
+(
+    faPatchVectorField,
+    edgeNormalFixedValueFaPatchVectorField
+);
+
+} // End namespace Foam
+
+
+// ************************************************************************* //

src/finiteArea/fields/faPatchFields/derived/edgeNormalFixedValue/edgeNormalFixedValueFaPatchVectorField.H

@@ -0,0 +1,173 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::edgeNormalFixedValueFaPatchVectorField
+
+Description
+    Edge normal fixed value vector field finite area boundary condition
+
+    Describes a surface normal vector boundary condition by its magnitude.
+    Note: The value is positive for outward-pointing vectors
+
+SourceFiles
+    edgeNormalFixedValueFaPatchVectorField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef edgeNormalFixedValueFaPatchVectorField_H
+#define edgeNormalFixedValueFaPatchVectorField_H
+
+#include "faPatchFields.H"
+#include "fixedValueFaPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                 Class edgeNormalFixedValueFaPatch Declaration
+\*---------------------------------------------------------------------------*/
+
+class edgeNormalFixedValueFaPatchVectorField
+:
+    public fixedValueFaPatchVectorField
+{
+    // Private data
+
+        //- Surface-normal velocity value
+        scalarField refValue_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("edgeNormalFixedValue");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        edgeNormalFixedValueFaPatchVectorField
+        (
+            const faPatch&,
+            const DimensionedField<vector, areaMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        edgeNormalFixedValueFaPatchVectorField
+        (
+            const faPatch&,
+            const DimensionedField<vector, areaMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        //  edgeNormalFixedValueFaPatchVectorField
+        //  onto a new patch
+        edgeNormalFixedValueFaPatchVectorField
+        (
+            const edgeNormalFixedValueFaPatchVectorField&,
+            const faPatch&,
+            const DimensionedField<vector, areaMesh>&,
+            const faPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        edgeNormalFixedValueFaPatchVectorField
+        (
+            const edgeNormalFixedValueFaPatchVectorField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<faPatchVectorField> clone() const
+        {
+            return tmp<faPatchVectorField>
+            (
+                new edgeNormalFixedValueFaPatchVectorField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        edgeNormalFixedValueFaPatchVectorField
+        (
+            const edgeNormalFixedValueFaPatchVectorField&,
+            const DimensionedField<vector, areaMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<faPatchVectorField> clone
+        (
+            const DimensionedField<vector, areaMesh>& iF
+        ) const
+        {
+            return tmp<faPatchVectorField>
+            (
+                new edgeNormalFixedValueFaPatchVectorField
+                (
+                    *this,
+                    iF
+                )
+            );
+        }
+
+
+
+    // Member functions
+
+        // Mapping functions
+
+            //- Map (and resize as needed) from self given a mapping object
+            virtual void autoMap
+            (
+                const faPatchFieldMapper&
+            );
+
+            //- Reverse map the given faPatchField onto this faPatchField
+            virtual void rmap
+            (
+                const faPatchVectorField&,
+                const labelList&
+            );
+
+
+    // Member functions
+
+        //- Update the coefficients associated with the patch field
+        virtual void updateCoeffs();
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteArea/fields/faPatchFields/derived/timeVaryingUniformFixedValue/timeVaryingUniformFixedValueFaPatchField.C

@@ -0,0 +1,137 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "timeVaryingUniformFixedValueFaPatchField.H"
+#include "foamTime.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::timeVaryingUniformFixedValueFaPatchField<Type>::
+timeVaryingUniformFixedValueFaPatchField
+(
+    const faPatch& p,
+    const DimensionedField<Type, areaMesh>& iF
+)
+:
+    fixedValueFaPatchField<Type>(p, iF),
+    timeSeries_()
+{}
+
+
+template<class Type>
+Foam::timeVaryingUniformFixedValueFaPatchField<Type>::
+timeVaryingUniformFixedValueFaPatchField
+(
+    const faPatch& p,
+    const DimensionedField<Type, areaMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFaPatchField<Type>(p, iF),
+    timeSeries_(dict)
+{
+   if (dict.found("value"))
+   {
+       faPatchField<Type>::operator==(Field<Type>("value", dict, p.size()));
+   }
+   else
+   {
+       updateCoeffs();
+   }
+}
+
+
+template<class Type>
+Foam::timeVaryingUniformFixedValueFaPatchField<Type>::
+timeVaryingUniformFixedValueFaPatchField
+(
+    const timeVaryingUniformFixedValueFaPatchField<Type>& ptf,
+    const faPatch& p,
+    const DimensionedField<Type, areaMesh>& iF,
+    const faPatchFieldMapper& mapper
+)
+:
+    fixedValueFaPatchField<Type>(ptf, p, iF, mapper),
+    timeSeries_(ptf.timeSeries_)
+{}
+
+
+template<class Type>
+Foam::timeVaryingUniformFixedValueFaPatchField<Type>::
+timeVaryingUniformFixedValueFaPatchField
+(
+    const timeVaryingUniformFixedValueFaPatchField<Type>& ptf
+)
+:
+    fixedValueFaPatchField<Type>(ptf),
+    timeSeries_(ptf.timeSeries_)
+{}
+
+
+template<class Type>
+Foam::timeVaryingUniformFixedValueFaPatchField<Type>::
+timeVaryingUniformFixedValueFaPatchField
+(
+    const timeVaryingUniformFixedValueFaPatchField<Type>& ptf,
+    const DimensionedField<Type, areaMesh>& iF
+)
+:
+    fixedValueFaPatchField<Type>(ptf, iF),
+    timeSeries_(ptf.timeSeries_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::timeVaryingUniformFixedValueFaPatchField<Type>::updateCoeffs()
+{
+    if (this->updated())
+    {
+        return;
+    }
+
+    faPatchField<Type>::operator==
+    (
+        timeSeries_(this->db().time().timeOutputValue())
+    );
+    fixedValueFaPatchField<Type>::updateCoeffs();
+}
+
+
+template<class Type>
+void Foam::timeVaryingUniformFixedValueFaPatchField<Type>::write
+(
+    Ostream& os
+) const
+{
+    faPatchField<Type>::write(os);
+    timeSeries_.write(os);
+    this->writeEntry("value", os);
+}
+
+
+// ************************************************************************* //

src/finiteArea/fields/faPatchFields/derived/timeVaryingUniformFixedValue/timeVaryingUniformFixedValueFaPatchField.H

@@ -0,0 +1,182 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::timeVaryingUniformFixedValueFaPatchField
+
+Description
+    A time-varying form of a uniform fixed value finite area
+    boundary condition.
+
+    Example of the boundary condition specification:
+    @verbatim
+    inlet
+    {
+        type            timeVaryingUniformFixedValue;
+        fileName        "$FOAM_CASE/time-series";
+        outOfBounds     clamp;           // (error|warn|clamp|repeat)
+    }
+    @endverbatim
+
+Note
+    This class is derived directly from a fixedValue patch rather than from
+    a uniformFixedValue patch.
+
+See Also
+    Foam::interpolationTable and Foam::fixedValueFaPatchField
+
+SourceFiles
+    timeVaryingUniformFixedValueFaPatchField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef timeVaryingUniformFixedValueFaPatchField_H
+#define timeVaryingUniformFixedValueFaPatchField_H
+
+#include "fixedValueFaPatchField.H"
+#include "interpolationTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+           Class timeVaryingUniformFixedValueFaPatch Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class timeVaryingUniformFixedValueFaPatchField
+:
+    public fixedValueFaPatchField<Type>
+{
+    // Private data
+
+        //- The time series being used, including the bounding treatment
+        interpolationTable<Type> timeSeries_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("timeVaryingUniformFixedValue");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        timeVaryingUniformFixedValueFaPatchField
+        (
+            const faPatch&,
+            const DimensionedField<Type, areaMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        timeVaryingUniformFixedValueFaPatchField
+        (
+            const faPatch&,
+            const DimensionedField<Type, areaMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given patch field onto a new patch
+        timeVaryingUniformFixedValueFaPatchField
+        (
+            const timeVaryingUniformFixedValueFaPatchField<Type>&,
+            const faPatch&,
+            const DimensionedField<Type, areaMesh>&,
+            const faPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        timeVaryingUniformFixedValueFaPatchField
+        (
+            const timeVaryingUniformFixedValueFaPatchField<Type>&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<faPatchField<Type> > clone() const
+        {
+            return tmp<faPatchField<Type> >
+            (
+                new timeVaryingUniformFixedValueFaPatchField<Type>(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        timeVaryingUniformFixedValueFaPatchField
+        (
+            const timeVaryingUniformFixedValueFaPatchField<Type>&,
+            const DimensionedField<Type, areaMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<faPatchField<Type> > clone
+        (
+            const DimensionedField<Type, areaMesh>& iF
+        ) const
+        {
+            return tmp<faPatchField<Type> >
+            (
+                new timeVaryingUniformFixedValueFaPatchField<Type>(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        // Access
+
+            //- Return the time series used
+            const interpolationTable<Type>& timeSeries() const
+            {
+                return timeSeries_;
+            }
+
+
+        // Evaluation functions
+
+            //- Update the coefficients associated with the patch field
+            virtual void updateCoeffs();
+
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "timeVaryingUniformFixedValueFaPatchField.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteArea/fields/faPatchFields/derived/timeVaryingUniformFixedValue/timeVaryingUniformFixedValueFaPatchFields.C

@@ -21,28 +21,23 @@ License
     You should have received a copy of the GNU General Public License
     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
 
-Description
-    Prints out a description of the streams
-
 \*---------------------------------------------------------------------------*/
 
-#include "IPstream.H"
-#include "OPstream.H"
+#include "timeVaryingUniformFixedValueFaPatchFields.H"
+#include "addToRunTimeSelectionTable.H"
+#include "areaFields.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-void Foam::IPstream::print(Ostream& os) const
+namespace Foam
 {
-    os  << "Reading from processor " << fromProcNo_
-        << " to processor " << myProcNo() << Foam::endl;
-}
 
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-void Foam::OPstream::print(Ostream& os) const
-{
-    os  << "Writing from processor " << toProcNo_
-        << " to processor " << myProcNo() << Foam::endl;
-}
+makeFaPatchFields(timeVaryingUniformFixedValue);
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
 
 // ************************************************************************* //

src/finiteArea/fields/faPatchFields/derived/timeVaryingUniformFixedValue/timeVaryingUniformFixedValueFaPatchFields.H

@@ -0,0 +1,49 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef timeVaryingUniformFixedValueFaPatchFields_H
+#define timeVaryingUniformFixedValueFaPatchFields_H
+
+#include "timeVaryingUniformFixedValueFaPatchField.H"
+#include "fieldTypes.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makeFaPatchTypeFieldTypedefs(timeVaryingUniformFixedValue)
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteArea/fields/faPatchFields/derived/timeVaryingUniformFixedValue/timeVaryingUniformFixedValueFaPatchFieldsFwd.H

@@ -0,0 +1,51 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef timeVaryingUniformFixedValueFaPatchFieldsFwd_H
+#define timeVaryingUniformFixedValueFaPatchFieldsFwd_H
+
+#include "faPatchField.H"
+#include "fieldTypes.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type> class timeVaryingUniformFixedValueFaPatchField;
+
+makeFaPatchTypeFieldTypedefs(timeVaryingUniformFixedValue);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteArea/fields/faPatchFields/faPatchField/faPatchField.H

@@ -468,21 +468,7 @@ public:
 #endif
 
 
-#define makeFaPatchTypeFieldTypeName(typePatchTypeField)                   \
-                                                                           \
-defineNamedTemplateTypeNameAndDebug(typePatchTypeField, 0);
-
-#define makeFaPatchFieldsTypeName(typePatchField)                          \
-                                                                           \
-makeFaPatchTypeFieldTypeName(typePatchField##FaPatchScalarField);          \
-makeFaPatchTypeFieldTypeName(typePatchField##FaPatchVectorField);          \
-makeFaPatchTypeFieldTypeName(typePatchField##FaPatchSphericalTensorField); \
-makeFaPatchTypeFieldTypeName(typePatchField##FaPatchSymmTensorField);      \
-makeFaPatchTypeFieldTypeName(typePatchField##FaPatchTensorField);
-
-#define makeFaPatchTypeField(PatchTypeField, typePatchTypeField)           \
-                                                                           \
-defineNamedTemplateTypeNameAndDebug(typePatchTypeField, 0);                \
+#define addToFaPatchFieldRunTimeSelection(PatchTypeField, typePatchTypeField) \
                                                                               \
 addToRunTimeSelectionTable                                                    \
 (                                                                             \
@@ -502,17 +488,58 @@ addToRunTimeSelectionTable                                                 \
 );
 
 
+#define makeFaPatchTypeFieldTypeName(typePatchTypeField)                      \
+                                                                              \
+defineNamedTemplateTypeNameAndDebug(typePatchTypeField, 0);
+
+
+#define makeFaPatchFieldsTypeName(typePatchField)                             \
+                                                                              \
+makeFaPatchTypeFieldTypeName(typePatchField##FaPatchScalarField);             \
+makeFaPatchTypeFieldTypeName(typePatchField##FaPatchVectorField);             \
+makeFaPatchTypeFieldTypeName(typePatchField##FaPatchSphericalTensorField);    \
+makeFaPatchTypeFieldTypeName(typePatchField##FaPatchSymmTensorField);         \
+makeFaPatchTypeFieldTypeName(typePatchField##FaPatchTensorField);
+
+
+#define makeFaPatchTypeField(PatchTypeField, typePatchTypeField)              \
+                                                                              \
+defineTypeNameAndDebug(typePatchTypeField, 0);                                \
+                                                                              \
+addToFaPatchFieldRunTimeSelection                                             \
+(                                                                             \
+    PatchTypeField, typePatchTypeField                                        \
+);
+
+#define makeTemplateFaPatchTypeField(PatchTypeField, typePatchTypeField)      \
+                                                                              \
+defineNamedTemplateTypeNameAndDebug(typePatchTypeField, 0);                   \
+                                                                              \
+addToFaPatchFieldRunTimeSelection                                             \
+(                                                                             \
+    PatchTypeField, typePatchTypeField                                        \
+);
+
+
 #define makeFaPatchFields(type)                                               \
                                                                               \
-makeFaPatchTypeField(faPatchScalarField, type##FaPatchScalarField);        \
-makeFaPatchTypeField(faPatchVectorField, type##FaPatchVectorField);        \
-makeFaPatchTypeField                                                       \
+makeTemplateFaPatchTypeField(faPatchScalarField, type##FaPatchScalarField);   \
+makeTemplateFaPatchTypeField(faPatchVectorField, type##FaPatchVectorField);   \
+makeTemplateFaPatchTypeField                                                  \
 (                                                                             \
     faPatchSphericalTensorField,                                              \
     type##FaPatchSphericalTensorField                                         \
 );                                                                            \
-makeFaPatchTypeField(faPatchSymmTensorField, type##FaPatchSymmTensorField);\
-makeFaPatchTypeField(faPatchTensorField, type##FaPatchTensorField);
+makeTemplateFaPatchTypeField                                                  \
+(                                                                             \
+    faPatchSymmTensorField,                                                   \
+    type##FaPatchSymmTensorField                                              \
+);                                                                            \
+makeTemplateFaPatchTypeField                                                  \
+(                                                                             \
+    faPatchTensorField,                                                       \
+    type##FaPatchTensorField                                                  \
+);
 
 
 #define makeFaPatchTypeFieldTypedefs(type)                                    \

src/finiteArea/finiteArea/divSchemes/gaussFaDivScheme/gaussFaDivScheme.C

@@ -62,7 +62,7 @@ gaussDivScheme<Type>::facDiv
         (
             this->mesh_.Le() & this->tinterpScheme_().interpolate(vf)
         )
-      - this->mesh_.faceCurvatures()*(this->mesh_.faceAreaNormals()&vf)
+        // Removed for consistencty.  Matthias Rauter, 6/Dec/2016
     );
 
     tDiv().rename("div(" + vf.name() + ')');

src/finiteArea/finiteArea/fac/fac.H

@@ -43,6 +43,8 @@ Description
 #include "facAverage.H"
 #include "facLnGrad.H"
 #include "facDdt.H"
+#include "facNGrad.H"
+#include "facNDiv.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/finiteArea/finiteArea/fac/facDiv.C

@@ -28,6 +28,7 @@ License
 #include "facEdgeIntegrate.H"
 #include "faDivScheme.H"
 #include "faConvectionScheme.H"
+#include "transformField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -48,7 +49,17 @@ div
     const GeometricField<Type, faePatchField, edgeMesh>& ssf
 )
 {
-    return fac::edgeIntegrate(ssf);
+    const areaVectorField& n = ssf.mesh().faceAreaNormals();
+
+    tmp<GeometricField<Type, faPatchField, areaMesh> > tDiv =
+        fac::edgeIntegrate(ssf);
+
+    GeometricField<Type, faPatchField, areaMesh>& Div = tDiv();
+
+    Div.internalField() = transform(tensor::I - sqr(n), Div.internalField());
+    Div.correctBoundaryConditions();
+
+    return tDiv;
 }
 
 
@@ -79,10 +90,34 @@ div
     const word& name
 )
 {
-    return fa::divScheme<Type>::New
+    const areaVectorField& n = vf.mesh().faceAreaNormals();
+
+    tmp
+    <
+        GeometricField
+        <
+            typename innerProduct<vector, Type>::type,
+            faPatchField,
+            areaMesh
+        >
+    > tDiv
+    (
+        fa::divScheme<Type>::New
         (
             vf.mesh(), vf.mesh().schemesDict().divScheme(name)
-    )().facDiv(vf);
+        )().facDiv(vf)
+    );
+    GeometricField
+    <
+        typename innerProduct<vector, Type>::type,
+        faPatchField,
+        areaMesh
+    >& Div = tDiv();
+
+    Div.internalField() = transform(tensor::I - sqr(n), Div.internalField());
+    Div.correctBoundaryConditions();
+
+    return tDiv;
 }
 
 
@@ -159,12 +194,24 @@ div
     const word& name
 )
 {
-    return fa::convectionScheme<Type>::New
+    const areaVectorField& n = vf.mesh().faceAreaNormals();
+
+    tmp<GeometricField<Type, faPatchField, areaMesh> > tDiv
+    (
+        fa::convectionScheme<Type>::New
         (
             vf.mesh(),
             flux,
             vf.mesh().schemesDict().divScheme(name)
-    )().facDiv(flux, vf);
+        )().facDiv(flux, vf)
+    );
+    GeometricField<Type, faPatchField, areaMesh>& Div = tDiv();
+
+    Div.internalField() = transform(tensor::I - sqr(n), Div.internalField());
+    Div.correctBoundaryConditions();
+
+    return tDiv;
+
 }
 
 

src/finiteArea/finiteArea/fac/facGrad.C

@@ -55,7 +55,18 @@ grad
     const GeometricField<Type, faePatchField, edgeMesh>& ssf
 )
 {
-    return fac::edgeIntegrate(ssf.mesh().Sf() * ssf);
+    const areaVectorField &n = ssf.mesh().faceAreaNormals();
+    typedef typename outerProduct<vector,Type>::type GradType;
+
+    tmp<GeometricField<GradType, faPatchField, areaMesh> > tgGrad =
+        fac::edgeIntegrate(ssf.mesh().Sf()*ssf);
+
+    GeometricField<GradType, faPatchField, areaMesh>& gGrad = tgGrad();
+
+    gGrad -= (gGrad & n)*n;
+    gGrad.correctBoundaryConditions();
+
+    return tgGrad;
 }
 
 template<class Type>
@@ -95,11 +106,22 @@ grad
     const word& name
 )
 {
-    return fa::gradScheme<Type>::New
+    const areaVectorField &n = vf.mesh().faceAreaNormals();
+    typedef typename outerProduct<vector,Type>::type GradType;
+
+    tmp<GeometricField<GradType, faPatchField, areaMesh> > tgGrad =
+        fa::gradScheme<Type>::New
         (
             vf.mesh(),
             vf.mesh().schemesDict().gradScheme(name)
         )().grad(vf);
+
+    GeometricField<GradType, faPatchField, areaMesh>& gGrad = tgGrad();
+
+    gGrad -= (gGrad & n)*n;
+    gGrad.correctBoundaryConditions();
+
+    return tgGrad;
 }
 
 

src/finiteArea/finiteArea/fac/facNDiv.C

@@ -0,0 +1,333 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     3.2
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "facNDiv.H"
+#include "faMesh.H"
+#include "facEdgeIntegrate.H"
+#include "faDivScheme.H"
+#include "faConvectionScheme.H"
+#include "transformField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace fac
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+tmp<GeometricField<Type, faPatchField, areaMesh> >
+ndiv
+(
+    const GeometricField<Type, faePatchField, edgeMesh>& ssf
+)
+{
+    const areaVectorField &n = ssf.mesh().faceAreaNormals();
+
+    tmp<GeometricField<Type, faPatchField, areaMesh> > v = fac::edgeIntegrate(ssf);
+
+    //v.internalField() = transform(n*n, v.internalField());
+    v.internalField() = (v.internalField()&n)*n;
+
+    return v;
+}
+
+
+template<class Type>
+tmp<GeometricField<Type, faPatchField, areaMesh> >
+ndiv
+(
+    const tmp<GeometricField<Type, faePatchField, edgeMesh> >& tssf
+)
+{
+    tmp<GeometricField<Type, faPatchField, areaMesh> > Div(fac::ndiv(tssf()));
+    tssf.clear();
+    return Div;
+}
+
+
+template<class Type>
+tmp
+<
+    GeometricField
+    <
+        typename innerProduct<vector, Type>::type, faPatchField, areaMesh
+    >
+>
+ndiv
+(
+    const GeometricField<Type, faPatchField, areaMesh>& vf,
+    const word& name
+)
+{
+    const areaVectorField &n = vf.mesh().faceAreaNormals();
+
+    tmp<GeometricField<Type, faPatchField, areaMesh> > tDiv
+    (
+        fa::divScheme<Type>::New
+            (
+                vf.mesh(), vf.mesh().schemesDict().divScheme(name)
+            )().facDiv(vf)
+    );
+
+    GeometricField<Type, faPatchField, areaMesh>& Div = tDiv();
+
+    //Div.internalField() = transform(n*n, Div.internalField());
+    Div.internalField() = (Div.internalField()&n)*n;
+    return tDiv;
+}
+
+
+template<class Type>
+tmp
+<
+    GeometricField
+    <
+        typename innerProduct<vector, Type>::type, faPatchField, areaMesh
+    >
+>
+ndiv
+(
+    const tmp<GeometricField<Type, faPatchField, areaMesh> >& tvvf,
+    const word& name
+)
+{
+    typedef typename innerProduct<vector, Type>::type DivType;
+    tmp<GeometricField<DivType, faPatchField, areaMesh> > Div
+    (
+        fac::ndiv(tvvf(), name)
+    );
+    tvvf.clear();
+    return Div;
+}
+
+template<class Type>
+tmp
+<
+    GeometricField
+    <
+        typename innerProduct<vector, Type>::type, faPatchField, areaMesh
+    >
+>
+ndiv
+(
+    const GeometricField<Type, faPatchField, areaMesh>& vf
+)
+{
+    return fac::ndiv(vf, "div("+vf.name()+')');
+}
+
+
+template<class Type>
+tmp
+<
+    GeometricField
+    <
+        typename innerProduct<vector, Type>::type, faPatchField, areaMesh
+    >
+>
+ndiv
+(
+    const tmp<GeometricField<Type, faPatchField, areaMesh> >& tvvf
+)
+{
+    typedef typename innerProduct<vector, Type>::type DivType;
+    tmp<GeometricField<DivType, faPatchField, areaMesh> > Div
+    (
+        fac::ndiv(tvvf())
+    );
+
+    tvvf.clear();
+    return Div;
+}
+
+
+template<class Type>
+tmp<GeometricField<Type, faPatchField, areaMesh> >
+ndiv
+(
+    const edgeScalarField& flux,
+    const GeometricField<Type, faPatchField, areaMesh>& vf,
+    const word& name
+)
+{
+    const areaVectorField &n = vf.mesh().faceAreaNormals();
+
+    tmp<GeometricField<Type, faPatchField, areaMesh> > tDiv
+    (
+        fa::convectionScheme<Type>::New
+            (
+                vf.mesh(),
+                flux,
+                vf.mesh().schemesDict().divScheme(name)
+            )().facDiv(flux, vf)
+    );
+
+    GeometricField<Type, faPatchField, areaMesh>& Div = tDiv();
+
+    //Div.internalField() = transform(n*n, Div.internalField());
+    Div.internalField() = (Div.internalField()&n)*n;
+
+    return tDiv;
+
+}
+
+
+template<class Type>
+tmp<GeometricField<Type, faPatchField, areaMesh> >
+ndiv
+(
+    const tmp<edgeScalarField>& tflux,
+    const GeometricField<Type, faPatchField, areaMesh>& vf,
+    const word& name
+)
+{
+    tmp<GeometricField<Type, faPatchField, areaMesh> > Div
+    (
+        fac::ndiv(tflux(), vf, name)
+    );
+    tflux.clear();
+    return Div;
+}
+
+
+template<class Type>
+tmp<GeometricField<Type, faPatchField, areaMesh> >
+ndiv
+(
+    const edgeScalarField& flux,
+    const tmp<GeometricField<Type, faPatchField, areaMesh> >& tvf,
+    const word& name
+)
+{
+    tmp<GeometricField<Type, faPatchField, areaMesh> > Div
+    (
+        fac::ndiv(flux, tvf(), name)
+    );
+    tvf.clear();
+    return Div;
+}
+
+
+template<class Type>
+tmp<GeometricField<Type, faPatchField, areaMesh> >
+ndiv
+(
+    const tmp<edgeScalarField>& tflux,
+    const tmp<GeometricField<Type, faPatchField, areaMesh> >& tvf,
+    const word& name
+)
+{
+    tmp<GeometricField<Type, faPatchField, areaMesh> > Div
+    (
+        fac::ndiv(tflux(), tvf(), name)
+    );
+    tflux.clear();
+    tvf.clear();
+    return Div;
+}
+
+
+template<class Type>
+tmp<GeometricField<Type, faPatchField, areaMesh> >
+ndiv
+(
+    const edgeScalarField& flux,
+    const GeometricField<Type, faPatchField, areaMesh>& vf
+)
+{
+    return fac::ndiv
+    (
+        flux, vf, "div("+flux.name()+','+vf.name()+')'
+    );
+}
+
+
+template<class Type>
+tmp<GeometricField<Type, faPatchField, areaMesh> >
+ndiv
+(
+    const tmp<edgeScalarField>& tflux,
+    const GeometricField<Type, faPatchField, areaMesh>& vf
+)
+{
+    tmp<GeometricField<Type, faPatchField, areaMesh> > Div
+    (
+        fac::ndiv(tflux(), vf)
+    );
+    tflux.clear();
+    return Div;
+}
+
+
+template<class Type>
+tmp<GeometricField<Type, faPatchField, areaMesh> >
+ndiv
+(
+    const edgeScalarField& flux,
+    const tmp<GeometricField<Type, faPatchField, areaMesh> >& tvf
+)
+{
+    tmp<GeometricField<Type, faPatchField, areaMesh> > Div
+    (
+        fac::ndiv(flux, tvf())
+    );
+    tvf.clear();
+    return Div;
+}
+
+
+template<class Type>
+tmp<GeometricField<Type, faPatchField, areaMesh> >
+ndiv
+(
+    const tmp<edgeScalarField>& tflux,
+    const tmp<GeometricField<Type, faPatchField, areaMesh> >& tvf
+)
+{
+    tmp<GeometricField<Type, faPatchField, areaMesh> > Div
+    (
+        fac::ndiv(tflux(), tvf())
+    );
+    tflux.clear();
+    tvf.clear();
+    return Div;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fac
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/finiteArea/finiteArea/fac/facNDiv.H

@@ -0,0 +1,187 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     3.2
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Namespace
+    fac
+
+Description
+    Calculate the divergence of the given field.
+
+SourceFiles
+    facDiv.C
+
+\*---------------------------------------------------------------------------*/
+
+
+#ifndef facNDiv_H
+#define facNDiv_H
+
+#include "areaFieldsFwd.H"
+#include "edgeFieldsFwd.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                     Namespace fac functions Declaration
+\*---------------------------------------------------------------------------*/
+
+namespace fac
+{
+    template<class Type>
+    tmp<GeometricField<Type, faPatchField, areaMesh> > ndiv
+    (
+        const GeometricField<Type, faePatchField, edgeMesh>&
+    );
+
+    template<class Type>
+    tmp<GeometricField<Type, faPatchField, areaMesh> > ndiv
+    (
+        const tmp<GeometricField<Type, faePatchField, edgeMesh> >&
+    );
+
+
+    template<class Type>
+    tmp
+    <
+        GeometricField
+        <typename innerProduct<vector, Type>::type, faPatchField, areaMesh>
+    > ndiv
+    (
+        const GeometricField<Type, faPatchField, areaMesh>&,
+        const word& name
+    );
+
+    template<class Type>
+    tmp
+    <
+        GeometricField
+        <typename innerProduct<vector, Type>::type, faPatchField, areaMesh>
+    > ndiv
+    (
+        const tmp<GeometricField<Type, faPatchField, areaMesh> >&,
+        const word& name
+    );
+
+
+    template<class Type>
+    tmp
+    <
+        GeometricField
+        <typename innerProduct<vector, Type>::type, faPatchField, areaMesh>
+    > ndiv
+    (
+        const GeometricField<Type, faPatchField, areaMesh>&
+    );
+
+    template<class Type>
+    tmp
+    <
+        GeometricField
+        <typename innerProduct<vector, Type>::type, faPatchField, areaMesh>
+    > ndiv
+    (
+        const tmp<GeometricField<Type, faPatchField, areaMesh> >&
+    );
+
+
+    template<class Type>
+    tmp<GeometricField<Type, faPatchField, areaMesh> > ndiv
+    (
+        const edgeScalarField&,
+        const GeometricField<Type, faPatchField, areaMesh>&,
+        const word& name
+    );
+
+    template<class Type>
+    tmp<GeometricField<Type, faPatchField, areaMesh> > ndiv
+    (
+        const tmp<edgeScalarField>&,
+        const GeometricField<Type, faPatchField, areaMesh>&,
+        const word& name
+    );
+
+    template<class Type>
+    tmp<GeometricField<Type, faPatchField, areaMesh> > ndiv
+    (
+        const edgeScalarField&,
+        const tmp<GeometricField<Type, faPatchField, areaMesh> >&,
+        const word& name
+    );
+
+    template<class Type>
+    tmp<GeometricField<Type, faPatchField, areaMesh> > ndiv
+    (
+        const tmp<edgeScalarField>&,
+        const tmp<GeometricField<Type, faPatchField, areaMesh> >&,
+        const word& name
+    );
+
+
+    template<class Type>
+    tmp<GeometricField<Type, faPatchField, areaMesh> > ndiv
+    (
+        const edgeScalarField&,
+        const GeometricField<Type, faPatchField, areaMesh>&
+    );
+
+    template<class Type>
+    tmp<GeometricField<Type, faPatchField, areaMesh> > ndiv
+    (
+        const tmp<edgeScalarField>&,
+        const GeometricField<Type, faPatchField, areaMesh>&
+    );
+
+    template<class Type>
+    tmp<GeometricField<Type, faPatchField, areaMesh> > ndiv
+    (
+        const edgeScalarField&,
+        const tmp<GeometricField<Type, faPatchField, areaMesh> >&
+    );
+
+    template<class Type>
+    tmp<GeometricField<Type, faPatchField, areaMesh> > ndiv
+    (
+        const tmp<edgeScalarField>&,
+        const tmp<GeometricField<Type, faPatchField, areaMesh> >&
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "facNDiv.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteArea/finiteArea/fac/facNGrad.C

@@ -0,0 +1,202 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     3.2
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "facNGrad.H"
+#include "areaFields.H"
+#include "edgeFields.H"
+#include "facEdgeIntegrate.H"
+#include "faMesh.H"
+#include "faGradScheme.H"
+#include "transformField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace fac
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+tmp
+<
+    GeometricField
+    <
+        typename outerProduct<vector, Type>::type, faPatchField, areaMesh
+    >
+>
+ngrad
+(
+    const GeometricField<Type, faePatchField, edgeMesh>& ssf
+)
+{
+    const areaVectorField &n = ssf.mesh().faceAreaNormals();
+    typedef typename outerProduct<vector,Type>::type GradType;
+
+    tmp<GeometricField<GradType, faPatchField, areaMesh> > tgGrad = fac::edgeIntegrate(ssf.mesh().Sf() * ssf);
+
+    GeometricField<GradType, faPatchField, areaMesh>& grad = tgGrad();
+
+    grad = (grad&n)*n;
+
+    return tgGrad;
+}
+
+template<class Type>
+tmp
+<
+    GeometricField
+    <
+        typename outerProduct<vector,Type>::type, faPatchField, areaMesh
+    >
+>
+ngrad
+(
+    const tmp<GeometricField<Type, faePatchField, edgeMesh> >& tssf
+)
+{
+    typedef typename outerProduct<vector, Type>::type GradType;
+    tmp<GeometricField<GradType, faPatchField, areaMesh> > Grad
+    (
+        fac::ngrad(tssf())
+    );
+    tssf.clear();
+    return Grad;
+}
+
+
+template<class Type>
+tmp
+<
+    GeometricField
+    <
+        typename outerProduct<vector,Type>::type, faPatchField, areaMesh
+    >
+>
+ngrad
+(
+    const GeometricField<Type, faPatchField, areaMesh>& vf,
+    const word& name
+)
+{
+    const areaVectorField &n = vf.mesh().faceAreaNormals();
+    typedef typename outerProduct<vector,Type>::type GradType;
+
+    tmp<GeometricField<GradType, faPatchField, areaMesh> > tgGrad = fa::gradScheme<Type>::New
+        (
+            vf.mesh(),
+            vf.mesh().schemesDict().gradScheme(name)
+        )().grad(vf);
+
+    GeometricField<GradType, faPatchField, areaMesh>& grad = tgGrad();
+
+    grad = (grad&n)*n;
+
+    return tgGrad;
+}
+
+
+template<class Type>
+tmp
+<
+    GeometricField
+    <
+        typename outerProduct<vector,Type>::type, faPatchField, areaMesh
+    >
+>
+ngrad
+(
+    const tmp<GeometricField<Type, faPatchField, areaMesh> >& tvf,
+    const word& name
+)
+{
+    tmp
+    <
+        GeometricField
+        <
+            typename outerProduct<vector, Type>::type, faPatchField, areaMesh
+        >
+    > tGrad
+    (
+        fac::ngrad(tvf(), name)
+    );
+    tvf.clear();
+    return tGrad;
+}
+
+
+template<class Type>
+tmp
+<
+    GeometricField
+    <
+        typename outerProduct<vector,Type>::type, faPatchField, areaMesh
+    >
+>
+ngrad
+(
+    const GeometricField<Type, faPatchField, areaMesh>& vf
+)
+{
+    return fac::ngrad(vf, "grad(" + vf.name() + ')');
+}
+
+
+template<class Type>
+tmp
+<
+    GeometricField
+    <
+        typename outerProduct<vector,Type>::type, faPatchField, areaMesh
+    >
+>
+ngrad
+(
+    const tmp<GeometricField<Type, faPatchField, areaMesh> >& tvf
+)
+{
+    typedef typename outerProduct<vector, Type>::type GradType;
+    tmp<GeometricField<GradType, faPatchField, areaMesh> > Grad
+    (
+        fac::ngrad(tvf())
+    );
+    tvf.clear();
+    return Grad;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fac
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/finiteArea/finiteArea/fac/facNGrad.H

@@ -0,0 +1,131 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     3.2
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Namespace
+    fac
+
+Description
+    Calculate the gradient normal to the surface of the given field.
+
+SourceFiles
+    facGrad.C
+
+\*---------------------------------------------------------------------------*/
+
+
+#ifndef facNGrad_H
+#define facNGrad_H
+
+#include "areaFieldsFwd.H"
+#include "edgeFieldsFwd.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                     Namespace fac functions Declaration
+\*---------------------------------------------------------------------------*/
+
+namespace fac
+{
+    template<class Type>
+    tmp
+    <
+        GeometricField
+        <typename outerProduct<vector, Type>::type, faPatchField, areaMesh>
+    > ngrad
+    (
+        const GeometricField<Type, faePatchField, edgeMesh>&
+    );
+
+    template<class Type>
+    tmp
+    <
+        GeometricField
+        <typename outerProduct<vector, Type>::type, faPatchField, areaMesh>
+    > ngrad
+    (
+        const tmp<GeometricField<Type, faePatchField, edgeMesh> >&
+    );
+
+    template<class Type>
+    tmp
+    <
+        GeometricField
+        <typename outerProduct<vector, Type>::type, faPatchField, areaMesh>
+    >ngrad
+    (
+        const GeometricField<Type, faPatchField, areaMesh>&,
+        const word& name
+    );
+
+    template<class Type>
+    tmp
+    <
+        GeometricField
+        <typename outerProduct<vector, Type>::type, faPatchField, areaMesh>
+    >ngrad
+    (
+        const tmp<GeometricField<Type, faPatchField, areaMesh> >&,
+        const word& name
+    );
+
+    template<class Type>
+    tmp
+    <
+        GeometricField
+        <typename outerProduct<vector, Type>::type, faPatchField, areaMesh>
+    >ngrad
+    (
+        const GeometricField<Type, faPatchField, areaMesh>&
+    );
+
+    template<class Type>
+    tmp
+    <
+        GeometricField
+        <typename outerProduct<vector, Type>::type, faPatchField, areaMesh>
+    >ngrad
+    (
+        const tmp<GeometricField<Type, faPatchField, areaMesh> >&
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "facNGrad.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteArea/finiteArea/fam/fam.H

@@ -43,6 +43,7 @@ SourceFiles
 #include "famDiv.H"
 #include "famLaplacian.H"
 #include "famSup.H"
+#include "famNDiv.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/finiteArea/finiteArea/fam/famDiv.C

@@ -49,14 +49,36 @@ div
     const word& name
 )
 {
-    return fa::convectionScheme<Type>::New
+    const areaVectorField& n = vf.mesh().faceAreaNormals();
+
+    tmp<faMatrix<Type> > tM
+    (
+        fa::convectionScheme<Type>::New
         (
             vf.mesh(),
             flux,
             vf.mesh().schemesDict().divScheme(name)
-    )().famDiv(flux, vf);
+        )().famDiv(flux, vf)
+    );
+    faMatrix<Type>& M = tM();
+
+    GeometricField<Type, faPatchField, areaMesh> v
+    (
+        fa::convectionScheme<Type>::New
+        (
+            vf.mesh(),
+            flux,
+            vf.mesh().schemesDict().divScheme(name)
+        )().facDiv(flux, vf)
+    );
+
+    //HJ Check if the product is from left or right.  HJ, 6/Dec/2016
+    M -= (v & n)*n;
+
+    return tM;
 }
 
+
 template<class Type>
 tmp<faMatrix<Type> >
 div

src/finiteArea/finiteArea/fam/famDiv.H

@@ -29,6 +29,7 @@ Description
 
 SourceFiles
     famDiv.C
+    vectorFamDiv.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -89,6 +90,9 @@ namespace fam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+// Template specialisation
+#include "vectorFamDiv.H"
+
 #ifdef NoRepository
 #   include "famDiv.C"
 #endif

src/finiteArea/finiteArea/fam/famNDiv.C

@@ -0,0 +1,108 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     3.2
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "famNDiv.H"
+#include "faMesh.H"
+#include "faMatrix.H"
+#include "faConvectionScheme.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace fam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+tmp<faMatrix<Type> >
+ndiv
+(
+    const edgeScalarField& flux,
+    GeometricField<Type, faPatchField, areaMesh>& vf,
+    const word& name
+)
+{
+    return fa::convectionScheme<Type>::New
+    (
+        vf.mesh(),
+        flux,
+        vf.mesh().schemesDict().divScheme(name)
+    )().famDiv(flux, vf);//TODO calculate normal
+}
+
+template<class Type>
+tmp<faMatrix<Type> >
+ndiv
+(
+    const tmp<edgeScalarField>& tflux,
+    GeometricField<Type, faPatchField, areaMesh>& vf,
+    const word& name
+)
+{
+    tmp<faMatrix<Type> > Div(fam::ndiv(tflux(), vf, name));
+    tflux.clear();
+    return Div;
+}
+
+
+template<class Type>
+tmp<faMatrix<Type> >
+ndiv
+(
+    const edgeScalarField& flux,
+    GeometricField<Type, faPatchField, areaMesh>& vf
+)
+{
+    return fam::ndiv(flux, vf, "div("+flux.name()+','+vf.name()+')');
+}
+
+template<class Type>
+tmp<faMatrix<Type> >
+ndiv
+(
+    const tmp<edgeScalarField>& tflux,
+    GeometricField<Type, faPatchField, areaMesh>& vf
+)
+{
+    tmp<faMatrix<Type> > Div(fam::ndiv(tflux(), vf));
+    tflux.clear();
+    return Div;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/finiteArea/finiteArea/fam/famNDiv.H

@@ -0,0 +1,100 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     3.2
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Namespace
+    fam
+
+Description
+    Calculate the matrix for the divergence of the given field and flux.
+
+SourceFiles
+    famDiv.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef famNDiv_H
+#define famNDiv_H
+
+#include "areaFieldsFwd.H"
+#include "edgeFieldsFwd.H"
+#include "edgeInterpolationScheme.H"
+#include "faMatrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                     Namespace fam functions Declaration
+\*---------------------------------------------------------------------------*/
+
+namespace fam
+{
+    template<class Type>
+    tmp<faMatrix<Type> > ndiv
+    (
+        const edgeScalarField&,
+        GeometricField<Type, faPatchField, areaMesh>&,
+        const word& name
+    );
+
+    template<class Type>
+    tmp<faMatrix<Type> > ndiv
+    (
+        const tmp<edgeScalarField>&,
+        GeometricField<Type, faPatchField, areaMesh>&,
+        const word& name
+    );
+
+    template<class Type>
+    tmp<faMatrix<Type> > ndiv
+    (
+        const edgeScalarField&,
+        GeometricField<Type, faPatchField, areaMesh>&
+    );
+
+    template<class Type>
+    tmp<faMatrix<Type> > ndiv
+    (
+        const tmp<edgeScalarField>&,
+        GeometricField<Type, faPatchField, areaMesh>&
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "famNDiv.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteArea/finiteArea/fam/vectorFamDiv.C

@@ -0,0 +1,69 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "vectorFamDiv.H"
+#include "faMesh.H"
+#include "faMatrix.H"
+#include "faConvectionScheme.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace fam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<>
+tmp<faMatrix<scalar> >
+div
+(
+    const edgeScalarField& flux,
+    const GeometricField<scalar, faPatchField, areaMesh>& vf,
+    const word& name
+)
+{
+    return fa::convectionScheme<scalar>::New
+    (
+        vf.mesh(),
+        flux,
+        vf.mesh().schemesDict().divScheme(name)
+    )().famDiv(flux, vf);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/finiteArea/finiteArea/fam/vectorFamDiv.H

@@ -0,0 +1,74 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Typedef
+    Foam::fam
+
+Description
+    Specialisation of fam div for a flux.
+
+SourceFiles
+    vectorFamDiv.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef vectorFamDiv_H
+#define vectorFamDiv_H
+
+#include "famDiv.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+namespace fam
+{
+
+/*---------------------------------------------------------------------------*\
+                     Namespace fam functions Declaration
+\*---------------------------------------------------------------------------*/
+
+template<>
+tmp<faMatrix<scalar> >
+div
+(
+    const edgeScalarField& flux,
+    const GeometricField<scalar, faPatchField, areaMesh>& vf,
+    const word& name
+);
+
+    
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteArea/finiteArea/gradSchemes/gaussFaGrad/gaussFaGrad.C

@@ -66,12 +66,7 @@ gaussGrad<Type>::grad
 
     GeometricField<GradType, faPatchField, areaMesh>& gGrad = tgGrad();
 
-    gGrad -= vsf*fac::edgeIntegrate(vsf.mesh().Le());
-
-    // Remove component of gradient normal to surface (area)
-    const areaVectorField& n = vsf.mesh().faceAreaNormals();
-
-    gGrad -= n*(n & gGrad);
+    // Removed for consistencty.  Matthias Rauter, 6/Dec/2016
     gGrad.correctBoundaryConditions();
 
     gGrad.rename("grad(" + vsf.name() + ')');
@@ -96,15 +91,9 @@ void gaussGrad<Type>::correctBoundaryConditions
         if (!vsf.boundaryField()[patchI].coupled())
         {
             vectorField m =
-                vsf.mesh().Le().boundaryField()[patchI]
-                /vsf.mesh().magLe().boundaryField()[patchI];
+                vsf.mesh().Le().boundaryField()[patchI]/
+                vsf.mesh().magLe().boundaryField()[patchI];
 
-            // Zeljko Tukovic
-//             gGrad.boundaryField()[patchI] =
-//                 m*vsf.boundaryField()[patchI].snGrad();
-
-            //HJ Not sure: should this be a complete correction or just the
-            //   tangential part?  HJ, 24/Jul/2009
             gGrad.boundaryField()[patchI] += m*
             (
                 vsf.boundaryField()[patchI].snGrad()

src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaGrad.C

@@ -138,11 +138,8 @@ leastSquaresFaGrad<Type>::grad
         }
     }
 
-
     // Remove component of gradient normal to surface (area)
-    const areaVectorField& n = vsf.mesh().faceAreaNormals();
-
-    lsGrad -= n*(n & lsGrad);
+    // Removed for consistencty.  Matthias Rauter, 6/Dec/2016
     lsGrad.correctBoundaryConditions();
 
     gaussGrad<Type>::correctBoundaryConditions(vsf, lsGrad);

src/finiteVolume/cfdTools/general/MRF/MRFZone.C

@@ -62,7 +62,7 @@ void Foam::MRFZone::setMRFFaces()
     if (cellZoneID_ != -1)
     {
         const labelList& cellLabels = mesh_.cellZones()[cellZoneID_];
-        forAll(cellLabels, i)
+        forAll (cellLabels, i)
         {
             zoneCell[cellLabels[i]] = true;
         }
@@ -84,15 +84,15 @@ void Foam::MRFZone::setMRFFaces()
 
     labelHashSet excludedPatches(excludedPatchLabels_);
 
-    forAll(patches, patchI)
+    forAll (patches, patchI)
     {
         const polyPatch& pp = patches[patchI];
 
         if (pp.isWall() && !excludedPatches.found(patchI))
         {
-            forAll(pp, i)
+            forAll (pp, i)
             {
-                label faceI = pp.start()+i;
+                label faceI = pp.start() + i;
 
                 if (zoneCell[own[faceI]])
                 {
@@ -103,9 +103,9 @@ void Foam::MRFZone::setMRFFaces()
         }
         else if (!isA<emptyPolyPatch>(pp))
         {
-            forAll(pp, i)
+            forAll (pp, i)
             {
-                label faceI = pp.start()+i;
+                label faceI = pp.start() + i;
 
                 if (zoneCell[own[faceI]])
                 {
@@ -138,11 +138,11 @@ void Foam::MRFZone::setMRFFaces()
     labelList nIncludedFaces(patches.size(), 0);
     labelList nExcludedFaces(patches.size(), 0);
 
-    forAll(patches, patchi)
+    forAll (patches, patchi)
     {
         const polyPatch& pp = patches[patchi];
 
-        forAll(pp, patchFacei)
+        forAll (pp, patchFacei)
         {
             label faceI = pp.start()+patchFacei;
 
@@ -159,7 +159,7 @@ void Foam::MRFZone::setMRFFaces()
 
     includedFaces_.setSize(patches.size());
     excludedFaces_.setSize(patches.size());
-    forAll(nIncludedFaces, patchi)
+    forAll (nIncludedFaces, patchi)
     {
         includedFaces_[patchi].setSize(nIncludedFaces[patchi]);
         excludedFaces_[patchi].setSize(nExcludedFaces[patchi]);
@@ -167,11 +167,11 @@ void Foam::MRFZone::setMRFFaces()
     nIncludedFaces = 0;
     nExcludedFaces = 0;
 
-    forAll(patches, patchi)
+    forAll (patches, patchi)
     {
         const polyPatch& pp = patches[patchi];
 
-        forAll(pp, patchFacei)
+        forAll (pp, patchFacei)
         {
             label faceI = pp.start()+patchFacei;
 
@@ -201,9 +201,9 @@ void Foam::MRFZone::setMRFFaces()
         internalFaces.write();
 
         faceSet MRFFaces(mesh_, "includedFaces", 100);
-        forAll(includedFaces_, patchi)
+        forAll (includedFaces_, patchi)
         {
-            forAll(includedFaces_[patchi], i)
+            forAll (includedFaces_[patchi], i)
             {
                 label patchFacei = includedFaces_[patchi][i];
                 MRFFaces.insert(patches[patchi].start()+patchFacei);
@@ -215,12 +215,12 @@ void Foam::MRFZone::setMRFFaces()
         MRFFaces.write();
 
         faceSet excludedFaces(mesh_, "excludedFaces", 100);
-        forAll(excludedFaces_, patchi)
+        forAll (excludedFaces_, patchi)
         {
-            forAll(excludedFaces_[patchi], i)
+            forAll (excludedFaces_[patchi], i)
             {
                 label patchFacei = excludedFaces_[patchi][i];
-                excludedFaces.insert(patches[patchi].start()+patchFacei);
+                excludedFaces.insert(patches[patchi].start() + patchFacei);
             }
         }
         Pout<< "Writing " << excludedFaces.size()
@@ -231,12 +231,30 @@ void Foam::MRFZone::setMRFFaces()
 }
 
 
+Foam::vector Foam::MRFZone::Omega() const
+{
+    if (rampTime_ < SMALL)
+    {
+        return omega_.value()*axis_.value();
+    }
+    else
+    {
+        // Ramping
+        Info<< "MRF Ramp: " << Foam::min(mesh_.time().value()/rampTime_, 1)
+            << endl;
+
+        return Foam::min(mesh_.time().value()/rampTime_, 1)*
+            omega_.value()*axis_.value();
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::MRFZone::MRFZone(const fvMesh& mesh, Istream& is)
 :
-    mesh_(mesh),
     name_(is),
+    mesh_(mesh),
     dict_(is),
     cellZoneID_(mesh_.cellZones().findZoneID(name_)),
     excludedPatchNames_
@@ -246,7 +264,7 @@ Foam::MRFZone::MRFZone(const fvMesh& mesh, Istream& is)
     origin_(dict_.lookup("origin")),
     axis_(dict_.lookup("axis")),
     omega_(dict_.lookup("omega")),
-    Omega_("Omega", omega_*axis_)
+    rampTime_(dict_.lookupOrDefault<scalar>("rampTime", 0))
 {
     if (dict_.found("patches"))
     {
@@ -260,12 +278,19 @@ Foam::MRFZone::MRFZone(const fvMesh& mesh, Istream& is)
 
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
+    if (mag(axis_.value()) < SMALL)
+    {
+        FatalErrorIn("MRFZone(const fvMesh&, Istream&)")
+            << "Axis vector has zero magnitude: " << axis_
+            << ".  This is not allowed"
+            << abort(FatalError);
+    }
+
     axis_ = axis_/mag(axis_);
-    Omega_ = omega_*axis_;
 
     excludedPatchLabels_.setSize(excludedPatchNames_.size());
 
-    forAll(excludedPatchNames_, i)
+    forAll (excludedPatchNames_, i)
     {
         excludedPatchLabels_[i] = patches.findPatchID(excludedPatchNames_[i]);
 
@@ -309,12 +334,12 @@ void Foam::MRFZone::addCoriolis(fvVectorMatrix& UEqn) const
     const scalarField& V = mesh_.V();
     vectorField& Usource = UEqn.source();
     const vectorField& U = UEqn.psi();
-    const vector& Omega = Omega_.value();
+    const vector rotVel = Omega();
 
-    forAll(cells, i)
+    forAll (cells, i)
     {
         label celli = cells[i];
-        Usource[celli] -= V[celli]*(Omega ^ U[celli]);
+        Usource[celli] -= V[celli]*(rotVel ^ U[celli]);
     }
 }
 
@@ -334,12 +359,12 @@ void Foam::MRFZone::addCoriolis
     const scalarField& V = mesh_.V();
     vectorField& Usource = UEqn.source();
     const vectorField& U = UEqn.psi();
-    const vector& Omega = Omega_.value();
+    const vector rotVel = Omega();
 
-    forAll(cells, i)
+    forAll (cells, i)
     {
         label celli = cells[i];
-        Usource[celli] -= V[celli]*rho[celli]*(Omega ^ U[celli]);
+        Usource[celli] -= V[celli]*rho[celli]*(rotVel ^ U[celli]);
     }
 }
 
@@ -349,20 +374,20 @@ void Foam::MRFZone::relativeVelocity(volVectorField& U) const
     const volVectorField& C = mesh_.C();
 
     const vector& origin = origin_.value();
-    const vector& Omega = Omega_.value();
+    const vector rotVel = Omega();
 
     const labelList& cells = mesh_.cellZones()[cellZoneID_];
 
-    forAll(cells, i)
+    forAll (cells, i)
     {
         label celli = cells[i];
-        U[celli] -= (Omega ^ (C[celli] - origin));
+        U[celli] -= (rotVel ^ (C[celli] - origin));
     }
 
     // Included patches
-    forAll(includedFaces_, patchi)
+    forAll (includedFaces_, patchi)
     {
-        forAll(includedFaces_[patchi], i)
+        forAll (includedFaces_[patchi], i)
         {
             label patchFacei = includedFaces_[patchi][i];
             U.boundaryField()[patchi][patchFacei] = vector::zero;
@@ -370,13 +395,13 @@ void Foam::MRFZone::relativeVelocity(volVectorField& U) const
     }
 
     // Excluded patches
-    forAll(excludedFaces_, patchi)
+    forAll (excludedFaces_, patchi)
     {
-        forAll(excludedFaces_[patchi], i)
+        forAll (excludedFaces_[patchi], i)
         {
             label patchFacei = excludedFaces_[patchi][i];
             U.boundaryField()[patchi][patchFacei] -=
-                (Omega ^ (C.boundaryField()[patchi][patchFacei] - origin));
+                (rotVel ^ (C.boundaryField()[patchi][patchFacei] - origin));
         }
     }
 }
@@ -387,35 +412,35 @@ void Foam::MRFZone::absoluteVelocity(volVectorField& U) const
     const volVectorField& C = mesh_.C();
 
     const vector& origin = origin_.value();
-    const vector& Omega = Omega_.value();
+    const vector rotVel = Omega();
 
     const labelList& cells = mesh_.cellZones()[cellZoneID_];
 
-    forAll(cells, i)
+    forAll (cells, i)
     {
         label celli = cells[i];
-        U[celli] += (Omega ^ (C[celli] - origin));
+        U[celli] += (rotVel ^ (C[celli] - origin));
     }
 
     // Included patches
-    forAll(includedFaces_, patchi)
+    forAll (includedFaces_, patchi)
     {
-        forAll(includedFaces_[patchi], i)
+        forAll (includedFaces_[patchi], i)
         {
             label patchFacei = includedFaces_[patchi][i];
             U.boundaryField()[patchi][patchFacei] =
-                (Omega ^ (C.boundaryField()[patchi][patchFacei] - origin));
+                (rotVel ^ (C.boundaryField()[patchi][patchFacei] - origin));
         }
     }
 
     // Excluded patches
-    forAll(excludedFaces_, patchi)
+    forAll (excludedFaces_, patchi)
     {
-        forAll(excludedFaces_[patchi], i)
+        forAll (excludedFaces_[patchi], i)
         {
             label patchFacei = excludedFaces_[patchi][i];
             U.boundaryField()[patchi][patchFacei] +=
-                (Omega ^ (C.boundaryField()[patchi][patchFacei] - origin));
+                (rotVel ^ (C.boundaryField()[patchi][patchFacei] - origin));
         }
     }
 }
@@ -460,36 +485,36 @@ void Foam::MRFZone::meshPhi
     const surfaceVectorField& Cf = mesh_.Cf();
 
     const vector& origin = origin_.value();
-    const vector& Omega = Omega_.value();
+    const vector rotVel = Omega();
 
     // Internal faces
-    forAll(internalFaces_, i)
+    forAll (internalFaces_, i)
     {
         label facei = internalFaces_[i];
-        phi[facei] = (Omega ^ (Cf[facei] - origin));
+        phi[facei] = (rotVel ^ (Cf[facei] - origin));
     }
 
     // Included patches
-    forAll(includedFaces_, patchi)
+    forAll (includedFaces_, patchi)
     {
-        forAll(includedFaces_[patchi], i)
+        forAll (includedFaces_[patchi], i)
         {
             label patchFacei = includedFaces_[patchi][i];
 
             phi.boundaryField()[patchi][patchFacei] =
-                (Omega ^ (Cf.boundaryField()[patchi][patchFacei] - origin));
+                (rotVel ^ (Cf.boundaryField()[patchi][patchFacei] - origin));
         }
     }
 
     // Excluded patches
-    forAll(excludedFaces_, patchi)
+    forAll (excludedFaces_, patchi)
     {
-        forAll(excludedFaces_[patchi], i)
+        forAll (excludedFaces_[patchi], i)
         {
             label patchFacei = excludedFaces_[patchi][i];
 
             phi.boundaryField()[patchi][patchFacei] =
-                (Omega ^ (Cf.boundaryField()[patchi][patchFacei] - origin));
+                (rotVel ^ (Cf.boundaryField()[patchi][patchFacei] - origin));
         }
     }
 }
@@ -497,20 +522,20 @@ void Foam::MRFZone::meshPhi
 void Foam::MRFZone::correctBoundaryVelocity(volVectorField& U) const
 {
     const vector& origin = origin_.value();
-    const vector& Omega = Omega_.value();
+    const vector rotVel = Omega();
 
     // Included patches
-    forAll(includedFaces_, patchi)
+    forAll (includedFaces_, patchi)
     {
         const vectorField& patchC = mesh_.Cf().boundaryField()[patchi];
 
         vectorField pfld(U.boundaryField()[patchi]);
 
-        forAll(includedFaces_[patchi], i)
+        forAll (includedFaces_[patchi], i)
         {
             label patchFacei = includedFaces_[patchi][i];
 
-            pfld[patchFacei] = (Omega ^ (patchC[patchFacei] - origin));
+            pfld[patchFacei] = (rotVel ^ (patchC[patchFacei] - origin));
         }
 
         U.boundaryField()[patchi] == pfld;
@@ -529,12 +554,12 @@ void Foam::MRFZone::Su
     const volVectorField& C = mesh_.C();
 
     const vector& origin = origin_.value();
-    const vector& Omega = Omega_.value();
+    const vector rotVel = Omega();
 
-    forAll(cells, i)
+    forAll (cells, i)
     {
         source[cells[i]] =
-            (Omega ^ (C[cells[i]] - origin)) & gradPhi[cells[i]];
+            (rotVel ^ (C[cells[i]] - origin)) & gradPhi[cells[i]];
     }
 }
 
@@ -550,13 +575,13 @@ void Foam::MRFZone::Su
     const volVectorField& C = mesh_.C();
 
     const vector& origin = origin_.value();
-    const vector& Omega = Omega_.value();
+    const vector rotVel = Omega();
 
-    forAll(cells, i)
+    forAll (cells, i)
     {
         source[cells[i]] =
-            ((Omega ^ (C[cells[i]] - origin)) & gradPhi[cells[i]])
-          - (Omega ^ phi[cells[i]]);
+            ((rotVel ^ (C[cells[i]] - origin)) & gradPhi[cells[i]])
+          - (rotVel ^ phi[cells[i]]);
     }
 }
 

src/finiteVolume/cfdTools/general/MRF/MRFZone.H

@@ -29,7 +29,7 @@ Description
     obtained from a control dictionary constructed from the given stream.
 
     The rotation of the MRF region is defined by an origin and axis of
-    rotation and an angular speed.
+    rotation and an angular speed in rad/s.
 
 SourceFiles
     MRFZone.C
@@ -65,15 +65,22 @@ class MRFZone
 {
     // Private data
 
-        const fvMesh& mesh_;
-
+        //- Name of the MRF zone
         const word name_;
 
+        //- Mesh reference
+        const fvMesh& mesh_;
+
+        //- Reference to dictionary
         const dictionary dict_;
 
+        //- MRF cell zone ID: cells in MRF region
         label cellZoneID_;
 
+        //- List of patch names touching the MRF zone which do not move
         const wordList excludedPatchNames_;
+
+        //- List of patch labels touching the MRF zone which do not move
         labelList excludedPatchLabels_;
 
         //- Internal faces that are part of MRF
@@ -85,17 +92,34 @@ class MRFZone
         //- Excluded faces (per patch) that do not move with the MRF
         labelListList excludedFaces_;
 
+        //- Origin of rotation
         const dimensionedVector origin_;
+
+        //- Axis of rotation
         dimensionedVector axis_;
-        const dimensionedScalar omega_;
-        dimensionedVector Omega_;
+
+        //- Rotational velocity in rad/s
+        dimensionedScalar omega_;
+
+        //- Ramping time scale
+        scalar rampTime_;
 
 
     // Private Member Functions
 
+        //- Disallow default bitwise copy construct
+        MRFZone(const MRFZone&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const MRFZone&);
+
+
         //- Divide faces in frame according to patch
         void setMRFFaces();
 
+        //- Return rotational vector
+        vector Omega() const;
+
         //- Make the given absolute mass/vol flux relative within MRF region
         template<class RhoFieldType>
         void relativeRhoFlux
@@ -112,12 +136,6 @@ class MRFZone
             surfaceScalarField& phi
         ) const;
 
-        //- Disallow default bitwise copy construct
-        MRFZone(const MRFZone&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const MRFZone&);
-
 
 public:
 
@@ -170,7 +188,11 @@ public:
         void addCoriolis(fvVectorMatrix& UEqn) const;
 
         //- Add the Coriolis force contribution to the momentum equation
-        void addCoriolis(const volScalarField& rho, fvVectorMatrix& UEqn) const;
+        void addCoriolis
+        (
+            const volScalarField& rho,
+            fvVectorMatrix& UEqn
+        ) const;
 
         //- Make the given absolute velocity relative within the MRF region
         void relativeVelocity(volVectorField& U) const;

src/finiteVolume/cfdTools/general/MRF/MRFZoneTemplates.C

@@ -41,19 +41,19 @@ void Foam::MRFZone::relativeRhoFlux
     const surfaceVectorField& Sf = mesh_.Sf();
 
     const vector& origin = origin_.value();
-    const vector& Omega = Omega_.value();
+    const vector rotVel = Omega();
 
     // Internal faces
-    forAll(internalFaces_, i)
+    forAll (internalFaces_, i)
     {
         label facei = internalFaces_[i];
-        phi[facei] -= rho[facei]*(Omega ^ (Cf[facei] - origin)) & Sf[facei];
+        phi[facei] -= rho[facei]*(rotVel ^ (Cf[facei] - origin)) & Sf[facei];
     }
 
     // Included patches
-    forAll(includedFaces_, patchi)
+    forAll (includedFaces_, patchi)
     {
-        forAll(includedFaces_[patchi], i)
+        forAll (includedFaces_[patchi], i)
         {
             label patchFacei = includedFaces_[patchi][i];
 
@@ -62,15 +62,15 @@ void Foam::MRFZone::relativeRhoFlux
     }
 
     // Excluded patches
-    forAll(excludedFaces_, patchi)
+    forAll (excludedFaces_, patchi)
     {
-        forAll(excludedFaces_[patchi], i)
+        forAll (excludedFaces_[patchi], i)
         {
             label patchFacei = excludedFaces_[patchi][i];
 
             phi.boundaryField()[patchi][patchFacei] -=
                 rho.boundaryField()[patchi][patchFacei]
-               *(Omega ^ (Cf.boundaryField()[patchi][patchFacei] - origin))
+               *(rotVel ^ (Cf.boundaryField()[patchi][patchFacei] - origin))
               & Sf.boundaryField()[patchi][patchFacei];
         }
     }
@@ -88,37 +88,37 @@ void Foam::MRFZone::absoluteRhoFlux
     const surfaceVectorField& Sf = mesh_.Sf();
 
     const vector& origin = origin_.value();
-    const vector& Omega = Omega_.value();
+    const vector rotVel = Omega();
 
     // Internal faces
-    forAll(internalFaces_, i)
+    forAll (internalFaces_, i)
     {
         label facei = internalFaces_[i];
-        phi[facei] += (Omega ^ (Cf[facei] - origin)) & Sf[facei];
+        phi[facei] += (rotVel ^ (Cf[facei] - origin)) & Sf[facei];
     }
 
     // Included patches
-    forAll(includedFaces_, patchi)
+    forAll (includedFaces_, patchi)
     {
-        forAll(includedFaces_[patchi], i)
+        forAll (includedFaces_[patchi], i)
         {
             label patchFacei = includedFaces_[patchi][i];
 
             phi.boundaryField()[patchi][patchFacei] +=
-                (Omega ^ (Cf.boundaryField()[patchi][patchFacei] - origin))
+                (rotVel ^ (Cf.boundaryField()[patchi][patchFacei] - origin))
               & Sf.boundaryField()[patchi][patchFacei];
         }
     }
 
     // Excluded patches
-    forAll(excludedFaces_, patchi)
+    forAll (excludedFaces_, patchi)
     {
-        forAll(excludedFaces_[patchi], i)
+        forAll (excludedFaces_[patchi], i)
         {
             label patchFacei = excludedFaces_[patchi][i];
 
             phi.boundaryField()[patchi][patchFacei] +=
-                (Omega ^ (Cf.boundaryField()[patchi][patchFacei] - origin))
+                (rotVel ^ (Cf.boundaryField()[patchi][patchFacei] - origin))
               & Sf.boundaryField()[patchi][patchFacei];
         }
     }

src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C

@@ -62,7 +62,7 @@ void Foam::setRefCell
                          ")",
                         dict
                     )   << "Illegal master cellID " << refCelli
-                        << ". Should be 0.." << field.mesh().nCells()
+                        << ". Should be 0." << field.mesh().nCells()
                         << exit(FatalIOError);
                 }
             }

src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchField.C

@@ -31,27 +31,28 @@ License
 #include "coeffFields.H"
 #include "demandDrivenData.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
 
 template<class Type>
-processorFvPatchField<Type>::processorFvPatchField
+Foam::processorFvPatchField<Type>::processorFvPatchField
 (
     const fvPatch& p,
     const DimensionedField<Type, volMesh>& iF
 )
 :
     coupledFvPatchField<Type>(p, iF),
-    procPatch_(refCast<const processorFvPatch>(p))
+    procPatch_(refCast<const processorFvPatch>(p)),
+    outstandingSendRequest_(-1),
+    outstandingRecvRequest_(-1),
+    sendBuf_(0),
+    receiveBuf_(0),
+    scalarSendBuf_(0),
+    scalarReceiveBuf_(0)
 {}
 
 
 template<class Type>
-processorFvPatchField<Type>::processorFvPatchField
+Foam::processorFvPatchField<Type>::processorFvPatchField
 (
     const fvPatch& p,
     const DimensionedField<Type, volMesh>& iF,
@@ -59,13 +60,18 @@ processorFvPatchField<Type>::processorFvPatchField
 )
 :
     coupledFvPatchField<Type>(p, iF, f),
-    procPatch_(refCast<const processorFvPatch>(p))
+    procPatch_(refCast<const processorFvPatch>(p)),
+    outstandingSendRequest_(-1),
+    outstandingRecvRequest_(-1),
+    sendBuf_(0),
+    receiveBuf_(0),
+    scalarSendBuf_(0),
+    scalarReceiveBuf_(0)
 {}
 
 
-// Construct by mapping given processorFvPatchField<Type>
 template<class Type>
-processorFvPatchField<Type>::processorFvPatchField
+Foam::processorFvPatchField<Type>::processorFvPatchField
 (
     const processorFvPatchField<Type>& ptf,
     const fvPatch& p,
@@ -74,9 +80,15 @@ processorFvPatchField<Type>::processorFvPatchField
 )
 :
     coupledFvPatchField<Type>(ptf, p, iF, mapper),
-    procPatch_(refCast<const processorFvPatch>(p))
+    procPatch_(refCast<const processorFvPatch>(p)),
+    outstandingSendRequest_(-1),
+    outstandingRecvRequest_(-1),
+    sendBuf_(0),
+    receiveBuf_(0),
+    scalarSendBuf_(0),
+    scalarReceiveBuf_(0)
 {
-    if (!isType<processorFvPatch>(this->patch()))
+    if (!isA<processorFvPatch>(this->patch()))
     {
         FatalErrorIn
         (
@@ -98,7 +110,7 @@ processorFvPatchField<Type>::processorFvPatchField
 
 
 template<class Type>
-processorFvPatchField<Type>::processorFvPatchField
+Foam::processorFvPatchField<Type>::processorFvPatchField
 (
     const fvPatch& p,
     const DimensionedField<Type, volMesh>& iF,
@@ -106,9 +118,15 @@ processorFvPatchField<Type>::processorFvPatchField
 )
 :
     coupledFvPatchField<Type>(p, iF, dict),
-    procPatch_(refCast<const processorFvPatch>(p))
+    procPatch_(refCast<const processorFvPatch>(p)),
+    outstandingSendRequest_(-1),
+    outstandingRecvRequest_(-1),
+    sendBuf_(0),
+    receiveBuf_(0),
+    scalarSendBuf_(0),
+    scalarReceiveBuf_(0)
 {
-    if (!isType<processorFvPatch>(p))
+    if (!isA<processorFvPatch>(p))
     {
         FatalIOErrorIn
         (
@@ -130,42 +148,79 @@ processorFvPatchField<Type>::processorFvPatchField
 
 
 template<class Type>
-processorFvPatchField<Type>::processorFvPatchField
+Foam::processorFvPatchField<Type>::processorFvPatchField
 (
     const processorFvPatchField<Type>& ptf
 )
 :
     processorLduInterfaceField(),
     coupledFvPatchField<Type>(ptf),
-    procPatch_(refCast<const processorFvPatch>(ptf.patch()))
+    procPatch_(refCast<const processorFvPatch>(ptf.patch())),
+    outstandingSendRequest_(-1),
+    outstandingRecvRequest_(-1),
+    sendBuf_(0),
+    receiveBuf_(0),
+    scalarSendBuf_(0),
+    scalarReceiveBuf_(0)
 {}
 
 
 template<class Type>
-processorFvPatchField<Type>::processorFvPatchField
+Foam::processorFvPatchField<Type>::processorFvPatchField
 (
     const processorFvPatchField<Type>& ptf,
     const DimensionedField<Type, volMesh>& iF
 )
 :
     coupledFvPatchField<Type>(ptf, iF),
-    procPatch_(refCast<const processorFvPatch>(ptf.patch()))
-{}
+    procPatch_(refCast<const processorFvPatch>(ptf.patch())),
+    outstandingSendRequest_(-1),
+    outstandingRecvRequest_(-1),
+    sendBuf_(0),
+    receiveBuf_(0),
+    scalarSendBuf_(0),
+    scalarReceiveBuf_(0)
+{
+    if (debug && !ptf.ready())
+    {
+        FatalErrorIn
+        (
+            "processorFvPatchField<Type>::processorFvPatchField\n"
+            "(\n"
+            "    const processorFvPatchField<Type>& ptf,\n"
+            "    const DimensionedField<Type, volMesh>& iF\n"
+            ")"
+        )   << "On patch " << procPatch_.name() << " outstanding request."
+            << abort(FatalError);
+    }
+}
 
 
 // * * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * //
 
 template<class Type>
-processorFvPatchField<Type>::~processorFvPatchField()
+Foam::processorFvPatchField<Type>::~processorFvPatchField()
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type>
-tmp<Field<Type> > processorFvPatchField<Type>::patchNeighbourField() const
+Foam::tmp<Foam::Field<Type> >
+Foam::processorFvPatchField<Type>::patchNeighbourField() const
 {
-    // Warning: returning own patch field, which on update stores
+    if (debug && !this->ready())
+    {
+        FatalErrorIn
+        (
+            "tmp<Field<Type> >"
+            "processorFvPatchField<Type>::patchNeighbourField() const"
+        )   << "On patch " << procPatch_.name()
+            << " outstanding request."
+            << abort(FatalError);
+    }
+
+    // Warning: returning own patch field, which only after update stores
     // actual neighbour data
     // HJ, 14/May/2009
     return *this;
@@ -173,27 +228,79 @@ tmp<Field<Type> > processorFvPatchField<Type>::patchNeighbourField() const
 
 
 template<class Type>
-void processorFvPatchField<Type>::initEvaluate
+void Foam::processorFvPatchField<Type>::initEvaluate
 (
     const Pstream::commsTypes commsType
 )
 {
     if (Pstream::parRun())
     {
-        procPatch_.compressedSend(commsType, this->patchInternalField()());
+        // Collect data into send buffer
+        sendBuf_ = this->patchInternalField();
+
+        if (commsType == Pstream::nonBlocking)
+        {
+            // Fast path. Receive into *this
+
+            outstandingRecvRequest_ = Pstream::nRequests();
+
+            IPstream::read
+            (
+                Pstream::nonBlocking,
+                procPatch_.neighbProcNo(),
+                reinterpret_cast<char*>(this->begin()),
+                this->byteSize(),
+                procPatch_.tag(),
+                procPatch_.comm()
+            );
+
+            outstandingSendRequest_ = Pstream::nRequests();
+
+            OPstream::write
+            (
+                Pstream::nonBlocking,
+                procPatch_.neighbProcNo(),
+                reinterpret_cast<const char*>(sendBuf_.begin()),
+                this->byteSize(),
+                procPatch_.tag(),
+                procPatch_.comm()
+            );
+        }
+        else
+        {
+            procPatch_.send(commsType, sendBuf_);
+        }
     }
 }
 
 
 template<class Type>
-void processorFvPatchField<Type>::evaluate
+void Foam::processorFvPatchField<Type>::evaluate
 (
     const Pstream::commsTypes commsType
 )
 {
     if (Pstream::parRun())
     {
-        procPatch_.compressedReceive<Type>(commsType, *this);
+        if (commsType == Pstream::nonBlocking)
+        {
+            // Fast path. Received into *this
+
+            if
+            (
+                outstandingRecvRequest_ >= 0
+             && outstandingRecvRequest_ < Pstream::nRequests()
+            )
+            {
+                Pstream::waitRequest(outstandingRecvRequest_);
+            }
+            outstandingSendRequest_ = -1;
+            outstandingRecvRequest_ = -1;
+        }
+        else
+        {
+            procPatch_.receive<Type>(commsType, *this);
+        }
 
         if (doTransform())
         {
@@ -204,14 +311,15 @@ void processorFvPatchField<Type>::evaluate
 
 
 template<class Type>
-tmp<Field<Type> > processorFvPatchField<Type>::snGrad() const
+Foam::tmp<Foam::Field<Type> >
+Foam::processorFvPatchField<Type>::snGrad() const
 {
     return this->patch().deltaCoeffs()*(*this - this->patchInternalField());
 }
 
 
 template<class Type>
-void processorFvPatchField<Type>::initInterfaceMatrixUpdate
+void Foam::processorFvPatchField<Type>::initInterfaceMatrixUpdate
 (
     const scalarField& psiInternal,
     scalarField&,
@@ -222,16 +330,68 @@ void processorFvPatchField<Type>::initInterfaceMatrixUpdate
     const bool switchToLhs
 ) const
 {
-    procPatch_.compressedSend
+    scalarSendBuf_ = this->patch().patchInternalField(psiInternal);
+
+    if (commsType == Pstream::nonBlocking)
+    {
+        // Fast path.
+        if (debug && !this->ready())
+        {
+            FatalErrorIn
             (
-        commsType,
-        this->patch().patchInternalField(psiInternal)()
+                "void processorFvPatchField<Type>::initInterfaceMatrixUpdate\n"
+                "(\n"
+                "    const scalarField& psiInternal,\n"
+                "    scalarField&,\n"
+                "    const lduMatrix&,\n"
+                "    const scalarField&,\n"
+                "    const direction,\n"
+                "    const Pstream::commsTypes commsType,\n"
+                "    const bool switchToLhs\n"
+                ") const"
+            )   << "On patch " << procPatch_.name()
+                << " outstanding request."
+                << abort(FatalError);
+        }
+
+        // Check buffer size
+        scalarReceiveBuf_.setSize(this->size());
+
+        outstandingRecvRequest_ = Pstream::nRequests();
+
+        IPstream::read
+        (
+            Pstream::nonBlocking,
+            procPatch_.neighbProcNo(),
+            reinterpret_cast<char*>(scalarReceiveBuf_.begin()),
+            scalarReceiveBuf_.byteSize(),
+            procPatch_.tag(),
+            procPatch_.comm()
         );
+
+        outstandingSendRequest_ = Pstream::nRequests();
+
+        OPstream::write
+        (
+            Pstream::nonBlocking,
+            procPatch_.neighbProcNo(),
+            reinterpret_cast<const char*>(scalarSendBuf_.begin()),
+            scalarSendBuf_.byteSize(),
+            procPatch_.tag(),
+            procPatch_.comm()
+        );
+    }
+    else
+    {
+        procPatch_.send(commsType, scalarSendBuf_);
+    }
+
+    const_cast<processorFvPatchField<Type>&>(*this).updatedMatrix() = false;
 }
 
 
 template<class Type>
-void processorFvPatchField<Type>::updateInterfaceMatrix
+void Foam::processorFvPatchField<Type>::updateInterfaceMatrix
 (
     const scalarField&,
     scalarField& result,
@@ -242,36 +402,65 @@ void processorFvPatchField<Type>::updateInterfaceMatrix
     const bool switchToLhs
 ) const
 {
-    scalarField pnf
+    if (this->updatedMatrix())
+    {
+        return;
+    }
+
+    if (commsType == Pstream::nonBlocking)
+    {
+        // Fast path.
+        if
         (
-        procPatch_.compressedReceive<scalar>(commsType, this->size())()
-    );
+            outstandingRecvRequest_ >= 0
+         && outstandingRecvRequest_ < Pstream::nRequests()
+        )
+        {
+            Pstream::waitRequest(outstandingRecvRequest_);
+        }
+
+        // Recv finished so assume sending finished as well.
+        outstandingSendRequest_ = -1;
+        outstandingRecvRequest_ = -1;
+    }
+    else
+    {
+        // Check size
+        scalarReceiveBuf_.setSize(this->size());
+
+        procPatch_.receive<scalar>(commsType, scalarReceiveBuf_);
+    }
+
+    // The data is now in scalarReceiveBuf_ for both cases
 
     // Transform according to the transformation tensor
-    transformCoupleField(pnf, cmpt);
+    transformCoupleField(scalarReceiveBuf_, cmpt);
 
     // Multiply the field by coefficients and add into the result
+
     const unallocLabelList& faceCells = this->patch().faceCells();
 
     if (switchToLhs)
     {
         forAll(faceCells, elemI)
         {
-            result[faceCells[elemI]] += coeffs[elemI]*pnf[elemI];
+            result[faceCells[elemI]] += coeffs[elemI]*scalarReceiveBuf_[elemI];
         }
     }
     else
     {
         forAll(faceCells, elemI)
         {
-            result[faceCells[elemI]] -= coeffs[elemI]*pnf[elemI];
+            result[faceCells[elemI]] -= coeffs[elemI]*scalarReceiveBuf_[elemI];
         }
     }
+
+    const_cast<processorFvPatchField<Type>&>(*this).updatedMatrix() = true;
 }
 
 
 template<class Type>
-void processorFvPatchField<Type>::initInterfaceMatrixUpdate
+void Foam::processorFvPatchField<Type>::initInterfaceMatrixUpdate
 (
     const Field<Type>& psiInternal,
     Field<Type>&,
@@ -281,16 +470,71 @@ void processorFvPatchField<Type>::initInterfaceMatrixUpdate
     const bool switchToLhs
 ) const
 {
-    procPatch_.compressedSend
+    sendBuf_ = this->patch().patchInternalField(psiInternal);
+
+    if (commsType == Pstream::nonBlocking)
+    {
+        // Fast path.
+        if (debug && !this->ready())
+        {
+            FatalErrorIn
+            (
+                "void processorFvPatchField<Type>::initInterfaceMatrixUpdate\n"
+                "(\n"
+                "    const Field<Type>& psiInternal,\n"
+                "    Field<Type>&,\n"
+                "    const BlockLduMatrix<Type>&,\n"
+                "    const CoeffField<Type>&,\n"
+                "    const Pstream::commsTypes commsType,\n"
+                "    const bool switchToLhs\n"
+                ") const"
+            )   << "On patch " << procPatch_.name()
+                << " outstanding request."
+                << abort(FatalError);
+        }
+
+        // Check buffer size
+        receiveBuf_.setSize(this->size());
+
+        outstandingRecvRequest_ = Pstream::nRequests();
+
+        IPstream::read
+        (
+            Pstream::nonBlocking,
+            procPatch_.neighbProcNo(),
+            reinterpret_cast<char*>(receiveBuf_.begin()),
+            receiveBuf_.byteSize(),
+            procPatch_.tag(),
+            procPatch_.comm()
+        );
+
+        outstandingSendRequest_ = Pstream::nRequests();
+
+        OPstream::write
+        (
+            Pstream::nonBlocking,
+            procPatch_.neighbProcNo(),
+            reinterpret_cast<const char*>(sendBuf_.begin()),
+            sendBuf_.byteSize(),
+            procPatch_.tag(),
+            procPatch_.comm()
+        );
+    }
+    else
+    {
+        procPatch_.send
         (
             commsType,
             this->patch().patchInternalField(psiInternal)()
         );
+    }
+
+    const_cast<processorFvPatchField<Type>&>(*this).updatedMatrix() = false;
 }
 
 
 template<class Type>
-void processorFvPatchField<Type>::updateInterfaceMatrix
+void Foam::processorFvPatchField<Type>::updateInterfaceMatrix
 (
     const Field<Type>& psiInternal,
     Field<Type>& result,
@@ -300,14 +544,40 @@ void processorFvPatchField<Type>::updateInterfaceMatrix
     const bool switchToLhs
 ) const
 {
-    Field<Type> pnf
-    (
-        procPatch_.compressedReceive<Type>(commsType, this->size())()
-    );
+    if (this->updatedMatrix())
+    {
+        return;
+    }
 
-    // Multiply neighbour field with coeffs and re-use pnf for result
+    if (commsType == Pstream::nonBlocking)
+    {
+        // Fast path.
+        if
+        (
+            outstandingRecvRequest_ >= 0
+         && outstandingRecvRequest_ < Pstream::nRequests()
+        )
+        {
+            Pstream::waitRequest(outstandingRecvRequest_);
+        }
+
+        // Recv finished so assume sending finished as well.
+        outstandingSendRequest_ = -1;
+        outstandingRecvRequest_ = -1;
+    }
+    else
+    {
+        // Check size
+        receiveBuf_.setSize(this->size());
+
+        procPatch_.receive<Type>(commsType, receiveBuf_);
+    }
+
+    // The data is now in receiveBuf_ for both cases
+
+    // Multiply neighbour field with coeffs and re-use buffer for result
     // of multiplication
-    multiply(pnf, coeffs, pnf);
+    multiply(receiveBuf_, coeffs, receiveBuf_);
 
     const unallocLabelList& faceCells = this->patch().faceCells();
 
@@ -315,21 +585,56 @@ void processorFvPatchField<Type>::updateInterfaceMatrix
     {
         forAll(faceCells, elemI)
         {
-            result[faceCells[elemI]] += pnf[elemI];
+            result[faceCells[elemI]] += receiveBuf_[elemI];
         }
     }
     else
     {
         forAll(faceCells, elemI)
         {
-            result[faceCells[elemI]] -= pnf[elemI];
+            result[faceCells[elemI]] -= receiveBuf_[elemI];
         }
     }
+
+    const_cast<processorFvPatchField<Type>&>(*this).updatedMatrix() = true;
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+template<class Type>
+bool Foam::processorFvPatchField<Type>::ready() const
+{
+    if
+    (
+        outstandingSendRequest_ >= 0
+     && outstandingSendRequest_ < Pstream::nRequests()
+    )
+    {
+        bool finished = Pstream::finishedRequest(outstandingSendRequest_);
+
+        if (!finished)
+        {
+            return false;
+        }
+    }
+    outstandingSendRequest_ = -1;
+
+    if
+    (
+        outstandingRecvRequest_ >= 0
+     && outstandingRecvRequest_ < Pstream::nRequests()
+    )
+    {
+        bool finished = Pstream::finishedRequest(outstandingRecvRequest_);
+
+        if (!finished)
+        {
+            return false;
+        }
+    }
+    outstandingRecvRequest_ = -1;
+
+    return true;
+}
 
-} // End namespace Foam
 
 // ************************************************************************* //

src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchField.H

@@ -69,6 +69,25 @@ class processorFvPatchField
         //- Local reference cast into the processor patch
         const processorFvPatch& procPatch_;
 
+        // Sending and receiving
+
+            //- Outstanding request
+            mutable label outstandingSendRequest_;
+
+            //- Outstanding request
+            mutable label outstandingRecvRequest_;
+
+            //- Send buffer.
+            mutable Field<Type> sendBuf_;
+
+            //- Receive buffer.
+            mutable Field<Type> receiveBuf_;
+
+            //- Scalar send buffer
+            mutable Field<scalar> scalarSendBuf_;
+
+            //- Scalar receive buffer
+            mutable Field<scalar> scalarReceiveBuf_;
 
 public:
 
@@ -143,7 +162,7 @@ public:
 
 
     //- Destructor
-    ~processorFvPatchField();
+    virtual ~processorFvPatchField();
 
 
     // Member functions
@@ -249,7 +268,7 @@ public:
             ) const;
 
 
-        //- Processor coupled interface functions
+        // Processor coupled interface functions
 
             //- Return processor number
             virtual int myProcNo() const
@@ -263,6 +282,19 @@ public:
                 return procPatch_.neighbProcNo();
             }
 
+
+            // Communication support
+
+                //- Is all data available
+                virtual bool ready() const;
+
+                //- Return communicator used for parallel communication
+                virtual int comm() const
+                {
+                    return procPatch_.comm();
+                }
+
+
             //- Does the patch field perform the transfromation
             virtual bool doTransform() const
             {

src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchScalarField.C

@@ -32,26 +32,6 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<>
-void processorFvPatchField<scalar>::initInterfaceMatrixUpdate
-(
-    const scalarField& psiInternal,
-    scalarField&,
-    const lduMatrix&,
-    const scalarField&,
-    const direction,
-    const Pstream::commsTypes commsType,
-    const bool switchToLhs
-) const
-{
-    procPatch_.compressedSend
-    (
-        commsType,
-        patch().patchInternalField(psiInternal)()
-    );
-}
-
-
 template<>
 void processorFvPatchField<scalar>::updateInterfaceMatrix
 (
@@ -64,82 +44,60 @@ void processorFvPatchField<scalar>::updateInterfaceMatrix
     const bool switchToLhs
 ) const
 {
-    scalarField pnf
-    (
-        procPatch_.compressedReceive<scalar>(commsType, this->size())()
-    );
+    if (this->updatedMatrix())
+    {
+        return;
+    }
 
-    const unallocLabelList& faceCells = patch().faceCells();
+    if (commsType == Pstream::nonBlocking)
+    {
+        // Fast path.
+        if
+        (
+            outstandingRecvRequest_ >= 0
+         && outstandingRecvRequest_ < Pstream::nRequests()
+        )
+        {
+            Pstream::waitRequest(outstandingRecvRequest_);
+        }
+
+        // Recv finished so assume sending finished as well.
+        outstandingSendRequest_ = -1;
+        outstandingRecvRequest_ = -1;
+    }
+    else
+    {
+        // Check size
+        scalarReceiveBuf_.setSize(this->size());
+
+        procPatch_.receive<scalar>(commsType, scalarReceiveBuf_);
+    }
+
+    // The data is now in scalarReceiveBuf_ for both cases
+
+    // Transform according to the transformation tensor
+    // No transform for scalar.  HJ, 29/Nov/2016
+    
+    // Multiply the field by coefficients and add into the result
+
+    const unallocLabelList& faceCells = this->patch().faceCells();
 
     if (switchToLhs)
     {
-        forAll (faceCells, facei)
+        forAll(faceCells, elemI)
         {
-            result[faceCells[facei]] += coeffs[facei]*pnf[facei];
+            result[faceCells[elemI]] += coeffs[elemI]*scalarReceiveBuf_[elemI];
         }
     }
     else
     {
-        forAll (faceCells, facei)
+        forAll(faceCells, elemI)
         {
-            result[faceCells[facei]] -= coeffs[facei]*pnf[facei];
+            result[faceCells[elemI]] -= coeffs[elemI]*scalarReceiveBuf_[elemI];
         }
     }
-}
 
-
-template<>
-void processorFvPatchField<scalar>::initInterfaceMatrixUpdate
-(
-    const scalarField& psiInternal,
-    scalarField&,
-    const BlockLduMatrix<scalar>&,
-    const CoeffField<scalar>&,
-    const Pstream::commsTypes commsType,
-    const bool switchToLhs
-) const
-{
-    procPatch_.compressedSend
-    (
-        commsType,
-        patch().patchInternalField(psiInternal)()
-    );
-}
-
-
-template<>
-void processorFvPatchField<scalar>::updateInterfaceMatrix
-(
-    const scalarField&,
-    scalarField& result,
-    const BlockLduMatrix<scalar>&,
-    const CoeffField<scalar>& coeffs,
-    const Pstream::commsTypes commsType,
-    const bool switchToLhs
-) const
-{
-    scalarField pnf
-    (
-        procPatch_.compressedReceive<scalar>(commsType, this->size())()
-    );
-
-    const unallocLabelList& faceCells = patch().faceCells();
-    const scalarField& scalarCoeffs = coeffs.asScalar();
-
-    if (switchToLhs)
-    {
-        forAll (faceCells, facei)
-        {
-            result[faceCells[facei]] += scalarCoeffs[facei]*pnf[facei];
-        }
-    }
-    else
-    {
-        forAll (faceCells, facei)
-        {
-            result[faceCells[facei]] -= scalarCoeffs[facei]*pnf[facei];
-        }
-    }
+    const_cast<processorFvPatchField<scalar>&>(*this).updatedMatrix() = true;
 }
 
 

src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchScalarField.H

@@ -36,18 +36,8 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<>
-void processorFvPatchField<scalar>::initInterfaceMatrixUpdate
-(
-    const scalarField&,
-    scalarField&,
-    const lduMatrix&,
-    const scalarField&,
-    const direction,
-    const Pstream::commsTypes commsType,
-    const bool switchToLhs
-) const;
-
+// Removed unnecessary specialisation for initInterfaceMatrixUpdate
+// HJ, 29/Nov/2016
 
 template<>
 void processorFvPatchField<scalar>::updateInterfaceMatrix
@@ -62,29 +52,11 @@ void processorFvPatchField<scalar>::updateInterfaceMatrix
 ) const;
 
 
-template<>
-void processorFvPatchField<scalar>::initInterfaceMatrixUpdate
-(
-    const scalarField& psiInternal,
-    scalarField&,
-    const BlockLduMatrix<scalar>&,
-    const CoeffField<scalar>&,
-    const Pstream::commsTypes commsType,
-    const bool switchToLhs
-) const;
+// Removed unnecessary specialisation for initInterfaceMatrixUpdate
+// and updateInterfaceMatrix for block matrices
+// HJ, 29/Nov/2016
 
 
-template<>
-void processorFvPatchField<scalar>::updateInterfaceMatrix
-(
-    const scalarField&,
-    scalarField& result,
-    const BlockLduMatrix<scalar>&,
-    const CoeffField<scalar>& coeffs,
-    const Pstream::commsTypes commsType,
-    const bool switchToLhs
-) const;
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/finiteVolume/fields/fvPatchFields/constraint/wedge/wedgeFvPatchField.C

@@ -141,7 +141,10 @@ tmp<Field<Type> > wedgeFvPatchField<Type>::snGrad() const
     Field<Type> pif = this->patchInternalField();
     return
     (
-        transform(refCast<const wedgeFvPatch>(this->patch()).cellT(), pif) - pif
+        transform
+        (
+            refCast<const wedgeFvPatch>(this->patch()).cellT(), pif
+        ) - pif
     )*(0.5*this->patch().deltaCoeffs());
 }
 

src/finiteVolume/fields/fvPatchFields/constraint/wedge/wedgeFvPatchField.H

@@ -44,7 +44,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                         Class wedgeFvPatch Declaration
+                     Class wedgeFvPatchField Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class Type>

src/finiteVolume/fields/fvPatchFields/derived/pressureDirectedInletOutletVelocity/pressureDirectedInletOutletVelocityFvPatchVectorField.C

@@ -46,7 +46,7 @@ pressureDirectedInletOutletVelocityFvPatchVectorField
     mixedFvPatchVectorField(p, iF),
     phiName_("phi"),
     rhoName_("rho"),
-    inletDir_(p.size())
+    inletDir_(p.size(), vector(1, 0, 0))
 {
     refValue() = *this;
     refGrad() = vector::zero;
@@ -84,6 +84,30 @@ pressureDirectedInletOutletVelocityFvPatchVectorField
     inletDir_("inletDirection", dict, p.size())
 {
     fvPatchVectorField::operator=(vectorField("value", dict, p.size()));
+
+    if (!inletDir_.empty())
+    {
+        if (min(mag(inletDir_)) < SMALL)
+        {
+            FatalErrorIn
+            (
+                "pressureDirectedInletOutletVelocityFvPatchVectorField::\n"
+                "pressureDirectedInletOutletVelocityFvPatchVectorField\n"
+                "(\n"
+                "    const fvPatch& p,\n"
+                "    const DimensionedField<vector, volMesh>& iF,\n"
+                "    const dictionary& dict\n"
+                ")"
+            )    << "Badly defined inlet direction for field "
+                 << this->dimensionedInternalField().name()
+                 << " and patch " << this->patch().name()
+                 << abort(FatalError);
+        }
+    }       
+
+    // Normalise to obtain the flow direction
+    inletDir_ /= (mag(inletDir_) + SMALL);
+
     refValue() = *this;
     refGrad() = vector::zero;
     valueFraction() = 0.0;

src/finiteVolume/fields/fvPatchFields/derived/pressureDirectedInletVelocity/pressureDirectedInletVelocityFvPatchVectorField.C

@@ -46,7 +46,7 @@ pressureDirectedInletVelocityFvPatchVectorField
     fixedValueFvPatchVectorField(p, iF),
     phiName_("phi"),
     rhoName_("rho"),
-    inletDir_(p.size())
+    inletDir_(p.size(), vector(1, 0, 0))
 {}
 
 
@@ -80,6 +80,29 @@ pressureDirectedInletVelocityFvPatchVectorField
     inletDir_("inletDirection", dict, p.size())
 {
     fvPatchVectorField::operator=(vectorField("value", dict, p.size()));
+
+    if (!inletDir_.empty())
+    {
+        if (min(mag(inletDir_)) < SMALL)
+        {
+            FatalErrorIn
+            (
+                "pressureDirectedInletVelocityFvPatchVectorField::\n"
+                "pressureDirectedInletVelocityFvPatchVectorField\n"
+                "(\n"
+                "    const fvPatch& p,\n"
+                "    const DimensionedField<vector, volMesh>& iF,\n"
+                "    const dictionary& dict\n"
+                ")"
+            )    << "Badly defined inlet direction for field "
+                 << this->dimensionedInternalField().name()
+                 << " and patch " << this->patch().name()
+                 << abort(FatalError);
+        }
+    }       
+
+    // Normalise to obtain the flow direction
+    inletDir_ /= (mag(inletDir_) + SMALL);
 }
 
 

src/finiteVolume/fields/fvPatchFields/derived/pressureInletOutletVelocity/pressureInletOutletVelocityFvPatchVectorField.C

@@ -192,7 +192,6 @@ void pressureInletOutletVelocityFvPatchVectorField::updateCoeffs()
     }
 
     directionMixedFvPatchVectorField::updateCoeffs();
-    directionMixedFvPatchVectorField::evaluate();
 }
 
 

src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C

@@ -698,7 +698,10 @@ void Foam::fvMatrix<Type>::relax()
 {
     if (psi_.mesh().solutionDict().relaxEquation(psi_.name()))
     {
-        relax(psi_.mesh().solutionDict().equationRelaxationFactor(psi_.name()));
+        relax
+        (
+            psi_.mesh().solutionDict().equationRelaxationFactor(psi_.name())
+        );
     }
     else
     {

src/finiteVolume/fvMesh/fvMesh.C

@@ -556,6 +556,12 @@ void Foam::fvMesh::updateMesh(const mapPolyMesh& mpm)
     // This is a temporary solution
     surfaceInterpolation::movePoints();
 
+    // Note: deltaCoeffs cannot be left on lazy evaluation on mesh motion
+    // because tangled comms will occur when they are accessed from
+    // individual boundary conditions
+    // HJ, VV and IG, 25/Oct/2016
+    deltaCoeffs();
+
     // Function object update moved to polyMesh
     // HJ, 29/Aug/2010
 }
@@ -577,6 +583,12 @@ void Foam::fvMesh::syncUpdateMesh()
     // This is a temporary solution
     surfaceInterpolation::movePoints();
 
+    // Note: deltaCoeffs cannot be left on lazy evaluation on mesh motion
+    // because tangled comms will occur when they are accessed from
+    // individual boundary conditions
+    // HJ, VV and IG, 25/Oct/2016
+    deltaCoeffs();
+
     // Function object update moved to polyMesh
     // HJ, 29/Aug/2010
 }
@@ -681,6 +693,12 @@ Foam::tmp<Foam::scalarField> Foam::fvMesh::movePoints(const pointField& p)
     // Function object update moved to polyMesh
     // HJ, 29/Aug/2010
 
+    // Note: deltaCoeffs cannot be left on lazy evaluation on mesh motion
+    // because tangled comms will occur when they are accessed from
+    // individual boundary conditions
+    // HJ, VV and IG, 25/Oct/2016
+    deltaCoeffs();
+
     return tsweptVols;
 }
 

src/finiteVolume/fvMesh/fvMesh.H

@@ -123,6 +123,7 @@ class fvMesh
 
        // Make geometric data
 
+            //- Make face areas
             void makeSf() const;
 
             //- Make face area magnitudes
@@ -216,8 +217,7 @@ public:
         );
 
 
-    // Destructor
-
+    //- Destructor
     virtual ~fvMesh();
 
 
@@ -285,6 +285,15 @@ public:
             }
 
 
+            // Communication support
+
+                //- Return communicator used for parallel communication
+                label comm() const
+                {
+                    return polyMesh::comm();
+                }
+
+
             //- Return cell volumes
             const DimensionedField<scalar, volMesh>& V() const;
 

src/finiteVolume/fvMesh/fvPatches/constraint/cyclicGgi/cyclicGgiFvPatch.H

@@ -90,8 +90,7 @@ public:
         {}
 
 
-    // Destructor
-
+    //- Destructor
     virtual ~cyclicGgiFvPatch()
     {}
 

src/finiteVolume/fvMesh/fvPatches/constraint/ggi/ggiFvPatch.C

@@ -258,6 +258,12 @@ bool Foam::ggiFvPatch::localParallel() const
 }
 
 
+const Foam::mapDistribute& Foam::ggiFvPatch::map() const
+{
+    return ggiPolyPatch_.map();
+}
+
+
 const Foam::scalarListList& Foam::ggiFvPatch::weights() const
 {
     if (ggiPolyPatch_.master())
@@ -271,6 +277,12 @@ const Foam::scalarListList& Foam::ggiFvPatch::weights() const
 }
 
 
+void Foam::ggiFvPatch::expandAddrToZone(labelField& lf) const
+{
+    lf = ggiPolyPatch_.fastExpand(lf);
+}
+
+
 Foam::tmp<Foam::labelField> Foam::ggiFvPatch::interfaceInternalField
 (
     const unallocLabelList& internalData
@@ -306,6 +318,7 @@ void Foam::ggiFvPatch::initInternalFieldTransfer
     const unallocLabelList& iF
 ) const
 {
+    // Label transfer is local without global reduction
     labelTransferBuffer_ = patchInternalField(iF);
 }
 

src/finiteVolume/fvMesh/fvPatches/constraint/ggi/ggiFvPatch.H

@@ -92,13 +92,27 @@ public:
         {}
 
 
-    // Destructor
-
+    //- Destructor
     virtual ~ggiFvPatch();
 
 
     // Member functions
 
+        // Communication support
+
+            //- Return communicator used for parallel communication
+            virtual int comm() const
+            {
+                return ggiPolyPatch_.comm();
+            }
+
+            //- Return message tag used for sending
+            virtual int tag() const
+            {
+                return ggiPolyPatch_.tag();
+            }
+
+
         // Access
 
             //- Return shadow patch
@@ -171,6 +185,9 @@ public:
             //- Is the patch localised on a single processor
             virtual bool localParallel() const;
 
+            //- Return mapDistribute
+            virtual const mapDistribute& map() const;
+
             //- Return weights.  Master side returns own weights and
             //  slave side returns weights from master
             virtual const scalarListList& weights() const;
@@ -187,6 +204,11 @@ public:
                 return coupledFvPatch::reverseT();
             }
 
+            //- Expand addressing to zone
+            //  Used in optimised AMG coarsening
+            virtual void expandAddrToZone(labelField&) const;
+
+
             //- Return the values of the given internal data adjacent to
             //  the interface as a field
             virtual tmp<labelField> interfaceInternalField

src/finiteVolume/fvMesh/fvPatches/constraint/mixingPlane/mixingPlaneFvPatch.H

@@ -95,13 +95,27 @@ public:
         {}
 
 
-    // Destructor
-
+    //- Destructor
     virtual ~mixingPlaneFvPatch();
 
 
     // Member functions
 
+        // Communication support
+
+            //- Return communicator used for parallel communication
+            virtual int comm() const
+            {
+                return mixingPlanePolyPatch_.comm();
+            }
+
+            //- Return message tag used for sending
+            virtual int tag() const
+            {
+                return mixingPlanePolyPatch_.tag();
+            }
+
+
         // Access
 
             //- Return shadow patch

src/finiteVolume/fvMesh/fvPatches/constraint/processor/processorFvPatch.H

@@ -86,8 +86,7 @@ public:
         {}
 
 
-    // Destructor
-
+    //- Destructor
     virtual ~processorFvPatch()
     {}
 
@@ -119,6 +118,21 @@ public:
             }
         }
 
+        // Communication support
+
+            //- Return communicator used for parallel communication
+            virtual int comm() const
+            {
+                return procPolyPatch_.comm();
+            }
+
+            //- Return message tag used for sending
+            virtual int tag() const
+            {
+                return procPolyPatch_.tag();
+            }
+
+
         //- Return face transformation tensor
         virtual const tensorField& forwardT() const
         {

src/finiteVolume/fvMesh/fvPatches/constraint/regionCouple/regionCoupleFvPatch.C

@@ -34,6 +34,7 @@ Author
 #include "fvMesh.H"
 #include "fvBoundaryMesh.H"
 #include "foamTime.H"
+#include "mapDistribute.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -302,6 +303,12 @@ bool Foam::regionCoupleFvPatch::localParallel() const
 }
 
 
+const Foam::mapDistribute& Foam::regionCoupleFvPatch::map() const
+{
+    return rcPolyPatch_.map();
+}
+
+
 const Foam::scalarListList& Foam::regionCoupleFvPatch::weights() const
 {
     if (rcPolyPatch_.master())
@@ -315,6 +322,16 @@ const Foam::scalarListList& Foam::regionCoupleFvPatch::weights() const
 }
 
 
+void Foam::regionCoupleFvPatch::expandAddrToZone(labelField& lf) const
+{
+    // Missing code.  Activate for AMG solvers across regionCoupleFvPatch
+    notImplemented
+    (
+        "void regionCoupleFvPatch::expandAddrToZone(labelField& lf) const"
+    );
+}
+
+
 Foam::tmp<Foam::labelField> Foam::regionCoupleFvPatch::interfaceInternalField
 (
     const unallocLabelList& internalData

src/finiteVolume/fvMesh/fvPatches/constraint/regionCouple/regionCoupleFvPatch.H

@@ -95,8 +95,7 @@ public:
         {}
 
 
-    // Destructor
-
+    //- Destructor
     virtual ~regionCoupleFvPatch();
 
 
@@ -123,6 +122,21 @@ public:
             virtual tmp<vectorField> delta() const;
 
 
+        // Communication support
+
+            //- Return communicator used for parallel communication
+            virtual int comm() const
+            {
+                return rcPolyPatch_.comm();
+            }
+
+            //- Return message tag used for sending
+            virtual int tag() const
+            {
+                return rcPolyPatch_.tag();
+            }
+
+
         // Interpolation
 
             //- Interpolate face field
@@ -180,6 +194,9 @@ public:
             //- Is the patch localised on a single processor
             virtual bool localParallel() const;
 
+            //- Return mapDistribute
+            const mapDistribute& map() const;
+
             //- Return weights.  Master side returns own weights and
             //  slave side returns weights from master
             virtual const scalarListList& weights() const;
@@ -196,6 +213,11 @@ public:
                 return coupledFvPatch::reverseT();
             }
 
+            //- Expand addressing to zone
+            //  Used in optimised AMG coarsening
+            virtual void expandAddrToZone(labelField&) const;
+
+
             //- Return the values of the given internal data adjacent to
             //  the interface as a field
             virtual tmp<labelField> interfaceInternalField

src/foam/Make/files

@@ -141,13 +141,13 @@ $(StringStreams)/StringStreamsPrint.C
 
 Pstreams = $(Streams)/Pstreams
 $(Pstreams)/Pstream.C
+$(Pstreams)/PstreamReduceOps.C
 $(Pstreams)/PstreamCommsStruct.C
 $(Pstreams)/PstreamGlobals.C
 $(Pstreams)/IPstream.C
 $(Pstreams)/OPstream.C
 $(Pstreams)/IPread.C
 $(Pstreams)/OPwrite.C
-$(Pstreams)/PstreamsPrint.C
 
 dictionary = db/dictionary
 $(dictionary)/dictionary.C

src/foam/containers/Lists/List/List.H

@@ -99,6 +99,7 @@ public:
         //- Return a null list
         inline static const List<T>& null();
 
+
     // Constructors
 
         //- Null constructor.

src/foam/db/IOstreams/Pstreams/IPread.C

@@ -38,6 +38,8 @@ Foam::IPstream::IPstream
     const commsTypes commsType,
     const int fromProcNo,
     const label bufSize,
+    const int tag,
+    const label comm,
     streamFormat format,
     versionNumber version
 )
@@ -45,6 +47,8 @@ Foam::IPstream::IPstream
     Pstream(commsType, bufSize),
     Istream(format, version),
     fromProcNo_(fromProcNo),
+    tag_(tag),
+    comm_(comm),
     messageSize_(0)
 {
     setOpened();
@@ -52,17 +56,31 @@ Foam::IPstream::IPstream
 
     MPI_Status status;
 
-    // If the buffer size is not specified, probe the incomming message
+    // If the buffer size is not specified, probe the incoming message
     // and set it
     if (!bufSize)
     {
-        MPI_Probe(procID(fromProcNo_), msgType(), MPI_COMM_WORLD, &status);
+        MPI_Probe
+        (
+            fromProcNo_,
+            tag_,
+            PstreamGlobals::MPICommunicators_[comm_],
+            &status
+        );
         MPI_Get_count(&status, MPI_BYTE, &messageSize_);
 
         buf_.setSize(messageSize_);
     }
 
-    messageSize_ = read(commsType, fromProcNo_, buf_.begin(), buf_.size());
+    messageSize_ = IPstream::read
+    (
+        commsType,
+        fromProcNo_,
+        buf_.begin(),
+        buf_.size(),
+        tag_,
+        comm_
+    );
 
     if (!messageSize_)
     {
@@ -83,9 +101,31 @@ Foam::label Foam::IPstream::read
     const commsTypes commsType,
     const int fromProcNo,
     char* buf,
-    const std::streamsize bufSize
+    const std::streamsize bufSize,
+    const int tag,
+    const label comm
 )
 {
+    if (debug)
+    {
+        Pout<< "IPstream::read : starting read from:" << fromProcNo
+            << " tag:" << tag << " comm:" << comm
+            << " wanted size:" << label(bufSize)
+            << " commsType:" << Pstream::commsTypeNames[commsType]
+            << Foam::endl;
+    }
+
+    if (Pstream::warnComm != -1 && comm != Pstream::warnComm)
+    {
+        Pout<< "IPstream::read : starting read from:" << fromProcNo
+            << " tag:" << tag << " comm:" << comm
+            << " wanted size:" << label(bufSize)
+            << " commsType:" << Pstream::commsTypeNames[commsType]
+            << " warnComm:" << Pstream::warnComm
+            << Foam::endl;
+        error::printStack(Pout);
+    }
+
     if (commsType == blocking || commsType == scheduled)
     {
         MPI_Status status;
@@ -96,10 +136,10 @@ Foam::label Foam::IPstream::read
             (
                 buf,
                 bufSize,
-                MPI_PACKED,
-                procID(fromProcNo),
-                msgType(),
-                MPI_COMM_WORLD,
+                MPI_BYTE,
+                fromProcNo,
+                tag,
+                PstreamGlobals::MPICommunicators_[comm],
                 &status
             )
         )
@@ -144,10 +184,10 @@ Foam::label Foam::IPstream::read
             (
                 buf,
                 bufSize,
-                MPI_PACKED,
-                procID(fromProcNo),
-                msgType(),
-                MPI_COMM_WORLD,
+                MPI_BYTE,
+                fromProcNo,
+                tag,
+                PstreamGlobals::MPICommunicators_[comm],
                 &request
             )
         )
@@ -162,8 +202,18 @@ Foam::label Foam::IPstream::read
             return 0;
         }
 
-        PstreamGlobals::IPstream_outstandingRequests_.append(request);
+        if (debug)
+        {
+            Pout<< "IPstream::read : started read from:" << fromProcNo
+                << " tag:" << tag << " read size:" << label(bufSize)
+                << " commsType:" << Pstream::commsTypeNames[commsType]
+                << " request:" << PstreamGlobals::outstandingRequests_.size()
+                << Foam::endl;
+        }
 
+        PstreamGlobals::outstandingRequests_.append(request);
+
+        // Assume the message is completely received.
         return 1;
     }
     else
@@ -180,56 +230,4 @@ Foam::label Foam::IPstream::read
 }
 
 
-void Foam::IPstream::waitRequests()
-{
-    if (PstreamGlobals::IPstream_outstandingRequests_.size())
-    {
-        if
-        (
-            MPI_Waitall
-            (
-                PstreamGlobals::IPstream_outstandingRequests_.size(),
-                PstreamGlobals::IPstream_outstandingRequests_.begin(),
-                MPI_STATUSES_IGNORE
-            )
-        )
-        {
-            FatalErrorIn
-            (
-                "IPstream::waitRequests()"
-            )   << "MPI_Waitall returned with error" << endl;
-        }
-
-        PstreamGlobals::IPstream_outstandingRequests_.clear();
-    }
-}
-
-
-bool Foam::IPstream::finishedRequest(const label i)
-{
-    if (i >= PstreamGlobals::IPstream_outstandingRequests_.size())
-    {
-        FatalErrorIn
-        (
-            "IPstream::finishedRequest(const label)"
-        )   << "There are "
-            << PstreamGlobals::IPstream_outstandingRequests_.size()
-            << " outstanding send requests and you are asking for i=" << i
-            << nl
-            << "Maybe you are mixing blocking/non-blocking comms?"
-            << Foam::abort(FatalError);
-    }
-
-    int flag;
-    MPI_Test
-    (
-        &PstreamGlobals::IPstream_outstandingRequests_[i],
-        &flag,
-        MPI_STATUS_IGNORE
-    );
-
-    return flag != 0;
-}
-
-
 // ************************************************************************* //

src/foam/db/IOstreams/Pstreams/IPreadToken.C

@@ -1,187 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | foam-extend: Open Source CFD
-   \\    /   O peration     | Version:     4.0
-    \\  /    A nd           | Web:         http://www.foam-extend.org
-     \\/     M anipulation  | For copyright notice see file Copyright
--------------------------------------------------------------------------------
-License
-    This file is part of foam-extend.
-
-    foam-extend 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.
-
-    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "IPstream.H"
-#include "token.H"
-#include <cctype>
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-Istream& IPstream::read(token& t)
-{
-    // Return the put back token if it exists
-    if (Istream::getBack(t))
-    {
-        return *this;
-    }
-
-    char c;
-
-    // return on error
-    if (!read(c))
-    {
-        t.setBad();
-        return *this;
-    }
-
-    // Set the line number of this token to the current stream line number
-    t.lineNumber() = lineNumber();
-
-    // Analyse input starting with this character.
-    switch (c)
-    {
-        // Punctuation
-        case token::END_STATEMENT :
-        case token::BEGIN_LIST :
-        case token::END_LIST :
-        case token::BEGIN_SQR :
-        case token::END_SQR :
-        case token::BEGIN_BLOCK :
-        case token::END_BLOCK :
-        case token::COLON :
-        case token::COMMA :
-        case token::ASSIGN :
-        case token::ADD :
-        case token::SUBTRACT :
-        case token::MULTIPLY :
-        case token::DIVIDE :
-        {
-            t = token::punctuationToken(c);
-            return *this;
-        }
-
-        // Word
-        case token::WORD :
-        {
-            word* wPtr = new word;
-            if (read(*wPtr))
-            {
-                if (token::compound::isCompound(*wPtr))
-                {
-                    t = token::compound::New(*wPtr, *this).ptr();
-                    delete wPtr;
-                }
-                else
-                {
-                    t = wPtr;
-                }
-            }
-            else
-            {
-                delete wPtr;
-                t.setBad();
-            }
-            return *this;
-        }
-
-        // String
-        case token::STRING :
-        {
-            string* sPtr = new string;
-            if (read(*sPtr))
-            {
-                t = sPtr;
-            }
-            else
-            {
-                delete sPtr;
-                t.setBad();
-            }
-            return *this;
-        }
-
-        // Label
-        case token::LABEL :
-        {
-            label l;
-            if (read(l))
-            {
-                t = l;
-            }
-            else
-            {
-                t.setBad();
-            }
-            return *this;
-        }
-
-        // floatScalar
-        case token::FLOAT_SCALAR :
-        {
-            floatScalar s;
-            if (read(s))
-            {
-                t = s;
-            }
-            else
-            {
-                t.setBad();
-            }
-            return *this;
-        }
-
-        // doubleScalar
-        case token::DOUBLE_SCALAR :
-        {
-            doubleScalar s;
-            if (read(s))
-            {
-                t = s;
-            }
-            else
-            {
-                t.setBad();
-            }
-            return *this;
-        }
-
-        // Character (returned as a single character word) or error
-        default:
-        {
-            if (isalpha(c))
-            {
-                t = word(c);
-                return *this;
-            }
-
-            setBad();
-            t.setBad();
-
-            return *this;
-        }
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// ************************************************************************* //

src/foam/db/IOstreams/Pstreams/IPstream.C

@@ -324,4 +324,13 @@ Foam::Istream& Foam::IPstream::rewind()
 }
 
 
+void Foam::IPstream::print(Ostream& os) const
+{
+    os  << "Reading from processor " << fromProcNo_
+        << " using communicator " << comm_
+        <<  " and tag " << tag_
+        << Foam::endl;
+}
+
+
 // ************************************************************************* //

src/foam/db/IOstreams/Pstreams/IPstream.H

@@ -58,6 +58,12 @@ class IPstream
         //- ID of sending processor
         int fromProcNo_;
 
+        //- Message tag
+        const int tag_;
+
+        //- Communicator
+        const label comm_;
+
         //- Message size
         label messageSize_;
 
@@ -86,12 +92,14 @@ public:
             const commsTypes commsType,
             const int fromProcNo,
             const label bufSize = 0,
+            const int tag = Pstream::msgType(),
+            const label comm = Pstream::worldComm,
             streamFormat format = BINARY,
             versionNumber version = currentVersion
         );
 
 
-    // Destructor
+    //- Destructor
     ~IPstream();
 
 
@@ -115,15 +123,11 @@ public:
                 const commsTypes commsType,
                 const int fromProcNo,
                 char* buf,
-                const std::streamsize bufSize
+                const std::streamsize bufSize,
+                const int tag = Pstream::msgType(),
+                const label communicator = 0
             );
 
-            //- Non-blocking receives: wait until all have finished.
-            static void waitRequests();
-
-            //- Non-blocking receives: has request i finished?
-            static bool finishedRequest(const label i);
-
             //- Return next token from stream
             Istream& read(token&);
 

src/foam/db/IOstreams/Pstreams/OPstream.C

@@ -92,13 +92,17 @@ Foam::OPstream::OPstream
     const commsTypes commsType,
     const int toProcNo,
     const label bufSize,
+    const int tag,
+    const label comm,
     streamFormat format,
     versionNumber version
 )
 :
     Pstream(commsType, bufSize),
     Ostream(format, version),
-    toProcNo_(toProcNo)
+    toProcNo_(toProcNo),
+    tag_(tag),
+    comm_(comm)
 {
     setOpened();
     setGood();
@@ -233,4 +237,12 @@ Foam::Ostream& Foam::OPstream::write(const char* data, std::streamsize count)
 }
 
 
+void Foam::OPstream::print(Ostream& os) const
+{
+    os  << "Writing from processor " << toProcNo_
+        << " to processor " << myProcNo() << " in communicator " << comm_
+        << " and tag " << tag_ << Foam::endl;
+}
+
+
 // ************************************************************************* //

src/foam/db/IOstreams/Pstreams/OPstream.H

@@ -58,6 +58,12 @@ class OPstream
          // ID of receiving processor
         int toProcNo_;
 
+        //- Message tag
+        const int tag_;
+
+        //- Communicator
+        const label comm_;
+
 
     // Private member functions
 
@@ -88,13 +94,14 @@ public:
             const commsTypes commsType,
             const int toProcNo,
             const label bufSize = 0,
+            const int tag = Pstream::msgType(),
+            const label comm = Pstream::worldComm,
             streamFormat format = BINARY,
             versionNumber version = currentVersion
         );
 
 
-    // Destructor
-
+    //- Destructor
     ~OPstream();
 
 
@@ -117,15 +124,11 @@ public:
                 const commsTypes commsType,
                 const int toProcNo,
                 const char* buf,
-                const std::streamsize bufSize
+                const std::streamsize bufSize,
+                const int tag = Pstream::msgType(),
+                const label comm = Pstream::worldComm
             );
 
-            //- Non-blocking writes: wait until all have finished.
-            static void waitRequests();
-
-            //- Non-blocking writes: has request i finished?
-            static bool finishedRequest(const label i);
-
             //- Write next token to stream
             Ostream& write(const token&);
 

src/foam/db/IOstreams/Pstreams/OPwrite.C

@@ -42,7 +42,9 @@ Foam::OPstream::~OPstream()
             commsType_,
             toProcNo_,
             buf_.begin(),
-            bufPosition_
+            bufPosition_,
+            tag_,
+            comm_
         )
     )
     {
@@ -60,9 +62,34 @@ bool Foam::OPstream::write
     const commsTypes commsType,
     const int toProcNo,
     const char* buf,
-    const std::streamsize bufSize
+    const std::streamsize bufSize,
+    const int tag,
+    const label comm
 )
 {
+    if (debug)
+    {
+        Pout<< "OPstream::write : starting write to:" << toProcNo
+            << " tag:" << tag
+            << " comm:" << comm << " size:" << label(bufSize)
+            << " commsType:" << Pstream::commsTypeNames[commsType]
+            << Foam::endl;
+    }
+    if (Pstream::warnComm != -1 && comm != Pstream::warnComm)
+    {
+        Pout<< "OPstream::write : starting write to:" << toProcNo
+            << " tag:" << tag
+            << " comm:" << comm << " size:" << label(bufSize)
+            << " commsType:" << Pstream::commsTypeNames[commsType]
+            << " warnComm:" << Pstream::warnComm
+            << Foam::endl;
+
+        error::printStack(Pout);
+    }
+
+
+    PstreamGlobals::checkCommunicator(comm, toProcNo);
+
     bool transferFailed = true;
 
     if (commsType == blocking)
@@ -71,11 +98,19 @@ bool Foam::OPstream::write
         (
             const_cast<char*>(buf),
             bufSize,
-            MPI_PACKED,
-            procID(toProcNo),
-            msgType(),
-            MPI_COMM_WORLD
+            MPI_BYTE,
+            toProcNo,   //procID(toProcNo),
+            tag,
+            PstreamGlobals::MPICommunicators_[comm] // MPI_COMM_WORLD
         );
+
+        if (debug)
+        {
+            Pout<< "OPstream::write : finished write to:" << toProcNo
+                << " tag:" << tag << " size:" << label(bufSize)
+                << " commsType:" << Pstream::commsTypeNames[commsType]
+                << Foam::endl;
+        }
     }
     else if (commsType == scheduled)
     {
@@ -83,11 +118,19 @@ bool Foam::OPstream::write
         (
             const_cast<char*>(buf),
             bufSize,
-            MPI_PACKED,
-            procID(toProcNo),
-            msgType(),
-            MPI_COMM_WORLD
+            MPI_BYTE,
+            toProcNo,   //procID(toProcNo),
+            tag,
+            PstreamGlobals::MPICommunicators_[comm] // MPI_COMM_WORLD
         );
+
+        if (debug)
+        {
+            Pout<< "OPstream::write : finished write to:" << toProcNo
+                << " tag:" << tag << " size:" << label(bufSize)
+                << " commsType:" << Pstream::commsTypeNames[commsType]
+                << Foam::endl;
+        }
     }
     else if (commsType == nonBlocking)
     {
@@ -97,14 +140,23 @@ bool Foam::OPstream::write
         (
             const_cast<char*>(buf),
             bufSize,
-            MPI_PACKED,
-            procID(toProcNo),
-            msgType(),
-            MPI_COMM_WORLD,
+            MPI_BYTE,
+            toProcNo,   //procID(toProcNo),
+            tag,
+            PstreamGlobals::MPICommunicators_[comm],// MPI_COMM_WORLD,
             &request
         );
 
-        PstreamGlobals::OPstream_outstandingRequests_.append(request);
+        if (debug)
+        {
+            Pout<< "OPstream::write : started write to:" << toProcNo
+                << " tag:" << tag << " size:" << label(bufSize)
+                << " commsType:" << Pstream::commsTypeNames[commsType]
+                << " request:" << PstreamGlobals::outstandingRequests_.size()
+                << Foam::endl;
+        }
+
+        PstreamGlobals::outstandingRequests_.append(request);
     }
     else
     {
@@ -120,56 +172,4 @@ bool Foam::OPstream::write
 }
 
 
-void Foam::OPstream::waitRequests()
-{
-    if (PstreamGlobals::OPstream_outstandingRequests_.size())
-    {
-        if
-        (
-            MPI_Waitall
-            (
-                PstreamGlobals::OPstream_outstandingRequests_.size(),
-                PstreamGlobals::OPstream_outstandingRequests_.begin(),
-                MPI_STATUSES_IGNORE
-            )
-        )
-        {
-            FatalErrorIn
-            (
-                "OPstream::waitRequests()"
-            )   << "MPI_Waitall returned with error" << Foam::endl;
-        }
-
-        PstreamGlobals::OPstream_outstandingRequests_.clear();
-    }
-}
-
-
-bool Foam::OPstream::finishedRequest(const label i)
-{
-    if (i >= PstreamGlobals::OPstream_outstandingRequests_.size())
-    {
-        FatalErrorIn
-        (
-            "OPstream::finishedRequest(const label)"
-        )   << "There are "
-            << PstreamGlobals::OPstream_outstandingRequests_.size()
-            << " outstanding send requests and you are asking for i=" << i
-            << nl
-            << "Maybe you are mixing blocking/non-blocking comms?"
-            << Foam::abort(FatalError);
-    }
-
-    int flag;
-    MPI_Test
-    (
-        &PstreamGlobals::OPstream_outstandingRequests_[i],
-        &flag,
-        MPI_STATUS_IGNORE
-    );
-
-    return flag != 0;
-}
-
-
 // ************************************************************************* //

src/foam/db/IOstreams/Pstreams/Pstream.H

@@ -29,11 +29,11 @@ Description
 
 SourceFiles
     Pstream.C
-    PstreamsPrint.C
     PstreamCommsStruct.C
     gatherScatter.C
     combineGatherScatter.C
     gatherScatterList.C
+    PstreamExchange.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -47,6 +47,8 @@ SourceFiles
 #include "NamedEnum.H"
 #include "dynamicLabelList.H"
 #include "optimisationSwitch.H"
+#include "ListOps.H"
+#include "LIFOStack.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -59,7 +61,6 @@ namespace Foam
 
 class Pstream
 {
-
 public:
 
     //- Types of communications
@@ -72,6 +73,7 @@ public:
 
     static const NamedEnum<commsTypes, 3> commsTypeNames;
 
+
     // Public classes
 
         //- Structure for communicating between processors
@@ -85,7 +87,7 @@ public:
                 //- procIDs of processors directly below me
                 labelList below_;
 
-                //- procIDs of all processors below (so not just directly below)
+                //- procIDs of all processors below (not just directly below)
                 labelList allBelow_;
 
                 //- procIDs of all processors not below. (inverse set of
@@ -165,7 +167,7 @@ public:
             template<class T>
             void operator()(T& x, const T& y) const
             {
-                forAll(y, i)
+                forAll (y, i)
                 {
                     if (y[i].size())
                     {
@@ -183,32 +185,38 @@ private:
         //- Is this a parallel run?
         static bool parRun_;
 
-        //- My processor index
-        static int myProcNo_;
-
-        //- Process IDs
-        static List<int> procIDs_;
-
-        //- Default message type
+        //- Default message type info
         static const int msgType_;
 
+        //- Stack of free comms
+        static LIFOStack<label> freeComms_;
+
+        //- My processor index
+        static DynamicList<int> myProcNo_;
+
+        //- Process IDs
+        static DynamicList<List<int> > procIDs_;
+
+        //- List of parent communicators
+        static DynamicList<label> parentCommunicator_;
+
         //- Structure for linear communications
-        static List<commsStruct> linearCommunication_;
+        static DynamicList<List<commsStruct> > linearCommunication_;
 
         //- Structure for tree communications
-        static List<commsStruct> treeCommunication_;
+        static DynamicList<List<commsStruct> > treeCommunication_;
 
 
     // Private member functions
 
         //- Set data for parallel running
-        static void setParRun();
+        static void setParRun(const label nProcs);
 
         //- Calculate linear communication schedule
-        static void calcLinearComm(const label nProcs);
+        static List<commsStruct> calcLinearComm(const label nProcs);
 
         //- Calculate tree communication schedule
-        static void calcTreeComm(const label nProcs);
+        static List<commsStruct> calcTreeComm(const label nProcs);
 
         //- Helper function for tree communication schedule determination
         //  Collects all processorIDs below a processor
@@ -223,6 +231,19 @@ private:
         //  Pstream::init()
         static void initCommunicationSchedule();
 
+        //- Allocate a communicator with index
+        static void allocatePstreamCommunicator
+        (
+            const label parentIndex,
+            const label index
+        );
+
+        //- Free a communicator
+        static void freePstreamCommunicator
+        (
+            const label index
+        );
+
 
 protected:
 
@@ -252,11 +273,6 @@ public:
 
     // Static data
 
-        //- Should compact transfer be used in which floats replace doubles
-        //  reducing the bandwidth requirement at the expense of some loss
-        //  in accuracy
-        static const debug::optimisationSwitch floatTransfer;
-
         //- Number of processors at which the sum algorithm changes from linear
         //  to tree
         static const debug::optimisationSwitch nProcsSimpleSum;
@@ -273,6 +289,15 @@ public:
             );
         }
 
+        //- Number of polling cycles in processor updates
+        static const debug::optimisationSwitch nPollProcInterfaces;
+
+        //- Default communicator (all processors)
+        static label worldComm;
+
+        //- Debugging: warn for use of any communicator differing from warnComm
+        static label warnComm;
+
 
     // Constructors
 
@@ -295,6 +320,79 @@ public:
 
     // Member functions
 
+        //- Allocate a new communicator
+        static label allocateCommunicator
+        (
+            const label parent,
+            const labelList& subRanks,
+            const bool doPstream = true
+        );
+
+        //- Free a previously allocated communicator
+        static void freeCommunicator
+        (
+            const label communicator,
+            const bool doPstream = true
+        );
+
+        //- Free all communicators
+        static void freeCommunicators(const bool doPstream);
+
+
+        //- Helper class for allocating/freeing communicators
+        class communicator
+        {
+            //- Communicator identifier
+            label comm_;
+
+            //- Disallow copy and assignment
+            communicator(const communicator&);
+            void operator=(const communicator&);
+
+        public:
+
+            //- Constructo from components
+            communicator
+            (
+                const label parent,
+                const labelList& subRanks,
+                const bool doPstream
+            )
+            :
+                comm_(allocateCommunicator(parent, subRanks, doPstream))
+            {}
+
+            //- Destructor
+            ~communicator()
+            {
+                freeCommunicator(comm_);
+            }
+
+            //- Cast to label
+            operator label() const
+            {
+                return comm_;
+            }
+        };
+
+
+        //- Return physical processor number (i.e. processor number in
+        //  worldComm) given communicator and processor
+        static int baseProcNo(const label myComm, const int procID);
+
+        //- Return processor number in communicator (given physical processor
+        //  number) (= reverse of baseProcNo)
+        static label procNo(const label comm, const int baseProcID);
+
+        //- Return processor number in communicator (given processor number
+        //  and communicator)
+        static label procNo
+        (
+            const label myComm,
+            const label currentComm,
+            const int currentProcID
+        );
+
         //- Add the valid option this type of communications library
         //  adds/requires on the command line
         static void addValidParOptions(HashTable<string>& validParOptions);
@@ -303,44 +401,78 @@ public:
         //  Spawns slave processes and initialises inter-communication
         static bool init(int& argc, char**& argv);
 
+
+        // Non-blocking comms
+
+            //- Get number of outstanding requests
+            static label nRequests();
+
+            //- Truncate number of outstanding requests
+            static void resetRequests(const label sz);
+
+            //- Wait until all requests (from start onwards) have finished.
+            static void waitRequests(const label start = 0);
+
+            //- Wait until request i has finished.
+            static void waitRequest(const label i);
+
+            //- Non-blocking comms: has request i finished?
+            static bool finishedRequest(const label i);
+
+            static int allocateTag(const char*);
+
+            static int allocateTag(const word&);
+
+            static void freeTag(const char*, const int tag);
+
+            static void freeTag(const word&, const int tag);
+
+
         //- Is this a parallel run?
         static bool& parRun()
         {
             return parRun_;
         }
 
-        //- Number of processes in parallel run
-        static label nProcs()
+        //- Number of processes in parallel run for a given communicator
+        static label nProcs(const label communicator = 0)
         {
-            return procIDs_.size();
+            return procIDs_[communicator].size();
         }
 
-        //- Am I the master process
-        static bool master()
-        {
-            return myProcNo_ == masterNo();
-        }
-
-        //- Process index of the master
+        //- Process index of the master for the global communicator
         static int masterNo()
         {
             return 0;
         }
 
-        //- Number of this process (starting from masterNo() = 0)
-        static int myProcNo()
+        //- Am I the master process
+        static bool master(const label communicator = 0)
         {
-            return myProcNo_;
+            return myProcNo_[communicator] == masterNo();
+        }
+
+
+        //- Number of this process (starting from masterNo() = 0)
+        static int myProcNo(const label communicator = 0)
+        {
+            return myProcNo_[communicator];
+        }
+
+        //- Return parent communicator
+        static label parent(const label communicator)
+        {
+            return parentCommunicator_(communicator);
         }
 
         //- Process IDs
-        static const List<int>& procIDs()
+        static const List<int>& procIDs(const int communicator)
         {
-            return procIDs_;
+            return procIDs_[communicator];
         }
 
         //- Process ID of given process index
-        static int procID(int procNo)
+        static List<int>& procID(const int procNo)
         {
             return procIDs_[procNo];
         }
@@ -352,21 +484,27 @@ public:
         }
 
         //- Process index of last slave
-        static int lastSlave()
+        static int lastSlave(const label communicator = 0)
         {
-            return nProcs() - 1;
+            return nProcs(communicator) - 1;
         }
 
         //- Communication schedule for linear all-to-master (proc 0)
-        static const List<commsStruct>& linearCommunication()
+        static const List<commsStruct>& linearCommunication
+        (
+            const label communicator = 0
+        )
         {
-            return linearCommunication_;
+            return linearCommunication_[communicator];
         }
 
         //- Communication schedule for tree all-to-master (proc 0)
-        static const List<commsStruct>& treeCommunication()
+        static const List<commsStruct>& treeCommunication
+        (
+            const label communicator = 0
+        )
         {
-            return treeCommunication_;
+            return treeCommunication_[communicator];
         }
 
         //- Message tag of standard messages
@@ -405,21 +543,40 @@ public:
             (
                 const List<commsStruct>& comms,
                 T& Value,
-                const BinaryOp& bop
+                const BinaryOp& bop,
+                const int tag,
+                const label comm
             );
 
             //- Like above but switches between linear/tree communication
             template <class T, class BinaryOp>
-            static void gather(T& Value, const BinaryOp& bop);
+            static void gather
+            (
+                T& Value,
+                const BinaryOp& bop,
+                const int tag = Pstream::msgType(),
+                const label comm = Pstream::worldComm
+            );
 
             //- Scatter data. Distribute without modification.
             //  Reverse of gather
             template <class T>
-            static void scatter(const List<commsStruct>& comms, T& Value);
+            static void scatter
+            (
+                const List<commsStruct>& comms,
+                T& Value,
+                const int tag,
+                const label comm
+            );
 
             //- Like above but switches between linear/tree communication
             template <class T>
-            static void scatter(T& Value);
+            static void scatter
+            (
+                T& Value,
+                const int tag = Pstream::msgType(),
+                const label comm = Pstream::worldComm
+            );
 
 
         // Combine variants. Inplace combine values from processors.
@@ -430,24 +587,39 @@ public:
             (
                 const List<commsStruct>& comms,
                 T& Value,
-                const CombineOp& cop
+                const CombineOp& cop,
+                const int tag,
+                const label comm
             );
 
             //- Like above but switches between linear/tree communication
             template <class T, class CombineOp>
-            static void combineGather(T& Value, const CombineOp& cop);
+            static void combineGather
+            (
+                T& Value,
+                const CombineOp& cop,
+                const int tag = Pstream::msgType(),
+                const label comm = Pstream::worldComm
+            );
 
             //- Scatter data. Reverse of combineGather
             template <class T>
             static void combineScatter
             (
                 const List<commsStruct>& comms,
-                T& Value
+                T& Value,
+                const int tag,
+                const label comm
             );
 
             //- Like above but switches between linear/tree communication
             template <class T>
-            static void combineScatter(T& Value);
+            static void combineScatter
+            (
+                T& Value,
+                const int tag = Pstream::msgType(),
+                const label comm = Pstream::worldComm
+            );
 
         // Combine variants working on whole List at a time.
 
@@ -456,7 +628,9 @@ public:
             (
                 const List<commsStruct>& comms,
                 List<T>& Value,
-                const CombineOp& cop
+                const CombineOp& cop,
+                const int tag,
+                const label comm
             );
 
             //- Like above but switches between linear/tree communication
@@ -464,7 +638,9 @@ public:
             static void listCombineGather
             (
                 List<T>& Value,
-                const CombineOp& cop
+                const CombineOp& cop,
+                const int tag = Pstream::msgType(),
+                const label comm = Pstream::worldComm
             );
 
             //- Scatter data. Reverse of combineGather
@@ -472,12 +648,19 @@ public:
             static void listCombineScatter
             (
                 const List<commsStruct>& comms,
-                List<T>& Value
+                List<T>& Value,
+                const int tag,
+                const label comm
             );
 
             //- Like above but switches between linear/tree communication
             template <class T>
-            static void listCombineScatter(List<T>& Value);
+            static void listCombineScatter
+            (
+                List<T>& Value,
+                const int tag = Pstream::msgType(),
+                const label comm = Pstream::worldComm
+            );
 
         // Combine variants working on whole map at a time. Container needs to
         // have iterators and find() defined.
@@ -487,7 +670,9 @@ public:
             (
                 const List<commsStruct>& comms,
                 Container& Values,
-                const CombineOp& cop
+                const CombineOp& cop,
+                const int tag,
+                const label comm
             );
 
             //- Like above but switches between linear/tree communication
@@ -495,7 +680,9 @@ public:
             static void mapCombineGather
             (
                 Container& Values,
-                const CombineOp& cop
+                const CombineOp& cop,
+                const int tag = Pstream::msgType(),
+                const label comm = Pstream::worldComm
             );
 
             //- Scatter data. Reverse of combineGather
@@ -503,13 +690,19 @@ public:
             static void mapCombineScatter
             (
                 const List<commsStruct>& comms,
-                Container& Values
+                Container& Values,
+                const int tag,
+                const label comm
             );
 
             //- Like above but switches between linear/tree communication
             template <class Container>
-            static void mapCombineScatter(Container& Values);
-
+            static void mapCombineScatter
+            (
+                Container& Values,
+                const int tag = Pstream::msgType(),
+                const label comm = Pstream::worldComm
+            );
 
 
         // Gather/scatter keeping the individual processor data separate.
@@ -521,24 +714,56 @@ public:
             static void gatherList
             (
                 const List<commsStruct>& comms,
-                List<T>& Values
+                List<T>& Values,
+                const int tag,
+                const label comm
             );
 
             //- Like above but switches between linear/tree communication
             template <class T>
-            static void gatherList(List<T>& Values);
+            static void gatherList
+            (
+                List<T>& Values,
+                const int tag = Pstream::msgType(),
+                const label comm = Pstream::worldComm
+            );
 
             //- Scatter data. Reverse of gatherList
             template <class T>
             static void scatterList
             (
                 const List<commsStruct>& comms,
-                List<T>& Values
+                List<T>& Values,
+                const int tag,
+                const label comm
             );
 
             //- Like above but switches between linear/tree communication
             template <class T>
-            static void scatterList(List<T>& Values);
+            static void scatterList
+            (
+                List<T>& Values,
+                const int tag = Pstream::msgType(),
+                const label comm = Pstream::worldComm
+            );
+
+
+        // Exchange
+
+            //- Exchange data. Sends sendData, receives into recvData, sets
+            //  sizes (not bytes). sizes[p0][p1] is what processor p0 has
+            //  sent to p1. Continuous data only.
+            //  If block=true will wait for all transfers to finish.
+            template<class Container, class T>
+            static void exchange
+            (
+                const List<Container>&,
+                List<Container>&,
+                labelListList& sizes,
+                const int tag = Pstream::msgType(),
+                const label comm = Pstream::worldComm,
+                const bool block = true
+            );
 };
 
 
@@ -561,6 +786,7 @@ Ostream& operator<<(Ostream&, const Pstream::commsStruct&);
 #   include "gatherScatter.C"
 #   include "combineGatherScatter.C"
 #   include "gatherScatterList.C"
+#   include "PstreamExchange.C"
 #endif
 
 

src/foam/db/IOstreams/Pstreams/PstreamCombineReduceOps.H

@@ -30,7 +30,6 @@ Description
     combined using the given combination function and the result is
     broadcast to all nodes
 
-
 \*---------------------------------------------------------------------------*/
 
 #ifndef PstreamCombineReduceOps_H
@@ -51,26 +50,62 @@ void combineReduce
 (
     const List<Pstream::commsStruct>& comms,
     T& Value,
-    const CombineOp& cop
+    const CombineOp& cop,
+    const int tag,
+    const label comm
 )
 {
-    Pstream::combineGather(comms, Value, cop);
-    Pstream::combineScatter(comms, Value);
+    Pstream::combineGather(comms, Value, cop, tag, comm);
+    Pstream::combineScatter(comms, Value, tag, comm);
 }
 
 
 template <class T, class CombineOp>
-void combineReduce(T& Value, const CombineOp& cop)
+void combineReduce
+(
+    T& Value,
+    const CombineOp& cop,
+    const int tag = Pstream::msgType(),
+    const label comm = Pstream::worldComm
+)
 {
     if (Pstream::nProcs() < Pstream::nProcsSimpleSum())
     {
-        Pstream::combineGather(Pstream::linearCommunication(), Value, cop);
-        Pstream::combineScatter(Pstream::linearCommunication(), Value);
+        Pstream::combineGather
+        (
+            Pstream::linearCommunication(comm),
+            Value,
+            cop,
+            tag,
+            comm
+        );
+
+        Pstream::combineScatter
+        (
+            Pstream::linearCommunication(comm),
+            Value,
+            tag,
+            comm
+        );
     }
     else
     {
-        Pstream::combineGather(Pstream::treeCommunication(), Value, cop);
-        Pstream::combineScatter(Pstream::treeCommunication(), Value);
+        Pstream::combineGather
+        (
+            Pstream::treeCommunication(comm),
+            Value,
+            cop,
+            tag,
+            comm
+        );
+
+        Pstream::combineScatter
+        (
+            Pstream::treeCommunication(comm),
+            Value,
+            tag,
+            comm
+        );
     }
 }
 

src/foam/db/IOstreams/Pstreams/PstreamCommsStruct.C

@@ -68,13 +68,13 @@ Foam::Pstream::commsStruct::commsStruct
 {
     boolList inBelow(nProcs, false);
 
-    forAll(allBelow, belowI)
+    forAll (allBelow, belowI)
     {
         inBelow[allBelow[belowI]] = true;
     }
 
     label notI = 0;
-    forAll(inBelow, procI)
+    forAll (inBelow, procI)
     {
         if ((procI != myProcID) && !inBelow[procI])
         {

src/foam/db/IOstreams/Pstreams/PstreamExchange.C

@@ -0,0 +1,159 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Description
+    Pstream exchange data.
+
+\*---------------------------------------------------------------------------*/
+
+#include "Pstream.H"
+#include "contiguous.H"
+#include "PstreamCombineReduceOps.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+//template<template<class> class ListType, class T>
+template<class Container, class T>
+void Pstream::exchange
+(
+    const List<Container>& sendBufs,
+    List<Container>& recvBufs,
+    labelListList& sizes,
+    const int tag,
+    const label comm,
+    const bool block
+)
+{
+    if (!contiguous<T>())
+    {
+        FatalErrorIn
+        (
+            "Pstream::exchange(..)"
+        )   << "Continuous data only." << Foam::abort(FatalError);
+    }
+
+    if (sendBufs.size() != Pstream::nProcs(comm))
+    {
+        FatalErrorIn
+        (
+            "Pstream::exchange(..)"
+        )   << "Size of list:" << sendBufs.size()
+            << " does not equal the number of processors:"
+            << Pstream::nProcs(comm)
+            << Foam::abort(FatalError);
+    }
+
+    sizes.setSize(Pstream::nProcs(comm));
+    labelList& nsTransPs = sizes[Pstream::myProcNo(comm)];
+    nsTransPs.setSize(Pstream::nProcs(comm));
+
+    forAll (sendBufs, procI)
+    {
+        nsTransPs[procI] = sendBufs[procI].size();
+    }
+
+    // Send sizes across. Note: blocks.
+    combineReduce(sizes, Pstream::listEq(), tag, comm);
+
+    if (Pstream::nProcs(comm) > 1)
+    {
+        label startOfRequests = Pstream::nRequests();
+
+        // Set up receives
+        // ~~~~~~~~~~~~~~~
+
+        recvBufs.setSize(sendBufs.size());
+        forAll (sizes, procI)
+        {
+            label nRecv = sizes[procI][Pstream::myProcNo(comm)];
+
+            if (procI != Pstream::myProcNo(comm) && nRecv > 0)
+            {
+                recvBufs[procI].setSize(nRecv);
+                IPstream::read
+                (
+                    Pstream::nonBlocking,
+                    procI,
+                    reinterpret_cast<char*>(recvBufs[procI].begin()),
+                    nRecv*sizeof(T),
+                    tag,
+                    comm
+                );
+            }
+        }
+
+
+        // Set up sends
+        // ~~~~~~~~~~~~
+
+        forAll (sendBufs, procI)
+        {
+            if (procI != Pstream::myProcNo(comm) && sendBufs[procI].size() > 0)
+            {
+                if
+                (
+                   !OPstream::write
+                    (
+                        Pstream::nonBlocking,
+                        procI,
+                        reinterpret_cast<const char*>(sendBufs[procI].begin()),
+                        sendBufs[procI].size()*sizeof(T),
+                        tag,
+                        comm
+                    )
+                )
+                {
+                    FatalErrorIn("Pstream::exchange(..)")
+                        << "Cannot send outgoing message. "
+                        << "to:" << procI << " nBytes:"
+                        << label(sendBufs[procI].size()*sizeof(T))
+                        << Foam::abort(FatalError);
+                }
+            }
+        }
+
+
+        // Wait for all to finish
+        // ~~~~~~~~~~~~~~~~~~~~~~
+
+        if (block)
+        {
+            Pstream::waitRequests(startOfRequests);
+        }
+    }
+
+    // Do myself
+    recvBufs[Pstream::myProcNo(comm)] = sendBufs[Pstream::myProcNo(comm)];
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/foam/db/IOstreams/Pstreams/PstreamGlobals.C

@@ -33,10 +33,50 @@ namespace Foam
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 // Outstanding non-blocking operations.
-//! @cond fileScope
-DynamicList<MPI_Request> PstreamGlobals::IPstream_outstandingRequests_;
-DynamicList<MPI_Request> PstreamGlobals::OPstream_outstandingRequests_;
-//! @endcond
+//! \cond fileScope
+DynamicList<MPI_Request> PstreamGlobals::outstandingRequests_;
+//! \endcond
+
+// Max outstanding message tag operations.
+//! \cond fileScope
+int PstreamGlobals::nTags_ = 0;
+//! \endcond
+
+// Free'd message tags
+//! \cond fileScope
+DynamicList<int> PstreamGlobals::freedTags_;
+//! \endcond
+
+
+// Allocated communicators.
+//! \cond fileScope
+DynamicList<MPI_Comm> PstreamGlobals::MPICommunicators_;
+DynamicList<MPI_Group> PstreamGlobals::MPIGroups_;
+//! \endcond
+
+void PstreamGlobals::checkCommunicator
+(
+    const label comm,
+    const label otherProcNo
+)
+{
+    if
+    (
+        comm < 0
+     || comm >= PstreamGlobals::MPICommunicators_.size()
+    )
+    {
+        FatalErrorIn
+        (
+            "PstreamGlobals::checkCommunicator(const label, const label)"
+        )   << "otherProcNo:" << otherProcNo << " : illegal communicator "
+            << comm << endl
+            << "Communicator should be within range 0.."
+            << PstreamGlobals::MPICommunicators_.size() - 1
+            << abort(FatalError);
+    }
+}
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/foam/db/IOstreams/Pstreams/PstreamGlobals.H

@@ -52,8 +52,18 @@ namespace Foam
 namespace PstreamGlobals
 {
 
-extern DynamicList<MPI_Request> IPstream_outstandingRequests_;
-extern DynamicList<MPI_Request> OPstream_outstandingRequests_;
+extern DynamicList<MPI_Request> outstandingRequests_;
+
+extern int nTags_;
+
+extern DynamicList<int> freedTags_;
+
+
+// Current communicators. First element will be MPI_COMM_WORLD
+extern DynamicList<MPI_Comm> MPICommunicators_;
+extern DynamicList<MPI_Group> MPIGroups_;
+
+void checkCommunicator(const label, const label procNo);
 
 };
 

src/foam/db/IOstreams/Pstreams/PstreamReduceOps.C

@@ -0,0 +1,732 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "mpi.h"
+
+#include "label.H"
+#include "Pstream.H"
+#include "PstreamReduceOps.H"
+#include "allReduce.H"
+
+// Check type of label for use in MPI calls
+#if defined(WM_INT)
+#   define MPI_LABEL MPI_INT
+#elif defined(WM_LONG)
+#   define MPI_LABEL MPI_LONG
+#elif defined(WM_LLONG)
+#   define MPI_LABEL MPI_LONG_LONG
+#endif
+
+// Check type of scalar for use in MPI calls
+#if defined(WM_SP)
+#   define MPI_SCALAR MPI_FLOAT
+#elif defined(WM_DP)
+#   define MPI_SCALAR MPI_DOUBLE
+#elif defined(WM_LDP)
+#   define MPI_SCALAR MPI_LONG_DOUBLE
+#endif
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+// Optimizing template specialisations using MPI_REDUCE
+void Foam::reduce
+(
+    bool& Value,
+    const andOp<bool>& bop,
+    const int tag,
+    const label comm
+)
+{
+    if (Pstream::warnComm != -1 && comm != Pstream::warnComm)
+    {
+        Pout<< "** reducing:" << Value << " with comm:" << comm
+            << " warnComm:" << Pstream::warnComm
+            << endl;
+        error::printStack(Pout);
+    }
+
+#   ifdef FULLDEBUG
+    // Check for processors that are not in the communicator
+    if (Pstream::myProcNo(comm) == -1)
+    {
+        FatalErrorIn
+        (
+            "void Foam::reduce\n"
+            "(\n"
+            "    bool& Value,\n"
+            "    const andOp<bool>& bop,\n"
+            "    const int tag,\n"
+            "    const label comm\n"
+            ")"
+        )   << "Reduce called on the processor which is not a member "
+            << "of comm.  This is not allowed"
+            << abort(FatalError);
+    }
+#   endif
+
+    // Note: C++ bool is a type separate from C and cannot be cast
+    // For safety and compatibility with compilers, convert bool to int
+    // to comply with MPI types.  HJ, 23/Sep/2016
+
+    int intBool = 0;
+
+    if (Value)
+    {
+        intBool = 1;
+    }
+
+    allReduce(intBool, 1, MPI_INT, MPI_LAND, bop, tag, comm);
+
+    if (intBool > 0)
+    {
+        Value = true;
+    }
+    else
+    {
+        Value = false;
+    }
+}
+
+
+void Foam::reduce
+(
+    bool& Value,
+    const orOp<bool>& bop,
+    const int tag,
+    const label comm
+)
+{
+    if (Pstream::warnComm != -1 && comm != Pstream::warnComm)
+    {
+        Pout<< "** reducing:" << Value << " with comm:" << comm
+            << " warnComm:" << Pstream::warnComm
+            << endl;
+        error::printStack(Pout);
+    }
+
+#   ifdef FULLDEBUG
+    // Check for processors that are not in the communicator
+    if (Pstream::myProcNo(comm) == -1)
+    {
+        FatalErrorIn
+        (
+            "void Foam::reduce\n"
+            "(\n"
+            "    bool& Value,\n"
+            "    const orOp<bool>& bop,\n"
+            "    const int tag,\n"
+            "    const label comm\n"
+            ")"
+        )   << "Reduce called on the processor which is not a member "
+            << "of comm.  This is not allowed"
+            << abort(FatalError);
+    }
+#   endif
+
+    // Note: C++ bool is a type separate from C and cannot be cast
+    // For safety and compatibility with compilers, convert bool to int
+    // to comply with MPI types.  HJ, 23/Sep/2016
+
+    int intBool = 0;
+
+    if (Value)
+    {
+        intBool = 1;
+    }
+
+    allReduce(intBool, 1, MPI_INT, MPI_LOR, bop, tag, comm);
+
+    if (intBool > 0)
+    {
+        Value = true;
+    }
+    else
+    {
+        Value = false;
+    }
+}
+
+
+void Foam::reduce
+(
+    label& Value,
+    const minOp<label>& bop,
+    const int tag,
+    const label comm
+)
+{
+    if (Pstream::warnComm != -1 && comm != Pstream::warnComm)
+    {
+        Pout<< "** reducing:" << Value << " with comm:" << comm
+            << " warnComm:" << Pstream::warnComm
+            << endl;
+        error::printStack(Pout);
+    }
+
+#   ifdef FULLDEBUG
+    // Check for processors that are not in the communicator
+    if (Pstream::myProcNo(comm) == -1)
+    {
+        FatalErrorIn
+        (
+            "void Foam::reduce\n"
+            "(\n"
+            "    label& Value,\n"
+            "    const minOp<label>& bop,\n"
+            "    const int tag,\n"
+            "    const label comm\n"
+            ")"
+        )   << "Reduce called on the processor which is not a member "
+            << "of comm.  This is not allowed"
+            << abort(FatalError);
+    }
+#   endif
+
+    allReduce(Value, 1, MPI_LABEL, MPI_MIN, bop, tag, comm);
+}
+
+
+void Foam::reduce
+(
+    label& Value,
+    const maxOp<label>& bop,
+    const int tag,
+    const label comm
+)
+{
+    if (Pstream::warnComm != -1 && comm != Pstream::warnComm)
+    {
+        Pout<< "** reducing:" << Value << " with comm:" << comm
+            << " warnComm:" << Pstream::warnComm
+            << endl;
+        error::printStack(Pout);
+    }
+
+#   ifdef FULLDEBUG
+    // Check for processors that are not in the communicator
+    if (Pstream::myProcNo(comm) == -1)
+    {
+        FatalErrorIn
+        (
+            "void Foam::reduce\n"
+            "(\n"
+            "    label& Value,\n"
+            "    const maxOp<label>& bop,\n"
+            "    const int tag,\n"
+            "    const label comm\n"
+            ")"
+        )   << "Reduce called on the processor which is not a member "
+            << "of comm.  This is not allowed"
+            << abort(FatalError);
+    }
+#   endif
+
+    allReduce(Value, 1, MPI_LABEL, MPI_MAX, bop, tag, comm);
+}
+
+
+void Foam::reduce
+(
+    label& Value,
+    const sumOp<label>& bop,
+    const int tag,
+    const label comm
+)
+{
+    if (Pstream::warnComm != -1 && comm != Pstream::warnComm)
+    {
+        Pout<< "** reducing:" << Value << " with comm:" << comm
+            << " warnComm:" << Pstream::warnComm
+            << endl;
+        error::printStack(Pout);
+    }
+
+#   ifdef FULLDEBUG
+    // Check for processors that are not in the communicator
+    if (Pstream::myProcNo(comm) == -1)
+    {
+        FatalErrorIn
+        (
+            "void Foam::reduce\n"
+            "(\n"
+            "    label& Value,\n"
+            "    const sumOp<label>& bop,\n"
+            "    const int tag,\n"
+            "    const label comm\n"
+            ")"
+        )   << "Reduce called on the processor which is not a member "
+            << "of comm.  This is not allowed"
+            << abort(FatalError);
+    }
+#   endif
+
+    allReduce(Value, 1, MPI_LABEL, MPI_SUM, bop, tag, comm);
+}
+
+
+void Foam::reduce
+(
+    scalar& Value,
+    const minOp<scalar>& bop,
+    const int tag,
+    const label comm
+)
+{
+    if (Pstream::warnComm != -1 && comm != Pstream::warnComm)
+    {
+        Pout<< "** reducing:" << Value << " with comm:" << comm
+            << " warnComm:" << Pstream::warnComm
+            << endl;
+        error::printStack(Pout);
+    }
+
+#   ifdef FULLDEBUG
+    // Check for processors that are not in the communicator
+    if (Pstream::myProcNo(comm) == -1)
+    {
+        FatalErrorIn
+        (
+            "void Foam::reduce\n"
+            "(\n"
+            "    scalar& Value,\n"
+            "    const minOp<scalar>& bop,\n"
+            "    const int tag,\n"
+            "    const label comm\n"
+            ")"
+        )   << "Reduce called on the processor which is not a member "
+            << "of comm.  This is not allowed"
+            << abort(FatalError);
+    }
+#   endif
+
+    allReduce(Value, 1, MPI_SCALAR, MPI_MIN, bop, tag, comm);
+}
+
+
+void Foam::reduce
+(
+    scalar& Value,
+    const maxOp<scalar>& bop,
+    const int tag,
+    const label comm
+)
+{
+    if (Pstream::warnComm != -1 && comm != Pstream::warnComm)
+    {
+        Pout<< "** reducing:" << Value << " with comm:" << comm
+            << " warnComm:" << Pstream::warnComm
+            << endl;
+        error::printStack(Pout);
+    }
+
+#   ifdef FULLDEBUG
+    // Check for processors that are not in the communicator
+    if (Pstream::myProcNo(comm) == -1)
+    {
+        FatalErrorIn
+        (
+            "void Foam::reduce\n"
+            "(\n"
+            "    scalar& Value,\n"
+            "    const maxOp<scalar>& bop,\n"
+            "    const int tag,\n"
+            "    const label comm\n"
+            ")"
+        )   << "Reduce called on the processor which is not a member "
+            << "of comm.  This is not allowed"
+            << abort(FatalError);
+    }
+#   endif
+
+    allReduce(Value, 1, MPI_SCALAR, MPI_MAX, bop, tag, comm);
+}
+
+
+void Foam::reduce
+(
+    scalar& Value,
+    const sumOp<scalar>& bop,
+    const int tag,
+    const label comm
+)
+{
+    if (Pstream::warnComm != -1 && comm != Pstream::warnComm)
+    {
+        Pout<< "** reducing:" << Value << " with comm:" << comm
+            << " warnComm:" << Pstream::warnComm
+            << endl;
+        error::printStack(Pout);
+    }
+
+#   ifdef FULLDEBUG
+    // Check for processors that are not in the communicator
+    if (Pstream::myProcNo(comm) == -1)
+    {
+        FatalErrorIn
+        (
+            "void Foam::reduce\n"
+            "(\n"
+            "    scalar& Value,\n"
+            "    const sumOp<scalar>& bop,\n"
+            "    const int tag,\n"
+            "    const label comm\n"
+            ")"
+        )   << "Reduce called on the processor which is not a member "
+            << "of comm.  This is not allowed"
+            << abort(FatalError);
+    }
+#   endif
+
+    allReduce(Value, 1, MPI_SCALAR, MPI_SUM, bop, tag, comm);
+}
+
+
+void Foam::reduce
+(
+    List<label>& Value,
+    const minOp<List<label> >& bop,
+    const int tag,
+    const label comm
+)
+{
+    if (!Pstream::parRun())
+    {
+        return;
+    }
+
+    if (Pstream::warnComm != -1 && comm != Pstream::warnComm)
+    {
+        Pout<< "** reducing:" << Value << " with comm:" << comm
+            << " warnComm:" << Pstream::warnComm
+            << endl;
+        error::printStack(Pout);
+    }
+
+#   ifdef FULLDEBUG
+    // Check for processors that are not in the communicator
+    if (Pstream::myProcNo(comm) == -1)
+    {
+        FatalErrorIn
+        (
+            "void Foam::reduce\n"
+            "(\n"
+            "    List<label>& Value,\n"
+            "    const minOp<List<label> >& bop,\n"
+            "    const int tag,\n"
+            "    const label comm\n"
+            ")"
+        )   << "Reduce called on the processor which is not a member "
+            << "of comm.  This is not allowed"
+            << abort(FatalError);
+    }
+#   endif
+
+    // Make a copy of the Value in the send buffer so that Value can be
+    // used for receive.  HJ, 8/Oct/2016
+    labelList send(Value);
+
+    int MPISize = Value.size();
+
+    MPI_Allreduce
+    (
+        send.begin(),
+        Value.begin(),
+        MPISize,
+        MPI_LABEL,
+        MPI_MIN,
+        PstreamGlobals::MPICommunicators_[comm]
+    );
+}
+
+
+void Foam::reduce
+(
+    List<label>& Value,
+    const maxOp<List<label> >& bop,
+    const int tag,
+    const label comm
+)
+{
+    if (!Pstream::parRun())
+    {
+        return;
+    }
+
+    if (Pstream::warnComm != -1 && comm != Pstream::warnComm)
+    {
+        Pout<< "** reducing:" << Value << " with comm:" << comm
+            << " warnComm:" << Pstream::warnComm
+            << endl;
+        error::printStack(Pout);
+    }
+
+#   ifdef FULLDEBUG
+    // Check for processors that are not in the communicator
+    if (Pstream::myProcNo(comm) == -1)
+    {
+        FatalErrorIn
+        (
+            "void Foam::reduce\n"
+            "(\n"
+            "    List<label>& Value,\n"
+            "    const maxOp<List<label> >& bop,\n"
+            "    const int tag,\n"
+            "    const label comm\n"
+            ")"
+        )   << "Reduce called on the processor which is not a member "
+            << "of comm.  This is not allowed"
+            << abort(FatalError);
+    }
+#   endif
+
+    // Make a copy of the Value in the send buffer so that Value can be
+    // used for receive.  HJ, 8/Oct/2016
+    labelList send(Value);
+
+    int MPISize = Value.size();
+
+    MPI_Allreduce
+    (
+        send.begin(),
+        Value.begin(),
+        MPISize,
+        MPI_LABEL,
+        MPI_MAX,
+        PstreamGlobals::MPICommunicators_[comm]
+    );
+}
+
+
+void Foam::reduce
+(
+    List<label>& Value,
+    const sumOp<List<label> >& bop,
+    const int tag,
+    const label comm
+)
+{
+    if (!Pstream::parRun())
+    {
+        return;
+    }
+
+    if (Pstream::warnComm != -1 && comm != Pstream::warnComm)
+    {
+        Pout<< "** reducing:" << Value << " with comm:" << comm
+            << " warnComm:" << Pstream::warnComm
+            << endl;
+        error::printStack(Pout);
+    }
+
+#   ifdef FULLDEBUG
+    // Check for processors that are not in the communicator
+    if (Pstream::myProcNo(comm) == -1)
+    {
+        FatalErrorIn
+        (
+            "void Foam::reduce\n"
+            "(\n"
+            "    List<label>& Value,\n"
+            "    const sumOp<List<label> >& bop,\n"
+            "    const int tag,\n"
+            "    const label comm\n"
+            ")"
+        )   << "Reduce called on the processor which is not a member "
+            << "of comm.  This is not allowed"
+            << abort(FatalError);
+    }
+#   endif
+
+    // Make a copy of the Value in the send buffer so that Value can be
+    // used for receive.  HJ, 8/Oct/2016
+    labelList send(Value);
+
+    int MPISize = Value.size();
+
+    MPI_Allreduce
+    (
+        send.begin(),
+        Value.begin(),
+        MPISize,
+        MPI_LABEL,
+        MPI_SUM,
+        PstreamGlobals::MPICommunicators_[comm]
+    );
+}
+
+
+void Foam::reduce
+(
+    vector2D& Value,
+    const sumOp<vector2D>& bop,
+    const int tag,
+    const label comm
+)
+{
+    if (Pstream::warnComm != -1 && comm != Pstream::warnComm)
+    {
+        Pout<< "** reducing:" << Value << " with comm:" << comm
+            << " warnComm:" << Pstream::warnComm
+            << endl;
+        error::printStack(Pout);
+    }
+
+#   ifdef FULLDEBUG
+    // Check for processors that are not in the communicator
+    if (Pstream::myProcNo(comm) == -1)
+    {
+        FatalErrorIn
+        (
+            "void Foam::reduce\n"
+            "(\n"
+            "    vector2D& Value,\n"
+            "    const sumOp<vector2D>& bop,\n"
+            "    const int tag,\n"
+            "    const label comm\n"
+            ")"
+        )   << "Reduce called on the processor which is not a member "
+            << "of comm.  This is not allowed"
+            << abort(FatalError);
+    }
+#   endif
+
+    allReduce(Value, 2, MPI_SCALAR, MPI_SUM, bop, tag, comm);
+}
+
+
+void Foam::sumReduce
+(
+    scalar& Value,
+    label& Count,
+    const int tag,
+    const label comm
+)
+{
+    if (Pstream::warnComm != -1 && comm != Pstream::warnComm)
+    {
+        Pout<< "** reducing:" << Value << " with comm:" << comm
+            << " warnComm:" << Pstream::warnComm
+            << endl;
+        error::printStack(Pout);
+    }
+
+#   ifdef FULLDEBUG
+    // Check for processors that are not in the communicator
+    if (Pstream::myProcNo(comm) == -1)
+    {
+        FatalErrorIn
+        (
+            "void Foam::sumReduce\n"
+            "(\n"
+            "    scalar& Value,\n"
+            "    label& Count,\n"
+            "    const int tag,\n"
+            "    const label comm\n"
+            ")"
+        )   << "Reduce called on the processor which is not a member "
+            << "of comm.  This is not allowed"
+            << abort(FatalError);
+    }
+#   endif
+
+    vector2D twoScalars(Value, scalar(Count));
+    reduce(twoScalars, sumOp<vector2D>(), tag, comm);
+
+    Value = twoScalars.x();
+    Count = twoScalars.y();
+}
+
+
+void Foam::reduce
+(
+    scalar& Value,
+    const sumOp<scalar>& bop,
+    const int tag,
+    const label comm,
+    label& requestID
+)
+{
+    if (!Pstream::parRun())
+    {
+        return;
+    }
+
+#   ifdef FULLDEBUG
+    // Check for processors that are not in the communicator
+    if (Pstream::myProcNo(comm) == -1)
+    {
+        FatalErrorIn
+        (
+            "void Foam::reduce\n"
+            "(\n"
+            "    scalar& Value,\n"
+            "    const sumOp<scalar>& bop,\n"
+            "    const int tag,\n"
+            "    const label comm,\n"
+            "    label& requestID\n"
+            ")"
+        )   << "Reduce called on the processor which is not a member "
+            << "of comm.  This is not allowed"
+            << abort(FatalError);
+    }
+#   endif
+
+#ifdef MPIX_COMM_TYPE_SHARED
+    // Assume mpich2 with non-blocking collectives extensions. Once mpi3
+    // is available this will change.
+    MPI_Request request;
+    scalar v = Value;
+    MPIX_Ireduce
+    (
+        &v,
+        &Value,
+        1,
+        MPI_SCALAR,
+        MPI_SUM,
+        0,              // root
+        PstreamGlobals::MPICommunicators_[comm],
+        &request
+    );
+
+    requestID = PstreamGlobals::outstandingRequests_.size();
+    PstreamGlobals::outstandingRequests_.append(request);
+
+    if (debug)
+    {
+        Pout<< "Pstream::allocateRequest for non-blocking reduce"
+            << " : request:" << requestID
+            << endl;
+    }
+#else
+    // Non-blocking not yet implemented in mpi
+    reduce(Value, bop, tag, comm);
+
+    requestID = -1;
+#endif
+}
+
+
+// ************************************************************************* //

src/foam/db/IOstreams/Pstreams/PstreamReduceOps.H

@@ -28,6 +28,7 @@ License
 
 #include "Pstream.H"
 #include "ops.H"
+#include "vector2D.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -42,11 +43,20 @@ void reduce
 (
     const List<Pstream::commsStruct>& comms,
     T& Value,
-    const BinaryOp& bop
+    const BinaryOp& bop,
+    const int tag,
+    const label comm
 )
 {
-    Pstream::gather(comms, Value, bop);
-    Pstream::scatter(comms, Value);
+    if (Pstream::warnComm != -1 && comm != Pstream::warnComm)
+    {
+        Pout<< "** reducing:" << Value << " with comm:" << comm
+            << endl;
+        error::printStack(Pout);
+    }
+
+    Pstream::gather(comms, Value, bop, tag, comm);
+    Pstream::scatter(comms, Value, tag, comm);
 }
 
 
@@ -55,16 +65,18 @@ template <class T, class BinaryOp>
 void reduce
 (
     T& Value,
-    const BinaryOp& bop
+    const BinaryOp& bop,
+    const int tag = Pstream::msgType(),
+    const label comm = Pstream::worldComm
 )
 {
-    if (Pstream::nProcs() < Pstream::nProcsSimpleSum())
+    if (Pstream::nProcs(comm) < Pstream::nProcsSimpleSum())
     {
-        reduce(Pstream::linearCommunication(), Value, bop);
+        reduce(Pstream::linearCommunication(comm), Value, bop, tag, comm);
     }
     else
     {
-        reduce(Pstream::treeCommunication(), Value, bop);
+        reduce(Pstream::treeCommunication(comm), Value, bop, tag, comm);
     }
 }
 
@@ -74,26 +86,185 @@ template <class T, class BinaryOp>
 T returnReduce
 (
     const T& Value,
-    const BinaryOp& bop
+    const BinaryOp& bop,
+    const int tag = Pstream::msgType(),
+    const label comm = Pstream::worldComm
 )
 {
     T WorkValue(Value);
 
-    if (Pstream::nProcs() < Pstream::nProcsSimpleSum())
+    if (Pstream::nProcs(comm) < Pstream::nProcsSimpleSum())
     {
-        reduce(Pstream::linearCommunication(), WorkValue, bop);
+        reduce(Pstream::linearCommunication(comm), WorkValue, bop, tag, comm);
     }
     else
     {
-        reduce(Pstream::treeCommunication(), WorkValue, bop);
+        reduce(Pstream::treeCommunication(comm), WorkValue, bop, tag, comm);
     }
 
     return WorkValue;
 }
 
 
-// Insist there is a specialisation for the reduction of a scalar
-void reduce(scalar& Value, const sumOp<scalar>& bop);
+// Reduce with sum of both value and count (for averaging)
+template<class T>
+void sumReduce
+(
+    T& Value,
+    label& Count,
+    const int tag = Pstream::msgType(),
+    const label comm = Pstream::worldComm
+)
+{
+    reduce(Value, sumOp<T>(), tag, comm);
+    reduce(Count, sumOp<label>(), tag, comm);
+}
+
+
+// Non-blocking version of reduce. Sets request.
+template<class T, class BinaryOp>
+void reduce
+(
+    T& Value,
+    const BinaryOp& bop,
+    const int tag,
+    const label comm,
+    label& request
+)
+{
+    notImplemented
+    (
+        "reduce(T&, const BinaryOp&, const int, const label, label&)"
+    );
+}
+
+
+// Insist there are specialisations for the common reductions of bool
+// See implementation notes on casting of a C++ bool into MPI_INT in
+// Pstream.C  HJ, 23/Sep/2016
+
+void reduce
+(
+    bool& Value,
+    const andOp<bool>& bop,
+    const int tag = Pstream::msgType(),
+    const label comm = Pstream::worldComm
+);
+
+void reduce
+(
+    bool& Value,
+    const orOp<bool>& bop,
+    const int tag = Pstream::msgType(),
+    const label comm = Pstream::worldComm
+);
+
+
+// Insist there are specialisations for the common reductions of labels
+
+void reduce
+(
+    label& Value,
+    const sumOp<label>& bop,
+    const int tag = Pstream::msgType(),
+    const label comm = Pstream::worldComm
+);
+
+void reduce
+(
+    label& Value,
+    const minOp<label>& bop,
+    const int tag = Pstream::msgType(),
+    const label comm = Pstream::worldComm
+);
+
+void reduce
+(
+    label& Value,
+    const maxOp<label>& bop,
+    const int tag = Pstream::msgType(),
+    const label comm = Pstream::worldComm
+);
+
+// Insist there are specialisations for the common reductions of scalars
+
+void reduce
+(
+    scalar& Value,
+    const sumOp<scalar>& bop,
+    const int tag = Pstream::msgType(),
+    const label comm = Pstream::worldComm
+);
+
+void reduce
+(
+    scalar& Value,
+    const minOp<scalar>& bop,
+    const int tag = Pstream::msgType(),
+    const label comm = Pstream::worldComm
+);
+
+void reduce
+(
+    scalar& Value,
+    const maxOp<scalar>& bop,
+    const int tag = Pstream::msgType(),
+    const label comm = Pstream::worldComm
+);
+
+void reduce
+(
+    vector2D& Value,
+    const sumOp<vector2D>& bop,
+    const int tag = Pstream::msgType(),
+    const label comm = Pstream::worldComm
+);
+
+void sumReduce
+(
+    scalar& Value,
+    label& Count,
+    const int tag = Pstream::msgType(),
+    const label comm = Pstream::worldComm
+);
+
+void reduce
+(
+    scalar& Value,
+    const sumOp<scalar>& bop,
+    const int tag,
+    const label comm,
+    label& request
+);
+
+
+// Insist there are specialisations for the common reductions of
+// lists of labels.  Note: template function specialisation must be the
+// exact match on argument types.  HJ, 8/Oct/2016
+
+void reduce
+(
+    List<label>& Value,
+    const sumOp<List<label> >& bop,
+    const int tag = Pstream::msgType(),
+    const label comm = Pstream::worldComm
+);
+
+void reduce
+(
+    List<label>& Value,
+    const minOp<List<label> >& bop,
+    const int tag = Pstream::msgType(),
+    const label comm = Pstream::worldComm
+);
+
+void reduce
+(
+    List<label>& Value,
+    const maxOp<List<label> >& bop,
+    const int tag = Pstream::msgType(),
+    const label comm = Pstream::worldComm
+);
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/foam/db/IOstreams/Pstreams/allReduce.H

@@ -0,0 +1,72 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Description
+    Various functions to wrap MPI_Allreduce
+
+SourceFiles
+    allReduceTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef allReduce_H
+#define allReduce_H
+
+#include "mpi.h"
+
+#include "Pstream.H"
+#include "PstreamGlobals.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+template<class Type, class BinaryOp>
+void allReduce
+(
+    Type& Value,
+    int count,
+    MPI_Datatype MPIType,
+    MPI_Op op,
+    const BinaryOp& bop,
+    const int tag,
+    const int communicator
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "allReduceTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/foam/db/IOstreams/Pstreams/allReduceTemplates.C

@@ -0,0 +1,88 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "allReduce.H"
+
+// * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
+
+template<class Type, class BinaryOp>
+void Foam::allReduce
+(
+    Type& Value,
+    int MPICount,
+    MPI_Datatype MPIType,
+    MPI_Op MPIOp,
+    const BinaryOp& bop,
+    const int tag,
+    const label comm
+)
+{
+    if (!Pstream::parRun())
+    {
+        return;
+    }
+
+    // Removed send-received loop: use Allreduce instead.
+    // HJ, 8/Oct/2016
+
+#   ifdef FULLDEBUG
+    // Check for processors that are not in the communicator
+    if (Pstream::myProcNo(comm) == -1)
+    {
+        FatalErrorIn
+        (
+            "void Foam::allReduce\n"
+            "(\n"
+            "    Type& Value,\n"
+            "    int MPICount,\n"
+            "    MPI_Datatype MPIType,\n"
+            "    MPI_Op MPIOp,\n"
+            "    const BinaryOp& bop,\n"
+            "    const int tag,\n"
+            "    const label comm\n"
+            ")"
+        )   << "Reduce called on the processor which is not a member "
+            << "of comm.  This is not allowed"
+            << abort(FatalError);
+    }
+#   endif
+
+    Type sum;
+
+    MPI_Allreduce
+    (
+        &Value,
+        &sum,
+        MPICount,
+        MPIType,
+        MPIOp,
+        PstreamGlobals::MPICommunicators_[comm]
+    );
+
+    Value = sum;
+}
+
+
+// ************************************************************************* //

src/foam/db/IOstreams/Pstreams/combineGatherScatter.C

@@ -50,16 +50,18 @@ void Pstream::combineGather
 (
     const List<Pstream::commsStruct>& comms,
     T& Value,
-    const CombineOp& cop
+    const CombineOp& cop,
+    const int tag,
+    const label comm
 )
 {
-    if (Pstream::parRun())
+    if (Pstream::nProcs(comm) > 1)
     {
         // Get my communication order
-        const commsStruct& myComm = comms[Pstream::myProcNo()];
+        const commsStruct& myComm = comms[Pstream::myProcNo(comm)];
 
         // Receive from my downstairs neighbours
-        forAll(myComm.below(), belowI)
+        forAll (myComm.below(), belowI)
         {
             label belowID = myComm.below()[belowI];
 
@@ -71,7 +73,9 @@ void Pstream::combineGather
                     Pstream::scheduled,
                     belowID,
                     reinterpret_cast<char*>(&value),
-                    sizeof(T)
+                    sizeof(T),
+                    tag,
+                    comm
                 );
 
                 if (debug > 1)
@@ -84,7 +88,7 @@ void Pstream::combineGather
             }
             else
             {
-                IPstream fromBelow(Pstream::scheduled, belowID);
+                IPstream fromBelow(Pstream::scheduled, belowID, 0, tag, comm);
                 T value(fromBelow);
 
                 if (debug > 1)
@@ -113,12 +117,21 @@ void Pstream::combineGather
                     Pstream::scheduled,
                     myComm.above(),
                     reinterpret_cast<const char*>(&Value),
-                    sizeof(T)
+                    sizeof(T),
+                    tag,
+                    comm
                 );
             }
             else
             {
-                OPstream toAbove(Pstream::scheduled, myComm.above());
+                OPstream toAbove
+                (
+                    Pstream::scheduled,
+                    myComm.above(),
+                    0,
+                    tag,
+                    comm
+                );
                 toAbove << Value;
             }
         }
@@ -127,26 +140,52 @@ void Pstream::combineGather
 
 
 template <class T, class CombineOp>
-void Pstream::combineGather(T& Value, const CombineOp& cop)
+void Pstream::combineGather
+(
+    T& Value,
+    const CombineOp& cop,
+    const int tag,
+    const label comm
+)
 {
-    if (Pstream::nProcs() < Pstream::nProcsSimpleSum())
+    if (Pstream::nProcs(comm) < Pstream::nProcsSimpleSum())
     {
-        combineGather(Pstream::linearCommunication(), Value, cop);
+        combineGather
+        (
+            Pstream::linearCommunication(comm),
+            Value,
+            cop,
+            tag,
+            comm
+        );
     }
     else
     {
-        combineGather(Pstream::treeCommunication(), Value, cop);
+        combineGather
+        (
+            Pstream::treeCommunication(comm),
+            Value,
+            cop,
+            tag,
+            comm
+        );
     }
 }
 
 
 template <class T>
-void Pstream::combineScatter(const List<Pstream::commsStruct>& comms, T& Value)
+void Pstream::combineScatter
+(
+    const List<Pstream::commsStruct>& comms,
+    T& Value,
+    const int tag,
+    const label comm
+)
 {
-    if (Pstream::parRun())
+    if (Pstream::nProcs(comm) > 1)
     {
         // Get my communication order
-        const Pstream::commsStruct& myComm = comms[Pstream::myProcNo()];
+        const Pstream::commsStruct& myComm = comms[Pstream::myProcNo(comm)];
 
         // Receive from up
         if (myComm.above() != -1)
@@ -158,12 +197,21 @@ void Pstream::combineScatter(const List<Pstream::commsStruct>& comms, T& Value)
                     Pstream::scheduled,
                     myComm.above(),
                     reinterpret_cast<char*>(&Value),
-                    sizeof(T)
+                    sizeof(T),
+                    tag,
+                    comm
                 );
             }
             else
             {
-                IPstream fromAbove(Pstream::scheduled, myComm.above());
+                IPstream fromAbove
+                (
+                    Pstream::scheduled,
+                    myComm.above(),
+                    0,
+                    tag,
+                    comm
+                );
                 Value = T(fromAbove);
             }
 
@@ -174,8 +222,11 @@ void Pstream::combineScatter(const List<Pstream::commsStruct>& comms, T& Value)
             }
         }
 
-        // Send to my downstairs neighbours
-        forAll(myComm.below(), belowI)
+        // Send to my downstairs neighbours.  Note reverse order (compared to
+        // receiving). This is to make sure to send to the critical path
+        // (only when using a tree schedule!) first.
+        // This is ESI Comms optimisation, v16.06.  HJ, 19/Sep/2016
+        forAllReverse (myComm.below(), belowI)
         {
             label belowID = myComm.below()[belowI];
 
@@ -191,12 +242,14 @@ void Pstream::combineScatter(const List<Pstream::commsStruct>& comms, T& Value)
                     Pstream::scheduled,
                     belowID,
                     reinterpret_cast<const char*>(&Value),
-                    sizeof(T)
+                    sizeof(T),
+                    tag,
+                    comm
                 );
             }
             else
             {
-                OPstream toBelow(Pstream::scheduled, belowID);
+                OPstream toBelow(Pstream::scheduled, belowID, 0, tag, comm);
                 toBelow << Value;
             }
         }
@@ -205,15 +258,20 @@ void Pstream::combineScatter(const List<Pstream::commsStruct>& comms, T& Value)
 
 
 template <class T>
-void Pstream::combineScatter(T& Value)
+void Pstream::combineScatter
+(
+    T& Value,
+    const int tag,
+    const label comm
+)
 {
-    if (Pstream::nProcs() < Pstream::nProcsSimpleSum())
+    if (Pstream::nProcs(comm) < Pstream::nProcsSimpleSum())
     {
-        combineScatter(Pstream::linearCommunication(), Value);
+        combineScatter(Pstream::linearCommunication(comm), Value, tag, comm);
     }
     else
     {
-        combineScatter(Pstream::treeCommunication(), Value);
+        combineScatter(Pstream::treeCommunication(comm), Value, tag, comm);
     }
 }
 
@@ -227,16 +285,18 @@ void Pstream::listCombineGather
 (
     const List<Pstream::commsStruct>& comms,
     List<T>& Values,
-    const CombineOp& cop
+    const CombineOp& cop,
+    const int tag,
+    const label comm
 )
 {
-    if (Pstream::parRun())
+    if (Pstream::nProcs(comm) > 1)
     {
         // Get my communication order
-        const commsStruct& myComm = comms[Pstream::myProcNo()];
+        const commsStruct& myComm = comms[Pstream::myProcNo(comm)];
 
         // Receive from my downstairs neighbours
-        forAll(myComm.below(), belowI)
+        forAll (myComm.below(), belowI)
         {
             label belowID = myComm.below()[belowI];
 
@@ -249,7 +309,9 @@ void Pstream::listCombineGather
                     Pstream::scheduled,
                     belowID,
                     reinterpret_cast<char*>(receivedValues.begin()),
-                    receivedValues.byteSize()
+                    receivedValues.byteSize(),
+                    tag,
+                    comm
                 );
 
                 if (debug > 1)
@@ -258,14 +320,14 @@ void Pstream::listCombineGather
                         << belowID << " data:" << receivedValues << endl;
                 }
 
-                forAll(Values, i)
+                forAll (Values, i)
                 {
                     cop(Values[i], receivedValues[i]);
                 }
             }
             else
             {
-                IPstream fromBelow(Pstream::scheduled, belowID);
+                IPstream fromBelow(Pstream::scheduled, belowID, 0, tag, comm);
                 List<T> receivedValues(fromBelow);
 
                 if (debug > 1)
@@ -274,7 +336,7 @@ void Pstream::listCombineGather
                         << belowID << " data:" << receivedValues << endl;
                 }
 
-                forAll(Values, i)
+                forAll (Values, i)
                 {
                     cop(Values[i], receivedValues[i]);
                 }
@@ -297,12 +359,21 @@ void Pstream::listCombineGather
                     Pstream::scheduled,
                     myComm.above(),
                     reinterpret_cast<const char*>(Values.begin()),
-                    Values.byteSize()
+                    Values.byteSize(),
+                    tag,
+                    comm
                 );
             }
             else
             {
-                OPstream toAbove(Pstream::scheduled, myComm.above());
+                OPstream toAbove
+                (
+                    Pstream::scheduled,
+                    myComm.above(),
+                    0,
+                    tag,
+                    comm
+                );
                 toAbove << Values;
             }
         }
@@ -311,15 +382,35 @@ void Pstream::listCombineGather
 
 
 template <class T, class CombineOp>
-void Pstream::listCombineGather(List<T>& Values, const CombineOp& cop)
+void Pstream::listCombineGather
+(
+    List<T>& Values,
+    const CombineOp& cop,
+    const int tag,
+    const label comm
+)
 {
-    if (Pstream::nProcs() < Pstream::nProcsSimpleSum())
+    if (Pstream::nProcs(comm) < Pstream::nProcsSimpleSum())
     {
-        listCombineGather(Pstream::linearCommunication(), Values, cop);
+        listCombineGather
+        (
+            Pstream::linearCommunication(comm),
+            Values,
+            cop,
+            tag,
+            comm
+        );
     }
     else
     {
-        listCombineGather(Pstream::treeCommunication(), Values, cop);
+        listCombineGather
+        (
+            Pstream::treeCommunication(comm),
+            Values,
+            cop,
+            tag,
+            comm
+        );
     }
 }
 
@@ -328,13 +419,15 @@ template <class T>
 void Pstream::listCombineScatter
 (
     const List<Pstream::commsStruct>& comms,
-    List<T>& Values
+    List<T>& Values,
+    const int tag,
+    const label comm
 )
 {
-    if (Pstream::parRun())
+    if (Pstream::nProcs(comm) > 1)
     {
         // Get my communication order
-        const Pstream::commsStruct& myComm = comms[Pstream::myProcNo()];
+        const Pstream::commsStruct& myComm = comms[Pstream::myProcNo(comm)];
 
         // Receive from up
         if (myComm.above() != -1)
@@ -346,12 +439,21 @@ void Pstream::listCombineScatter
                     Pstream::scheduled,
                     myComm.above(),
                     reinterpret_cast<char*>(Values.begin()),
-                    Values.byteSize()
+                    Values.byteSize(),
+                    tag,
+                    comm
                 );
             }
             else
             {
-                IPstream fromAbove(Pstream::scheduled, myComm.above());
+                IPstream fromAbove
+                (
+                    Pstream::scheduled,
+                    myComm.above(),
+                    0,
+                    tag,
+                    comm
+                );
                 fromAbove >> Values;
             }
 
@@ -362,8 +464,11 @@ void Pstream::listCombineScatter
             }
         }
 
-        // Send to my downstairs neighbours
-        forAll(myComm.below(), belowI)
+        // Send to my downstairs neighbours.  Note reverse order (compared to
+        // receiving). This is to make sure to send to the critical path
+        // (only when using a tree schedule!) first.
+        // This is ESI Comms optimisation, v16.06.  HJ, 19/Sep/2016
+        forAllReverse (myComm.below(), belowI)
         {
             label belowID = myComm.below()[belowI];
 
@@ -379,12 +484,14 @@ void Pstream::listCombineScatter
                     Pstream::scheduled,
                     belowID,
                     reinterpret_cast<const char*>(Values.begin()),
-                    Values.byteSize()
+                    Values.byteSize(),
+                    tag,
+                    comm
                 );
             }
             else
             {
-                OPstream toBelow(Pstream::scheduled, belowID);
+                OPstream toBelow(Pstream::scheduled, belowID, 0, tag, comm);
                 toBelow << Values;
             }
         }
@@ -393,44 +500,60 @@ void Pstream::listCombineScatter
 
 
 template <class T>
-void Pstream::listCombineScatter(List<T>& Values)
+void Pstream::listCombineScatter
+(
+    List<T>& Values,
+    const int tag,
+    const label comm
+)
 {
-    if (Pstream::nProcs() < Pstream::nProcsSimpleSum())
+    if (Pstream::nProcs(comm) < Pstream::nProcsSimpleSum())
     {
-        listCombineScatter(Pstream::linearCommunication(), Values);
+        listCombineScatter
+        (
+            Pstream::linearCommunication(),
+            Values,
+            tag,
+            comm
+        );
     }
     else
     {
-        listCombineScatter(Pstream::treeCommunication(), Values);
+        listCombineScatter
+        (
+            Pstream::treeCommunication(),
+            Values,
+            tag,
+            comm
+        );
     }
 }
 
 
-
-
 // Same thing but for sparse list (map)
 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-
 template <class Container, class CombineOp>
 void Pstream::mapCombineGather
 (
     const List<Pstream::commsStruct>& comms,
     Container& Values,
-    const CombineOp& cop
+    const CombineOp& cop,
+    const int tag,
+    const label comm
 )
 {
-    if (Pstream::parRun())
+    if (Pstream::nProcs(comm) > 1)
     {
         // Get my communication order
-        const commsStruct& myComm = comms[Pstream::myProcNo()];
+        const commsStruct& myComm = comms[Pstream::myProcNo(comm)];
 
         // Receive from my downstairs neighbours
-        forAll(myComm.below(), belowI)
+        forAll (myComm.below(), belowI)
         {
             label belowID = myComm.below()[belowI];
 
-            IPstream fromBelow(Pstream::scheduled, belowID);
+            IPstream fromBelow(Pstream::scheduled, belowID, 0, tag, comm);
             Container receivedValues(fromBelow);
 
             if (debug > 1)
@@ -470,7 +593,7 @@ void Pstream::mapCombineGather
                     << " data:" << Values << endl;
             }
 
-            OPstream toAbove(Pstream::scheduled, myComm.above());
+            OPstream toAbove(Pstream::scheduled, myComm.above(), 0, tag, comm);
             toAbove << Values;
         }
     }
@@ -478,15 +601,35 @@ void Pstream::mapCombineGather
 
 
 template <class Container, class CombineOp>
-void Pstream::mapCombineGather(Container& Values, const CombineOp& cop)
+void Pstream::mapCombineGather
+(
+    Container& Values,
+    const CombineOp& cop,
+    const int tag,
+    const label comm
+)
 {
-    if (Pstream::nProcs() < Pstream::nProcsSimpleSum())
+    if (Pstream::nProcs(comm) < Pstream::nProcsSimpleSum())
     {
-        mapCombineGather(Pstream::linearCommunication(), Values, cop);
+        mapCombineGather
+        (
+            Pstream::linearCommunication(),
+            Values,
+            cop,
+            tag,
+            comm
+        );
     }
     else
     {
-        mapCombineGather(Pstream::treeCommunication(), Values, cop);
+        mapCombineGather
+        (
+            Pstream::treeCommunication(),
+            Values,
+            cop,
+            tag,
+            comm
+        );
     }
 }
 
@@ -495,18 +638,27 @@ template <class Container>
 void Pstream::mapCombineScatter
 (
     const List<Pstream::commsStruct>& comms,
-    Container& Values
+    Container& Values,
+    const int tag,
+    const label comm
 )
 {
-    if (Pstream::parRun())
+    if (Pstream::nProcs(comm) > 1)
     {
         // Get my communication order
-        const Pstream::commsStruct& myComm = comms[Pstream::myProcNo()];
+        const Pstream::commsStruct& myComm = comms[Pstream::myProcNo(comm)];
 
         // Receive from up
         if (myComm.above() != -1)
         {
-            IPstream fromAbove(Pstream::scheduled, myComm.above());
+            IPstream fromAbove
+            (
+                Pstream::scheduled,
+                myComm.above(),
+                0,
+                tag,
+                comm
+            );
             fromAbove >> Values;
 
             if (debug > 1)
@@ -516,8 +668,11 @@ void Pstream::mapCombineScatter
             }
         }
 
-        // Send to my downstairs neighbours
-        forAll(myComm.below(), belowI)
+        // Send to my downstairs neighbours.  Note reverse order (compared to
+        // receiving). This is to make sure to send to the critical path
+        // (only when using a tree schedule!) first.
+        // This is ESI Comms optimisation, v16.06.  HJ, 19/Sep/2016
+        forAllReverse (myComm.below(), belowI)
         {
             label belowID = myComm.below()[belowI];
 
@@ -526,7 +681,7 @@ void Pstream::mapCombineScatter
                 Pout<< " sending to " << belowID << " data:" << Values << endl;
             }
 
-            OPstream toBelow(Pstream::scheduled, belowID);
+            OPstream toBelow(Pstream::scheduled, belowID, 0, tag, comm);
             toBelow << Values;
         }
     }
@@ -534,15 +689,32 @@ void Pstream::mapCombineScatter
 
 
 template <class Container>
-void Pstream::mapCombineScatter(Container& Values)
+void Pstream::mapCombineScatter
+(
+    Container& Values,
+    const int tag,
+    const label comm
+)
 {
-    if (Pstream::nProcs() < Pstream::nProcsSimpleSum())
+    if (Pstream::nProcs(comm) < Pstream::nProcsSimpleSum())
     {
-        mapCombineScatter(Pstream::linearCommunication(), Values);
+        mapCombineScatter
+        (
+            Pstream::linearCommunication(),
+            Values,
+            tag,
+            comm
+        );
     }
     else
     {
-        mapCombineScatter(Pstream::treeCommunication(), Values);
+        mapCombineScatter
+        (
+            Pstream::treeCommunication(),
+            Values,
+            tag,
+            comm
+        );
     }
 }
 

src/foam/db/IOstreams/Pstreams/gatherScatter.C

@@ -45,16 +45,25 @@ void Pstream::gather
 (
     const List<Pstream::commsStruct>& comms,
     T& Value,
-    const BinaryOp& bop
+    const BinaryOp& bop,
+    const int tag,
+    const label comm
 )
 {
-    if (Pstream::parRun())
+    // Return if not active in this comm
+    // HJ, 12/Sep/2016
+    if (Pstream::myProcNo(comm) == -1)
+    {
+        return;
+    }
+
+    if (Pstream::nProcs(comm) > 1)
     {
         // Get my communication order
-        const commsStruct& myComm = comms[Pstream::myProcNo()];
+        const commsStruct& myComm = comms[Pstream::myProcNo(comm)];
 
         // Receive from my downstairs neighbours
-        forAll(myComm.below(), belowI)
+        forAll (myComm.below(), belowI)
         {
             T value;
 
@@ -65,12 +74,21 @@ void Pstream::gather
                     Pstream::scheduled,
                     myComm.below()[belowI],
                     reinterpret_cast<char*>(&value),
-                    sizeof(T)
+                    sizeof(T),
+                    tag,
+                    comm
                 );
             }
             else
             {
-                IPstream fromBelow(Pstream::scheduled, myComm.below()[belowI]);
+                IPstream fromBelow
+                (
+                    Pstream::scheduled,
+                    myComm.below()[belowI],
+                    0,
+                    tag,
+                    comm
+                );
                 fromBelow >> value;
             }
 
@@ -87,12 +105,21 @@ void Pstream::gather
                     Pstream::scheduled,
                     myComm.above(),
                     reinterpret_cast<const char*>(&Value),
-                    sizeof(T)
+                    sizeof(T),
+                    tag,
+                    comm
                 );
             }
             else
             {
-                OPstream toAbove(Pstream::scheduled, myComm.above());
+                OPstream toAbove
+                (
+                    Pstream::scheduled,
+                    myComm.above(),
+                    0,
+                    tag,
+                    comm
+                );
                 toAbove << Value;
             }
         }
@@ -101,26 +128,45 @@ void Pstream::gather
 
 
 template <class T, class BinaryOp>
-void Pstream::gather(T& Value, const BinaryOp& bop)
+void Pstream::gather
+(
+    T& Value,
+    const BinaryOp& bop,
+    const int tag,
+    const label comm
+)
 {
-    if (Pstream::nProcs() < Pstream::nProcsSimpleSum())
+    if (Pstream::nProcs(comm) < Pstream::nProcsSimpleSum())
     {
-        gather(Pstream::linearCommunication(), Value, bop);
+        gather(Pstream::linearCommunication(comm), Value, bop, tag, comm);
     }
     else
     {
-        gather(Pstream::treeCommunication(), Value, bop);
+        gather(Pstream::treeCommunication(comm), Value, bop, tag, comm);
     }
 }
 
 
 template <class T>
-void Pstream::scatter(const List<Pstream::commsStruct>& comms, T& Value)
+void Pstream::scatter
+(
+    const List<Pstream::commsStruct>& comms,
+    T& Value,
+    const int tag,
+    const label comm
+)
 {
-    if (Pstream::parRun())
+    // Return if not active in this comm
+    // HJ, 12/Sep/2016
+    if (Pstream::myProcNo(comm) == -1)
+    {
+        return;
+    }
+
+    if (Pstream::nProcs(comm) > 1)
     {
         // Get my communication order
-        const commsStruct& myComm = comms[Pstream::myProcNo()];
+        const commsStruct& myComm = comms[Pstream::myProcNo(comm)];
 
         // Receive from up
         if (myComm.above() != -1)
@@ -132,18 +178,30 @@ void Pstream::scatter(const List<Pstream::commsStruct>& comms, T& Value)
                     Pstream::scheduled,
                     myComm.above(),
                     reinterpret_cast<char*>(&Value),
-                    sizeof(T)
+                    sizeof(T),
+                    tag,
+                    comm
                 );
             }
             else
             {
-                IPstream fromAbove(Pstream::scheduled, myComm.above());
+                IPstream fromAbove
+                (
+                    Pstream::scheduled,
+                    myComm.above(),
+                    0,
+                    tag,
+                    comm
+                );
                 fromAbove >> Value;
             }
         }
 
-        // Send to my downstairs neighbours
-        forAll(myComm.below(), belowI)
+        // Send to my downstairs neighbours.  Note reverse order (compared to
+        // receiving). This is to make sure to send to the critical path
+        // (only when using a tree schedule!) first.
+        // This is ESI Comms optimisation, v16.06.  HJ, 19/Sep/2016
+        forAllReverse (myComm.below(), belowI)
         {
             if (contiguous<T>())
             {
@@ -152,12 +210,21 @@ void Pstream::scatter(const List<Pstream::commsStruct>& comms, T& Value)
                     Pstream::scheduled,
                     myComm.below()[belowI],
                     reinterpret_cast<const char*>(&Value),
-                    sizeof(T)
+                    sizeof(T),
+                    tag,
+                    comm
                 );
             }
             else
             {
-                OPstream toBelow(Pstream::scheduled,myComm.below()[belowI]);
+                OPstream toBelow
+                (
+                    Pstream::scheduled,
+                    myComm.below()[belowI],
+                    0,
+                    tag,
+                    comm
+                );
                 toBelow << Value;
             }
         }
@@ -166,15 +233,15 @@ void Pstream::scatter(const List<Pstream::commsStruct>& comms, T& Value)
 
 
 template <class T>
-void Pstream::scatter(T& Value)
+void Pstream::scatter(T& Value, const int tag, const label comm)
 {
-    if (Pstream::nProcs() < Pstream::nProcsSimpleSum())
+    if (Pstream::nProcs(comm) < Pstream::nProcsSimpleSum())
     {
-        scatter(Pstream::linearCommunication(), Value);
+        scatter(Pstream::linearCommunication(comm), Value, tag, comm);
     }
     else
     {
-        scatter(Pstream::treeCommunication(), Value);
+        scatter(Pstream::treeCommunication(comm), Value, tag, comm);
     }
 }
 

src/foam/db/IOstreams/Pstreams/gatherScatterList.C

@@ -49,28 +49,42 @@ template <class T>
 void Pstream::gatherList
 (
     const List<Pstream::commsStruct>& comms,
-    List<T>& Values
+    List<T>& Values,
+    const int tag,
+    const label comm
 )
 {
-    if (Pstream::parRun())
+    // Return if not active in this comm
+    // HJ, 12/Sep/2016
+    if (Pstream::myProcNo(comm) == -1)
     {
-        if (Values.size() != Pstream::nProcs())
+        return;
+    }
+
+    if (Pstream::nProcs(comm) > 1)
+    {
+        if (Values.size() != Pstream::nProcs(comm))
         {
             FatalErrorIn
             (
-                "Pstream::gatherList(const List<Pstream::commsStruct>&"
-                ", List<T>)"
+                "void Pstream::gatherList\n"
+                "(\n"
+                "    const List<Pstream::commsStruct>& comms,\n"
+                "    List<T>& Values,\n"
+                "    const int tag,\n"
+                "    const label comm\n"
+                ")"
             )   << "Size of list:" << Values.size()
                 << " does not equal the number of processors:"
-                << Pstream::nProcs()
+                << Pstream::nProcs(comm)
                 << Foam::abort(FatalError);
         }
 
         // Get my communication order
-        const commsStruct& myComm = comms[Pstream::myProcNo()];
+        const commsStruct& myComm = comms[Pstream::myProcNo(comm)];
 
         // Receive from my downstairs neighbours
-        forAll(myComm.below(), belowI)
+        forAll (myComm.below(), belowI)
         {
             label belowID = myComm.below()[belowI];
             const labelList& belowLeaves = comms[belowID].allBelow();
@@ -84,19 +98,21 @@ void Pstream::gatherList
                     Pstream::scheduled,
                     belowID,
                     reinterpret_cast<char*>(receivedValues.begin()),
-                    receivedValues.byteSize()
+                    receivedValues.byteSize(),
+                    tag,
+                    comm
                 );
 
                 Values[belowID] = receivedValues[0];
 
-                forAll(belowLeaves, leafI)
+                forAll (belowLeaves, leafI)
                 {
                     Values[belowLeaves[leafI]] = receivedValues[leafI + 1];
                 }
             }
             else
             {
-                IPstream fromBelow(Pstream::scheduled, belowID);
+                IPstream fromBelow(Pstream::scheduled, belowID, 0, tag, comm);
                 fromBelow >> Values[belowID];
 
                 if (debug > 1)
@@ -107,7 +123,7 @@ void Pstream::gatherList
                 }
 
                 // Receive from all other processors below belowID
-                forAll(belowLeaves, leafI)
+                forAll (belowLeaves, leafI)
                 {
                     label leafID = belowLeaves[leafI];
                     fromBelow >> Values[leafID];
@@ -132,16 +148,16 @@ void Pstream::gatherList
             if (debug > 1)
             {
                 Pout<< " sending to " << myComm.above()
-                    << " data from: " << Pstream::myProcNo()
-                    << " data: " << Values[Pstream::myProcNo()] << endl;
+                    << " data from: " << Pstream::myProcNo(comm)
+                    << " data: " << Values[Pstream::myProcNo(comm)] << endl;
             }
 
             if (contiguous<T>())
             {
                 List<T> sendingValues(belowLeaves.size() + 1);
-                sendingValues[0] = Values[Pstream::myProcNo()];
+                sendingValues[0] = Values[Pstream::myProcNo(comm)];
 
-                forAll(belowLeaves, leafI)
+                forAll (belowLeaves, leafI)
                 {
                     sendingValues[leafI + 1] = Values[belowLeaves[leafI]];
                 }
@@ -151,15 +167,24 @@ void Pstream::gatherList
                     Pstream::scheduled,
                     myComm.above(),
                     reinterpret_cast<const char*>(sendingValues.begin()),
-                    sendingValues.byteSize()
+                    sendingValues.byteSize(),
+                    tag,
+                    comm
                 );
             }
             else
             {
-                OPstream toAbove(Pstream::scheduled, myComm.above());
-                toAbove << Values[Pstream::myProcNo()];
+                OPstream toAbove
+                (
+                    Pstream::scheduled,
+                    myComm.above(),
+                    0,
+                    tag,
+                    comm
+                );
+                toAbove << Values[Pstream::myProcNo(comm)];
 
-                forAll(belowLeaves, leafI)
+                forAll (belowLeaves, leafI)
                 {
                     label leafID = belowLeaves[leafI];
 
@@ -178,15 +203,15 @@ void Pstream::gatherList
 
 
 template <class T>
-void Pstream::gatherList(List<T>& Values)
+void Pstream::gatherList(List<T>& Values, const int tag, const label comm)
 {
-    if (Pstream::nProcs() < Pstream::nProcsSimpleSum())
+    if (Pstream::nProcs(comm) < Pstream::nProcsSimpleSum())
     {
-        gatherList(Pstream::linearCommunication(), Values);
+        gatherList(Pstream::linearCommunication(comm), Values, tag, comm);
     }
     else
     {
-        gatherList(Pstream::treeCommunication(), Values);
+        gatherList(Pstream::treeCommunication(comm), Values, tag, comm);
     }
 }
 
@@ -195,25 +220,39 @@ template <class T>
 void Pstream::scatterList
 (
     const List<Pstream::commsStruct>& comms,
-    List<T>& Values
+    List<T>& Values,
+    const int tag,
+    const label comm
 )
 {
-    if (Pstream::parRun())
+    // Return if not active in this comm
+    // HJ, 12/Sep/2016
+    if (Pstream::myProcNo(comm) == -1)
     {
-        if (Values.size() != Pstream::nProcs())
+        return;
+    }
+
+    if (Pstream::nProcs(comm) > 1)
+    {
+        if (Values.size() != Pstream::nProcs(comm))
         {
             FatalErrorIn
             (
-                "Pstream::scatterList(const List<Pstream::commsStruct>&"
-                ", List<T>)"
+                "void Pstream::scatterList\n"
+                "(\n"
+                "    const List<Pstream::commsStruct>& comms,\n"
+                "    List<T>& Values,\n"
+                "    const int tag,\n"
+                "    const label comm\n"
+                ")"
             )   << "Size of list:" << Values.size()
                 << " does not equal the number of processors:"
-                << Pstream::nProcs()
+                << Pstream::nProcs(comm)
                 << Foam::abort(FatalError);
         }
 
         // Get my communication order
-        const commsStruct& myComm = comms[Pstream::myProcNo()];
+        const commsStruct& myComm = comms[Pstream::myProcNo(comm)];
 
         // Receive from up
         if (myComm.above() != -1)
@@ -229,24 +268,33 @@ void Pstream::scatterList
                     Pstream::scheduled,
                     myComm.above(),
                     reinterpret_cast<char*>(receivedValues.begin()),
-                    receivedValues.byteSize()
+                    receivedValues.byteSize(),
+                    tag,
+                    comm
                 );
 
-                forAll(notBelowLeaves, leafI)
+                forAll (notBelowLeaves, leafI)
                 {
                     Values[notBelowLeaves[leafI]] = receivedValues[leafI];
                 }
             }
             else
             {
-                IPstream fromAbove(Pstream::scheduled, myComm.above());
+                IPstream fromAbove
+                (
+                    Pstream::scheduled,
+                    myComm.above(),
+                    0,
+                    tag,
+                    comm
+                );
 
-                forAll(notBelowLeaves, leafI)
+                forAll (notBelowLeaves, leafI)
                 {
                     label leafID = notBelowLeaves[leafI];
                     fromAbove >> Values[leafID];
 
-                    if (debug)
+                    if (debug > 1)
                     {
                         Pout<< " received through "
                             << myComm.above() << " data for:" << leafID
@@ -256,8 +304,11 @@ void Pstream::scatterList
             }
         }
 
-        // Send to my downstairs neighbours
-        forAll(myComm.below(), belowI)
+        // Send to my downstairs neighbours.  Note reverse order (compared to
+        // receiving). This is to make sure to send to the critical path
+        // (only when using a tree schedule!) first.
+        // This is ESI Comms optimisation, v16.06.  HJ, 19/Sep/2016
+        forAllReverse (myComm.below(), belowI)
         {
             label belowID = myComm.below()[belowI];
             const labelList& notBelowLeaves = comms[belowID].allNotBelow();
@@ -266,7 +317,7 @@ void Pstream::scatterList
             {
                 List<T> sendingValues(notBelowLeaves.size());
 
-                forAll(notBelowLeaves, leafI)
+                forAll (notBelowLeaves, leafI)
                 {
                     sendingValues[leafI] = Values[notBelowLeaves[leafI]];
                 }
@@ -276,20 +327,22 @@ void Pstream::scatterList
                     Pstream::scheduled,
                     belowID,
                     reinterpret_cast<const char*>(sendingValues.begin()),
-                    sendingValues.byteSize()
+                    sendingValues.byteSize(),
+                    tag,
+                    comm
                 );
             }
             else
             {
-                OPstream toBelow(Pstream::scheduled, belowID);
+                OPstream toBelow(Pstream::scheduled, belowID, 0, tag, comm);
 
                 // Send data destined for all other processors below belowID
-                forAll(notBelowLeaves, leafI)
+                forAll (notBelowLeaves, leafI)
                 {
                     label leafID = notBelowLeaves[leafI];
                     toBelow << Values[leafID];
 
-                    if (debug)
+                    if (debug > 1)
                     {
                         Pout<< " sent through "
                             << belowID << " data for:" << leafID
@@ -303,15 +356,15 @@ void Pstream::scatterList
 
 
 template <class T>
-void Pstream::scatterList(List<T>& Values)
+void Pstream::scatterList(List<T>& Values, const int tag, const label comm)
 {
-    if (Pstream::nProcs() < Pstream::nProcsSimpleSum())
+    if (Pstream::nProcs(comm) < Pstream::nProcsSimpleSum())
     {
-        scatterList(Pstream::linearCommunication(), Values);
+        scatterList(Pstream::linearCommunication(comm), Values, tag, comm);
     }
     else
     {
-        scatterList(Pstream::treeCommunication(), Values);
+        scatterList(Pstream::treeCommunication(comm), Values, tag, comm);
     }
 }
 

src/foam/db/regIOobject/regIOobjectRead.C

@@ -28,7 +28,6 @@ License
 #include "objectRegistry.H"
 #include "PstreamReduceOps.H"
 
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::Istream& Foam::regIOobject::readStream()

src/foam/dimensionedTypes/dimensionedType/dimensionedType.C

@@ -132,45 +132,6 @@ dimensioned<Type>::dimensioned
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template <class Type>
-const word& dimensioned<Type>::name() const
-{
-    return name_;
-}
-
-template <class Type>
-word& dimensioned<Type>::name()
-{
-    return name_;
-}
-
-
-template <class Type>
-const dimensionSet& dimensioned<Type>::dimensions() const
-{
-    return dimensions_;
-}
-
-template <class Type>
-dimensionSet& dimensioned<Type>::dimensions()
-{
-    return dimensions_;
-}
-
-
-template <class Type>
-const Type& dimensioned<Type>::value() const
-{
-    return value_;
-}
-
-template <class Type>
-Type& dimensioned<Type>::value()
-{
-    return value_;
-}
-
-
 template <class Type>
 dimensioned<typename dimensioned<Type>::cmptType> dimensioned<Type>::component
 (

src/foam/dimensionedTypes/dimensionedType/dimensionedType.H

@@ -130,22 +130,41 @@ public:
     // Member functions
 
         //- Return const reference to name.
-        const word& name() const;
+        inline const word& name() const
+        {
+            return name_;
+        }
 
         //- Return non-const reference to name.
-        word& name();
+        inline word& name()
+        {
+            return name_;
+        }
+
 
         //- Return const reference to dimensions.
-        const dimensionSet& dimensions() const;
+        const dimensionSet& dimensions() const
+        {
+            return dimensions_;
+        }
 
         //- Return non-const reference to dimensions.
-        dimensionSet& dimensions();
+        inline dimensionSet& dimensions()
+        {
+            return dimensions_;
+        }
 
         //- Return const reference to value.
-        const Type& value() const;
+        inline const Type& value() const
+        {
+            return value_;
+        }
 
         //- Return non-const reference to value.
-        Type& value();
+        inline Type& value()
+        {
+            return value_;
+        }
 
         //- Return a component as a dimensioned<cmptType>
         dimensioned<cmptType> component(const direction) const;

src/foam/interpolations/GGIInterpolation/GGIInterpolationWeights.C

@@ -418,7 +418,21 @@ void GGIInterpolation<MasterPatch, SlavePatch>::calcAddressing() const
                         neighbPointsInUV
                     );
 
-                if (intersectionArea > VSMALL) // Or > areaErrorTol_ ???
+                scalar intersectionTestArea =
+                    Foam::max
+                    (
+                        VSMALL,
+                        areaErrorTol_()*
+                        Foam::max
+                        (
+                            surfaceAreaMasterPointsInUV,
+                            surfaceAreaNeighbPointsInUV
+                        )
+                    );
+
+                // Fix: previously checked for VSMALL.
+                // HJ, 19/Sep2/106
+                if (intersectionArea > intersectionTestArea)
                 {
                     // We compute the GGI weights based on this
                     // intersection area, and on the individual face
@@ -452,17 +466,17 @@ void GGIInterpolation<MasterPatch, SlavePatch>::calcAddressing() const
                 }
                 else
                 {
-                    WarningIn
-                    (
-                        "GGIInterpolation<MasterPatch, SlavePatch>::"
-                        "calcAddressing()"
-                    )   << "polygonIntersection is returning a "
-                        << "zero surface area between " << nl
-                        << "     Master face: " << faceMi
-                        << " and Neighbour face: " << curCMN[neighbI]
-                        << " intersection area = " << intersectionArea << nl
-                        << "Please check the two quick-check algorithms for "
-                        << "GGIInterpolation.  Something is  missing." << endl;
+//                     WarningIn
+//                     (
+//                         "GGIInterpolation<MasterPatch, SlavePatch>::"
+//                         "calcAddressing()"
+//                     )   << "polygonIntersection is returning a "
+//                         << "zero surface area between " << nl
+//                         << "     Master face: " << faceMi
+//                         << " and Neighbour face: " << curCMN[neighbI]
+//                         << " intersection area = " << intersectionArea << nl
+//                         << "Please check the two quick-check algorithms for "
+//                         << "GGIInterpolation.  Something is  missing." << endl;
                 }
             }
         }

src/foam/matrices/blockLduMatrix/BlockAMG/BlockAMGCycle.H

@@ -128,8 +128,7 @@ public:
         BlockAMGCycle(autoPtr<BlockAMGLevel<Type> > levelPtr);
 
 
-    // Destructor
-
+    //- Destructor
     virtual ~BlockAMGCycle();
 
 

src/foam/matrices/blockLduMatrix/BlockAMG/BlockAMGInterfaceFields/GGIBlockAMGInterfaceField/GGIBlockAMGInterfaceField.C

@@ -100,7 +100,6 @@ Foam::GGIBlockAMGInterfaceField<Type>::GGIBlockAMGInterfaceField
     BlockAMGInterfaceField<Type>(AMGCp, fineInterfaceField),
     ggiInterface_(refCast<const ggiAMGInterface>(AMGCp)),
     doTransform_(false),
-    rank_(),
     fieldTransferBuffer_()
 {
     // If the interface based on a patch this must be taken care specially of
@@ -113,7 +112,6 @@ Foam::GGIBlockAMGInterfaceField<Type>::GGIBlockAMGInterfaceField
             );
 
         doTransform_ = p.doTransform();
-        rank_ = p.rank();
     }
     else if (isA<ggiLduInterfaceField>(fineInterfaceField))
     {
@@ -121,7 +119,6 @@ Foam::GGIBlockAMGInterfaceField<Type>::GGIBlockAMGInterfaceField
             refCast<const ggiLduInterfaceField >(fineInterfaceField);
 
         doTransform_ = p.doTransform();
-        rank_ = p.rank();
     }
     else
     {

src/foam/matrices/blockLduMatrix/BlockAMG/BlockAMGInterfaceFields/GGIBlockAMGInterfaceField/GGIBlockAMGInterfaceField.H

@@ -65,9 +65,6 @@ class GGIBlockAMGInterfaceField
         //- Is the transform required
         bool doTransform_;
 
-        //- Rank of component for transformation.
-        int rank_;
-
         //- Field transfer buffer
         mutable Field<Type> fieldTransferBuffer_;
 
@@ -199,11 +196,6 @@ public:
                 return ggiInterface_.reverseT();
             }
 
-            //- Return rank of component for transform
-            virtual int rank() const
-            {
-                return rank_;
-            }
 
         // Transfer buffer access
 
@@ -212,7 +204,6 @@ public:
             {
                 return fieldTransferBuffer_;
             }
-
 };