diff --git a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C
index 1d7cd88a29..5f429e1214 100644
--- a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C
+++ b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C
@@ -32,7 +32,7 @@ License
 #include "coupledPolyPatch.H"
 #include "sampledSurface.H"
 #include "mergePoints.H"
-#include "indirectPrimitivePatch.H"
+#include "uindirectPrimitivePatch.H"
 #include "PatchTools.H"
 #include "addToRunTimeSelectionTable.H"
 
@@ -130,115 +130,225 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::obr() const
 
 void Foam::functionObjects::fieldValues::surfaceFieldValue::setFaceZoneFaces()
 {
-    const label zoneId = mesh_.faceZones().findZoneID(regionName_);
+    // Indices for all matches, already sorted
+    const labelList zoneIds
+    (
+        mesh_.faceZones().indices(selectionNames_)
+    );
 
-    if (zoneId < 0)
+    // Total number of faces selected
+    label numFaces = 0;
+    for (const label zoneId : zoneIds)
     {
-        FatalErrorInFunction
-            << type() << " " << name() << ": "
-            << regionTypeNames_[regionType_] << '(' << regionName_ << "):" << nl
-            << "    Unknown face zone name: " << regionName_
-            << ". Valid face zones are: " << mesh_.faceZones().names()
-            << nl << exit(FatalError);
+        numFaces += mesh_.faceZones()[zoneId].size();
     }
 
-    const faceZone& fZone = mesh_.faceZones()[zoneId];
-
-    DynamicList<label> faceIds(fZone.size());
-    DynamicList<label> facePatchIds(fZone.size());
-    DynamicList<bool> faceFlip(fZone.size());
-
-    forAll(fZone, i)
+    if (zoneIds.empty())
     {
-        const label facei = fZone[i];
+        FatalErrorInFunction
+            << type() << ' ' << name() << ": "
+            << regionTypeNames_[regionType_] << '(' << regionName_ << "):" << nl
+            << "    No matching face zone(s): "
+            << flatOutput(selectionNames_)  << nl
+            << "    Known face zones: "
+            << flatOutput(mesh_.faceZones().names()) << nl
+            << exit(FatalError);
+    }
 
-        label faceId = -1;
-        label facePatchId = -1;
-        if (mesh_.isInternalFace(facei))
+    // Could also check this
+    #if 0
+    if (!returnReduce(bool(numFaces), orOp<bool>()))
+    {
+        WarningInFunction
+            << type() << ' ' << name() << ": "
+            << regionTypeNames_[regionType_] << '(' << regionName_ << "):" << nl
+            << "    The faceZone specification: "
+            << flatOutput(selectionNames_) << nl
+            << "    resulted in 0 faces" << nl
+            << exit(FatalError);
+    }
+    #endif
+
+    faceId_.resize(numFaces);
+    facePatchId_.resize(numFaces);
+    faceFlip_.resize(numFaces);
+
+    numFaces = 0;
+
+    for (const label zoneId : zoneIds)
+    {
+        const faceZone& fZone = mesh_.faceZones()[zoneId];
+
+        forAll(fZone, i)
         {
-            faceId = facei;
-            facePatchId = -1;
-        }
-        else
-        {
-            facePatchId = mesh_.boundaryMesh().whichPatch(facei);
-            const polyPatch& pp = mesh_.boundaryMesh()[facePatchId];
-            if (isA<coupledPolyPatch>(pp))
+            const label meshFacei = fZone[i];
+            const bool isFlip = fZone.flipMap()[i];
+
+            // Internal faces
+            label faceId = meshFacei;
+            label facePatchId = -1;
+
+            // Boundary faces
+            if (!mesh_.isInternalFace(meshFacei))
             {
-                if (refCast<const coupledPolyPatch>(pp).owner())
+                facePatchId = mesh_.boundaryMesh().whichPatch(meshFacei);
+                const polyPatch& pp = mesh_.boundaryMesh()[facePatchId];
+
+                if (isA<coupledPolyPatch>(pp))
                 {
-                    faceId = pp.whichFace(facei);
+                    if (refCast<const coupledPolyPatch>(pp).owner())
+                    {
+                        faceId = pp.whichFace(meshFacei);
+                    }
+                    else
+                    {
+                        faceId = -1;
+                    }
+                }
+                else if (!isA<emptyPolyPatch>(pp))
+                {
+                    faceId = meshFacei - pp.start();
                 }
                 else
                 {
                     faceId = -1;
+                    facePatchId = -1;
                 }
             }
-            else if (!isA<emptyPolyPatch>(pp))
-            {
-                faceId = facei - pp.start();
-            }
-            else
-            {
-                faceId = -1;
-                facePatchId = -1;
-            }
-        }
 
-        if (faceId >= 0)
-        {
-            faceIds.append(faceId);
-            facePatchIds.append(facePatchId);
-            faceFlip.append(fZone.flipMap()[i] ? true : false);
+            if (faceId >= 0)
+            {
+                faceId_[numFaces] = faceId;
+                facePatchId_[numFaces] = facePatchId;
+                faceFlip_[numFaces] = isFlip;
+
+                ++numFaces;
+            }
         }
     }
 
-    faceId_.transfer(faceIds);
-    facePatchId_.transfer(facePatchIds);
-    faceFlip_.transfer(faceFlip);
+    // Shrink to size used
+    faceId_.resize(numFaces);
+    facePatchId_.resize(numFaces);
+    faceFlip_.resize(numFaces);
     nFaces_ = returnReduce(faceId_.size(), sumOp<label>());
-
-    if (debug)
-    {
-        Pout<< "Original face zone size = " << fZone.size()
-            << ", new size = " << faceId_.size() << endl;
-    }
 }
 
 
 void Foam::functionObjects::fieldValues::surfaceFieldValue::setPatchFaces()
 {
-    const label patchid = mesh_.boundaryMesh().findPatchID(regionName_);
+    // Patch indices for all matches
+    labelList patchIds;
 
-    if (patchid < 0)
+    // Total number of faces selected
+    label numFaces = 0;
+
+    labelList selected
+    (
+        mesh_.boundaryMesh().patchSet
+        (
+            selectionNames_,
+            false  // warnNotFound - we do that ourselves
+        ).sortedToc()
+    );
+
+    DynamicList<label> bad;
+    for (const label patchi : selected)
+    {
+        const polyPatch& pp = mesh_.boundaryMesh()[patchi];
+
+        if (isA<emptyPolyPatch>(pp))
+        {
+            bad.append(patchi);
+        }
+        else
+        {
+            numFaces += pp.size();
+        }
+    }
+
+    if (bad.size())
+    {
+        label nGood = (selected.size() - bad.size());
+
+        auto& os = (nGood > 0 ? WarningInFunction : FatalErrorInFunction);
+
+        os  << "Cannot sample an empty patch" << nl;
+
+        for (const label patchi : bad)
+        {
+            os  << "    "
+                << mesh_.boundaryMesh()[patchi].name() << nl;
+        }
+
+        if (nGood)
+        {
+            os  << "No non-empty patches selected" << endl
+                << exit(FatalError);
+        }
+        else
+        {
+            os  << "Selected " << nGood << " non-empty patches" << nl;
+        }
+
+        patchIds.resize(nGood);
+        nGood = 0;
+        for (const label patchi : selected)
+        {
+            if (!bad.found(patchi))
+            {
+                patchIds[nGood] = patchi;
+                ++nGood;
+            }
+        }
+    }
+    else
+    {
+        patchIds = std::move(selected);
+    }
+
+    if (patchIds.empty())
     {
         FatalErrorInFunction
-            << type() << " " << name() << ": "
+            << type() << ' ' << name() << ": "
             << regionTypeNames_[regionType_] << '(' << regionName_ << "):" << nl
-            << "    Unknown patch name: " << regionName_
-            << ". Valid patch names are: "
+            << "    No matching patch name(s): "
+            << flatOutput(selectionNames_)  << nl
+            << "    Known patch names:" << nl
             << mesh_.boundaryMesh().names() << nl
             << exit(FatalError);
     }
 
-    const polyPatch& pp = mesh_.boundaryMesh()[patchid];
-
-    label nFaces = pp.size();
-    if (isA<emptyPolyPatch>(pp))
+    // Could also check this
+    #if 0
+    if (!returnReduce(bool(numFaces), orOp<bool>()))
     {
-        nFaces = 0;
+        WarningInFunction
+            << type() << ' ' << name() << ": "
+            << regionTypeNames_[regionType_] << '(' << regionName_ << "):" << nl
+            << "    The patch specification: "
+            << flatOutput(selectionNames_) << nl
+            << "    resulted in 0 faces" << nl
+            << exit(FatalError);
     }
+    #endif
 
-    faceId_.setSize(nFaces);
-    facePatchId_.setSize(nFaces);
-    faceFlip_.setSize(nFaces);
+    faceId_.resize(numFaces);
+    facePatchId_.resize(numFaces);
+    faceFlip_.resize(numFaces);
     nFaces_ = returnReduce(faceId_.size(), sumOp<label>());
 
-    forAll(faceId_, facei)
+    numFaces = 0;
+    for (const label patchi : patchIds)
     {
-        faceId_[facei] = facei;
-        facePatchId_[facei] = patchid;
-        faceFlip_[facei] = false;
+        const polyPatch& pp = mesh_.boundaryMesh()[patchi];
+        const label len = pp.size();
+
+        SubList<label>(faceId_, len, numFaces) = identity(len);
+        SubList<label>(facePatchId_, len, numFaces) = patchi;
+        SubList<bool>(faceFlip_, len, numFaces) = false;
+
+        numFaces += len;
     }
 }
 
