surfaceScalarField Ef
        (
            IOobject
            (
                "Ef",
                mesh.time().timeName(),
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            fvc::interpolate(j0Re) & mesh.Sf() / (fvc::interpolate(sigma))
        );

        Info<< nl << "Calculate applied current distribution" << endl;
        
        // TODO read 'correctNonOrthogonal' from the coil subDict
        while (control.correctNonOrthogonal())
        {
            V0Re.storePrevIter();
            fvScalarMatrix V0ReEqn
            (
                fvm::laplacian // check sign
                (
                    dimensionedScalar("1", dimless, 1),
                    V0Re
                )
            ==
               fvc::div(Ef)
            );

            // TODO ADD REFERENCE LOOKUP
            V0ReEqn.setReference(0, 0.0);
            V0ReEqn.solve(control.subDict(Region::COIL, V0Re.name()));

            if (control.finalNonOrthogonalIter())
            {
                Ef -= V0ReEqn.flux();
            }
            else
            {
                V0Re.relax();
            }
        }

I have been solving the potentialFlow equation with little bit different variables, However what happends in the solver is that when I run it once it goes 1 times through the loop after which it tries to solve and it gets a floatingPoint error, when I run it again everything works and it goes through the loops and converges. - I have looked at removing

fv::div(Ef)

from the loop interely which still results in the issues, which I think narrows it down to V0Re, I check the initilization and nothing changed. So I am bit clueless on what is going wrong.