From a9cb40b55b3aeec9cb60feeec02d45ae44b132c3 Mon Sep 17 00:00:00 2001
From: Henry Weller <http://openfoam.org>
Date: Mon, 1 Jan 2018 16:02:57 +0000
Subject: [PATCH] reactingEulerFoam::phasePair: Provide more convenient method
 to "loop" over pair

Checking a pair contains a particular phase and adding a contribution from the
"other" phase can now be written:

            if (pair.contains(phase))
            {
                const phaseModel& otherPhase = pair.other(phase);

                phiHbyAs[phasei] +=
                    fvc::interpolate(rAUs[phasei]*K)
                   *MRF.absolute(otherPhase.phi());

                HbyAs[phasei] += rAUs[phasei]*K*otherPhase.U();
            }

which previously would have been written as a loop over the pair and excluding
self reference:

            const phaseModel* phase1 = &pair.phase1();
            const phaseModel* phase2 = &pair.phase2();

            forAllConstIter(phasePair, pair, iter)
            {
                if (phase1 == &phase)
                {
                    phiHbyAs[phasei] +=
                        fvc::interpolate(rAUs[phasei]*K)
                       *MRF.absolute(phase2->phi());

                    HbyAs[phasei] += rAUs[phasei]*K*phase2->U();
                }

                Swap(phase1, phase2);
            }
---
 .../phasePair/phasePair/phasePair.H           | 16 +++++++---
 .../phasePair/phasePair/phasePairI.H          | 32 ++++++++++++++++++-
 .../reactingMultiphaseEulerFoam/pU/pEqn.H     | 18 ++++-------
 3 files changed, 49 insertions(+), 17 deletions(-)

diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phasePair/phasePair/phasePair.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phasePair/phasePair/phasePair.H
index 4ee9eefbc..26ef837e8 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phasePair/phasePair/phasePair.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phasePair/phasePair/phasePair.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2014-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -150,13 +150,21 @@ public:
 
         // Access
 
-            // Phase 1
+            //- Return phase 1
             inline const phaseModel& phase1() const;
 
-            // Phase 2
+            //- Return phase 2
             inline const phaseModel& phase2() const;
 
-            // Gravitation acceleration
+            //- Return true if this phasePair contains the given phase
+            inline bool contains(const phaseModel& phase) const;
+
+            //- Return the other phase relative to the given phase
+            //  Generates a FatalError if this phasePair does not contain
+            //  the given phase
+            inline const phaseModel& other(const phaseModel& phase) const;
+
+            //- Return gravitation acceleration
             inline const uniformDimensionedVectorField& g() const;
 };
 
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phasePair/phasePair/phasePairI.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phasePair/phasePair/phasePairI.H
index c461b7d91..6eaf1a612 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phasePair/phasePair/phasePairI.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phasePair/phasePair/phasePairI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2014-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,6 +37,36 @@ inline const Foam::phaseModel& Foam::phasePair::phase2() const
 }
 
 
+inline bool Foam::phasePair::contains(const phaseModel& phase) const
+{
+    return &phase1_ == &phase || & phase2_ == &phase;
+}
+
+
+inline const Foam::phaseModel& Foam::phasePair::other
+(
+    const phaseModel& phase
+) const
+{
+    if (&phase1_ == &phase)
+    {
+        return phase2_;
+    }
+    else if (&phase2_ == &phase)
+    {
+        return phase1_;
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "this phasePair does not contain phase " << phase.name()
+            << exit(FatalError);
+
+        return phase;
+    }
+}
+
+
 inline const Foam::uniformDimensionedVectorField& Foam::phasePair::g() const
 {
     return g_;
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H
index 696385d7b..0fdad237b 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H
@@ -227,21 +227,15 @@ while (pimple.correct())
 
             const phasePair& pair(fluid.phasePairs()[KdIter.key()]);
 
-            const phaseModel* phase1 = &pair.phase1();
-            const phaseModel* phase2 = &pair.phase2();
-
-            forAllConstIter(phasePair, pair, iter)
+            if (pair.contains(phase))
             {
-                if (phase1 == &phase)
-                {
-                    phiHbyAs[phasei] +=
-                        fvc::interpolate(rAUs[phasei]*K)
-                       *MRF.absolute(phase2->phi());
+                const phaseModel& otherPhase = pair.other(phase);
 
-                    HbyAs[phasei] += rAUs[phasei]*K*phase2->U();
-                }
+                phiHbyAs[phasei] +=
+                    fvc::interpolate(rAUs[phasei]*K)
+                   *MRF.absolute(otherPhase.phi());
 
-                Swap(phase1, phase2);
+                HbyAs[phasei] += rAUs[phasei]*K*otherPhase.U();
             }
         }