From 3e8b97fef603d3ab1469097f8216cc58678d0d0e Mon Sep 17 00:00:00 2001
From: Will Bainbridge <http://cfd.direct>
Date: Wed, 30 Nov 2022 12:31:33 +0000
Subject: [PATCH] functionObjects: fieldAverage: Corrected initialisation
 behaviour

The fieldAverage can now average fields that do not exist at
construction time, and it also supports restart on cases in which
the mesh topology is changing.
---
 .../fields/volFields/volFieldsFwd.H           |   3 +
 .../field/fieldAverage/fieldAverage.C         | 343 ++++++++++------
 .../field/fieldAverage/fieldAverage.H         |  60 +--
 .../fieldAverageItem/fieldAverageItem.H       |   8 +-
 .../fieldAverage/fieldAverageTemplates.C      | 375 ++++++++++--------
 5 files changed, 469 insertions(+), 320 deletions(-)

diff --git a/src/finiteVolume/fields/volFields/volFieldsFwd.H b/src/finiteVolume/fields/volFields/volFieldsFwd.H
index 664f4a1eeb..51a57fb3f2 100644
--- a/src/finiteVolume/fields/volFields/volFieldsFwd.H
+++ b/src/finiteVolume/fields/volFields/volFieldsFwd.H
@@ -54,6 +54,9 @@ class GeometricField;
 template<class Type>
 using VolField = GeometricField<Type, fvPatchField, volMesh>;
 
+template<class Type>
+using VolInternalField = typename VolField<Type>::Internal;
+
 typedef GeometricField<label, fvPatchField, volMesh> volLabelField;
 typedef GeometricField<scalar, fvPatchField, volMesh> volScalarField;
 typedef GeometricField<vector, fvPatchField, volMesh> volVectorField;
diff --git a/src/functionObjects/field/fieldAverage/fieldAverage.C b/src/functionObjects/field/fieldAverage/fieldAverage.C
index b7d77c9d92..2e49b21a46 100644
--- a/src/functionObjects/field/fieldAverage/fieldAverage.C
+++ b/src/functionObjects/field/fieldAverage/fieldAverage.C
@@ -56,8 +56,43 @@ const Foam::NamedEnum
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
-void Foam::functionObjects::fieldAverage::resetFields()
+void Foam::functionObjects::fieldAverage::initialise()
 {
+    // Initialise any unset times
+    forAll(totalTime_, fieldi)
+    {
+        if (totalTime_[fieldi] < 0)
+        {
+            totalTime_[fieldi] = obr_.time().deltaTValue();
+        }
+    }
+
+    // Initialise mean fields
+    forAll(faItems_, fieldi)
+    {
+        initialiseMeanField<scalar>(fieldi);
+        initialiseMeanField<vector>(fieldi);
+        initialiseMeanField<sphericalTensor>(fieldi);
+        initialiseMeanField<symmTensor>(fieldi);
+        initialiseMeanField<tensor>(fieldi);
+    }
+    forAll(faItems_, fieldi)
+    {
+        initialisePrime2MeanField<scalar, scalar>(fieldi);
+        initialisePrime2MeanField<vector, symmTensor>(fieldi);
+    }
+}
+
+
+void Foam::functionObjects::fieldAverage::restart()
+{
+    Log << "    Restarting averaging at time " << obr_.time().name() << nl;
+
+    // Clear the times
+    totalIter_ = 1;
+    totalTime_ = -1;
+
+    // Clear mean fields
     forAll(faItems_, i)
     {
         if (faItems_[i].mean())
@@ -76,70 +111,19 @@ void Foam::functionObjects::fieldAverage::resetFields()
             }
         }
     }
