/*
  Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
  Copyright 2017 Statoil ASA.

  This file is part of the Open Porous Media project (OPM).

  OPM is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 3 of the License, or
  (at your option) any later version.

  OPM is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
*/


#include <opm/simulators/wells/MSWellHelpers.hpp>
#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/MSW/Valve.hpp>

namespace Opm
{


    template <typename TypeTag>
    MultisegmentWell<TypeTag>::
    MultisegmentWell(const Well& well, const int time_step,
                     const ModelParameters& param,
                     const RateConverterType& rate_converter,
                     const int pvtRegionIdx,
                     const int num_components,
                     const int num_phases,
                     const int index_of_well,
                     const int first_perf_index,
                     const std::vector<PerforationData>& perf_data)
        : Base(well, time_step, param, rate_converter, pvtRegionIdx, num_components, num_phases, index_of_well, first_perf_index, perf_data)
    , segment_perforations_(numberOfSegments())
    , segment_inlets_(numberOfSegments())
    , cell_perforation_depth_diffs_(number_of_perforations_, 0.0)
    , cell_perforation_pressure_diffs_(number_of_perforations_, 0.0)
    , perforation_segment_depth_diffs_(number_of_perforations_, 0.0)
    , segment_fluid_initial_(numberOfSegments(), std::vector<double>(num_components_, 0.0))
    , segment_densities_(numberOfSegments(), 0.0)
    , segment_viscosities_(numberOfSegments(), 0.0)
    , segment_mass_rates_(numberOfSegments(), 0.0)
    , segment_depth_diffs_(numberOfSegments(), 0.0)
    , upwinding_segments_(numberOfSegments(), 0)
    , segment_phase_fractions_(numberOfSegments(), std::vector<EvalWell>(num_components_, 0.0)) // number of phase here?
    , segment_phase_viscosities_(numberOfSegments(), std::vector<EvalWell>(num_components_, 0.0)) // number of phase here?
    {
        // not handling solvent or polymer for now with multisegment well
        if (has_solvent) {
            OPM_THROW(std::runtime_error, "solvent is not supported by multisegment well yet");
        }

        if (has_polymer) {
            OPM_THROW(std::runtime_error, "polymer is not supported by multisegment well yet");
        }

        if (Base::has_energy) {
            OPM_THROW(std::runtime_error, "energy is not supported by multisegment well yet");
        }

        if (Base::has_foam) {
            OPM_THROW(std::runtime_error, "foam is not supported by multisegment well yet");
        }

        if (Base::has_brine) {
            OPM_THROW(std::runtime_error, "brine is not supported by multisegment well yet");
        }
        // since we decide to use the WellSegments from the well parser. we can reuse a lot from it.
        // for other facilities needed but not available from parser, we need to process them here

        // initialize the segment_perforations_ and update perforation_segment_depth_diffs_
        const WellConnections& completion_set = well_ecl_.getConnections();
        // index of the perforation within wells struct
        // there might be some perforations not active, which causes the number of the perforations in
        // well_ecl_ and wells struct different
        // the current implementation is a temporary solution for now, it should be corrected from the parser
        // side
        int i_perf_wells = 0;
        perf_depth_.resize(number_of_perforations_, 0.);
        for (size_t perf = 0; perf < completion_set.size(); ++perf) {
            const Connection& connection = completion_set.get(perf);
            if (connection.state() == Connection::State::OPEN) {
                const int segment_index = segmentNumberToIndex(connection.segment());
                segment_perforations_[segment_index].push_back(i_perf_wells);
                perf_depth_[i_perf_wells] = connection.depth();
                const double segment_depth = segmentSet()[segment_index].depth();
                perforation_segment_depth_diffs_[i_perf_wells] = perf_depth_[i_perf_wells] - segment_depth;
                i_perf_wells++;
            }
        }

        // initialize the segment_inlets_
        for (int seg = 0; seg < numberOfSegments(); ++seg) {
            const Segment& segment = segmentSet()[seg];
            const int segment_number = segment.segmentNumber();
            const int outlet_segment_number = segment.outletSegment();
            if (outlet_segment_number > 0) {
                const int segment_index = segmentNumberToIndex(segment_number);
                const int outlet_segment_index = segmentNumberToIndex(outlet_segment_number);
                segment_inlets_[outlet_segment_index].push_back(segment_index);
            }
        }

        // calculating the depth difference between the segment and its oulet_segments
        // for the top segment, we will make its zero unless we find other purpose to use this value
        for (int seg = 1; seg < numberOfSegments(); ++seg) {
            const double segment_depth = segmentSet()[seg].depth();
            const int outlet_segment_number = segmentSet()[seg].outletSegment();
            const Segment& outlet_segment = segmentSet()[segmentNumberToIndex(outlet_segment_number)];
            const double outlet_depth = outlet_segment.depth();
            segment_depth_diffs_[seg] = segment_depth - outlet_depth;
        }
    }





    template <typename TypeTag>
    void
    MultisegmentWell<TypeTag>::
    init(const PhaseUsage* phase_usage_arg,
         const std::vector<double>& depth_arg,
         const double gravity_arg,
         const int num_cells)
    {
        Base::init(phase_usage_arg, depth_arg, gravity_arg, num_cells);

        // TODO: for StandardWell, we need to update the perf depth here using depth_arg.
        // for MultisegmentWell, it is much more complicated.
        // It can be specified directly, it can be calculated from the segment depth,
        // it can also use the cell center, which is the same for StandardWell.
        // For the last case, should we update the depth with the depth_arg? For the
        // future, it can be a source of wrong result with Multisegment well.
        // An indicator from the opm-parser should indicate what kind of depth we should use here.

        // \Note: we do not update the depth here. And it looks like for now, we only have the option to use
        // specified perforation depth
        initMatrixAndVectors(num_cells);

        // calcuate the depth difference between the perforations and the perforated grid block
        for (int perf = 0; perf < number_of_perforations_; ++perf) {
            const int cell_idx = well_cells_[perf];
            cell_perforation_depth_diffs_[perf] = depth_arg[cell_idx] - perf_depth_[perf];
        }
    }





    template <typename TypeTag>
    void
    MultisegmentWell<TypeTag>::
    initMatrixAndVectors(const int num_cells) const
    {
        duneB_.setBuildMode( OffDiagMatWell::row_wise );
        duneC_.setBuildMode( OffDiagMatWell::row_wise );
        duneD_.setBuildMode( DiagMatWell::row_wise );

        // set the size and patterns for all the matrices and vectors
        // [A C^T    [x    =  [ res
        //  B D] x_well]      res_well]

        // calculatiing the NNZ for duneD_
        // NNZ = number_of_segments + 2 * (number_of_inlets / number_of_outlets)
        {
            int nnz_d = numberOfSegments();
            for (const std::vector<int>& inlets : segment_inlets_) {
                nnz_d += 2 * inlets.size();
            }
            duneD_.setSize(numberOfSegments(), numberOfSegments(), nnz_d);
        }
        duneB_.setSize(numberOfSegments(), num_cells, number_of_perforations_);
        duneC_.setSize(numberOfSegments(), num_cells, number_of_perforations_);

        // we need to add the off diagonal ones
        for (auto row = duneD_.createbegin(), end = duneD_.createend(); row != end; ++row) {
            // the number of the row corrspnds to the segment now
            const int seg = row.index();
            // adding the item related to outlet relation
            const Segment& segment = segmentSet()[seg];
            const int outlet_segment_number = segment.outletSegment();
            if (outlet_segment_number > 0) { // if there is a outlet_segment
                const int outlet_segment_index = segmentNumberToIndex(outlet_segment_number);
                row.insert(outlet_segment_index);
            }

            // Add nonzeros for diagonal
            row.insert(seg);

            // insert the item related to its inlets
            for (const int& inlet : segment_inlets_[seg]) {
                row.insert(inlet);
            }
        }

        // make the C matrix
        for (auto row = duneC_.createbegin(), end = duneC_.createend(); row != end; ++row) {
            // the number of the row corresponds to the segment number now.
            for (const int& perf : segment_perforations_[row.index()]) {
                const int cell_idx = well_cells_[perf];
                row.insert(cell_idx);
            }
        }

        // make the B^T matrix
        for (auto row = duneB_.createbegin(), end = duneB_.createend(); row != end; ++row) {
            // the number of the row corresponds to the segment number now.
            for (const int& perf : segment_perforations_[row.index()]) {
                const int cell_idx = well_cells_[perf];
                row.insert(cell_idx);
            }
        }

        resWell_.resize( numberOfSegments() );

        primary_variables_.resize(numberOfSegments());
        primary_variables_evaluation_.resize(numberOfSegments());
    }





    template <typename TypeTag>
    void
    MultisegmentWell<TypeTag>::
    initPrimaryVariablesEvaluation() const
    {
        for (int seg = 0; seg < numberOfSegments(); ++seg) {
            for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) {
                primary_variables_evaluation_[seg][eq_idx] = 0.0;
                primary_variables_evaluation_[seg][eq_idx].setValue(primary_variables_[seg][eq_idx]);
                primary_variables_evaluation_[seg][eq_idx].setDerivative(eq_idx + numEq, 1.0);
            }
        }
    }





    template <typename TypeTag>
    void
    MultisegmentWell<TypeTag>::
    assembleWellEq(const Simulator& ebosSimulator,
                   const std::vector<Scalar>& B_avg,
                   const double dt,
                   WellState& well_state,
                   Opm::DeferredLogger& deferred_logger)
    {

        const bool use_inner_iterations = param_.use_inner_iterations_ms_wells_;
        if (use_inner_iterations) {
            this->iterateWellEquations(ebosSimulator, B_avg, dt, well_state, deferred_logger);
        }

        const auto& summary_state = ebosSimulator.vanguard().summaryState();
        const auto inj_controls = well_ecl_.isInjector() ? well_ecl_.injectionControls(summary_state) : Well::InjectionControls(0);
        const auto prod_controls = well_ecl_.isProducer() ? well_ecl_.productionControls(summary_state) : Well::ProductionControls(0);
        assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, deferred_logger);
    }





    template <typename TypeTag>
    void
    MultisegmentWell<TypeTag>::
    updateWellStateWithTarget(const Simulator& ebos_simulator,
                              WellState& well_state,
                              Opm::DeferredLogger&  deferred_logger) const
    {
        // segRates and  segPressure are used to initilize the primaryvariables for MSW wells
        // first initialize wellRates and then use it to compute segRates
        // When THP is supported for MSW wells this code and its fried in the standard model
        // can be merge.

        const auto& well = well_ecl_;
        const int well_index = index_of_well_;
        const int top_segment_index = well_state.topSegmentIndex(index_of_well_);
        const auto& pu = phaseUsage();
        const int np = well_state.numPhases();
        const auto& summaryState = ebos_simulator.vanguard().summaryState();

        if (wellIsStopped_) {
            for (int p = 0; p<np; ++p) {
                well_state.wellRates()[well_index*np + p] = 0.0;
            }
            return;
        }

        if (well.isInjector() )
        {
            const auto& controls = well.injectionControls(summaryState);

            InjectorType injectorType = controls.injector_type;
            int phasePos;
            switch (injectorType) {
            case InjectorType::WATER:
            {
                phasePos = pu.phase_pos[BlackoilPhases::Aqua];
                break;
            }
            case InjectorType::OIL:
            {
                phasePos = pu.phase_pos[BlackoilPhases::Liquid];
                break;
            }
            case InjectorType::GAS:
            {
                phasePos = pu.phase_pos[BlackoilPhases::Vapour];
                break;
            }
            default:
                throw("Expected WATER, OIL or GAS as type for injectors " + well.name());
            }

            const Opm::Well::InjectorCMode& current = well_state.currentInjectionControls()[well_index];

            switch(current) {
            case Well::InjectorCMode::RATE:
            {
                well_state.wellRates()[well_index*np + phasePos] = controls.surface_rate;
                break;
            }

            case Well::InjectorCMode::RESV:
            {
                std::vector<double> convert_coeff(number_of_phases_, 1.0);
                Base::rateConverter_.calcCoeff(/*fipreg*/ 0, Base::pvtRegionIdx_, convert_coeff);
                const double coeff = convert_coeff[phasePos];
                well_state.wellRates()[well_index*np + phasePos] = controls.reservoir_rate/coeff;
                break;
            }

            case Well::InjectorCMode::THP:
            {
                std::vector<double> rates(3, 0.0);
                for (int p = 0; p<np; ++p) {
                    rates[p] = well_state.wellRates()[well_index*np + p];
                }
                double bhp = calculateBhpFromThp(rates, well, summaryState, deferred_logger);
                well_state.bhp()[well_index] = bhp;
                break;
            }
            case Well::InjectorCMode::BHP:
            {
                well_state.segPress()[top_segment_index] = controls.bhp_limit;
                break;
            }
            case Well::InjectorCMode::GRUP:
            {
                //do nothing at the moment
                break;
            }
            case Well::InjectorCMode::CMODE_UNDEFINED:
            {
                OPM_DEFLOG_THROW(std::runtime_error, "Well control must be specified for well "  + name(), deferred_logger );
            }

            }
        }
        //Producer
        else
        {
            const Well::ProducerCMode& current = well_state.currentProductionControls()[well_index];
            const auto& controls = well.productionControls(summaryState);

            switch (current) {
            case Well::ProducerCMode::ORAT:
            {
                double current_rate = -well_state.wellRates()[ well_index*np + pu.phase_pos[Oil] ];

                if (current_rate == 0.0)
                    break;

                for (int p = 0; p<np; ++p) {
                    well_state.wellRates()[well_index*np + p] *= controls.oil_rate/current_rate;
                }
                break;
            }
            case Well::ProducerCMode::WRAT:
            {
                double current_rate = -well_state.wellRates()[ well_index*np + pu.phase_pos[Water] ];

                if (current_rate == 0.0)
                    break;

                for (int p = 0; p<np; ++p) {
                    well_state.wellRates()[well_index*np + p] *= controls.water_rate/current_rate;
                }
                break;
            }
            case Well::ProducerCMode::GRAT:
            {
                double current_rate = -well_state.wellRates()[ well_index*np + pu.phase_pos[Gas] ];

                if (current_rate == 0.0)
                    break;

                for (int p = 0; p<np; ++p) {
                    well_state.wellRates()[well_index*np + p] *= controls.gas_rate/current_rate;
                }
                break;

            }
            case Well::ProducerCMode::LRAT:
            {
                double current_rate = -well_state.wellRates()[ well_index*np + pu.phase_pos[Water] ]
                        - well_state.wellRates()[ well_index*np + pu.phase_pos[Oil] ];

                if (current_rate == 0.0)
                    break;

                for (int p = 0; p<np; ++p) {
                    well_state.wellRates()[well_index*np + p] *= controls.liquid_rate/current_rate;
                }
                break;
            }
            case Well::ProducerCMode::CRAT:
            {
                OPM_DEFLOG_THROW(std::runtime_error, "CRAT control not supported " << name(), deferred_logger);
            }
            case Well::ProducerCMode::RESV:
            {
                std::vector<double> convert_coeff(number_of_phases_, 1.0);
                Base::rateConverter_.calcCoeff(/*fipreg*/ 0, Base::pvtRegionIdx_, convert_coeff);
                double total_res_rate = 0.0;
                for (int p = 0; p<np; ++p) {
                    total_res_rate -= well_state.wellRates()[well_index*np + p] * convert_coeff[p];
                }
                if (total_res_rate == 0.0)
                    break;

                if (controls.prediction_mode) {
                    for (int p = 0; p<np; ++p) {
                        well_state.wellRates()[well_index*np + p] *= controls.resv_rate/total_res_rate;
                    }
                } else {
                    std::vector<double> hrates(number_of_phases_,0.);
                    if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
                        hrates[pu.phase_pos[Water]] = controls.water_rate;
                    }
                    if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
                        hrates[pu.phase_pos[Oil]] = controls.oil_rate;
                    }
                    if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
                        hrates[pu.phase_pos[Gas]] = controls.gas_rate;
                    }
                    std::vector<double> hrates_resv(number_of_phases_,0.);
                    Base::rateConverter_.calcReservoirVoidageRates(/*fipreg*/ 0, Base::pvtRegionIdx_, hrates, hrates_resv);
                    double target = std::accumulate(hrates_resv.begin(), hrates_resv.end(), 0.0);
                    for (int p = 0; p<np; ++p) {
                        well_state.wellRates()[well_index*np + p] *= target/total_res_rate;
                    }

                }
                break;
            }
            case Well::ProducerCMode::BHP:
            {
                well_state.segPress()[top_segment_index] = controls.bhp_limit;
                break;
            }
            case Well::ProducerCMode::THP:
            {
                std::vector<double> rates(3, 0.0);
                for (int p = 0; p<np; ++p) {
                    rates[p] = well_state.wellRates()[well_index*np + p];
                }
                double bhp = calculateBhpFromThp(rates, well, summaryState, deferred_logger);
                well_state.bhp()[well_index] = bhp;
                break;
            }
            case Well::ProducerCMode::GRUP:
            {
                //do nothing at the moment
                break;
            }
            case Well::ProducerCMode::CMODE_UNDEFINED:
            {
                OPM_DEFLOG_THROW(std::runtime_error, "Well control must be specified for well "  + name(), deferred_logger );
            }
            case Well::ProducerCMode::NONE:
            {
                OPM_DEFLOG_THROW(std::runtime_error, "Well control must be specified for well "  + name() , deferred_logger);
            }

            }
        }

        // compute the segment rates based on the wellRates
        initSegmentRatesWithWellRates(well_state);
    }





    template <typename TypeTag>
    void
    MultisegmentWell<TypeTag>::
    initSegmentRatesWithWellRates(WellState& well_state) const
    {
        for (int phase = 0; phase < number_of_phases_; ++phase) {
            const double perf_phaserate = well_state.wellRates()[number_of_phases_ * index_of_well_ + phase] / number_of_perforations_;
            for (int perf = 0; perf < number_of_perforations_; ++perf) {
                well_state.perfPhaseRates()[number_of_phases_ * (first_perf_ + perf) + phase] = perf_phaserate;
            }
        }

        const std::vector<double> perforation_rates(well_state.perfPhaseRates().begin() + number_of_phases_ * first_perf_,
                                                    well_state.perfPhaseRates().begin() +
                                                    number_of_phases_ * (first_perf_ + number_of_perforations_) );
        std::vector<double> segment_rates;
        WellState::calculateSegmentRates(segment_inlets_, segment_perforations_, perforation_rates, number_of_phases_,
                                         0, segment_rates);
        const int top_segment_index = well_state.topSegmentIndex(index_of_well_);
        std::copy(segment_rates.begin(), segment_rates.end(),
                  well_state.segRates().begin() + number_of_phases_ * top_segment_index );
        // we need to check the top segment rates should be same with the well rates
    }





    template <typename TypeTag>
    ConvergenceReport
    MultisegmentWell<TypeTag>::
    getWellConvergence(const WellState& well_state, const std::vector<double>& B_avg, Opm::DeferredLogger& deferred_logger, const bool relax_tolerance) const
    {
        assert(int(B_avg.size()) == num_components_);

        // checking if any residual is NaN or too large. The two large one is only handled for the well flux
        std::vector<std::vector<double>> abs_residual(numberOfSegments(), std::vector<double>(numWellEq, 0.0));
        for (int seg = 0; seg < numberOfSegments(); ++seg) {
            for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) {
                abs_residual[seg][eq_idx] = std::abs(resWell_[seg][eq_idx]);
            }
        }

        std::vector<double> maximum_residual(numWellEq, 0.0);

        ConvergenceReport report;
        // TODO: the following is a little complicated, maybe can be simplified in some way?
        for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) {
            for (int seg = 0; seg < numberOfSegments(); ++seg) {
                if (eq_idx < num_components_) { // phase or component mass equations
                    const double flux_residual = B_avg[eq_idx] * abs_residual[seg][eq_idx];
                    if (flux_residual > maximum_residual[eq_idx]) {
                        maximum_residual[eq_idx] = flux_residual;
                    }
                } else { // pressure or control equation
                    // for the top segment (seg == 0), it is control equation, will be checked later separately
                    if (seg > 0) {
                        // Pressure equation
                        const double pressure_residual = abs_residual[seg][eq_idx];
                        if (pressure_residual > maximum_residual[eq_idx]) {
                            maximum_residual[eq_idx] = pressure_residual;
                        }
                    }
                }
            }
        }

        using CR = ConvergenceReport;
        for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) {
            if (eq_idx < num_components_) { // phase or component mass equations
                const double flux_residual = maximum_residual[eq_idx];
                // TODO: the report can not handle the segment number yet.

                if (std::isnan(flux_residual)) {
                    report.setWellFailed({CR::WellFailure::Type::MassBalance, CR::Severity::NotANumber, eq_idx, name()});
                } else if (flux_residual > param_.max_residual_allowed_) {
                    report.setWellFailed({CR::WellFailure::Type::MassBalance, CR::Severity::TooLarge, eq_idx, name()});
                } else if (!relax_tolerance && flux_residual > param_.tolerance_wells_) {
                    report.setWellFailed({CR::WellFailure::Type::MassBalance, CR::Severity::Normal, eq_idx, name()});
                } else if (flux_residual > param_.relaxed_inner_tolerance_flow_ms_well_) {
                    report.setWellFailed({CR::WellFailure::Type::MassBalance, CR::Severity::Normal, eq_idx, name()});
                }
            } else { // pressure equation
                const double pressure_residual = maximum_residual[eq_idx];
                const int dummy_component = -1;
                if (std::isnan(pressure_residual)) {
                    report.setWellFailed({CR::WellFailure::Type::Pressure, CR::Severity::NotANumber, dummy_component, name()});
                } else if (std::isinf(pressure_residual)) {
                    report.setWellFailed({CR::WellFailure::Type::Pressure, CR::Severity::TooLarge, dummy_component, name()});
                } else if (!relax_tolerance && pressure_residual > param_.tolerance_pressure_ms_wells_) {
                    report.setWellFailed({CR::WellFailure::Type::Pressure, CR::Severity::Normal, dummy_component, name()});
                } else if (pressure_residual > param_.relaxed_inner_tolerance_pressure_ms_well_) {
                    report.setWellFailed({CR::WellFailure::Type::Pressure, CR::Severity::Normal, dummy_component, name()});
                }
            }
        }

        checkConvergenceControlEq(well_state, report, deferred_logger);

        return report;
    }





    template <typename TypeTag>
    void
    MultisegmentWell<TypeTag>::
    apply(const BVector& x, BVector& Ax) const
    {
        BVectorWell Bx(duneB_.N());

        duneB_.mv(x, Bx);

        // invDBx = duneD^-1 * Bx_
        const BVectorWell invDBx = mswellhelpers::applyUMFPack(duneD_, duneDSolver_, Bx);

        // Ax = Ax - duneC_^T * invDBx
        duneC_.mmtv(invDBx,Ax);
    }





    template <typename TypeTag>
    void
    MultisegmentWell<TypeTag>::
    apply(BVector& r) const
    {
        // invDrw_ = duneD^-1 * resWell_
        const BVectorWell invDrw = mswellhelpers::applyUMFPack(duneD_, duneDSolver_, resWell_);
        // r = r - duneC_^T * invDrw
        duneC_.mmtv(invDrw, r);
    }



