[Onurb]

Hello,

I am modifying the icoFoam solver in the following way:

fvScalarMatrix CEqn 
(
    fvm::ddt(C) + fvm::div(phi,C) - fvm::laplacian (D,C) == fvm::Sp(-CTE_1,C)
);
CEqn.solve();

fvVectorMatrix UEqn
(
    fvm::ddt(U) + fvm::div(phi,U) - fvm::laplacian(nu,U) == fvm::Sp(-CTE_2,U) 
);

if (piso.momentumPredictor())
{
    solve(UEqn == -fvc::grad(p));
}

while (piso.correct())
{
    volScalarField rAU(1.0/UEqn.A());

    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));

    surfaceScalarField phiHbyA
    (
        "phiHbyA", fvc::flux(HbyA) + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
    );

    adjustPhi(phiHbyA, U, p);

    constrainPressure(p, U, phiHbyA, rAU);

    while (piso.correctNonOrthogonal())
    {
        fvScalarMatrix pEqn
        (
            fvm::laplacian(rAU, p) == fvc::div(phiHbyA) - CTE_3*C
        );

        pEqn.setReference(pRefCell, pRefValue);

        pEqn.solve();

        if (piso.finalNonOrthogonalIter())
        {
            phi = phiHbyA + pEqn.flux();
        }
    }
    #include "continuityErrs.H"            

    U = HbyA - rAU*fvc::grad(p);
    U.correctBoundaryConditions();
}

The solver converges until the source term in the continuity equation became nonzero. If I set CTE_3*C = 0, everything works fine. Could someone verify if is there any mistake in the implementation?