{
    fvScalarMatrix testSource(testEtaMesh, testEtaMesh.dimensions()*dimVol/dimTime);

    testSource.source() = 1.0;
    
    fvScalarMatrix QiEqn
    (
        fvm::ddt(testEtaMesh)
//      + mvConvection->fvmDiv(phi, he)
//      + fvc::ddt(rho, K) + fvc::div(phi, K)
//      + (
//            he.name() == "ea"
//          ? fvc::div
//            (
//                fvc::absolute(phi/fvc::interpolate(rho), U),
//                p,
//                "div(phiv,p)"
//            )
//          : -dpdt
//        )
      - fvm::laplacian(DQi, testEtaMesh)
      ==
      testSource
    );

    QiEqn.relax();

    QiEqn.solve();
}
