ENH: support 'probes' style of ensemble output for sampledSets (#2389)

- with the special setFormat "probes", all of the sampled sets are
  treated more similarly to probes, with an ensemble output to raw
  probed format.

  This is of course less useful when the number of sampled points
  becomes very large.

src/sampling/probes/probes.C

@@ -50,6 +50,119 @@ namespace Foam
 }
 
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::probes::createProbeFiles(const wordList& fieldNames)
+{
+    // Open new output streams
+
+    bool needsNewFiles = false;
+    for (const word& fieldName : fieldNames)
+    {
+        if (!probeFilePtrs_.found(fieldName))
+        {
+            needsNewFiles = true;
+            break;
+        }
+    }
+
+    if (needsNewFiles && Pstream::master())
+    {
+        DebugInfo
+            << "Probing fields: " << fieldNames << nl
+            << "Probing locations: " << *this << nl
+            << endl;
+
+        // Put in undecomposed case
+        // (Note: gives problems for distributed data running)
+
+        fileName probeSubDir = name();
+        if (mesh_.name() != polyMesh::defaultRegion)
+        {
+            probeSubDir = probeSubDir/mesh_.name();
+        }
+
+        fileName probeDir
+        (
+            mesh_.time().globalPath()
+          / functionObject::outputPrefix
+          / probeSubDir
+          / mesh_.time().timeName()
+        );
+        probeDir.clean();  // Remove unneeded ".."
+
+        // Create directory if needed
+        Foam::mkDir(probeDir);
+
+        for (const word& fieldName : fieldNames)
+        {
+            if (probeFilePtrs_.found(fieldName))
+            {
+                // Safety
+                continue;
+            }
+
+            auto osPtr = autoPtr<OFstream>::New(probeDir/fieldName);
+            auto& os = *osPtr;
+
+            probeFilePtrs_.insert(fieldName, osPtr);
+
+            DebugInfo<< "open probe stream: " << os.name() << endl;
+
+            const unsigned int width(IOstream::defaultPrecision() + 7);
+
+            forAll(*this, probei)
+            {
+                os  << "# Probe " << probei << ' ' << operator[](probei);
+
+                if (processor_[probei] == -1)
+                {
+                    os  << "  # Not Found";
+                }
+                // Only for patchProbes
+                else if (probei < patchIDList_.size())
+                {
+                    const label patchi = patchIDList_[probei];
+                    if (patchi != -1)
+                    {
+                        const polyBoundaryMesh& bm = mesh_.boundaryMesh();
+                        if
+                        (
+                            patchi < bm.nNonProcessor()
+                         || processor_[probei] == Pstream::myProcNo()
+                        )
+                        {
+                            os  << " at patch " << bm[patchi].name();
+                        }
+                        os  << " with a distance of "
+                            << mag(operator[](probei)-oldPoints_[probei])
+                            << " m to the original point "
+                            << oldPoints_[probei];
+                    }
+                }
+
+                os  << nl;
+            }
+
+            os  << '#' << setw(IOstream::defaultPrecision() + 6)
+                << "Probe";
+
+            forAll(*this, probei)
+            {
+                if (includeOutOfBounds_ || processor_[probei] != -1)
+                {
+                    os  << ' ' << setw(width) << probei;
+                }
+            }
+            os  << nl;
+
+            os  << '#' << setw(IOstream::defaultPrecision() + 6)
+                << "Time" << endl;
+        }
+    }
+}
+
+
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
 void Foam::probes::findElements(const fvMesh& mesh)
@@ -243,99 +356,9 @@ Foam::label Foam::probes::prepare(unsigned request)
             }
         }
 
