diff --git a/src/fvAgglomerationMethods/pairPatchAgglomeration/pairPatchAgglomeration.C b/src/fvAgglomerationMethods/pairPatchAgglomeration/pairPatchAgglomeration.C
index 9a7462127c..9467906f15 100644
--- a/src/fvAgglomerationMethods/pairPatchAgglomeration/pairPatchAgglomeration.C
+++ b/src/fvAgglomerationMethods/pairPatchAgglomeration/pairPatchAgglomeration.C
@@ -44,9 +44,7 @@ bool Foam::pairPatchAgglomeration::continueAgglomerating
 {
     // Check the need for further agglomeration on all processors
     label localnCoarseFaces = nCoarseFaces;
-//     reduce(localnCoarseFaces, sumOp<label>());
     bool contAgg = localnCoarseFaces >= nFacesInCoarsestLevel_;
-    //reduce(contAgg, andOp<bool>());
     return contAgg;
 }
 
diff --git a/src/thermophysicalModels/radiationModels/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C b/src/thermophysicalModels/radiationModels/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C
index 5817279076..52b3192c6d 100644
--- a/src/thermophysicalModels/radiationModels/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C
+++ b/src/thermophysicalModels/radiationModels/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -189,51 +189,40 @@ updateCoeffs()
     }
 
     scalarField& Iw = *this;
-    const vectorField n(patch().Sf()/patch().magSf());
+    const vectorField n(patch().nf());
 
     radiativeIntensityRay& ray =
         const_cast<radiativeIntensityRay&>(dom.IRay(rayId));
 
-    ray.Qr().boundaryField()[patchI] += Iw*(n & ray.dAve());
+    const scalarField nAve(n & ray.dAve());
 
-    scalarField temissivity = emissivity();
+    ray.Qr().boundaryField()[patchI] += Iw*nAve;
+
+    const scalarField temissivity = emissivity();
+
+    scalarField& Qem = ray.Qem().boundaryField()[patchI];
+    scalarField& Qin = ray.Qin().boundaryField()[patchI];
+
+    const vector& myRayId = dom.IRay(rayId).d();
+
+    const scalarField& Ir = dom.Qin();
 
     forAll(Iw, faceI)
     {
-        scalar Ir = 0.0;
-
-        for (label rayI=0; rayI < dom.nRay(); rayI++)
-        {
-            const vector& d = dom.IRay(rayI).d();
-
-            const scalarField& IFace =
-                dom.IRay(rayI).ILambda(lambdaId).boundaryField()[patchI];
-
-            if ((-n[faceI] & d) < 0.0)
-            {
-                // q into the wall
-                const vector& dAve = dom.IRay(rayI).dAve();
-                Ir += IFace[faceI]*mag(n[faceI] & dAve);
-            }
-        }
-
-        const vector& d = dom.IRay(rayId).d();
-
-        if ((-n[faceI] & d) > 0.0)
+        if ((-n[faceI] & myRayId) > 0.0)
         {
             // direction out of the wall
             refGrad()[faceI] = 0.0;
             valueFraction()[faceI] = 1.0;
             refValue()[faceI] =
                 (
-                    Ir*(scalar(1.0) - temissivity[faceI])
+                    Ir[faceI]*(scalar(1.0) - temissivity[faceI])
                   + temissivity[faceI]*physicoChemical::sigma.value()
                   * pow4(Tp[faceI])
                 )/pi;
 
-             // Emmited heat flux from this ray direction
-            ray.Qem().boundaryField()[patchI][faceI] =
-                refValue()[faceI]*(n[faceI] & ray.dAve());
+            // Emmited heat flux from this ray direction
+            Qem[faceI] = refValue()[faceI]*nAve[faceI];
         }
         else
         {
@@ -242,9 +231,8 @@ updateCoeffs()
             refGrad()[faceI] = 0.0;
             refValue()[faceI] = 0.0; //not used
 
-             // Incident heat flux on this ray direction
-            ray.Qin().boundaryField()[patchI][faceI] =
-                Iw[faceI]*(n[faceI] & ray.dAve());
+            // Incident heat flux on this ray direction
+            Qin[faceI] = Iw[faceI]*nAve[faceI];
         }
     }
 
diff --git a/src/thermophysicalModels/radiationModels/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.H b/src/thermophysicalModels/radiationModels/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.H
index f1a38077e0..8f8b1aef17 100644
--- a/src/thermophysicalModels/radiationModels/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.H
+++ b/src/thermophysicalModels/radiationModels/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/thermophysicalModels/radiationModels/radiationModel/fvDOM/fvDOM/fvDOM.C b/src/thermophysicalModels/radiationModels/radiationModel/fvDOM/fvDOM/fvDOM.C
index 10961db032..b147a9f7ed 100644
--- a/src/thermophysicalModels/radiationModels/radiationModel/fvDOM/fvDOM/fvDOM.C
+++ b/src/thermophysicalModels/radiationModels/radiationModel/fvDOM/fvDOM/fvDOM.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -27,6 +27,7 @@ License
 #include "absorptionEmissionModel.H"
 #include "scatterModel.H"
 #include "constants.H"
