fvOptions: Changed the source, constraint and correct functions to const

Most fvOptions change the state of the fields and equations they are applied to
but do not change internal state so it makes more sense that the interface is
const, consistent with MeshObjects.  For the few fvOptions which do maintain a
changing state the member data is now mutable.

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H

@@ -50,7 +50,7 @@
     radiationModel& rad = radiation[i];
 
     IOMRFZoneList& MRF = MRFfluid[i];
-    fv::options& fvOptions = fluidFvOptions[i];
+    const fv::options& fvOptions = fluidFvOptions[i];
 
     const dimensionedScalar initialMass
     (

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H

@@ -28,6 +28,6 @@ volScalarField& e = thermo.he();
 
 const volScalarField& betav = betavSolid[i];
 
-fv::options& fvOptions = solidHeatSources[i];
+const fv::options& fvOptions = solidHeatSources[i];
 
 solidNoLoopControl& pimple = pimples.solid(i);

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.C

@@ -48,7 +48,7 @@ namespace Foam
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::fv::VoFSolidificationMeltingSource::update()
+void Foam::fv::VoFSolidificationMeltingSource::update() const
 {
     if (curTimeIndex_ == mesh_.time().timeIndex())
     {
@@ -156,7 +156,7 @@ void Foam::fv::VoFSolidificationMeltingSource::addSup
 (
     fvMatrix<scalar>& eqn,
     const label fieldi
-)
+) const
 {
     apply(geometricOneField(), eqn);
 }
@@ -167,7 +167,7 @@ void Foam::fv::VoFSolidificationMeltingSource::addSup
     const volScalarField& rho,
     fvMatrix<scalar>& eqn,
     const label fieldi
-)
+) const
 {
     apply(rho, eqn);
 }
@@ -177,7 +177,7 @@ void Foam::fv::VoFSolidificationMeltingSource::addSup
 (
     fvMatrix<vector>& eqn,
     const label fieldi
-)
+) const
 {
     if (debug)
     {
@@ -209,7 +209,7 @@ void Foam::fv::VoFSolidificationMeltingSource::addSup
     const volScalarField& rho,
     fvMatrix<vector>& eqn,
     const label fieldi
-)
+) const
 {
     // Momentum source uses a Boussinesq approximation - redirect
     addSup(eqn, fieldi);

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFSolidificationMeltingSource/VoFSolidificationMeltingSource.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2017-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2017-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -117,10 +117,10 @@ class VoFSolidificationMeltingSource
         scalar q_;
 
         //- Solid phase fraction
-        volScalarField alphaSolid_;
+        mutable volScalarField alphaSolid_;
 
         //- Current time index (used for updating)
-        label curTimeIndex_;
+        mutable label curTimeIndex_;
 
 
     // Private Member Functions
@@ -129,11 +129,11 @@ class VoFSolidificationMeltingSource
         word alphaSolidName() const;
 
         //- Update the model
-        void update();
+        void update() const;
 
         //- Helper function to apply to the energy equation
         template<class RhoFieldType>
-        void apply(const RhoFieldType& rho, fvMatrix<scalar>& eqn);
+        void apply(const RhoFieldType& rho, fvMatrix<scalar>& eqn) const;
 
 
 public:
@@ -165,10 +165,18 @@ public:
         // Add explicit and implicit contributions
 
             //- Add explicit contribution to enthalpy equation
-            virtual void addSup(fvMatrix<scalar>& eqn, const label fieldi);
+            virtual void addSup
+            (
+                fvMatrix<scalar>& eqn,
+                const label fieldi
+            ) const;
 
             //- Add implicit contribution to momentum equation
-            virtual void addSup(fvMatrix<vector>& eqn, const label fieldi);
+            virtual void addSup
+            (
+                fvMatrix<vector>& eqn,
+                const label fieldi
+            ) const;
 
 
         // Add explicit and implicit contributions to compressible equation
@@ -179,7 +187,7 @@ public:
                 const volScalarField& rho,
                 fvMatrix<scalar>& eqn,
                 const label fieldi
-            );
+            ) const;
 
             //- Add implicit contribution to compressible momentum equation
             virtual void addSup
@@ -187,7 +195,7 @@ public:
                 const volScalarField& rho,
                 fvMatrix<vector>& eqn,
                 const label fieldi
-            );
+            ) const;
 
 
         // IO

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFSolidificationMeltingSource/VoFSolidificationMeltingSourceTemplates.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2017-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2017-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -33,7 +33,7 @@ void Foam::fv::VoFSolidificationMeltingSource::apply
 (
     const RhoFieldType& rho,
     fvMatrix<scalar>& eqn
-)
+) const
 {
     if (debug)
     {

applications/solvers/multiphase/multiphaseEulerFoam/multiphaseCompressibleMomentumTransportModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C

@@ -477,7 +477,7 @@ void Foam::RASModels::kineticTheoryModel::correct()
         // 'thermal' conductivity (Table 3.3, p. 49)
         kappa_ = conductivityModel_->kappa(alpha, Theta_, gs0_, rho, da, e_);
 
-        fv::options& fvOptions(fv::options::New(mesh_));
+        const fv::options& fvOptions(fv::options::New(mesh_));
 
         // Construct the granular temperature equation (Eq. 3.20, p. 44)
         // NB. note that there are two typos in Eq. 3.20:

applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam/createFieldRefs.H

@@ -1,4 +1,4 @@
 surfaceScalarField& phi = fluid.phi();
 
 const IOMRFZoneList& MRF = fluid.MRF();
-fv::options& fvOptions = fluid.fvOptions();
+const fv::options& fvOptions = fluid.fvOptions();

applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/IATE/IATE.C

@@ -144,7 +144,7 @@ void Foam::diameterModels::IATE::correct()
         R += sources_[j].R(alphaAv, kappai_);
     }
 
-    fv::options& fvOptions(fv::options::New(phase().mesh()));
+    const fv::options& fvOptions(fv::options::New(phase().mesh()));
 
     // Construct the interfacial curvature equation
     fvScalarMatrix kappaiEqn

applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.H

@@ -436,7 +436,7 @@ public:
             inline const IOMRFZoneList& MRF() const;
 
             //- Access the fvOptions
-            inline fv::options& fvOptions() const;
+            inline const fv::options& fvOptions() const;
 
 
         // Sub-model lookup

applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystemI.H

@@ -161,7 +161,7 @@ inline const Foam::IOMRFZoneList& Foam::phaseSystem::MRF() const
 }
 
 
-inline Foam::fv::options& Foam::phaseSystem::fvOptions() const
+inline const Foam::fv::options& Foam::phaseSystem::fvOptions() const
 {
     return fv::options::New(mesh_);
 }

etc/codeTemplates/dynamicCode/codedFvOptionTemplate.C

@@ -131,7 +131,7 @@ ${typeName}FvOption${SourceType}::
 void ${typeName}FvOption${SourceType}::correct
 (
     GeometricField<${TemplateType}, fvPatchField, volMesh>& fld
-)
+) const
 {
     if (${verbose:-false})
     {
@@ -148,7 +148,7 @@ void ${typeName}FvOption${SourceType}::addSup
 (
     fvMatrix<${TemplateType}>& eqn,
     const label fieldi
-)
+) const
 {
     if (${verbose:-false})
     {
@@ -166,7 +166,7 @@ void ${typeName}FvOption${SourceType}::addSup
     const volScalarField& rho,
     fvMatrix<${TemplateType}>& eqn,
     const label fieldi
-)
+) const
 {
     if (${verbose:-false})
     {
@@ -183,7 +183,7 @@ void ${typeName}FvOption${SourceType}::constrain
 (
     fvMatrix<${TemplateType}>& eqn,
     const label fieldi
-)
+) const
 {
     if (${verbose:-false})
     {

etc/codeTemplates/dynamicCode/codedFvOptionTemplate.H

@@ -144,14 +144,14 @@ public:
             virtual void correct
             (
                 GeometricField<${TemplateType}, fvPatchField, volMesh>&
-            );
+            ) const;
 
             //- Explicit and implicit matrix contributions
             virtual void addSup
             (
                 fvMatrix<${TemplateType}>& eqn,
                 const label fieldi
-            );
+            ) const;
 
             //- Explicit and implicit matrix contributions for compressible
             //  equations
@@ -160,14 +160,14 @@ public:
                 const volScalarField& rho,
                 fvMatrix<${TemplateType}>& eqn,
                 const label fieldi
-            );
+            ) const;
 
             //- Set value
             virtual void constrain
             (
                 fvMatrix<${TemplateType}>& eqn,
                 const label fieldi
-            );
+            ) const;
 };
 
 

src/MomentumTransportModels/momentumTransportModels/Base/kOmegaSST/kOmegaSSTBase.C

@@ -429,7 +429,7 @@ void kOmegaSST<MomentumTransportModel, BasicMomentumTransportModel>::correct()
     const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
     const volVectorField& U = this->U_;
     volScalarField& nut = this->nut_;
-    fv::options& fvOptions(fv::options::New(this->mesh_));
+    const fv::options& fvOptions(fv::options::New(this->mesh_));
 
     BasicMomentumTransportModel::correct();
 

src/MomentumTransportModels/momentumTransportModels/LES/DeardorffDiffStress/DeardorffDiffStress.C

@@ -163,7 +163,7 @@ void DeardorffDiffStress<BasicMomentumTransportModel>::correct()
     const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
     const volVectorField& U = this->U_;
     volSymmTensorField& R = this->R_;
-    fv::options& fvOptions(fv::options::New(this->mesh_));
+    const fv::options& fvOptions(fv::options::New(this->mesh_));
 
     ReynoldsStress<LESModel<BasicMomentumTransportModel>>::correct();
 

src/MomentumTransportModels/momentumTransportModels/LES/SpalartAllmarasDES/SpalartAllmarasDES.C

@@ -454,7 +454,7 @@ void SpalartAllmarasDES<BasicMomentumTransportModel>::correct()
     const rhoField& rho = this->rho_;
     const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
     const volVectorField& U = this->U_;
-    fv::options& fvOptions(fv::options::New(this->mesh_));
+    const fv::options& fvOptions(fv::options::New(this->mesh_));
 
     LESeddyViscosity<BasicMomentumTransportModel>::correct();
 

src/MomentumTransportModels/momentumTransportModels/LES/dynamicKEqn/dynamicKEqn.C

@@ -233,7 +233,7 @@ void dynamicKEqn<BasicMomentumTransportModel>::correct()
     const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
     const volVectorField& U = this->U_;
     volScalarField& nut = this->nut_;
-    fv::options& fvOptions(fv::options::New(this->mesh_));
+    const fv::options& fvOptions(fv::options::New(this->mesh_));
 
     LESeddyViscosity<BasicMomentumTransportModel>::correct();
 

src/MomentumTransportModels/momentumTransportModels/LES/dynamicLagrangian/dynamicLagrangian.C

@@ -159,7 +159,7 @@ void dynamicLagrangian<BasicMomentumTransportModel>::correct()
     const rhoField& rho = this->rho_;
     const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
     const volVectorField& U = this->U_;
-    fv::options& fvOptions(fv::options::New(this->mesh_));
+    const fv::options& fvOptions(fv::options::New(this->mesh_));
 
     LESeddyViscosity<BasicMomentumTransportModel>::correct();
 

src/MomentumTransportModels/momentumTransportModels/LES/kEqn/kEqn.C

@@ -159,7 +159,7 @@ void kEqn<BasicMomentumTransportModel>::correct()
     const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
     const volVectorField& U = this->U_;
     volScalarField& nut = this->nut_;
-    fv::options& fvOptions(fv::options::New(this->mesh_));
+    const fv::options& fvOptions(fv::options::New(this->mesh_));
 
     LESeddyViscosity<BasicMomentumTransportModel>::correct();
 

src/MomentumTransportModels/momentumTransportModels/RAS/LRR/LRR.C

@@ -272,7 +272,7 @@ void LRR<BasicMomentumTransportModel>::correct()
     const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
     const volVectorField& U = this->U_;
     volSymmTensorField& R = this->R_;
-    fv::options& fvOptions(fv::options::New(this->mesh_));
+    const fv::options& fvOptions(fv::options::New(this->mesh_));
 
     ReynoldsStress<RASModel<BasicMomentumTransportModel>>::correct();
 

src/MomentumTransportModels/momentumTransportModels/RAS/LaunderSharmaKE/LaunderSharmaKE.C