#if HAVE_CUDA || HAVE_OPENCL
    template<typename TypeTag>
    void
    MultisegmentWell<TypeTag>::
    addWellContribution(WellContributions& wellContribs) const
    {
        unsigned int Nb = duneB_.M();       // number of blockrows in matrix A
        unsigned int Mb = duneB_.N();       // number of blockrows in duneB_, duneC_ and duneD_
        unsigned int BnumBlocks = duneB_.nonzeroes();
        unsigned int DnumBlocks = duneD_.nonzeroes();

        // duneC
        std::vector<unsigned int> Ccols;
        std::vector<double> Cvals;
        Ccols.reserve(BnumBlocks);
        Cvals.reserve(BnumBlocks * numEq * numWellEq);
        for (auto rowC = duneC_.begin(); rowC != duneC_.end(); ++rowC) {
            for (auto colC = rowC->begin(), endC = rowC->end(); colC != endC; ++colC) {
                Ccols.emplace_back(colC.index());
                for (int i = 0; i < numWellEq; ++i) {
                    for (int j = 0; j < numEq; ++j) {
                        Cvals.emplace_back((*colC)[i][j]);
                    }
                }
            }
        }

        // duneD
        Dune::UMFPack<DiagMatWell> umfpackMatrix(duneD_, 0);
        double *Dvals = umfpackMatrix.getInternalMatrix().getValues();
        int *Dcols = umfpackMatrix.getInternalMatrix().getColStart();
        int *Drows = umfpackMatrix.getInternalMatrix().getRowIndex();

        // duneB
        std::vector<unsigned int> Bcols;
        std::vector<unsigned int> Brows;
        std::vector<double> Bvals;
        Bcols.reserve(BnumBlocks);
        Brows.reserve(Mb+1);
        Bvals.reserve(BnumBlocks * numEq * numWellEq);
        Brows.emplace_back(0);
        unsigned int sumBlocks = 0;
        for (auto rowB = duneB_.begin(); rowB != duneB_.end(); ++rowB) {
            int sizeRow = 0;
            for (auto colB = rowB->begin(), endB = rowB->end(); colB != endB; ++colB) {
                Bcols.emplace_back(colB.index());
                for (int i = 0; i < numWellEq; ++i) {
                    for (int j = 0; j < numEq; ++j) {
                        Bvals.emplace_back((*colB)[i][j]);
                    }
                }
                sizeRow++;
            }
            sumBlocks += sizeRow;
            Brows.emplace_back(sumBlocks);
        }

        wellContribs.addMultisegmentWellContribution(numEq, numWellEq, Nb, Mb, BnumBlocks, Bvals, Bcols, Brows, DnumBlocks, Dvals, Dcols, Drows, Cvals);
    }
