diff --git a/src/faOptions/corrections/limitVelocity/limitVelocity.H b/src/faOptions/corrections/limitVelocity/limitVelocity.H
index e5da8a400a..c4abb1c8d5 100644
--- a/src/faOptions/corrections/limitVelocity/limitVelocity.H
+++ b/src/faOptions/corrections/limitVelocity/limitVelocity.H
@@ -65,8 +65,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef limitVelocity_H
-#define limitVelocity_H
+#ifndef Foam_fa_limitVelocity_H
+#define Foam_fa_limitVelocity_H
 
 #include "faceSetOption.H"
 
@@ -83,7 +83,7 @@ namespace fa
 
 class limitVelocity
 :
-    public faceSetOption
+    public fa::faceSetOption
 {
 protected:
 
@@ -96,17 +96,6 @@ protected:
         scalar max_;
 
 
-private:
-
-    // Private Member Functions
-
-        //- No copy construct
-        limitVelocity(const limitVelocity&) = delete;
-
-        //- No copy assignment
-        void operator=(const limitVelocity&) = delete;
-
-
 public:
 
     //- Runtime type information
@@ -124,6 +113,12 @@ public:
             const fvPatch& patch
         );
 
+        //- No copy construct
+        limitVelocity(const limitVelocity&) = delete;
+
+        //- No copy assignment
+        void operator=(const limitVelocity&) = delete;
+
 
     //- Destructor
     virtual ~limitVelocity() = default;
diff --git a/src/faOptions/faOption/faOptions.C b/src/faOptions/faOption/faOptions.C
index a97e165597..7c2401be57 100644
--- a/src/faOptions/faOption/faOptions.C
+++ b/src/faOptions/faOption/faOptions.C
@@ -103,26 +103,19 @@ Foam::fa::options& Foam::fa::options::New(const fvPatch& p)
 {
     const fvMesh& mesh = p.boundaryMesh().mesh();
 
-    if (mesh.thisDb().foundObject<options>(typeName))
-    {
-        return const_cast<options&>
-        (
-            mesh.lookupObject<options>(typeName)
-        );
-    }
-    else
-    {
-        if (debug)
-        {
-            InfoInFunction
-                << "Constructing " << typeName
-                << " for region " << mesh.name() << endl;
-        }
+    options* ptr = mesh.thisDb().getObjectPtr<options>(typeName);
 
-        options* objectPtr = new options(p);
-        regIOobject::store(objectPtr);
-        return *objectPtr;
+    if (!ptr)
+    {
+        DebugInFunction
+            << "Constructing " << typeName
+            << " for region " << mesh.name() << endl;
+
+        ptr = new options(p);
+        regIOobject::store(ptr);
     }
+
+    return *ptr;
 }
 
 
@@ -133,10 +126,8 @@ bool Foam::fa::options::read()
         optionList::read(*this);
         return true;
     }
-    else
-    {
-        return false;
-    }
+
+    return false;
 }
 
 
diff --git a/src/faOptions/faOption/faOptions.H b/src/faOptions/faOption/faOptions.H
index a80969885c..165c340164 100644
--- a/src/faOptions/faOption/faOptions.H
+++ b/src/faOptions/faOption/faOptions.H
@@ -86,8 +86,7 @@ public:
 
 
     //- Destructor
-    virtual ~options()
-    {}
+    virtual ~options() = default;
 
 
     // Member Functions
diff --git a/src/faOptions/faceSetOption/faceSetOption.C b/src/faOptions/faceSetOption/faceSetOption.C
index aa33a78507..c28d589ed5 100644
--- a/src/faOptions/faceSetOption/faceSetOption.C
+++ b/src/faOptions/faceSetOption/faceSetOption.C
@@ -26,6 +26,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "faceSetOption.H"
+#include "faceSet.H"
 #include "areaFields.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -46,7 +47,10 @@ const Foam::Enum
 Foam::fa::faceSetOption::selectionModeTypeNames_
 ({
     { selectionModeType::smAll, "all" },
-    { selectionModeType::smVolFaceZone, "volFaceZone" }
+    { selectionModeType::smFaceSet, "faceSet" },
+    { selectionModeType::smFaceZone, "faceZone" },
+    { selectionModeType::smPatch, "patch" },
+    { selectionModeType::smFaceZone, "volFaceZone" }  // Compat?
 });
 
 