@@ -246,7 +246,7 @@ void LaunderSharmaKE<BasicMomentumTransportModel>::correct()
     const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
     const volVectorField& U = this->U_;
     volScalarField& nut = this->nut_;
-    fv::options& fvOptions(fv::options::New(this->mesh_));
+    const fv::options& fvOptions(fv::options::New(this->mesh_));
 
     eddyViscosity<RASModel<BasicMomentumTransportModel>>::correct();
 

src/MomentumTransportModels/momentumTransportModels/RAS/RNGkEpsilon/RNGkEpsilon.C

@@ -248,7 +248,7 @@ void RNGkEpsilon<BasicMomentumTransportModel>::correct()
     const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
     const volVectorField& U = this->U_;
     volScalarField& nut = this->nut_;
-    fv::options& fvOptions(fv::options::New(this->mesh_));
+    const fv::options& fvOptions(fv::options::New(this->mesh_));
 
     eddyViscosity<RASModel<BasicMomentumTransportModel>>::correct();
 

src/MomentumTransportModels/momentumTransportModels/RAS/SSG/SSG.C

@@ -282,7 +282,7 @@ void SSG<BasicMomentumTransportModel>::correct()
     const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
     const volVectorField& U = this->U_;
     volSymmTensorField& R = this->R_;
-    fv::options& fvOptions(fv::options::New(this->mesh_));
+    const fv::options& fvOptions(fv::options::New(this->mesh_));
 
     ReynoldsStress<RASModel<BasicMomentumTransportModel>>::correct();
 

src/MomentumTransportModels/momentumTransportModels/RAS/SpalartAllmaras/SpalartAllmaras.C

@@ -355,7 +355,7 @@ void SpalartAllmaras<BasicMomentumTransportModel>::correct()
     const alphaField& alpha = this->alpha_;
     const rhoField& rho = this->rho_;
     const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
-    fv::options& fvOptions(fv::options::New(this->mesh_));
+    const fv::options& fvOptions(fv::options::New(this->mesh_));
 
     eddyViscosity<RASModel<BasicMomentumTransportModel>>::correct();
 

src/MomentumTransportModels/momentumTransportModels/RAS/kEpsilon/kEpsilon.C

@@ -227,7 +227,7 @@ void kEpsilon<BasicMomentumTransportModel>::correct()
     const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
     const volVectorField& U = this->U_;
     volScalarField& nut = this->nut_;
-    fv::options& fvOptions(fv::options::New(this->mesh_));
+    const fv::options& fvOptions(fv::options::New(this->mesh_));
 
     eddyViscosity<RASModel<BasicMomentumTransportModel>>::correct();
 

src/MomentumTransportModels/momentumTransportModels/RAS/kOmega/kOmega.C

@@ -216,7 +216,7 @@ void kOmega<BasicMomentumTransportModel>::correct()
     const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
     const volVectorField& U = this->U_;
     volScalarField& nut = this->nut_;
-    fv::options& fvOptions(fv::options::New(this->mesh_));
+    const fv::options& fvOptions(fv::options::New(this->mesh_));
 
     eddyViscosity<RASModel<BasicMomentumTransportModel>>::correct();
 

src/MomentumTransportModels/momentumTransportModels/RAS/kOmegaSSTLM/kOmegaSSTLM.C

@@ -497,7 +497,7 @@ void kOmegaSSTLM<BasicMomentumTransportModel>::correctReThetatGammaInt()
     const tmp<volScalarField> tnu = this->nu();
     const volScalarField::Internal& nu = tnu()();
     const volScalarField::Internal& y = this->y_();
-    fv::options& fvOptions(fv::options::New(this->mesh_));
+    const fv::options& fvOptions(fv::options::New(this->mesh_));
 
     // Fields derived from the velocity gradient
     tmp<volTensorField> tgradU = fvc::grad(U);

src/MomentumTransportModels/momentumTransportModels/RAS/realizableKE/realizableKE.C

@@ -256,7 +256,7 @@ void realizableKE<BasicMomentumTransportModel>::correct()
     const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
     const volVectorField& U = this->U_;
     volScalarField& nut = this->nut_;
-    fv::options& fvOptions(fv::options::New(this->mesh_));
+    const fv::options& fvOptions(fv::options::New(this->mesh_));
 
     eddyViscosity<RASModel<BasicMomentumTransportModel>>::correct();
 

src/MomentumTransportModels/momentumTransportModels/RAS/v2f/v2f.C

@@ -282,7 +282,7 @@ void v2f<BasicMomentumTransportModel>::correct()
     const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
     const volVectorField& U = this->U_;
     volScalarField& nut = this->nut_;
-    fv::options& fvOptions(fv::options::New(this->mesh_));
+    const fv::options& fvOptions(fv::options::New(this->mesh_));
 
     eddyViscosity<RASModel<BasicMomentumTransportModel>>::correct();
 

src/MomentumTransportModels/momentumTransportModels/laminar/Maxwell/Maxwell.C

@@ -322,7 +322,7 @@ void Maxwell<BasicMomentumTransportModel>::correct()
     const rhoField& rho = this->rho_;
     const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
     const volVectorField& U = this->U_;
-    fv::options& fvOptions(fv::options::New(this->mesh_));
+    const fv::options& fvOptions(fv::options::New(this->mesh_));
 
     laminarModel<BasicMomentumTransportModel>::correct();
 

src/MomentumTransportModels/phaseCompressible/RAS/mixtureKEpsilon/mixtureKEpsilon.C

@@ -596,7 +596,7 @@ void mixtureKEpsilon<BasicMomentumTransportModel>::correct()
     volScalarField& km = km_();
     volScalarField& epsilonm = epsilonm_();
 
-    fv::options& fvOptions(fv::options::New(this->mesh_));
+    const fv::options& fvOptions(fv::options::New(this->mesh_));
 
     eddyViscosity<RASModel<BasicMomentumTransportModel>>::correct();
 

src/atmosphericModels/kEpsilonLopesdaCosta/kEpsilonLopesdaCosta.C

@@ -98,7 +98,7 @@ template<class BasicMomentumTransportModel>
 void kEpsilonLopesdaCosta<BasicMomentumTransportModel>::
 setPorosityCoefficients()
 {
-    fv::options::optionList& fvOptions(fv::options::New(this->mesh_));
+    const fv::options::optionList& fvOptions(fv::options::New(this->mesh_));
 
     forAll(fvOptions, i)
     {
@@ -398,7 +398,7 @@ void kEpsilonLopesdaCosta<BasicMomentumTransportModel>::correct()
     const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
     const volVectorField& U = this->U_;
     volScalarField& nut = this->nut_;
-    fv::options& fvOptions(fv::options::New(this->mesh_));
+    const fv::options& fvOptions(fv::options::New(this->mesh_));
 
     eddyViscosity<RASModel<BasicMomentumTransportModel>>::correct();
 

src/finiteVolume/cfdTools/general/fvOptions/fvOption.C

@@ -104,7 +104,7 @@ Foam::fv::option::~option()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-bool Foam::fv::option::isActive()
+bool Foam::fv::option::isActive() const
 {
     return active_;
 }
@@ -134,7 +134,7 @@ void Foam::fv::option::addSup
 (
     fvMatrix<scalar>& eqn,
     const label fieldi
-)
+) const
 {}
 
 
@@ -142,7 +142,7 @@ void Foam::fv::option::addSup
 (
     fvMatrix<vector>& eqn,
     const label fieldi
-)
+) const
 {}
 
 
@@ -150,7 +150,7 @@ void Foam::fv::option::addSup
 (
     fvMatrix<sphericalTensor>& eqn,
     const label fieldi
-)
+) const
 {}
 
 
@@ -158,7 +158,7 @@ void Foam::fv::option::addSup
 (
     fvMatrix<symmTensor>& eqn,
     const label fieldi
-)
+) const
 {}
 
 
@@ -166,7 +166,7 @@ void Foam::fv::option::addSup
 (
     fvMatrix<tensor>& eqn,
     const label fieldi
-)
+) const
 {}
 
 
@@ -175,7 +175,7 @@ void Foam::fv::option::addSup
     const volScalarField& rho,
     fvMatrix<scalar>& eqn,
     const label fieldi
-)
+) const
 {}
 
 
@@ -184,7 +184,7 @@ void Foam::fv::option::addSup
     const volScalarField& rho,
     fvMatrix<vector>& eqn,
     const label fieldi
-)
+) const
 {}
 
 
@@ -193,7 +193,7 @@ void Foam::fv::option::addSup
     const volScalarField& rho,
     fvMatrix<sphericalTensor>& eqn,
     const label fieldi
-)
+) const
 {}
 
 
@@ -202,7 +202,7 @@ void Foam::fv::option::addSup
     const volScalarField& rho,
     fvMatrix<symmTensor>& eqn,
     const label fieldi
-)
+) const
 {}
 
 
@@ -211,7 +211,7 @@ void Foam::fv::option::addSup
     const volScalarField& rho,
     fvMatrix<tensor>& eqn,
     const label fieldi
-)
+) const
 {}
 
 
@@ -221,7 +221,7 @@ void Foam::fv::option::addSup
     const volScalarField& rho,
     fvMatrix<scalar>& eqn,
     const label fieldi
-)
+) const
 {
     addSup(alpha*rho, eqn, fieldi);
 }
@@ -233,7 +233,7 @@ void Foam::fv::option::addSup
     const volScalarField& rho,
     fvMatrix<vector>& eqn,
     const label fieldi
-)
+) const
 {
     addSup(alpha*rho, eqn, fieldi);
 }
@@ -245,7 +245,7 @@ void Foam::fv::option::addSup
     const volScalarField& rho,
     fvMatrix<sphericalTensor>& eqn,
     const label fieldi
-)
+) const
 {
     addSup(alpha*rho, eqn, fieldi);
 }
@@ -257,7 +257,7 @@ void Foam::fv::option::addSup
     const volScalarField& rho,
     fvMatrix<symmTensor>& eqn,
     const label fieldi
-)
+) const
 {
     addSup(alpha*rho, eqn, fieldi);
 }
@@ -269,17 +269,25 @@ void Foam::fv::option::addSup
     const volScalarField& rho,
     fvMatrix<tensor>& eqn,
     const label fieldi
-)
+) const
 {
     addSup(alpha*rho, eqn, fieldi);
 }
 
 
-void Foam::fv::option::constrain(fvMatrix<scalar>& eqn, const label fieldi)
+void Foam::fv::option::constrain
+(
+    fvMatrix<scalar>& eqn,
+    const label fieldi
+) const
 {}
 
 
-void Foam::fv::option::constrain(fvMatrix<vector>& eqn, const label fieldi)
+void Foam::fv::option::constrain
+(
+    fvMatrix<vector>& eqn,
+    const label fieldi
+) const
 {}
 
 
@@ -287,7 +295,7 @@ void Foam::fv::option::constrain
 (
     fvMatrix<sphericalTensor>& eqn,
     const label fieldi
-)
+) const
 {}
 
 
@@ -295,31 +303,34 @@ void Foam::fv::option::constrain
 (
     fvMatrix<symmTensor>& eqn,
     const label fieldi
-)
+) const
 {}
 
 
-void Foam::fv::option::constrain(fvMatrix<tensor>& eqn, const label fieldi)
+void Foam::fv::option::constrain
+(
+    fvMatrix<tensor>& eqn, const label fieldi
+) const
 {}
 
 
-void Foam::fv::option::correct(volScalarField& field)
+void Foam::fv::option::correct(volScalarField& field) const
 {}
 
 
-void Foam::fv::option::correct(volVectorField& field)
+void Foam::fv::option::correct(volVectorField& field) const
 {}
 
 
-void Foam::fv::option::correct(volSphericalTensorField& field)
+void Foam::fv::option::correct(volSphericalTensorField& field) const
 {}
 
 
-void Foam::fv::option::correct(volSymmTensorField& field)
+void Foam::fv::option::correct(volSymmTensorField& field) const
 {}
 
 
-void Foam::fv::option::correct(volTensorField& field)
+void Foam::fv::option::correct(volTensorField& field) const
 {}
 
 

src/finiteVolume/cfdTools/general/fvOptions/fvOption.H

@@ -48,7 +48,6 @@ SourceFiles
 #include "volFieldsFwd.H"
 #include "dictionary.H"
 #include "Switch.H"
