diff --git a/applications/solvers/combustion/PDRFoam/XiEqns b/applications/solvers/combustion/PDRFoam/XiEqns
deleted file mode 100644
index 8de3bffe4b1c689d3ed8ebc4e914a8448add523f..0000000000000000000000000000000000000000
--- a/applications/solvers/combustion/PDRFoam/XiEqns
+++ /dev/null
@@ -1,91 +0,0 @@
-    // Calculate Xi according to the selected flame wrinkling model
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    // Calculate coefficients for Gulder's flame speed correlation
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    volScalarField up = uPrimeCoef*sqrt((2.0/3.0)*k);
-    volScalarField epsilonPlus = pow(uPrimeCoef, 3)*epsilon;
-
-    volScalarField tauEta = sqrt(mag(thermo->muu()/(rhou*epsilonPlus)));
-    volScalarField Reta = up/
-    (
-        sqrt(epsilonPlus*tauEta)
-      + dimensionedScalar("1e-8", up.dimensions(), 1e-8)
-    );
-
-    else if (XiModel == "algebraic")
-    {
-        // Simple algebraic model for Xi based on Gulders correlation
-        // with a linear correction function to give a plausible profile for Xi
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-        volScalarField GEta = GEtaCoef/tauEta;
-        volScalarField XiEqEta = 1.0 + XiCoef*sqrt(up/(Su + SuMin))*Reta;
-
-        volScalarField R =
-            GEta*XiEqEta/(XiEqEta - 0.999) + GIn*XiIn/(XiIn - 0.999);
-
-        volScalarField XiEqStar = R/(R - GEta - GIn);
-
-        //- Calculate the unweighted b
-        //volScalarField Rrho = thermo->rhou()/thermo->rhob();
-        //volScalarField bbar = b/(b + (Rrho*(1.0 - b)));
-
-        Xi == 1.0 + (1.0 + (2*XiShapeCoef)*(0.5 - b))*(XiEqStar - 1.0);
-    }
-    else if (XiModel == "transport")
-    {
-        // Calculate Xi transport coefficients based on Gulders correlation
-        // and DNS data for the rate of generation
-        // with a linear correction function to give a plausible profile for Xi
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-        volScalarField GEta = GEtaCoef/tauEta;
-        volScalarField XiEqEta = 1.0 + XiCoef*sqrt(up/(Su + SuMin))*Reta;
-
-        volScalarField R =
-            GEta*XiEqEta/(XiEqEta - 0.999) + GIn*XiIn/(XiIn - 0.999);
-
-        volScalarField XiEqStar = R/(R - GEta - GIn);
-
-        volScalarField XiEq =
-            1.0 + (1.0 + (2*XiShapeCoef)*(0.5 - b))*(XiEqStar - 1.0);
-
-        volScalarField G = R*(XiEq - 1.0)/XiEq;
-
-    // Calculate Xi flux
-    // ~~~~~~~~~~~~~~~~~
-    surfaceScalarField phiXi =
-        phiSt
-      + (
-          - fvc::interpolate(fvc::laplacian(Dbf, b)/mgb)*nf
-          + fvc::interpolate(rho)*fvc::interpolate(Su*(1.0/Xi - Xi))*nf
-        );
-
-
-        // Solve for the flame wrinkling
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-        solve
-        (
-            betav*fvm::ddt(rho, Xi)
-          + mvConvection->fvmDiv(phi, Xi)
-          + fvm::div(phiXi, Xi, "div(phiXi,Xi)")
-          - fvm::Sp(fvc::div(phiXi), Xi)
-         ==
-            betav*rho*R
-          - fvm::Sp(betav*rho*(R - G), Xi)
-        );
-
-
-        // Correct boundedness of Xi
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~
-        Xi.max(1.0);
-        Xi = min(Xi, 2.0*XiEq);
-        Info<< "max(Xi) = " << max(Xi).value() << endl;
-        Info<< "max(XiEq) = " << max(XiEq).value() << endl;
-    }
-    else
-    {
-        FatalError
-            << args.executable() << " : Unknown Xi model " << XiModel
-            << abort(FatalError);
-    }