Info<< "Reading velocity field U\n" << endl;
volVectorField U
(
    IOobject
    (
        "U",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);

// Initialise the velocity internal field to zero
U = dimensionedVector("0", U.dimensions(), Zero);

surfaceScalarField phi
(
    IOobject
    (
        "phi",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    fvc::flux(U)
);

if (args.optionFound("initialiseUBCs"))
{
    U.correctBoundaryConditions();
    phi = fvc::flux(U);
}


// 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
args.optionReadIfPresent("pName", pName);

// Infer the pressure BCs from the velocity BCs
wordList pBCTypes
(
    U.boundaryField().size(),
    fixedValueFvPatchScalarField::typeName
);

forAll(U.boundaryField(), patchi)
{
    if (U.boundaryField()[patchi].fixesValue())
    {
        pBCTypes[patchi] = zeroGradientFvPatchScalarField::typeName;
    }
}

// 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
(
    IOobject
    (
        pName,
        runTime.timeName(),
        mesh,
        IOobject::READ_IF_PRESENT,
        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;
autoPtr<volScalarField> PhiPtr;

IOobject io
(
    "Phi",
    runTime.timeName(),
    mesh,
    IOobject::MUST_READ,
    IOobject::NO_WRITE
);

if (io.typeHeaderOk<volScalarField>())
{
    PhiPtr.reset(new volScalarField(io, mesh));
}
else
{
    // Cannot just use p.boundaryField().types() since does not initialise
    // complex boundary types. Instead re-clone them from p.

    io.readOpt() = IOobject::NO_READ;
    PhiPtr.reset
    (
        new volScalarField
        (
            io,
            mesh,
            dimensionedScalar("Phi", dimLength*dimVelocity, 0),
            p.boundaryField().types()
        )
    );

    const volScalarField::GeometricBoundaryField& bp = p.boundaryField();
    volScalarField::GeometricBoundaryField& bPhi = PhiPtr().boundaryField();

    forAll(bp, patchI)
    {
        bPhi.set(patchI, bp[patchI].clone(PhiPtr().dimensionedInternalField()));
    }
}
volScalarField& Phi = PhiPtr();

label PhiRefCell = 0;
scalar PhiRefValue = 0;
setRefCell
(
    Phi,
    potentialFlow.dict(),
    PhiRefCell,
    PhiRefValue
);
mesh.setFluxRequired(Phi.name());