-#include "runTimeSelectionTables.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -93,7 +92,7 @@ protected:
         wordList fieldNames_;
 
         //- Applied flag list - corresponds to each fieldNames_ entry
-        List<bool> applied_;
+        mutable List<bool> applied_;
 
 
 public:
@@ -203,7 +202,7 @@ public:
             inline bool active() const;
 
             //- Set the applied flag to true for field index fieldi
-            inline void setApplied(const label fieldi);
+            inline void setApplied(const label fieldi) const;
 
 
         // Edit
@@ -215,7 +214,7 @@ public:
         // Checks
 
             //- Is the source active?
-            virtual bool isActive();
+            virtual bool isActive() const;
 
             //- Return index of field name if found in fieldNames list
             virtual label applyToField(const word& fieldName) const;
@@ -232,31 +231,31 @@ public:
                 (
                     fvMatrix<scalar>& eqn,
                     const label fieldi
-                );
+                ) const;
 
                 virtual void addSup
                 (
                     fvMatrix<vector>& eqn,
                     const label fieldi
-                );
+                ) const;
 
                 virtual void addSup
                 (
                     fvMatrix<symmTensor>& eqn,
                     const label fieldi
-                );
+                ) const;
 
                 virtual void addSup
                 (
                     fvMatrix<sphericalTensor>& eqn,
                     const label fieldi
-                );
+                ) const;
 
                 virtual void addSup
                 (
                     fvMatrix<tensor>& eqn,
                     const label fieldi
-                );
+                ) const;
 
 
             // Explicit and implicit sources for compressible equations
@@ -266,35 +265,35 @@ public:
                     const volScalarField& rho,
                     fvMatrix<scalar>& eqn,
                     const label fieldi
-                );
+                ) const;
 
                 virtual void addSup
                 (
                     const volScalarField& rho,
                     fvMatrix<vector>& eqn,
                     const label fieldi
-                );
+                ) const;
 
                 virtual void addSup
                 (
                     const volScalarField& rho,
                     fvMatrix<symmTensor>& eqn,
                     const label fieldi
-                );
+                ) const;
 
                 virtual void addSup
                 (
                     const volScalarField& rho,
                     fvMatrix<sphericalTensor>& eqn,
                     const label fieldi
-                );
+                ) const;
 
                 virtual void addSup
                 (
                     const volScalarField& rho,
                     fvMatrix<tensor>& eqn,
                     const label fieldi
-                );
+                ) const;
 
 
             // Explicit and implicit sources for phase equations
@@ -305,7 +304,7 @@ public:
                     const volScalarField& rho,
                     fvMatrix<scalar>& eqn,
                     const label fieldi
-                );
+                ) const;
 
                 virtual void addSup
                 (
@@ -313,7 +312,7 @@ public:
                     const volScalarField& rho,
                     fvMatrix<vector>& eqn,
                     const label fieldi
-                );
+                ) const;
 
                 virtual void addSup
                 (
@@ -321,7 +320,7 @@ public:
                     const volScalarField& rho,
                     fvMatrix<symmTensor>& eqn,
                     const label fieldi
-                );
+                ) const;
 
                 virtual void addSup
                 (
@@ -329,7 +328,7 @@ public:
                     const volScalarField& rho,
                     fvMatrix<sphericalTensor>& eqn,
                     const label fieldi
-                );
+                ) const;
 
                 virtual void addSup
                 (
@@ -337,7 +336,7 @@ public:
                     const volScalarField& rho,
                     fvMatrix<tensor>& eqn,
                     const label fieldi
-                );
+                ) const;
 
 
             // Constraints
@@ -346,40 +345,40 @@ public:
                 (
                     fvMatrix<scalar>& eqn,
                     const label fieldi
-                );
+                ) const;
 
                 virtual void constrain
                 (
                     fvMatrix<vector>& eqn,
                     const label fieldi
-                );
+                ) const;
 
                 virtual void constrain
                 (
                     fvMatrix<sphericalTensor>& eqn,
                     const label fieldi
-                );
+                ) const;
 
                 virtual void constrain
                 (
                     fvMatrix<symmTensor>& eqn,
                     const label fieldi
-                );
+                ) const;
 
                 virtual void constrain
                 (
                     fvMatrix<tensor>& eqn,
                     const label fieldi
-                );
+                ) const;
 
 
             // Correction
 
-                virtual void correct(volScalarField& field);
+                virtual void correct(volScalarField& field) const;
-                virtual void correct(volVectorField& field);
+                virtual void correct(volVectorField& field) const;
-                virtual void correct(volSphericalTensorField& field);
+                virtual void correct(volSphericalTensorField& field) const;
-                virtual void correct(volSymmTensorField& field);
+                virtual void correct(volSymmTensorField& field) const;
-                virtual void correct(volTensorField& field);
+                virtual void correct(volTensorField& field) const;
 
 
         // Mesh changes

src/finiteVolume/cfdTools/general/fvOptions/fvOptionI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -49,7 +49,7 @@ inline bool Foam::fv::option::active() const
 }
 
 
-inline void Foam::fv::option::setApplied(const label fieldi)
+inline void Foam::fv::option::setApplied(const label fieldi) const
 {
     applied_[fieldi] = true;
 }

src/finiteVolume/cfdTools/general/fvOptions/fvOptionList.H

@@ -92,7 +92,7 @@ protected:
             GeometricField<Type, fvPatchField, volMesh>& field,
             const word& fieldName,
             const dimensionSet& ds
-        );
+        ) const;
 
 
 private:
@@ -142,7 +142,7 @@ public:
             tmp<fvMatrix<Type>> operator()
             (
                 GeometricField<Type, fvPatchField, volMesh>& field
-            );
+            ) const;
 
             //- Return source for equation with specified name
             template<class Type>
@@ -150,7 +150,7 @@ public:
             (
                 GeometricField<Type, fvPatchField, volMesh>& field,
                 const word& fieldName
-            );
+            ) const;
 
             //- Return source for equation
             template<class Type>
@@ -158,7 +158,7 @@ public:
             (
                 const volScalarField& rho,
                 GeometricField<Type, fvPatchField, volMesh>& field
-            );
+            ) const;
 
             //- Return source for equation with specified name
             template<class Type>
@@ -167,7 +167,7 @@ public:
                 const volScalarField& rho,
                 GeometricField<Type, fvPatchField, volMesh>& field,
                 const word& fieldName
-            );
+            ) const;
 
             //- Return source for equation
             template<class Type>
@@ -176,7 +176,7 @@ public:
                 const volScalarField& alpha,
                 const volScalarField& rho,
                 GeometricField<Type, fvPatchField, volMesh>& field
-            );
+            ) const;
 
             //- Return source for equation with specified name
             template<class Type>
@@ -186,7 +186,7 @@ public:
                 const volScalarField& rho,
                 GeometricField<Type, fvPatchField, volMesh>& field,
                 const word& fieldName
-            );
+            ) const;
 
             //- Return source for equation
             template<class Type>
@@ -195,7 +195,7 @@ public:
                 const volScalarField& alpha,
                 const geometricOneField& rho,
                 GeometricField<Type, fvPatchField, volMesh>& field
-            );
+            ) const;
 
             //- Return source for equation
             template<class Type>
@@ -204,7 +204,7 @@ public:
                 const geometricOneField& alpha,
                 const volScalarField& rho,
                 GeometricField<Type, fvPatchField, volMesh>& field
-            );
+            ) const;
 
             //- Return source for equation
             template<class Type>
@@ -213,14 +213,14 @@ public:
                 const geometricOneField& alpha,
                 const geometricOneField& rho,
                 GeometricField<Type, fvPatchField, volMesh>& field
-            );
+            ) const;
 
             //- Return source for equation with second time derivative
             template<class Type>
             tmp<fvMatrix<Type>> d2dt2
             (
                 GeometricField<Type, fvPatchField, volMesh>& field
-            );
+            ) const;
 
             //- Return source for equation with second time derivative
             template<class Type>
@@ -228,21 +228,24 @@ public:
             (
                 GeometricField<Type, fvPatchField, volMesh>& field,
                 const word& fieldName
-            );
+            ) const;
 
 
         // Constraints
 
             //- Apply constraints to equation
             template<class Type>
-            void constrain(fvMatrix<Type>& eqn);
+            void constrain(fvMatrix<Type>& eqn) const;
 
 
         // Correction
 
             //- Apply correction to field
             template<class Type>
-            void correct(GeometricField<Type, fvPatchField, volMesh>& field);
+            void correct
+            (
+                GeometricField<Type, fvPatchField, volMesh>& field
+            ) const;
 
 
         // Mesh changes

src/finiteVolume/cfdTools/general/fvOptions/fvOptionListTemplates.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -31,7 +31,7 @@ Foam::tmp<Foam::fvMatrix<Type>> Foam::fv::optionList::source
     GeometricField<Type, fvPatchField, volMesh>& field,
     const word& fieldName,
     const dimensionSet& ds
-)
+) const
 {
     checkApplied();
 
@@ -40,7 +40,7 @@ Foam::tmp<Foam::fvMatrix<Type>> Foam::fv::optionList::source
 
     forAll(*this, i)
     {
-        option& source = this->operator[](i);
+        const option& source = this->operator[](i);
 
         label fieldi = source.applyToField(fieldName);
 
@@ -69,7 +69,7 @@ template<class Type>
 Foam::tmp<Foam::fvMatrix<Type>> Foam::fv::optionList::operator()
 (
     GeometricField<Type, fvPatchField, volMesh>& field
-)
+) const
 {
     return this->operator()(field, field.name());
 }
@@ -80,7 +80,7 @@ Foam::tmp<Foam::fvMatrix<Type>> Foam::fv::optionList::operator()
 (
     GeometricField<Type, fvPatchField, volMesh>& field,
     const word& fieldName
-)
+) const
 {
     return source(field, fieldName, field.dimensions()/dimTime*dimVolume);
 }
@@ -91,7 +91,7 @@ Foam::tmp<Foam::fvMatrix<Type>> Foam::fv::optionList::operator()
 (
     const volScalarField& rho,
     GeometricField<Type, fvPatchField, volMesh>& field
-)
+) const
 {
     return this->operator()(rho, field, field.name());
 }
@@ -103,7 +103,7 @@ Foam::tmp<Foam::fvMatrix<Type>> Foam::fv::optionList::operator()
     const volScalarField& rho,
     GeometricField<Type, fvPatchField, volMesh>& field,
     const word& fieldName
-)
+) const
 {
     checkApplied();
 
@@ -117,7 +117,7 @@ Foam::tmp<Foam::fvMatrix<Type>> Foam::fv::optionList::operator()
 
     forAll(*this, i)
     {
-        option& source = this->operator[](i);
+        const option& source = this->operator[](i);
 
         label fieldi = source.applyToField(fieldName);
 
@@ -148,7 +148,7 @@ Foam::tmp<Foam::fvMatrix<Type>> Foam::fv::optionList::operator()
     const volScalarField& alpha,
     const volScalarField& rho,
     GeometricField<Type, fvPatchField, volMesh>& field
-)
+) const
 {
     return this->operator()(alpha, rho, field, field.name());
 }
@@ -161,7 +161,7 @@ Foam::tmp<Foam::fvMatrix<Type>> Foam::fv::optionList::operator()
     const volScalarField& rho,
     GeometricField<Type, fvPatchField, volMesh>& field,
     const word& fieldName
-)
+) const
 {
     checkApplied();
 
@@ -176,7 +176,7 @@ Foam::tmp<Foam::fvMatrix<Type>> Foam::fv::optionList::operator()
 
     forAll(*this, i)
     {
-        option& source = this->operator[](i);
+        const option& source = this->operator[](i);
 
         label fieldi = source.applyToField(fieldName);
 
@@ -207,7 +207,7 @@ Foam::tmp<Foam::fvMatrix<Type>> Foam::fv::optionList::operator()
     const geometricOneField& alpha,
     const geometricOneField& rho,
     GeometricField<Type, fvPatchField, volMesh>& field
-)
+) const
 {
     return this->operator()(field, field.name());
 }
@@ -219,7 +219,7 @@ Foam::tmp<Foam::fvMatrix<Type>> Foam::fv::optionList::operator()
     const volScalarField& alpha,
     const geometricOneField& rho,
     GeometricField<Type, fvPatchField, volMesh>& field
