Overhauled fvDOM and added specular reflective BCs with tutorial case

src/thermophysicalModels/radiation/Make/files

@@ -33,6 +33,7 @@ derivedFvPatchFields/MarshakRadiation/MarshakRadiationMixedFvPatchScalarField.C
 derivedFvPatchFields/MarshakRadiationFixedT/MarshakRadiationFixedTMixedFvPatchScalarField.C
 derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C
 derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.C
+derivedFvPatchFields/wideBandSpecularRadiation/wideBandSpecularRadiationMixedFvPatchScalarField.C
 derivedFvPatchFields/greyDiffusiveViewFactor/greyDiffusiveViewFactorFixedValueFvPatchScalarField.C
 
 LIB = $(FOAM_LIBBIN)/libradiation

src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C

@@ -29,13 +29,64 @@ License
 #include "volFields.H"
 
 #include "fvDOM.H"
-#include "radiationConstants.H"
+
-#include "mathematicalConstants.H"
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace radiation
+{
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void
+greyDiffusiveRadiationMixedFvPatchScalarField::calcSumOutgoingAngles() const
+{
+    if (sumOutgoingAnglesPtr_)
+    {
+        FatalErrorIn
+        (
+            "wideBandDiffusiveRadiationMixedFvPatchScalarField"
+            "::calcSumOutgoingAngles()"
+        )   << "sumOutgoingAnglesPtr_ already calculated"
+            << abort(FatalError);
+    }
+
+    sumOutgoingAnglesPtr_ = new scalarField(this->size(), 0.);
+    scalarField& sumOutgoingAngles = *sumOutgoingAnglesPtr_;
+
+
+    // Get access to radiation model, and recast as fvDOM
+    const radiationModel& radiation =
+        db().lookupObject<radiationModel>("radiationProperties");
+
+    const fvDOM& dom = dynamic_cast<const fvDOM&>(radiation);
+
+
+    for(label rayI = 0; rayI < dom.nRay(); rayI++)
+    {
+        // Calculate cosine of angle between face and ray
+        scalarField outgoingAngles = this->patch().nf() & dom.IRay(rayI).dAve();
+
+        // For outgoing rays, outgoingAngles will be negative
+        sumOutgoingAngles += neg(outgoingAngles)*(-outgoingAngles);
+    }
+}
+
+const Foam::scalarField&
+greyDiffusiveRadiationMixedFvPatchScalarField::sumOutgoingAngles() const
+{
+    if (!sumOutgoingAnglesPtr_)
+    {
+        calcSumOutgoingAngles();
+    }
+
+    return *sumOutgoingAnglesPtr_;
+}
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::radiation::greyDiffusiveRadiationMixedFvPatchScalarField::
+greyDiffusiveRadiationMixedFvPatchScalarField::
 greyDiffusiveRadiationMixedFvPatchScalarField
 (
     const fvPatch& p,
@@ -44,7 +95,8 @@ greyDiffusiveRadiationMixedFvPatchScalarField
 :
     mixedFvPatchScalarField(p, iF),
     TName_("undefinedT"),
-    emissivity_(0.0)
+    emissivity_(0.0),
+    sumOutgoingAnglesPtr_(NULL)
 {
     refValue() = 0.0;
     refGrad() = 0.0;
@@ -52,7 +104,7 @@ greyDiffusiveRadiationMixedFvPatchScalarField
 }
 
 
-Foam::radiation::greyDiffusiveRadiationMixedFvPatchScalarField::
+greyDiffusiveRadiationMixedFvPatchScalarField::
 greyDiffusiveRadiationMixedFvPatchScalarField
 (
     const greyDiffusiveRadiationMixedFvPatchScalarField& ptf,
@@ -63,11 +115,12 @@ greyDiffusiveRadiationMixedFvPatchScalarField
 :
     mixedFvPatchScalarField(ptf, p, iF, mapper),
     TName_(ptf.TName_),
-    emissivity_(ptf.emissivity_)
+    emissivity_(ptf.emissivity_),
+    sumOutgoingAnglesPtr_(NULL)
 {}
 
 
-Foam::radiation::greyDiffusiveRadiationMixedFvPatchScalarField::
+greyDiffusiveRadiationMixedFvPatchScalarField::
 greyDiffusiveRadiationMixedFvPatchScalarField
 (
     const fvPatch& p,
@@ -77,7 +130,8 @@ greyDiffusiveRadiationMixedFvPatchScalarField
 :
     mixedFvPatchScalarField(p, iF),
     TName_(dict.lookup("T")),
-    emissivity_(readScalar(dict.lookup("emissivity")))
+    emissivity_(readScalar(dict.lookup("emissivity"))),
+    sumOutgoingAnglesPtr_(NULL)
 {
     if (dict.found("refValue"))
     {
@@ -110,7 +164,7 @@ greyDiffusiveRadiationMixedFvPatchScalarField
 }
 
 
-Foam::radiation::greyDiffusiveRadiationMixedFvPatchScalarField::
+greyDiffusiveRadiationMixedFvPatchScalarField::
 greyDiffusiveRadiationMixedFvPatchScalarField
 (
     const greyDiffusiveRadiationMixedFvPatchScalarField& ptf
@@ -118,11 +172,12 @@ greyDiffusiveRadiationMixedFvPatchScalarField
 :
     mixedFvPatchScalarField(ptf),
     TName_(ptf.TName_),
-    emissivity_(ptf.emissivity_)
+    emissivity_(ptf.emissivity_),
+    sumOutgoingAnglesPtr_(NULL)
 {}
 
 
-Foam::radiation::greyDiffusiveRadiationMixedFvPatchScalarField::
+greyDiffusiveRadiationMixedFvPatchScalarField::
 greyDiffusiveRadiationMixedFvPatchScalarField
 (
     const greyDiffusiveRadiationMixedFvPatchScalarField& ptf,
@@ -131,100 +186,88 @@ greyDiffusiveRadiationMixedFvPatchScalarField
 :
     mixedFvPatchScalarField(ptf, iF),
     TName_(ptf.TName_),
-    emissivity_(ptf.emissivity_)
+    emissivity_(ptf.emissivity_),
+    sumOutgoingAnglesPtr_(NULL)
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::radiation::greyDiffusiveRadiationMixedFvPatchScalarField::
+void greyDiffusiveRadiationMixedFvPatchScalarField::updateCoeffs()
-updateCoeffs()
 {
     if (this->updated())
     {
         return;
     }
 
-    const scalarField& Tp =
-        lookupPatchField<volScalarField, scalar>(TName_);
 
+
+
+    const label patchI = this->patch().index();
+
+    // Access radiation model
     const radiationModel& radiation =
         db().lookupObject<radiationModel>("radiationProperties");
 
-    const fvDOM& dom(refCast<const fvDOM>(radiation));
+    const fvDOM& dom = dynamic_cast<const fvDOM&>(radiation);
 
+    if (dom.nLambda() == 0)
+    {
+        FatalErrorIn
+        (
+            ""
+            "wideBandDiffusiveRadiationMixedFvPatchScalarField::updateCoeffs"
+        )   << " a non-grey boundary condition is used with a grey "
+            << "absorption model" << nl << exit(FatalError);
+    }
+
+    // Get rayId and lambda Id for this ray
     label rayId = -1;
     label lambdaId = -1;
     dom.setRayIdLambdaId(dimensionedInternalField().name(), rayId, lambdaId);
 
-    const label patchI = patch().index();
+    // Make shortcut to ray belonging to this field
+    const radiativeIntensityRay& ray = dom.IRay(rayId);
 
-    if (dom.nLambda() != 1)
+    // Access incoming radiation for this patch for this band
-    {
+    const scalarField& Qin = dom.Qin(lambdaId)[patchI];
-        FatalErrorIn
+
+    // Access black body radiation for this patch for this band
+    const scalarField& Eb =
+        dom.blackBody().bLambda(lambdaId).boundaryField()[patchI];
+
+    // Get face normals
+    vectorField nHat = patch().nf();
+
+    // Calculate cos of incoming angle of current ray with every face
+    scalarField incomingAngle = this->patch().nf() & ray.dAve();
+
+    // Set to zeroGradient (=0; incomingAngle > 0) for faces with incoming rays
+    // and to fixedValue (=1; incomingAngle < 0) for outgoing rays
+    this->valueFraction() = neg(incomingAngle);
+
+    // Set intensity value for fixedValue part (see reference in header file)
+    this->refValue() =
+        emissivity_*Eb + ((1 - emissivity_)*Qin)/sumOutgoingAngles();
+
+    // Update boundary field now, so values for incoming and outgoing rays
+    // are in balance
+    scalarField::operator=
     (
-            "Foam::radiation::"
+        this->valueFraction()*this->refValue()
-            "greyDiffusiveRadiationMixedFvPatchScalarField::updateCoeffs"
+      +
-        )   << " a grey boundary condition is used with a non-grey "
+        (1.0 - this->valueFraction())*
-            << "absorption model" << nl << exit(FatalError);
-    }
-
-    scalarField& Iw = *this;
-    vectorField n = patch().Sf()/patch().magSf();
-
-    radiativeIntensityRay& ray =
-        const_cast<radiativeIntensityRay&>(dom.IRay(rayId));
-
-    ray.Qr().boundaryField()[patchI] += Iw*(n & ray.dAve());
-
-    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)
-        {
-            // direction out of the wall
-            refGrad()[faceI] = 0.0;
-            valueFraction()[faceI] = 1.0;
-            refValue()[faceI] =
         (
-                    Ir*(1.0 - emissivity_)
+            this->patchInternalField()
-                  + emissivity_*radiation::sigmaSB.value()*pow4(Tp[faceI])
+          + this->refGrad()/this->patch().deltaCoeffs()
         )
-               /mathematicalConstant::pi;
+    );
-
-        }
-        else
-        {
-            // direction into the wall
-            valueFraction()[faceI] = 0.0;
-            refGrad()[faceI] = 0.0;
-            refValue()[faceI] = 0.0; //not used
-        }
-    }
 
     mixedFvPatchScalarField::updateCoeffs();
 }
 
 
-void Foam::radiation::greyDiffusiveRadiationMixedFvPatchScalarField::write
+void greyDiffusiveRadiationMixedFvPatchScalarField::write
 (
     Ostream& os
 ) const
@@ -237,17 +280,13 @@ void Foam::radiation::greyDiffusiveRadiationMixedFvPatchScalarField::write
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-namespace Foam
-{
-namespace radiation
-{
 makePatchTypeField
 (
     fvPatchScalarField,
     greyDiffusiveRadiationMixedFvPatchScalarField
 );
-}
-}
 
+} // End namespace radiation
+} // End namespace Foam
 
 // ************************************************************************* //

src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.H

@@ -25,7 +25,9 @@ Class
     Foam::greyDiffusiveRadiationMixedFvPatchScalarField
 
 Description
-    Radiation temperature specified
+    Radiation boundary will diffusive reflection and black body emission.
+    See Modest "Radiative heat transfer", Capter 16.6 "The finite volume method"
+    Eq.16.63
 
 SourceFiles
     greyDiffusiveRadiationMixedFvPatchScalarField.C
@@ -59,6 +61,15 @@ class greyDiffusiveRadiationMixedFvPatchScalarField
         //- Emissivity
         scalar emissivity_;
 
+        //- Sum of cosine factors of all outgoing rays (lazy evaluated)
+        mutable scalarField *sumOutgoingAnglesPtr_;
+
+
+    // Private member functions
+
+        //- Calculate sum of cosine factors of all outgoing rays
+        void calcSumOutgoingAngles() const;
+
 
 public:
 
@@ -157,6 +168,9 @@ public:
                 return emissivity_;
             }
 
+            //- Return sum of cosine factors of all outgoing rays
+            const scalarField& sumOutgoingAngles() const;
+
 
         // Evaluation functions
 

src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.C

@@ -29,13 +29,64 @@ License
 #include "volFields.H"
 
 #include "fvDOM.H"
-#include "wideBandAbsorptionEmission.H"
+
-#include "radiationConstants.H"
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-#include "mathematicalConstants.H"
+
+namespace Foam
+{
+namespace radiation
+{
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void
+wideBandDiffusiveRadiationMixedFvPatchScalarField::calcSumOutgoingAngles() const
+{
+    if (sumOutgoingAnglesPtr_)
+    {
+        FatalErrorIn
+        (
+            "wideBandDiffusiveRadiationMixedFvPatchScalarField"
+            "::calcSumOutgoingAngles()"
+        )   << "sumOutgoingAnglesPtr_ already calculated"
+            << abort(FatalError);
+    }
+
+    sumOutgoingAnglesPtr_ = new scalarField(this->size(), 0.);
+    scalarField& sumOutgoingAngles = *sumOutgoingAnglesPtr_;
+
+
+    // Get access to radiation model, and recast as fvDOM
+    const radiationModel& radiation =
+        db().lookupObject<radiationModel>("radiationProperties");
+
+    const fvDOM& dom = dynamic_cast<const fvDOM&>(radiation);
+
+
+    for(label rayI = 0; rayI < dom.nRay(); rayI++)
+    {
+        // Calculate cosine of angle between face and ray
+        scalarField outgoingAngles = this->patch().nf() & dom.IRay(rayI).dAve();
+
+        // For outgoing rays, outgoingAngles will be negative
+        sumOutgoingAngles += neg(outgoingAngles)*(-outgoingAngles);
+    }
+}
+
+const Foam::scalarField&
+wideBandDiffusiveRadiationMixedFvPatchScalarField::sumOutgoingAngles() const
+{
+    if (!sumOutgoingAnglesPtr_)
+    {
+        calcSumOutgoingAngles();
+    }
+
+    return *sumOutgoingAnglesPtr_;
+}
+
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::radiation::wideBandDiffusiveRadiationMixedFvPatchScalarField::
+wideBandDiffusiveRadiationMixedFvPatchScalarField::
 wideBandDiffusiveRadiationMixedFvPatchScalarField
 (
     const fvPatch& p,
@@ -44,7 +95,8 @@ wideBandDiffusiveRadiationMixedFvPatchScalarField
 :
     mixedFvPatchScalarField(p, iF),
     TName_("undefinedT"),
-    emissivity_(0.0)
+    emissivity_(0.0),
+    sumOutgoingAnglesPtr_(NULL)
 {
     refValue() = 0.0;
     refGrad() = 0.0;
@@ -52,7 +104,7 @@ wideBandDiffusiveRadiationMixedFvPatchScalarField
 }
 
 
-Foam::radiation::wideBandDiffusiveRadiationMixedFvPatchScalarField::
+wideBandDiffusiveRadiationMixedFvPatchScalarField::
 wideBandDiffusiveRadiationMixedFvPatchScalarField
 (
     const wideBandDiffusiveRadiationMixedFvPatchScalarField& ptf,
@@ -63,11 +115,12 @@ wideBandDiffusiveRadiationMixedFvPatchScalarField
 :
     mixedFvPatchScalarField(ptf, p, iF, mapper),
     TName_(ptf.TName_),
-    emissivity_(ptf.emissivity_)
+    emissivity_(ptf.emissivity_),
+    sumOutgoingAnglesPtr_(NULL)
 {}
 
 
-Foam::radiation::wideBandDiffusiveRadiationMixedFvPatchScalarField::
+wideBandDiffusiveRadiationMixedFvPatchScalarField::
 wideBandDiffusiveRadiationMixedFvPatchScalarField
 (
     const fvPatch& p,
@@ -77,7 +130,8 @@ wideBandDiffusiveRadiationMixedFvPatchScalarField
 :
     mixedFvPatchScalarField(p, iF),
     TName_(dict.lookup("T")),
-    emissivity_(readScalar(dict.lookup("emissivity")))
+    emissivity_(readScalar(dict.lookup("emissivity"))),
+    sumOutgoingAnglesPtr_(NULL)
 {
     if (dict.found("value"))
     {
@@ -95,7 +149,7 @@ wideBandDiffusiveRadiationMixedFvPatchScalarField
     {
         // No value given. Restart as fixedValue b.c.
 
-        // Bugfix: Do not initialize from temperautre because it is unavailable
+        // Bugfix: Do not initialize from temperature because it is unavailable
         // when running, e.g. decomposePar and loading radiation as
         // shared library. Initialize to zero instead.
         // 26 Mar 2014 - DC
@@ -109,7 +163,7 @@ wideBandDiffusiveRadiationMixedFvPatchScalarField
 }
 
 
-Foam::radiation::wideBandDiffusiveRadiationMixedFvPatchScalarField::
+wideBandDiffusiveRadiationMixedFvPatchScalarField::
 wideBandDiffusiveRadiationMixedFvPatchScalarField
 (
     const wideBandDiffusiveRadiationMixedFvPatchScalarField& ptf
@@ -117,11 +171,12 @@ wideBandDiffusiveRadiationMixedFvPatchScalarField
 :
     mixedFvPatchScalarField(ptf),
     TName_(ptf.TName_),
-    emissivity_(ptf.emissivity_)
+    emissivity_(ptf.emissivity_),
+    sumOutgoingAnglesPtr_(NULL)
 {}
 
 
-Foam::radiation::wideBandDiffusiveRadiationMixedFvPatchScalarField::
+wideBandDiffusiveRadiationMixedFvPatchScalarField::
 wideBandDiffusiveRadiationMixedFvPatchScalarField
 (
     const wideBandDiffusiveRadiationMixedFvPatchScalarField& ptf,
@@ -130,97 +185,85 @@ wideBandDiffusiveRadiationMixedFvPatchScalarField
 :
     mixedFvPatchScalarField(ptf, iF),
     TName_(ptf.TName_),
-    emissivity_(ptf.emissivity_)
+    emissivity_(ptf.emissivity_),
+    sumOutgoingAnglesPtr_(NULL)
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::radiation::wideBandDiffusiveRadiationMixedFvPatchScalarField::
+void wideBandDiffusiveRadiationMixedFvPatchScalarField::updateCoeffs()
-updateCoeffs()
 {
     if (this->updated())
     {
         return;
     }
 
+    const label patchI = this->patch().index();
+
+    // Access radiation model
     const radiationModel& radiation =
         db().lookupObject<radiationModel>("radiationProperties");
 
-    const fvDOM& dom(refCast<const fvDOM>(radiation));
+    const fvDOM& dom = dynamic_cast<const fvDOM&>(radiation);
-
-    label rayId = -1;
-    label lambdaId = -1;
-    dom.setRayIdLambdaId(dimensionedInternalField().name(), rayId, lambdaId);
-
-    const label patchI = patch().index();
 
     if (dom.nLambda() == 0)
     {
         FatalErrorIn
         (
-            "Foam::radiation::"
+            ""
             "wideBandDiffusiveRadiationMixedFvPatchScalarField::updateCoeffs"
         )   << " a non-grey boundary condition is used with a grey "
             << "absorption model" << nl << exit(FatalError);
     }
 
-    scalarField& Iw = *this;
+    // Get rayId and lambda Id for this ray
-    vectorField n = patch().Sf()/patch().magSf();
+    label rayId = -1;
+    label lambdaId = -1;
+    dom.setRayIdLambdaId(dimensionedInternalField().name(), rayId, lambdaId);
 
-    radiativeIntensityRay& ray =
+    // Make shortcut to ray belonging to this field
-        const_cast<radiativeIntensityRay&>(dom.IRay(rayId));
+    const radiativeIntensityRay& ray = dom.IRay(rayId);
 
-    ray.Qr().boundaryField()[patchI] += Iw*(n & ray.dAve());
+    // Access incoming radiation for this patch for this band
+    const scalarField& Qin = dom.Qin(lambdaId)[patchI];
 
-    const scalarField Eb =
+    // Access black body radiation for this patch for this band
+    const scalarField& Eb =
         dom.blackBody().bLambda(lambdaId).boundaryField()[patchI];
 
-    forAll(Iw, faceI)
+    // Get face normals
-    {
+    vectorField nHat = patch().nf();
-        scalar Ir = 0.0;
-        for (label rayI=0; rayI < dom.nRay(); rayI++)
-        {
-            const vector& d = dom.IRay(rayI).d();
 
-            const scalarField& IFace =
+    // Calculate cos of incoming angle of current ray with every face
-                dom.IRay(rayI).ILambda(lambdaId).boundaryField()[patchI];
+    scalarField incomingAngle = this->patch().nf() & ray.dAve();
 
-            if ((-n[faceI] & d) < 0.0) // qin into the wall
+    // Set to zeroGradient (=0; incomingAngle > 0) for faces with incoming rays
-            {
+    // and to fixedValue (=1; incomingAngle < 0) for outgoing rays
-                const vector& dAve = dom.IRay(rayI).dAve();
+    this->valueFraction() = neg(incomingAngle);
-                Ir = Ir + IFace[faceI]*mag(n[faceI] & dAve);
-            }
-        }
 
-        const vector& d = dom.IRay(rayId).d();
+    // Set intensity value for fixedValue part (see reference in header file)
+    this->refValue() =
+        emissivity_*Eb + ((1 - emissivity_)*Qin)/sumOutgoingAngles();
 
-        if ((-n[faceI] & d) > 0.0)
+    // Update boundary field now, so values for incoming and outgoing rays
-        {
+    // are in balance
-            // direction out of the wall
+    scalarField::operator=
-            refGrad()[faceI] = 0.0;
-            valueFraction()[faceI] = 1.0;
-            refValue()[faceI] =
     (
-                    Ir*(1.0 - emissivity_)
+        this->valueFraction()*this->refValue()
-                  + emissivity_*Eb[faceI]
+      +
+        (1.0 - this->valueFraction())*
+        (
+            this->patchInternalField()
+          + this->refGrad()/this->patch().deltaCoeffs()
         )
-               /mathematicalConstant::pi;
+    );
-        }
-        else
-        {
-            // direction into the wall
-            valueFraction()[faceI] = 0.0;
-            refGrad()[faceI] = 0.0;
-            refValue()[faceI] = 0.0; //not used
-        }
-    }
 
     mixedFvPatchScalarField::updateCoeffs();
 }
 
 
-void Foam::radiation::wideBandDiffusiveRadiationMixedFvPatchScalarField::write
+void wideBandDiffusiveRadiationMixedFvPatchScalarField::write
 (
     Ostream& os
 ) const
@@ -234,17 +277,14 @@ void Foam::radiation::wideBandDiffusiveRadiationMixedFvPatchScalarField::write
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-namespace Foam
-{
-namespace radiation
-{
 makePatchTypeField
 (
     fvPatchScalarField,
     wideBandDiffusiveRadiationMixedFvPatchScalarField
 );
-}
-}
 
 
+} // End namespace radiation
+} // End namespace Foam
+
 // ************************************************************************* //

src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandDiffusiveRadiation/wideBandDiffusiveRadiationMixedFvPatchScalarField.H

@@ -25,7 +25,9 @@ Class
     Foam::wideBandDiffusiveRadiationMixedFvPatchScalarField
 
 Description
-    Radiation temperature specified
+    Radiation boundary will diffusive reflection and black body emission.
+    See Modest "Radiative heat transfer", Capter 16.6 "The finite volume method"
+    Eq.16.63
 
 SourceFiles
     wideBandDiffusiveRadiationMixedFvPatchScalarField.C
@@ -59,6 +61,15 @@ class wideBandDiffusiveRadiationMixedFvPatchScalarField
         //- Emissivity
         scalar emissivity_;
 
+        //- Sum of cosine factors of all outgoing rays (lazy evaluated)
+        mutable scalarField *sumOutgoingAnglesPtr_;
+
+
+    // Private member functions
+
+        //- Calculate sum of cosine factors of all outgoing rays
+        void calcSumOutgoingAngles() const;
+
 
 public:
 
@@ -157,6 +168,9 @@ public:
                 return emissivity_;
             }
 
+            //- Return sum of cosine factors of all outgoing rays
+            const scalarField& sumOutgoingAngles() const;
+
 
         // Evaluation functions
 

src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandSpecularRadiation/wideBandSpecularRadiationMixedFvPatchScalarField.C

@@ -0,0 +1,416 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     |
+    \\  /    A nd           | For copyright notice see file Copyright
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation, either version 3 of the License, or (at your
+    option) any later version.
+
+    foam-extend is distributed in the hope that it will be useful, but
+    WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+    General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "radiationModel.H"
+#include "fvDOM.H"
+#include "wideBandSpecularRadiationMixedFvPatchScalarField.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace radiation
+{
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void wideBandSpecularRadiationMixedFvPatchScalarField::
+calcReceivedRayIDs() const
+{
+    if (receivedRayIDPtr_)
+    {
+        FatalErrorIn
+        (
+            "wideBandSpecularRadiationMixedFvPatchScalarField"
+            "::calcreceivedRayIDs()"
+        )   << "receivedRayIDPtr already calculated"
+            << abort(FatalError);
+    }
+
+    receivedRayIDPtr_ = new labelListList(this->size());
+    labelListList& receivedRayID = *receivedRayIDPtr_;
+
+    // Get access to radiation model, and recast as fvDOM
+    const radiationModel& radiation =
+        this->db().lookupObject<radiationModel>("radiationProperties");
+
+    const fvDOM& dom(refCast<const fvDOM>(radiation));
+
+    // Get rayId and lambda Id for this ray
+    label rayId = -1;
+    label lambdaId = -1;
+    dom.setRayIdLambdaId(dimensionedInternalField().name(), rayId, lambdaId);
+
+    // Get index of active ray
+    word rayAndBand = this->dimensionedInternalField().name();
+    rayAndBand =
+        rayAndBand.substr(rayAndBand.find("_")+1, rayAndBand.size()-1);
+
+    // Get all ray directions
+    List<vector> dAve(dom.nRay());
+    forAll (dAve, rayI)
+    {
+        dAve[rayI] = dom.IRay(rayI).dAve()/mag(dom.IRay(rayI).dAve());
+    }
+
+    // Get face normal vectors
+    vectorField nHat = this->patch().nf();
+
+
+    // For each face, and for each reflected ray, try to find another
+    // ray that is better suited to pick it up.
+    forAll(receivedRayID, faceI)
+    {
+        // Access the list of received rays for this face
+        labelList& receivedRayIDs = receivedRayID[faceI];
+
+        // Check whether ray is going into surface
+        // -> no reflection
+        if ((dAve[rayId] & nHat[faceI]) > 0)
+        {
+            receivedRayIDs.setSize(0);
+            continue;
+        }
+
+
+        // For each face, initialize list of picked up
+        // rays to maximum possible size
+        receivedRayIDs.setSize(dom.nRay() - 1);
+
+        // Count actual number of receiving rays
+        label nReceiving = 0;
+
+        forAll(dAve, incomingRayI)
+        {
+            // Calculate reflected ray direction for rayI
+            vector dReflected =
+                dAve[incomingRayI]
+              - 2*(dAve[incomingRayI] & nHat[faceI])*nHat[faceI];
+
+            // Get dot product with this ray
+            scalar dotProductThisRay = dAve[rayId] & dReflected;
+
+            // Assume this ray is closest to the reflected ray
+            bool closest = true;
+
+            // And look through all other rays to find a better suited one
+            forAll(dAve, receivingRayI)
+            {
+                if (receivingRayI == rayId)
+                {
+                    // Do not compare this ray with itself
+                    continue;
+                }
+
+                scalar dotProductOtherRay = dAve[receivingRayI] & dReflected;
+                if (dotProductThisRay < dotProductOtherRay)
+                {
+                    // If another ray is closer, stop search
+                    closest = false;
+                    break;
+                }
+            }
+
+            if (closest)
+            {
+                // Could not find better suited ray, so this ray needs to
+                // pick it up. Add incoming ray to list for this face
+                receivedRayIDs[nReceiving] = incomingRayI;
+                nReceiving ++;
+            }
+
+        }
+
+        // Resize list of picked up rays for this face
+        receivedRayIDs.setSize(nReceiving);
+    }
+
+    // Sanity check on last ray
+    if (rayId == (dom.nRay() - 1))
+    {
+
+        forAll(nHat, faceI)
+        {
+            label incomingRays = 0;
+            label pickedUpRays = 0;
+
+            for (label rayI = 0; rayI < dom.nRay(); rayI++)
+            {
+                // Is this ray going into the wall?
+                if ((dAve[rayI] &  nHat[faceI]) > 0)
+                {
+                    incomingRays++;
+                }
+
+                // How many rays are picked up by this ray?
+                const wideBandSpecularRadiationMixedFvPatchScalarField& rayBC =
+                    refCast
+                    <
+                        const wideBandSpecularRadiationMixedFvPatchScalarField
+                    >
+                    (
+                        dom.IRayLambda
+                        (
+                            rayI,
+                            lambdaId
+                        ).boundaryField()[patch().index()]
+                    );
+
+                const labelListList& receivedRaysList = rayBC.receivedRayIDs();
+
+                pickedUpRays += receivedRaysList[faceI].size();
+            }
+
+            if (incomingRays != pickedUpRays)
+            {
+                FatalErrorIn
+                (
+                    "wideBandSpecularRadiationMixedFvPatchScalarField"
+                    "::calcreceivedRayIDs()"
+                )   << "Sanity checked failed on patch " << patch().name()
+                    << " face " << faceI << nl
+                    << "number of rays with direction into wall: "
+                    << incomingRays
+                    << ", number of reflected rays picked up by outgoing rays: "
+                    << pickedUpRays << nl
+                    << abort(FatalError);
+            }
+        }
+    }
+
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+wideBandSpecularRadiationMixedFvPatchScalarField::
+wideBandSpecularRadiationMixedFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    mixedFvPatchField<scalar>(p, iF),
+    receivedRayIDPtr_(NULL)
+{
+    this->refValue() = 0;
+    this->refGrad() = 0;
+    this->valueFraction() = 0;
+}
+
+
+wideBandSpecularRadiationMixedFvPatchScalarField::
+wideBandSpecularRadiationMixedFvPatchScalarField
+(
+    const wideBandSpecularRadiationMixedFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    mixedFvPatchField<scalar>(ptf, p, iF, mapper),
+    receivedRayIDPtr_(NULL)
+{}
+
+
+wideBandSpecularRadiationMixedFvPatchScalarField::
+wideBandSpecularRadiationMixedFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    mixedFvPatchField<scalar>(p, iF),
+    receivedRayIDPtr_(NULL)
+{
+    this->refValue() = 0;
+    this->refGrad() = 0;
+    this->valueFraction() = 0;
+
+    if (dict.found("patchType"))
+    {
+        word& pType = const_cast<word &>(this->patchType());
+        pType = word(dict.lookup("patchType"));
+}
+
+    if (dict.found("value"))
+    {
+        fvPatchField<scalar>::operator=
+        (
+            scalarField("value", dict, p.size())
+        );
+    }
+    else
+    {
+        fvPatchScalarField::operator=(this->refValue());
+    }
+}
+
+
+wideBandSpecularRadiationMixedFvPatchScalarField::
+wideBandSpecularRadiationMixedFvPatchScalarField
+(
+    const wideBandSpecularRadiationMixedFvPatchScalarField& ptf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    mixedFvPatchField<scalar>(ptf, iF),
+    receivedRayIDPtr_(NULL)
+{}
+
+
+wideBandSpecularRadiationMixedFvPatchScalarField::
+wideBandSpecularRadiationMixedFvPatchScalarField
+(
+    const wideBandSpecularRadiationMixedFvPatchScalarField& ptf
+)
+:
+    mixedFvPatchField<scalar>(ptf),
+    receivedRayIDPtr_(NULL)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::labelListList&
+wideBandSpecularRadiationMixedFvPatchScalarField::receivedRayIDs() const
+{
+    if (!receivedRayIDPtr_)
+    {
+        calcReceivedRayIDs();
+    }
+
+    return *receivedRayIDPtr_;
+}
+
+
+// Evaluate the field on the patch
+void wideBandSpecularRadiationMixedFvPatchScalarField::updateCoeffs()
+{
+    if (this->updated())
+    {
+        return;
+    }
+
+    vectorField nHat = this->patch().nf();
+
+    // Get access to radiation model, and recast as fvDOM
+    const radiationModel& radiation =
+        this->db().lookupObject<radiationModel>("radiationProperties");
+
+    const fvDOM& dom(refCast<const fvDOM>(radiation));
+
+    // Get rayId and lambda Id for this ray
+    label rayId = -1;
+    label lambdaId = -1;
+    dom.setRayIdLambdaId(dimensionedInternalField().name(), rayId, lambdaId);
+
+    const labelListList& receivedRayIDs = this->receivedRayIDs();
+
+    scalarField IValue(this->size(), 0);
+
+    labelList faceCellIDs = this->patch().faceCells();
+
+    // Loop over all faces and add values from all received rays
+    forAll(receivedRayIDs, faceI)
+    {
+        if (receivedRayIDs[faceI].size() == 0)
+        {
+            // Ray goes into face -> act as zeroGradient
+            this->valueFraction()[faceI] = 0;
+        }
+        else
+        {
+            // Ray goes out of face -> act as fixedValue
+            this->valueFraction()[faceI] = 1;
+
+            // Pick up all reflected rays
+            forAll(receivedRayIDs[faceI], receivedRayI)
+            {
+                // Get ray from object registry
+                IValue[faceI] +=
+                    dom.IRayLambda
+                    (
+                        receivedRayIDs[faceI][receivedRayI],
+                        lambdaId
+                    ).internalField()[faceCellIDs[faceI]];
+            }
+        }
+    }
+
+    // Set value for patch
+    this->refValue() = IValue;
+
+    scalarField::operator=
+    (
+        this->valueFraction()*this->refValue()
+      +
+        (1.0 - this->valueFraction())*
+        (
+            this->patchInternalField()
+          + this->refGrad()/this->patch().deltaCoeffs()
+        )
+    );
+
+    mixedFvPatchField<scalar>::updateCoeffs();
+}
+
+
+void wideBandSpecularRadiationMixedFvPatchScalarField::write(Ostream& os) const
+{
+    fvPatchScalarField::write(os);
+    this->writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+void wideBandSpecularRadiationMixedFvPatchScalarField::operator=
+(
+    const fvPatchScalarField& ptf
+)
+{
+    fvPatchScalarField::operator=
+    (
+        this->valueFraction()*this->refValue()
+        + (1 - this->valueFraction())*ptf
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeField
+(
+   fvPatchScalarField,
+   wideBandSpecularRadiationMixedFvPatchScalarField
+);
+
+
+} // End namespace radiation
+} // End namespace Foam
+
+// ************************************************************************* //

src/thermophysicalModels/radiation/derivedFvPatchFields/wideBandSpecularRadiation/wideBandSpecularRadiationMixedFvPatchScalarField.H

@@ -0,0 +1,162 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     |
+    \\  /    A nd           | For copyright notice see file Copyright
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation, either version 3 of the License, or (at your
+    option) any later version.
+
+    foam-extend is distributed in the hope that it will be useful, but
+    WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+    General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::wideBandSpecularRadiationMixedFvPatchScalarField
+
+Description
+    Foam::wideBandSpecularRadiationMixedFvPatchScalarField
+
+SourceFiles
+    wideBandSpecularRadiationMixedFvPatchScalarField.C
+
+Author
+    Dominik Christ, Wikki Ltd. All rights reserved.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef wideBandSpecularRadiationMixedFvPatchScalarField_H
+#define wideBandSpecularRadiationMixedFvPatchScalarField_H
+
+#include "mixedFvPatchField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace radiation
+{
+
+/*---------------------------------------------------------------------------*\
+            Class wideBandSpecularRadiationMixedFvPatchScalarField Declaration
+\*---------------------------------------------------------------------------*/
+
+class wideBandSpecularRadiationMixedFvPatchScalarField
+:
+    public mixedFvPatchField<scalar>
+{
+
+    //- For each face, ray index of ray in direction of reflection
+    mutable labelListList *receivedRayIDPtr_;
+
+    //- Calculate IDs of rays in direction of reflection
+    void calcReceivedRayIDs() const;
+
+public:
+
+    //- Runtime type information
+    TypeName("wideBandSpecularRadiation");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        wideBandSpecularRadiationMixedFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        wideBandSpecularRadiationMixedFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given wideBandSpecularRadiationMixedFvPatchScalarField onto a new patch
+        wideBandSpecularRadiationMixedFvPatchScalarField
+        (
+            const wideBandSpecularRadiationMixedFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        wideBandSpecularRadiationMixedFvPatchScalarField
+        (
+            const wideBandSpecularRadiationMixedFvPatchScalarField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchField<scalar> > clone() const
+        {
+            return tmp<fvPatchField<scalar> >
+            (
+                new wideBandSpecularRadiationMixedFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        wideBandSpecularRadiationMixedFvPatchScalarField
+        (
+            const wideBandSpecularRadiationMixedFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchField<scalar> > clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchField<scalar> >
+            (
+                new wideBandSpecularRadiationMixedFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        // Access
+
+            //- Return ray IDs of reflected rays
+            const labelListList& receivedRayIDs() const;
+
+        //- Update the coefficients associated with the patch field
+        virtual void updateCoeffs();
+
+        //- Write
+        virtual void write(Ostream&) const;
+
+
+    // Member operators
+
+        virtual void operator=(const fvPatchScalarField& pvf);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace radiation
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/radiation/radiationModel/fvDOM/blackBodyEmission/blackBodyEmission.C

@@ -202,7 +202,7 @@ Foam::radiation::blackBodyEmission::EbDeltaLambdaT
     const Vector2D<scalar>& band
 ) const
 {
-    tmp<volScalarField> Eb
+    tmp<volScalarField> tEb
     (
         new volScalarField
         (
@@ -217,30 +217,34 @@ Foam::radiation::blackBodyEmission::EbDeltaLambdaT
             radiation::sigmaSB*pow4(T)
         )
     );
+    volScalarField& Eb = tEb();
 
+    if (magSqr(band - Vector2D<scalar>::one) > SMALL) // Multiple bands?
+    {
+        forAll(Eb.internalField(), cellI)
+        {
+            scalar percentileUpper = fLambdaT(band[1]*T.internalField()[cellI]);
+            scalar percentileLower = fLambdaT(band[0]*T.internalField()[cellI]);
+            Eb.internalField()[cellI] *= (percentileUpper - percentileLower);
+        }
 
-    if (band == Vector2D<scalar>::one)
+        forAll(Eb.boundaryField(), patchI)
         {
-        return Eb;
+            forAll(Eb.boundaryField()[patchI], faceI)
-    }
-    else
             {
-        forAll(T, i)
+                scalar percentileUpper =
-        {
+                    fLambdaT(band[1]*T.boundaryField()[patchI][faceI]);
-            scalar T1 = fLambdaT(band[1]*T[i]);
+                scalar percentileLower =
-            scalar T2 = fLambdaT(band[0]*T[i]);
+                    fLambdaT(band[0]*T.boundaryField()[patchI][faceI]);
-            dimensionedScalar fLambdaDelta
+                Eb.boundaryField()[patchI][faceI] *=
-            (
+                    (percentileUpper - percentileLower);
-                "fLambdaDelta",
-                dimless,
-                T1 - T2
-            );
-            Eb()[i] = Eb()[i]*fLambdaDelta.value();
             }
-        return Eb;
         }
     }
 
+    return tEb;
+}
+
 
 void Foam::radiation::blackBodyEmission::correct
 (

src/thermophysicalModels/radiation/radiationModel/fvDOM/fvDOM/fvDOM.C

@@ -25,10 +25,12 @@ License
 
 #include "fvDOM.H"
 #include "addToRunTimeSelectionTable.H"
-
+#include "fvm.H"
 #include "mathematicalConstants.H"
 #include "radiationConstants.H"
 
+using namespace Foam::mathematicalConstant;
+
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
@@ -57,66 +59,27 @@ Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
         IOobject
         (
             "G",
-            mesh_.time().timeName(),
+            mesh().time().timeName(),
             T.db(),
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        mesh_,
+        mesh(),
         dimensionedScalar("G", dimMass/pow3(dimTime), 0.0)
     ),
-    Qr_
-    (
-        IOobject
-        (
-            "Qr",
-            mesh_.time().timeName(),
-            T.db(),
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("Qr", dimMass/pow3(dimTime), 0.0)
-    ),
     a_
     (
         IOobject
         (
             "a",
-            mesh_.time().timeName(),
+            mesh().time().timeName(),
             T.db(),
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        mesh_,
+        mesh(),
         dimensionedScalar("a", dimless/dimLength, 0.0)
     ),
-    e_
-    (
-        IOobject
-        (
-            "e",
-            mesh_.time().timeName(),
-            T.db(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("a", dimless/dimLength, 0.0)
-    ),
-    E_
-    (
-        IOobject
-        (
-            "E",
-            mesh_.time().timeName(),
-            T.db(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("E", dimMass/dimLength/pow3(dimTime), 0.0)
-    ),
     nTheta_(readLabel(coeffs_.lookup("nTheta"))),
     nPhi_(readLabel(coeffs_.lookup("nPhi"))),
     nRay_(0),
@@ -124,15 +87,20 @@ Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
     aLambda_(nLambda_),
     blackBody_(nLambda_, T),
     IRay_(0),
+    Qem_(nLambda_),
+    Qin_(nLambda_),
     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)
 {
-    if (mesh_.nSolutionD() == 3)    //3D
+    if (mesh().nSolutionD() == 3)    //3D
     {
         nRay_ = 4*nPhi_*nTheta_;
         IRay_.setSize(nRay_);
-        scalar deltaPhi = mathematicalConstant::pi/(2.0*nPhi_);
+        scalar deltaPhi = pi/(2.0*nPhi_);
-        scalar deltaTheta = mathematicalConstant::pi/nTheta_;
+        scalar deltaTheta = pi/nTheta_;
         label i = 0;
         for (label n = 1; n <= nTheta_; n++)
         {
@@ -146,14 +114,15 @@ Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
                     new radiativeIntensityRay
                     (
                         *this,
-                        mesh_,
+                        mesh(),
                         phii,
                         thetai,
                         deltaPhi,
                         deltaTheta,
                         nLambda_,
                         absorptionEmission_,
-                        blackBody_
+                        blackBody_,
+                        i
                     )
                 );
                 i++;
@@ -162,13 +131,13 @@ Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
     }
     else
     {
-        if (mesh_.nSolutionD() == 2)    //2D (X & Y)
+        if (mesh().nSolutionD() == 2)    //2D (X & Y)
         {
-            scalar thetai = mathematicalConstant::piByTwo;
+            scalar thetai = piByTwo;
-            scalar deltaTheta = mathematicalConstant::pi;
+            scalar deltaTheta = pi;
             nRay_ = 4*nPhi_;
             IRay_.setSize(nRay_);
-            scalar deltaPhi = mathematicalConstant::pi/(2.0*nPhi_);
+            scalar deltaPhi = pi/(2.0*nPhi_);
             label i = 0;
             for (label m = 1; m <= 4*nPhi_; m++)
             {
@@ -179,14 +148,15 @@ Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
                     new radiativeIntensityRay
                     (
                         *this,
-                        mesh_,
+                        mesh(),
                         phii,
                         thetai,
                         deltaPhi,
                         deltaTheta,
                         nLambda_,
                         absorptionEmission_,
-                        blackBody_
+                        blackBody_,
+                        i
                     )
                 );
                 i++;
@@ -194,11 +164,11 @@ Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
         }
         else    //1D (X)
         {
-            scalar thetai = mathematicalConstant::piByTwo;
+            scalar thetai = piByTwo;
-            scalar deltaTheta = mathematicalConstant::pi;
+            scalar deltaTheta = pi;
             nRay_ = 2;
             IRay_.setSize(nRay_);
-            scalar deltaPhi = mathematicalConstant::pi;
+            scalar deltaPhi = pi;
             label i = 0;
             for (label m = 1; m <= 2; m++)
             {
@@ -209,19 +179,19 @@ Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
                     new radiativeIntensityRay
                     (
                         *this,
-                        mesh_,
+                        mesh(),
                         phii,
                         thetai,
                         deltaPhi,
                         deltaTheta,
                         nLambda_,
                         absorptionEmission_,
-                        blackBody_
+                        blackBody_,
+                        i
                     )
                 );
                 i++;
             }
-
         }
     }
 
@@ -237,7 +207,7 @@ Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
                 IOobject
                 (
                     "aLambda_" + Foam::name(lambdaI) ,
-                    mesh_.time().timeName(),
+                    mesh().time().timeName(),
                     T.db(),
                     IOobject::NO_READ,
                     IOobject::NO_WRITE
@@ -245,6 +215,30 @@ Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
                 a_
             )
         );
+
+        Qin_.set
+        (
+            lambdaI,
+            volScalarField::GeometricBoundaryField
+            (
+                mesh().boundary(),
+                mesh().V(),           // Dummy internal field,
+                calculatedFvPatchScalarField::typeName
+            )
+        );
+        Qin_[lambdaI] = 0;
+
+        Qem_.set
+        (
+            lambdaI,
+            volScalarField::GeometricBoundaryField
+            (
+                mesh().boundary(),
+                mesh().V(),           // Dummy internal field,
+                calculatedFvPatchScalarField::typeName
+            )
+        );
+        Qem_[lambdaI] = 0;
     }
 
     Info<< "fvDOM : Allocated " << IRay_.size()
@@ -254,6 +248,43 @@ Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
         Info<< '\t' << IRay_[i].I().name()
             << '\t' << IRay_[i].dAve() << 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());
+                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;
 }
 
@@ -266,12 +297,108 @@ Foam::radiation::fvDOM::~fvDOM()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+Foam::tmp<Foam::volScalarField::GeometricBoundaryField>
+Foam::radiation::fvDOM::Qin() const
+{
+    tmp<volScalarField::GeometricBoundaryField> tQin
+    (
+        new volScalarField::GeometricBoundaryField
+        (
+            mesh().boundary(),
+            mesh().V(),           // Dummy internal field,
+            calculatedFvPatchScalarField::typeName
+        )
+    );
+    volScalarField::GeometricBoundaryField& sumQin = tQin();
+
+    sumQin = 0;
+
+    forAll(Qin_, lambdaI)
+    {
+        sumQin += Qin(lambdaI);
+    }
+
+    return tQin;
+}
+
+
+Foam::tmp<Foam::volScalarField::GeometricBoundaryField>
+Foam::radiation::fvDOM::Qem() const
+{
+    tmp<volScalarField::GeometricBoundaryField> tsumQem
+    (
+        new volScalarField::GeometricBoundaryField
+        (
+            mesh().boundary(),
+            mesh().V(),           // Dummy internal field,
+            calculatedFvPatchScalarField::typeName
+        )
+    );
+    volScalarField::GeometricBoundaryField& sumQem = tsumQem();
+
+    sumQem = 0;
+
+    forAll(Qem_, lambdaI)
+    {
+        sumQem += Qem(lambdaI);
+    }
+
+    return tsumQem;
+}
+
+Foam::tmp<Foam::volScalarField::GeometricBoundaryField>
+Foam::radiation::fvDOM::Qr() const
+{
+    tmp<volScalarField::GeometricBoundaryField> tQr
+    (
+        new volScalarField::GeometricBoundaryField
+        (
+            mesh().boundary(),
+            mesh().V(),           // Dummy internal field,
+            calculatedFvPatchScalarField::typeName
+        )
+    );
+    volScalarField::GeometricBoundaryField& Qr = tQr();
+
+    Qr = 0;
+
+    forAll(Qem_, lambdaI)
+    {
+        Qr += Qin(lambdaI);
+        Qr -= Qem(lambdaI);
+    }
+
+    return tQr;
+}
+
+
+Foam::tmp<Foam::volScalarField::GeometricBoundaryField>
+Foam::radiation::fvDOM::Qr(scalar lambdaI) const
+{
+    tmp<volScalarField::GeometricBoundaryField> tQr
+    (
+        new volScalarField::GeometricBoundaryField
+        (
+            mesh().boundary(),
+            mesh().V(),           // Dummy internal field,
+            calculatedFvPatchScalarField::typeName
+        )
+    );
+    volScalarField::GeometricBoundaryField& Qr = tQr();
+
+    Qr = 0;
+    Qr += Qin(lambdaI);
+    Qr -= Qem(lambdaI);
+
+    return tQr;
+}
+
+
 bool Foam::radiation::fvDOM::read()
 {
     if (radiationModel::read())
     {
         // Only reading solution parameters - not changing ray geometry
-
         coeffs_.readIfPresent("convergence", convergence_);
         coeffs_.readIfPresent("maxIter", maxIter_);
 
@@ -290,19 +417,96 @@ void Foam::radiation::fvDOM::calculate()
 
     updateBlackBodyEmission();
 
-    scalar maxResidual = 0.0;
+    scalar maxResidual = 0;
     label radIter = 0;
     do
     {
-        radIter++;
+        Info << "Radiation solver iter: " << radIter << endl;
+
+        Info << "Updating Radiation BCs..." << flush;
         forAll(IRay_, rayI)
         {
-            maxResidual = 0.0;
+            IRay_[rayI].updateBCs();
-            scalar maxBandResidual = IRay_[rayI].correct();
+        }
-            maxResidual = max(maxBandResidual, maxResidual);
+        Info << "done." << endl;
+
+        // For debug purposes, recalculate Qin and Qem to see if balances match
+        if (debug)
+        {
+            // Update radiation balances
+            forAll(aLambda_, lambdaI)
+            {
+                Qem_[lambdaI] = 0;
+                Qin_[lambdaI] = 0;
+
+                forAll(Qem_[lambdaI], patchI)
+                {
+                    // Loop over all rays
+                    forAll(IRay_, rayI)
+                    {
+                        const fvPatchScalarField& curPatch =
+                            IRay_[rayI].ILambda
+                            (
+                                lambdaI
+                            ).boundaryField()[patchI];
+
+                        scalarField incomingAngle =
+                            curPatch.patch().nf() & IRay_[rayI].dAve();
+
+                        Qin_[lambdaI][patchI] +=
+                            pos(incomingAngle)*curPatch*incomingAngle;
+
+                        Qem_[lambdaI][patchI] +=
+                            neg(incomingAngle)*curPatch*(-incomingAngle);
                     }
 
-        Info << "Radiation solver iter: " << radIter << endl;
+                    Info << "Patch " << Qem_[lambdaI][patchI].patch().name()
+                        << " band " << lambdaI
+                        << ": Radiation incoming "
+                        << sum(Qin_[lambdaI][patchI])
+                        << "  outgoing "
+                        << sum(Qem_[lambdaI][patchI])
+                        << endl;
+                }
+            }
+        }
+
+        // Solve ray transport equations
+        forAll(IRay_, rayI)
+        {
+            maxResidual = 0;
+            scalar maxRayResidual = IRay_[rayI].correct();
+            maxResidual = max(maxRayResidual, maxResidual);
+        }
+
+        // Update radiation balances
+        forAll(aLambda_, lambdaI)
+        {
+            Qem_[lambdaI] = 0;
+            Qin_[lambdaI] = 0;
+
+            forAll(Qem_[lambdaI], patchI)
+            {
+                // Loop over all rays
+                forAll(IRay_, rayI)
+                {
+                    const fvPatchScalarField& curPatch =
+                        IRay_[rayI].ILambda(lambdaI).boundaryField()[patchI];
+
+                    scalarField incomingAngle =
+                        curPatch.patch().nf() & IRay_[rayI].dAve();
+
+                    Qin_[lambdaI][patchI] +=
+                        pos(incomingAngle)*curPatch*incomingAngle;
+
+                    Qem_[lambdaI][patchI] +=
+                        neg(incomingAngle)*curPatch*(-incomingAngle);
+                }
+            }
+        }
+
+        radIter++;
+        Info << "Max residual: " << maxResidual << endl;
 
     } while(maxResidual > convergence_ && radIter < maxIter_);
 
@@ -319,8 +523,8 @@ Foam::tmp<Foam::volScalarField> Foam::radiation::fvDOM::Rp() const
             IOobject
             (
                 "Rp",
-                mesh_.time().timeName(),
+                mesh().time().timeName(),
-                mesh_,
+                mesh(),
                 IOobject::NO_READ,
                 IOobject::NO_WRITE,
                 false
@@ -340,7 +544,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;
 }
@@ -358,14 +562,11 @@ void Foam::radiation::fvDOM::updateBlackBodyEmission()
 void Foam::radiation::fvDOM::updateG()
 {
     G_ = dimensionedScalar("zero",dimMass/pow3(dimTime), 0.0);
-    Qr_ = dimensionedScalar("zero",dimMass/pow3(dimTime), 0.0);
 
     forAll(IRay_, rayI)
     {
         IRay_[rayI].addIntensity();
         G_ += IRay_[rayI].I()*IRay_[rayI].omega();
-        //Qr_ += IRay_[rayI].Qr();
-        Qr_.boundaryField() += IRay_[rayI].Qr().boundaryField();
     }
 }
 

src/thermophysicalModels/radiation/radiationModel/fvDOM/fvDOM/fvDOM.H

@@ -80,18 +80,9 @@ class fvDOM
         //- Incident radiation  [W/m2]
         volScalarField G_;
 
-        //- Total radiative heat flux [W/m2]
-        volScalarField Qr_;
-
         //- Total absorption coefficient [1/m]
         volScalarField a_;
 
-        //- Total emission coefficient [1/m]
-        volScalarField e_;
-
-        //- Emission contribution [Kg/m/s^3]
-        volScalarField E_;
-
         //- Number of solid angles in theta
         label nTheta_;
 
@@ -113,12 +104,27 @@ class fvDOM
         //- List of pointers to radiative intensity rays
         PtrList<radiativeIntensityRay> IRay_;
 
+        //- Emmited radiative heat flux [W/m2], per band
+        PtrList<volScalarField::GeometricBoundaryField> Qem_;
+
+        //- Incidet radiative heat flux [W/m2], per band
+        PtrList<volScalarField::GeometricBoundaryField> Qin_;
+
         //- Convergence criterion
         scalar convergence_;
 
         //- 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
 
@@ -150,7 +156,17 @@ public:
 
     // Member functions
 
-        // Edit
+        //- Total incident radiative heat flux field
+        tmp<volScalarField::GeometricBoundaryField> Qin() const;
+
+        //- Total emitted radiative heat flux field
+        tmp<volScalarField::GeometricBoundaryField> Qem() const;
+
+        //- Total radiative heat flux as difference of Qem and Qin
+        tmp<volScalarField::GeometricBoundaryField> Qr() const;
+
+        //- Band-wise radiative heat flux as difference of Qem and Qin
+        tmp<volScalarField::GeometricBoundaryField> Qr(scalar lambda) const;
 
         //- Solve radiation equation(s)
         void calculate();
@@ -210,11 +226,33 @@ public:
             //- Const access to incident radiation field
             inline const volScalarField& G() const;
 
-            //- Const access to total radiative heat flux field
+            //- Const access to band wise emitted radiative heat flux field
-            inline const volScalarField& Qr() const;
+            inline const volScalarField::GeometricBoundaryField& Qem
+            (
+                scalar lambda
+            ) const;
+
+            //- Const access to band wise incident radiative heat flux field
+            inline const volScalarField::GeometricBoundaryField& Qin
+            (
+                scalar lambda
+            ) const;
 
             //- 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;
 };
 
 
@@ -222,6 +260,7 @@ public:
 
 #include "fvDOMI.H"
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace radiation

src/thermophysicalModels/radiation/radiationModel/fvDOM/fvDOM/fvDOMI.H

@@ -29,7 +29,6 @@ Foam::radiation::fvDOM::IRay(const label rayI) const
     return  IRay_[rayI];
 }
 
-
 inline const Foam::volScalarField&
 Foam::radiation::fvDOM::IRayLambda
 (
@@ -86,12 +85,34 @@ inline const Foam::volScalarField& Foam::radiation::fvDOM::G() const
 }
 
 
-inline const Foam::volScalarField& Foam::radiation::fvDOM::Qr() const
+inline const Foam::volScalarField::GeometricBoundaryField&
+Foam::radiation::fvDOM::Qin(scalar lambdaI) const
 {
-    return Qr_;
+    return Qin_[lambdaI];
 }
 
 
+/*
+inline Foam::volScalarField::GeometricBoundaryField&
+Foam::radiation::fvDOM::Qin(scalar lambdaI)
+{
+    return Qin_[lambdaI];
+}
+*/
+inline const Foam::volScalarField::GeometricBoundaryField&
+Foam::radiation::fvDOM::Qem(scalar lambdaI) const
+{
+    return Qem_[lambdaI];
+}
+
+
+/*
+inline Foam::volScalarField::GeometricBoundaryField&
+Foam::radiation::fvDOM::Qem(scalar lambdaI)
+{
+    return Qem_[lambdaI];
+}
+*/
 inline const Foam::radiation::blackBodyEmission&
 Foam::radiation::fvDOM::blackBody() const
 {
@@ -99,4 +120,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_;
+}
+
+
 // ************************************************************************* //

src/thermophysicalModels/radiation/radiationModel/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C

@@ -26,6 +26,9 @@ License
 #include "radiativeIntensityRay.H"
 #include "fvm.H"
 #include "fvDOM.H"
+#include "mathematicalConstants.H"
+
+using namespace Foam::mathematicalConstant;
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -47,7 +50,8 @@ Foam::radiation::radiativeIntensityRay::radiativeIntensityRay
     const scalar deltaTheta,
     const label nLambda,
     const absorptionEmissionModel& absorptionEmission,
-    const blackBodyEmission& blackBody
+    const blackBodyEmission& blackBody,
+    const label rayId
 )
 :
     dom_(dom),
@@ -67,26 +71,14 @@ Foam::radiation::radiativeIntensityRay::radiativeIntensityRay
         mesh_,
         dimensionedScalar("I", dimMass/pow3(dimTime), 0.0)
     ),
-    Qr_
-    (
-        IOobject
-        (
-            "Qr" + name(rayId),
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("Qr", dimMass/pow3(dimTime), 0.0)
-    ),
     d_(vector::zero),
     dAve_(vector::zero),
     theta_(theta),
     phi_(phi),
     omega_(0.0),
     nLambda_(nLambda),
-    ILambda_(nLambda)
+    ILambda_(nLambda),
+    myRayId_(rayId)
 {
     scalar sinTheta = Foam::sin(theta);
     scalar cosTheta = Foam::cos(theta);
@@ -108,6 +100,7 @@ Foam::radiation::radiativeIntensityRay::radiativeIntensityRay
         0.5*deltaPhi*Foam::sin(2.0*theta)*Foam::sin(deltaTheta)
     );
 
+//    dAve_ /= mag(dAve_);
 
     autoPtr<volScalarField> IDefaultPtr;
 
@@ -164,7 +157,6 @@ Foam::radiation::radiativeIntensityRay::radiativeIntensityRay
             );
         }
     }
-    rayId++;
 }
 
 
@@ -178,24 +170,32 @@ Foam::radiation::radiativeIntensityRay::~radiativeIntensityRay()
 
 Foam::scalar Foam::radiation::radiativeIntensityRay::correct()
 {
-    // reset boundary heat flux to zero
-    //Qr_ = dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0);
-    Qr_.boundaryField() = 0.0;
-
     scalar maxResidual = -GREAT;
 
     forAll(ILambda_, lambdaI)
     {
         const volScalarField& k = dom_.aLambda(lambdaI);
 
-        surfaceScalarField Ji = dAve_ & mesh_.Sf();
+        tmp<fvScalarMatrix> divJiILambda;
+
+        // Retrieve advection term from cache or create from scratch
+        if (!dom_.cacheDiv())
+        {
+            const surfaceScalarField Ji(dAve_ & mesh_.Sf());
+
+            divJiILambda = fvm::div(Ji, ILambda_[lambdaI], "div(Ji,Ii_h)");
+        }
+        else
+        {
+            divJiILambda = dom_.fvRayDiv(myRayId_, lambdaI);
+        }
 
         fvScalarMatrix IiEq
         (
-            fvm::div(Ji, ILambda_[lambdaI], "div(Ji,Ii_h)")
+            divJiILambda()
           + fvm::Sp(k*omega_, ILambda_[lambdaI])
         ==
-            1.0/Foam::mathematicalConstant::pi*omega_
+            1.0/pi*omega_
           * (
                 k*blackBody_.bLambda(lambdaI)
               + absorptionEmission_.ECont(lambdaI)/4
@@ -209,7 +209,6 @@ Foam::scalar Foam::radiation::radiativeIntensityRay::correct()
             IiEq,
             mesh_.solutionDict().solver("Ii")
         ).initialResidual();
-
         maxResidual = max(eqnResidual, maxResidual);
     }
 
@@ -217,6 +216,15 @@ Foam::scalar Foam::radiation::radiativeIntensityRay::correct()
 }
 
 
+void Foam::radiation::radiativeIntensityRay::updateBCs()
+{
+    forAll(ILambda_, lambdaI)
+    {
+        ILambda_[lambdaI].boundaryField().updateCoeffs();
+    }
+}
+
+
 void Foam::radiation::radiativeIntensityRay::addIntensity()
 {
     I_ = dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0);

src/thermophysicalModels/radiation/radiationModel/fvDOM/radiativeIntensityRay/radiativeIntensityRay.H

@@ -38,6 +38,7 @@ SourceFiles
 #include "absorptionEmissionModel.H"
 #include "blackBodyEmission.H"
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -78,9 +79,6 @@ private:
         //- Total radiative intensity / [W/m2]
         volScalarField I_;
 
-        //- Total radiative heat flux on boundary
-        volScalarField Qr_;
-
         //- Direction
         vector d_;
 
@@ -105,6 +103,9 @@ private:
         //- Global ray id - incremented in constructor
         static label rayId;
 
+        //- My ray Id
+        label myRayId_;
+
 
     // Private member functions
 
@@ -130,7 +131,8 @@ public:
             const scalar deltaTheta,
             const label lambda,
             const absorptionEmissionModel& absEmmModel_,
-            const blackBodyEmission& blackBody
+            const blackBodyEmission& blackBody,
+            const label rayId
         );
 
 
@@ -155,6 +157,9 @@ public:
                 const scalar lambda
             );
 
+            //- Add radiative intensities from all the bands
+            void updateBCs();
+
             //- Add radiative intensities from all the bands
             void addIntensity();
 
@@ -164,12 +169,6 @@ public:
             //- Return intensity
             inline const volScalarField& I() const;
 
-            //- Return const access to the boundary heat flux
-            inline const volScalarField& Qr() const;
-
-            //- Return non-const access to the boundary heat flux
-            inline volScalarField& Qr();
-
             //- Return direction
             inline const vector& d() const;
 
@@ -188,6 +187,9 @@ public:
             //- Return the solid angle
             inline scalar omega() const;
 
+            //- Return the radiative intensity for a given wavelength
+            inline volScalarField& ILambda(const label lambdaI);
+
             //- Return the radiative intensity for a given wavelength
             inline const volScalarField& ILambda(const label lambdaI) const;
 };

src/thermophysicalModels/radiation/radiationModel/fvDOM/radiativeIntensityRay/radiativeIntensityRayI.H

@@ -29,18 +29,6 @@ inline const Foam::volScalarField& Foam::radiation::radiativeIntensityRay::I() c
 }
 
 
-inline const Foam::volScalarField& Foam::radiation::radiativeIntensityRay::Qr() const
-{
-    return Qr_;
-}
-
-
-inline Foam::volScalarField& Foam::radiation::radiativeIntensityRay::Qr()
-{
-    return Qr_;
-}
-
-
 inline const Foam::vector& Foam::radiation::radiativeIntensityRay::d() const
 {
     return d_;
@@ -76,8 +64,17 @@ inline Foam::scalar Foam::radiation::radiativeIntensityRay::omega() const
     return omega_;
 }
 
+inline Foam::volScalarField&
+Foam::radiation::radiativeIntensityRay::ILambda
+(
+    const label lambdaI
+)
+{
+    return ILambda_[lambdaI];
+}
 
-inline const Foam::volScalarField& Foam::radiation::radiativeIntensityRay::ILambda
+inline const Foam::volScalarField&
+Foam::radiation::radiativeIntensityRay::ILambda
 (
     const label lambdaI
 ) const

src/thermophysicalModels/radiation/radiationModel/radiationModel/radiationModel.H

@@ -168,6 +168,12 @@ public:
 
         // Access
 
+            //- Const access to mesh
+            const fvMesh& mesh() const
+            {
+                return mesh_;
+            };
+
             //- Source term component (for power of T^4)
             virtual tmp<volScalarField> Rp() const = 0;
 

tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOMSpecularCeiling/0/G

@@ -0,0 +1,56 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | foam-extend: Open Source CFD                    |
+|  \\    /   O peration     | Version:     3.1                                |
+|   \\  /    A nd           | Web:         http://www.extend-project.de       |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      G;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 0 -3 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    floor
+    {
+        type            MarshakRadiation;
+        T               T;
+        emissivity      1;
+        value           uniform 0;
+    }
+
+    fixedWalls
+    {
+        type            MarshakRadiation;
+        T               T;
+        emissivity      1;
+        value           uniform 0;
+    }
+
+    ceiling
+    {
+        type            MarshakRadiation;
+        T               T;
+        emissivity      1;
+        value           uniform 0;
+    }
+
+    box
+    {
+        type            MarshakRadiation;
+        T               T;
+        emissivity      1;
+        value           uniform 0;
+    }
+}
+
+// ************************************************************************* //

tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOMSpecularCeiling/0/IDefault

@@ -0,0 +1,39 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | foam-extend: Open Source CFD                    |
+|  \\    /   O peration     | Version:     3.1                                |
+|   \\  /    A nd           | Web:         http://www.extend-project.de       |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      IDefault;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 0 -3 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    ".*"
+    {
+        type            greyDiffusiveRadiation;
+        T               T;
+        emissivity      0.5;
+        value           uniform 0;
+    }
+    "ceiling"
+    {
+        type            wideBandSpecularRadiation;
+        T               T;
+        emissivity      0.5;
+        value           uniform 0;
+    }
+}
+
+// ************************************************************************* //

tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOMSpecularCeiling/0/T

@@ -0,0 +1,47 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | foam-extend: Open Source CFD                    |
+|  \\    /   O peration     | Version:     3.1                                |
+|   \\  /    A nd           | Web:         http://www.extend-project.de       |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      T;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 1 0 0 0];
+
+internalField   uniform 300;
+
+boundaryField
+{
+    floor
+    {
+        type            fixedValue;
+        value           uniform 300.0;
+    }
+
+    ceiling
+    {
+        type            fixedValue;
+        value           uniform 300.0;
+    }
+
+    fixedWalls
+    {
+        type            zeroGradient;
+    }
+
+    box
+    {
+        type            fixedValue;
+        value           uniform 500.0;
+    }
+}
+
+// ************************************************************************* //

tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOMSpecularCeiling/0/U

@@ -0,0 +1,48 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | foam-extend: Open Source CFD                    |
+|  \\    /   O peration     | Version:     3.1                                |
+|   \\  /    A nd           | Web:         http://www.extend-project.de       |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    floor
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+
+    ceiling
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+
+    fixedWalls
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+
+    box
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+}
+
+// ************************************************************************* //

tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOMSpecularCeiling/0/alphat

@@ -0,0 +1,50 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | foam-extend: Open Source CFD                    |
+|  \\    /   O peration     | Version:     3.1                                |
+|   \\  /    A nd           | Web:         http://www.extend-project.de       |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       volScalarField;
+    location    "0";
+    object      alphat;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    box
+    {
+        type            alphatWallFunction;
+        Prt             0.85;
+        value           uniform 0;
+    }
+    floor
+    {
+        type            alphatWallFunction;
+        Prt             0.85;
+        value           uniform 0;
+    }
+    ceiling
+    {
+        type            alphatWallFunction;
+        Prt             0.85;
+        value           uniform 0;
+    }
+    fixedWalls
+    {
+        type            alphatWallFunction;
+        Prt             0.85;
+        value           uniform 0;
+    }
+}
+
+// ************************************************************************* //

tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOMSpecularCeiling/0/epsilon

@@ -0,0 +1,58 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | foam-extend: Open Source CFD                    |
+|  \\    /   O peration     | Version:     3.1                                |
+|   \\  /    A nd           | Web:         http://www.extend-project.de       |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       volScalarField;
+    location    "0";
+    object      epsilon;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -3 0 0 0 0];
+
+internalField   uniform 0.01;
+
+boundaryField
+{
+    box
+    {
+        type            compressible::epsilonWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 0.01;
+    }
+    floor
+    {
+        type            compressible::epsilonWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 0.01;
+    }
+    ceiling
+    {
+        type            compressible::epsilonWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 0.01;
+    }
+    fixedWalls
+    {
+        type            compressible::epsilonWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 0.01;
+    }
+}
+
+// ************************************************************************* //

tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOMSpecularCeiling/0/k

@@ -0,0 +1,46 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | foam-extend: Open Source CFD                    |
+|  \\    /   O peration     | Version:     3.1                                |
+|   \\  /    A nd           | Web:         http://www.extend-project.de       |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       volScalarField;
+    location    "0";
+    object      k;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 0.1;
+
+boundaryField
+{
+    box
+    {
+        type            compressible::kqRWallFunction;
+        value           uniform 0.1;
+    }
+    floor
+    {
+        type            compressible::kqRWallFunction;
+        value           uniform 0.1;
+    }
+    ceiling
+    {
+        type            compressible::kqRWallFunction;
+        value           uniform 0.1;
+    }
+    fixedWalls
+    {
+        type            compressible::kqRWallFunction;
+        value           uniform 0.1;
+    }
+}
+
+// ************************************************************************* //

tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOMSpecularCeiling/0/mut

@@ -0,0 +1,58 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | foam-extend: Open Source CFD                    |
+|  \\    /   O peration     | Version:     3.1                                |
+|   \\  /    A nd           | Web:         http://www.extend-project.de       |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       volScalarField;
+    location    "0";
+    object      mut;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    box
+    {
+        type            mutWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 0;
+    }
+    floor
+    {
+        type            mutWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 0;
+    }
+    ceiling
+    {
+        type            mutWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 0;
+    }
+    fixedWalls
+    {
+        type            mutWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 0;
+    }
+}
+
+// ************************************************************************* //

tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOMSpecularCeiling/0/p

@@ -0,0 +1,48 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | foam-extend: Open Source CFD                    |
+|  \\    /   O peration     | Version:     3.1                                |
+|   \\  /    A nd           | Web:         http://www.extend-project.de       |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 100000;
+
+boundaryField
+{
+    floor
+    {
+        type            buoyantPressure;
+        value           uniform 100000;
+    }
+
+    ceiling
+    {
+        type            buoyantPressure;
+        value           uniform 100000;
+    }
+
+    fixedWalls
+    {
+        type            buoyantPressure;
+        value           uniform 100000;
+    }
+
+    box
+    {
+        type            buoyantPressure;
+        value           uniform 100000;
+    }
+}
+
+// ************************************************************************* //

tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOMSpecularCeiling/constant/RASProperties

@@ -0,0 +1,24 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | foam-extend: Open Source CFD                    |
+|  \\    /   O peration     | Version:     3.1                                |
+|   \\  /    A nd           | Web:         http://www.extend-project.de       |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      RASProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+RASModel        kEpsilon;
+
+turbulence      on;
+
+printCoeffs     on;
+
+// ************************************************************************* //

tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOMSpecularCeiling/constant/g

@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | foam-extend: Open Source CFD                    |
+|  \\    /   O peration     | Version:     3.1                                |
+|   \\  /    A nd           | Web:         http://www.extend-project.de       |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       uniformDimensionedVectorField;
+    location    "constant";
+    object      g;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -2 0 0 0 0];
+value           ( 0 0 -9.81 );
+
+// ************************************************************************* //

tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOMSpecularCeiling/constant/polyMesh/blockMeshDict

@@ -0,0 +1,186 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | foam-extend: Open Source CFD                    |
+|  \\    /   O peration     | Version:     3.1                                |
+|   \\  /    A nd           | Web:         http://www.extend-project.de       |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 1;
+
+vertices
+(
+    ( 0.0  0.0  0.0)
+    ( 0.5  0.0  0.0)
+    ( 1.5  0.0  0.0)
+    (10.0  0.0  0.0)
+    ( 0.0  0.5  0.0)
+    ( 0.5  0.5  0.0)
+    ( 1.5  0.5  0.0)
+    (10.0  0.5  0.0)
+    ( 0.0  1.5  0.0)
+    ( 0.5  1.5  0.0)
+    ( 1.5  1.5  0.0)
+    (10.0  1.5  0.0)
+    ( 0.0  6.0  0.0)
+    ( 0.5  6.0  0.0)
+    ( 1.5  6.0  0.0)
+    (10.0  6.0  0.0)
+
+    ( 0.0  0.0  0.5)
+    ( 0.5  0.0  0.5)
+    ( 1.5  0.0  0.5)
+    (10.0  0.0  0.5)
+    ( 0.0  0.5  0.5)
+    ( 0.5  0.5  0.5)
+    ( 1.5  0.5  0.5)
+    (10.0  0.5  0.5)
+    ( 0.0  1.5  0.5)
+    ( 0.5  1.5  0.5)
+    ( 1.5  1.5  0.5)
+    (10.0  1.5  0.5)
+    ( 0.0  6.0  0.5)
+    ( 0.5  6.0  0.5)
+    ( 1.5  6.0  0.5)
+    (10.0  6.0  0.5)
+
+    ( 0.0  0.0  2.0)
+    ( 0.5  0.0  2.0)
+    ( 1.5  0.0  2.0)
+    (10.0  0.0  2.0)
+    ( 0.0  0.5  2.0)
+    ( 0.5  0.5  2.0)
+    ( 1.5  0.5  2.0)
+    (10.0  0.5  2.0)
+    ( 0.0  1.5  2.0)
+    ( 0.5  1.5  2.0)
+    ( 1.5  1.5  2.0)
+    (10.0  1.5  2.0)
+    ( 0.0  6.0  2.0)
+    ( 0.5  6.0  2.0)
+    ( 1.5  6.0  2.0)
+    (10.0  6.0  2.0)
+);
+
+blocks
+(
+    hex ( 0  1  5  4 16 17 21 20) ( 5  5  5) simpleGrading (1 1 1)
+    hex ( 1  2  6  5 17 18 22 21) (10  5  5) simpleGrading (1 1 1)
+    hex ( 2  3  7  6 18 19 23 22) (80  5  5) simpleGrading (1 1 1)
+    hex ( 4  5  9  8 20 21 25 24) ( 5 10  5) simpleGrading (1 1 1)
+    hex ( 6  7 11 10 22 23 27 26) (80 10  5) simpleGrading (1 1 1)
+    hex ( 8  9 13 12 24 25 29 28) ( 5 40  5) simpleGrading (1 1 1)
+    hex ( 9 10 14 13 25 26 30 29) (10 40  5) simpleGrading (1 1 1)
+    hex (10 11 15 14 26 27 31 30) (80 40  5) simpleGrading (1 1 1)
+
+    hex (16 17 21 20 32 33 37 36) ( 5  5 15) simpleGrading (1 1 1)
+    hex (17 18 22 21 33 34 38 37) (10  5 15) simpleGrading (1 1 1)
+    hex (18 19 23 22 34 35 39 38) (80  5 15) simpleGrading (1 1 1)
+    hex (20 21 25 24 36 37 41 40) ( 5 10 15) simpleGrading (1 1 1)
+
+    hex (21 22 26 25 37 38 42 41) (10 10 15) simpleGrading (1 1 1)
+
+    hex (22 23 27 26 38 39 43 42) (80 10 15) simpleGrading (1 1 1)
+    hex (24 25 29 28 40 41 45 44) ( 5 40 15) simpleGrading (1 1 1)
+    hex (25 26 30 29 41 42 46 45) (10 40 15) simpleGrading (1 1 1)
+    hex (26 27 31 30 42 43 47 46) (80 40 15) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+    box
+    {
+        type wall;
+        faces
+        (
+            ( 6 22 21  5)
+            (10 26 22  6)
+            ( 9 25 26 10)
+            ( 5 21 25  9)
+            (22 26 25 21)
+        );
+    }
+    floor
+    {
+        type wall;
+        faces
+        (
+            ( 1  5  4  0)
+            ( 2  6  5  1)
+            ( 3  7  6  2)
+            ( 5  9  8  4)
+            ( 7 11 10  6)
+            ( 9 13 12  8)
+            (10 14 13  9)
+            (11 15 14 10)
+        );
+    }
+    ceiling
+    {
+        type wall;
+        faces
+        (
+            (33 37 36 32)
+            (34 38 37 33)
+            (35 39 38 34)
+            (37 41 40 36)
+            (38 42 41 37)
+            (39 43 42 38)
+            (41 45 44 40)
+            (42 46 45 41)
+            (43 47 46 42)
+        );
+    }
+    fixedWalls
+    {
+        type wall;
+        faces
+        (
+            ( 1 17 16  0)
+            ( 2 18 17  1)
+            ( 3 19 18  2)
+            (17 33 32 16)
+            (18 34 33 17)
+            (19 35 34 18)
+
+            ( 7 23 19  3)
+            (11 27 23  7)
+            (15 31 27 11)
+            (23 39 35 19)
+            (27 43 39 23)
+            (31 47 43 27)
+
+            (14 30 31 15)
+            (13 29 30 14)
+            (12 28 29 13)
+            (30 46 47 31)
+            (29 45 46 30)
+            (28 44 45 29)
+
+            ( 8 24 28 12)
+            ( 4 20 24  8)
+            ( 0 16 20  4)
+            (24 40 44 28)
+            (20 36 40 24)
+            (16 32 36 20)
+        );
+    }
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //

tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOMSpecularCeiling/constant/radiationProperties

@@ -0,0 +1,58 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | foam-extend: Open Source CFD                    |
+|  \\    /   O peration     | Version:     3.1                                |
+|   \\  /    A nd           | Web:         http://www.extend-project.de       |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      radiationProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+radiation       on;
+
+radiationModel  fvDOM;
+
+noRadiation
+{
+}
+
+P1Coeffs
+{
+}
+
+fvDOMCoeffs
+{
+    nPhi    4;          // azimuthal angles in PI/2 on X-Y.(from Y to X)
+    nTheta  4;          // polar angles in PI (from Z to X-Y plane)
+    convergence 1e-3;   // convergence criteria for radiation iteration
+    maxIter 10;         // maximum number of iterations
+}
+
+// Number of flow iterations per radiation iteration
+solverFreq 10;
+
+absorptionEmissionModel constantAbsorptionEmission;
+
+constantAbsorptionEmissionCoeffs
+{
+    a               a [ 0 -1 0 0 0 0 0 ] 0.01;
+    e               e [ 0 -1 0 0 0 0 0 ] 0;
+    E               E [ 1 -1 -3 0 0 0 0 ] 0;
+}
+
+scatterModel    constantScatter;
+
+constantScatterCoeffs
+{
+    sigma           sigma [ 0 -1 0 0 0 0 0 ] 0;
+    C               C [ 0 0 0 0 0 0 0 ] 0;
+}
+
+// ************************************************************************* //

tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOMSpecularCeiling/constant/thermophysicalProperties

@@ -0,0 +1,24 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | foam-extend: Open Source CFD                    |
+|  \\    /   O peration     | Version:     3.1                                |
+|   \\  /    A nd           | Web:         http://www.extend-project.de       |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType      hPsiThermo<pureMixture<constTransport<specieThermo<hConstThermo<perfectGas>>>>>;
+
+mixture         air 1 28.9 1000 0 1.8e-05 0.7;
+
+pRef            100000;
+
+// ************************************************************************* //

tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOMSpecularCeiling/system/controlDict

@@ -0,0 +1,48 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | foam-extend: Open Source CFD                    |
+|  \\    /   O peration     | Version:     3.1                                |
+|   \\  /    A nd           | Web:         http://www.extend-project.de       |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     buoyantSimpleRadiationFoam;
+
+startFrom       latestTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         1000;
+
+deltaT          1;
+
+writeControl    timeStep;
+
+writeInterval   100;
+
+purgeWrite      0;
+
+writeFormat     binary;
+
+writePrecision  6;
+
+writeCompression uncompressed;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable yes;
+
+// ************************************************************************* //

tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOMSpecularCeiling/system/fvSchemes

@@ -0,0 +1,69 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | foam-extend: Open Source CFD                    |
+|  \\    /   O peration     | Version:     3.1                                |
+|   \\  /    A nd           | Web:         http://www.extend-project.de       |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         steadyState;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+}
+
+divSchemes
+{
+    default         none;
+    div(phi,U)      Gauss upwind;
+    div(phi,h)      Gauss upwind;
+    div(phi,k)      Gauss upwind;
+    div(phi,epsilon) Gauss upwind;
+    div(phi,R)      Gauss upwind;
+    div(R)          Gauss linear;
+    div(Ji,Ii_h)    Gauss linearUpwind Gauss linear; //Gauss upwind;
+    div((muEff*dev2(grad(U).T()))) Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         none;
+    laplacian(muEff,U) Gauss linear corrected;
+    laplacian((rho*(1|A(U))),p) Gauss linear corrected;
+    laplacian(alphaEff,h) Gauss linear corrected;
+    laplacian(DkEff,k) Gauss linear corrected;
+    laplacian(DepsilonEff,epsilon) Gauss linear corrected;
+    laplacian(DREff,R) Gauss linear corrected;
+    laplacian(gammaRad,G) Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+fluxRequired
+{
+    default         no;
+    p;
+}
+
+// ************************************************************************* //

tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOMSpecularCeiling/system/fvSolution

@@ -0,0 +1,91 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | foam-extend: Open Source CFD                    |
+|  \\    /   O peration     | Version:     3.1                                |
+|   \\  /    A nd           | Web:         http://www.extend-project.de       |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    p
+    {
+        solver          GAMG;
+        tolerance       1e-06;
+        relTol          0.01;
+        smoother        GaussSeidel;
+        cacheAgglomeration true;
+        nCellsInCoarsestLevel 10;
+        agglomerator    faceAreaPair;
+        mergeLevels     1;
+    }
+
+    U
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-05;
+        relTol          0.1;
+    }
+
+    h
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-05;
+        relTol          0.1;
+    }
+
+    k
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-05;
+        relTol          0.001;
+    }
+
+    epsilon
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-05;
+        relTol          0.1;
+    }
+
+    Ii
+    {
+        solver          BiCGStab;
+        preconditioner  DILU;
+        tolerance       1e-05;
+        relTol          0;
+    }
+}
+
+SIMPLE
+{
+    nNonOrthogonalCorrectors 0;
+    pRefCell       0;
+    pRefValue      0;
+}
+
+relaxationFactors
+{
+    rho             1.0;
+    p               0.3;
+    U               0.7;
+    h               0.7;
+    k               0.7;
+    epsilon         0.7;
+    "ILambda.*"     1;
+}
+
+// ************************************************************************* //