[reverseila]

The equations are from here

Here is the code:

// Including necessary class definition and libraries ( Header files )
#include "fvCFD.H"
#include "simpleControl.H"

// handler for printing stuff
#define print(X) Info << X << endl

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
    // Creating objects from the above defined classes
    // These .H files are basically code snippets ( "inclusion" files )
    #include "setRootCaseLists.H"
    #include "createTime.H"
    #include "createMesh.H"

    simpleControl simple(mesh);

    #include "createFields.H"

    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    Info<< "\nCalculating scalar transport\n" << endl;

    // Trigger the objects to store the old value of theirs at each iteration
    c.oldTime();
    p.oldTime();
    pSt.oldTime();
    // pD.oldTime();
    w.oldTime();
    q.oldTime();

    // Define constant coefficients before the loop to only be calculated once
    dimensionedScalar coeff1 = (Pmi/Pma);
    dimensionedScalar coeff2 = 1e-9 * coeff1 * (VT/De);
    dimensionedScalar coeff3 = -2 * F / CapSt;

    while (simple.loop(runTime))
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        #include "CourantNo.H"

        // Potential Equation
        fvScalarMatrix pEqn
            (
             - fvm::laplacian(c, p)
             ==
             coeff2 * fvc::ddt(q)
            );

        pEqn.relax();
        pEqn.solve();

        while (simple.correctNonOrthogonal())
        {
            // Concentration Equation
            fvScalarMatrix cEqn
                (
                 fvm::ddt(c)
                 + coeff_conv * fvm::div(phi, c)
                 - coeff_diff_1 * fvm::laplacian(De, c) // spacer *2 : *0
                 - coeff_diff_2 * fvm::laplacian(De, c) // electrode *0 : *1
                 ==
                 coeff1 * fvc::ddt(w)
                );

            cEqn.relax();
            cEqn.solve();

        }

        q = src_coeff_p * (-c) * exp(Mu) * sinh(pD/VT);
        pSt = coeff3 * q;
        pD = max(-1 * src_coeff_p * p + EV * src_coeff_c - pSt);
        //w = c * exp(Mu) * cosh(pD/VT);
        w = sqrt(q * q + 4 * src_coeff_c * c * c);

        runTime.write();
    }

    Info<< endl;
    Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
        << "  ClockTime = " << runTime.elapsedClockTime() << " s"
        << nl << endl;

    Info<< "End\n" << endl;

    return 0;

I'm not sure why my source terms are zero