-)
+) const
 {
     volScalarField one
     (
@@ -246,7 +246,7 @@ Foam::tmp<Foam::fvMatrix<Type>> Foam::fv::optionList::operator()
     const geometricOneField& alpha,
     const volScalarField& rho,
     GeometricField<Type, fvPatchField, volMesh>& field
-)
+) const
 {
     return this->operator()(rho, field, field.name());
 }
@@ -256,7 +256,7 @@ template<class Type>
 Foam::tmp<Foam::fvMatrix<Type>> Foam::fv::optionList::d2dt2
 (
     GeometricField<Type, fvPatchField, volMesh>& field
-)
+) const
 {
     return this->d2dt2(field, field.name());
 }
@@ -267,20 +267,20 @@ Foam::tmp<Foam::fvMatrix<Type>> Foam::fv::optionList::d2dt2
 (
     GeometricField<Type, fvPatchField, volMesh>& field,
     const word& fieldName
-)
+) const
 {
     return source(field, fieldName, field.dimensions()/sqr(dimTime)*dimVolume);
 }
 
 
 template<class Type>
-void Foam::fv::optionList::constrain(fvMatrix<Type>& eqn)
+void Foam::fv::optionList::constrain(fvMatrix<Type>& eqn) const
 {
     checkApplied();
 
     forAll(*this, i)
     {
-        option& source = this->operator[](i);
+        const option& source = this->operator[](i);
 
         label fieldi = source.applyToField(eqn.psi().name());
 
@@ -307,13 +307,13 @@ template<class Type>
 void Foam::fv::optionList::correct
 (
     GeometricField<Type, fvPatchField, volMesh>& field
-)
+) const
 {
     const word& fieldName = field.name();
 
     forAll(*this, i)
     {
-        option& source = this->operator[](i);
+        const option& source = this->operator[](i);
 
         label fieldi = source.applyToField(fieldName);
 

src/finiteVolume/cfdTools/general/fvOptions/fvOptions.H

@@ -75,17 +75,7 @@ public:
         //- Disallow default bitwise copy construction
         options(const options&) = delete;
 
-        // using MeshObject<fvMesh, UpdateableMeshObject, options>::New;
+        using MeshObject<fvMesh, UpdateableMeshObject, options>::New;
-
-        //- Construct fvOptions and register to datbase if not present
-        //  otherwise lookup and return
-        static options& New(const fvMesh& mesh)
-        {
-            return const_cast<options&>
-            (
-                MeshObject<fvMesh, UpdateableMeshObject, options>::New(mesh)
-            );
-        }
 
 
     //- Destructor

src/finiteVolume/cfdTools/general/include/createFvOptions.H

@@ -1,4 +1,4 @@
-fv::options& fvOptions(fv::options::New(mesh));
+const fv::options& fvOptions(fv::options::New(mesh));
 
 if (!fvOptions.optionList::size())
 {

src/fvOptions/cellSetOption/cellSetOption.C

@@ -215,7 +215,7 @@ Foam::fv::cellSetOption::~cellSetOption()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-bool Foam::fv::cellSetOption::isActive()
+bool Foam::fv::cellSetOption::isActive() const
 {
     if (option::isActive() && inTimeLimits(mesh_.time().value()))
     {

src/fvOptions/cellSetOption/cellSetOption.H

@@ -185,7 +185,7 @@ public:
         // Checks
 
             //- Is the source active?
-            virtual bool isActive();
+            virtual bool isActive() const;
 
 
         // Mesh changes

src/fvOptions/constraints/fixedTemperatureConstraint/fixedTemperatureConstraint.C

@@ -115,7 +115,7 @@ void Foam::fv::fixedTemperatureConstraint::constrain
 (
     fvMatrix<scalar>& eqn,
     const label
-)
+) const
 {
     const basicThermo& thermo =
         mesh_.lookupObject<basicThermo>

src/fvOptions/constraints/fixedTemperatureConstraint/fixedTemperatureConstraint.H

@@ -140,7 +140,11 @@ public:
     // Member Functions
 
         //- Constrain energy equation to fix the temperature
-        virtual void constrain(fvMatrix<scalar>& eqn, const label fieldi);
+        virtual void constrain
+        (
+            fvMatrix<scalar>& eqn,
+            const label fieldi
+        ) const;
 
         //- Read dictionary
         virtual bool read(const dictionary& dict);

src/fvOptions/constraints/fixedValueConstraint/FixedValueConstraint.C

@@ -81,7 +81,7 @@ void Foam::fv::FixedValueConstraint<Type>::constrain
 (
     fvMatrix<Type>& eqn,
     const label fieldi
-)
+) const
 {
     DebugInfo
         << "FixedValueConstraint<"

src/fvOptions/constraints/fixedValueConstraint/FixedValueConstraint.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2016-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -105,7 +105,7 @@ public:
         virtual bool read(const dictionary& dict);
 
         //- Set value on field
-        virtual void constrain(fvMatrix<Type>& eqn, const label fieldi);
+        virtual void constrain(fvMatrix<Type>& eqn, const label fieldi) const;
 };
 
 

src/fvOptions/corrections/limitTemperature/limitTemperature.C

@@ -92,7 +92,7 @@ bool Foam::fv::limitTemperature::read(const dictionary& dict)
 }
 
 
-void Foam::fv::limitTemperature::correct(volScalarField& he)
+void Foam::fv::limitTemperature::correct(volScalarField& he) const
 {
     const basicThermo& thermo =
         mesh_.lookupObject<basicThermo>

src/fvOptions/corrections/limitTemperature/limitTemperature.H

@@ -114,7 +114,7 @@ public:
         virtual bool read(const dictionary& dict);
 
         //- Correct the energy field
-        virtual void correct(volScalarField& he);
+        virtual void correct(volScalarField& he) const;
 
 
     // Member Operators

src/fvOptions/corrections/limitVelocity/limitVelocity.C

@@ -80,7 +80,7 @@ bool Foam::fv::limitVelocity::read(const dictionary& dict)
 }
 
 
-void Foam::fv::limitVelocity::correct(volVectorField& U)
+void Foam::fv::limitVelocity::correct(volVectorField& U) const
 {
     const scalar maxSqrU = sqr(max_);
 

src/fvOptions/corrections/limitVelocity/limitVelocity.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2016-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -108,7 +108,7 @@ public:
         virtual bool read(const dictionary& dict);
 
         //- Correct the energy field
-        virtual void correct(volVectorField& U);
+        virtual void correct(volVectorField& U) const;
 
 
     // Member Operators

src/fvOptions/sources/derived/accelerationSource/accelerationSource.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2018-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -64,7 +64,7 @@ void Foam::fv::accelerationSource::addSup
 (
     fvMatrix<vector>& eqn,
     const label fieldi
-)
+) const
 {
     add(geometricOneField(), eqn, fieldi);
 }
@@ -75,7 +75,7 @@ void Foam::fv::accelerationSource::addSup
     const volScalarField& rho,
     fvMatrix<vector>& eqn,
     const label fieldi
-)
+) const
 {
     add(rho, eqn, fieldi);
 }
@@ -87,7 +87,7 @@ void Foam::fv::accelerationSource::addSup
     const volScalarField& rho,
     fvMatrix<vector>& eqn,
     const label fieldi
-)
+) const
 {
     add((alpha*rho)(), eqn, fieldi);
 }

src/fvOptions/sources/derived/accelerationSource/accelerationSource.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2018-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2018-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -89,7 +89,7 @@ class accelerationSource
             const AlphaRhoFieldType& rho,
             fvMatrix<vector>& eqn,
             const label fieldi
-        );
+        ) const;
 
 
 public:
@@ -124,7 +124,7 @@ public:
             (
                 fvMatrix<vector>& eqn,
                 const label fieldi
-            );
+            ) const;
 
             //- Source term to compressible momentum equation
             virtual void addSup
@@ -132,7 +132,7 @@ public:
                 const volScalarField& rho,
                 fvMatrix<vector>& eqn,
                 const label fieldi
-            );
+            ) const;
 
             //- Source term to phase momentum equation
             virtual void addSup
@@ -141,7 +141,7 @@ public:
                 const volScalarField& rho,
                 fvMatrix<vector>& eqn,
                 const label fieldi
-            );
+            ) const;
 
 
         // IO

src/fvOptions/sources/derived/accelerationSource/accelerationSourceTemplates.C

@@ -31,7 +31,7 @@ void Foam::fv::accelerationSource::add
     const AlphaRhoFieldType& alphaRho,
     fvMatrix<vector>& eqn,
     const label fieldi
-)
+) const
 {
     const DimensionedField<scalar, volMesh>& V = mesh_.V();
 

src/fvOptions/sources/derived/actuationDiskSource/actuationDiskSource.C

@@ -113,7 +113,7 @@ void Foam::fv::actuationDiskSource::addSup
 (
     fvMatrix<vector>& eqn,
     const label fieldi
-)
+) const
 {
     const scalarField& cellsV = mesh_.V();
     vectorField& Usource = eqn.source();
@@ -138,7 +138,7 @@ void Foam::fv::actuationDiskSource::addSup
     const volScalarField& rho,
     fvMatrix<vector>& eqn,
     const label fieldi
-)
+) const
 {
     const scalarField& cellsV = mesh_.V();
     vectorField& Usource = eqn.source();

src/fvOptions/sources/derived/actuationDiskSource/actuationDiskSource.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -187,7 +187,7 @@ public:
             (
                 fvMatrix<vector>& eqn,
                 const label fieldi
-            );
+            ) const;
 
             //- Source term to compressible momentum equation
             virtual void addSup
@@ -195,7 +195,7 @@ public:
                 const volScalarField& rho,
                 fvMatrix<vector>& eqn,
                 const label fieldi
-            );
+            ) const;
 
 
         // IO

src/fvOptions/sources/derived/buoyancyEnergy/buoyancyEnergy.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2015-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -77,7 +77,7 @@ void Foam::fv::buoyancyEnergy::addSup
     const volScalarField& rho,
     fvMatrix<scalar>& eqn,
     const label fieldi
-)
+) const
 {
     const uniformDimensionedVectorField& g =
         mesh_.lookupObject<uniformDimensionedVectorField>("g");

src/fvOptions/sources/derived/buoyancyEnergy/buoyancyEnergy.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2015-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -97,7 +97,7 @@ public:
                 const volScalarField& rho,
                 fvMatrix<scalar>& eqn,
                 const label fieldi
-            );
+            ) const;
 
 
         // IO

src/fvOptions/sources/derived/buoyancyForce/buoyancyForce.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2015-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -86,7 +86,7 @@ void Foam::fv::buoyancyForce::addSup
 (
     fvMatrix<vector>& eqn,
     const label fieldi
-)
+) const
 {
     eqn += g_;
 }
@@ -97,7 +97,7 @@ void Foam::fv::buoyancyForce::addSup
     const volScalarField& rho,
     fvMatrix<vector>& eqn,
     const label fieldi
-)
+) const
 {
     eqn += rho*g_;
 }

src/fvOptions/sources/derived/buoyancyForce/buoyancyForce.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2015-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -95,7 +95,7 @@ public:
             (
                 fvMatrix<vector>& eqn,
                 const label fieldi
-            );
+            ) const;
 
             //- Add explicit contribution to compressible momentum equation
             virtual void addSup
@@ -103,7 +103,7 @@ public:
                 const volScalarField& rho,
                 fvMatrix<vector>& eqn,
                 const label fieldi
-            );
+            ) const;
 
 
         // IO

src/fvOptions/sources/derived/damping/isotropicDamping/isotropicDamping.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2019-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -46,7 +46,7 @@ void Foam::fv::isotropicDamping::add
 (
     const volScalarField::Internal& forceCoeff,
     fvMatrix<vector>& eqn
-)
+) const
 {
     eqn -= fvm::Sp(forceCoeff, eqn.psi());
     eqn += forceCoeff*value_;
@@ -76,7 +76,7 @@ void Foam::fv::isotropicDamping::addSup
 (
     fvMatrix<vector>& eqn,
     const label fieldi
-)
+) const
 {
     add(this->forceCoeff(), eqn);
 }
@@ -87,7 +87,7 @@ void Foam::fv::isotropicDamping::addSup
     const volScalarField& rho,
     fvMatrix<vector>& eqn,
     const label fieldi
-)
+) const
 {
     add(rho*forceCoeff(), eqn);
 }
@@ -99,7 +99,7 @@ void Foam::fv::isotropicDamping::addSup
     const volScalarField& rho,
     fvMatrix<vector>& eqn,
     const label fieldi
