From c685e86f915c66d62e9790ba28e9dd8fa5b06c47 Mon Sep 17 00:00:00 2001
From: Andreas Lauser <and@poware.org>
Date: Sun, 19 Nov 2017 13:03:26 +0100
Subject: [PATCH] add an energy "phase"

This is quite a hack: Even though energy is not a "phase" and it is
also not considered in MaxNumPhases and pu.num_phases because this
would break a lot of assumptions in old code, it is nevertheless
assigned an "canonical index" that can be translated "active index"
via PhaseUsage::phase_pos[]. This awkwardness is needed because much
of the legacy OPM code conflates the concepts of "fluid phase" and
"conserved quantity" and fixing that issue would basically mean an
almost complete rewrite of much of the legacy code. That said, the
same statement applies to polymer and solvent, but these are currently
handled as even more second-class citizens because they are not even
given a canonical index and also cannot be translated into an active
one.
---
 opm/core/props/BlackoilPhases.hpp     |  7 ++++---
 opm/core/props/phaseUsageFromDeck.hpp | 27 ++++++++++++++++++++++++---
 2 files changed, 28 insertions(+), 6 deletions(-)

diff --git a/opm/core/props/BlackoilPhases.hpp b/opm/core/props/BlackoilPhases.hpp
index 5bf979ea..0cc025bc 100644
--- a/opm/core/props/BlackoilPhases.hpp
+++ b/opm/core/props/BlackoilPhases.hpp
@@ -29,17 +29,18 @@ namespace Opm
     public:
         static const int MaxNumPhases = 3;
         // enum ComponentIndex { Water = 0, Oil = 1, Gas = 2 };
-        enum PhaseIndex { Aqua = 0, Liquid = 1, Vapour = 2 };
+        enum PhaseIndex { Aqua = 0, Liquid = 1, Vapour = 2, Energy = 3 };
 
     };
 
     struct PhaseUsage : public BlackoilPhases
     {
         int num_phases;
-        int phase_used[MaxNumPhases];
-        int phase_pos[MaxNumPhases];
+        int phase_used[MaxNumPhases + 1];
+        int phase_pos[MaxNumPhases + 1];
         bool has_solvent;
         bool has_polymer;
+        bool has_energy;
     };
 
     /// Check or assign presence of a formed, free phase.  Limited to
diff --git a/opm/core/props/phaseUsageFromDeck.hpp b/opm/core/props/phaseUsageFromDeck.hpp
index 65d9572c..ded3edb5 100644
--- a/opm/core/props/phaseUsageFromDeck.hpp
+++ b/opm/core/props/phaseUsageFromDeck.hpp
@@ -50,9 +50,16 @@ namespace Opm
             pu.phase_used[BlackoilPhases::Vapour] = 1;
         }
         pu.num_phases = 0;
-        for (int i = 0; i < BlackoilPhases::MaxNumPhases; ++i) {
-            pu.phase_pos[i] = pu.num_phases;
-            pu.num_phases += pu.phase_used[i];
+        int numActivePhases = 0;
+        for (int phaseIdx = 0; phaseIdx < BlackoilPhases::MaxNumPhases; ++phaseIdx) {
+            if (!pu.phase_used[numActivePhases]) {
+                pu.phase_pos[phaseIdx] = -1;
+            }
+            else {
+                pu.phase_pos[phaseIdx] = numActivePhases;
+                ++ numActivePhases;
+                pu.num_phases = numActivePhases;
+            }
         }
 
         // Only 2 or 3 phase systems handled.
@@ -78,6 +85,17 @@ namespace Opm
             pu.has_polymer = true;
         }
 
+        // Add energy info
+        pu.has_energy = phase.active(Phase::ENERGY);
+        if (pu.has_energy) {
+            // this is quite a hack: even though energy is not considered as in
+            // MaxNumPhases and pu.num_phases because this would break a lot of
+            // assumptions in old code, it is nevertheless an index to be translated
+            // to. polymer and solvent are even larger hacks because not even this can be
+            // done for them.
+            pu.phase_pos[BlackoilPhases::Energy] = numActivePhases;
+            ++ numActivePhases;
+        }
         return pu;
     }
 
@@ -136,6 +154,9 @@ namespace Opm
             pu.has_polymer = true;
         }
 
+        // Add energy info
+        pu.has_energy = phase.active(Phase::ENERGY);
+
         return pu;
     }