-        if (!(request & ACTION_WRITE))
+        if ((request & ACTION_WRITE) && !currentFields.empty())
         {
-            // No writing - can return now
-            return nFields;
-        }
-        else if (currentFields.empty())
-        {
-            // No new fields - can return now
-            return nFields;
-        }
-
-
-        // Have new fields - open streams for them
-
-        // Put in undecomposed case
-        // (Note: gives problems for distributed data running)
-
-        fileName probeSubDir = name();
-        if (mesh_.name() != polyMesh::defaultRegion)
-        {
-            probeSubDir = probeSubDir/mesh_.name();
-        }
-
-        fileName probeDir
-        (
-            mesh_.time().globalPath()
-          / functionObject::outputPrefix
-          / probeSubDir
-          / mesh_.time().timeName()
-        );
-        probeDir.clean();  // Remove unneeded ".."
-
-        // Create directory if needed
-        mkDir(probeDir);
-
-        for (const word& fieldName : currentFields.sortedToc())
-        {
-            auto osPtr = autoPtr<OFstream>::New(probeDir/fieldName);
-            auto& os = *osPtr;
-
-            probeFilePtrs_.insert(fieldName, osPtr);
-
-            DebugInfo<< "open probe stream: " << os.name() << endl;
-
-            const unsigned int w = IOstream::defaultPrecision() + 7;
-
-            forAll(*this, probei)
-            {
-                os  << "# Probe " << probei << ' ' << operator[](probei);
-
-                if (processor_[probei] == -1)
-                {
-                    os  << "  # Not Found";
-                }
-                // Only for patchProbes
-                else if (probei < patchIDList_.size())
-                {
-                    const label patchi = patchIDList_[probei];
-                    if (patchi != -1)
-                    {
-                        const polyBoundaryMesh& bm = mesh_.boundaryMesh();
-                        if
-                        (
-                            patchi < bm.nNonProcessor()
-                         || processor_[probei] == Pstream::myProcNo()
-                        )
-                        {
-                            os  << " at patch " << bm[patchi].name();
-                        }
-                        os  << " with a distance of "
-                            << mag(operator[](probei)-oldPoints_[probei])
-                            << " m to the original point "
-                            << oldPoints_[probei];
-                    }
-                }
-
-                os  << nl;
-            }
-
-            os  << '#' << setw(IOstream::defaultPrecision() + 6)
-                << "Probe";
-
-            forAll(*this, probei)
-            {
-                if (includeOutOfBounds_ || processor_[probei] != -1)
-                {
-                    os  << ' ' << setw(w) << probei;
-                }
-            }
-            os  << nl;
-
-            os  << '#' << setw(IOstream::defaultPrecision() + 6)
-                << "Time" << endl;
+            createProbeFiles(currentFields.sortedToc());
         }
     }
 

src/sampling/probes/probes.H

@@ -234,6 +234,9 @@ private:
 
     // Private Member Functions
 
+        //- Create new streams as required
+        void createProbeFiles(const wordList& fieldNames);
+
         //- Write field values
         template<class Type>
         void writeValues

src/sampling/probes/probesTemplates.C