-}
-
-
-void Foam::functionObjects::fieldAverage::initialise()
-{
-    if (!totalIter_.size())
-    {
-        totalIter_.setSize(faItems_.size(), 1);
-    }
-
-    if (!totalTime_.size())
-    {
-        totalTime_.setSize(faItems_.size(), obr_.time().deltaTValue());
-    }
-    else
-    {
-        // Check if totalTime_ has been set otherwise initialise
-        forAll(totalTime_, fieldi)
-        {
-            if (totalTime_[fieldi] < 0)
-            {
-                totalTime_[fieldi] = obr_.time().deltaTValue();
-            }
-        }
-    }
-
-    resetFields();
-
-    // Add mean fields to the field lists
-    forAll(faItems_, fieldi)
-    {
-        addMeanField<scalar>(fieldi);
-        addMeanField<vector>(fieldi);
-        addMeanField<sphericalTensor>(fieldi);
-        addMeanField<symmTensor>(fieldi);
-        addMeanField<tensor>(fieldi);
-    }
-
-    // Add prime-squared mean fields to the field lists
-    forAll(faItems_, fieldi)
-    {
-        addPrime2MeanField<scalar, scalar>(fieldi);
-        addPrime2MeanField<vector, symmTensor>(fieldi);
-    }
-
-    // ensure first averaging works unconditionally
-    prevTimeIndex_ = -1;
-}
-
-
-void Foam::functionObjects::fieldAverage::restart()
-{
-    Log << "    Restarting averaging at time " << obr_.time().name()
-        << nl << endl;
-
-    totalIter_.clear();
-    totalTime_.clear();
 
+    // Re-create any mean fields
     initialise();
+
+    // Ensure first averaging works unconditionally
+    prevTimeIndex_ = -1;
 }
 
 
 void Foam::functionObjects::fieldAverage::calcAverages()
 {
+    Log << type() << " " << name() << ":" << nl;
+
     const label currentTimeIndex = obr_.time().timeIndex();
     const scalar currentTime = obr_.time().value();
 
@@ -147,19 +131,20 @@ void Foam::functionObjects::fieldAverage::calcAverages()
     {
         return;
     }
-    else
-    {
-        prevTimeIndex_ = currentTimeIndex;
-    }
+
+    prevTimeIndex_ = currentTimeIndex;
 
     if (periodicRestart_ && currentTime > restartPeriod_*periodIndex_)
     {
         restart();
         periodIndex_++;
     }
+    else
+    {
+        initialise();
+    }
 
-    Log << type() << " " << name() << nl
-        << "    Calculating averages" << nl;
+    Log << "    Calculating averages" << nl;
 
     addMeanSqrToPrime2Mean<scalar, scalar>();
     addMeanSqrToPrime2Mean<vector, symmTensor>();
@@ -185,7 +170,8 @@ void Foam::functionObjects::fieldAverage::calcAverages()
 
 void Foam::functionObjects::fieldAverage::writeAverages() const
 {
-    Log << "    Writing average fields" << endl;
+    Log << type() << " " << name() << ":" << nl
+        << "    Writing average fields" << endl;
 
     writeFields<scalar>();
     writeFields<vector>();
@@ -193,12 +179,6 @@ void Foam::functionObjects::fieldAverage::writeAverages() const
     writeFields<symmTensor>();
     writeFields<tensor>();
 
-    Log << endl;
-}
-
-
-void Foam::functionObjects::fieldAverage::writeAveragingProperties() const
-{
     timeIOdictionary propsDict
     (
         IOobject
@@ -227,19 +207,49 @@ void Foam::functionObjects::fieldAverage::writeAveragingProperties() const
 }
 
 
-void Foam::functionObjects::fieldAverage::readAveragingProperties()
+void Foam::functionObjects::fieldAverage::read
+(
+    const dictionary& dict,
+    const bool construct
+)
 {
-    if ((restartOnRestart_ || restartOnOutput_) && log)
+    dict.readIfPresent("restartOnRestart", restartOnRestart_);
+    dict.readIfPresent("restartOnOutput", restartOnOutput_);
+    dict.readIfPresent("periodicRestart", periodicRestart_);
+
+    if (periodicRestart_)
     {
-        Info<< "    Starting averaging at time " << obr_.time().name()
-            << nl;
+        dict.lookup("restartPeriod") >> restartPeriod_;
     }
-    else
+
+    mean_ = dict.lookupOrDefault<Switch>("mean", true);
+    prime2Mean_ = dict.lookupOrDefault<Switch>("prime2Mean", false);
+    base_ = baseTypeNames_[dict.lookupOrDefault<word>("base", "time")];
+    window_ = dict.lookupOrDefault<scalar>("window", -1);
+    windowName_ = dict.lookupOrDefault<word>("windowName", "");
+
+    if (construct)
     {
+        // First read of a run. Look for properties dict and read total
+        // iter/time values from it if available.
+
+        faItems_.clear();
+        totalIter_.clear();
+        totalTime_.clear();
+
+        faItems_ =
+            PtrList<fieldAverageItem>
+            (
+                dict.lookup("fields"),
+                fieldAverageItem::iNew(*this)
+            );
+        totalIter_.setSize(faItems_.size(), 1);
+        totalTime_.setSize(faItems_.size(), -1);
+
         typeIOobject<timeIOdictionary> propsDictHeader
         (
             name() + "Properties",
-            obr_.time().startTime().name(),
+            obr_.time().name(),
             "uniform",
             obr_,
             IOobject::MUST_READ_IF_MODIFIED,
@@ -247,40 +257,120 @@ void Foam::functionObjects::fieldAverage::readAveragingProperties()
             false
         );
 
-        if (!propsDictHeader.headerOk())
+        dictionary propsDict;
+        if
+        (
+            !restartOnRestart_
+         && !restartOnOutput_
+         && propsDictHeader.headerOk()
+        )
         {
-            Log << "    Starting averaging at time "
-                << obr_.time().name() << nl;
-
-            return;
+            propsDict = timeIOdictionary(propsDictHeader);
         }
 
-        timeIOdictionary propsDict(propsDictHeader);
-
-        Log << "    Restarting averaging for fields:" << nl;
-
-        totalIter_.setSize(faItems_.size(), 1);
-
-        // Initialise totalTime with negative values
-        // to indicate that it has not been set
-        totalTime_.setSize(faItems_.size(), -1);
+        bool first = true;
+        forAll(faItems_, fieldi)
+        {
+            const word& fieldName = faItems_[fieldi].fieldName();
+            if (!propsDict.found(fieldName))
+            {
+                if (first)
+                {
+                    Log << "    Starting averaging for fields:" << nl;
+                    first = false;
+                }
+                Log << "        " << fieldName << nl;
+            }
+        }
 
+        first = true;
         forAll(faItems_, fieldi)
         {
             const word& fieldName = faItems_[fieldi].fieldName();
             if (propsDict.found(fieldName))
             {
-                dictionary fieldDict(propsDict.subDict(fieldName));
-
-                totalIter_[fieldi] = fieldDict.lookup<label>("totalIter");
-                totalTime_[fieldi] = fieldDict.lookup<scalar>("totalTime");
-
+                if (first)
+                {
+                    Log << "    Restarting averaging for fields:" << nl;
+                    first = false;
+                }
+                const dictionary& fieldPropsDict = propsDict.subDict(fieldName);
+                totalIter_[fieldi] = fieldPropsDict.lookup<label>("totalIter");
+                totalTime_[fieldi] = fieldPropsDict.lookup<scalar>("totalTime");
                 Log << "        " << fieldName
                     << " iters = " << totalIter_[fieldi]
                     << " time = " << totalTime_[fieldi] << nl;
             }
         }
     }
+    else
+    {
+        // Re-read during a run. Read the items and copy the per-field total
+        // iter/times from the old items to the new.
+
+        PtrList<fieldAverageItem> faItems0;
+        List<label> totalIter0;
+        List<scalar> totalTime0;
+        faItems0.transfer(faItems_);
+        totalIter0.transfer(totalIter_);
+        totalTime0.transfer(totalTime_);
+
+        faItems_ =
+            PtrList<fieldAverageItem>
+            (
+                dict.lookup("fields"),
+                fieldAverageItem::iNew(*this)
+            );
+        totalIter_.resize(faItems_.size(), 1);
+        totalTime_.resize(faItems_.size(), -1);
+
+        // Map from field to old-field index
+        labelList fieldiFieldi0s(faItems_.size(), -1);
+        forAll(faItems_, fieldi)
+        {
+            const word& fieldName = faItems_[fieldi].fieldName();
+            forAll(faItems0, fieldi0)
+            {
+                if (faItems0[fieldi0].fieldName() == fieldName)
+                {
+                    fieldiFieldi0s[fieldi] = fieldi0;
+                    break;
+                }
+            }
+        }
+
+        bool first = true;
+        forAll(faItems_, fieldi)
+        {
+            if (fieldiFieldi0s[fieldi] == -1)
+            {
+                if (first)
+                {
+                    Log << "    Starting averaging for fields:" << nl;
+                    first = true;
+                }
+                Log << "        " << faItems_[fieldi].fieldName() << nl;
+            }
+        }
+
+        first = true;
+        forAll(faItems_, fieldi)
+        {
+            if (fieldiFieldi0s[fieldi] != -1)
+            {
+                if (first)
+                {
+                    Log << "    Continuing averaging for fields:" << nl;
+                    first = true;
+                }
+                totalIter_[fieldi] = totalIter0[fieldiFieldi0s[fieldi]];
+                totalTime_[fieldi] = totalTime0[fieldiFieldi0s[fieldi]];
+                Log << "        " << faItems_[fieldi].fieldName()
+                    << " iters = " << totalIter_[fieldi]
+                    << " time = " << totalTime_[fieldi] << nl;
+            }
+        }
+    }
 }
 
 
@@ -307,7 +397,28 @@ Foam::functionObjects::fieldAverage::fieldAverage
     totalTime_(),
     periodIndex_(1)
 {
-    read(dict);
+    fvMeshFunctionObject::read(dict);
+
+    Log << type() << " " << name << ":" << nl;
+
+    read(dict, true);
+
+    // Read any available mean fields
+    forAll(faItems_, fieldi)
+    {
+        readMeanField<scalar>(fieldi);
+        readMeanField<vector>(fieldi);
+        readMeanField<sphericalTensor>(fieldi);
+        readMeanField<symmTensor>(fieldi);
+        readMeanField<tensor>(fieldi);
+    }
+    forAll(faItems_, fieldi)
+    {
+        readPrime2MeanField<scalar, scalar>(fieldi);
+        readPrime2MeanField<vector, symmTensor>(fieldi);
+    }
+
+    Log << endl;
 }
 
 
@@ -323,36 +434,9 @@ bool Foam::functionObjects::fieldAverage::read(const dictionary& dict)
 {
     fvMeshFunctionObject::read(dict);
 
-
     Log << type() << " " << name() << ":" << nl;
 
-    dict.readIfPresent("restartOnRestart", restartOnRestart_);
-    dict.readIfPresent("restartOnOutput", restartOnOutput_);
-    dict.readIfPresent("periodicRestart", periodicRestart_);
-
-    mean_ = dict.lookupOrDefault<Switch>("mean", true);
-    prime2Mean_ = dict.lookupOrDefault<Switch>("prime2Mean", false);
-    base_ = baseTypeNames_
-    [
-        dict.lookupOrDefault<word>("base", "time")
-    ];
-    window_ = dict.lookupOrDefault<scalar>("window", -1);
-    windowName_ = dict.lookupOrDefault<word>("windowName", "");
-
-    faItems_ = PtrList<fieldAverageItem>
-    (
-        dict.lookup("fields"),
-        fieldAverageItem::iNew(*this)
-    );
-
-    if (periodicRestart_)
-    {
-        dict.lookup("restartPeriod") >> restartPeriod_;
-    }
-
-    readAveragingProperties();
-
-    initialise();
+    read(dict, false);
 
     Log << endl;
 
@@ -383,19 +467,16 @@ bool Foam::functionObjects::fieldAverage::execute()
 
         return false;
     }
-    else
-    {
-        calcAverages();
 
-        return true;
-    }
+    calcAverages();
+
+    return true;
 }
 
 
 bool Foam::functionObjects::fieldAverage::write()
 {
     writeAverages();
-    writeAveragingProperties();
 
     if (restartOnOutput_)
     {
diff --git a/src/functionObjects/field/fieldAverage/fieldAverage.H b/src/functionObjects/field/fieldAverage/fieldAverage.H
index be695a25e9..ce6e9d2c33 100644
--- a/src/functionObjects/field/fieldAverage/fieldAverage.H
+++ b/src/functionObjects/field/fieldAverage/fieldAverage.H
@@ -127,7 +127,6 @@ class fieldAverage
 :
     public fvMeshFunctionObject
 {
-
 public:
 
         //- Enumeration defining the averaging base type
@@ -179,26 +178,20 @@ protected:
         //  calculated and output
         PtrList<fieldAverageItem> faItems_;
 
-        // Counters
+        //- Iteration steps counter
+        List<label> totalIter_;
 
-            //- Iteration steps counter
-            List<label> totalIter_;
+        //- Total time counter
+        List<scalar> totalTime_;
 
-            //- Total time counter
-            List<scalar> totalTime_;
-
-            //- Index for periodic restart
-            label periodIndex_;
+        //- Index for periodic restart
+        label periodIndex_;
 
 
     // Protected Member Functions
 
         // Initialisation routines
 
-            //- Checkout fields (causes deletion) from the database
-            //  and reset lists
-            void resetFields();
-
             //- Reset lists (clear existing values) and initialise averaging.
             //  Check requested field averages are valid, populate field lists
             void initialise();
@@ -206,21 +199,37 @@ protected:
             //- Restart averaging for restartOnOutput
             void restart();
 
-            //- Add mean average field to database
+            //- Read the mean average field
             template<class Type>
-            void addMeanFieldType(const label fieldi);
+            void readMeanFieldType(const label fieldi);
 
-            //- Add mean average field to database
+            //- Read the mean average field
             template<class Type>
-            void addMeanField(const label fieldi);
+            void readMeanField(const label fieldi);
 
-            //- Add prime-squared average field to database
-            template<class Type1, class Type2>
-            void addPrime2MeanFieldType(const label fieldi);
+            //- Initialise the mean average field
+            template<class Type>
+            void initialiseMeanFieldType(const label fieldi);
 
-            //- Add prime-squared average field to database
+            //- Initialise the mean average field
+            template<class Type>
+            void initialiseMeanField(const label fieldi);
+
+            //- Read the prime-squared average field
             template<class Type1, class Type2>
-            void addPrime2MeanField(const label fieldi);
+            void readPrime2MeanFieldType(const label fieldi);
+
+            //- Read the prime-squared average field
+            template<class Type1, class Type2>
+            void readPrime2MeanField(const label fieldi);
+
+            //- Initialise the prime-squared average field
+            template<class Type1, class Type2>
+            void initialisePrime2MeanFieldType(const label fieldi);
+
+            //- Initialise the prime-squared average field
+            template<class Type1, class Type2>
+            void initialisePrime2MeanField(const label fieldi);
 
 
         // Calculation functions
@@ -266,11 +275,8 @@ protected:
             template<class Type>
             void writeFields() const;
 
-            //- Write averaging properties - steps and time
-            void writeAveragingProperties() const;
-
-            //- Read averaging properties - steps and time
-            void readAveragingProperties();
+            //- Read
+            void read(const dictionary& dict, const bool construct);
 
 
 public:
diff --git a/src/functionObjects/field/fieldAverage/fieldAverageItem/fieldAverageItem.H b/src/functionObjects/field/fieldAverage/fieldAverageItem/fieldAverageItem.H
index d7aea4552a..eb0ea98823 100644
--- a/src/functionObjects/field/fieldAverage/fieldAverageItem/fieldAverageItem.H
+++ b/src/functionObjects/field/fieldAverage/fieldAverageItem/fieldAverageItem.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -32,15 +32,9 @@ Description
     {
         mean            on;   // (default = on)
         prime2Mean      on;   // (default = off)
-        base            time; // time or iteration (default = time)
-        window          200;  // optional averaging window
-        windowName      w1;   // optional window name (default = "")
     }
     \endverbatim
 
-    The averaging window corresponds to the averaging interval (iters or time)
-    If not specified, the averaging is over 'all iters/time'
-
 SourceFiles
     fieldAverageItem.C
     fieldAverageItemIO.C
diff --git a/src/functionObjects/field/fieldAverage/fieldAverageTemplates.C b/src/functionObjects/field/fieldAverage/fieldAverageTemplates.C
index c6a8b25b2b..d94af60829 100644
--- a/src/functionObjects/field/fieldAverage/fieldAverageTemplates.C
+++ b/src/functionObjects/field/fieldAverage/fieldAverageTemplates.C
@@ -23,6 +23,7 @@ License
 
 \*---------------------------------------------------------------------------*/
 
+#include "fieldAverage.H"
 #include "fieldAverageItem.H"
 #include "volFields.H"
 #include "surfaceFields.H"
@@ -30,18 +31,71 @@ License
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 template<class Type>
-void Foam::functionObjects::fieldAverage::addMeanFieldType(const label fieldi)
+void Foam::functionObjects::fieldAverage::readMeanFieldType(const label fieldi)
+{
+    const word& meanFieldName = faItems_[fieldi].meanFieldName();
+
+    IOobject meanFieldIo
+    (
+        meanFieldName,
+        obr_.time().name(),
+        obr_,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    );
+
+    if
+    (
+        meanFieldIo.headerOk()
+     && meanFieldIo.headerClassName() == Type::typeName
+    )
+    {
+        if (obr_.found(meanFieldName))
+        {
+            Log << "    Cannot read average field " << meanFieldName
+                << " since an object with that name already exists."
+                << " Disabling averaging for field." << endl;
+
+            faItems_[fieldi].mean() = false;
+        }
+        else
+        {
+            Log << "    Reading field " << meanFieldName << endl;
+
+            obr_.store(new Type(meanFieldIo, mesh_));
+        }
+    }
+}
+
+
+template<class Type>
+void Foam::functionObjects::fieldAverage::readMeanField(const label fieldi)
+{
+    if (faItems_[fieldi].mean())
+    {
+        readMeanFieldType<VolField<Type>>(fieldi);
+        readMeanFieldType<VolInternalField<Type>>(fieldi);
+        readMeanFieldType<SurfaceField<Type>>(fieldi);
+    }
+}
+
+
+template<class Type>
+void Foam::functionObjects::fieldAverage::initialiseMeanFieldType
+(
+    const label fieldi
+)
 {
     const word& fieldName = faItems_[fieldi].fieldName();
     const word& meanFieldName = faItems_[fieldi].meanFieldName();
 
-    Log << "    Reading/initialising field " << meanFieldName << endl;
-
     if (obr_.foundObject<Type>(meanFieldName))
-    {}
+    {
+        // Do nothing ...
+    }
     else if (obr_.found(meanFieldName))
     {
-        Log << "    Cannot allocate average field " << meanFieldName
+        Log << "    Cannot initialise average field " << meanFieldName
             << " since an object with that name already exists."
             << " Disabling averaging for field." << endl;
 
@@ -49,9 +103,10 @@ void Foam::functionObjects::fieldAverage::addMeanFieldType(const label fieldi)
     }
     else
     {
+        Log << "    Initialising field " << meanFieldName << endl;
+
         const Type& baseField = obr_.lookupObject<Type>(fieldName);
 
-        // Store on registry
         obr_.store
         (
             new Type
@@ -59,12 +114,8 @@ void Foam::functionObjects::fieldAverage::addMeanFieldType(const label fieldi)
                 IOobject
                 (
                     meanFieldName,
-                    obr_.time().timeName(obr_.time().startTime().value()),
-                    obr_,
-                    restartOnOutput_
-                  ? IOobject::NO_READ
-                  : IOobject::READ_IF_PRESENT,
-                    IOobject::NO_WRITE
+                    obr_.time().name(),
+                    obr_
                 ),
                 1*baseField
             )
@@ -74,38 +125,100 @@ void Foam::functionObjects::fieldAverage::addMeanFieldType(const label fieldi)
 
 
 template<class Type>
-void Foam::functionObjects::fieldAverage::addMeanField(const label fieldi)
+void Foam::functionObjects::fieldAverage::initialiseMeanField
+(
+    const label fieldi
+)
 {
     if (faItems_[fieldi].mean())
     {
-        typedef GeometricField<Type, fvPatchField, volMesh>
-            VolFieldType;
-
-        typedef typename VolFieldType::Internal InternalType;
-
-        typedef GeometricField<Type, fvsPatchField, surfaceMesh>
-            SurfaceFieldType;
-
         const word& fieldName = faItems_[fieldi].fieldName();
 
-        if (obr_.foundObject<VolFieldType>(fieldName))
+        if (obr_.foundObject<VolField<Type>>(fieldName))
         {
-            addMeanFieldType<VolFieldType>(fieldi);
+            initialiseMeanFieldType<VolField<Type>>(fieldi);
         }
-        else if (obr_.foundObject<InternalType>(fieldName))
+        else if (obr_.foundObject<VolInternalField<Type>>(fieldName))
         {
-            addMeanFieldType<InternalType>(fieldi);
+            initialiseMeanFieldType<VolInternalField<Type>>(fieldi);
         }
-        else if (obr_.foundObject<SurfaceFieldType>(fieldName))
+        else if (obr_.foundObject<SurfaceField<Type>>(fieldName))
         {
-            addMeanFieldType<SurfaceFieldType>(fieldi);
+            initialiseMeanFieldType<SurfaceField<Type>>(fieldi);
         }
     }
 }
 
 
 template<class Type1, class Type2>
-void Foam::functionObjects::fieldAverage::addPrime2MeanFieldType
+void Foam::functionObjects::fieldAverage::readPrime2MeanFieldType
+(
+    const label fieldi
+)
+{
+    const word& prime2MeanFieldName = faItems_[fieldi].prime2MeanFieldName();
+
+    IOobject prime2MeanFieldIo
+    (
+        prime2MeanFieldName,
+        obr_.time().name(),
+        obr_,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    );
+
+    if
+    (
+        prime2MeanFieldIo.headerOk()
+     && prime2MeanFieldIo.headerClassName() == Type2::typeName
+    )
+    {
+        if (obr_.found(prime2MeanFieldName))
+        {
+            Log << "    Cannot read average field " << prime2MeanFieldName
+                << " since an object with that name already exists."
+                << " Disabling averaging for field." << endl;
+
+            faItems_[fieldi].mean() = false;
+        }
+        else
+        {
+            Log << "    Reading field " << prime2MeanFieldName << endl;
+
+            obr_.store(new Type2(prime2MeanFieldIo, mesh_));
+        }
+    }
+}
+
+
+template<class Type1, class Type2>
+void Foam::functionObjects::fieldAverage::readPrime2MeanField
+(
+    const label fieldi
+)
+{
+    if (faItems_[fieldi].prime2Mean())
+    {
+        if (!faItems_[fieldi].mean())
+        {
+            FatalErrorInFunction
+                << "To calculate the prime-squared average, the "
+                << "mean average must also be selected for field "
+                << faItems_[fieldi].fieldName() << nl << exit(FatalError);
+        }
+
+        readPrime2MeanFieldType
+        <VolField<Type1>, VolField<Type2>>(fieldi);
+        readPrime2MeanFieldType
+        <VolInternalField<Type1>, VolInternalField<Type2>>(fieldi);
+        readPrime2MeanFieldType
+        <SurfaceField<Type1>, SurfaceField<Type2>>(fieldi);
+    }
+}
+
+
+template<class Type1, class Type2>
+void Foam::functionObjects::fieldAverage::initialisePrime2MeanFieldType
 (
     const label fieldi
 )
@@ -114,13 +227,13 @@ void Foam::functionObjects::fieldAverage::addPrime2MeanFieldType
     const word& meanFieldName = faItems_[fieldi].meanFieldName();
     const word& prime2MeanFieldName = faItems_[fieldi].prime2MeanFieldName();
 
-    Log << "    Reading/initialising field " << prime2MeanFieldName << nl;
-
     if (obr_.foundObject<Type2>(prime2MeanFieldName))
-    {}
+    {
+        // Do nothing ...
+    }
     else if (obr_.found(prime2MeanFieldName))
     {
-        Log << "    Cannot allocate average field " << prime2MeanFieldName
+        Log << "    Cannot initialise average field " << prime2MeanFieldName
             << " since an object with that name already exists."
             << " Disabling averaging for field." << nl;
 
@@ -128,10 +241,11 @@ void Foam::functionObjects::fieldAverage::addPrime2MeanFieldType
     }
     else
     {
+        Log << "    Initialising field " << prime2MeanFieldName << nl;
+
         const Type1& baseField = obr_.lookupObject<Type1>(fieldName);
         const Type1& meanField = obr_.lookupObject<Type1>(meanFieldName);
 
-        // Store on registry
         obr_.store
         (
             new Type2
@@ -139,12 +253,8 @@ void Foam::functionObjects::fieldAverage::addPrime2MeanFieldType
                 IOobject
                 (
                     prime2MeanFieldName,
-                    obr_.time().timeName(obr_.time().startTime().value()),
-                    obr_,
-                    restartOnOutput_
-                  ? IOobject::NO_READ
-                  : IOobject::READ_IF_PRESENT,
-                    IOobject::NO_WRITE
+                    obr_.time().name(),
+                    obr_
                 ),
                 sqr(baseField) - sqr(meanField)
             )
@@ -154,16 +264,11 @@ void Foam::functionObjects::fieldAverage::addPrime2MeanFieldType
 
 
 template<class Type1, class Type2>
-void Foam::functionObjects::fieldAverage::addPrime2MeanField(const label fieldi)
+void Foam::functionObjects::fieldAverage::initialisePrime2MeanField
+(
+    const label fieldi
+)
 {
-    typedef GeometricField<Type1, fvPatchField, volMesh> VolFieldType1;
-    typedef typename VolFieldType1::Internal InternalType1;
-    typedef GeometricField<Type1, fvsPatchField, surfaceMesh> SurfaceFieldType1;
-
-    typedef GeometricField<Type2, fvPatchField, volMesh> VolFieldType2;
-    typedef typename VolFieldType2::Internal InternalType2;
-    typedef GeometricField<Type2, fvsPatchField, surfaceMesh> SurfaceFieldType2;
-
     if (faItems_[fieldi].prime2Mean())
     {
         const word& fieldName = faItems_[fieldi].fieldName();
@@ -176,20 +281,20 @@ void Foam::functionObjects::fieldAverage::addPrime2MeanField(const label fieldi)
                 << fieldName << nl << exit(FatalError);
         }
 
-        if (obr_.foundObject<VolFieldType1>(fieldName))
+        if (obr_.foundObject<VolField<Type1>>(fieldName))
         {
-            addPrime2MeanFieldType<VolFieldType1, VolFieldType2>(fieldi);
+            initialisePrime2MeanFieldType
+            <VolField<Type1>, VolField<Type2>>(fieldi);
         }
-        else if (obr_.foundObject<InternalType1>(fieldName))
+        else if (obr_.foundObject<VolInternalField<Type1>>(fieldName))
         {
-            addPrime2MeanFieldType<InternalType1, InternalType2>(fieldi);
+            initialisePrime2MeanFieldType
+            <VolInternalField<Type1>, VolInternalField<Type2>>(fieldi);
         }
-        else if (obr_.foundObject<SurfaceFieldType1>(fieldName))
+        else if (obr_.foundObject<SurfaceField<Type1>>(fieldName))
         {
-            addPrime2MeanFieldType<SurfaceFieldType1, SurfaceFieldType2>
-            (
-                fieldi
-            );
+            initialisePrime2MeanFieldType
+            <SurfaceField<Type1>, SurfaceField<Type2>>(fieldi);
         }
     }
 }
@@ -203,63 +308,56 @@ void Foam::functionObjects::fieldAverage::calculateMeanFieldType
 {
     const word& fieldName = faItems_[fieldi].fieldName();
 
-    if (obr_.foundObject<Type>(fieldName))
+    const Type& baseField = obr_.lookupObject<Type>(fieldName);
+
+    Type& meanField =
+        obr_.lookupObjectRef<Type>(faItems_[fieldi].meanFieldName());
+
+    scalar dt = obr_.time().deltaTValue();
+    scalar Dt = totalTime_[fieldi];
+
+    if (iterBase())
     {
-        const Type& baseField = obr_.lookupObject<Type>(fieldName);
-
-        Type& meanField =
-            obr_.lookupObjectRef<Type>(faItems_[fieldi].meanFieldName());
-
-        scalar dt = obr_.time().deltaTValue();
-        scalar Dt = totalTime_[fieldi];
-
-        if (iterBase())
-        {
-            dt = 1;
-            Dt = scalar(totalIter_[fieldi]);
-        }
-
-        scalar beta = dt/Dt;
-
-        if (window() > 0)
-        {
-            const scalar w = window();
-
-            if (Dt - dt >= w)
-            {
-                beta = dt/w;
-            }
-        }
-
-        meanField = (1 - beta)*meanField + beta*baseField;
+        dt = 1;
+        Dt = scalar(totalIter_[fieldi]);
     }
+
+    scalar beta = dt/Dt;
+
+    if (window() > 0)
+    {
+        const scalar w = window();
+
+        if (Dt - dt >= w)
+        {
+            beta = dt/w;
+        }
+    }
+
+    meanField = (1 - beta)*meanField + beta*baseField;
 }
 
 
 template<class Type>
 void Foam::functionObjects::fieldAverage::calculateMeanFields() const
 {
-    typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
-    typedef typename VolFieldType::Internal InternalType;
-    typedef GeometricField<Type, fvsPatchField, surfaceMesh> SurfaceFieldType;
-
     forAll(faItems_, fieldi)
     {
         if (faItems_[fieldi].mean())
         {
             const word& fieldName = faItems_[fieldi].fieldName();
 
-            if (obr_.foundObject<VolFieldType>(fieldName))
+            if (obr_.foundObject<VolField<Type>>(fieldName))
             {
-                calculateMeanFieldType<VolFieldType>(fieldi);
+                calculateMeanFieldType<VolField<Type>>(fieldi);
             }
-            else if (obr_.foundObject<InternalType>(fieldName))
+            else if (obr_.foundObject<VolInternalField<Type>>(fieldName))
             {
-                calculateMeanFieldType<InternalType>(fieldi);
+                calculateMeanFieldType<VolInternalField<Type>>(fieldi);
             }
-            else if (obr_.foundObject<SurfaceFieldType>(fieldName))
+            else if (obr_.foundObject<SurfaceField<Type>>(fieldName))
             {
-                calculateMeanFieldType<SurfaceFieldType>(fieldi);
+                calculateMeanFieldType<SurfaceField<Type>>(fieldi);
             }
         }
     }
@@ -312,41 +410,26 @@ void Foam::functionObjects::fieldAverage::calculatePrime2MeanFieldType
 template<class Type1, class Type2>
 void Foam::functionObjects::fieldAverage::calculatePrime2MeanFields() const
 {
-    typedef GeometricField<Type1, fvPatchField, volMesh> VolFieldType1;
-    typedef typename VolFieldType1::Internal InternalType1;
-    typedef GeometricField<Type1, fvsPatchField, surfaceMesh> SurfaceFieldType1;
-
-    typedef GeometricField<Type2, fvPatchField, volMesh> VolFieldType2;
-    typedef typename VolFieldType2::Internal InternalType2;
-    typedef GeometricField<Type2, fvsPatchField, surfaceMesh> SurfaceFieldType2;
-
     forAll(faItems_, fieldi)
     {
         if (faItems_[fieldi].prime2Mean())
         {
             const word& fieldName = faItems_[fieldi].fieldName();
 
-            if (obr_.foundObject<VolFieldType1>(fieldName))
-            {
-                calculatePrime2MeanFieldType<VolFieldType1, VolFieldType2>
-                (
-                    fieldi
-                );
-            }
-            else if (obr_.foundObject<InternalType1>(fieldName))
-            {
-                calculatePrime2MeanFieldType<InternalType1, InternalType2>
-                (
-                    fieldi
-                );
-            }
-            else if (obr_.foundObject<SurfaceFieldType1>(fieldName))
+            if (obr_.foundObject<VolField<Type1>>(fieldName))
             {
                 calculatePrime2MeanFieldType
-                <
-                    SurfaceFieldType1,
-                    SurfaceFieldType2
-                >(fieldi);
+                <VolField<Type1>, VolField<Type2>>(fieldi);
+            }
+            else if (obr_.foundObject<VolInternalField<Type1>>(fieldName))
+            {
+                calculatePrime2MeanFieldType
+                <VolInternalField<Type1>, VolInternalField<Type2>>(fieldi);
+            }
+            else if (obr_.foundObject<SurfaceField<Type1>>(fieldName))
+            {
+                calculatePrime2MeanFieldType
+                <SurfaceField<Type1>, SurfaceField<Type2>>(fieldi);
             }
         }
     }
@@ -372,40 +455,26 @@ void Foam::functionObjects::fieldAverage::addMeanSqrToPrime2MeanType
 template<class Type1, class Type2>
 void Foam::functionObjects::fieldAverage::addMeanSqrToPrime2Mean() const
 {
-    typedef GeometricField<Type1, fvPatchField, volMesh> VolFieldType1;
-    typedef typename VolFieldType1::Internal InternalType1;
-    typedef GeometricField<Type1, fvsPatchField, surfaceMesh> SurfaceFieldType1;
-
-    typedef GeometricField<Type2, fvPatchField, volMesh> VolFieldType2;
-    typedef typename VolFieldType2::Internal InternalType2;
-    typedef GeometricField<Type2, fvsPatchField, surfaceMesh> SurfaceFieldType2;
-
     forAll(faItems_, fieldi)
     {
         if (faItems_[fieldi].prime2Mean())
         {
             const word& fieldName = faItems_[fieldi].fieldName();
 
-            if (obr_.foundObject<VolFieldType1>(fieldName))
+            if (obr_.foundObject<VolField<Type1>>(fieldName))
             {
-                addMeanSqrToPrime2MeanType<VolFieldType1, VolFieldType2>
-                (
-                    fieldi
-                );
+                addMeanSqrToPrime2MeanType
+                <VolField<Type1>, VolField<Type2>>(fieldi);
             }
-            else if (obr_.foundObject<InternalType1>(fieldName))
+            else if (obr_.foundObject<VolInternalField<Type1>>(fieldName))
             {
-                addMeanSqrToPrime2MeanType<InternalType1, InternalType2>
-                (
-                    fieldi
-                );
+                addMeanSqrToPrime2MeanType
+                <VolInternalField<Type1>, VolInternalField<Type2>>(fieldi);
             }
-            else if (obr_.foundObject<SurfaceFieldType1>(fieldName))
+            else if (obr_.foundObject<SurfaceField<Type1>>(fieldName))
             {
-                addMeanSqrToPrime2MeanType<SurfaceFieldType1, SurfaceFieldType2>
-                (
-                    fieldi
-                );
+                addMeanSqrToPrime2MeanType
+                <SurfaceField<Type1>, SurfaceField<Type2>>(fieldi);
             }
         }
     }
@@ -429,25 +498,21 @@ void Foam::functionObjects::fieldAverage::writeFieldType
 template<class Type>
 void Foam::functionObjects::fieldAverage::writeFields() const
 {
-    typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
-    typedef typename VolFieldType::Internal InternalType;
-    typedef GeometricField<Type, fvsPatchField, surfaceMesh> SurfaceFieldType;
-
     forAll(faItems_, fieldi)
     {
         if (faItems_[fieldi].mean())
         {
             const word& fieldName = faItems_[fieldi].meanFieldName();
-            writeFieldType<VolFieldType>(fieldName);
-            writeFieldType<InternalType>(fieldName);
-            writeFieldType<SurfaceFieldType>(fieldName);
+            writeFieldType<VolField<Type>>(fieldName);
+            writeFieldType<VolInternalField<Type>>(fieldName);
+            writeFieldType<SurfaceField<Type>>(fieldName);
         }
         if (faItems_[fieldi].prime2Mean())
         {
             const word& fieldName = faItems_[fieldi].prime2MeanFieldName();
-            writeFieldType<VolFieldType>(fieldName);
-            writeFieldType<InternalType>(fieldName);
-            writeFieldType<SurfaceFieldType>(fieldName);
+            writeFieldType<VolField<Type>>(fieldName);
+            writeFieldType<VolInternalField<Type>>(fieldName);
+            writeFieldType<SurfaceField<Type>>(fieldName);
         }
     }
 }