@@ -252,24 +362,25 @@ void Foam::functionObjects::fieldValues::surfaceFieldValue::combineMeshGeometry
     List<faceList> allFaces(Pstream::nProcs());
     List<pointField> allPoints(Pstream::nProcs());
 
-    labelList globalFacesIs(faceId_);
-    forAll(globalFacesIs, i)
     {
-        if (facePatchId_[i] != -1)
+        IndirectList<face> selectedFaces(mesh_.faces(), labelList(faceId_));
+        labelList& meshFaceIds = selectedFaces.addressing();
+
+        forAll(meshFaceIds, i)
         {
             const label patchi = facePatchId_[i];
-            globalFacesIs[i] += mesh_.boundaryMesh()[patchi].start();
+            if (patchi != -1)
+            {
+                meshFaceIds[i] += mesh_.boundaryMesh()[patchi].start();
+            }
         }
-    }
 
-    // Add local faces and points to the all* lists
-    indirectPrimitivePatch pp
-    (
-        IndirectList<face>(mesh_.faces(), globalFacesIs),
-        mesh_.points()
-    );
-    allFaces[Pstream::myProcNo()] = pp.localFaces();
-    allPoints[Pstream::myProcNo()] = pp.localPoints();
+        // Add local faces and points to the all* lists
+        uindirectPrimitivePatch pp(selectedFaces, mesh_.points());
+
+        allFaces[Pstream::myProcNo()] = pp.localFaces();
+        allPoints[Pstream::myProcNo()] = pp.localPoints();
+    }
 
     Pstream::gatherList(allFaces);
     Pstream::gatherList(allPoints);
