From 72d0b8123bf4fb3b67b74e98a85767ccacdf674e Mon Sep 17 00:00:00 2001
From: Markus Blatt <markus@dr-blatt.de>
Date: Wed, 8 Jun 2016 11:07:34 +0200
Subject: [PATCH] Fixes the parallel Euclidian product for multiple phases.

This is used to compute the Euclidian product for the saturations.
Thes are ordered in an interleaved manner (all saturations for cell
with index 0, the all for index 1, ...). Up to now the implementation
assumed a different ordering: blockwise (all saturations for phase 0 first,
then all saturations phase 1, ...).
With this commit the computation uses the right assumption.
---
 opm/autodiff/BlackoilModelBase_impl.hpp | 38 ++++++++++++++++++++-----
 1 file changed, 31 insertions(+), 7 deletions(-)

diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp
index 3f2b0e5a0..434d9f2ca 100644
--- a/opm/autodiff/BlackoilModelBase_impl.hpp
+++ b/opm/autodiff/BlackoilModelBase_impl.hpp
@@ -1121,6 +1121,10 @@ namespace detail {
         }
 
         /// \brief Compute the Euclidian norm of a vector
+        /// \warning In the case that num_components is greater than 1
+        ///          an interleaved ordering is assumed. E.g. for each cell
+        ///          all phases of that cell are stored consecutively. First
+        ///          the ones for cell 0, then the ones for cell 1, ... .
         /// \param it              begin iterator for the given vector
         /// \param end             end iterator for the given vector
         /// \param num_components  number of components (i.e. phases) in the vector
@@ -1136,19 +1140,39 @@ namespace detail {
             {
                 const ParallelISTLInformation& info =
                     boost::any_cast<const ParallelISTLInformation&>(pinfo);
-                int size_per_component = (end - it) / num_components;
+                typedef typename Iterator::value_type Scalar;
+                Scalar product = 0.0;
+                int size_per_component = (end - it);
+                size_per_component /= num_components; // two lines to supresse unused warning.
                 assert((end - it) == num_components * size_per_component);
-                double component_product = 0.0;
-                for( int i = 0; i < num_components; ++i )
+
+                if( num_components == 1 )
                 {
                     auto component_container =
-                        boost::make_iterator_range(it + i * size_per_component,
-                                                   it + (i + 1) * size_per_component);
+                        boost::make_iterator_range(it, end);
                     info.computeReduction(component_container,
                                            Opm::Reduction::makeInnerProductFunctor<double>(),
-                                           component_product);
+                                           product);
                 }
-                return component_product;
+                else
+                {
+                    auto& maskContainer = info.getOwnerMask();
+                    auto mask = maskContainer.begin();
+                    assert(maskContainer.size() == size_per_component);
+
+                    for(int cell = 0; cell < size_per_component; ++cell, ++mask)
+                    {
+                        Scalar cell_product = (*it) * (*it);
+                        ++it;
+                        for(int component=1; component < num_components;
+                            ++component, ++it)
+                        {
+                            cell_product += (*it) * (*it);
+                        }
+                        product += cell_product * (*mask);
+                    }
+                }
+                return product;
             }
             else
 #endif