Merge branch 'master' into suppress-warnings

Conflicts:
	cmake/Modules/UseWarnings.cmake

CMakeLists_files.cmake

@@ -28,6 +28,7 @@
 list (APPEND MAIN_SOURCE_FILES
 	opm/autodiff/BlackoilPropsAd.cpp
 	opm/autodiff/BlackoilPropsAdInterface.cpp
+	opm/autodiff/NewtonIterationBlackoilSimple.cpp
 	opm/autodiff/GridHelpers.cpp
 	opm/autodiff/ImpesTPFAAD.cpp
 	opm/autodiff/SimulatorCompressibleAd.cpp
@@ -103,6 +104,9 @@ list (APPEND PUBLIC_HEADER_FILES
 	opm/autodiff/ImpesTPFAAD.hpp
 	opm/autodiff/FullyImplicitBlackoilSolver.hpp
 	opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
+	opm/autodiff/NewtonIterationBlackoilInterface.hpp
+	opm/autodiff/NewtonIterationBlackoilSimple.hpp
+	opm/autodiff/LinearisedBlackoilResidual.hpp
 	opm/autodiff/SimulatorCompressibleAd.hpp
 	opm/autodiff/SimulatorFullyImplicitBlackoil.hpp
 	opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp

cmake/Modules/FindSuiteSparse.cmake

@@ -232,14 +232,6 @@ if (UMFPACK_LIBRARY)
 	  set (UMFPACK_EXTRA_LIBS "-NOTFOUND")
 	endif (CHOLMOD_LIBRARIES)
   endif (HAVE_UMFPACK_WITHOUT_CHOLMOD)
-  # test if umfpack is underlinked (CentOS 5.9), i.e. doesn't specify
-  # that it depends on amd. in that case, force amd to be linked
-  if (UMFPACK_EXTRA_LIBS AND (CMAKE_CXX_PLATFORM_ID STREQUAL "Linux") AND CMAKE_COMPILER_IS_GNUCC)
-	try_compile_umfpack (HAVE_UMFPACK_NOT_UNDERLINKED "-Wl,--as-needed" ${UMFPACK_EXTRA_LIBS})
-	if (NOT HAVE_UMFPACK_NOT_UNDERLINKED)
-	  list (APPEND UMFPACK_LINKER_FLAGS "-Wl,--no-as-needed")
-	endif (NOT HAVE_UMFPACK_NOT_UNDERLINKED)
-  endif (UMFPACK_EXTRA_LIBS AND (CMAKE_CXX_PLATFORM_ID STREQUAL "Linux") AND CMAKE_COMPILER_IS_GNUCC)
   list (APPEND UMFPACK_LIBRARIES ${UMFPACK_EXTRA_LIBS})
   list (REVERSE UMFPACK_LIBRARIES)
   list (REMOVE_DUPLICATES UMFPACK_LIBRARIES)

cmake/Modules/FindUG.cmake

@@ -0,0 +1,112 @@
+#
+# This module first tests for UG and then sets the necessary flags
+# and config.h defines. If UG is found UG_FOUND will be true.
+#
+
+# this function is required in order not to pollute the global
+# namespace with the macros defined in ug-config*.cmake
+function(opmFindUg)
+  if(NOT UG_ROOT)
+    # check whether UG is in /usr/local
+    if(EXISTS "/usr/local/include/ug")
+      set(UG_ROOT "/usr/local")
+
+      # check whether UG is in /usr
+    elseif(EXISTS "/usr/include/ug")
+      set(UG_ROOT "/usr")
+      
+      # oops
+    else()
+      message(STATUS "Could not find UG. It seems to be not installed.")
+      return()
+    endif()
+  endif()
+
+  if(UG_ROOT AND NOT UG_DIR)
+    # define the directory where the config file resides
+    if(EXISTS "${UG_ROOT}/lib/cmake/ug/ug-config.cmake")
+      set(UG_DIR ${UG_ROOT}/lib/cmake/ug)
+    elseif(EXISTS "${UG_ROOT}/lib64/cmake/ug/ug-config.cmake")
+      set(UG_DIR ${UG_ROOT}/lib64/cmake/ug)
+    else()
+      message(WARNING "Could not find file ug-config.cmake relative to given UG_ROOT")
+      return()
+    endif()
+  endif()
+
+  # include the config mode files kindly provided by UG...
+  include(${UG_DIR}/ug-config-version.cmake)
+  include(${UG_DIR}/ug-config.cmake)
+
+  set(UG_FOUND "1")
+  if(NOT UG_FOR_DUNE STREQUAL "yes")
+    set(UG_FOUND "0")
+    message(WARNING "UG was not configured for DUNE. Did pass --enable-dune to its configure?")
+    return()
+  endif()
+
+  set(HAVE_UG ${UG_FOUND})
+
+  # parse version
+  string(REGEX REPLACE "([0-9]*)\\.[0-9]*\\..*" "\\1" UG_VERSION_MAJOR "${PACKAGE_VERSION}")
+  string(REGEX REPLACE "[0-9]*\\.([0-9]*)\\..*" "\\1" UG_VERSION_MINOR "${PACKAGE_VERSION}")
+  string(REGEX REPLACE "[0-9]*\\.[0-9]*\\.([0-9]*).*" "\\1" UG_VERSION_REVISION "${PACKAGE_VERSION}")
+
+  string(REGEX REPLACE ".*-patch([0-9]*)" "\\1" TMP "${PACKAGE_VERSION}")
+  if(TMP STREQUAL "${PACKAGE_VERSION}")
+    set(UG_VERSION_PATCHLEVEL "")
+  else()
+    set(UG_VERSION_PATCHLEVEL "${TMP}")
+  endif()
+
+  # Adjust compiler/linker arguments
+  set(UG_LIBRARY_DIR "${libdir}")
+  
+  foreach (UG_RAW_LIB "-lugS2" "-lugS3" "-ldevS")
+    string(REGEX REPLACE "-l(.*)" "\\1" UG_LIB "${UG_RAW_LIB}")
+    set(UG_LIB_FILE  "${UG_LIBRARY_DIR}/lib${UG_LIB}.a")
+    if (EXISTS "${UG_LIB_FILE}")
+      set(UG_LIBS "${UG_LIBS}" ${UG_LIB_FILE})
+    else()
+      set(UG_LIBS "${UG_LIBS}" ${UG_LIB})
+    endif()
+  endforeach()
+
+  set(UG_LIBRARIES "${UG_LIBS}")
+
+  # export all variables which need to be seen globally
+  set(UG_FOUND "${UG_FOUND}" PARENT_SCOPE)
+  set(HAVE_UG "${HAVE_UG}" PARENT_SCOPE)
+  set(UG_INCLUDE_DIRS "${UG_INCLUDES}" PARENT_SCOPE)
+  set(UG_LIBRARIES "${UG_LIBRARIES}" PARENT_SCOPE)
+  set(UG_VERSION_MAJOR "${UG_VERSION_MAJOR}" PARENT_SCOPE)
+  set(UG_VERSION_MINOR "${UG_VERSION_MINOR}" PARENT_SCOPE)
+  set(UG_VERSION_REVISION "${UG_VERSION_REVISION}" PARENT_SCOPE)
+  set(UG_VERSION_PATCHLEVEL "${UG_VERSION_PATCHLEVEL}" PARENT_SCOPE)
+
+  set(UG_DEFINITIONS "${UG_COMPILE_FLAGS}" PARENT_SCOPE)
+endfunction()
+
+if (NOT HAVE_UG)
+  opmFindUg()
+  
+  set(HAVE_UG "${HAVE_UG}" CACHE BOOL "UG library is available")
+  set(UG_INCLUDE_DIRS "${UG_INCLUDE_DIRS}" CACHE STRING "Directory containing the headers of the UG library")
+  set(UG_LIBRARIES "${UG_LIBRARIES}" CACHE STRING "The libraries which need to be linked to be able to use the UG library")
+  set(UG_DEFINITIONS "${UG_DEFINITIONS}" CACHE STRING "The compiler flags for the UG library")
+  set(UG_VERSION_MAJOR "${UG_VERSION_MAJOR}" CACHE INT "Major version of the UG release")
+  set(UG_VERSION_MINOR "${UG_VERSION_MINOR}" CACHE INT "Minor version of the UG release")
+  set(UG_VERSION_REVISION "${UG_VERSION_REVISION}" CACHE INT "Revision of the UG release")
+  set(UG_VERSION_PATCHLEVEL "${UG_VERSION_PATCHLEVEL}" CACHE INT "Patchlevel of the UG release")
+
+  mark_as_advanced(HAVE_UG)
+  mark_as_advanced(UG_INCLUDE_DIRS)
+  mark_as_advanced(UG_LIBRARIES)
+  mark_as_advanced(UG_DEFINITIONS)
+  mark_as_advanced(UG_VERSION_MAJOR)
+  mark_as_advanced(UG_VERSION_MINOR)
+  mark_as_advanced(UG_VERSION_REVISION)
+  mark_as_advanced(UG_VERSION_PATCHLEVEL)
+else()
+  set(UG_FOUND "0")
+endif()

cmake/Modules/Finddune-grid.cmake

@@ -23,7 +23,8 @@ find_opm_package (
   dune-common REQUIRED;
   dune-geometry REQUIRED;
   MPI;
-  ALUGrid
+  ALUGrid;
+  UG
   "
   # header to search for
   "dune/grid/onedgrid.hh"

cmake/Modules/UseOnlyNeeded.cmake

@@ -23,9 +23,11 @@ function (prepend var_name value)
   endif (NOT ("${_var_pre}" STREQUAL "${value}"))
 endfunction (prepend var_name value)
 
+option (ONLY_NEEDED_LIBRARIES "Instruct the linker to not use libraries which are unused" OFF)
+
 # only ELF shared objects can be underlinked, and only GNU will accept
 # these parameters; otherwise just leave it to the defaults
-if ((CMAKE_CXX_PLATFORM_ID STREQUAL "Linux") AND CMAKE_COMPILER_IS_GNUCC)
+if ((CMAKE_CXX_PLATFORM_ID STREQUAL "Linux") AND CMAKE_COMPILER_IS_GNUCC AND ONLY_NEEDED_LIBRARIES)
   # these are the modules whose probes will turn up incompatible
   # flags on some systems
   set (_maybe_underlinked
@@ -45,4 +47,4 @@ if ((CMAKE_CXX_PLATFORM_ID STREQUAL "Linux") AND CMAKE_COMPILER_IS_GNUCC)
 	prepend (CMAKE_MODULE_LINKER_FLAGS "-Wl,--as-needed")
 	prepend (CMAKE_SHARED_LINKER_FLAGS "-Wl,--as-needed")
   endif (NOT _underlinked)
-endif ((CMAKE_CXX_PLATFORM_ID STREQUAL "Linux")  AND CMAKE_COMPILER_IS_GNUCC)
+endif ((CMAKE_CXX_PLATFORM_ID STREQUAL "Linux")  AND CMAKE_COMPILER_IS_GNUCC AND ONLY_NEEDED_LIBRARIES)

cmake/Modules/UseWarnings.cmake

@@ -13,7 +13,7 @@ if (CXX_COMPAT_GCC)
   endif (_warn_flag)
 endif ()
 
-option(SILENCE_DUNE_WARNINGS "Disable warnings from DUNE?" OFF)
+option(SILENCE_EXTERNAL_WARNINGS "Disable some warnings from external packages (requires GCC 4.6 or newer)" OFF)
 if(SILENCE_EXTERNAL_WARNINGS AND CXX_COMPAT_GCC)
   file(WRITE ${CMAKE_BINARY_DIR}/disable_warning_pragmas.h "
 #pragma GCC diagnostic push

cmake/Modules/opm-core-prereqs.cmake

@@ -5,6 +5,8 @@
 set (opm-core_CONFIG_VAR
 	HAVE_ERT
 	HAVE_SUITESPARSE_UMFPACK_H
+	HAVE_DUNE_ISTL
+	HAVE_MPI
 	)
 
 # dependencies
@@ -27,6 +29,8 @@ set (opm-core_DEPS
 	"TinyXML"
 	# Ensembles-based Reservoir Tools (ERT)
 	"ERT"
+	# Look for MPI support
+	"MPI"
 	# DUNE dependency
 	"dune-common"
 	"dune-istl"

cmake/Scripts/configure

@@ -36,6 +36,7 @@ Optional Features:
   --config-cache          Reuse build configuration cache from a previous run
 
 Optional Packages:
+  --with-ug=PATH          use the UG libraries from a specified location
   --with-alugrid=PATH     use the ALUGrid library from a specified location
   --with-metis=PATH       use the METIS graph partitioning library from a specified location
   --with-boost=PATH       use Boost library from a specified location
@@ -256,6 +257,7 @@ for OPT in "$@"; do
             superlu     |\
             SuiteSparse |\
             TinyXML     |\
+            ug          |\
             opm         |\
             opm-*       |\
             dune        |\

examples/sim_fibo_ad.cpp

@@ -37,6 +37,7 @@
 #include <opm/core/props/rock/RockCompressibility.hpp>
 
 #include <opm/core/linalg/LinearSolverFactory.hpp>
+#include <opm/autodiff/NewtonIterationBlackoilSimple.hpp>
 
 #include <opm/core/simulator/BlackoilState.hpp>
 #include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
@@ -141,6 +142,7 @@ try
 
     // Linear solver.
     LinearSolverFactory linsolver(param);
+    NewtonIterationBlackoilSimple fis_solver(linsolver);
 
     // Write parameters used for later reference.
     bool output = param.getDefault("output", true);
@@ -210,7 +212,7 @@ try
                                                  *new_props,
                                                  rock_comp->isActive() ? rock_comp.get() : 0,
                                                  wells,
-                                                 linsolver,
+                                                 fis_solver,
                                                  grav);
         SimulatorReport episodeReport = simulator.run(simtimer, state, well_state);
 

examples/sim_fibo_ad_cp.cpp

@@ -59,6 +59,7 @@
 #include <opm/core/props/rock/RockCompressibility.hpp>
 
 #include <opm/core/linalg/LinearSolverFactory.hpp>
+#include <opm/autodiff/NewtonIterationBlackoilSimple.hpp>
 
 #include <opm/core/simulator/BlackoilState.hpp>
 #include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
@@ -195,6 +196,7 @@ try
 
     // Linear solver.
     LinearSolverFactory linsolver(param);
+    NewtonIterationBlackoilSimple fis_solver(linsolver);
 
     // Write parameters used for later reference.
     bool output = param.getDefault("output", true);
@@ -269,7 +271,7 @@ try
                                                                *new_props,
                                                                rock_comp->isActive() ? rock_comp.get() : 0,
                                                                wells,
-                                                               linsolver,
+                                                               fis_solver,
                                                                grav);
         SimulatorReport episodeReport = simulator.run(simtimer, state, well_state);
 

examples/test_implicit_ad.cpp

@@ -24,6 +24,7 @@
 #include <opm/autodiff/GeoProps.hpp>
 #include <opm/autodiff/BlackoilPropsAd.hpp>
 #include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
+#include <opm/autodiff/NewtonIterationBlackoilSimple.hpp>
 
 #include <opm/core/grid.h>
 #include <opm/core/wells.h>
@@ -103,8 +104,9 @@ try
     Opm::DerivedGeology geo(*g, props, grav);
 
     Opm::LinearSolverFactory linsolver(param);
+    Opm::NewtonIterationBlackoilSimple fis_solver(linsolver);
 
-    Opm::FullyImplicitBlackoilSolver<UnstructuredGrid> solver(*g, props, geo, 0, *wells, linsolver);
+    Opm::FullyImplicitBlackoilSolver<UnstructuredGrid> solver(*g, props, geo, 0, *wells, fis_solver);
 
     Opm::BlackoilState state;
     initStateBasic(*g, props0, param, 0.0, state);

opm/autodiff/FullyImplicitBlackoilSolver.hpp

@@ -23,6 +23,8 @@
 #include <opm/autodiff/AutoDiffBlock.hpp>
 #include <opm/autodiff/AutoDiffHelpers.hpp>
 #include <opm/autodiff/BlackoilPropsAdInterface.hpp>
+#include <opm/autodiff/LinearisedBlackoilResidual.hpp>
+#include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
 
 struct UnstructuredGrid;
 struct Wells;
@@ -31,7 +33,7 @@ namespace Opm {
 
     class DerivedGeology;
     class RockCompressibility;
-    class LinearSolverInterface;
+    class NewtonIterationBlackoilInterface;
     class BlackoilState;
     class WellStateFullyImplicitBlackoil;
 
@@ -65,7 +67,7 @@ namespace Opm {
                                     const DerivedGeology&           geo  ,
                                     const RockCompressibility*      rock_comp_props,
                                     const Wells&                    wells,
-                                    const LinearSolverInterface&    linsolver);
+                                    const NewtonIterationBlackoilInterface& linsolver);
 
         /// Take a single forward step, modifiying
         ///   state.pressure()
@@ -126,7 +128,7 @@ namespace Opm {
         const DerivedGeology&           geo_;
         const RockCompressibility*      rock_comp_props_;
         const Wells&                    wells_;
-        const LinearSolverInterface&    linsolver_;
+        const NewtonIterationBlackoilInterface&    linsolver_;
         // For each canonical phase -> true if active
         const std::vector<bool>         active_;
         // Size = # active faces. Maps active -> canonical phase indices.
@@ -140,14 +142,7 @@ namespace Opm {
         std::vector<PhasePresence> phaseCondition_;
         V well_perforation_pressure_diffs_; // Diff to bhp for each well perforation.
 
-        // The mass_balance vector has one element for each active phase,
+        LinearisedBlackoilResidual residual_;
-        // each of which has size equal to the number of cells.
-        // The well_eq has size equal to the number of wells.
-        struct {
-            std::vector<ADB> mass_balance;
-            ADB well_flux_eq;
-            ADB well_eq;
-        } residual_;
 
         // Private methods.
         SolutionState

opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp

@@ -203,7 +203,7 @@ namespace {
                                 const DerivedGeology&           geo  ,
                                 const RockCompressibility*      rock_comp_props,
                                 const Wells&                    wells,
-                                const LinearSolverInterface&    linsolver)
+                                const NewtonIterationBlackoilInterface&    linsolver)
         : grid_  (grid)
         , fluid_ (fluid)
         , geo_   (geo)
@@ -662,13 +662,13 @@ namespace {
             // std::cout << "===== rq_[" << phase << "].mflux = \n" << std::endl;
             // std::cout << rq_[phase].mflux;
 
-            residual_.mass_balance[ phaseIdx ] =
+            residual_.material_balance_eq[ phaseIdx ] =
                 pvdt*(rq_[phaseIdx].accum[1] - rq_[phaseIdx].accum[0])
                 + ops_.div*rq_[phaseIdx].mflux;
 
 
             // DUMP(ops_.div*rq_[phase].mflux);
-            // DUMP(residual_.mass_balance[phase]);
+            // DUMP(residual_.material_balance_eq[phase]);
         }
 
         // -------- Extra (optional) rs and rv contributions to the mass balance equations --------
@@ -682,16 +682,16 @@ namespace {
                                                 rq_[po].head.value());
             const ADB rs_face = upwindOil.select(state.rs);
 
-            residual_.mass_balance[ Gas ] += ops_.div * (rs_face * rq_[po].mflux);
+            residual_.material_balance_eq[ Gas ] += ops_.div * (rs_face * rq_[po].mflux);
 
             const int pg = fluid_.phaseUsage().phase_pos[ Gas ];
             const UpwindSelector<double> upwindGas(grid_, ops_,
                                                 rq_[pg].head.value());
             const ADB rv_face = upwindGas.select(state.rv);
 
-            residual_.mass_balance[ Oil ] += ops_.div * (rv_face * rq_[pg].mflux);
+            residual_.material_balance_eq[ Oil ] += ops_.div * (rv_face * rq_[pg].mflux);
 
-            // DUMP(residual_.mass_balance[ Gas ]);
+            // DUMP(residual_.material_balance_eq[ Gas ]);
 
         }
 
@@ -870,7 +870,7 @@ namespace {
 
         // Add well contributions to mass balance equations
         for (int phase = 0; phase < np; ++phase) {
-            residual_.mass_balance[phase] -= superset(cq_s[phase],well_cells,nc);
+            residual_.material_balance_eq[phase] -= superset(cq_s[phase],well_cells,nc);
         }
 
         // Add WELL EQUATIONS
@@ -1154,7 +1154,7 @@ namespace {
             // const ADB well_contrib = superset(perf_flux*perf_b, well_cells, nc);
             well_contribs[phase] = superset(perf_flux*perf_b, well_cells, nc);
             // DUMP(well_contribs[phase]);
-            residual_.mass_balance[phase] += well_contribs[phase];
+            residual_.material_balance_eq[phase] += well_contribs[phase];
         }
         if (active_[Gas] && active_[Oil]) {
             const int oilpos = pu.phase_pos[Oil];
@@ -1164,8 +1164,8 @@ namespace {
             well_rates_all += superset(wops_.p2w * (well_perf_rates[oilpos]*rs_perf), Span(nw, 1, gaspos*nw), nw*np);
             well_rates_all += superset(wops_.p2w * (well_perf_rates[gaspos]*rv_perf), Span(nw, 1, oilpos*nw), nw*np);
             // DUMP(well_contribs[gaspos] + well_contribs[oilpos]*state.rs);
-            residual_.mass_balance[gaspos] += well_contribs[oilpos]*state.rs;
+            residual_.material_balance_eq[gaspos] += well_contribs[oilpos]*state.rs;
-            residual_.mass_balance[oilpos] += well_contribs[gaspos]*state.rv;
+            residual_.material_balance_eq[oilpos] += well_contribs[gaspos]*state.rv;
         }
 
         // Set the well flux equation
@@ -1180,40 +1180,7 @@ namespace {
     template<class T>
     V FullyImplicitBlackoilSolver<T>::solveJacobianSystem() const
     {
-        const int np = fluid_.numPhases();
+        return linsolver_.computeNewtonIncrement(residual_);
-        ADB mass_res = residual_.mass_balance[0];
-        for (int phase = 1; phase < np; ++phase) {
-            mass_res = vertcat(mass_res, residual_.mass_balance[phase]);
-        }
-        const ADB well_res = vertcat(residual_.well_flux_eq, residual_.well_eq);
-        const ADB total_residual = collapseJacs(vertcat(mass_res, well_res));
-        // DUMP(total_residual);
-
-        const Eigen::SparseMatrix<double, Eigen::RowMajor> matr = total_residual.derivative()[0];
-
-        V dx(V::Zero(total_residual.size()));
-        Opm::LinearSolverInterface::LinearSolverReport rep
-            = linsolver_.solve(matr.rows(), matr.nonZeros(),
-                               matr.outerIndexPtr(), matr.innerIndexPtr(), matr.valuePtr(),
-                               total_residual.value().data(), dx.data());
-        /*
-        std::ofstream filestream("matrix.out");
-        filestream << matr;
-        filestream.close();
-        std::ofstream filestream2("sol.out");
-        filestream2 << dx;
-        filestream2.close();
-        std::ofstream filestream3("r.out");
-        filestream3 << total_residual.value();
-        filestream3.close(); */
-
-
-        if (!rep.converged) {
-            OPM_THROW(std::runtime_error,
-                      "FullyImplicitBlackoilSolver::solveJacobianSystem(): "
-                      "Linear solver convergence failure.");
-        }
-        return dx;
     }
 
 
@@ -1630,10 +1597,9 @@ namespace {
     FullyImplicitBlackoilSolver<T>::residualNorm() const
     {
         double globalNorm = 0;
-        std::vector<ADB>::const_iterator quantityIt = residual_.mass_balance.begin();
+        std::vector<ADB>::const_iterator quantityIt = residual_.material_balance_eq.begin();
-        const std::vector<ADB>::const_iterator endQuantityIt = residual_.mass_balance.end();
+        const std::vector<ADB>::const_iterator endQuantityIt = residual_.material_balance_eq.end();
-        for (; quantityIt != endQuantityIt; ++quantityIt)
+        for (; quantityIt != endQuantityIt; ++quantityIt) {
-        {
             const double quantityResid = (*quantityIt).value().matrix().norm();
             if (!std::isfinite(quantityResid)) {
                 OPM_THROW(Opm::NumericalProblem,

opm/autodiff/LinearisedBlackoilResidual.hpp

@@ -0,0 +1,70 @@
+/*
+  Copyright 2014 SINTEF ICT, Applied Mathematics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_LINEARISEDBLACKOILRESIDUAL_HEADER_INCLUDED
+#define OPM_LINEARISEDBLACKOILRESIDUAL_HEADER_INCLUDED
+
+#include <opm/autodiff/AutoDiffBlock.hpp>
+
+namespace Opm
+{
+
+    /// Residual structure of the fully implicit solver.
+    /// All equations are given as AD types, with multiple
+    /// jacobian blocks corresponding to the primary unknowns. The
+    /// primary unknowns are for a three-phase simulation, in order:
+    ///    p    (pressure)
+    ///    sw   (water saturation)
+    ///    xvar (gas saturation, gas-oil ratio or oil-gas ratio)
+    ///    qs   (well outflows by well and phase)
+    ///    bhp  (bottom hole pressures)
+    /// In the above, the xvar variable will have a different
+    /// meaning from cell to cell, corresponding to the state in
+    /// that cell (saturated, undersaturated oil or undersaturated
+    /// gas). In a two-phase simulation, either sw or xvar is not
+    /// used, depending on which phase is missing.
+    ///
+    /// Note: this class is strongly coupled to the class
+    /// FullyImplicitBlackoilSolver, and is separated from that
+    /// class to facilitate the development of linear solver
+    /// strategies outside that class.
+    struct LinearisedBlackoilResidual {
+        /// A type alias for the automatic differentiation type.
+        typedef AutoDiffBlock<double> ADB;
+        /// The material_balance_eq vector has one element for each
+        /// active phase, each of which has size equal to the number
+        /// of cells. Each material balance equation is given in terms
+        /// of surface volumes (in SI units, that is standard m^3).
+        std::vector<ADB> material_balance_eq;
+        /// The well_flux_eq has size equal to the number of wells
+        /// times the number of phases. It contains the well flow
+        /// equations, relating the total well flows to
+        /// bottom-hole pressures and reservoir conditions.
+        ADB well_flux_eq;
+        /// The well_eq has size equal to the number of wells. It
+        /// contains the well control equations, that is for each
+        /// well either a rate specification or bottom hole
+        /// pressure specification.
+        ADB well_eq;
+    };
+
+} // namespace Opm
+
+
+#endif // OPM_LINEARISEDBLACKOILRESIDUAL_HEADER_INCLUDED

opm/autodiff/NewtonIterationBlackoilInterface.hpp

@@ -0,0 +1,46 @@
+/*
+  Copyright 2014 SINTEF ICT, Applied Mathematics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_NEWTONITERATIONBLACKOILINTERFACE_HEADER_INCLUDED
+#define OPM_NEWTONITERATIONBLACKOILINTERFACE_HEADER_INCLUDED
+
+#include <opm/autodiff/LinearisedBlackoilResidual.hpp>
+
+namespace Opm
+{
+
+    /// Interface class for (linear) solvers for the fully implicit black-oil system.
+    class NewtonIterationBlackoilInterface
+    {
+    public:
+        /// Return type for linearSolve(). A simple, non-ad vector type.
+        typedef LinearisedBlackoilResidual::ADB::V SolutionVector;
+
+        /// Solve the linear system Ax = b, with A being the
+        /// combined derivative matrix of the residual and b
+        /// being the residual itself.
+        /// \param[in] residual   residual object containing A and b.
+        /// \return               the solution x
+        virtual SolutionVector computeNewtonIncrement(const LinearisedBlackoilResidual& residual) const = 0;
+    };
+
+} // namespace Opm
+
+
+#endif // OPM_NEWTONITERATIONBLACKOILINTERFACE_HEADER_INCLUDED

opm/autodiff/NewtonIterationBlackoilSimple.cpp

@@ -0,0 +1,69 @@
+/*
+  Copyright 2014 SINTEF ICT, Applied Mathematics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <config.h>
+
+#include <opm/autodiff/NewtonIterationBlackoilSimple.hpp>
+#include <opm/autodiff/AutoDiffHelpers.hpp>
+#include <opm/core/utility/ErrorMacros.hpp>
+
+namespace Opm
+{
+
+    /// Construct a system solver.
+    /// \param[in] linsolver   linear solver to use
+    NewtonIterationBlackoilSimple::NewtonIterationBlackoilSimple(const LinearSolverInterface& linsolver)
+        : linsolver_(linsolver)
+    {
+    }
+
+    /// Solve the linear system Ax = b, with A being the
+    /// combined derivative matrix of the residual and b
+    /// being the residual itself.
+    /// \param[in] residual   residual object containing A and b.
+    /// \return               the solution x
+    NewtonIterationBlackoilSimple::SolutionVector
+    NewtonIterationBlackoilSimple::computeNewtonIncrement(const LinearisedBlackoilResidual& residual) const
+    {
+        typedef LinearisedBlackoilResidual::ADB ADB;
+        const int np = residual.material_balance_eq.size();
+        ADB mass_res = residual.material_balance_eq[0];
+        for (int phase = 1; phase < np; ++phase) {
+            mass_res = vertcat(mass_res, residual.material_balance_eq[phase]);
+        }
+        const ADB well_res = vertcat(residual.well_flux_eq, residual.well_eq);
+        const ADB total_residual = collapseJacs(vertcat(mass_res, well_res));
+
+        const Eigen::SparseMatrix<double, Eigen::RowMajor> matr = total_residual.derivative()[0];
+
+        SolutionVector dx(SolutionVector::Zero(total_residual.size()));
+        Opm::LinearSolverInterface::LinearSolverReport rep
+            = linsolver_.solve(matr.rows(), matr.nonZeros(),
+                               matr.outerIndexPtr(), matr.innerIndexPtr(), matr.valuePtr(),
+                               total_residual.value().data(), dx.data());
+        if (!rep.converged) {
+            OPM_THROW(std::runtime_error,
+                      "FullyImplicitBlackoilSolver::solveJacobianSystem(): "
+                      "Linear solver convergence failure.");
+        }
+        return dx;
+    }
+
+} // namespace Opm
+

opm/autodiff/NewtonIterationBlackoilSimple.hpp

@@ -0,0 +1,55 @@
+/*
+  Copyright 2014 SINTEF ICT, Applied Mathematics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_NEWTONITERATIONBLACKOILSIMPLE_HEADER_INCLUDED
+#define OPM_NEWTONITERATIONBLACKOILSIMPLE_HEADER_INCLUDED
+
+
+#include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
+#include <opm/core/linalg/LinearSolverInterface.hpp>
+
+namespace Opm
+{
+
+    /// This class solves the fully implicit black-oil system by
+    /// simply concatenating the Jacobian matrices and passing the
+    /// resulting system to a linear solver. The linear solver used
+    /// can be passed in as a constructor argument.
+    class NewtonIterationBlackoilSimple : public NewtonIterationBlackoilInterface
+    {
+    public:
+        /// Construct a system solver.
+        /// \param[in] linsolver   linear solver to use
+        NewtonIterationBlackoilSimple(const LinearSolverInterface& linsolver);
+
+        /// Solve the system of linear equations Ax = b, with A being the
+        /// combined derivative matrix of the residual and b
+        /// being the residual itself.
+        /// \param[in] residual   residual object containing A and b.
+        /// \return               the solution x
+        virtual SolutionVector computeNewtonIncrement(const LinearisedBlackoilResidual& residual) const;
+
+    private:
+        const LinearSolverInterface& linsolver_;
+    };
+
+} // namespace Opm
+
+
+#endif // OPM_NEWTONITERATIONBLACKOILSIMPLE_HEADER_INCLUDED

opm/autodiff/SimulatorFullyImplicitBlackoil.hpp

@@ -33,7 +33,7 @@ namespace Opm
     class BlackoilPropsAdInterface;
     class RockCompressibility;
     class WellsManager;
-    class LinearSolverInterface;
+    class NewtonIterationBlackoilInterface;
     class SimulatorTimer;
     class BlackoilState;
     class WellStateFullyImplicitBlackoil;
@@ -73,7 +73,7 @@ namespace Opm
                                        BlackoilPropsAdInterface& props,
                                        const RockCompressibility* rock_comp_props,
                                        WellsManager& wells_manager,
-                                       LinearSolverInterface& linsolver,
+                                       NewtonIterationBlackoilInterface& linsolver,
                                        const double* gravity);
 
         /// Run the simulation.

opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp

@@ -65,7 +65,7 @@ namespace Opm
              BlackoilPropsAdInterface& props,
              const RockCompressibility* rock_comp_props,
              WellsManager& wells_manager,
-             LinearSolverInterface& linsolver,
+             NewtonIterationBlackoilInterface& linsolver,
              const double* gravity);
 
         SimulatorReport run(SimulatorTimer& timer,
@@ -106,7 +106,7 @@ namespace Opm
                                                                    BlackoilPropsAdInterface& props,
                                                                    const RockCompressibility* rock_comp_props,
                                                                    WellsManager& wells_manager,
-                                                                   LinearSolverInterface& linsolver,
+                                                                   NewtonIterationBlackoilInterface& linsolver,
                                                                    const double* gravity)
 
     {
@@ -191,7 +191,7 @@ namespace Opm
                                                BlackoilPropsAdInterface& props,
                                                const RockCompressibility* rock_comp_props,
                                                WellsManager& wells_manager,
-                                               LinearSolverInterface& linsolver,
+                                               NewtonIterationBlackoilInterface& linsolver,
                                                const double* gravity)
         : grid_(grid),
           props_(props),