@@ -54,15 +58,44 @@ Foam::fa::faceSetOption::selectionModeTypeNames_
 
 void Foam::fa::faceSetOption::setSelection(const dictionary& dict)
 {
+    selectionNames_.clear();
+
     switch (selectionMode_)
     {
         case smAll:
         {
             break;
         }
-        case smVolFaceZone:
+        case smFaceSet:
         {
-            dict.readEntry("faceZone", zoneName_);
+            selectionNames_.resize(1);
+            dict.readEntry("faceSet", selectionNames_.first());
+            break;
+        }
+        case smFaceZone:
+        {
+            if
+            (
+                !dict.readIfPresent("faceZones", selectionNames_)
+             || selectionNames_.empty()
+            )
+            {
+                selectionNames_.resize(1);
+                dict.readEntry("faceZone", selectionNames_.first());
+            }
+            break;
+        }
+        case smPatch:
+        {
+            if
+            (
+                !dict.readIfPresent("patches", selectionNames_)
+             || selectionNames_.empty()
+            )
+            {
+                selectionNames_.resize(1);
+                dict.readEntry("patch", selectionNames_.first());
+            }
             break;
         }
         default:
@@ -82,21 +115,22 @@ void Foam::fa::faceSetOption::setArea()
 {
     // Set area information
 
-    scalar sumArea = 0.0;
+    scalar sumArea = 0;
     for (const label facei : faces_)
     {
         sumArea += regionMesh().S()[facei];
     }
     reduce(sumArea, sumOp<scalar>());
 
-    const scalar AOld = A_;
+    const scalar old(A_);
     A_ = sumArea;
 
-    // Convert both areas to representation using current writeprecision
-    word AOldName(Time::timeName(AOld, IOstream::defaultPrecision()));
-    word AName(Time::timeName(A_, IOstream::defaultPrecision()));
-
-    if (AName != AOldName)
+    // Compare area values, stringified using current write precision
+    if
+    (
+        Time::timeName(old, IOstream::defaultPrecision())
+     != Time::timeName(A_, IOstream::defaultPrecision())
+    )
     {
         Info<< indent
             << "- selected " << returnReduce(faces_.size(), sumOp<label>())
@@ -109,37 +143,32 @@ void Foam::fa::faceSetOption::setFaceSelection()
 {
     switch (selectionMode_)
     {
-        case smVolFaceZone:
+        case smAll:
+        {
+            Info<< indent << "- selecting all faces" << endl;
+            faces_ = identity(regionMesh().nFaces());
+
+            break;
+        }
+
+        case smFaceSet:
         {
             Info<< indent
-                << "- selecting faces using volume-mesh faceZone "
-                << zoneName_ << nl;
+                << "- selecting face subset using volume-mesh faceSet "
+                << zoneName() << nl;
 
-            // Also handles groups, multiple zones (as wordRe match) ...
-            labelList zoneIDs = mesh_.faceZones().indices(zoneName_);
-
-            if (zoneIDs.empty())
-            {
-                FatalErrorInFunction
-                    << "No matching faceZones: " << zoneName_ << nl
-                    << "Valid zones : "
-                    << flatOutput(mesh_.faceZones().names()) << nl
-                    << "Valid groups: "
-                    << flatOutput(mesh_.faceZones().groupNames())
-                    << nl
-                    << exit(FatalError);
-            }
-
-            const bitSet isZoneFace(mesh_.faceZones().selection(zoneIDs));
+            const faceSet subset(mesh_, zoneName());
 
             const labelUList& faceLabels = regionMesh().faceLabels();
 
             faces_.resize_nocopy(faceLabels.size());
 
             label nUsed = 0;
-            for (const label facei : faceLabels)
+            forAll(faceLabels, facei)
             {
-                if (isZoneFace[facei])
+                const label meshFacei = faceLabels[facei];
+
+                if (subset.test(meshFacei))
                 {
                     faces_[nUsed] = facei;
                     ++nUsed;
@@ -149,13 +178,93 @@ void Foam::fa::faceSetOption::setFaceSelection()
             break;
         }
 
-        case smAll:
+        case smFaceZone:
         {
-            Info<< indent << "- selecting all faces" << endl;
-            faces_ = identity(regionMesh().nFaces());
+            Info<< indent
+                << "- selecting face subset using volume-mesh faceZones "
+                << flatOutput(selectionNames_) << nl;
 
+            const auto& zones = mesh_.faceZones();
+
+            // Also handles groups, multiple zones etc ...
+            labelList zoneIDs = zones.indices(selectionNames_);
+
+            if (zoneIDs.empty())
+            {
+                FatalErrorInFunction
+                    << "No matching faceZones: "
+                    << flatOutput(selectionNames_) << nl
+                    << "Valid zones : "
+                    << flatOutput(zones.names()) << nl
+                    << "Valid groups: "
+                    << flatOutput(zones.groupNames())
+                    << nl
+                    << exit(FatalError);
+            }
+
+            const bitSet subset(mesh_.faceZones().selection(zoneIDs));
+
+            const labelUList& faceLabels = regionMesh().faceLabels();
+
+            faces_.resize_nocopy(faceLabels.size());
+
+            label nUsed = 0;
+            forAll(faceLabels, facei)
+            {
+                const label meshFacei = faceLabels[facei];
+
+                if (subset.test(meshFacei))
+                {
+                    faces_[nUsed] = facei;
+                    ++nUsed;
+                }
+            }
+            faces_.resize(nUsed);
             break;
         }
+
+        case smPatch:
+        {
+            Info<< indent
+                << "- selecting face subset using volume-mesh patches "
+                << flatOutput(selectionNames_) << nl;
+
+            const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
+
+            // Also handles groups, multiple patches etc ...
+            labelList patchIDs = pbm.indices(selectionNames_);
+
+            if (patchIDs.empty())
+            {
+                FatalErrorInFunction
+                    << "No matching patches: "
+                    << flatOutput(selectionNames_) << nl
+                    << "Valid patches : "
+                    << flatOutput(pbm.names()) << nl
+                    << "Valid groups: "
+                    << flatOutput(pbm.groupNames()) << nl
+                    << exit(FatalError);
+            }
+
+            const List<labelPair>& patchFaces = regionMesh().whichPatchFaces();
+
+            faces_.resize_nocopy(patchFaces.size());
+
+            label nUsed = 0;
+            forAll(patchFaces, facei)
+            {
+                const label patchi = patchFaces[facei].first();
+
+                if (patchIDs.found(patchi))
+                {
+                    faces_[nUsed] = facei;
+                    ++nUsed;
+                }
+            }
+            faces_.resize(nUsed);
+            break;
+        }
+
         default:
         {
             FatalErrorInFunction
@@ -166,6 +275,12 @@ void Foam::fa::faceSetOption::setFaceSelection()
                 << exit(FatalError);
         }
     }
+
+    if (smAll != selectionMode_ && returnReduceAnd(faces_.empty()))
+    {
+        WarningInFunction
+            << "No faces selected!" << endl;
+    }
 }
 
 
@@ -183,7 +298,7 @@ Foam::fa::faceSetOption::faceSetOption
     timeStart_(-1),
     duration_(0),
     selectionMode_(selectionModeTypeNames_.get("selectionMode", coeffs_)),
-    zoneName_(),
+    selectionNames_(),
     A_(0)
 {
     if (isActive())
@@ -229,6 +344,8 @@ bool Foam::fa::faceSetOption::read(const dictionary& dict)
 {
     if (fa::option::read(dict))
     {
+        timeStart_ = -1;
+
         if (coeffs_.readIfPresent("timeStart", timeStart_))
         {
             coeffs_.readEntry("duration", duration_);
diff --git a/src/faOptions/faceSetOption/faceSetOption.H b/src/faOptions/faceSetOption/faceSetOption.H
index c67407270c..16416dc4fe 100644
--- a/src/faOptions/faceSetOption/faceSetOption.H
+++ b/src/faOptions/faceSetOption/faceSetOption.H
@@ -39,18 +39,23 @@ Usage
         ...
 
         // Mandatory entries (unmodifiable)
-        selectionMode     all;
+        selectionMode   all;
 
         // Optional entries (runtime modifiable)
-        timeStart         1.0;
+        timeStart       1.0;
 
         // Conditional mandatory entries (runtime modifiable)
 
-            // when timeStart entry is present
-            duration          1.4;
+        // when timeStart entry is present
+        duration        1.4;
 
-            // when selectionMode=volFaceZone
-            faceZone          <faceZoneName>;
+        // when selectionMode=faceZone
+        faceZones       (<name> ...);
+        //or: faceZone  <name>;
+
+        // when selectionMode=patch
+        patches         (<name> ...)
+        //or: patch     <name>;
     }
     \endverbatim
 
@@ -61,13 +66,21 @@ Usage
       timeStart     | Start time of faOption             | scalar | no  | -1
       duration      | Duration of faOption execution <!--
                           --> starting from timeStart    | scalar | cndtnl | 0
-      faceZone      | Name of operand faceZone           | word  | cndtnl | -
+      faceSet       | Name of operand faceSet            | word   | cndtnl | -
+      faceZone      | Name of operand faceZone(s)        | wordRe | cndtnl | -
+      faceZones     | Names of operand faceZones         | wordRes| cndtnl | -
+      patch         | Name of operand poly patch(s)      | wordRe | cndtnl | -
+      patches       | Names of operand poly patches      | wordRes| cndtnl | -
     \endtable
 
     Options for the \c selectionMode entry:
     \verbatim
       all       | Use all faces in the computational domain
-      faceZone  | Use a given faceZone
+      faceSet   | Use subset corresponding to specified (volume) faceSet
+      faceZones | Use subset corresponding to specified (volume) faceZones
+      faceZone  | Use subset corresponding to specified (volume) faceZones
+      patch     | Use subset corresponding of specified volume patches
+      patches   | Use subset corresponding of specified volume patches
     \endverbatim
 
     The inherited entries are elaborated in:
@@ -78,6 +91,7 @@ Note
 
 SourceFiles
     faceSetOption.C
+    faceSetOptionI.H
 
 \*---------------------------------------------------------------------------*/
 
@@ -111,8 +125,10 @@ public:
         //- Enumeration for selection mode types
         enum selectionModeType
         {
-            smAll,
-            smVolFaceZone
+            smAll,          //!< "all" finite-area faces
+            smFaceSet,      //!< "faceSet" : subset with (volume) face set
+            smFaceZone,     //!< "faceZone" : subset with (volume) zone faces
+            smPatch,        //!< "patch" : subset with (volume) patches
         };
 
         //- List of selection mode type names
@@ -132,8 +148,8 @@ protected:
         //- Face selection mode
         selectionModeType selectionMode_;
 
-        //- Name of zone for (volume) "faceZone" selection
-        wordRe zoneName_;
+        //- Face selection names (for set, zone or patch selections)
+        wordRes selectionNames_;
 
         //- Set of faces to apply source to
         labelList faces_;
@@ -153,6 +169,10 @@ protected:
         //- Recalculate the area
         void setArea();
 
+        //- Zero all non-selected locations within field
+        template<class Type>
+        inline void subsetFilter(List<Type>& field) const;
+
 
 public:
 
@@ -189,20 +209,23 @@ public:
             //- Return true if within time limits
             inline bool inTimeLimits(const scalar timeValue) const;
 
-            //- Return const access to the face selection mode
-            inline selectionModeType selectionMode() const noexcept;
-
             //- True if sub-selection should be used
             inline bool useSubMesh() const noexcept;
 
-            //- Return const access to the name of (volume) face zone
-            //- for "faceZone" selection mode
-            const wordRe& zoneName() const noexcept { return zoneName_; }
+            //- Return the face selection mode
+            inline selectionModeType selectionMode() const noexcept;
+
+            //- Return const access to the selection names
+            //- (set, zone or patch selection)
+            inline const wordRes& selectionNames() const noexcept;
+
+            //- Return const access to the first set/zone/patch name
+            inline const wordRe& zoneName() const;
 
             //- Return const access to the total face area
             inline scalar A() const noexcept;
 
-            //- Return const access to the face selection
+            //- Return const access to the local finite-area face selection
             inline const labelList& faces() const noexcept;
 
 
diff --git a/src/faOptions/faceSetOption/faceSetOptionI.H b/src/faOptions/faceSetOption/faceSetOptionI.H
index 76afe7b5d5..40cde9f860 100644
--- a/src/faOptions/faceSetOption/faceSetOptionI.H
+++ b/src/faOptions/faceSetOption/faceSetOptionI.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -25,6 +25,25 @@ License
 
 \*---------------------------------------------------------------------------*/
 
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+inline void Foam::fa::faceSetOption::subsetFilter(List<Type>& field) const
+{
+    if (selectionMode_ != selectionModeType::smAll)
+    {
+        List<Type> filtered(field.size(), Zero);
+
+        for (const label facei : faces_)
+        {
+            filtered[facei] = field[facei];
+        }
+
+        field.swap(filtered);
+    }
+}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 inline Foam::scalar Foam::fa::faceSetOption::timeStart() const noexcept
@@ -60,6 +79,19 @@ Foam::fa::faceSetOption::selectionMode() const noexcept
 }
 
 
+inline const Foam::wordRes&
+Foam::fa::faceSetOption::selectionNames() const noexcept
+{
+    return selectionNames_;
+}
+
+
+inline const Foam::wordRe& Foam::fa::faceSetOption::zoneName() const
+{
+    return (selectionNames_.empty() ? wordRe::null : selectionNames_.first());
+}
+
+
 inline bool Foam::fa::faceSetOption::useSubMesh() const noexcept
 {
     return selectionMode_ != selectionModeType::smAll;
diff --git a/src/faOptions/sources/derived/contactHeatFluxSource/contactHeatFluxSource.C b/src/faOptions/sources/derived/contactHeatFluxSource/contactHeatFluxSource.C
index 5865449c5a..0e4f52da14 100644
--- a/src/faOptions/sources/derived/contactHeatFluxSource/contactHeatFluxSource.C
+++ b/src/faOptions/sources/derived/contactHeatFluxSource/contactHeatFluxSource.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -59,22 +59,8 @@ Foam::fa::contactHeatFluxSource::contactHeatFluxSource
     TName_(dict.getOrDefault<word>("T", "T")),
     TprimaryName_(dict.get<word>("Tprimary")),
     Tp_(mesh().lookupObject<volScalarField>(TprimaryName_)),
-    Tw1_
-    (
-        IOobject
-        (
-            "Tw1_" + sourceName,
-            mesh().time().timeName(),
-            mesh(),
-            IOobject::READ_IF_PRESENT,
-            IOobject::NO_WRITE
-        ),
-        regionMesh(),
-        dimensionedScalar(dimTemperature, Zero),
-        zeroGradientFaPatchScalarField::typeName
-    ),
-    thicknessLayers_(Zero),
-    kappaLayers_(Zero),
+    thicknessLayers_(),
+    kappaLayers_(),
     contactRes_(0),
     curTimeIndex_(-1)
 {
@@ -88,57 +74,31 @@ Foam::fa::contactHeatFluxSource::contactHeatFluxSource
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::tmp<Foam::areaScalarField> Foam::fa::contactHeatFluxSource::htc() const
+Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::areaMesh>>
+Foam::fa::contactHeatFluxSource::htc() const
 {
-    IOobject io
+    auto thtc = DimensionedField<scalar, areaMesh>::New
     (
-        "thtc",
-        mesh().time().timeName(),
-        mesh(),
-        IOobject::NO_READ,
-        IOobject::NO_WRITE
+        "htc_" + option::name(),
+        regionMesh(),
+        dimensionedScalar(dimPower/dimArea/dimTemperature, Zero)
     );
+    auto& htc = thtc.ref();
 
-    tmp<areaScalarField> thtc
-    (
-        new areaScalarField
-        (
-            io,
-            regionMesh(),
-            dimensionedScalar(dimPower/dimArea/dimTemperature, Zero)
-        )
-    );
-
-    areaScalarField& htc = thtc.ref();
-
-    const volScalarField::Boundary& vfb = Tp_.boundaryField();
-
-    htc.primitiveFieldRef() =
+    htc.field() =
         temperatureCoupledBase::kappa
         (
-            vsm().mapInternalToSurface<scalar>(vfb)()
+            vsm().mapInternalToSurface<scalar>(Tp_)()
         )*patch().deltaCoeffs();
 
     if (contactRes_ != 0)
     {
-        tmp<areaScalarField> tcontact
-        (
-            new areaScalarField
-            (
-                io,
-                regionMesh(),
-                dimensionedScalar
-                (
-                    "contact",
-                    dimPower/dimArea/dimTemperature,
-                    contactRes_
-                )
-            )
-        );
-        areaScalarField& contact = tcontact.ref();
-        htc.primitiveFieldRef() += contact.primitiveField();
+        htc.field() += contactRes_;
     }
 
+    // Zero htc for non-mapped faces
+    faceSetOption::subsetFilter(htc.field());
+
     return thtc;
 }
 
@@ -153,19 +113,26 @@ void Foam::fa::contactHeatFluxSource::addSup
 {
     if (isActive())
     {
-        DebugInfo<< name() << ": applying source to " << eqn.psi().name()
-                 << endl;
+        DebugInfo
+            << name() << ": applying source to "
+            << eqn.psi().name() << endl;
 
         if (curTimeIndex_ != mesh().time().timeIndex())
         {
-            const volScalarField::Boundary& vfb = Tp_.boundaryField();
+            tmp<DimensionedField<scalar, areaMesh>> htcw(htc());
 
-            Tw1_.primitiveFieldRef() =
-                this->vsm().mapInternalToSurface<scalar>(vfb);
+            // Wall temperature - mapped from primary field to finite-area
+            auto Twall = DimensionedField<scalar, areaMesh>::New
+            (
+                "Tw_" + option::name(),
+                regionMesh(),
+                dimensionedScalar(dimTemperature, Zero)
+            );
 
-            tmp<areaScalarField> htcw = htc();
+            Twall.ref().field() =
+                this->vsm().mapInternalToSurface<scalar>(Tp_);
 
-            eqn += -fam::Sp(htcw(), eqn.psi()) + htcw()*Tw1_;
+            eqn += -fam::Sp(htcw(), eqn.psi()) + htcw()*Twall;
 
             curTimeIndex_ = mesh().time().timeIndex();
         }
@@ -179,17 +146,21 @@ bool Foam::fa::contactHeatFluxSource::read(const dictionary& dict)
     {
         coeffs_.readIfPresent("T", TName_);
 
+        contactRes_ = 0;
+
         if (dict.readIfPresent("thicknessLayers", thicknessLayers_))
         {
             dict.readEntry("kappaLayers", kappaLayers_);
 
-            if (thicknessLayers_.size() > 0)
+            // Calculate effective thermal resistance by harmonic averaging
+
+            forAll(thicknessLayers_, iLayer)
+            {
+                contactRes_ += thicknessLayers_[iLayer]/kappaLayers_[iLayer];
+            }
+
+            if (thicknessLayers_.size())
             {
-                // Calculate effective thermal resistance by harmonic averaging
-                forAll(thicknessLayers_, iLayer)
-                {
-                    contactRes_ += thicknessLayers_[iLayer]/kappaLayers_[iLayer];
-                }
                 contactRes_ = scalar(1)/contactRes_;
             }
         }
diff --git a/src/faOptions/sources/derived/contactHeatFluxSource/contactHeatFluxSource.H b/src/faOptions/sources/derived/contactHeatFluxSource/contactHeatFluxSource.H
index fbc4c6c0a0..3ceada68bf 100644
--- a/src/faOptions/sources/derived/contactHeatFluxSource/contactHeatFluxSource.H
+++ b/src/faOptions/sources/derived/contactHeatFluxSource/contactHeatFluxSource.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -76,8 +76,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef fa_contactHeatFluxSource_H
-#define fa_contactHeatFluxSource_H
+#ifndef Foam_fa_contactHeatFluxSource_H
+#define Foam_fa_contactHeatFluxSource_H
 
 #include "faOption.H"
 #include "Function1.H"
@@ -112,9 +112,6 @@ class contactHeatFluxSource
         //- Primary region temperature
         const volScalarField& Tp_;
 
-        //- Temperature - wall [K]
-        areaScalarField Tw1_;
-
         //- Thickness of layers
         scalarList thicknessLayers_;
 
@@ -130,8 +127,8 @@ class contactHeatFluxSource
 
     // Private Member Functions
 
-        //- Return htc from the primary region
-        tmp<areaScalarField> htc() const;
+        //- Return htc [W/m2/K] coupling to the primary region
+        tmp<DimensionedField<scalar, areaMesh>> htc() const;
 
 
 public:
diff --git a/src/faOptions/sources/derived/externalFileSource/externalFileSource.C b/src/faOptions/sources/derived/externalFileSource/externalFileSource.C
index ed557286e6..e885726d93 100644
--- a/src/faOptions/sources/derived/externalFileSource/externalFileSource.C
+++ b/src/faOptions/sources/derived/externalFileSource/externalFileSource.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -64,13 +64,14 @@ Foam::fa::externalFileSource::externalFileSource
             mesh_.time().timeName(),
             mesh_,
             IOobject::NO_READ,
-            IOobject::NO_WRITE
+            IOobject::NO_WRITE,
+            false  // Do not register
         ),
         regionMesh(),
-        dimensionedScalar("pExt", dimPressure, Zero),
-        zeroGradientFaPatchScalarField::typeName
+        dimensionedScalar(dimPressure, Zero)
     ),
-    value_
+    curTimeIndex_(-1),
+    mapping_
     (
         new PatchFunction1Types::MappedFile<scalar>
         (
@@ -80,8 +81,7 @@ Foam::fa::externalFileSource::externalFileSource
             tableName_,          // field table name
             true                 // face values
         )
-    ),
-    curTimeIndex_(-1)
+    )
 {
     fieldNames_.resize(1, fieldName_);
 
@@ -93,6 +93,17 @@ Foam::fa::externalFileSource::externalFileSource
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+void Foam::fa::externalFileSource::updateMapping()
+{
+    const scalar t = mesh().time().value();
+
+    pExt_.field() = mapping_->value(t);
+
+    // Zero pressure for non-mapped faces
+    faceSetOption::subsetFilter(pExt_.field());
+}
+
+
 void Foam::fa::externalFileSource::addSup
 (
     const areaScalarField& solidMass,
@@ -100,16 +111,18 @@ void Foam::fa::externalFileSource::addSup
     const label fieldi
 )
 {
-    const scalar t = mesh().time().value();
-
-    if (isActive() && t > timeStart() && t < (timeStart() + duration()))
+    if (isActive())
     {
-        DebugInfo<< name() << ": applying source to " << eqn.psi().name()<<endl;
+        DebugInfo
+            << name() << ": applying source to "
+            << eqn.psi().name() << endl;
 
         if (curTimeIndex_ != mesh().time().timeIndex())
         {
-            pExt_.field() = value_->value(t);
-            eqn += pExt_/solidMass;
+            updateMapping();
+
+            eqn += pExt_/solidMass.internalField();
+
             curTimeIndex_ = mesh().time().timeIndex();
         }
     }
diff --git a/src/faOptions/sources/derived/externalFileSource/externalFileSource.H b/src/faOptions/sources/derived/externalFileSource/externalFileSource.H
index 974f23925a..c5c2447106 100644
--- a/src/faOptions/sources/derived/externalFileSource/externalFileSource.H
+++ b/src/faOptions/sources/derived/externalFileSource/externalFileSource.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -31,7 +31,7 @@ Group
 
 Description
     Applies sources on a specified field within a specified region
-    by using an external table file for compressible flows.
+    by using an external table file.
 
 Usage
     Minimal example by using \c constant/faOptions:
@@ -68,8 +68,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef fa_externalFileSource_H
-#define fa_externalFileSource_H
+#ifndef Foam_fa_externalFileSource_H
+#define Foam_fa_externalFileSource_H
 
 #include "faOption.H"
 #include "areaFields.H"
@@ -100,14 +100,20 @@ class externalFileSource
         word tableName_;
 
         //- External pressure field
-        areaScalarField pExt_;
-
-        //- Mapped data from file
-        autoPtr<PatchFunction1Types::MappedFile<scalar>> value_;
+        DimensionedField<scalar, areaMesh> pExt_;
 
         //- Current time index (used for updating)
         label curTimeIndex_;
 
+        //- Mapped data from file
+        autoPtr<PatchFunction1Types::MappedFile<scalar>> mapping_;
+
+
+    // Private Member Functions
+
+        //- Update the pExt_ mapping
+        void updateMapping();
+
 
 public:
 
diff --git a/src/faOptions/sources/derived/externalHeatFluxSource/externalHeatFluxSource.C b/src/faOptions/sources/derived/externalHeatFluxSource/externalHeatFluxSource.C
index e29fa73367..c8cb3d0c19 100644
--- a/src/faOptions/sources/derived/externalHeatFluxSource/externalHeatFluxSource.C
+++ b/src/faOptions/sources/derived/externalHeatFluxSource/externalHeatFluxSource.C
@@ -71,9 +71,9 @@ Foam::fa::externalHeatFluxSource::externalHeatFluxSource
     fa::faceSetOption(sourceName, modelType, dict, patch),
     mode_(operationModeNames.get("mode", dict)),
     TName_(dict.getOrDefault<word>("T", "T")),
-    Q_(0),
-    q_(0),
-    h_(0),
+    Q_(nullptr),
+    q_(nullptr),
+    h_(nullptr),
     Ta_(nullptr),
     emissivity_(dict.getOrDefault<scalar>("emissivity", 0))
 {
@@ -97,73 +97,98 @@ void Foam::fa::externalHeatFluxSource::addSup
 {
     if (isActive())
     {
-        DebugInfo<< name() << ": applying source to "
+        DebugInfo
+            << name() << ": applying source to "
             << eqn.psi().name() << endl;
 
-        IOobject io
-        (
-            "Q",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE,
-            false
-        );
+        scalar qflux = 0;
 
-        auto tQ = tmp<areaScalarField>::New
-        (
-            io,
-            regionMesh(),
-            dimensionedScalar("q", dimPower/sqr(dimLength), 0),
-            zeroGradientFaPatchScalarField::typeName
-        );
-        areaScalarField& Q = tQ.ref();
+        const scalar timeVal = mesh_.time().timeOutputValue();
 
         switch (mode_)
         {
             case fixedPower:
             {
-                Q.primitiveFieldRef() = Q_/regionMesh().S().field();
+                // From [W] to [W/m2]
+                qflux = Q_->value(timeVal)/(faceSetOption::A() + VSMALL);
+                break;
+            }
+
+            case fixedHeatFlux:
+            {
+                qflux = q_->value(timeVal);
+                break;
+            }
+
+            default:
+            {
+                break;
+            }
+        }
+
+        switch (mode_)
+        {
+            case fixedPower:
+            case fixedHeatFlux:
+            {
+                auto tQ = DimensionedField<scalar, areaMesh>::New
+                (
+                    "Q",
+                    regionMesh(),
+                    dimensionedScalar(dimPower/sqr(dimLength), Zero)
+                );
+                auto& Q = tQ.ref();
+
+                if (faceSetOption::useSubMesh())
+                {
+                    UIndirectList<scalar>(Q.field(), faceSetOption::faces())
+                        = qflux;
+                }
+                else
+                {
+                    Q.field() = qflux;
+                }
+
                 eqn += Q;
 
                 break;
             }
-            case fixedHeatFlux:
-            {
-                Q.primitiveFieldRef() = q_;
-                eqn += Q;
-                break;
-            }
+
             case fixedHeatTransferCoeff:
             {
                 const dimensionedScalar Ta
                 (
                     "Ta",
                     dimTemperature,
-                    Ta_->value(mesh_.time().timeOutputValue())
+                    Ta_->value(timeVal)
                 );
 
-                areaScalarField hp
+                auto thp = DimensionedField<scalar, areaMesh>::New
                 (
-                    io,
+                    "h",
                     regionMesh(),
                     dimensionedScalar
                     (
                         "h",
                         dimPower/sqr(dimLength)/dimTemperature,
-                        h_
+                        h_->value(timeVal)
                     )
                 );
+                auto& hp = thp.ref();
 
-                const areaScalarField hpTa(hp*Ta);
+                DimensionedField<scalar, areaMesh> hpTa(hp*Ta);
 
                 if (emissivity_ > 0)
                 {
                     hp -= emissivity_*sigma.value()*pow3(eqn.psi());
                 }
 
-                eqn -= fam::SuSp(hp, eqn.psi()) - hpTa;
+                // Zero htc for non-mapped faces
+                faceSetOption::subsetFilter(hp.field());
+                faceSetOption::subsetFilter(hpTa.field());
 
+                eqn -= fam::SuSp(hp, eqn.psi()) - hpTa;
+                break;
             }
         }
     }
@@ -183,17 +208,17 @@ bool Foam::fa::externalHeatFluxSource::read(const dictionary& dict)
         {
             case fixedPower:
             {
-                dict.readEntry("Q", Q_);
+                Q_ = Function1<scalar>::New("Q", dict, &mesh_);
                 break;
             }
             case fixedHeatFlux:
             {
-                dict.readEntry("q", q_);
+                Q_ = Function1<scalar>::New("q", dict, &mesh_);
                 break;
             }
             case fixedHeatTransferCoeff:
             {
-                dict.readEntry("h", h_);
+                h_ = Function1<scalar>::New("h", dict, &mesh_);
                 Ta_ = Function1<scalar>::New("Ta", dict, &mesh_);
                 break;
             }
@@ -202,7 +227,7 @@ bool Foam::fa::externalHeatFluxSource::read(const dictionary& dict)
         return true;
     }
 
-     return false;
+    return false;
 }
 
 
diff --git a/src/faOptions/sources/derived/externalHeatFluxSource/externalHeatFluxSource.H b/src/faOptions/sources/derived/externalHeatFluxSource/externalHeatFluxSource.H
index 26e8385276..a7adc670dc 100644
--- a/src/faOptions/sources/derived/externalHeatFluxSource/externalHeatFluxSource.H
+++ b/src/faOptions/sources/derived/externalHeatFluxSource/externalHeatFluxSource.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -89,9 +89,9 @@ Usage
       T         | Name of operand temperature field    | word  | no   | T
       emissivity | Surface emissivity for radiative flux to ambient <!--
                  -->                                   | scalar | no  | 0
-      Q         | Fixed heat power [W]                 | scalar | cndtnl | -
-      q         | Fixed heat flux [W/m2]               | scalar | cndtnl | -
-      h         | Heat transfer coefficient [W/m^2/K]  | scalar | cndtnl | -
+      Q         | Fixed heat power [W]                 | Function1 | cndtnl | -
+      q         | Fixed heat flux [W/m2]               | Function1 | cndtnl | -
+      h         | Heat transfer coefficient [W/m^2/K]  | Function1 | cndtnl | -
       Ta        | Ambient temperature [K]           | Function1 | cndtnl | -
     \endtable
 
@@ -114,8 +114,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef fa_externalHeatFluxSource_H
-#define fa_externalHeatFluxSource_H
+#ifndef Foam_fa_externalHeatFluxSource_H
+#define Foam_fa_externalHeatFluxSource_H
 
 #include "faOption.H"
 #include "Function1.H"
@@ -164,13 +164,13 @@ private:
         word TName_;
 
         //- Heat power [W]
-        scalar Q_;
+        autoPtr<Function1<scalar>> Q_;
 
         //- Heat flux [W/m2]
-        scalar q_;
+        autoPtr<Function1<scalar>> q_;
 
         //- Heat transfer coefficient [W/m2K]
-        scalar h_;
+        autoPtr<Function1<scalar>> h_;
 
         //- Ambient temperature [K]
         autoPtr<Function1<scalar>> Ta_;
@@ -187,7 +187,7 @@ public:
 
     // Constructors
 
-         //- Construct from explicit source name and mesh
+        //- Construct from explicit source name and mesh
         externalHeatFluxSource
         (
             const word& sourceName,
diff --git a/src/faOptions/sources/derived/jouleHeatingSource/jouleHeatingSource.C b/src/faOptions/sources/derived/jouleHeatingSource/jouleHeatingSource.C
index 68a02faa3c..9b946c3566 100644
--- a/src/faOptions/sources/derived/jouleHeatingSource/jouleHeatingSource.C
+++ b/src/faOptions/sources/derived/jouleHeatingSource/jouleHeatingSource.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -59,8 +59,8 @@ Foam::fa::jouleHeatingSource::jouleHeatingSource
         IOobject
         (
             typeName + ":V_" + regionName_,
-            mesh().time().timeName(),
-            mesh(),
+            regionMesh().thisDb().time().timeName(),
+            regionMesh().thisDb(),
             IOobject::MUST_READ,
             IOobject::AUTO_WRITE
         ),
@@ -105,8 +105,9 @@ void Foam::fa::jouleHeatingSource::addSup
 {
     if (isActive())
     {
-        DebugInfo<< name() << ": applying source to " << eqn.psi().name()
-                 << endl;
+        DebugInfo
+            << name() << ": applying source to "
+            << eqn.psi().name() << endl;
 
         if (curTimeIndex_ != mesh().time().timeIndex())
         {
@@ -142,6 +143,14 @@ void Foam::fa::jouleHeatingSource::addSup
         // Add the Joule heating contribution
         areaVectorField gradV("gradV", fac::grad(V_));
 
+        if (debug > 1 && mesh().time().outputTime())
+        {
+            areaScalarField qgradV("gradVSource", (gradV & gradV));
+            qgradV.write();
+        }
+
+        tmp<areaScalarField> tsource;
+
         if (anisotropicElectricalConductivity_)
         {
             const auto& sigma =
@@ -150,7 +159,7 @@ void Foam::fa::jouleHeatingSource::addSup
                     typeName + ":sigma_" + regionName_
                 );
 
-            eqn += (h*sigma & gradV) & gradV;
+            tsource = (h*sigma & gradV) & gradV;
         }
         else
         {
@@ -160,14 +169,13 @@ void Foam::fa::jouleHeatingSource::addSup
                     typeName + ":sigma_" + regionName_
                 );
 
-            eqn += (h*sigma*gradV) & gradV;
-
-            if (mesh().time().outputTime() && debug)
-            {
-                areaScalarField qgradV("gradVSource", (gradV & gradV));
-                qgradV.write();
-            }
+            tsource = (h*sigma*gradV) & gradV;
         }
+
+        // Apply subMesh filter
+        faceSetOption::subsetFilter(tsource.ref().primitiveFieldRef());
+
+        eqn += tsource;
     }
 }
 
diff --git a/src/finiteVolume/cfdTools/general/fvOptions/fvOptions.C b/src/finiteVolume/cfdTools/general/fvOptions/fvOptions.C
index 4361f20a47..a341b20543 100644
--- a/src/finiteVolume/cfdTools/general/fvOptions/fvOptions.C
+++ b/src/finiteVolume/cfdTools/general/fvOptions/fvOptions.C
@@ -102,23 +102,19 @@ Foam::fv::options::options
 
 Foam::fv::options& Foam::fv::options::New(const fvMesh& mesh)
 {
-    if (mesh.thisDb().foundObject<options>(typeName))
-    {
-        return const_cast<options&>
-        (
-            mesh.lookupObject<options>(typeName)
-        );
-    }
-    else
+    options* ptr = mesh.thisDb().getObjectPtr<options>(typeName);
+
+    if (!ptr)
     {
         DebugInFunction
             << "Constructing " << typeName
             << " for region " << mesh.name() << nl;
 
-        options* objectPtr = new options(mesh);
-        regIOobject::store(objectPtr);
-        return *objectPtr;
+        ptr = new options(mesh);
+        regIOobject::store(ptr);
     }
+
+    return *ptr;
 }
 
 
diff --git a/src/fvOptions/cellSetOption/cellSetOption.C b/src/fvOptions/cellSetOption/cellSetOption.C
index 39ff229933..a5844021ef 100644
--- a/src/fvOptions/cellSetOption/cellSetOption.C
+++ b/src/fvOptions/cellSetOption/cellSetOption.C
@@ -48,11 +48,11 @@ const Foam::Enum
 >
 Foam::fv::cellSetOption::selectionModeTypeNames_
 ({
+    { selectionModeType::smAll, "all" },
     { selectionModeType::smGeometric, "geometric" },
     { selectionModeType::smPoints, "points" },
     { selectionModeType::smCellSet, "cellSet" },
     { selectionModeType::smCellZone, "cellZone" },
-    { selectionModeType::smAll, "all" },
 });
 
 
@@ -60,8 +60,14 @@ Foam::fv::cellSetOption::selectionModeTypeNames_
 
 void Foam::fv::cellSetOption::setSelection(const dictionary& dict)
 {
+    selectionNames_.clear();
+
     switch (selectionMode_)
     {
+        case smAll:
+        {
+            break;
+        }
         case smGeometric:
         {
             geometricSelection_ = dict.subDict("selection");
@@ -74,16 +80,21 @@ void Foam::fv::cellSetOption::setSelection(const dictionary& dict)
         }
         case smCellSet:
         {
-            dict.readEntry("cellSet", zoneName_);
+            selectionNames_.resize(1);
+            dict.readEntry("cellSet", selectionNames_.first());
             break;
         }
         case smCellZone:
         {
-            dict.readEntry("cellZone", zoneName_);
-            break;
-        }
-        case smAll:
-        {
+            if
+            (
+                !dict.readIfPresent("cellZones", selectionNames_)
+             || selectionNames_.empty()
+            )
+            {
+                selectionNames_.resize(1);
+                dict.readEntry("cellZone", selectionNames_.first());
+            }
             break;
         }
         default:
@@ -110,14 +121,15 @@ void Foam::fv::cellSetOption::setVol()
     }
     reduce(sumVol, sumOp<scalar>());
 
-    const scalar VOld = V_;
+    const scalar old(V_);
     V_ = sumVol;
 
-    // Convert both volumes to representation using current writeprecision
-    word VOldName(Time::timeName(VOld, IOstream::defaultPrecision()));
-    word VName(Time::timeName(V_, IOstream::defaultPrecision()));
-
-    if (VName != VOldName)
+    // Compare volume values, stringified using current write precision
+    if
+    (
+        Time::timeName(old, IOstream::defaultPrecision())
+     != Time::timeName(V_, IOstream::defaultPrecision())
+    )
     {
         Info<< indent
             << "- selected " << returnReduce(cells_.size(), sumOp<label>())
@@ -130,6 +142,13 @@ void Foam::fv::cellSetOption::setCellSelection()
 {
     switch (selectionMode_)
     {
+        case smAll:
+        {
+            Info<< indent << "- selecting all cells" << endl;
+
+            cells_ = identity(mesh_.nCells());
+            break;
+        }
         case smGeometric:
         {
             Info<< indent << "- selecting cells geometrically" << endl;
@@ -174,48 +193,47 @@ void Foam::fv::cellSetOption::setCellSelection()
         case smCellSet:
         {
             Info<< indent
-                << "- selecting cells using cellSet " << zoneName_ << endl;
+                << "- selecting cells using cellSet "
+                << zoneName() << endl;
 
-            cells_ = cellSet(mesh_, zoneName_).sortedToc();
+            cells_ = cellSet(mesh_, zoneName()).sortedToc();
             break;
         }
         case smCellZone:
         {
             Info<< indent
-                << "- selecting cells using cellZone " << zoneName_ << endl;
+                << "- selecting cells using cellZones "
+                << flatOutput(selectionNames_) << nl;
 
-            // Also handles groups, multiple zones (as wordRe match) ...
-            labelList zoneIDs = mesh_.cellZones().indices(zoneName_);
+            const auto& zones = mesh_.cellZones();
+
+            // Also handles groups, multiple zones etc ...
+            labelList zoneIDs = zones.indices(selectionNames_);
 
             if (zoneIDs.empty())
             {
                 FatalErrorInFunction
-                    << "No matching cellZones: " << zoneName_ << nl
+                    << "No matching cellZones: "
+                    << flatOutput(selectionNames_) << nl
                     << "Valid zones : "
-                    << flatOutput(mesh_.cellZones().names()) << nl
+                    << flatOutput(zones.names()) << nl
                     << "Valid groups: "
-                    << flatOutput(mesh_.cellZones().groupNames())
+                    << flatOutput(zones.groupNames())
                     << nl
                     << exit(FatalError);
             }
 
             if (zoneIDs.size() == 1)
             {
-                cells_ = mesh_.cellZones()[zoneIDs.first()];
+                cells_ = zones[zoneIDs.first()];
+                // TBD: Foam::sort(cells_);
             }
             else
             {
-                cells_ = mesh_.cellZones().selection(zoneIDs).sortedToc();
+                cells_ = zones.selection(zoneIDs).sortedToc();
             }
             break;
         }
-        case smAll:
-        {
-            Info<< indent << "- selecting all cells" << endl;
-
-            cells_ = identity(mesh_.nCells());
-            break;
-        }
         default:
         {
             FatalErrorInFunction
@@ -227,11 +245,7 @@ void Foam::fv::cellSetOption::setCellSelection()
         }
     }
 
-    if
-    (
-        smAll != selectionMode_
-     && returnReduce(cells_.empty(), andOp<bool>())
-    )
+    if (smAll != selectionMode_ && returnReduceAnd(cells_.empty()))
     {
         WarningInFunction
             << "No cells selected!" << endl;
@@ -253,7 +267,7 @@ Foam::fv::cellSetOption::cellSetOption
     timeStart_(-1),
     duration_(0),
     selectionMode_(selectionModeTypeNames_.get("selectionMode", coeffs_)),
-    zoneName_(),
+    selectionNames_(),
     points_(),
     geometricSelection_(),
     V_(0)
@@ -307,6 +321,8 @@ bool Foam::fv::cellSetOption::read(const dictionary& dict)
 {
     if (fv::option::read(dict))
     {
+        timeStart_ = -1;
+
         if (coeffs_.readIfPresent("timeStart", timeStart_))
         {
             coeffs_.readEntry("duration", duration_);
diff --git a/src/fvOptions/cellSetOption/cellSetOption.H b/src/fvOptions/cellSetOption/cellSetOption.H
index 01d943a3b1..32fe7f1fdd 100644
--- a/src/fvOptions/cellSetOption/cellSetOption.H
+++ b/src/fvOptions/cellSetOption/cellSetOption.H
@@ -56,6 +56,8 @@ Usage
             // when selectionMode=cellZone
             cellZone          <name>;
 
+            //OR: cellZones   (<name> ...);
+
             // when selectionMode=points
             points            (<point1> <point2> ... <pointN>);
 
@@ -92,7 +94,8 @@ Usage
       duration      | Duration of fvOption execution <!--
                           --> starting from timeStart    | scalar | cndtnl | 0
       cellSet       | Name of operand cellSet            | word  | cndtnl | -
-      cellZone      | Name of operand cellZone           | word  | cndtnl | -
+      cellZone      | Name of operand cellZone           | wordRe | cndtnl | -
+      cellZones     | Name of operand cellZones          | wordRes  | cndtnl | -
       points        | Set of points in global coordinate <!--
                             --> system              | vectorList | cndtnl | -
       selection     | Dictionary of geometric selections | dict | cndtnl | -
@@ -179,8 +182,8 @@ protected:
         //- Cell selection mode
         selectionModeType selectionMode_;
 
-        //- Name of set/zone for "cellSet" and "cellZone" selectionMode
-        wordRe zoneName_;
+        //- Face selection names (for set or zone selections)
+        wordRes selectionNames_;
 
         //- List of points for "points" selectionMode
         List<point> points_;
@@ -242,15 +245,18 @@ public:
             //- True if within time limits
             inline bool inTimeLimits(const scalar timeValue) const;
 
-            //- Return the cell selection mode
-            inline selectionModeType selectionMode() const noexcept;
-
             //- True if sub-selection should be used
             inline bool useSubMesh() const noexcept;
 
-            //- Return const access to the name of cell set/zone
-            //- for "cellSet" / "cellZone" selection modes
-            const wordRe& zoneName() const noexcept { return zoneName_; }
+            //- Return the cell selection mode
+            inline selectionModeType selectionMode() const noexcept;
+
+            //- Return const access to the selection names
+            //- (set or zone selection)
+            inline const wordRes& selectionNames() const noexcept;
+
+            //- Return const access to the first set/zone name
+            inline const wordRe& zoneName() const;
 
             //- Return const access to the total cell volume
             inline scalar V() const noexcept;
@@ -284,7 +290,7 @@ public:
 
         //- The name of the cell set/zone [as a word]
         //- for "cellSet" / "cellZone" selection modes)
-        const word& cellSetName() const noexcept { return zoneName_; }
+        const word& cellSetName() const { return zoneName(); }
 };
 
 
diff --git a/src/fvOptions/cellSetOption/cellSetOptionI.H b/src/fvOptions/cellSetOption/cellSetOptionI.H
index 1009439479..c246a4c340 100644
--- a/src/fvOptions/cellSetOption/cellSetOptionI.H
+++ b/src/fvOptions/cellSetOption/cellSetOptionI.H
@@ -54,6 +54,12 @@ inline bool Foam::fv::cellSetOption::inTimeLimits(const scalar timeValue) const
 }
 
 
+inline bool Foam::fv::cellSetOption::useSubMesh() const noexcept
+{
+    return selectionMode_ != selectionModeType::smAll;
+}
+
+
 inline Foam::fv::cellSetOption::selectionModeType
 Foam::fv::cellSetOption::selectionMode() const noexcept
 {
@@ -61,9 +67,16 @@ Foam::fv::cellSetOption::selectionMode() const noexcept
 }
 
 
-inline bool Foam::fv::cellSetOption::useSubMesh() const noexcept
+inline const Foam::wordRes&
+Foam::fv::cellSetOption::selectionNames() const noexcept
 {
-    return selectionMode_ != selectionModeType::smAll;
+    return selectionNames_;
+}
+
+
+inline const Foam::wordRe& Foam::fv::cellSetOption::zoneName() const
+{
+    return (selectionNames_.empty() ? wordRe::null : selectionNames_.first());
 }