@@ -283,53 +394,41 @@ void Foam::functionObjects::fieldValues::surfaceFieldValue::combineMeshGeometry
         nPoints += allPoints[proci].size();
     }
 
-    faces.setSize(nFaces);
-    points.setSize(nPoints);
+    faces.resize(nFaces);
+    points.resize(nPoints);
 
     nFaces = 0;
     nPoints = 0;
 
-    // My own data first
+    // My data first
     {
-        const faceList& fcs = allFaces[Pstream::myProcNo()];
-        for (const face& f : fcs)
+        for (const face& f : allFaces[Pstream::myProcNo()])
         {
-            face& newF = faces[nFaces++];
-            newF.setSize(f.size());
-            forAll(f, fp)
-            {
-                newF[fp] = f[fp] + nPoints;
-            }
+            faces[nFaces++] = offsetOp<face>()(f, nPoints);
         }
 
-        const pointField& pts = allPoints[Pstream::myProcNo()];
-        for (const point& pt : pts)
+        for (const point& p : allPoints[Pstream::myProcNo()])
         {
-            points[nPoints++] = pt;
+            points[nPoints++] = p;
         }
     }
 
     // Other proc data follows
     forAll(allFaces, proci)
     {
-        if (proci != Pstream::myProcNo())
+        if (proci == Pstream::myProcNo())
         {
-            const faceList& fcs = allFaces[proci];
-            for (const face& f : fcs)
-            {
-                face& newF = faces[nFaces++];
-                newF.setSize(f.size());
-                forAll(f, fp)
-                {
-                    newF[fp] = f[fp] + nPoints;
-                }
-            }
+            continue;
+        }
 
-            const pointField& pts = allPoints[proci];
-            for (const point& pt : pts)
-            {
-                points[nPoints++] = pt;
-            }
+        for (const face& f : allFaces[proci])
+        {
+            faces[nFaces++] = offsetOp<face>()(f, nPoints);
+        }
+
+        for (const point& p : allPoints[proci])
+        {
+            points[nPoints++] = p;
         }
     }
 