@@ -136,16 +136,16 @@ void Foam::probes::writeValues
 {
     if (Pstream::master())
     {
-        const unsigned int w = IOstream::defaultPrecision() + 7;
+        const unsigned int width(IOstream::defaultPrecision() + 7);
         OFstream& os = *probeFilePtrs_[fieldName];
 
-        os  << setw(w) << timeValue;
+        os  << setw(width) << timeValue;
 
         forAll(values, probei)
         {
             if (includeOutOfBounds_ || processor_[probei] != -1)
             {
-                os  << ' ' << setw(w) << values[probei];
+                os  << ' ' << setw(width) << values[probei];
             }
         }
         os  << endl;

src/sampling/sampledSet/sampledSets/sampledSets.C

@@ -32,6 +32,7 @@ License
 #include "globalIndex.H"
 #include "volFields.H"
 #include "mapPolyMesh.H"
+#include "IOmanip.H"
 #include "IOobjectList.H"
 #include "UIndirectList.H"
 #include "ListOps.H"
@@ -76,6 +77,81 @@ Foam::autoPtr<Foam::coordSetWriter> Foam::sampledSets::newWriter
 }
 
 
+Foam::OFstream* Foam::sampledSets::createProbeFile(const word& fieldName)
+{
+    // Open new output stream
+
+    OFstream* osptr = probeFilePtrs_.lookup(fieldName, nullptr);
+
+    if (!osptr && Pstream::master())
+    {
+        // Put in undecomposed case
+        // (Note: gives problems for distributed data running)
+
+        fileName probeSubDir = name();
+        if (mesh_.name() != polyMesh::defaultRegion)
+        {
+            probeSubDir = probeSubDir/mesh_.name();
+        }
+
+        fileName probeDir
+        (
+            mesh_.time().globalPath()
+          / functionObject::outputPrefix
+          / probeSubDir
+          / mesh_.time().timeName()
+        );
+        probeDir.clean();  // Remove unneeded ".."
+
+        // Create directory if needed
+        Foam::mkDir(probeDir);
+
+        probeFilePtrs_.insert
+        (
+            fieldName,
+            autoPtr<OFstream>::New(probeDir/fieldName)
+        );
+        osptr = probeFilePtrs_.lookup(fieldName, nullptr);
+
+        if (osptr)
+        {
+            auto& os = *osptr;
+
+            DebugInfo<< "open probe stream: " << os.name() << endl;
+
+            const unsigned int width(IOstream::defaultPrecision() + 7);
+
+            label nPoints = 0;
+            forAll(*this, seti)
+            {
+                const coordSet& s = gatheredSets_[seti];
+
+                const pointField& pts = static_cast<const pointField&>(s);
+
+                for (const point& p : pts)
+                {
+                    os  << "# Probe " << nPoints++ << ' ' << p << nl;
+                }
+            }
+
+            os  << '#' << setw(IOstream::defaultPrecision() + 6)
+                << "Probe";
+
+            for (label probei = 0; probei < nPoints; ++probei)
+            {
+                os  << ' ' << setw(width) << probei;
+            }
+            os  << nl;
+
+            os  << '#' << setw(IOstream::defaultPrecision() + 6)
+                << "Time" << endl;
+        }
+    }
+
+    return osptr;
+}
+
+
 void Foam::sampledSets::gatherAllSets()
 {
     // Any writer references will become invalid
@@ -101,7 +177,7 @@ void Foam::sampledSets::gatherAllSets()
 }
 
 
-Foam::IOobjectList Foam::sampledSets::preCheckFields()
+Foam::IOobjectList Foam::sampledSets::preCheckFields(unsigned request)
 {
     wordList allFields;    // Just needed for warnings
     HashTable<wordHashSet> selected;
@@ -146,7 +222,7 @@ Foam::IOobjectList Foam::sampledSets::preCheckFields()
         }
     }
 