#endif


    template <typename TypeTag>
    void
    MultisegmentWell<TypeTag>::
    recoverWellSolutionAndUpdateWellState(const BVector& x,
                                          WellState& well_state,
                                          Opm::DeferredLogger& deferred_logger) const
    {
        BVectorWell xw(1);
        recoverSolutionWell(x, xw);
        updateWellState(xw, well_state, deferred_logger);
    }





    template <typename TypeTag>
    void
    MultisegmentWell<TypeTag>::
    computeWellPotentials(const Simulator& ebosSimulator,
                          const std::vector<Scalar>& B_avg,
                          const WellState& well_state,
                          std::vector<double>& well_potentials,
                          Opm::DeferredLogger& deferred_logger)
    {
        const int np = number_of_phases_;
        well_potentials.resize(np, 0.0);

        // Stopped wells have zero potential.
        if (this->wellIsStopped()) {
            return;
        }

        // If the well is pressure controlled the potential equals the rate.
/*         {
            bool pressure_controlled_well = false;
            if (this->isInjector()) {
                const Opm::Well::InjectorCMode& current = well_state.currentInjectionControls()[index_of_well_];
                if (current == Well::InjectorCMode::BHP || current == Well::InjectorCMode::THP) {
                    pressure_controlled_well = true;
                }
            } else {
                const Opm::Well::ProducerCMode& current = well_state.currentProductionControls()[index_of_well_];
                if (current == Well::ProducerCMode::BHP || current == Well::ProducerCMode::THP) {
                    pressure_controlled_well = true;
                }
            }
            if (pressure_controlled_well) {
                for (int compIdx = 0; compIdx < num_components_; ++compIdx) {
                    const EvalWell rate = this->getQs(compIdx);
                    well_potentials[ebosCompIdxToFlowCompIdx(compIdx)] = rate.value();
                }
                return;
            }
        } */

        // creating a copy of the well itself, to avoid messing up the explicit informations
        // during this copy, the only information not copied properly is the well controls
        MultisegmentWell<TypeTag> well(*this);
        well.debug_cost_counter_ = 0;

        well.updatePrimaryVariables(well_state, deferred_logger);

        // initialize the primary variables in Evaluation, which is used in computePerfRate for computeWellPotentials
        // TODO: for computeWellPotentials, no derivative is required actually
        well.initPrimaryVariablesEvaluation();

        // does the well have a THP related constraint?
        const auto& summaryState = ebosSimulator.vanguard().summaryState();
        const Well::ProducerCMode& current_control = well_state.currentProductionControls()[this->index_of_well_];
        if ( !well.Base::wellHasTHPConstraints(summaryState) || current_control == Well::ProducerCMode::BHP) {
            well.computeWellRatesAtBhpLimit(ebosSimulator, B_avg, well_potentials, deferred_logger);
        } else {
            well_potentials = well.computeWellPotentialWithTHP(ebosSimulator, B_avg, deferred_logger);
        }
        deferred_logger.debug("Cost in iterations of finding well potential for well "
                              + name() + ": " + std::to_string(well.debug_cost_counter_));
    }




    template<typename TypeTag>
    void
    MultisegmentWell<TypeTag>::
    computeWellRatesAtBhpLimit(const Simulator& ebosSimulator,
                               const std::vector<Scalar>& B_avg,
                               std::vector<double>& well_flux,
                               Opm::DeferredLogger& deferred_logger) const
    {
        if (well_ecl_.isInjector()) {
            const auto controls = well_ecl_.injectionControls(ebosSimulator.vanguard().summaryState());
            computeWellRatesWithBhp(ebosSimulator, B_avg, controls.bhp_limit, well_flux, deferred_logger);
        } else {
            const auto controls = well_ecl_.productionControls(ebosSimulator.vanguard().summaryState());
            computeWellRatesWithBhp(ebosSimulator, B_avg, controls.bhp_limit, well_flux, deferred_logger);
        }
    }




    template<typename TypeTag>
    void
    MultisegmentWell<TypeTag>::
    computeWellRatesWithBhp(const Simulator& ebosSimulator,
                            const std::vector<Scalar>& B_avg,
                            const Scalar bhp,
                            std::vector<double>& well_flux,
                            Opm::DeferredLogger& deferred_logger) const
    {
        // creating a copy of the well itself, to avoid messing up the explicit informations
        // during this copy, the only information not copied properly is the well controls
        MultisegmentWell<TypeTag> well_copy(*this);
        well_copy.debug_cost_counter_ = 0;

        // store a copy of the well state, we don't want to update the real well state
        WellState well_state_copy = ebosSimulator.problem().wellModel().wellState();

        // Get the current controls.
        const auto& summary_state = ebosSimulator.vanguard().summaryState();
        auto inj_controls = well_copy.well_ecl_.isInjector()
            ? well_copy.well_ecl_.injectionControls(summary_state)
            : Well::InjectionControls(0);
        auto prod_controls = well_copy.well_ecl_.isProducer()
            ? well_copy.well_ecl_.productionControls(summary_state) :
            Well::ProductionControls(0);

        //  Set current control to bhp, and bhp value in state, modify bhp limit in control object.
        if (well_copy.well_ecl_.isInjector()) {
            inj_controls.bhp_limit = bhp;
            well_state_copy.currentInjectionControls()[index_of_well_] = Well::InjectorCMode::BHP;
        } else {
            prod_controls.bhp_limit = bhp;
            well_state_copy.currentProductionControls()[index_of_well_] = Well::ProducerCMode::BHP;
        }
        well_state_copy.bhp()[well_copy.index_of_well_] = bhp;
        well_copy.calculateExplicitQuantities(ebosSimulator, well_state_copy, deferred_logger);
        const double dt = ebosSimulator.timeStepSize();
        // iterate to get a solution at the given bhp.
        well_copy.iterateWellEqWithControl(ebosSimulator, B_avg, dt, inj_controls, prod_controls, well_state_copy,
                                           deferred_logger);

        // compute the potential and store in the flux vector.
        well_flux.clear();
        const int np = number_of_phases_;
        well_flux.resize(np, 0.0);
        for (int compIdx = 0; compIdx < num_components_; ++compIdx) {
            const EvalWell rate = well_copy.getQs(compIdx);
            well_flux[ebosCompIdxToFlowCompIdx(compIdx)] = rate.value();
        }
        debug_cost_counter_ += well_copy.debug_cost_counter_;
    }



    template<typename TypeTag>
    std::vector<double>
    MultisegmentWell<TypeTag>::
    computeWellPotentialWithTHP(const Simulator& ebos_simulator,
                                const std::vector<Scalar>& B_avg,
                                Opm::DeferredLogger& deferred_logger) const
    {
        std::vector<double> potentials(number_of_phases_, 0.0);
        const auto& summary_state = ebos_simulator.vanguard().summaryState();

        const auto& well = well_ecl_;
        if (well.isInjector()){
            auto bhp_at_thp_limit = computeBhpAtThpLimitInj(ebos_simulator, B_avg, summary_state, deferred_logger);
            if (bhp_at_thp_limit) {
                const auto& controls = well_ecl_.injectionControls(summary_state);
                const double bhp = std::min(*bhp_at_thp_limit, controls.bhp_limit);
                computeWellRatesWithBhp(ebos_simulator, B_avg, bhp, potentials, deferred_logger);
                deferred_logger.debug("Converged thp based potential calculation for well "
                                      + name() + ", at bhp = " + std::to_string(bhp));
            } else {
                deferred_logger.warning("FAILURE_GETTING_CONVERGED_POTENTIAL",
                                        "Failed in getting converged thp based potential calculation for well "
                                        + name() + ". Instead the bhp based value is used");
                const auto& controls = well_ecl_.injectionControls(summary_state);
                const double bhp = controls.bhp_limit;
                computeWellRatesWithBhp(ebos_simulator, B_avg, bhp, potentials, deferred_logger);
            }
        } else {
            auto bhp_at_thp_limit = computeBhpAtThpLimitProd(ebos_simulator, B_avg, summary_state, deferred_logger);
            if (bhp_at_thp_limit) {
                const auto& controls = well_ecl_.productionControls(summary_state);
                const double bhp = std::max(*bhp_at_thp_limit, controls.bhp_limit);
                computeWellRatesWithBhp(ebos_simulator, B_avg, bhp, potentials, deferred_logger);
                deferred_logger.debug("Converged thp based potential calculation for well "
                                      + name() + ", at bhp = " + std::to_string(bhp));
            } else {
                deferred_logger.warning("FAILURE_GETTING_CONVERGED_POTENTIAL",
                                        "Failed in getting converged thp based potential calculation for well "
                                        + name() + ". Instead the bhp based value is used");
                const auto& controls = well_ecl_.productionControls(summary_state);
                const double bhp = controls.bhp_limit;
                computeWellRatesWithBhp(ebos_simulator, B_avg, bhp, potentials, deferred_logger);
            }
        }

        return potentials;
    }



    template <typename TypeTag>
    void
    MultisegmentWell<TypeTag>::
    updatePrimaryVariables(const WellState& well_state, Opm::DeferredLogger& /* deferred_logger */) const
    {
        // TODO: to test using rate conversion coefficients to see if it will be better than
        // this default one

        const Well& well = Base::wellEcl();

        // the index of the top segment in the WellState
        const int top_segment_index = well_state.topSegmentIndex(index_of_well_);
        const std::vector<double>& segment_rates = well_state.segRates();
        const PhaseUsage& pu = phaseUsage();

        for (int seg = 0; seg < numberOfSegments(); ++seg) {
            // calculate the total rate for each segment
            double total_seg_rate = 0.0;
            const int seg_index = top_segment_index + seg;
            // the segment pressure
            primary_variables_[seg][SPres] = well_state.segPress()[seg_index];
            // TODO: under what kind of circustances, the following will be wrong?
            // the definition of g makes the gas phase is always the last phase
            for (int p = 0; p < number_of_phases_; p++) {
                total_seg_rate += scalingFactor(p) * segment_rates[number_of_phases_ * seg_index + p];
            }

            primary_variables_[seg][GTotal] = total_seg_rate;
            if (std::abs(total_seg_rate) > 0.) {
                if (has_water) {
                    const int water_pos = pu.phase_pos[Water];
                    primary_variables_[seg][WFrac] = scalingFactor(water_pos) * segment_rates[number_of_phases_ * seg_index + water_pos] / total_seg_rate;
                }
                if (has_gas) {
                    const int gas_pos = pu.phase_pos[Gas];
                    primary_variables_[seg][GFrac] = scalingFactor(gas_pos) * segment_rates[number_of_phases_ * seg_index + gas_pos] / total_seg_rate;
                }
            } else { // total_seg_rate == 0
                if (this->isInjector()) {
                    // only single phase injection handled
                    auto phase = well.getInjectionProperties().injectorType;

                    if (has_water) {
                        if (phase == InjectorType::WATER) {
                            primary_variables_[seg][WFrac] = 1.0;
                        } else {
                            primary_variables_[seg][WFrac] = 0.0;
                        }
                    }

                    if (has_gas) {
                        if (phase == InjectorType::GAS) {
                            primary_variables_[seg][GFrac] = 1.0;
                        } else {
                            primary_variables_[seg][GFrac] = 0.0;
                        }
                    }

                } else if (this->isProducer()) { // producers
                    if (has_water) {
                        primary_variables_[seg][WFrac] = 1.0 / number_of_phases_;
                    }

                    if (has_gas) {
                        primary_variables_[seg][GFrac] = 1.0 / number_of_phases_;
                    }
                }
            }
        }
    }





    template <typename TypeTag>
    void
    MultisegmentWell<TypeTag>::
    recoverSolutionWell(const BVector& x, BVectorWell& xw) const
    {
        BVectorWell resWell = resWell_;
        // resWell = resWell - B * x
        duneB_.mmv(x, resWell);
        // xw = D^-1 * resWell
        xw = mswellhelpers::applyUMFPack(duneD_, duneDSolver_, resWell);
    }





    template <typename TypeTag>
    void
    MultisegmentWell<TypeTag>::
    solveEqAndUpdateWellState(WellState& well_state, Opm::DeferredLogger& deferred_logger)
    {
        // We assemble the well equations, then we check the convergence,
        // which is why we do not put the assembleWellEq here.
        const BVectorWell dx_well = mswellhelpers::applyUMFPack(duneD_, duneDSolver_, resWell_);

        updateWellState(dx_well, well_state, deferred_logger);
    }





    template <typename TypeTag>
    void
    MultisegmentWell<TypeTag>::
    computePerfCellPressDiffs(const Simulator& ebosSimulator)
    {
        for (int perf = 0; perf < number_of_perforations_; ++perf) {

            std::vector<double> kr(number_of_phases_, 0.0);
            std::vector<double> density(number_of_phases_, 0.0);

            const int cell_idx = well_cells_[perf];
            const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
            const auto& fs = intQuants.fluidState();

            double sum_kr = 0.;

            const PhaseUsage& pu = phaseUsage();
            if (pu.phase_used[Water]) {
                const int water_pos = pu.phase_pos[Water];
                kr[water_pos] = intQuants.relativePermeability(FluidSystem::waterPhaseIdx).value();
                sum_kr += kr[water_pos];
                density[water_pos] = fs.density(FluidSystem::waterPhaseIdx).value();
            }

            if (pu.phase_used[Oil]) {
                const int oil_pos = pu.phase_pos[Oil];
                kr[oil_pos] = intQuants.relativePermeability(FluidSystem::oilPhaseIdx).value();
                sum_kr += kr[oil_pos];
                density[oil_pos] = fs.density(FluidSystem::oilPhaseIdx).value();
            }

            if (pu.phase_used[Gas]) {
                const int gas_pos = pu.phase_pos[Gas];
                kr[gas_pos] = intQuants.relativePermeability(FluidSystem::gasPhaseIdx).value();
                sum_kr += kr[gas_pos];
                density[gas_pos] = fs.density(FluidSystem::gasPhaseIdx).value();
            }

            assert(sum_kr != 0.);

            // calculate the average density
            double average_density = 0.;
            for (int p = 0; p < number_of_phases_; ++p) {
                average_density += kr[p] * density[p];
            }
            average_density /= sum_kr;

            cell_perforation_pressure_diffs_[perf] = gravity_ * average_density * cell_perforation_depth_diffs_[perf];
        }
    }





    template <typename TypeTag>
    void
    MultisegmentWell<TypeTag>::
    computeInitialSegmentFluids(const Simulator& ebos_simulator)
    {
        for (int seg = 0; seg < numberOfSegments(); ++seg) {
            // TODO: trying to reduce the times for the surfaceVolumeFraction calculation
            const double surface_volume = getSegmentSurfaceVolume(ebos_simulator, seg).value();
            for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) {
                segment_fluid_initial_[seg][comp_idx] = surface_volume * surfaceVolumeFraction(seg, comp_idx).value();
            }
        }
    }





    template <typename TypeTag>
    void
    MultisegmentWell<TypeTag>::
    updateWellState(const BVectorWell& dwells,
                    WellState& well_state,
                    Opm::DeferredLogger& deferred_logger,
                    const double relaxation_factor) const
    {
        const double dFLimit = param_.dwell_fraction_max_;
        const double max_pressure_change = param_.max_pressure_change_ms_wells_;
        const std::vector<std::array<double, numWellEq> > old_primary_variables = primary_variables_;

        for (int seg = 0; seg < numberOfSegments(); ++seg) {
            if (has_water) {
                const int sign = dwells[seg][WFrac] > 0. ? 1 : -1;
                const double dx_limited = sign * std::min(std::abs(dwells[seg][WFrac]) * relaxation_factor, dFLimit);
                primary_variables_[seg][WFrac] = old_primary_variables[seg][WFrac] - dx_limited;
            }

            if (has_gas) {
                const int sign = dwells[seg][GFrac] > 0. ? 1 : -1;
                const double dx_limited = sign * std::min(std::abs(dwells[seg][GFrac]) * relaxation_factor, dFLimit);
                primary_variables_[seg][GFrac] = old_primary_variables[seg][GFrac] - dx_limited;
            }

            // handling the overshooting or undershooting of the fractions
            processFractions(seg);

            // update the segment pressure
            {
                const int sign = dwells[seg][SPres] > 0.? 1 : -1;
                const double dx_limited = sign * std::min(std::abs(dwells[seg][SPres]) * relaxation_factor, max_pressure_change);
                primary_variables_[seg][SPres] = std::max( old_primary_variables[seg][SPres] - dx_limited, 1e5);
            }

            // update the total rate // TODO: should we have a limitation of the total rate change?
            {
                primary_variables_[seg][GTotal] = old_primary_variables[seg][GTotal] - relaxation_factor * dwells[seg][GTotal];

                // make sure that no injector produce and no producer inject
                if (seg == 0) {
                    if (this->isInjector()) {
                        primary_variables_[seg][GTotal] = std::max( primary_variables_[seg][GTotal], 0.0);
                    } else {
                        primary_variables_[seg][GTotal] = std::min( primary_variables_[seg][GTotal], 0.0);
                    }
                }
            }

        }

        updateWellStateFromPrimaryVariables(well_state, deferred_logger);
        Base::calculateReservoirRates(well_state);
    }





    template <typename TypeTag>
    void
    MultisegmentWell<TypeTag>::
    calculateExplicitQuantities(const Simulator& ebosSimulator,
                                const WellState& well_state,
                                Opm::DeferredLogger& deferred_logger)
    {
        updatePrimaryVariables(well_state, deferred_logger);
        initPrimaryVariablesEvaluation();
        computePerfCellPressDiffs(ebosSimulator);
        computeInitialSegmentFluids(ebosSimulator);
    }





    template<typename TypeTag>
    void
    MultisegmentWell<TypeTag>::
    addWellContributions(SparseMatrixAdapter& /* jacobian */) const
    {
        OPM_THROW(std::runtime_error, "addWellContributions is not supported by multisegment well yet");
    }





    template <typename TypeTag>
    const WellSegments&
    MultisegmentWell<TypeTag>::
    segmentSet() const
    {
        return well_ecl_.getSegments();
    }





    template <typename TypeTag>
    int
    MultisegmentWell<TypeTag>::
    numberOfSegments() const
    {
        return segmentSet().size();
    }





    template <typename TypeTag>
    int
    MultisegmentWell<TypeTag>::
    numberOfPerforations() const
    {
        return segmentSet().number_of_perforations_;
    }





    template <typename TypeTag>
    WellSegments::CompPressureDrop
    MultisegmentWell<TypeTag>::
    compPressureDrop() const
    {
        return segmentSet().compPressureDrop();
    }







    template <typename TypeTag>
    int
    MultisegmentWell<TypeTag>::
    segmentNumberToIndex(const int segment_number) const
    {
        return segmentSet().segmentNumberToIndex(segment_number);
    }





    template <typename TypeTag>
    typename MultisegmentWell<TypeTag>::EvalWell
    MultisegmentWell<TypeTag>::
    volumeFraction(const int seg, const unsigned compIdx) const
    {

        if (has_water && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx)) {
            return primary_variables_evaluation_[seg][WFrac];
        }

        if (has_gas && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx)) {
            return primary_variables_evaluation_[seg][GFrac];
        }

        // Oil fraction
        EvalWell oil_fraction = 1.0;
        if (has_water) {
            oil_fraction -= primary_variables_evaluation_[seg][WFrac];
        }

        if (has_gas) {
            oil_fraction -= primary_variables_evaluation_[seg][GFrac];
        }
        /* if (has_solvent) {
            oil_fraction -= primary_variables_evaluation_[seg][SFrac];
        } */
        return oil_fraction;
    }




    template <typename TypeTag>
    typename MultisegmentWell<TypeTag>::EvalWell
    MultisegmentWell<TypeTag>::
    volumeFractionScaled(const int seg, const int comp_idx) const
    {
        // For reservoir rate control, the distr in well control is used for the
        // rate conversion coefficients. For the injection well, only the distr of the injection
        // phase is not zero.
        const double scale = scalingFactor(ebosCompIdxToFlowCompIdx(comp_idx));
        if (scale > 0.) {
            return volumeFraction(seg, comp_idx) / scale;
        }

        return volumeFraction(seg, comp_idx);
    }





    template <typename TypeTag>
    typename MultisegmentWell<TypeTag>::EvalWell
    MultisegmentWell<TypeTag>::
    surfaceVolumeFraction(const int seg, const int comp_idx) const
    {
        EvalWell sum_volume_fraction_scaled = 0.;
        for (int idx = 0; idx < num_components_; ++idx) {
            sum_volume_fraction_scaled += volumeFractionScaled(seg, idx);
        }

        assert(sum_volume_fraction_scaled.value() != 0.);

        return volumeFractionScaled(seg, comp_idx) / sum_volume_fraction_scaled;
    }





    template <typename TypeTag>
    void
    MultisegmentWell<TypeTag>::
    computePerfRatePressure(const IntensiveQuantities& int_quants,
                            const std::vector<EvalWell>& mob_perfcells,
                            const int seg,
                            const int perf,
                            const EvalWell& segment_pressure,
                            const bool& allow_cf,
                            std::vector<EvalWell>& cq_s,
                            EvalWell& perf_press,
                            double& perf_dis_gas_rate,
                            double& perf_vap_oil_rate,
                            Opm::DeferredLogger& deferred_logger) const

    {
        std::vector<EvalWell> cmix_s(num_components_, 0.0);

        // the composition of the components inside wellbore
        for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) {
            cmix_s[comp_idx] = surfaceVolumeFraction(seg, comp_idx);
        }

        const auto& fs = int_quants.fluidState();

        const EvalWell pressure_cell = extendEval(fs.pressure(FluidSystem::oilPhaseIdx));
        const EvalWell rs = extendEval(fs.Rs());
        const EvalWell rv = extendEval(fs.Rv());

        // not using number_of_phases_ because of solvent
        std::vector<EvalWell> b_perfcells(num_components_, 0.0);

        for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
            if (!FluidSystem::phaseIsActive(phaseIdx)) {
                continue;
            }

            const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
            b_perfcells[compIdx] = extendEval(fs.invB(phaseIdx));
        }

        // pressure difference between the segment and the perforation
        const EvalWell perf_seg_press_diff = gravity_ * segment_densities_[seg] * perforation_segment_depth_diffs_[perf];
        // pressure difference between the perforation and the grid cell
        const double cell_perf_press_diff = cell_perforation_pressure_diffs_[perf];

        perf_press = pressure_cell - cell_perf_press_diff;
        // Pressure drawdown (also used to determine direction of flow)
        // TODO: not 100% sure about the sign of the seg_perf_press_diff
        const EvalWell drawdown = perf_press - (segment_pressure + perf_seg_press_diff);

        // producing perforations
        if ( drawdown > 0.0) {
            // Do nothing is crossflow is not allowed
            if (!allow_cf && this->isInjector()) {
                return;
            }

            // compute component volumetric rates at standard conditions
            for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) {
                const EvalWell cq_p = - well_index_[perf] * (mob_perfcells[comp_idx] * drawdown);
                cq_s[comp_idx] = b_perfcells[comp_idx] * cq_p;
            }

            if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
                const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
                const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
                const EvalWell cq_s_oil = cq_s[oilCompIdx];
                const EvalWell cq_s_gas = cq_s[gasCompIdx];
                cq_s[gasCompIdx] += rs * cq_s_oil;
                cq_s[oilCompIdx] += rv * cq_s_gas;
            }
        } else { // injecting perforations
            // Do nothing if crossflow is not allowed
            if (!allow_cf && this->isProducer()) {
                return;
            }

            // for injecting perforations, we use total mobility
            EvalWell total_mob = mob_perfcells[0];
            for (int comp_idx = 1; comp_idx < num_components_; ++comp_idx) {
                total_mob += mob_perfcells[comp_idx];
            }

            // injection perforations total volume rates
            const EvalWell cqt_i = - well_index_[perf] * (total_mob * drawdown);

            // compute volume ratio between connection and at standard conditions
            EvalWell volume_ratio = 0.0;
            if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
                const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
                volume_ratio += cmix_s[waterCompIdx] / b_perfcells[waterCompIdx];
            }

            if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
                const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
                const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);

                // Incorporate RS/RV factors if both oil and gas active
                // TODO: not sure we use rs rv from the perforation cells when handling injecting perforations
                // basically, for injecting perforations, the wellbore is the upstreaming side.
                const EvalWell d = 1.0 - rv * rs;

                if (d.value() == 0.0) {
                    OPM_DEFLOG_THROW(Opm::NumericalIssue, "Zero d value obtained for well " << name() << " during flux calcuation"
                                                  << " with rs " << rs << " and rv " << rv, deferred_logger);
                }

                const EvalWell tmp_oil = (cmix_s[oilCompIdx] - rv * cmix_s[gasCompIdx]) / d;
                volume_ratio += tmp_oil / b_perfcells[oilCompIdx];

                const EvalWell tmp_gas = (cmix_s[gasCompIdx] - rs * cmix_s[oilCompIdx]) / d;
                volume_ratio += tmp_gas / b_perfcells[gasCompIdx];
            } else { // not having gas and oil at the same time
                if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
                    const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
                    volume_ratio += cmix_s[oilCompIdx] / b_perfcells[oilCompIdx];
                }
                if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
                    const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
                    volume_ratio += cmix_s[gasCompIdx] / b_perfcells[gasCompIdx];
                }
            }
            // injecting connections total volumerates at standard conditions
            EvalWell cqt_is = cqt_i / volume_ratio;
            for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) {
                cq_s[comp_idx] = cmix_s[comp_idx] * cqt_is;
            }
        } // end for injection perforations

        // calculating the perforation solution gas rate and solution oil rates
        if (this->isProducer()) {
            if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
                const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
                const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
                // TODO: the formulations here remain to be tested with cases with strong crossflow through production wells
                // s means standard condition, r means reservoir condition
                // q_os = q_or * b_o + rv * q_gr * b_g
                // q_gs = q_gr * g_g + rs * q_or * b_o
                // d = 1.0 - rs * rv
                // q_or = 1 / (b_o * d) * (q_os - rv * q_gs)
                // q_gr = 1 / (b_g * d) * (q_gs - rs * q_os)

                const double d = 1.0 - rv.value() * rs.value();
                // vaporized oil into gas
                // rv * q_gr * b_g = rv * (q_gs - rs * q_os) / d
                perf_vap_oil_rate = rv.value() * (cq_s[gasCompIdx].value() - rs.value() * cq_s[oilCompIdx].value()) / d;
                // dissolved of gas in oil
                // rs * q_or * b_o = rs * (q_os - rv * q_gs) / d
                perf_dis_gas_rate = rs.value() * (cq_s[oilCompIdx].value() - rv.value() * cq_s[gasCompIdx].value()) / d;
            }
        }
    }





    template <typename TypeTag>
    typename MultisegmentWell<TypeTag>::EvalWell
    MultisegmentWell<TypeTag>::
    extendEval(const Eval& in) const
    {
        EvalWell out = 0.0;
        out.setValue(in.value());
        for(int eq_idx = 0; eq_idx < numEq;++eq_idx) {
            out.setDerivative(eq_idx, in.derivative(eq_idx));
        }
        return out;
    }





    template <typename TypeTag>
    void
    MultisegmentWell<TypeTag>::
    computeSegmentFluidProperties(const Simulator& ebosSimulator)
    {
        // TODO: the concept of phases and components are rather confusing in this function.
        // needs to be addressed sooner or later.

        // get the temperature for later use. It is only useful when we are not handling
        // thermal related simulation
        // basically, it is a single value for all the segments

        EvalWell temperature;
        EvalWell saltConcentration;
        // not sure how to handle the pvt region related to segment
        // for the current approach, we use the pvt region of the first perforated cell
        // although there are some text indicating using the pvt region of the lowest
        // perforated cell
        // TODO: later to investigate how to handle the pvt region
        int pvt_region_index;
        {
            // using the first perforated cell
            const int cell_idx = well_cells_[0];
            const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
            const auto& fs = intQuants.fluidState();
            temperature.setValue(fs.temperature(FluidSystem::oilPhaseIdx).value());
            saltConcentration = extendEval(fs.saltConcentration());
            pvt_region_index = fs.pvtRegionIndex();
        }

        std::vector<double> surf_dens(num_components_);
        // Surface density.
        for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
            if (!FluidSystem::phaseIsActive(phaseIdx)) {
                continue;
            }

            const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
            surf_dens[compIdx] = FluidSystem::referenceDensity( phaseIdx, pvt_region_index );
        }

        for (int seg = 0; seg < numberOfSegments(); ++seg) {
            // the compostion of the components inside wellbore under surface condition
            std::vector<EvalWell> mix_s(num_components_, 0.0);
            for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) {
                mix_s[comp_idx] = surfaceVolumeFraction(seg, comp_idx);
            }

            std::vector<EvalWell> b(num_components_, 0.0);
            std::vector<EvalWell> visc(num_components_, 0.0);

            const EvalWell seg_pressure = getSegmentPressure(seg);
            if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
                const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
                b[waterCompIdx] =
                    FluidSystem::waterPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, saltConcentration);
                visc[waterCompIdx] =
                    FluidSystem::waterPvt().viscosity(pvt_region_index, temperature, seg_pressure, saltConcentration);
            }

            EvalWell rv(0.0);
            // gas phase
            if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
                const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
                if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
                    const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
                    const EvalWell rvmax = FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvt_region_index, temperature, seg_pressure);
                    if (mix_s[oilCompIdx] > 0.0) {
                        if (mix_s[gasCompIdx] > 0.0) {
                            rv = mix_s[oilCompIdx] / mix_s[gasCompIdx];
                        }

                        if (rv > rvmax) {
                            rv = rvmax;
                        }
                        b[gasCompIdx] =
                            FluidSystem::gasPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, rv);
                        visc[gasCompIdx] =
                            FluidSystem::gasPvt().viscosity(pvt_region_index, temperature, seg_pressure, rv);
                    } else { // no oil exists
                        b[gasCompIdx] =
                            FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure);
                        visc[gasCompIdx] =
                            FluidSystem::gasPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure);
                    }
                } else { // no Liquid phase
                    // it is the same with zero mix_s[Oil]
                    b[gasCompIdx] =
                        FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure);
                    visc[gasCompIdx] =
                        FluidSystem::gasPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure);
                }
            }

            EvalWell rs(0.0);
            // oil phase
            if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
                const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
                if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
                    const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
                    const EvalWell rsmax = FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvt_region_index, temperature, seg_pressure);
                    if (mix_s[gasCompIdx] > 0.0) {
                        if (mix_s[oilCompIdx] > 0.0) {
                            rs = mix_s[gasCompIdx] / mix_s[oilCompIdx];
                        }

                        if (rs > rsmax) {
                            rs = rsmax;
                        }
                        b[oilCompIdx] =
                            FluidSystem::oilPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, rs);
                        visc[oilCompIdx] =
                            FluidSystem::oilPvt().viscosity(pvt_region_index, temperature, seg_pressure, rs);
                    } else { // no oil exists
                        b[oilCompIdx] =
                            FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure);
                        visc[oilCompIdx] =
                            FluidSystem::oilPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure);
                    }
                } else { // no Liquid phase
                    // it is the same with zero mix_s[Oil]
                    b[oilCompIdx] =
                        FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure);
                    visc[oilCompIdx] =
                        FluidSystem::oilPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure);
                }
            }

            segment_phase_viscosities_[seg] = visc;

            std::vector<EvalWell> mix(mix_s);
            if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
                const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
                const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);

                const EvalWell d = 1.0 - rs * rv;

                if (rs != 0.0) { // rs > 0.0?
                    mix[gasCompIdx] = (mix_s[gasCompIdx] - mix_s[oilCompIdx] * rs) / d;
                }
                if (rv != 0.0) { // rv > 0.0?
                    mix[oilCompIdx] = (mix_s[oilCompIdx] - mix_s[gasCompIdx] * rv) / d;
                }
            }

            EvalWell volrat(0.0);
            for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) {
                volrat += mix[comp_idx] / b[comp_idx];
            }

            segment_viscosities_[seg] = 0.;
            // calculate the average viscosity
            for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) {
                const EvalWell fraction =  mix[comp_idx] / b[comp_idx] / volrat;
                // TODO: a little more work needs to be done to handle the negative fractions here
                segment_phase_fractions_[seg][comp_idx] = fraction; // >= 0.0 ? fraction : 0.0;
                segment_viscosities_[seg] += visc[comp_idx] * segment_phase_fractions_[seg][comp_idx];
            }

            EvalWell density(0.0);
            for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) {
                density += surf_dens[comp_idx] * mix_s[comp_idx];
            }
            segment_densities_[seg] = density / volrat;

            // calculate the mass rates
            segment_mass_rates_[seg] = 0.;
            for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) {
                const EvalWell rate = getSegmentRateUpwinding(seg, comp_idx);
                segment_mass_rates_[seg] += rate * surf_dens[comp_idx];
            }
        }
    }





    template <typename TypeTag>
    typename MultisegmentWell<TypeTag>::EvalWell
    MultisegmentWell<TypeTag>::
    getSegmentPressure(const int seg) const
    {
        return primary_variables_evaluation_[seg][SPres];
    }





    template <typename TypeTag>
    typename MultisegmentWell<TypeTag>::EvalWell
    MultisegmentWell<TypeTag>::
    getBhp() const
    {
        return getSegmentPressure(0);
    }





    template <typename TypeTag>
    typename MultisegmentWell<TypeTag>::EvalWell
    MultisegmentWell<TypeTag>::
    getSegmentRate(const int seg,
                   const int comp_idx) const
    {
        return primary_variables_evaluation_[seg][GTotal] * volumeFractionScaled(seg, comp_idx);
    }





    template <typename TypeTag>
    typename MultisegmentWell<TypeTag>::EvalWell
    MultisegmentWell<TypeTag>::
    getQs(const int comp_idx) const
    {
        return getSegmentRate(0, comp_idx);
    }





    template <typename TypeTag>
    typename MultisegmentWell<TypeTag>::EvalWell
    MultisegmentWell<TypeTag>::
    getSegmentRateUpwinding(const int seg,
                            const size_t comp_idx) const
    {
        const int seg_upwind = upwinding_segments_[seg];
        // the result will contain the derivative with resepct to GTotal in segment seg,
        // and the derivatives with respect to WFrac GFrac in segment seg_upwind.
        // the derivative with respect to SPres should be zero.
        if (seg == 0 && this->isInjector()) {
            const Well& well = Base::wellEcl();
            auto phase = well.getInjectionProperties().injectorType;

            if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)
                    && Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx) == comp_idx
                    && phase == InjectorType::WATER)
                return primary_variables_evaluation_[seg][GTotal] / scalingFactor(ebosCompIdxToFlowCompIdx(comp_idx));


            if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)
                    && Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx) == comp_idx
                    && phase == InjectorType::OIL)
                return primary_variables_evaluation_[seg][GTotal] / scalingFactor(ebosCompIdxToFlowCompIdx(comp_idx));

            if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)
                    && Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx) == comp_idx
                    && phase == InjectorType::GAS)
                return primary_variables_evaluation_[seg][GTotal] / scalingFactor(ebosCompIdxToFlowCompIdx(comp_idx));

            return 0.0;
        }

        const EvalWell segment_rate = primary_variables_evaluation_[seg][GTotal] * volumeFractionScaled(seg_upwind, comp_idx);

        assert(segment_rate.derivative(SPres + numEq) == 0.);

        return segment_rate;
    }





    template <typename TypeTag>
    typename MultisegmentWell<TypeTag>::EvalWell
    MultisegmentWell<TypeTag>::
    getSegmentGTotal(const int seg) const
    {
        return primary_variables_evaluation_[seg][GTotal];
    }





    template <typename TypeTag>
    typename MultisegmentWell<TypeTag>::EvalWell
    MultisegmentWell<TypeTag>::
    getWQTotal() const
    {
        return getSegmentGTotal(0);
    }





    template <typename TypeTag>
    void
    MultisegmentWell<TypeTag>::
    getMobility(const Simulator& ebosSimulator,
                const int perf,
                std::vector<EvalWell>& mob) const
    {
        // TODO: most of this function, if not the whole function, can be moved to the base class
        const int cell_idx = well_cells_[perf];
        assert (int(mob.size()) == num_components_);
        const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
        const auto& materialLawManager = ebosSimulator.problem().materialLawManager();

        // either use mobility of the perforation cell or calcualte its own
        // based on passing the saturation table index
        const int satid = saturation_table_number_[perf] - 1;
        const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx);
        if( satid == satid_elem ) { // the same saturation number is used. i.e. just use the mobilty from the cell

            for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
                if (!FluidSystem::phaseIsActive(phaseIdx)) {
                    continue;
                }

                const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
                mob[activeCompIdx] = extendEval(intQuants.mobility(phaseIdx));
            }
            // if (has_solvent) {
            //     mob[contiSolventEqIdx] = extendEval(intQuants.solventMobility());
            // }
        } else {

            const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx);
            Eval relativePerms[3] = { 0.0, 0.0, 0.0 };
            MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState());

            // reset the satnumvalue back to original
            materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx);

            // compute the mobility
            for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
                if (!FluidSystem::phaseIsActive(phaseIdx)) {
                    continue;
                }

                const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
                mob[activeCompIdx] = extendEval(relativePerms[phaseIdx] / intQuants.fluidState().viscosity(phaseIdx));
            }

        }
    }




    template <typename TypeTag>
    void
    MultisegmentWell<TypeTag>::
    assembleControlEq(const WellState& well_state,
                      const Opm::Schedule& schedule,
                      const SummaryState& summaryState,
                      const Well::InjectionControls& inj_controls,
                      const Well::ProductionControls& prod_controls,
                      Opm::DeferredLogger& deferred_logger)
    {

        EvalWell control_eq(0.0);

        const auto& well = well_ecl_;

        auto getRates = [&]() {
            std::vector<EvalWell> rates(3, 0.0);
            if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
                rates[Water] = getQs(Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx));
            }
            if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
                rates[Oil] = getQs(Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx));
            }
            if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
                rates[Gas] = getQs(Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx));
            }
            return rates;
        };

        if (wellIsStopped_) {
            control_eq = getWQTotal();
        } else if (this->isInjector() ) {
            // Find scaling factor to get injection rate,
            const InjectorType injectorType = inj_controls.injector_type;
            double scaling = 1.0;
            const auto& pu = phaseUsage();
            switch (injectorType) {
            case InjectorType::WATER:
            {
                scaling = scalingFactor(pu.phase_pos[BlackoilPhases::Aqua]);
                break;
            }
            case InjectorType::OIL:
            {
                scaling = scalingFactor(pu.phase_pos[BlackoilPhases::Liquid]);
                break;
            }
            case InjectorType::GAS:
            {
                scaling = scalingFactor(pu.phase_pos[BlackoilPhases::Vapour]);
                break;
            }
            default:
                throw("Expected WATER, OIL or GAS as type for injectors " + well.name());
            }
            const EvalWell injection_rate = getWQTotal() / scaling;
            // Setup function for evaluation of BHP from THP (used only if needed).
            auto bhp_from_thp = [&]() {
                const auto rates = getRates();
                return calculateBhpFromThp(rates, well, summaryState, deferred_logger);
            };
            // Call generic implementation.
            Base::assembleControlEqInj(well_state, schedule, summaryState, inj_controls, getBhp(), injection_rate, bhp_from_thp, control_eq, deferred_logger);
        } else {
            // Find rates.
            const auto rates = getRates();
            // Setup function for evaluation of BHP from THP (used only if needed).
            auto bhp_from_thp = [&]() {
                return calculateBhpFromThp(rates, well, summaryState, deferred_logger);
            };
            // Call generic implementation.
            Base::assembleControlEqProd(well_state, schedule, summaryState, prod_controls, getBhp(), rates, bhp_from_thp, control_eq, deferred_logger);
        }

        // using control_eq to update the matrix and residuals
        resWell_[0][SPres] = control_eq.value();
        for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
            duneD_[0][0][SPres][pv_idx] = control_eq.derivative(pv_idx + numEq);
        }
    }


    template<typename TypeTag>
    void
    MultisegmentWell<TypeTag>::
    updateThp(WellState& well_state, Opm::DeferredLogger& deferred_logger) const
    {
        // When there is no vaild VFP table provided, we set the thp to be zero.
        if (!this->isVFPActive(deferred_logger) || this->wellIsStopped()) {
            well_state.thp()[index_of_well_] = 0.;
            return;
        }

        // the well is under other control types, we calculate the thp based on bhp and rates
        std::vector<double> rates(3, 0.0);

        const Opm::PhaseUsage& pu = phaseUsage();
        if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
            rates[ Water ] = well_state.wellRates()[index_of_well_ * number_of_phases_ + pu.phase_pos[ Water ] ];
        }
        if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
            rates[ Oil ] = well_state.wellRates()[index_of_well_ * number_of_phases_ + pu.phase_pos[ Oil ] ];
        }
        if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
            rates[ Gas ] = well_state.wellRates()[index_of_well_ * number_of_phases_ + pu.phase_pos[ Gas ] ];
        }

        const double bhp = well_state.bhp()[index_of_well_];

        well_state.thp()[index_of_well_] = calculateThpFromBhp(rates, bhp, deferred_logger);

    }



    template<typename TypeTag>
    double
    MultisegmentWell<TypeTag>::
    calculateThpFromBhp(const std::vector<double>& rates,
                        const double bhp,
                        Opm::DeferredLogger& deferred_logger) const
    {
        assert(int(rates.size()) == 3); // the vfp related only supports three phases now.

        const double aqua = rates[Water];
        const double liquid = rates[Oil];
        const double vapour = rates[Gas];

        // pick the density in the top segment
        const double rho = segment_densities_[0].value();

        double thp = 0.0;
        if (this->isInjector()) {
            const int table_id = well_ecl_.vfp_table_number();
            const double vfp_ref_depth = vfp_properties_->getInj()->getTable(table_id)->getDatumDepth();
            const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_);

            thp = vfp_properties_->getInj()->thp(table_id, aqua, liquid, vapour, bhp + dp);
        }
        else if (this->isProducer()) {
            const int table_id = well_ecl_.vfp_table_number();
            const double alq = well_ecl_.alq_value();
            const double vfp_ref_depth = vfp_properties_->getProd()->getTable(table_id)->getDatumDepth();
            const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_);

            thp = vfp_properties_->getProd()->thp(table_id, aqua, liquid, vapour, bhp + dp, alq);
        }
        else {
            OPM_DEFLOG_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well", deferred_logger);
        }

        return thp;
    }


    template<typename TypeTag>
    template<class ValueType>
    ValueType
    MultisegmentWell<TypeTag>::
    calculateBhpFromThp(const std::vector<ValueType>& rates,
                        const Well& well,
                        const SummaryState& summaryState,
                        Opm::DeferredLogger& deferred_logger) const
    {
        // TODO: when well is under THP control, the BHP is dependent on the rates,
        // the well rates is also dependent on the BHP, so it might need to do some iteration.
        // However, when group control is involved, change of the rates might impacts other wells
        // so iterations on a higher level will be required. Some investigation might be needed when
        // we face problems under THP control.

        assert(int(rates.size()) == 3); // the vfp related only supports three phases now.

        const ValueType aqua = rates[Water];
        const ValueType liquid = rates[Oil];
        const ValueType vapour = rates[Gas];

        // pick the density in the top layer
        // TODO: it is possible it should be a Evaluation
        const double rho = segment_densities_[0].value();

        if (well.isInjector() )
        {
            const auto& controls = well.injectionControls(summaryState);
            const double vfp_ref_depth = vfp_properties_->getInj()->getTable(controls.vfp_table_number)->getDatumDepth();
            const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_);
            return vfp_properties_->getInj()->bhp(controls.vfp_table_number, aqua, liquid, vapour, controls.thp_limit) - dp;
         }
         else if (well.isProducer()) {
             const auto& controls = well.productionControls(summaryState);
             const double vfp_ref_depth = vfp_properties_->getProd()->getTable(controls.vfp_table_number)->getDatumDepth();
             const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_);
             return vfp_properties_->getProd()->bhp(controls.vfp_table_number, aqua, liquid, vapour, controls.thp_limit, controls.alq_value) - dp;
         }
         else {
             OPM_DEFLOG_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well", deferred_logger);
         }
    }



    template <typename TypeTag>
    void
    MultisegmentWell<TypeTag>::
    assemblePressureEq(const int seg, WellState& well_state) const
    {
        assert(seg != 0); // not top segment

        // for top segment, the well control equation will be used.
        EvalWell pressure_equation = getSegmentPressure(seg);

        // we need to handle the pressure difference between the two segments
        // we only consider the hydrostatic pressure loss first
        // TODO: we might be able to add member variables to store these values, then we update well state
        // after converged
        const auto hydro_pressure_drop = getHydroPressureLoss(seg);
        well_state.segPressDropHydroStatic()[seg] = hydro_pressure_drop.value();
        pressure_equation -= hydro_pressure_drop;

        if (frictionalPressureLossConsidered()) {
            const auto friction_pressure_drop = getFrictionPressureLoss(seg);
            pressure_equation -= friction_pressure_drop;
            well_state.segPressDropFriction()[seg] = friction_pressure_drop.value();
        }

        resWell_[seg][SPres] = pressure_equation.value();
        const int seg_upwind = upwinding_segments_[seg];
        duneD_[seg][seg][SPres][SPres] += pressure_equation.derivative(SPres + numEq);
        duneD_[seg][seg][SPres][GTotal] += pressure_equation.derivative(GTotal + numEq);
        if (has_water) {
            duneD_[seg][seg_upwind][SPres][WFrac] += pressure_equation.derivative(WFrac + numEq);
        }
        if (has_gas) {
            duneD_[seg][seg_upwind][SPres][GFrac] += pressure_equation.derivative(GFrac + numEq);
        }

        // contribution from the outlet segment
        const int outlet_segment_index = segmentNumberToIndex(segmentSet()[seg].outletSegment());
        const EvalWell outlet_pressure = getSegmentPressure(outlet_segment_index);

        resWell_[seg][SPres] -= outlet_pressure.value();
        for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
            duneD_[seg][outlet_segment_index][SPres][pv_idx] = -outlet_pressure.derivative(pv_idx + numEq);
        }

        if (accelerationalPressureLossConsidered()) {
            handleAccelerationPressureLoss(seg, well_state);
        }
    }





    template <typename TypeTag>
    typename MultisegmentWell<TypeTag>::EvalWell
    MultisegmentWell<TypeTag>::
    getHydroPressureLoss(const int seg) const
    {
        return segment_densities_[seg] * gravity_ * segment_depth_diffs_[seg];
    }





    template <typename TypeTag>
    typename MultisegmentWell<TypeTag>::EvalWell
    MultisegmentWell<TypeTag>::
    getFrictionPressureLoss(const int seg) const
    {
        const EvalWell mass_rate = segment_mass_rates_[seg];
        const int seg_upwind = upwinding_segments_[seg];
        EvalWell density = segment_densities_[seg_upwind];
        EvalWell visc = segment_viscosities_[seg_upwind];
        // WARNING
        // We disregard the derivatives from the upwind density to make sure derivatives
        // wrt. to different segments dont get mixed.
        if (seg != seg_upwind) {
            density.clearDerivatives();
            visc.clearDerivatives();
        }
        const int outlet_segment_index = segmentNumberToIndex(segmentSet()[seg].outletSegment());
        const double length = segmentSet()[seg].totalLength() - segmentSet()[outlet_segment_index].totalLength();
        assert(length > 0.);
        const double roughness = segmentSet()[seg].roughness();
        const double area = segmentSet()[seg].crossArea();
        const double diameter = segmentSet()[seg].internalDiameter();

        const double sign = mass_rate < 0. ? 1.0 : - 1.0;

        return sign * mswellhelpers::frictionPressureLoss(length, diameter, area, roughness, density, mass_rate, visc);
    }





    template <typename TypeTag>
    void
    MultisegmentWell<TypeTag>::
    handleAccelerationPressureLoss(const int seg, WellState& well_state) const
    {
        const double area = segmentSet()[seg].crossArea();
        const EvalWell mass_rate = segment_mass_rates_[seg];
        const int seg_upwind = upwinding_segments_[seg];
        EvalWell density = segment_densities_[seg_upwind];
        // WARNING
        // We disregard the derivatives from the upwind density to make sure derivatives
        // wrt. to different segments dont get mixed.
        if (seg != seg_upwind) {
            density.clearDerivatives();
        }

        EvalWell accelerationPressureLoss = mswellhelpers::velocityHead(area, mass_rate, density);
        // handling the velocity head of intlet segments
        for (const int inlet : segment_inlets_[seg]) {
            const int seg_upwind_inlet = upwinding_segments_[inlet];
            const double inlet_area = segmentSet()[inlet].crossArea();
            EvalWell inlet_density = segment_densities_[seg_upwind_inlet];
            // WARNING
            // We disregard the derivatives from the upwind density to make sure derivatives
            // wrt. to different segments dont get mixed.
            if (inlet != seg_upwind_inlet) {
                inlet_density.clearDerivatives();
            }
            const EvalWell inlet_mass_rate = segment_mass_rates_[inlet];
            accelerationPressureLoss -= mswellhelpers::velocityHead(std::max(inlet_area, area), inlet_mass_rate, inlet_density);
        }
        
        // We change the sign of the accelerationPressureLoss for injectors. 
        // Is this correct? Testing indicates that this is what the reference simulator does
        const double sign = mass_rate < 0. ? 1.0 : - 1.0;
        accelerationPressureLoss *= sign;

        well_state.segPressDropAcceleration()[seg] = accelerationPressureLoss.value();

        resWell_[seg][SPres] -= accelerationPressureLoss.value();
        duneD_[seg][seg][SPres][SPres] -= accelerationPressureLoss.derivative(SPres + numEq);
        duneD_[seg][seg][SPres][GTotal] -= accelerationPressureLoss.derivative(GTotal + numEq);
        if (has_water) {
            duneD_[seg][seg_upwind][SPres][WFrac] -= accelerationPressureLoss.derivative(WFrac + numEq);
        }
        if (has_gas) {
            duneD_[seg][seg_upwind][SPres][GFrac] -= accelerationPressureLoss.derivative(GFrac + numEq);
        }
    }





    template <typename TypeTag>
    void
    MultisegmentWell<TypeTag>::
    processFractions(const int seg) const
    {
        const PhaseUsage& pu = phaseUsage();

        std::vector<double> fractions(number_of_phases_, 0.0);

        assert( FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) );
        const int oil_pos = pu.phase_pos[Oil];
        fractions[oil_pos] = 1.0;

        if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) {
            const int water_pos = pu.phase_pos[Water];
            fractions[water_pos] = primary_variables_[seg][WFrac];
            fractions[oil_pos] -= fractions[water_pos];
        }

        if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) {
            const int gas_pos = pu.phase_pos[Gas];
            fractions[gas_pos] = primary_variables_[seg][GFrac];
            fractions[oil_pos] -= fractions[gas_pos];
        }

        if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
            const int water_pos = pu.phase_pos[Water];
            if (fractions[water_pos] < 0.0) {
                if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) {
                    fractions[pu.phase_pos[Gas]] /= (1.0 - fractions[water_pos]);
                }
                fractions[oil_pos] /= (1.0 - fractions[water_pos]);
                fractions[water_pos] = 0.0;
            }
        }

        if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
            const int gas_pos = pu.phase_pos[Gas];
            if (fractions[gas_pos] < 0.0) {
                if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) {
                    fractions[pu.phase_pos[Water]] /= (1.0 - fractions[gas_pos]);
                }
                fractions[oil_pos] /= (1.0 - fractions[gas_pos]);
                fractions[gas_pos] = 0.0;
            }
        }

        if (fractions[oil_pos] < 0.0) {
            if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) {
                fractions[pu.phase_pos[Water]] /= (1.0 - fractions[oil_pos]);
            }
            if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) {
                fractions[pu.phase_pos[Gas]] /= (1.0 - fractions[oil_pos]);
            }
            fractions[oil_pos] = 0.0;
        }

        if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) {
            primary_variables_[seg][WFrac] = fractions[pu.phase_pos[Water]];
        }

        if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) {
            primary_variables_[seg][GFrac] = fractions[pu.phase_pos[Gas]];
        }
    }





    template <typename TypeTag>
    void
    MultisegmentWell<TypeTag>::
    checkWellOperability(const Simulator& /* ebos_simulator */,
                         const WellState& /* well_state */,
                         Opm::DeferredLogger& deferred_logger)
    {
        const bool checkOperability = EWOMS_GET_PARAM(TypeTag, bool, EnableWellOperabilityCheck);
        if (!checkOperability) {
            return;
        }

        // focusing on PRODUCER for now
        if (this->isInjector()) {
            return;
        }

        if (!this->underPredictionMode() ) {
            return;
        }

        const std::string msg = "Support of well operability checking for multisegment wells is not implemented "
                                "yet, checkWellOperability() for " + name() + " will do nothing";
        deferred_logger.warning("NO_OPERATABILITY_CHECKING_MS_WELLS", msg);
    }





    template <typename TypeTag>
    void
    MultisegmentWell<TypeTag>::
    updateWellStateFromPrimaryVariables(WellState& well_state, Opm::DeferredLogger& deferred_logger) const
    {
        const PhaseUsage& pu = phaseUsage();
        assert( FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) );
        const int oil_pos = pu.phase_pos[Oil];

        for (int seg = 0; seg < numberOfSegments(); ++seg) {
            std::vector<double> fractions(number_of_phases_, 0.0);
            fractions[oil_pos] = 1.0;

            if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) {
                const int water_pos = pu.phase_pos[Water];
                fractions[water_pos] = primary_variables_[seg][WFrac];
                fractions[oil_pos] -= fractions[water_pos];
            }

            if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) {
                const int gas_pos = pu.phase_pos[Gas];
                fractions[gas_pos] = primary_variables_[seg][GFrac];
                fractions[oil_pos] -= fractions[gas_pos];
            }

            // convert the fractions to be Q_p / G_total to calculate the phase rates
            for (int p = 0; p < number_of_phases_; ++p) {
                const double scale = scalingFactor(p);
                // for injection wells, there should only one non-zero scaling factor
                if (scale > 0.) {
                    fractions[p] /= scale;
                } else {
                    // this should only happens to injection wells
                    fractions[p] = 0.;
                }
            }

            // calculate the phase rates based on the primary variables
            const double g_total = primary_variables_[seg][GTotal];
            const int top_segment_index = well_state.topSegmentIndex(index_of_well_);
            for (int p = 0; p < number_of_phases_; ++p) {
                const double phase_rate = g_total * fractions[p];
                well_state.segRates()[(seg + top_segment_index) * number_of_phases_ + p] = phase_rate;
                if (seg == 0) { // top segment
                    well_state.wellRates()[index_of_well_ * number_of_phases_ + p] = phase_rate;
                }
            }

            // update the segment pressure
            well_state.segPress()[seg + top_segment_index] = primary_variables_[seg][SPres];
            if (seg == 0) { // top segment
                well_state.bhp()[index_of_well_] = well_state.segPress()[seg + top_segment_index];
            }
        }
        updateThp(well_state, deferred_logger);
    }




    template <typename TypeTag>
    bool
    MultisegmentWell<TypeTag>::
    frictionalPressureLossConsidered() const
    {
        // HF- and HFA needs to consider frictional pressure loss
        return (segmentSet().compPressureDrop() != WellSegments::CompPressureDrop::H__);
    }





    template <typename TypeTag>
    bool
    MultisegmentWell<TypeTag>::
    accelerationalPressureLossConsidered() const
    {
        return (segmentSet().compPressureDrop() == WellSegments::CompPressureDrop::HFA);
    }





    template<typename TypeTag>
    bool
    MultisegmentWell<TypeTag>::
    iterateWellEqWithControl(const Simulator& ebosSimulator,
                             const std::vector<Scalar>& B_avg,
                             const double dt,
                             const Well::InjectionControls& inj_controls,
                             const Well::ProductionControls& prod_controls,
                             WellState& well_state,
                             Opm::DeferredLogger& deferred_logger)
    {
        const int max_iter_number = param_.max_inner_iter_ms_wells_;
        const WellState well_state0 = well_state;
        const std::vector<Scalar> residuals0 = getWellResiduals(B_avg);
        std::vector<std::vector<Scalar> > residual_history;
        std::vector<double> measure_history;
        int it = 0;
        // relaxation factor
        double relaxation_factor = 1.;
        const double min_relaxation_factor = 0.6;
        bool converged = false;
        int stagnate_count = 0;
        bool relax_convergence = false;
        for (; it < max_iter_number; ++it, ++debug_cost_counter_) {

            assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, deferred_logger);

            const BVectorWell dx_well = mswellhelpers::applyUMFPack(duneD_, duneDSolver_, resWell_);

            if (it > param_.strict_inner_iter_ms_wells_)
                relax_convergence = true;

            const auto report = getWellConvergence(well_state, B_avg, deferred_logger, relax_convergence);
            if (report.converged()) {
                converged = true;
                break;
            }

            residual_history.push_back(getWellResiduals(B_avg));
            measure_history.push_back(getResidualMeasureValue(well_state, residual_history[it], deferred_logger) );

            bool is_oscillate = false;
            bool is_stagnate = false;

            detectOscillations(measure_history, it, is_oscillate, is_stagnate);
            // TODO: maybe we should have more sophiscated strategy to recover the relaxation factor,
            // for example, to recover it to be bigger

            if (is_oscillate || is_stagnate) {
                // HACK!
                std::ostringstream sstr;
                if (relaxation_factor == min_relaxation_factor) {
                    // Still stagnating, terminate iterations if 5 iterations pass.
                    ++stagnate_count;
                    if (stagnate_count == 6) {
                        sstr << " well " << name() << " observes severe stagnation and/or oscillation. We relax the tolerance and check for convergence. \n";
                        const auto reportStag = getWellConvergence(well_state, B_avg, deferred_logger, true);
                        if (reportStag.converged()) {
                            converged = true;
                            sstr << " well " << name() << " manages to get converged with relaxed tolerances in " << it << " inner iterations";
                            deferred_logger.debug(sstr.str());
                            return converged;
                        }
                    }
                }

                // a factor value to reduce the relaxation_factor
                const double reduction_mutliplier = 0.9;
                relaxation_factor = std::max(relaxation_factor * reduction_mutliplier, min_relaxation_factor);

                // debug output
                if (is_stagnate) {
                    sstr << " well " << name() << " observes stagnation in inner iteration " << it << "\n";

                }
                if (is_oscillate) {
                    sstr << " well " << name() << " observes oscillation in inner iteration " << it << "\n";
                }
                sstr << " relaxation_factor is " << relaxation_factor << " now\n";
                deferred_logger.debug(sstr.str());
            }
            updateWellState(dx_well, well_state, deferred_logger, relaxation_factor);
            initPrimaryVariablesEvaluation();
        }

        // TODO: we should decide whether to keep the updated well_state, or recover to use the old well_state
        if (converged) {
            std::ostringstream sstr;
            sstr << "     Well " << name() << " converged in " << it << " inner iterations.";
            if (relax_convergence)
                sstr << "      (A relaxed tolerance was used after "<< param_.strict_inner_iter_ms_wells_ << " iterations)";
            deferred_logger.debug(sstr.str());
        } else {
            std::ostringstream sstr;
            sstr << "     Well " << name() << " did not converge in " << it << " inner iterations.";
#define EXTRA_DEBUG_MSW 0
#if EXTRA_DEBUG_MSW
            sstr << "***** Outputting the residual history for well " << name() << " during inner iterations:";
            for (int i = 0; i < it; ++i) {
                const auto& residual = residual_history[i];
                sstr << " residual at " << i << "th iteration ";
                for (const auto& res : residual) {
                    sstr << " " << res;
                }
                sstr << " " << measure_history[i] << " \n";
            }
#endif
            deferred_logger.debug(sstr.str());
        }

        return converged;
    }





    template<typename TypeTag>
    void
    MultisegmentWell<TypeTag>::
    assembleWellEqWithoutIteration(const Simulator& ebosSimulator,
                                   const double dt,
                                   const Well::InjectionControls& inj_controls,
                                   const Well::ProductionControls& prod_controls,
                                   WellState& well_state,
                                   Opm::DeferredLogger& deferred_logger)
    {

        // update the upwinding segments
        updateUpwindingSegments();

        // calculate the fluid properties needed.
        computeSegmentFluidProperties(ebosSimulator);

        // clear all entries
        duneB_ = 0.0;
        duneC_ = 0.0;

        duneD_ = 0.0;
        resWell_ = 0.0;

        duneDSolver_.reset();

        well_state.wellVaporizedOilRates()[index_of_well_] = 0.;
        well_state.wellDissolvedGasRates()[index_of_well_] = 0.;

        // for the black oil cases, there will be four equations,
        // the first three of them are the mass balance equations, the last one is the pressure equations.
        //
        // but for the top segment, the pressure equation will be the well control equation, and the other three will be the same.

        const bool allow_cf = getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator);

        const int nseg = numberOfSegments();

        for (int seg = 0; seg < nseg; ++seg) {
            // calculating the accumulation term
            // TODO: without considering the efficiencty factor for now
            {
                const EvalWell segment_surface_volume = getSegmentSurfaceVolume(ebosSimulator, seg);

                // Add a regularization_factor to increase the accumulation term
                // This will make the system less stiff and help convergence for 
                // difficult cases
                const Scalar regularization_factor =  param_.regularization_factor_ms_wells_;
                // for each component
                for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) {
                    const EvalWell accumulation_term = regularization_factor * (segment_surface_volume * surfaceVolumeFraction(seg, comp_idx)
                                                     - segment_fluid_initial_[seg][comp_idx]) / dt;

                    resWell_[seg][comp_idx] += accumulation_term.value();
                    for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
                        duneD_[seg][seg][comp_idx][pv_idx] += accumulation_term.derivative(pv_idx + numEq);
                    }
                }
            }
            // considering the contributions due to flowing out from the segment
            {
                for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) {
                    const EvalWell segment_rate = getSegmentRateUpwinding(seg, comp_idx) * well_efficiency_factor_;

                    const int seg_upwind = upwinding_segments_[seg];
                    // segment_rate contains the derivatives with respect to GTotal in seg,
                    // and WFrac and GFrac in seg_upwind
                    resWell_[seg][comp_idx] -= segment_rate.value();
                    duneD_[seg][seg][comp_idx][GTotal] -= segment_rate.derivative(GTotal + numEq);
                    if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
                        duneD_[seg][seg_upwind][comp_idx][WFrac] -= segment_rate.derivative(WFrac + numEq);
                    }
                    if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
                        duneD_[seg][seg_upwind][comp_idx][GFrac] -= segment_rate.derivative(GFrac + numEq);
                    }
                    // pressure derivative should be zero
                }
            }

            // considering the contributions from the inlet segments
            {
                for (const int inlet : segment_inlets_[seg]) {
                    for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) {
                        const EvalWell inlet_rate = getSegmentRateUpwinding(inlet, comp_idx) * well_efficiency_factor_;

                        const int inlet_upwind = upwinding_segments_[inlet];
                        // inlet_rate contains the derivatives with respect to GTotal in inlet,
                        // and WFrac and GFrac in inlet_upwind
                        resWell_[seg][comp_idx] += inlet_rate.value();
                        duneD_[seg][inlet][comp_idx][GTotal] += inlet_rate.derivative(GTotal + numEq);
                        if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
                            duneD_[seg][inlet_upwind][comp_idx][WFrac] += inlet_rate.derivative(WFrac + numEq);
                        }
                        if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
                            duneD_[seg][inlet_upwind][comp_idx][GFrac] += inlet_rate.derivative(GFrac + numEq);
                        }
                        // pressure derivative should be zero
                    }
                }
            }

            // calculating the perforation rate for each perforation that belongs to this segment
            const EvalWell seg_pressure = getSegmentPressure(seg);
            for (const int perf : segment_perforations_[seg]) {
                const int cell_idx = well_cells_[perf];
                const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
                std::vector<EvalWell> mob(num_components_, 0.0);
                getMobility(ebosSimulator, perf, mob);
                std::vector<EvalWell> cq_s(num_components_, 0.0);
                EvalWell perf_press;
                double perf_dis_gas_rate = 0.;
                double perf_vap_oil_rate = 0.;

                computePerfRatePressure(int_quants, mob, seg, perf, seg_pressure, allow_cf, cq_s, perf_press, perf_dis_gas_rate, perf_vap_oil_rate, deferred_logger);

                // updating the solution gas rate and solution oil rate
                if (this->isProducer()) {
                    well_state.wellDissolvedGasRates()[index_of_well_] += perf_dis_gas_rate;
                    well_state.wellVaporizedOilRates()[index_of_well_] += perf_vap_oil_rate;
                }

                // store the perf pressure and rates
                const int rate_start_offset = (first_perf_ + perf) * number_of_phases_;
                for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) {
                    well_state.perfPhaseRates()[rate_start_offset + ebosCompIdxToFlowCompIdx(comp_idx)] = cq_s[comp_idx].value();
                }
                well_state.perfPress()[first_perf_ + perf] = perf_press.value();

                for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) {
                    // the cq_s entering mass balance equations need to consider the efficiency factors.
                    const EvalWell cq_s_effective = cq_s[comp_idx] * well_efficiency_factor_;

                    connectionRates_[perf][comp_idx] = Base::restrictEval(cq_s_effective);

                    // subtract sum of phase fluxes in the well equations.
                    resWell_[seg][comp_idx] += cq_s_effective.value();

                    // assemble the jacobians
                    for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {

                        // also need to consider the efficiency factor when manipulating the jacobians.
                        duneC_[seg][cell_idx][pv_idx][comp_idx] -= cq_s_effective.derivative(pv_idx + numEq); // intput in transformed matrix

                        // the index name for the D should be eq_idx / pv_idx
                        duneD_[seg][seg][comp_idx][pv_idx] += cq_s_effective.derivative(pv_idx + numEq);
                    }

                    for (int pv_idx = 0; pv_idx < numEq; ++pv_idx) {
                        // also need to consider the efficiency factor when manipulating the jacobians.
                        duneB_[seg][cell_idx][comp_idx][pv_idx] += cq_s_effective.derivative(pv_idx);
                    }
                }
            }

            // the fourth dequation, the pressure drop equation
            if (seg == 0) { // top segment, pressure equation is the control equation
                const auto& summaryState = ebosSimulator.vanguard().summaryState();
                const Opm::Schedule& schedule = ebosSimulator.vanguard().schedule();
                assembleControlEq(well_state, schedule, summaryState, inj_controls, prod_controls, deferred_logger);
            } else {
                // TODO: maybe the following should go to the function assemblePressureEq()
		switch(segmentSet()[seg].segmentType()) {
		    case Segment::SegmentType::SICD :
		        assembleSICDPressureEq(seg, well_state);
			break;
		    case Segment::SegmentType::VALVE :
		        assembleValvePressureEq(seg, well_state);
			break;
		    default :
		        assemblePressureEq(seg, well_state);
		}
            }

            well_state.segPressDrop()[seg] = well_state.segPressDropHydroStatic()[seg] +
                                             well_state.segPressDropFriction()[seg] +
                                             well_state.segPressDropAcceleration()[seg];
        }
    }



    template<typename TypeTag>
    bool
    MultisegmentWell<TypeTag>::
    openCrossFlowAvoidSingularity(const Simulator& ebos_simulator) const
    {
        return !getAllowCrossFlow() && allDrawDownWrongDirection(ebos_simulator);
    }


    template<typename TypeTag>
    bool
    MultisegmentWell<TypeTag>::
    allDrawDownWrongDirection(const Simulator& ebos_simulator) const
    {
        bool all_drawdown_wrong_direction = true;
        const int nseg = numberOfSegments();

        for (int seg = 0; seg < nseg; ++seg) {
            const EvalWell segment_pressure = getSegmentPressure(seg);
            for (const int perf : segment_perforations_[seg]) {

                const int cell_idx = well_cells_[perf];
                const auto& intQuants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
                const auto& fs = intQuants.fluidState();

                // pressure difference between the segment and the perforation
                const EvalWell perf_seg_press_diff = gravity_ * segment_densities_[seg] * perforation_segment_depth_diffs_[perf];
                // pressure difference between the perforation and the grid cell
                const double cell_perf_press_diff = cell_perforation_pressure_diffs_[perf];

                const double pressure_cell = (fs.pressure(FluidSystem::oilPhaseIdx)).value();
                const double perf_press = pressure_cell - cell_perf_press_diff;
                // Pressure drawdown (also used to determine direction of flow)
                // TODO: not 100% sure about the sign of the seg_perf_press_diff
                const EvalWell drawdown = perf_press - (segment_pressure + perf_seg_press_diff);

                // for now, if there is one perforation can produce/inject in the correct
                // direction, we consider this well can still produce/inject.
                // TODO: it can be more complicated than this to cause wrong-signed rates
                if ( (drawdown < 0. && this->isInjector()) ||
                     (drawdown > 0. && this->isProducer()) )  {
                    all_drawdown_wrong_direction = false;
                    break;
                }
            }
        }

        return all_drawdown_wrong_direction;
    }


    template<typename TypeTag>
    void
    MultisegmentWell<TypeTag>::
    wellTestingPhysical(const Simulator& /* simulator */, const std::vector<double>& /* B_avg */,
                        const double /* simulation_time */, const int /* report_step */,
                        WellState& /* well_state */, WellTestState& /* welltest_state */, Opm::DeferredLogger& deferred_logger)
    {
        const std::string msg = "Support of well testing for physical limits for multisegment wells is not "
                                "implemented yet, wellTestingPhysical() for " + name() + " will do nothing";
        deferred_logger.warning("NO_WELLTESTPHYSICAL_CHECKING_MS_WELLS", msg);
    }





    template<typename TypeTag>
    void
    MultisegmentWell<TypeTag>::
    updateWaterThroughput(const double dt OPM_UNUSED, WellState& well_state OPM_UNUSED) const
    {
    }





    template<typename TypeTag>
    typename MultisegmentWell<TypeTag>::EvalWell
    MultisegmentWell<TypeTag>::
    getSegmentSurfaceVolume(const Simulator& ebos_simulator, const int seg_idx) const
    {
        EvalWell temperature;
        EvalWell saltConcentration;
        int pvt_region_index;
        {
            // using the pvt region of first perforated cell
            // TODO: it should be a member of the WellInterface, initialized properly
            const int cell_idx = well_cells_[0];
            const auto& intQuants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
            const auto& fs = intQuants.fluidState();
            temperature.setValue(fs.temperature(FluidSystem::oilPhaseIdx).value());
            saltConcentration = extendEval(fs.saltConcentration());
            pvt_region_index = fs.pvtRegionIndex();
        }

        const EvalWell seg_pressure = getSegmentPressure(seg_idx);

        std::vector<EvalWell> mix_s(num_components_, 0.0);
        for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) {
            mix_s[comp_idx] = surfaceVolumeFraction(seg_idx, comp_idx);
        }

        std::vector<EvalWell> b(num_components_, 0.);
        if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
            const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
            b[waterCompIdx] =
                FluidSystem::waterPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, saltConcentration);
        }

        EvalWell rv(0.0);
        // gas phase
        if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
            const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
            if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
                const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
                EvalWell rvmax = FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvt_region_index, temperature, seg_pressure);
                if (rvmax < 0.0) { // negative rvmax can happen if the seg_pressure is outside the range of the table
                    rvmax = 0.0;
                }
                if (mix_s[oilCompIdx] > 0.0) {
                    if (mix_s[gasCompIdx] > 0.0) {
                        rv = mix_s[oilCompIdx] / mix_s[gasCompIdx];
                    }

                    if (rv > rvmax) {
                        rv = rvmax;
                    }
                    b[gasCompIdx] =
                        FluidSystem::gasPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, rv);
                } else { // no oil exists
                    b[gasCompIdx] =
                        FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure);
                }
            } else { // no Liquid phase
                // it is the same with zero mix_s[Oil]
                b[gasCompIdx] =
                    FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure);
            }
        }

        EvalWell rs(0.0);
        // oil phase
        if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
            const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
            if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
                const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
                EvalWell rsmax = FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvt_region_index, temperature, seg_pressure);
                if (rsmax < 0.0) { // negative rsmax can happen if the seg_pressure is outside the range of the table
                    rsmax = 0.0;
                }
                if (mix_s[gasCompIdx] > 0.0) {
                    if (mix_s[oilCompIdx] > 0.0) {
                        rs = mix_s[gasCompIdx] / mix_s[oilCompIdx];
                    }
                    // std::cout << " rs " << rs.value() << " rsmax " << rsmax.value() << std::endl;

                    if (rs > rsmax) {
                        rs = rsmax;
                    }
                    b[oilCompIdx] =
                        FluidSystem::oilPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, rs);
                } else { // no oil exists
                    b[oilCompIdx] =
                        FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure);
                }
            } else { // no gas phase
                // it is the same with zero mix_s[Gas]
                b[oilCompIdx] =
                    FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure);
            }
        }

        std::vector<EvalWell> mix(mix_s);
        if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
            const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
            const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);

            const EvalWell d = 1.0 - rs * rv;
            if (d <= 0.0 || d > 1.0) {
                OPM_THROW(Opm::NumericalIssue, "Problematic d value " << d << " obtained for well " << name()
                                               << " during convertion to surface volume with rs " << rs
                                               << ", rv " << rv << " and pressure " << seg_pressure
                                               << " obtaining d " << d);
            }

            if (rs > 0.0) { // rs > 0.0?
                mix[gasCompIdx] = (mix_s[gasCompIdx] - mix_s[oilCompIdx] * rs) / d;
            }
            if (rv > 0.0) { // rv > 0.0?
                mix[oilCompIdx] = (mix_s[oilCompIdx] - mix_s[gasCompIdx] * rv) / d;
            }
        }

        EvalWell vol_ratio(0.0);
        for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) {
            vol_ratio += mix[comp_idx] / b[comp_idx];
        }

        // We increase the segment volume with a factor 10 to stabilize the system.
        const double volume = segmentSet()[seg_idx].volume();

        return volume / vol_ratio;
    }





    template<typename TypeTag>
    std::vector<typename MultisegmentWell<TypeTag>::Scalar>
    MultisegmentWell<TypeTag>::
    getWellResiduals(const std::vector<Scalar>& B_avg) const
    {
        assert(int(B_avg.size() ) == num_components_);
        std::vector<Scalar> residuals(numWellEq + 1, 0.0);

        for (int seg = 0; seg < numberOfSegments(); ++seg) {
            for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) {
                double residual = 0.;
                if (eq_idx < num_components_) {
                    residual = std::abs(resWell_[seg][eq_idx]) * B_avg[eq_idx];
                } else {
                    if (seg > 0) {
                        residual = std::abs(resWell_[seg][eq_idx]);
                    }
                }
                if (std::isnan(residual) || std::isinf(residual)) {
                    OPM_THROW(Opm::NumericalIssue, "nan or inf value for residal get for well " << name()
                                                    << " segment " << seg << " eq_idx " << eq_idx);
                }

                if (residual > residuals[eq_idx]) {
                    residuals[eq_idx] = residual;
                }
            }
        }

        // handling the control equation residual
        {
            const double control_residual = std::abs(resWell_[0][numWellEq - 1]);
            if (std::isnan(control_residual) || std::isinf(control_residual)) {
               OPM_THROW(Opm::NumericalIssue, "nan or inf value for control residal get for well " << name());
            }
            residuals[numWellEq] = control_residual;
        }

        return residuals;
    }





    /// Detect oscillation or stagnation based on the residual measure history
    template<typename TypeTag>
    void
    MultisegmentWell<TypeTag>::
    detectOscillations(const std::vector<double>& measure_history,
                       const int it, bool& oscillate, bool& stagnate) const
    {
        if ( it < 2 ) {
            oscillate = false;
            stagnate = false;
            return;
        }

        stagnate = true;
        const double F0 = measure_history[it];
        const double F1 = measure_history[it - 1];
        const double F2 = measure_history[it - 2];
        const double d1 = std::abs((F0 - F2) / F0);
        const double d2 = std::abs((F0 - F1) / F0);

        const double oscillaton_rel_tol = 0.2;
        oscillate = (d1 < oscillaton_rel_tol) && (oscillaton_rel_tol < d2);

        const double stagnation_rel_tol = 1.e-2;
        stagnate = std::abs((F1 - F2) / F2) <= stagnation_rel_tol;
    }





    template<typename TypeTag>
    double
    MultisegmentWell<TypeTag>::
    getResidualMeasureValue(const WellState& well_state,
                            const std::vector<double>& residuals,
                            DeferredLogger& deferred_logger) const
    {
        assert(int(residuals.size()) == numWellEq + 1);

        const double rate_tolerance = param_.tolerance_wells_;
        int count = 0;
        double sum = 0;
        for (int eq_idx = 0; eq_idx < numWellEq - 1; ++eq_idx) {
            if (residuals[eq_idx] > rate_tolerance) {
                sum += residuals[eq_idx] / rate_tolerance;
                ++count;
            }
        }

        const double pressure_tolerance = param_.tolerance_pressure_ms_wells_;
        if (residuals[SPres] > pressure_tolerance) {
            sum += residuals[SPres] / pressure_tolerance;
            ++count;
        }

        const double control_tolerance = getControlTolerance(well_state, deferred_logger);
        if (residuals[SPres + 1] > control_tolerance) {
            sum += residuals[SPres + 1] / control_tolerance;
            ++count;
        }

        // if (count == 0), it should be converged.
        assert(count != 0);

        return sum;
    }





    template<typename TypeTag>
    double
    MultisegmentWell<TypeTag>::
    getControlTolerance(const WellState& well_state,
                        DeferredLogger& deferred_logger) const
    {
        double control_tolerance = 0.;

        const int well_index = index_of_well_;
        if (this->isInjector() )
        {
            const Opm::Well::InjectorCMode& current = well_state.currentInjectionControls()[well_index];
            switch(current) {
            case Well::InjectorCMode::THP:
                control_tolerance = param_.tolerance_pressure_ms_wells_;
                break;
            case Well::InjectorCMode::BHP:
                control_tolerance = param_.tolerance_wells_;
                break;
            case Well::InjectorCMode::RATE:
            case Well::InjectorCMode::RESV:
                control_tolerance = param_.tolerance_wells_;
                break;
            case Well::InjectorCMode::GRUP:
                control_tolerance = param_.tolerance_wells_;
                break;
            default:
                OPM_DEFLOG_THROW(std::runtime_error, "Unknown well control control types for well " << name(), deferred_logger);
            }
        }

        if (this->isProducer() )
        {
            const Well::ProducerCMode& current = well_state.currentProductionControls()[well_index];
            switch(current) {
            case Well::ProducerCMode::THP:
                control_tolerance = param_.tolerance_pressure_ms_wells_; // 0.1 bar
                break;
            case Well::ProducerCMode::BHP:
                control_tolerance = param_.tolerance_wells_; // 0.01 bar
                break;
            case Well::ProducerCMode::ORAT:
            case Well::ProducerCMode::WRAT:
            case Well::ProducerCMode::GRAT:
            case Well::ProducerCMode::LRAT:
            case Well::ProducerCMode::RESV:
            case Well::ProducerCMode::CRAT:
                control_tolerance = param_.tolerance_wells_; // smaller tolerance for rate control
                break;
            case Well::ProducerCMode::GRUP:
                control_tolerance = param_.tolerance_wells_; // smaller tolerance for rate control
                break;
            default:
                OPM_DEFLOG_THROW(std::runtime_error, "Unknown well control control types for well " << name(), deferred_logger);
            }
        }

        return control_tolerance;
    }





    template<typename TypeTag>
    void
    MultisegmentWell<TypeTag>::
    checkConvergenceControlEq(const WellState& well_state,
                              ConvergenceReport& report,
                              DeferredLogger& deferred_logger) const
    {
        double control_tolerance = 0.;
        using CR = ConvergenceReport;
        CR::WellFailure::Type ctrltype = CR::WellFailure::Type::Invalid;

        const int well_index = index_of_well_;
        if (this->isInjector() )
        {
            const Opm::Well::InjectorCMode& current = well_state.currentInjectionControls()[well_index];
            switch(current) {
            case Well::InjectorCMode::THP:
                ctrltype = CR::WellFailure::Type::ControlTHP;
                control_tolerance = param_.tolerance_pressure_ms_wells_;
                break;
            case Well::InjectorCMode::BHP:
                ctrltype = CR::WellFailure::Type::ControlBHP;
                control_tolerance = param_.tolerance_pressure_ms_wells_;
                break;
            case Well::InjectorCMode::RATE:
            case Well::InjectorCMode::RESV:
                ctrltype = CR::WellFailure::Type::ControlRate;
                control_tolerance = param_.tolerance_wells_;
                break;
            case Well::InjectorCMode::GRUP:
                ctrltype = CR::WellFailure::Type::ControlRate;
                control_tolerance = param_.tolerance_wells_;
                break;
            default:
                OPM_DEFLOG_THROW(std::runtime_error, "Unknown well control control types for well " << name(), deferred_logger);
            }
        }

        if (this->isProducer() )
        {
            const Well::ProducerCMode& current = well_state.currentProductionControls()[well_index];
            switch(current) {
            case Well::ProducerCMode::THP:
                ctrltype = CR::WellFailure::Type::ControlTHP;
                control_tolerance = param_.tolerance_pressure_ms_wells_;
                break;
            case Well::ProducerCMode::BHP:
                ctrltype = CR::WellFailure::Type::ControlBHP;
                control_tolerance = param_.tolerance_pressure_ms_wells_;
                break;
            case Well::ProducerCMode::ORAT:
            case Well::ProducerCMode::WRAT:
            case Well::ProducerCMode::GRAT:
            case Well::ProducerCMode::LRAT:
            case Well::ProducerCMode::RESV:
            case Well::ProducerCMode::CRAT:
                ctrltype = CR::WellFailure::Type::ControlRate;
                control_tolerance = param_.tolerance_wells_;
                break;
            case Well::ProducerCMode::GRUP:
                ctrltype = CR::WellFailure::Type::ControlRate;
                control_tolerance = param_.tolerance_wells_;
                break;
            default:
                OPM_DEFLOG_THROW(std::runtime_error, "Unknown well control control types for well " << name(), deferred_logger);
            }
        }

        const double well_control_residual = std::abs(resWell_[0][SPres]);
        const int dummy_component = -1;
        const double max_residual_allowed = param_.max_residual_allowed_;
        if (std::isnan(well_control_residual)) {
            report.setWellFailed({ctrltype, CR::Severity::NotANumber, dummy_component, name()});
        } else if (well_control_residual > max_residual_allowed * 10.) {
            report.setWellFailed({ctrltype, CR::Severity::TooLarge, dummy_component, name()});
        } else if ( well_control_residual > control_tolerance) {
            report.setWellFailed({ctrltype, CR::Severity::Normal, dummy_component, name()});
        }
    }






    template<typename TypeTag>
    void
    MultisegmentWell<TypeTag>::
    updateUpwindingSegments()
    {
        for (int seg = 0; seg < numberOfSegments(); ++seg) {
            // special treatment is needed for segment 0
            if (seg == 0) {
                // we are not supposed to have injecting producers and producing injectors
                assert( ! (this->isProducer() && primary_variables_evaluation_[seg][GTotal] > 0.) );
                assert( ! (this->isInjector() && primary_variables_evaluation_[seg][GTotal] < 0.) );
                upwinding_segments_[seg] = seg;
                continue;
            }

            // for other normal segments
            if (primary_variables_evaluation_[seg][GTotal] <= 0.) {
                upwinding_segments_[seg] = seg;
            } else {
                const int outlet_segment_index = segmentNumberToIndex(segmentSet()[seg].outletSegment());
                upwinding_segments_[seg] = outlet_segment_index;
            }
        }
    }












    template<typename TypeTag>
    void
    MultisegmentWell<TypeTag>::
    assembleSICDPressureEq(const int seg, WellState& well_state) const
    {
        // TODO: upwinding needs to be taken care of
        // top segment can not be a spiral ICD device
        assert(seg != 0);

        // the pressure equation is something like
        // p_seg - deltaP - p_outlet = 0.
        // the major part is how to calculate the deltaP

        EvalWell pressure_equation = getSegmentPressure(seg);

        const auto sicd_pressure_drop = pressureDropSpiralICD(seg);
        pressure_equation = pressure_equation - sicd_pressure_drop;
        well_state.segPressDropFriction()[seg] = sicd_pressure_drop.value();

        const int seg_upwind = upwinding_segments_[seg];
        resWell_[seg][SPres] = pressure_equation.value();
        duneD_[seg][seg][SPres][SPres] += pressure_equation.derivative(SPres + numEq);
        duneD_[seg][seg][SPres][GTotal] += pressure_equation.derivative(GTotal + numEq);
        if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
            duneD_[seg][seg_upwind][SPres][WFrac] += pressure_equation.derivative(WFrac + numEq);
        }
        if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
            duneD_[seg][seg_upwind][SPres][GFrac] += pressure_equation.derivative(GFrac + numEq);
        }

        // contribution from the outlet segment
        const int outlet_segment_index = segmentNumberToIndex(segmentSet()[seg].outletSegment());
        const EvalWell outlet_pressure = getSegmentPressure(outlet_segment_index);

        resWell_[seg][SPres] -= outlet_pressure.value();
        for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
            duneD_[seg][outlet_segment_index][SPres][pv_idx] = -outlet_pressure.derivative(pv_idx + numEq);

        }
    }




    template<typename TypeTag>
    void
    MultisegmentWell<TypeTag>::
    assembleValvePressureEq(const int seg, WellState& well_state) const
    {
        // top segment can not be a spiral ICD device
        assert(seg != 0);

        // const Valve& valve = *segmentSet()[seg].Valve();

        // the pressure equation is something like
        // p_seg - deltaP - p_outlet = 0.
        // the major part is how to calculate the deltaP

        EvalWell pressure_equation = getSegmentPressure(seg);

        const int seg_upwind = upwinding_segments_[seg];

        const auto valve_pressure_drop = pressureDropValve(seg);
        pressure_equation = pressure_equation - valve_pressure_drop;
        well_state.segPressDropFriction()[seg] = valve_pressure_drop.value();

        resWell_[seg][SPres] = pressure_equation.value();
        duneD_[seg][seg][SPres][SPres] += pressure_equation.derivative(SPres + numEq);
        duneD_[seg][seg][SPres][GTotal] += pressure_equation.derivative(GTotal + numEq);
        if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
            duneD_[seg][seg_upwind][SPres][WFrac] += pressure_equation.derivative(WFrac + numEq);
        }
        if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
            duneD_[seg][seg_upwind][SPres][GFrac] += pressure_equation.derivative(GFrac + numEq);
        }

        // contribution from the outlet segment
        const int outlet_segment_index = segmentNumberToIndex(segmentSet()[seg].outletSegment());
        const EvalWell outlet_pressure = getSegmentPressure(outlet_segment_index);

        resWell_[seg][SPres] -= outlet_pressure.value();
        for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
            duneD_[seg][outlet_segment_index][SPres][pv_idx] = -outlet_pressure.derivative(pv_idx + numEq);
        }
    }





    template<typename TypeTag>
    std::optional<double>
    MultisegmentWell<TypeTag>::
    computeBhpAtThpLimitProd(const Simulator& ebos_simulator,
                             const std::vector<Scalar>& B_avg,
                             const SummaryState& summary_state,
                             DeferredLogger& deferred_logger) const
    {
        // Given a VFP function returning bhp as a function of phase
        // rates and thp:
        //     fbhp(rates, thp),
        // a function extracting the particular flow rate used for VFP
        // lookups:
        //     flo(rates)
        // and the inflow function (assuming the reservoir is fixed):
        //     frates(bhp)
        // we want to solve the equation:
        //     fbhp(frates(bhp, thplimit)) - bhp = 0
        // for bhp.
        //
        // This may result in 0, 1 or 2 solutions. If two solutions,
        // the one corresponding to the lowest bhp (and therefore
        // highest rate) should be returned.

        // Make the fbhp() function.
        const auto& controls = well_ecl_.productionControls(summary_state);
        const auto& table = *(vfp_properties_->getProd()->getTable(controls.vfp_table_number));
        const double vfp_ref_depth = table.getDatumDepth();
        const double rho = segment_densities_[0].value(); // Use the density at the top perforation.
        const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_);
        auto fbhp = [this, &controls, dp](const std::vector<double>& rates) {
            assert(rates.size() == 3);
            return this->vfp_properties_->getProd()
            ->bhp(controls.vfp_table_number, rates[Water], rates[Oil], rates[Gas], controls.thp_limit, controls.alq_value) - dp;
        };

        // Make the flo() function.
        auto flo_type = table.getFloType();
        auto flo = [flo_type](const std::vector<double>& rates) {
            return detail::getFlo(rates[Water], rates[Oil], rates[Gas], flo_type);
        };

        // Make the frates() function.
        auto frates = [this, &ebos_simulator, &B_avg, &deferred_logger](const double bhp) {
            // Not solving the well equations here, which means we are
            // calculating at the current Fg/Fw values of the
            // well. This does not matter unless the well is
            // crossflowing, and then it is likely still a good
            // approximation.
            std::vector<double> rates(3);
            computeWellRatesWithBhp(ebos_simulator, B_avg, bhp, rates, deferred_logger);
            return rates;
        };

        // Find the bhp-point where production becomes nonzero.
        double bhp_max = 0.0;
        {
            auto fflo = [&flo, &frates](double bhp) { return flo(frates(bhp)); };
            double low = controls.bhp_limit;
            double high = maxPerfPress(ebos_simulator) + 1.0 * unit::barsa;
            double f_low = fflo(low);
            double f_high = fflo(high);
            deferred_logger.debug("computeBhpAtThpLimitProd(): well = " + name() +
                                  "  low = " + std::to_string(low) +
                                  "  high = " + std::to_string(high) +
                                  "  f(low) = " + std::to_string(f_low) +
                                  "  f(high) = " + std::to_string(f_high));
            int adjustments = 0;
            const int max_adjustments = 10;
            const double adjust_amount = 5.0 * unit::barsa;
            while (f_low * f_high > 0.0 && adjustments < max_adjustments) {
                // Same sign, adjust high to see if we can flip it.
                high += adjust_amount;
                f_high = fflo(high);
                ++adjustments;
            }
            if (f_low * f_high > 0.0) {
                if (f_low > 0.0) {
                    // Even at the BHP limit, we are injecting.
                    // There will be no solution here, return an
                    // empty optional.
                    deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE_INOPERABLE",
                                            "Robust bhp(thp) solve failed due to inoperability for well " + name());
                    return std::optional<double>();
                } else {
                    // Still producing, even at high bhp.
                    assert(f_high < 0.0);
                    bhp_max = high;
                }
            } else {
                // Bisect to find a bhp point where we produce, but
                // not a large amount ('eps' below).
                const double eps = 0.1 * std::fabs(table.getFloAxis().front());
                const int maxit = 50;
                int it = 0;
                while (std::fabs(f_low) > eps && it < maxit) {
                    const double curr = 0.5*(low + high);
                    const double f_curr = fflo(curr);
                    if (f_curr * f_low > 0.0) {
                        low = curr;
                        f_low = f_curr;
                    } else {
                        high = curr;
                        f_high = f_curr;
                    }
                    ++it;
                }
                bhp_max = low;
            }
            deferred_logger.debug("computeBhpAtThpLimitProd(): well = " + name() +
                                  "  low = " + std::to_string(low) +
                                  "  high = " + std::to_string(high) +
                                  "  f(low) = " + std::to_string(f_low) +
                                  "  f(high) = " + std::to_string(f_high) +
                                  "  bhp_max = " + std::to_string(bhp_max));
        }

        // Define the equation we want to solve.
        auto eq = [&fbhp, &frates](double bhp) {
            return fbhp(frates(bhp)) - bhp;
        };

        // Find appropriate brackets for the solution.
        double low = controls.bhp_limit;
        double high = bhp_max;
        {
            double eq_high = eq(high);
            double eq_low = eq(low);
            const double eq_bhplimit = eq_low;
            deferred_logger.debug("computeBhpAtThpLimitProd(): well = " + name() +
                                  "  low = " + std::to_string(low) +
                                  "  high = " + std::to_string(high) +
                                  "  eq(low) = " + std::to_string(eq_low) +
                                  "  eq(high) = " + std::to_string(eq_high));
            if (eq_low * eq_high > 0.0) {
                // Failed to bracket the zero.
                // If this is due to having two solutions, bisect until bracketed.
                double abs_low = std::fabs(eq_low);
                double abs_high = std::fabs(eq_high);
                int bracket_attempts = 0;
                const int max_bracket_attempts = 20;
                double interval = high - low;
                const double min_interval = 1.0 * unit::barsa;
                while (eq_low * eq_high > 0.0 && bracket_attempts < max_bracket_attempts && interval > min_interval) {
                    if (abs_high < abs_low) {
                        low = 0.5 * (low + high);
                        eq_low = eq(low);
                        abs_low = std::fabs(eq_low);
                    } else {
                        high = 0.5 * (low + high);
                        eq_high = eq(high);
                        abs_high = std::fabs(eq_high);
                    }
                    ++bracket_attempts;
                }
                if (eq_low * eq_high > 0.0) {
                    // Still failed bracketing!
                    const double limit = 3.0 * unit::barsa;
                    if (std::min(abs_low, abs_high) < limit) {
                        // Return the least bad solution if less off than 3 bar.
                        deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE_BRACKETING_FAILURE",
                                                "Robust bhp(thp) not solved precisely for well " + name());
                        return abs_low < abs_high ? low : high;
                    } else {
                        // Return failure.
                        deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE_BRACKETING_FAILURE",
                                                "Robust bhp(thp) solve failed due to bracketing failure for well " + name());
                        return std::optional<double>();
                    }
                }
            }
            // We have a bracket!
            // Now, see if (bhplimit, low) is a bracket in addition to (low, high).
            // If so, that is the bracket we shall use, choosing the solution with the
            // highest flow.
            if (eq_low * eq_bhplimit <= 0.0) {
                high = low;
                low = controls.bhp_limit;
            }
        }

        // Solve for the proper solution in the given interval.
        const int max_iteration = 100;
        const double bhp_tolerance = 0.01 * unit::barsa;
        int iteration = 0;
        try {
            const double solved_bhp = RegulaFalsiBisection<ThrowOnError>::
                solve(eq, low, high, max_iteration, bhp_tolerance, iteration);
            return solved_bhp;
        }
        catch (...) {
            deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE",
                                    "Robust bhp(thp) solve failed for well " + name());
            return std::optional<double>();
	}
    }




    template<typename TypeTag>
    std::optional<double>
    MultisegmentWell<TypeTag>::
    computeBhpAtThpLimitInj(const Simulator& ebos_simulator,
                            const std::vector<Scalar>& B_avg,
                            const SummaryState& summary_state,
                            DeferredLogger& deferred_logger) const
    {
        // Given a VFP function returning bhp as a function of phase
        // rates and thp:
        //     fbhp(rates, thp),
        // a function extracting the particular flow rate used for VFP
        // lookups:
        //     flo(rates)
        // and the inflow function (assuming the reservoir is fixed):
        //     frates(bhp)
        // we want to solve the equation:
        //     fbhp(frates(bhp, thplimit)) - bhp = 0
        // for bhp.
        //
        // This may result in 0, 1 or 2 solutions. If two solutions,
        // the one corresponding to the lowest bhp (and therefore
        // highest rate) is returned.
        //
        // In order to detect these situations, we will find piecewise
        // linear approximations both to the inverse of the frates
        // function and to the fbhp function.
        //
        // We first take the FLO sample points of the VFP curve, and
        // find the corresponding bhp values by solving the equation:
        //     flo(frates(bhp)) - flo_sample = 0
        // for bhp, for each flo_sample. The resulting (flo_sample,
        // bhp_sample) values give a piecewise linear approximation to
        // the true inverse inflow function, at the same flo values as
        // the VFP data.
        //
        // Then we extract a piecewise linear approximation from the
        // multilinear fbhp() by evaluating it at the flo_sample
        // points, with fractions given by the frates(bhp_sample)
        // values.
        //
        // When we have both piecewise linear curves defined on the
        // same flo_sample points, it is easy to distinguish between
        // the 0, 1 or 2 solution cases, and obtain the right interval
        // in which to solve for the solution we want (with highest
        // flow in case of 2 solutions).

        // Make the fbhp() function.
        const auto& controls = well_ecl_.injectionControls(summary_state);
        const auto& table = *(vfp_properties_->getInj()->getTable(controls.vfp_table_number));
        const double vfp_ref_depth = table.getDatumDepth();
        const double rho = segment_densities_[0].value(); // Use the density at the top perforation.
        const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_);
        auto fbhp = [this, &controls, dp](const std::vector<double>& rates) {
            assert(rates.size() == 3);
            return this->vfp_properties_->getInj()
                    ->bhp(controls.vfp_table_number, rates[Water], rates[Oil], rates[Gas], controls.thp_limit) - dp;
        };

        // Make the flo() function.
        auto flo_type = table.getFloType();
        auto flo = [flo_type](const std::vector<double>& rates) {
            return detail::getFlo(rates[Water], rates[Oil], rates[Gas], flo_type);
        };

        // Make the frates() function.
        auto frates = [this, &ebos_simulator, &B_avg, &deferred_logger](const double bhp) {
            // Not solving the well equations here, which means we are
            // calculating at the current Fg/Fw values of the
            // well. This does not matter unless the well is
            // crossflowing, and then it is likely still a good
            // approximation.
            std::vector<double> rates(3);
            computeWellRatesWithBhp(ebos_simulator, B_avg, bhp, rates, deferred_logger);
            return rates;
        };

        // Get the flo samples, add extra samples at low rates and bhp
        // limit point if necessary.
        std::vector<double> flo_samples = table.getFloAxis();
        if (flo_samples[0] > 0.0) {
            const double f0 = flo_samples[0];
            flo_samples.insert(flo_samples.begin(), { f0/20.0, f0/10.0, f0/5.0, f0/2.0 });
        }
        const double flo_bhp_limit = flo(frates(controls.bhp_limit));
        if (flo_samples.back() < flo_bhp_limit) {
            flo_samples.push_back(flo_bhp_limit);
        }

        // Find bhp values for inflow relation corresponding to flo samples.
        std::vector<double> bhp_samples;
        for (double flo_sample : flo_samples) {
            if (flo_sample > flo_bhp_limit) {
                // We would have to go over the bhp limit to obtain a
                // flow of this magnitude. We associate all such flows
                // with simply the bhp limit. The first one
                // encountered is considered valid, the rest not. They
                // are therefore skipped.
                bhp_samples.push_back(controls.bhp_limit);
                break;
            }
            auto eq = [&flo, &frates, flo_sample](double bhp) {
                return flo(frates(bhp)) - flo_sample;
            };
            // TODO: replace hardcoded low/high limits.
            const double low = 10.0 * unit::barsa;
            const double high = 800.0 * unit::barsa;
            const int max_iteration = 100;
            const double flo_tolerance = 0.05 * std::fabs(flo_samples.back());
            int iteration = 0;
            try {
                const double solved_bhp = RegulaFalsiBisection<WarnAndContinueOnError>::
                        solve(eq, low, high, max_iteration, flo_tolerance, iteration);
                bhp_samples.push_back(solved_bhp);
            }
            catch (...) {
                // Use previous value (or max value if at start) if we failed.
                bhp_samples.push_back(bhp_samples.empty() ? low : bhp_samples.back());
                deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE_EXTRACT_SAMPLES",
                                        "Robust bhp(thp) solve failed extracting bhp values at flo samples for well " + name());
            }
        }

        // Find bhp values for VFP relation corresponding to flo samples.
        const int num_samples = bhp_samples.size(); // Note that this can be smaller than flo_samples.size()
        std::vector<double> fbhp_samples(num_samples);
        for (int ii = 0; ii < num_samples; ++ii) {
            fbhp_samples[ii] = fbhp(frates(bhp_samples[ii]));
        }