-)
+) const
 {
     add(alpha()*rho()*this->forceCoeff(), eqn);
 }

src/fvOptions/sources/derived/damping/isotropicDamping/isotropicDamping.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2019-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -129,7 +129,7 @@ class isotropicDamping
         (
             const volScalarField::Internal& forceCoeff,
             fvMatrix<vector>& eqn
-        );
+        ) const;
 
 
 public:
@@ -164,7 +164,7 @@ public:
             (
                 fvMatrix<vector>& eqn,
                 const label fieldi
-            );
+            ) const;
 
             //- Source term to compressible momentum equation
             virtual void addSup
@@ -172,7 +172,7 @@ public:
                 const volScalarField& rho,
                 fvMatrix<vector>& eqn,
                 const label fieldi
-            );
+            ) const;
 
             //- Source term to phase momentum equation
             virtual void addSup
@@ -181,7 +181,7 @@ public:
                 const volScalarField& rho,
                 fvMatrix<vector>& eqn,
                 const label fieldi
-            );
+            ) const;
 
 
         // IO

src/fvOptions/sources/derived/damping/verticalDamping/verticalDamping.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2017-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2017-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -46,7 +46,7 @@ void Foam::fv::verticalDamping::add
 (
     const volVectorField& alphaRhoU,
     fvMatrix<vector>& eqn
-)
+) const
 {
     const uniformDimensionedVectorField& g =
         mesh_.lookupObject<uniformDimensionedVectorField>("g");
@@ -79,7 +79,7 @@ void Foam::fv::verticalDamping::addSup
 (
     fvMatrix<vector>& eqn,
     const label fieldi
-)
+) const
 {
     add(eqn.psi(), eqn);
 }
@@ -90,7 +90,7 @@ void Foam::fv::verticalDamping::addSup
     const volScalarField& rho,
     fvMatrix<vector>& eqn,
     const label fieldi
-)
+) const
 {
     add(rho*eqn.psi(), eqn);
 }
@@ -102,7 +102,7 @@ void Foam::fv::verticalDamping::addSup
     const volScalarField& rho,
     fvMatrix<vector>& eqn,
     const label fieldi
-)
+) const
 {
     add(alpha*rho*eqn.psi(), eqn);
 }

src/fvOptions/sources/derived/damping/verticalDamping/verticalDamping.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2017-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2017-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -133,7 +133,7 @@ class verticalDamping
         (
             const volVectorField& alphaRhoU,
             fvMatrix<vector>& eqn
-        );
+        ) const;
 
 
 public:
@@ -168,7 +168,7 @@ public:
             (
                 fvMatrix<vector>& eqn,
                 const label fieldi
-            );
+            ) const;
 
             //- Source term to compressible momentum equation
             virtual void addSup
@@ -176,7 +176,7 @@ public:
                 const volScalarField& rho,
                 fvMatrix<vector>& eqn,
                 const label fieldi
-            );
+            ) const;
 
             //- Source term to phase momentum equation
             virtual void addSup
@@ -185,7 +185,7 @@ public:
                 const volScalarField& rho,
                 fvMatrix<vector>& eqn,
                 const label fieldi
-            );
+            ) const;
 };
 
 

src/fvOptions/sources/derived/effectivenessHeatExchangerSource/effectivenessHeatExchangerSource.C

@@ -121,7 +121,7 @@ void Foam::fv::effectivenessHeatExchangerSource::initialise()
 void Foam::fv::effectivenessHeatExchangerSource::calculateTotalArea
 (
     scalar& area
-)
+) const
 {
     area = 0;
     forAll(faceId_, i)
@@ -198,7 +198,7 @@ void Foam::fv::effectivenessHeatExchangerSource::addSup
     const volScalarField& rho,
     fvMatrix<scalar>& eqn,
     const label
-)
+) const
 {
     const basicThermo& thermo =
         mesh_.lookupObject<basicThermo>(basicThermo::dictName);

src/fvOptions/sources/derived/effectivenessHeatExchangerSource/effectivenessHeatExchangerSource.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2013-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -201,7 +201,7 @@ private:
         void initialise();
 
         //- Calculate total area of faceZone across processors
-        void calculateTotalArea(scalar& area);
+        void calculateTotalArea(scalar& area) const;
 
 
 public:
@@ -242,7 +242,7 @@ public:
             (
                 fvMatrix<scalar>& eqn,
                 const label fieldi
-            )
+            ) const
             {
                 NotImplemented;
             }
@@ -256,7 +256,7 @@ public:
                 const volScalarField& rho,
                 fvMatrix<scalar>& eqn,
                 const label fieldi
-            );
+            ) const;
 
 
         // IO

src/fvOptions/sources/derived/explicitPorositySource/explicitPorositySource.C

@@ -80,7 +80,7 @@ void Foam::fv::explicitPorositySource::addSup
 (
     fvMatrix<vector>& eqn,
     const label fieldi
-)
+) const
 {
     fvMatrix<vector> porosityEqn(eqn.psi(), eqn.dimensions());
     porosityPtr_->addResistance(porosityEqn);
@@ -93,7 +93,7 @@ void Foam::fv::explicitPorositySource::addSup
     const volScalarField& rho,
     fvMatrix<vector>& eqn,
     const label fieldi
-)
+) const
 {
     fvMatrix<vector> porosityEqn(eqn.psi(), eqn.dimensions());
     porosityPtr_->addResistance(porosityEqn);
@@ -107,7 +107,7 @@ void Foam::fv::explicitPorositySource::addSup
     const volScalarField& rho,
     fvMatrix<vector>& eqn,
     const label fieldi
-)
+) const
 {
     fvMatrix<vector> porosityEqn(eqn.psi(), eqn.dimensions());
     porosityPtr_->addResistance(porosityEqn);

src/fvOptions/sources/derived/explicitPorositySource/explicitPorositySource.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2012-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -91,7 +91,7 @@ protected:
     // Protected data
 
         //- Run-time selectable porosity model
-        autoPtr<porosityModel> porosityPtr_;
+        mutable autoPtr<porosityModel> porosityPtr_;
 
 
 public:
@@ -134,7 +134,7 @@ public:
             (
                 fvMatrix<vector>& eqn,
                 const label fieldi
-            );
+            ) const;
 
             //- Add implicit contribution to compressible momentum equation
             virtual void addSup
@@ -142,7 +142,7 @@ public:
                 const volScalarField& rho,
                 fvMatrix<vector>& eqn,
                 const label fieldi
-            );
+            ) const;
 
             //- Add implicit contribution to phase momentum equation
             virtual void addSup
@@ -151,7 +151,7 @@ public:
                 const volScalarField& rho,
                 fvMatrix<vector>& eqn,
                 const label fieldi
-            );
+            ) const;
 
 
         // IO

src/fvOptions/sources/derived/meanVelocityForce/meanVelocityForce.C

@@ -147,7 +147,7 @@ Foam::scalar Foam::fv::meanVelocityForce::magUbarAve
 }
 
 
-void Foam::fv::meanVelocityForce::correct(volVectorField& U)
+void Foam::fv::meanVelocityForce::correct(volVectorField& U) const
 {
     const scalarField& rAU = rAPtr_();
 
@@ -195,7 +195,7 @@ void Foam::fv::meanVelocityForce::addSup
 (
     fvMatrix<vector>& eqn,
     const label fieldi
-)
+) const
 {
     volVectorField::Internal Su
     (
@@ -224,7 +224,7 @@ void Foam::fv::meanVelocityForce::addSup
     const volScalarField& rho,
     fvMatrix<vector>& eqn,
     const label fieldi
-)
+) const
 {
     this->addSup(eqn, fieldi);
 }
@@ -234,7 +234,7 @@ void Foam::fv::meanVelocityForce::constrain
 (
     fvMatrix<vector>& eqn,
     const label
-)
+) const
 {
     if (rAPtr_.empty())
     {

src/fvOptions/sources/derived/meanVelocityForce/meanVelocityForce.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -77,10 +77,10 @@ protected:
         vector Ubar_;
 
         //- Pressure gradient before correction
-        scalar gradP0_;
+        mutable scalar gradP0_;
 
         //- Change in pressure gradient
-        scalar dGradP_;
+        mutable scalar dGradP_;
 
         //- Flow direction
         vector flowDir_;
@@ -89,7 +89,7 @@ protected:
         scalar relaxation_;
 
         //- Matrix 1/A coefficients field pointer
-        autoPtr<volScalarField> rAPtr_;
+        mutable autoPtr<volScalarField> rAPtr_;
 
 
     // Protected Member Functions
@@ -102,7 +102,7 @@ protected:
         void writeProps(const scalar gradP) const;
 
         //- Correct driving force for a constant mass flow rate
-        void update(fvMatrix<vector>& eqn);
+        void update(fvMatrix<vector>& eqn) const;
 
 
 public:
@@ -131,14 +131,14 @@ public:
         // Evaluate
 
             //- Correct the pressure gradient
-            virtual void correct(volVectorField& U);
+            virtual void correct(volVectorField& U) const;
 
             //- Add explicit contribution to momentum equation
             virtual void addSup
             (
                 fvMatrix<vector>& eqn,
                 const label fieldi
-            );
+            ) const;
 
             //- Add explicit contribution to compressible momentum equation
             virtual void addSup
@@ -146,14 +146,14 @@ public:
                 const volScalarField& rho,
                 fvMatrix<vector>& eqn,
                 const label fieldi
-            );
+            ) const;
 
             //- Set 1/A coefficient
             virtual void constrain
             (
                 fvMatrix<vector>& eqn,
                 const label fieldi
-            );
+            ) const;
 
 
         // IO

src/fvOptions/sources/derived/phaseLimitStabilization/PhaseLimitStabilization.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2017-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2017-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -58,7 +58,7 @@ void Foam::fv::PhaseLimitStabilization<Type>::addSup
     const volScalarField& rho,
     fvMatrix<Type>& eqn,
     const label fieldi
-)
+) const
 {
     const GeometricField<Type, fvPatchField, volMesh>& psi = eqn.psi();
 

src/fvOptions/sources/derived/phaseLimitStabilization/PhaseLimitStabilization.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2017-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2017-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -122,7 +122,7 @@ public:
             const volScalarField& rho,
             fvMatrix<Type>& eqn,
             const label fieldi
-        );
+        ) const;
 
         //- Read dictionary
         virtual bool read(const dictionary& dict);

src/fvOptions/sources/derived/radialActuationDiskSource/radialActuationDiskSource.C

@@ -66,7 +66,7 @@ void Foam::fv::radialActuationDiskSource::addSup
 (
     fvMatrix<vector>& eqn,
     const label fieldi
-)
+) const
 {
     const scalarField& cellsV = mesh_.V();
     vectorField& Usource = eqn.source();
@@ -91,7 +91,7 @@ void Foam::fv::radialActuationDiskSource::addSup
     const volScalarField& rho,
     fvMatrix<vector>& eqn,
     const label fieldi
-)
+) const
 {
     const scalarField& cellsV = mesh_.V();
     vectorField& Usource = eqn.source();

src/fvOptions/sources/derived/radialActuationDiskSource/radialActuationDiskSource.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -143,7 +143,7 @@ public:
         (
             fvMatrix<vector>& eqn,
             const label fieldi
-        );
+        ) const;
 
         //- Source term to compressible momentum equation
         virtual void addSup
@@ -151,7 +151,7 @@ public:
             const volScalarField& rho,
             fvMatrix<vector>& eqn,
             const label fieldi
-        );
+        ) const;
 
 
     // IO

src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSource.C

@@ -513,7 +513,7 @@ void Foam::fv::rotorDiskSource::addSup
 (
     fvMatrix<vector>& eqn,
     const label fieldi
-)
+) const
 {
     volVectorField force
     (
@@ -554,7 +554,7 @@ void Foam::fv::rotorDiskSource::addSup
     const volScalarField& rho,
     fvMatrix<vector>& eqn,
     const label fieldi
-)
+) const
 {
     volVectorField force
     (

src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSource.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -158,7 +158,7 @@ protected:
     // Protected data
 
         //- Reference density for incompressible case
-        scalar rhoRef_;
+        mutable scalar rhoRef_;
 
         //- Rotational speed [rad/s]
         //  Positive anti-clockwise when looking along -ve lift direction
@@ -303,7 +303,7 @@ public:
             (
                 fvMatrix<vector>& eqn,
                 const label fieldi
-            );
+            ) const;
 
             //- Source term to compressible momentum equation
             virtual void addSup
@@ -311,7 +311,7 @@ public:
                 const volScalarField& rho,
                 fvMatrix<vector>& eqn,
                 const label fieldi
-            );
+            ) const;
 
 
         // IO

