diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.C
index 593ed919ab..86499a7827 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.C
@@ -46,20 +46,23 @@ PopulationBalancePhaseSystem
         const diameterModels::populationBalanceModel& popBal =
             populationBalances_[popBali];
 
-        forAll(popBal.velocityGroups(), velGrp1i)
+        forAllConstIter
+        (
+            HashTable<const diameterModels::velocityGroup*>,
+            popBal.velocityGroupPtrs(),
+            iter1
+        )
         {
-            const diameterModels::velocityGroup& velGrp1 =
-                popBal.velocityGroups()[velGrp1i];
+            const diameterModels::velocityGroup& velGrp1 = *iter1();
 
-            for
+            forAllConstIter
             (
-                label velGrp2i = velGrp1i + 1;
-                velGrp2i < popBal.velocityGroups().size();
-                ++ velGrp2i
+                HashTable<const diameterModels::velocityGroup*>,
+                popBal.velocityGroupPtrs(),
+                iter2
             )
             {
-                const diameterModels::velocityGroup& velGrp2 =
-                    popBal.velocityGroups()[velGrp2i];
+                const diameterModels::velocityGroup& velGrp2 = *iter2();
 
                 const phaseInterface interface
                 (
@@ -67,28 +70,36 @@ PopulationBalancePhaseSystem
                     velGrp2.phase()
                 );
 
-                this->template validateMassTransfer
-                    <diameterModels::populationBalanceModel>(interface);
-
-                dmdtfs_.insert
+                if
                 (
-                    interface,
-                    new volScalarField
+                    &velGrp1 != &velGrp2
+                    &&
+                    !dmdtfs_.found(interface)
+                )
+                {
+                    this->template validateMassTransfer
+                        <diameterModels::populationBalanceModel>(interface);
+
+                    dmdtfs_.insert
                     (
-                        IOobject
+                        interface,
+                        new volScalarField
                         (
-                            IOobject::groupName
+                            IOobject
                             (
-                                "populationBalance:dmdtf",
-                                interface.name()
+                                IOobject::groupName
+                                (
+                                    "populationBalance:dmdtf",
+                                    interface.name()
+                                ),
+                                this->mesh().time().timeName(),
+                                this->mesh()
                             ),
-                            this->mesh().time().timeName(),
-                            this->mesh()
-                        ),
-                        this->mesh(),
-                        dimensionedScalar(dimDensity/dimTime, 0)
-                    )
-                );
+                            this->mesh(),
+                            dimensionedScalar(dimDensity/dimTime, 0)
+                        )
+                    );
+                }
             }
         }
     }
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/shapeModels/fractal/fractal.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/shapeModels/fractal/fractal.C
index e56019fe73..b10e1f363c 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/shapeModels/fractal/fractal.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/shapeModels/fractal/fractal.C
@@ -26,9 +26,9 @@ License
 #include "fractal.H"
 #include "addToRunTimeSelectionTable.H"
 #include "sinteringModel.H"
