diff --git a/applications/solvers/combustion/dieselEngineFoam/spraySummary.H b/applications/solvers/combustion/dieselEngineFoam/spraySummary.H
index 553cb52940..907784b438 100644
--- a/applications/solvers/combustion/dieselEngineFoam/spraySummary.H
+++ b/applications/solvers/combustion/dieselEngineFoam/spraySummary.H
@@ -1,19 +1,20 @@
         label Nparcels = dieselSpray.size();
         reduce(Nparcels, sumOp<label>());
 
-        Info<< "\nNumber of parcels in system | "
+        Info<< "\nNumber of parcels in system.... | "
             << Nparcels << endl
-            << "Injected liquid mass....... | "
+            << "Injected liquid mass........... | "
             << 1e6*dieselSpray.injectedMass(runTime.value()) << " mg" << endl
-            << "Liquid Mass in system...... | "
+            << "Liquid Mass in system.......... | "
             << 1e6*dieselSpray.liquidMass() << " mg" << endl
-            << "SMD, Dmax.................. | "
+            << "SMD, Dmax...................... | "
             << dieselSpray.smd()*1e6 << " mu, "
             << dieselSpray.maxD()*1e6 << " mu"
             << endl;
 
-        scalar evapMass = 
-            dieselSpray.injectedMass(runTime.value()) - dieselSpray.liquidMass();
+        scalar evapMass =
+            dieselSpray.injectedMass(runTime.value())
+          - dieselSpray.liquidMass();
 
         scalar gasMass = fvc::domainIntegrate(rho).value();
 
@@ -26,6 +27,6 @@
 
         scalar addedMass = gasMass - gasMass0;
 
-        Info<< "Added gas mass = " << 1e6*addedMass << " mg" << nl
-            << "Evaporation Continuity Error| "
+        Info<< "Added gas mass................. | " << 1e6*addedMass << " mg"
+            << nl << "Evaporation Continuity Error... | "
             << 1e6*(addedMass - evapMass) << " mg" << endl;
diff --git a/applications/utilities/postProcessing/patch/patchIntegrate/patchIntegrate.C b/applications/utilities/postProcessing/patch/patchIntegrate/patchIntegrate.C
index 69ffc17f91..dc92aade57 100644
--- a/applications/utilities/postProcessing/patch/patchIntegrate/patchIntegrate.C
+++ b/applications/utilities/postProcessing/patch/patchIntegrate/patchIntegrate.C
@@ -31,6 +31,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
+#include "cyclicPolyPatch.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 // Main program:
@@ -75,26 +76,54 @@ int main(int argc, char *argv[])
             }
 
             // Give patch area
-            Info<< "    Patch area = " << gSum(mesh.Sf().boundaryField()[patchi]) << endl;
-
-            if (fieldHeader.headerClassName() == "volScalarField")
+            if (isType<cyclicPolyPatch>(mesh.boundaryMesh()[patchi]))
             {
-                Info<< "    Reading volScalarField " << fieldName << endl;
-                volScalarField field(fieldHeader, mesh);
+                Info<< "    Cyclic patch area: " << nl;
+                label nFaces = mesh.boundaryMesh()[patchi].size();
+                vector sum1 = vector::zero;
+                vector sum2 = vector::zero;
+                for (label i=0; i<nFaces/2; i++)
+                {
+                    sum1 += mesh.Sf().boundaryField()[patchi][i];
+                    sum2 += mesh.Sf().boundaryField()[patchi][i+nFaces/2];
+                }
+                reduce(sum1, sumOp<vector>());
+                reduce(sum2, sumOp<vector>());
+                Info<< "    - half 1 = " << sum1 << ", " << mag(sum1) << nl
+                    << "    - half 2 = " << sum2 << ", " << mag(sum2) << nl
+                    << "    - total  = " << (sum1 + sum2) << ", "
+                    << mag(sum1 + sum2) << endl;;
+            }
+            else
+            {
+                Info<< "    Patch area = "
+                    << gSum(mesh.Sf().boundaryField()[patchi]) << endl;
+            }
 
+            // Read field and calc integral
+            if (fieldHeader.headerClassName() == volScalarField::typeName)
+            {
+                Info<< "    Reading " << volScalarField::typeName << " "
+                    << fieldName << endl;
+
+                volScalarField field(fieldHeader, mesh);
                 vector sumField = gSum
                 (
                     mesh.Sf().boundaryField()[patchi]
-                  * field.boundaryField()[patchi]
+                   *field.boundaryField()[patchi]
                 );
 
                 Info<< "    Integral of " << fieldName << " over patch "
                     << patchName << '[' << patchi << ']' << " = "
                     << sumField << nl;
             }