src/fvOptions/sources/derived/rotorDiskSource/trimModel/fixed/fixedTrim.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2012-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -90,7 +90,7 @@ void Foam::fixedTrim::correct
 (
     const vectorField& U,
     vectorField& force
-)
+) const
 {}
 
 
@@ -98,7 +98,8 @@ void Foam::fixedTrim::correct
 (
     const volScalarField rho,
     const vectorField& U,
-    vectorField& force)
+    vectorField& force
+) const
 {}
 
 

src/fvOptions/sources/derived/rotorDiskSource/trimModel/fixed/fixedTrim.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2012-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -84,7 +84,7 @@ public:
         (
             const vectorField& U,
             vectorField& force
-        );
+        ) const;
 
         //- Correct the model for compressible flow
         virtual void correct
@@ -92,7 +92,7 @@ public:
             const volScalarField rho,
             const vectorField& U,
             vectorField& force
-        );
+        ) const;
 };
 
 

src/fvOptions/sources/derived/rotorDiskSource/trimModel/targetCoeff/targetCoeffTrim.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2012-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -101,7 +101,7 @@ void Foam::targetCoeffTrim::correctTrim
     const RhoFieldType& rho,
     const vectorField& U,
     vectorField& force
-)
+) const
 {
     if (rotor_.mesh().time().timeIndex() % calcFrequency_ == 0)
     {
@@ -276,7 +276,7 @@ void Foam::targetCoeffTrim::correct
 (
     const vectorField& U,
     vectorField& force
-)
+) const
 {
     correctTrim(geometricOneField(), U, force);
 }
@@ -287,7 +287,7 @@ void Foam::targetCoeffTrim::correct
     const volScalarField rho,
     const vectorField& U,
     vectorField& force
-)
+) const
 {
     correctTrim(rho, U, force);
 }

src/fvOptions/sources/derived/rotorDiskSource/trimModel/targetCoeff/targetCoeffTrim.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2012-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -104,7 +104,7 @@ protected:
         vector target_;
 
         //- Pitch angles (collective, roll, pitch) [rad]
-        vector theta_;
+        mutable vector theta_;
 
         //- Maximum number of iterations in trim routine
         label nIter_;
@@ -141,7 +141,7 @@ protected:
             const RhoFieldType& rho,
             const vectorField& U,
             vectorField& force
-        );
+        ) const;
 
 
 public:
@@ -169,7 +169,7 @@ public:
         (
             const vectorField& U,
             vectorField& force
-        );
+        ) const;
 
         //- Correct the model for compressible flow
         virtual void correct
@@ -177,7 +177,7 @@ public:
             const volScalarField rho,
             const vectorField& U,
             vectorField& force
-        );
+        ) const;
 };
 
 

src/fvOptions/sources/derived/rotorDiskSource/trimModel/trimModel/trimModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2012-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -124,7 +124,7 @@ public:
         (
             const vectorField& U,
             vectorField& force
-        ) = 0;
+        ) const = 0;
 
         //- Correct the model for compressible flow
         virtual void correct
@@ -132,7 +132,7 @@ public:
             const volScalarField rho,
             const vectorField& U,
             vectorField& force
-        ) = 0;
+        ) const = 0;
 };
 
 

src/fvOptions/sources/derived/solidEquilibriumEnergySource/solidEquilibriumEnergySource.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2019-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -123,7 +123,7 @@ void Foam::fv::solidEquilibriumEnergySource::addSup
     const volScalarField& rho,
     fvMatrix<scalar>& eqn,
     const label fieldi
-)
+) const
 {
     const volScalarField alphahe(thermo().alphahe());
 
@@ -147,7 +147,7 @@ void Foam::fv::solidEquilibriumEnergySource::addSup
     const volScalarField& rho,
     fvMatrix<scalar>& eqn,
     const label fieldi
-)
+) const
 {
     const volScalarField alphahe(alpha*thermo().alphahe());
 

src/fvOptions/sources/derived/solidEquilibriumEnergySource/solidEquilibriumEnergySource.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2019-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -125,7 +125,7 @@ public:
                 const volScalarField&,
                 fvMatrix<scalar>&,
                 const label
-            );
+            ) const;
 
             //- Explicit and implicit sources for phase equations
             virtual void addSup
@@ -134,7 +134,7 @@ public:
                 const volScalarField&,
                 fvMatrix<scalar>&,
                 const label
-            );
+            ) const;
 
 
         // IO

src/fvOptions/sources/derived/solidificationMeltingSource/solidificationMeltingSource.C

@@ -131,7 +131,10 @@ Foam::vector Foam::fv::solidificationMeltingSource::g() const
 }
 
 
-void Foam::fv::solidificationMeltingSource::update(const volScalarField& Cp)
+void Foam::fv::solidificationMeltingSource::update
+(
+    const volScalarField& Cp
+) const
 {
     if (curTimeIndex_ == mesh_.time().timeIndex())
     {
@@ -264,7 +267,7 @@ void Foam::fv::solidificationMeltingSource::addSup
 (
     fvMatrix<scalar>& eqn,
     const label fieldi
-)
+) const
 {
     apply(geometricOneField(), eqn);
 }
@@ -275,7 +278,7 @@ void Foam::fv::solidificationMeltingSource::addSup
     const volScalarField& rho,
     fvMatrix<scalar>& eqn,
     const label fieldi
-)
+) const
 {
     apply(rho, eqn);
 }
@@ -285,7 +288,7 @@ void Foam::fv::solidificationMeltingSource::addSup
 (
     fvMatrix<vector>& eqn,
     const label fieldi
-)
+) const
 {
     if (debug)
     {
@@ -325,7 +328,7 @@ void Foam::fv::solidificationMeltingSource::addSup
     const volScalarField& rho,
     fvMatrix<vector>& eqn,
     const label fieldi
-)
+) const
 {
     // Momentum source uses a Boussinesq approximation - redirect
     addSup(eqn, fieldi);

src/fvOptions/sources/derived/solidificationMeltingSource/solidificationMeltingSource.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2014-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -203,13 +203,13 @@ private:
         scalar beta_;
 
         //- Phase fraction indicator field
-        volScalarField alpha1_;
+        mutable volScalarField alpha1_;
 
         //- Current time index (used for updating)
-        label curTimeIndex_;
+        mutable label curTimeIndex_;
 
         //- Temperature change cached for source calculation when alpha1 updated
-        scalarField deltaT_;
+        mutable scalarField deltaT_;
 
 
     // Private Member Functions
@@ -221,11 +221,11 @@ private:
         vector g() const;
 
         //- Update the model
-        void update(const volScalarField& Cp);
+        void update(const volScalarField& Cp) const;
 
         //- Helper function to apply to the energy equation
         template<class RhoFieldType>
-        void apply(const RhoFieldType& rho, fvMatrix<scalar>& eqn);
+        void apply(const RhoFieldType& rho, fvMatrix<scalar>& eqn) const;
 
 
 public:
@@ -257,10 +257,18 @@ public:
         // Add explicit and implicit contributions
 
             //- Add explicit contribution to enthalpy equation
-            virtual void addSup(fvMatrix<scalar>& eqn, const label fieldi);
+            virtual void addSup
+            (
+                fvMatrix<scalar>& eqn,
+                const label fieldi
+            ) const;
 
             //- Add implicit contribution to momentum equation
-            virtual void addSup(fvMatrix<vector>& eqn, const label fieldi);
+            virtual void addSup
+            (
+                fvMatrix<vector>& eqn,
+                const label fieldi
+            ) const;
 
 
         // Add explicit and implicit contributions to compressible equation
@@ -271,7 +279,7 @@ public:
                 const volScalarField& rho,
                 fvMatrix<scalar>& eqn,
                 const label fieldi
-            );
+            ) const;
 
             //- Add implicit contribution to compressible momentum equation
             virtual void addSup
@@ -279,7 +287,7 @@ public:
                 const volScalarField& rho,
                 fvMatrix<vector>& eqn,
                 const label fieldi
-            );
+            ) const;
 
 
         // IO

src/fvOptions/sources/derived/solidificationMeltingSource/solidificationMeltingSourceTemplates.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2014-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -33,7 +33,7 @@ void Foam::fv::solidificationMeltingSource::apply
 (
     const RhoFieldType& rho,
     fvMatrix<scalar>& eqn
-)
+) const
 {
     if (debug)
     {

src/fvOptions/sources/derived/tabulatedAccelerationSource/tabulatedAccelerationSource.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2015-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -77,7 +77,7 @@ void Foam::fv::tabulatedAccelerationSource::addSup
 (
     fvMatrix<vector>& eqn,
     const label fieldi
-)
+) const
 {
     addSup<geometricOneField>(geometricOneField(), eqn, fieldi);
 }
@@ -88,7 +88,7 @@ void Foam::fv::tabulatedAccelerationSource::addSup
     const volScalarField& rho,
     fvMatrix<vector>& eqn,
     const label fieldi
-)
+) const
 {
     addSup<volScalarField>(rho, eqn, fieldi);
 }

src/fvOptions/sources/derived/tabulatedAccelerationSource/tabulatedAccelerationSource.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2015-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -79,6 +79,7 @@ protected:
 
         dimensionedVector g0_;
 
+
 private:
 
     // Private Member Functions
@@ -90,7 +91,7 @@ private:
             const RhoFieldType& rho,
             fvMatrix<vector>& eqn,
             const label fieldi
-        );
+        ) const;
 
 
 public:
@@ -129,7 +130,7 @@ public:
         (
             fvMatrix<vector>& eqn,
             const label fieldi
-        );
+        ) const;
 
         //- Source term to compressible momentum equation
         virtual void addSup
@@ -137,7 +138,7 @@ public:
             const volScalarField& rho,
             fvMatrix<vector>& eqn,
             const label fieldi
-        );
+        ) const;
 
         //- Read dictionary
         virtual bool read(const dictionary& dict);

src/fvOptions/sources/derived/tabulatedAccelerationSource/tabulatedAccelerationSourceTemplates.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2015-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -36,7 +36,7 @@ void Foam::fv::tabulatedAccelerationSource::addSup
     const RhoFieldType& rho,
     fvMatrix<vector>& eqn,
     const label fieldi
