diff --git a/opm/simulators/wells/StandardWellConnections.cpp b/opm/simulators/wells/StandardWellConnections.cpp index 1cc5f0ac3..a17ea877a 100644 --- a/opm/simulators/wells/StandardWellConnections.cpp +++ b/opm/simulators/wells/StandardWellConnections.cpp @@ -22,24 +22,32 @@ #include #include +#include + +#include +#include +#include + +#include #include #include +#include + namespace Opm { -template -StandardWellConnections:: +template +StandardWellConnections:: StandardWellConnections(const WellInterfaceGeneric& well) - : perf_densities_(well_.numPerfs()) - , perf_pressure_diffs_(well_.numPerfs()) + : perf_densities_(well.numPerfs()) + , perf_pressure_diffs_(well.numPerfs()) , well_(well) { } -template -void -StandardWellConnections:: +template +void StandardWellConnections:: computeConnectionPressureDelta() { // Algorithm: @@ -74,6 +82,225 @@ computeConnectionPressureDelta() well_.parallelWellInfo().partialSumPerfValues(beg, end); } -template class StandardWellConnections; +template +void StandardWellConnections:: +computeConnectionDensities(const std::vector& perfComponentRates, + const std::vector& b_perf, + const std::vector& rsmax_perf, + const std::vector& rvmax_perf, + const std::vector& rvwmax_perf, + const std::vector& surf_dens_perf, + DeferredLogger& deferred_logger) +{ + // Verify that we have consistent input. + const int nperf = well_.numPerfs(); + const int num_comp = well_.numComponents(); + + // 1. Compute the flow (in surface volume units for each + // component) exiting up the wellbore from each perforation, + // taking into account flow from lower in the well, and + // in/out-flow at each perforation. + std::vector q_out_perf((nperf)*num_comp, 0.0); + + // Step 1 depends on the order of the perforations. Hence we need to + // do the modifications globally. + // Create and get the global perforation information and do this sequentially + // on each process + + const auto& factory = well_.parallelWellInfo().getGlobalPerfContainerFactory(); + auto global_q_out_perf = factory.createGlobal(q_out_perf, num_comp); + auto global_perf_comp_rates = factory.createGlobal(perfComponentRates, num_comp); + + // TODO: investigate whether we should use the following techniques to calcuate the composition of flows in the wellbore + // Iterate over well perforations from bottom to top. + for (int perf = factory.numGlobalPerfs() - 1; perf >= 0; --perf) { + for (int component = 0; component < num_comp; ++component) { + auto index = perf * num_comp + component; + if (perf == factory.numGlobalPerfs() - 1) { + // This is the bottom perforation. No flow from below. + global_q_out_perf[index] = 0.0; + } else { + // Set equal to flow from below. + global_q_out_perf[index] = global_q_out_perf[index + num_comp]; + } + // Subtract outflow through perforation. + global_q_out_perf[index] -= global_perf_comp_rates[index]; + } + } + + // Copy the data back to local view + factory.copyGlobalToLocal(global_q_out_perf, q_out_perf, num_comp); + + // 2. Compute the component mix at each perforation as the + // absolute values of the surface rates divided by their sum. + // Then compute volume ratios (formation factors) for each perforation. + // Finally compute densities for the segments associated with each perforation. + std::vector mix(num_comp,0.0); + std::vector x(num_comp); + std::vector surf_dens(num_comp); + + for (int perf = 0; perf < nperf; ++perf) { + // Find component mix. + const Scalar tot_surf_rate = std::accumulate(q_out_perf.begin() + num_comp*perf, + q_out_perf.begin() + num_comp*(perf+1), 0.0); + if (tot_surf_rate != 0.0) { + for (int component = 0; component < num_comp; ++component) { + mix[component] = std::fabs(q_out_perf[perf*num_comp + component]/tot_surf_rate); + } + } else if (num_comp == 1) { + mix[num_comp-1] = 1.0; + } else { + std::fill(mix.begin(), mix.end(), 0.0); + // No flow => use well specified fractions for mix. + if (well_.isInjector()) { + switch (well_.wellEcl().injectorType()) { + case InjectorType::WATER: + mix[FluidSystem::waterCompIdx] = 1.0; + break; + case InjectorType::GAS: + mix[FluidSystem::gasCompIdx] = 1.0; + break; + case InjectorType::OIL: + mix[FluidSystem::oilCompIdx] = 1.0; + break; + case InjectorType::MULTI: + // Not supported. + // deferred_logger.warning("MULTI_PHASE_INJECTOR_NOT_SUPPORTED", + // "Multi phase injectors are not supported, requested for well " + name()); + break; + } + } else { + assert(well_.isProducer()); + // For the frist perforation without flow we use the preferred phase to decide the mix initialization. + if (perf == 0) { // + switch (well_.wellEcl().getPreferredPhase()) { + case Phase::OIL: + mix[FluidSystem::oilCompIdx] = 1.0; + break; + case Phase::GAS: + mix[FluidSystem::gasCompIdx] = 1.0; + break; + case Phase::WATER: + mix[FluidSystem::waterCompIdx] = 1.0; + break; + default: + // No others supported. + break; + } + // For the rest of the perforation without flow we use mix from the above perforation. + } else { + mix = x; + } + + } + } + // Compute volume ratio. + x = mix; + + // Subtract dissolved gas from oil phase and vapporized oil from gas phase and vaporized water from gas phase + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { + const unsigned gaspos = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + const unsigned oilpos = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); + Scalar rs = 0.0; + Scalar rv = 0.0; + if (!rsmax_perf.empty() && mix[oilpos] > 1e-12) { + rs = std::min(mix[gaspos]/mix[oilpos], rsmax_perf[perf]); + } + if (!rvmax_perf.empty() && mix[gaspos] > 1e-12) { + rv = std::min(mix[oilpos]/mix[gaspos], rvmax_perf[perf]); + } + const Scalar d = 1.0 - rs*rv; + if (d <= 0.0) { + std::ostringstream sstr; + sstr << "Problematic d value " << d << " obtained for well " << well_.name() + << " during ccomputeConnectionDensities with rs " << rs + << ", rv " << rv + << " obtaining d " << d + << " Continue as if no dissolution (rs = 0) and vaporization (rv = 0) " + << " for this connection."; + deferred_logger.debug(sstr.str()); + } else { + if (rs > 0.0) { + // Subtract gas in oil from gas mixture + x[gaspos] = (mix[gaspos] - mix[oilpos]*rs)/d; + } + if (rv > 0.0) { + // Subtract oil in gas from oil mixture + x[oilpos] = (mix[oilpos] - mix[gaspos]*rv)/d; + } + } + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + //matrix system: (mix[oilpos] = q_os, x[oilpos] = bo*q_or, etc...) + //┌ ┐ ┌ ┐ ┌ ┐ + //│mix[oilpos] │ | 1 Rv 0 | |x[oilpos] | + //│mix[gaspos] │ = │ Rs 1 0 │ │x[gaspos] │ + //│mix[waterpos]│ │ 0 Rvw 1 │ │x[waterpos │ + //└ ┘ └ ┘ └ ┘ + const unsigned waterpos = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); + Scalar rvw = 0.0; + if (!rvwmax_perf.empty() && mix[gaspos] > 1e-12) { + rvw = std::min(mix[waterpos]/mix[gaspos], rvwmax_perf[perf]); + } + if (rvw > 0.0) { + // Subtract water in gas from water mixture + x[waterpos] = mix[waterpos] - x[gaspos] * rvw; + } + } + } else if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + //no oil + const unsigned gaspos = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + const unsigned waterpos = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); + Scalar rvw = 0.0; + if (!rvwmax_perf.empty() && mix[gaspos] > 1e-12) { + rvw = std::min(mix[waterpos]/mix[gaspos], rvwmax_perf[perf]); + } + if (rvw > 0.0) { + // Subtract water in gas from water mixture + x[waterpos] = mix[waterpos] - mix[gaspos] * rvw; + } + } + + Scalar volrat = 0.0; + for (int component = 0; component < num_comp; ++component) { + volrat += x[component] / b_perf[perf*num_comp+ component]; + } + for (int component = 0; component < num_comp; ++component) { + surf_dens[component] = surf_dens_perf[perf*num_comp+ component]; + } + + // Compute segment density. + perf_densities_[perf] = std::inner_product(surf_dens.begin(), surf_dens.end(), mix.begin(), 0.0) / volrat; + } +} + +#define INSTANCE(...) \ +template class StandardWellConnections,__VA_ARGS__,double>; + +// One phase +INSTANCE(BlackOilOnePhaseIndices<0u,0u,0u,0u,false,false,0u,1u,0u>) +INSTANCE(BlackOilOnePhaseIndices<0u,0u,0u,1u,false,false,0u,1u,0u>) +INSTANCE(BlackOilOnePhaseIndices<0u,0u,0u,0u,false,false,0u,1u,5u>) + +// Two phase +INSTANCE(BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,0u,0u>) +INSTANCE(BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,1u,0u>) +INSTANCE(BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,2u,0u>) +INSTANCE(BlackOilTwoPhaseIndices<0u,0u,1u,0u,false,false,0u,2u,0u>) +INSTANCE(BlackOilTwoPhaseIndices<0u,0u,1u,0u,false,true,0u,2u,0u>) +INSTANCE(BlackOilTwoPhaseIndices<0u,0u,0u,1u,false,false,0u,1u,0u>) +INSTANCE(BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,true,0u,0u,0u>) +INSTANCE(BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,true,0u,2u,0u>) +INSTANCE(BlackOilTwoPhaseIndices<0u,0u,2u,0u,false,false,0u,2u,0u>) + +// Blackoil +INSTANCE(BlackOilIndices<0u,0u,0u,0u,false,false,0u,0u>) +INSTANCE(BlackOilIndices<1u,0u,0u,0u,false,false,0u,0u>) +INSTANCE(BlackOilIndices<0u,1u,0u,0u,false,false,0u,0u>) +INSTANCE(BlackOilIndices<0u,0u,1u,0u,false,false,0u,0u>) +INSTANCE(BlackOilIndices<0u,0u,0u,1u,false,false,0u,0u>) +INSTANCE(BlackOilIndices<0u,0u,0u,0u,true,false,0u,0u>) +INSTANCE(BlackOilIndices<0u,0u,0u,0u,false,true,0u,0u>) +INSTANCE(BlackOilIndices<0u,0u,0u,1u,false,true,0u,0u>) +INSTANCE(BlackOilIndices<0u,0u,0u,1u,false,false,1u,0u>) } diff --git a/opm/simulators/wells/StandardWellConnections.hpp b/opm/simulators/wells/StandardWellConnections.hpp index 0567706f9..f08738136 100644 --- a/opm/simulators/wells/StandardWellConnections.hpp +++ b/opm/simulators/wells/StandardWellConnections.hpp @@ -28,9 +28,10 @@ namespace Opm { +class DeferredLogger; class WellInterfaceGeneric; -template +template class StandardWellConnections { public: @@ -38,6 +39,16 @@ public: void computeConnectionPressureDelta(); + // TODO: not total sure whether it is a good idea to put this function here + // the major reason to put here is to avoid the usage of Wells struct + void computeConnectionDensities(const std::vector& perfComponentRates, + const std::vector& b_perf, + const std::vector& rsmax_perf, + const std::vector& rvmax_perf, + const std::vector& rvwmax_perf, + const std::vector& surf_dens_perf, + DeferredLogger& deferred_logger); + Scalar getRho() const { return this->perf_densities_.empty() ? 0.0 : perf_densities_[0]; diff --git a/opm/simulators/wells/StandardWellEval.cpp b/opm/simulators/wells/StandardWellEval.cpp index eb5e57555..5da334d0f 100644 --- a/opm/simulators/wells/StandardWellEval.cpp +++ b/opm/simulators/wells/StandardWellEval.cpp @@ -155,199 +155,6 @@ getWellConvergence(const WellState& well_state, return report; } -template -void -StandardWellEval:: -computeConnectionDensities(const std::vector& perfComponentRates, - const std::vector& b_perf, - const std::vector& rsmax_perf, - const std::vector& rvmax_perf, - const std::vector& rvwmax_perf, - const std::vector& surf_dens_perf, - DeferredLogger& deferred_logger) -{ - // Verify that we have consistent input. - const int nperf = baseif_.numPerfs(); - const int num_comp = baseif_.numComponents(); - - // 1. Compute the flow (in surface volume units for each - // component) exiting up the wellbore from each perforation, - // taking into account flow from lower in the well, and - // in/out-flow at each perforation. - std::vector q_out_perf((nperf)*num_comp, 0.0); - - // Step 1 depends on the order of the perforations. Hence we need to - // do the modifications globally. - // Create and get the global perforation information and do this sequentially - // on each process - - const auto& factory = baseif_.parallelWellInfo().getGlobalPerfContainerFactory(); - auto global_q_out_perf = factory.createGlobal(q_out_perf, num_comp); - auto global_perf_comp_rates = factory.createGlobal(perfComponentRates, num_comp); - - // TODO: investigate whether we should use the following techniques to calcuate the composition of flows in the wellbore - // Iterate over well perforations from bottom to top. - for (int perf = factory.numGlobalPerfs() - 1; perf >= 0; --perf) { - for (int component = 0; component < num_comp; ++component) { - auto index = perf * num_comp + component; - if (perf == factory.numGlobalPerfs() - 1) { - // This is the bottom perforation. No flow from below. - global_q_out_perf[index] = 0.0; - } else { - // Set equal to flow from below. - global_q_out_perf[index] = global_q_out_perf[index + num_comp]; - } - // Subtract outflow through perforation. - global_q_out_perf[index] -= global_perf_comp_rates[index]; - } - } - - // Copy the data back to local view - factory.copyGlobalToLocal(global_q_out_perf, q_out_perf, num_comp); - - // 2. Compute the component mix at each perforation as the - // absolute values of the surface rates divided by their sum. - // Then compute volume ratios (formation factors) for each perforation. - // Finally compute densities for the segments associated with each perforation. - std::vector mix(num_comp,0.0); - std::vector x(num_comp); - std::vector surf_dens(num_comp); - - for (int perf = 0; perf < nperf; ++perf) { - // Find component mix. - const double tot_surf_rate = std::accumulate(q_out_perf.begin() + num_comp*perf, - q_out_perf.begin() + num_comp*(perf+1), 0.0); - if (tot_surf_rate != 0.0) { - for (int component = 0; component < num_comp; ++component) { - mix[component] = std::fabs(q_out_perf[perf*num_comp + component]/tot_surf_rate); - } - } else if (num_comp == 1) { - mix[num_comp-1] = 1.0; - } else { - std::fill(mix.begin(), mix.end(), 0.0); - // No flow => use well specified fractions for mix. - if (baseif_.isInjector()) { - switch (baseif_.wellEcl().injectorType()) { - case InjectorType::WATER: - mix[FluidSystem::waterCompIdx] = 1.0; - break; - case InjectorType::GAS: - mix[FluidSystem::gasCompIdx] = 1.0; - break; - case InjectorType::OIL: - mix[FluidSystem::oilCompIdx] = 1.0; - break; - case InjectorType::MULTI: - // Not supported. - // deferred_logger.warning("MULTI_PHASE_INJECTOR_NOT_SUPPORTED", - // "Multi phase injectors are not supported, requested for well " + name()); - break; - } - } else { - assert(baseif_.isProducer()); - // For the frist perforation without flow we use the preferred phase to decide the mix initialization. - if (perf == 0) { // - switch (baseif_.wellEcl().getPreferredPhase()) { - case Phase::OIL: - mix[FluidSystem::oilCompIdx] = 1.0; - break; - case Phase::GAS: - mix[FluidSystem::gasCompIdx] = 1.0; - break; - case Phase::WATER: - mix[FluidSystem::waterCompIdx] = 1.0; - break; - default: - // No others supported. - break; - } - // For the rest of the perforation without flow we use mix from the above perforation. - } else { - mix = x; - } - - } - } - // Compute volume ratio. - x = mix; - - // Subtract dissolved gas from oil phase and vapporized oil from gas phase and vaporized water from gas phase - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { - const unsigned gaspos = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); - const unsigned oilpos = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); - double rs = 0.0; - double rv = 0.0; - if (!rsmax_perf.empty() && mix[oilpos] > 1e-12) { - rs = std::min(mix[gaspos]/mix[oilpos], rsmax_perf[perf]); - } - if (!rvmax_perf.empty() && mix[gaspos] > 1e-12) { - rv = std::min(mix[oilpos]/mix[gaspos], rvmax_perf[perf]); - } - const double d = 1.0 - rs*rv; - if (d <= 0.0) { - std::ostringstream sstr; - sstr << "Problematic d value " << d << " obtained for well " << baseif_.name() - << " during ccomputeConnectionDensities with rs " << rs - << ", rv " << rv - << " obtaining d " << d - << " Continue as if no dissolution (rs = 0) and vaporization (rv = 0) " - << " for this connection."; - deferred_logger.debug(sstr.str()); - } else { - if (rs > 0.0) { - // Subtract gas in oil from gas mixture - x[gaspos] = (mix[gaspos] - mix[oilpos]*rs)/d; - } - if (rv > 0.0) { - // Subtract oil in gas from oil mixture - x[oilpos] = (mix[oilpos] - mix[gaspos]*rv)/d; - } - } - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - //matrix system: (mix[oilpos] = q_os, x[oilpos] = bo*q_or, etc...) - //┌ ┐ ┌ ┐ ┌ ┐ - //│mix[oilpos] │ | 1 Rv 0 | |x[oilpos] | - //│mix[gaspos] │ = │ Rs 1 0 │ │x[gaspos] │ - //│mix[waterpos]│ │ 0 Rvw 1 │ │x[waterpos │ - //└ ┘ └ ┘ └ ┘ - const unsigned waterpos = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); - double rvw = 0.0; - if (!rvwmax_perf.empty() && mix[gaspos] > 1e-12) { - rvw = std::min(mix[waterpos]/mix[gaspos], rvwmax_perf[perf]); - } - if (rvw > 0.0) { - // Subtract water in gas from water mixture - x[waterpos] = mix[waterpos] - x[gaspos] * rvw; - } - } - } else if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - //no oil - const unsigned gaspos = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); - const unsigned waterpos = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); - double rvw = 0.0; - if (!rvwmax_perf.empty() && mix[gaspos] > 1e-12) { - rvw = std::min(mix[waterpos]/mix[gaspos], rvwmax_perf[perf]); - } - if (rvw > 0.0) { - // Subtract water in gas from water mixture - x[waterpos] = mix[waterpos] - mix[gaspos] * rvw; - } - } - - double volrat = 0.0; - for (int component = 0; component < num_comp; ++component) { - volrat += x[component] / b_perf[perf*num_comp+ component]; - } - for (int component = 0; component < num_comp; ++component) { - surf_dens[component] = surf_dens_perf[perf*num_comp+ component]; - } - - // Compute segment density. - this->connections_.perf_densities_[perf] = std::inner_product(surf_dens.begin(), surf_dens.end(), mix.begin(), 0.0) / volrat; - } -} - - template void StandardWellEval:: diff --git a/opm/simulators/wells/StandardWellEval.hpp b/opm/simulators/wells/StandardWellEval.hpp index 47570068e..97afaa480 100644 --- a/opm/simulators/wells/StandardWellEval.hpp +++ b/opm/simulators/wells/StandardWellEval.hpp @@ -78,16 +78,6 @@ protected: // computing the accumulation term for later use in well mass equations void computeAccumWell(); - // TODO: not total sure whether it is a good idea to put this function here - // the major reason to put here is to avoid the usage of Wells struct - void computeConnectionDensities(const std::vector& perfComponentRates, - const std::vector& b_perf, - const std::vector& rsmax_perf, - const std::vector& rvmax_perf, - const std::vector& rvwmax_perf, - const std::vector& surf_dens_perf, - DeferredLogger& deferred_logger); - ConvergenceReport getWellConvergence(const WellState& well_state, const std::vector& B_avg, const double maxResidualAllowed, @@ -111,7 +101,7 @@ protected: std::vector F0_; StandardWellEquations linSys_; //!< Linear equation system - StandardWellConnections connections_; //!< Connection level values + StandardWellConnections connections_; //!< Connection level values }; } diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index e1ff6ca52..f1e264c54 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -1625,8 +1625,7 @@ namespace Opm } } - this->computeConnectionDensities(perfRates, b_perf, rsmax_perf, rvmax_perf, rvwmax_perf, surf_dens_perf, deferred_logger); - + this->connections_.computeConnectionDensities(perfRates, b_perf, rsmax_perf, rvmax_perf, rvwmax_perf, surf_dens_perf, deferred_logger); this->connections_.computeConnectionPressureDelta(); }