@@ -383,11 +482,7 @@ combineSurfaceGeometry
             PatchTools::gatherAndMerge
             (
                 mergeDim,
-                primitivePatch
-                (
-                    SubList<face>(s.faces(), s.faces().size()),
-                    s.points()
-                ),
+                primitivePatch(SubList<face>(s.faces()), s.points()),
                 points,
                 faces,
                 pointsMap
@@ -413,11 +508,7 @@ combineSurfaceGeometry
             PatchTools::gatherAndMerge
             (
                 mergeDim,
-                primitivePatch
-                (
-                    SubList<face>(s.faces(), s.faces().size()),
-                    s.points()
-                ),
+                primitivePatch(SubList<face>(s.faces()), s.points()),
                 points,
                 faces,
                 pointsMap
@@ -522,7 +613,7 @@ bool Foam::functionObjects::fieldValues::surfaceFieldValue::update()
     if (nFaces_ == 0)
     {
         FatalErrorInFunction
-            << type() << " " << name() << ": "
+            << type() << ' ' << name() << ": "
             << regionTypeNames_[regionType_] << '(' << regionName_ << "):" << nl
             << "    Region has no faces" << exit(FatalError);
     }
@@ -548,7 +639,7 @@ void Foam::functionObjects::fieldValues::surfaceFieldValue::writeFileHeader
     if (canWriteHeader() && (operation_ != opNone))
     {
         writeCommented(os, "Region type : ");
-        os << regionTypeNames_[regionType_] << " " << regionName_ << endl;
+        os << regionTypeNames_[regionType_] << ' ' << regionName_ << nl;
 
         writeHeaderValue(os, "Faces", nFaces_);
         writeHeaderValue(os, "Area", totalArea_);
@@ -570,7 +661,7 @@ void Foam::functionObjects::fieldValues::surfaceFieldValue::writeFileHeader
         for (const word& fieldName : fields_)
         {
             os  << tab << operationTypeNames_[operation_]
-                << "(" << fieldName << ")";
+                << '(' << fieldName << ')';
         }
 
         os  << endl;
@@ -821,9 +912,10 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::surfaceFieldValue
             true  // Failsafe behaviour
         )
     ),
-    weightFieldName_("none"),
     needsUpdate_(true),
     writeArea_(false),
+    selectionNames_(),
+    weightFieldName_("none"),
     totalArea_(0),
     nFaces_(0),
     faceId_(),
@@ -854,9 +946,10 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::surfaceFieldValue
             true  // Failsafe behaviour
         )
     ),