+#include "fvm.H"
 
 using namespace Foam::constant;
 using namespace Foam::constant::mathematical;
@@ -172,10 +173,39 @@ void Foam::radiation::fvDOM::initialise()
     Info<< "fvDOM : Allocated " << IRay_.size()
         << " rays with average orientation:" << nl;
 
-    forAll(IRay_, i)
+    if (cacheDiv_)
     {
-        Info<< '\t' << IRay_[i].I().name()
-            << '\t' << IRay_[i].dAve() << nl;
+        Info<< "Caching div fvMatrix..."<< endl;
+        for (label lambdaI = 0; lambdaI < nLambda_; lambdaI++)
+        {
+            fvRayDiv_[lambdaI].setSize(nRay_);
+
+            forAll(IRay_, rayId)
+            {
+                const surfaceScalarField Ji = IRay_[rayId].dAve() & mesh_.Sf();
+                const volScalarField& iRayLambdaI =
+                    IRay_[rayId].ILambda(lambdaI);
+
+                fvRayDiv_[lambdaI].set
+                (
+                    rayId,
+                    new fvScalarMatrix
+                    (
+                        fvm::div(Ji, iRayLambdaI, "div(Ji,Ii_h)")
+                    )
+                );
+            }
+        }
+    }
+
+    forAll(IRay_, rayId)
+    {
+        if (omegaMax_ <  IRay_[rayId].omega())
+        {
+            omegaMax_ = IRay_[rayId].omega();
+        }
+        Info<< '\t' << IRay_[rayId].I().name() << " : " << "omega : "
+            << '\t' << IRay_[rayId].omega() << nl;
     }
 
     Info<< endl;
@@ -194,7 +224,7 @@ Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
             "G",
             mesh_.time().timeName(),
             mesh_,
-            IOobject::READ_IF_PRESENT,
+            IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
         mesh_,
@@ -247,7 +277,7 @@ Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
             mesh_.time().timeName(),
             mesh_,
             IOobject::NO_READ,
-            IOobject::NO_WRITE
+            IOobject::AUTO_WRITE
         ),
         mesh_,
         dimensionedScalar("a", dimless/dimLength, 0.0)
@@ -260,7 +290,10 @@ Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
     blackBody_(nLambda_, T),
     IRay_(0),
     convergence_(coeffs_.lookupOrDefault<scalar>("convergence", 0.0)),
-    maxIter_(coeffs_.lookupOrDefault<label>("maxIter", 50))
+    maxIter_(coeffs_.lookupOrDefault<label>("maxIter", 50)),
+    fvRayDiv_(nLambda_),
+    cacheDiv_(coeffs_.lookupOrDefault<bool>("cacheDiv", false)),
+    omegaMax_(0)
 {
     initialise();
 }
@@ -346,7 +379,10 @@ Foam::radiation::fvDOM::fvDOM
     blackBody_(nLambda_, T),
     IRay_(0),
     convergence_(coeffs_.lookupOrDefault<scalar>("convergence", 0.0)),
-    maxIter_(coeffs_.lookupOrDefault<label>("maxIter", 50))
+    maxIter_(coeffs_.lookupOrDefault<label>("maxIter", 50)),
+    fvRayDiv_(nLambda_),
+    cacheDiv_(coeffs_.lookupOrDefault<bool>("cacheDiv", false)),
+    omegaMax_(0)
 {
     initialise();
 }
@@ -364,7 +400,6 @@ bool Foam::radiation::fvDOM::read()
     if (radiationModel::read())
     {
         // Only reading solution parameters - not changing ray geometry
-
         coeffs_.readIfPresent("convergence", convergence_);
         coeffs_.readIfPresent("maxIter", maxIter_);
 
@@ -383,19 +418,30 @@ void Foam::radiation::fvDOM::calculate()
 
     updateBlackBodyEmission();
 
+    // Set rays convergence false
+    List<bool> rayIdConv(nRay_, false);
+
     scalar maxResidual = 0.0;
     label radIter = 0;
     do
     {
+        Info<< "Radiation solver iter: " << radIter << endl;
+
         radIter++;
+        maxResidual = 0.0;
         forAll(IRay_, rayI)
         {
-            maxResidual = 0.0;
-            scalar maxBandResidual = IRay_[rayI].correct();
-            maxResidual = max(maxBandResidual, maxResidual);
-        }
+            if (!rayIdConv[rayI])
+            {
+                scalar maxBandResidual = IRay_[rayI].correct();
+                maxResidual = max(maxBandResidual, maxResidual);
 
-        Info<< "Radiation solver iter: " << radIter << endl;
+                if (maxBandResidual < convergence_)
+                {
+                    rayIdConv[rayI] = true;
+                }
+            }
+        }
 
     } while (maxResidual > convergence_ && radIter < maxIter_);
 
