void restoreOPM_XWELKeyword(const std::string& restart_filename, int reportstep, bool unified, WellState& wellstate)
    {
        const char * keyword = "OPM_XWEL";
        const char* filename = restart_filename.c_str();
        ecl_file_type* file_type = ecl_file_open(filename, 0);

        if (file_type != NULL) {

            bool block_selected = unified ? ecl_file_select_rstblock_report_step(file_type , reportstep) : true;

            if (block_selected) {
                ecl_kw_type* xwel = ecl_file_iget_named_kw(file_type , keyword, 0);
                const double* xwel_data = ecl_kw_get_double_ptr(xwel);
                std::copy_n(xwel_data + wellstate.getRestartTemperatureOffset(), wellstate.temperature().size(), wellstate.temperature().begin());
                std::copy_n(xwel_data + wellstate.getRestartBhpOffset(), wellstate.bhp().size(), wellstate.bhp().begin());
                std::copy_n(xwel_data + wellstate.getRestartPerfPressOffset(), wellstate.perfPress().size(), wellstate.perfPress().begin());
                std::copy_n(xwel_data + wellstate.getRestartPerfRatesOffset(), wellstate.perfRates().size(), wellstate.perfRates().begin());
                std::copy_n(xwel_data + wellstate.getRestartWellRatesOffset(), wellstate.wellRates().size(), wellstate.wellRates().begin());
            } else {
                std::string error_str = "Restart file " +  restart_filename + " does not contain data for report step " + std::to_string(reportstep) + "!\n";
                throw std::runtime_error(error_str);
            }
            ecl_file_close(file_type);
        } else {
            std::string error_str = "Restart file " + restart_filename + " not found!\n";
            throw std::runtime_error(error_str);
        }
    }

    // Solve with no rock compressibility (linear eqn).
    void IncompTpfa::solveIncomp(const double dt,
                                 SimulatorState& state,
                                 WellState& well_state)
    {
        // Set up properties.
        computePerSolveDynamicData(dt, state, well_state);

        // Assemble.
        UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
        int ok = ifs_tpfa_assemble(gg, &forces_, &trans_[0], &gpress_omegaweighted_[0], h_);
        if (!ok) {
            OPM_THROW(std::runtime_error, "Failed assembling pressure system.");
        }

        // Solve.
        linsolver_.solve(h_->A, h_->b, h_->x);

        // Obtain solution.
        assert(int(state.pressure().size()) == grid_.number_of_cells);
        assert(int(state.faceflux().size()) == grid_.number_of_faces);
        ifs_tpfa_solution soln = { NULL, NULL, NULL, NULL };
        soln.cell_press = &state.pressure()[0];
        soln.face_flux  = &state.faceflux()[0];
        if (wells_ != NULL) {
            assert(int(well_state.bhp().size()) == wells_->number_of_wells);
            assert(int(well_state.perfRates().size()) == wells_->well_connpos[ wells_->number_of_wells ]);
            soln.well_flux = &well_state.perfRates()[0];
            soln.well_press = &well_state.bhp()[0];
        }
        ifs_tpfa_press_flux(gg, &forces_, &trans_[0], h_, &soln);
    }

    void SimulatorBase<Implementation>::computeWellPotentials(const Wells* wells,
                                                              const BlackoilState& x,
                                                              const WellState& xw,
                                                              std::vector<double>& well_potentials)
    {
        const int nw = wells->number_of_wells;
        const int np = wells->number_of_phases;

        well_potentials.clear();
        well_potentials.resize(nw*np,0.0);       
        for (int w = 0; w < nw; ++w) {
            for (int perf = wells->well_connpos[w]; perf < wells->well_connpos[w + 1]; ++perf) {
                const double well_cell_pressure = x.pressure()[wells->well_cells[perf]];
                const double drawdown_used = well_cell_pressure - xw.perfPress()[perf];
                const WellControls* ctrl = wells->ctrls[w];
                const int nwc = well_controls_get_num(ctrl);
                //Loop over all controls until we find a BHP control
                //that specifies what we need...
                double bhp = 0.0;
                for (int ctrl_index=0; ctrl_index < nwc; ++ctrl_index) {
                    if (well_controls_iget_type(ctrl, ctrl_index) == BHP) {
                        bhp = well_controls_iget_target(ctrl, ctrl_index);
                    }
                    // TODO: do something for thp;
                }
                // Calculate the pressure difference in the well perforation
                const double dp = xw.perfPress()[perf] - xw.bhp()[w];
                const double drawdown_maximum = well_cell_pressure - (bhp + dp);

                for (int phase = 0; phase < np; ++phase) {
                    well_potentials[w*np + phase] += (drawdown_maximum / drawdown_used * xw.perfPhaseRates()[perf*np + phase]);
                }
            }
        }
    }

    void
    BlackoilMultiSegmentModel<Grid>::
    prepareStep(const SimulatorTimerInterface& timer,
                const ReservoirState& reservoir_state,
                const WellState& well_state)
    {
        const double dt = timer.currentStepLength();
        pvdt_ = geo_.poreVolume() / dt;
        if (active_[Gas]) {
            updatePrimalVariableFromState(reservoir_state);
        }

        const int nw = wellsMultiSegment().size();

        if ( !msWellOps().has_multisegment_wells ) {
            wellModel().segVDt() = V::Zero(nw);
            return;
        }

        const int nseg_total = well_state.numSegments();
        std::vector<double> segment_volume;
        segment_volume.reserve(nseg_total);
        for (int w = 0; w < nw; ++w) {
            WellMultiSegmentConstPtr well = wellsMultiSegment()[w];
            const std::vector<double>& segment_volume_well = well->segmentVolume();
            segment_volume.insert(segment_volume.end(), segment_volume_well.begin(), segment_volume_well.end());
        }
        assert(int(segment_volume.size()) == nseg_total);
        wellModel().segVDt() = Eigen::Map<V>(segment_volume.data(), nseg_total) / dt;
    }

    void
    BlackoilOutputWriter::
    writeTimeStepSerial(const SimulatorTimerInterface& timer,
                        const SimulationDataContainer& state,
                        const WellState& wellState,
                        const data::Solution& simProps,
                        bool substep)
    {
        // Matlab output
        if( matlabWriter_ ) {
            matlabWriter_->writeTimeStep( timer, state, wellState, substep );
        }

        // ECL output
        if ( eclWriter_ )
        {
            const auto& initConfig = eclipseState_.getInitConfig();
            if (initConfig.restartRequested() && ((initConfig.getRestartStep()) == (timer.currentStepNum()))) {
                std::cout << "Skipping restart write in start of step " << timer.currentStepNum() << std::endl;
            } else {
                data::Solution combined_sol = simToSolution(state, phaseUsage_); // Get "normal" data (SWAT, PRESSURE, ...)
                combined_sol.insert(simProps.begin(), simProps.end());           // ... insert "extra" data (KR, VISC, ...)
                eclWriter_->writeTimeStep(timer.reportStepNum(),
                                          substep,
                                          timer.simulationTimeElapsed(),
                                          combined_sol,
                                          wellState.report(phaseUsage_));
            }
        }

        // write backup file
        if( backupfile_.is_open() )
        {
            int reportStep      = timer.reportStepNum();
            int currentTimeStep = timer.currentStepNum();
            if( (reportStep == currentTimeStep || // true for SimulatorTimer
                 currentTimeStep == 0 || // true for AdaptiveSimulatorTimer at reportStep
                 timer.done() ) // true for AdaptiveSimulatorTimer at reportStep
               && lastBackupReportStep_ != reportStep ) // only backup report step once
            {
                // store report step
                lastBackupReportStep_ = reportStep;
                // write resport step number
                backupfile_.write( (const char *) &reportStep, sizeof(int) );

                try {
                    backupfile_ << state;

                    const WellStateFullyImplicitBlackoil& boWellState = static_cast< const WellStateFullyImplicitBlackoil& > (wellState);
                    backupfile_ << boWellState;
                }
                catch ( const std::bad_cast& e )
                {
                }

                backupfile_ << std::flush;
            }
        } // end backup
    }

 /// Compute two-phase transport source terms from well terms.
 /// Note: Unlike the incompressible version of this function,
 ///       this version computes surface volume injection rates,
 ///       production rates are still total reservoir volumes.
 /// \param[in]  props         Fluid and rock properties.
 /// \param[in]  wells         Wells data structure.
 /// \param[in]  well_state    Well pressures and fluxes.
 /// \param[out] transport_src The transport source terms. They are to be interpreted depending on sign:
 ///                           (+) positive  inflow of first (water) phase (surface volume),
 ///                           (-) negative  total outflow of both phases (reservoir volume).
 void computeTransportSource(const BlackoilPropertiesInterface& props,
                             const Wells* wells,
                             const WellState& well_state,
                             std::vector<double>& transport_src)
 {
     int nc = props.numCells();
     transport_src.clear();
     transport_src.resize(nc, 0.0);
     // Well contributions.
     if (wells) {
         const int nw = wells->number_of_wells;
         const int np = wells->number_of_phases;
         if (np != 2) {
             OPM_THROW(std::runtime_error, "computeTransportSource() requires a 2 phase case.");
         }
         std::vector<double> A(np*np);
         for (int w = 0; w < nw; ++w) {
             const double* comp_frac = wells->comp_frac + np*w;
             for (int perf = wells->well_connpos[w]; perf < wells->well_connpos[w + 1]; ++perf) {
                 const int perf_cell = wells->well_cells[perf];
                 double perf_rate = well_state.perfRates()[perf];
                 if (perf_rate > 0.0) {
                     // perf_rate is a total inflow reservoir rate, we want a surface water rate.
                     if (wells->type[w] != INJECTOR) {
                         std::cout << "**** Warning: crossflow in well "
                                   << w << " perf " << perf - wells->well_connpos[w]
                                   << " ignored. Reservoir rate was "
                                   << perf_rate/Opm::unit::day << " m^3/day." << std::endl;
                         perf_rate = 0.0;
                     } else {
                         assert(std::fabs(comp_frac[0] + comp_frac[1] - 1.0) < 1e-6);
                         perf_rate *= comp_frac[0]; // Water reservoir volume rate.
                         props.matrix(1, &well_state.perfPress()[perf], comp_frac, &perf_cell, &A[0], 0);
                         perf_rate *= A[0];         // Water surface volume rate.
                     }
                 }
                 transport_src[perf_cell] += perf_rate;
             }
         }
     }
 }

        SimulatorReport
        assemble(const ReservoirState& reservoir_state,
                 WellState& well_state,
                 const bool initial_assembly)
        {
            SimulatorReport report;

            report += Base::assemble(reservoir_state, well_state, initial_assembly);
            if (initial_assembly) {
            }

            // Compute pressure residual.
            ADB pressure_residual = ADB::constant(V::Zero(residual_.material_balance_eq[0].size()));
            for (int phase = 0; phase < numPhases(); ++phase) {
                pressure_residual += residual_.material_balance_eq[phase] * scaling_[phase];
            }
            residual_.material_balance_eq[0] = pressure_residual; // HACK

            // Compute total reservoir volume flux.
            const int n = sd_.rq[0].mflux.size();
            V flux = V::Zero(n);
            for (int phase = 0; phase < numPhases(); ++phase) {
                UpwindSelector<double> upwind(grid_, ops_, sd_.rq[phase].dh.value());
                flux += sd_.rq[phase].mflux.value() / upwind.select(sd_.rq[phase].b.value());
            }

            // Storing the fluxes in the assemble() method is a bit of
            // a hack, but alternatives either require a more
            // significant redesign of the base class or are
            // inefficient.
            ReservoirState& s = const_cast<ReservoirState&>(reservoir_state);
            s.faceflux().resize(n);
            std::copy_n(flux.data(), n, s.faceflux().begin());
            if (asImpl().localWellsActive()) {
                const V& wflux = asImpl().wellModel().getStoredWellPerforationFluxes();
                assert(int(well_state.perfRates().size()) == wflux.size());
                std::copy_n(wflux.data(), wflux.size(), well_state.perfRates().begin());
            }

            return report;
        }

    /// Compute per-iteration dynamic properties.
    void IncompTpfa::computePerIterationDynamicData(const double /*dt*/,
                                                    const SimulatorState& state,
                                                    const WellState& well_state)
    {
        // These are the variables that get computed by this function:
        //
        // std::vector<double> porevol_
        // std::vector<double> rock_comp_
        // std::vector<double> pressures_

        computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol_);
        if (rock_comp_props_ && rock_comp_props_->isActive()) {
            for (int cell = 0; cell < grid_.number_of_cells; ++cell) {
                rock_comp_[cell] = rock_comp_props_->rockComp(state.pressure()[cell]);
            }
        }
        if (wells_) {
            std::copy(state.pressure().begin(), state.pressure().end(), pressures_.begin());
            std::copy(well_state.bhp().begin(), well_state.bhp().end(), pressures_.begin() + grid_.number_of_cells);
        }
    }

    /// Compute the output.
    void IncompTpfa::computeResults(SimulatorState& state,
                                    WellState& well_state) const
    {
        // Make sure h_ contains the direct-solution matrix
        // and right hand side (not jacobian and residual).
        // TODO: optimize by only adjusting b and diagonal of A.
        UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
        ifs_tpfa_assemble(gg, &forces_, &trans_[0], &gpress_omegaweighted_[0], h_);


        // Make sure h_->x contains the direct solution vector.
        assert(int(state.pressure().size()) == grid_.number_of_cells);
        assert(int(state.faceflux().size()) == grid_.number_of_faces);
        std::copy(state.pressure().begin(), state.pressure().end(), h_->x);
        std::copy(well_state.bhp().begin(), well_state.bhp().end(), h_->x + grid_.number_of_cells);

        // Obtain solution.
        ifs_tpfa_solution soln = { NULL, NULL, NULL, NULL };
        soln.cell_press = &state.pressure()[0];
        soln.face_flux  = &state.faceflux()[0];
        if (wells_ != NULL) {
            assert(int(well_state.bhp().size()) == wells_->number_of_wells);
            assert(int(well_state.perfRates().size()) == wells_->well_connpos[ wells_->number_of_wells ]);
            soln.well_flux = &well_state.perfRates()[0];
            soln.well_press = &well_state.bhp()[0];
        }
        ifs_tpfa_press_flux(gg, &forces_, &trans_[0], h_, &soln); // TODO: Check what parts of h_ are used here.
    }

 /// Compute the residual and Jacobian.
 void CompressibleTpfa::assemble(const double dt,
                                 const BlackoilState& state,
                                 const WellState& well_state)
 {
     const double* cell_press = &state.pressure()[0];
     const double* well_bhp = well_state.bhp().empty() ? NULL : &well_state.bhp()[0];
     const double* z = &state.surfacevol()[0];
     UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
     CompletionData completion_data;
     completion_data.wdp = ! wellperf_wdp_.empty() ? &wellperf_wdp_[0] : 0;
     completion_data.A = ! wellperf_A_.empty() ? &wellperf_A_[0] : 0;
     completion_data.phasemob = ! wellperf_phasemob_.empty() ? &wellperf_phasemob_[0] : 0;
     cfs_tpfa_res_wells wells_tmp;
     wells_tmp.W = const_cast<Wells*>(wells_);
     wells_tmp.data = &completion_data;
     cfs_tpfa_res_forces forces;
     forces.wells = &wells_tmp;
     forces.src = NULL; // Check if it is legal to leave it as NULL.
     compr_quantities_gen cq;
     cq.nphases = props_.numPhases();
     cq.Ac = &cell_A_[0];
     cq.dAc = &cell_dA_[0];
     cq.Af = &face_A_[0];
     cq.phasemobf = &face_phasemob_[0];
     cq.voldiscr = &cell_voldisc_[0];
     int was_adjusted = 0;
     if (! (rock_comp_props_ && rock_comp_props_->isActive())) {
         was_adjusted =
             cfs_tpfa_res_assemble(gg, dt, &forces, z, &cq, &trans_[0],
                                   &face_gravcap_[0], cell_press, well_bhp,
                                   &porevol_[0], h_);
     } else {
         was_adjusted =
             cfs_tpfa_res_comprock_assemble(gg, dt, &forces, z, &cq, &trans_[0],
                                            &face_gravcap_[0], cell_press, well_bhp,
                                            &porevol_[0], &initial_porevol_[0],
                                            &rock_comp_[0], h_);
     }
     singular_ = (was_adjusted == 1);
 }

    IterationReport
    BlackoilMultiSegmentModel<Grid>::solveWellEq(const std::vector<ADB>& mob_perfcells,
                                                 const std::vector<ADB>& b_perfcells,
                                                 SolutionState& state,
                                                 WellState& well_state)
    {
        IterationReport iter_report = Base::solveWellEq(mob_perfcells, b_perfcells, state, well_state);

        if (iter_report.converged) {
            // We must now update the state.segp and state.segqs members,
            // that the base version does not know about.
            const int np = numPhases();
            const int nseg_total =well_state.numSegments();
            {
                // We will set the segp primary variable to the new ones,
                // but we do not change the derivatives here.
                ADB::V new_segp = Eigen::Map<ADB::V>(well_state.segPress().data(), nseg_total);
                // Avoiding the copy below would require a value setter method
                // in AutoDiffBlock.
                std::vector<ADB::M> old_segp_derivs = state.segp.derivative();
                state.segp = ADB::function(std::move(new_segp), std::move(old_segp_derivs));
            }
            {
                // Need to reshuffle well rates, from phase running fastest
                // to wells running fastest.
                // The transpose() below switches the ordering.
                const DataBlock segrates = Eigen::Map<const DataBlock>(well_state.segPhaseRates().data(), nseg_total, np).transpose();
                ADB::V new_segqs = Eigen::Map<const V>(segrates.data(), nseg_total * np);
                std::vector<ADB::M> old_segqs_derivs = state.segqs.derivative();
                state.segqs = ADB::function(std::move(new_segqs), std::move(old_segqs_derivs));
            }

            // This is also called by the base version, but since we have updated
            // state.segp we must call it again.
            asImpl().computeWellConnectionPressures(state, well_state);
        }

        return iter_report;
    }

 /// Compute per-iteration dynamic properties for wells.
 void CompressibleTpfa::computeWellDynamicData(const double /*dt*/,
                                               const BlackoilState& /*state*/,
                                               const WellState& well_state)
 {
     // These are the variables that get computed by this function:
     //
     // std::vector<double> wellperf_A_;
     // std::vector<double> wellperf_phasemob_;
     const int np = props_.numPhases();
     const int nw = (wells_ != 0) ? wells_->number_of_wells : 0;
     const int nperf = (wells_ != 0) ? wells_->well_connpos[nw] : 0;
     wellperf_A_.resize(nperf*np*np);
     wellperf_phasemob_.resize(nperf*np);
     // The A matrix is set equal to the perforation grid cells'
     // matrix for producers, computed from bhp and injection
     // component fractions from
     // The mobilities are set equal to the perforation grid cells'
     // mobilities for producers.
     std::vector<double> mu(np);
     for (int w = 0; w < nw; ++w) {
         bool producer = (wells_->type[w] == PRODUCER);
         const double* comp_frac = &wells_->comp_frac[np*w];
         for (int j = wells_->well_connpos[w]; j < wells_->well_connpos[w+1]; ++j) {
             const int c = wells_->well_cells[j];
             double* wpA = &wellperf_A_[np*np*j];
             double* wpM = &wellperf_phasemob_[np*j];
             if (producer) {
                 const double* cA = &cell_A_[np*np*c];
                 std::copy(cA, cA + np*np, wpA);
                 const double* cM = &cell_phasemob_[np*c];
                 std::copy(cM, cM + np, wpM);
             } else {
                 const double bhp = well_state.bhp()[w];
                 double perf_p = bhp + wellperf_wdp_[j];
                 // Hack warning: comp_frac is used as a component
                 // surface-volume variable in calls to matrix() and
                 // viscosity(), but as a saturation in the call to
                 // relperm(). This is probably ok as long as injectors
                 // only inject pure fluids.
                 props_.matrix(1, &perf_p, comp_frac, &c, wpA, NULL);
                 props_.viscosity(1, &perf_p, comp_frac, &c, &mu[0], NULL);
                 assert(std::fabs(std::accumulate(comp_frac, comp_frac + np, 0.0) - 1.0) < 1e-6);
                 props_.relperm  (1, comp_frac, &c, wpM , NULL);
                 for (int phase = 0; phase < np; ++phase) {
                     wpM[phase] /= mu[phase];
                 }
             }
         }
     }
 }

 void SimulatorBase<Implementation>::computeWellPotentials(const Wells* wells,
                                                           const WellState& xw,
                                                           std::vector<double>& well_potentials)
 {
     const int nw = wells->number_of_wells;
     const int np = wells->number_of_phases;
     well_potentials.clear();
     well_potentials.resize(nw*np,0.0);
     for (int w = 0; w < nw; ++w) {
         for (int perf = wells->well_connpos[w]; perf < wells->well_connpos[w + 1]; ++perf) {
             for (int phase = 0; phase < np; ++phase) {
                 well_potentials[w*np + phase] += xw.wellPotentials()[perf*np + phase];
             }
         }
     }
 }

    /// Compute the output.
    void CompressibleTpfa::computeResults(BlackoilState& state,
                                          WellState& well_state) const
    {
        UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
        CompletionData completion_data;
        completion_data.wdp = ! wellperf_wdp_.empty() ? const_cast<double*>(&wellperf_wdp_[0]) : 0;
        completion_data.A = ! wellperf_A_.empty() ? const_cast<double*>(&wellperf_A_[0]) : 0;
        completion_data.phasemob = ! wellperf_phasemob_.empty() ? const_cast<double*>(&wellperf_phasemob_[0]) : 0;
        cfs_tpfa_res_wells wells_tmp;
        wells_tmp.W = const_cast<Wells*>(wells_);
        wells_tmp.data = &completion_data;
        cfs_tpfa_res_forces forces;
        forces.wells = &wells_tmp;
        forces.src = NULL;

        double* wpress = ! well_state.bhp      ().empty() ? & well_state.bhp      ()[0] : 0;
        double* wflux  = ! well_state.perfRates().empty() ? & well_state.perfRates()[0] : 0;

        cfs_tpfa_res_flux(gg,
                          &forces,
                          props_.numPhases(),
                          &trans_[0],
                          &cell_phasemob_[0],
                          &face_phasemob_[0],
                          &face_gravcap_[0],
                          &state.pressure()[0],
                          wpress,
                          &state.faceflux()[0],
                          wflux);
        cfs_tpfa_res_fpress(gg,
                            props_.numPhases(),
                            &htrans_[0],
                            &face_phasemob_[0],
                            &face_gravcap_[0],
                            h_,
                            &state.pressure()[0],
                            &state.faceflux()[0],
                            &state.facepressure()[0]);

        // Compute well perforation pressures (not done by the C code).
        if (wells_ != 0) {
            const int nw = wells_->number_of_wells;
            for (int w = 0; w < nw; ++w) {
                for (int j = wells_->well_connpos[w]; j < wells_->well_connpos[w+1]; ++j) {
                    const double bhp = well_state.bhp()[w];
                    well_state.perfPress()[j] = bhp + wellperf_wdp_[j];
                }
            }
        }
    }

    void SimulatorBase<Implementation>::computeRESV(const std::size_t               step,
                                                    const Wells*                    wells,
                                                    const BlackoilState&            x,
                                                    WellState& xw)
    {
        typedef SimFIBODetails::WellMap WellMap;

        const std::vector<WellConstPtr>& w_ecl = eclipse_state_->getSchedule()->getWells(step);
        const WellMap& wmap = SimFIBODetails::mapWells(w_ecl);

        const std::vector<int>& resv_wells = SimFIBODetails::resvWells(wells, step, wmap);

        const std::size_t number_resv_wells        = resv_wells.size();
        std::size_t       global_number_resv_wells = number_resv_wells;
#if HAVE_MPI
        if ( solver_.parallelInformation().type() == typeid(ParallelISTLInformation) )
        {
            const auto& info =
                boost::any_cast<const ParallelISTLInformation&>(solver_.parallelInformation());
            global_number_resv_wells = info.communicator().sum(global_number_resv_wells);
            if ( global_number_resv_wells )
            {
                // At least one process has resv wells. Therefore rate converter needs
                // to calculate averages over regions that might cross process
                // borders. This needs to be done by all processes and therefore
                // outside of the next if statement.
                rateConverter_.defineState(x, boost::any_cast<const ParallelISTLInformation&>(solver_.parallelInformation()));
            }
        }
        else
#endif
        {
            if ( global_number_resv_wells )
            {
                rateConverter_.defineState(x);
            }
        }

        if (! resv_wells.empty()) {
            const PhaseUsage&                    pu = props_.phaseUsage();
            const std::vector<double>::size_type np = props_.numPhases();

            std::vector<double> distr (np);
            std::vector<double> hrates(np);
            std::vector<double> prates(np);

            for (std::vector<int>::const_iterator
                     rp = resv_wells.begin(), e = resv_wells.end();
                 rp != e; ++rp)
            {
                WellControls* ctrl = wells->ctrls[*rp];
                const bool is_producer = wells->type[*rp] == PRODUCER;

                // RESV control mode, all wells
                {
                    const int rctrl = SimFIBODetails::resv_control(ctrl);

                    if (0 <= rctrl) {
                        const std::vector<double>::size_type off = (*rp) * np;

                        if (is_producer) {
                            // Convert to positive rates to avoid issues
                            // in coefficient calculations.
                            std::transform(xw.wellRates().begin() + (off + 0*np),
                                           xw.wellRates().begin() + (off + 1*np),
                                           prates.begin(), std::negate<double>());
                        } else {
                            std::copy(xw.wellRates().begin() + (off + 0*np),
                                      xw.wellRates().begin() + (off + 1*np),
                                      prates.begin());
                        }

                        const int fipreg = 0; // Hack.  Ignore FIP regions.
                        rateConverter_.calcCoeff(prates, fipreg, distr);

                        well_controls_iset_distr(ctrl, rctrl, & distr[0]);
                    }
                }

                // RESV control, WCONHIST wells.  A bit of duplicate
                // work, regrettably.
                if (is_producer && wells->name[*rp] != 0) {
                    WellMap::const_iterator i = wmap.find(wells->name[*rp]);

                    if (i != wmap.end()) {
                        WellConstPtr wp = i->second;

                        const WellProductionProperties& p =
                            wp->getProductionProperties(step);

                        if (! p.predictionMode) {
                            // History matching (WCONHIST/RESV)
                            SimFIBODetails::historyRates(pu, p, hrates);

                            const int fipreg = 0; // Hack.  Ignore FIP regions.
                            rateConverter_.calcCoeff(hrates, fipreg, distr);

                            // WCONHIST/RESV target is sum of all
                            // observed phase rates translated to
                            // reservoir conditions.  Recall sign
//.........这里部分代码省略.........

    void
    StandardWellsSolvent::
    computePropertiesForWellConnectionPressures(const SolutionState& state,
                                                const WellState& xw,
                                                std::vector<double>& b_perf,
                                                std::vector<double>& rsmax_perf,
                                                std::vector<double>& rvmax_perf,
                                                std::vector<double>& surf_dens_perf)
    {
        // 1. Compute properties required by computeConnectionPressureDelta().
        //    Note that some of the complexity of this part is due to the function
        //    taking std::vector<double> arguments, and not Eigen objects.
        const int nperf = wells().well_connpos[wells().number_of_wells];
        const int nw = wells().number_of_wells;

        // Compute the average pressure in each well block
        const Vector perf_press = Eigen::Map<const V>(xw.perfPress().data(), nperf);
        Vector avg_press = perf_press*0;
        for (int w = 0; w < nw; ++w) {
            for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) {
                const double p_above = perf == wells().well_connpos[w] ? state.bhp.value()[w] : perf_press[perf - 1];
                const double p_avg = (perf_press[perf] + p_above)/2;
                avg_press[perf] = p_avg;
            }
        }

        const std::vector<int>& well_cells = wellOps().well_cells;

        // Use cell values for the temperature as the wells don't knows its temperature yet.
        const ADB perf_temp = subset(state.temperature, well_cells);

        // Compute b, rsmax, rvmax values for perforations.
        // Evaluate the properties using average well block pressures
        // and cell values for rs, rv, phase condition and temperature.
        const ADB avg_press_ad = ADB::constant(avg_press);
        std::vector<PhasePresence> perf_cond(nperf);
        for (int perf = 0; perf < nperf; ++perf) {
            perf_cond[perf] = (*phase_condition_)[well_cells[perf]];
        }

        const PhaseUsage& pu = fluid_->phaseUsage();
        DataBlock b(nperf, pu.num_phases);

        const Vector bw = fluid_->bWat(avg_press_ad, perf_temp, well_cells).value();
        if (pu.phase_used[BlackoilPhases::Aqua]) {
            b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw;
        }

        assert((*active_)[Oil]);
        assert((*active_)[Gas]);
        const ADB perf_rv = subset(state.rv, well_cells);
        const ADB perf_rs = subset(state.rs, well_cells);
        const Vector perf_so =  subset(state.saturation[pu.phase_pos[Oil]].value(), well_cells);
        if (pu.phase_used[BlackoilPhases::Liquid]) {
            const Vector bo = fluid_->bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).value();
            //const V bo_eff = subset(rq_[pu.phase_pos[Oil] ].b , well_cells).value();
            b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo;
            // const Vector rssat = fluidRsSat(avg_press, perf_so, well_cells);
            const Vector rssat = fluid_->rsSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value();
            rsmax_perf.assign(rssat.data(), rssat.data() + nperf);
        } else {
            rsmax_perf.assign(0.0, nperf);
        }
        V surf_dens_copy = superset(fluid_->surfaceDensity(0, well_cells), Span(nperf, pu.num_phases, 0), nperf*pu.num_phases);
        for (int phase = 1; phase < pu.num_phases; ++phase) {
            if ( phase == pu.phase_pos[BlackoilPhases::Vapour]) {
                continue; // the gas surface density is added after the solvent is accounted for.
            }
            surf_dens_copy += superset(fluid_->surfaceDensity(phase, well_cells), Span(nperf, pu.num_phases, phase), nperf*pu.num_phases);
        }

        if (pu.phase_used[BlackoilPhases::Vapour]) {
            // Unclear wether the effective or the pure values should be used for the wells
            // the current usage of unmodified properties values gives best match.
            //V bg_eff = subset(rq_[pu.phase_pos[Gas]].b,well_cells).value();
            Vector bg = fluid_->bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value();
            Vector rhog = fluid_->surfaceDensity(pu.phase_pos[BlackoilPhases::Vapour], well_cells);
            // to handle solvent related
            if (has_solvent_) {

                const Vector bs = solvent_props_->bSolvent(avg_press_ad,well_cells).value();
                //const V bs_eff = subset(rq_[solvent_pos_].b,well_cells).value();

                // number of cells
                const int nc = state.pressure.size();

                const ADB zero = ADB::constant(Vector::Zero(nc));
                const ADB& ss = state.solvent_saturation;
                const ADB& sg = ((*active_)[ Gas ]
                                 ? state.saturation[ pu.phase_pos[ Gas ] ]
                                 : zero);

                Selector<double> zero_selector(ss.value() + sg.value(), Selector<double>::Zero);
                Vector F_solvent = subset(zero_selector.select(ss, ss / (ss + sg)),well_cells).value();

                Vector injectedSolventFraction = Eigen::Map<const Vector>(&xw.solventFraction()[0], nperf);

                Vector isProducer = Vector::Zero(nperf);
                Vector ones = Vector::Constant(nperf,1.0);
                for (int w = 0; w < nw; ++w) {
//.........这里部分代码省略.........

    /// Solve pressure equation, by Newton iterations.
    void CompressibleTpfa::solve(const double dt,
                                 BlackoilState& state,
                                 WellState& well_state)
    {
        const int nc = grid_.number_of_cells;
        const int nw = (wells_ != 0) ? wells_->number_of_wells : 0;

        // Set up dynamic data.
        computePerSolveDynamicData(dt, state, well_state);
        computePerIterationDynamicData(dt, state, well_state);

        // Assemble J and F.
        assemble(dt, state, well_state);

        double inc_norm = 0.0;
        int iter = 0;
        double res_norm = residualNorm();
        std::cout << "\nIteration         Residual        Change in p\n"
                  << std::setw(9) << iter
                  << std::setw(18) << res_norm
                  << std::setw(18) << '*' << std::endl;
        while ((iter < maxiter_) && (res_norm > residual_tol_)) {
            // Solve for increment in Newton method:
            //   incr = x_{n+1} - x_{n} = -J^{-1}F
            // (J is Jacobian matrix, F is residual)
            solveIncrement();
            ++iter;

            // Update pressure vars with increment.
            for (int c = 0; c < nc; ++c) {
                state.pressure()[c] += pressure_increment_[c];
            }
            for (int w = 0; w < nw; ++w) {
                well_state.bhp()[w] += pressure_increment_[nc + w];
            }

            // Stop iterating if increment is small.
            inc_norm = incrementNorm();
            if (inc_norm <= change_tol_) {
                std::cout << std::setw(9) << iter
                          << std::setw(18) << '*'
                          << std::setw(18) << inc_norm << std::endl;
                break;
            }

            // Set up dynamic data.
            computePerIterationDynamicData(dt, state, well_state);

            // Assemble J and F.
            assemble(dt, state, well_state);

            // Update residual norm.
            res_norm = residualNorm();

            std::cout << std::setw(9) << iter
                      << std::setw(18) << res_norm
                      << std::setw(18) << inc_norm << std::endl;
        }

        if ((iter == maxiter_) && (res_norm > residual_tol_) && (inc_norm > change_tol_)) {
            OPM_THROW(std::runtime_error, "CompressibleTpfa::solve() failed to converge in " << maxiter_ << " iterations.");
        }

        std::cout << "Solved pressure in " << iter << " iterations." << std::endl;

        // Compute fluxes and face pressures.
        computeResults(state, well_state);
    }