A WIP version BlackoilMultiSegmentModel

and also a Simulator Class and example for multisegment wells.

CMakeLists_files.cmake

@@ -78,6 +78,7 @@ list (APPEND TEST_DATA_FILES
 list (APPEND EXAMPLE_SOURCE_FILES
 	examples/find_zero.cpp
 	examples/flow.cpp
+	examples/flow_multisegment.cpp
 	examples/flow_solvent.cpp
 	examples/sim_2p_incomp_ad.cpp
 	examples/sim_simple.cpp
@@ -115,6 +116,8 @@ list (APPEND PUBLIC_HEADER_FILES
 	opm/autodiff/BlackoilSolventModel.hpp
 	opm/autodiff/BlackoilSolventModel_impl.hpp
 	opm/autodiff/BlackoilSolventState.hpp
+	opm/autodiff/BlackoilMultiSegmentModel.hpp
+	opm/autodiff/BlackoilMultiSegmentModel_impl.hpp
 	opm/autodiff/fastSparseProduct.hpp
 	opm/autodiff/DuneMatrix.hpp
 	opm/autodiff/ExtractParallelGridInformationToISTL.hpp
@@ -138,6 +141,8 @@ list (APPEND PUBLIC_HEADER_FILES
 	opm/autodiff/SimulatorFullyImplicitBlackoil.hpp
 	opm/autodiff/SimulatorFullyImplicitBlackoilSolvent.hpp
 	opm/autodiff/SimulatorFullyImplicitBlackoilSolvent_impl.hpp
+	opm/autodiff/SimulatorFullyImplicitBlackoilMultiSegment.hpp
+	opm/autodiff/SimulatorFullyImplicitBlackoilMultiSegment_impl.hpp
 	opm/autodiff/SimulatorIncompTwophaseAd.hpp
 	opm/autodiff/TransportSolverTwophaseAd.hpp
 	opm/autodiff/WellDensitySegmented.hpp

examples/flow_multisegment.cpp

@@ -0,0 +1,418 @@
+/*
+  Copyright 2013 SINTEF ICT, Applied Mathematics.
+  Copyright 2014 Dr. Blatt - HPC-Simulation-Software & Services
+  Copyright 2015 IRIS AS
+
+  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/>.
+*/
+
+
+#if HAVE_CONFIG_H
+#include "config.h"
+#endif // HAVE_CONFIG_H
+
+
+#include <dune/common/version.hh>
+
+#include <opm/core/utility/platform_dependent/disable_warnings.h>
+
+#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 3)
+#include <dune/common/parallel/mpihelper.hh>
+#else
+#include <dune/common/mpihelper.hh>
+#endif
+
+#if HAVE_DUNE_CORNERPOINT && WANT_DUNE_CORNERPOINTGRID
+#define USE_DUNE_CORNERPOINTGRID 1
+#include <dune/grid/CpGrid.hpp>
+#include <dune/grid/common/GridAdapter.hpp>
+#else
+#undef USE_DUNE_CORNERPOINTGRID
+#endif
+
+#include <opm/core/utility/platform_dependent/reenable_warnings.h>
+
+#include <opm/core/pressure/FlowBCManager.hpp>
+
+#include <opm/core/grid.h>
+#include <opm/core/grid/cornerpoint_grid.h>
+#include <opm/core/grid/GridManager.hpp>
+#include <opm/autodiff/GridHelpers.hpp>
+#include <opm/autodiff/createGlobalCellArray.hpp>
+
+#include <opm/core/wells.h>
+#include <opm/core/wells/WellsManager.hpp>
+#include <opm/core/utility/ErrorMacros.hpp>
+#include <opm/core/simulator/initState.hpp>
+#include <opm/core/simulator/initStateEquil.hpp>
+#include <opm/core/simulator/SimulatorReport.hpp>
+#include <opm/core/simulator/SimulatorTimer.hpp>
+#include <opm/core/utility/miscUtilities.hpp>
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
+#include <opm/core/utility/thresholdPressures.hpp> // Note: the GridHelpers must be included before this (to make overloads available). \TODO: Fix.
+
+#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
+#include <opm/core/props/BlackoilPropertiesBasic.hpp>
+#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
+#include <opm/core/props/rock/RockCompressibility.hpp>
+
+#include <opm/core/linalg/LinearSolverFactory.hpp>
+#include <opm/autodiff/NewtonIterationBlackoilSimple.hpp>
+#include <opm/autodiff/NewtonIterationBlackoilCPR.hpp>
+#include <opm/autodiff/NewtonIterationBlackoilInterleaved.hpp>
+
+#include <opm/core/simulator/BlackoilState.hpp>
+#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
+
+#include <opm/autodiff/SimulatorFullyImplicitBlackoilMultiSegment.hpp>
+#include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
+#include <opm/autodiff/RedistributeDataHandles.hpp>
+
+#include <opm/core/utility/share_obj.hpp>
+
+#include <opm/parser/eclipse/OpmLog/OpmLog.hpp>
+#include <opm/parser/eclipse/OpmLog/StreamLog.hpp>
+#include <opm/parser/eclipse/OpmLog/CounterLog.hpp>
+#include <opm/parser/eclipse/Deck/Deck.hpp>
+#include <opm/parser/eclipse/Parser/Parser.hpp>
+#include <opm/parser/eclipse/Parser/ParseMode.hpp>
+#include <opm/parser/eclipse/EclipseState/checkDeck.hpp>
+#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
+
+#include <boost/filesystem.hpp>
+#include <boost/algorithm/string.hpp>
+
+#include <memory>
+#include <algorithm>
+#include <cstdlib>
+#include <iostream>
+#include <vector>
+#include <numeric>
+#include <cstdlib>
+
+namespace
+{
+    void warnIfUnusedParams(const Opm::parameter::ParameterGroup& param)
+    {
+        if (param.anyUnused()) {
+            std::cout << "--------------------   Unused parameters:   --------------------\n";
+            param.displayUsage();
+            std::cout << "----------------------------------------------------------------" << std::endl;
+        }
+    }
+} // anon namespace
+
+
+
+// ----------------- Main program -----------------
+int
+main(int argc, char** argv)
+try
+{
+    using namespace Opm;
+#if USE_DUNE_CORNERPOINTGRID
+    // Must ensure an instance of the helper is created to initialise MPI.
+    const Dune::MPIHelper& mpi_helper = Dune::MPIHelper::instance(argc, argv);
+    const int mpi_rank = mpi_helper.rank();
+    const int mpi_size = mpi_helper.size();
+#else
+    // default values for serial run
+    const int mpi_rank = 0;
+    const int mpi_size = 1;
+#endif
+
+    // Write parameters used for later reference. (only if rank is zero)
+    const bool output_cout = ( mpi_rank == 0 );
+
+    if(output_cout)
+    {
+        std::cout << "**********************************************************************\n";
+        std::cout << "*                                                                    *\n";
+        std::cout << "*                   This is Flow (version 2015.04)                   *\n";
+        std::cout << "*                                                                    *\n";
+        std::cout << "* Flow is a simulator for fully implicit three-phase black-oil flow, *\n";
+        std::cout << "*            and is part of OPM. For more information see:           *\n";
+        std::cout << "*                       http://opm-project.org                       *\n";
+        std::cout << "*                                                                    *\n";
+        std::cout << "**********************************************************************\n\n";
+    }
+
+    // Read parameters, see if a deck was specified on the command line.
+    if ( output_cout )
+    {
+        std::cout << "---------------    Reading parameters     ---------------" << std::endl;
+    }
+
+    parameter::ParameterGroup param(argc, argv, false, output_cout);
+    if( !output_cout )
+    {
+        param.disableOutput();
+    }
+
+    if (!param.unhandledArguments().empty()) {
+        if (param.unhandledArguments().size() != 1) {
+            std::cerr << "You can only specify a single input deck on the command line.\n";
+            return EXIT_FAILURE;
+        } else {
+            param.insertParameter("deck_filename", param.unhandledArguments()[0]);
+        }
+    }
+
+    // We must have an input deck. Grid and props will be read from that.
+    if (!param.has("deck_filename")) {
+        std::cerr << "This program must be run with an input deck.\n"
+            "Specify the deck filename either\n"
+            "    a) as a command line argument by itself\n"
+            "    b) as a command line parameter with the syntax deck_filename=<path to your deck>, or\n"
+            "    c) as a parameter in a parameter file (.param or .xml) passed to the program.\n";
+        return EXIT_FAILURE;
+    }
+
+    // bool check_well_controls = false;
+    // int max_well_control_iterations = 0;
+    double gravity[3] = { 0.0 };
+    std::string deck_filename = param.get<std::string>("deck_filename");
+
+    // Write parameters used for later reference. (only if rank is zero)
+    bool output = ( mpi_rank == 0 ) && param.getDefault("output", true);
+    std::string output_dir;
+    if (output) {
+        // Create output directory if needed.
+        output_dir =
+            param.getDefault("output_dir", std::string("output"));
+        boost::filesystem::path fpath(output_dir);
+        try {
+            create_directories(fpath);
+        }
+        catch (...) {
+            std::cerr << "Creating directories failed: " << fpath << std::endl;
+            return EXIT_FAILURE;
+        }
+        // Write simulation parameters.
+        param.writeParam(output_dir + "/simulation.param");
+    }
+
+    std::string logFile = output_dir + "/LOGFILE.txt";
+    Opm::ParserPtr parser(new Opm::Parser());
+    {
+        std::shared_ptr<Opm::StreamLog> streamLog = std::make_shared<Opm::StreamLog>(logFile , Opm::Log::DefaultMessageTypes);
+        std::shared_ptr<Opm::CounterLog> counterLog = std::make_shared<Opm::CounterLog>(Opm::Log::DefaultMessageTypes);
+
+        Opm::OpmLog::addBackend( "STREAM" , streamLog );
+        Opm::OpmLog::addBackend( "COUNTER" , counterLog );
+    }
+
+    Opm::ParseMode parseMode({{ ParseMode::PARSE_RANDOM_SLASH , InputError::IGNORE }});
+    Opm::DeckConstPtr deck;
+    std::shared_ptr<EclipseState> eclipseState;
+    try {
+        deck = parser->parseFile(deck_filename, parseMode);
+        Opm::checkDeck(deck);
+        eclipseState.reset(new Opm::EclipseState(deck , parseMode));
+    }
+    catch (const std::invalid_argument& e) {
+        std::cerr << "Failed to create valid ECLIPSESTATE object. See logfile: " << logFile << std::endl;
+        std::cerr << "Exception caught: " << e.what() << std::endl;
+        return EXIT_FAILURE;
+    }
+
+    std::vector<double> porv = eclipseState->getDoubleGridProperty("PORV")->getData();
+#if USE_DUNE_CORNERPOINTGRID
+    // Dune::CpGrid as grid manager
+    typedef Dune::CpGrid  Grid;
+    // Grid init
+    Grid grid;
+    grid.processEclipseFormat(deck, false, false, false, porv);
+#else
+    // UnstructuredGrid as grid manager
+    typedef UnstructuredGrid  Grid;
+    GridManager gridManager( eclipseState->getEclipseGrid(), porv );
+    const Grid& grid = *(gridManager.c_grid());
+#endif
+
+    // Possibly override IOConfig setting (from deck) for how often RESTART files should get written to disk (every N report step)
+    if (param.has("output_interval")) {
+        int output_interval = param.get<int>("output_interval");
+        IOConfigPtr ioConfig = eclipseState->getIOConfig();
+        ioConfig->overrideRestartWriteInterval((size_t)output_interval);
+    }
+
+    const PhaseUsage pu = Opm::phaseUsageFromDeck(deck);
+
+    std::vector<int> compressedToCartesianIdx;
+    Opm::createGlobalCellArray(grid, compressedToCartesianIdx);
+
+    typedef BlackoilPropsAdFromDeck::MaterialLawManager MaterialLawManager;
+    auto materialLawManager = std::make_shared<MaterialLawManager>();
+    materialLawManager->initFromDeck(deck, eclipseState, compressedToCartesianIdx);
+
+    // Rock and fluid init
+    BlackoilPropertiesFromDeck props( deck, eclipseState, materialLawManager,
+                                      Opm::UgGridHelpers::numCells(grid),
+                                      Opm::UgGridHelpers::globalCell(grid),
+                                      Opm::UgGridHelpers::cartDims(grid),
+                                      param);
+
+    BlackoilPropsAdFromDeck new_props( deck, eclipseState, materialLawManager, grid );
+    // check_well_controls = param.getDefault("check_well_controls", false);
+    // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
+    // Rock compressibility.
+
+    RockCompressibility rock_comp(deck, eclipseState);
+
+    // Gravity.
+    gravity[2] = deck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity;
+
+    BlackoilState state;
+    // Init state variables (saturation and pressure).
+    if (param.has("init_saturation")) {
+        initStateBasic(Opm::UgGridHelpers::numCells(grid),
+                       Opm::UgGridHelpers::globalCell(grid),
+                       Opm::UgGridHelpers::cartDims(grid),
+                       Opm::UgGridHelpers::numFaces(grid),
+                       Opm::UgGridHelpers::faceCells(grid),
+                       Opm::UgGridHelpers::beginFaceCentroids(grid),
+                       Opm::UgGridHelpers::beginCellCentroids(grid),
+                       Opm::UgGridHelpers::dimensions(grid),
+                       props, param, gravity[2], state);
+
+        initBlackoilSurfvol(Opm::UgGridHelpers::numCells(grid), props, state);
+
+        enum { Oil = BlackoilPhases::Liquid, Gas = BlackoilPhases::Vapour };
+        if (pu.phase_used[Oil] && pu.phase_used[Gas]) {
+            const int numPhases = props.numPhases();
+            const int numCells  = Opm::UgGridHelpers::numCells(grid);
+            for (int c = 0; c < numCells; ++c) {
+                state.gasoilratio()[c] = state.surfacevol()[c*numPhases + pu.phase_pos[Gas]]
+                    / state.surfacevol()[c*numPhases + pu.phase_pos[Oil]];
+            }
+        }
+    } else if (deck->hasKeyword("EQUIL") && props.numPhases() == 3) {
+        state.init(Opm::UgGridHelpers::numCells(grid),
+                   Opm::UgGridHelpers::numFaces(grid),
+                   props.numPhases());
+        const double grav = param.getDefault("gravity", unit::gravity);
+        initStateEquil(grid, props, deck, eclipseState, grav, state);
+        state.faceflux().resize(Opm::UgGridHelpers::numFaces(grid), 0.0);
+    } else {
+        initBlackoilStateFromDeck(Opm::UgGridHelpers::numCells(grid),
+                                  Opm::UgGridHelpers::globalCell(grid),
+                                  Opm::UgGridHelpers::numFaces(grid),
+                                  Opm::UgGridHelpers::faceCells(grid),
+                                  Opm::UgGridHelpers::beginFaceCentroids(grid),
+                                  Opm::UgGridHelpers::beginCellCentroids(grid),
+                                  Opm::UgGridHelpers::dimensions(grid),
+                                  props, deck, gravity[2], state);
+    }
+
+
+    // The capillary pressure is scaled in new_props to match the scaled capillary pressure in props.
+    if (deck->hasKeyword("SWATINIT")) {
+        const int numCells = Opm::UgGridHelpers::numCells(grid);
+        std::vector<int> cells(numCells);
+        for (int c = 0; c < numCells; ++c) { cells[c] = c; }
+        std::vector<double> pc = state.saturation();
+        props.capPress(numCells, state.saturation().data(), cells.data(), pc.data(),NULL);
+        new_props.setSwatInitScaling(state.saturation(),pc);
+    }
+
+    bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0);
+    const double *grav = use_gravity ? &gravity[0] : 0;
+
+    const bool use_local_perm = param.getDefault("use_local_perm", true);
+
+    DerivedGeology geoprops(grid, new_props, eclipseState, use_local_perm, grav);
+    boost::any parallel_information;
+
+    // At this point all properties and state variables are correctly initialized
+    // If there are more than one processors involved, we now repartition the grid
+    // and initilialize new properties and states for it.
+    if( mpi_size > 1 )
+    {
+        Opm::distributeGridAndData( grid, eclipseState, state, new_props, geoprops, parallel_information, use_local_perm );
+    }
+
+    // create output writer after grid is distributed, otherwise the parallel output
+    // won't work correctly since we need to create a mapping from the distributed to
+    // the global view
+    Opm::BlackoilOutputWriter outputWriter(grid, param, eclipseState, pu );
+
+    // Solver for Newton iterations.
+    std::unique_ptr<NewtonIterationBlackoilInterface> fis_solver;
+    if (param.getDefault("use_interleaved", true)) {
+        fis_solver.reset(new NewtonIterationBlackoilInterleaved(param, parallel_information));
+    } else if (param.getDefault("use_cpr", true)) {
+        fis_solver.reset(new NewtonIterationBlackoilCPR(param, parallel_information));
+    } else {
+        fis_solver.reset(new NewtonIterationBlackoilSimple(param, parallel_information));
+    }
+
+    Opm::ScheduleConstPtr schedule = eclipseState->getSchedule();
+    Opm::TimeMapConstPtr timeMap(schedule->getTimeMap());
+    SimulatorTimer simtimer;
+
+    // initialize variables
+    simtimer.init(timeMap);
+
+    std::vector<double> threshold_pressures = thresholdPressures(parseMode, eclipseState, grid);
+
+    SimulatorFullyImplicitBlackoilMultiSegment< Grid >  simulator(param,
+                                                      grid,
+                                                      geoprops,
+                                                      new_props,
+                                                      rock_comp.isActive() ? &rock_comp : 0,
+                                                      *fis_solver,
+                                                      grav,
+                                                      deck->hasKeyword("DISGAS"),
+                                                      deck->hasKeyword("VAPOIL"),
+                                                      eclipseState,
+                                                      outputWriter,
+                                                      threshold_pressures);
+
+    if (!schedule->initOnly()){
+        if( output_cout )
+        {
+            std::cout << "\n\n================ Starting main simulation loop ===============\n"
+                      << std::flush;
+        }
+
+        SimulatorReport fullReport = simulator.run(simtimer, state);
+
+        if( output_cout )
+        {
+            std::cout << "\n\n================    End of simulation     ===============\n\n";
+            fullReport.reportFullyImplicit(std::cout);
+        }
+
+        if (output) {
+            std::string filename = output_dir + "/walltime.txt";
+            std::fstream tot_os(filename.c_str(),std::fstream::trunc | std::fstream::out);
+            fullReport.reportParam(tot_os);
+            warnIfUnusedParams(param);
+        }
+    } else {
+        outputWriter.writeInit( simtimer );
+        if ( output_cout )
+        {
+            std::cout << "\n\n================ Simulation turned off ===============\n" << std::flush;
+        }
+    }
+}
+catch (const std::exception &e) {
+    std::cerr << "Program threw an exception: " << e.what() << "\n";
+    return EXIT_FAILURE;
+}
+

opm/autodiff/BlackoilModel.hpp

@@ -27,7 +27,6 @@
 #include <opm/core/simulator/BlackoilState.hpp>
 #include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
 #include <opm/autodiff/BlackoilModelParameters.hpp>
-#include <opm/autodiff/WellStateMultiSegment.hpp>
 
 namespace Opm {
 

opm/autodiff/BlackoilMultiSegmentModel.hpp

@@ -0,0 +1,567 @@
+/*
+  Copyright 2013, 2015 SINTEF ICT, Applied Mathematics.
+  Copyright 2014, 2015 Statoil ASA.
+  Copyright 2014, 2015 Dr. Markus Blatt - HPC-Simulation-Software & Services
+  Copyright 2015 NTNU
+
+  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_BLACKOILMULTISEGMENTMODEL_HEADER_INCLUDED
+#define OPM_BLACKOILMULTISEGMENTMODEL_HEADER_INCLUDED
+
+#include <opm/core/simulator/BlackoilState.hpp>
+#include <opm/autodiff/BlackoilModelBase.hpp>
+#include <opm/autodiff/BlackoilModelParameters.hpp>
+#include <opm/autodiff/WellStateMultiSegment.hpp>
+#include <opm/autodiff/WellMultiSegment.hpp>
+
+namespace Opm {
+
+    struct BlackoilMultiSegmentSolutionState : public DefaultBlackoilSolutionState
+    {
+        explicit BlackoilMultiSegmentSolutionState(const int np)
+            : DefaultBlackoilSolutionState(np),
+              pseg( ADB::null())
+        {
+        }
+        ADB pseg; // the segment pressures
+    };
+
+    /// A model implementation for three-phase black oil with support
+    /// for multi-segment wells.
+    ///
+    /// It uses automatic differentiation via the class AutoDiffBlock
+    /// to simplify assembly of the jacobian matrix.
+    /// \tparam  Grid            UnstructuredGrid or CpGrid.
+    /// \tparam  Implementation  Provides concrete state types.
+    template<class Grid>
+    class BlackoilMultiSegmentModel : public BlackoilModelBase<Grid, BlackoilMultiSegmentModel<Grid>>
+    {
+    public:
+
+        typedef BlackoilModelBase<Grid, BlackoilMultiSegmentModel<Grid> > Base; // base class
+        typedef typename Base::ReservoirState ReservoirState;
+        typedef BlackoilMultiSegmentSolutionState SolutionState;
+
+        // ---------  Public methods  ---------
+
+        /// Construct the model. It will retain references to the
+        /// arguments of this functions, and they are expected to
+        /// remain in scope for the lifetime of the solver.
+        /// \param[in] param            parameters
+        /// \param[in] grid             grid data structure
+        /// \param[in] fluid            fluid properties
+        /// \param[in] geo              rock properties
+        /// \param[in] rock_comp_props  if non-null, rock compressibility properties
+        /// \param[in] wells            well structure
+        /// \param[in] vfp_properties   Vertical flow performance tables
+        /// \param[in] linsolver        linear solver
+        /// \param[in] eclState         eclipse state
+        /// \param[in] has_disgas       turn on dissolved gas
+        /// \param[in] has_vapoil       turn on vaporized oil feature
+        /// \param[in] terminal_output  request output to cout/cerr
+        BlackoilMultiSegmentModel(const typename Base::ModelParameters&  param,
+                          const Grid&                     grid ,
+                          const BlackoilPropsAdInterface& fluid,
+                          const DerivedGeology&           geo  ,
+                          const RockCompressibility*      rock_comp_props,
+                          const Wells*                    wells,
+                          const NewtonIterationBlackoilInterface& linsolver,
+                          Opm::EclipseStateConstPtr eclState,
+                          const bool has_disgas,
+                          const bool has_vapoil,
+                          const bool terminal_output);
+
+        /// \brief Set threshold pressures that prevent or reduce flow.
+        /// This prevents flow across faces if the potential
+        /// difference is less than the threshold. If the potential
+        /// difference is greater, the threshold value is subtracted
+        /// before calculating flow. This is treated symmetrically, so
+        /// flow is prevented or reduced in both directions equally.
+        /// \param[in]  threshold_pressures_by_face   array of size equal to the number of faces
+        ///                                   of the grid passed in the constructor.
+        /* void setThresholdPressures(const std::vector<double>& threshold_pressures_by_face);
+        oafafd
+
+        /// Called once before each time step.
+        /// \param[in] dt                     time step size
+        /// \param[in, out] reservoir_state   reservoir state variables
+        /// \param[in, out] well_state        well state variables
+        void prepareStep(const double dt,
+                         ReservoirState& reservoir_state,
+                         WellState& well_state);
+
+        /// Called once after each time step.
+        /// In this class, this function does nothing.
+        /// \param[in] dt                     time step size
+        /// \param[in, out] reservoir_state   reservoir state variables
+        /// \param[in, out] well_state        well state variables
+        void afterStep(const double dt,
+                       ReservoirState& reservoir_state,
+                       WellState& well_state);
+
+        /// Assemble the residual and Jacobian of the nonlinear system.
+        /// \param[in]      reservoir_state   reservoir state variables
+        /// \param[in, out] well_state        well state variables
+        /// \param[in]      initial_assembly  pass true if this is the first call to assemble() in this timestep
+        void assemble(const ReservoirState& reservoir_state,
+                      WellState& well_state,
+                      const bool initial_assembly);
+
+        void assemble(const ReservoirState& reservoir_state,
+                      WellStateMultiSegment& well_state,
+                      const bool initial_assembly);
+
+        /// \brief Compute the residual norms of the mass balance for each phase,
+        /// the well flux, and the well equation.
+        /// \return a vector that contains for each phase the norm of the mass balance
+        /// and afterwards the norm of the residual of the well flux and the well equation.
+        std::vector<double> computeResidualNorms() const;
+
+        /// The size (number of unknowns) of the nonlinear system of equations.
+        int sizeNonLinear() const;
+
+        /// Number of linear iterations used in last call to solveJacobianSystem().
+        int linearIterationsLastSolve() const;
+
+        /// Solve the Jacobian system Jx = r where J is the Jacobian and
+        /// r is the residual.
+        V solveJacobianSystem() const;
+
+        /// Apply an update to the primary variables, chopped if appropriate.
+        /// \param[in]      dx                updates to apply to primary variables
+        /// \param[in, out] reservoir_state   reservoir state variables
+        /// \param[in, out] well_state        well state variables
+        void updateState(const V& dx,
+                         ReservoirState& reservoir_state,
+                         WellState& well_state);
+
+        /// Return true if output to cout is wanted.
+        bool terminalOutputEnabled() const;
+
+        /// Compute convergence based on total mass balance (tol_mb) and maximum
+        /// residual mass balance (tol_cnv).
+        /// \param[in]   dt          timestep length
+        /// \param[in]   iteration   current iteration number
+        bool getConvergence(const double dt, const int iteration);
+
+        /// The number of active phases in the model.
+        int numPhases() const;
+
+        /// Update the scaling factors for mass balance equations
+        void updateEquationsScaling();
+
+    protected:
+
+        // ---------  Types and enums  ---------
+
+        typedef Eigen::Array<double,
+                             Eigen::Dynamic,
+                             Eigen::Dynamic,
+                             Eigen::RowMajor> DataBlock;
+
+        struct ReservoirResidualQuant {
+            ReservoirResidualQuant();
+            std::vector<ADB> accum; // Accumulations
+            ADB              mflux; // Mass flux (surface conditions)
+            ADB              b;     // Reciprocal FVF
+            ADB              dh;    // Pressure drop across int. interfaces
+            ADB              mob;   // Phase mobility (per cell)
+        };
+
+        struct WellOps {
+            WellOps(const Wells* wells);
+            Eigen::SparseMatrix<double> w2p;              // well -> perf (scatter)
+            Eigen::SparseMatrix<double> p2w;              // perf -> well (gather)
+        };
+
+        // ---------  Data members  ---------
+
+        const Grid&         grid_;
+        const BlackoilPropsAdInterface& fluid_;
+        const DerivedGeology&           geo_;
+        const RockCompressibility*      rock_comp_props_;
+        const Wells*                    wells_;
+        // FOR TEMPORARY
+        // SHOUlD BE A REFERENCE
+        const std::vector<WellMultiSegment> wells_multi_segment_;
+        VFPProperties                   vfp_properties_;
+        const NewtonIterationBlackoilInterface&    linsolver_;
+        // For each canonical phase -> true if active
+        const std::vector<bool>         active_;
+        // Size = # active phases. Maps active -> canonical phase indices.
+        const std::vector<int>          canph_;
+        const std::vector<int>          cells_;  // All grid cells
+        HelperOps                       ops_;
+        const WellOps                   wops_;
+        const bool has_disgas_;
+        const bool has_vapoil_;
+
+        ModelParameters                 param_;
+        bool use_threshold_pressure_;
+        bool wells_active_;
+        V threshold_pressures_by_interior_face_;
+
+        std::vector<ReservoirResidualQuant> rq_;
+        std::vector<PhasePresence> phaseCondition_;
+        V isRs_;
+        V isRv_;
+        V isSg_;
+
+        // For the non-segmented well, it should be the density with AVG or SEG way.
+        // while usually SEG way
+        V well_perforation_densities_; //Density of each well perforation
+
+        // ADB version, when using AVG way, the calculation of the density and hydrostatic head
+        // is implicit
+        ADB well_perforation_densities_adb_;
+
+        // Diff to the pressure of the related segment.
+        // When the well is a usual well, the bhp will be the pressure of the top segment
+        // For mutlti-segmented wells, only AVG is allowed.
+        // For non-segmented wells, typically SEG is used. AVG way might not have been
+        // implemented yet.
+
+        // Diff to bhp for each well perforation. only for usual wells.
+        // For segmented wells, they are zeros.
+        V well_perforation_pressure_diffs_; // Diff to bhp for each well perforation.
+
+        // ADB version. Eventually, only ADB version will be kept.
+        ADB well_perforation_pressure_diffs_adb_;
+
+        // Pressure correction due to the different depth of the perforation
+        // and the cell center of the grid block
+        // For the non-segmented wells, since the perforation are forced to be
+        // at the center of the grid cell, it should be ZERO.
+        // It should only apply to the mutli-segmented wells.
+        V well_perforation_pressure_cell_diffs_;
+        ADB well_perforation_pressure_cell_diffs_adb_;
+
+        // Pressure correction due to the depth differennce between segment depth and perforation depth.
+        // TODO: It should be able to be merge as a part of the perforation_pressure_diffs_.
+        ADB well_perforations_segment_pressure_diffs_;
+
+        // the average of the fluid densities in the grid block
+        // which is used to calculate the hydrostatic head correction due to the depth difference of the perforation
+        // and the cell center of the grid block
+        V well_perforation_cell_densities_;
+        ADB well_perforation_cell_densities_adb_;
+
+        V well_perforatoin_cell_pressure_diffs_;
+
+        LinearisedBlackoilResidual residual_;
+
+        /// \brief Whether we print something to std::cout
+        bool terminal_output_;
+
+        std::vector<int>         primalVariable_;
+        V pvdt_;
+
+        // ---------  Protected methods  ---------
+
+        /// Access the most-derived class used for
+        /// static polymorphism (CRTP).
+        Implementation& asImpl()
+        {
+            return static_cast<Implementation&>(*this);
+        }
+
+        /// Access the most-derived class used for
+        /// static polymorphism (CRTP).
+        const Implementation& asImpl() const
+        {
+            return static_cast<const Implementation&>(*this);
+        }
+
+        // return true if wells are available in the reservoir
+        bool wellsActive() const { return wells_active_; }
+        // return true if wells are available on this process
+        bool localWellsActive() const { return wells_ ? (wells_->number_of_wells > 0 ) : false; }
+
+        // return wells object
+        const Wells& wells () const { assert( bool(wells_ != 0) ); return *wells_; }
+        const std::vector<WellMultiSegment>& wellsMultiSegment() const { return wells_multi_segment_; }
+
+        void
+        makeConstantState(SolutionState& state) const;
+
+        SolutionState
+        variableState(const ReservoirState& x,
+                      const WellState& xw) const;
+
+        SolutionState
+        variableState(const ReservoirState& x,
+                      const WellStateMultiSegment& xw) const;
+
+        std::vector<V>
+        variableStateInitials(const ReservoirState& x,
+                              const WellState& xw) const;
+        std::vector<V>
+        variableStateInitials(const ReservoirState& x,
+                              const WellStateMultiSegment& xw) const;
+        void
+        variableReservoirStateInitials(const ReservoirState& x,
+                                       std::vector<V>& vars0) const;
+        void
+        variableWellStateInitials(const WellState& xw,
+                                  std::vector<V>& vars0) const;
+
+        void variableWellStateInitials(const WellStateMultiSegment& xw,
+                                       std::vector<V>& vars0) const;
+
+        void
+        variableWellState(const WellStateMultiSegment& xw,
+                                  std::vector<V>& vars0) const;
+
+        std::vector<int>
+        variableStateIndices() const;
+
+        std::vector<int>
+        variableWellStateIndices() const;
+
+        SolutionState
+        variableStateExtractVars(const ReservoirState& x,
+                                 const std::vector<int>& indices,
+                                 std::vector<ADB>& vars) const;
+
+        void
+        variableStateExtractWellsVars(const std::vector<int>& indices,
+                                      std::vector<ADB>& vars,
+                                      SolutionState& state) const;
+
+        void
+        computeAccum(const SolutionState& state,
+                     const int            aix  );
+
+        void computeWellConnectionPressures(const SolutionState& state,
+                                            const WellState& xw);
+
+        void computeWellConnectionPressures(const SolutionState& state,
+                                            const WellStateMultiSegment& xw);
+
+        void
+        assembleMassBalanceEq(const SolutionState& state);
+
+        void
+        solveWellEq(const std::vector<ADB>& mob_perfcells,
+                    const std::vector<ADB>& b_perfcells,
+                    SolutionState& state,
+                    WellState& well_state);
+
+        void
+        computeWellFlux(const SolutionState& state,
+                        const std::vector<ADB>& mob_perfcells,
+                        const std::vector<ADB>& b_perfcells,
+                        V& aliveWells,
+                        std::vector<ADB>& cq_s);
+
+        void
+        computeWellFlux(const MultiSegmentBlackoilSolutionState& state,
+                        const std::vector<ADB>& mob_perfcells,
+                        const std::vector<ADB>& b_perfcells,
+                        V& aliveWells,
+                        std::vector<ADB>& cq_s);
+        void
+        updatePerfPhaseRatesAndPressures(const std::vector<ADB>& cq_s,
+                                         const SolutionState& state,
+                                         WellState& xw);
+
+        void
+        updatePerfPhaseRatesAndPressures(const std::vector<ADB>& cq_s,
+                                         const MultiSegmentBlackoilSolutionState& state,
+                                         WellStateMultiSegment& xw);
+
+        void
+        addWellFluxEq(const std::vector<ADB>& cq_s,
+                      const SolutionState& state);
+
+        void
+        addWellFluxEq(const std::vector<ADB>& cq_s,
+                      const MultiSegmentBlackoilSolutionState& state);
+
+        void
+        addWellContributionToMassBalanceEq(const std::vector<ADB>& cq_s,
+                                           const SolutionState& state,
+                                           const WellState& xw);
+
+        void
+        addWellControlEq(const SolutionState& state,
+                         const WellState& xw,
+                         const V& aliveWells);
+
+        void
+        addWellControlEq(const MultiSegmentBlackoilSolutionState& state,
+                         const WellStateMultiSegment& xw,
+                         const V& aliveWells);
+
+
+        void updateWellControls(WellState& xw) const;
+        void updateWellControls(WellStateMultiSegment& xw) const;
+
+        void updateWellState(const V& dwells,
+                             WellState& well_state);
+
+        bool getWellConvergence(const int iteration);
+
+        bool isVFPActive() const;
+
+        std::vector<ADB>
+        computePressures(const ADB& po,
+                         const ADB& sw,
+                         const ADB& so,
+                         const ADB& sg) const;
+
+        V
+        computeGasPressure(const V& po,
+                           const V& sw,
+                           const V& so,
+                           const V& sg) const;
+
+        std::vector<ADB>
+        computeRelPerm(const SolutionState& state) const;
+
+        void
+        computeMassFlux(const int               actph ,
+                        const V&                transi,
+                        const ADB&              kr    ,
+                        const ADB&              p     ,
+                        const SolutionState&    state );
+
+        void applyThresholdPressures(ADB& dp);
+
+        ADB
+        fluidViscosity(const int               phase,
+                       const ADB&              p    ,
+                       const ADB&              temp ,
+                       const ADB&              rs   ,
+                       const ADB&              rv   ,
+                       const std::vector<PhasePresence>& cond) const;
+
+        ADB
+        fluidReciprocFVF(const int               phase,
+                         const ADB&              p    ,
+                         const ADB&              temp ,
+                         const ADB&              rs   ,
+                         const ADB&              rv   ,
+                         const std::vector<PhasePresence>& cond) const;
+
+        ADB
+        fluidDensity(const int  phase,
+                     const ADB& b,
+                     const ADB& rs,
+                     const ADB& rv) const;
+
+        V
+        fluidRsSat(const V&                p,
+                   const V&                so,
+                   const std::vector<int>& cells) const;
+
+        ADB
+        fluidRsSat(const ADB&              p,
+                   const ADB&              so,
+                   const std::vector<int>& cells) const;
+
+        V
+        fluidRvSat(const V&                p,
+                   const V&                so,
+                   const std::vector<int>& cells) const;
+
+        ADB
+        fluidRvSat(const ADB&              p,
+                   const ADB&              so,
+                   const std::vector<int>& cells) const;
+
+        ADB
+        poroMult(const ADB& p) const;
+
+        ADB
+        transMult(const ADB& p) const;
+
+        const std::vector<PhasePresence>
+        phaseCondition() const {return phaseCondition_;}
+
+        void
+        classifyCondition(const ReservoirState& state);
+
+
+        /// update the primal variable for Sg, Rv or Rs. The Gas phase must
+        /// be active to call this method.
+        void
+        updatePrimalVariableFromState(const ReservoirState& state);
+
+        /// Update the phaseCondition_ member based on the primalVariable_ member.
+        /// Also updates isRs_, isRv_ and isSg_;
+        void
+        updatePhaseCondFromPrimalVariable();
+
+        /// \brief Compute the reduction within the convergence check.
+        /// \param[in] B     A matrix with MaxNumPhases columns and the same number rows
+        ///                  as the number of cells of the grid. B.col(i) contains the values
+        ///                  for phase i.
+        /// \param[in] tempV A matrix with MaxNumPhases columns and the same number rows
+        ///                  as the number of cells of the grid. tempV.col(i) contains the
+        ///                   values
+        ///                  for phase i.
+        /// \param[in] R     A matrix with MaxNumPhases columns and the same number rows
+        ///                  as the number of cells of the grid. B.col(i) contains the values
+        ///                  for phase i.
+        /// \param[out] R_sum An array of size MaxNumPhases where entry i contains the sum
+        ///                   of R for the phase i.
+        /// \param[out] maxCoeff An array of size MaxNumPhases where entry i contains the
+        ///                   maximum of tempV for the phase i.
+        /// \param[out] B_avg An array of size MaxNumPhases where entry i contains the average
+        ///                   of B for the phase i.
+        /// \param[out] maxNormWell The maximum of the well equations for each phase.
+        /// \param[in]  nc    The number of cells of the local grid.
+        /// \param[in]  nw    The number of wells on the local grid.
+        /// \return The total pore volume over all cells.
+        double
+        convergenceReduction(const Eigen::Array<double, Eigen::Dynamic, MaxNumPhases>& B,
+                             const Eigen::Array<double, Eigen::Dynamic, MaxNumPhases>& tempV,
+                             const Eigen::Array<double, Eigen::Dynamic, MaxNumPhases>& R,
+                             std::array<double,MaxNumPhases>& R_sum,
+                             std::array<double,MaxNumPhases>& maxCoeff,
+                             std::array<double,MaxNumPhases>& B_avg,
+                             std::vector<double>& maxNormWell,
+                             int nc,
+                             int nw) const;
+
+        double dpMaxRel() const { return param_.dp_max_rel_; }
+        double dsMax() const { return param_.ds_max_; }
+        double drMaxRel() const { return param_.dr_max_rel_; }
+        double maxResidualAllowed() const { return param_.max_residual_allowed_; } */
+
+    };
+
+    /// Providing types by template specialisation of ModelTraits for BlackoilMultiSegmentModel.
+    template <class Grid>
+    struct ModelTraits< BlackoilMultiSegmentModel<Grid> >
+    {
+        typedef BlackoilState ReservoirState;
+        typedef WellStateMultiSegment WellState;
+        typedef BlackoilModelParameters ModelParameters;
+        typedef BlackoilMultiSegmentSolutionState SolutionState;
+    };
+
+
+
+
+} // namespace Opm
+
+#include "BlackoilMultiSegmentModel_impl.hpp"
+
+#endif // OPM_BLACKOILMULTISEGMENTMODEL_HEADER_INCLUDED

opm/autodiff/SimulatorBase_impl.hpp

@@ -20,7 +20,6 @@
 */
 
 #include <algorithm>
-#include <opm/autodiff/WellMultiSegment.hpp>
 
 namespace Opm
 {
@@ -138,24 +137,6 @@ namespace Opm
             WellState well_state;
             well_state.init(wells, state, prev_well_state);
 
-            const std::vector<WellConstPtr>& wells_ecl = eclipse_state_->getSchedule()->getWells(timer.currentStepNum());
-
-            std::vector<WellMutliSegmentPtr> wells_multisegment(wells_ecl.size());
-
-            // wells_multisegment.resize(wells_ecl.size());
-
-            for (size_t i = 0; i < wells_multisegment.size(); ++i) {
-                wells_multisegment[i].reset(new WellMultiSegment(wells_ecl[i], timer.currentStepNum(), wells));
-            }
-
-            // for DEBUGGING OUTPUT
-            std::cout << " the number of the wells from EclipseState " << wells_ecl.size() << std::endl;
-            for (size_t i = 0; i < wells_ecl.size(); ++i) {
-                std::cout << " well name " << wells_ecl[i]->name() << std::endl;
-                std::cout << " segment wells " << wells_ecl[i]->isMultiSegment() << std::endl;
-            }
-            std::cin.ignore();
-
             // give the polymer and surfactant simulators the chance to do their stuff
             asImpl().handleAdditionalWellInflow(timer, wells_manager, well_state, wells);
 

opm/autodiff/SimulatorFullyImplicitBlackoilMultiSegment.hpp

@@ -0,0 +1,94 @@
+/*
+  Copyright 2013 SINTEF ICT, Applied Mathematics.
+  Copyright 2015 Andreas Lauser
+
+  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_SIMULATORFULLYIMPLICITBLACKOILMULTISEGMENT_HEADER_INCLUDED
+#define OPM_SIMULATORFULLYIMPLICITBLACKOILMULTISEGMENT_HEADER_INCLUDED
+
+#include "SimulatorBase.hpp"
+
+#include "NewtonSolver.hpp"
+
+#include <opm/autodiff/BlackoilMultiSegmentModel.hpp>
+#include <opm/autodiff/WellStateMultiSegment.hpp>
+
+namespace Opm {
+
+template <class GridT>
+class SimulatorFullyImplicitBlackoilMultiSegment;
+
+template <class GridT>
+struct SimulatorTraits<SimulatorFullyImplicitBlackoilMultiSegment<GridT> >
+{
+    typedef WellStateMultiSegment WellState;
+    typedef BlackoilState ReservoirState;
+    typedef BlackoilOutputWriter OutputWriter;
+    typedef GridT Grid;
+    typedef BlackoilMultiSegmentModel<Grid> Model;
+    typedef NewtonSolver<Model> Solver;
+};
+
+/// a simulator for the blackoil model
+template <class GridT>
+class SimulatorFullyImplicitBlackoilMultiSegment
+    : public SimulatorBase<SimulatorFullyImplicitBlackoilMultiSegment<GridT> >
+{
+    typedef SimulatorBase<SimulatorFullyImplicitBlackoilMultiSegment<GridT> > Base;
+    typedef SimulatorFullyImplicitBlackoilMultiSegment<GridT> ThisType;
+    typedef SimulatorTraits<ThisType> Traits;
+    typedef typename Traits::ReservoirState ReservoirState;
+    typedef typename Traits::WellState WellState;
+public:
+    // forward the constructor to the base class
+    SimulatorFullyImplicitBlackoilMultiSegment(const parameter::ParameterGroup& param,
+                                   const GridT& grid,
+                                   const DerivedGeology& geo,
+                                   BlackoilPropsAdInterface& props,
+                                   const RockCompressibility* rock_comp_props,
+                                   NewtonIterationBlackoilInterface& linsolver,
+                                   const double* gravity,
+                                   const bool disgas,
+                                   const bool vapoil,
+                                   std::shared_ptr<EclipseState> eclipse_state,
+                                   BlackoilOutputWriter& output_writer,
+                                   const std::vector<double>& threshold_pressures_by_face)
+    : Base(param, grid, geo, props, rock_comp_props, linsolver, gravity, disgas, vapoil,
+           eclipse_state, output_writer, threshold_pressures_by_face)
+    {}
+
+    SimulatorReport run(SimulatorTimer& timer,
+                        ReservoirState& state);
+
+protected:
+    using Base::output_writer_;
+    using Base::param_;
+    using Base::solver_;
+    using Base::terminal_output_;
+    using Base::eclipse_state_;
+    using Base::grid_;
+    using Base::props_;
+    using Base::is_parallel_run_;
+    using Base::allcells_;
+};
+
+} // namespace Opm
+
+#include "SimulatorFullyImplicitBlackoilMultiSegment_impl.hpp"
+
+#endif // OPM_SIMULATORFULLYIMPLICITBLACKOILMULTISEGMENT_HEADER_INCLUDED

opm/autodiff/SimulatorFullyImplicitBlackoilMultiSegment_impl.hpp

@@ -0,0 +1,175 @@
+/*
+  Copyright 2013 SINTEF ICT, Applied Mathematics.
+  Copyright 2014 IRIS AS
+  Copyright 2015 Andreas Lauser
+
+  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/>.
+*/
+
+
+namespace Opm
+{
+
+    template <class GridT>
+    SimulatorReport SimulatorFullyImplicitBlackoilMultiSegment<GridT>::run(SimulatorTimer& timer,
+                                                                           ReservoirState& state)
+    {
+        WellState prev_well_state;
+
+        // Create timers and file for writing timing info.
+        Opm::time::StopWatch solver_timer;
+        double stime = 0.0;
+        Opm::time::StopWatch step_timer;
+        Opm::time::StopWatch total_timer;
+        total_timer.start();
+        std::string tstep_filename = output_writer_.outputDirectory() + "/step_timing.txt";
+        std::ofstream tstep_os(tstep_filename.c_str());
+
+        // adaptive time stepping
+        std::unique_ptr< AdaptiveTimeStepping > adaptiveTimeStepping;
+        if( param_.getDefault("timestep.adaptive", true ) )
+        {
+            adaptiveTimeStepping.reset( new AdaptiveTimeStepping( param_, solver_.parallelInformation(), terminal_output_ ) );
+        }
+
+        // init output writer
+        output_writer_.writeInit( timer );
+
+        std::string restorefilename = param_.getDefault("restorefile", std::string("") );
+        if( ! restorefilename.empty() )
+        {
+            // -1 means that we'll take the last report step that was written
+            const int desiredRestoreStep = param_.getDefault("restorestep", int(-1) );
+            output_writer_.restore( timer, state, prev_well_state, restorefilename, desiredRestoreStep );
+        }
+
+        unsigned int totalNewtonIterations = 0;
+        unsigned int totalLinearIterations = 0;
+
+        // Main simulation loop.
+        while (!timer.done()) {
+            // Report timestep.
+            step_timer.start();
+            if ( terminal_output_ )
+            {
+                timer.report(std::cout);
+            }
+
+            // Create wells and well state.
+            WellsManager wells_manager(eclipse_state_,
+                                       timer.currentStepNum(),
+                                       Opm::UgGridHelpers::numCells(grid_),
+                                       Opm::UgGridHelpers::globalCell(grid_),
+                                       Opm::UgGridHelpers::cartDims(grid_),
+                                       Opm::UgGridHelpers::dimensions(grid_),
+                                       Opm::UgGridHelpers::cell2Faces(grid_),
+                                       Opm::UgGridHelpers::beginFaceCentroids(grid_),
+                                       props_.permeability(),
+                                       is_parallel_run_);
+            const Wells* wells = wells_manager.c_wells();
+            WellState well_state;
+            // well_state.init(wells, state, prev_well_state);
+
+            const std::vector<WellConstPtr>& wells_ecl = eclipse_state_->getSchedule()->getWells(timer.currentStepNum());
+            std::vector<WellMutliSegmentPtr> wells_multisegment(wells_ecl.size());
+            // wells_multisegment.resize(wells_ecl.size());
+            for (size_t i = 0; i < wells_multisegment.size(); ++i) {
+                wells_multisegment[i].reset(new WellMultiSegment(wells_ecl[i], timer.currentStepNum(), wells));
+            }
+            // for DEBUGGING OUTPUT
+            std::cout << " the number of the wells from EclipseState " << wells_ecl.size() << std::endl;
+            for (size_t i = 0; i < wells_ecl.size(); ++i) {
+                std::cout << " well name " << wells_ecl[i]->name() << std::endl;
+                std::cout << " segment wells " << wells_ecl[i]->isMultiSegment() << std::endl;
+            }
+            std::cin.ignore();
+
+            well_state.init(wells_multisegment, state, prev_well_state);
+
+            // give the polymer and surfactant simulators the chance to do their stuff
+            Base::asImpl().handleAdditionalWellInflow(timer, wells_manager, well_state, wells);
+
+            // write simulation state at the report stage
+            output_writer_.writeTimeStep( timer, state, well_state );
+
+            // Max oil saturation (for VPPARS), hysteresis update.
+            props_.updateSatOilMax(state.saturation());
+            props_.updateSatHyst(state.saturation(), allcells_);
+
+            // Compute reservoir volumes for RESV controls.
+            Base::asImpl().computeRESV(timer.currentStepNum(), wells, state, well_state);
+
+            // Run a multiple steps of the solver depending on the time step control.
+            solver_timer.start();
+
+            auto solver = Base::asImpl().createSolver(wells);
+
+            // If sub stepping is enabled allow the solver to sub cycle
+            // in case the report steps are too large for the solver to converge
+            //
+            // \Note: The report steps are met in any case
+            // \Note: The sub stepping will require a copy of the state variables
+            if( adaptiveTimeStepping ) {
+                adaptiveTimeStepping->step( timer, *solver, state, well_state,  output_writer_ );
+            }
+            else {
+                // solve for complete report step
+                solver->step(timer.currentStepLength(), state, well_state);
+            }
+
+            // take time that was used to solve system for this reportStep
+            solver_timer.stop();
+
+            // accumulate the number of Newton and Linear Iterations
+            totalNewtonIterations += solver->newtonIterations();
+            totalLinearIterations += solver->linearIterations();
+
+            // Report timing.
+            const double st = solver_timer.secsSinceStart();
+
+            if ( terminal_output_ )
+            {
+                std::cout << "Fully implicit solver took: " << st << " seconds." << std::endl;
+            }
+
+            stime += st;
+            if ( output_writer_.output() ) {
+                SimulatorReport step_report;
+                step_report.pressure_time = st;
+                step_report.total_time =  step_timer.secsSinceStart();
+                step_report.reportParam(tstep_os);
+            }
+
+            // Increment timer, remember well state.
+            ++timer;
+            prev_well_state = well_state;
+        }
+
+        // Write final simulation state.
+        output_writer_.writeTimeStep( timer, state, prev_well_state );
+
+        // Stop timer and create timing report
+        total_timer.stop();
+        SimulatorReport report;
+        report.pressure_time = stime;
+        report.transport_time = 0.0;
+        report.total_time = total_timer.secsSinceStart();
+        report.total_newton_iterations = totalNewtonIterations;
+        report.total_linear_iterations = totalLinearIterations;
+        return report;
+    }
+
+} // namespace Opm

opm/autodiff/WellStateMultiSegment.hpp

@@ -23,10 +23,10 @@
 
 #include <opm/core/wells.h>
 #include <opm/core/well_controls.h>
-#include <opm/core/simulator/WellState.hpp>
 #include <opm/core/utility/ErrorMacros.hpp>
 #include <opm/autodiff/AutoDiffBlock.hpp>
 #include <opm/autodiff/WellMultiSegment.hpp>
+#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
 #include <vector>
 #include <cassert>
 #include <string>
@@ -42,7 +42,7 @@ namespace Opm
     /// since we are avoiding to use the old wells structure
     /// it makes it might be a good idea not to relate this State to the WellState
     class WellStateMultiSegment
-         : public WellState
+         : public WellStateFullyImplicitBlackoil
     {
     public: