From d23270c98fb4ce615ad63af1ecab61aa3acfcd51 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Fri, 18 Nov 2016 12:43:53 +0100 Subject: [PATCH] Fix the 2p simulator Only tested for oil+water case The blockmatrix and vectors are hardcoded to be 3 and a trivial equation is used for the Gas phase. --- opm/autodiff/BlackoilModelEbos.hpp | 278 +++++++++--------- .../SimulatorFullyImplicitBlackoilEbos.hpp | 2 +- ...mulatorFullyImplicitBlackoilOutputEbos.hpp | 2 +- opm/autodiff/StandardWellsDense.hpp | 184 ++++++++---- .../WellStateFullyImplicitBlackoilDense.hpp | 30 +- 5 files changed, 293 insertions(+), 203 deletions(-) diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index ae1d5981a..911a96522 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -622,99 +622,102 @@ namespace Opm { double& so = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Oil ]]; so -= step * dso; - // const double drmaxrel = drMaxRel(); - // Update rs and rv - if (has_disgas_) { - double& rs = reservoir_state.gasoilratio()[cell_idx]; - rs -= drs; - rs = std::max(rs, 0.0); - - } - if (has_vapoil_) { - double& rv = reservoir_state.rv()[cell_idx]; - rv -= drv; - rv = std::max(rv, 0.0); - } - - // Sg is used as primal variable for water only cells. - const double epsilon = 1e-4; //std::sqrt(std::numeric_limits::epsilon()); - double& sw = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Water ]]; - double& sg = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Gas ]]; - double& rs = reservoir_state.gasoilratio()[cell_idx]; - double& rv = reservoir_state.rv()[cell_idx]; - - - // phase translation sg <-> rs - const HydroCarbonState hydroCarbonState = reservoir_state.hydroCarbonState()[cell_idx]; - const auto& intQuants = *(ebosSimulator_.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); - const auto& fs = intQuants.fluidState(); - switch (hydroCarbonState) { - case HydroCarbonState::GasAndOil: { - - if (sw > (1.0 - epsilon)) // water only i.e. do nothing - break; - - if (sg <= 0.0 && has_disgas_) { - reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::OilOnly; // sg --> rs - sg = 0; - so = 1.0 - sw - sg; - const double& rsSat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]); + // phase for when oil and gas + if (active_[Gas] && active_[Oil] ) { + // const double drmaxrel = drMaxRel(); + // Update rs and rv + if (has_disgas_) { double& rs = reservoir_state.gasoilratio()[cell_idx]; - rs = rsSat*(1-epsilon); - } else if (so <= 0.0 && has_vapoil_) { - reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasOnly; // sg --> rv - so = 0; - sg = 1.0 - sw - so; + rs -= drs; + rs = std::max(rs, 0.0); + + } + if (has_vapoil_) { double& rv = reservoir_state.rv()[cell_idx]; - // use gas pressure? + rv -= drv; + rv = std::max(rv, 0.0); + } + + // Sg is used as primal variable for water only cells. + const double epsilon = 1e-4; //std::sqrt(std::numeric_limits::epsilon()); + double& sw = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Water ]]; + double& sg = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Gas ]]; + double& rs = reservoir_state.gasoilratio()[cell_idx]; + double& rv = reservoir_state.rv()[cell_idx]; + + + // phase translation sg <-> rs + const HydroCarbonState hydroCarbonState = reservoir_state.hydroCarbonState()[cell_idx]; + const auto& intQuants = *(ebosSimulator_.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); + const auto& fs = intQuants.fluidState(); + switch (hydroCarbonState) { + case HydroCarbonState::GasAndOil: { + + if (sw > (1.0 - epsilon)) // water only i.e. do nothing + break; + + if (sg <= 0.0 && has_disgas_) { + reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::OilOnly; // sg --> rs + sg = 0; + so = 1.0 - sw - sg; + const double& rsSat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]); + double& rs = reservoir_state.gasoilratio()[cell_idx]; + rs = rsSat*(1-epsilon); + } else if (so <= 0.0 && has_vapoil_) { + reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasOnly; // sg --> rv + so = 0; + sg = 1.0 - sw - so; + double& rv = reservoir_state.rv()[cell_idx]; + // use gas pressure? + const double& rvSat = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]); + rv = rvSat*(1-epsilon); + } + break; + } + case HydroCarbonState::OilOnly: { + if (sw > (1.0 - epsilon)) { + // water only change to Sg + rs = 0; + rv = 0; + reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasAndOil; + //std::cout << "watonly rv -> sg" << cell_idx << std::endl; + break; + } + + + + + const double& rsSat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]); + if (rs > ( rsSat * (1+epsilon) ) ) { + reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasAndOil; + sg = epsilon; + so -= epsilon; + rs = rsSat; + } + break; + } + case HydroCarbonState::GasOnly: { + if (sw > (1.0 - epsilon)) { + // water only change to Sg + rs = 0; + rv = 0; + reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasAndOil; + //std::cout << "watonly rv -> sg" << cell_idx << std::endl; + break; + } const double& rvSat = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]); - rv = rvSat*(1-epsilon); - } - break; - } - case HydroCarbonState::OilOnly: { - if (sw > (1.0 - epsilon)) { - // water only change to Sg - rs = 0; - rv = 0; - reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasAndOil; - //std::cout << "watonly rv -> sg" << cell_idx << std::endl; + if (rv > rvSat * (1+epsilon) ) { + reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasAndOil; + so = epsilon; + rv = rvSat; + sg -= epsilon; + } break; } - - - - const double& rsSat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]); - if (rs > ( rsSat * (1+epsilon) ) ) { - reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasAndOil; - sg = epsilon; - so -= epsilon; - rs = rsSat; + default: + OPM_THROW(std::logic_error, "Unknown primary variable enum value in cell " << cell_idx << ": " << hydroCarbonState); } - break; - } - case HydroCarbonState::GasOnly: { - if (sw > (1.0 - epsilon)) { - // water only change to Sg - rs = 0; - rv = 0; - reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasAndOil; - //std::cout << "watonly rv -> sg" << cell_idx << std::endl; - break; - } - const double& rvSat = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]); - if (rv > rvSat * (1+epsilon) ) { - reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasAndOil; - so = epsilon; - rv = rvSat; - sg -= epsilon; - } - break; - } - - default: - OPM_THROW(std::logic_error, "Unknown primary variable enum value in cell " << cell_idx << ": " << hydroCarbonState); } } @@ -1237,56 +1240,61 @@ namespace Opm { cellPv[BlackoilIndices::waterSaturationIdx] = saturations[cellIdx*numPhases + pu.phase_pos[Water]]; // set switching variable and interpretation - if( reservoirState.hydroCarbonState()[cellIdx] == HydroCarbonState::OilOnly && has_disgas_ ) - { - cellPv[BlackoilIndices::compositionSwitchIdx] = rs[cellIdx]; + if (active_[Gas] ) { + if( reservoirState.hydroCarbonState()[cellIdx] == HydroCarbonState::OilOnly && has_disgas_ ) + { + cellPv[BlackoilIndices::compositionSwitchIdx] = rs[cellIdx]; + cellPv[BlackoilIndices::pressureSwitchIdx] = oilPressure[cellIdx]; + cellPv.setPrimaryVarsMeaning( PrimaryVariables::Sw_po_Rs ); + } + else if( reservoirState.hydroCarbonState()[cellIdx] == HydroCarbonState::GasOnly && has_vapoil_ ) + { + // this case (-> gas only with vaporized oil in the gas) is + // relatively expensive as it requires to compute the capillary + // pressure in order to get the gas phase pressure. (the reason why + // ebos uses the gas pressure here is that it makes the common case + // of the primary variable switching code fast because to determine + // whether the oil phase appears one needs to compute the Rv value + // for the saturated gas phase and if this is not available as a + // primary variable, it needs to be computed.) luckily for here, the + // gas-only case is not too common, so the performance impact of this + // is limited. + typedef Opm::SimpleModularFluidState SatOnlyFluidState; + SatOnlyFluidState fluidState; + fluidState.setSaturation(FluidSystem::waterPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Water]]); + fluidState.setSaturation(FluidSystem::oilPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Oil]]); + fluidState.setSaturation(FluidSystem::gasPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Gas]]); + + double pC[/*numPhases=*/3] = { 0.0, 0.0, 0.0 }; + const MaterialLawParams& matParams = simulator.problem().materialLawParams(cellIdx); + MaterialLaw::capillaryPressures(pC, matParams, fluidState); + double pg = oilPressure[cellIdx] + (pC[FluidSystem::gasPhaseIdx] - pC[FluidSystem::oilPhaseIdx]); + + cellPv[BlackoilIndices::compositionSwitchIdx] = rv[cellIdx]; + cellPv[BlackoilIndices::pressureSwitchIdx] = pg; + cellPv.setPrimaryVarsMeaning( PrimaryVariables::Sw_pg_Rv ); + } + else + { + assert( reservoirState.hydroCarbonState()[cellIdx] == HydroCarbonState::GasAndOil); + cellPv[BlackoilIndices::compositionSwitchIdx] = saturations[cellIdx*numPhases + pu.phase_pos[Gas]]; + cellPv[BlackoilIndices::pressureSwitchIdx] = oilPressure[ cellIdx ]; + cellPv.setPrimaryVarsMeaning( PrimaryVariables::Sw_po_Sg ); + } + } else { + // for oil-water case oil pressure should be used as primary variable cellPv[BlackoilIndices::pressureSwitchIdx] = oilPressure[cellIdx]; - cellPv.setPrimaryVarsMeaning( PrimaryVariables::Sw_po_Rs ); - } - else if( reservoirState.hydroCarbonState()[cellIdx] == HydroCarbonState::GasOnly && has_vapoil_ ) - { - // this case (-> gas only with vaporized oil in the gas) is - // relatively expensive as it requires to compute the capillary - // pressure in order to get the gas phase pressure. (the reason why - // ebos uses the gas pressure here is that it makes the common case - // of the primary variable switching code fast because to determine - // whether the oil phase appears one needs to compute the Rv value - // for the saturated gas phase and if this is not available as a - // primary variable, it needs to be computed.) luckily for here, the - // gas-only case is not too common, so the performance impact of this - // is limited. - typedef Opm::SimpleModularFluidState SatOnlyFluidState; - SatOnlyFluidState fluidState; - fluidState.setSaturation(FluidSystem::waterPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Water]]); - fluidState.setSaturation(FluidSystem::oilPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Oil]]); - fluidState.setSaturation(FluidSystem::gasPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Gas]]); - - double pC[/*numPhases=*/3] = { 0.0, 0.0, 0.0 }; - const MaterialLawParams& matParams = simulator.problem().materialLawParams(cellIdx); - MaterialLaw::capillaryPressures(pC, matParams, fluidState); - double pg = oilPressure[cellIdx] + (pC[FluidSystem::gasPhaseIdx] - pC[FluidSystem::oilPhaseIdx]); - - cellPv[BlackoilIndices::compositionSwitchIdx] = rv[cellIdx]; - cellPv[BlackoilIndices::pressureSwitchIdx] = pg; - cellPv.setPrimaryVarsMeaning( PrimaryVariables::Sw_pg_Rv ); - } - else - { - assert( reservoirState.hydroCarbonState()[cellIdx] == HydroCarbonState::GasAndOil); - cellPv[BlackoilIndices::compositionSwitchIdx] = saturations[cellIdx*numPhases + pu.phase_pos[Gas]]; - cellPv[BlackoilIndices::pressureSwitchIdx] = oilPressure[ cellIdx ]; - cellPv.setPrimaryVarsMeaning( PrimaryVariables::Sw_po_Sg ); } } diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp index 353b7196c..106ed4bca 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp @@ -233,7 +233,7 @@ public: const Wells* wells = wells_manager.c_wells(); WellState well_state; - well_state.init(wells, state, prev_well_state); + well_state.init(wells, state, prev_well_state, props_.phaseUsage()); // give the polymer and surfactant simulators the chance to do their stuff handleAdditionalWellInflow(timer, wells_manager, well_state, wells); diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilOutputEbos.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilOutputEbos.hpp index edc912949..8b7f0a051 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilOutputEbos.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilOutputEbos.hpp @@ -254,7 +254,7 @@ namespace Opm std::unordered_set()); const Wells* wells = wellsmanager.c_wells(); - wellstate.resize(wells, simulatorstate); //Resize for restart step + wellstate.resize(wells, simulatorstate, phaseusage ); //Resize for restart step auto restarted = Opm::init_from_restart_file( eclipseState_, Opm::UgGridHelpers::numCells(grid) ); diff --git a/opm/autodiff/StandardWellsDense.hpp b/opm/autodiff/StandardWellsDense.hpp index ea5161d90..765f1f493 100644 --- a/opm/autodiff/StandardWellsDense.hpp +++ b/opm/autodiff/StandardWellsDense.hpp @@ -55,6 +55,12 @@ namespace Opm { +enum WellVariablePositions { + XvarWell = 0, + WFrac = 1, + GFrac = 2 +}; + /// Class for handling the standard well model. template @@ -62,15 +68,17 @@ namespace Opm { public: // --------- Types --------- - typedef DenseAd::Evaluation EvalWell; typedef WellStateFullyImplicitBlackoilDense WellState; typedef BlackoilModelParameters ModelParameters; typedef double Scalar; - typedef Dune::FieldVector VectorBlockType; - typedef Dune::FieldMatrix MatrixBlockType; + static const int blocksize = 3; + typedef Dune::FieldVector VectorBlockType; + typedef Dune::FieldMatrix MatrixBlockType; typedef Dune::BCRSMatrix Mat; typedef Dune::BlockVector BVector; + typedef DenseAd::Evaluation EvalWell; + // --------- Public methods --------- StandardWellsDense(const Wells* wells_arg, @@ -106,6 +114,10 @@ namespace Opm { return; } + // assumes the gas fractions are stored after water fractions + // WellVariablePositions needs to be changed for 2p runs + assert (np == 3 || np == 2 && !pu.phase_used[Gas] ); + fluid_ = fluid_arg; active_ = active_arg; vfp_properties_ = vfp_properties_arg; @@ -235,10 +247,16 @@ namespace Opm { for (int p2 = 0; p2 < np; ++p2) { if (!only_wells) { ebosJac[cell_idx][cell_idx][flowPhaseToEbosCompIdx(p1)][flowToEbosPvIdx(p2)] -= cq_s[p1].derivative(p2); - duneB_[w][cell_idx][flowToEbosPvIdx(p2)][flowPhaseToEbosCompIdx(p1)] -= cq_s[p1].derivative(p2+3); // intput in transformed matrix + duneB_[w][cell_idx][flowToEbosPvIdx(p2)][flowPhaseToEbosCompIdx(p1)] -= cq_s[p1].derivative(p2+blocksize); // intput in transformed matrix duneC_[w][cell_idx][flowPhaseToEbosCompIdx(p1)][flowToEbosPvIdx(p2)] -= cq_s[p1].derivative(p2); } - invDuneD_[w][w][flowPhaseToEbosCompIdx(p1)][flowToEbosPvIdx(p2)] -= cq_s[p1].derivative(p2+3); + invDuneD_[w][w][flowPhaseToEbosCompIdx(p1)][flowToEbosPvIdx(p2)] -= cq_s[p1].derivative(p2+blocksize); + } + + // add trivial equation for 2p cases (Only support water + oil) + if (np == 2) { + assert((*active_)[ Gas ]); + invDuneD_[w][w][flowPhaseToEbosCompIdx(Gas)][flowToEbosPvIdx(Gas)] = 1.0; } // Store the perforation phase flux for later usage. @@ -253,7 +271,7 @@ namespace Opm { EvalWell resWell_loc = (wellVolumeFraction(w, p1) - F0_[w + nw*p1]) * volume / dt; resWell_loc += getQs(w, p1); for (int p2 = 0; p2 < np; ++p2) { - invDuneD_[w][w][flowPhaseToEbosCompIdx(p1)][flowToEbosPvIdx(p2)] += resWell_loc.derivative(p2+3); + invDuneD_[w][w][flowPhaseToEbosCompIdx(p1)][flowToEbosPvIdx(p2)] += resWell_loc.derivative(p2+blocksize); } resWell_[w][flowPhaseToEbosCompIdx(p1)] += resWell_loc.value(); } @@ -464,11 +482,11 @@ namespace Opm { } - typedef DenseAd::Evaluation Eval; + typedef DenseAd::Evaluation Eval; EvalWell extendEval(Eval in) const { EvalWell out = 0.0; out.setValue(in.value()); - for(int i = 0;i<3;++i) { + for(int i = 0; i < blocksize;++i) { out.setDerivative(i, in.derivative(flowToEbosPvIdx(i))); } return out; @@ -482,7 +500,7 @@ namespace Opm { for (int w = 0; w < nw; ++w) { wellVariables_[w + nw*phaseIdx] = 0.0; wellVariables_[w + nw*phaseIdx].setValue(xw.wellSolutions()[w + nw* phaseIdx]); - wellVariables_[w + nw*phaseIdx].setDerivative(np + phaseIdx, 1.0); + wellVariables_[w + nw*phaseIdx].setDerivative(blocksize + phaseIdx, 1.0); } } } @@ -819,10 +837,14 @@ namespace Opm { const PhaseUsage& pu = fluid_->phaseUsage(); const int np = fluid_->numPhases(); b_perf.resize(nperf*np); - rsmax_perf.resize(nperf); - rvmax_perf.resize(nperf); surf_dens_perf.resize(nperf*np); + //rs and rv are only used if both oil and gas is present + if (pu.phase_used[BlackoilPhases::Vapour] && pu.phase_pos[BlackoilPhases::Liquid]) { + rsmax_perf.resize(nperf); + rvmax_perf.resize(nperf); + } + // Compute the average pressure in each well block for (int w = 0; w < nw; ++w) { for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) { @@ -915,34 +937,64 @@ namespace Opm { // update the second and third well variable (The flux fractions) std::vector F(np,0.0); - const int sign2 = dwells[w][flowPhaseToEbosCompIdx(1)] > 0 ? 1: -1; - const double dx2_limited = sign2 * std::min(std::abs(dwells[w][flowPhaseToEbosCompIdx(1)]),dFLimit); - well_state.wellSolutions()[nw + w] = xvar_well_old[nw + w] - dx2_limited; - const int sign3 = dwells[w][flowPhaseToEbosCompIdx(2)] > 0 ? 1: -1; - const double dx3_limited = sign3 * std::min(std::abs(dwells[w][flowPhaseToEbosCompIdx(2)]),dFLimit); - well_state.wellSolutions()[2*nw + w] = xvar_well_old[2*nw + w] - dx3_limited; - F[Water] = well_state.wellSolutions()[nw + w]; - F[Gas] = well_state.wellSolutions()[2*nw + w]; - F[Oil] = 1.0 - F[Water] - F[Gas]; - - if (F[Water] < 0.0) { - F[Gas] /= (1.0 - F[Water]); - F[Oil] /= (1.0 - F[Water]); - F[Water] = 0.0; + if ((*active_)[ Water ]) { + const int sign2 = dwells[w][flowPhaseToEbosCompIdx(WFrac)] > 0 ? 1: -1; + const double dx2_limited = sign2 * std::min(std::abs(dwells[w][flowPhaseToEbosCompIdx(WFrac)]),dFLimit); + well_state.wellSolutions()[WFrac*nw + w] = xvar_well_old[WFrac*nw + w] - dx2_limited; } - if (F[Gas] < 0.0) { - F[Water] /= (1.0 - F[Gas]); - F[Oil] /= (1.0 - F[Gas]); - F[Gas] = 0.0; + + if ((*active_)[ Gas ]) { + const int sign3 = dwells[w][flowPhaseToEbosCompIdx(GFrac)] > 0 ? 1: -1; + const double dx3_limited = sign3 * std::min(std::abs(dwells[w][flowPhaseToEbosCompIdx(GFrac)]),dFLimit); + well_state.wellSolutions()[GFrac*nw + w] = xvar_well_old[GFrac*nw + w] - dx3_limited; + } + + assert((*active_)[ Oil ]); + F[Oil] = 1.0; + if ((*active_)[ Water ]) { + F[Water] = well_state.wellSolutions()[WFrac*nw + w]; + F[Oil] -= F[Water]; + } + + if ((*active_)[ Gas ]) { + F[Gas] = well_state.wellSolutions()[GFrac*nw + w]; + F[Oil] -= F[Gas]; + } + + if ((*active_)[ Water ]) { + if (F[Water] < 0.0) { + if ((*active_)[ Gas ]) { + F[Gas] /= (1.0 - F[Water]); + } + F[Oil] /= (1.0 - F[Water]); + F[Water] = 0.0; + } + } + if ((*active_)[ Gas ]) { + if (F[Gas] < 0.0) { + if ((*active_)[ Water ]) { + F[Water] /= (1.0 - F[Gas]); + } + F[Oil] /= (1.0 - F[Gas]); + F[Gas] = 0.0; + } } if (F[Oil] < 0.0) { - F[Water] /= (1.0 - F[Oil]); - F[Gas] /= (1.0 - F[Oil]); + if ((*active_)[ Water ]) { + F[Water] /= (1.0 - F[Oil]); + } + if ((*active_)[ Gas ]) { + F[Gas] /= (1.0 - F[Oil]); + } F[Oil] = 0.0; } - well_state.wellSolutions()[nw + w] = F[Water]; - well_state.wellSolutions()[2*nw + w] = F[Gas]; - //std::cout << wells().name[w] << " "<< F[Water] << " " << F[Oil] << " " << F[Gas] << std::endl; + + if ((*active_)[ Water ]) { + well_state.wellSolutions()[WFrac*nw + w] = F[Water]; + } + if ((*active_)[ Gas ]) { + well_state.wellSolutions()[GFrac*nw + w] = F[Gas]; + } // The interpretation of the first well variable depends on the well control const WellControls* wc = wells().ctrls[w]; @@ -969,18 +1021,18 @@ namespace Opm { case THP: // The BHP and THP both uses the total rate as first well variable. case BHP: { - well_state.wellSolutions()[w] = xvar_well_old[w] - dwells[w][flowPhaseToEbosCompIdx(0)]; + well_state.wellSolutions()[nw*XvarWell + w] = xvar_well_old[nw*XvarWell + w] - dwells[w][flowPhaseToEbosCompIdx(XvarWell)]; switch (wells().type[w]) { case INJECTOR: for (int p = 0; p < np; ++p) { const double comp_frac = wells().comp_frac[np*w + p]; - well_state.wellRates()[w*np + p] = comp_frac * well_state.wellSolutions()[w]; + well_state.wellRates()[w*np + p] = comp_frac * well_state.wellSolutions()[nw*XvarWell + w]; } break; case PRODUCER: for (int p = 0; p < np; ++p) { - well_state.wellRates()[w*np + p] = well_state.wellSolutions()[w] * F[p]; + well_state.wellRates()[w*np + p] = well_state.wellSolutions()[nw*XvarWell + w] * F[p]; } break; } @@ -1038,10 +1090,10 @@ namespace Opm { case SURFACE_RATE: // Both rate controls use bhp as first well variable case RESERVOIR_RATE: { - const int sign1 = dwells[w][flowPhaseToEbosCompIdx(0)] > 0 ? 1: -1; - const double dx1_limited = sign1 * std::min(std::abs(dwells[w][flowPhaseToEbosCompIdx(0)]),std::abs(xvar_well_old[w])*dBHPLimit); - well_state.wellSolutions()[w] = std::max(xvar_well_old[w] - dx1_limited,1e5); - well_state.bhp()[w] = well_state.wellSolutions()[w]; + const int sign1 = dwells[w][flowPhaseToEbosCompIdx(XvarWell)] > 0 ? 1: -1; + const double dx1_limited = sign1 * std::min(std::abs(dwells[w][flowPhaseToEbosCompIdx(XvarWell)]),std::abs(xvar_well_old[nw*XvarWell + w])*dBHPLimit); + well_state.wellSolutions()[nw*XvarWell + w] = std::max(xvar_well_old[nw*XvarWell + w] - dx1_limited,1e5); + well_state.bhp()[w] = well_state.wellSolutions()[nw*XvarWell + w]; if (well_controls_iget_type(wc, current) == SURFACE_RATE) { if (wells().type[w]==PRODUCER) { @@ -1235,14 +1287,14 @@ namespace Opm { case BHP: { const WellType& well_type = wells().type[w]; - xw.wellSolutions()[w] = 0.0; + xw.wellSolutions()[nw*XvarWell + w] = 0.0; if (well_type == INJECTOR) { for (int p = 0; p < np; ++p) { - xw.wellSolutions()[w] += xw.wellRates()[np*w + p] * wells().comp_frac[np*w + p]; + xw.wellSolutions()[nw*XvarWell + w] += xw.wellRates()[np*w + p] * wells().comp_frac[np*w + p]; } } else { for (int p = 0; p < np; ++p) { - xw.wellSolutions()[w] += g[p] * xw.wellRates()[np*w + p]; + xw.wellSolutions()[nw*XvarWell + w] += g[p] * xw.wellRates()[np*w + p]; } } } @@ -1252,7 +1304,7 @@ namespace Opm { case RESERVOIR_RATE: // Intentional fall-through case SURFACE_RATE: { - xw.wellSolutions()[w] = xw.bhp()[w]; + xw.wellSolutions()[nw*XvarWell + w] = xw.bhp()[w]; } break; } @@ -1262,11 +1314,20 @@ namespace Opm { tot_well_rate += g[p] * xw.wellRates()[np*w + p]; } if(std::abs(tot_well_rate) > 0) { - xw.wellSolutions()[nw + w] = g[Water] * xw.wellRates()[np*w + Water] / tot_well_rate; - xw.wellSolutions()[2*nw + w] = g[Gas] * xw.wellRates()[np*w + Gas] / tot_well_rate ; + if ((*active_)[ Water ]) { + xw.wellSolutions()[WFrac*nw + w] = g[Water] * xw.wellRates()[np*w + Water] / tot_well_rate; + } + if ((*active_)[ Gas ]) { + xw.wellSolutions()[GFrac*nw + w] = g[Gas] * xw.wellRates()[np*w + Gas] / tot_well_rate ; + } } else { - xw.wellSolutions()[nw + w] = wells().comp_frac[np*w + Water]; - xw.wellSolutions()[2 * nw + w] = wells().comp_frac[np*w + Gas]; + if ((*active_)[ Water ]) { + xw.wellSolutions()[WFrac*nw + w] = wells().comp_frac[np*w + Water]; + } + + if ((*active_)[ Gas ]) { + xw.wellSolutions()[GFrac*nw + w] = wells().comp_frac[np*w + Gas]; + } } } } @@ -1482,13 +1543,15 @@ namespace Opm { return bhp; } - return wellVariables_[wellIdx]; + const int nw = wells().number_of_wells; + return wellVariables_[nw*XvarWell + wellIdx]; } EvalWell getQs(const int wellIdx, const int phaseIdx) const { EvalWell qs = 0.0; const WellControls* wc = wells().ctrls[wellIdx]; const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; const double target_rate = well_controls_get_current_target(wc); if (wells().type[wellIdx] == INJECTOR) { @@ -1497,7 +1560,7 @@ namespace Opm { return qs; if (well_controls_get_current_type(wc) == BHP || well_controls_get_current_type(wc) == THP) { - return wellVariables_[wellIdx]; + return wellVariables_[nw*XvarWell + wellIdx]; } qs.setValue(target_rate); return qs; @@ -1505,7 +1568,7 @@ namespace Opm { // Producers if (well_controls_get_current_type(wc) == BHP || well_controls_get_current_type(wc) == THP ) { - return wellVariables_[wellIdx] * wellVolumeFractionScaled(wellIdx,phaseIdx); + return wellVariables_[nw*XvarWell + wellIdx] * wellVolumeFractionScaled(wellIdx,phaseIdx); } if (well_controls_get_current_type(wc) == SURFACE_RATE) { const double comp_frac = wells().comp_frac[np*wellIdx + phaseIdx]; @@ -1530,18 +1593,25 @@ namespace Opm { } EvalWell wellVolumeFraction(const int wellIdx, const int phaseIdx) const { - assert(wells().number_of_phases == 3); const int nw = wells().number_of_wells; if (phaseIdx == Water) { - return wellVariables_[nw + wellIdx]; + return wellVariables_[WFrac * nw + wellIdx]; } if (phaseIdx == Gas) { - return wellVariables_[2*nw + wellIdx]; + return wellVariables_[GFrac * nw + wellIdx]; } - // Oil - return 1.0 - wellVariables_[nw + wellIdx] - wellVariables_[2 * nw + wellIdx]; + // Oil fraction + EvalWell well_fraction = 1.0; + if ((*active_)[Water]) { + well_fraction -= wellVariables_[WFrac * nw + wellIdx]; + } + + if ((*active_)[Gas]) { + well_fraction -= wellVariables_[GFrac * nw + wellIdx]; + } + return well_fraction; } EvalWell wellVolumeFractionScaled(const int wellIdx, const int phaseIdx) const { diff --git a/opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp b/opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp index 565e899c6..0718ac931 100644 --- a/opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp +++ b/opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp @@ -62,7 +62,7 @@ namespace Opm /// to give useful initial values to the bhp(), wellRates() /// and perfPhaseRates() fields, depending on controls template - void init(const Wells* wells, const State& state, const PrevState& prevState) + void init(const Wells* wells, const State& state, const PrevState& prevState, const PhaseUsage& pu) { // call init on base class BaseType :: init(wells, state, prevState); @@ -115,19 +115,31 @@ namespace Opm } break; } -#warning "HACK: that assert() was probably there for a reason!" - //assert(np == 3); + double total_rates = 0.0; for (int p = 0; p < np; ++p) { total_rates += g[p] * wellRates()[np*w + p]; } + const int waterpos = pu.phase_pos[Water]; + const int gaspos = pu.phase_pos[Gas]; + + assert (np == 3 || np == 2 && !pu.phase_used[Gas] ); + // assumes the gas fractions are stored after water fractions if(std::abs(total_rates) > 0) { - wellSolutions()[nw + w] = g[Water] * wellRates()[np*w + Water] / total_rates; //wells->comp_frac[np*w + Water]; // Water; - wellSolutions()[2*nw + w] = g[Gas] * wellRates()[np*w + Gas] / total_rates ; //wells->comp_frac[np*w + Gas]; //Gas + if( pu.phase_used[Water] ) { + wellSolutions()[nw + w] = g[Water] * wellRates()[np*w + waterpos] / total_rates; + } + if( pu.phase_used[Gas] ) { + wellSolutions()[2*nw + w] = g[Gas] * wellRates()[np*w + gaspos] / total_rates ; + } } else { - wellSolutions()[nw + w] = wells->comp_frac[np*w + Water]; - wellSolutions()[2*nw + w] = wells->comp_frac[np*w + Gas]; + if( pu.phase_used[Water] ) { + wellSolutions()[nw + w] = wells->comp_frac[np*w + waterpos]; + } + if( pu.phase_used[Gas] ) { + wellSolutions()[2*nw + w] = wells->comp_frac[np*w + gaspos]; + } } } @@ -159,9 +171,9 @@ namespace Opm } template - void resize(const Wells* wells, const State& state) { + void resize(const Wells* wells, const State& state, const PhaseUsage& pu ) { const WellStateFullyImplicitBlackoilDense dummy_state{}; // Init with an empty previous state only resizes - init(wells, state, dummy_state) ; + init(wells, state, dummy_state, pu) ; }