//mixture fraction variance equation
//cf) mfVar should not be ZERO

Info<<"Calculating epsilon/k"<<endl;
volScalarField epsk(turbulence->epsilon() / turbulence->k());

Info<<"Calculating mixture fraction gradient"<<endl;
volVectorField Gradmf = fvc::grad(mf);

fvScalarMatrix mfVarEqn
(
    fvm::div(phi, mfVar)
  - fvm::laplacian((1/Sc) * turbulence->mut(), mfVar)
 ==
    2 * ((1/Sc) * turbulence->mut()) * (Gradmf & Gradmf)
  - fvm::Sp(2 * rho * epsk, mfVar)
);

mfVarEqn.relax();

mfVarEqn.solve();

mfVar.writeMinMax(Info);

mfVar.max(0.0);