@@ -433,7 +479,7 @@ Foam::radiation::fvDOM::Ru() const
     const DimensionedField<scalar, volMesh> E =
         absorptionEmission_->ECont()().dimensionedInternalField();
     const DimensionedField<scalar, volMesh> a =
-        a_.dimensionedInternalField(); //absorptionEmission_->aCont()()
+        a_.dimensionedInternalField();
 
     return  a*G - E;
 }
diff --git a/src/thermophysicalModels/radiationModels/radiationModel/fvDOM/fvDOM/fvDOM.H b/src/thermophysicalModels/radiationModels/radiationModel/fvDOM/fvDOM/fvDOM.H
index 48efae9ee6..3d9950db4c 100644
--- a/src/thermophysicalModels/radiationModels/radiationModel/fvDOM/fvDOM/fvDOM.H
+++ b/src/thermophysicalModels/radiationModels/radiationModel/fvDOM/fvDOM/fvDOM.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -30,6 +30,7 @@ Description
     directions in a participating media, not including scatter.
 
     Available absorption models:
+        constantAbsorptionEmission
         greyMeanAbsoprtionEmission
         wideBandAbsorptionEmission
 
@@ -37,9 +38,15 @@ Description
     \verbatim
         fvDOMCoeffs
         {
-            nPhi    1;          // azimuthal angles in PI/2 on X-Y.(from Y to X)
-            nTheta  2;          // polar angles in PI (from Z to X-Y plane)
-            convergence 1e-4;   // convergence criteria for radiation iteration
+            nPhi        4;          // azimuthal angles in PI/2 on X-Y.
+                                    //(from Y to X)
+            nTheta      0;          // polar angles in PI (from Z to X-Y plane)
+            convergence 1e-3;       // convergence criteria for radiation
+                                    //iteration
+            maxIter     4;          // maximum number of iterations
+            cacheDiv    true;       // cache the div of the RTE equation.
+            //NOTE: Caching div is "only" accurate if the upwind scheme is used
+            //in div(Ji,Ii_h)
         }
 
         solverFreq   1; // Number of flow iterations per radiation iteration
@@ -61,6 +68,7 @@ SourceFiles
 
 #include "radiativeIntensityRay.H"
 #include "radiationModel.H"
+#include "fvMatrices.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -79,6 +87,7 @@ class fvDOM
 {
     // Private data
 
+
         //- Incident radiation  [W/m2]
         volScalarField G_;
 
@@ -121,6 +130,15 @@ class fvDOM
         //- Maximum number of iterations
         scalar maxIter_;
 
+        //- List of cached fvMatrices for rays
+        List<PtrList<fvScalarMatrix> >fvRayDiv_;
+
+        //- Cache convection div matrix
+        bool cacheDiv_;
+
+        //- Maximum omega weight
+        scalar omegaMax_;
+
 
     // Private Member Functions
 
@@ -229,6 +247,19 @@ public:
 
             //- Const access to black body
             inline const blackBodyEmission& blackBody() const;
+
+            //- Const access to cached fvMatrix
+            inline const fvScalarMatrix& fvRayDiv
+            (
+                const label lambdaI,
+                const label rayId
+            ) const;
+
+            //- Caching div(Ji, Ilamda)
+            inline bool cacheDiv() const;
+
+            //- Return omegaMax
+            inline scalar omegaMax() const;
 };
 
 
diff --git a/src/thermophysicalModels/radiationModels/radiationModel/fvDOM/fvDOM/fvDOMI.H b/src/thermophysicalModels/radiationModels/radiationModel/fvDOM/fvDOM/fvDOMI.H
index de6e44f523..8bc9c6425e 100644
--- a/src/thermophysicalModels/radiationModels/radiationModel/fvDOM/fvDOM/fvDOMI.H
+++ b/src/thermophysicalModels/radiationModels/radiationModel/fvDOM/fvDOM/fvDOMI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -102,6 +102,7 @@ inline const Foam::volScalarField& Foam::radiation::fvDOM::Qem() const
     return Qem_;
 }
 
+
 inline const Foam::radiation::blackBodyEmission&
 Foam::radiation::fvDOM::blackBody() const
 {
@@ -109,4 +110,26 @@ Foam::radiation::fvDOM::blackBody() const
 }
 
 
