solarHeatLoad.C

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.
 
    OpenFOAM 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.
 
    OpenFOAM 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 OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
\*---------------------------------------------------------------------------*/
 
#include "solarHeatLoad.H"
#include "surfaceFields.H"
#include "vectorList.H"
#include "addToRunTimeSelectionTable.H"
#include "boundaryRadiationProperties.H"
#include "uniformDimensionedFields.H"
#include "cyclicAMIPolyPatch.H"
#include "mappedPatchBase.H"
#include "wallPolyPatch.H"
 
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
namespace Foam
{
    namespace radiation
    {
        defineTypeNameAndDebug(solarHeatLoad, 1);
        addToRadiationRunTimeSelectionTables(solarHeatLoad);
    }
}
 
const Foam::word Foam::radiation::solarHeatLoad::viewFactorWalls
    = "viewFactorWall";
 
// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
bool Foam::radiation::solarHeatLoad::updateHitFaces()
{
    if (hitFaces_.empty())
    {
         hitFaces_.reset(new faceHeatShading(mesh_, solarCalc_.direction()));
         return true;
    }
    else
    {
        switch (solarCalc_.sunDirectionModel())
        {
            case solarCalculator::mSunDirConstant:
            {
                return false;
                break;
            }
            case solarCalculator::mSunDirTraking:
            {
                label updateIndex = label
                (
                    mesh_.time().value()/solarCalc_.sunTrackingUpdateInterval()
                );
 
                if (updateIndex > updateTimeIndex_)
                {
                    Info << "Updating Sun position..." << endl;
                    updateTimeIndex_ = updateIndex;
                    solarCalc_.correctSunDirection();
                    hitFaces_->direction() = solarCalc_.direction();
                    hitFaces_->correct();
//Info << "now the sun direction is : " << solarCalc_.direction();
//add write for debug
//                  if (debug)
//                  {
//                      hitFaces_->write();
//                  }
                    return true;
                    break;
                }
            }
        }
    }
 
    return false;
}
 
 
void Foam::radiation::solarHeatLoad::updateAbsorptivity
(
    const labelHashSet& includePatches
)
{
    const boundaryRadiationProperties& boundaryRadiation =
        boundaryRadiationProperties::New(mesh_);
 
    forAllConstIter(labelHashSet, includePatches, iter)
    {
        const label patchID = iter.key();
        absorptivity_[patchID].setSize(nBands_);
        for (label bandI = 0; bandI < nBands_; bandI++)
        {
            absorptivity_[patchID][bandI] =
                boundaryRadiation.absorptivity(patchID, bandI);
        }
    }
}
 
 
void Foam::radiation::solarHeatLoad::updateDirectHitRadiation
(
    const labelList& hitFacesId,
    const labelHashSet& includeMappedPatchBasePatches
)
{
    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
    const scalarField& V = mesh_.V();
 
    forAll (hitFacesId, i)
    {
        const label faceI = hitFacesId[i];
        label patchID = patches.whichPatch(faceI);
        const polyPatch& pp = patches[patchID];
        const label localFaceI = faceI - pp.start();
        const vector qPrim =
              solarCalc_.directSolarRad()
            * solarCalc_.direction();
//Info << "cell hit face " << i << " from solar direction is  " << solarCalc_.direction()<< endl;
//Info << "cell hit face " << i << " 's solar load  is  " << solarCalc_.directSolarRad()<< endl;
//Info << "cell hit face " << i <<" heat from sun is " << qPrim<< endl<< endl;
 
//        if (includeMappedPatchBasePatches[patchID])
        {
            const vectorField n = pp.faceNormals();
 
            for (label bandI = 0; bandI < nBands_; bandI++)
            {
 
//Info <<endl<< "for patch name " << pp.name() << endl;
//Info << "cell hit face [" << i <<"] spectralDistribution is " << spectralDistribution_[bandI]<< endl;
//Info << "cell hit face [" << i <<"] absorptivity is " <<absorptivity_[patchID][bandI]()[localFaceI]<< endl;
//Info << "For Qr, face normals is " << n[localFaceI] << endl;
//Info << "For Qr, qPrm & face normals ==  " << (qPrim & n[localFaceI]) << endl<<endl;
                Qr_.boundaryField()[patchID][localFaceI] +=
                    (qPrim & n[localFaceI])
                  * spectralDistribution_[bandI]
                  * absorptivity_[patchID][bandI]()[localFaceI];
            }
        }
 //       else
        {
            const vectorField& sf = mesh_.Sf().boundaryField()[patchID];
            const label cellI = pp.faceCells()[localFaceI];
 
            for (label bandI = 0; bandI < nBands_; bandI++)
            {
//Info<<endl << "face normals is " << sf[localFaceI] << endl;
//Info << "qPrm & face normals == is " << (qPrim & sf[localFaceI]) << endl<<endl;
                Ru_[cellI] +=
                    (qPrim & sf[localFaceI])
                  * spectralDistribution_[bandI]
                  * absorptivity_[patchID][bandI]()[localFaceI]
                  / V[cellI];
 
                //Ru_[cellI]=0;
            }
        }
    }
//  Info << "All direct radiation heat flux is " << Qr_ << endl;
}
 