-    weightFieldName_("none"),
     needsUpdate_(true),
     writeArea_(false),
+    selectionNames_(),
+    weightFieldName_("none"),
     totalArea_(0),
     nFaces_(0),
     faceId_(),
@@ -876,18 +969,54 @@ bool Foam::functionObjects::fieldValues::surfaceFieldValue::read
 {
     fieldValue::read(dict);
 
-    weightFieldName_ = "none";
     needsUpdate_ = true;
     writeArea_ = dict.getOrDefault("writeArea", false);
+    weightFieldName_ = "none";
     totalArea_ = 0;
     nFaces_ = 0;
     faceId_.clear();
     facePatchId_.clear();
     faceFlip_.clear();
-    sampledPtr_.clear();
-    surfaceWriterPtr_.clear();
+    sampledPtr_.reset(nullptr);
+    surfaceWriterPtr_.reset(nullptr);
+
+    // Can have "name" (word) and/or "names" (wordRes)
+    //
+    // If "names" exists AND contains a literal (non-regex) that can be used
+    // as a suitable value for "name", the "name" entry becomes optional.
+
+    regionName_.clear();
+    selectionNames_.clear();
+
+    {
+        dict.readIfPresent("names", selectionNames_);
+
+        for (const auto& item : selectionNames_)
+        {
+            if (item.isLiteral())
+            {
+                regionName_ = item;
+                break;
+            }
+        }
+
+        // Mandatory if we didn't pick up a value from selectionNames_
+        dict.readEntry
+        (
+            "name",
+            regionName_,
+            keyType::LITERAL,
+            regionName_.empty()
+        );
+
+        // Ensure there is always content for selectionNames_
+        if (selectionNames_.empty())
+        {
+            selectionNames_.resize(1);
+            selectionNames_.first() = regionName_;
+        }
+    }
 
-    dict.readEntry("name", regionName_);
 
     // Create sampled surface, but leave 'expired' (ie, no update) since it
     // may depend on fields or data that do not yet exist
@@ -901,7 +1030,7 @@ bool Foam::functionObjects::fieldValues::surfaceFieldValue::read
         );
     }
 
-    Info<< type() << " " << name() << ":" << nl
+    Info<< type() << ' ' << name() << ':' << nl
         << "    operation     = ";
 
     if (postOperation_ != postOpNone)
diff --git a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.H b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.H
index 1dde638e55..5687e9f4a4 100644
--- a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.H
+++ b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.H
@@ -31,17 +31,17 @@ Group
     grpFieldFunctionObjects
 
 Description
-    Provides a 'face regionType' variant of the \c fieldValues function object.
+    A \c face regionType variant of the \c fieldValues function object.
 
     Given a list of user-specified fields and a selection of mesh (or general
     surface) faces, a number of operations can be performed, such as sums,
     averages and integrations.
 
     For example, to calculate the volumetric or mass flux across a patch,
-    apply the 'sum' operator to the flux field (typically \c phi).
+    apply the \c sum operator to the flux field (typically \c phi).
 
 Usage