-)
+) const
 {
     Vector<vector> acceleration(motion_.acceleration());
 

src/fvOptions/sources/derived/volumeFractionSource/volumeFractionSource.C

@@ -125,7 +125,7 @@ void Foam::fv::volumeFractionSource::addDivSup
 (
     fvMatrix<Type>& eqn,
     const label fieldi
-)
+) const
 {
     const surfaceScalarField& phi =
         mesh().lookupObject<surfaceScalarField>(phiName_);
@@ -140,7 +140,7 @@ void Foam::fv::volumeFractionSource::addUDivSup
 (
     fvMatrix<vector>& eqn,
     const label fieldi
-)
+) const
 {
     const surfaceScalarField& phi =
         mesh().lookupObject<surfaceScalarField>(phiName_);
@@ -157,7 +157,7 @@ void Foam::fv::volumeFractionSource::addRhoDivSup
 (
     fvMatrix<scalar>& eqn,
     const label fieldi
-)
+) const
 {
     const surfaceScalarField& phi =
         mesh().lookupObject<surfaceScalarField>(phiName_);
@@ -174,7 +174,7 @@ void Foam::fv::volumeFractionSource::addLaplacianSup
     const AlphaFieldType& alpha,
     fvMatrix<Type>& eqn,
     const label fieldi
-)
+) const
 {
     const volScalarField B(1 - this->alpha());
 
@@ -220,7 +220,7 @@ void Foam::fv::volumeFractionSource::addSup
 (
     fvMatrix<scalar>& eqn,
     const label fieldi
-)
+) const
 {
     if (fieldNames_[fieldi] == rhoName_)
     {
@@ -238,7 +238,7 @@ void Foam::fv::volumeFractionSource::addSup
 (
     fvMatrix<vector>& eqn,
     const label fieldi
-)
+) const
 {
     if (fieldNames_[fieldi] == UName_)
     {
@@ -256,7 +256,7 @@ void Foam::fv::volumeFractionSource::addSup
 (
     fvMatrix<sphericalTensor>& eqn,
     const label fieldi
-)
+) const
 {
     addDivSup(eqn, fieldi);
     addLaplacianSup(geometricOneField(), eqn, fieldi);
@@ -267,7 +267,7 @@ void Foam::fv::volumeFractionSource::addSup
 (
     fvMatrix<symmTensor>& eqn,
     const label fieldi
-)
+) const
 {
     addDivSup(eqn, fieldi);
     addLaplacianSup(geometricOneField(), eqn, fieldi);
@@ -278,7 +278,7 @@ void Foam::fv::volumeFractionSource::addSup
 (
     fvMatrix<tensor>& eqn,
     const label fieldi
-)
+) const
 {
     addDivSup(eqn, fieldi);
     addLaplacianSup(geometricOneField(), eqn, fieldi);
@@ -290,7 +290,7 @@ void Foam::fv::volumeFractionSource::addSup
     const volScalarField& rho,
     fvMatrix<scalar>& eqn,
     const label fieldi
-)
+) const
 {
     if (fieldNames_[fieldi] == rhoName_)
     {
@@ -309,7 +309,7 @@ void Foam::fv::volumeFractionSource::addSup
     const volScalarField& rho,
     fvMatrix<vector>& eqn,
     const label fieldi
-)
+) const
 {
     if (fieldNames_[fieldi] == UName_)
     {
@@ -328,7 +328,7 @@ void Foam::fv::volumeFractionSource::addSup
     const volScalarField& rho,
     fvMatrix<sphericalTensor>& eqn,
     const label fieldi
-)
+) const
 {
     addDivSup(eqn, fieldi);
     addLaplacianSup(geometricOneField(), eqn, fieldi);
@@ -340,7 +340,7 @@ void Foam::fv::volumeFractionSource::addSup
     const volScalarField& rho,
     fvMatrix<symmTensor>& eqn,
     const label fieldi
-)
+) const
 {
     addDivSup(eqn, fieldi);
     addLaplacianSup(geometricOneField(), eqn, fieldi);
@@ -352,7 +352,7 @@ void Foam::fv::volumeFractionSource::addSup
     const volScalarField& rho,
     fvMatrix<tensor>& eqn,
     const label fieldi
-)
+) const
 {
     addDivSup(eqn, fieldi);
     addLaplacianSup(geometricOneField(), eqn, fieldi);
@@ -365,7 +365,7 @@ void Foam::fv::volumeFractionSource::addSup
     const volScalarField& rho,
     fvMatrix<scalar>& eqn,
     const label fieldi
-)
+) const
 {
     if (fieldNames_[fieldi] == rhoName_)
     {
@@ -385,7 +385,7 @@ void Foam::fv::volumeFractionSource::addSup
     const volScalarField& rho,
     fvMatrix<vector>& eqn,
     const label fieldi
-)
+) const
 {
     if (fieldNames_[fieldi] == UName_)
     {
@@ -405,7 +405,7 @@ void Foam::fv::volumeFractionSource::addSup
     const volScalarField& rho,
     fvMatrix<sphericalTensor>& eqn,
     const label fieldi
-)
+) const
 {
     addDivSup(eqn, fieldi);
     addLaplacianSup(alpha, eqn, fieldi);
@@ -418,7 +418,7 @@ void Foam::fv::volumeFractionSource::addSup
     const volScalarField& rho,
     fvMatrix<symmTensor>& eqn,
     const label fieldi
-)
+) const
 {
     addDivSup(eqn, fieldi);
     addLaplacianSup(alpha, eqn, fieldi);
@@ -431,7 +431,7 @@ void Foam::fv::volumeFractionSource::addSup
     const volScalarField& rho,
     fvMatrix<tensor>& eqn,
     const label fieldi
-)
+) const
 {
     addDivSup(eqn, fieldi);
     addLaplacianSup(alpha, eqn, fieldi);

src/fvOptions/sources/derived/volumeFractionSource/volumeFractionSource.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2019-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -115,11 +115,16 @@ class volumeFractionSource
 
         //- Add time-derivative terms to an equation
         template <class Type>
-        void addDdtSup(fvMatrix<Type>&, const label);
+        void addDdtSup(fvMatrix<Type>&, const label) const;
 
         //- Add time-derivative terms to a compressible equation
         template <class Type>
-        void addDdtSup(const volScalarField&, fvMatrix<Type>&, const label);
+        void addDdtSup
+        (
+            const volScalarField&,
+            fvMatrix<Type>&,
+            const label
+        ) const;
 
         //- Add time-derivative terms to a phase-compressible equation
         template <class Type>
@@ -129,17 +134,17 @@ class volumeFractionSource
             const volScalarField&,
             fvMatrix<Type>&,
             const label
-        );
+        ) const;
 
         //- Add divergence terms to an equation
         template <class Type>
-        void addDivSup(fvMatrix<Type>&, const label);
+        void addDivSup(fvMatrix<Type>&, const label) const;
 
         //- Add divergence terms to the momentum equation
-        void addUDivSup(fvMatrix<vector>&, const label);
+        void addUDivSup(fvMatrix<vector>&, const label) const;
 
         //- Add divergence terms to the continuity equation
-        void addRhoDivSup(fvMatrix<scalar>&, const label);
+        void addRhoDivSup(fvMatrix<scalar>&, const label) const;
 
         //- Add Laplacian terms to an equation
         template <class Type, class AlphaFieldType>
@@ -148,7 +153,7 @@ class volumeFractionSource
             const AlphaFieldType& alpha,
             fvMatrix<Type>& eqn,
             const label fieldi
-        );
+        ) const;
 
 
 public:
@@ -182,15 +187,19 @@ public:
 
             // Explicit and implicit sources
 
-                virtual void addSup(fvMatrix<scalar>&, const label);
+                virtual void addSup(fvMatrix<scalar>&, const label) const;
 
-                virtual void addSup(fvMatrix<vector>&, const label);
+                virtual void addSup(fvMatrix<vector>&, const label) const;
 
-                virtual void addSup(fvMatrix<symmTensor>&, const label);
+                virtual void addSup(fvMatrix<symmTensor>&, const label) const;
 
-                virtual void addSup(fvMatrix<sphericalTensor>&, const label);
+                virtual void addSup
+                (
+                    fvMatrix<sphericalTensor>&,
+                    const label
+                ) const;
 
-                virtual void addSup(fvMatrix<tensor>&, const label);
+                virtual void addSup(fvMatrix<tensor>&, const label) const;
 
 
             // Explicit and implicit sources for compressible equations
@@ -200,35 +209,35 @@ public:
                     const volScalarField&,
                     fvMatrix<scalar>&,
                     const label
-                );
+                ) const;
 
                 virtual void addSup
                 (
                     const volScalarField&,
                     fvMatrix<vector>&,
                     const label
-                );
+                ) const;
 
                 virtual void addSup
                 (
                     const volScalarField&,
                     fvMatrix<symmTensor>&,
                     const label
-                );
+                ) const;
 
                 virtual void addSup
                 (
                     const volScalarField&,
                     fvMatrix<sphericalTensor>&,
                     const label
-                );
+                ) const;
 
                 virtual void addSup
                 (
                     const volScalarField&,
                     fvMatrix<tensor>&,
                     const label
-                );
+                ) const;
 
 
             // Explicit and implicit sources for phase equations
@@ -239,7 +248,7 @@ public:
                     const volScalarField&,
                     fvMatrix<scalar>&,
                     const label
-                );
+                ) const;
 
                 virtual void addSup
                 (
@@ -247,7 +256,7 @@ public:
                     const volScalarField&,
                     fvMatrix<vector>&,
                     const label
-                );
+                ) const;
 
                 virtual void addSup
                 (
@@ -255,7 +264,7 @@ public:
                     const volScalarField&,
                     fvMatrix<symmTensor>&,
                     const label
-                );
+                ) const;
 
                 virtual void addSup
                 (
@@ -263,7 +272,7 @@ public:
                     const volScalarField&,
                     fvMatrix<sphericalTensor>&,
                     const label
-                );
+                ) const;
 
                 virtual void addSup
                 (
@@ -271,7 +280,7 @@ public:
                     const volScalarField&,
                     fvMatrix<tensor>&,
                     const label
-                );
+                ) const;
 
 
         // IO

src/fvOptions/sources/general/codedSource/CodedSource.C

@@ -158,7 +158,7 @@ template<class Type>
 void Foam::fv::CodedSource<Type>::correct
 (
     GeometricField<Type, fvPatchField, volMesh>& field
-)
+) const
 {
     if (debug)
     {
@@ -176,7 +176,7 @@ void Foam::fv::CodedSource<Type>::addSup
 (
     fvMatrix<Type>& eqn,
     const label fieldi
-)
+) const
 {
     if (debug)
     {
@@ -195,7 +195,7 @@ void Foam::fv::CodedSource<Type>::addSup
     const volScalarField& rho,
     fvMatrix<Type>& eqn,
     const label fieldi
-)
+) const
 {
     if (debug)
     {
@@ -213,7 +213,7 @@ void Foam::fv::CodedSource<Type>::constrain
 (
     fvMatrix<Type>& eqn,
     const label fieldi
-)
+) const
 {
     if (debug)
     {

src/fvOptions/sources/general/codedSource/CodedSource.H

@@ -191,20 +191,21 @@ public:
         //- Dynamically compiled fvOption
         option& redirectFvOption() const;
 
+
         // Evaluation
 
             //- Correct field
             virtual void correct
             (
                 GeometricField<Type, fvPatchField, volMesh>&
-            );
+            ) const;
 
             //- Explicit and implicit matrix contributions
             virtual void addSup
             (
                 fvMatrix<Type>& eqn,
                 const label fieldi
-            );
+            ) const;
 
             //- Explicit and implicit matrix contributions
             //  to compressible equation
@@ -213,14 +214,14 @@ public:
                 const volScalarField& rho,
                 fvMatrix<Type>& eqn,
                 const label fieldi
-            );
+            ) const;
 
             //- Set value
             virtual void constrain
             (
                 fvMatrix<Type>& eqn,
                 const label fieldi
-            );
+            ) const;
 
 
         // IO

src/fvOptions/sources/general/semiImplicitSource/semiImplicitSource.C

@@ -64,7 +64,7 @@ void Foam::fv::semiImplicitSource::addSupType
 (
     fvMatrix<Type>& eqn,
     const label fieldi
-)
+) const
 {
     if (debug)
     {
@@ -132,7 +132,7 @@ void Foam::fv::semiImplicitSource::addSupType
     const volScalarField& rho,
     fvMatrix<Type>& eqn,
     const label fieldi
-)
+) const
 {
     if (debug)
     {
@@ -174,7 +174,7 @@ void Foam::fv::semiImplicitSource::addSup
 (
     fvMatrix<scalar>& eqn,
     const label fieldi
-)
+) const
 {
     addSupType(eqn, fieldi);
 }
@@ -184,7 +184,7 @@ void Foam::fv::semiImplicitSource::addSup
 (
     fvMatrix<vector>& eqn,
     const label fieldi
-)
+) const
 {
     addSupType(eqn, fieldi);
 }
@@ -194,7 +194,7 @@ void Foam::fv::semiImplicitSource::addSup
 (
     fvMatrix<symmTensor>& eqn,
     const label fieldi
-)
+) const
 {
     addSupType(eqn, fieldi);
 }
@@ -204,7 +204,7 @@ void Foam::fv::semiImplicitSource::addSup
 (
     fvMatrix<sphericalTensor>& eqn,
     const label fieldi
-)
+) const
 {
     addSupType(eqn, fieldi);
 }
@@ -214,7 +214,7 @@ void Foam::fv::semiImplicitSource::addSup
 (
     fvMatrix<tensor>& eqn,
     const label fieldi
-)
+) const
 {
     addSupType(eqn, fieldi);
 }
@@ -225,7 +225,7 @@ void Foam::fv::semiImplicitSource::addSup
     const volScalarField& rho,
     fvMatrix<scalar>& eqn,
     const label fieldi
-)
+) const
 {
     addSupType(eqn, fieldi);
 }
@@ -236,7 +236,7 @@ void Foam::fv::semiImplicitSource::addSup
     const volScalarField& rho,
     fvMatrix<vector>& eqn,
     const label fieldi
-)
+) const
 {
     addSupType(eqn, fieldi);
 }
