From 0b81cd8f6fc78c9fa5af1004a2405e1d92871f30 Mon Sep 17 00:00:00 2001
From: Markus Blatt <markus@dr-blatt.de>
Date: Mon, 19 Jan 2015 20:22:34 +0100
Subject: [PATCH] Reduces code duplication and reorganizes code when checking
 convergence.

The computations made to check the convergence are the same for all existing
phases. Therefore this patch uses loops over the phase indices when cmputing
them,

In the convergence check there are several reductions (maxCoeff(), sum())
that will trigger communication in a parallel run. This patch seperates the
reductions from the other computations. The idea is to one reduction for the
reductions that need to done as global communication is expensive.
---
 .../FullyImplicitBlackoilSolver_impl.hpp      | 121 ++++++++----------
 1 file changed, 51 insertions(+), 70 deletions(-)

diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
index 4c8fed7b8..ca4e1ce5b 100644
--- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
+++ b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
@@ -1900,79 +1900,60 @@ namespace {
 
         const std::vector<PhasePresence> cond = phaseCondition();
 
-        double CNVW = 0.;
-        double CNVO = 0.;
-        double CNVG = 0.;
+        double CNV[MaxNumPhases] = {0., 0., 0.};
+        double R_sum[MaxNumPhases] = {0., 0., 0.};
+        double B_avg[MaxNumPhases] = {0., 0., 0.};
+        double maxCoeff[MaxNumPhases] = {0., 0., 0.};
+        double mass_balance_residual[MaxNumPhases] = {0., 0., 0.};
+        Eigen::Array<V::Scalar, Eigen::Dynamic, MaxNumPhases> B;
+        Eigen::Array<V::Scalar, Eigen::Dynamic, MaxNumPhases> R;
+        Eigen::Array<V::Scalar, Eigen::Dynamic, MaxNumPhases> tempV;
 
-        double RW_sum = 0.;
-        double RO_sum = 0.;
-        double RG_sum = 0.;
-
-        double BW_avg = 0.;
-        double BO_avg = 0.;
-        double BG_avg = 0.;
-
-        if (active_[Water]) {
-            const int pos = pu.phase_pos[Water];
-            const ADB& tempBW = rq_[pos].b;
-            V BW = 1./tempBW.value();
-            V RW = residual_.material_balance_eq[Water].value();
-            BW_avg = BW.sum()/nc;
-            const V tempV = RW.abs()/pv;
-
-            CNVW = BW_avg * dt * tempV.maxCoeff();
-            RW_sum = RW.sum();
+        for(int idx=0; idx<MaxNumPhases; ++idx)
+        {
+            if (active_[idx]) {
+                const int pos = pu.phase_pos[idx];
+                const ADB& tempB = rq_[pos].b;
+                B.col(pu.phase_pos[idx])  = 1./tempB.value();
+                R.col(pu.phase_pos[idx]) = residual_.material_balance_eq[pu.phase_pos[idx]].value();
+                tempV.col(pu.phase_pos[idx]) = R.col(pu.phase_pos[idx]).abs()/pv;
+            }
+        }
+        // Do the global reductions
+        for(int idx=0; idx<MaxNumPhases; ++idx)
+        {
+            if (active_[idx]) {
+                B_avg[pu.phase_pos[idx]] = B.sum()/nc;
+                maxCoeff[pu.phase_pos[idx]]=tempV.col(pu.phase_pos[idx]).maxCoeff();
+                R_sum[pu.phase_pos[idx]] = R.col(pu.phase_pos[idx]).sum();
+            }else{
+                R_sum[pu.phase_pos[idx]] = B_avg[pu.phase_pos[idx]] = 0.;
+            }
         }
 
-        if (active_[Oil]) {
-            // Omit the disgas here. We should add it.
-            const int pos = pu.phase_pos[Oil];
-            const ADB& tempBO = rq_[pos].b;
-            V BO = 1./tempBO.value();
-            V RO = residual_.material_balance_eq[Oil].value();
-            BO_avg = BO.sum()/nc;
-            const V tempV = RO.abs()/pv;
-
-            CNVO = BO_avg * dt * tempV.maxCoeff();
-            RO_sum = RO.sum();
+        bool converged_MB = true;
+        bool converged_CNV = true;
+        // Finish computation
+        for(int idx=0; idx<MaxNumPhases; ++idx)
+        {
+            CNV[idx] = B_avg[pu.phase_pos[idx]] * dt * maxCoeff[pu.phase_pos[idx]];
+            mass_balance_residual[idx] = std::abs(B_avg[pu.phase_pos[idx]]*R_sum[pu.phase_pos[idx]]) * dt / pvSum;
+            converged_MB = converged_MB && (mass_balance_residual[idx] < tol_mb);
+            converged_CNV = converged_CNV && (CNV[idx] < tol_cnv);
         }
 
-        if (active_[Gas]) {
-            // Omit the vapoil here. We should add it.
-            const int pos = pu.phase_pos[Gas];
-            const ADB& tempBG = rq_[pos].b;
-            V BG = 1./tempBG.value();
-            V RG = residual_.material_balance_eq[Gas].value();
-            BG_avg = BG.sum()/nc;
-            const V tempV = RG.abs()/pv;
-
-            CNVG = BG_avg * dt * tempV.maxCoeff();
-            RG_sum = RG.sum();
-        }
-
-        const double mass_balance_residual_water = std::abs(BW_avg*RW_sum) * dt / pvSum;
-        const double mass_balance_residual_oil   = std::abs(BO_avg*RO_sum) * dt / pvSum;
-        const double mass_balance_residual_gas   = std::abs(BG_avg*RG_sum) * dt / pvSum;
-
-        bool converged_MB = (mass_balance_residual_water < tol_mb)
-                         && (mass_balance_residual_oil   < tol_mb)
-                         && (mass_balance_residual_gas   < tol_mb);
-
-        bool converged_CNV = (CNVW < tol_cnv) && (CNVO < tol_cnv) && (CNVG < tol_cnv);
-
-        const double residualWellFlux = infinityNorm( residual_.well_flux_eq );
-        const double residualWell     = infinityNorm( residual_.well_eq );
-
+        double residualWellFlux = infinityNorm(residual_.well_flux_eq);
+        double residualWell = infinityNorm(residual_.well_eq);
         bool converged_Well = (residualWellFlux < 1./Opm::unit::day) && (residualWell < Opm::unit::barsa);
         bool converged = converged_MB && converged_CNV && converged_Well;
 
         // if one of the residuals is NaN, throw exception, so that the solver can be restarted
-        if( std::isnan(mass_balance_residual_water) || mass_balance_residual_water > maxResidualAllowed() ||
-            std::isnan(mass_balance_residual_oil)   || mass_balance_residual_oil   > maxResidualAllowed() ||
-            std::isnan(mass_balance_residual_gas)   || mass_balance_residual_gas   > maxResidualAllowed() ||
-            std::isnan(CNVW) || CNVW > maxResidualAllowed() ||
-            std::isnan(CNVO) || CNVO > maxResidualAllowed() ||
-            std::isnan(CNVG) || CNVG > maxResidualAllowed() ||
+        if( std::isnan(mass_balance_residual[Water]) || mass_balance_residual[Water] > maxResidualAllowed() ||
+            std::isnan(mass_balance_residual[Oil])   || mass_balance_residual[Oil]   > maxResidualAllowed() ||
+            std::isnan(mass_balance_residual[Gas])   || mass_balance_residual[Gas]   > maxResidualAllowed() ||
+            std::isnan(CNV[Water]) || CNV[Water] > maxResidualAllowed() ||
+            std::isnan(CNV[Oil]) || CNV[Oil] > maxResidualAllowed() ||
+            std::isnan(CNV[Gas]) || CNV[Gas] > maxResidualAllowed() ||
             std::isnan(residualWellFlux) || residualWellFlux > maxResidualAllowed() ||
             std::isnan(residualWell)     || residualWell     > maxResidualAllowed() )
         {
@@ -1985,12 +1966,12 @@ namespace {
         const std::streamsize oprec = std::cout.precision(3);
         const std::ios::fmtflags oflags = std::cout.setf(std::ios::scientific);
         std::cout << std::setw(4) << iteration
-                  << std::setw(11) << mass_balance_residual_water
-                  << std::setw(11) << mass_balance_residual_oil
-                  << std::setw(11) << mass_balance_residual_gas
-                  << std::setw(11) << CNVW
-                  << std::setw(11) << CNVO
-                  << std::setw(11) << CNVG
+                  << std::setw(11) << mass_balance_residual[Water]
+                  << std::setw(11) << mass_balance_residual[Oil]
+                  << std::setw(11) << mass_balance_residual[Gas]
+                  << std::setw(11) << CNV[Water]
+                  << std::setw(11) << CNV[Oil]
+                  << std::setw(11) << CNV[Gas]
                   << std::setw(11) << residualWellFlux
                   << std::setw(11) << residualWell
                   << std::endl;