/////////////////////////////////////////////////////////////////////////////////
//
//  Copyright (C) 2017     Statoil ASA
//
//  ResInsight 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.
//
//  ResInsight 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 at <http://www.gnu.org/licenses/gpl.html>
//  for more details.
//
/////////////////////////////////////////////////////////////////////////////////

#include "RigAccWellFlowCalculator.h"

#include "RiaWellFlowDefines.h"

#include "RigActiveCellInfo.h"
#include "RigFlowDiagResults.h"
#include "RigMainGrid.h"
#include "RigSimWellData.h"
#include "RigSimulationWellCoordsAndMD.h"
#include "RigWellResultPoint.h"

#include <cmath>

//==================================================================================================
///
///
//==================================================================================================

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
size_t RigEclCellIndexCalculator::resultCellIndex( size_t gridIndex, size_t gridCellIndex ) const
{
    size_t reservoirCellIndex = m_mainGrid->reservoirCellIndexByGridAndGridLocalCellIndex( gridIndex, gridCellIndex );

    return m_activeCellInfo->cellResultIndex( reservoirCellIndex );
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RigEclCellIndexCalculator::isCellVisible( size_t gridIndex, size_t gridCellIndex ) const
{
    if ( !m_cellVisibilities ) return true;

    size_t reservoirCellIndex = m_mainGrid->reservoirCellIndexByGridAndGridLocalCellIndex( gridIndex, gridCellIndex );

    return ( *m_cellVisibilities )[reservoirCellIndex];
}

//==================================================================================================
///
///
//==================================================================================================

//--------------------------------------------------------------------------------------------------
///
/// The pipeBranchesWellResultPoints are describing the lines between the points, starting with the first line
//  and is thus expected to be one less than the number of centerline points
//--------------------------------------------------------------------------------------------------
RigAccWellFlowCalculator::RigAccWellFlowCalculator( const std::vector<std::vector<cvf::Vec3d>>&          pipeBranchesCLCoords,
                                                    const std::vector<std::vector<RigWellResultPoint>>&  pipeBranchesWellResultPoints,
                                                    const std::map<QString, const std::vector<double>*>& tracerCellFractionValues,
                                                    const RigEclCellIndexCalculator&                     cellIndexCalculator,
                                                    double                                               smallContribThreshold,
                                                    bool                                                 isProducer )
    : m_pipeBranchesCLCoords( pipeBranchesCLCoords )
    , m_pipeBranchesWellResultPoints( pipeBranchesWellResultPoints )
    , m_tracerCellFractionValues( &tracerCellFractionValues )
    , m_cellIndexCalculator( cellIndexCalculator )
    , m_smallContributionsThreshold( smallContribThreshold )
    , m_isProducer( isProducer )
    , m_useTotalWellPhaseRateOnly( false )
{
    CVF_ASSERT( !pipeBranchesWellResultPoints.empty() );

    m_connectionFlowPrBranch.resize( m_pipeBranchesWellResultPoints.size() );
    m_pseudoLengthFlowPrBranch.resize( m_pipeBranchesWellResultPoints.size() );

    for ( const auto& it : ( *m_tracerCellFractionValues ) )
        m_tracerNames.push_back( it.first );

    m_tracerNames.push_back( RiaDefines::reservoirTracerName() );

    initializePipeBranchesMeasuredDepths();
    calculateFlowData();
    sortTracers();
    groupSmallContributions();
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RigAccWellFlowCalculator::RigAccWellFlowCalculator( const std::vector<std::vector<cvf::Vec3d>>&         pipeBranchesCLCoords,
                                                    const std::vector<std::vector<RigWellResultPoint>>& pipeBranchesWellResultPoints,
                                                    double                                              smallContribThreshold )
    : m_pipeBranchesCLCoords( pipeBranchesCLCoords )
    , m_pipeBranchesWellResultPoints( pipeBranchesWellResultPoints )
    , m_tracerCellFractionValues( nullptr )
    , m_cellIndexCalculator( RigEclCellIndexCalculator( nullptr, nullptr, nullptr ) )
    , m_smallContributionsThreshold( smallContribThreshold )
    , m_isProducer( true )
    , m_useTotalWellPhaseRateOnly( false )
{
    CVF_ASSERT( !pipeBranchesWellResultPoints.empty() );

    m_connectionFlowPrBranch.resize( m_pipeBranchesWellResultPoints.size() );
    m_pseudoLengthFlowPrBranch.resize( m_pipeBranchesWellResultPoints.size() );

    if ( !m_useTotalWellPhaseRateOnly )
    {
        m_tracerNames.push_back( RIG_FLOW_OIL_NAME );
        m_tracerNames.push_back( RIG_FLOW_GAS_NAME );
        m_tracerNames.push_back( RIG_FLOW_WATER_NAME );
    }
    else
    {
        m_tracerNames.push_back( RIG_FLOW_TOTAL_NAME );
    }

    initializePipeBranchesMeasuredDepths();
    calculateFlowData();

    if ( !m_useTotalWellPhaseRateOnly ) sortTracers();
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RigAccWellFlowCalculator::RigAccWellFlowCalculator( const std::vector<cvf::Vec3d>&         pipeBranchCLCoords,
                                                    const std::vector<RigWellResultPoint>& pipeBranchesWellResultPoints,
                                                    const std::vector<double>&             pipeBranchMeasuredDepths,
                                                    bool                                   totalFlowOnly )
    : m_tracerCellFractionValues( nullptr )
    , m_cellIndexCalculator( RigEclCellIndexCalculator( nullptr, nullptr, nullptr ) )
    , m_smallContributionsThreshold( 0.0 )
    , m_isProducer( true )
    , m_useTotalWellPhaseRateOnly( totalFlowOnly )
{
    CVF_ASSERT( !pipeBranchesWellResultPoints.empty() );

    m_pipeBranchesCLCoords.push_back( pipeBranchCLCoords );
    m_pipeBranchesWellResultPoints.push_back( pipeBranchesWellResultPoints );
    m_pipeBranchesMeasuredDepths.push_back( pipeBranchMeasuredDepths );

    m_connectionFlowPrBranch.resize( m_pipeBranchesWellResultPoints.size() );
    m_pseudoLengthFlowPrBranch.resize( m_pipeBranchesWellResultPoints.size() );

    if ( !m_useTotalWellPhaseRateOnly )
    {
        m_tracerNames.push_back( RIG_FLOW_OIL_NAME );
        m_tracerNames.push_back( RIG_FLOW_GAS_NAME );
        m_tracerNames.push_back( RIG_FLOW_WATER_NAME );
    }
    else
    {
        m_tracerNames.push_back( RIG_FLOW_TOTAL_NAME );
    }

    initializePipeBranchesMeasuredDepths();
    calculateFlowData();

    if ( !m_useTotalWellPhaseRateOnly ) sortTracers();
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigAccWellFlowCalculator::initializePipeBranchesMeasuredDepths()
{
    for ( const auto& branchClPoints : m_pipeBranchesCLCoords )
    {
        RigSimulationWellCoordsAndMD mdCalculator( branchClPoints );
        m_pipeBranchesMeasuredDepths.push_back( mdCalculator.measuredDepths() );
    }
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const std::vector<double>& RigAccWellFlowCalculator::connectionNumbersFromTop( size_t branchIdx ) const
{
    return m_connectionFlowPrBranch[branchIdx].depthValuesFromTop;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const std::vector<double>& RigAccWellFlowCalculator::accumulatedTracerFlowPrConnection( const QString& tracerName, size_t branchIdx ) const
{
    auto flowPrTracerIt = m_connectionFlowPrBranch[branchIdx].accFlowPrTracer.find( tracerName );
    if ( flowPrTracerIt != m_connectionFlowPrBranch[branchIdx].accFlowPrTracer.end() )
    {
        return flowPrTracerIt->second;
    }
    else
    {
        CVF_ASSERT( false );
        static std::vector<double> dummy;
        return dummy;
    }
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const std::vector<double>& RigAccWellFlowCalculator::tracerFlowPrConnection( const QString& tracerName, size_t branchIdx ) const
{
    auto flowPrTracerIt = m_connectionFlowPrBranch[branchIdx].flowPrTracer.find( tracerName );
    if ( flowPrTracerIt != m_connectionFlowPrBranch[branchIdx].flowPrTracer.end() )
    {
        return flowPrTracerIt->second;
    }
    else
    {
        CVF_ASSERT( false );
        static std::vector<double> dummy;
        return dummy;
    }
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const std::vector<double>& RigAccWellFlowCalculator::pseudoLengthFromTop( size_t branchIdx ) const
{
    return m_pseudoLengthFlowPrBranch[branchIdx].depthValuesFromTop;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const std::vector<double>& RigAccWellFlowCalculator::trueVerticalDepth( size_t branchIdx ) const
{
    return m_pseudoLengthFlowPrBranch[branchIdx].trueVerticalDepth;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const std::vector<double>& RigAccWellFlowCalculator::accumulatedTracerFlowPrPseudoLength( const QString& tracerName, size_t branchIdx ) const
{
    auto flowPrTracerIt = m_pseudoLengthFlowPrBranch[branchIdx].accFlowPrTracer.find( tracerName );
    if ( flowPrTracerIt != m_pseudoLengthFlowPrBranch[branchIdx].accFlowPrTracer.end() )
    {
        return flowPrTracerIt->second;
    }
    else
    {
        CVF_ASSERT( false );
        static std::vector<double> dummy;
        return dummy;
    }
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const std::vector<double>& RigAccWellFlowCalculator::tracerFlowPrPseudoLength( const QString& tracerName, size_t branchIdx ) const
{
    auto flowPrTracerIt = m_pseudoLengthFlowPrBranch[branchIdx].flowPrTracer.find( tracerName );
    if ( flowPrTracerIt != m_pseudoLengthFlowPrBranch[branchIdx].flowPrTracer.end() )
    {
        return flowPrTracerIt->second;
    }
    else
    {
        CVF_ASSERT( false );
        static std::vector<double> dummy;
        return dummy;
    }
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<std::pair<QString, double>> RigAccWellFlowCalculator::totalWellFlowPrTracer() const
{
    std::vector<QString>                    tracerNames = this->tracerNames();
    std::vector<std::pair<QString, double>> tracerWithValues;

    for ( const QString& tracerName : tracerNames )
    {
        const std::vector<double>& accFlow = this->accumulatedTracerFlowPrConnection( tracerName, 0 );
        tracerWithValues.push_back( std::make_pair( tracerName, accFlow.back() ) );
    }

    return tracerWithValues;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<std::pair<QString, double>> RigAccWellFlowCalculator::totalTracerFractions() const
{
    std::vector<std::pair<QString, double>> totalFlows = totalWellFlowPrTracer();

    float sumTracerFlows = 0.0f;
    for ( const auto& tracerVal : totalFlows )
    {
        sumTracerFlows += tracerVal.second;
    }

    if ( sumTracerFlows == 0.0 ) totalFlows.clear();

    for ( auto& tracerPair : totalFlows )
    {
        tracerPair.second = tracerPair.second / sumTracerFlows;
    }

    return totalFlows;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RigAccWellFlowCalculator::isWellFlowConsistent() const
{
    bool isConsistent = true;
    for ( const std::vector<RigWellResultPoint>& branch : m_pipeBranchesWellResultPoints )
    {
        for ( const RigWellResultPoint& wrp : branch )
        {
            isConsistent = isFlowRateConsistent( wrp.flowRate() );

            if ( !isConsistent ) break;
        }
        if ( !isConsistent ) break;
    }
    return isConsistent;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------

std::vector<double> RigAccWellFlowCalculator::calculateAccumulatedFractions( const std::vector<double>& accumulatedFlowPrTracer ) const
{
    double totalFlow = 0.0;
    for ( double tracerFlow : accumulatedFlowPrTracer )
    {
        totalFlow += tracerFlow;
    }

    std::vector<double> flowFractionsPrTracer( accumulatedFlowPrTracer.size(), 0.0 );

    if ( totalFlow == 0.0 || !isFlowRateConsistent( totalFlow ) ) // If we have no accumulated flow, we set all the flow
                                                                  // associated to the last tracer, which is the reservoir
    {
        flowFractionsPrTracer.back() = 1.0;
        return flowFractionsPrTracer;
    }

    for ( size_t tIdx = 0; tIdx < accumulatedFlowPrTracer.size(); ++tIdx )
    {
        double tracerFlow           = accumulatedFlowPrTracer[tIdx];
        flowFractionsPrTracer[tIdx] = tracerFlow / totalFlow;
    }

    return flowFractionsPrTracer;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RigAccWellFlowCalculator::isConnectionFlowConsistent( const RigWellResultPoint& wellCell ) const
{
    if ( !m_tracerCellFractionValues ) return true; // No flow diagnostics.

    return isFlowRateConsistent( wellCell.flowRate() );
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RigAccWellFlowCalculator::isFlowRateConsistent( double flowRate ) const
{
    if ( !m_tracerCellFractionValues ) return true; // No flow diagnostics.

    return ( flowRate >= 0.0 && m_isProducer ) || ( flowRate <= 0.0 && !m_isProducer );
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigAccWellFlowCalculator::calculateFlowData()
{
    for ( size_t branchIndex = 0; branchIndex < m_connectionFlowPrBranch.size(); branchIndex++ )
    {
        const auto branch = m_connectionFlowPrBranch[branchIndex];
        if ( branch.depthValuesFromTop.empty() )
        {
            calculateAccumulatedFlowPrConnection( branchIndex, 0 );
        }
    }

    for ( size_t branchIndex = 0; branchIndex < m_pseudoLengthFlowPrBranch.size(); branchIndex++ )
    {
        const auto branch = m_pseudoLengthFlowPrBranch[branchIndex];
        if ( branch.depthValuesFromTop.empty() )
        {
            calculateFlowPrPseudoLength( branchIndex, 0.0 );
        }
    }
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigAccWellFlowCalculator::calculateAccumulatedFlowPrConnection( size_t branchIdx, size_t startConnectionNumberFromTop )
{
    const std::vector<RigWellResultPoint>& branchCells = m_pipeBranchesWellResultPoints[branchIdx];

    std::vector<size_t> resPointUniqueIndexFromBottom = wrpToUniqueWrpIndexFromBottom( branchCells );

    size_t prevConnIndx = -1;
    int    clSegIdx     = static_cast<int>( branchCells.size() ) - 1;

    std::vector<double> accFlowPrTracer( m_tracerNames.size(), 0.0 );

    while ( clSegIdx >= 0 )
    {
        // Skip point if referring to the same cell as the previous centerline segment did
        {
            if ( resPointUniqueIndexFromBottom[clSegIdx] == prevConnIndx )
            {
                --clSegIdx;
                continue;
            }

            prevConnIndx = resPointUniqueIndexFromBottom[clSegIdx];
        }

        // Accumulate the connection-cell's fraction flows

        const RigWellResultPoint& wellCell = branchCells[clSegIdx];

        std::vector<double> flowPrTracer = calculateWellCellFlowPrTracer( wellCell, accFlowPrTracer );

        addDownStreamBranchFlow( &accFlowPrTracer, flowPrTracer );

        if ( !isConnectionFlowConsistent( wellCell ) )
        {
            // Associate all the flow with the reservoir tracer for inconsistent flow direction
            flowPrTracer        = std::vector<double>( flowPrTracer.size(), 0.0 );
            flowPrTracer.back() = wellCell.flowRate();
        }

        // Add the total accumulated (fraction) flows from any branches connected to this cell

        size_t connNumFromTop = connectionIndexFromTop( resPointUniqueIndexFromBottom, clSegIdx ) + startConnectionNumberFromTop;

        std::vector<size_t> downStreamBranchIndices = findDownStreamBranchIdxs( branchCells[clSegIdx] );
        for ( size_t dsBidx : downStreamBranchIndices )
        {
            BranchFlow& downStreamBranchFlow = m_connectionFlowPrBranch[dsBidx];
            if ( dsBidx != branchIdx && downStreamBranchFlow.depthValuesFromTop.size() == 0 ) // Not this branch or
                                                                                              // already calculated
            {
                calculateAccumulatedFlowPrConnection( dsBidx, connNumFromTop );
                std::vector<double> accBranchFlowPrTracer = accumulatedDsBranchFlowPrTracer( downStreamBranchFlow );
                addDownStreamBranchFlow( &accFlowPrTracer, accBranchFlowPrTracer );
                if ( m_pipeBranchesWellResultPoints[dsBidx].size() <= 3 )
                {
                    // Short branch. Will not be visible. Show branch flow as addition to this connections direct flow
                    addDownStreamBranchFlow( &flowPrTracer, accBranchFlowPrTracer );
                }
            }
        }

        // Push back the accumulated result into the storage

        BranchFlow& branchFlow = m_connectionFlowPrBranch[branchIdx];

        storeFlowOnDepth( &branchFlow, connNumFromTop, accFlowPrTracer, flowPrTracer );

        --clSegIdx;
    }
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigAccWellFlowCalculator::calculateFlowPrPseudoLength( size_t branchIdx, double startPseudoLengthFromTop )
{
    const std::vector<RigWellResultPoint>& branchCells    = m_pipeBranchesWellResultPoints[branchIdx];
    const std::vector<cvf::Vec3d>&         branchClPoints = m_pipeBranchesCLCoords[branchIdx];
    const std::vector<double>&             branchMDs      = m_pipeBranchesMeasuredDepths[branchIdx];

    int clSegIdx = static_cast<int>( branchCells.size() ) - 1;

    std::vector<double> accFlowPrTracer( m_tracerNames.size(), 0.0 );

    BranchFlow& branchFlow = m_pseudoLengthFlowPrBranch[branchIdx];

    RigWellResultPoint previousResultPoint;

    while ( clSegIdx >= 0 )
    {
        int cellBottomPointIndex = -1;
        int cellUpperPointIndex  = -1;
        int currentSegmentIndex  = -1;

        // Find the complete cell span
        {
            cellBottomPointIndex = clSegIdx + 1;

            previousResultPoint = branchCells[clSegIdx];
            --clSegIdx;
            while ( clSegIdx >= 0 && previousResultPoint.isEqual( branchCells[clSegIdx] ) )
            {
                --clSegIdx;
            }

            cellUpperPointIndex = clSegIdx + 1;
            currentSegmentIndex = cellUpperPointIndex;
        }
        const RigWellResultPoint& wellCell                 = branchCells[currentSegmentIndex];
        std::vector<double>       flowPrTracerToAccumulate = calculateWellCellFlowPrTracer( wellCell, accFlowPrTracer );

        double pseudoLengthFromTop_lower = branchMDs[cellBottomPointIndex] + startPseudoLengthFromTop;
        double tvd_lower                 = -branchClPoints[cellBottomPointIndex][2];

        // Push back the new start-of-cell flow, with the previously accumulated result into the storage

        std::vector<double> flowPrTracer;
        if ( !isConnectionFlowConsistent( wellCell ) )
        {
            // Associate all the flow with the reservoir tracer for inconsistent flow direction
            flowPrTracer        = std::vector<double>( flowPrTracerToAccumulate.size(), 0.0 );
            flowPrTracer.back() = wellCell.flowRate();
        }
        else
        {
            flowPrTracer = flowPrTracerToAccumulate;
        }

        storeFlowOnDepthWTvd( &branchFlow, pseudoLengthFromTop_lower, tvd_lower, accFlowPrTracer, flowPrTracer );

        // Accumulate the connection-cell's fraction flows

        addDownStreamBranchFlow( &accFlowPrTracer, flowPrTracerToAccumulate );

        double pseudoLengthFromTop_upper = branchMDs[cellUpperPointIndex] + startPseudoLengthFromTop;
        double tvd_upper                 = -branchClPoints[cellUpperPointIndex][2];

        // Push back the accumulated result into the storage

        storeFlowOnDepthWTvd( &branchFlow, pseudoLengthFromTop_upper, tvd_upper, accFlowPrTracer, flowPrTracer );

        // Add the total accumulated (fraction) flows from any branches connected to this cell

        std::vector<size_t> downStreamBranchIndices = findDownStreamBranchIdxs( branchCells[cellUpperPointIndex] );
        for ( size_t dsBidx : downStreamBranchIndices )
        {
            BranchFlow& downStreamBranchFlow = m_pseudoLengthFlowPrBranch[dsBidx];
            if ( dsBidx != branchIdx && downStreamBranchFlow.depthValuesFromTop.size() == 0 ) // Not this branch or
                                                                                              // already calculated
            {
                calculateFlowPrPseudoLength( dsBidx, pseudoLengthFromTop_upper );
                std::vector<double> accBranchFlowPrTracer = accumulatedDsBranchFlowPrTracer( downStreamBranchFlow );
                addDownStreamBranchFlow( &accFlowPrTracer, accBranchFlowPrTracer );
                if ( m_pipeBranchesWellResultPoints[dsBidx].size() <= 3 )
                {
                    // Short branch. Will not be visible. Show branch flow as addition to this connections direct flow
                    addDownStreamBranchFlow( &flowPrTracer, accBranchFlowPrTracer );
                }
            }
        }

        // Push back the accumulated result after adding the branch result into the storage

        if ( downStreamBranchIndices.size() )
            storeFlowOnDepthWTvd( &branchFlow, pseudoLengthFromTop_upper, tvd_upper, accFlowPrTracer, flowPrTracer );
    }
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigAccWellFlowCalculator::addDownStreamBranchFlow( std::vector<double>*       accFlowPrTracer,
                                                        const std::vector<double>& accBranchFlowPrTracer ) const
{
    double totalThisBranchFlow = 0.0;
    for ( double tracerFlow : *accFlowPrTracer )
    {
        totalThisBranchFlow += tracerFlow;
    }

    double totalDsBranchFlow = 0.0;
    for ( double tracerFlow : accBranchFlowPrTracer )
    {
        totalDsBranchFlow += tracerFlow;
    }

    bool isAccumulationConsistent = isFlowRateConsistent( totalThisBranchFlow ); // If inconsistent, is it always only
                                                                                 // the Reservoir tracer that has the
                                                                                 // flow ?
    bool isBranchConsistent = isFlowRateConsistent( totalDsBranchFlow );

    if ( isAccumulationConsistent == isBranchConsistent )
    {
        for ( size_t tracerIdx = 0; tracerIdx < ( *accFlowPrTracer ).size(); ++tracerIdx )
        {
            ( *accFlowPrTracer )[tracerIdx] += accBranchFlowPrTracer[tracerIdx];
        }
        return;
    }

    double totalAccFlow = totalThisBranchFlow + totalDsBranchFlow;

    if ( !isFlowRateConsistent( totalAccFlow ) )
    {
        // Reset the accumulated values, as everything must be moved to the "Reservoir" tracer.
        for ( double& val : ( *accFlowPrTracer ) )
            val = 0.0;

        // Put all flow into the Reservoir tracer
        accFlowPrTracer->back() = totalThisBranchFlow + totalDsBranchFlow;

        return;
    }

    // We will end up with a consistent accumulated flow, and need to keep the accumulated distribution in this branch
    // or to use the ds branch distribution

    std::vector<double> accFractionsPrTracer;

    if ( !isAccumulationConsistent && isBranchConsistent )
    {
        accFractionsPrTracer = calculateAccumulatedFractions( accBranchFlowPrTracer );
    }
    else if ( isAccumulationConsistent && !isBranchConsistent )
    {
        accFractionsPrTracer = calculateAccumulatedFractions( *accFlowPrTracer );
    }

    // Set the accumulated values to the totalFlow times the tracer fraction selected.

    for ( size_t tIdx = 0; tIdx < accFlowPrTracer->size(); ++tIdx )
    {
        ( *accFlowPrTracer )[tIdx] = accFractionsPrTracer[tIdx] * ( totalAccFlow );
    }
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigAccWellFlowCalculator::storeFlowOnDepth( BranchFlow*                branchFlow,
                                                 double                     depthValue,
                                                 const std::vector<double>& accFlowPrTracer,
                                                 const std::vector<double>& flowPrTracer )
{
    size_t tracerIdx = 0;
    for ( const auto& tracerName : m_tracerNames )
    {
        branchFlow->accFlowPrTracer[tracerName].push_back( accFlowPrTracer[tracerIdx] );
        branchFlow->flowPrTracer[tracerName].push_back( flowPrTracer[tracerIdx] );
        tracerIdx++;
    }

    branchFlow->depthValuesFromTop.push_back( depthValue );
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigAccWellFlowCalculator::storeFlowOnDepthWTvd( BranchFlow*                branchFlow,
                                                     double                     depthValue,
                                                     double                     trueVerticalDepth,
                                                     const std::vector<double>& accFlowPrTracer,
                                                     const std::vector<double>& flowPrTracer )
{
    size_t tracerIdx = 0;
    for ( const auto& tracerName : m_tracerNames )
    {
        branchFlow->accFlowPrTracer[tracerName].push_back( accFlowPrTracer[tracerIdx] );
        branchFlow->flowPrTracer[tracerName].push_back( flowPrTracer[tracerIdx] );
        tracerIdx++;
    }

    branchFlow->depthValuesFromTop.push_back( depthValue );
    branchFlow->trueVerticalDepth.push_back( trueVerticalDepth );
}

std::vector<double> RigAccWellFlowCalculator::accumulatedDsBranchFlowPrTracer( const BranchFlow& downStreamBranchFlow ) const
{
    std::vector<double> accBranchFlowPrTracer( m_tracerNames.size(), 0.0 );

    size_t tracerIdx = 0;
    for ( const auto& tracerName : m_tracerNames )
    {
        const auto trNameAccFlowsPair = downStreamBranchFlow.accFlowPrTracer.find( tracerName );
        if ( trNameAccFlowsPair != downStreamBranchFlow.accFlowPrTracer.end() )
        {
            accBranchFlowPrTracer[tracerIdx] = trNameAccFlowsPair->second.back();
        }
        tracerIdx++;
    }

    return accBranchFlowPrTracer;
}

//--------------------------------------------------------------------------------------------------
/// Calculate the flow pr tracer. If inconsistent flow, keep the existing fractions constant
//--------------------------------------------------------------------------------------------------
std::vector<double> RigAccWellFlowCalculator::calculateWellCellFlowPrTracer( const RigWellResultPoint& wellCell,
                                                                             const std::vector<double>& currentAccumulatedFlowPrTracer ) const
{
    std::vector<double> flowPrTracer( m_tracerNames.size(), 0.0 );

    if ( wellCell.isCell() && !m_cellIndexCalculator.isCellVisible( wellCell.gridIndex(), wellCell.cellIndex() ) ) return flowPrTracer;

    if ( !isConnectionFlowConsistent( wellCell ) )
    {
        double flowRate = wellCell.flowRate();
        flowPrTracer    = calculateAccumulatedFractions( currentAccumulatedFlowPrTracer );
        for ( double& accFraction : flowPrTracer )
        {
            accFraction *= flowRate;
        }

        return flowPrTracer;
    }

    if ( m_tracerCellFractionValues )
    {
        if ( wellCell.isCell() && wellCell.isOpen() )
        {
            size_t resCellIndex              = m_cellIndexCalculator.resultCellIndex( wellCell.gridIndex(), wellCell.cellIndex() );
            size_t tracerIdx                 = 0;
            double totalTracerFractionInCell = 0.0;
            for ( const auto& tracerFractionValsPair : ( *m_tracerCellFractionValues ) )
            {
                const std::vector<double>* fractionVals = tracerFractionValsPair.second;
                if ( fractionVals )
                {
                    double cellTracerFraction = ( *fractionVals )[resCellIndex];
                    if ( cellTracerFraction != HUGE_VAL && cellTracerFraction == cellTracerFraction )
                    {
                        double tracerFlow       = cellTracerFraction * wellCell.flowRate();
                        flowPrTracer[tracerIdx] = tracerFlow;

                        totalTracerFractionInCell += cellTracerFraction;
                    }
                }
                tracerIdx++;
            }

            double reservoirFraction   = 1.0 - totalTracerFractionInCell;
            double reservoirTracerFlow = reservoirFraction * wellCell.flowRate();
            flowPrTracer[tracerIdx]    = reservoirTracerFlow;
        }
    }
    else
    {
        if ( !m_useTotalWellPhaseRateOnly )
        {
            flowPrTracer[0] = wellCell.oilRate();
            flowPrTracer[1] = wellCell.gasRate();
            flowPrTracer[2] = wellCell.waterRate();
        }
        else
        {
            flowPrTracer[0] = wellCell.flowRate();
        }
    }

    return flowPrTracer;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<size_t> RigAccWellFlowCalculator::wrpToUniqueWrpIndexFromBottom( const std::vector<RigWellResultPoint>& branchCells ) const
{
    std::vector<size_t> resPointToConnectionIndexFromBottom;
    resPointToConnectionIndexFromBottom.resize( branchCells.size(), -1 );

    size_t connIdxFromBottom = 0;
    int    clSegIdx          = static_cast<int>( branchCells.size() ) - 1;

    if ( clSegIdx < 0 ) return resPointToConnectionIndexFromBottom;

    size_t prevGridIdx     = branchCells[clSegIdx].gridIndex();
    size_t prevGridCellIdx = branchCells[clSegIdx].cellIndex();
    int    prevErtSegId    = branchCells[clSegIdx].segmentId();
    int    prevErtBranchId = branchCells[clSegIdx].branchId();

    while ( clSegIdx >= 0 )
    {
        if ( branchCells[clSegIdx].isValid() &&
             ( branchCells[clSegIdx].gridIndex() != prevGridIdx || branchCells[clSegIdx].cellIndex() != prevGridCellIdx ||
               branchCells[clSegIdx].segmentId() != prevErtSegId || branchCells[clSegIdx].branchId() != prevErtBranchId ) )
        {
            ++connIdxFromBottom;

            prevGridIdx     = branchCells[clSegIdx].gridIndex();
            prevGridCellIdx = branchCells[clSegIdx].cellIndex();
            prevErtSegId    = branchCells[clSegIdx].segmentId();
            prevErtBranchId = branchCells[clSegIdx].branchId();
        }

        resPointToConnectionIndexFromBottom[clSegIdx] = connIdxFromBottom;

        --clSegIdx;
    }

    return resPointToConnectionIndexFromBottom;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
size_t RigAccWellFlowCalculator::connectionIndexFromTop( const std::vector<size_t>& resPointToConnectionIndexFromBottom, size_t clSegIdx )
{
    return resPointToConnectionIndexFromBottom.front() - resPointToConnectionIndexFromBottom[clSegIdx];
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<size_t> RigAccWellFlowCalculator::findDownStreamBranchIdxs( const RigWellResultPoint& connectionPoint ) const
{
    std::vector<size_t> downStreamBranchIdxs;

    for ( size_t bIdx = 0; bIdx < m_pipeBranchesWellResultPoints.size(); ++bIdx )
    {
        if ( m_pipeBranchesWellResultPoints[bIdx][0].gridIndex() == connectionPoint.gridIndex() &&
             m_pipeBranchesWellResultPoints[bIdx][0].cellIndex() == connectionPoint.cellIndex() &&
             m_pipeBranchesWellResultPoints[bIdx][0].branchId() == connectionPoint.branchId() &&
             m_pipeBranchesWellResultPoints[bIdx][0].segmentId() == connectionPoint.segmentId() )
        {
            downStreamBranchIdxs.push_back( bIdx );
        }
    }
    return downStreamBranchIdxs;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------

void RigAccWellFlowCalculator::sortTracers()
{
    std::multimap<double, QString> sortedTracers;
    for ( const QString& tracerName : m_tracerNames )
    {
        const std::vector<double>& mainBranchAccFlow = accumulatedTracerFlowPrConnection( tracerName, 0 );

        double totalFlow = 0.0;

        if ( mainBranchAccFlow.size() ) totalFlow = -fabs( mainBranchAccFlow.back() ); // Based on size in reverse order (biggest to least)

        sortedTracers.insert( { totalFlow, tracerName } );
    }

    m_tracerNames.clear();
    for ( const auto& tracerPair : sortedTracers )
    {
        m_tracerNames.push_back( tracerPair.second );
    }
}

//--------------------------------------------------------------------------------------------------
/// Concatenate small tracers into an "Other" group
//--------------------------------------------------------------------------------------------------
void RigAccWellFlowCalculator::groupSmallContributions()
{
    if ( !( m_smallContributionsThreshold > 0.0 ) ) return;

    // Find the tracers we need to group

    std::vector<QString> tracersToGroup;
    {
        bool hasConsistentWellFlow = isWellFlowConsistent();

        std::vector<std::pair<QString, double>> totalTracerFractions = this->totalTracerFractions();

        for ( const auto& tracerPair : totalTracerFractions )
        {
            if ( fabs( tracerPair.second ) <= m_smallContributionsThreshold &&
                 ( hasConsistentWellFlow || tracerPair.first != RiaDefines::reservoirTracerName() ) ) // Do not group
                                                                                                      // the Reservoir
                                                                                                      // tracer if the
                                                                                                      // well flow is
                                                                                                      // inconsistent,
                                                                                                      // because cross
                                                                                                      // flow is shown
                                                                                                      // as the
                                                                                                      // reservoir
                                                                                                      // fraction
            {
                tracersToGroup.push_back( tracerPair.first );
            }
        }
    }

    if ( tracersToGroup.size() < 2 ) return; // Must at least group two ...

    // Concatenate the values for each branch, erasing the tracers being grouped, replaced with the concatenated values

    for ( BranchFlow& brRes : m_connectionFlowPrBranch )
    {
        groupSmallTracers( &brRes.accFlowPrTracer, tracersToGroup );
        groupSmallTracers( &brRes.flowPrTracer, tracersToGroup );
    }

    for ( BranchFlow& brRes : m_pseudoLengthFlowPrBranch )
    {
        groupSmallTracers( &brRes.accFlowPrTracer, tracersToGroup );
        groupSmallTracers( &brRes.flowPrTracer, tracersToGroup );
    }

    // Remove the grouped tracer names from the tracerName list, and replace with the "Others" name

    std::vector<QString> filteredTracernames;
    for ( const QString& tracerName : m_tracerNames )
    {
        bool isDeleted = false;
        for ( const QString& deletedTracerName : tracersToGroup )
        {
            if ( tracerName == deletedTracerName )
            {
                isDeleted = true;
                break;
            }
        }

        if ( !isDeleted ) filteredTracernames.push_back( tracerName );
    }

    m_tracerNames.swap( filteredTracernames );
    m_tracerNames.push_back( RIG_TINY_TRACER_GROUP_NAME );
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigAccWellFlowCalculator::groupSmallTracers( std::map<QString, std::vector<double>>* branchFlowSet,
                                                  const std::vector<QString>&             tracersToGroup )
{
    if ( branchFlowSet->empty() ) return;

    size_t              depthCount = branchFlowSet->begin()->second.size();
    std::vector<double> groupedAccFlowValues( depthCount, 0.0 );

    for ( const QString& tracername : tracersToGroup )
    {
        auto it = branchFlowSet->find( tracername );

        if ( it != branchFlowSet->end() )
        {
            const std::vector<double>& tracerVals = it->second;
            for ( size_t cIdx = 0; cIdx < groupedAccFlowValues.size(); ++cIdx )
            {
                groupedAccFlowValues[cIdx] += tracerVals[cIdx];
            }
        }

        branchFlowSet->erase( it );
    }

    ( *branchFlowSet )[RIG_TINY_TRACER_GROUP_NAME] = groupedAccFlowValues;
}