-#include "fvm.H"
-#include "fvcDdt.H"
-#include "fvcDiv.H"
+#include "fvmDdt.H"
+#include "fvmDiv.H"
+#include "fvmSup.H"
 #include "mixedFvPatchField.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -122,7 +122,7 @@ Foam::diameterModels::shapeModels::fractal::fractal
 
         if
         (
-            isA<mixedFvPatchScalarField>(kappa_.boundaryFieldRef()[patchi])
+            isA<const mixedFvPatchScalarField>(kappa_.boundaryField()[patchi])
         )
         {
             mixedFvPatchScalarField& kappa =
@@ -205,12 +205,12 @@ void Foam::diameterModels::shapeModels::fractal::correct()
 
     fvScalarMatrix kappaEqn
     (
-        fvm::ddt(alpha, fi, kappa_) + fvm::div(fAlphaPhi, kappa_)
+        fvm::ddt(alpha, fi, kappa_)
+      + fvm::div(fAlphaPhi, kappa_)
       ==
       - sinteringModel_->R()
       + Su_
-      - fvm::Sp(popBal.Sp(fi.i()())*fi, kappa_)
-
+      - fvm::Sp(popBal.Sp(fi.i())*fi, kappa_)
       - correction
         (
             fvm::Sp
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/sizeGroup.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/sizeGroup.C
index 7db66a5957..0f5a8b3477 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/sizeGroup.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/sizeGroup.C
@@ -75,7 +75,7 @@ Foam::diameterModels::sizeGroup::sizeGroup
 
         if
         (
-            isA<mixedFvPatchScalarField>(this->boundaryFieldRef()[patchi])
+            isA<const mixedFvPatchScalarField>(this->boundaryField()[patchi])
         )
         {
             mixedFvPatchScalarField& f =
@@ -108,7 +108,7 @@ Foam::diameterModels::sizeGroup::clone() const
 }
 
 
-const Foam::autoPtr<Foam::label>& Foam::diameterModels::sizeGroup::i() const
+const Foam::label& Foam::diameterModels::sizeGroup::i() const
 {
     if (!i_.valid())
     {
@@ -127,7 +127,7 @@ const Foam::autoPtr<Foam::label>& Foam::diameterModels::sizeGroup::i() const
         }
     }
 
-    return i_;
+    return i_();
 }
 
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/sizeGroup.H b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/sizeGroup.H
index 1fd317d8e4..43121086eb 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/sizeGroup.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/sizeGroup.H
@@ -192,8 +192,8 @@ public:
         //- Return const-reference to diameterModel of the phase
         inline const autoPtr<shapeModel>& shapeModelPtr() const;
 
-        //- Return label of the sizeGroup within the population balance
-        const autoPtr<label>& i() const;
+        //- Return index of the size group within the population balance
+        const label& i() const;
 
         //- Return representative surface area of the sizeGroup
         const tmp<volScalarField> a() const;
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/driftModels/constantDrift/constantDrift.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/driftModels/constantDrift/constantDrift.C
index de960a56b0..562dfe4f31 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/driftModels/constantDrift/constantDrift.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/driftModels/constantDrift/constantDrift.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2017-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2017-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -63,7 +63,7 @@ void Foam::diameterModels::driftModels::constantDrift::addToDriftRate
     const label i
 )
 {
-    driftRate += rate_;
+    driftRate += popBal_.sizeGroups()[i].phase()*rate_;
 }
 
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/driftModels/densityChange/densityChange.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/driftModels/densityChange/densityChange.C
index e3da300864..d01d7344f0 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/driftModels/densityChange/densityChange.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/driftModels/densityChange/densityChange.C
@@ -70,11 +70,11 @@ void Foam::diameterModels::driftModels::densityChangeDrift::addToDriftRate
     const volScalarField& rho = phase.thermo().rho();
 
     driftRate -=
-        fi.x()/(rho*max(alpha, phase.residualAlpha()))
+        fi.x()
        *(
             fvc::ddt(alpha, rho) + fvc::div(phase.alphaRhoPhi())
           - fvc::Sp(fvc::ddt(alpha) + fvc::div(phase.alphaPhi()), rho)
-        );
+        )/rho;
 }
 
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/driftModels/phaseChange/phaseChange.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/driftModels/phaseChange/phaseChange.C
index 940096ce33..6e73374ba7 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/driftModels/phaseChange/phaseChange.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/driftModels/phaseChange/phaseChange.C
@@ -97,15 +97,20 @@ void Foam::diameterModels::driftModels::phaseChange::precompute()
 
     forAll(interfaces_, k)
     {
-        forAll(popBal_.velocityGroups(), j)
+        forAllConstIter
+        (
+            HashTable<const velocityGroup*>,
+            popBal_.velocityGroupPtrs(),
+            iter
+        )
         {
-            const velocityGroup& vgj = popBal_.velocityGroups()[j];
+            const velocityGroup& velGrp = *iter();
 
-            if (interfaces_[k].contains(vgj.phase()))
+            if (interfaces_[k].contains(velGrp.phase()))
             {
-                forAll(vgj.sizeGroups(), i)
+                forAll(velGrp.sizeGroups(), i)
                 {
-                    const sizeGroup& fi = vgj.sizeGroups()[i];
+                    const sizeGroup& fi = velGrp.sizeGroups()[i];
 
                     if (numberWeighted_)
                     {
@@ -128,11 +133,11 @@ void Foam::diameterModels::driftModels::phaseChange::addToDriftRate
     const label i
 )
 {
-    const velocityGroup& vg = popBal_.sizeGroups()[i].VelocityGroup();
+    const velocityGroup& velGrp = popBal_.sizeGroups()[i].VelocityGroup();
 
     forAll(interfaces_, k)
     {
-        if (interfaces_[k].contains(vg.phase()))
+        if (interfaces_[k].contains(velGrp.phase()))
         {
             const volScalarField& dmidtf =
                 popBal_.mesh().lookupObject<volScalarField>
@@ -149,7 +154,7 @@ void Foam::diameterModels::driftModels::phaseChange::addToDriftRate
                 );
 
             const scalar dmidtfSign =
-                interfaces_[k].index(vg.phase()) == 0 ? +1 : -1;
+                interfaces_[k].index(velGrp.phase()) == 0 ? +1 : -1;
 
             const sizeGroup& fi = popBal_.sizeGroups()[i];
 
@@ -163,7 +168,7 @@ void Foam::diameterModels::driftModels::phaseChange::addToDriftRate
                 dDriftRate.ref() *= fi.a();
             }
 
-            driftRate += dDriftRate;
+            driftRate += velGrp.phase()*dDriftRate;
         }
     }
 }
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C
index e442a878eb..e4f7116dbe 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C
@@ -31,9 +31,12 @@ License
 #include "driftModel.H"
 #include "nucleationModel.H"
 #include "surfaceTensionModel.H"
-#include "fvm.H"
+#include "fvmDdt.H"
+#include "fvmDiv.H"
+#include "fvmSup.H"
 #include "fvcDdt.H"
 #include "fvcDiv.H"
+#include "fvcSup.H"
 #include "shapeModel.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -60,12 +63,26 @@ void Foam::diameterModels::populationBalanceModel::registerVelocityGroups()
 
             if (velGroup.popBalName() == this->name())
             {
-                velocityGroups_.resize(velocityGroups_.size() + 1);
+                velocityGroupPtrs_.insert(velGroup.phase().name(), &velGroup);
 
-                velocityGroups_.set
+                dilatationErrors_.insert
                 (
-                    velocityGroups_.size() - 1,
-                    &const_cast<velocityGroup&>(velGroup)
+                    velGroup.phase().name(),
+                    volScalarField
+                    (
+                        IOobject
+                        (
+                            IOobject::groupName
+                            (
+                                "dilatationError",
+                                velGroup.phase().name()
+                            ),
+                            fluid_.time().timeName(),
+                            mesh_
+                        ),
+                        mesh_,
+                        dimensionedScalar(inv(dimTime), 0)
+                    )
                 );
 
                 forAll(velGroup.sizeGroups(), i)
@@ -446,22 +463,21 @@ void Foam::diameterModels::populationBalanceModel::drift
         const sizeGroup& fe = sizeGroups()[i+1];
         volScalarField& Sue = Sui_;
 
-        Sp_[i] += 1/(fe.x() - fp.x())*pos(driftRate_())*driftRate_()*fp.phase();
+        Sp_[i] += 1/(fe.x() - fp.x())*pos(driftRate_())*driftRate_();
 
         Sue =
-            fe.x()/(fp.x()*(fe.x() - fp.x()))
-           *pos(driftRate_())*driftRate_()*fp*fp.phase();
+            fe.x()/(fp.x()*(fe.x() - fp.x()))*pos(driftRate_())*driftRate_()*fp;
 
         Su_[i+1] += Sue;
 
-        const phaseInterface interfaceij(fp.phase(), fe.phase());
+        const phaseInterface interfacepe(fp.phase(), fe.phase());
 
-        if (pDmdt_.found(interfaceij))
+        if (pDmdt_.found(interfacepe))
         {
             const scalar dmdtSign =
-                interfaceij.index(fp.phase()) == 0 ? +1 : -1;
+                interfacepe.index(fp.phase()) == 0 ? +1 : -1;
 
-            *pDmdt_[interfaceij] -= dmdtSign*Sue*fp.phase().rho();
+            *pDmdt_[interfacepe] -= dmdtSign*Sue*fp.phase().rho();
         }
 
         sizeGroups_[i+1].shapeModelPtr()->addDrift(Sue, fp, model);
@@ -469,7 +485,7 @@ void Foam::diameterModels::populationBalanceModel::drift
 
     if (i == sizeGroups().size() - 1)
     {
-        Sp_[i] -= pos(driftRate_())*driftRate_()*fp.phase()/fp.x();
+        Sp_[i] -= pos(driftRate_())*driftRate_()/fp.x();
     }
 
     if (i > 0)
@@ -477,22 +493,21 @@ void Foam::diameterModels::populationBalanceModel::drift
         const sizeGroup& fw = sizeGroups()[i-1];
         volScalarField& Suw = Sui_;
 
-        Sp_[i] += 1/(fw.x() - fp.x())*neg(driftRate_())*driftRate_()*fp.phase();
+        Sp_[i] += 1/(fw.x() - fp.x())*neg(driftRate_())*driftRate_();
 
         Suw =
-            fw.x()/(fp.x()*(fw.x() - fp.x()))
-           *neg(driftRate_())*driftRate_()*fp*fp.phase();
+            fw.x()/(fp.x()*(fw.x() - fp.x()))*neg(driftRate_())*driftRate_()*fp;
 
         Su_[i-1] += Suw;
 
-        const phaseInterface interfaceih(fp.phase(), fw.phase());
+        const phaseInterface interfacepw(fp.phase(), fw.phase());
 
-        if (pDmdt_.found(interfaceih))
+        if (pDmdt_.found(interfacepw))
         {
             const scalar dmdtSign =
-                interfaceih.index(fp.phase()) == 0 ? +1 : -1;
+                interfacepw.index(fp.phase()) == 0 ? +1 : -1;
 
-            *pDmdt_[interfaceih] -= dmdtSign*Suw*fp.phase().rho();
+            *pDmdt_[interfacepw] -= dmdtSign*Suw*fp.phase().rho();
         }
 
         sizeGroups_[i-1].shapeModelPtr()->addDrift(Suw, fp, model);
@@ -500,7 +515,7 @@ void Foam::diameterModels::populationBalanceModel::drift
 
     if (i == 0)
     {
-        Sp_[i] -= neg(driftRate_())*driftRate_()*fp.phase()/fp.x();
+        Sp_[i] -= neg(driftRate_())*driftRate_()/fp.x();
     }
 }
 
@@ -609,13 +624,43 @@ void Foam::diameterModels::populationBalanceModel::sources()
 }
 
 
+void Foam::diameterModels::populationBalanceModel::correctDilatationError()
+{
+    forAllIter
+    (
+        HashTable<volScalarField>,
+        dilatationErrors_,
+        iter
+    )
+    {
+        volScalarField& dilatationError = iter();
+        const word& phaseName = iter.key();
+        const phaseModel& phase = fluid_.phases()[phaseName];
+        const velocityGroup& velGroup = *velocityGroupPtrs_[phaseName];
+        const volScalarField& alpha = phase;
+        const volScalarField& rho = phase.thermo().rho();
+
+        dilatationError =
+            fvc::ddt(alpha) + fvc::div(phase.alphaPhi())
+          - (fluid_.fvModels().source(alpha, rho) & rho)/rho;
+
+        forAll(velGroup.sizeGroups(), i)
+        {
+            const sizeGroup& fi = velGroup.sizeGroups()[i];
+
+            dilatationError -= Su_[fi.i()] - fvc::Sp(Sp_[fi.i()], fi);
+        }
+    }
+}
+
+
 void Foam::diameterModels::populationBalanceModel::calcAlphas()
 {
     alphas_() = Zero;
 
-    forAll(velocityGroups_, v)
+    forAllConstIter(HashTable<const velocityGroup*>, velocityGroupPtrs_, iter)
     {
-        const phaseModel& phase = velocityGroups_[v].phase();
+        const phaseModel& phase = iter()->phase();
 
         alphas_() += max(phase, phase.residualAlpha());
     }
@@ -637,9 +682,9 @@ Foam::diameterModels::populationBalanceModel::calcDsm()
 
     volScalarField& invDsm = tInvDsm.ref();
 
-    forAll(velocityGroups_, v)
+    forAllConstIter(HashTable<const velocityGroup*>, velocityGroupPtrs_, iter)
     {
-        const phaseModel& phase = velocityGroups_[v].phase();
+        const phaseModel& phase = iter()->phase();
 
         invDsm += max(phase, phase.residualAlpha())/(phase.d()*alphas_());
     }
@@ -652,9 +697,9 @@ void Foam::diameterModels::populationBalanceModel::calcVelocity()
 {
     U_() = Zero;
 
-    forAll(velocityGroups_, v)
+    forAllConstIter(HashTable<const velocityGroup*>, velocityGroupPtrs_, iter)
     {
-        const phaseModel& phase = velocityGroups_[v].phase();
+        const phaseModel& phase = iter()->phase();
 
         U_() += max(phase, phase.residualAlpha())*phase.U()/alphas_();
     }
@@ -704,7 +749,6 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel
             IOobject::groupName("alpha", dict_.lookup("continuousPhase"))
         )
     ),
-    velocityGroups_(),
     sizeGroups_(),
     v_(),
     delta_(),
@@ -890,7 +934,7 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel
         );
     }
 
