From 27c62303ad4b8b90ece4fe41089926d59aec3e06 Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Mon, 7 Jan 2019 09:20:51 +0100
Subject: [PATCH] ENH: for-range, forAllIters() ... in applications/solvers

- reduced clutter when iterating over containers
---
 .../multiphaseMixtureThermo.C                 | 201 +++++------
 .../compressibleMultiphaseInterFoam/pEqn.H    |  19 +-
 .../icoReactingMultiphaseInterFoam/YEqns.H    |   4 +-
 .../laserDTRM/laserDTRM.C                     |  22 +-
 .../interfaceCompositionModel.H               |   1 -
 .../MassTransferPhaseSystem.C                 |  21 +-
 .../MultiComponentPhaseModel.C                |   2 +-
 .../multiphaseSystem/multiphaseSystem.C       |  79 ++---
 .../phasesSystem/phaseSystem/phaseSystem.C    | 329 +++++++++---------
 .../phaseSystem/phaseSystemTemplates.H        |  25 +-
 .../twoPhaseMixtureEThermo.C                  |   8 +-
 .../twoPhaseMixtureEThermo.H                  |   6 +-
 .../multiphaseEulerFoam/CourantNo.H           |   4 +-
 .../multiphase/multiphaseEulerFoam/DDtU.H     |   4 +-
 .../multiphaseEulerFoam/MRFCorrectBCs.H       |   4 +-
 .../multiphase/multiphaseEulerFoam/UEqns.H    |   9 +-
 .../multiphaseEulerFoam/createFields.H        |   3 +-
 .../multiphaseEulerFoam/createMRFZones.H      |   4 +-
 .../multiphaseSystem/multiphaseSystem.C       | 236 ++++++-------
 .../multiphaseSystem/phaseModel/phaseModel.C  |   6 +-
 .../multiphase/multiphaseEulerFoam/pEqn.H     |  46 +--
 .../multiphaseEulerFoam/zonePhaseVolumes.H    |   6 +-
 .../multiphaseMixture/multiphaseMixture.C     | 100 +++---
 .../interfaceCompositionModels/Henry/Henry.C  |  11 +-
 .../interfaceCompositionModels/Henry/Henry.H  |   2 +-
 .../InterfaceCompositionModel.C               |  22 +-
 .../InterfaceCompositionModel.H               |   6 +-
 .../Raoult/Raoult.C                           |  20 +-
 .../Raoult/Raoult.H                           |   2 +-
 .../interfaceCompositionModel.C               |   6 -
 .../interfaceCompositionModel.H               |   4 +-
 .../blendingMethods/hyperbolic/hyperbolic.C   |  19 +-
 .../blendingMethods/hyperbolic/hyperbolic.H   |   2 +-
 .../blendingMethods/linear/linear.C           |  42 +--
 .../blendingMethods/linear/linear.H           |   2 +-
 .../HeatAndMassTransferPhaseSystem.C          |  12 +-
 .../HeatTransferPhaseSystem.C                 |  12 +-
 ...terfaceCompositionPhaseChangePhaseSystem.C |  22 +-
 .../MomentumTransferPhaseSystem.C             |  12 +-
 .../ThermalPhaseChangePhaseSystem.C           |  16 +-
 .../IsothermalPhaseModel.C                    |   7 -
 .../IsothermalPhaseModel.H                    |   2 +-
 .../phaseSystems/phaseSystem/phaseSystem.C    |  54 +--
 .../phaseSystem/phaseSystemTemplates.C        |  14 +-
 .../multiphaseSystem/multiphaseSystem.C       |  94 +++--
 .../multiphaseSystem/multiphaseSystem.H       |   4 +-
 .../blendingMethods/hyperbolic/hyperbolic.C   |  19 +-
 .../blendingMethods/hyperbolic/hyperbolic.H   |   2 +-
 .../blendingMethods/linear/linear.C           |  42 +--
 .../blendingMethods/linear/linear.H           |   2 +-
 .../twoPhaseSystem/twoPhaseSystem.C           |   6 +-
 51 files changed, 687 insertions(+), 910 deletions(-)

diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C
index 78c75b9348e..b16588a2768 100644
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C
@@ -51,9 +51,9 @@ void Foam::multiphaseMixtureThermo::calcAlphas()
     scalar level = 0.0;
     alphas_ == 0.0;
 
-    forAllIter(PtrDictionary<phaseModel>, phases_, phase)
+    for (const phaseModel& phase : phases_)
     {
-        alphas_ += level*phase();
+        alphas_ += level * phase;
         level += 1.0;
     }
 }
