diff --git a/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.C b/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.C
index d1cdcad666..164534ff57 100644
--- a/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.C
+++ b/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.C
@@ -190,31 +190,6 @@ void Foam::radiation::fvDOM::initialise()
     Info<< "fvDOM : Allocated " << IRay_.size()
         << " rays with average orientation:" << nl;
 
-    if (cacheDiv_)
-    {
-        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())
@@ -308,8 +283,6 @@ Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
     IRay_(0),
     convergence_(coeffs_.lookupOrDefault<scalar>("convergence", 0.0)),
     maxIter_(coeffs_.lookupOrDefault<label>("maxIter", 50)),
-    fvRayDiv_(nLambda_),
-    cacheDiv_(coeffs_.lookupOrDefault<bool>("cacheDiv", false)),
     omegaMax_(0)
 {
     initialise();
@@ -397,8 +370,6 @@ Foam::radiation::fvDOM::fvDOM
     IRay_(0),
     convergence_(coeffs_.lookupOrDefault<scalar>("convergence", 0.0)),
     maxIter_(coeffs_.lookupOrDefault<label>("maxIter", 50)),
-    fvRayDiv_(nLambda_),
-    cacheDiv_(coeffs_.lookupOrDefault<bool>("cacheDiv", false)),
     omegaMax_(0)
 {
     initialise();
diff --git a/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.H b/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.H
index 08823fa28b..1154b6cafc 100644
--- a/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.H
+++ b/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.H
@@ -44,9 +44,6 @@ Description
             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
@@ -130,12 +127,6 @@ class fvDOM
         //- Maximum number of iterations
         label maxIter_;
 
-        //- List of cached fvMatrices for rays
-        List<PtrList<fvScalarMatrix>>fvRayDiv_;
-
-        //- Cache convection div matrix
-        bool cacheDiv_;
-
         //- Maximum omega weight
         scalar omegaMax_;
 
@@ -248,16 +239,6 @@ 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/radiation/radiationModels/fvDOM/fvDOM/fvDOMI.H b/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOMI.H
index 8e4c201b31..d2255b2de0 100644
--- a/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOMI.H
+++ b/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOMI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -110,22 +110,6 @@ 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/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C b/src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C
index 4afda38d4a..d776541384 100644
--- a/src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C
+++ b/src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C
@@ -214,50 +214,28 @@ Foam::scalar Foam::radiation::radiativeIntensityRay::correct()
     {
         const volScalarField& k = dom_.aLambda(lambdaI);
 
-        tmp<fvScalarMatrix> IiEq;
+        const surfaceScalarField Ji(dAve_ & mesh_.Sf());
 
-        if (!dom_.cacheDiv())
-        {
-            const surfaceScalarField Ji(dAve_ & mesh_.Sf());
+        fvScalarMatrix IiEq
+        (
+            fvm::div(Ji, ILambda_[lambdaI], "div(Ji,Ii_h)")
+          + fvm::Sp(k*omega_, ILambda_[lambdaI])
+        ==
+            1.0/constant::mathematical::pi*omega_
+           *(
+                // Remove aDisp from k
+                (k - absorptionEmission_.aDisp(lambdaI))
+               *blackBody_.bLambda(lambdaI)
 
-            IiEq =
-            (
-                fvm::div(Ji, ILambda_[lambdaI], "div(Ji,Ii_h)")
-              + fvm::Sp(k*omega_, ILambda_[lambdaI])
-            ==
-                1.0/constant::mathematical::pi*omega_
-               *(
-                    // Remove aDisp from k
-                    (k - absorptionEmission_.aDisp(lambdaI))
-                   *blackBody_.bLambda(lambdaI)
+              + absorptionEmission_.E(lambdaI)/4
+            )
+        );
 
-                  + absorptionEmission_.E(lambdaI)/4
-                )
-            );
-        }
-        else
-        {
-            IiEq =
-            (
-               dom_.fvRayDiv(myRayId_, lambdaI)
-             + fvm::Sp(k*omega_, ILambda_[lambdaI])
-           ==
-               1.0/constant::mathematical::pi*omega_
-             * (
-                   // Remove aDisp from k
-                   (k - absorptionEmission_.aDisp(lambdaI))
-                  *blackBody_.bLambda(lambdaI)
-
-                 + absorptionEmission_.E(lambdaI)/4
-               )
-            );
-        }
-
-        IiEq.ref().relax();
+        IiEq.relax();
 
         const solverPerformance ILambdaSol = solve
         (
-            IiEq.ref(),
+            IiEq,
             mesh_.solver("Ii")
         );
 
diff --git a/tutorials/combustion/fireFoam/les/flameSpreadWaterSuppressionPanel/constant/radiationProperties b/tutorials/combustion/fireFoam/les/flameSpreadWaterSuppressionPanel/constant/radiationProperties
index 242c5aace4..c8eb0182a8 100644
--- a/tutorials/combustion/fireFoam/les/flameSpreadWaterSuppressionPanel/constant/radiationProperties
+++ b/tutorials/combustion/fireFoam/les/flameSpreadWaterSuppressionPanel/constant/radiationProperties
@@ -26,7 +26,6 @@ fvDOMCoeffs
     nTheta  2;           // polar angles in PI (from Z to X-Y plane)
     convergence 1e-2;   // convergence criteria for radiation iteration
     maxIter 3;          // maximum number of iterations
-    cacheDiv    true;   // cache the div of the RTE equation.
 }
 
 // Number of flow iterations per radiation iteration
diff --git a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/radiationProperties b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/radiationProperties
index b779d05d38..53e91a9c76 100644
--- a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/radiationProperties
+++ b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/radiationProperties
@@ -26,10 +26,6 @@ fvDOMCoeffs
     nTheta  6;          // polar angles in PI (from Z to X-Y plane)
     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
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/radiationProperties b/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/radiationProperties
index a25c81382a..d8da2d50ca 100644
--- a/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/radiationProperties
+++ b/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/radiationProperties
@@ -25,10 +25,6 @@ fvDOMCoeffs
     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
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/radiationProperties b/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/radiationProperties
index 621656b060..d7a5e81b93 100644
--- a/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/radiationProperties
+++ b/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/radiationProperties
@@ -26,10 +26,6 @@ fvDOMCoeffs
     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
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/hotRadiationRoomFvDOM/constant/radiationProperties b/tutorials/heatTransfer/buoyantSimpleFoam/hotRadiationRoomFvDOM/constant/radiationProperties
index 76bbe562b8..0c47e3668c 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/hotRadiationRoomFvDOM/constant/radiationProperties
+++ b/tutorials/heatTransfer/buoyantSimpleFoam/hotRadiationRoomFvDOM/constant/radiationProperties
@@ -25,7 +25,6 @@ fvDOMCoeffs
     nTheta      5;          // polar angles in PI (from Z to X-Y plane)
     convergence 1e-3;   // convergence criteria for radiation iteration
     maxIter     10;         // maximum number of iterations
-    cacheDiv    false;
 }
 
 // Number of flow iterations per radiation iteration