From b987e4b3245ccfbf4894575c6e24ea450e9bc1fa Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Tue, 9 May 2017 08:21:51 +0200 Subject: [PATCH 1/6] Implement solvent model in flow_ebos 1) Extends the well model to account for solvent surface volumes 2) Add solvent to updateState 3) Add solvent to well and field output The solvent parts is encapsled in if (has_solvent_) and should not effect the standard runs. --- opm/autodiff/BlackoilModelEbos.hpp | 98 +++++++++-- .../SimulatorFullyImplicitBlackoilEbos.hpp | 4 +- .../SimulatorFullyImplicitBlackoilOutput.hpp | 2 +- opm/autodiff/StandardWellsDense.hpp | 21 ++- opm/autodiff/StandardWellsDense_impl.hpp | 161 ++++++++++++++++-- .../WellStateFullyImplicitBlackoilDense.hpp | 47 ++++- 6 files changed, 298 insertions(+), 35 deletions(-) diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 81a1dc5ab..3905786f2 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -84,6 +84,7 @@ NEW_TYPE_TAG(EclFlowProblem, INHERITS_FROM(BlackOilModel, EclBaseProblem)); SET_BOOL_PROP(EclFlowProblem, DisableWells, true); SET_BOOL_PROP(EclFlowProblem, EnableDebuggingChecks, false); SET_BOOL_PROP(EclFlowProblem, ExportGlobalTransmissibility, true); +SET_BOOL_PROP(EclFlowProblem, EnableSolvent, false); // SWATINIT is done by the flow part of flow_ebos. this can be removed once the legacy // code for fluid and satfunc handling gets fully retired. @@ -119,6 +120,7 @@ namespace Opm { typedef double Scalar; static const int numEq = BlackoilIndices::numEq; + static const int solventCompIdx = 3; //TODO get this from ebos typedef Dune::FieldVector VectorBlockType; typedef Dune::FieldMatrix MatrixBlockType; typedef Dune::BCRSMatrix Mat; @@ -162,6 +164,7 @@ namespace Opm { , active_(detail::activePhases(fluid.phaseUsage())) , has_disgas_(FluidSystem::enableDissolvedGas()) , has_vapoil_(FluidSystem::enableVaporizedOil()) + , has_solvent_(GET_PROP_VALUE(TypeTag, EnableSolvent)) , param_( param ) , well_model_ (well_model) , terminal_output_ (terminal_output) @@ -624,6 +627,10 @@ namespace Opm { } dso -= dsg; + // solvent + const double dss = has_solvent_ ? dx[cell_idx][BlackoilIndices::solventSaturationIdx] : 0.0; + dso -= dss; + // Appleyard chop process. maxVal = std::max(std::abs(dsg),maxVal); double step = dsMax()/maxVal; @@ -638,6 +645,12 @@ namespace Opm { double& sg = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Gas ]]; sg -= step * dsg; } + + if (has_solvent_) { + double& ss = reservoir_state.getCellData( reservoir_state.SSOL )[cell_idx]; + ss -= step * dss; + } + double& so = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Oil ]]; so -= step * dso; @@ -682,12 +695,20 @@ namespace Opm { if (sg <= 0.0 && has_disgas_) { reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::OilOnly; // sg --> rs sg = 0; - so = 1.0 - sw - sg; + so = 1.0 - sw; + if (has_solvent_) { + double& ss = reservoir_state.getCellData( reservoir_state.SSOL )[cell_idx]; + so -= ss; + } rs *= (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; + sg = 1.0 - sw; + if (has_solvent_) { + double& ss = reservoir_state.getCellData( reservoir_state.SSOL )[cell_idx]; + sg -= ss; + } rv *= (1-epsilon); } break; @@ -862,6 +883,7 @@ namespace Opm { const double tol_wells = param_.tolerance_wells_; const int nc = Opm::AutoDiffGrid::numCells(grid_); + const int np = numPhases(); const int numComp = numComponents(); Vector R_sum(numComp); @@ -895,16 +917,25 @@ namespace Opm { const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0); const auto& fs = intQuants.fluidState(); - for ( int compIdx = 0; compIdx < numComp; ++compIdx ) + for ( int phaseIdx = 0; phaseIdx < np; ++phaseIdx ) { - Vector& R2_idx = R2[ compIdx ]; - Vector& B_idx = B[ compIdx ]; - const int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(compIdx); - const int ebosCompIdx = flowPhaseToEbosCompIdx(compIdx); + Vector& R2_idx = R2[ phaseIdx ]; + Vector& B_idx = B[ phaseIdx ]; + const int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phaseIdx); + const int ebosCompIdx = flowPhaseToEbosCompIdx(phaseIdx); B_idx [cell_idx] = 1.0 / fs.invB(ebosPhaseIdx).value(); R2_idx[cell_idx] = ebosResid[cell_idx][ebosCompIdx]; } + + if (has_solvent_ ) { + Vector& R2_idx = R2[ solventCompIdx ]; + Vector& B_idx = B[ solventCompIdx ]; + B_idx [cell_idx] = 1.0 / intQuants.solventInverseFormationVolumeFactor().value(); + const int ebosCompIdx = flowPhaseToEbosCompIdx(solventCompIdx); + R2_idx[cell_idx] = ebosResid[cell_idx][ebosCompIdx]; + } + } Vector pv_vector; @@ -969,6 +1000,9 @@ namespace Opm { const std::string& phaseName = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(phaseIdx)); key[ phaseIdx ] = std::toupper( phaseName.front() ); } + if (has_solvent_) { + key[ solventCompIdx ] = "S"; + } for (int compIdx = 0; compIdx < numComp; ++compIdx) { msg += " MB(" + key[ compIdx ] + ") "; @@ -1029,7 +1063,11 @@ namespace Opm { if (numPhases() == 2) { return 2; } - return FluidSystem::numComponents; + int numComp = FluidSystem::numComponents; + if (has_solvent_) + numComp ++; + + return numComp; } /// Wrapper required due to not following generic API @@ -1273,6 +1311,11 @@ namespace Opm { VectorType& pcSwMdc_ow = simData.getCellData( "PCSWMDC_OW" ); VectorType& krnSwMdc_ow = simData.getCellData( "KRNSWMDC_OW" ); + if (has_solvent_) { + simData.registerCellData( "SSOL", 1 ); + } + VectorType& ssol = has_solvent_ ? simData.getCellData( "SSOL" ) : zero; + std::vector failed_cells_pb; std::vector failed_cells_pd; const auto& gridView = ebosSimulator().gridView(); @@ -1356,6 +1399,11 @@ namespace Opm { krOil[cellIdx] = intQuants.relativePermeability(FluidSystem::oilPhaseIdx).value(); } + if (has_solvent_) + { + ssol[cellIdx] = intQuants.solventSaturation().value(); + } + // hack to make the intial output of rs and rv Ecl compatible. // For cells with swat == 1 Ecl outputs; rs = rsSat and rv=rvSat, in all but the initial step // where it outputs rs and rv values calculated by the initialization. To be compatible we overwrite @@ -1464,6 +1512,7 @@ namespace Opm { const std::vector cells_; // All grid cells const bool has_disgas_; const bool has_vapoil_; + const bool has_solvent_; ModelParameters param_; SimulatorReport failureReport_; @@ -1524,6 +1573,10 @@ namespace Opm { // set water saturation cellPv[BlackoilIndices::waterSaturationIdx] = saturations[cellIdx*numPhases + pu.phase_pos[Water]]; + if (has_solvent_) { + cellPv[BlackoilIndices::solventSaturationIdx] = reservoirState.getCellData( reservoirState.SSOL )[cellIdx]; + } + // set switching variable and interpretation if (active_[Gas] ) { if( reservoirState.hydroCarbonState()[cellIdx] == HydroCarbonState::OilOnly && has_disgas_ ) @@ -1599,19 +1652,18 @@ namespace Opm { int flowToEbosPvIdx( const int flowPv ) const { - assert(flowPv < 3); - const int flowToEbos[ 3 ] = { + const int flowToEbos[ 4 ] = { BlackoilIndices::pressureSwitchIdx, BlackoilIndices::waterSaturationIdx, - BlackoilIndices::compositionSwitchIdx + BlackoilIndices::compositionSwitchIdx, + BlackoilIndices::solventSaturationIdx }; return flowToEbos[ flowPv ]; } int flowPhaseToEbosCompIdx( const int phaseIdx ) const { - assert(phaseIdx < 3); - const int phaseToComp[ 3 ] = { FluidSystem::waterCompIdx, FluidSystem::oilCompIdx, FluidSystem::gasCompIdx }; + const int phaseToComp[ 4 ] = { FluidSystem::waterCompIdx, FluidSystem::oilCompIdx, FluidSystem::gasCompIdx, solventCompIdx }; return phaseToComp[ phaseIdx ]; } @@ -1644,6 +1696,13 @@ namespace Opm { cellRes[ flowPhaseToEbosCompIdx( flowPhaseIdx ) ] /= refDens; cellRes[ flowPhaseToEbosCompIdx( flowPhaseIdx ) ] *= cellVolume; } + if (has_solvent_) { + // no need to store refDens for all cells? + const auto& intQuants = ebosSimulator_.model().cachedIntensiveQuantities(cellIdx, /*timeIdx=*/0); + const auto& refDens = intQuants->solventRefDensity(); + cellRes[ solventCompIdx ] /= refDens; + cellRes[ solventCompIdx ] *= cellVolume; + } } for( auto row = ebosJac.begin(); row != endrow; ++row ) @@ -1669,13 +1728,24 @@ namespace Opm { (*col)[ebosCompIdx][flowToEbosPvIdx(pvIdx)] *= cellVolume; } } + if (has_solvent_) { + // TODO store refDens pr pvtRegion? + const auto& intQuants = ebosSimulator_.model().cachedIntensiveQuantities(rowIdx, /*timeIdx=*/0); + const auto& refDens = intQuants->solventRefDensity(); + for( int pvIdx=0; pvIdx < numEq; ++pvIdx ) + { + (*col)[solventCompIdx][flowToEbosPvIdx(pvIdx)] /= refDens; + (*col)[solventCompIdx][flowToEbosPvIdx(pvIdx)] *= cellVolume; + } + } } } } int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const { - const int flowToEbos[ 3 ] = { FluidSystem::waterPhaseIdx, FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx }; + assert(phaseIdx < 3); + const int flowToEbos[ 3 ] = { FluidSystem::waterPhaseIdx, FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx}; return flowToEbos[ phaseIdx ]; } diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp index 387a055b6..ffa5ad684 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp @@ -252,7 +252,9 @@ public: // Run a multiple steps of the solver depending on the time step control. solver_timer.start(); - WellModel well_model(wells, &(wells_manager.wellCollection()), model_param_, terminal_output_); + const auto& wells_ecl = eclState().getSchedule().getWells(timer.currentStepNum()); + const WellModel well_model(wells, &(wells_manager.wellCollection()), wells_ecl, model_param_, terminal_output_, timer.currentStepNum()); + auto solver = createSolver(well_model); std::vector> currentFluidInPlace; diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp index 871bc9e98..49fefbbc0 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp @@ -62,7 +62,7 @@ namespace Opm { class SimulationDataContainer; - class BlackoilState; + class BlackoilSolventState; void outputStateVtk(const UnstructuredGrid& grid, const Opm::SimulationDataContainer& state, diff --git a/opm/autodiff/StandardWellsDense.hpp b/opm/autodiff/StandardWellsDense.hpp index 5c2906a49..6c862d6d0 100644 --- a/opm/autodiff/StandardWellsDense.hpp +++ b/opm/autodiff/StandardWellsDense.hpp @@ -62,7 +62,8 @@ namespace Opm { enum WellVariablePositions { XvarWell = 0, WFrac = 1, - GFrac = 2 + GFrac = 2, + SFrac = 3 }; @@ -85,6 +86,7 @@ enum WellVariablePositions { typedef double Scalar; static const int numEq = BlackoilIndices::numEq; static const int numWellEq = numEq; //number of wellEq is the same as numEq in the model + static const int solventCompIdx = 3; //TODO get this from ebos typedef Dune::FieldVector VectorBlockType; typedef Dune::FieldMatrix MatrixBlockType; typedef Dune::BCRSMatrix Mat; @@ -102,8 +104,10 @@ enum WellVariablePositions { // --------- Public methods --------- StandardWellsDense(const Wells* wells_arg, WellCollection* well_collection, + const std::vector< const Well* >& wells_ecl, const ModelParameters& param, - const bool terminal_output); + const bool terminal_output, + const int current_index); void init(const PhaseUsage phase_usage_arg, const std::vector& active_arg, @@ -118,10 +122,14 @@ enum WellVariablePositions { /// The number of components in the model. int numComponents() const { - if (phase_usage_.num_phases == 2) { + if (numPhases() == 2) { return 2; } - return FluidSystem::numComponents; + int numComp = FluidSystem::numComponents; + if (has_solvent_) + numComp ++; + + return numComp; } @@ -279,12 +287,15 @@ enum WellVariablePositions { protected: bool wells_active_; const Wells* wells_; + const std::vector< const Well* > wells_ecl_; // Well collection is used to enforce the group control WellCollection* well_collection_; ModelParameters param_; bool terminal_output_; + bool has_solvent_; + int current_timeIdx_; PhaseUsage phase_usage_; std::vector active_; @@ -385,6 +396,8 @@ enum WellVariablePositions { const int well_index, const double initial_bhp, // bhp from BHP constraints const std::vector& initial_potential) const; + + double wsolvent(const int well_index) const; }; diff --git a/opm/autodiff/StandardWellsDense_impl.hpp b/opm/autodiff/StandardWellsDense_impl.hpp index d1af3d5dd..751f09832 100644 --- a/opm/autodiff/StandardWellsDense_impl.hpp +++ b/opm/autodiff/StandardWellsDense_impl.hpp @@ -7,13 +7,18 @@ namespace Opm { StandardWellsDense:: StandardWellsDense(const Wells* wells_arg, WellCollection* well_collection, + const std::vector< const Well* >& wells_ecl, const ModelParameters& param, - const bool terminal_output) + const bool terminal_output, + const int current_timeIdx) : wells_active_(wells_arg!=nullptr) , wells_(wells_arg) + , wells_ecl_(wells_ecl) , well_collection_(well_collection) , param_(param) , terminal_output_(terminal_output) + , has_solvent_(GET_PROP_VALUE(TypeTag, EnableSolvent)) + , current_timeIdx_(current_timeIdx) , well_perforation_efficiency_factors_((wells_!=nullptr ? wells_->well_connpos[wells_->number_of_wells] : 0), 1.0) , well_perforation_densities_( wells_ ? wells_arg->well_connpos[wells_arg->number_of_wells] : 0) , well_perforation_pressure_diffs_( wells_ ? wells_arg->well_connpos[wells_arg->number_of_wells] : 0) @@ -168,6 +173,7 @@ namespace Opm { { const int nw = wells().number_of_wells; const int numComp = numComponents(); + const int np = numPhases(); // clear all entries duneB_ = 0.0; @@ -223,7 +229,11 @@ namespace Opm { } // Store the perforation phase flux for later usage. - well_state.perfPhaseRates()[perf*numComp + componentIdx] = cq_s[componentIdx].value(); + if (componentIdx == solventCompIdx) {// if (flowPhaseToEbosCompIdx(componentIdx) == Solvent) + well_state.perfRateSolvent()[perf] = cq_s[componentIdx].value(); + } else { + well_state.perfPhaseRates()[perf*np + componentIdx] = cq_s[componentIdx].value(); + } } // Store the perforation pressure for later usage. @@ -267,6 +277,9 @@ namespace Opm { int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase); mob[phase] = extendEval(intQuants.mobility(ebosPhaseIdx)); } + if (has_solvent_) { + mob[solventCompIdx] = extendEval(intQuants.solventMobility()); + } } else { const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx); @@ -281,6 +294,11 @@ namespace Opm { int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase); mob[phase] = extendEval(relativePerms[ebosPhaseIdx] / intQuants.fluidState().viscosity(ebosPhaseIdx)); } + + // this may not work if viscosity and relperms has been modified? + if (has_solvent_) { + OPM_THROW(std::runtime_error, "individual mobility for wells does not work in combination with solvent"); + } } } @@ -443,8 +461,7 @@ namespace Opm { StandardWellsDense:: flowPhaseToEbosCompIdx( const int phaseIdx ) const { - assert(phaseIdx < 3); - const int phaseToComp[ 3 ] = { FluidSystem::waterCompIdx, FluidSystem::oilCompIdx, FluidSystem::gasCompIdx }; + const int phaseToComp[ 4 ] = { FluidSystem::waterCompIdx, FluidSystem::oilCompIdx, FluidSystem::gasCompIdx, solventCompIdx }; return phaseToComp[ phaseIdx ]; } @@ -457,11 +474,11 @@ namespace Opm { StandardWellsDense:: flowToEbosPvIdx( const int flowPv ) const { - assert(flowPv < 3); - const int flowToEbos[ 3 ] = { + const int flowToEbos[ 4 ] = { BlackoilIndices::pressureSwitchIdx, BlackoilIndices::waterSaturationIdx, - BlackoilIndices::compositionSwitchIdx + BlackoilIndices::compositionSwitchIdx, + BlackoilIndices::solventSaturationIdx }; return flowToEbos[ flowPv ]; } @@ -754,6 +771,9 @@ namespace Opm { int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase); b_perfcells_dense[phase] = extendEval(fs.invB(ebosPhaseIdx)); } + if (has_solvent_) { + b_perfcells_dense[solventCompIdx] = extendEval(intQuants.solventInverseFormationVolumeFactor()); + } // Pressure drawdown (also used to determine direction of flow) EvalWell well_pressure = bhp + cdp; @@ -803,6 +823,10 @@ namespace Opm { volumeRatio += cmix_s[watpos] / b_perfcells_dense[watpos]; } + if (has_solvent_) { + volumeRatio += cmix_s[solventCompIdx] / b_perfcells_dense[solventCompIdx]; + } + if (active_[Oil] && active_[Gas]) { const int oilpos = pu.phase_pos[Oil]; const int gaspos = pu.phase_pos[Gas]; @@ -989,6 +1013,10 @@ namespace Opm { B += 1 / fs.invB(ebosPhaseIdx).value(); } + if (has_solvent_) { + auto& B = B_avg[ solventCompIdx ]; + B += 1 / intQuants.solventInverseFormationVolumeFactor().value(); + } } // compute global average @@ -1181,6 +1209,12 @@ namespace Opm { for (int p = 0; p < pu.num_phases; ++p) { surf_dens_perf[numComp*perf + p] = FluidSystem::referenceDensity( flowPhaseToEbosPhaseIdx( p ), fs.pvtRegionIndex()); } + + #warning HACK use cell values for solvent injector + if (has_solvent_) { + b_perf[numComp*perf + solventCompIdx] = intQuants.solventInverseFormationVolumeFactor().value(); + surf_dens_perf[numComp*perf + solventCompIdx] = intQuants.solventRefDensity(); + } } } } @@ -1220,6 +1254,12 @@ namespace Opm { well_state.wellSolutions()[GFrac*nw + w] = xvar_well_old[GFrac*nw + w] - dx3_limited; } + if (has_solvent_) { + const int sign4 = dwells[w][flowPhaseToEbosCompIdx(SFrac)] > 0 ? 1: -1; + const double dx4_limited = sign4 * std::min(std::abs(dwells[w][flowPhaseToEbosCompIdx(SFrac)]),dFLimit); + well_state.wellSolutions()[SFrac*nw + w] = xvar_well_old[SFrac*nw + w] - dx4_limited; + } + assert(active_[ Oil ]); F[Oil] = 1.0; if (active_[ Water ]) { @@ -1232,11 +1272,20 @@ namespace Opm { F[Oil] -= F[Gas]; } + double F_solvent = 0.0; + if (has_solvent_) { + F_solvent = well_state.wellSolutions()[SFrac*nw + w]; + F[Oil] -= F_solvent; + } + if (active_[ Water ]) { if (F[Water] < 0.0) { if (active_[ Gas ]) { F[Gas] /= (1.0 - F[Water]); } + if (has_solvent_) { + F_solvent /= (1.0 - F[Water]); + } F[Oil] /= (1.0 - F[Water]); F[Water] = 0.0; } @@ -1246,6 +1295,9 @@ namespace Opm { if (active_[ Water ]) { F[Water] /= (1.0 - F[Gas]); } + if (has_solvent_) { + F_solvent /= (1.0 - F[Gas]); + } F[Oil] /= (1.0 - F[Gas]); F[Gas] = 0.0; } @@ -1257,6 +1309,9 @@ namespace Opm { if (active_[ Gas ]) { F[Gas] /= (1.0 - F[Oil]); } + if (has_solvent_) { + F_solvent /= (1.0 - F[Oil]); + } F[Oil] = 0.0; } @@ -1266,6 +1321,14 @@ namespace Opm { if (active_[ Gas ]) { well_state.wellSolutions()[GFrac*nw + w] = F[Gas]; } + if(has_solvent_) { + well_state.wellSolutions()[SFrac*nw + w] = F_solvent; + } + + # warning F_solvent is added to F_gas. This means that well_rate[Gas] also contains solvent. More testing is needed to make sure this is correct for output, wellControls, well groups, THP etc. + if (has_solvent_){ + F[Gas] += F_solvent; + } // The interpretation of the first well variable depends on the well control const WellControls* wc = wells().ctrls[w]; @@ -1699,9 +1762,21 @@ namespace Opm { const double grav) { // Compute densities + const int nperf = depth_perf.size(); + const int numComponent = b_perf.size() / nperf; + const int np = wells().number_of_phases; + std::vector perfRates(b_perf.size(),0.0); + for (int perf = 0; perf < nperf; ++perf) { + for (int phase = 0; phase < np; ++phase) { + perfRates[perf*numComponent + phase] = xw.perfPhaseRates()[perf*np + phase]; + } + if(has_solvent_) { + perfRates[perf*numComponent + solventCompIdx] = xw.perfRateSolvent()[perf]; + } + } well_perforation_densities_ = WellDensitySegmented::computeConnectionDensities( - wells(), phase_usage_, xw.perfPhaseRates(), + wells(), phase_usage_, perfRates, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); // Compute pressure deltas @@ -2065,13 +2140,27 @@ namespace Opm { EvalWell qs = 0.0; const WellControls* wc = wells().ctrls[wellIdx]; const int np = wells().number_of_phases; - assert(compIdx < np); + assert(compIdx < numComponents()); const int nw = wells().number_of_wells; + const auto pu = phase_usage_; const double target_rate = well_controls_get_current_target(wc); // TODO: the formulation for the injectors decides it only work with single phase // surface rate injection control. Improvement will be required. if (wells().type[wellIdx] == INJECTOR) { + if (has_solvent_ ) { + if (well_controls_get_current_type(wc) == BHP || well_controls_get_current_type(wc) == THP) { + OPM_THROW(std::runtime_error,"BHP controlled solvent injector is unsupported. Check well " + << wells().name [wellIdx] ); + } + if (compIdx == pu.phase_pos[ Gas ]) { //gas + qs.setValue(target_rate * (1.0 - wsolvent(wellIdx))); + return qs; + } else if (compIdx == solventCompIdx) { // solvent + qs.setValue(wsolvent(wellIdx) * target_rate); + return qs; + } + } const double comp_frac = wells().comp_frac[np*wellIdx + compIdx]; if (comp_frac == 0.0) { return qs; @@ -2175,6 +2264,10 @@ namespace Opm { return wellVariables_[GFrac * nw + wellIdx]; } + if (compIdx == solventCompIdx) { + return wellVariables_[SFrac * nw + wellIdx]; + } + // Oil fraction EvalWell well_fraction = 1.0; if (active_[Water]) { @@ -2184,6 +2277,9 @@ namespace Opm { if (active_[Gas]) { well_fraction -= wellVariables_[GFrac * nw + wellIdx]; } + if (has_solvent_) { + well_fraction -= wellVariables_[SFrac * nw + wellIdx]; + } return well_fraction; } @@ -2208,8 +2304,7 @@ namespace Opm { return wellVolumeFraction(wellIdx, compIdx); } } - assert(compIdx < 3); - std::vector g = {1,1,0.01}; + std::vector g = {1,1,0.01,0.01}; return (wellVolumeFraction(wellIdx, compIdx) / g[compIdx]); } @@ -2601,7 +2696,10 @@ namespace Opm { xw.wellSolutions()[WFrac*nw + well_index] = g[Water] * xw.wellRates()[np*well_index + Water] / tot_well_rate; } if (active_[ Gas ]) { - xw.wellSolutions()[GFrac*nw + well_index] = g[Gas] * xw.wellRates()[np*well_index + Gas] / tot_well_rate ; + xw.wellSolutions()[GFrac*nw + well_index] = g[Gas] * (1.0 - wsolvent(well_index)) * xw.wellRates()[np*well_index + Gas] / tot_well_rate ; + } + if (has_solvent_) { + xw.wellSolutions()[SFrac*nw + well_index] = g[Gas] * wsolvent(well_index) * xw.wellRates()[np*well_index + Gas] / tot_well_rate ; } } else { const WellType& well_type = wells().type[well_index]; @@ -2617,7 +2715,10 @@ namespace Opm { if (active_[Gas]) { if (distr[Gas] > 0.0) { - xw.wellSolutions()[GFrac * nw + well_index] = 1.0; + xw.wellSolutions()[GFrac * nw + well_index] = 1.0 - wsolvent(well_index); + if (has_solvent_) { + xw.wellSolutions()[SFrac * nw + well_index] = wsolvent(well_index); + } } else { xw.wellSolutions()[GFrac * nw + well_index] = 0.0; } @@ -2879,4 +2980,38 @@ namespace Opm { return potentials; } + template + double + StandardWellsDense:: + wsolvent(const int well_index) const { + + if (!has_solvent_) { + return 0.0; + } + + // loop over all wells until we find the well with the matching name + for (const auto& well : wells_ecl_) { + if (well->getStatus( current_timeIdx_ ) == WellCommon::SHUT) { + continue; + } + + WellInjectionProperties injection = well->getInjectionProperties(current_timeIdx_); + if (injection.injectorType == WellInjector::GAS) { + + double solventFraction = well->getSolventFraction(current_timeIdx_); + + // Look until we find the correct well + if (well->name() == wells().name[well_index]) { + return solventFraction; + } + } + } + // we didn't find it return 0; + // or should we throw? + return 0.0; + } + + + + } // namespace Opm diff --git a/opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp b/opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp index de4b945f2..fc1ef92c2 100644 --- a/opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp +++ b/opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp @@ -88,9 +88,14 @@ namespace Opm if (nw == 0) { return; } + const int nperf = wells_->well_connpos[nw]; + perfRateSolvent_.clear(); + perfRateSolvent_.resize(nperf, 0.0); + const int np = wells_->number_of_phases; + const int numComp = pu.has_solvent? np+1:np; well_solutions_.clear(); - well_solutions_.resize(nw * np, 0.0); + well_solutions_.resize(nw * numComp, 0.0); std::vector g = {1.0,1.0,0.01}; for (int w = 0; w < nw; ++w) { WellControls* wc = wells_->ctrls[w]; @@ -136,8 +141,12 @@ namespace Opm 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 ; + wellSolutions()[2*nw + w] = g[Gas] * (wellRates()[np*w + gaspos] - solventWellRate(w)) / total_rates ; } + if( pu.has_solvent) { + wellSolutions()[3*nw + w] = g[Gas] * solventWellRate(w) / total_rates; + } + } else { @@ -147,6 +156,10 @@ namespace Opm if( pu.phase_used[Gas] ) { wellSolutions()[2*nw + w] = wells_->comp_frac[np*w + gaspos]; } + if(pu.has_solvent) { + wellSolutions()[3*nw + w] = 0; + } + } } } @@ -163,14 +176,44 @@ namespace Opm std::vector& wellSolutions() { return well_solutions_; } const std::vector& wellSolutions() const { return well_solutions_; } + /// One rate pr well connection. + std::vector& perfRateSolvent() { return perfRateSolvent_; } + const std::vector& perfRateSolvent() const { return perfRateSolvent_; } + + /// One rate pr well + const double solventWellRate(const int w) const { + double solvent_well_rate = 0.0; + for (int perf = wells_->well_connpos[w]; perf < wells_->well_connpos[w+1]; ++perf ) { + solvent_well_rate += perfRateSolvent_[perf]; + } + return solvent_well_rate; + } + + data::Wells report(const PhaseUsage& pu) const override { data::Wells res = BaseType::report(pu); + const int nw = WellState::numWells(); + // If there are now wells numPhases throws a floating point + // exception. + if (nw == 0) { + return res; + } + if(pu.has_solvent) { + // add solvent component + for( int w = 0; w < nw; ++w ) { + using rt = data::Rates::opt; + res.at( wells_->name[ w ]).rates.set( rt::solvent, solventWellRate(w) ); + } + } + return res; } private: std::vector well_solutions_; + std::vector perfRateSolvent_; + }; } // namespace Opm From 50c1a1404a17a13befaa0fa3fcaa08b06e6939a0 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Wed, 10 May 2017 11:16:28 +0200 Subject: [PATCH 2/6] Remove BlackoilSolventState The solvent saturation is added to BlackoilState and the BlackoilSolventState is thus redundant. --- CMakeLists_files.cmake | 2 - opm/autodiff/BlackoilSolventModel.hpp | 3 +- opm/autodiff/BlackoilSolventModel_impl.hpp | 2 +- opm/autodiff/BlackoilSolventState.cpp | 34 -------------- opm/autodiff/BlackoilSolventState.hpp | 46 ------------------- opm/autodiff/Compat.cpp | 3 +- .../SimulatorFullyImplicitBlackoilOutput.hpp | 2 +- .../SimulatorFullyImplicitBlackoilSolvent.hpp | 3 +- 8 files changed, 5 insertions(+), 90 deletions(-) delete mode 100644 opm/autodiff/BlackoilSolventState.cpp delete mode 100644 opm/autodiff/BlackoilSolventState.hpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 2392fa88d..fcabf8bb6 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -48,7 +48,6 @@ list (APPEND MAIN_SOURCE_FILES opm/autodiff/VFPInjProperties.cpp opm/autodiff/WellMultiSegment.cpp opm/autodiff/MultisegmentWells.cpp - opm/autodiff/BlackoilSolventState.cpp opm/autodiff/MissingFeatures.cpp opm/polymer/PolymerState.cpp opm/polymer/PolymerBlackoilState.cpp @@ -169,7 +168,6 @@ list (APPEND PUBLIC_HEADER_FILES opm/autodiff/BlackoilSequentialModel.hpp opm/autodiff/BlackoilSolventModel.hpp opm/autodiff/BlackoilSolventModel_impl.hpp - opm/autodiff/BlackoilSolventState.hpp opm/autodiff/BlackoilMultiSegmentModel.hpp opm/autodiff/BlackoilMultiSegmentModel_impl.hpp opm/autodiff/BlackoilTransportModel.hpp diff --git a/opm/autodiff/BlackoilSolventModel.hpp b/opm/autodiff/BlackoilSolventModel.hpp index 2c56f4fbe..87606e24a 100644 --- a/opm/autodiff/BlackoilSolventModel.hpp +++ b/opm/autodiff/BlackoilSolventModel.hpp @@ -22,7 +22,6 @@ #include #include -#include #include #include #include @@ -255,7 +254,7 @@ namespace Opm { template struct ModelTraits< BlackoilSolventModel > { - typedef BlackoilSolventState ReservoirState; + typedef BlackoilState ReservoirState; typedef WellStateFullyImplicitBlackoilSolvent WellState; typedef BlackoilModelParameters ModelParameters; typedef BlackoilSolventSolutionState SolutionState; diff --git a/opm/autodiff/BlackoilSolventModel_impl.hpp b/opm/autodiff/BlackoilSolventModel_impl.hpp index 807bb6fbf..422c67eda 100644 --- a/opm/autodiff/BlackoilSolventModel_impl.hpp +++ b/opm/autodiff/BlackoilSolventModel_impl.hpp @@ -135,7 +135,7 @@ namespace Opm { // Initial solvent concentration. if (has_solvent_) { - const auto& solvent_saturation = x.getCellData( BlackoilSolventState::SSOL ); + const auto& solvent_saturation = x.getCellData( BlackoilState::SSOL ); const int nc = solvent_saturation.size(); const V ss = Eigen::Map(solvent_saturation.data() , nc); diff --git a/opm/autodiff/BlackoilSolventState.cpp b/opm/autodiff/BlackoilSolventState.cpp deleted file mode 100644 index 90e8bccd7..000000000 --- a/opm/autodiff/BlackoilSolventState.cpp +++ /dev/null @@ -1,34 +0,0 @@ -/* - Copyright 2015 IRIS AS, Applied Mathematics. - - 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 . -*/ - - -#include - -namespace Opm -{ - const std::string BlackoilSolventState::SSOL = "SSOL"; - - BlackoilSolventState::BlackoilSolventState( int number_of_cells, int number_of_faces, int number_of_phases) - : BlackoilState( number_of_cells , number_of_faces , number_of_phases) - { - registerCellData( SSOL , 1 ); - } - -} // namespace Opm - diff --git a/opm/autodiff/BlackoilSolventState.hpp b/opm/autodiff/BlackoilSolventState.hpp deleted file mode 100644 index 0103e0904..000000000 --- a/opm/autodiff/BlackoilSolventState.hpp +++ /dev/null @@ -1,46 +0,0 @@ -/* - Copyright 2015 IRIS AS, Applied Mathematics. - - 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 . -*/ - -#ifndef OPM_BLACKOILSOLVENTSTATE_HEADER_INCLUDED -#define OPM_BLACKOILSOLVENTSTATE_HEADER_INCLUDED - -#include - - #include - - -namespace Opm -{ - - /// Simulator state for blackoil simulator with solvent. - /// We use the Blackoil state parameters. - class BlackoilSolventState : public BlackoilState - { - public: - static const std::string SSOL; - - BlackoilSolventState( int number_of_cells, int number_of_faces, int number_of_phases); - }; - -} // namespace Opm - - - - -#endif // OPM_BLACKOILSOLVENTSTATE_HEADER_INCLUDED diff --git a/opm/autodiff/Compat.cpp b/opm/autodiff/Compat.cpp index 6a5735394..7dbc4746d 100644 --- a/opm/autodiff/Compat.cpp +++ b/opm/autodiff/Compat.cpp @@ -29,7 +29,6 @@ #include #include #include -#include #include #include #include @@ -111,7 +110,7 @@ data::Solution simToSolution( const SimulationDataContainer& reservoir, } if (phases.has_solvent) { - sol.insert( "SSOL", UnitSystem::measure::identity, reservoir.getCellData( BlackoilSolventState::SSOL ) , data::TargetType::RESTART_SOLUTION ); + sol.insert( "SSOL", UnitSystem::measure::identity, reservoir.getCellData( BlackoilState::SSOL ) , data::TargetType::RESTART_SOLUTION ); } return sol; diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp index 49fefbbc0..871bc9e98 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp @@ -62,7 +62,7 @@ namespace Opm { class SimulationDataContainer; - class BlackoilSolventState; + class BlackoilState; void outputStateVtk(const UnstructuredGrid& grid, const Opm::SimulationDataContainer& state, diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilSolvent.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilSolvent.hpp index 4bc7ee446..45c1ac6cc 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilSolvent.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilSolvent.hpp @@ -22,7 +22,6 @@ #include #include -#include #include #include @@ -85,7 +84,7 @@ namespace Opm struct SimulatorTraits > { typedef WellStateFullyImplicitBlackoilSolvent WellState; - typedef BlackoilSolventState ReservoirState; + typedef BlackoilState ReservoirState; typedef BlackoilOutputWriter OutputWriter; typedef GridT Grid; typedef BlackoilSolventModel Model; From 441a8895ac5d6381470e5ea040ad568dada02c67 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Tue, 23 May 2017 15:46:11 +0200 Subject: [PATCH 3/6] Improvments for solvent model - add dss to appleyard chopping - support for bhp injectors with solvent - copy perfSolventRates between the time steps. - fix bug in well access indicies when numComponents ~= numPhases --- opm/autodiff/BlackoilModelEbos.hpp | 2 + opm/autodiff/StandardWellsDense_impl.hpp | 37 +++++++++------ .../WellStateFullyImplicitBlackoilDense.hpp | 46 +++++++++++++++++-- 3 files changed, 67 insertions(+), 18 deletions(-) diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 3905786f2..29109fab7 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -633,6 +633,8 @@ namespace Opm { // Appleyard chop process. maxVal = std::max(std::abs(dsg),maxVal); + maxVal = std::max(std::abs(dss),maxVal); + double step = dsMax()/maxVal; step = std::min(step, 1.0); diff --git a/opm/autodiff/StandardWellsDense_impl.hpp b/opm/autodiff/StandardWellsDense_impl.hpp index 751f09832..3517518c5 100644 --- a/opm/autodiff/StandardWellsDense_impl.hpp +++ b/opm/autodiff/StandardWellsDense_impl.hpp @@ -1157,14 +1157,14 @@ namespace Opm { if (pu.phase_used[BlackoilPhases::Vapour]) { const int gaspos = pu.phase_pos[BlackoilPhases::Vapour] + perf * numComp; - const int gaspos_well = pu.phase_pos[BlackoilPhases::Vapour] + w * numComp; + const int gaspos_well = pu.phase_pos[BlackoilPhases::Vapour] + w * pu.num_phases; if (pu.phase_used[BlackoilPhases::Liquid]) { - const int oilpos_well = pu.phase_pos[BlackoilPhases::Liquid] + w * numComp; + const int oilpos_well = pu.phase_pos[BlackoilPhases::Liquid] + w * pu.num_phases; const double oilrate = std::abs(xw.wellRates()[oilpos_well]); //in order to handle negative rates in producers rvmax_perf[perf] = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg); if (oilrate > 0) { - const double gasrate = std::abs(xw.wellRates()[gaspos_well]); + const double gasrate = std::abs(xw.wellRates()[gaspos_well]) - xw.solventWellRate(w); double rv = 0.0; if (gasrate > 0) { rv = oilrate / gasrate; @@ -1184,11 +1184,11 @@ namespace Opm { if (pu.phase_used[BlackoilPhases::Liquid]) { const int oilpos = pu.phase_pos[BlackoilPhases::Liquid] + perf * numComp; - const int oilpos_well = pu.phase_pos[BlackoilPhases::Liquid] + w * numComp; + const int oilpos_well = pu.phase_pos[BlackoilPhases::Liquid] + w * pu.num_phases; if (pu.phase_used[BlackoilPhases::Vapour]) { rsmax_perf[perf] = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), temperature, p_avg); - const int gaspos_well = pu.phase_pos[BlackoilPhases::Vapour] + w * numComp; - const double gasrate = std::abs(xw.wellRates()[gaspos_well]); + const int gaspos_well = pu.phase_pos[BlackoilPhases::Vapour] + w * pu.num_phases; + const double gasrate = std::abs(xw.wellRates()[gaspos_well]) - xw.solventWellRate(w); if (gasrate > 0) { const double oilrate = std::abs(xw.wellRates()[oilpos_well]); double rs = 0.0; @@ -2149,17 +2149,24 @@ namespace Opm { // surface rate injection control. Improvement will be required. if (wells().type[wellIdx] == INJECTOR) { if (has_solvent_ ) { + double comp_frac = 0.0; + if (compIdx == solventCompIdx) { // solvent + comp_frac = wells().comp_frac[np*wellIdx + pu.phase_pos[ Gas ]] * wsolvent(wellIdx); + } else if (compIdx == pu.phase_pos[ Gas ]) { + comp_frac = wells().comp_frac[np*wellIdx + compIdx] * (1.0 - wsolvent(wellIdx)); + } else { + comp_frac = wells().comp_frac[np*wellIdx + compIdx]; + } + if (comp_frac == 0.0) { + return qs; //zero + } + if (well_controls_get_current_type(wc) == BHP || well_controls_get_current_type(wc) == THP) { - OPM_THROW(std::runtime_error,"BHP controlled solvent injector is unsupported. Check well " - << wells().name [wellIdx] ); - } - if (compIdx == pu.phase_pos[ Gas ]) { //gas - qs.setValue(target_rate * (1.0 - wsolvent(wellIdx))); - return qs; - } else if (compIdx == solventCompIdx) { // solvent - qs.setValue(wsolvent(wellIdx) * target_rate); - return qs; + return comp_frac * wellVariables_[nw*XvarWell + wellIdx]; } + + qs.setValue(comp_frac * target_rate); + return qs; } const double comp_frac = wells().comp_frac[np*wellIdx + compIdx]; if (comp_frac == 0.0) { diff --git a/opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp b/opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp index fc1ef92c2..f28d79d8b 100644 --- a/opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp +++ b/opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp @@ -67,6 +67,48 @@ namespace Opm // call init on base class BaseType :: init(wells, state, prevState); + + const int nw = wells->number_of_wells; + if (nw == 0) { + return; + } + const int nperf = wells->well_connpos[nw]; + perfRateSolvent_.clear(); + perfRateSolvent_.resize(nperf, 0.0); + + if(pu.has_solvent) { + + // intialize wells that have been there before + // order may change so the mapping is based on the well name + if( ! prevState.wellMap().empty() ) + { + typedef typename WellMapType :: const_iterator const_iterator; + const_iterator end = prevState.wellMap().end(); + for (int w = 0; w < nw; ++w) { + std::string name( wells->name[ w ] ); + const_iterator it = prevState.wellMap().find( name ); + if( it != end ) + { + const int newIndex = w; + + // perfSolventRates + int oldPerf_idx = (*it).second[ 1 ]; + const int num_perf_old_well = (*it).second[ 2 ]; + const int num_perf_this_well = wells->well_connpos[newIndex + 1] - wells->well_connpos[newIndex]; + if( num_perf_old_well == num_perf_this_well ) + { + for (int perf = wells->well_connpos[ newIndex ]; + perf < wells->well_connpos[ newIndex + 1]; ++perf, ++oldPerf_idx ) + { + perfRateSolvent()[ perf ] = prevState.perfRateSolvent()[ oldPerf_idx ]; + } + } + } + } + } + } + + // TODO: the reason to keep this is to avoid getting defaulted value BHP // some facilities needed from opm-parser or opm-core // It is a little tricky, since sometimes before applying group control, the only @@ -88,9 +130,7 @@ namespace Opm if (nw == 0) { return; } - const int nperf = wells_->well_connpos[nw]; - perfRateSolvent_.clear(); - perfRateSolvent_.resize(nperf, 0.0); + const int np = wells_->number_of_phases; const int numComp = pu.has_solvent? np+1:np; From 683ba7e9298fbd95a75e086b95259eb463dfe07b Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Tue, 23 May 2017 16:08:08 +0200 Subject: [PATCH 4/6] Fix fallout from rebase --- opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp index ffa5ad684..277d3cac9 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp @@ -253,7 +253,7 @@ public: solver_timer.start(); const auto& wells_ecl = eclState().getSchedule().getWells(timer.currentStepNum()); - const WellModel well_model(wells, &(wells_manager.wellCollection()), wells_ecl, model_param_, terminal_output_, timer.currentStepNum()); + WellModel well_model(wells, &(wells_manager.wellCollection()), wells_ecl, model_param_, terminal_output_, timer.currentStepNum()); auto solver = createSolver(well_model); From 8b75e2eedcfd6bc0d497f2600d9b7802c19a24a1 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Wed, 24 May 2017 10:53:48 +0200 Subject: [PATCH 5/6] Minor fixes solvent model 1) Fix GRAT controlled wells 2) Fix seg_fault in computeWellPotentials for solvent simulations 3) Fix mem_issus for RESV combinded with solvent --- opm/autodiff/StandardWellsDense_impl.hpp | 23 +++++++++++++++++++---- 1 file changed, 19 insertions(+), 4 deletions(-) diff --git a/opm/autodiff/StandardWellsDense_impl.hpp b/opm/autodiff/StandardWellsDense_impl.hpp index 3517518c5..1bd934594 100644 --- a/opm/autodiff/StandardWellsDense_impl.hpp +++ b/opm/autodiff/StandardWellsDense_impl.hpp @@ -2213,17 +2213,27 @@ namespace Opm { assert(phase_under_control >= 0); + EvalWell wellVolumeFractionScaledPhaseUnderControl = wellVolumeFractionScaled(wellIdx, phase_under_control); + if (has_solvent_ && phase_under_control == Gas) { + // for GRAT controlled wells solvent is included in the target + wellVolumeFractionScaledPhaseUnderControl += wellVolumeFractionScaled(wellIdx, solventCompIdx); + } + if (compIdx == phase_under_control) { + if (has_solvent_ && phase_under_control == Gas) { + qs.setValue(target_rate * wellVolumeFractionScaled(wellIdx, Gas).value() / wellVolumeFractionScaledPhaseUnderControl.value() ); + return qs; + } qs.setValue(target_rate); return qs; } // TODO: not sure why the single phase under control will have near zero fraction const double eps = 1e-6; - if (wellVolumeFractionScaled(wellIdx, phase_under_control) < eps) { + if (wellVolumeFractionScaledPhaseUnderControl < eps) { return qs; } - return (target_rate * wellVolumeFractionScaled(wellIdx, compIdx) / wellVolumeFractionScaled(wellIdx, phase_under_control)); + return (target_rate * wellVolumeFractionScaled(wellIdx, compIdx) / wellVolumeFractionScaledPhaseUnderControl); } // when it is a combined two phase rate limit, such like LRAT @@ -2301,6 +2311,10 @@ namespace Opm { { const WellControls* wc = wells().ctrls[wellIdx]; if (well_controls_get_current_type(wc) == RESERVOIR_RATE) { + + if (has_solvent_ && compIdx == solventCompIdx) { + return wellVolumeFraction(wellIdx, compIdx); + } const double* distr = well_controls_get_current_distr(wc); assert(compIdx < 3); if (distr[compIdx] > 0.) { @@ -2780,6 +2794,7 @@ namespace Opm { std::vector& well_flux) const { const int np = wells().number_of_phases; + const int numComp = numComponents(); well_flux.resize(np, 0.0); const bool allow_cf = allow_cross_flow(well_index, ebosSimulator); @@ -2787,8 +2802,8 @@ namespace Opm { const int cell_index = wells().well_cells[perf]; const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_index, /*timeIdx=*/ 0)); // flux for each perforation - std::vector cq_s(np, 0.0); - std::vector mob(np, 0.0); + std::vector cq_s(numComp, 0.0); + std::vector mob(numComp, 0.0); getMobility(ebosSimulator, perf, cell_index, mob); computeWellFlux(well_index, wells().WI[perf], intQuants, mob, bhp, wellPerforationPressureDiffs()[perf], allow_cf, cq_s); From f671af6cd61074cb6543c3610f4880882359a489 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Tue, 30 May 2017 14:33:17 +0200 Subject: [PATCH 6/6] Clean-up of the solvent implementation --- opm/autodiff/BlackoilModelEbos.hpp | 3 +-- opm/autodiff/StandardWellsDense.hpp | 3 ++- opm/autodiff/StandardWellsDense_impl.hpp | 7 ++++--- opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp | 6 +++--- 4 files changed, 10 insertions(+), 9 deletions(-) diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 29109fab7..797d67f81 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -934,8 +934,7 @@ namespace Opm { Vector& R2_idx = R2[ solventCompIdx ]; Vector& B_idx = B[ solventCompIdx ]; B_idx [cell_idx] = 1.0 / intQuants.solventInverseFormationVolumeFactor().value(); - const int ebosCompIdx = flowPhaseToEbosCompIdx(solventCompIdx); - R2_idx[cell_idx] = ebosResid[cell_idx][ebosCompIdx]; + R2_idx[cell_idx] = ebosResid[cell_idx][solventCompIdx]; } } diff --git a/opm/autodiff/StandardWellsDense.hpp b/opm/autodiff/StandardWellsDense.hpp index 6c862d6d0..8ba09e20e 100644 --- a/opm/autodiff/StandardWellsDense.hpp +++ b/opm/autodiff/StandardWellsDense.hpp @@ -126,8 +126,9 @@ enum WellVariablePositions { return 2; } int numComp = FluidSystem::numComponents; - if (has_solvent_) + if (has_solvent_) { numComp ++; + } return numComp; } diff --git a/opm/autodiff/StandardWellsDense_impl.hpp b/opm/autodiff/StandardWellsDense_impl.hpp index 1bd934594..11781c98b 100644 --- a/opm/autodiff/StandardWellsDense_impl.hpp +++ b/opm/autodiff/StandardWellsDense_impl.hpp @@ -1210,7 +1210,7 @@ namespace Opm { surf_dens_perf[numComp*perf + p] = FluidSystem::referenceDensity( flowPhaseToEbosPhaseIdx( p ), fs.pvtRegionIndex()); } - #warning HACK use cell values for solvent injector + // We use cell values for solvent injector if (has_solvent_) { b_perf[numComp*perf + solventCompIdx] = intQuants.solventInverseFormationVolumeFactor().value(); surf_dens_perf[numComp*perf + solventCompIdx] = intQuants.solventRefDensity(); @@ -1325,7 +1325,8 @@ namespace Opm { well_state.wellSolutions()[SFrac*nw + w] = F_solvent; } - # warning F_solvent is added to F_gas. This means that well_rate[Gas] also contains solvent. More testing is needed to make sure this is correct for output, wellControls, well groups, THP etc. + // F_solvent is added to F_gas. This means that well_rate[Gas] also contains solvent. + // More testing is needed to make sure this is correct for well groups and THP. if (has_solvent_){ F[Gas] += F_solvent; } @@ -3029,7 +3030,7 @@ namespace Opm { } } // we didn't find it return 0; - // or should we throw? + assert(false); return 0.0; } diff --git a/opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp b/opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp index f28d79d8b..840c7ed33 100644 --- a/opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp +++ b/opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp @@ -76,7 +76,7 @@ namespace Opm perfRateSolvent_.clear(); perfRateSolvent_.resize(nperf, 0.0); - if(pu.has_solvent) { + if (pu.has_solvent) { // intialize wells that have been there before // order may change so the mapping is based on the well name @@ -196,7 +196,7 @@ namespace Opm if( pu.phase_used[Gas] ) { wellSolutions()[2*nw + w] = wells_->comp_frac[np*w + gaspos]; } - if(pu.has_solvent) { + if (pu.has_solvent) { wellSolutions()[3*nw + w] = 0; } @@ -238,7 +238,7 @@ namespace Opm if (nw == 0) { return res; } - if(pu.has_solvent) { + if (pu.has_solvent) { // add solvent component for( int w = 0; w < nw; ++w ) { using rt = data::Rates::opt;