-            else if (fieldHeader.headerClassName() == "surfaceScalarField")
+            else if
+            (
+                fieldHeader.headerClassName() == surfaceScalarField::typeName
+            )
             {
-                Info<< "    Reading surfaceScalarField " << fieldName << endl;
+                Info<< "    Reading " << surfaceScalarField::typeName << " "
+                    << fieldName << endl;
 
                 surfaceScalarField field(fieldHeader, mesh);
                 scalar sumField = gSum(field.boundaryField()[patchi]);
@@ -106,8 +135,10 @@ int main(int argc, char *argv[])
             else
             {
                 FatalError
-                    << "Only possible to integrate volScalarFields "
-                    << "and surfaceScalarFields" << nl << exit(FatalError);
+                    << "Only possible to integrate "
+                    << volScalarField::typeName << "s "
+                    << "and " << surfaceScalarField::typeName << "s"
+                    << nl << exit(FatalError);
             }
         }
         else
diff --git a/etc/settings.csh b/etc/settings.csh
index ee6a08c94a..73b14eb38c 100644
--- a/etc/settings.csh
+++ b/etc/settings.csh
@@ -121,7 +121,7 @@ unset MPI_ARCH_PATH
 
 switch ("$WM_MPLIB")
 case OPENMPI:
-    set mpi_version=openmpi-1.2.8
+    set mpi_version=openmpi-1.3
     setenv MPI_HOME $WM_THIRD_PARTY_DIR/$mpi_version
     setenv MPI_ARCH_PATH $MPI_HOME/platforms/$WM_OPTIONS
 
diff --git a/etc/settings.sh b/etc/settings.sh
index 2621e1530c..ff56af0b95 100644
--- a/etc/settings.sh
+++ b/etc/settings.sh
@@ -151,7 +151,7 @@ unset MPI_ARCH_PATH
 
 case "$WM_MPLIB" in
 OPENMPI)
-    mpi_version=openmpi-1.2.8
+    mpi_version=openmpi-1.3
     export MPI_HOME=$WM_THIRD_PARTY_DIR/$mpi_version
     export MPI_ARCH_PATH=$MPI_HOME/platforms/$WM_OPTIONS
 
diff --git a/src/finiteVolume/cfdTools/general/bound/bound.C b/src/finiteVolume/cfdTools/general/bound/bound.C
index ccceae463b..b629033754 100644
--- a/src/finiteVolume/cfdTools/general/bound/bound.C
+++ b/src/finiteVolume/cfdTools/general/bound/bound.C
@@ -53,7 +53,6 @@ void Foam::bound(volScalarField& vsf, const dimensionedScalar& vsf0)
             vsf0.value()
         );
 
-        vsf.correctBoundaryConditions();
         vsf.boundaryField() = max(vsf.boundaryField(), vsf0.value());
     }
 }
diff --git a/src/finiteVolume/cfdTools/general/include/setDeltaT.H b/src/finiteVolume/cfdTools/general/include/setDeltaT.H
index 5b75787589..58cd25f233 100644
--- a/src/finiteVolume/cfdTools/general/include/setDeltaT.H
+++ b/src/finiteVolume/cfdTools/general/include/setDeltaT.H
@@ -27,7 +27,7 @@ Global
 
 Description
     Reset the timestep to maintain a constant maximum courant Number.
-    Reduction of time-step is imediate but increase is damped to avoid
+    Reduction of time-step is immediate, but increase is damped to avoid
     unstable oscillations.
 
 \*---------------------------------------------------------------------------*/
@@ -45,7 +45,7 @@ if (adjustTimeStep)
             maxDeltaT
         )
     );
-    
+
     Info<< "deltaT = " <<  runTime.deltaT().value() << endl;
 }
 
