driftFluxFoam: separate advection and diffusion for predictor-corrector MULES to reduce operator-splitting error

applications/solvers/multiphase/driftFluxFoam/alphaEqn.H

@@ -19,16 +19,7 @@
             ).fvmDiv(phi, alpha1)
         );
 
-        solve
-        (
-            alpha1Eqn
-          - fv::gaussLaplacianScheme<scalar, scalar>
-            (
-                mesh,
-                linear<scalar>(mesh),
-                fv::uncorrectedSnGrad<scalar>(mesh)
-            ).fvmLaplacian(fvc::interpolate(turbulence->nut()), alpha1)
-        );
+        solve(alpha1Eqn);
 
         Info<< "Phase-1 volume fraction = "
             << alpha1.weightedAverage(mesh.Vsc()).value()
@@ -42,6 +33,7 @@
         if (alphaApplyPrevCorr && tphiAlphaCorr0.valid())
         {
             Info<< "Applying the previous iteration correction flux" << endl;
+
             #ifdef LTSSOLVE
             MULES::LTScorrect
             (

applications/solvers/multiphase/driftFluxFoam/alphaEqnSubCycle.H

@@ -50,7 +50,6 @@
 
     // Apply the diffusion term separately to allow implicit solution
     // and boundedness of the explicit advection
-    if (!MULESCorr)
     {
         fvScalarMatrix alpha1Eqn
         (