-    Minimal example by using \c system/controlDict.functions:
+    A minimal example:
     \verbatim
     surfaceFieldValuePatch1
     {
@@ -56,6 +56,8 @@ Usage
         name            <patch>;
 
         // Optional entries (runtime modifiable)
+        names           (<patch-name> <patch-regex>);
+
         postOperation   none;
         weightField     alpha1;
         scaleFactor     1.0;
@@ -79,6 +81,8 @@ Usage
         name            <faceZone>;
 
         // Optional entries (runtime modifiable)
+        names           (<zone-name> <zone-regex>);
+
         postOperation   none;
         weightField     alpha1;
         scaleFactor     1.0;
@@ -92,12 +96,13 @@ Usage
 
     where the entries mean:
     \table
-      Property     | Description                        | Type | Req'd | Dflt
+      Property     | Description                        | Type | Reqd | Dflt
       type         | Type name: surfaceFieldValue       | word |  yes  | -
-      libs         | Library name: fieldFunctionObjects | word |  yes  | -
-      fields       | Names of operand fields            | wordList | yes | -
+      libs         | Libraries: fieldFunctionObjects    | wordList |  yes  | -
       regionType   | Face regionType: see below         | word |  yes  | -
-      name         | Name for regionType                | word |  yes  | -
+      fields       | Names of operand fields            | wordList | yes | -
+      name         | Name of the regionType             | word |  yes  | -
+      names        | Extended selection                 | word/regex list | no | -
       operation    | Operation type: see below          | word |  yes  | -
       postOperation | Post-operation type: see below    | word |  no   | none
       weightField  | Name of field to apply weighting   | word |  no   | none
@@ -112,9 +117,9 @@ Usage
 
     Options for the \c regionType entry:
     \plaintable
-      faceZone              | The \b name entry to specify the faceZone
-      patch                 | The \b name entry to specify the patch
-      functionObjectSurface | The \b name entry to specify a polySurface
+      faceZone      | The \b name entry specifies a faceZone. Supports \b names
+      patch         | The \b name entry specifies a patch. Supports \b names
+      functionObjectSurface | The \b name entry specifies a polySurface
       sampledSurface        | A \b sampledSurfaceDict sub-dictionary and \b name
     \endplaintable
 
@@ -156,12 +161,16 @@ Usage
     Usage by the \c postProcess utility is not available.
 
 Note
+    - Some types (eg, patch or faceZone) support the selection of multiple
+      entries, which can be specified as list of words or regular expressions
+      by \b names entry.<br>
+      If the \b names enty exists \em and contains a literal that can be used
+      as a suitable value for \b name, the \b name entry becomes optional
+      instead of being mandatory.
     - The values reported by the \c areaNormalAverage and \c areaNormalIntegrate
       operations are written as the first component of a field with the same
       rank as the input field.
-    - Faces on empty patches get ignored.
-    - If the field is a volField the \c faceZone
-      can only consist of boundary faces.
+    - Faces on empty patches are ignored.
     - Using \c functionObjectSurface:
       - The keyword %subRegion should not be used to select surfaces.
         Instead specify the regionType 'functionObjectSurface' and provide
@@ -242,8 +251,8 @@ public:
         //- Region type enumeration
         enum regionTypes
         {
-            stFaceZone = 0x01,   //!< Calculate on a faceZone
-            stPatch    = 0x02,   //!< Calculate on a patch
+            stFaceZone = 0x01,   //!< Calculate with faceZone(s)
+            stPatch    = 0x02,   //!< Calculate with patch(es)
             stObject   = 0x11,   //!< Calculate with function object surface
             stSampled  = 0x12    //!< Sample onto surface and calculate
         };
@@ -385,15 +394,18 @@ protected:
         //- Optional post-evaluation operation
         postOperationType postOperation_;
 
-        //- Weight field name - optional
-        word weightFieldName_;
-
         //- Track if the surface needs an update
         bool needsUpdate_;
 
         //- Optionally write the area of the surfaceFieldValue
         bool writeArea_;
 
+        //- Extended selections
+        wordRes selectionNames_;
+
+        //- Weight field name - optional
+        word weightFieldName_;
+
         //- Total area of the surfaceFieldValue
         scalar totalArea_;
 
@@ -401,17 +413,20 @@ protected:
         label nFaces_;
 
 
-        // If operating on mesh faces (faceZone, patch)
+    // If operating on mesh faces (faceZone, patch)
 
-            //- Local list of face IDs
-            labelList faceId_;
+        //- Local list of face IDs
+        labelList faceId_;
 
-            //- Local list of patch ID per face
-            labelList facePatchId_;
+        //- Local list of patch ID per face
+        labelList facePatchId_;
 
-            //- List representing the face flip map
-            //  (false: use as-is, true: negate)
-            boolList faceFlip_;
+        //- List representing the face flip map
+        //  (false: use as-is, true: negate)
+        boolList faceFlip_;
+
+
+    // Demand-driven
 
         //- The sampledSurface (when operating on sampledSurface)
         autoPtr<sampledSurface> sampledPtr_;
diff --git a/tutorials/lagrangian/reactingParcelFoam/verticalChannel/Allrun b/tutorials/lagrangian/reactingParcelFoam/verticalChannel/Allrun
index 9ec0596bab..a5d9a2cb05 100755
--- a/tutorials/lagrangian/reactingParcelFoam/verticalChannel/Allrun
+++ b/tutorials/lagrangian/reactingParcelFoam/verticalChannel/Allrun
@@ -3,15 +3,12 @@ cd "${0%/*}" || exit                                # Run from this directory
 . ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
 #------------------------------------------------------------------------------
 
-# create mesh
-runApplication blockMesh
-
 restore0Dir
 
-# initialise with potentialFoam solution
+runApplication blockMesh
+
 runApplication potentialFoam
 
-# run the solver
 runApplication $(getApplication)
 
 #------------------------------------------------------------------------------
diff --git a/tutorials/lagrangian/reactingParcelFoam/verticalChannel/Allrun-parallel b/tutorials/lagrangian/reactingParcelFoam/verticalChannel/Allrun-parallel
new file mode 100755
index 0000000000..8bb64db819
--- /dev/null
+++ b/tutorials/lagrangian/reactingParcelFoam/verticalChannel/Allrun-parallel
@@ -0,0 +1,16 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
+#------------------------------------------------------------------------------
+
+restore0Dir
+
+runApplication blockMesh
+
+runApplication potentialFoam
+
+runApplication decomposePar
+
+runParallel $(getApplication)
+
+#------------------------------------------------------------------------------
diff --git a/tutorials/lagrangian/reactingParcelFoam/verticalChannel/system/controlDict b/tutorials/lagrangian/reactingParcelFoam/verticalChannel/system/controlDict
index 9ac52b9bae..3442d41a3a 100644
--- a/tutorials/lagrangian/reactingParcelFoam/verticalChannel/system/controlDict
+++ b/tutorials/lagrangian/reactingParcelFoam/verticalChannel/system/controlDict
@@ -53,7 +53,27 @@ maxDeltaT       1e-03;
 
 functions
 {
-    surfaceFieldValue1
+    avgInlets
+    {
+        type            surfaceFieldValue;
+        libs            (fieldFunctionObjects);
+        enabled         yes;
+        writeControl    writeTime;
+        log             yes;
+        writeFields     off;
+        surfaceFormat   vtk;
+        regionType      patch;
+        // name            inlet;
+        names           (inlet "inlet.*");
+        operation       weightedAverage;
+        weightField     phi;
+        fields
+        (
+            T
+        );
+    }
+
+    avgOutlets
     {
         type            surfaceFieldValue;
         libs            (fieldFunctionObjects);
diff --git a/tutorials/lagrangian/reactingParcelFoam/verticalChannel/system/decomposeParDict b/tutorials/lagrangian/reactingParcelFoam/verticalChannel/system/decomposeParDict
new file mode 100644
index 0000000000..e12fff0f4b
--- /dev/null
+++ b/tutorials/lagrangian/reactingParcelFoam/verticalChannel/system/decomposeParDict
@@ -0,0 +1,30 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2006                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains 4;
+
+method          scotch;
+
+coeffs
+{
+    n           (2 2 1);
+}
+
+distributed     no;
+
+roots           ( );
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/system/controlDict b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/system/controlDict
index b12d88a3d2..ae2bb7accb 100644
--- a/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/system/controlDict
+++ b/tutorials/lagrangian/simpleReactingParcelFoam/verticalChannel/system/controlDict
@@ -48,7 +48,7 @@ runTimeModifiable yes;
 
 functions
 {
-    surfaceFieldValue1
+    avgOutlets
     {
         type            surfaceFieldValue;
         libs            (fieldFunctionObjects);