@@ -247,7 +247,7 @@ void Foam::fv::semiImplicitSource::addSup
     const volScalarField& rho,
     fvMatrix<symmTensor>& eqn,
     const label fieldi
-)
+) const
 {
     addSupType(eqn, fieldi);
 }
@@ -258,7 +258,7 @@ void Foam::fv::semiImplicitSource::addSup
     const volScalarField& rho,
     fvMatrix<sphericalTensor>& eqn,
     const label fieldi
-)
+) const
 {
     addSupType(eqn, fieldi);
 }
@@ -269,7 +269,7 @@ void Foam::fv::semiImplicitSource::addSup
     const volScalarField& rho,
     fvMatrix<tensor>& eqn,
     const label fieldi
-)
+) const
 {
     addSupType(eqn, fieldi);
 }
@@ -281,7 +281,7 @@ void Foam::fv::semiImplicitSource::addSup
     const volScalarField& rho,
     fvMatrix<scalar>& eqn,
     const label fieldi
-)
+) const
 {
     addSupType(eqn, fieldi);
 }
@@ -293,7 +293,7 @@ void Foam::fv::semiImplicitSource::addSup
     const volScalarField& rho,
     fvMatrix<vector>& eqn,
     const label fieldi
-)
+) const
 {
     addSupType(eqn, fieldi);
 }
@@ -305,7 +305,7 @@ void Foam::fv::semiImplicitSource::addSup
     const volScalarField& rho,
     fvMatrix<symmTensor>& eqn,
     const label fieldi
-)
+) const
 {
     addSupType(eqn, fieldi);
 }
@@ -317,7 +317,7 @@ void Foam::fv::semiImplicitSource::addSup
     const volScalarField& rho,
     fvMatrix<sphericalTensor>& eqn,
     const label fieldi
-)
+) const
 {
     addSupType(eqn, fieldi);
 }
@@ -329,7 +329,7 @@ void Foam::fv::semiImplicitSource::addSup
     const volScalarField& rho,
     fvMatrix<tensor>& eqn,
     const label fieldi
-)
+) const
 {
     addSupType(eqn, fieldi);
 }
@@ -393,4 +393,5 @@ bool Foam::fv::semiImplicitSource::read(const dictionary& dict)
     }
 }
 
+
 // ************************************************************************* //

src/fvOptions/sources/general/semiImplicitSource/semiImplicitSource.H

@@ -147,11 +147,16 @@ private:
 
         //- Add divergence terms to an equation
         template <class Type>
-        void addSupType(fvMatrix<Type>&, const label);
+        void addSupType(fvMatrix<Type>&, const label) const;
 
         //- Add divergence terms to an equation
         template <class Type>
-        void addSupType(const volScalarField&, fvMatrix<Type>&, const label);
+        void addSupType
+        (
+            const volScalarField&,
+            fvMatrix<Type>&,
+            const label
+        ) const;
 
 
 public:
@@ -185,31 +190,31 @@ public:
                 (
                     fvMatrix<scalar>& eqn,
                     const label fieldi
-                );
+                ) const;
 
                 virtual void addSup
                 (
                     fvMatrix<vector>& eqn,
                     const label fieldi
-                );
+                ) const;
 
                 virtual void addSup
                 (
                     fvMatrix<symmTensor>& eqn,
                     const label fieldi
-                );
+                ) const;
 
                 virtual void addSup
                 (
                     fvMatrix<sphericalTensor>& eqn,
                     const label fieldi
-                );
+                ) const;
 
                 virtual void addSup
                 (
                     fvMatrix<tensor>& eqn,
                     const label fieldi
-                );
+                ) const;
 
 
             // Explicit and implicit sources for compressible equations
@@ -219,35 +224,35 @@ public:
                     const volScalarField& rho,
                     fvMatrix<scalar>& eqn,
                     const label fieldi
-                );
+                ) const;
 
                 virtual void addSup
                 (
                     const volScalarField& rho,
                     fvMatrix<vector>& eqn,
                     const label fieldi
-                );
+                ) const;
 
                 virtual void addSup
                 (
                     const volScalarField& rho,
                     fvMatrix<symmTensor>& eqn,
                     const label fieldi
-                );
+                ) const;
 
                 virtual void addSup
                 (
                     const volScalarField& rho,
                     fvMatrix<sphericalTensor>& eqn,
                     const label fieldi
-                );
+                ) const;
 
                 virtual void addSup
                 (
                     const volScalarField& rho,
                     fvMatrix<tensor>& eqn,
                     const label fieldi
-                );
+                ) const;
 
 
             // Explicit and implicit sources for phase equations
@@ -258,7 +263,7 @@ public:
                     const volScalarField& rho,
                     fvMatrix<scalar>& eqn,
                     const label fieldi
-                );
+                ) const;
 
                 virtual void addSup
                 (
@@ -266,7 +271,7 @@ public:
                     const volScalarField& rho,
                     fvMatrix<vector>& eqn,
                     const label fieldi
-                );
+                ) const;
 
                 virtual void addSup
                 (
@@ -274,7 +279,7 @@ public:
                     const volScalarField& rho,
                     fvMatrix<symmTensor>& eqn,
                     const label fieldi
-                );
+                ) const;
 
                 virtual void addSup
                 (
@@ -282,7 +287,7 @@ public:
                     const volScalarField& rho,
                     fvMatrix<sphericalTensor>& eqn,
                     const label fieldi
-                );
+                ) const;
 
                 virtual void addSup
                 (
@@ -290,7 +295,7 @@ public:
                     const volScalarField& rho,
                     fvMatrix<tensor>& eqn,
                     const label fieldi
-                );
+                ) const;
 
 
         // IO

src/fvOptions/sources/interRegion/interRegionExplicitPorositySource/interRegionExplicitPorositySource.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2012-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -47,7 +47,7 @@ namespace fv
 
 // * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
 
-void Foam::fv::interRegionExplicitPorositySource::initialise()
+void Foam::fv::interRegionExplicitPorositySource::initialise() const
 {
     if (!firstIter_)
     {
@@ -135,7 +135,7 @@ void Foam::fv::interRegionExplicitPorositySource::addSup
 (
     fvMatrix<vector>& eqn,
     const label fieldi
-)
+) const
 {
     initialise();
 
@@ -189,7 +189,7 @@ void Foam::fv::interRegionExplicitPorositySource::addSup
     const volScalarField& rho,
     fvMatrix<vector>& eqn,
     const label fieldi
-)
+) const
 {
     initialise();
 

src/fvOptions/sources/interRegion/interRegionExplicitPorositySource/interRegionExplicitPorositySource.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2012-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -85,10 +85,10 @@ protected:
     // Protected data
 
         //- Run-time selectable porosity model
-        autoPtr<porosityModel> porosityPtr_;
+        mutable autoPtr<porosityModel> porosityPtr_;
 
         //- First iteration
-        bool firstIter_;
+        mutable bool firstIter_;
 
         //- Velocity field name, default = U
         word UName_;
@@ -101,7 +101,7 @@ protected:
     // Protected Member Functions
 
         //- Initialise
-        void initialise();
+        void initialise() const;
 
 
 public:
@@ -142,7 +142,7 @@ public:
             (
                 fvMatrix<vector>& eqn,
                 const label fieldi
-            );
+            ) const;
 
 
         // Add explicit and implicit contributions to compressible equation
@@ -153,7 +153,7 @@ public:
                 const volScalarField& rho,
                 fvMatrix<vector>& eqn,
                 const label fieldi
-            );
+            ) const;
 
 
         // IO

src/fvOptions/sources/interRegion/interRegionHeatTransfer/constantHeatTransfer/constantHeatTransfer.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -104,7 +104,7 @@ Foam::fv::constantHeatTransfer::~constantHeatTransfer()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::fv::constantHeatTransfer::calculateHtc()
+void Foam::fv::constantHeatTransfer::calculateHtc() const
 {}
 
 

src/fvOptions/sources/interRegion/interRegionHeatTransfer/constantHeatTransfer/constantHeatTransfer.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -85,7 +85,7 @@ public:
     // Public Functions
 
         //- Calculate the heat transfer coefficient
-        virtual void calculateHtc();
+        virtual void calculateHtc() const;
 
 
         // IO

src/fvOptions/sources/interRegion/interRegionHeatTransfer/interRegionHeatTransferModel/interRegionHeatTransferModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -43,7 +43,7 @@ namespace fv
 
 // * * * * * * * * * * * *  Protected member functions * * * * * * * * * * * //
 
-void Foam::fv::interRegionHeatTransferModel::setNbrModel()
+void Foam::fv::interRegionHeatTransferModel::setNbrModel() const
 {
     if (!firstIter_)
     {
@@ -84,7 +84,7 @@ void Foam::fv::interRegionHeatTransferModel::setNbrModel()
 }
 
 
-void Foam::fv::interRegionHeatTransferModel::correct()
+void Foam::fv::interRegionHeatTransferModel::correct() const
 {
     if (master_)
     {
@@ -167,7 +167,7 @@ void Foam::fv::interRegionHeatTransferModel::addSup
 (
     fvMatrix<scalar>& eqn,
     const label fieldi
-)
+) const
 {
     setNbrModel();
 
@@ -254,7 +254,7 @@ void Foam::fv::interRegionHeatTransferModel::addSup
     const volScalarField& rho,
     fvMatrix<scalar>& eqn,
     const label fieldi
-)
+) const
 {
     addSup(eqn, fieldi);
 }

src/fvOptions/sources/interRegion/interRegionHeatTransfer/interRegionHeatTransferModel/interRegionHeatTransferModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -64,16 +64,16 @@ protected:
         word nbrModelName_;
 
         //- Pointer to neighbour interRegionHeatTransferModel
-        interRegionHeatTransferModel* nbrModel_;
+        mutable interRegionHeatTransferModel* nbrModel_;
 
         //- First iteration
-        bool firstIter_;
+        mutable bool firstIter_;
 
         //- Time index - used for updating htc
-        label timeIndex_;
+        mutable label timeIndex_;
 
         //- Heat transfer coefficient [W/m^2/k] times area/volume [1/m]
-        volScalarField htc_;
+        mutable volScalarField htc_;
 
         //- Flag to activate semi-implicit coupling
         bool semiImplicit_;
@@ -88,13 +88,13 @@ protected:
     // Protected member functions
 
         //- Set the neighbour interRegionHeatTransferModel
-        void setNbrModel();
+        void setNbrModel() const;
 
         //- Inherit correct from interRegionOption
         using interRegionOption::correct;
 
         //- Correct to calculate the inter-region heat transfer coefficient
-        void correct();
+        void correct() const;
 
         //- Interpolate field with nbrModel specified
         template<class Type>
@@ -175,7 +175,7 @@ public:
             (
                 fvMatrix<scalar>& eqn,
                 const label fieldi
-            );
+            ) const;
 
             //- Source term to compressible energy equation
             virtual void addSup
@@ -183,10 +183,10 @@ public:
                 const volScalarField& rho,
                 fvMatrix<scalar>& eqn,
                 const label fieldi
-            );
+            ) const;
 
             //- Calculate heat transfer coefficient
-            virtual void calculateHtc() = 0;
+            virtual void calculateHtc() const = 0;
 
 
         // IO

src/fvOptions/sources/interRegion/interRegionHeatTransfer/tabulatedHeatTransfer/tabulatedHeatTransfer.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -46,7 +46,7 @@ namespace fv
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 const Foam::interpolation2DTable<Foam::scalar>&
-Foam::fv::tabulatedHeatTransfer::hTable()
+Foam::fv::tabulatedHeatTransfer::hTable() const
 {
     if (!hTable_.valid())
     {
@@ -57,7 +57,7 @@ Foam::fv::tabulatedHeatTransfer::hTable()
 }
 
 
-const Foam::volScalarField& Foam::fv::tabulatedHeatTransfer::AoV()
+const Foam::volScalarField& Foam::fv::tabulatedHeatTransfer::AoV() const
 {
     if (!AoV_.valid())
     {
@@ -109,7 +109,7 @@ Foam::fv::tabulatedHeatTransfer::~tabulatedHeatTransfer()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::fv::tabulatedHeatTransfer::calculateHtc()
+void Foam::fv::tabulatedHeatTransfer::calculateHtc() const
 {
     const fvMesh& nbrMesh = mesh_.time().lookupObject<fvMesh>(nbrRegionName());