diff --git a/applications/solvers/compressible/sonicDyMFoam/createFields.H b/applications/solvers/compressible/sonicDyMFoam/createFields.H
index fac70b85e..6d3406829 100644
--- a/applications/solvers/compressible/sonicDyMFoam/createFields.H
+++ b/applications/solvers/compressible/sonicDyMFoam/createFields.H
@@ -6,6 +6,7 @@
     );
     basicPsiThermo& thermo = pThermo();
 
+    volScalarField& T = const_cast<volScalarField&>(thermo.T());
     volScalarField& p = thermo.p();
     volScalarField& e = thermo.e();
     const volScalarField& psi = thermo.psi();
diff --git a/applications/solvers/compressible/sonicDyMFoam/pEqn.H b/applications/solvers/compressible/sonicDyMFoam/pEqn.H
new file mode 100644
index 000000000..4be700755
--- /dev/null
+++ b/applications/solvers/compressible/sonicDyMFoam/pEqn.H
@@ -0,0 +1,69 @@
+{
+    U = UEqn.H()/UEqn.A();
+
+#   include "limitU.H"
+
+    surfaceScalarField phid
+    (
+        "phid",
+        fvc::interpolate(psi)*
+        (
+            (fvc::interpolate(U) & mesh.Sf())
+          - fvc::meshPhi(rho, U)
+        )
+    );
+
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        // Store pressure for under-relaxation
+        p.storePrevIter();
+
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi, p)
+          + fvm::div(phid, p)
+          - fvm::laplacian(rho/UEqn.A(), p)
+        );
+
+        if
+        (
+//              oCorr == nOuterCorr - 1
+            corr == nCorr - 1
+         && nonOrth == nNonOrthCorr
+        )
+        {
+            pEqn.solve
+            (
+                mesh.solutionDict().solver(p.name() + "Final")
+            );
+        }
+        else
+        {
+            pEqn.solve(mesh.solutionDict().solver(p.name()));
+        }
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            phi = pEqn.flux();
+        }
+
+
+        // Bound the pressure
+        if (min(p) < pMin || max(p) > pMax)
+        {
+            p.max(pMin);
+            p.min(pMax);
+            p.correctBoundaryConditions();
+        }
+
+        // Relax the pressure
+        p.relax();
+    }
+
+#   include "compressibleContinuityErrs.H"
+
+    U -= fvc::grad(p)/UEqn.A();
+    U.correctBoundaryConditions();
+
+#   include "limitU.H"
+}
diff --git a/applications/solvers/compressible/sonicDyMFoam/sonicDyMFoam.C b/applications/solvers/compressible/sonicDyMFoam/sonicDyMFoam.C
index f14dd3367..377f63da7 100644
--- a/applications/solvers/compressible/sonicDyMFoam/sonicDyMFoam.C
+++ b/applications/solvers/compressible/sonicDyMFoam/sonicDyMFoam.C
@@ -77,10 +77,19 @@ int main(int argc, char *argv[])
         }
 
         // Mesh motion update
-//         if (correctPhi && meshChanged)
-//         {
+        if (meshChanged)
+        {
+            T.correctBoundaryConditions();
+            p.correctBoundaryConditions();
+            e.correctBoundaryConditions();
+            thermo.correct();
+            rho = thermo.rho();
+
+            if (correctPhi)
+            {
 // #           include "correctPhi.H"
-//         }
+            }
+        }
 
         if (meshChanged)
         {
@@ -103,73 +112,7 @@ int main(int argc, char *argv[])
             // --- PISO loop
             for (int corr = 0; corr < nCorr; corr++)
             {
-                U = UEqn.H()/UEqn.A();
-
-#               include "limitU.H"
-
-                surfaceScalarField phid
-                (
-                    "phid",
-                    fvc::interpolate(psi)*
-                    (
-                        (fvc::interpolate(U) & mesh.Sf())
-                      - fvc::meshPhi(rho, U)
-                    )
-                );
-
-                for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
-                {
-                    // Store pressure for under-relaxation
-                    p.storePrevIter();
-
-                    fvScalarMatrix pEqn
-                    (
-                        fvm::ddt(psi, p)
-                      + fvm::div(phid, p)
-                      - fvm::laplacian(rho/UEqn.A(), p)
-                    );
-
-                    if
-                    (
-//                         oCorr == nOuterCorr - 1
-                        corr == nCorr - 1
-                     && nonOrth == nNonOrthCorr
-                    )
-                    {
-                        pEqn.solve
-                        (
-                            mesh.solutionDict().solver(p.name() + "Final")
-                        );
-                    }
-                    else
-                    {
-                        pEqn.solve(mesh.solutionDict().solver(p.name()));
-                    }
-
-                    if (nonOrth == nNonOrthCorr)
-                    {
-                        phi = pEqn.flux();
-                    }
-
-
-                    // Bound the pressure
-                    if (min(p) < pMin || max(p) > pMax)
-                    {
-                        p.max(pMin);
-                        p.min(pMax);
-                        p.correctBoundaryConditions();
-                    }
-
-                    // Relax the pressure
-                    p.relax();
-                }
-
-#               include "compressibleContinuityErrs.H"
-
-                U -= fvc::grad(p)/UEqn.A();
-                U.correctBoundaryConditions();
-
-#               include "limitU.H"
+#               include "pEqn.H"
             }
 
             // Recalculate density