// #define EXTRA_THP_DEBUGGING
#ifdef EXTRA_THP_DEBUGGING
        std::string dbgmsg;
        dbgmsg += "flo: ";
        for (int ii = 0; ii < num_samples; ++ii) {
            dbgmsg += "  " + std::to_string(flo_samples[ii]);
        }
        dbgmsg += "\nbhp: ";
        for (int ii = 0; ii < num_samples; ++ii) {
            dbgmsg += "  " + std::to_string(bhp_samples[ii]);
        }
        dbgmsg += "\nfbhp: ";
        for (int ii = 0; ii < num_samples; ++ii) {
            dbgmsg += "  " + std::to_string(fbhp_samples[ii]);
        }
        OpmLog::debug(dbgmsg);
#endif // EXTRA_THP_DEBUGGING

        // Look for sign changes for the (fbhp_samples - bhp_samples) piecewise linear curve.
        // We only look at the valid
        int sign_change_index = -1;
        for (int ii = 0; ii < num_samples - 1; ++ii) {
            const double curr = fbhp_samples[ii] - bhp_samples[ii];
            const double next = fbhp_samples[ii + 1] - bhp_samples[ii + 1];
            if (curr * next < 0.0) {
                // Sign change in the [ii, ii + 1] interval.
                sign_change_index = ii; // May overwrite, thereby choosing the highest-flo solution.
            }
        }

        // Handle the no solution case.
        if (sign_change_index == -1) {
            return std::optional<double>();
        }

        // Solve for the proper solution in the given interval.
        auto eq = [&fbhp, &frates](double bhp) {
            return fbhp(frates(bhp)) - bhp;
        };
        // TODO: replace hardcoded low/high limits.
        const double low = bhp_samples[sign_change_index + 1];
        const double high = bhp_samples[sign_change_index];
        const int max_iteration = 100;
        const double bhp_tolerance = 0.01 * unit::barsa;
        int iteration = 0;
        if (low == high) {
            // We are in the high flow regime where the bhp_samples
            // are all equal to the bhp_limit.
            assert(low == controls.bhp_limit);
            deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE",
                                    "Robust bhp(thp) solve failed for well " + name());
            return std::optional<double>();
        }
        try {
            const double solved_bhp = RegulaFalsiBisection<WarnAndContinueOnError>::
                    solve(eq, low, high, max_iteration, bhp_tolerance, iteration);
#ifdef EXTRA_THP_DEBUGGING
            OpmLog::debug("*****    " + name() + "    solved_bhp = " + std::to_string(solved_bhp)
                          + "    flo_bhp_limit = " + std::to_string(flo_bhp_limit));
#endif // EXTRA_THP_DEBUGGING
            return solved_bhp;
        }
        catch (...) {
            deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE",
                                    "Robust bhp(thp) solve failed for well " + name());
            return std::optional<double>();
        }

    }





    template<typename TypeTag>
    double
    MultisegmentWell<TypeTag>::
    maxPerfPress(const Simulator& ebos_simulator) const
    {
        double max_pressure = 0.0;
        const int nseg = numberOfSegments();
        for (int seg = 0; seg < nseg; ++seg) {
            for (const int perf : segment_perforations_[seg]) {
                const int cell_idx = well_cells_[perf];
                const auto& int_quants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
                const auto& fs = int_quants.fluidState();
                double pressure_cell = fs.pressure(FluidSystem::oilPhaseIdx).value();
                max_pressure = std::max(max_pressure, pressure_cell);
            }
        }
        return max_pressure;
    }





    template<typename TypeTag>
    typename MultisegmentWell<TypeTag>::EvalWell
    MultisegmentWell<TypeTag>::
    pressureDropSpiralICD(const int seg) const
    {
        const SICD& sicd = segmentSet()[seg].spiralICD();

        const int seg_upwind = upwinding_segments_[seg];
        const std::vector<EvalWell>& phase_fractions = segment_phase_fractions_[seg_upwind];
        const std::vector<EvalWell>& phase_viscosities = segment_phase_viscosities_[seg_upwind];

        EvalWell water_fraction = 0.;
        EvalWell water_viscosity = 0.;
        if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
            const int water_pos = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
            water_fraction = phase_fractions[water_pos];
            water_viscosity = phase_viscosities[water_pos];
        }

        EvalWell oil_fraction = 0.;
        EvalWell oil_viscosity = 0.;
        if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
            const int oil_pos = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
            oil_fraction = phase_fractions[oil_pos];
            oil_viscosity = phase_viscosities[oil_pos];
        }

        EvalWell gas_fraction = 0.;
        EvalWell gas_viscosity = 0.;
        if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
            const int gas_pos = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
            gas_fraction = phase_fractions[gas_pos];
            gas_viscosity = phase_viscosities[gas_pos];
        }

        EvalWell density = segment_densities_[seg_upwind];
        // WARNING
        // We disregard the derivatives from the upwind density to make sure derivatives
        // wrt. to different segments dont get mixed.
        if (seg != seg_upwind) {
            water_fraction.clearDerivatives();
            water_viscosity.clearDerivatives();
            oil_fraction.clearDerivatives();
            oil_viscosity.clearDerivatives();
            gas_fraction.clearDerivatives();
            gas_viscosity.clearDerivatives();
            density.clearDerivatives();
        }

        const EvalWell liquid_emulsion_viscosity = mswellhelpers::emulsionViscosity(water_fraction, water_viscosity,
                                                     oil_fraction, oil_viscosity, sicd);
        const EvalWell mixture_viscosity = (water_fraction + oil_fraction) * liquid_emulsion_viscosity + gas_fraction * gas_viscosity;

        const EvalWell reservoir_rate = segment_mass_rates_[seg] / density;

        const EvalWell reservoir_rate_icd = reservoir_rate * sicd.scalingFactor();

        const double viscosity_cali = sicd.viscosityCalibration();

        using MathTool = MathToolbox<EvalWell>;

        const double density_cali = sicd.densityCalibration();
        const EvalWell temp_value1 = MathTool::pow(density / density_cali, 0.75);
        const EvalWell temp_value2 = MathTool::pow(mixture_viscosity / viscosity_cali, 0.25);

        // formulation before 2016, base_strength is used
        // const double base_strength = sicd.strength() / density_cali;
        // formulation since 2016, strength is used instead
        const double strength = sicd.strength();

        const double sign = reservoir_rate_icd <= 0. ? 1.0 : -1.0;

        return sign * temp_value1 * temp_value2 * strength * reservoir_rate_icd * reservoir_rate_icd;
    }




    template<typename TypeTag>
    typename MultisegmentWell<TypeTag>::EvalWell
    MultisegmentWell<TypeTag>::
    pressureDropValve(const int seg) const
    {
        const Valve& valve = segmentSet()[seg].valve();

        const EvalWell& mass_rate = segment_mass_rates_[seg];
        const int seg_upwind = upwinding_segments_[seg];
        EvalWell visc = segment_viscosities_[seg_upwind];
        EvalWell density = segment_densities_[seg_upwind];
        // WARNING
        // We disregard the derivatives from the upwind density to make sure derivatives
        // wrt. to different segments dont get mixed.
        if (seg != seg_upwind) {
            visc.clearDerivatives();
            density.clearDerivatives();
        }

        const double additional_length = valve.pipeAdditionalLength();
        const double roughness = valve.pipeRoughness();
        const double diameter = valve.pipeDiameter();
        const double area = valve.pipeCrossArea();

        const EvalWell friction_pressure_loss =
            mswellhelpers::frictionPressureLoss(additional_length, diameter, area, roughness, density, mass_rate, visc);

        const double area_con = valve.conCrossArea();
        const double cv = valve.conFlowCoefficient();

        const EvalWell constriction_pressure_loss =
            mswellhelpers::valveContrictionPressureLoss(mass_rate, density, area_con, cv);

        const double sign = mass_rate <= 0. ? 1.0 : -1.0;
        return sign * (friction_pressure_loss + constriction_pressure_loss);
    }
    }