+inline const Foam::fvScalarMatrix& Foam::radiation::fvDOM::fvRayDiv
+(
+    const label rayId,
+    const label lambdaI
+) const
+{
+    return fvRayDiv_[lambdaI][rayId];
+}
+
+
+inline bool Foam::radiation::fvDOM::cacheDiv() const
+{
+    return cacheDiv_;
+}
+
+
+inline Foam::scalar Foam::radiation::fvDOM::omegaMax() const
+{
+    return omegaMax_;
+}
+
+
 // ************************************************************************* //
diff --git a/src/thermophysicalModels/radiationModels/radiationModel/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C b/src/thermophysicalModels/radiationModels/radiationModel/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C
index 947e3e5030..28644026a4 100644
--- a/src/thermophysicalModels/radiationModels/radiationModel/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C
+++ b/src/thermophysicalModels/radiationModels/radiationModel/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -113,7 +113,8 @@ Foam::radiation::radiativeIntensityRay::radiativeIntensityRay
     phi_(phi),
     omega_(0.0),
     nLambda_(nLambda),
-    ILambda_(nLambda)
+    ILambda_(nLambda),
+    myRayId_(rayId)
 {
     scalar sinTheta = Foam::sin(theta);
     scalar cosTheta = Foam::cos(theta);
@@ -135,7 +136,6 @@ Foam::radiation::radiativeIntensityRay::radiativeIntensityRay
         0.5*deltaPhi*Foam::sin(2.0*theta)*Foam::sin(deltaTheta)
     );
 
-
     autoPtr<volScalarField> IDefaultPtr;
 
     forAll(ILambda_, lambdaI)
@@ -213,29 +213,51 @@ Foam::scalar Foam::radiation::radiativeIntensityRay::correct()
     {
         const volScalarField& k = dom_.aLambda(lambdaI);
 
-        const surfaceScalarField Ji(dAve_ & mesh_.Sf());
+        tmp<fvScalarMatrix> IiEq;
 
-        fvScalarMatrix IiEq
+        if (!dom_.cacheDiv())
+        {
+            const surfaceScalarField Ji(dAve_ & mesh_.Sf());
+
+            IiEq =
+            (
+                fvm::div(Ji, ILambda_[lambdaI], "div(Ji,Ii_h)")
+              + fvm::Sp(k*omega_, ILambda_[lambdaI])
+            ==
+                1.0/constant::mathematical::pi*omega_
+              * (
+                    k*blackBody_.bLambda(lambdaI)
+                  + absorptionEmission_.ECont(lambdaI)/4
+                )
+            );
+        }
+        else
+        {
+            IiEq =
+            (
+               dom_.fvRayDiv(myRayId_, lambdaI)
+             + fvm::Sp(k*omega_, ILambda_[lambdaI])
+           ==
+               1.0/constant::mathematical::pi*omega_
+             * (
+                   k*blackBody_.bLambda(lambdaI)
+                 + absorptionEmission_.ECont(lambdaI)/4
+               )
+            );
+        }
+
+        IiEq().relax();
+
+        const solverPerformance ILambdaSol = solve
         (
-            fvm::div(Ji, ILambda_[lambdaI], "div(Ji,Ii_h)")
-          + fvm::Sp(k*omega_, ILambda_[lambdaI])
-         ==
-            1.0/constant::mathematical::pi*omega_
-           *(
-                k*blackBody_.bLambda(lambdaI)
-              + absorptionEmission_.ECont(lambdaI)/4
-            )
+            IiEq(),
+            mesh_.solver("Ii")
         );
 
-        IiEq.relax();
+        const scalar initialRes =
+            ILambdaSol.initialResidual()*omega_/dom_.omegaMax();
 
-        scalar eqnResidual = solve
-        (
-            IiEq,
-            mesh_.solver("Ii")
-        ).initialResidual();
-
-        maxResidual = max(eqnResidual, maxResidual);
+        maxResidual = max(initialRes, maxResidual);
     }
 
     return maxResidual;
diff --git a/src/thermophysicalModels/radiationModels/radiationModel/fvDOM/radiativeIntensityRay/radiativeIntensityRay.H b/src/thermophysicalModels/radiationModels/radiationModel/fvDOM/radiativeIntensityRay/radiativeIntensityRay.H
index 83eed35d99..25f363044e 100644
--- a/src/thermophysicalModels/radiationModels/radiationModel/fvDOM/radiativeIntensityRay/radiativeIntensityRay.H
+++ b/src/thermophysicalModels/radiationModels/radiationModel/fvDOM/radiativeIntensityRay/radiativeIntensityRay.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -38,6 +38,7 @@ SourceFiles
 #include "absorptionEmissionModel.H"
 #include "blackBodyEmission.H"
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -111,6 +112,9 @@ private:
         //- Global ray id - incremented in constructor
         static label rayId;
 