-    if (missed.size())
+    if (missed.size() && (request != ACTION_NONE))
     {
         WarningInFunction
             << nl
@@ -189,11 +265,29 @@ Foam::IOobjectList Foam::sampledSets::preCheckFields()
 
     const label nFields = selectedFieldNames_.size();
 
-    forAll(writers_, seti)
+    if (writeAsProbes_)
     {
-        coordSetWriter& writer = writers_[seti];
+        // Close streams for fields that no longer exist
+        forAllIters(probeFilePtrs_, iter)
+        {
+            if (!selectedFieldNames_.found(iter.key()))
+            {
+                DebugInfo
+                    << "close probe stream: "
+                    << iter()->name() << endl;
 
-        writer.nFields(nFields);
+                probeFilePtrs_.remove(iter);
+            }
+        }
+    }
+    else if ((request & ACTION_WRITE) != 0)
+    {
+        forAll(writers_, seti)
+        {
+            coordSetWriter& writer = writers_[seti];
+
+            writer.nFields(nFields);
+        }
     }
 
     return objects;
@@ -213,7 +307,7 @@ void Foam::sampledSets::initDict(const dictionary& dict, const bool initial)
     if (eptr && eptr->isDict())
     {
         PtrList<sampledSet> sampSets(eptr->dict().size());
-        if (initial)
+        if (initial && !writeAsProbes_)
         {
             writers_.resize(sampSets.size());
         }
@@ -247,7 +341,7 @@ void Foam::sampledSets::initDict(const dictionary& dict, const bool initial)
             sampSets.set(seti, sampSet);
 
             // Define writer, but do not attached
-            if (initial)
+            if (initial && !writeAsProbes_)
             {
                 writers_.set
                 (
@@ -263,7 +357,7 @@ void Foam::sampledSets::initDict(const dictionary& dict, const bool initial)
         }
 
         sampSets.resize(seti);
-        if (initial)
+        if (initial && !writeAsProbes_)
         {
             writers_.resize(seti);
         }
@@ -283,7 +377,7 @@ void Foam::sampledSets::initDict(const dictionary& dict, const bool initial)
         );
 
         PtrList<sampledSet> sampSets(input.size());
-        if (initial)
+        if (initial && !writeAsProbes_)
         {
             writers_.resize(sampSets.size());
         }
@@ -305,7 +399,7 @@ void Foam::sampledSets::initDict(const dictionary& dict, const bool initial)
             sampSets.set(seti, sampSet);
 
             // Define writer, but do not attached
-            if (initial)
+            if (initial && !writeAsProbes_)
             {
                 writers_.set
                 (
@@ -321,7 +415,7 @@ void Foam::sampledSets::initDict(const dictionary& dict, const bool initial)
         }
 
         sampSets.resize(seti);
-        if (initial)
+        if (initial && !writeAsProbes_)
         {
             writers_.resize(seti);
         }
@@ -351,6 +445,7 @@ Foam::sampledSets::sampledSets
     verbose_(false),
     onExecute_(false),
     needsCorrect_(false),
+    writeAsProbes_(false),
     outputPath_
     (
         time_.globalPath()/functionObject::outputPrefix/name
@@ -359,8 +454,9 @@ Foam::sampledSets::sampledSets
     samplePointScheme_(),
     writeFormat_(),
     writeFormatOptions_(dict.subOrEmptyDict("formatOptions")),
-    writers_(),
     selectedFieldNames_(),
+    writers_(),
+    probeFilePtrs_(),
     gatheredSets_(),
     gatheredSorting_(),
     globalIndices_()
@@ -369,7 +465,6 @@ Foam::sampledSets::sampledSets
     {
         outputPath_ /= mesh_.name();
     }
-
     outputPath_.clean();  // Remove unneeded ".."
 
     read(dict);
@@ -391,6 +486,7 @@ Foam::sampledSets::sampledSets
     verbose_(false),
     onExecute_(false),
     needsCorrect_(false),
+    writeAsProbes_(false),
     outputPath_
     (
         time_.globalPath()/functionObject::outputPrefix/name
@@ -399,8 +495,9 @@ Foam::sampledSets::sampledSets
     samplePointScheme_(),
     writeFormat_(),
     writeFormatOptions_(dict.subOrEmptyDict("formatOptions")),
-    writers_(),
     selectedFieldNames_(),
+    writers_(),
+    probeFilePtrs_(),
     gatheredSets_(),
     gatheredSorting_(),
     globalIndices_()
@@ -456,6 +553,15 @@ bool Foam::sampledSets::read(const dictionary& dict)
     {
         dict.readEntry("setFormat", writeFormat_);
     }
+
+    // Hard-coded handling of ensemble 'probes' writer
+    writeAsProbes_ = ("probes" == writeFormat_);
+    if (!writeAsProbes_)
+    {
+        // Close all streams
+        probeFilePtrs_.clear();
+    }
+
     // const dictionary formatOptions(dict.subOrEmptyDict("formatOptions"));
     // Writer type and format options
     // const word writerType =
@@ -472,17 +578,28 @@ bool Foam::sampledSets::read(const dictionary& dict)
         fieldSelection_.uniq();
 
         // Report
-        forAll(*this, seti)
+        if (writeAsProbes_)
         {
-            const sampledSet& s = (*this)[seti];
+            Info<< "Sampled set as probes ensemble:" << nl;
 
-            if (!seti)
+            forAll(*this, seti)
             {
-                Info<< "Sampled set:" << nl;
+                const sampledSet& s = (*this)[seti];
+                Info<< "  " << s.name();
             }
+            Info<< nl;
+        }
+        else
+        {
+            Info<< "Sampled set:" << nl;
 
-            Info<< "    " << s.name() << " -> " << writers_[seti].type()
-                << nl;
+            forAll(*this, seti)
+            {
+                const sampledSet& s = (*this)[seti];
+
+                Info<< "    " << s.name() << " -> "
+                    << writers_[seti].type() << nl;
+            }
         }
 
         Info<< endl;
@@ -504,6 +621,11 @@ bool Foam::sampledSets::read(const dictionary& dict)
         Pout<< ')' << endl;
     }
 
+    if (writeAsProbes_)
+    {
+        (void) preCheckFields(ACTION_NONE);
+    }
+
     // FUTURE:
     // Ensure all sets and merge information are expired
     // expire(true);
@@ -532,7 +654,7 @@ bool Foam::sampledSets::performAction(unsigned request)
     // Determine availability of fields.
     // Count number of fields (only seems to be needed for VTK legacy)
 
-    IOobjectList objects = preCheckFields();
+    IOobjectList objects = preCheckFields(request);
 
     const label nFields = selectedFieldNames_.size();
 
@@ -543,34 +665,40 @@ bool Foam::sampledSets::performAction(unsigned request)
     }
 
     // Update writers
-
-    forAll(*this, seti)
+    if (!writeAsProbes_)
     {
-        const coordSet& s = gatheredSets_[seti];
-
-        if (request & ACTION_WRITE)
+        forAll(*this, seti)
         {
-            coordSetWriter& writer = writers_[seti];
+            const coordSet& s = gatheredSets_[seti];
 
-            if (writer.needsUpdate())
+            if ((request & ACTION_WRITE) != 0)
             {
-                writer.setCoordinates(s);
-            }
+                coordSetWriter& writer = writers_[seti];
 
-            if (writer.buffering())
-            {
-                writer.open
-                (
-                    outputPath_
-                  / word(s.name() + coordSetWriter::suffix(selectedFieldNames_))
-                );
-            }
-            else
-            {
-                writer.open(outputPath_/s.name());
-            }
+                if (writer.needsUpdate())
+                {
+                    writer.setCoordinates(s);
+                }
 
-            writer.beginTime(mesh_.time());
+                if (writer.buffering())
+                {
+                    writer.open
+                    (
+                        outputPath_
+                      / word
+                        (
+                            s.name()
+                          + coordSetWriter::suffix(selectedFieldNames_)
+                        )
+                    );
+                }
+                else
+                {
+                    writer.open(outputPath_/s.name());
+                }
+
+                writer.beginTime(mesh_.time());
+            }
         }
     }
 
@@ -584,19 +712,22 @@ bool Foam::sampledSets::performAction(unsigned request)
 
 
     // Finish this time step
-    forAll(writers_, seti)
+    if (!writeAsProbes_)
     {
-        // Write geometry if no fields were written so that we still
-        // can have something to look at
-
-        if (request & ACTION_WRITE)
+        forAll(writers_, seti)
         {
-            /// if (!writers_[seti].wroteData())
-            /// {
-            ///     writers_[seti].write();
-            /// }
+            // Write geometry if no fields were written so that we still
+            // can have something to look at
 
-            writers_[seti].endTime();
+            if ((request & ACTION_WRITE) != 0)
+            {
+                /// if (!writers_[seti].wroteData())
+                /// {
+                ///     writers_[seti].write();
+                /// }
+
+                writers_[seti].endTime();
+            }
         }
     }
 

src/sampling/sampledSet/sampledSets/sampledSets.H

@@ -93,6 +93,10 @@ Description
         formatOptions | dictionary of format options        | no  |
     \endtable
 
+Note
+    Special setFormat \c probes can be used to output ensemble results
+    in a format similar to the probes function object.
+
 SourceFiles
     sampledSets.C
     sampledSetsImpl.C
@@ -103,6 +107,8 @@ SourceFiles
 #define Foam_sampledSets_H
 
 #include "fvMeshFunctionObject.H"
+#include "HashPtrTable.H"
+#include "OFstream.H"
 #include "sampledSet.H"
 #include "meshSearch.H"
 #include "coordSet.H"
@@ -157,6 +163,9 @@ class sampledSets
         //- Correct meshSearch and update sets
         bool needsCorrect_;
 
+        //- Output in raw format like 'probes' does
+        bool writeAsProbes_;
+
         //- Output path
         fileName outputPath_;
 
@@ -181,11 +190,14 @@ class sampledSets
 
     // Output control
 
+        //- Current list of field names selected for sampling
+        DynamicList<word> selectedFieldNames_;
+
         //- The coordSet writers (one per sampled set)
         PtrList<coordSetWriter> writers_;
 
-        //- Current list of field names selected for sampling
-        DynamicList<word> selectedFieldNames_;
+        //- Current open files (non-empty on master only) when writing as probes
+        HashPtrTable<OFstream> probeFilePtrs_;
 
 
     // Merging
@@ -203,6 +215,9 @@ class sampledSets
 
     // Private Member Functions
 
+        //- Get or create new stream as required (on master)
+        OFstream* createProbeFile(const word& fieldName);
+
         //- A new coordSet writer, with per-set formatOptions
         static autoPtr<coordSetWriter> newWriter
         (
@@ -218,7 +233,7 @@ class sampledSets
         //  Returns empty IOobjectList if loadFromFiles_ is not active
         //
         //  Adjusts selectedFieldNames_
-        IOobjectList preCheckFields();
+        IOobjectList preCheckFields(unsigned request);
 
         //- Setup the sets (optional writers)
         void initDict(const dictionary& dict, const bool initial);

src/sampling/sampledSet/sampledSets/sampledSetsImpl.C

@@ -107,6 +107,7 @@ void Foam::sampledSets::performAction
 )
 {
     const word& fieldName = fld.name();
+    const scalar timeValue = fld.time().timeOutputValue();
 
     // The interpolator for this field
     autoPtr<interpolation<Type>> interpPtr;
@@ -116,6 +117,19 @@ void Foam::sampledSets::performAction
         interpPtr.reset(interpolation<Type>::New(samplePointScheme_, fld));
     }
 
+    const unsigned int width(IOstream::defaultPrecision() + 7);
+    OFstream* osptr = nullptr;
+
+    if (writeAsProbes_ && (request & ACTION_WRITE))
+    {
+        osptr = createProbeFile(fieldName);
+
+        if (Pstream::master() && osptr)
+        {
+            (*osptr) << setw(width) << timeValue;
+        }
+    }
+
     // Ensemble min/max/avg values
     Type avgEnsemble = Zero;
     label sizeEnsemble = 0;
@@ -215,12 +229,31 @@ void Foam::sampledSets::performAction
                 << "    max: " << limits.max() << nl << nl;
         }
 
-        if (request & ACTION_WRITE)
+        if ((request & ACTION_WRITE) != 0)
         {
-            writeCoordSet<Type>(writers_[seti], values, fieldName);
+            if (writeAsProbes_)
+            {
+                if (osptr)
+                {
+                    for (const Type& val : values)
+                    {
+                        (*osptr) << ' ' << setw(width) << val;
+                    }
+                }
+            }
+            else
+            {
+                writeCoordSet<Type>(writers_[seti], values, fieldName);
+            }
         }
     }
 
+    // Finish probes write
+    if (Pstream::master() && osptr)
+    {
+        (*osptr) << endl;
+    }
+
     if (sizeEnsemble)
     {
         avgEnsemble /= sizeEnsemble;