void Foam::radiation::solarHeatLoad::updateSkyDiffusiveRadiation
(
    const labelHashSet& includePatches,
    const labelHashSet& includeMappedPatchBasePatches
)
{
    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
    const scalarField& V = mesh_.V();
 
    switch(solarCalc_.sunLoadModel())
    {
        case solarCalculator::mSunLoadFairWeatherConditions:
        case solarCalculator::mSunLoadTheoreticalMaximum:
        {
            forAllConstIter(labelHashSet, includePatches, iter)
            {
                const label patchID = iter.key();
                const polyPatch& pp = patches[patchID];
                const scalarField& sf = mesh_.magSf().boundaryField()[patchID];
 
                const vectorField n = pp.faceNormals();
                const labelList& cellIds = pp.faceCells();
 
                forAll (n, faceI)
                {
                    const scalar cosEpsilon(verticalDir_ & -n[faceI]);
 
                    scalar Ed(0.0);
                    scalar Er(0.0);
                    const scalar cosTheta(solarCalc_.direction() & -n[faceI]);
 
                    {
                        // Above the horizon
                        if (cosEpsilon == 0.0)
                        {
                            // Vertical walls
                            scalar Y(0);
 
                            if (cosTheta > -0.2)
                            {
                                Y = 0.55+0.437*cosTheta + 0.313*sqr(cosTheta);
                            }
                            else
                            {
                                Y = 0.45;
                            }
 
                            Ed = solarCalc_.C()*Y*solarCalc_.directSolarRad();
                        }
                        else
                        {
                            //Other than vertical walls
                            Ed =
                                solarCalc_.C()
                              * solarCalc_.directSolarRad()
                              * (1.0 + cosEpsilon)/2.0;
                        }
 
                        // Ground reflected
                        Er =
                            solarCalc_.directSolarRad()
                         * (solarCalc_.C() + Foam::sin(solarCalc_.beta()))
                         * solarCalc_.groundReflectivity()
                         * (1.0 - cosEpsilon)/2.0;
                    }
 
                    const label cellI = cellIds[faceI];
                 //   if (includeMappedPatchBasePatches[patchID])
//本来这部分includeMappedPatchBasePatches判断内部的计算的Q辐射量是为了耦合边界使用的
                    {
                        for (label bandI = 0; bandI < nBands_; bandI++)
                        {
                // 这个是计算所有的天空漫反射，Qr_是包含天空漫反射和地面反射的辐射量
                            QdiffRad_.boundaryField()[patchID][faceI] +=
                                (Ed + Er)
                              * spectralDistribution_[bandI]
                              * absorptivity_[patchID][bandI]()[faceI];
                        }
                    }
 
 
                  //  else
                    //这里是方程计算会用到的Ru
                    {
                        for (label bandI = 0; bandI < nBands_; bandI++)
                        {
                            Ru_[cellI] +=
                                  (Ed + Er)
                                * spectralDistribution_[bandI]
                                * absorptivity_[patchID][bandI]()[faceI]
                                * sf[faceI]/V[cellI];
 
                            //Ru_[cellI]=0;
                        }
                    }
                }
            }
//          Info << "theor max  Sky reflective  radiation heat flux is " << QdiffRad_ << endl;
        }
        break;
 
        case solarCalculator::mSunLoadConstant:
        {
            forAllConstIter(labelHashSet, includePatches, iter)
            {
                const label patchID = iter.key();
                const polyPatch& pp = patches[patchID];
                const scalarField& sf = mesh_.magSf().boundaryField()[patchID];
 
                const labelList& cellIds = pp.faceCells();
                forAll (pp, faceI)
                {
                    const label cellI = cellIds[faceI];
                    //if (includeMappedPatchBasePatches[patchID])
                    {
                        for (label bandI = 0; bandI < nBands_; bandI++)
                        {
 
                            QdiffRad_.boundaryField()[patchID][faceI] +=
                                solarCalc_.diffuseSolarRad()
                              * spectralDistribution_[bandI]
                              * absorptivity_[patchID][bandI]()[faceI];
                        }
                    }
                    //else
                    {
                        for (label bandI = 0; bandI < nBands_; bandI++)
                        {
                            Ru_[cellI] +=
                                (
                                    spectralDistribution_[bandI]
                                  * absorptivity_[patchID][bandI]()[faceI]
                                  * solarCalc_.diffuseSolarRad()
                                )*sf[faceI]/V[cellI];
 
                        //  Ru_[cellI]=0;
                        }
                    }
                }
            }
//          Info << "constant model, Sky reflective  radiation heat flux is " << QdiffRad_ << endl;
            break;
        }
    }
}
 
 
void Foam::radiation::solarHeatLoad::initialise(const dictionary& coeffs)
{
 
    if (coeffs.found("gridUp"))
    {
         coeffs.lookup("gridUp") >> verticalDir_;
         verticalDir_ /= mag(verticalDir_);
    }
    else if (mesh_.foundObject<uniformDimensionedVectorField>("g"))
    {
        const uniformDimensionedVectorField& g =
            mesh_.lookupObject<uniformDimensionedVectorField>("g");
        verticalDir_ = (-g/mag(g)).value();
    }
 
    includePatches_ = mesh_.boundaryMesh().findIndices(viewFactorWalls);
 
    coeffs.lookup("useVFbeamToDiffuse") >> useVFbeamToDiffuse_;
 
    coeffs.lookup("spectralDistribution") >> spectralDistribution_;
    spectralDistribution_ =
        spectralDistribution_/sum(spectralDistribution_);
 
    nBands_ = spectralDistribution_.size();
 
    if (useVFbeamToDiffuse_)
    {
        map_.reset
        (
            new IOmapDistribute
            (
                IOobject
                (
                    "mapDist",
                    mesh_.facesInstance(),
                    mesh_,
                    IOobject::MUST_READ,
                    IOobject::NO_WRITE
                )
            )
        );
    }
 
    if (coeffs.found("solidCoupled"))
    {
         coeffs.lookup("solidCoupled") >> solidCoupled_;
    }
 
    if (coeffs.found("wallCoupled"))
    {
         coeffs.lookup("wallCoupled") >> wallCoupled_;
    }
 
    if (coeffs.found("updateAbsorptivity"))
    {
        coeffs.lookup("updateAbsorptivity") >> updateAbsorptivity_;
    }
}
 
 
void Foam::radiation::solarHeatLoad::calculateQdiff
(
    const labelHashSet& includePatches,
    const labelHashSet& includeMappedPatchBasePatches
)
{
 
    scalarListIOList FmyProc
    (
        IOobject
        (
            "F",
            mesh_.facesInstance(),
            mesh_,
            IOobject::MUST_READ,
            IOobject::NO_WRITE,
            false
        )
    );
 
 
    if (finalAgglom_.size() > 0 && coarseMesh_.empty())
    {
        coarseMesh_.reset
        (
            new singleCellFvMesh
            (
                IOobject
                (
                    "coarse." + mesh_.name(),
                    mesh_.polyMesh::instance(),
                    mesh_.time(),
                    IOobject::NO_READ,
                    IOobject::NO_WRITE
                ),
                mesh_,
                finalAgglom_
            )
        );
    }
 
    label nLocalVFCoarseFaces = 0;
    forAll(includePatches_, i)
    {
        const label patchI = includePatches_[i];
        nLocalVFCoarseFaces +=
            coarseMesh_->boundaryMesh()[patchI].size();
    }
 
    label totalFVNCoarseFaces = nLocalVFCoarseFaces;
    reduce(totalFVNCoarseFaces, sumOp<label>());
 
    // Calculate weighted absorptivity on coarse patches
    List<scalar> localCoarseRave(nLocalVFCoarseFaces);
    List<scalar> localCoarsePartialArea(nLocalVFCoarseFaces);
    List<vector> localCoarseNorm(nLocalVFCoarseFaces);
 
    scalarField compactCoarseRave(map_->constructSize(), 0.0);
    scalarField compactCoarsePartialArea(map_->constructSize(), 0.0);
    vectorList compactCoarseNorm(map_->constructSize(), vector::zero);
 
    const boundaryRadiationProperties& boundaryRadiation =
        boundaryRadiationProperties::New(mesh_);
 
    coarseToFine_.setSize(includePatches_.size());
 
    const labelList& hitFacesId = hitFaces_->rayStartFaces();
 
    label startI = 0;
    label compactI = 0;
    forAll (includePatches_, i)
    {
        const label patchID = includePatches_[i];
        const polyPatch& pp = mesh_.boundaryMesh()[patchID];
 
        const polyPatch& cpp = coarseMesh_->boundaryMesh()[patchID];
 
        const labelList& agglom = finalAgglom_[patchID];
        //if (pp.size() > 0)
        if (agglom.size() > 0)
        {
            label nAgglom = max(agglom) + 1;
            coarseToFine_[i] = invertOneToMany(nAgglom, agglom);
        }
 
        scalarField r(pp.size(), 0.0);
        for (label bandI = 0; bandI < nBands_; bandI++)
        {
            const tmp<scalarField> tr =
               spectralDistribution_[bandI]
              *boundaryRadiation.reflectivity(patchID, bandI);
            r += tr();
        }
 
        scalarList Rave(cpp.size(), 0.0);
        scalarList area(cpp.size(), 0.0);
 
        const scalarField& sf = mesh_.magSf().boundaryField()[patchID];
 
        const labelList& coarsePatchFace =
            coarseMesh_->patchFaceMap()[patchID];
 
        forAll (cpp, coarseI)
        {
            const label coarseFaceID = coarsePatchFace[coarseI];
            const labelList& fineFaces =
                coarseToFine_[i][coarseFaceID];
 
            UIndirectList<scalar> fineSf
            (
                sf,
                fineFaces
            );
            scalar fineArea = sum(fineSf());
 
            scalar fullArea = 0.0;
            forAll(fineFaces, j)
            {
                label faceI = fineFaces[j];
                label globalFaceI = faceI + pp.start();
 
                if (findIndex(hitFacesId, globalFaceI) != -1)
                {
                    fullArea += sf[faceI];
                }
                Rave[coarseI] += (r[faceI]*sf[faceI])/fineArea;
            }
            localCoarsePartialArea[compactI++] = fullArea/fineArea;
        }
 
        SubList<scalar>
        (
            localCoarseRave,
            Rave.size(),
            startI
        ).assign(Rave);
 
 
        const vectorList coarseNSf = cpp.faceNormals();
        SubList<vector>
        (
            localCoarseNorm,
            cpp.size(),
            startI
        ).assign(coarseNSf);
        startI += cpp.size();
    }
 
 
     SubList<scalar>(compactCoarsePartialArea, nLocalVFCoarseFaces).assign
     (
         localCoarsePartialArea
     );
 
     SubList<scalar>(compactCoarseRave, nLocalVFCoarseFaces).assign
     (
         localCoarseRave
     );
 
     SubList<vector>(compactCoarseNorm, nLocalVFCoarseFaces).assign
     (
         localCoarseNorm
     );
 
     map_->distribute(compactCoarsePartialArea);
     map_->distribute(compactCoarseRave);
     map_->distribute(compactCoarseNorm);
 
 
    // Calculate coarse hitFaces and Sun direct hit heat fluxes
    scalarList localqDiffusive(nLocalVFCoarseFaces, 0.0);
 
    label locaFaceI = 0;
    forAll (includePatches_, i)
    {
        const label patchID = includePatches_[i];
        const polyPatch& pp = coarseMesh_->boundaryMesh()[patchID];
        const polyPatch& ppf = mesh_.boundaryMesh()[patchID];
 
        const labelList& coarsePatchFace =
            coarseMesh_->patchFaceMap()[patchID];
        const scalarField& sf = mesh_.magSf().boundaryField()[patchID];
 
        scalarField a(ppf.size(), 0.0);
        for (label bandI = 0; bandI < nBands_; bandI++)
        {
            const tmp<scalarField> ta =
                spectralDistribution_[bandI]
              * absorptivity_[patchID][bandI]();
            a += ta();
        }
 
        forAll (pp, coarseI)
        {
            const label coarseFaceID = coarsePatchFace[coarseI];
            const labelList& fineFaces = coarseToFine_[i][coarseFaceID];
            UIndirectList<scalar> fineSf
            (
                sf,
                fineFaces
            );
            scalar fineArea = sum(fineSf());
 
            scalar aAve = 0.0;
            forAll(fineFaces, j)
            {
                label faceI = fineFaces[j];
                aAve += (a[faceI]*sf[faceI])/fineArea;
            }
 
            const scalarList& vf = FmyProc[locaFaceI];
 
            const labelList& compactFaces = visibleFaceFaces_[locaFaceI];
 
 
            forAll (compactFaces, j)
            {
                label compactI = compactFaces[j];
 
                localqDiffusive[locaFaceI] +=
                      compactCoarsePartialArea[compactI]
                    * aAve
                    * (solarCalc_.directSolarRad()*solarCalc_.direction())
                    & compactCoarseNorm[compactI]
                    * vf[j]
                    * compactCoarseRave[compactI];
 
            }
            locaFaceI++;
        }
    }
 
    // Fill QsecondRad_
    label compactId = 0;
    forAll (includePatches_, i)
    {
        const label patchID = includePatches_[i];
        const polyPatch& pp = coarseMesh_->boundaryMesh()[patchID];
 
        if (pp.size() > 0)
        {
            scalarField& Qrp = QrefecFormOtherWallRad_.boundaryField()[patchID];
 
            const labelList& coarsePatchFace =
                coarseMesh_->patchFaceMap()[patchID];
 
            forAll(pp, coarseI)
            {
                const label coarseFaceID = coarsePatchFace[coarseI];
                const labelList& fineFaces = coarseToFine_[i][coarseFaceID];
                forAll(fineFaces, k)
                {
                    label faceI = fineFaces[k];
                    Qrp[faceI] = localqDiffusive[compactId];
                }
                compactId ++;
            }
        }
    }
 
    const scalarField& V = mesh_.V();
    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
    forAllConstIter(labelHashSet, includePatches, iter)
    {
        const label patchID = iter.key();
        const scalarField& qSecond = QrefecFormOtherWallRad_.boundaryField()[patchID];
//Info << "other wall  reflective  radiation heat flux is " << QrefecFormOtherWallRad_ << endl;
        //if (includeMappedPatchBasePatches[patchID])
        {
            QdiffRad_.boundaryField()[patchID] += qSecond;//add the groundReflectivity radiation to the toatal reflective radiation 
        }
        //else
        {
            const polyPatch& pp = patches[patchID];
            const labelList& cellIds = pp.faceCells();
            const scalarField& sf = mesh_.magSf().boundaryField()[patchID];
            forAll (pp, faceI)
            {
                const label cellI = cellIds[faceI];
                Ru_[cellI] += qSecond[faceI]*sf[faceI]/V[cellI];
 
                //Ru_[cellI]=0;
            }
        }
    }
//  Info << "add other wall reflective heat, total diffuse Solar radiation heat flux is " << QdiffRad_ << endl;
}
 
 
// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
Foam::radiation::solarHeatLoad::solarHeatLoad(const volScalarField& T)
:
    radiationModel(typeName, T),
    finalAgglom_
    (
        IOobject
        (
            "finalAgglom",
            mesh_.facesInstance(),
            mesh_,
            IOobject::READ_IF_PRESENT,
            IOobject::NO_WRITE,
            false
        )
    ),
    coarseMesh_(),
    Qr_
    (
        IOobject
        (
            "Qr",
            mesh_.time().timeName(),
            mesh_,
            IOobject::READ_IF_PRESENT,
            IOobject::AUTO_WRITE
        ),
        mesh_,
        dimensionedScalar("Qr", dimMass/pow3(dimTime), 0.0)
    ),
   QrefecFormOtherWallRad_ 
    (
        IOobject
        (
            "Qrefec",
            mesh_.time().timeName(),
            mesh_,
            IOobject::READ_IF_PRESENT,
            IOobject::NO_WRITE
        ),
        mesh_,
        dimensionedScalar("Qrefec", dimMass/pow3(dimTime), 0.0)
    ),
    QdiffRad_
    (
        IOobject
        (
            "QdiffRad",
            mesh_.time().timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        mesh_,
        dimensionedScalar("QdiffRad", dimMass/pow3(dimTime), 0.0)
    ),
    QtotalRad_
    (
        IOobject
        (
            "QtotalRad_",
            mesh_.time().timeName(),
            mesh_,
            IOobject::READ_IF_PRESENT,
            IOobject::AUTO_WRITE
        ),
        mesh_,
        dimensionedScalar("QtotalRad_", dimMass/pow3(dimTime), 0.0)
    ),
    hitFaces_(),
    Ru_
    (
        IOobject
        (
            "Ru",
            mesh_.time().timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        mesh_,
        dimensionedScalar("Ru", dimMass/dimLength/pow3(dimTime), 0.0)
    ),
    solarCalc_(this->subDict(typeName + "Coeffs"), mesh_),
    verticalDir_(vector::zero),
    useVFbeamToDiffuse_(false),
    includePatches_(mesh_.boundary().size(), -1),
    coarseToFine_(),
    nBands_(1),
    spectralDistribution_(nBands_),
    map_(),
    visibleFaceFaces_
    (
        IOobject
        (
            "visibleFaceFaces",
            mesh_.facesInstance(),
            mesh_,
            IOobject::READ_IF_PRESENT,
            IOobject::NO_WRITE,
            false
        )
    ),
    solidCoupled_(true),
    absorptivity_(mesh_.boundaryMesh().size()),
    updateAbsorptivity_(false),
    firstIter_(true),
    updateTimeIndex_(0)
{
    initialise(coeffs_);
}
 
 
Foam::radiation::solarHeatLoad::solarHeatLoad
(
    const dictionary& dict,
    const volScalarField& T
)
:
    radiationModel(typeName, dict, T),
    finalAgglom_
    (
        IOobject
        (
            "finalAgglom",
            mesh_.facesInstance(),
            mesh_,
            IOobject::READ_IF_PRESENT,
            IOobject::NO_WRITE,
            false
        )
    ),
    coarseMesh_(),
    Qr_
    (
        IOobject
        (
            "Qr",
            mesh_.time().timeName(),
            mesh_,
            IOobject::READ_IF_PRESENT,
            IOobject::AUTO_WRITE
        ),
        mesh_,
        dimensionedScalar("Qr", dimMass/pow3(dimTime), 0.0)
    ),
    QrefecFormOtherWallRad_
    (
        IOobject
        (
            "Qrefec",
            mesh_.time().timeName(),
            mesh_,
            IOobject::READ_IF_PRESENT,
            IOobject::NO_WRITE
        ),
        mesh_,
        dimensionedScalar("Qrefec", dimMass/pow3(dimTime), 0.0)
    ),
    QdiffRad_
    (
        IOobject
        (
            "QdiffRad",
            mesh_.time().timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        mesh_,
        dimensionedScalar("QdiffRad", dimMass/pow3(dimTime), 0.0)
    ),
    QtotalRad_
    (
        IOobject
        (
            "QtotalRad_",
            mesh_.time().timeName(),
            mesh_,
            IOobject::READ_IF_PRESENT,
            IOobject::AUTO_WRITE
        ),
        mesh_,
        dimensionedScalar("QtotalRad_", dimMass/pow3(dimTime), 0.0)
    ),
    hitFaces_(),
    Ru_
    (
        IOobject
        (
            "Ru",
            mesh_.time().timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        mesh_,
        dimensionedScalar("Ru", dimMass/dimLength/pow3(dimTime), 0.0)
    ),
    solarCalc_(coeffs_, mesh_),
    verticalDir_(vector::zero),
    useVFbeamToDiffuse_(false),
    includePatches_(mesh_.boundary().size(), -1),
    coarseToFine_(),
    nBands_(1),
    spectralDistribution_(nBands_),
    map_(),
    visibleFaceFaces_
    (
        IOobject
        (
            "visibleFaceFaces",
            mesh_.facesInstance(),
            mesh_,
            IOobject::READ_IF_PRESENT,
            IOobject::NO_WRITE,
            false
        )
    ),
    solidCoupled_(true),
    wallCoupled_(false),
    absorptivity_(mesh_.boundaryMesh().size()),
    updateAbsorptivity_(false),
    firstIter_(true),
    updateTimeIndex_(0)
{
    initialise(coeffs_);
}
 
 
Foam::radiation::solarHeatLoad::solarHeatLoad
(
    const dictionary& dict,
    const volScalarField& T,
    const word radWallFieldName
)
:
    radiationModel("none", T),
    finalAgglom_
    (
        IOobject
        (
            "finalAgglom",
            mesh_.facesInstance(),
            mesh_,
            IOobject::READ_IF_PRESENT,
            IOobject::NO_WRITE,
            false
        )
    ),
    coarseMesh_(),
    Qr_
    (
        IOobject
        (
            radWallFieldName,
            mesh_.time().timeName(),
            mesh_,
            IOobject::READ_IF_PRESENT,
            IOobject::AUTO_WRITE
        ),
        mesh_,
        dimensionedScalar("Qr", dimMass/pow3(dimTime), 0.0)
    ),
    QrefecFormOtherWallRad_
    (
        IOobject
        (
            "Qrefec",
            mesh_.time().timeName(),
            mesh_,
            IOobject::READ_IF_PRESENT,
            IOobject::NO_WRITE
        ),
        mesh_,
        dimensionedScalar("Qrefec", dimMass/pow3(dimTime), 0.0)
    ),
    QdiffRad_
    (
        IOobject
        (
            "QdiffRad",
            mesh_.time().timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        mesh_,
        dimensionedScalar("QdiffRad", dimMass/pow3(dimTime), 0.0)
    ),
    QtotalRad_
    (
        IOobject
        (
            "QtotalRad_",
            mesh_.time().timeName(),
            mesh_,
            IOobject::READ_IF_PRESENT,
            IOobject::AUTO_WRITE
        ),
        mesh_,
        dimensionedScalar("QtotalRad_", dimMass/pow3(dimTime), 0.0)
    ),
    hitFaces_(),
    Ru_
    (
        IOobject
        (
            "Ru",
            mesh_.time().timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        mesh_,
        dimensionedScalar("Ru", dimMass/dimLength/pow3(dimTime), 0.0)
    ),
    solarCalc_(dict, mesh_),
    verticalDir_(vector::zero),
    useVFbeamToDiffuse_(false),
    includePatches_(mesh_.boundary().size(), -1),
    coarseToFine_(),
    nBands_(1),
    spectralDistribution_(nBands_),
    map_(),
    visibleFaceFaces_
    (
        IOobject
        (
            "visibleFaceFaces",
            mesh_.facesInstance(),
            mesh_,
            IOobject::READ_IF_PRESENT,
            IOobject::NO_WRITE,
            false
        )
    ),
    solidCoupled_(true),
    wallCoupled_(false),
    absorptivity_(mesh_.boundaryMesh().size()),
    updateAbsorptivity_(false),
    firstIter_(true)
{
    initialise(dict);
}
 
 
// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
Foam::radiation::solarHeatLoad::~solarHeatLoad()
{}
 
 
// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
Foam::label  Foam::radiation::solarHeatLoad::nBands() const
{
    return nBands_;
}
 
 
bool Foam::radiation::solarHeatLoad::read()
{
    if (radiationModel::read())
    {
        return true;
    }
    else
    {
        return false;
    }
}
 
 
void Foam::radiation::solarHeatLoad::calculate()
{
    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
    labelHashSet includePatches;
    forAll(patches, patchI)
    {
        const polyPatch& pp = patches[patchI];
        if (!pp.coupled() && !isA<cyclicAMIPolyPatch>(pp))
        {
            includePatches.insert(patchI);
        }
    }
    //Info << "all IncludePatches the number is " << includePatches.size() << endl;
 
    labelHashSet includeWallGetRadiationPatches;
    forAll(patches, patchI)
    {
        const polyPatch& pp = patches[patchI];
    //  Info << "boundary type is " << pp.type() << endl;
        if  (isA<wallPolyPatch>(pp))
        {
           includeWallGetRadiationPatches.insert(patchI);
        }
    }
 
//  Info << "includeWallGetRadiationPatches the number is " <<includeWallGetRadiationPatches.size() << endl;
 
    if (updateAbsorptivity_ || firstIter_)
    {
        updateAbsorptivity(includePatches);
    }
 
    bool facesChanged = updateHitFaces();
 
//  Info << "faceChanged is = " << facesChanged <<endl;
 
    if (facesChanged)
    {
        // Reset Ru and Qr
        Ru_ = dimensionedScalar("Ru", dimMass/dimLength/pow3(dimTime), 0.0);
        Qr_.boundaryField() = 0.0;
        QdiffRad_.boundaryField() = 0.0;
 
        // Add direct hit radation
        const labelList& hitFacesId = hitFaces_->rayStartFaces();
        updateDirectHitRadiation(hitFacesId,includeWallGetRadiationPatches);
 
        // Add sky diffusive radiation
        updateSkyDiffusiveRadiation
        (
            includePatches,
            includeWallGetRadiationPatches
        );
 
        // Add indirect diffusive radiation
        if (useVFbeamToDiffuse_)
        {
            calculateQdiff(includePatches,includeWallGetRadiationPatches);
        }
 
        firstIter_ = false;
 
//Info << "direct Qr is = " << Qr_ <<endl;
        QtotalRad_ = Qr_ + QdiffRad_;
 
    }
 
    if (debug)
    {
        if (mesh_.time().outputTime())
        {
            Ru_.write();
        }
    }
}
 
 
Foam::tmp<Foam::volScalarField> Foam::radiation::solarHeatLoad::Rp() const
{
    return tmp<volScalarField>
    (
        new volScalarField
        (
            IOobject
            (
                "Rp",
                mesh_.time().timeName(),
                mesh_,
                IOobject::NO_READ,
                IOobject::NO_WRITE,
                false
            ),
            mesh_,
            dimensionedScalar
            (
                "zero",
                dimMass/pow3(dimTime)/dimLength/pow4(dimTemperature),
                0.0
            )
        )
    );
}
 
 
Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
Foam::radiation::solarHeatLoad::Ru() const
{
    return Ru_;
}
 
// ************************************************************************* //