diff --git a/applications/solvers/incompressible/channelFoam/channelFoam.C b/applications/solvers/incompressible/channelFoam/channelFoam.C
index 740370230..2bf2783b2 100644
--- a/applications/solvers/incompressible/channelFoam/channelFoam.C
+++ b/applications/solvers/incompressible/channelFoam/channelFoam.C
@@ -26,7 +26,8 @@ Application
 
 Description
     Incompressible LES solver for flow in a channel.
-    Consistent formulation without time-step and relaxation dependence by Jasak
+    Consistent formulation without time-step and relaxation dependence by
+    Jasak and Tukovic.
 
 Author
     Hrvoje Jasak, Wikki Ltd.  All rights reserved
@@ -66,6 +67,9 @@ int main(int argc, char *argv[])
 
         sgsModel->correct();
 
+        // Time derivative matrix
+        fvVectorMatrix ddtUEqn(fvm::ddt(U));
+
         // Convection-diffusion matrix
         fvVectorMatrix HUEqn
         (
@@ -75,24 +79,23 @@ int main(int argc, char *argv[])
             flowDirection*gradP
         );
 
-        // Time derivative matrix
-        fvVectorMatrix ddtUEqn(fvm::ddt(U));
-
         if (piso.momentumPredictor())
         {
             solve(ddtUEqn + HUEqn == -fvc::grad(p));
         }
 
-        // Prepare clean Ap without time derivative contribution
-        // HJ, 26/Oct/2015
-        volScalarField aU = HUEqn.A();
+        // Prepare clean 1/a_p without time derivative contribution
+        volScalarField rAU = 1.0/HUEqn.A();
 
         // --- PISO loop
 
         while (piso.correct())
         {
-            U = HUEqn.H()/aU;
-            phi = (fvc::interpolate(U) & mesh.Sf());
+            // Calculate U from convection-diffusion matrix
+            U = rAU*HUEqn.H();
+
+            // Consistently calculate flux
+            piso.calcTransientConsistentFlux(phi, U, rAU, ddtUEqn);
 
             adjustPhi(phi, U, p);
 
@@ -100,7 +103,14 @@ int main(int argc, char *argv[])
             {
                 fvScalarMatrix pEqn
                 (
-                    fvm::laplacian(1/aU, p) == fvc::div(phi)
+                    fvm::laplacian
+                    (
+                        fvc::interpolate(rAU)/piso.aCoeff(),
+                        p,
+                        "laplacian(rAU," + p.name() + ')'
+                    )
+                 ==
+                    fvc::div(phi)
                 );
 
                 pEqn.setReference(pRefCell, pRefValue);
@@ -117,13 +127,8 @@ int main(int argc, char *argv[])
 
 #           include "continuityErrs.H"
 
-            // Note: cannot call H(U) here because the velocity is not complete
-            // HJ, 22/Jan/2016
-            U = 1.0/(aU + ddtUEqn.A())*
-                (
-                    U*aU - fvc::grad(p) + ddtUEqn.H()
-                );
-            U.correctBoundaryConditions();
+            // Consistently reconstruct velocity after pressure equation
+            piso.reconstructTransientVelocity(U, ddtUEqn, rAU, p, phi);
         }
 
         // Correct driving force for a constant mass flow rate
@@ -135,9 +140,9 @@ int main(int argc, char *argv[])
         // Calculate the pressure gradient increment needed to
         // adjust the average flow-rate to the correct value
         dimensionedScalar gragPplus =
-            (magUbar - magUbarStar)*aU.weightedAverage(mesh.V());
+            (magUbar - magUbarStar)/rAU.weightedAverage(mesh.V());
 
-        U += gragPplus/aU*flowDirection;
+        U += gragPplus*rAU*flowDirection;
 
         gradP += gragPplus;