WIP: Updates to field average windowing methods

src/functionObjects/field/fieldAverage/fieldAverage.C

@@ -207,7 +207,12 @@ void Foam::functionObjects::fieldAverage::writeAveragingProperties()
 
 void Foam::functionObjects::fieldAverage::readAveragingProperties()
 {
-    if (restartOnRestart_ || restartOnOutput_)
+    if
+    (
+        time_.timeIndex() == time_.startTimeIndex()
+     || restartOnRestart_
+     || restartOnOutput_
+    )
     {
         Info<< "    Starting averaging at time "
             << obr().time().timeOutputValue()
@@ -280,12 +285,6 @@ Foam::functionObjects::fieldAverage::fieldAverage
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::functionObjects::fieldAverage::~fieldAverage()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 bool Foam::functionObjects::fieldAverage::read(const dictionary& dict)

src/functionObjects/field/fieldAverage/fieldAverage.H

@@ -299,7 +299,7 @@ public:
 
 
     //- Destructor
-    virtual ~fieldAverage();
+    virtual ~fieldAverage() = default;
 
 
     // Member Functions

src/functionObjects/field/fieldAverage/fieldAverageItem/fieldAverageItem.C

@@ -123,15 +123,20 @@ void Foam::functionObjects::fieldAverageItem::addToWindow
     const scalar deltaT
 )
 {
-    windowTimes_.push(deltaT);
+//    windowTimes_.push(deltaT);
+    windowTimes_.push(0);
     windowFieldNames_.push(fieldName);
 }
 
 
 void Foam::functionObjects::fieldAverageItem::evolve(const objectRegistry& obr)
 {
+    const scalar dt = obr.time().deltaTValue();
+
     totalIter_++;
-    totalTime_ += obr.time().deltaTValue();
+    totalTime_ += dt;
+
+    // Advance time previous fields have existed
     forAllIters(windowTimes_, timeIter)
     {
         timeIter() += obr.time().deltaTValue();
@@ -142,7 +147,7 @@ void Foam::functionObjects::fieldAverageItem::evolve(const objectRegistry& obr)
 
     while (removeItem && windowTimes_.size())
     {
-        removeItem = !(inWindow(windowTimes_.first()));
+        removeItem = !(inWindow(dt, windowTimes_.first()));
 
         if (removeItem)
         {

src/functionObjects/field/fieldAverage/fieldAverageItem/fieldAverageItem.H

@@ -255,7 +255,7 @@ public:
             inline word windowFieldName(const word& prefix) const;
 
             //- Return true if time is inside window (including boundaries)
-            inline bool inWindow(const scalar t) const;
+            inline bool inWindow(const scalar dt, const scalar t) const;
 
             //- Return true if we wish to store window fields
             inline bool storeWindowFields() const;

src/functionObjects/field/fieldAverage/fieldAverageItem/fieldAverageItemI.H

@@ -190,17 +190,21 @@ Foam::word Foam::functionObjects::fieldAverageItem::windowFieldName
 }
 
 
-bool Foam::functionObjects::fieldAverageItem::inWindow(const scalar t) const
+bool Foam::functionObjects::fieldAverageItem::inWindow
+(
+    const scalar dt,
+    const scalar t
+) const
 {
     switch (base_)
     {
         case baseType::ITER:
         {
-            return round(t) <= round(window_) + 1;
+            return windowTimes_.size() <= round(window_); //  + 1;
         }
         case baseType::TIME:
         {
-            return t <= window_;
+            return t - dt <= window_;
         }
         default:
         {

src/functionObjects/field/fieldAverage/fieldAverageItem/fieldAverageItemTemplates.C

@@ -85,8 +85,6 @@ bool Foam::functionObjects::fieldAverageItem::calculateMeanField
                     //       needs to do 1 field lookup
 
                     label n = windowTimes_.size();
-                    const Type& lastField =
-                        obr.lookupObject<Type>(windowFieldNames_.first());
 
                     if (n <= round(window_))
                     {
@@ -95,6 +93,8 @@ bool Foam::functionObjects::fieldAverageItem::calculateMeanField
                     }
                     else
                     {
+                        const Type& lastField =
+                            obr.lookupObject<Type>(windowFieldNames_.first());
                         meanField += (baseField - lastField)/scalar(n - 1);
                     }
 
@@ -102,40 +102,56 @@ bool Foam::functionObjects::fieldAverageItem::calculateMeanField
                 }
                 case baseType::TIME:
                 {
-                    // Assuming non-uniform time step
-                    // Note: looks up all window fields from the registry
-
-                    meanField = 0*baseField;
-                    FIFOStack<scalar>::const_iterator timeIter =
-                        windowTimes_.begin();
-                    FIFOStack<word>::const_iterator nameIter =
-                        windowFieldNames_.begin();
-
-                    const Type* wOld = nullptr;
-
-                    for
-                    (
-                        ;
-                        timeIter != windowTimes_.end();
-                        ++timeIter, ++nameIter
-                    )
+                    if (windowTimes_.size() < 2)
                     {
-                        const word& fieldName = nameIter();
-                        const scalar dt = timeIter();
-                        const Type* w = obr.lookupObjectPtr<Type>(fieldName);
-
-                        meanField += dt*(*w);
-
-                        if (wOld)
+                        meanField = 1*baseField;
+                    }
+                    else
+                    {
+                        meanField = 0*baseField;
+                        FIFOStack<scalar>::const_iterator timeIter =
+                            windowTimes_.begin();
+                        FIFOStack<word>::const_iterator nameIter =
+                            windowFieldNames_.begin();
+
+                        const Type* w0Ptr =
+                            obr.lookupObjectPtr<Type>(nameIter());
+
+                        ++nameIter;
+                        scalar t0 = timeIter();
+                        ++timeIter;
+
+                        for
+                        (
+                            ;
+                            timeIter != windowTimes_.end();
+                            ++timeIter, ++nameIter
+                        )
                         {
-                            meanField -= dt*(*wOld);
+                            const Type* wPtr =
+                                obr.lookupObjectPtr<Type>(nameIter());
+                            const scalar t = timeIter();
+                            const Type& w = *wPtr;
+                            const Type& w0 = *w0Ptr;
+
+                            if (t0 > window_)
+                            {
+                                scalar tStar = max(0, window_ - t);
+                                meanField +=
+                                    ((w0 - w)/(t0 - t)*tStar + 2*w)*tStar;
+                            }
+                            else
+                            {
+                                meanField += (w + w0)*(t0 - t);
+                            }
+
+                            w0Ptr = wPtr;
+                            t0 = t;
                         }
 
-                        wOld = w;
+                        meanField /= 2*min(windowTimes_.first(), window_);
                     }
 
-                    meanField /= windowTimes_.first();
-
                     break;
                 }
                 default:
@@ -221,58 +237,107 @@ bool Foam::functionObjects::fieldAverageItem::calculatePrime2MeanField
         }
         case windowType::EXACT:
         {
-            // Not storing old time mean fields - treat all as TIME (integrated)
-            prime2MeanField = 0*prime2MeanField;
             FIFOStack<scalar>::const_iterator timeIter =
                 windowTimes_.begin();
             FIFOStack<word>::const_iterator nameIter =
                 windowFieldNames_.begin();
 
-            switch (base_)
+            if (windowFieldNames_.size() < 2)
             {
-                case baseType::ITER:
-                {
-                    // ITER method stores an additional entry compared to TIME
-                    ++timeIter;
-                    ++nameIter;
-
-                    if (timeIter == windowTimes_.end()) return false;
-
-                    break;
-                }
-                default:
-                {}
+                prime2MeanField = sqr(baseField - meanField);
             }
+            else
+            {
+                switch (base_)
+                {
+                    case baseType::ITER:
+                    {
+                        prime2MeanField = 0*prime2MeanField;
+                        label n = min(windowFieldNames_.size(), window_);
+                        DebugVar(windowFieldNames_.size());
+                        DebugVar(window_);
+                        if (windowFieldNames_.size() > window_)
+                        {
+                            ++nameIter;
+                            ++timeIter;
+                        }
 
+                        for
+                        (
+                            ;
+                            nameIter != windowFieldNames_.end();
+                            ++nameIter
+                        )
+                        {
+                            const Type1* wPtr =
+                                obr.lookupObjectPtr<Type1>(nameIter());
+                            const Type1& w = *wPtr;
 
-            scalar windowLength = timeIter();
+                            prime2MeanField += sqr(w - meanField);
+                        }
 
-            const Type1* wOld = nullptr;
+                        prime2MeanField /= n;
 
-            for
-            (
-                ;
-                timeIter != windowTimes_.end();
-                ++timeIter, ++nameIter
-            )
-            {
-                const word& fieldName = nameIter();
-                const scalar dt = timeIter();
-                const Type1* w = obr.lookupObjectPtr<Type1>(fieldName);
+                        break;
+                    }
+                    case baseType::TIME:
+                    {
+                        prime2MeanField = 0*prime2MeanField;
+                        const Type1* w0Ptr =
+                            obr.lookupObjectPtr<Type1>(nameIter());
+
+                        ++nameIter;
+                        scalar t0 = timeIter();
+                        ++timeIter;
+                        for
+                        (
+                            ;
+                            timeIter != windowTimes_.end();
+                            ++timeIter, ++nameIter
+                        )
+                        {
+                            const Type1* w1Ptr =
+                                obr.lookupObjectPtr<Type1>(nameIter());
+                            const scalar t1 = timeIter();
+
+                            const Type1& w0 = *w0Ptr;
+                            const Type1& w1 = *w1Ptr;
+
+                            if (t0 > window_)
+                            {
+                                scalar tStar = max(0, window_ - t1);
+                                prime2MeanField +=
+                                    tStar
+                                   *sqr
+                                    (
+                                        0.5*((w1 - w0)/(t1 - t0)*tStar + 2*w1)
+                                      - meanField
+                                    );
+                            }
+                            else
+                            {
+                                prime2MeanField +=
+                                    (t0 - t1)*(sqr(0.5*(w0 + w1) - meanField));
+                            }
+
+                            w0Ptr = w1Ptr;
+                            t0 = t1;
+                        }
 
-                prime2MeanField += dt*(sqr((*w) - meanField));
+                        prime2MeanField /= min(windowTimes_.first(), window_);
 
-                if (wOld)
-                {
-                    prime2MeanField -= dt*(sqr((*wOld) - meanField));
+                        break;
+                    }
+                    default:
+                    {
+                        FatalErrorInFunction
+                            << "Unhandled baseType enumeration "
+                            << baseTypeNames_[base_]
+                            << abort(FatalError);
+                    }
                 }
-
-                wOld = w;
             }
 
-            prime2MeanField /= windowLength;
-
-
             break;
         }
         default: