diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/UEqns.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/UEqns.H
index c513c24266..6e85586838 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/UEqns.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/UEqns.H
@@ -60,22 +60,16 @@ PtrList<fvVectorMatrix> UEqns(phases.size());
 
             const phasePair& pair(fluid.phasePairs()[VmIter.key()]);
 
-            const phaseModel* phase1 = &pair.phase1();
-            const phaseModel* phase2 = &pair.phase2();
-
-            forAllConstIter(phasePair, pair, iter)
+            if (pair.contains(phase))
             {
-                if (phase1 == &phase)
-                {
-                    UEqns[phasei] +=
-                        Vm
-                       *(
-                            UgradUs[phasei]
-                          - (UgradUs[phase2->index()] & phase2->U())
-                        );
-                }
+                const phaseModel& otherPhase = pair.other(phase);
 
-                Swap(phase1, phase2);
+                UEqns[phasei] +=
+                    Vm
+                   *(
+                        UgradUs[phasei]
+                      - (UgradUs[otherPhase.index()] & otherPhase.U())
+                    );
             }
         }
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/pEqn.H
index e920119b9a..b211c71f02 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/pEqn.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/pEqn.H
@@ -194,18 +194,10 @@ while (pimple.correct())
 
             const phasePair& pair(fluid.phasePairs()[KdfIter.key()]);
 
-            const phaseModel* phase1 = &pair.phase1();
-            const phaseModel* phase2 = &pair.phase2();
-
-            forAllConstIter(phasePair, pair, iter)
+            if (pair.contains(phase))
             {
-                if (phase1 == &phase)
-                {
-                    phiHbyAs[phasei] +=
-                        rAUfs[phasei]*Kf*MRF.absolute(phase2->phi());
-                }
-
-                Swap(phase1, phase2);
+                phiHbyAs[phasei] +=
+                    rAUfs[phasei]*Kf*MRF.absolute(pair.other(phase).phi());
             }
         }
 
@@ -221,22 +213,14 @@ while (pimple.correct())
 
             const phasePair& pair(fluid.phasePairs()[VmfIter.key()]);
 
-            const phaseModel* phase1 = &pair.phase1();
-            const phaseModel* phase2 = &pair.phase2();
-
-            forAllConstIter(phasePair, pair, iter)
+            if (pair.contains(phase))
             {
-                if (phase1 == &phase)
-                {
-                    phiHbyAs[phasei] +=
-                        rAUfs[phasei]
-                       *(
-                            Vmf*byDt(MRF.absolute(phase1->phi()().oldTime()))
-                          + Vmf*ddtPhis[phase2->index()]
-                        );
-                }
-
-                Swap(phase1, phase2);
+                phiHbyAs[phasei] +=
+                    rAUfs[phasei]
+                   *(
+                        Vmf*byDt(MRF.absolute(phase.phi().oldTime()))
+                      + Vmf*ddtPhis[pair.other(phase).index()]
+                    );
             }
         }