+        //- My ray Id
+        label myRayId_;
+
 
     // Private Member Functions
 
@@ -209,6 +213,7 @@ public:
 
             //- Return the radiative intensity for a given wavelength
             inline const volScalarField& ILambda(const label lambdaI) const;
+
 };
 
 
diff --git a/src/thermophysicalModels/radiationModels/submodels/absorptionEmissionModel/absorptionEmissionModel/absorptionEmissionModel.C b/src/thermophysicalModels/radiationModels/submodels/absorptionEmissionModel/absorptionEmissionModel/absorptionEmissionModel.C
index 982c7e3137..563b46727f 100644
--- a/src/thermophysicalModels/radiationModels/submodels/absorptionEmissionModel/absorptionEmissionModel/absorptionEmissionModel.C
+++ b/src/thermophysicalModels/radiationModels/submodels/absorptionEmissionModel/absorptionEmissionModel/absorptionEmissionModel.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -252,8 +252,8 @@ void Foam::radiation::absorptionEmissionModel::correct
     PtrList<volScalarField>& aj
 ) const
 {
-    a.internalField() = this->a();
-    aj[0].internalField() =  a.internalField();
+    a = this->a();
+    aj[0] =  a;
 }
 
 
diff --git a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/0/G b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/0/G
new file mode 100644
index 0000000000..fd53ea0248
--- /dev/null
+++ b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/0/G
@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      G;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 0 -3 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    ".*"
+    {
+        type            MarshakRadiation;
+        T               T;
+        emissivityMode  lookup;
+        emissivity      uniform 1.0;
+        value           uniform 0;
+    }
+
+    "(region0_to.*)"
+    {
+        type            MarshakRadiation;
+        T               T;
+        emissivityMode  solidRadiation;
+        value           uniform 0;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/radiationProperties b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/radiationProperties
index 89ee652071..132a150419 100644
--- a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/radiationProperties
+++ b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/radiationProperties
@@ -24,22 +24,24 @@ fvDOMCoeffs
 {
     nPhi    3;          // azimuthal angles in PI/2 on X-Y.(from Y to X)
     nTheta  6;          // polar angles in PI (from Z to X-Y plane)
-    convergence 1e-4;   // convergence criteria for radiation iteration
-    maxIter 2;         // maximum number of iterations
+    convergence 0.05;   // convergence criteria for radiation iteration
+    maxIter 3;          // maximum number of iterations
+    cacheDiv    true;   // cache the div of the RTE equation.
+
+// NOTE: Caching div is "only" accurate if the upwind scheme is used in
+// div(Ji,Ii_h)
 }
 
 // Number of flow iterations per radiation iteration
 solverFreq 10;
 
 absorptionEmissionModel constantAbsorptionEmission;
-//absorptionEmissionModel greyMeanAbsorptionEmission;
 
 constantAbsorptionEmissionCoeffs
 {
-    absorptivity    absorptivity [ 0 -1 0 0 0 0 0 ] 0.01;
-    emissivity      emissivity [ 0 -1 0 0 0 0 0 ] 0.01;
-    E               E [ 1 -1 -3 0 0 0 0 ] 0;
-}
+    absorptivity    absorptivity    [ m^-1 ]       0.01;
+    emissivity      emissivity      [ m^-1 ]       0.01;
+    E               E               [ kg m^-1 s^-3 ]  0;
 
 greyMeanAbsorptionEmissionCoeffs
 {
@@ -49,7 +51,7 @@ greyMeanAbsorptionEmissionCoeffs
 
     CO2
     {
-        Tcommon         300;   //Common Temp
+        Tcommon         200;   //Common Temp
         invTemp         true;   //Is the polynomio using inverse temperature.
         Tlow            200;   //Low Temp
         Thigh           2500;  //High Temp
@@ -77,7 +79,7 @@ greyMeanAbsorptionEmissionCoeffs
 
     H2O
     {
-        Tcommon         300;
+        Tcommon         200;
         invTemp         true;
         Tlow            200;
         Thigh           2500;
@@ -104,7 +106,7 @@ greyMeanAbsorptionEmissionCoeffs
 
     C3H8//CH4
     {
-        Tcommon         300;
+        Tcommon         200;
         Tlow            200;
         Thigh           2000;
         invTemp         false;
@@ -131,7 +133,7 @@ greyMeanAbsorptionEmissionCoeffs
 
     O2
     {
-        Tcommon         300;
+        Tcommon         200;
         invTemp         true;
         Tlow            200;
         Thigh           2500;
@@ -159,7 +161,7 @@ greyMeanAbsorptionEmissionCoeffs
 
     N2
     {
-        Tcommon         300;
+        Tcommon         200;
         invTemp         true;
         Tlow            200;
         Thigh           2500;
diff --git a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/system/controlDict b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/system/controlDict
index 3efab07aeb..eb03206ed1 100644
--- a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/system/controlDict
+++ b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/system/controlDict
@@ -16,7 +16,7 @@ FoamFile
 
 application     fireFoam;
 
-startFrom       startTime;
+startFrom       latestTime;
 
 startTime       0;
 
@@ -28,11 +28,11 @@ deltaT          0.03;
 
 writeControl    adjustableRunTime;
 
-writeInterval   1;
+writeInterval   0.5;
 
 purgeWrite      0;
 
-writeFormat     ascii;
+writeFormat     binary;
 
 writePrecision  12;
 
diff --git a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/system/fvSolution b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/system/fvSolution
index 3a3bd03c73..7345991f71 100644
--- a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/system/fvSolution
+++ b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/system/fvSolution
@@ -34,7 +34,7 @@ solvers
     p_rgh
     {
         solver              GAMG;
-        tolerance           1e-5;
+        tolerance           1e-7;
         relTol              0.01;
         smoother            GaussSeidel;
         cacheAgglomeration  true;
@@ -45,28 +45,25 @@ solvers
 
     p_rghFinal
     {
-        solver              GAMG;
-        tolerance           1e-6;
+        $p_rgh;
+        tolerance           1e-7;
         relTol              0;
-        smoother            GaussSeidel;
-        cacheAgglomeration  true;
-        nCellsInCoarsestLevel 10;
-        agglomerator     faceAreaPair;
-        mergeLevels      1;
     };
 
 
-    "(U|k)"
+    "(U|Yi|k|h|omega)"
     {
         solver          PBiCG;
         preconditioner  DILU;
-        tolerance       1e-6;
-        relTol          0.01;
+        tolerance       1e-7;
+        relTol          0.1;
+        nSweeps         1;
     };
 
-    "(U|k)Final"
+    "(U|Yi|k|h|omega)Final"
     {
         $U;
+        tolerance       1e-7;
         relTol          0;
     };
 
@@ -84,18 +81,21 @@ solvers
         solver              GAMG;
         tolerance           1e-4;
         relTol              0;
-        smoother            DILU;
+        smoother            symGaussSeidel;
         cacheAgglomeration  true;
         nCellsInCoarsestLevel   10;
         agglomerator    faceAreaPair;
-        mergeLevels     1;
+        mergeLevels         1;
+        maxIter             1;
+        nPreSweeps          0;
+        nPostSweeps         1;
     }
 
-    G
+    "(G)Final"
     {
         solver          PCG;
         preconditioner  DIC;
-        tolerance       1e-06;
+        tolerance       1e-04;
         relTol          0;
     }
 
@@ -104,8 +104,8 @@ solvers
 PIMPLE
 {
     momentumPredictor yes;
-    nOuterCorrectors  2;
-    nCorrectors       1;
+    nOuterCorrectors  1;
+    nCorrectors       3;
     nNonOrthogonalCorrectors 0;
 }
 
@@ -117,7 +117,7 @@ relaxationFactors
     equations
     {
         "(U|k).*"                   1;
-        "(C3H8|O2|H2O|CO2|h).*"    1;
+        "(C3H8|O2|H2O|CO2|h).*"     1;
     }
 }
 
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/radiationProperties b/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/radiationProperties
index 9ff932edc9..62b677bd77 100644
--- a/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/radiationProperties
+++ b/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/radiationProperties
@@ -19,10 +19,14 @@ radiationModel  fvDOM;
 
 fvDOMCoeffs
 {
-    nPhi    4;          // azimuthal angles in PI/2 on X-Y.(from Y to X)
-    nTheta  0;          // polar angles in PI (from Z to X-Y plane)
-    convergence 1e-3;   // convergence criteria for radiation iteration
-    maxIter 1;          // maximum number of iterations
+    nPhi        4;          // azimuthal angles in PI/2 on X-Y.(from Y to X)
+    nTheta      0;          // polar angles in PI (from Z to X-Y plane)
+    convergence 1e-2;       // convergence criteria for radiation iteration
+    maxIter     4;          // maximum number of iterations
+    cacheDiv    true;       // cache the div of the RTE equation.
+
+// NOTE: Caching div is "only" accurate if the upwind scheme is used in
+// div(Ji,Ii_h)
 }
 
 // Number of flow iterations per radiation iteration
@@ -32,9 +36,9 @@ absorptionEmissionModel greyMeanAbsorptionEmission;
 
 constantAbsorptionEmissionCoeffs
 {
-   absorptivity    absorptivity    [ 0 -1 0 0 0 0 0 ] 0.1;
-   emissivity      emissivity      [ 0 -1 0 0 0 0 0 ] 0.1;
-   E               E   [ 1 -1 -3 0 0 0 0 ] 0;
+    absorptivity    absorptivity    [ m^-1 ]       0.01;
+    emissivity      emissivity      [ m^-1 ]       0.01;
+    E               E               [ kg m^-1 s^-3 ]  0;
 }
 
 greyMeanAbsorptionEmissionCoeffs
@@ -45,7 +49,7 @@ greyMeanAbsorptionEmissionCoeffs
 
     CO2
     {
-        Tcommon         300;   //Common Temp
+        Tcommon         200;   //Common Temp
         invTemp         true;   //Is the polynomio using inverse temperature.
         Tlow            200;   //Low Temp
         Thigh           2500;  //High Temp
@@ -73,7 +77,7 @@ greyMeanAbsorptionEmissionCoeffs
 
     H2O
     {
-        Tcommon         300;
+        Tcommon         200;
         invTemp         true;
         Tlow            200;
         Thigh           2500;
@@ -100,7 +104,7 @@ greyMeanAbsorptionEmissionCoeffs
 
     CH4
     {
-        Tcommon         300;
+        Tcommon         200;
         Tlow            200;
         Thigh           2500;
         invTemp         false;
@@ -127,7 +131,7 @@ greyMeanAbsorptionEmissionCoeffs
 
     O2
     {
-        Tcommon         300;
+        Tcommon         200;
         invTemp         true;
         Tlow            200;
         Thigh           2500;
@@ -155,7 +159,7 @@ greyMeanAbsorptionEmissionCoeffs
 
     N2
     {
-        Tcommon         300;
+        Tcommon         200;
         invTemp         true;
         Tlow            200;
         Thigh           2500;
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire2D/system/controlDict b/tutorials/combustion/fireFoam/les/smallPoolFire2D/system/controlDict
index 8420fc7598..689f2a866c 100644
--- a/tutorials/combustion/fireFoam/les/smallPoolFire2D/system/controlDict
+++ b/tutorials/combustion/fireFoam/les/smallPoolFire2D/system/controlDict
@@ -48,8 +48,6 @@ runTimeModifiable yes;
 
 adjustTimeStep  yes;
 
-maxCo           0.25;
-
-maxDeltaT       0.1;
+maxCo           0.5;
 
 // ************************************************************************* //
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire2D/system/fvSolution b/tutorials/combustion/fireFoam/les/smallPoolFire2D/system/fvSolution
index d7bebcf0bb..fba625916d 100644
--- a/tutorials/combustion/fireFoam/les/smallPoolFire2D/system/fvSolution
+++ b/tutorials/combustion/fireFoam/les/smallPoolFire2D/system/fvSolution
@@ -72,21 +72,24 @@ solvers
 
     Ii
     {
-        solver              GAMG;
+         solver              GAMG;
         tolerance           1e-4;
-        relTol              0;
-        smoother            DILU;
+        relTol              0.1;
+        smoother            symGaussSeidel;
         cacheAgglomeration  true;
         nCellsInCoarsestLevel   10;
-        agglomerator    faceAreaPair;
-        mergeLevels     1;
+        agglomerator      faceAreaPair;
+        mergeLevels         1;
+        maxIter             3;
+        nPreSweeps          0;
+        nPostSweeps         1;
     }
 
     G
     {
         solver          PCG;
         preconditioner  DIC;
-        tolerance       1e-06;
+        tolerance       1e-04;
         relTol          0;
     }
 
@@ -108,7 +111,7 @@ relaxationFactors
     equations
     {
         "(U|k).*"                   1;
-        "(C3H8|O2|H2O|CO2|h).*"    0.9;
+        "(C3H8|O2|H2O|CO2|h).*"     1;
     }
 }
 
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire3D/0/IDefault b/tutorials/combustion/fireFoam/les/smallPoolFire3D/0/IDefault
index e70bf3dec0..1da0ef50b9 100644
--- a/tutorials/combustion/fireFoam/les/smallPoolFire3D/0/IDefault
+++ b/tutorials/combustion/fireFoam/les/smallPoolFire3D/0/IDefault
@@ -25,7 +25,7 @@ boundaryField
         type            greyDiffusiveRadiation;
         T               T;
         emissivityMode  lookup;
-        emissivity      uniform 1.0;
+        emissivity      uniform 1;
         value           uniform 0;
     }
 }
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/radiationProperties b/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/radiationProperties
index 513c11e7db..63e3de884c 100644
--- a/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/radiationProperties
+++ b/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/radiationProperties
@@ -20,10 +20,14 @@ radiationModel  fvDOM;
 
 fvDOMCoeffs
 {
-    nPhi    2;          // azimuthal angles in PI/2 on X-Y.(from Y to X)
-    nTheta  2;          // polar angles in PI (from Z to X-Y plane)
-    convergence 1e-3;   // convergence criteria for radiation iteration
-    maxIter 10;         // maximum number of iterations
+    nPhi        2;          // azimuthal angles in PI/2 on X-Y.(from Y to X)
+    nTheta      2;          // polar angles in PI (from Z to X-Y plane)
+    convergence 1e-1;       // convergence criteria for radiation iteration
+    maxIter     1;          // maximum number of iterations
+    cacheDiv    true;       // cache the div of the RTE equation.
+
+//  NOTE: Caching div is "only" accurate if the upwind scheme is used in
+//  div(Ji,Ii_h)
 }
 
 // Number of flow iterations per radiation iteration
@@ -33,9 +37,9 @@ absorptionEmissionModel greyMeanAbsorptionEmission;
 
 constantAbsorptionEmissionCoeffs
 {
-    absorptivity    absorptivity    [ 0 -1 0 0 0 0 0 ] 0.01;
-    emissivity      emissivity      [ 0 -1 0 0 0 0 0 ] 0.01;
-    E               E   [ 1 -1 -3 0 0 0 0 ] 0;
+    absorptivity    absorptivity    [ m^-1 ]       0.01;
+    emissivity      emissivity      [ m^-1 ]       0.01;
+    E               E               [ kg m^-1 s^-3 ]  0;
 }
 
 greyMeanAbsorptionEmissionCoeffs
@@ -46,7 +50,7 @@ greyMeanAbsorptionEmissionCoeffs
 
     CO2
     {
-        Tcommon         300;   //Common Temp
+        Tcommon         200;   //Common Temp
         invTemp         true;   //Is the polynomio using inverse temperature.
         Tlow            200;   //Low Temp
         Thigh           2500;  //High Temp
@@ -74,7 +78,7 @@ greyMeanAbsorptionEmissionCoeffs
 
     H2O
     {
-        Tcommon         300;
+        Tcommon         200;
         invTemp         true;
         Tlow            200;
         Thigh           2500;
@@ -101,7 +105,7 @@ greyMeanAbsorptionEmissionCoeffs
 
     CH4
     {
-        Tcommon         300;
+        Tcommon         200;
         Tlow            200;
         Thigh           2500;
         invTemp         false;
@@ -128,7 +132,7 @@ greyMeanAbsorptionEmissionCoeffs
 
     O2
     {
-        Tcommon         300;
+        Tcommon         200;
         invTemp         true;
         Tlow            200;
         Thigh           2500;
@@ -156,7 +160,7 @@ greyMeanAbsorptionEmissionCoeffs
 
     N2
     {
-        Tcommon         300;
+        Tcommon         200;
         invTemp         true;
         Tlow            200;
         Thigh           2500;
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire3D/system/controlDict b/tutorials/combustion/fireFoam/les/smallPoolFire3D/system/controlDict
index 6a45d3a30f..6f9c9cf1e6 100644
--- a/tutorials/combustion/fireFoam/les/smallPoolFire3D/system/controlDict
+++ b/tutorials/combustion/fireFoam/les/smallPoolFire3D/system/controlDict
@@ -22,13 +22,13 @@ startTime       0.0;
 
 stopAt          endTime;
 
-endTime         2.75;
+endTime         0.75;
 
 deltaT          0.001;
 
 writeControl    adjustableRunTime;
 
-writeInterval   0.05;
+writeInterval   0.5;
 
 purgeWrite      0;
 
@@ -48,7 +48,7 @@ runTimeModifiable yes;
 
 adjustTimeStep  yes;
 
-maxCo           0.2;
+maxCo           0.5;
 
 maxDeltaT       0.1;
 
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire3D/system/fvSolution b/tutorials/combustion/fireFoam/les/smallPoolFire3D/system/fvSolution
index 829f7a4a7a..11cfaf95bd 100644
--- a/tutorials/combustion/fireFoam/les/smallPoolFire3D/system/fvSolution
+++ b/tutorials/combustion/fireFoam/les/smallPoolFire3D/system/fvSolution
@@ -66,12 +66,14 @@ solvers
         solver              GAMG;
         tolerance           1e-4;
         relTol              0;
-        smoother            DILU;
+        smoother            symGaussSeidel;
         cacheAgglomeration  true;
         nCellsInCoarsestLevel   10;
-        agglomerator    faceAreaPair;
-        mergeLevels     1;
-        maxIter         20;
+        agglomerator        faceAreaPair;
+        mergeLevels         1;
+        maxIter             1;
+        nPreSweeps          0;
+        nPostSweeps         1;
     }
 
     G
@@ -88,7 +90,7 @@ PIMPLE
 {
     momentumPredictor yes;
     nOuterCorrectors  1;
-    nCorrectors       2;
+    nCorrectors       3;
     nNonOrthogonalCorrectors 0;
 }