diff --git a/src/finiteVolume/cfdTools/general/include/setInitialDeltaT.H b/src/finiteVolume/cfdTools/general/include/setInitialDeltaT.H
index 142371758e..c3d9bdad3b 100644
--- a/src/finiteVolume/cfdTools/general/include/setInitialDeltaT.H
+++ b/src/finiteVolume/cfdTools/general/include/setInitialDeltaT.H
@@ -33,7 +33,7 @@ Description
 
 if (adjustTimeStep)
 {
-    if (CoNum > SMALL)
+    if ((runTime.timeIndex() == 0) && (CoNum > SMALL))
     {
         runTime.setDeltaT
         (
diff --git a/src/lagrangian/dieselSpray/injector/unitInjector/unitInjector.C b/src/lagrangian/dieselSpray/injector/unitInjector/unitInjector.C
index caceb9d7d2..6b1b07dc7a 100644
--- a/src/lagrangian/dieselSpray/injector/unitInjector/unitInjector.C
+++ b/src/lagrangian/dieselSpray/injector/unitInjector/unitInjector.C
@@ -30,19 +30,19 @@ License
 #include "mathematicalConstants.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
 namespace Foam
 {
+    defineTypeNameAndDebug(unitInjector, 0);
 
-defineTypeNameAndDebug(unitInjector, 0);
-
-addToRunTimeSelectionTable
-(
-    injectorType,
-    unitInjector,
-    dictionary
-);
+    addToRunTimeSelectionTable
+    (
+        injectorType,
+        unitInjector,
+        dictionary
+    );
 }
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -74,26 +74,36 @@ Foam::unitInjector::unitInjector
     // check if time entries for soi and eoi match
     if (mag(massFlowRateProfile_[0][0]-TProfile_[0][0]) > SMALL)
     {
-        FatalError << "unitInjector::unitInjector(const time& t, const dictionary dict) " << endl
-            << " start-times do not match for TemperatureProfile and massFlowRateProfile."
-            << abort(FatalError);
+        FatalErrorIn
+        (
+            "unitInjector::unitInjector(const time& t, const dictionary dict)"
+        )<< "start-times do not match for TemperatureProfile and "
+         << " massFlowRateProfile." << nl << exit (FatalError);
     }
 
-    if (mag(massFlowRateProfile_[massFlowRateProfile_.size()-1][0]-TProfile_[TProfile_.size()-1][0]) > SMALL)
+    if
+    (
+        mag(massFlowRateProfile_[massFlowRateProfile_.size()-1][0]
+      - TProfile_[TProfile_.size()-1][0])
+      > SMALL
+    )
     {
-        FatalError << "unitInjector::unitInjector(const time& t, const dictionary dict) " << endl
-            << " end-times do not match for TemperatureProfile and massFlowRateProfile."
-            << abort(FatalError);
+        FatalErrorIn
+        (
+            "unitInjector::unitInjector(const time& t, const dictionary dict)"
+        )<< "end-times do not match for TemperatureProfile and "
+         << "massFlowRateProfile." << nl << exit(FatalError);
     }
 
     // convert CA to real time
     forAll(massFlowRateProfile_, i)
     {
-        massFlowRateProfile_[i][0] = t.userTimeToTime(massFlowRateProfile_[i][0]);
+        massFlowRateProfile_[i][0] =
+            t.userTimeToTime(massFlowRateProfile_[i][0]);
         velocityProfile_[i][0] = massFlowRateProfile_[i][0];
         injectionPressureProfile_[i][0] = massFlowRateProfile_[i][0];
     }
-    
+
     forAll(TProfile_, i)
     {
         TProfile_[i][0] = t.userTimeToTime(TProfile_[i][0]);
@@ -105,14 +115,14 @@ Foam::unitInjector::unitInjector
     {
         // correct the massFlowRateProfile to match the injected mass
         massFlowRateProfile_[i][1] *= mass_/integratedMFR;
-        
+
         CdProfile_[i][0] = massFlowRateProfile_[i][0];
         CdProfile_[i][1] = Cd_;
     }
 
     // Normalize the direction vector
     direction_ /= mag(direction_);
-    
+
     setTangentialVectors();
 
     // check molar fractions
@@ -124,9 +134,9 @@ Foam::unitInjector::unitInjector
 
     if (mag(Xsum - 1.0) > SMALL)
     {
-        Info << "Warning!!!\n unitInjector::unitInjector(const time& t, Istream& is)"
-            << "X does not add up to 1.0, correcting molar fractions."
-            << endl;
+        WarningIn("unitInjector::unitInjector(const time& t, Istream& is)")
+            << "X does not sum to 1.0, correcting molar fractions."
+            << nl << endl;
         forAll(X_, i)
         {
             X_[i] /= Xsum;
@@ -169,18 +179,18 @@ Foam::label Foam::unitInjector::nParcelsToInject
     const scalar time1
 ) const
 {
-
     scalar mInj = mass_*(fractionOfInjection(time1)-fractionOfInjection(time0));
     label nParcels = label(mInj/averageParcelMass_ + 0.49);
-    
     return nParcels;
 }
 
+
 const Foam::vector Foam::unitInjector::position(const label n) const
 {
     return position_;
 }
 
+
 Foam::vector Foam::unitInjector::position
 (
     const label n,
@@ -212,7 +222,7 @@ Foam::vector Foam::unitInjector::position
         scalar iAngle = 2.0*mathematicalConstant::pi*rndGen.scalar01();
 
         return
-        ( 
+        (
             position_
           + iRadius
           * (
@@ -220,22 +230,25 @@ Foam::vector Foam::unitInjector::position
             + tangentialInjectionVector2_*sin(iAngle)
           )
         );
-        
+
     }
 
     return position_;
 }
 
+
 Foam::label Foam::unitInjector::nHoles() const
 {
     return 1;
 }
 
+
 Foam::scalar Foam::unitInjector::d() const
 {
     return d_;
 }
 
+
 const Foam::vector& Foam::unitInjector::direction
 (
     const label i,
@@ -245,6 +258,7 @@ const Foam::vector& Foam::unitInjector::direction
     return direction_;
 }
 
+
 Foam::scalar Foam::unitInjector::mass
 (
     const scalar time0,
@@ -255,7 +269,7 @@ Foam::scalar Foam::unitInjector::mass
 {
     scalar mInj = mass_*(fractionOfInjection(time1)-fractionOfInjection(time0));
 
-    // correct mass if calculation is 2D 
+    // correct mass if calculation is 2D
     if (twoD)
     {
         mInj *= 0.5*angleOfWedge/mathematicalConstant::pi;
@@ -264,82 +278,80 @@ Foam::scalar Foam::unitInjector::mass
     return mInj;
 }
 
+
 Foam::scalar Foam::unitInjector::mass() const
 {
     return mass_;
 }
 
+
 const Foam::scalarField& Foam::unitInjector::X() const
 {
     return X_;
 }
 
+
 Foam::List<Foam::unitInjector::pair> Foam::unitInjector::T() const
 {
     return TProfile_;
 }
 
+
 Foam::scalar Foam::unitInjector::T(const scalar time) const
 {
     return getTableValue(TProfile_, time);
 }
 
+
 Foam::scalar Foam::unitInjector::tsoi() const
 {
     return massFlowRateProfile_[0][0];
 }
 
+
 Foam::scalar Foam::unitInjector::teoi() const
 {
     return massFlowRateProfile_[massFlowRateProfile_.size()-1][0];
 }
 
-Foam::scalar Foam::unitInjector::massFlowRate
-(
-    const scalar time
-) const
+
+Foam::scalar Foam::unitInjector::massFlowRate(const scalar time) const
 {
     return getTableValue(massFlowRateProfile_, time);
 }
 
-Foam::scalar Foam::unitInjector::injectionPressure
-(
-    const scalar time
-) const
+
+Foam::scalar Foam::unitInjector::injectionPressure(const scalar time) const
 {
     return getTableValue(injectionPressureProfile_, time);
 }
 
-Foam::scalar Foam::unitInjector::velocity
-(
-    const scalar time
-) const
+
+Foam::scalar Foam::unitInjector::velocity(const scalar time) const
 {
     return getTableValue(velocityProfile_, time);
 }
 
+
 Foam::List<Foam::unitInjector::pair> Foam::unitInjector::CdProfile() const
 {
     return CdProfile_;
 }
 
-Foam::scalar Foam::unitInjector::Cd
-(
-    const scalar time
-) const
+
+Foam::scalar Foam::unitInjector::Cd(const scalar time) const
 {
     return Cd_;
 }
 
+
 Foam::scalar Foam::unitInjector::fractionOfInjection(const scalar time) const
 {
     return integrateTable(massFlowRateProfile_, time)/mass_;
 }
 
-Foam::scalar Foam::unitInjector::injectedMass
-(
-    const scalar t
-) const
+
+Foam::scalar Foam::unitInjector::injectedMass(const scalar t) const
 {
     return mass_*fractionOfInjection(t);
 }
@@ -351,7 +363,6 @@ void Foam::unitInjector::correctProfiles
     const scalar referencePressure
 )
 {
-
     scalar A = 0.25*mathematicalConstant::pi*pow(d_, 2.0);
     scalar pDummy = 1.0e+5;
 
@@ -365,14 +376,17 @@ void Foam::unitInjector::correctProfiles
     }
 }
 
-Foam::vector Foam::unitInjector::tan1(const label n) const
+
+Foam::vector Foam::unitInjector::tan1(const label) const
 {
     return tangentialInjectionVector1_;
 }
 
-Foam::vector Foam::unitInjector::tan2(const label n) const
+
+Foam::vector Foam::unitInjector::tan2(const label) const
 {
     return tangentialInjectionVector2_;
 }
 
+
 // ************************************************************************* //
diff --git a/src/lagrangian/dieselSpray/spray/sprayInject.C b/src/lagrangian/dieselSpray/spray/sprayInject.C
index 104debda87..394f361c68 100644
--- a/src/lagrangian/dieselSpray/spray/sprayInject.C
+++ b/src/lagrangian/dieselSpray/spray/sprayInject.C
@@ -63,14 +63,10 @@ void spray::inject()
         {
             Np = max(1, Np);
             scalar mp = mass/Np/nHoles;
-        
-            // constT is only larger than zero for the first 
+
+            // constT is only larger than zero for the first
             // part of the injection
-            scalar constT = max
-            (
-                0.0,
-                it->tsoi() - time0
-            );
+            scalar constT = max(0.0, it->tsoi() - time0);
 
             // deltaT is the duration of injection during this timestep
             scalar deltaT = min
@@ -103,9 +99,10 @@ void spray::inject()
                         axisOfWedgeNormal_,
                         rndGen_
                     );
-                
+
                     scalar diameter = injection().d0(i, toi);
-                    vector direction = injection().direction(i, n, toi, diameter);
+                    vector direction =
+                        injection().direction(i, n, toi, diameter);
                     vector U = injection().velocity(i, toi)*direction;
 
                     scalar symComponent = direction & axisOfSymmetry_;
@@ -117,13 +114,13 @@ void spray::inject()
                     scalar ddev = breakup().yDot0();
 
                     label injectorCell = mesh_.findCell(injectionPosition);
-                
+
 #                   include "findInjectorCell.H"
-   
+
                     if (injectorCell >= 0)
                     {
                         scalar liquidCore = 1.0;
-                
+
                         // construct the parcel that is to be injected
 
                         parcel* pPtr = new parcel
@@ -148,19 +145,16 @@ void spray::inject()
                             fuels_->components()
                         );
 
-                        injectedLiquidKE_ += 0.5*pPtr->m()*pow(mag(U), 2.0);
-                    
+                        injectedLiquidKE_ += 0.5*pPtr->m()*magSqr(U);
+
                         scalar dt = time - toi;
 
                         pPtr->stepFraction() =
                             (runTime_.deltaT().value() - dt)
                            /runTime_.deltaT().value();
 
-                        bool keepParcel = pPtr->move
-                        (
-                            *this
-                        );
-  
+                        bool keepParcel = pPtr->move(*this);
+
                         if (keepParcel)
                         {
                             addParticle(pPtr);
diff --git a/src/lagrangian/dieselSpray/spray/sprayOps.C b/src/lagrangian/dieselSpray/spray/sprayOps.C
index 04694a4ddd..72e44a6bbe 100644
--- a/src/lagrangian/dieselSpray/spray/sprayOps.C
+++ b/src/lagrangian/dieselSpray/spray/sprayOps.C
@@ -65,7 +65,7 @@ void spray::evolve()
     inject();
     atomizationLoop();
     breakupLoop();
-            
+
     UInterpolator_.clear();
     rhoInterpolator_.clear();
     pInterpolator_.clear();
@@ -89,12 +89,7 @@ void spray::move()
 
 void spray::breakupLoop()
 {
-    for
-    (
-        spray::iterator elmnt = begin();
-        elmnt != end();
-        ++elmnt
-    )
+    forAllIter(spray::iterator, *this, elmnt)
     {
         // interpolate...
         vector velocity = UInterpolator().interpolate
@@ -128,12 +123,7 @@ void spray::breakupLoop()
 
 void spray::atomizationLoop()
 {
-    for
-    (
-        spray::iterator elmnt = begin();
-        elmnt != end();
-        ++elmnt
-    )
+    forAllIter(spray::iterator, *this, elmnt)
     {
         // interpolate...
         vector velocity = UInterpolator().interpolate
diff --git a/src/lagrangian/dieselSpray/spraySubModels/breakupModel/reitzKHRT/reitzKHRT.C b/src/lagrangian/dieselSpray/spraySubModels/breakupModel/reitzKHRT/reitzKHRT.C
index 050db9f356..85de5e5e97 100644
--- a/src/lagrangian/dieselSpray/spraySubModels/breakupModel/reitzKHRT/reitzKHRT.C
+++ b/src/lagrangian/dieselSpray/spraySubModels/breakupModel/reitzKHRT/reitzKHRT.C
@@ -108,31 +108,19 @@ void reitzKHRT::breakupParcel
 
     // frequency of the fastest growing KH-wave
     scalar omegaKH =
-    (
-        0.34 + 0.38*pow(weGas, 1.5)
-    )/
-    (
-        (1 + ohnesorge)*(1 + 1.4*pow(taylor, 0.6))
-    )*sqrt(sigma/(rhoLiquid*pow(r, 3)) );
+        (0.34 + 0.38*pow(weGas, 1.5))
+       /((1 + ohnesorge)*(1 + 1.4*pow(taylor, 0.6)))
+       *sqrt(sigma/(rhoLiquid*pow(r, 3)));
 
-    // ... and the corresponding KH wave-length.
+    // corresponding KH wave-length.
     scalar lambdaKH =
-    9.02*r*
-    (
-        1.0 + 0.45*sqrt(ohnesorge)
-    )*
-    (
-        1.0 + 0.4*pow(taylor, 0.7)
-    )/
-    pow
-    (
-        (
-            1.0 + 0.865*pow(weGas, 1.67)
-        ),
-        0.6
-    );
+        9.02
+       *r
+       *(1.0 + 0.45*sqrt(ohnesorge))
+       *(1.0 + 0.4*pow(taylor, 0.7))
+       /pow(1.0 + 0.865*pow(weGas, 1.67), 0.6);
 
-    // the characteristic Kelvin-Helmholtz breakup time
+    // characteristic Kelvin-Helmholtz breakup time
     scalar tauKH = 3.726*b1_*r/(omegaKH*lambdaKH);
 
     // stable KH diameter
@@ -140,7 +128,6 @@ void reitzKHRT::breakupParcel
 
     // the frequency of the fastest growing RT wavelength.
     scalar helpVariable = mag(gt*(rhoLiquid - rhoGas));
-
     scalar omegaRT = sqrt
     (
         2.0*pow(helpVariable, 1.5)
@@ -148,12 +135,9 @@ void reitzKHRT::breakupParcel
     );
 
     // RT wave number
-    scalar KRT = sqrt
-    (
-        helpVariable/(3.0*sigma + VSMALL)
-    );
+    scalar KRT = sqrt(helpVariable/(3.0*sigma + VSMALL));
 
-    // the wavelength of the fastest growing Raleigh-Taylor frequency
+    // wavelength of the fastest growing RT frequency
     scalar lambdaRT = 2.0*mathematicalConstant::pi*cRT_/(KRT + VSMALL);
 
     // if lambdaRT < diameter, then RT waves are growing on the surface
@@ -163,7 +147,7 @@ void reitzKHRT::breakupParcel
         p.ct() += deltaT;
     }
 
-    // the characteristic RT breakup time
+    // characteristic RT breakup time
     scalar tauRT = cTau_/(omegaRT + VSMALL);
 
     // check if we have RT breakup
@@ -178,7 +162,7 @@ void reitzKHRT::breakupParcel
     // otherwise check for KH breakup
     else if (dc < p.d())
     {
-        // no breakup below weber = 12
+        // no breakup below Weber = 12
         if (weGas > weberLimit_)
         {
 
@@ -188,25 +172,24 @@ void reitzKHRT::breakupParcel
             // reduce the diameter according to the rate-equation
             p.d() = (fraction*dc + p.d())/(1.0 + fraction);
 
-            scalar dc3 = pow(dc, 3.0);
-            scalar ms = rhoLiquid*Np*dc3*mathematicalConstant::pi/6.0;
+            scalar ms = rhoLiquid*Np*pow3(dc)*mathematicalConstant::pi/6.0;
             p.ms() += ms;
 
-            label nParcels = spray_.injectors()[injector].properties()->nParcelsToInject
-            (
-                spray_.injectors()[injector].properties()->tsoi(),
-                spray_.injectors()[injector].properties()->teoi()
-            );
+            // Total number of parcels for the whole injection event
+            label nParcels =
+                spray_.injectors()[injector].properties()->nParcelsToInject
+                (
+                    spray_.injectors()[injector].properties()->tsoi(),
+                    spray_.injectors()[injector].properties()->teoi()
+                );
 
-            scalar averageParcelMass = spray_.injectors()[injector].properties()->mass()/nParcels;
+            scalar averageParcelMass =
+                spray_.injectors()[injector].properties()->mass()/nParcels;
 
-            if
-            (
-                (p.ms()/averageParcelMass > msLimit_)
-            )
+            if (p.ms()/averageParcelMass > msLimit_)
             {
                 // set the initial ms value to -GREAT. This prevents
-                // new droplets from being formed from the childDroplet
+                // new droplets from being formed from the child droplet
                 // from the KH instability
 
                 // mass of stripped child parcel
diff --git a/src/postProcessing/functionObjects/field/fieldMinMax/fieldMinMax.C b/src/postProcessing/functionObjects/field/fieldMinMax/fieldMinMax.C
index 172a1325c8..051f9f9a6d 100644
--- a/src/postProcessing/functionObjects/field/fieldMinMax/fieldMinMax.C
+++ b/src/postProcessing/functionObjects/field/fieldMinMax/fieldMinMax.C
@@ -38,6 +38,18 @@ namespace Foam
 }
 
 
+template<>
+const char* Foam::NamedEnum<Foam::fieldMinMax::modeType, 2>::names[] =
+{
+    "magnitude",
+    "component"
+};
+
+
+const Foam::NamedEnum<Foam::fieldMinMax::modeType, 2>
+Foam::fieldMinMax::modeTypeNames_;
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::fieldMinMax::fieldMinMax
@@ -52,6 +64,7 @@ Foam::fieldMinMax::fieldMinMax
     obr_(obr),
     active_(true),
     log_(false),
+    mode_(mdMag),
     fieldSet_(),
     fieldMinMaxFilePtr_(NULL)
 {
@@ -85,6 +98,7 @@ void Foam::fieldMinMax::read(const dictionary& dict)
     {
         log_ = dict.lookupOrDefault<Switch>("log", false);
 
+        mode_ = modeTypeNames_[dict.lookup("mode")];
         dict.lookup("fields") >> fieldSet_;
     }
 }
@@ -176,16 +190,13 @@ void Foam::fieldMinMax::calcMinMaxFields<Foam::scalar>
 {
     if (obr_.foundObject<volScalarField>(fieldName))
     {
-        const volScalarField& field =
-            obr_.lookupObject<volScalarField>(fieldName);
-        scalar minValue = min(field).value();
-        scalar maxValue = max(field).value();
-
-        reduce(minValue, minOp<scalar>());
-        reduce(maxValue, maxOp<scalar>());
-
         if (Pstream::master())
         {
+            const volScalarField& field =
+                obr_.lookupObject<volScalarField>(fieldName);
+            scalar minValue = min(field).value();
+            scalar maxValue = max(field).value();
+
             fieldMinMaxFilePtr_() << obr_.time().value() << tab
                 << fieldName << tab << minValue << tab << maxValue << endl;
 
diff --git a/src/postProcessing/functionObjects/field/fieldMinMax/fieldMinMax.H b/src/postProcessing/functionObjects/field/fieldMinMax/fieldMinMax.H
index 1f2ac449c6..17cfcf54e3 100644
--- a/src/postProcessing/functionObjects/field/fieldMinMax/fieldMinMax.H
+++ b/src/postProcessing/functionObjects/field/fieldMinMax/fieldMinMax.H
@@ -49,6 +49,7 @@ SourceFiles
 #include "OFstream.H"
 #include "Switch.H"
 #include "pointFieldFwd.H"
+#include "NamedEnum.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -66,9 +67,18 @@ class mapPolyMesh;
 
 class fieldMinMax
 {
+public:
+
+    enum modeType
+    {
+        mdMag,
+        mdCmpt
+    };
+
+
 protected:
 
-    // Private data
+    // Protected data
 
         //- Name of this set of forces,
         //  Also used as the name of the probes directory.
@@ -82,11 +92,17 @@ protected:
         //- Switch to send output to Info as well as to file
         Switch log_;
 
-        //- Patches to integrate forces over
+        //- Mode type names
+        static const NamedEnum<modeType, 2> modeTypeNames_;
+
+        //- Mode for min/max - only applicable for ranks > 0
+        modeType mode_;
+
+        //- Fields to assess min/max
         wordList fieldSet_;
 
 
-        //- Forces/moment file ptr
+        //- Min/max file ptr
         autoPtr<OFstream> fieldMinMaxFilePtr_;
 
 
diff --git a/src/postProcessing/functionObjects/field/fieldMinMax/fieldMinMaxTemplates.C b/src/postProcessing/functionObjects/field/fieldMinMax/fieldMinMaxTemplates.C
index 354668d3f4..06d6bd585e 100644
--- a/src/postProcessing/functionObjects/field/fieldMinMax/fieldMinMaxTemplates.C
+++ b/src/postProcessing/functionObjects/field/fieldMinMax/fieldMinMaxTemplates.C
@@ -38,24 +38,60 @@ void Foam::fieldMinMax::calcMinMaxFields(const word& fieldName)
 
     if (obr_.foundObject<fieldType>(fieldName))
     {
-        const fieldType& field = obr_.lookupObject<fieldType>(fieldName);
-        scalar minValue = min(mag(field)).value();
-        scalar maxValue = max(mag(field)).value();
-
-        reduce(minValue, minOp<scalar>());
-        reduce(maxValue, maxOp<scalar>());
-
         if (Pstream::master())
         {
-            fieldMinMaxFilePtr_() << obr_.time().value() << tab
-                << fieldName << tab << minValue << tab << maxValue << endl;
-
-            if (log_)
+            const fieldType& field = obr_.lookupObject<fieldType>(fieldName);
+            switch (mode_)
             {
-                Info<< "fieldMinMax output:" << nl
-                    << "    min(mag(" << fieldName << ")) = " << minValue << nl
-                    << "    max(mag(" << fieldName << ")) = " << maxValue << nl
-                    << endl;
+                case mdMag:
+                {
+                    scalar minValue = min(mag(field)).value();
+                    scalar maxValue = max(mag(field)).value();
+
+                    fieldMinMaxFilePtr_() << obr_.time().value() << tab
+                        << fieldName << tab << minValue << tab << maxValue
+                        << endl;
+
+                    if (log_)
+                    {
+                        Info<< "fieldMinMax output:" << nl
+                            << "    min(mag(" << fieldName << ")) = "
+                            << minValue << nl
+                            << "    max(mag(" << fieldName << ")) = "
+                            << maxValue << nl
+                            << endl;
+                    }
+                    break;
+                }
+                case mdCmpt:
+                {
+                    Type minValue = min(field).value();
+                    Type maxValue = max(field).value();
+
+                    fieldMinMaxFilePtr_() << obr_.time().value() << tab
+                        << fieldName << tab << minValue << tab << maxValue
+                        << endl;
+
+                    if (log_)
+                    {
+                        Info<< "fieldMinMax output:" << nl
+                            << "    cmptMin(" << fieldName << ") = "
+                            << minValue << nl
+                            << "    cmptMax(" << fieldName << ") = "
+                            << maxValue << nl
+                            << endl;
+                    }
+                    break;
+                }
+                default:
+                {
+                    FatalErrorIn
+                    (
+                        "Foam::fieldMinMax::calcMinMaxFields"
+                        "(const word& fieldName)"
+                    )<< "Unknown min/max mode: " << modeTypeNames_[mode_]
+                     << exit(FatalError);
+                }
             }
         }
     }
diff --git a/tutorials/dieselFoam/aachenBomb/0/alphat b/tutorials/dieselFoam/aachenBomb/0/alphat
new file mode 100644
index 0000000000..1681c54f00
--- /dev/null
+++ b/tutorials/dieselFoam/aachenBomb/0/alphat
@@ -0,0 +1,32 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       volScalarField;
+    location    "0";
+    object      alphat;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    walls
+    {
+        type            alphatWallFunction;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/dieselFoam/aachenBomb/0/epsilon b/tutorials/dieselFoam/aachenBomb/0/epsilon
index de3d91a3ed..c14d3d7e02 100644
--- a/tutorials/dieselFoam/aachenBomb/0/epsilon
+++ b/tutorials/dieselFoam/aachenBomb/0/epsilon
@@ -1,40 +1,32 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.5                                   |
-|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
 FoamFile
 {
     version     2.0;
-    format      ascii;
+    format      binary;
     class       volScalarField;
+    location    "0";
     object      epsilon;
 }
-// ************************************************************************* //
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-dimensions          [ 0 2 -3 0 0 0 0 ];
+dimensions      [0 2 -3 0 0 0 0];
 
-internalField       uniform 90.0;
+internalField   uniform 90;
 
 boundaryField
 {
     walls
     {
-        type                zeroGradient;
+        type            epsilonWallFunction;
+        value           uniform 90;
     }
-
-    front
-    {
-        type                wedge;
-    }
-
-    back
-    {
-        type                wedge;
-    }
-
 }
 
+
 // ************************************************************************* //
diff --git a/tutorials/dieselFoam/aachenBomb/0/k b/tutorials/dieselFoam/aachenBomb/0/k
index 1ca9a40d69..e6e58aecfa 100644
--- a/tutorials/dieselFoam/aachenBomb/0/k
+++ b/tutorials/dieselFoam/aachenBomb/0/k
@@ -1,40 +1,32 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.5                                   |
-|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
 FoamFile
 {
     version     2.0;
-    format      ascii;
+    format      binary;
     class       volScalarField;
+    location    "0";
     object      k;
 }
-// ************************************************************************* //
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-dimensions          [ 0 2 -2 0 0 0 0 ];
+dimensions      [0 2 -2 0 0 0 0];
 
-internalField       uniform 1.0;
+internalField   uniform 1;
 
 boundaryField
 {
     walls
     {
-        type                zeroGradient;
+        type            kQRWallFunction;
+        value           uniform 1;
     }
-
-    front
-    {
-        type                wedge;
-    }
-
-    back
-    {
-        type                wedge;
-    }
-
 }
 
+
 // ************************************************************************* //
diff --git a/tutorials/dieselFoam/aachenBomb/0/mut b/tutorials/dieselFoam/aachenBomb/0/mut
new file mode 100644
index 0000000000..2750ce32ec
--- /dev/null
+++ b/tutorials/dieselFoam/aachenBomb/0/mut
@@ -0,0 +1,32 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       volScalarField;
+    location    "0";
+    object      mut;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    walls
+    {
+        type            mutWallFunction;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/dieselFoam/aachenBomb/constant/polyMesh/boundary b/tutorials/dieselFoam/aachenBomb/constant/polyMesh/boundary
new file mode 100644
index 0000000000..5e28555623
--- /dev/null
+++ b/tutorials/dieselFoam/aachenBomb/constant/polyMesh/boundary
@@ -0,0 +1,28 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       polyBoundaryMesh;
+    location    "constant/polyMesh";
+    object      boundary;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+1
+(
+    walls
+    {
+        type            wall;
+        nFaces          19762;
+        startFace       494419;
+    }
+)
+
+// ************************************************************************* //
diff --git a/tutorials/dieselFoam/aachenBomb/constant/turbulenceProperties b/tutorials/dieselFoam/aachenBomb/constant/turbulenceProperties
new file mode 100644
index 0000000000..3d7edb19ec
--- /dev/null
+++ b/tutorials/dieselFoam/aachenBomb/constant/turbulenceProperties
@@ -0,0 +1,19 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.5                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType      RASModel;
+
+// ************************************************************************* //
diff --git a/tutorials/dieselFoam/aachenBomb/system/fvSolution b/tutorials/dieselFoam/aachenBomb/system/fvSolution
index 4f4a45349f..3183443d6e 100644
--- a/tutorials/dieselFoam/aachenBomb/system/fvSolution
+++ b/tutorials/dieselFoam/aachenBomb/system/fvSolution
@@ -16,44 +16,51 @@ FoamFile
 
 solvers
 {
-    rho PCG
+    rho
     {
+        solver           PCG;
         preconditioner   DIC;
         tolerance        1e-06;
         relTol           0;
     };
-    U PBiCG
+    U
     {
+        solver           PBiCG;
         preconditioner   DILU;
         tolerance        1e-06;
         relTol           0;
     };
-    p PCG
+    p
     {
+        solver           PCG;
         preconditioner   DIC;
         tolerance        1e-09;
         relTol           0;
     };
-    Yi PBiCG
+    Yi
     {
+        solver           PBiCG;
         preconditioner   DILU;
         tolerance        1e-06;
         relTol           0;
     };
-    h PBiCG
+    h
     {
+        solver           PBiCG;
         preconditioner   DILU;
         tolerance        1e-06;
         relTol           0;
     };
-    k PBiCG
+    k
     {
+        solver           PBiCG;
         preconditioner   DILU;
         tolerance        1e-06;
         relTol           0;
     };
-    epsilon PBiCG
+    epsilon
     {
+        solver           PBiCG;
         preconditioner   DILU;
         tolerance        1e-06;
         relTol           0;