Merge branch 'master' of ssh://hunt/home/hunt2/OpenFOAM/OpenFOAM-dev

applications/solvers/compressible/rhoSimpleFoam/hEqn.H

@@ -5,7 +5,7 @@
       - fvm::Sp(fvc::div(phi), h)
       - fvm::laplacian(turbulence->alphaEff(), h)
      ==
-        fvc::div(phi/fvc::interpolate(rho)*fvc::interpolate(p, "div(U,p)"))
+        fvc::div(phi/fvc::interpolate(rho), p, "div(U,p)")
       - p*fvc::div(phi/fvc::interpolate(rho))
     );
 

applications/solvers/compressible/rhoSimpleFoam/pEqn.H

-volScalarField AU = UEqn().A();
-U = UEqn().H()/AU;
+rho = thermo->rho();
+
+volScalarField rUA = 1.0/UEqn().A();
+U = rUA*UEqn().H();
 UEqn.clear();
-phi = fvc::interpolate(rho*U) & mesh.Sf();
-bool closedVolume = adjustPhi(phi, U, p);
 
-for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+bool closedVolume = false;
+
+if (transonic)
 {
-    fvScalarMatrix pEqn
+    surfaceScalarField phid
     (
-        fvm::laplacian(rho/AU, p) == fvc::div(phi)
+        "phid",
+        fvc::interpolate(thermo->psi())*(fvc::interpolate(U) & mesh.Sf())
     );
 
-    pEqn.setReference(pRefCell, pRefValue);
-    // retain the residual from the first iteration
-    if (nonOrth == 0)
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        eqnResidual = pEqn.solve().initialResidual();
-        maxResidual = max(eqnResidual, maxResidual);
-    }
-    else
-    {
-        pEqn.solve();
+        fvScalarMatrix pEqn0
+        (
+            fvm::div(phid, p)
+          - fvm::laplacian(rho*rUA, p)
+        );
+        fvScalarMatrix pEqn = pEqn0;
+        pEqn.relax(mesh.relaxationFactor("pEqn"));
+
+        pEqn.setReference(pRefCell, pRefValue);
+
+        // retain the residual from the first iteration
+        if (nonOrth == 0)
+        {
+            eqnResidual = pEqn.solve().initialResidual();
+            maxResidual = max(eqnResidual, maxResidual);
+        }
+        else
+        {
+            pEqn.solve();
+        }
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            phi == pEqn0.flux();
+        }
     }
+}
+else
+{
+    phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
+    //phi = fvc::interpolate(rho*U) & mesh.Sf();
+    closedVolume = adjustPhi(phi, U, p);
 
-    if (nonOrth == nNonOrthCorr)
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        phi -= pEqn.flux();
+        fvScalarMatrix pEqn
+        (
+            fvm::laplacian(rho*rUA, p) == fvc::div(phi)
+        );
+
+        pEqn.setReference(pRefCell, pRefValue);
+
+        // retain the residual from the first iteration
+        if (nonOrth == 0)
+        {
+            eqnResidual = pEqn.solve().initialResidual();
+            maxResidual = max(eqnResidual, maxResidual);
+        }
+        else
+        {
+            pEqn.solve();
+        }
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            phi -= pEqn.flux();
+        }
     }
 }
 
+
 #include "incompressible/continuityErrs.H"
 
 // Explicitly relax pressure for momentum corrector
 p.relax();
 
-U -= fvc::grad(p)/AU;
+rho = thermo->rho();
+rho.relax();
+Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;
+
+U -= rUA*fvc::grad(p);
 U.correctBoundaryConditions();
 
 bound(p, pMin);
@@ -46,7 +98,3 @@ if (closedVolume)
     p += (initialMass - fvc::domainIntegrate(thermo->psi()*p))
         /fvc::domainIntegrate(thermo->psi());
 }
-
-rho = thermo->rho();
-rho.relax();
-Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;