{
    p.boundaryField().updateCoeffs();

    volScalarField A
    (
  	IOobject
  	(
    	 "A",
    	 runTime.timeName(),
    	 mesh,
    	 IOobject::NO_READ,
    	 IOobject::NO_WRITE
  	),
	UrelEqn().A()
    );
    A.write();

    volVectorField H
    (
  	IOobject
  	(
    	 "H",
    	 runTime.timeName(),
    	 mesh,
    	 IOobject::NO_READ,
    	 IOobject::NO_WRITE
  	),
	UrelEqn().H()
    );
    H.write();

    volScalarField H1
    (
  	IOobject
  	(
    	 "H1",
    	 runTime.timeName(),
    	 mesh,
    	 IOobject::NO_READ,
    	 IOobject::NO_WRITE
  	),
	UrelEqn().H1()
    );
    H1.write();

    volScalarField rAUrel(1.0/UrelEqn().A());
    Urel = rAUrel*UrelEqn().H();
    UrelEqn.clear();

    phi = fvc::interpolate(Urel, "interpolate(HbyA)") & mesh.Sf();
    adjustPhi(phi, Urel, p);

    // Non-orthogonal pressure corrector loop
    while (simple.correctNonOrthogonal())
    {
        fvScalarMatrix pEqn
        (
            fvm::laplacian(rAUrel, p) == fvc::div(phi)
        );

        pEqn.setReference(pRefCell, pRefValue);

        pEqn.solve();

        if (simple.finalNonOrthogonalIter())
        {
            phi -= pEqn.flux();
        }
    }

    #include "continuityErrs.H"

    // Explicitly relax pressure for momentum corrector
    p.relax();

    // Momentum corrector
    Urel -= rAUrel*fvc::grad(p);
    Urel.correctBoundaryConditions();
    sources.correct(Urel);
}