diff --git a/applications/solvers/basic/potentialFoam/createControls.H b/applications/solvers/basic/potentialFoam/createControls.H
deleted file mode 100644
index 7015cae02e8c7400d4c59a211637015c558d6b63..0000000000000000000000000000000000000000
--- a/applications/solvers/basic/potentialFoam/createControls.H
+++ /dev/null
@@ -1,9 +0,0 @@
-const dictionary& potentialFlow
-(
-    mesh.solutionDict().subDict("potentialFlow")
-);
-
-const int nNonOrthCorr
-(
-    potentialFlow.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0)
-);
diff --git a/applications/solvers/basic/potentialFoam/createFields.H b/applications/solvers/basic/potentialFoam/createFields.H
index ce912510ecbf992542d9ba8d3206ec4bb199d6f0..36c4dcc08afe5a408899504a0e3168ca0842d8ed 100644
--- a/applications/solvers/basic/potentialFoam/createFields.H
+++ b/applications/solvers/basic/potentialFoam/createFields.H
@@ -12,6 +12,7 @@ volVectorField U
     mesh
 );
 
+// Initialise the velocity internal field to zero
 U = dimensionedVector("0", U.dimensions(), vector::zero);
 
 surfaceScalarField phi
@@ -34,7 +35,9 @@ if (args.optionFound("initialiseUBCs"))
 }
 
 
-// Default name for the pressure field
+// Construct a pressure field
+// If it is available read it otherwise construct from the velocity BCs
+// converting fixed-value BCs to zero-gradient and vice versa.
 word pName("p");
 
 // Update name of the pressure field from the command-line option
@@ -55,6 +58,9 @@ forAll(U.boundaryField(), patchi)
     }
 }
 
+// Note that registerObject is false for the pressure field. The pressure
+// field in this solver doesn't have a physical value during the solution.
+// It shouldn't be looked up and used by sub models or boundary conditions.
 Info<< "Constructing pressure field " << pName << nl << endl;
 volScalarField p
 (
@@ -64,13 +70,28 @@ volScalarField p
         runTime.timeName(),
         mesh,
         IOobject::READ_IF_PRESENT,
-        IOobject::NO_WRITE
+        IOobject::NO_WRITE,
+        false
     ),
     mesh,
     dimensionedScalar(pName, sqr(dimVelocity), 0),
     pBCTypes
 );
 
+label pRefCell = 0;
+scalar pRefValue = 0.0;
+if (args.optionFound("writep"))
+{
+    setRefCell
+    (
+        p,
+        potentialFlow.dict(),
+        pRefCell,
+        pRefValue
+    );
+}
+
+
 Info<< "Constructing velocity potential field Phi\n" << endl;
 volScalarField Phi
 (
diff --git a/applications/solvers/basic/potentialFoam/potentialFoam.C b/applications/solvers/basic/potentialFoam/potentialFoam.C
index 5c56133788200eed506481d321c1bf5fd55f84fd..d6dc4b78a55711201c7767804259d4d132e0073d 100644
--- a/applications/solvers/basic/potentialFoam/potentialFoam.C
+++ b/applications/solvers/basic/potentialFoam/potentialFoam.C
@@ -24,13 +24,64 @@ License
 Application
     potentialFoam
 
-Description
-    Potential flow solver which solves for the velocity potential
-    from which the flux-field is obtained and velocity field by reconstructing
-    the flux.
+Group
+    grpBasicSolvers
 
-    This application is particularly useful to generate starting fields for
-    Navier-Stokes codes.
+Description
+    Potential flow solver.
+
+    \heading Solver details
+    The potential flow solution is typically employed to generate initial fields
+    for full Navier-Stokes codes.  The flow is evolved using the equation:
+
+    \f[
+        \laplacian \Phi = \div(\vec{U})
+    \f]
+
+    Where:
+    \vartable
+        \Phi      | Velocity potential [m2/s]
+        \vec{U}   | Velocity [m/s]
+    \endvartable
+
+    The corresponding pressure field could be calculated from the divergence
+    of the Euler equation:
+
+    \f[
+        \laplacian p + \div(\div(\vec{U}\otimes\vec{U})) = 0
+    \f]
+
+    but this generates excessive pressure variation in regions of large
+    velocity gradient normal to the flow direction.  A better option is to
+    calculate the pressure field corresponding to velocity variation along the
+    stream-lines:
+
+    \f[
+        \laplacian p + \div(\vec{F}\cdot\div(\vec{U}\otimes\vec{U})) = 0
+    \f]
+    where the flow direction tensor \f$\vec{F}\f$ is obtained from
+    \f[
+        \vec{F} = \hat{\vec{U}}\otimes\hat{\vec{U}}
+    \f]
+
+    \heading Required fields
+    \plaintable
+        U         | Velocity [m/s]
+    \endplaintable
+
+    \heading Optional fields
+    \plaintable
+        p         | Kinematic pressure [m2/s2]
+        Phi       | Velocity potential [m2/s]
+                  | Generated from p (if present) or U if not present
+    \endplaintable
+
+    \heading Options
+    \plaintable
+        -writep   | write the Euler pressure
+        -writePhi | Write the final velocity potential
+        -initialiseUBCs | Update the velocity boundaries before solving for Phi
+    \endplaintable
 
 \*---------------------------------------------------------------------------*/
 
@@ -58,19 +109,13 @@ int main(int argc, char *argv[])
     argList::addBoolOption
     (
         "writePhi",
-        "Write the velocity potential field"
+        "Write the final velocity potential field"
     );
 
     argList::addBoolOption
     (
         "writep",
-        "Calculate and write the pressure field"
-    );
-
-    argList::addBoolOption
-    (
-        "withFunctionObjects",
-        "execute functionObjects"
+        "Calculate and write the Euler pressure field"
     );
 
     #include "setRootCase.H"
@@ -137,7 +182,7 @@ int main(int argc, char *argv[])
         Phi.write();
     }
 
-    // Calculate the pressure field
+    // Calculate the pressure field from the Euler equation
     if (args.optionFound("writep"))
     {
         Info<< nl << "Calculating approximate pressure field" << endl;