-    if (velocityGroups_.size() > 1)
+    if (velocityGroupPtrs_.size() > 1)
     {
         alphas_.set
         (
@@ -1060,6 +1104,8 @@ void Foam::diameterModels::populationBalanceModel::solve()
                 sources();
             }
 
+            correctDilatationError();
+
             maxInitialResidual = 0;
 
             forAll(sizeGroups(), i)
@@ -1068,15 +1114,18 @@ void Foam::diameterModels::populationBalanceModel::solve()
                 const phaseModel& phase = fi.phase();
                 const volScalarField& alpha = phase;
                 const volScalarField& rho = phase.thermo().rho();
+                const volScalarField& dilatationError =
+                    dilatationErrors_[phase.name()];
 
                 fvScalarMatrix sizeGroupEqn
                 (
-                    fvm::ddt(alpha, fi) + fvm::div(phase.alphaPhi(), fi)
+                    fvm::ddt(alpha, fi)
+                  + fvm::div(phase.alphaPhi(), fi)
+                  - fvm::Sp(dilatationError, fi)
                 ==
-                    Su_[i]
+                    fvc::Su(Su_[i], fi)
                   - fvm::Sp(Sp_[i], fi)
                   + fluid_.fvModels().source(alpha, rho, fi)/rho
-
                   - correction
                     (
                         fvm::Sp
@@ -1121,7 +1170,7 @@ void Foam::diameterModels::populationBalanceModel::solve()
 
 void Foam::diameterModels::populationBalanceModel::correct()
 {
-    if (velocityGroups_.size() > 1)
+    if (velocityGroupPtrs_.size() > 1)
     {
         calcAlphas();
         dsm_() = calcDsm();
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.H b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.H
index 17bca4f326..62ea1f1b83 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.H
@@ -257,10 +257,10 @@ private:
         //- Continuous phase
         const phaseModel& continuousPhase_;
 
-        //- velocityGroups belonging to this populationBalance
-        UPtrList<velocityGroup> velocityGroups_;
+        //- Velocity groups belonging to this populationBalance
+        HashTable<const velocityGroup*> velocityGroupPtrs_;
 
-        //- sizeGroups belonging to this populationBalance
+        //- Size groups belonging to this populationBalance
         UPtrList<sizeGroup> sizeGroups_;
 
         //- sizeGroup boundaries
@@ -275,6 +275,9 @@ private:
         //- Implicitly treated sources
         PtrList<volScalarField> Sp_;
 
+        //- Dilatation errors
+        HashTable<volScalarField> dilatationErrors_;
+
         //- Field for caching sources
         volScalarField Sui_;
 
@@ -357,7 +360,7 @@ private:
 
         void sources();
 
-        void dmdt();
+        void correctDilatationError();
 
         void calcAlphas();
 
@@ -454,15 +457,18 @@ public:
         //- Return continuous phase
         inline const phaseModel& continuousPhase() const;
 
-        //- Return the velocityGroups belonging to this populationBalance
-        inline const UPtrList<velocityGroup>& velocityGroups() const;
+        //- Return the velocity groups belonging to this populationBalance
+        inline const HashTable<const velocityGroup*>& velocityGroupPtrs() const;
 
-        //- Return the sizeGroups belonging to this populationBalance
+        //- Return the size groups belonging to this populationBalance
         inline const UPtrList<sizeGroup>& sizeGroups() const;
 
         //- Return implicit source terms
         inline const volScalarField& Sp(const label i) const;
 
+        //- Return dilatation obtained from source terms
+        inline const HashTable<volScalarField>& sourceDilatation() const;
+
         //- Return coalescence relevant size group pairs
         inline const List<labelPair>& coalescencePairs() const;
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModelI.H b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModelI.H
index f38dd139e8..48b5c8be46 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModelI.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModelI.H
@@ -82,10 +82,10 @@ Foam::diameterModels::populationBalanceModel::continuousPhase() const
 }
 
 
-inline const Foam::UPtrList<Foam::diameterModels::velocityGroup>&
-Foam::diameterModels::populationBalanceModel::velocityGroups() const
+inline const Foam::HashTable<const Foam::diameterModels::velocityGroup*>&
+Foam::diameterModels::populationBalanceModel::velocityGroupPtrs() const
 {
-    return velocityGroups_;
+    return velocityGroupPtrs_;
 }
 
 
@@ -120,13 +120,13 @@ Foam::diameterModels::populationBalanceModel::binaryBreakupPairs() const
 inline const Foam::volScalarField&
 Foam::diameterModels::populationBalanceModel::alphas() const
 {
-    if (velocityGroups_.size() > 1)
+    if (velocityGroupPtrs_.size() > 1)
     {
         return alphas_();
     }
     else
     {
-        return velocityGroups_.first().phase();
+        return velocityGroupPtrs_[velocityGroupPtrs_.begin().key()]->phase();
     }
 }
 
@@ -134,13 +134,14 @@ Foam::diameterModels::populationBalanceModel::alphas() const
 inline const Foam::volVectorField&
 Foam::diameterModels::populationBalanceModel::U() const
 {
-    if (velocityGroups_.size() > 1)
+    if (velocityGroupPtrs_.size() > 1)
     {
         return U_();
     }
     else
     {
-        return velocityGroups_.first().phase().U();
+        return
+            velocityGroupPtrs_[velocityGroupPtrs_.begin().key()]->phase().U();
     }
 }
 
diff --git a/test/multiphase/multiphaseEulerFoam/populationBalance/drift/constant/fvModels b/test/multiphase/multiphaseEulerFoam/populationBalance/drift/constant/fvModels
index 941923760a..cc8a2e57e6 100644
--- a/test/multiphase/multiphaseEulerFoam/populationBalance/drift/constant/fvModels
+++ b/test/multiphase/multiphaseEulerFoam/populationBalance/drift/constant/fvModels
@@ -14,51 +14,6 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-airSource
-{
-    type            massSource;
-
-    selectionMode   all;
-
-    massFlowRate    0.02186682;
-
-    phase           air;
-    rho             thermo:rho.air;
-
-    fieldValues
-    {
-        f0.air      0;
-        f1.air      0;
-        f2.air      0;
-        f3.air      0;
-        f4.air      0;
-        f5.air      0;
-        f6.air      0;
-        f7.air      0;
-        f8.air      0;
-        f9.air      0;
-        f10.air     0;
-        f11.air     0;
-        f12.air     0;
-        f13.air     0;
-        f14.air     0;
-        f15.air     0;
-        f16.air     0;
-        f17.air     0;
-        f18.air     0;
-        f19.air     0;
-        f20.air     0;
-        f21.air     0;
-        f22.air     0;
-        f23.air     0;
-        f24.air     0;
-        f25.air     0;
-        f26.air     0;
-        f27.air     0;
-        f28.air     0;
-    }
-}
-
 waterSink
 {
     type            massSource;
@@ -71,7 +26,9 @@ waterSink
     rho             thermo:rho.water;
 
     fieldValues
-    {}
+    {
+        U.water     (0 0 0);
+    }
 }
 
 // ************************************************************************* //
diff --git a/test/multiphase/multiphaseEulerFoam/populationBalance/drift/constant/physicalProperties.air b/test/multiphase/multiphaseEulerFoam/populationBalance/drift/constant/physicalProperties.air
index 2c9bc0082d..ec67d27ce9 100644
--- a/test/multiphase/multiphaseEulerFoam/populationBalance/drift/constant/physicalProperties.air
+++ b/test/multiphase/multiphaseEulerFoam/populationBalance/drift/constant/physicalProperties.air
@@ -20,7 +20,7 @@ thermoType
     mixture         pureMixture;
     transport       const;
     thermo          hConst;
-    equationOfState rhoConst;
+    equationOfState perfectGas;
     specie          specie;
     energy          sensibleInternalEnergy;
 }
@@ -31,10 +31,6 @@ mixture
     {
         molWeight   28.9;
     }
-    equationOfState
-    {
-        rho         1;
-    }
     thermodynamics
     {
         Cp          1007;
diff --git a/test/multiphase/multiphaseEulerFoam/populationBalance/drift/system/fvSolution b/test/multiphase/multiphaseEulerFoam/populationBalance/drift/system/fvSolution
index 5f62647fd7..bab15a15b5 100644
--- a/test/multiphase/multiphaseEulerFoam/populationBalance/drift/system/fvSolution
+++ b/test/multiphase/multiphaseEulerFoam/populationBalance/drift/system/fvSolution
@@ -60,7 +60,7 @@ PIMPLE
     pRefCell        0;
     pRefValue       1e5;
 
-    flow            no;
+    flow            yes;
     thermophysics   yes;
     models          yes;
 }
diff --git a/test/multiphase/multiphaseEulerFoam/populationBalance/isothermalGrowth/validation/createGraphs b/test/multiphase/multiphaseEulerFoam/populationBalance/isothermalGrowth/validation/createGraphs
index 7f97180f8d..843ccf6778 100755
--- a/test/multiphase/multiphaseEulerFoam/populationBalance/isothermalGrowth/validation/createGraphs
+++ b/test/multiphase/multiphaseEulerFoam/populationBalance/isothermalGrowth/validation/createGraphs
@@ -12,7 +12,6 @@ gnuplot<<EOF
 
     set xlabel 'd (mm)'
     set ylabel 'y (m)'
-    set yrange [0.05:]
     set key bottom center
 
     latestTime = system("foamListTimes -case .. -latestTime")
diff --git a/test/multiphase/multiphaseEulerFoam/populationBalance/negativeDrift/constant/fvModels b/test/multiphase/multiphaseEulerFoam/populationBalance/negativeDrift/constant/fvModels
index c165b0f115..81d711b964 100644
--- a/test/multiphase/multiphaseEulerFoam/populationBalance/negativeDrift/constant/fvModels
+++ b/test/multiphase/multiphaseEulerFoam/populationBalance/negativeDrift/constant/fvModels
@@ -14,52 +14,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-airSource
-{
-    type            massSource;
-
-    selectionMode   all;
-
-    massFlowRate    -0.030785654;
-
-    phase           air;
-    rho             thermo:rho.air;
-
-    fieldValues
-    {
-        f0.air      0;
-        f1.air      0;
-        f2.air      0;
-        f3.air      0;
-        f4.air      0;
-        f5.air      0;
-        f6.air      0;
-        f7.air      0;
-        f8.air      0;
-        f9.air      0;
-        f10.air     0;
-        f11.air     0;
-        f12.air     0;
-        f13.air     0;
-        f14.air     0;
-        f15.air     0;
-        f16.air     0;
-        f17.air     0;
-        f18.air     0;
-        f19.air     0;
-        f20.air     0;
-        f21.air     0;
-        f22.air     0;
-        f23.air     0;
-        f24.air     0;
-        f25.air     0;
-        f26.air     0;
-        f27.air     0;
-        f28.air     0;
-    }
-}
-
-waterSink
+waterSource
 {
     type            massSource;
 
@@ -71,7 +26,9 @@ waterSink
     rho             thermo:rho.water;
 
     fieldValues
-    {}
+    {
+        U.water     (0 0 0);
+    }
 }
 
 // ************************************************************************* //
diff --git a/test/multiphase/multiphaseEulerFoam/populationBalance/negativeDrift/constant/physicalProperties.air b/test/multiphase/multiphaseEulerFoam/populationBalance/negativeDrift/constant/physicalProperties.air
index 2c9bc0082d..ec67d27ce9 100644
--- a/test/multiphase/multiphaseEulerFoam/populationBalance/negativeDrift/constant/physicalProperties.air
+++ b/test/multiphase/multiphaseEulerFoam/populationBalance/negativeDrift/constant/physicalProperties.air
@@ -20,7 +20,7 @@ thermoType
     mixture         pureMixture;
     transport       const;
     thermo          hConst;
-    equationOfState rhoConst;
+    equationOfState perfectGas;
     specie          specie;
     energy          sensibleInternalEnergy;
 }
@@ -31,10 +31,6 @@ mixture
     {
         molWeight   28.9;
     }
-    equationOfState
-    {
-        rho         1;
-    }
     thermodynamics
     {
         Cp          1007;
diff --git a/test/multiphase/multiphaseEulerFoam/populationBalance/negativeDrift/system/fvSolution b/test/multiphase/multiphaseEulerFoam/populationBalance/negativeDrift/system/fvSolution
index 5f62647fd7..bab15a15b5 100644
--- a/test/multiphase/multiphaseEulerFoam/populationBalance/negativeDrift/system/fvSolution
+++ b/test/multiphase/multiphaseEulerFoam/populationBalance/negativeDrift/system/fvSolution
@@ -60,7 +60,7 @@ PIMPLE
     pRefCell        0;
     pRefValue       1e5;
 
-    flow            no;
+    flow            yes;
     thermophysics   yes;
     models          yes;
 }
diff --git a/tutorials/multiphase/multiphaseEulerFoam/RAS/wallBoilingPolydisperse/system/fvSolution b/tutorials/multiphase/multiphaseEulerFoam/RAS/wallBoilingPolydisperse/system/fvSolution
index 772b122868..abc735062b 100644
--- a/tutorials/multiphase/multiphaseEulerFoam/RAS/wallBoilingPolydisperse/system/fvSolution
+++ b/tutorials/multiphase/multiphaseEulerFoam/RAS/wallBoilingPolydisperse/system/fvSolution
@@ -28,7 +28,7 @@ solvers
         tolerance               1e-4;
         scale                   true;
         solveOnFinalIterOnly    true;
-        sourceUpdateInterval    1;
+        sourceUpdateInterval    20;
     }
 
     p_rgh
diff --git a/tutorials/multiphase/multiphaseEulerFoam/RAS/wallBoilingPolydisperseTwoGroups/system/fvSolution b/tutorials/multiphase/multiphaseEulerFoam/RAS/wallBoilingPolydisperseTwoGroups/system/fvSolution
index 772b122868..abc735062b 100644
--- a/tutorials/multiphase/multiphaseEulerFoam/RAS/wallBoilingPolydisperseTwoGroups/system/fvSolution
+++ b/tutorials/multiphase/multiphaseEulerFoam/RAS/wallBoilingPolydisperseTwoGroups/system/fvSolution
@@ -28,7 +28,7 @@ solvers
         tolerance               1e-4;
         scale                   true;
         solveOnFinalIterOnly    true;
-        sourceUpdateInterval    1;
+        sourceUpdateInterval    20;
     }
 
     p_rgh