src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -92,7 +92,6 @@ Usage
     spectralDistribution is the energy spectral distribution of the collimated
     external beam.
 
-
 SourceFiles
     fvDOM.C
 
@@ -121,8 +120,7 @@ class fvDOM
 :
     public radiationModel
 {
-    // Private data
-
+    // Private Data
 
         //- Incident radiation  [W/m2]
         volScalarField G_;
@@ -184,6 +182,9 @@ class fvDOM
         //- Spectral energy distribution for the external beam
         scalarList spectralDistribution_;
 
+        //- Time-dependent spectral distributions
+        autoPtr<TimeFunction1<scalarField>> spectralDistributions_;
+
         //- Solar calculator
         autoPtr<solarCalculator> solarCalculator_;
 
@@ -222,10 +223,10 @@ public:
 
 
     //- Destructor
-    virtual ~fvDOM();
+    virtual ~fvDOM() = default;
 
 
-    // Member functions
+    // Member Functions
 
         // Edit
 

src/thermophysicalModels/radiation/radiationModels/solarLoad/solarLoad.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2015 OpenFOAM Foundation
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -37,7 +37,6 @@ License
 #include "wallPolyPatch.H"
 #include "constants.H"
 
-
 using namespace Foam::constant;
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -83,7 +82,7 @@ void Foam::radiation::solarLoad::updateReflectedRays
     {
         if (includePatches[patchID])
         {
-            for (label bandI = 0; bandI < nBands_; bandI++)
+            for (label bandI = 0; bandI < nBands_; ++bandI)
             {
                 qrBf[patchID] +=
                 reflectedFaces_->qreflective(bandI).boundaryField()[patchID];
@@ -94,7 +93,7 @@ void Foam::radiation::solarLoad::updateReflectedRays
             const scalarField& sf = mesh_.magSf().boundaryField()[patchID];
             const labelUList& cellIs = patches[patchID].faceCells();
 
-            for (label bandI = 0; bandI < nBands_; bandI++)
+            for (label bandI = 0; bandI < nBands_; ++bandI)
             {
                 forAll(cellIs, i)
                 {
@@ -109,9 +108,6 @@ void Foam::radiation::solarLoad::updateReflectedRays
             }
         }
     }
-
-
-
 }
 
 
@@ -133,7 +129,7 @@ bool Foam::radiation::solarLoad::updateHitFaces()
             }
             case solarCalculator::mSunDirTracking:
             {
-                label updateIndex = label
+                const label updateIndex = label
                 (
                     mesh_.time().value()/solarCalc_.sunTrackingUpdateInterval()
                 );
@@ -167,7 +163,7 @@ void Foam::radiation::solarLoad::updateAbsorptivity
     for (const label patchID : includePatches)
     {
         absorptivity_[patchID].setSize(nBands_);
-        for (label bandI = 0; bandI < nBands_; bandI++)
+        for (label bandI = 0; bandI < nBands_; ++bandI)
         {
             absorptivity_[patchID][bandI] =
                 boundaryRadiation.absorptivity(patchID, bandI);
@@ -189,7 +185,7 @@ void Foam::radiation::solarLoad::updateDirectHitRadiation
     // Reset qr and qrPrimary
     qrBf = 0.0;
 
-    for (label bandI = 0; bandI < nBands_; bandI++)
+    for (label bandI = 0; bandI < nBands_; ++bandI)
     {
         volScalarField::Boundary& qprimaryBf =
             qprimaryRad_[bandI].boundaryFieldRef();
@@ -199,7 +195,7 @@ void Foam::radiation::solarLoad::updateDirectHitRadiation
         forAll(hitFacesId, i)
         {
             const label faceI = hitFacesId[i];
-            label patchID = patches.whichPatch(faceI);
+            const label patchID = patches.whichPatch(faceI);
             const polyPatch& pp = patches[patchID];
             const label localFaceI = faceI - pp.start();
             const vector qPrim =
@@ -303,7 +299,7 @@ void Foam::radiation::solarLoad::updateSkyDiffusiveRadiation
                     const label cellI = cellIds[faceI];
                     if (includeMappedPatchBasePatches[patchID])
                     {
-                        for (label bandI = 0; bandI < nBands_; bandI++)
+                        for (label bandI = 0; bandI < nBands_; ++bandI)
                         {
                             qrBf[patchID][faceI] +=
                                 (Ed + Er)
@@ -313,7 +309,7 @@ void Foam::radiation::solarLoad::updateSkyDiffusiveRadiation
                     }
                     else
                     {
-                        for (label bandI = 0; bandI < nBands_; bandI++)
+                        for (label bandI = 0; bandI < nBands_; ++bandI)
                         {
                             Ru_[cellI] +=
                                 (Ed + Er)
@@ -328,6 +324,7 @@ void Foam::radiation::solarLoad::updateSkyDiffusiveRadiation
         break;
 
         case solarCalculator::mSunLoadConstant:
+        case solarCalculator::mSunLoadTimeDependent:
         {
             for (const label patchID : includePatches)
             {
@@ -340,7 +337,7 @@ void Foam::radiation::solarLoad::updateSkyDiffusiveRadiation
                     const label cellI = cellIds[faceI];
                     if (includeMappedPatchBasePatches[patchID])
                     {
-                        for (label bandI = 0; bandI < nBands_; bandI++)
+                        for (label bandI = 0; bandI < nBands_; ++bandI)
                         {
                             qrBf[patchID][faceI] +=
                                 solarCalc_.diffuseSolarRad()
@@ -350,7 +347,7 @@ void Foam::radiation::solarLoad::updateSkyDiffusiveRadiation
                     }
                     else
                     {
-                        for (label bandI = 0; bandI < nBands_; bandI++)
+                        for (label bandI = 0; bandI < nBands_; ++bandI)
                         {
                             Ru_[cellI] +=
                                 (
@@ -370,28 +367,24 @@ void Foam::radiation::solarLoad::updateSkyDiffusiveRadiation
 
 void Foam::radiation::solarLoad::initialise(const dictionary& coeffs)
 {
-    coeffs.readEntry("spectralDistribution", spectralDistribution_);
-
-    nBands_ = spectralDistribution_.size();
-
-    qprimaryRad_.setSize(nBands_);
+    spectralDistributions_.reset
+    (
+        new TimeFunction1<scalarField>
+        (
+            mesh_.time(), "spectralDistribution", coeffs
+        )
+    );
 
-    if (coeffs.readIfPresent("gridUp", verticalDir_))
-    {
-         verticalDir_.normalise();
-    }
-    else
-    {
-        const uniformDimensionedVectorField& g =
-            meshObjects::gravity::New(mesh_.time());
-        verticalDir_ = (-g/mag(g)).value();
-    }
+    spectralDistribution_ =
+        spectralDistributions_->value(mesh_.time().value());
 
-    coeffs.readEntry("useReflectedRays", useReflectedRays_);
+    nBands_ = spectralDistribution_.size();
 
     spectralDistribution_ =
         spectralDistribution_/sum(spectralDistribution_);
 
+    qprimaryRad_.setSize(nBands_);
+
     forAll(qprimaryRad_, bandI)
     {
         qprimaryRad_.set
@@ -413,9 +406,21 @@ void Foam::radiation::solarLoad::initialise(const dictionary& coeffs)
         );
     }
 
+    if (coeffs.readIfPresent("gridUp", verticalDir_))
+    {
+         verticalDir_.normalise();
+    }
+    else
+    {
+        const uniformDimensionedVectorField& g =
+            meshObjects::gravity::New(mesh_.time());
+        verticalDir_ = (-g/mag(g)).value();
+    }
+
     coeffs.readIfPresent("solidCoupled", solidCoupled_);
     coeffs.readIfPresent("wallCoupled", wallCoupled_);
     coeffs.readIfPresent("updateAbsorptivity", updateAbsorptivity_);
+    coeffs.readEntry("useReflectedRays", useReflectedRays_);
 }
 
 /*
@@ -721,6 +726,7 @@ void Foam::radiation::solarLoad::calculateQdiff
 Foam::radiation::solarLoad::solarLoad(const volScalarField& T)
 :
     radiationModel(typeName, T),
+    solarCalc_(coeffs_, mesh_),
     dict_(coeffs_),
     qr_
     (
@@ -750,17 +756,18 @@ Foam::radiation::solarLoad::solarLoad(const volScalarField& T)
         mesh_,
         dimensionedScalar(dimMass/dimLength/pow3(dimTime), Zero)
     ),
-    solarCalc_(coeffs_, mesh_),
-    verticalDir_(Zero),
-    useReflectedRays_(false),
+    absorptivity_(mesh_.boundaryMesh().size()),
     spectralDistribution_(),
-    nBands_(0),
+    spectralDistributions_(),
     qprimaryRad_(0),
+    verticalDir_(Zero),
+    nBands_(0),
+    updateTimeIndex_(0),
     solidCoupled_(true),
-    absorptivity_(mesh_.boundaryMesh().size()),
+    wallCoupled_(false),
     updateAbsorptivity_(false),
-    firstIter_(true),
-    updateTimeIndex_(0)
+    useReflectedRays_(false),
+    firstIter_(true)
 {
     initialise(coeffs_);
 }
@@ -773,6 +780,7 @@ Foam::radiation::solarLoad::solarLoad
 )
 :
     radiationModel(typeName, dict, T),
+    solarCalc_(dict, mesh_),
     dict_(dict),
     qr_
     (
@@ -802,18 +810,18 @@ Foam::radiation::solarLoad::solarLoad
         mesh_,
         dimensionedScalar(dimMass/dimLength/pow3(dimTime), Zero)
     ),
-    solarCalc_(dict, mesh_),
-    verticalDir_(Zero),
-    useReflectedRays_(false),
+    absorptivity_(mesh_.boundaryMesh().size()),
     spectralDistribution_(),
-    nBands_(0),
+    spectralDistributions_(),
     qprimaryRad_(0),
+    verticalDir_(Zero),
+    nBands_(0),
+    updateTimeIndex_(0),
     solidCoupled_(true),
     wallCoupled_(false),
-    absorptivity_(mesh_.boundaryMesh().size()),
     updateAbsorptivity_(false),
-    firstIter_(true),
-    updateTimeIndex_(0)
+    useReflectedRays_(false),
+    firstIter_(true)
 {
     initialise(dict);
 }
@@ -821,12 +829,6 @@ Foam::radiation::solarLoad::solarLoad
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::label  Foam::radiation::solarLoad::nBands() const
-{
-    return nBands_;
-}
-
-
 bool Foam::radiation::solarLoad::read()
 {
     if (radiationModel::read())
@@ -870,13 +872,25 @@ void Foam::radiation::solarLoad::calculate()
         updateAbsorptivity(includePatches);
     }
 
-    bool facesChanged = updateHitFaces();
+    const bool facesChanged = updateHitFaces();
 
-    if (facesChanged)
+    const bool timeDependentLoad =
+        solarCalc_.sunLoadModel() == solarCalculator::mSunLoadTimeDependent;
+
+    if (facesChanged || timeDependentLoad)
     {
         // Reset Ru
         Ru_ = dimensionedScalar("Ru", dimMass/dimLength/pow3(dimTime), Zero);
 
+        solarCalc_.correctDirectSolarRad();
+        solarCalc_.correctDiffuseSolarRad();
+
+        spectralDistribution_ =
+            spectralDistributions_->value(mesh_.time().value());
+
+        spectralDistribution_ =
+            spectralDistribution_/sum(spectralDistribution_);
+
         // Add direct hit radiation
         const labelList& hitFacesId = hitFaces_->rayStartFaces();
         updateDirectHitRadiation(hitFacesId, includeMappedPatchBasePatches);

src/thermophysicalModels/radiation/radiationModels/solarLoad/solarLoad.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2015 OpenFOAM Foundation
-    Copyright (C) 2018 OpenCFD Ltd.
+    Copyright (C) 2018-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -31,7 +31,7 @@ Group
     grpRadiationModels
 
 Description
-    The solar load radiation model includes Sun primary hits, their
+    The \c solarLoad radiation model includes Sun primary hits, their
     reflective fluxes and diffusive sky radiative fluxes.
 
     The primary hit rays are calculated using a face shading algorithm.
@@ -43,14 +43,52 @@ Description
     Fair Weather Conditions Method from the ASHRAE Handbook.
 
     By default the energy is included in cells adjacent to the patches into
-    the energy Equation (wallCoupled = false). On coupled patches the flux is
+    the energy equation (\c wallCoupled=false). On coupled patches the flux is
     by default added to the wall and considered into the solid
-    (solidCoupled = true).
-
-    The solarLoad model can be used in conjuntion with fvDOM and viewFactor
-    radiation models. The flag useSolarLoad must be true on the rediation
-    dictionary.
-
+    (\c solidCoupled=true).
+
+    The \c solarLoad model can be used in conjuntion with \c fvDOM and
+    \c viewFactor radiation models. The flag \c useSolarLoad must be
+    \c true on the \c radiationProperties dictionary.
+
+Usage
+    Minimal examples by using \c constant/radiationProperties:
+
+    \verbatim
+    solarLoadCoeffs
+    {
+        // Mandatory entries
+        useReflectedRays            true;
+        spectralDistribution        (1 5 1 2);
+
+        // Optional entries
+        solidCoupled                true;
+        wallCoupled                 false;
+        updateAbsorptivity          true;
+
+        // Mandatory/Optional (inherited) entries
+        ...
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property     | Description                      | Type  | Reqd   | Deflt
+      useReflectedRays | Flag to use reflected rays   | bool  | yes    | -
+      spectralDistribution | Spectral distribution for the integrated <!--
+                   --> solar heat flux | TimeFunction1\<scalarField\> | yes | -
+      solidCoupled     | Flag to couple solids through mapped         <!--
+                   --> boundary patch using qr        | bool  | no     | true
+      wallCoupled      | Flag to couple wall patches using qr         <!--
+                   -->                                | bool  | no     | false
+      updateAbsorptivity | Flag to enable absorptivity updates        <!--
+                   -->                                | bool  | no     | false
+    \endtable
+
+    The inherited entries are elaborated in:
+     - \link radiationModel.H \endlink
+     - \link solarCalculator.H \endlink
+     - \link TimeFunction1.H \endlink
 
 SourceFiles
     solarLoad.C
@@ -65,6 +103,7 @@ SourceFiles
 #include "faceShading.H"
 #include "faceReflecting.H"
 #include "solarCalculator.H"
+#include "TimeFunction1.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -81,10 +120,10 @@ class solarLoad
 :
     public radiationModel
 {
+    // Private Data
 
-private:
-
-    // Private data
+        //- Solar calculator
+        solarCalculator solarCalc_;
 
         //- Dictionary
         dictionary dict_;
@@ -102,46 +141,46 @@ private:
         //- and wallCoupled false
         DimensionedField<scalar, volMesh> Ru_;
 
-        //- Solar calculator
-        solarCalculator solarCalc_;
-
-        //- Vertical direction (Default is g vector)
-        vector verticalDir_;
-
-        //- Include reflected rays from specular surfaces
-        bool useReflectedRays_;
+        //- Absorptivity list
+        List<List<tmp<scalarField>>> absorptivity_;
 
         //- Spectral distribution for the integrated solar heat flux
         scalarList spectralDistribution_;
 
-        //-Number of bands
-        label nBands_;
+        //- Time-dependent spectral distributions
+        autoPtr<TimeFunction1<scalarField>> spectralDistributions_;
 
         //- Primary solar radiative heat flux per band [W/m2]
         PtrList<volScalarField> qprimaryRad_;
 
-        //- Couple solids through mapped boundary patch using qr (default:true)
+        //- Vertical direction (default is g vector)
+        vector verticalDir_;
+
+        //- Number of bands
+        label nBands_;
+
+        //- Update Sun position index
+        label updateTimeIndex_;
+
+        //- Couple solids through mapped boundary patch using qr
         bool solidCoupled_;
 
-        //- Couple wall patches using qr (default:false)
+        //- Couple wall patches using qr
         bool wallCoupled_;
 
-        //- Absorptivity list
-        List<List<tmp<scalarField>>> absorptivity_;
-
         //- Update absorptivity
         bool updateAbsorptivity_;
 
+        //- Include reflected rays from specular surfaces
+        bool useReflectedRays_;
+
         //- First iteration
         bool firstIter_;
 
-        //- Update Sun position index
-        label updateTimeIndex_;
-
 
     // Private Member Functions
 
-        //- Initialise
+        //- Initialise model parameters
         void initialise(const dictionary&);
 
         //- Update direct hit faces radiation
@@ -166,6 +205,7 @@ private:
         //- Update absorptivity
         void updateAbsorptivity(const labelHashSet& includePatches);
 
+
         //- No copy construct
         solarLoad(const solarLoad&) = delete;
 
@@ -192,15 +232,18 @@ public:
     virtual ~solarLoad() = default;
 
 
-    // Member functions
+    // Member Functions
 
-        // Edit
+        // Evaluation
 
-            //- Solve
+            //- Read radiationProperties dictionary
+            bool read();
+
+            //- Solve radiation equations
             void calculate();
 
-            //- Read radiation properties dictionary
-            bool read();
+
+        // Access
 
             //- Source term component (for power of T^4)
             virtual tmp<volScalarField> Rp() const;
@@ -208,17 +251,17 @@ public:
             //- Source term component (constant)
             virtual tmp<DimensionedField<scalar, volMesh>> Ru() const;
 
-
-    // Access
-
-        //- Number of bands
-        label nBands() const;
-
-        //- Primary solar heat flux
-        const volScalarField& qprimaryRad(const label bandI) const
-        {
-            return qprimaryRad_[bandI];
-        }
+            //- Return const access to the number of bands
+            label nBands() const noexcept
+            {
+                return nBands_;
+            }
+
+            //- Return const access to the primary solar heat flux
+            const volScalarField& qprimaryRad(const label bandI) const
+            {
+                return qprimaryRad_[bandI];
+            }
 };
 
 

src/thermophysicalModels/radiation/submodels/solarCalculator/solarCalculator.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2015-2020 OpenCFD Ltd.
+    Copyright (C) 2015-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -46,8 +46,12 @@ const Foam::Enum
 >
 Foam::solarCalculator::sunDirectionModelTypeNames_
 ({
+    { sunDirModel::mSunDirConstant, "constant" },
+    { sunDirModel::mSunDirTracking, "tracking" },
+
+    // old long names (v2012 and earlier)
     { sunDirModel::mSunDirConstant, "sunDirConstant" },
-    { sunDirModel::mSunDirTracking, "sunDirTracking" },
+    { sunDirModel::mSunDirTracking, "sunDirTracking" }
 });
 
 
@@ -55,14 +59,20 @@ const Foam::Enum
 <
     Foam::solarCalculator::sunLModel
 >
-Foam::solarCalculator::sunLoadModelTypeNames_
+Foam::solarCalculator::sunLModelTypeNames_
 ({
+    { sunLModel::mSunLoadConstant, "constant" },
+    { sunLModel::mSunLoadTimeDependent, "timeDependent" },
+    { sunLModel::mSunLoadFairWeatherConditions, "fairWeather" },
+    { sunLModel::mSunLoadTheoreticalMaximum, "theoreticalMaximum" },
+
+    // old long names (v2012 and earlier)
     { sunLModel::mSunLoadConstant, "sunLoadConstant" },
     {
         sunLModel::mSunLoadFairWeatherConditions,
         "sunLoadFairWeatherConditions"
     },
-    { sunLModel::mSunLoadTheoreticalMaximum, "sunLoadTheoreticalMaximum" },
+    { sunLModel::mSunLoadTheoreticalMaximum, "sunLoadTheoreticalMaximum" }
 });
 
 
@@ -70,41 +80,34 @@ Foam::solarCalculator::sunLoadModelTypeNames_
 
 void Foam::solarCalculator::calculateBetaTheta()
 {
-    scalar runTime = 0.0;
-    switch (sunDirectionModel_)
+    scalar runTime = 0;
+
+    if (sunDirectionModel_ == mSunDirTracking)
     {
-        case mSunDirTracking:
-        {
-            runTime = mesh_.time().value();
-            break;
-        }
-        case mSunDirConstant:
-        {
-            break;
-        }
+        runTime = mesh_.time().value();
     }
 
-    scalar LSM = 15.0*(dict_.get<scalar>("localStandardMeridian"));
+    const scalar LSM = 15.0*(dict_.get<scalar>("localStandardMeridian"));
 
-    scalar D = dict_.get<scalar>("startDay") + runTime/86400.0;
-    scalar M = 6.24004 + 0.0172*D;
-    scalar EOT = -7.659*sin(M) + 9.863*sin(2*M + 3.5932);
+    const scalar D = dict_.get<scalar>("startDay") + runTime/86400.0;
+    const scalar M = 6.24004 + 0.0172*D;
+    const scalar EOT = -7.659*sin(M) + 9.863*sin(2*M + 3.5932);
 
     dict_.readEntry("startTime", startTime_);
 
-    scalar LST =  startTime_ + runTime/3600.0;
+    const scalar LST =  startTime_ + runTime/3600.0;
 
-    scalar LON = dict_.get<scalar>("longitude");
+    const scalar LON = dict_.get<scalar>("longitude");
 
-    scalar AST = LST + EOT/60.0 + (LON - LSM)/15;
+    const scalar AST = LST + EOT/60.0 + (LON - LSM)/15;
 
-    scalar delta = 23.45*sin(degToRad((360*(284 + D))/365));
+    const scalar delta = 23.45*sin(degToRad((360*(284 + D))/365));
 
-    scalar H = degToRad(15*(AST - 12));
+    const scalar H = degToRad(15*(AST - 12));
 
-    scalar L = degToRad(dict_.get<scalar>("latitude"));
+    const scalar L = degToRad(dict_.get<scalar>("latitude"));
 
-    scalar deltaRad = degToRad(delta);
+    const scalar deltaRad = degToRad(delta);
     beta_ = max(asin(cos(L)*cos(deltaRad)*cos(H) + sin(L)*sin(deltaRad)), 1e-3);
     theta_ = acos((sin(beta_)*sin(L) - sin(deltaRad))/(cos(beta_)*cos(L)));
 
@@ -150,7 +153,7 @@ void Foam::solarCalculator::calculateSunDirection()
 }
 
 
-void Foam::solarCalculator::init()
+void Foam::solarCalculator::initialise()
 {
     switch (sunDirectionModel_)
     {
@@ -165,7 +168,6 @@ void Foam::solarCalculator::init()
                 calculateBetaTheta();
                 calculateSunDirection();
             }
-
             break;
         }
         case mSunDirTracking:
@@ -197,6 +199,32 @@ void Foam::solarCalculator::init()
             dict_.readEntry("diffuseSolarRad", diffuseSolarRad_);
             break;
         }
+        case mSunLoadTimeDependent:
+        {
+            directSolarRads_.reset
+            (
+                new TimeFunction1<scalar>
+                (
+                    mesh_.time(),
+                    "directSolarRad",
+                    dict_
+                )
+            );
+
+            diffuseSolarRads_.reset
+            (
+                new TimeFunction1<scalar>
+                (
+                    mesh_.time(),
+                    "diffuseSolarRad",
+                    dict_
+                )
+            );
+
+            directSolarRad_ = directSolarRads_->value(mesh_.time().value());
+            diffuseSolarRad_ = diffuseSolarRads_->value(mesh_.time().value());
+            break;
+        }
         case mSunLoadFairWeatherConditions:
         {
             dict_.readIfPresent
@@ -207,7 +235,8 @@ void Foam::solarCalculator::init()
 
             dict_.readEntry("A", A_);
             dict_.readEntry("B", B_);
-
+            dict_.readEntry("C", C_);
+            dict_.readEntry("groundReflectivity", groundReflectivity_);
             if (!dict_.readIfPresent("beta", beta_))
             {
                 calculateBetaTheta();
@@ -216,17 +245,16 @@ void Foam::solarCalculator::init()
             directSolarRad_ =
                 (1.0 - 0.75*pow(skyCloudCoverFraction_, 3.0))
               * A_/exp(B_/sin(beta_));
-
-            dict_.readEntry("groundReflectivity", groundReflectivity_);
             break;
         }
         case mSunLoadTheoreticalMaximum:
         {
             dict_.readEntry("Setrn", Setrn_);
             dict_.readEntry("SunPrime", SunPrime_);
-            directSolarRad_ = Setrn_*SunPrime_;
-
             dict_.readEntry("groundReflectivity", groundReflectivity_);
+            dict_.readEntry("C", C_);
+
+            directSolarRad_ = Setrn_*SunPrime_;
             break;
         }
     }
@@ -243,26 +271,32 @@ Foam::solarCalculator::solarCalculator
 :
     mesh_(mesh),
     dict_(dict),
-    direction_(Zero),
-    directSolarRad_(0.0),
-    diffuseSolarRad_(0.0),
-    groundReflectivity_(0.0),
-    A_(0.0),
-    B_(0.0),
-    beta_(0.0),
-    theta_(0.0),
-    skyCloudCoverFraction_(0.0),
-    Setrn_(0.0),
-    SunPrime_(0.0),
-    C_(dict.get<scalar>("C")),
     sunDirectionModel_
     (
         sunDirectionModelTypeNames_.get("sunDirectionModel", dict)
     ),
-    sunLoadModel_(sunLoadModelTypeNames_.get("sunLoadModel", dict)),
-    coord_()
+    sunLoadModel_(sunLModelTypeNames_.get("sunLoadModel", dict)),
+    direction_(Zero),
+    sunTrackingUpdateInterval_(0),
+    startTime_(0),
+    gridUp_(Zero),
+    eastDir_(Zero),
+    coord_(),
+    directSolarRad_(0),
+    diffuseSolarRad_(0),
+    directSolarRads_(),
+    diffuseSolarRads_(),
+    skyCloudCoverFraction_(0),
+    groundReflectivity_(0),
+    A_(0),
+    B_(0),
+    beta_(0),
+    theta_(0),
+    C_(0.058),
+    Setrn_(0),
+    SunPrime_(0)
 {
-    init();
+    initialise();
 }
 
 
@@ -270,19 +304,29 @@ Foam::solarCalculator::solarCalculator
 
 void Foam::solarCalculator::correctSunDirection()
 {
-    switch (sunDirectionModel_)
+    if (sunDirectionModel_ == mSunDirTracking)
     {
-        case mSunDirConstant:
-        {
-            break;
-        }
-        case mSunDirTracking:
-        {
-            calculateBetaTheta();
-            calculateSunDirection();
-            directSolarRad_ = A_/exp(B_/sin(max(beta_, ROOTVSMALL)));
-            break;
-        }
+        calculateBetaTheta();
+        calculateSunDirection();
+        directSolarRad_ = A_/exp(B_/sin(max(beta_, ROOTVSMALL)));
+    }
+}
+
+
+void Foam::solarCalculator::correctDirectSolarRad()
+{
+    if (sunLoadModel_ == mSunLoadTimeDependent)
+    {
+        directSolarRad_ = directSolarRads_->value(mesh_.time().value());
+    }
+}
+
+
+void Foam::solarCalculator::correctDiffuseSolarRad()
+{
+    if (sunLoadModel_ == mSunLoadTimeDependent)
+    {
+        diffuseSolarRad_ = diffuseSolarRads_->value(mesh_.time().value());
     }
 }
 

src/thermophysicalModels/radiation/submodels/solarCalculator/solarCalculator.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2015-2017 OpenCFD Ltd.
+    Copyright (C) 2015-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,58 +27,188 @@ Class
     Foam::solarCalculator
 
 Description
-    The solar calculator model provides information about the Sun direction
-    and Sun load model. The available models are:
-
-    For the Sun direction:
-    1) SunDirConstant : the direction is given in 'sunDirection'
-    2) SunDirTracking : the direction is calculated from the following
-       parameters:
-            localStandardMeridian : GMT (Local Zone Meridian) in hours
-            startDay :  day from 1 to 365)
-            startTime:  in hours
-            longitude:  in degrees
-            latitude:   in degrees
-            gridUp:     grid orientation upwards
-            gridEast    grid orientation eastwards
-
-        This model should be use in transient calculations.
-        The keyword 'sunTrackingUpdateInterval' (in hours) specifies on which
-        interval is the Sun direction updated.
-
-
-    Solar Load models available:
-    1) SunLoadConstant: direct and diffusive heat fluxes are provided by the
-       entries 'directSolarRad' and 'diffuseSolarRad'
-
-    2) SunLoadFairWeatherConditions: The solar fluxes are calculated following
-       the Fair Weather Conditions Method from the ASHRAE Handbook. The entries
-       are:
-            skyCloudCoverFraction: Fraction of sky covered by clouds (0-1)
-            A :     Apparent solar irradiation at air mass m = 0
-            B :     Atmospheric extinction coefficient
-            beta:   Solar altitude (in degrees) above the horizontal. This
-            can be read or calculated providing the respective parameters
-            for Sun position explained above.
-            groundReflectivity : ground reflectivity
-
-            In this model the flux is calculated as:
-
-                directSolarRad =
-                    (1 - 0.75*skyCloudCoverFraction^3)*A/exp(B/sin(beta));
-
-    3) SunLoadTheoreticalMaximum: The entries are:
-            Setrn
-            SunPrime:
-            groundReflectivity : ground reflectivity
-
-            In this model the flux is calculated as:
-
-                directSolarRad = Setrn*SunPrime;
-
-    The diffuse on vertical/horizontal walls and ground-reflected radiation are
+    A solar calculator model providing models
+    for the solar direction and solar loads.
+
+    Available models for the solar direction:
+    - \c constant: Constant sunbeam direction.
+    - \c tracking: Transient model calculating sunbeam direction
+    based on a given set of parameters.
+
+    Available models for the solar load:
+    - \c constant: Constant solar load.
+    - \c timeDependent: Time-dependent solar load.
+    - \c fairWeather: Solar fluxes are calculated following
+       the "Fair Weather Conditions Method from the ASHRAE Handbook".
+    - \c theoreticalMaximum: Theoretically maximum solar load.
+
+Usage
+    Minimal examples by using \c constant/radiationProperties:
+
+    \c sunDirectionModel - Option-1:
+    \verbatim
+    solarLoadCoeffs
+    {
+        sunDirectionModel           constant;
+        sunDirection                (1 0 -1);
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property     | Description                  | Type   | Reqd | Deflt
+      sunDirection | Sunbeam direction            | vector | no   | calculated
+    \endtable
+
+
+    \c sunDirectionModel - Option-2:
+    \verbatim
+    solarLoadCoeffs
+    {
+        sunDirectionModel           tracking;
+        sunTrackingUpdateInterval   800;
+        localStandardMeridian       9;
+        startDay                    204;
+        startTime                   15;
+        longitude                   139.74;
+        latitude                    35.658;
+        gridUp                      (0 0 1);
+        gridEast                    (1 0 0);
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property     | Description                       | Type   | Reqd | Deflt
+      sunTrackingUpdateInterval  | Interval to update the Sun direction  <!--
+                   --> [decimal hours]                 | scalar | yes  | -
+      localStandardMeridian | GMT (Local Zone Meridian) [decimal hours]<!--
+                   -->                                 | scalar | yes  | -
+      startDay     | Day from 1 to 365                 | scalar | yes  | -
+      startTime    | Start time for the Sun position [decimal hours]   <!--
+                   -->                                 | scalar | yes  | -
+      longitude    | Geographic coordinate specifying the east–west <!--
+                   --> position of a point on the surface of a planetary <!--
+                   --> body [degree]                   | scalar | yes  | -
+      latitude     | Geographic coordinate specifying the north–south  <!--
+                   --> position of a point on the surface of a planetary <!--
+                   --> body [degree]                   | scalar | yes  | -
+      gridUp       | Grid orientation upwards          | vector | yes  | -
+      gridEast     | Grid orientation eastwards        | vector | yes  | -
+    \endtable
+
+
+    \c sunLoadModel - Option-1:
+    \verbatim
+    solarLoadCoeffs
+    {
+        sunLoadModel                constant;
+        directSolarRad              100;
+        diffuseSolarRad             0;
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property     | Description                       | Type   | Reqd | Deflt
+      directSolarRad | Direct solar irradiation [W/m2] | scalar | yes  | -
+      diffuseSolarRad | Diffuse solar irradiation on vertical surfaces <!--
+                      --> [W/m2]                       | scalar | yes  | -
+    \endtable
+
+
+    \c sunLoadModel - Option-2:
+    \verbatim
+    solarLoadCoeffs
+    {
+        sunLoadModel                timeDependent;
+        directSolarRad              <TimeFunction1<scalar>>;
+        diffuseSolarRad             <TimeFunction1<scalar>>;
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property     | Description                       | Type    | Reqd | Deflt
+      directSolarRad  | Time-series of direct solar irradiation        <!--
+                   --> [W/m2]          | TimeFunction1\<scalar\> | yes  | -
+      diffuseSolarRad | Time-series of diffuse solar irradiation on    <!--
+                   --> vertical surfaces [W/m2]                        <!--
+                   -->                  | TimeFunction1\<scalar\> | yes | -
+    \endtable
+
+
+    \c sunLoadModel - Option-3:
+    \verbatim
+    solarLoadCoeffs
+    {
+        sunLoadModel                fairWeather;
+        skyCloudCoverFraction       0.25;
+        groundReflectivity          1.0;
+        A                           0.1;
+        B                           0.2;
+        C                           0.058;
+        beta                        0.15;
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property  | Description                       | Type   | Reqd | Deflt
+      A         | Apparent solar irradiation at air mass m = 0      <!--
+                -->                                 | scalar | yes  | -
+      B         | Atmospheric extinction coefficient                <!--
+                -->                                 | scalar | yes  | -
+      C         | Solar diffusivity constant        | scalar | yes  | -
+      groundReflectivity | Ground reflectivity      | scalar | yes  | -
+      skyCloudCoverFraction | Fraction of sky covered by clouds [0,1]  <!--
+                -->                                 | scalar | no   | 0
+      beta      | Solar altitude (in degrees) above the horizontal  <!--
+                -->                                 | scalar | no   | calculated
+    \endtable
+
+    In this model the flux is calculated as:
+    \verbatim
+    directSolarRad = (1 - 0.75*skyCloudCoverFraction^3)*A/exp(B/sin(beta));
+    \endverbatim
+
+
+    \c sunLoadModel - Option-4:
+    \verbatim
+    solarLoadCoeffs
+    {
+        sunLoadModel                theoreticalMaximum;
+        Setrn                       1.0;
+        SunPrime                    4.0;
+        groundReflectivity          1.0;
+        C                           0.058;
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property     | Description                       | Type   | Reqd | Deflt
+      Setrn        | Parameter in maximum theoretical direct solar    <!--
+                   --> model                           | scalar | yes  | -
+      SunPrime     | Parameter in maximum theoretical direct solar    <!--
+                   --> model                           | scalar | yes  | -
+      groundReflectivity | Ground reflectivity         | scalar | yes  | -
+      C            | Solar diffusivity constant        | scalar | yes  | -
+    \endtable
+
+    In this model the flux is calculated as:
+    \verbatim
+    directSolarRad = Setrn*SunPrime;
+    \endverbatim
+
+Note
+  - The \c sunDirectionModel:tracking can only be used
+  in transient calculations.
+  - The keyword \c sunTrackingUpdateInterval (in hours) specifies on which
+    interval is the Sun direction updated.
+  - The diffuse on vertical/horizontal walls and ground-reflected radiation are
     calculated following the ASHRAE Handbook.
-
+  - The range of \c skyCloudCoverFraction is [0,1].
 
 SourceFiles
     solarCalculator.C
@@ -93,6 +223,7 @@ SourceFiles
 #include "DynamicField.H"
 #include "HashSet.H"
 #include "coordinateSystem.H"
+#include "TimeFunction1.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -107,37 +238,34 @@ class solarCalculator
 {
 public:
 
-    //  Public enumeration
+    //  Public Enumeration
 
-        //- Sun direction models
+        //- Options for the Sun direction models
         enum sunDirModel
         {
             mSunDirConstant,
             mSunDirTracking
         };
 
-        //- Direct sun load models
+        //- Names for sunDirModel
+        static const Enum<sunDirModel> sunDirectionModelTypeNames_;
+
+        //- Options for the Sun load models
         enum sunLModel
         {
             mSunLoadConstant,
+            mSunLoadTimeDependent,
             mSunLoadFairWeatherConditions,
             mSunLoadTheoreticalMaximum
         };
 
-
-protected:
-
-         //- Sun direction models
-        static const Enum<sunDirModel> sunDirectionModelTypeNames_;
-
-        //- Sun load models
-        static const Enum<sunLModel> sunLoadModelTypeNames_;
+        //- Names for sunLModel
+        static const Enum<sunLModel> sunLModelTypeNames_;
 
 
 private:
 
-    // Private data
-
+    // Private Data
 
         //- Reference to mesh
         const fvMesh& mesh_;
@@ -145,56 +273,78 @@ private:
         //- Dictionary
         dictionary dict_;
 
-        //- Direction
-        vector direction_;
+        //- Sun direction model
+        sunDirModel sunDirectionModel_;
 
-        //- Direct solar irradiation
-        scalar directSolarRad_;
+        //- Sun load model
+        sunLModel sunLoadModel_;
 
-        //- Diffuse solar irradiation on vertical surfaces
-        scalar diffuseSolarRad_;
 
-        //- Ground reflectivity
-        scalar groundReflectivity_;
+        // sunDirectionModel = constant
 
-        //- Fair weather direct solar load model parameters
-        scalar A_;
-        scalar B_;
-        scalar beta_;
-        scalar theta_;
+            //- Sunbeam direction
+            vector direction_;
 
-        //- Sky cloud cover fraction [0-1]
-        scalar skyCloudCoverFraction_;
 
+        // sunDirectionModel = tracking
 
-        //- Maximum theoretical direct solar load model parameters
-        scalar Setrn_;
-        scalar SunPrime_;
+            //- Interval to update the Sun direction [decimal hours]
+            scalar sunTrackingUpdateInterval_;
 
+            //- Start time for the Sun position [decimal hours]
+            scalar startTime_;
 
-        //- Diffusive solar load model parameters
-        scalar C_;
+            //- Up grid orientation
+            vector gridUp_;
 
-        //- Sun direction model
-        sunDirModel sunDirectionModel_;
+            //- East grid orientation
+            vector eastDir_;
 
-        //- Sun load model
-        sunLModel sunLoadModel_;
+            //- Grid coordinate system
+            autoPtr<coordinateSystem> coord_;
+
+
+        // sunLoadModel = constant
+
+            //- Direct solar irradiation
+            scalar directSolarRad_;
+
+            //- Diffuse solar irradiation on vertical surfaces
+            scalar diffuseSolarRad_;
+
+
+        // sunLoadModel = timeDependent
+
+            //- Time-series of direct solar irradiation
+            autoPtr<TimeFunction1<scalar>> directSolarRads_;
+
+            //- Time-series of diffuse solar irradiation on vertical surfaces
+            autoPtr<TimeFunction1<scalar>> diffuseSolarRads_;
+
+
+        // sunLoadModel = fairWeather
+
+            //- Sky cloud cover fraction [0-1]
+            scalar skyCloudCoverFraction_;
 
-        //- Grid coordinate system
-        autoPtr<coordinateSystem> coord_;
+            //- Ground reflectivity
+            scalar groundReflectivity_;
 
-        //- East grid orientation
-        vector eastDir_;
+            //- Fair weather direct solar load model parameters
+            scalar A_;
+            scalar B_;
+            scalar beta_;
+            scalar theta_;
 
-        //- Up grid orientation
-        vector gridUp_;
+            //- Diffusive solar load model parameter
+            scalar C_;
 
-        //- Interval in decimal hours to update Sun direction for SunDirTracking
-        scalar sunTrackingUpdateInterval_;
 
-        //- Start time for the Sun position (decimal hours)
-        scalar startTime_;
+        // sunLoadModel = theoreticalMaximum
+
+            //- Maximum theoretical direct solar load model parameters
+            scalar Setrn_;
+            scalar SunPrime_;
 
 
         //- No copy construct
@@ -204,15 +354,15 @@ private:
         void operator=(const solarCalculator&) = delete;
 
 
-    // Private members
+    // Private Member Functions
 
-        //- Init
-        void init();
+        //- Initialise model parameters
+        void initialise();
 
         //- Calculate beta and theta angles
         void calculateBetaTheta();
 
-        //- Calculate Sun direction
+        //- Calculate the Sun direction
         void calculateSunDirection();
 
 
@@ -234,101 +384,109 @@ public:
 
     // Member Functions
 
+        // Evaluation
+
+            //- Correct the Sun direction
+            void correctSunDirection();
+
+            //- Correct direct solar irradiation
+            void correctDirectSolarRad();
+
+            //- Correct diffuse solar irradiation
+            void correctDiffuseSolarRad();
+
+
         // Access
 
-            //- const access to direction
-            const vector& direction() const
+            //- Return const access to the Sun direction model
+            const sunDirModel& sunDirectionModel() const noexcept
             {
-                return direction_;
+                return sunDirectionModel_;
+            }
+
+            //- Return const access to the Sun load model
+            const sunLModel& sunLoadModel() const noexcept
+            {
+                return sunLoadModel_;
             }
 
-            //- Non-const access to direction
+            //- Return non-const access to the Sun direction
             vector& direction()
             {
                 return direction_;
             }
 
-            //- Return direct solar irradiation
+            //- Return const access to the Sun direction
+            const vector& direction() const noexcept
+            {
+                return direction_;
+            }
+
+            //- Return non-const access to the direct solar irradiation
             scalar& directSolarRad()
             {
                 return directSolarRad_;
             }
 
-            //- Return const access to direct solar irradiation
-            const scalar& directSolarRad() const
+            //- Return const access to the direct solar irradiation
+            const scalar& directSolarRad() const noexcept
             {
                 return directSolarRad_;
             }
 
-            //- Return diffuse solar irradiation
+            //- Return non-const access to the diffuse solar irradiation
             scalar& diffuseSolarRad()
             {
                 return diffuseSolarRad_;
             }
 
-            //- Return diffuse solar irradiation
-            const scalar& diffuseSolarRad() const
+            //- Return const access to the diffuse solar irradiation
+            const scalar& diffuseSolarRad() const noexcept
             {
                 return diffuseSolarRad_;
             }
 
-            //- Return C constant
-            scalar C()
+            //- Return const access to the C constant
+            scalar C() const noexcept
             {
                 return C_;
             }
 
-            //- Return beta
-            scalar beta()
+            //- Return const access to beta
+            scalar beta() const noexcept
             {
                 return beta_;
             }
 
-            //- Return theta
-            scalar theta()
+            //- Return const access to theta
+            scalar theta() const noexcept
             {
                 return theta_;
             }
 
-            //- Return Sun direction model
-            sunDirModel sunDirectionModel() const
-            {
-                return sunDirectionModel_;
-            }
-
-            //- Return Sun load model
-            sunLModel sunLoadModel() const
-            {
-                return sunLoadModel_;
-            }
-
-            //- Return groundReflectivity
-            scalar groundReflectivity()
+            //- Return const access to the ground reflectivity
+            scalar groundReflectivity() const noexcept
             {
                 return groundReflectivity_;
             }
 
-            //- Return coordinateSystem
-            const coordinateSystem& coord()
+            //- Return const access to the coordinate system
+            const coordinateSystem& coord() const noexcept
             {
                 return *coord_;
             }
 
-            //- Return sunTrackingUpdateInterval
-            scalar sunTrackingUpdateInterval()
+            //- Return const access to sunTrackingUpdateInterval
+            scalar sunTrackingUpdateInterval() const noexcept
             {
                 return sunTrackingUpdateInterval_;
             }
 
-            //- Return startTime
-            scalar startTime()
+            //- Return const access to startTime
+            scalar startTime() const noexcept
             {
                 return startTime_;
             }
-
-
-        //- Recalculate
-        void correctSunDirection();
 };
 
 

src/transportModels/geometricVoF/advectionSchemes/isoAdvection/isoAdvectionTemplates.C

@@ -249,8 +249,13 @@ void Foam::isoAdvection::boundFlux
 
             // First try to pass surplus fluid on to neighbour cells that are
             // not filled and to which dVf < phi*dt
-            while (mag(alphaOvershoot) > aTol && nFacesToPassFluidThrough > 0)
+            for (label iter=0; iter<10; iter++)
             {
+                if(mag(alphaOvershoot) < aTol || nFacesToPassFluidThrough == 0)
+                {
+                    break;
+                }
+
                 DebugInfo
                     << "\n\nBounding cell " << celli
                     << " with alpha overshooting " << alphaOvershoot

src/transportModels/geometricVoF/reconstructionSchemes/plicSchemes/plicRDF/plicRDF.C

@@ -233,8 +233,7 @@ void Foam::reconstruction::plicRDF::setInitNormals(bool interpolate)
 
 void Foam::reconstruction::plicRDF::calcResidual
 (
-    Map<scalar>& normalResidual,
-    Map<scalar>& avgAngle
+    List<normalRes>& normalResidual
 )
 {
     exchangeFields_.setUpCommforZone(interfaceCell_,false);
@@ -246,7 +245,6 @@ void Foam::reconstruction::plicRDF::calcResidual
 
     const labelListList& stencil = exchangeFields_.getStencil();
 
-    normalResidual.clear();
 
     forAll(interfaceLabels_, i)
     {
@@ -260,13 +258,13 @@ void Foam::reconstruction::plicRDF::calcResidual
         scalar maxDiffNormal = GREAT;
         scalar weight= 0;
         const vector cellNormal = normal_[celli]/mag(normal_[celli]);
-
-        for (const label gblIdx : stencil[celli])
+        forAll(stencil[celli],j)
         {
+            const label gblIdx = stencil[celli][j];
             vector normal =
                 exchangeFields_.getValue(normal_, mapNormal, gblIdx);
 
-            if (mag(normal) != 0 && i != 0)
+            if (mag(normal) != 0 && j != 0)
             {
                 vector n = normal/mag(normal);
                 scalar cosAngle = max(min((cellNormal & n), 1), -1);
@@ -291,8 +289,9 @@ void Foam::reconstruction::plicRDF::calcResidual
         vector newCellNormal = normalised(interfaceNormal_[i]);
 
         scalar normalRes = (1 - (cellNormal & newCellNormal));
-        avgAngle.insert(celli, avgDiffNormal);
-        normalResidual.insert(celli, normalRes);
+        normalResidual[i].celli = celli;
+        normalResidual[i].normalResidual = normalRes;
+        normalResidual[i].avgAngle = avgDiffNormal;
     }
 }
 
@@ -399,14 +398,14 @@ void Foam::reconstruction::plicRDF::reconstruct(bool forceUpdate)
     // nextToInterface is update on setInitNormals
     const boolList& nextToInterface_ = RDF_.nextToInterface();
 
-    labelHashSet tooCoarse;
+    bitSet tooCoarse(mesh_.nCells(),false);
 
     for (int iter=0; iter<iteration_; ++iter)
     {
         forAll(interfaceLabels_, i)
         {
             const label celli = interfaceLabels_[i];
-            if (mag(interfaceNormal_[i]) == 0 || tooCoarse.found(celli))
+            if (mag(interfaceNormal_[i]) == 0 || tooCoarse.test(celli))
             {
                 continue;
             }
@@ -439,8 +438,7 @@ void Foam::reconstruction::plicRDF::reconstruct(bool forceUpdate)
 
         normal_.correctBoundaryConditions();
         centre_.correctBoundaryConditions();
-        Map<scalar> residual;
-        Map<scalar> avgAngle;
+        List<normalRes> normalResidual(interfaceLabels_.size());
 
         surfaceVectorField::Boundary nHatb(mesh_.Sf().boundaryField());
         nHatb *= 1/(mesh_.magSf().boundaryField());
@@ -456,43 +454,42 @@ void Foam::reconstruction::plicRDF::reconstruct(bool forceUpdate)
             );
             RDF_.updateContactAngle(alpha1_, U_, nHatb);
             gradSurf(RDF_);
-            calcResidual(residual, avgAngle);
+            calcResidual(normalResidual);
         }
 
-
         label resCounter = 0;
         scalar avgRes = 0;
         scalar avgNormRes = 0;
 
-        Map<scalar>::iterator resIter = residual.begin();
-        Map<scalar>::iterator avgAngleIter = avgAngle.begin();
-
-        while (resIter.found())
+        forAll(normalResidual,i)
         {
-            if (avgAngleIter() > 0.26 && iter > 0) // 15 deg
+
+            const label celli = normalResidual[i].celli;
+            const scalar normalRes= normalResidual[i].normalResidual;
+            const scalar avgA = normalResidual[i].avgAngle;
+
+            if (avgA > 0.26 && iter > 0) // 15 deg
             {
-                tooCoarse.set(resIter.key());
+                tooCoarse.set(celli);
             }
             else
             {
-                avgRes += resIter();
+                avgRes += normalRes;
                 scalar normRes = 0;
-                scalar discreteError = 0.01*sqr(avgAngleIter());
+                scalar discreteError = 0.01*sqr(avgA);
                 if (discreteError != 0)
                 {
-                    normRes= resIter()/max(discreteError, tol_);
+                    normRes= normalRes/max(discreteError, tol_);
                 }
                 else
                 {
-                    normRes= resIter()/tol_;
+                    normRes= normalRes/tol_;
                 }
                 avgNormRes += normRes;
                 resCounter++;
 
             }
 
-            ++resIter;
-            ++avgAngleIter;
         }
 
         reduce(avgRes,sumOp<scalar>());

src/transportModels/geometricVoF/reconstructionSchemes/plicSchemes/plicRDF/plicRDF.H

@@ -78,6 +78,14 @@ class plicRDF
 {
     // Private Data
 
+        //-  residuals storage
+        struct normalRes
+        {
+            label celli;
+            scalar normalResidual;
+            scalar avgAngle;
+        };
+
         //- Reference to mesh
         const fvMesh& mesh_;
 
@@ -129,8 +137,7 @@ class plicRDF
         //- compute the normal residuals
         void calcResidual
         (
-            Map<scalar>& normalResidual,
-            Map<scalar>& avgAngle
+            List<normalRes>& normalResidual
         );
 
         //- interpolation of the normals from the previous time step

tutorials/DNS/dnsFoam/boxTurb16/0/U → tutorials/DNS/dnsFoam/boxTurb16/0.orig/U

@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volVectorField;
-    location    "0";
     object      U;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

tutorials/DNS/dnsFoam/boxTurb16/0/p → tutorials/DNS/dnsFoam/boxTurb16/0.orig/p

@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    location    "1";
     object      p;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -21,27 +20,7 @@ internalField   uniform 0;
 
 boundaryField
 {
-    patch0_half0
-    {
-        type            cyclic;
-    }
-    patch1_half0
-    {
-        type            cyclic;
-    }
-    patch2_half0
-    {
-        type            cyclic;
-    }
-    patch2_half1
-    {
-        type            cyclic;
-    }
-    patch1_half1
-    {
-        type            cyclic;
-    }
-    patch0_half1
+    ".*"
     {
         type            cyclic;
     }

tutorials/DNS/dnsFoam/boxTurb16/Allclean

@@ -3,7 +3,6 @@ cd "${0%/*}" || exit                                # Run from this directory
 . ${WM_PROJECT_DIR:?}/bin/tools/CleanFunctions      # Tutorial clean functions
 #------------------------------------------------------------------------------
 
-cleanCase
-rm -f 0/enstrophy
+cleanCase0
 
 #------------------------------------------------------------------------------

tutorials/DNS/dnsFoam/boxTurb16/Allrun

@@ -3,9 +3,14 @@ cd "${0%/*}" || exit                                # Run from this directory
 . ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
 #------------------------------------------------------------------------------
 
+restore0Dir
+
 runApplication blockMesh
+
 runApplication boxTurb
+
 runApplication $(getApplication)
-runApplication -s enstrophy  postProcess -func enstrophy
+
+runApplication postProcess -func enstrophy
 
 #------------------------------------------------------------------------------

tutorials/DNS/dnsFoam/boxTurb16/constant/boxTurbDict

@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "constant";
     object      boxTurbDict;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

tutorials/DNS/dnsFoam/boxTurb16/constant/transportProperties

@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "constant";
     object      transportProperties;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

tutorials/DNS/dnsFoam/boxTurb16/constant/turbulenceProperties

@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "constant";
     object      turbulenceProperties;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

tutorials/DNS/dnsFoam/boxTurb16/system/controlDict

@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "system";
     object      controlDict;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

tutorials/DNS/dnsFoam/boxTurb16/system/fvSchemes

@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "system";
     object      fvSchemes;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

tutorials/DNS/dnsFoam/boxTurb16/system/fvSolution

@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "system";
     object      fvSolution;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

tutorials/electromagnetics/electrostaticFoam/chargedWire/0/phi

@@ -53,4 +53,5 @@ boundaryField
     }
 }
 
+
 // ************************************************************************* //

tutorials/electromagnetics/electrostaticFoam/chargedWire/0/rho

@@ -52,4 +52,5 @@ boundaryField
     }
 }
 
+
 // ************************************************************************* //