@@ -121,18 +121,18 @@ Foam::multiphaseMixtureThermo::multiphaseMixtureThermo
 
 void Foam::multiphaseMixtureThermo::correct()
 {
-    forAllIter(PtrDictionary<phaseModel>, phases_, phasei)
+    for (phaseModel& phase : phases_)
     {
-        phasei().correct();
+        phase.correct();
     }
 
-    PtrDictionary<phaseModel>::iterator phasei = phases_.begin();
+    auto phasei = phases_.cbegin();
 
     psi_ = phasei()*phasei().thermo().psi();
     mu_ = phasei()*phasei().thermo().mu();
     alpha_ = phasei()*phasei().thermo().alpha();
 
-    for (++phasei; phasei != phases_.end(); ++phasei)
+    for (++phasei; phasei != phases_.cend(); ++phasei)
     {
         psi_ += phasei()*phasei().thermo().psi();
         mu_ += phasei()*phasei().thermo().mu();
@@ -143,20 +143,20 @@ void Foam::multiphaseMixtureThermo::correct()
 
 void Foam::multiphaseMixtureThermo::correctRho(const volScalarField& dp)
 {
-    forAllIter(PtrDictionary<phaseModel>, phases_, phasei)
+    for (phaseModel& phase : phases_)
     {
-        phasei().thermo().rho() +=  phasei().thermo().psi()*dp;
+        phase.thermo().rho() += phase.thermo().psi()*dp;
     }
 }
 
 
 Foam::word Foam::multiphaseMixtureThermo::thermoName() const
 {
-    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+    auto phasei = phases_.cbegin();
 
     word name = phasei().thermo().thermoName();
 
-    for (++ phasei; phasei != phases_.end(); ++ phasei)
+    for (++phasei; phasei != phases_.cend(); ++phasei)
     {
         name += ',' + phasei().thermo().thermoName();
     }
@@ -167,27 +167,29 @@ Foam::word Foam::multiphaseMixtureThermo::thermoName() const
 
 bool Foam::multiphaseMixtureThermo::incompressible() const
 {
-    bool ico = true;
-
-    forAllConstIter(PtrDictionary<phaseModel>, phases_, phase)
+    for (const phaseModel& phase : phases_)
     {
-        ico &= phase().thermo().incompressible();
+        if (!phase.thermo().incompressible())
+        {
+            return false;
+        }
     }
 
-    return ico;
+    return true;
 }
 
 
 bool Foam::multiphaseMixtureThermo::isochoric() const
 {
-    bool iso = true;
-
-    forAllConstIter(PtrDictionary<phaseModel>, phases_, phase)
+    for (const phaseModel& phase : phases_)
     {
-        iso &= phase().thermo().incompressible();
+        if (!phase.thermo().isochoric())
+        {
+            return false;
+        }
     }
 
-    return iso;
+    return true;
 }
 
 
@@ -197,11 +199,11 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::he
     const volScalarField& T
 ) const
 {
-    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+    auto phasei = phases_.cbegin();
 
     tmp<volScalarField> the(phasei()*phasei().thermo().he(p, T));
 
-    for (++phasei; phasei != phases_.end(); ++phasei)
+    for (++phasei; phasei != phases_.cend(); ++phasei)
     {
         the.ref() += phasei()*phasei().thermo().he(p, T);
     }
@@ -217,14 +219,14 @@ Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::he
     const labelList& cells
 ) const
 {
-    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+    auto phasei = phases_.cbegin();
 
     tmp<scalarField> the
     (
         scalarField(phasei(), cells)*phasei().thermo().he(p, T, cells)
     );
 
-    for (++phasei; phasei != phases_.end(); ++phasei)
+    for (++phasei; phasei != phases_.cend(); ++phasei)
     {
         the.ref() +=
             scalarField(phasei(), cells)*phasei().thermo().he(p, T, cells);
@@ -241,14 +243,14 @@ Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::he
     const label patchi
 ) const
 {
-    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+    auto phasei = phases_.cbegin();
 
     tmp<scalarField> the
     (
         phasei().boundaryField()[patchi]*phasei().thermo().he(p, T, patchi)
     );
 
-    for (++phasei; phasei != phases_.end(); ++phasei)
+    for (++phasei; phasei != phases_.cend(); ++phasei)
     {
         the.ref() +=
             phasei().boundaryField()[patchi]*phasei().thermo().he(p, T, patchi);
@@ -260,11 +262,11 @@ Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::he
 
 Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::hc() const
 {
-    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+    auto phasei = phases_.cbegin();
 
     tmp<volScalarField> thc(phasei()*phasei().thermo().hc());
 
-    for (++phasei; phasei != phases_.end(); ++phasei)
+    for (++phasei; phasei != phases_.cend(); ++phasei)
     {
         thc.ref() += phasei()*phasei().thermo().hc();
     }
@@ -301,11 +303,11 @@ Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::THE
 
 Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::rho() const
 {
-    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+    auto phasei = phases_.cbegin();
 
     tmp<volScalarField> trho(phasei()*phasei().thermo().rho());
 
-    for (++phasei; phasei != phases_.end(); ++phasei)
+    for (++phasei; phasei != phases_.cend(); ++phasei)
     {
         trho.ref() += phasei()*phasei().thermo().rho();
     }
@@ -319,14 +321,14 @@ Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::rho
     const label patchi
 ) const
 {
-    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+    auto phasei = phases_.cbegin();
 
     tmp<scalarField> trho
     (
         phasei().boundaryField()[patchi]*phasei().thermo().rho(patchi)
     );
 
-    for (++phasei; phasei != phases_.end(); ++phasei)
+    for (++phasei; phasei != phases_.cend(); ++phasei)
     {
         trho.ref() +=
             phasei().boundaryField()[patchi]*phasei().thermo().rho(patchi);
@@ -338,11 +340,11 @@ Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::rho
 
 Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::Cp() const
 {
-    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+    auto phasei = phases_.cbegin();
 
     tmp<volScalarField> tCp(phasei()*phasei().thermo().Cp());
 
-    for (++phasei; phasei != phases_.end(); ++phasei)
+    for (++phasei; phasei != phases_.cend(); ++phasei)
     {
         tCp.ref() += phasei()*phasei().thermo().Cp();
     }
@@ -358,14 +360,14 @@ Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::Cp
     const label patchi
 ) const
 {
-    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+    auto phasei = phases_.cbegin();
 
     tmp<scalarField> tCp
     (
         phasei().boundaryField()[patchi]*phasei().thermo().Cp(p, T, patchi)
     );
 
-    for (++phasei; phasei != phases_.end(); ++phasei)
+    for (++phasei; phasei != phases_.cend(); ++phasei)
     {
         tCp.ref() +=
             phasei().boundaryField()[patchi]*phasei().thermo().Cp(p, T, patchi);
@@ -377,11 +379,11 @@ Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::Cp
 
 Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::Cv() const
 {
-    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+    auto phasei = phases_.cbegin();
 
     tmp<volScalarField> tCv(phasei()*phasei().thermo().Cv());
 
-    for (++phasei; phasei != phases_.end(); ++phasei)
+    for (++phasei; phasei != phases_.cend(); ++phasei)
     {
         tCv.ref() += phasei()*phasei().thermo().Cv();
     }
@@ -397,14 +399,14 @@ Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::Cv
     const label patchi
 ) const
 {
-    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+    auto phasei = phases_.cbegin();
 
     tmp<scalarField> tCv
     (
         phasei().boundaryField()[patchi]*phasei().thermo().Cv(p, T, patchi)
     );
 
-    for (++phasei; phasei != phases_.end(); ++phasei)
+    for (++phasei; phasei != phases_.cend(); ++phasei)
     {
         tCv.ref() +=
             phasei().boundaryField()[patchi]*phasei().thermo().Cv(p, T, patchi);
@@ -416,11 +418,11 @@ Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::Cv
 
 Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::gamma() const
 {
-    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+    auto phasei = phases_.cbegin();
 
     tmp<volScalarField> tgamma(phasei()*phasei().thermo().gamma());
 
-    for (++phasei; phasei != phases_.end(); ++phasei)
+    for (++phasei; phasei != phases_.cend(); ++phasei)
     {
         tgamma.ref() += phasei()*phasei().thermo().gamma();
     }
@@ -436,14 +438,14 @@ Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::gamma
     const label patchi
 ) const
 {
-    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+    auto phasei = phases_.cbegin();
 
     tmp<scalarField> tgamma
     (
         phasei().boundaryField()[patchi]*phasei().thermo().gamma(p, T, patchi)
     );
 
-    for (++phasei; phasei != phases_.end(); ++phasei)
+    for (++phasei; phasei != phases_.cend(); ++phasei)
     {
         tgamma.ref() +=
             phasei().boundaryField()[patchi]
@@ -456,11 +458,11 @@ Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::gamma
 
 Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::Cpv() const
 {
-    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+    auto phasei = phases_.cbegin();
 
     tmp<volScalarField> tCpv(phasei()*phasei().thermo().Cpv());
 
-    for (++phasei; phasei != phases_.end(); ++phasei)
+    for (++phasei; phasei != phases_.cend(); ++phasei)
     {
         tCpv.ref() += phasei()*phasei().thermo().Cpv();
     }
@@ -476,14 +478,14 @@ Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::Cpv
     const label patchi
 ) const
 {
-    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+    auto phasei = phases_.cbegin();
 
     tmp<scalarField> tCpv
     (
         phasei().boundaryField()[patchi]*phasei().thermo().Cpv(p, T, patchi)
     );
 
-    for (++phasei; phasei != phases_.end(); ++phasei)
+    for (++phasei; phasei != phases_.cend(); ++phasei)
     {
         tCpv.ref() +=
             phasei().boundaryField()[patchi]
@@ -496,11 +498,11 @@ Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::Cpv
 
 Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::CpByCpv() const
 {
-    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+    auto phasei = phases_.cbegin();
 
     tmp<volScalarField> tCpByCpv(phasei()*phasei().thermo().CpByCpv());
 
-    for (++phasei; phasei != phases_.end(); ++phasei)
+    for (++phasei; phasei != phases_.cend(); ++phasei)
     {
         tCpByCpv.ref() += phasei()*phasei().thermo().CpByCpv();
     }
@@ -516,14 +518,14 @@ Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::CpByCpv
     const label patchi
 ) const
 {
-    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+    auto phasei = phases_.cbegin();
 
     tmp<scalarField> tCpByCpv
     (
         phasei().boundaryField()[patchi]*phasei().thermo().CpByCpv(p, T, patchi)
     );
 
-    for (++phasei; phasei != phases_.end(); ++phasei)
+    for (++phasei; phasei != phases_.cend(); ++phasei)
     {
         tCpByCpv.ref() +=
             phasei().boundaryField()[patchi]
@@ -536,11 +538,11 @@ Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::CpByCpv
 
 Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::W() const
 {
-    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+    auto phasei = phases_.cbegin();
 
     tmp<volScalarField> tW(phasei()*phasei().thermo().W());
 
-    for (++phasei; phasei != phases_.end(); ++phasei)
+    for (++phasei; phasei != phases_.cend(); ++phasei)
     {
         tW.ref() += phasei()*phasei().thermo().W();
     }
@@ -566,11 +568,11 @@ Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::nu
 
 Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::kappa() const
 {
-    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+    auto phasei = phases_.cbegin();
 
     tmp<volScalarField> tkappa(phasei()*phasei().thermo().kappa());
 
-    for (++phasei; phasei != phases_.end(); ++phasei)
+    for (++phasei; phasei != phases_.cend(); ++phasei)
     {
         tkappa.ref() += phasei()*phasei().thermo().kappa();
     }
@@ -584,14 +586,14 @@ Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::kappa
     const label patchi
 ) const
 {
-    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+    auto phasei = phases_.cbegin();
 
     tmp<scalarField> tkappa
     (
         phasei().boundaryField()[patchi]*phasei().thermo().kappa(patchi)
     );
 
-    for (++phasei; phasei != phases_.end(); ++phasei)
+    for (++phasei; phasei != phases_.cend(); ++phasei)
     {
         tkappa.ref() +=
             phasei().boundaryField()[patchi]*phasei().thermo().kappa(patchi);
@@ -606,11 +608,11 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::kappaEff
     const volScalarField& alphat
 ) const
 {
-    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+    auto phasei = phases_.cbegin();
 
     tmp<volScalarField> tkappaEff(phasei()*phasei().thermo().kappaEff(alphat));
 
-    for (++phasei; phasei != phases_.end(); ++phasei)
+    for (++phasei; phasei != phases_.cend(); ++phasei)
     {
         tkappaEff.ref() += phasei()*phasei().thermo().kappaEff(alphat);
     }
@@ -625,7 +627,7 @@ Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::kappaEff
     const label patchi
 ) const
 {
-    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+    auto phasei = phases_.cbegin();
 
     tmp<scalarField> tkappaEff
     (
@@ -633,7 +635,7 @@ Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::kappaEff
        *phasei().thermo().kappaEff(alphat, patchi)
     );
 
-    for (++phasei; phasei != phases_.end(); ++phasei)
+    for (++phasei; phasei != phases_.cend(); ++phasei)
     {
         tkappaEff.ref() +=
             phasei().boundaryField()[patchi]
@@ -649,11 +651,11 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::alphaEff
     const volScalarField& alphat
 ) const
 {
-    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+    auto phasei = phases_.cbegin();
 
     tmp<volScalarField> talphaEff(phasei()*phasei().thermo().alphaEff(alphat));
 
-    for (++phasei; phasei != phases_.end(); ++phasei)
+    for (++phasei; phasei != phases_.cend(); ++phasei)
     {
         talphaEff.ref() += phasei()*phasei().thermo().alphaEff(alphat);
     }
@@ -668,7 +670,7 @@ Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::alphaEff
     const label patchi
 ) const
 {
-    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+    auto phasei = phases_.cbegin();
 
     tmp<scalarField> talphaEff
     (
@@ -676,7 +678,7 @@ Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::alphaEff
        *phasei().thermo().alphaEff(alphat, patchi)
     );
 
-    for (++phasei; phasei != phases_.end(); ++phasei)
+    for (++phasei; phasei != phases_.cend(); ++phasei)
     {
         talphaEff.ref() +=
             phasei().boundaryField()[patchi]
@@ -689,11 +691,11 @@ Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::alphaEff
 
 Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::rCv() const
 {
-    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
+    auto phasei = phases_.cbegin();
 
     tmp<volScalarField> trCv(phasei()/phasei().thermo().Cv());
 
-    for (++phasei; phasei != phases_.end(); ++phasei)
+    for (++phasei; phasei != phases_.cend(); ++phasei)
     {
         trCv.ref() += phasei()/phasei().thermo().Cv();
     }
@@ -723,21 +725,19 @@ Foam::multiphaseMixtureThermo::surfaceTensionForce() const
     surfaceScalarField& stf = tstf.ref();
     stf.setOriented();
 
-    forAllConstIter(PtrDictionary<phaseModel>, phases_, phase1)
+    forAllConstIters(phases_, phase1)
     {
-        const phaseModel& alpha1 = phase1();
+        const phaseModel& alpha1 = *phase1;
 
-        PtrDictionary<phaseModel>::const_iterator phase2 = phase1;
-        ++phase2;
+        auto phase2 = phase1;
 
-        for (; phase2 != phases_.end(); ++phase2)
+        for (++phase2; phase2 != phases_.cend(); ++phase2)
         {
-            const phaseModel& alpha2 = phase2();
+            const phaseModel& alpha2 = *phase2;
 
-            sigmaTable::const_iterator sigma =
-                sigmas_.find(interfacePair(alpha1, alpha2));
+            auto sigma = sigmas_.cfind(interfacePair(alpha1, alpha2));
 
-            if (sigma == sigmas_.end())
+            if (!sigma.found())
             {
                 FatalErrorInFunction
                     << "Cannot find interface " << interfacePair(alpha1, alpha2)
@@ -745,7 +745,7 @@ Foam::multiphaseMixtureThermo::surfaceTensionForce() const
                     << exit(FatalError);
             }
 
-            stf += dimensionedScalar("sigma", dimSigma_, sigma())
+            stf += dimensionedScalar("sigma", dimSigma_, *sigma)
                *fvc::interpolate(K(alpha1, alpha2))*
                 (
                     fvc::interpolate(alpha2)*fvc::snGrad(alpha1)
@@ -862,11 +862,10 @@ void Foam::multiphaseMixtureThermo::correctContactAngle
                /mesh_.magSf().boundaryField()[patchi]
             );
 
-            alphaContactAngleFvPatchScalarField::thetaPropsTable::
-                const_iterator tp =
-                acap.thetaProps().find(interfacePair(alpha1, alpha2));
+            const auto tp =
+                acap.thetaProps().cfind(interfacePair(alpha1, alpha2));
 
-            if (tp == acap.thetaProps().end())
+            if (!tp.found())
             {
                 FatalErrorInFunction
                     << "Cannot find interface " << interfacePair(alpha1, alpha2)
@@ -875,12 +874,12 @@ void Foam::multiphaseMixtureThermo::correctContactAngle
                     << exit(FatalError);
             }
 
-            bool matched = (tp.key().first() == alpha1.name());
+            const bool matched = (tp.key().first() == alpha1.name());
 
             const scalar theta0 = degToRad(tp().theta0(matched));
             scalarField theta(boundary[patchi].size(), theta0);
 
-            scalar uTheta = tp().uTheta();
+            const scalar uTheta = tp().uTheta();
 
             // Calculate the dynamic contact angle if required
             if (uTheta > SMALL)
@@ -972,10 +971,10 @@ Foam::multiphaseMixtureThermo::nearInterface() const
         )
     );
 
-    forAllConstIter(PtrDictionary<phaseModel>, phases_, phase)
+    for (const phaseModel& phase : phases_)
     {
         tnearInt.ref() =
-            max(tnearInt(), pos0(phase() - 0.01)*pos0(0.99 - phase()));
+            max(tnearInt(), pos0(phase - 0.01)*pos0(0.99 - phase));
     }
 
     return tnearInt;
@@ -987,22 +986,20 @@ void Foam::multiphaseMixtureThermo::solveAlphas
     const scalar cAlpha
 )
 {
-    static label nSolves=-1;
-    nSolves++;
+    static label nSolves(-1);
+    ++nSolves;
 
-    word alphaScheme("div(phi,alpha)");
-    word alpharScheme("div(phirb,alpha)");
+    const word alphaScheme("div(phi,alpha)");
+    const word alpharScheme("div(phirb,alpha)");
 
     surfaceScalarField phic(mag(phi_/mesh_.magSf()));
     phic = min(cAlpha*phic, max(phic));
 
     PtrList<surfaceScalarField> alphaPhiCorrs(phases_.size());
-    int phasei = 0;
 
-    forAllIter(PtrDictionary<phaseModel>, phases_, phase)
+    int phasei = 0;
+    for (phaseModel& alpha : phases_)
     {
-        phaseModel& alpha = phase();
-
         alphaPhiCorrs.set
         (
             phasei,
@@ -1020,10 +1017,8 @@ void Foam::multiphaseMixtureThermo::solveAlphas
 
         surfaceScalarField& alphaPhiCorr = alphaPhiCorrs[phasei];
 
-        forAllIter(PtrDictionary<phaseModel>, phases_, phase2)
+        for (phaseModel& alpha2 : phases_)
         {
-            phaseModel& alpha2 = phase2();
-
             if (&alpha2 == &alpha) continue;
 
             surfaceScalarField phir(phic*nHatf(alpha, alpha2));
@@ -1050,7 +1045,7 @@ void Foam::multiphaseMixtureThermo::solveAlphas
             true
         );
 
-        phasei++;
+        ++phasei;
     }
 
     MULES::limitSum(alphaPhiCorrs);
@@ -1075,10 +1070,8 @@ void Foam::multiphaseMixtureThermo::solveAlphas
 
     phasei = 0;
 
-    forAllIter(PtrDictionary<phaseModel>, phases_, phase)
+    for (phaseModel& alpha : phases_)
     {
-        phaseModel& alpha = phase();
-
         surfaceScalarField& alphaPhi = alphaPhiCorrs[phasei];
         alphaPhi += upwind<scalar>(mesh_, phi_).flux(alpha);
 
@@ -1124,10 +1117,8 @@ void Foam::multiphaseMixtureThermo::solveAlphas
             }
         }
 
-        forAllConstIter(PtrDictionary<phaseModel>, phases_, phase2)
+        for (const phaseModel& alpha2 : phases_)
         {
-            const phaseModel& alpha2 = phase2();
-
             if (&alpha2 == &alpha) continue;
 
             const scalarField& dgdt2 = alpha2.dgdt();
@@ -1165,7 +1156,7 @@ void Foam::multiphaseMixtureThermo::solveAlphas
 
         sumAlpha += alpha;
 
-        phasei++;
+        ++phasei;
     }
 
     Info<< "Phase-sum volume fraction, min, max = "
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H
index e2ccdb9a4eb..479bc2cd34f 100644
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H
@@ -39,7 +39,7 @@
             ).ptr()
         );
 
-        phasei++;
+        ++phasei;
     }
 
     // Cache p_rgh prior to solve for density update
@@ -73,7 +73,7 @@
                 p_rghEqnComp.ref() += hmm;
             }
 
-            phasei++;
+            ++phasei;
         }
 
         solve
@@ -86,18 +86,13 @@
         if (pimple.finalNonOrthogonalIter())
         {
             phasei = 0;
-            forAllIter
-            (
-                PtrDictionary<phaseModel>,
-                mixture.phases(),
-                phase
-            )
+            for (phaseModel& phase : mixture.phases())
             {
-                phase().dgdt() =
-                    pos0(phase())
-                  *(p_rghEqnComps[phasei] & p_rgh)/phase().thermo().rho();
+                phase.dgdt() =
+                    pos0(phase)
+                  *(p_rghEqnComps[phasei] & p_rgh)/phase.thermo().rho();
 
-                phasei++;
+                ++phasei;
             }
 
             phi = phiHbyA + p_rghEqnIncomp.flux();
diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/YEqns.H b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/YEqns.H
index 643cd115723..31d10788bc5 100644
--- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/YEqns.H
+++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/YEqns.H
@@ -1,7 +1,6 @@
 {
-    forAllIter(UPtrList<phaseModel>, fluid.phases(), iter)
+    for (phaseModel& phase : fluid.phases())
     {
-        phaseModel& phase = iter();
         PtrList<volScalarField>& Y = phase.Y();
 
         if (!Y.empty())
@@ -55,4 +54,3 @@
         }
     }
 }
-
diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/laserDTRM/laserDTRM.C b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/laserDTRM/laserDTRM.C
index d57eeefd485..e57cc8df7a9 100644
--- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/laserDTRM/laserDTRM.C
+++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/laserDTRM/laserDTRM.C
@@ -133,7 +133,7 @@ void Foam::radiation::laserDTRM::initialiseReflection()
     {
         dictTable modelDicts(lookup("reflectionModel"));
 
-        forAllConstIter(dictTable, modelDicts, iter)
+        forAllConstIters(modelDicts, iter)
         {
             const phasePairKey& key = iter.key();
 
@@ -142,7 +142,7 @@ void Foam::radiation::laserDTRM::initialiseReflection()
                 key,
                 reflectionModel::New
                 (
-                    *iter,
+                    iter.object(),
                     mesh_
                 )
             );
@@ -612,20 +612,23 @@ void Foam::radiation::laserDTRM::calculate()
         reflectionUPtr.resize(reflections_.size());
 
         label reflectionModelId(0);
-        forAllIter(reflectionModelTable, reflections_, iter1)
+        forAllIters(reflections_, iter1)
         {
             reflectionModel& model = iter1()();
 
             reflectionUPtr.set(reflectionModelId, &model);
 
-            const word alpha1Name = "alpha." + iter1.key().first();
-            const word alpha2Name = "alpha." + iter1.key().second();
-
             const volScalarField& alphaFrom =
-                mesh_.lookupObject<volScalarField>(alpha1Name);
+                mesh_.lookupObject<volScalarField>
+                (
+                    IOobject::groupName("alpha", iter1.key().first())
+                );
 
             const volScalarField& alphaTo =
-                mesh_.lookupObject<volScalarField>(alpha2Name);
+                mesh_.lookupObject<volScalarField>
+                (
+                    IOobject::groupName("alpha", iter1.key().second())
+                );
 
             const volVectorField nHatPhase(nHatfv(alphaFrom, alphaTo));
 
@@ -700,9 +703,8 @@ void Foam::radiation::laserDTRM::calculate()
         DynamicList<point>  positionsMyProc;
         DynamicList<point>  p0MyProc;
 
-        forAllIter(Cloud<DTRMParticle>, DTRMCloud_, iter)
+        for (const DTRMParticle& p : DTRMCloud_)
         {
-            DTRMParticle& p = iter();
             positionsMyProc.append(p.position());
             p0MyProc.append(p.p0());
         }
diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/interfaceCompositionModel/interfaceCompositionModel.H b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/interfaceCompositionModel/interfaceCompositionModel.H
index 953a02ac191..ecf52f3f6bd 100644
--- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/interfaceCompositionModel/interfaceCompositionModel.H
+++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/massTransferModels/interfaceCompositionModel/interfaceCompositionModel.H
@@ -41,7 +41,6 @@ SourceFiles
 
 #include "volFields.H"
 #include "dictionary.H"
-#include "hashedWordList.H"
 #include "runTimeSelectionTables.H"
 #include "Enum.H"
 
diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/MassTransferPhaseSystem/MassTransferPhaseSystem.C b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/MassTransferPhaseSystem/MassTransferPhaseSystem.C
index eec86cb2b4d..7f304e34a5e 100644
--- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/MassTransferPhaseSystem/MassTransferPhaseSystem.C
+++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/MassTransferPhaseSystem/MassTransferPhaseSystem.C
@@ -44,7 +44,7 @@ Foam::MassTransferPhaseSystem<BasePhaseSystem>::MassTransferPhaseSystem
 {
     this->generatePairsAndSubModels("massTransferModel", massTransferModels_);
 
-    forAllConstIter(massTransferModelTable, massTransferModels_, iterModel)
+    forAllConstIters(massTransferModels_, iterModel)
     {
         if (!dmdt_.found(iterModel()->pair()))
         {
@@ -184,22 +184,17 @@ Foam::MassTransferPhaseSystem<BasePhaseSystem>::heatTransfer
 
     fvScalarMatrix& eqn = tEqnPtr.ref();
 
-    forAllIter(phaseSystem::phaseModelTable,  this->phaseModels_, iteri)
+    forAllConstIters(this->phaseModels_, iteri)
     {
-        phaseModel& phasei = iteri()();
+        const phaseModel& phasei = iteri()();
 
-        phaseSystem::phaseModelTable::iterator iterk = iteri;
-        iterk++;
-        for
-        (
-            ;
-            iterk != this->phaseModels_.end();
-            ++iterk
-        )
+        auto iterk = iteri;
+
+        for (++iterk; iterk != this->phaseModels_.end(); ++iterk)
         {
             if (iteri()().name() != iterk()().name())
             {
-                phaseModel& phasek = iterk()();
+                const phaseModel& phasek = iterk()();
 
                 // Phase i to phase k
                 const phasePairKey keyik(phasei.name(), phasek.name(), true);
@@ -290,7 +285,7 @@ void Foam::MassTransferPhaseSystem<BasePhaseSystem>::massSpeciesTransfer
 )
 {
     // Fill the volumetric mass transfer for species
-    forAllIter(massTransferModelTable, massTransferModels_, iter)
+    forAllConstIters(massTransferModels_, iter)
     {
         if (iter()->transferSpecie() == speciesName)
         {
diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C
index f07ecf7ba5c..aee73b89e95 100644
--- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C
+++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C
@@ -63,7 +63,7 @@ MultiComponentPhaseModel
         ).ptr()
     );
 
-    if (thermoPtr_->composition().species().size() == 0)
+    if (thermoPtr_->composition().species().empty())
     {
         FatalErrorInFunction
             << " The selected thermo is pure. Use a multicomponent thermo."
diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/multiphaseSystem/multiphaseSystem.C
index 2440b43e021..0fd59e03019 100644
--- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/multiphaseSystem/multiphaseSystem.C
+++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/multiphaseSystem/multiphaseSystem.C
@@ -65,18 +65,18 @@ Foam::multiphaseSystem::multiphaseSystem
     Su_(phaseModels_.size()),
     Sp_(phaseModels_.size())
 {
-    label phaseI = 0;
+    label phasei = 0;
     phases_.setSize(phaseModels_.size());
-    forAllConstIter(HashTable<autoPtr<phaseModel>>, phaseModels_, iter)
+    forAllIters(phaseModels_, iter)
     {
-        phaseModel& pm = const_cast<phaseModel&>(iter()());
-        phases_.set(phaseI++, &pm);
+        phaseModel& pm = iter()();
+        phases_.set(phasei++, &pm);
     }
 
     // Initiate Su and Sp
-    forAllConstIter(HashTable<autoPtr<phaseModel>>, phaseModels_, iter)
+    forAllConstIters(phaseModels_, iter)
     {
-        phaseModel& pm = const_cast<phaseModel&>(iter()());
+        const phaseModel& pm = iter()();
 
         Su_.insert
         (
@@ -117,7 +117,7 @@ Foam::multiphaseSystem::multiphaseSystem
 
 void Foam::multiphaseSystem::calculateSuSp()
 {
-    forAllIter(phaseSystem::phasePairTable, totalPhasePairs_, iter)
+    forAllConstIters(totalPhasePairs_, iter)
     {
         const phasePair& pair = iter()();
 
@@ -297,10 +297,9 @@ void Foam::multiphaseSystem::solve()
 
     for (int acorr=0; acorr<nAlphaCorr; acorr++)
     {
-        int phasei = 0;
-        forAllIter(UPtrList<phaseModel>, phases_, iter)
+        label phasei = 0;
+        for (phaseModel& phase1 : phases_)
         {
-            phaseModel& phase1 = iter();
             const volScalarField& alpha1 = phase1;
 
             phiAlphaCorrs.set
@@ -320,15 +319,11 @@ void Foam::multiphaseSystem::solve()
 
             surfaceScalarField& phiAlphaCorr = phiAlphaCorrs[phasei];
 
-            forAllIter(UPtrList<phaseModel>, phases_, iter2)
+            for (phaseModel& phase2 : phases_)
             {
-                phaseModel& phase2 = iter2();
                 const volScalarField& alpha2 = phase2;
 
-                if (&phase2 == &phase1)
-                {
-                    continue;
-                }
+                if (&phase2 == &phase1) continue;
 
                 const phasePairKey key12(phase1.name(), phase2.name());
 
@@ -337,7 +332,7 @@ void Foam::multiphaseSystem::solve()
                     FatalErrorInFunction
                         << "Phase compression factor (cAlpha) not found for : "
                         << key12
-                    << exit(FatalError);
+                        << exit(FatalError);
                 }
                 scalar cAlpha = cAlphas_.find(key12)();
 
@@ -380,13 +375,12 @@ void Foam::multiphaseSystem::solve()
                 }
             }
 
-            phasei++;
+            ++phasei;
         }
 
-        // Set Su and Sp tp zero
-        forAllIter(UPtrList<phaseModel>, phases_, iter)
+        // Set Su and Sp to zero
+        for (const phaseModel& phase : phases_)
         {
-            phaseModel& phase = iter();
             Su_[phase.name()] = dimensionedScalar("Su", dimless/dimTime, Zero);
             Sp_[phase.name()] = dimensionedScalar("Sp", dimless/dimTime, Zero);
 
@@ -402,9 +396,8 @@ void Foam::multiphaseSystem::solve()
 
         // Limit phiAlphaCorr on each phase
         phasei = 0;
-        forAllIter(UPtrList<phaseModel>, phases_, iter)
+        for (phaseModel& phase : phases_)
         {
-            phaseModel& phase = iter();
             volScalarField& alpha1 = phase;
 
             surfaceScalarField& phiAlphaCorr = phiAlphaCorrs[phasei];
@@ -425,7 +418,7 @@ void Foam::multiphaseSystem::solve()
                 0,
                 true
             );
-            phasei ++;
+            ++phasei;
         }
 
         MULES::limitSum(phiAlphaCorrs);
@@ -443,9 +436,8 @@ void Foam::multiphaseSystem::solve()
         );
 
         phasei = 0;
-        forAllIter(UPtrList<phaseModel>, phases_, iter)
+        for (phaseModel& phase : phases_)
         {
-            phaseModel& phase = iter();
             volScalarField& alpha1 = phase;
 
             const volScalarField::Internal& Su = Su_[phase.name()];
@@ -508,13 +500,9 @@ void Foam::multiphaseSystem::solve()
                 }
 
                 phase.alphaPhi() /= nAlphaSubCycles;
-
             }
             else
             {
-                phaseModel& phase = iter();
-                volScalarField& alpha1 = phase;
-
                 MULES::explicitSolve
                 (
                     geometricOneField(),
@@ -530,7 +518,7 @@ void Foam::multiphaseSystem::solve()
                 phase.alphaPhi() = phiAlpha;
             }
 
-            phasei++;
+            ++phasei;
         }
 
         if (acorr == nAlphaCorr - 1)
@@ -550,15 +538,13 @@ void Foam::multiphaseSystem::solve()
             // Reset rhoPhi
             rhoPhi_ = dimensionedScalar("rhoPhi", dimMass/dimTime, Zero);
 
-            forAllIter(UPtrList<phaseModel>, phases_, iter)
+            for (phaseModel& phase : phases_)
             {
-                phaseModel& phase = iter();
                 volScalarField& alpha1 = phase;
                 sumAlpha += alpha1;
 
                 // Update rhoPhi
-                rhoPhi_ +=
-                    fvc::interpolate(phase.rho())*phase.alphaPhi();
+                rhoPhi_ += fvc::interpolate(phase.rho()) * phase.alphaPhi();
 
             }
 
@@ -570,9 +556,8 @@ void Foam::multiphaseSystem::solve()
 
             volScalarField sumCorr(1.0 - sumAlpha);
 
-            forAllIter(UPtrList<phaseModel>, phases_, iter)
+            for (phaseModel& phase : phases_)
             {
-                phaseModel& phase = iter();
                 volScalarField& alpha = phase;
                 alpha += alpha*sumCorr;
 
@@ -619,24 +604,16 @@ Foam::dimensionedScalar Foam::multiphaseSystem::ddtAlphaMax() const
 
 Foam::scalar Foam::multiphaseSystem::maxDiffNo() const
 {
-    phaseModelTable::const_iterator phaseModelIter = phaseModels_.begin();
-    tmp<surfaceScalarField> kapparhoCpbyDelta(phaseModelIter()->diffNo());
+    auto iter = phaseModels_.cbegin();
 
-    scalar DiNum =
-        max(kapparhoCpbyDelta.ref()).value()*mesh_.time().deltaT().value();
+    scalar maxVal = max(iter()->diffNo()).value();
 
-    ++phaseModelIter;
-    for (; phaseModelIter != phaseModels_.end(); ++ phaseModelIter)
+    for (++iter; iter != phaseModels_.cend(); ++iter)
     {
-        kapparhoCpbyDelta = phaseModelIter()->diffNo();
-        DiNum =
-            max
-            (
-                DiNum,
-                max(kapparhoCpbyDelta).value()*mesh_.time().deltaT().value()
-            );
+        maxVal = max(maxVal, max(iter()->diffNo()).value());
     }
-    return DiNum;
+
+    return maxVal * mesh_.time().deltaT().value();
 }
 
 
diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/phaseSystem.C b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/phaseSystem.C
index 21b29248faa..ea7aa28fd19 100644
--- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/phaseSystem.C
+++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/phaseSystem.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2017-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -58,15 +58,15 @@ Foam::phaseSystem::generatePhaseModels(const wordList& phaseNames) const
 {
     phaseModelTable phaseModels;
 
-    forAllConstIter(wordList, phaseNames, phaseNameIter)
+    for (const word& phaseName : phaseNames)
     {
         phaseModels.insert
         (
-            *phaseNameIter,
+            phaseName,
             phaseModel::New
             (
                 *this,
-                *phaseNameIter
+                phaseName
             )
         );
     }
@@ -80,21 +80,17 @@ Foam::tmp<Foam::surfaceScalarField> Foam::phaseSystem::generatePhi
     const phaseModelTable& phaseModels
 ) const
 {
-    phaseModelTable::const_iterator phaseModelIter = phaseModels.begin();
+    auto iter = phaseModels.cbegin();
 
     auto tmpPhi = tmp<surfaceScalarField>::New
     (
         "phi",
-        fvc::interpolate(phaseModelIter()())*phaseModelIter()->phi()
+        fvc::interpolate(iter()()) * iter()->phi()
     );
 
-    ++phaseModelIter;
-
-    for (; phaseModelIter != phaseModels.end(); ++phaseModelIter)
+    for (++iter; iter != phaseModels.cend(); ++iter)
     {
-        tmpPhi.ref() +=
-            fvc::interpolate(phaseModelIter()())
-           *phaseModelIter()->phi();
+        tmpPhi.ref() += fvc::interpolate(iter()()) * iter()->phi();
     }
 
     return tmpPhi;
@@ -103,7 +99,7 @@ Foam::tmp<Foam::surfaceScalarField> Foam::phaseSystem::generatePhi
 
 void Foam::phaseSystem::generatePairs(const dictTable& modelDicts)
 {
-    forAllConstIter(dictTable, modelDicts, iter)
+    forAllConstIters(modelDicts, iter)
     {
         const phasePairKey& key = iter.key();
 
@@ -150,10 +146,9 @@ void Foam::phaseSystem::generatePairs(const dictTable& modelDicts)
 
 void Foam::phaseSystem::generatePairsTable()
 {
-
-    forAllConstIter(phaseModelTable, phaseModels_, phaseIter1)
+    forAllConstIters(phaseModels_, phaseIter1)
     {
-        forAllConstIter(phaseModelTable, phaseModels_, phaseIter2)
+        forAllConstIters(phaseModels_, phaseIter2)
         {
             if (phaseIter1()->name() != phaseIter2()->name())
             {
@@ -316,17 +311,16 @@ Foam::tmp<Foam::scalarField> Foam::phaseSystem::he
 
 Foam::tmp<Foam::volScalarField> Foam::phaseSystem::hc() const
 {
-    phaseModelTable::const_iterator phaseModelIter = phaseModels_.begin();
+    auto iter = phaseModels_.cbegin();
 
     tmp<volScalarField> tAlphaHc
     (
-        phaseModelIter()()*phaseModelIter()->hc()
+        iter()() * iter()->hc()
     );
 
-    ++phaseModelIter;
-    for (; phaseModelIter != phaseModels_.end(); phaseModelIter++)
+    for (++iter; iter != phaseModels_.cend(); ++iter)
     {
-        tAlphaHc.ref() += phaseModelIter()()*phaseModelIter()->hc();
+        tAlphaHc.ref() += iter()() * iter()->hc();
     }
 
     return tAlphaHc;
@@ -361,17 +355,16 @@ Foam::tmp<Foam::scalarField> Foam::phaseSystem::THE
 
 Foam::tmp<Foam::volScalarField> Foam::phaseSystem::rho() const
 {
-    phaseModelTable::const_iterator phaseModelIter = phaseModels_.begin();
+    auto iter = phaseModels_.cbegin();
 
     tmp<volScalarField> tmpRho
     (
-        phaseModelIter()()*phaseModelIter()->rho()
+        iter()() * iter()->rho()
     );
 
-    ++phaseModelIter;
-    for (; phaseModelIter != phaseModels_.end(); phaseModelIter++)
+    for (++iter; iter != phaseModels_.cend(); ++iter)
     {
-        tmpRho.ref() += phaseModelIter()()*phaseModelIter()->rho();
+        tmpRho.ref() += iter()() * iter()->rho();
     }
 
     return tmpRho;
@@ -380,20 +373,21 @@ Foam::tmp<Foam::volScalarField> Foam::phaseSystem::rho() const
 
 Foam::tmp<Foam::scalarField> Foam::phaseSystem::rho(const label patchI) const
 {
-    phaseModelTable::const_iterator phaseModelIter = phaseModels_.begin();
+    auto iter = phaseModels_.cbegin();
 
     tmp<scalarField> tmpRho
     (
-        phaseModelIter()().boundaryField()[patchI]
-      * phaseModelIter()->rho()().boundaryField()[patchI]
+        iter()().boundaryField()[patchI]
+      * iter()->rho()().boundaryField()[patchI]
     );
 
-    ++phaseModelIter;
-    for (; phaseModelIter != phaseModels_.end(); phaseModelIter++)
+    for (++iter; iter != phaseModels_.cend(); ++iter)
     {
         tmpRho.ref() +=
-            phaseModelIter()().boundaryField()[patchI]
-          * phaseModelIter()->rho()().boundaryField()[patchI];
+        (
+            iter()().boundaryField()[patchI]
+          * iter()->rho()().boundaryField()[patchI]
+        );
     }
 
     return tmpRho;
@@ -402,17 +396,16 @@ Foam::tmp<Foam::scalarField> Foam::phaseSystem::rho(const label patchI) const
 
 Foam::tmp<Foam::volScalarField> Foam::phaseSystem::Cp() const
 {
-    phaseModelTable::const_iterator phaseModelIter = phaseModels_.begin();
+    auto iter = phaseModels_.cbegin();
 
     tmp<volScalarField> tmpCp
     (
-        phaseModelIter()()*phaseModelIter()->Cp()
+        iter()() * iter()->Cp()
     );
 
-    ++phaseModelIter;
-    for (; phaseModelIter != phaseModels_.end(); phaseModelIter++)
+    for (++iter; iter != phaseModels_.cend(); ++iter)
     {
-        tmpCp.ref() += phaseModelIter()()*phaseModelIter()->Cp();
+        tmpCp.ref() += iter()() * iter()->Cp();
     }
 
     return tmpCp;
@@ -426,17 +419,16 @@ Foam::tmp<Foam::scalarField> Foam::phaseSystem::Cp
     const label patchI
 ) const
 {
-    phaseModelTable::const_iterator phaseModelIter = phaseModels_.begin();
+    auto iter = phaseModels_.cbegin();
 
     tmp<scalarField> tmpCp
     (
-        phaseModelIter()()*phaseModelIter()->Cp(p, T, patchI)
+        iter()() * iter()->Cp(p, T, patchI)
     );
 
-    ++phaseModelIter;
-    for (; phaseModelIter != phaseModels_.end(); ++phaseModelIter)
+    for (++iter; iter != phaseModels_.cend(); ++iter)
     {
-        tmpCp.ref() += phaseModelIter()()*phaseModelIter()->Cp(p, T, patchI);
+        tmpCp.ref() += iter()() * iter()->Cp(p, T, patchI);
     }
 
     return tmpCp;
@@ -445,17 +437,16 @@ Foam::tmp<Foam::scalarField> Foam::phaseSystem::Cp
 
 Foam::tmp<Foam::volScalarField> Foam::phaseSystem::Cv() const
 {
-    phaseModelTable::const_iterator phaseModelIter = phaseModels_.begin();
+    auto iter = phaseModels_.cbegin();
 
     tmp<volScalarField> tmpCv
     (
-        phaseModelIter()()*phaseModelIter()->Cv()
+        iter()() * iter()->Cv()
     );
 
-    ++phaseModelIter;
-    for (; phaseModelIter != phaseModels_.end(); phaseModelIter++)
+    for (++iter; iter != phaseModels_.cend(); ++iter)
     {
-        tmpCv.ref() += phaseModelIter()()*phaseModelIter()->Cv();
+        tmpCv.ref() += iter()() * iter()->Cv();
     }
 
     return tmpCv;
@@ -469,17 +460,16 @@ Foam::tmp<Foam::scalarField> Foam::phaseSystem::Cv
     const label patchI
 ) const
 {
-    phaseModelTable::const_iterator phaseModelIter = phaseModels_.begin();
+    auto iter = phaseModels_.cbegin();
 
     tmp<scalarField> tmpCv
     (
-        phaseModelIter()()*phaseModelIter()->Cv(p, T, patchI)
+        iter()() * iter()->Cv(p, T, patchI)
     );
 
-    ++phaseModelIter;
-    for (; phaseModelIter != phaseModels_.end(); ++ phaseModelIter)
+    for (++iter; iter != phaseModels_.cend(); ++iter)
     {
-        tmpCv.ref() += phaseModelIter()()*phaseModelIter()->Cv(p, T, patchI);
+        tmpCv.ref() += iter()() * iter()->Cv(p, T, patchI);
     }
 
     return tmpCv;
@@ -488,23 +478,22 @@ Foam::tmp<Foam::scalarField> Foam::phaseSystem::Cv
 
 Foam::tmp<Foam::volScalarField> Foam::phaseSystem::gamma() const
 {
-    phaseModelTable::const_iterator phaseModelIter = phaseModels_.begin();
+    auto iter = phaseModels_.cbegin();
 
     tmp<volScalarField> tmpCp
     (
-        phaseModelIter()()*phaseModelIter()->Cp()
+        iter()() * iter()->Cp()
     );
 
     tmp<volScalarField> tmpCv
     (
-        phaseModelIter()()*phaseModelIter()->Cv()
+        iter()() * iter()->Cv()
     );
 
-    ++phaseModelIter;
-    for (; phaseModelIter != phaseModels_.end(); phaseModelIter++)
+    for (++iter; iter != phaseModels_.cend(); ++iter)
     {
-        tmpCp.ref() += phaseModelIter()()*phaseModelIter()->Cp();
-        tmpCv.ref() += phaseModelIter()()*phaseModelIter()->Cv();
+        tmpCp.ref() += iter()() * iter()->Cp();
+        tmpCv.ref() += iter()() * iter()->Cv();
     }
 
     return (tmpCp/tmpCv);
@@ -527,17 +516,16 @@ Foam::tmp<Foam::scalarField> Foam::phaseSystem::gamma
 
 Foam::tmp<Foam::volScalarField> Foam::phaseSystem::Cpv() const
 {
-    phaseModelTable::const_iterator phaseModelIter = phaseModels_.begin();
+    auto iter = phaseModels_.cbegin();
 
     tmp<volScalarField> tmpCpv
     (
-        phaseModelIter()()*phaseModelIter()->Cpv()
+        iter()() * iter()->Cpv()
     );
 
-    ++phaseModelIter;
-    for (; phaseModelIter != phaseModels_.end(); phaseModelIter++)
+    for (++iter; iter != phaseModels_.cend(); ++iter)
     {
-        tmpCpv.ref() += phaseModelIter()()*phaseModelIter()->Cpv();
+        tmpCpv.ref() += iter()() * iter()->Cpv();
     }
 
     return tmpCpv;
@@ -551,17 +539,16 @@ Foam::tmp<Foam::scalarField> Foam::phaseSystem::Cpv
     const label patchI
 ) const
 {
-    phaseModelTable::const_iterator phaseModelIter = phaseModels_.begin();
+    auto iter = phaseModels_.cbegin();
 
     tmp<scalarField> tmpCpv
     (
-        phaseModelIter()()*phaseModelIter()->Cpv(p, T, patchI)
+        iter()() * iter()->Cpv(p, T, patchI)
     );
 
-    ++phaseModelIter;
-    for (; phaseModelIter != phaseModels_.end(); ++ phaseModelIter)
+    for (++iter; iter != phaseModels_.cend(); ++iter)
     {
-        tmpCpv.ref() += phaseModelIter()()*phaseModelIter()->Cpv(p, T, patchI);
+        tmpCpv.ref() += iter()() * iter()->Cpv(p, T, patchI);
     }
 
     return tmpCpv;
@@ -570,17 +557,16 @@ Foam::tmp<Foam::scalarField> Foam::phaseSystem::Cpv
 
 Foam::tmp<Foam::volScalarField> Foam::phaseSystem::CpByCpv() const
 {
-    phaseModelTable::const_iterator phaseModelIter = phaseModels_.begin();
+    auto iter = phaseModels_.cbegin();
 
     tmp<volScalarField> tmpCpByCpv
     (
-        phaseModelIter()()*phaseModelIter()->CpByCpv()
+        iter()() * iter()->CpByCpv()
     );
 
-    ++phaseModelIter;
-    for (; phaseModelIter != phaseModels_.end(); phaseModelIter++)
+    for (++iter; iter != phaseModels_.cend(); ++iter)
     {
-        tmpCpByCpv.ref() += phaseModelIter()()*phaseModelIter()->CpByCpv();
+        tmpCpByCpv.ref() += iter()() * iter()->CpByCpv();
     }
 
     return tmpCpByCpv;
@@ -594,20 +580,21 @@ Foam::tmp<Foam::scalarField> Foam::phaseSystem::CpByCpv
     const label patchI
 ) const
 {
-    phaseModelTable::const_iterator phaseModelIter = phaseModels_.begin();
+    auto iter = phaseModels_.cbegin();
 
     tmp<scalarField> tmpCpv
     (
-        phaseModelIter()().boundaryField()[patchI]
-       *phaseModelIter()->CpByCpv(p, T, patchI)
+        iter()().boundaryField()[patchI]
+      * iter()->CpByCpv(p, T, patchI)
     );
 
-    ++phaseModelIter;
-    for (; phaseModelIter != phaseModels_.end(); ++ phaseModelIter)
+    for (++iter; iter != phaseModels_.cend(); ++iter)
     {
         tmpCpv.ref() +=
-            phaseModelIter()().boundaryField()[patchI]
-           *phaseModelIter()->CpByCpv(p, T, patchI);
+        (
+            iter()().boundaryField()[patchI]
+          * iter()->CpByCpv(p, T, patchI)
+        );
     }
 
     return tmpCpv;
@@ -623,17 +610,16 @@ Foam::tmp<Foam::volScalarField> Foam::phaseSystem::W() const
 
 Foam::tmp<Foam::volScalarField> Foam::phaseSystem::kappa() const
 {
-    phaseModelTable::const_iterator phaseModelIter = phaseModels_.begin();
+    auto iter = phaseModels_.cbegin();
 
     tmp<volScalarField> tmpkappa
     (
-        phaseModelIter()()*phaseModelIter()->kappa()
+        iter()() * iter()->kappa()
     );
 
-    ++phaseModelIter;
-    for (; phaseModelIter != phaseModels_.end(); ++ phaseModelIter)
+    for (++iter; iter != phaseModels_.cend(); ++iter)
     {
-        tmpkappa.ref() += phaseModelIter()()*phaseModelIter()->kappa();
+        tmpkappa.ref() += iter()() * iter()->kappa();
     }
 
     return tmpkappa;
@@ -642,20 +628,21 @@ Foam::tmp<Foam::volScalarField> Foam::phaseSystem::kappa() const
 
 Foam::tmp<Foam::scalarField> Foam::phaseSystem::kappa(const label patchI) const
 {
-    phaseModelTable::const_iterator phaseModelIter = phaseModels_.begin();
+    auto iter = phaseModels_.cbegin();
 
     tmp<scalarField> tmpKappa
     (
-        phaseModelIter()().boundaryField()[patchI]
-       *phaseModelIter()->kappa(patchI)
+        iter()().boundaryField()[patchI]
+      * iter()->kappa(patchI)
     );
 
-    ++phaseModelIter;
-    for (; phaseModelIter != phaseModels_.end(); ++ phaseModelIter)
+    for (++iter; iter != phaseModels_.cend(); ++iter)
     {
         tmpKappa.ref() +=
-            phaseModelIter()().boundaryField()[patchI]
-           *phaseModelIter()->kappa(patchI);
+        (
+            iter()().boundaryField()[patchI]
+          * iter()->kappa(patchI)
+        );
     }
 
     return tmpKappa;
@@ -688,17 +675,16 @@ Foam::tmp<Foam::volScalarField> Foam::phaseSystem::alphaEff
     const volScalarField& alphat
 ) const
 {
-    phaseModelTable::const_iterator phaseModelIter = phaseModels_.begin();
+    auto iter = phaseModels_.cbegin();
 
     tmp<volScalarField> tmpAlpha
     (
-        phaseModelIter()()*phaseModelIter()->alpha()
+        iter()() * iter()->alpha()
     );
 
-    ++phaseModelIter;
-    for (; phaseModelIter != phaseModels_.end(); ++ phaseModelIter)
+    for (++iter; iter != phaseModels_.cend(); ++iter)
     {
-        tmpAlpha.ref() += phaseModelIter()()*phaseModelIter()->alpha();
+        tmpAlpha.ref() += iter()() * iter()->alpha();
     }
 
     tmpAlpha.ref() += alphat;
@@ -713,20 +699,21 @@ Foam::tmp<Foam::scalarField> Foam::phaseSystem::alphaEff
     const label patchI
 ) const
 {
-    phaseModelTable::const_iterator phaseModelIter = phaseModels_.begin();
+    auto iter = phaseModels_.cbegin();
 
     tmp<scalarField> tmpAlpha
     (
-        phaseModelIter()().boundaryField()[patchI]
-       *phaseModelIter()->alpha(patchI)
+        iter()().boundaryField()[patchI]
+      * iter()->alpha(patchI)
     );
 
-    ++phaseModelIter;
-    for (; phaseModelIter != phaseModels_.end(); ++ phaseModelIter)
+    for (++iter; iter != phaseModels_.cend(); ++iter)
     {
         tmpAlpha.ref() +=
-            phaseModelIter()().boundaryField()[patchI]
-           *phaseModelIter()->alpha(patchI);
+        (
+            iter()().boundaryField()[patchI]
+          * iter()->alpha(patchI)
+        );
     }
 
     tmpAlpha.ref() += alphat;
@@ -743,17 +730,16 @@ const Foam::dimensionedScalar& Foam::phaseSystem::Prt() const
 
 Foam::tmp<Foam::volScalarField> Foam::phaseSystem::mu() const
 {
-    phaseModelTable::const_iterator phaseModelIter = phaseModels_.begin();
+    auto iter = phaseModels_.cbegin();
 
     tmp<volScalarField> tmpMu
     (
-        phaseModelIter()()*phaseModelIter()->mu()
+        iter()() * iter()->mu()
     );
 
-    ++phaseModelIter;
-    for (; phaseModelIter != phaseModels_.end(); ++ phaseModelIter)
+    for (++iter; iter != phaseModels_.cend(); ++iter)
     {
-        tmpMu.ref() += phaseModelIter()()*phaseModelIter()->mu();
+        tmpMu.ref() += iter()() * iter()->mu();
     }
 
     return tmpMu;
@@ -762,20 +748,21 @@ Foam::tmp<Foam::volScalarField> Foam::phaseSystem::mu() const
 
 Foam::tmp<Foam::scalarField> Foam::phaseSystem::mu(const label patchI) const
 {
-    phaseModelTable::const_iterator phaseModelIter = phaseModels_.begin();
+    auto iter = phaseModels_.cbegin();
 
     tmp<scalarField> tmpMu
     (
-        phaseModelIter()().boundaryField()[patchI]
-       *phaseModelIter()->mu(patchI)
+        iter()().boundaryField()[patchI]
+      * iter()->mu(patchI)
     );
 
-    ++phaseModelIter;
-    for (; phaseModelIter != phaseModels_.end(); ++ phaseModelIter)
+    for (++iter; iter != phaseModels_.cend(); ++iter)
     {
         tmpMu.ref() +=
-            phaseModelIter()().boundaryField()[patchI]
-           *phaseModelIter()->mu(patchI);
+        (
+            iter()().boundaryField()[patchI]
+          * iter()->mu(patchI)
+        );
     }
 
     return tmpMu;
@@ -784,17 +771,16 @@ Foam::tmp<Foam::scalarField> Foam::phaseSystem::mu(const label patchI) const
 
 Foam::tmp<Foam::volScalarField> Foam::phaseSystem::nu() const
 {
-    phaseModelTable::const_iterator phaseModelIter = phaseModels_.begin();
+    auto iter = phaseModels_.cbegin();
 
     tmp<volScalarField> tmpNu
     (
-        phaseModelIter()()*phaseModelIter()->nu()
+        iter()() * iter()->nu()
     );
 
-    ++phaseModelIter;
-    for (; phaseModelIter != phaseModels_.end(); ++ phaseModelIter)
+    for (++iter; iter != phaseModels_.cend(); ++iter)
     {
-        tmpNu.ref() += phaseModelIter()()*phaseModelIter()->nu();
+        tmpNu.ref() += iter()() * iter()->nu();
     }
 
     return tmpNu;
@@ -803,25 +789,27 @@ Foam::tmp<Foam::volScalarField> Foam::phaseSystem::nu() const
 
 Foam::tmp<Foam::scalarField> Foam::phaseSystem::nu(const label patchI) const
 {
-    phaseModelTable::const_iterator phaseModelIter = phaseModels_.begin();
+    auto iter = phaseModels_.cbegin();
 
     tmp<scalarField> tmpNu
     (
-        phaseModelIter()().boundaryField()[patchI]
-       *phaseModelIter()->nu(patchI)
+        iter()().boundaryField()[patchI]
+      * iter()->nu(patchI)
     );
 
-    ++phaseModelIter;
-    for (; phaseModelIter != phaseModels_.end(); ++ phaseModelIter)
+    for (++iter; iter != phaseModels_.cend(); ++iter)
     {
         tmpNu.ref() +=
-            phaseModelIter()().boundaryField()[patchI]
-           *phaseModelIter()->nu(patchI);
+        (
+            iter()().boundaryField()[patchI]
+          * iter()->nu(patchI)
+        );
     }
 
     return tmpNu;
 }
 
+
 const Foam::surfaceScalarField& Foam::phaseSystem::phi() const
 {
     return phi_;
@@ -848,18 +836,18 @@ Foam::surfaceScalarField& Foam::phaseSystem::rhoPhi()
 
 void Foam::phaseSystem::correct()
 {
-    forAllIter(phaseModelTable, phaseModels_, phaseModelIter)
+    forAllIters(phaseModels_, iter)
     {
-        phaseModelIter()->correct();
+        iter()->correct();
     }
 }
 
 
 void Foam::phaseSystem::correctTurbulence()
 {
-    forAllIter(phaseModelTable, phaseModels_, phaseModelIter)
+    forAllIters(phaseModels_, iter)
     {
-        phaseModelIter()->correctTurbulence();
+        iter()->correctTurbulence();
     }
 }
 
@@ -891,14 +879,15 @@ Foam::phaseSystem::phasePairTable& Foam::phaseSystem::totalPhasePairs()
 
 bool Foam::phaseSystem::incompressible() const
 {
-    bool incompressible = true;
-
-    forAllConstIter(phaseModelTable, phaseModels_, phaseModelIter)
+    forAllConstIters(phaseModels_, iter)
     {
-        incompressible *= phaseModelIter()->thermo().incompressible();
+        if (!iter()->thermo().incompressible())
+        {
+            return false;
+        }
     }
 
-    return incompressible;
+    return true;
 }
 
 
@@ -910,14 +899,15 @@ bool Foam::phaseSystem::incompressible(const word phaseName) const
 
 bool Foam::phaseSystem::isochoric() const
 {
-     bool isochoric = true;
-
-    forAllConstIter(phaseModelTable, phaseModels_, phaseModelIter)
+    forAllConstIters(phaseModels_, iter)
     {
-        isochoric *= phaseModelIter()->thermo().isochoric();
+        if (!iter()->thermo().isochoric())
+        {
+            return false;
+        }
     }
 
-    return isochoric;
+    return true;
 }
 
 
@@ -939,22 +929,21 @@ Foam::phaseSystem::surfaceTensionForce() const
             mesh_
         ),
         mesh_,
-        dimensionedScalar(dimensionSet(1, -2, -2, 0, 0), Zero)
+        dimensionedScalar(dimForce, Zero)
     );
 
     auto& stf = tstf.ref();
     stf.setOriented();
 
-    if (surfaceTensionModels_.size() > 0)
+    if (surfaceTensionModels_.size())
     {
-        forAllConstIter(phaseModelTable, phaseModels_, iter1)
+        forAllConstIters(phaseModels_, iter1)
         {
             const volScalarField& alpha1 = iter1()();
 
-            phaseModelTable::const_iterator iter2 = iter1;
-            ++iter2;
+            auto iter2 = iter1;
 
-            for (; iter2 != phaseModels_.end(); ++iter2)
+            for (++iter2; iter2 != phaseModels_.cend(); ++iter2)
             {
                 const volScalarField& alpha2 = iter2()();
 
@@ -990,14 +979,14 @@ Foam::tmp<Foam::volVectorField> Foam::phaseSystem::U() const
             mesh_
         ),
         mesh_,
-        dimensionedVector("U", dimVelocity, Zero)
+        dimensionedVector(dimVelocity, Zero)
     );
 
     auto& stf = tstf.ref();
 
-    forAllConstIter(phaseModelTable, phaseModels_, iter1)
+    forAllConstIters(phaseModels_, iter)
     {
-        stf += iter1()()*iter1()->U();
+        stf += iter()() * iter()->U();
     }
 
     return tstf;
@@ -1025,22 +1014,17 @@ void Foam::phaseSystem::addInterfacePorosity(fvVectorMatrix& UEqn)
     const scalarField& Vc = mesh_.V();
     scalarField& Udiag = UEqn.diag();
 
-    forAllIter(phaseModelTable,  phaseModels_, iteri)
+    forAllConstIters(phaseModels_, iteri)
     {
         const phaseModel& phasei = iteri();
 
-        phaseModelTable::iterator iterk = iteri;
-        ++iterk;
-        for
-        (
-            ;
-            iterk != phaseModels_.end();
-            ++iterk
-        )
+        auto iterk = iteri;
+
+        for (++iterk; iterk != phaseModels_.cend(); ++iterk)
         {
             if (iteri()().name() != iterk()().name())
             {
-                phaseModel& phasek = iterk()();
+                const phaseModel& phasek = iterk()();
 
                 // Phase i and k
                 const phasePairKey keyik
@@ -1092,7 +1076,7 @@ Foam::tmp<Foam::volScalarField> Foam::phaseSystem::nearInterface
 
 Foam::tmp<Foam::volScalarField> Foam::phaseSystem::nearInterface() const
 {
-    auto tnI = tmp<volScalarField>::New
+    auto tnearInt = tmp<volScalarField>::New
     (
         IOobject
         (
@@ -1104,20 +1088,19 @@ Foam::tmp<Foam::volScalarField> Foam::phaseSystem::nearInterface() const
         dimensionedScalar(dimless, Zero)
     );
 
-    auto& nI = tnI.ref();
+    auto& nearInt = tnearInt.ref();
 
-    forAllConstIter(phaseModelTable, phaseModels_, iter1)
+    forAllConstIters(phaseModels_, iter1)
     {
         const volScalarField& alpha1 = iter1()();
 
-        phaseModelTable::const_iterator iter2 = iter1;
-        ++iter2;
+        auto iter2 = iter1;
 
-        for (; iter2 != phaseModels_.end(); ++iter2)
+        for (++iter2; iter2 != phaseModels_.cend(); ++iter2)
         {
             const volScalarField& alpha2 = iter2()();
 
-            nI +=
+            nearInt +=
             (
                 pos(alpha1 - 0.1)*pos(0.9 - alpha1)
                *pos(alpha2 - 0.1)*pos(0.9 - alpha2)
@@ -1125,7 +1108,7 @@ Foam::tmp<Foam::volScalarField> Foam::phaseSystem::nearInterface() const
         }
     }
 
-    return tnI;
+    return tnearInt;
 }
 
 
@@ -1168,9 +1151,7 @@ bool Foam::phaseSystem::read()
 {
     if (regIOobject::read())
     {
-        bool readOK = true;
-
-        return readOK;
+        return true;
     }
 
     return false;
diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/phaseSystemTemplates.H b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/phaseSystemTemplates.H
index b58184d98eb..642c88c3e09 100644
--- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/phaseSystemTemplates.H
+++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/phasesSystem/phaseSystem/phaseSystemTemplates.H
@@ -23,7 +23,6 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
 template<class modelType>
@@ -38,7 +37,7 @@ void Foam::phaseSystem::createSubModels
     >& models
 )
 {
-    forAllConstIter(dictTable, modelDicts, iter)
+    forAllConstIters(modelDicts, iter)
     {
         const phasePairKey& key = iter.key();
 
@@ -47,8 +46,8 @@ void Foam::phaseSystem::createSubModels
             key,
             modelType::New
             (
-               *iter,
-               phasePairs_[key]
+                iter.object(),
+                phasePairs_[key]
             )
         );
     }
@@ -68,7 +67,7 @@ void Foam::phaseSystem::createSubModels
     >& models
 )
 {
-    forAllConstIter(dictTable, modelDicts, iter)
+    forAllConstIters(modelDicts, iter)
     {
         const phasePairKey& key = iter.key();
 
@@ -77,8 +76,8 @@ void Foam::phaseSystem::createSubModels
             key,
             modelType::New
             (
-               *iter,
-               mesh
+                iter.object(),
+                mesh
             )
         );
     }
@@ -142,18 +141,18 @@ void Foam::phaseSystem::generatePairsAndSubModels
         HashTable<autoPtr<modelType>, phasePairKey, phasePairKey::hash>
         modelTypeTable;
 
-    forAllConstIter(wordList, phaseNames_, phaseNameIter)
+    for (const word& phaseName : phaseNames_)
     {
         modelTypeTable tempModels;
         generatePairsAndSubModels
         (
-            IOobject::groupName(modelName, *phaseNameIter),
+            IOobject::groupName(modelName, phaseName),
             tempModels
         );
 
-        forAllConstIter(typename modelTypeTable, tempModels, tempModelIter)
+        forAllConstIters(tempModels, tempModelIter)
         {
-            const phasePairKey key(tempModelIter.key());
+            const phasePairKey& key = tempModelIter.key();
 
             if (!models.found(key))
             {
@@ -164,9 +163,9 @@ void Foam::phaseSystem::generatePairsAndSubModels
                 );
             }
 
-            models[tempModelIter.key()].insert
+            models[key].insert
             (
-                *phaseNameIter,
+                phaseName,
                 *tempModelIter
             );
         }
diff --git a/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/twoPhaseMixtureEThermo/twoPhaseMixtureEThermo.C b/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/twoPhaseMixtureEThermo/twoPhaseMixtureEThermo.C
index 0782d26f727..6e9751ca11a 100644
--- a/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/twoPhaseMixtureEThermo/twoPhaseMixtureEThermo.C
+++ b/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/twoPhaseMixtureEThermo/twoPhaseMixtureEThermo.C
@@ -115,12 +115,6 @@ Foam::twoPhaseMixtureEThermo::twoPhaseMixtureEThermo
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
-
-Foam::twoPhaseMixtureEThermo::~twoPhaseMixtureEThermo()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 void Foam::twoPhaseMixtureEThermo::correct()
@@ -182,7 +176,7 @@ Foam::tmp<Foam::scalarField> Foam::twoPhaseMixtureEThermo::he
 
     forAll(T, i)
     {
-        label celli = cells[i];
+        const label celli = cells[i];
         he[i] =
             (
                 (T[i] - TSat_.value())
diff --git a/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/twoPhaseMixtureEThermo/twoPhaseMixtureEThermo.H b/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/twoPhaseMixtureEThermo/twoPhaseMixtureEThermo.H
index 5e5937d21fd..07c3c612de4 100644
--- a/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/twoPhaseMixtureEThermo/twoPhaseMixtureEThermo.H
+++ b/applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/twoPhaseMixtureEThermo/twoPhaseMixtureEThermo.H
@@ -94,7 +94,7 @@ public:
 
 
     //- Destructor
-    virtual ~twoPhaseMixtureEThermo();
+    virtual ~twoPhaseMixtureEThermo() = default;
 
 
     // Member Functions
@@ -162,14 +162,14 @@ public:
         //- i.e. rho != f(p)
         bool incompressible() const
         {
-            return (true);
+            return true;
         }
 
         //- Return true if the equation of state is isochoric
         //- i.e. rho = const
         bool isochoric() const
         {
-            return (false);
+            return false;
         }
 
         //- Return rho of the mixture
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/CourantNo.H b/applications/solvers/multiphase/multiphaseEulerFoam/CourantNo.H
index 005d01ad844..5b9935b0c23 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/CourantNo.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/CourantNo.H
@@ -39,12 +39,12 @@ if (mesh.nInternalFaces())
         fvc::surfaceSum(mag(phi))().primitiveField()
     );
 
-    forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+    for (const phaseModel& phase : fluid.phases())
     {
         sumPhi = max
         (
             sumPhi,
-            fvc::surfaceSum(mag(iter().phi()))().primitiveField()
+            fvc::surfaceSum(mag(phase.phi()))().primitiveField()
         );
     }
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/DDtU.H b/applications/solvers/multiphase/multiphaseEulerFoam/DDtU.H
index 04b66a307a6..afbfa692847 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/DDtU.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/DDtU.H
@@ -1,7 +1,5 @@
-forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+for (phaseModel& phase : fluid.phases())
 {
-    phaseModel& phase = iter();
-
     phase.DDtU() =
         fvc::ddt(phase.U())
       + fvc::div(phase.phi(), phase.U())
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/MRFCorrectBCs.H b/applications/solvers/multiphase/multiphaseEulerFoam/MRFCorrectBCs.H
index e871719cd8a..14a7f1b4a92 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/MRFCorrectBCs.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/MRFCorrectBCs.H
@@ -1,6 +1,6 @@
-    forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+    for (phaseModel& phase : fluid.phases())
     {
-        MRF.correctBoundaryVelocity(iter().U());
+        MRF.correctBoundaryVelocity(phase.U());
     }
 
     MRF.correctBoundaryVelocity(U);
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/UEqns.H b/applications/solvers/multiphase/multiphaseEulerFoam/UEqns.H
index 6c6da05f7a5..228e5ce50e0 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/UEqns.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/UEqns.H
@@ -3,14 +3,13 @@
 PtrList<fvVectorMatrix> UEqns(fluid.phases().size());
 autoPtr<multiphaseSystem::dragCoeffFields> dragCoeffs(fluid.dragCoeffs());
 
-int phasei = 0;
-forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+label phasei = 0;
+for (phaseModel& phase : fluid.phases())
 {
-    phaseModel& phase = iter();
     const volScalarField& alpha = phase;
     volVectorField& U = phase.U();
 
-    volScalarField nuEff(turbulence->nut() + iter().nu());
+    volScalarField nuEff(turbulence->nut() + phase.nu());
 
     UEqns.set
     (
@@ -58,5 +57,5 @@ forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
     );
     //UEqns[phasei].relax();
 
-    phasei++;
+    ++phasei;
 }
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/createFields.H b/applications/solvers/multiphase/multiphaseEulerFoam/createFields.H
index fa5ff769815..87a539ba752 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/createFields.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/createFields.H
@@ -42,9 +42,8 @@ surfaceScalarField phi
 
 multiphaseSystem fluid(U, phi);
 
-forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+for (const phaseModel& phase : fluid.phases())
 {
-    phaseModel& phase = iter();
     const volScalarField& alpha = phase;
 
     U += alpha*phase.U();
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/createMRFZones.H b/applications/solvers/multiphase/multiphaseEulerFoam/createMRFZones.H
index ea04f3efa04..00b57d132c2 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/createMRFZones.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/createMRFZones.H
@@ -1,8 +1,8 @@
     IOMRFZoneList MRF(mesh);
 
-    forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+    for (phaseModel& phase : fluid.phases())
     {
-        MRF.correctBoundaryVelocity(iter().U());
+        MRF.correctBoundaryVelocity(phase.U());
     }
 
     MRF.correctBoundaryVelocity(U);
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
index 30e47c042dc..209ce8f8d2c 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
@@ -44,9 +44,9 @@ void Foam::multiphaseSystem::calcAlphas()
     scalar level = 0.0;
     alphas_ == 0.0;
 
-    forAllIter(PtrDictionary<phaseModel>, phases_, iter)
+    for (const phaseModel& phase : phases_)
     {
-        alphas_ += level*iter();
+        alphas_ += level * phase;
         level += 1.0;
     }
 }
@@ -57,9 +57,8 @@ void Foam::multiphaseSystem::solveAlphas()
     PtrList<surfaceScalarField> alphaPhiCorrs(phases_.size());
     int phasei = 0;
 
-    forAllIter(PtrDictionary<phaseModel>, phases_, iter)
+    for (phaseModel& phase : phases_)
     {
-        phaseModel& phase = iter();
         volScalarField& alpha1 = phase;
 
         alphaPhiCorrs.set
@@ -79,21 +78,17 @@ void Foam::multiphaseSystem::solveAlphas()
 
         surfaceScalarField& alphaPhiCorr = alphaPhiCorrs[phasei];
 
-        forAllIter(PtrDictionary<phaseModel>, phases_, iter2)
+        for (phaseModel& phase2 : phases_)
         {
-            phaseModel& phase2 = iter2();
             volScalarField& alpha2 = phase2;
 
             if (&phase2 == &phase) continue;
 
             surfaceScalarField phir(phase.phi() - phase2.phi());
 
-            scalarCoeffSymmTable::const_iterator cAlpha
-            (
-                cAlphas_.find(interfacePair(phase, phase2))
-            );
+            const auto cAlpha = cAlphas_.cfind(interfacePair(phase, phase2));
 
-            if (cAlpha != cAlphas_.end())
+            if (cAlpha.found())
             {
                 surfaceScalarField phic
                 (
@@ -103,7 +98,7 @@ void Foam::multiphaseSystem::solveAlphas()
                 phir += min(cAlpha()*phic, max(phic))*nHatf(phase, phase2);
             }
 
-            word phirScheme
+            const word phirScheme
             (
                 "div(phir," + alpha2.name() + ',' + alpha1.name() + ')'
             );
@@ -132,7 +127,7 @@ void Foam::multiphaseSystem::solveAlphas()
             true
         );
 
-        phasei++;
+        ++phasei;
     }
 
     MULES::limitSum(alphaPhiCorrs);
@@ -151,10 +146,8 @@ void Foam::multiphaseSystem::solveAlphas()
 
     phasei = 0;
 
-    forAllIter(PtrDictionary<phaseModel>, phases_, iter)
+    for (phaseModel& phase : phases_)
     {
-        phaseModel& phase = iter();
-
         surfaceScalarField& alphaPhi = alphaPhiCorrs[phasei];
         alphaPhi += upwind<scalar>(mesh_, phi_).flux(phase);
         phase.correctInflowOutflow(alphaPhi);
@@ -178,7 +171,7 @@ void Foam::multiphaseSystem::solveAlphas()
 
         sumAlpha += phase;
 
-        phasei++;
+        ++phasei;
     }
 
     Info<< "Phase-sum volume fraction, min, max = "
@@ -189,9 +182,8 @@ void Foam::multiphaseSystem::solveAlphas()
 
     // Correct the sum of the phase-fractions to avoid 'drift'
     volScalarField sumCorr(1.0 - sumAlpha);
-    forAllIter(PtrDictionary<phaseModel>, phases_, iter)
+    for (phaseModel& phase : phases_)
     {
-        phaseModel& phase = iter();
         volScalarField& alpha = phase;
         alpha += alpha*sumCorr;
     }
@@ -270,11 +262,10 @@ void Foam::multiphaseSystem::correctContactAngle
                /mesh_.magSf().boundaryField()[patchi]
             );
 
-            alphaContactAngleFvPatchScalarField::thetaPropsTable::
-                const_iterator tp =
-                acap.thetaProps().find(interfacePair(phase1, phase2));
+            const auto tp =
+                acap.thetaProps().cfind(interfacePair(phase1, phase2));
 
-            if (tp == acap.thetaProps().end())
+            if (!tp.found())
             {
                 FatalErrorInFunction
                     << "Cannot find interface " << interfacePair(phase1, phase2)
@@ -416,7 +407,7 @@ Foam::multiphaseSystem::multiphaseSystem
 
     interfaceDictTable dragModelsDict(lookup("drag"));
 
-    forAllConstIter(interfaceDictTable, dragModelsDict, iter)
+    forAllConstIters(dragModelsDict, iter)
     {
         dragModels_.insert
         (
@@ -430,39 +421,24 @@ Foam::multiphaseSystem::multiphaseSystem
         );
     }
 
-    forAllConstIter(PtrDictionary<phaseModel>, phases_, iter1)
+    for (const phaseModel& phase1 : phases_)
     {
-        const phaseModel& phase1 = iter1();
-
-        forAllConstIter(PtrDictionary<phaseModel>, phases_, iter2)
+        for (const phaseModel& phase2 : phases_)
         {
-            const phaseModel& phase2 = iter2();
-
-            if (&phase2 != &phase1)
+            if (&phase2 == &phase1)
             {
-                scalarCoeffSymmTable::const_iterator sigma
-                (
-                    sigmas_.find(interfacePair(phase1, phase2))
-                );
+                continue;
+            }
 
-                if (sigma != sigmas_.end())
-                {
-                    scalarCoeffSymmTable::const_iterator cAlpha
-                    (
-                        cAlphas_.find(interfacePair(phase1, phase2))
-                    );
-
-                    if (cAlpha == cAlphas_.end())
-                    {
-                        WarningInFunction
-                          << "Compression coefficient not specified for "
-                             "phase pair ("
-                          << phase1.name() << ' ' << phase2.name()
-                          << ") for which a surface tension "
-                             "coefficient is specified"
-                          << endl;
-                    }
-                }
+            const interfacePair key(phase1, phase2);
+
+            if (sigmas_.found(key) && !cAlphas_.found(key))
+            {
+                WarningInFunction
+                    << "Compression coefficient not specified for phase pair ("
+                    << phase1.name() << ' ' << phase2.name()
+                    << ") for which a surface tension coefficient is specified"
+                    << endl;
             }
         }
     }
@@ -473,12 +449,12 @@ Foam::multiphaseSystem::multiphaseSystem
 
 Foam::tmp<Foam::volScalarField> Foam::multiphaseSystem::rho() const
 {
-    PtrDictionary<phaseModel>::const_iterator iter = phases_.begin();
+    auto iter = phases_.cbegin();
 
     tmp<volScalarField> trho = iter()*iter().rho();
     volScalarField& rho = trho.ref();
 
-    for (++iter; iter != phases_.end(); ++iter)
+    for (++iter; iter != phases_.cend(); ++iter)
     {
         rho += iter()*iter().rho();
     }
@@ -490,12 +466,12 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseSystem::rho() const
 Foam::tmp<Foam::scalarField>
 Foam::multiphaseSystem::rho(const label patchi) const
 {
-    PtrDictionary<phaseModel>::const_iterator iter = phases_.begin();
+    auto iter = phases_.cbegin();
 
     tmp<scalarField> trho = iter().boundaryField()[patchi]*iter().rho().value();
     scalarField& rho = trho.ref();
 
-    for (++iter; iter != phases_.end(); ++iter)
+    for (++iter; iter != phases_.cend(); ++iter)
     {
         rho += iter().boundaryField()[patchi]*iter().rho().value();
     }
@@ -506,12 +482,12 @@ Foam::multiphaseSystem::rho(const label patchi) const
 
 Foam::tmp<Foam::volScalarField> Foam::multiphaseSystem::nu() const
 {
-    PtrDictionary<phaseModel>::const_iterator iter = phases_.begin();
+    auto iter = phases_.cbegin();
 
     tmp<volScalarField> tmu = iter()*(iter().rho()*iter().nu());
     volScalarField& mu = tmu.ref();
 
-    for (++iter; iter != phases_.end(); ++iter)
+    for (++iter; iter != phases_.cend(); ++iter)
     {
         mu += iter()*(iter().rho()*iter().nu());
     }
@@ -523,14 +499,14 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseSystem::nu() const
 Foam::tmp<Foam::scalarField>
 Foam::multiphaseSystem::nu(const label patchi) const
 {
-    PtrDictionary<phaseModel>::const_iterator iter = phases_.begin();
+    auto iter = phases_.cbegin();
 
     tmp<scalarField> tmu =
         iter().boundaryField()[patchi]
        *(iter().rho().value()*iter().nu().value());
     scalarField& mu = tmu.ref();
 
-    for (++iter; iter != phases_.end(); ++iter)
+    for (++iter; iter != phases_.cend(); ++iter)
     {
         mu +=
             iter().boundaryField()[patchi]
@@ -557,33 +533,30 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseSystem::Cvm
                 mesh_
             ),
             mesh_,
-            dimensionedScalar(dimensionSet(1, -3, 0, 0, 0), Zero)
+            dimensionedScalar(dimDensity, Zero)
         )
     );
 
-    forAllConstIter(PtrDictionary<phaseModel>, phases_, iter)
+    for (const phaseModel& phase2 : phases_)
     {
-        const phaseModel& phase2 = iter();
+        if (&phase2 == &phase)
+        {
+            continue;
+        }
+
+        auto iterCvm = Cvms_.cfind(interfacePair(phase, phase2));
 
-        if (&phase2 != &phase)
+        if (iterCvm.found())
         {
-            scalarCoeffTable::const_iterator Cvm
-            (
-                Cvms_.find(interfacePair(phase, phase2))
-            );
+            tCvm.ref() += iterCvm()*phase2.rho()*phase2;
+        }
+        else
+        {
+            iterCvm = Cvms_.cfind(interfacePair(phase2, phase));
 
-            if (Cvm != Cvms_.end())
+            if (iterCvm.found())
             {
-                tCvm.ref() += Cvm()*phase2.rho()*phase2;
-            }
-            else
-            {
-                Cvm = Cvms_.find(interfacePair(phase2, phase));
-
-                if (Cvm != Cvms_.end())
-                {
-                    tCvm.ref() += Cvm()*phase.rho()*phase2;
-                }
+                tCvm.ref() += iterCvm()*phase.rho()*phase2;
             }
         }
     }
@@ -612,29 +585,26 @@ Foam::tmp<Foam::volVectorField> Foam::multiphaseSystem::Svm
         )
     );
 
-    forAllConstIter(PtrDictionary<phaseModel>, phases_, iter)
+    for (const phaseModel& phase2 : phases_)
     {
-        const phaseModel& phase2 = iter();
+        if (&phase2 == &phase)
+        {
+            continue;
+        }
+
+        auto Cvm = Cvms_.cfind(interfacePair(phase, phase2));
 
-        if (&phase2 != &phase)
+        if (Cvm.found())
         {
-            scalarCoeffTable::const_iterator Cvm
-            (
-                Cvms_.find(interfacePair(phase, phase2))
-            );
+            tSvm.ref() += Cvm()*phase2.rho()*phase2*phase2.DDtU();
+        }
+        else
+        {
+            Cvm = Cvms_.cfind(interfacePair(phase2, phase));
 
-            if (Cvm != Cvms_.end())
+            if (Cvm.found())
             {
-                tSvm.ref() += Cvm()*phase2.rho()*phase2*phase2.DDtU();
-            }
-            else
-            {
-                Cvm = Cvms_.find(interfacePair(phase2, phase));
-
-                if (Cvm != Cvms_.end())
-                {
-                    tSvm.ref() += Cvm()*phase.rho()*phase2*phase2.DDtU();
-                }
+                tSvm.ref() += Cvm()*phase.rho()*phase2*phase2.DDtU();
             }
         }
     }
@@ -666,7 +636,7 @@ Foam::multiphaseSystem::dragCoeffs() const
 {
     auto dragCoeffsPtr = autoPtr<dragCoeffFields>::New();
 
-    forAllConstIter(dragModelTable, dragModels_, iter)
+    forAllConstIters(dragModels_, iter)
     {
         const dragModel& dm = *iter();
 
@@ -739,7 +709,7 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseSystem::dragCoeff
     for
     (
         ;
-        dmIter != dragModels_.end() && dcIter != dragCoeffs.end();
+        dmIter.good() && dcIter.good();
         ++dmIter, ++dcIter
     )
     {
@@ -778,27 +748,24 @@ Foam::tmp<Foam::surfaceScalarField> Foam::multiphaseSystem::surfaceTension
     );
     tSurfaceTension.ref().setOriented();
 
-    forAllConstIter(PtrDictionary<phaseModel>, phases_, iter)
+    for (const phaseModel& phase2 : phases_)
     {
-        const phaseModel& phase2 = iter();
-
-        if (&phase2 != &phase1)
+        if (&phase2 == &phase1)
         {
-            scalarCoeffSymmTable::const_iterator sigma
-            (
-                sigmas_.find(interfacePair(phase1, phase2))
-            );
+            continue;
+        }
 
-            if (sigma != sigmas_.end())
-            {
-                tSurfaceTension.ref() +=
-                    dimensionedScalar("sigma", dimSigma_, sigma())
-                   *fvc::interpolate(K(phase1, phase2))*
-                    (
-                        fvc::interpolate(phase2)*fvc::snGrad(phase1)
-                      - fvc::interpolate(phase1)*fvc::snGrad(phase2)
-                    );
-            }
+        const auto sigma = sigmas_.cfind(interfacePair(phase1, phase2));
+
+        if (sigma.found())
+        {
+            tSurfaceTension.ref() +=
+                dimensionedScalar("sigma", dimSigma_, *sigma)
+               *fvc::interpolate(K(phase1, phase2))*
+                (
+                    fvc::interpolate(phase2)*fvc::snGrad(phase1)
+                  - fvc::interpolate(phase1)*fvc::snGrad(phase2)
+                );
         }
     }
 
@@ -824,10 +791,10 @@ Foam::multiphaseSystem::nearInterface() const
         )
     );
 
-    forAllConstIter(PtrDictionary<phaseModel>, phases_, iter)
+    for (const phaseModel& phase : phases_)
     {
         tnearInt.ref() =
-            max(tnearInt(), pos0(iter() - 0.01)*pos0(0.99 - iter()));
+            max(tnearInt(), pos0(phase - 0.01)*pos0(0.99 - phase));
     }
 
     return tnearInt;
@@ -836,9 +803,9 @@ Foam::multiphaseSystem::nearInterface() const
 
 void Foam::multiphaseSystem::solve()
 {
-    forAllIter(PtrDictionary<phaseModel>, phases_, iter)
+    for (phaseModel& phase : phases_)
     {
-        iter().correct();
+        phase.correct();
     }
 
     const Time& runTime = mesh_.time();
@@ -853,10 +820,9 @@ void Foam::multiphaseSystem::solve()
         PtrList<volScalarField> alpha0s(phases_.size());
         PtrList<surfaceScalarField> alphaPhiSums(phases_.size());
 
-        int phasei = 0;
-        forAllIter(PtrDictionary<phaseModel>, phases_, iter)
+        label phasei = 0;
+        for (phaseModel& phase : phases_)
         {
-            phaseModel& phase = iter();
             volScalarField& alpha = phase;
 
             alpha0s.set
@@ -881,7 +847,7 @@ void Foam::multiphaseSystem::solve()
                 )
             );
 
-            phasei++;
+            ++phasei;
         }
 
         for
@@ -896,18 +862,18 @@ void Foam::multiphaseSystem::solve()
         {
             solveAlphas();
 
-            int phasei = 0;
-            forAllIter(PtrDictionary<phaseModel>, phases_, iter)
+            label phasei = 0;
+            for (const phaseModel& phase : phases_)
             {
-                alphaPhiSums[phasei] += iter().alphaPhi()/nAlphaSubCycles;
-                phasei++;
+                alphaPhiSums[phasei] += phase.alphaPhi()/nAlphaSubCycles;
+
+                ++phasei;
             }
         }
 
         phasei = 0;
-        forAllIter(PtrDictionary<phaseModel>, phases_, iter)
+        for (phaseModel& phase : phases_)
         {
-            phaseModel& phase = iter();
             volScalarField& alpha = phase;
 
             phase.alphaPhi() = alphaPhiSums[phasei];
@@ -920,7 +886,7 @@ void Foam::multiphaseSystem::solve()
             alpha.oldTime() = alpha0s[phasei];
             alpha.oldTime().timeIndex() = runTime.timeIndex();
 
-            phasei++;
+            ++phasei;
         }
     }
     else
@@ -939,9 +905,9 @@ bool Foam::multiphaseSystem::read()
         PtrList<entry> phaseData(lookup("phases"));
         label phasei = 0;
 
-        forAllIter(PtrDictionary<phaseModel>, phases_, iter)
+        for (phaseModel& phase : phases_)
         {
-            readOK &= iter().read(phaseData[phasei++].dict());
+            readOK &= phase.read(phaseData[phasei++].dict());
         }
 
         readEntry("sigmas", sigmas_);
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/phaseModel/phaseModel.C b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/phaseModel/phaseModel.C
index b94cb287e6a..ec1a6b20335 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/phaseModel/phaseModel.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/phaseModel/phaseModel.C
@@ -229,12 +229,8 @@ bool Foam::phaseModel::read(const dictionary& phaseDict)
 
         return true;
     }
-    // else
-    // {
-    //     return false;
-    // }
 
-    return true;
+    // return false;
 }
 
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
index 6c2612b2319..4d3b927e4e3 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
@@ -31,17 +31,15 @@
     PtrList<surfaceScalarField> rAlphaAUfs(fluid.phases().size());
 
     phasei = 0;
-    forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+    for (phaseModel& phase : fluid.phases())
     {
-        phaseModel& phase = iter();
-
         MRF.makeAbsolute(phase.phi().oldTime());
         MRF.makeAbsolute(phase.phi());
 
         HbyAs.set(phasei, new volVectorField(phase.U()));
         phiHbyAs.set(phasei, new surfaceScalarField(1.0*phase.phi()));
 
-        phasei++;
+        ++phasei;
     }
 
     surfaceScalarField phiHbyA
@@ -62,9 +60,8 @@
     surfaceScalarField ghSnGradRho(ghf*fvc::snGrad(rho)*mesh.magSf());
 
     phasei = 0;
-    forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+    for (phaseModel& phase : fluid.phases())
     {
-        phaseModel& phase = iter();
         const volScalarField& alpha = phase;
 
         alphafs.set(phasei, fvc::interpolate(alpha).ptr());
@@ -115,14 +112,13 @@
              - ghSnGradRho
             )/phase.rho();
 
-        multiphaseSystem::dragModelTable::const_iterator dmIter =
-            fluid.dragModels().begin();
-        multiphaseSystem::dragCoeffFields::const_iterator dcIter =
-            dragCoeffs().begin();
+        auto dmIter = fluid.dragModels().cbegin();
+        auto dcIter = dragCoeffs().cbegin();
+
         for
         (
             ;
-            dmIter != fluid.dragModels().end() && dcIter != dragCoeffs().end();
+            dmIter.good() && dcIter.good();
             ++dmIter, ++dcIter
         )
         {
@@ -164,7 +160,7 @@
 
         phiHbyA += alphafs[phasei]*phiHbyAs[phasei];
 
-        phasei++;
+        ++phasei;
     }
 
     surfaceScalarField rAUf
@@ -180,12 +176,11 @@
     );
 
     phasei = 0;
-    forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+    for (const phaseModel& phase : fluid.phases())
     {
-        phaseModel& phase = iter();
         rAUf += mag(alphafs[phasei]*rAlphaAUfs[phasei])/phase.rho();
 
-        phasei++;
+        ++phasei;
     }
 
     // Update the fixedFluxPressure BCs to ensure flux consistency
@@ -193,15 +188,13 @@
         surfaceScalarField::Boundary phib(phi.boundaryField());
         phib = 0;
         phasei = 0;
-        forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+        for (const phaseModel& phase : fluid.phases())
         {
-            phaseModel& phase = iter();
-
             phib +=
                 alphafs[phasei].boundaryField()
                *(mesh.Sf().boundaryField() & phase.U().boundaryField());
 
-            phasei++;
+            ++phasei;
         }
 
         setSnGrad<fixedFluxPressureFvPatchScalarField>
@@ -238,20 +231,20 @@
             surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAUf);
 
             phasei = 0;
-            phi = dimensionedScalar("phi", phi.dimensions(), 0);
-            forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
-            {
-                phaseModel& phase = iter();
+            phi = dimensionedScalar("phi", phi.dimensions(), Zero);
 
+            for (phaseModel& phase : fluid.phases())
+            {
                 phase.phi() =
                     phiHbyAs[phasei]
                   + rAlphaAUfs[phasei]*mSfGradp/phase.rho();
+
                 phi +=
                     alphafs[phasei]*phiHbyAs[phasei]
                   + mag(alphafs[phasei]*rAlphaAUfs[phasei])
                    *mSfGradp/phase.rho();
 
-                phasei++;
+                ++phasei;
             }
 
             // dgdt =
@@ -270,9 +263,8 @@
             U = dimensionedVector("U", dimVelocity, Zero);
 
             phasei = 0;
-            forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+            for (phaseModel& phase : fluid.phases())
             {
-                phaseModel& phase = iter();
                 const volScalarField& alpha = phase;
 
                 phase.U() =
@@ -293,7 +285,7 @@
 
                 U += alpha*phase.U();
 
-                phasei++;
+                ++phasei;
             }
         }
     }
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/zonePhaseVolumes.H b/applications/solvers/multiphase/multiphaseEulerFoam/zonePhaseVolumes.H
index 5f8f3fb19a9..467351ab18d 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/zonePhaseVolumes.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/zonePhaseVolumes.H
@@ -5,14 +5,12 @@
     {
         const labelList& cellLabels = mesh.cellZones()[czi];
 
-        forAllConstIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+        for (const volScalarField& alpha : fluid.phases())
         {
-            const volScalarField& alpha = iter();
             scalar phaseVolume = 0;
 
-            forAll(cellLabels, cli)
+            for (const label celli : cellLabels)
             {
-                label celli = cellLabels[cli];
                 phaseVolume += alpha[celli]*V[celli];
             }
 
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
index fa49d98ad21..a8c4feb2ffb 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
+++ b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
@@ -42,9 +42,9 @@ void Foam::multiphaseMixture::calcAlphas()
     scalar level = 0.0;
     alphas_ == 0.0;
 
-    forAllIter(PtrDictionary<phase>, phases_, iter)
+    for (const phase& ph : phases_)
     {
-        alphas_ += level*iter();
+        alphas_ += level * ph;
         level += 1.0;
     }
 }
@@ -135,12 +135,12 @@ Foam::multiphaseMixture::multiphaseMixture
 Foam::tmp<Foam::volScalarField>
 Foam::multiphaseMixture::rho() const
 {
-    PtrDictionary<phase>::const_iterator iter = phases_.begin();
+    auto iter = phases_.cbegin();
 
     tmp<volScalarField> trho = iter()*iter().rho();
     volScalarField& rho = trho.ref();
 
-    for (++iter; iter != phases_.end(); ++iter)
+    for (++iter; iter != phases_.cend(); ++iter)
     {
         rho += iter()*iter().rho();
     }
@@ -152,12 +152,12 @@ Foam::multiphaseMixture::rho() const
 Foam::tmp<Foam::scalarField>
 Foam::multiphaseMixture::rho(const label patchi) const
 {
-    PtrDictionary<phase>::const_iterator iter = phases_.begin();
+    auto iter = phases_.cbegin();
 
     tmp<scalarField> trho = iter().boundaryField()[patchi]*iter().rho().value();
     scalarField& rho = trho.ref();
 
-    for (++iter; iter != phases_.end(); ++iter)
+    for (++iter; iter != phases_.cend(); ++iter)
     {
         rho += iter().boundaryField()[patchi]*iter().rho().value();
     }
@@ -169,12 +169,12 @@ Foam::multiphaseMixture::rho(const label patchi) const
 Foam::tmp<Foam::volScalarField>
 Foam::multiphaseMixture::mu() const
 {
-    PtrDictionary<phase>::const_iterator iter = phases_.begin();
+    auto iter = phases_.cbegin();
 
     tmp<volScalarField> tmu = iter()*iter().rho()*iter().nu();
     volScalarField& mu = tmu.ref();
 
-    for (++iter; iter != phases_.end(); ++iter)
+    for (++iter; iter != phases_.cend(); ++iter)
     {
         mu += iter()*iter().rho()*iter().nu();
     }
@@ -186,20 +186,25 @@ Foam::multiphaseMixture::mu() const
 Foam::tmp<Foam::scalarField>
 Foam::multiphaseMixture::mu(const label patchi) const
 {
-    PtrDictionary<phase>::const_iterator iter = phases_.begin();
+    auto iter = phases_.cbegin();
 
     tmp<scalarField> tmu =
+    (
         iter().boundaryField()[patchi]
        *iter().rho().value()
-       *iter().nu(patchi);
+       *iter().nu(patchi)
+    );
+
     scalarField& mu = tmu.ref();
 
-    for (++iter; iter != phases_.end(); ++iter)
+    for (++iter; iter != phases_.cend(); ++iter)
     {
         mu +=
+        (
             iter().boundaryField()[patchi]
            *iter().rho().value()
-           *iter().nu(patchi);
+           *iter().nu(patchi)
+        );
     }
 
     return tmu;
@@ -209,13 +214,13 @@ Foam::multiphaseMixture::mu(const label patchi) const
 Foam::tmp<Foam::surfaceScalarField>
 Foam::multiphaseMixture::muf() const
 {
-    PtrDictionary<phase>::const_iterator iter = phases_.begin();
+    auto iter = phases_.cbegin();
 
     tmp<surfaceScalarField> tmuf =
         fvc::interpolate(iter())*iter().rho()*fvc::interpolate(iter().nu());
     surfaceScalarField& muf = tmuf.ref();
 
-    for (++iter; iter != phases_.end(); ++iter)
+    for (++iter; iter != phases_.cend(); ++iter)
     {
         muf +=
             fvc::interpolate(iter())*iter().rho()*fvc::interpolate(iter().nu());
@@ -267,21 +272,19 @@ Foam::multiphaseMixture::surfaceTensionForce() const
     surfaceScalarField& stf = tstf.ref();
     stf.setOriented();
 
-    forAllConstIter(PtrDictionary<phase>, phases_, iter1)
+    forAllConstIters(phases_, iter1)
     {
         const phase& alpha1 = iter1();
 
-        PtrDictionary<phase>::const_iterator iter2 = iter1;
-        ++iter2;
+        auto iter2 = iter1;
 
-        for (; iter2 != phases_.end(); ++iter2)
+        for (++iter2; iter2 != phases_.cend(); ++iter2)
         {
             const phase& alpha2 = iter2();
 
-            sigmaTable::const_iterator sigma =
-                sigmas_.find(interfacePair(alpha1, alpha2));
+            auto sigma = sigmas_.cfind(interfacePair(alpha1, alpha2));
 
-            if (sigma == sigmas_.end())
+            if (!sigma.found())
             {
                 FatalErrorInFunction
                     << "Cannot find interface " << interfacePair(alpha1, alpha2)
@@ -289,7 +292,7 @@ Foam::multiphaseMixture::surfaceTensionForce() const
                     << exit(FatalError);
             }
 
-            stf += dimensionedScalar("sigma", dimSigma_, sigma())
+            stf += dimensionedScalar("sigma", dimSigma_, *sigma)
                *fvc::interpolate(K(alpha1, alpha2))*
                 (
                     fvc::interpolate(alpha2)*fvc::snGrad(alpha1)
@@ -354,9 +357,9 @@ void Foam::multiphaseMixture::solve()
 
 void Foam::multiphaseMixture::correct()
 {
-    forAllIter(PtrDictionary<phase>, phases_, iter)
+    for (phase& ph : phases_)
     {
-        iter().correct();
+        ph.correct();
     }
 }
 
@@ -431,11 +434,10 @@ void Foam::multiphaseMixture::correctContactAngle
                /mesh_.magSf().boundaryField()[patchi]
             );
 
-            alphaContactAngleFvPatchScalarField::thetaPropsTable::
-                const_iterator tp =
-                acap.thetaProps().find(interfacePair(alpha1, alpha2));
+            const auto tp =
+                acap.thetaProps().cfind(interfacePair(alpha1, alpha2));
 
-            if (tp == acap.thetaProps().end())
+            if (!tp.found())
             {
                 FatalErrorInFunction
                     << "Cannot find interface " << interfacePair(alpha1, alpha2)
@@ -444,12 +446,12 @@ void Foam::multiphaseMixture::correctContactAngle
                     << exit(FatalError);
             }
 
-            bool matched = (tp.key().first() == alpha1.name());
+            const bool matched = (tp.key().first() == alpha1.name());
 
             const scalar theta0 = degToRad(tp().theta0(matched));
             scalarField theta(boundary[patchi].size(), theta0);
 
-            scalar uTheta = tp().uTheta();
+            const scalar uTheta = tp().uTheta();
 
             // Calculate the dynamic contact angle if required
             if (uTheta > SMALL)
@@ -541,10 +543,9 @@ Foam::multiphaseMixture::nearInterface() const
         )
     );
 
-    forAllConstIter(PtrDictionary<phase>, phases_, iter)
+    for (const phase& ph : phases_)
     {
-        tnearInt.ref() =
-            max(tnearInt(), pos0(iter() - 0.01)*pos0(0.99 - iter()));
+        tnearInt.ref() = max(tnearInt(), pos0(ph - 0.01)*pos0(0.99 - ph));
     }
 
     return tnearInt;
@@ -556,11 +557,11 @@ void Foam::multiphaseMixture::solveAlphas
     const scalar cAlpha
 )
 {
-    static label nSolves=-1;
-    nSolves++;
+    static label nSolves(-1);
+    ++nSolves;
 
-    word alphaScheme("div(phi,alpha)");
-    word alpharScheme("div(phirb,alpha)");
+    const word alphaScheme("div(phi,alpha)");
+    const word alpharScheme("div(phirb,alpha)");
 
     surfaceScalarField phic(mag(phi_/mesh_.magSf()));
     phic = min(cAlpha*phic, max(phic));
@@ -568,10 +569,8 @@ void Foam::multiphaseMixture::solveAlphas
     PtrList<surfaceScalarField> alphaPhiCorrs(phases_.size());
     int phasei = 0;
 
-    forAllIter(PtrDictionary<phase>, phases_, iter)
+    for (phase& alpha : phases_)
     {
-        phase& alpha = iter();
-
         alphaPhiCorrs.set
         (
             phasei,
@@ -589,10 +588,8 @@ void Foam::multiphaseMixture::solveAlphas
 
         surfaceScalarField& alphaPhiCorr = alphaPhiCorrs[phasei];
 
-        forAllIter(PtrDictionary<phase>, phases_, iter2)
+        for (phase& alpha2 : phases_)
         {
-            phase& alpha2 = iter2();
-
             if (&alpha2 == &alpha) continue;
 
             surfaceScalarField phir(phic*nHatf(alpha, alpha2));
@@ -619,12 +616,12 @@ void Foam::multiphaseMixture::solveAlphas
             true
         );
 
-        phasei++;
+        ++phasei;
     }
 
     MULES::limitSum(alphaPhiCorrs);
 
-    rhoPhi_ = dimensionedScalar(dimensionSet(1, 0, -1, 0, 0), Zero);
+    rhoPhi_ = dimensionedScalar(dimMass/dimTime, Zero);
 
     volScalarField sumAlpha
     (
@@ -640,10 +637,8 @@ void Foam::multiphaseMixture::solveAlphas
 
     phasei = 0;
 
-    forAllIter(PtrDictionary<phase>, phases_, iter)
+    for (phase& alpha : phases_)
     {
-        phase& alpha = iter();
-
         surfaceScalarField& alphaPhi = alphaPhiCorrs[phasei];
         alphaPhi += upwind<scalar>(mesh_, phi_).flux(alpha);
 
@@ -666,7 +661,7 @@ void Foam::multiphaseMixture::solveAlphas
 
         sumAlpha += alpha;
 
-        phasei++;
+        ++phasei;
     }
 
     Info<< "Phase-sum volume fraction, min, max = "
@@ -677,9 +672,8 @@ void Foam::multiphaseMixture::solveAlphas
 
     // Correct the sum of the phase-fractions to avoid 'drift'
     volScalarField sumCorr(1.0 - sumAlpha);
-    forAllIter(PtrDictionary<phase>, phases_, iter)
+    for (phase& alpha : phases_)
     {
-        phase& alpha = iter();
         alpha += alpha*sumCorr;
     }
 
@@ -696,9 +690,9 @@ bool Foam::multiphaseMixture::read()
         PtrList<entry> phaseData(lookup("phases"));
         label phasei = 0;
 
-        forAllIter(PtrDictionary<phase>, phases_, iter)
+        for (phase& ph : phases_)
         {
-            readOK &= iter().read(phaseData[phasei++].dict());
+            readOK &= ph.read(phaseData[phasei++].dict());
         }
 
         readEntry("sigmas", sigmas_);
diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/interfaceCompositionModels/Henry/Henry.C b/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/interfaceCompositionModels/Henry/Henry.C
index 671bfaba4c6..00dc6556bfc 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/interfaceCompositionModels/Henry/Henry.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/interfaceCompositionModels/Henry/Henry.C
@@ -57,13 +57,6 @@ Foam::interfaceCompositionModels::Henry<Thermo, OtherThermo>::Henry
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class Thermo, class OtherThermo>
-Foam::interfaceCompositionModels::Henry<Thermo, OtherThermo>::~Henry()
-{}
-
-
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 template<class Thermo, class OtherThermo>
@@ -74,9 +67,9 @@ void Foam::interfaceCompositionModels::Henry<Thermo, OtherThermo>::update
 {
     YSolvent_ = scalar(1);
 
-    forAllConstIter(hashedWordList, this->speciesNames_, iter)
+    for (const word& speciesName : this->speciesNames_)
     {
-        YSolvent_ -= Yf(*iter, Tf);
+        YSolvent_ -= Yf(speciesName, Tf);
     }
 }
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/interfaceCompositionModels/Henry/Henry.H b/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/interfaceCompositionModels/Henry/Henry.H
index f6e49da88eb..1bae9a089f8 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/interfaceCompositionModels/Henry/Henry.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/interfaceCompositionModels/Henry/Henry.H
@@ -85,7 +85,7 @@ public:
 
 
     //- Destructor
-    virtual ~Henry();
+    virtual ~Henry() = default;
 
 
     // Member Functions
diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/interfaceCompositionModels/InterfaceCompositionModel/InterfaceCompositionModel.C b/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/interfaceCompositionModels/InterfaceCompositionModel/InterfaceCompositionModel.C
index f58cdbc562f..cbc83409d37 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/interfaceCompositionModels/InterfaceCompositionModel/InterfaceCompositionModel.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/interfaceCompositionModels/InterfaceCompositionModel/InterfaceCompositionModel.C
@@ -93,14 +93,6 @@ Foam::InterfaceCompositionModel<Thermo, OtherThermo>::InterfaceCompositionModel
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class Thermo, class OtherThermo>
-Foam::InterfaceCompositionModel<Thermo, OtherThermo>::
-~InterfaceCompositionModel()
-{}
-
-
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 template<class Thermo, class OtherThermo>
@@ -153,7 +145,7 @@ Foam::InterfaceCompositionModel<Thermo, OtherThermo>::D
         )
     );
 
-    volScalarField& D(tmpD.ref());
+    volScalarField& D = tmpD.ref();
 
     forAll(p, celli)
     {
@@ -207,7 +199,7 @@ Foam::InterfaceCompositionModel<Thermo, OtherThermo>::L
         )
     );
 
-    volScalarField& L(tmpL.ref());
+    volScalarField& L = tmpL.ref();
 
     forAll(p, celli)
     {
@@ -229,18 +221,18 @@ void Foam::InterfaceCompositionModel<Thermo, OtherThermo>::addMDotL
     volScalarField& mDotLPrime
 ) const
 {
-    forAllConstIter(hashedWordList, this->speciesNames_, iter)
+    for (const word& speciesName : this->speciesNames_)
     {
         volScalarField rhoKDL
         (
             thermo_.rhoThermo::rho()
            *K
-           *D(*iter)
-           *L(*iter, Tf)
+           *D(speciesName)
+           *L(speciesName, Tf)
         );
 
-        mDotL += rhoKDL*dY(*iter, Tf);
-        mDotLPrime += rhoKDL*YfPrime(*iter, Tf);
+        mDotL += rhoKDL*dY(speciesName, Tf);
+        mDotLPrime += rhoKDL*YfPrime(speciesName, Tf);
     }
 }
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/interfaceCompositionModels/InterfaceCompositionModel/InterfaceCompositionModel.H b/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/interfaceCompositionModels/InterfaceCompositionModel/InterfaceCompositionModel.H
index 6aeea6e097d..55c9e3ee821 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/interfaceCompositionModels/InterfaceCompositionModel/InterfaceCompositionModel.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/interfaceCompositionModels/InterfaceCompositionModel/InterfaceCompositionModel.H
@@ -59,7 +59,7 @@ class InterfaceCompositionModel
 {
 protected:
 
-    // Private data
+    // Protected Data
 
         //- Thermo
         const Thermo& thermo_;
@@ -71,7 +71,7 @@ protected:
         const dimensionedScalar Le_;
 
 
-    // Private member functions
+    // Protected Member Functions
 
         //- Get a reference to the local thermo for a pure mixture
         template<class ThermoType>
@@ -105,7 +105,7 @@ public:
 
 
     //- Destructor
-    ~InterfaceCompositionModel();
+    ~InterfaceCompositionModel() = default;
 
 
     // Member Functions
diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/interfaceCompositionModels/Raoult/Raoult.C b/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/interfaceCompositionModels/Raoult/Raoult.C
index 7a4245275b1..be9678378ea 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/interfaceCompositionModels/Raoult/Raoult.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/interfaceCompositionModels/Raoult/Raoult.C
@@ -58,16 +58,16 @@ Foam::interfaceCompositionModels::Raoult<Thermo, OtherThermo>::Raoult
         dimensionedScalar(dimless/dimTemperature, Zero)
     )
 {
-    forAllConstIter(hashedWordList, this->speciesNames_, iter)
+    for (const word& speciesName : this->speciesNames_)
     {
         speciesModels_.insert
         (
-            *iter,
+            speciesName,
             autoPtr<interfaceCompositionModel>
             (
                 interfaceCompositionModel::New
                 (
-                    dict.subDict(*iter),
+                    dict.subDict(speciesName),
                     pair
                 )
             )
@@ -76,13 +76,6 @@ Foam::interfaceCompositionModels::Raoult<Thermo, OtherThermo>::Raoult
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class Thermo, class OtherThermo>
-Foam::interfaceCompositionModels::Raoult<Thermo, OtherThermo>::~Raoult()
-{}
-
-
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 template<class Thermo, class OtherThermo>
@@ -93,12 +86,7 @@ void Foam::interfaceCompositionModels::Raoult<Thermo, OtherThermo>::update
 {
     YNonVapour_ = scalar(1);
 
-    forAllIter
-    (
-        HashTable<autoPtr<interfaceCompositionModel>>,
-        speciesModels_,
-        iter
-    )
+    forAllIters(speciesModels_, iter)
     {
         iter()->update(Tf);
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/interfaceCompositionModels/Raoult/Raoult.H b/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/interfaceCompositionModels/Raoult/Raoult.H
index c295f8f6fb6..a55116c5959 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/interfaceCompositionModels/Raoult/Raoult.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/interfaceCompositionModels/Raoult/Raoult.H
@@ -86,7 +86,7 @@ public:
 
 
     //- Destructor
-    virtual ~Raoult();
+    virtual ~Raoult() = default;
 
 
     // Member Functions
diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/interfaceCompositionModels/interfaceCompositionModel/interfaceCompositionModel.C b/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/interfaceCompositionModels/interfaceCompositionModel/interfaceCompositionModel.C
index 583854f3a5f..e4bcf250986 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/interfaceCompositionModels/interfaceCompositionModel/interfaceCompositionModel.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/interfaceCompositionModels/interfaceCompositionModel/interfaceCompositionModel.C
@@ -49,12 +49,6 @@ Foam::interfaceCompositionModel::interfaceCompositionModel
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::interfaceCompositionModel::~interfaceCompositionModel()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 const Foam::hashedWordList& Foam::interfaceCompositionModel::species() const
diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/interfaceCompositionModels/interfaceCompositionModel/interfaceCompositionModel.H b/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/interfaceCompositionModels/interfaceCompositionModel/interfaceCompositionModel.H
index bfbdee52b69..26067ab4cab 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/interfaceCompositionModels/interfaceCompositionModel/interfaceCompositionModel.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/interfaceCompositionModels/interfaceCompositionModel/interfaceCompositionModel.H
@@ -99,7 +99,7 @@ public:
 
 
     //- Destructor
-    virtual ~interfaceCompositionModel();
+    virtual ~interfaceCompositionModel() = default;
 
 
     // Selectors
@@ -120,7 +120,7 @@ public:
         const hashedWordList& species() const;
 
         //- Returns whether the species is transported by the model and
-        //  provides the name of the diffused species
+        //- provides the name of the diffused species
         bool transports
         (
             word& speciesName
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/BlendedInterfacialModel/blendingMethods/hyperbolic/hyperbolic.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/BlendedInterfacialModel/blendingMethods/hyperbolic/hyperbolic.C
index 930e483d9ac..fd2823261ac 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/BlendedInterfacialModel/blendingMethods/hyperbolic/hyperbolic.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/BlendedInterfacialModel/blendingMethods/hyperbolic/hyperbolic.C
@@ -54,25 +54,22 @@ Foam::blendingMethods::hyperbolic::hyperbolic
     blendingMethod(dict),
     transitionAlphaScale_("transitionAlphaScale", dimless, dict)
 {
-    forAllConstIter(wordList, phaseNames, iter)
+    for (const word& phaseName : phaseNames)
     {
-        const word name(IOobject::groupName("minContinuousAlpha", *iter));
-
         minContinuousAlpha_.insert
         (
-            *iter,
-            dimensionedScalar(name, dimless, dict)
+            phaseName,
+            dimensionedScalar
+            (
+                IOobject::groupName("minContinuousAlpha", phaseName),
+                dimless,
+                dict
+            )
         );
     }
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::blendingMethods::hyperbolic::~hyperbolic()
-{}
-
-
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::volScalarField> Foam::blendingMethods::hyperbolic::f1
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/BlendedInterfacialModel/blendingMethods/hyperbolic/hyperbolic.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/BlendedInterfacialModel/blendingMethods/hyperbolic/hyperbolic.H
index 7f62628ca36..b8695016cb7 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/BlendedInterfacialModel/blendingMethods/hyperbolic/hyperbolic.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/BlendedInterfacialModel/blendingMethods/hyperbolic/hyperbolic.H
@@ -77,7 +77,7 @@ public:
 
 
     //- Destructor
-    ~hyperbolic();
+    ~hyperbolic() = default;
 
 
     // Member Functions
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/BlendedInterfacialModel/blendingMethods/linear/linear.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/BlendedInterfacialModel/blendingMethods/linear/linear.C
index 5fc159420df..4bdf362bb92 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/BlendedInterfacialModel/blendingMethods/linear/linear.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/BlendedInterfacialModel/blendingMethods/linear/linear.C
@@ -53,39 +53,39 @@ Foam::blendingMethods::linear::linear
 :
     blendingMethod(dict)
 {
-    forAllConstIter(wordList, phaseNames, iter)
+    for (const word& phaseName : phaseNames)
     {
-        const word nameFull
-        (
-            IOobject::groupName("minFullyContinuousAlpha", *iter)
-        );
-
         minFullyContinuousAlpha_.insert
         (
-            *iter,
-            dimensionedScalar(nameFull, dimless, dict)
-        );
-
-        const word namePart
-        (
-            IOobject::groupName("minPartlyContinuousAlpha", *iter)
+            phaseName,
+            dimensionedScalar
+            (
+                IOobject::groupName("minFullyContinuousAlpha", phaseName),
+                dimless,
+                dict
+            )
         );
 
         minPartlyContinuousAlpha_.insert
         (
-            *iter,
-            dimensionedScalar(namePart, dimless, dict)
+            phaseName,
+            dimensionedScalar
+            (
+                IOobject::groupName("minPartlyContinuousAlpha", phaseName),
+                dimless,
+                dict
+            )
         );
 
         if
         (
-            minFullyContinuousAlpha_[*iter]
-          < minPartlyContinuousAlpha_[*iter]
+            minFullyContinuousAlpha_[phaseName]
+          < minPartlyContinuousAlpha_[phaseName]
         )
         {
             FatalErrorInFunction
                 << "The supplied fully continuous volume fraction for "
-                << *iter
+                << phaseName
                 << " is less than the partly continuous value."
                 << endl << exit(FatalError);
         }
@@ -93,12 +93,6 @@ Foam::blendingMethods::linear::linear
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::blendingMethods::linear::~linear()
-{}
-
-
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::volScalarField> Foam::blendingMethods::linear::f1
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/BlendedInterfacialModel/blendingMethods/linear/linear.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/BlendedInterfacialModel/blendingMethods/linear/linear.H
index 13072bf5453..eb1affbbafb 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/BlendedInterfacialModel/blendingMethods/linear/linear.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/BlendedInterfacialModel/blendingMethods/linear/linear.H
@@ -77,7 +77,7 @@ public:
 
 
     //- Destructor
-    ~linear();
+    ~linear() = default;
 
 
     // Member Functions
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C
index b985ef2ad00..95cb33d1a2f 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C
@@ -359,16 +359,10 @@ bool Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::read()
 {
     if (BasePhaseSystem::read())
     {
-        bool readOK = true;
-
-        // Models ...
-
-        return readOK;
-    }
-    else
-    {
-        return false;
+        return true;
     }
+
+    return false;
 }
 
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.C
index 7f8b768c941..97a173bd4c4 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.C
@@ -164,16 +164,10 @@ bool Foam::HeatTransferPhaseSystem<BasePhaseSystem>::read()
 {
     if (BasePhaseSystem::read())
     {
-        bool readOK = true;
-
-        // Models ...
-
-        return readOK;
-    }
-    else
-    {
-        return false;
+        return true;
     }
+
+    return false;
 }
 
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.C
index 94c67fa820b..2a213d61f52 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.C
@@ -97,17 +97,13 @@ massTransfer() const
     }
 
     // Sum up the contribution from each interface composition model
-    forAllConstIters
-    (
-        interfaceCompositionModels_,
-        interfaceCompositionModelIter
-    )
+    forAllConstIters(interfaceCompositionModels_, modelIter)
     {
         const phasePair& pair =
-            *(this->phasePairs_[interfaceCompositionModelIter.key()]);
+            *(this->phasePairs_[modelIter.key()]);
 
         const interfaceCompositionModel& compositionModel =
-            *(interfaceCompositionModelIter.object());
+            *(modelIter.object());
 
         const phaseModel& phase = pair.phase1();
         const phaseModel& otherPhase = pair.phase2();
@@ -288,16 +284,10 @@ bool Foam::InterfaceCompositionPhaseChangePhaseSystem<BasePhaseSystem>::read()
 {
     if (BasePhaseSystem::read())
     {
-        bool readOK = true;
-
-        // Models ...
-
-        return readOK;
-    }
-    else
-    {
-        return false;
+        return true;
     }
+
+    return false;
 }
 
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C
index 619dbd16316..3e5a515619a 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C
@@ -592,16 +592,10 @@ bool Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::read()
 {
     if (BasePhaseSystem::read())
     {
-        bool readOK = true;
-
-        // Read models ...
-
-        return readOK;
-    }
-    else
-    {
-        return false;
+        return true;
     }
+
+    return false;
 }
 
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
index dfd3a346cc1..88c8725a0a3 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
@@ -211,7 +211,7 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::massTransfer() const
         const phaseModel& phase = pair.phase1();
         const phaseModel& otherPhase = pair.phase2();
 
-        const word name
+        const word thisName
         (
             IOobject::groupName(volatile_, phase.name())
         );
@@ -225,7 +225,7 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::massTransfer() const
         const volScalarField dmdt12(posPart(dmdt));
         const volScalarField dmdt21(negPart(dmdt));
 
-        *eqns[name] += fvm::Sp(dmdt21, eqns[name]->psi()) - dmdt21;
+        *eqns[thisName] += fvm::Sp(dmdt21, eqns[thisName]->psi()) - dmdt21;
         *eqns[otherName] += dmdt12 - fvm::Sp(dmdt12, eqns[otherName]->psi());
     }
 
@@ -485,16 +485,10 @@ bool Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::read()
 {
     if (BasePhaseSystem::read())
     {
-        bool readOK = true;
-
-        // Models ...
-
-        return readOK;
-    }
-    else
-    {
-        return false;
+        return true;
     }
+
+    return false;
 }
 
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/IsothermalPhaseModel/IsothermalPhaseModel.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/IsothermalPhaseModel/IsothermalPhaseModel.C
index 8677e02a475..cf84c8435ae 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/IsothermalPhaseModel/IsothermalPhaseModel.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/IsothermalPhaseModel/IsothermalPhaseModel.C
@@ -40,13 +40,6 @@ Foam::IsothermalPhaseModel<BasePhaseModel>::IsothermalPhaseModel
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class BasePhaseModel>
-Foam::IsothermalPhaseModel<BasePhaseModel>::~IsothermalPhaseModel()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class BasePhaseModel>
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/IsothermalPhaseModel/IsothermalPhaseModel.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/IsothermalPhaseModel/IsothermalPhaseModel.H
index f8bb151fe30..086b20d840c 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/IsothermalPhaseModel/IsothermalPhaseModel.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/IsothermalPhaseModel/IsothermalPhaseModel.H
@@ -66,7 +66,7 @@ public:
 
 
     //- Destructor
-    virtual ~IsothermalPhaseModel();
+    virtual ~IsothermalPhaseModel() = default;
 
 
     // Member Functions
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.C
index 7848257f810..f78aeef07fd 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.C
@@ -55,7 +55,7 @@ Foam::tmp<Foam::surfaceScalarField> Foam::phaseSystem::calcPhi
         )
     );
 
-    for (label phasei=1; phasei<phaseModels.size(); phasei++)
+    for (label phasei=1; phasei<phaseModels.size(); ++phasei)
     {
         tmpPhi.ref() +=
             fvc::interpolate(phaseModels[phasei])*phaseModels[phasei].phi();
@@ -70,7 +70,7 @@ void Foam::phaseSystem::generatePairs
     const dictTable& modelDicts
 )
 {
-    forAllConstIter(dictTable, modelDicts, iter)
+    forAllConstIters(modelDicts, iter)
     {
         const phasePairKey& key = iter.key();
 
@@ -188,14 +188,16 @@ Foam::phaseSystem::~phaseSystem()
 
 Foam::tmp<Foam::volScalarField> Foam::phaseSystem::rho() const
 {
+    auto phasei = phaseModels_.cbegin();
+
     tmp<volScalarField> tmpRho
     (
-        phaseModels_[0]*phaseModels_[0].rho()
+        phasei() * phasei().rho()
     );
 
-    for (label phasei=1; phasei<phaseModels_.size(); phasei++)
+    for (++phasei; phasei != phaseModels_.cend(); ++phasei)
     {
-        tmpRho.ref() += phaseModels_[phasei]*phaseModels_[phasei].rho();
+        tmpRho.ref() += phasei() * phasei().rho();
     }
 
     return tmpRho;
@@ -204,14 +206,16 @@ Foam::tmp<Foam::volScalarField> Foam::phaseSystem::rho() const
 
 Foam::tmp<Foam::volVectorField> Foam::phaseSystem::U() const
 {
+    auto phasei = phaseModels_.cbegin();
+
     tmp<volVectorField> tmpU
     (
-        phaseModels_[0]*phaseModels_[0].U()
+        phasei() * phasei().U()
     );
 
-    for (label phasei=1; phasei<phaseModels_.size(); phasei++)
+    for (++phasei; phasei != phaseModels_.cend(); ++phasei)
     {
-        tmpU.ref() += phaseModels_[phasei]*phaseModels_[phasei].U();
+        tmpU.ref() += phasei() * phasei().U();
     }
 
     return tmpU;
@@ -274,9 +278,9 @@ void Foam::phaseSystem::solve()
 
 void Foam::phaseSystem::correct()
 {
-    forAll(phaseModels_, phasei)
+    for (phaseModel& phase : phaseModels_)
     {
-        phaseModels_[phasei].correct();
+        phase.correct();
     }
 }
 
@@ -285,44 +289,44 @@ void Foam::phaseSystem::correctKinematics()
 {
     bool updateDpdt = false;
 
-    forAll(phaseModels_, phasei)
+    for (phaseModel& phase : phaseModels_)
     {
-        phaseModels_[phasei].correctKinematics();
+        phase.correctKinematics();
 
-        updateDpdt = updateDpdt || phaseModels_[phasei].thermo().dpdt();
+        updateDpdt = updateDpdt || phase.thermo().dpdt();
     }
 
     // Update the pressure time-derivative if required
     if (updateDpdt)
     {
-        dpdt_ = fvc::ddt(phaseModels_.begin()().thermo().p());
+        dpdt_ = fvc::ddt(phaseModels_.cbegin()().thermo().p());
     }
 }
 
 
 void Foam::phaseSystem::correctThermo()
 {
-    forAll(phaseModels_, phasei)
+    for (phaseModel& phase : phaseModels_)
     {
-        phaseModels_[phasei].correctThermo();
+        phase.correctThermo();
     }
 }
 
 
 void Foam::phaseSystem::correctTurbulence()
 {
-    forAll(phaseModels_, phasei)
+    for (phaseModel& phase : phaseModels_)
     {
-        phaseModels_[phasei].correctTurbulence();
+        phase.correctTurbulence();
     }
 }
 
 
 void Foam::phaseSystem::correctEnergyTransport()
 {
-    forAll(phaseModels_, phasei)
+    for (phaseModel& phase : phaseModels_)
     {
-        phaseModels_[phasei].correctEnergyTransport();
+        phase.correctEnergyTransport();
     }
 }
 
@@ -333,19 +337,17 @@ bool Foam::phaseSystem::read()
     {
         bool readOK = true;
 
-        forAll(phaseModels_, phasei)
+        for (phaseModel& phase : phaseModels_)
         {
-            readOK &= phaseModels_[phasei].read();
+            readOK &= phase.read();
         }
 
         // models ...
 
         return readOK;
     }
-    else
-    {
-        return false;
-    }
+
+    return false;
 }
 
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystemTemplates.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystemTemplates.C
index 768c905e1ff..824f98b33cb 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystemTemplates.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystemTemplates.C
@@ -39,7 +39,7 @@ void Foam::phaseSystem::createSubModels
     >& models
 )
 {
-    forAllConstIter(dictTable, modelDicts, iter)
+    forAllConstIters(modelDicts, iter)
     {
         const phasePairKey& key = iter.key();
 
@@ -104,7 +104,7 @@ void Foam::phaseSystem::generatePairsAndSubModels
 
     autoPtr<modelType> noModel;
 
-    forAllConstIter(typename modelTypeTable, tempModels, iter)
+    forAllConstIters(tempModels, iter)
     {
         if (!iter().valid())
         {
@@ -167,18 +167,18 @@ void Foam::phaseSystem::generatePairsAndSubModels
         HashTable<autoPtr<modelType>, phasePairKey, phasePairKey::hash>
         modelTypeTable;
 
-    forAll(phaseModels_, phasei)
+    for (const phaseModel& phase : phaseModels_)
     {
         modelTypeTable tempModels;
         generatePairsAndSubModels
         (
-            IOobject::groupName(modelName, phaseModels_[phasei].name()),
+            IOobject::groupName(modelName, phase.name()),
             tempModels
         );
 
-        forAllConstIter(typename modelTypeTable, tempModels, tempModelIter)
+        forAllConstIters(tempModels, tempModelIter)
         {
-            const phasePairKey key(tempModelIter.key());
+            const phasePairKey& key = tempModelIter.key();
 
             if (!models.found(key))
             {
@@ -191,7 +191,7 @@ void Foam::phaseSystem::generatePairsAndSubModels
 
             models[tempModelIter.key()].insert
             (
-                phaseModels_[phasei].name(),
+                phase.name(),
                 *tempModelIter
             );
         }
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
index 40f8edfd4b5..95b301284c7 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
@@ -58,9 +58,9 @@ void Foam::multiphaseSystem::calcAlphas()
     scalar level = 0.0;
     alphas_ == 0.0;
 
-    forAll(phases(), i)
+    for (const phaseModel& phase : phases())
     {
-        alphas_ += level*phases()[i];
+        alphas_ += level * phase;
         level += 1.0;
     }
 }
@@ -70,15 +70,16 @@ void Foam::multiphaseSystem::solveAlphas()
 {
     bool LTS = fv::localEulerDdt::enabled(mesh_);
 
-    forAll(phases(), phasei)
+    for (phaseModel& phase : phases())
     {
-        phases()[phasei].correctBoundaryConditions();
+        phase.correctBoundaryConditions();
     }
 
     PtrList<surfaceScalarField> alphaPhiCorrs(phases().size());
-    forAll(phases(), phasei)
+
+    label phasei = 0;
+    for (phaseModel& phase : phases())
     {
-        phaseModel& phase = phases()[phasei];
         volScalarField& alpha1 = phase;
 
         alphaPhiCorrs.set
@@ -98,21 +99,18 @@ void Foam::multiphaseSystem::solveAlphas()
 
         surfaceScalarField& alphaPhiCorr = alphaPhiCorrs[phasei];
 
-        forAll(phases(), phasej)
+        for (const phaseModel& phase2 : phases())
         {
-            phaseModel& phase2 = phases()[phasej];
-            volScalarField& alpha2 = phase2;
+            const volScalarField& alpha2 = phase2;
 
             if (&phase2 == &phase) continue;
 
             surfaceScalarField phir(phase.phi() - phase2.phi());
 
-            cAlphaTable::const_iterator cAlpha
-            (
-                cAlphas_.find(phasePairKey(phase.name(), phase2.name()))
-            );
+            const auto cAlpha =
+                cAlphas_.cfind(phasePairKey(phase.name(), phase2.name()));
 
-            if (cAlpha != cAlphas_.end())
+            if (cAlpha.found())
             {
                 surfaceScalarField phic
                 (
@@ -122,7 +120,7 @@ void Foam::multiphaseSystem::solveAlphas()
                 phir += min(cAlpha()*phic, max(phic))*nHatf(phase, phase2);
             }
 
-            word phirScheme
+            const word phirScheme
             (
                 "div(phir," + alpha2.name() + ',' + alpha1.name() + ')'
             );
@@ -171,6 +169,8 @@ void Foam::multiphaseSystem::solveAlphas()
                 true
             );
         }
+
+        ++phasei;
     }
 
     MULES::limitSum(alphaPhiCorrs);
@@ -244,9 +244,8 @@ void Foam::multiphaseSystem::solveAlphas()
             }
         }
 
-        forAll(phases(), phasej)
+        for (const phaseModel& phase2 : phases())
         {
-            const phaseModel& phase2 = phases()[phasej];
             const volScalarField& alpha2 = phase2;
 
             if (&phase2 == &phase) continue;
@@ -303,9 +302,9 @@ void Foam::multiphaseSystem::solveAlphas()
 
     // Correct the sum of the phase-fractions to avoid 'drift'
     volScalarField sumCorr(1.0 - sumAlpha);
-    forAll(phases(), phasei)
+    for (phaseModel& phase : phases())
     {
-        volScalarField& alpha = phases()[phasei];
+        volScalarField& alpha = phase;
         alpha += alpha*sumCorr;
     }
 }
@@ -381,12 +380,11 @@ void Foam::multiphaseSystem::correctContactAngle
                /mesh_.magSf().boundaryField()[patchi]
             );
 
-            alphaContactAngleFvPatchScalarField::thetaPropsTable::
-                const_iterator tp =
-                    acap.thetaProps()
-                   .find(phasePairKey(phase1.name(), phase2.name()));
+            const auto tp =
+                acap.thetaProps()
+                    .cfind(phasePairKey(phase1.name(), phase2.name()));
 
-            if (tp == acap.thetaProps().end())
+            if (!tp.found())
             {
                 FatalErrorInFunction
                     << "Cannot find interface "
@@ -396,7 +394,7 @@ void Foam::multiphaseSystem::correctContactAngle
                     << exit(FatalError);
             }
 
-            bool matched = (tp.key().first() == phase1.name());
+            const bool matched = (tp.key().first() == phase1.name());
 
             const scalar theta0 = degToRad(tp().theta0(matched));
             scalarField theta(boundary[patchi].size(), theta0);
@@ -506,20 +504,14 @@ Foam::multiphaseSystem::multiphaseSystem
         1e-8/cbrt(average(mesh_.V()))
     )
 {
-    forAll(phases(), phasei)
+    for (const phaseModel& phase : phases())
     {
-        volScalarField& alphai = phases()[phasei];
-        mesh_.setFluxRequired(alphai.name());
+        const volScalarField& alpha = phase;
+        mesh_.setFluxRequired(alpha.name());
     }
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::multiphaseSystem::~multiphaseSystem()
-{}
-
-
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::surfaceScalarField> Foam::multiphaseSystem::surfaceTension
@@ -544,17 +536,15 @@ Foam::tmp<Foam::surfaceScalarField> Foam::multiphaseSystem::surfaceTension
 
     tSurfaceTension.ref().setOriented();
 
-    forAll(phases(), phasej)
+    for (const phaseModel& phase2 : phases())
     {
-        const phaseModel& phase2 = phases()[phasej];
-
         if (&phase2 != &phase1)
         {
             phasePairKey key12(phase1.name(), phase2.name());
 
-            cAlphaTable::const_iterator cAlpha(cAlphas_.find(key12));
+            const auto cAlpha = cAlphas_.cfind(key12);
 
-            if (cAlpha != cAlphas_.end())
+            if (cAlpha.found())
             {
                 tSurfaceTension.ref() +=
                     fvc::interpolate(sigma(key12)*K(phase1, phase2))
@@ -588,12 +578,12 @@ Foam::multiphaseSystem::nearInterface() const
         )
     );
 
-    forAll(phases(), phasei)
+    for (const phaseModel& phase : phases())
     {
         tnearInt.ref() = max
         (
             tnearInt(),
-            pos0(phases()[phasei] - 0.01)*pos0(0.99 - phases()[phasei])
+            pos0(phase - 0.01)*pos0(0.99 - phase)
         );
     }
 
@@ -623,9 +613,9 @@ void Foam::multiphaseSystem::solve()
         PtrList<volScalarField> alpha0s(phases().size());
         PtrList<surfaceScalarField> alphaPhiSums(phases().size());
 
-        forAll(phases(), phasei)
+        label phasei = 0;
+        for (phaseModel& phase : phases())
         {
-            phaseModel& phase = phases()[phasei];
             volScalarField& alpha = phase;
 
             alpha0s.set
@@ -649,6 +639,8 @@ void Foam::multiphaseSystem::solve()
                     dimensionedScalar(dimensionSet(0, 3, -1, 0, 0), Zero)
                 )
             );
+
+            ++phasei;
         }
 
         for
@@ -663,15 +655,18 @@ void Foam::multiphaseSystem::solve()
         {
             solveAlphas();
 
-            forAll(phases(), phasei)
+            phasei = 0;
+            for (const phaseModel& phase : phases())
             {
-                alphaPhiSums[phasei] += phases()[phasei].alphaPhi();
+                alphaPhiSums[phasei] += phase.alphaPhi();
+
+                ++phasei;
             }
         }
 
-        forAll(phases(), phasei)
+        phasei = 0;
+        for (phaseModel& phase : phases())
         {
-            phaseModel& phase = phases()[phasei];
             volScalarField& alpha = phase;
 
             phase.alphaPhi() = alphaPhiSums[phasei]/nAlphaSubCycles;
@@ -683,6 +678,8 @@ void Foam::multiphaseSystem::solve()
             // Reset the old-time field value
             alpha.oldTime() = alpha0s[phasei];
             alpha.oldTime().timeIndex() = runTime.timeIndex();
+
+            ++phasei;
         }
     }
     else
@@ -690,9 +687,8 @@ void Foam::multiphaseSystem::solve()
         solveAlphas();
     }
 
-    forAll(phases(), phasei)
+    for (phaseModel& phase : phases())
     {
-        phaseModel& phase = phases()[phasei];
         phase.alphaRhoPhi() = fvc::interpolate(phase.rho())*phase.alphaPhi();
 
         // Ensure the phase-fractions are bounded
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H
index 0c43f6dd8f1..53e7e96f72e 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H
@@ -140,11 +140,11 @@ public:
     // Constructors
 
         //- Construct from fvMesh
-        multiphaseSystem(const fvMesh&);
+        multiphaseSystem(const fvMesh& mesh);
 
 
     //- Destructor
-    virtual ~multiphaseSystem();
+    virtual ~multiphaseSystem() = default;
 
 
     // Selectors
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/blendingMethods/hyperbolic/hyperbolic.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/blendingMethods/hyperbolic/hyperbolic.C
index 9305f85c101..e577b35d6be 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/blendingMethods/hyperbolic/hyperbolic.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/blendingMethods/hyperbolic/hyperbolic.C
@@ -54,25 +54,22 @@ Foam::blendingMethods::hyperbolic::hyperbolic
     blendingMethod(dict),
     transitionAlphaScale_("transitionAlphaScale", dimless, dict)
 {
-    forAllConstIter(wordList, phaseNames, iter)
+    for (const word& phaseName : phaseNames)
     {
-        const word name(IOobject::groupName("maxDispersedAlpha", *iter));
-
         maxDispersedAlpha_.insert
         (
-            *iter,
-            dimensionedScalar(name, dimless, dict)
+            phaseName,
+            dimensionedScalar
+            (
+                IOobject::groupName("maxDispersedAlpha", phaseName),
+                dimless,
+                dict
+            )
         );
     }
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::blendingMethods::hyperbolic::~hyperbolic()
-{}
-
-
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::volScalarField> Foam::blendingMethods::hyperbolic::f1
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/blendingMethods/hyperbolic/hyperbolic.H b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/blendingMethods/hyperbolic/hyperbolic.H
index f518d7f0835..0d21b7396e7 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/blendingMethods/hyperbolic/hyperbolic.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/blendingMethods/hyperbolic/hyperbolic.H
@@ -77,7 +77,7 @@ public:
 
 
     //- Destructor
-    ~hyperbolic();
+    ~hyperbolic() = default;
 
 
     // Member Functions
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/blendingMethods/linear/linear.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/blendingMethods/linear/linear.C
index 1880e2be15b..adbe2c1929f 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/blendingMethods/linear/linear.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/blendingMethods/linear/linear.C
@@ -53,39 +53,39 @@ Foam::blendingMethods::linear::linear
 :
     blendingMethod(dict)
 {
-    forAllConstIter(wordList, phaseNames, iter)
+    for (const word& phaseName : phaseNames)
     {
-        const word nameFull
-        (
-            IOobject::groupName("maxFullyDispersedAlpha", *iter)
-        );
-
         maxFullyDispersedAlpha_.insert
         (
-            *iter,
-            dimensionedScalar(nameFull, dimless, dict)
-        );
-
-        const word namePart
-        (
-            IOobject::groupName("maxPartlyDispersedAlpha", *iter)
+            phaseName,
+            dimensionedScalar
+            (
+                IOobject::groupName("maxFullyDispersedAlpha", phaseName),
+                dimless,
+                dict
+            )
         );
 
         maxPartlyDispersedAlpha_.insert
         (
-            *iter,
-            dimensionedScalar(namePart, dimless, dict)
+            phaseName,
+            dimensionedScalar
+            (
+                IOobject::groupName("maxPartlyDispersedAlpha", phaseName),
+                dimless,
+                dict
+            )
         );
 
         if
         (
-            maxFullyDispersedAlpha_[*iter]
-          > maxPartlyDispersedAlpha_[*iter]
+            maxFullyDispersedAlpha_[phaseName]
+          > maxPartlyDispersedAlpha_[phaseName]
         )
         {
             FatalErrorInFunction
                 << "The supplied fully dispersed volume fraction for "
-                << *iter
+                << phaseName
                 << " is greater than the partly dispersed value."
                 << endl << exit(FatalError);
         }
@@ -93,12 +93,6 @@ Foam::blendingMethods::linear::linear
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::blendingMethods::linear::~linear()
-{}
-
-
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::volScalarField> Foam::blendingMethods::linear::f1
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/blendingMethods/linear/linear.H b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/blendingMethods/linear/linear.H
index 54dcc73ecf3..2b6249a6f33 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/blendingMethods/linear/linear.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/blendingMethods/linear/linear.H
@@ -77,7 +77,7 @@ public:
 
 
     //- Destructor
-    ~linear();
+    ~linear() = default;
 
 
     // Member Functions
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
index b5796197add..b73a894a62c 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
@@ -561,10 +561,8 @@ bool Foam::twoPhaseSystem::read()
 
         return readOK;
     }
-    else
-    {
-        return false;
-    }
+
+    return false;
 }
 
 
-- 
GitLab