From 0068c175a7dc05ed5e79d04573d4093796bef91d Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Wed, 7 Jun 2017 09:29:31 +0200 Subject: [PATCH 1/6] Add polymer option to flow_ebos No extra equation is added for polymer in the well equation. Seperate executables are added for polymer: flow_ebos_polymer and solvent: flow_ebos_solvent Tested and verified on the test cases in polymer_test_suite This PR should not effect the performance and results of the blackoil simulator --- CMakeLists_files.cmake | 4 + examples/flow_ebos.cpp | 2 +- examples/flow_ebos_polymer.cpp | 42 +++ examples/flow_ebos_solvent.cpp | 42 +++ opm/autodiff/BlackoilModelEbos.hpp | 84 ++++- opm/autodiff/Compat.cpp | 11 + opm/autodiff/FlowMainEbos.hpp | 27 +- opm/autodiff/FlowMainPolymer.hpp | 4 - .../SimulatorFullyImplicitBlackoilEbos.hpp | 22 +- opm/autodiff/StandardWellsDense.hpp | 31 +- opm/autodiff/StandardWellsDense_impl.hpp | 310 +++++++++++++++--- 11 files changed, 486 insertions(+), 93 deletions(-) create mode 100644 examples/flow_ebos_polymer.cpp create mode 100644 examples/flow_ebos_solvent.cpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 499067468..47209bce9 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -112,6 +112,8 @@ list (APPEND EXAMPLE_SOURCE_FILES examples/flow_legacy.cpp examples/flow_sequential.cpp examples/flow_ebos.cpp + examples/flow_ebos_solvent.cpp + examples/flow_ebos_polymer.cpp examples/flow_multisegment.cpp examples/flow_solvent.cpp examples/sim_2p_incomp.cpp @@ -132,6 +134,8 @@ list (APPEND PROGRAM_SOURCE_FILES examples/sim_2p_incomp_ad.cpp examples/sim_2p_comp_reorder.cpp examples/flow_ebos.cpp + examples/flow_ebos_solvent.cpp + examples/flow_ebos_polymer.cpp examples/flow_legacy.cpp examples/flow_sequential.cpp examples/flow_solvent.cpp diff --git a/examples/flow_ebos.cpp b/examples/flow_ebos.cpp index 9ff12a482..90c44433f 100644 --- a/examples/flow_ebos.cpp +++ b/examples/flow_ebos.cpp @@ -34,6 +34,6 @@ // ----------------- Main program ----------------- int main(int argc, char** argv) { - Opm::FlowMainEbos mainfunc; + Opm::FlowMainEbos mainfunc; return mainfunc.execute(argc, argv); } diff --git a/examples/flow_ebos_polymer.cpp b/examples/flow_ebos_polymer.cpp new file mode 100644 index 000000000..d09a0d759 --- /dev/null +++ b/examples/flow_ebos_polymer.cpp @@ -0,0 +1,42 @@ +/* + Copyright 2017 IRIS AS + + 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 . +*/ + + +#if HAVE_CONFIG_H +#include "config.h" +#endif // HAVE_CONFIG_H + +#include +#include +#include +#include +#include + +namespace Ewoms { +namespace Properties { +NEW_TYPE_TAG(EclFlowPolymerProblem, INHERITS_FROM(EclFlowProblem)); +SET_BOOL_PROP(EclFlowPolymerProblem, EnablePolymer, true); +}} + +// ----------------- Main program ----------------- +int main(int argc, char** argv) +{ + Opm::FlowMainEbos mainfunc; + return mainfunc.execute(argc, argv); +} diff --git a/examples/flow_ebos_solvent.cpp b/examples/flow_ebos_solvent.cpp new file mode 100644 index 000000000..443c400f8 --- /dev/null +++ b/examples/flow_ebos_solvent.cpp @@ -0,0 +1,42 @@ +/* + Copyright 2017 IRIS AS + + 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 . +*/ + + +#if HAVE_CONFIG_H +#include "config.h" +#endif // HAVE_CONFIG_H + +#include +#include +#include +#include +#include + +namespace Ewoms { +namespace Properties { +NEW_TYPE_TAG(EclFlowSolventProblem, INHERITS_FROM(EclFlowProblem)); +SET_BOOL_PROP(EclFlowSolventProblem, EnableSolvent, true); +}} + +// ----------------- Main program ----------------- +int main(int argc, char** argv) +{ + Opm::FlowMainEbos mainfunc; + return mainfunc.execute(argc, argv); +} diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 9a7d45e24..4d48e28bc 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -83,7 +83,6 @@ 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. @@ -97,6 +96,7 @@ namespace Opm { /// where gas can be dissolved in oil and vice versa. It /// uses an industry-standard TPFA discretization with per-phase /// upwind weighting of mobilities. + template class BlackoilModelEbos { public: @@ -105,7 +105,6 @@ namespace Opm { typedef WellStateFullyImplicitBlackoilDense WellState; typedef BlackoilModelParameters ModelParameters; - typedef typename TTAG(EclFlowProblem) TypeTag; typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; @@ -118,7 +117,11 @@ namespace Opm { typedef double Scalar; static const int numEq = BlackoilIndices::numEq; - static const int solventCompIdx = 3; //TODO get this from ebos + static const int contiSolventEqIdx = BlackoilIndices::contiSolventEqIdx; + static const int contiPolymerEqIdx = BlackoilIndices::contiPolymerEqIdx; + static const int solventSaturationIdx = BlackoilIndices::solventSaturationIdx; + static const int polymerConcentrationIdx = BlackoilIndices::polymerConcentrationIdx; + typedef Dune::FieldVector VectorBlockType; typedef Dune::FieldMatrix MatrixBlockType; typedef Dune::BCRSMatrix Mat; @@ -149,7 +152,8 @@ namespace Opm { const ModelParameters& param, const StandardWellsDense& well_model, const NewtonIterationBlackoilInterface& linsolver, - const bool terminal_output) + const bool terminal_output + ) : ebosSimulator_(ebosSimulator) , grid_(ebosSimulator_.gridManager().grid()) , istlSolver_( dynamic_cast< const ISTLSolverType* > (&linsolver) ) @@ -161,8 +165,9 @@ namespace Opm { , has_disgas_(FluidSystem::enableDissolvedGas()) , has_vapoil_(FluidSystem::enableVaporizedOil()) , has_solvent_(GET_PROP_VALUE(TypeTag, EnableSolvent)) + , has_polymer_(GET_PROP_VALUE(TypeTag, EnablePolymer)) , param_( param ) - , well_model_ (well_model) + , well_model_ (well_model) , terminal_output_ (terminal_output) , rate_converter_(phaseUsage_, ebosSimulator_.problem().pvtRegionArray().empty()?nullptr:ebosSimulator_.problem().pvtRegionArray().data(), AutoDiffGrid::numCells(grid_), std::vector(AutoDiffGrid::numCells(grid_),0)) , current_relaxation_(1.0) @@ -615,6 +620,9 @@ namespace Opm { const double dss = has_solvent_ ? dx[cell_idx][BlackoilIndices::solventSaturationIdx] : 0.0; dso -= dss; + // polymer + const double dc = has_polymer_ ? dx[cell_idx][BlackoilIndices::polymerConcentrationIdx] : 0.0; + // Appleyard chop process. maxVal = std::max(std::abs(dsg),maxVal); maxVal = std::max(std::abs(dss),maxVal); @@ -637,6 +645,12 @@ namespace Opm { ss -= step * dss; } + if (has_polymer_) { + double& c = reservoir_state.getCellData( reservoir_state.POLYMER )[cell_idx]; + c -= step * dc; + c = std::max(c, 0.0); + } + double& so = reservoir_state.saturation()[cell_idx*np + pu.phase_pos[ Oil ]]; so -= step * dso; @@ -867,10 +881,16 @@ namespace Opm { } if ( has_solvent_ ) { - B_avg[ solventCompIdx ] += 1.0 / intQuants.solventInverseFormationVolumeFactor().value(); - const auto R2 = ebosResid[cell_idx][solventCompIdx]; - R_sum[ solventCompIdx ] += R2; - maxCoeff[ solventCompIdx ] = std::max( maxCoeff[ solventCompIdx ], std::abs( R2 ) / pvValue ); + B_avg[ contiSolventEqIdx ] += 1.0 / intQuants.solventInverseFormationVolumeFactor().value(); + const auto R2 = ebosResid[cell_idx][contiSolventEqIdx]; + R_sum[ contiSolventEqIdx ] += R2; + maxCoeff[ contiSolventEqIdx ] = std::max( maxCoeff[ contiSolventEqIdx ], std::abs( R2 ) / pvValue ); + } + if (has_polymer_ ) { + B_avg[ contiPolymerEqIdx ] += 1.0 / intQuants.solventInverseFormationVolumeFactor().value(); + const auto R2 = ebosResid[cell_idx][contiPolymerEqIdx]; + R_sum[ contiPolymerEqIdx ] += R2; + maxCoeff[ contiPolymerEqIdx ] = std::max( maxCoeff[ contiPolymerEqIdx ], std::abs( R2 ) / pvValue ); } } @@ -938,7 +958,11 @@ namespace Opm { key[ phaseIdx ] = std::toupper( phaseName.front() ); } if (has_solvent_) { - key[ solventCompIdx ] = "S"; + key[ solventSaturationIdx ] = "S"; + } + + if (has_polymer_) { + key[ polymerConcentrationIdx ] = "P"; } for (int compIdx = 0; compIdx < numComp; ++compIdx) { @@ -1004,6 +1028,9 @@ namespace Opm { if (has_solvent_) numComp ++; + if (has_polymer_) + numComp ++; + return numComp; } @@ -1450,6 +1477,7 @@ namespace Opm { const bool has_disgas_; const bool has_vapoil_; const bool has_solvent_; + const bool has_polymer_; ModelParameters param_; SimulatorReport failureReport_; @@ -1513,6 +1541,11 @@ namespace Opm { cellPv[BlackoilIndices::solventSaturationIdx] = reservoirState.getCellData( reservoirState.SSOL )[cellIdx]; } + if (has_polymer_) { + cellPv[BlackoilIndices::polymerConcentrationIdx] = reservoirState.getCellData( reservoirState.POLYMER )[cellIdx]; + } + + // set switching variable and interpretation if (active_[Gas] ) { if( reservoirState.hydroCarbonState()[cellIdx] == HydroCarbonState::OilOnly && has_disgas_ ) @@ -1588,18 +1621,23 @@ namespace Opm { int flowToEbosPvIdx( const int flowPv ) const { - const int flowToEbos[ 4 ] = { + const int flowToEbos[ 3 ] = { BlackoilIndices::pressureSwitchIdx, BlackoilIndices::waterSaturationIdx, - BlackoilIndices::compositionSwitchIdx, - BlackoilIndices::solventSaturationIdx + BlackoilIndices::compositionSwitchIdx }; + + if (flowPv > 2 ) + return flowPv; + return flowToEbos[ flowPv ]; } int flowPhaseToEbosCompIdx( const int phaseIdx ) const { - const int phaseToComp[ 4 ] = { FluidSystem::waterCompIdx, FluidSystem::oilCompIdx, FluidSystem::gasCompIdx, solventCompIdx }; + const int phaseToComp[ 3 ] = { FluidSystem::waterCompIdx, FluidSystem::oilCompIdx, FluidSystem::gasCompIdx}; + if (phaseIdx > 2 ) + return phaseIdx; return phaseToComp[ phaseIdx ]; } @@ -1636,8 +1674,11 @@ namespace Opm { // 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; + cellRes[ contiSolventEqIdx ] /= refDens; + cellRes[ contiSolventEqIdx ] *= cellVolume; + } + if (has_polymer_) { + cellRes[ contiPolymerEqIdx ] *= cellVolume; } } @@ -1670,10 +1711,17 @@ namespace Opm { const auto& refDens = intQuants->solventRefDensity(); for( int pvIdx=0; pvIdx < numEq; ++pvIdx ) { - (*col)[solventCompIdx][flowToEbosPvIdx(pvIdx)] /= refDens; - (*col)[solventCompIdx][flowToEbosPvIdx(pvIdx)] *= cellVolume; + (*col)[contiSolventEqIdx][flowToEbosPvIdx(pvIdx)] /= refDens; + (*col)[contiSolventEqIdx][flowToEbosPvIdx(pvIdx)] *= cellVolume; } } + if (has_polymer_) { + for( int pvIdx=0; pvIdx < numEq; ++pvIdx ) + { + (*col)[contiPolymerEqIdx][flowToEbosPvIdx(pvIdx)] *= cellVolume; + } + } + } } } diff --git a/opm/autodiff/Compat.cpp b/opm/autodiff/Compat.cpp index 7dbc4746d..53fd72cee 100644 --- a/opm/autodiff/Compat.cpp +++ b/opm/autodiff/Compat.cpp @@ -23,7 +23,9 @@ #include #include +#include +#include #include #include #include @@ -113,6 +115,15 @@ data::Solution simToSolution( const SimulationDataContainer& reservoir, sol.insert( "SSOL", UnitSystem::measure::identity, reservoir.getCellData( BlackoilState::SSOL ) , data::TargetType::RESTART_SOLUTION ); } + if (phases.has_polymer) { + if (reservoir.hasCellData( PolymerBlackoilState::CONCENTRATION )) { // compatibility with legacy polymer + sol.insert( "POLYMER", UnitSystem::measure::identity, reservoir.getCellData( PolymerBlackoilState::CONCENTRATION ) , data::TargetType::RESTART_SOLUTION ); + } else { + sol.insert( "POLYMER", UnitSystem::measure::identity, reservoir.getCellData( BlackoilState::POLYMER ) , data::TargetType::RESTART_SOLUTION ); + } + + } + return sol; } diff --git a/opm/autodiff/FlowMainEbos.hpp b/opm/autodiff/FlowMainEbos.hpp index e0a5c30be..9ba105843 100755 --- a/opm/autodiff/FlowMainEbos.hpp +++ b/opm/autodiff/FlowMainEbos.hpp @@ -61,6 +61,7 @@ namespace Opm { // The FlowMain class is the ebos based black-oil simulator. + template class FlowMainEbos { enum FileOutputValue{ @@ -73,7 +74,6 @@ namespace Opm }; public: - typedef TTAG(EclFlowProblem) TypeTag; typedef typename GET_PROP(TypeTag, MaterialLaw)::EclMaterialLawManager MaterialLawManager; typedef typename GET_PROP_TYPE(TypeTag, Simulator) EbosSimulator; typedef typename GET_PROP_TYPE(TypeTag, ElementMapper) ElementMapper; @@ -82,7 +82,7 @@ namespace Opm typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; - typedef Opm::SimulatorFullyImplicitBlackoilEbos Simulator; + typedef Opm::SimulatorFullyImplicitBlackoilEbos Simulator; typedef typename Simulator::ReservoirState ReservoirState; typedef typename Simulator::OutputWriter OutputWriter; @@ -534,6 +534,27 @@ namespace Opm } initHydroCarbonState(*state_, pu, Opm::UgGridHelpers::numCells(grid), deck().hasKeyword("DISGAS"), deck().hasKeyword("VAPOIL")); + + // Get initial polymer concentration from ebos + if (GET_PROP_VALUE(TypeTag, EnablePolymer)) { + auto& cpolymer = state_->getCellData( state_->POLYMER ); + const int numCells = Opm::UgGridHelpers::numCells(grid); + for (int c = 0; c < numCells; ++c) { + cpolymer[c] = ebosProblem().polymerConcentration(c); + } + } + // Get initial solvent saturation from ebos + if (GET_PROP_VALUE(TypeTag, EnableSolvent)) { + auto& solvent = state_->getCellData( state_->SSOL ); + auto& sat = state_->saturation(); + const int np = props.numPhases(); + const int numCells = Opm::UgGridHelpers::numCells(grid); + for (int c = 0; c < numCells; ++c) { + solvent[c] = ebosProblem().solventSaturation(c); + sat[c * np + pu.phase_pos[Water]]; + } + } + } // Extract messages from parser. @@ -667,7 +688,7 @@ namespace Opm // fis_solver_ void setupLinearSolver() { - typedef typename BlackoilModelEbos :: ISTLSolverType ISTLSolverType; + typedef typename BlackoilModelEbos :: ISTLSolverType ISTLSolverType; extractParallelGridInformationToISTL(grid(), parallel_information_); fis_solver_.reset( new ISTLSolverType( param_, parallel_information_ ) ); diff --git a/opm/autodiff/FlowMainPolymer.hpp b/opm/autodiff/FlowMainPolymer.hpp index 898c7bef0..092482ef6 100644 --- a/opm/autodiff/FlowMainPolymer.hpp +++ b/opm/autodiff/FlowMainPolymer.hpp @@ -89,10 +89,6 @@ namespace Opm } } - - - - // Setup linear solver. // Writes to: // fis_solver_ diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp index 5330d0591..90ff48d8d 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp @@ -42,14 +42,12 @@ namespace Opm { -class SimulatorFullyImplicitBlackoilEbos; -//class StandardWellsDense; /// a simulator for the blackoil model +template class SimulatorFullyImplicitBlackoilEbos { public: - typedef typename TTAG(EclFlowProblem) TypeTag; typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; @@ -58,11 +56,12 @@ public: typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw; + typedef Ewoms::BlackOilPolymerModule PolymerModule; typedef WellStateFullyImplicitBlackoilDense WellState; typedef BlackoilState ReservoirState; typedef BlackoilOutputWriter OutputWriter; - typedef BlackoilModelEbos Model; + typedef BlackoilModelEbos Model; typedef BlackoilModelParameters ModelParameters; typedef NonlinearSolver Solver; typedef StandardWellsDense WellModel; @@ -255,7 +254,8 @@ public: solver_timer.start(); const auto& wells_ecl = eclState().getSchedule().getWells(timer.currentStepNum()); - 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); @@ -413,11 +413,12 @@ public: { return ebosSimulator_.gridManager().grid(); } protected: - void handleAdditionalWellInflow(SimulatorTimer& /* timer */, + void handleAdditionalWellInflow(SimulatorTimer& timer, WellsManager& /* wells_manager */, WellState& /* well_state */, const Wells* /* wells */) - { } + { + } std::unique_ptr createSolver(WellModel& well_model) { @@ -439,7 +440,8 @@ protected: legacyDepth_, legacyPoreVolume_, rateConverter_.get(), - globalNumCells); + globalNumCells, + grid()); auto model = std::unique_ptr(new Model(ebosSimulator_, model_param_, well_model, @@ -852,6 +854,10 @@ protected: } } + + + + // Data. Simulator& ebosSimulator_; diff --git a/opm/autodiff/StandardWellsDense.hpp b/opm/autodiff/StandardWellsDense.hpp index 6b0429b6d..2557e5cde 100644 --- a/opm/autodiff/StandardWellsDense.hpp +++ b/opm/autodiff/StandardWellsDense.hpp @@ -57,6 +57,8 @@ #include +#include + namespace Opm { enum WellVariablePositions { @@ -85,8 +87,11 @@ 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 + static const int numWellEq = GET_PROP_VALUE(TypeTag, EnablePolymer)? 3:numEq; // //numEq; //number of wellEq is only for 3 for polymer + static const int contiSolventEqIdx = BlackoilIndices::contiSolventEqIdx; + static const int contiPolymerEqIdx = BlackoilIndices::contiPolymerEqIdx; + static const int solventSaturationIdx = BlackoilIndices::solventSaturationIdx; + static const int polymerConcentrationIdx = BlackoilIndices::polymerConcentrationIdx; typedef Dune::FieldVector VectorBlockType; typedef Dune::FieldMatrix MatrixBlockType; typedef Dune::BCRSMatrix Mat; @@ -94,8 +99,7 @@ enum WellVariablePositions { typedef DenseAd::Evaluation EvalWell; typedef DenseAd::Evaluation Eval; - - + typedef Ewoms::BlackOilPolymerModule PolymerModule; // For the conversion between the surface volume rate and resrevoir voidage rate using RateConverterType = RateConverter:: @@ -116,7 +120,8 @@ enum WellVariablePositions { const std::vector& depth_arg, const std::vector& pv_arg, const RateConverterType* rate_converter, - long int global_nc); + long int global_nc, + const auto& grid); /// The number of components in the model. @@ -294,6 +299,7 @@ enum WellVariablePositions { ModelParameters param_; bool terminal_output_; bool has_solvent_; + bool has_polymer_; int current_timeIdx_; PhaseUsage phase_usage_; @@ -314,6 +320,13 @@ enum WellVariablePositions { std::vector well_perforation_densities_; std::vector well_perforation_pressure_diffs_; + std::vector wpolymer_; + std::vector wsolvent_; + + std::vector wells_rep_radius_; + std::vector wells_perf_length_; + std::vector wells_bore_diameter_; + std::vector wellVariables_; std::vector F0_; @@ -397,6 +410,14 @@ enum WellVariablePositions { const std::vector& initial_potential) const; double wsolvent(const int well_index) const; + + double wpolymer(const int well_index) const; + + void setupCompressedToCartesian(const int* global_cell, int number_of_cells, std::map& cartesian_to_compressed ) const; + + void computeRepRadiusPerfLength(const auto& grid); + + }; diff --git a/opm/autodiff/StandardWellsDense_impl.hpp b/opm/autodiff/StandardWellsDense_impl.hpp index dfd06cbd0..854aeb868 100644 --- a/opm/autodiff/StandardWellsDense_impl.hpp +++ b/opm/autodiff/StandardWellsDense_impl.hpp @@ -18,6 +18,7 @@ namespace Opm { , param_(param) , terminal_output_(terminal_output) , has_solvent_(GET_PROP_VALUE(TypeTag, EnableSolvent)) + , has_polymer_(GET_PROP_VALUE(TypeTag, EnablePolymer)) , 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) @@ -47,7 +48,8 @@ namespace Opm { const std::vector& depth_arg, const std::vector& pv_arg, const RateConverterType* rate_converter, - long int global_nc) + long int global_nc, + const auto& grid) { // has to be set always for the convergence check! global_nc_ = global_nc; @@ -118,6 +120,14 @@ namespace Opm { // resize temporary class variables Cx_.resize( duneC_.N() ); invDrw_.resize( invDuneD_.N() ); + + if (has_polymer_) + { + if (PolymerModule::hasPlyshlog()) { + computeRepRadiusPerfLength(grid); + } + } + } @@ -197,6 +207,37 @@ namespace Opm { getMobility(ebosSimulator, perf, cell_idx, mob); computeWellFlux(w, wells().WI[perf], intQuants, mob, bhp, wellPerforationPressureDiffs()[perf], allow_cf, cq_s); + if (has_polymer_) { + if (PolymerModule::hasPlyshlog()) { + // compute the well water velocity based on the perforation rates. + double area = 2 * M_PI * wells_rep_radius_[perf] * wells_perf_length_[perf]; + const auto& materialLawManager = ebosSimulator.problem().materialLawManager(); + const auto& scaledDrainageInfo = + materialLawManager->oilWaterScaledEpsInfoDrainage(cell_idx); + const Scalar& Swcr = scaledDrainageInfo.Swcr; + const EvalWell poro = extendEval(intQuants.porosity()); + const EvalWell Sw = extendEval(intQuants.fluidState().saturation(flowPhaseToEbosPhaseIdx(Water))); + // guard against zero porosity and no water + const EvalWell denom = Opm::max( (area * poro * (Sw - Swcr)), 1e-12); + EvalWell waterVelocity = cq_s[ Water ] / denom * extendEval(intQuants.fluidState().invB(flowPhaseToEbosPhaseIdx(Water))); + EvalWell polymerConcentration = extendEval(intQuants.polymerConcentration()); + + if (PolymerModule::hasShrate()) { + // TODO Use the same conversion as for the reservoar equations. + // Need the "permeability" of the well? + // For now use the same formula as in legacy. + waterVelocity *= PolymerModule::shrate( intQuants.pvtRegionIndex() ) / wells_bore_diameter_[perf]; + } + EvalWell shearFactor = PolymerModule::computeShearFactor(polymerConcentration, + intQuants.pvtRegionIndex(), + waterVelocity); + + // modify the mobility with the shear factor and recompute the well fluxes. + mob[ Water ] /= shearFactor; + computeWellFlux(w, wells().WI[perf], intQuants, mob, bhp, wellPerforationPressureDiffs()[perf], allow_cf, cq_s); + } + } + for (int componentIdx = 0; componentIdx < numComp; ++componentIdx) { // the cq_s entering mass balance equations need to consider the efficiency factors. @@ -209,17 +250,16 @@ namespace Opm { } // subtract sum of phase fluxes in the well equations. - resWell_[w][componentIdx] -= cq_s[componentIdx].value(); + resWell_[w][flowPhaseToEbosCompIdx(componentIdx)] -= cq_s[componentIdx].value(); // assemble the jacobians for (int pvIdx = 0; pvIdx < numWellEq; ++pvIdx) { if (!only_wells) { // also need to consider the efficiency factor when manipulating the jacobians. ebosJac[cell_idx][cell_idx][flowPhaseToEbosCompIdx(componentIdx)][flowToEbosPvIdx(pvIdx)] -= cq_s_effective.derivative(pvIdx); - duneB_[w][cell_idx][pvIdx][flowPhaseToEbosCompIdx(componentIdx)] -= cq_s_effective.derivative(pvIdx+numEq); // intput in transformed matrix - duneC_[w][cell_idx][componentIdx][flowToEbosPvIdx(pvIdx)] -= cq_s_effective.derivative(pvIdx); + duneB_[w][cell_idx][flowToEbosPvIdx(pvIdx)][flowPhaseToEbosCompIdx(componentIdx)] -= cq_s_effective.derivative(pvIdx+numEq); // intput in transformed matrix } - invDuneD_[w][w][componentIdx][pvIdx] -= cq_s[componentIdx].derivative(pvIdx+numEq); + invDuneD_[w][w][flowPhaseToEbosCompIdx(componentIdx)][pvIdx] -= cq_s[componentIdx].derivative(pvIdx+numEq); } // add trivial equation for 2p cases (Only support water + oil) @@ -229,7 +269,7 @@ namespace Opm { } // Store the perforation phase flux for later usage. - if (componentIdx == solventCompIdx) {// if (flowPhaseToEbosCompIdx(componentIdx) == Solvent) + if (has_solvent_ && componentIdx == solventSaturationIdx) {// if (flowPhaseToEbosCompIdx(componentIdx) == Solvent) well_state.perfRateSolvent()[perf] = cq_s[componentIdx].value(); } else { well_state.perfPhaseRates()[perf*np + componentIdx] = cq_s[componentIdx].value(); @@ -455,13 +495,31 @@ namespace Opm { + template + int + StandardWellsDense:: + flowToEbosPvIdx( const int flowPv ) const + { + const int flowToEbos[ 3 ] = { + BlackoilIndices::pressureSwitchIdx, + BlackoilIndices::waterSaturationIdx, + BlackoilIndices::compositionSwitchIdx + }; + + if (flowPv > 2 ) + return flowPv; + + return flowToEbos[ flowPv ]; + } template int StandardWellsDense:: flowPhaseToEbosCompIdx( const int phaseIdx ) const { - const int phaseToComp[ 4 ] = { FluidSystem::waterCompIdx, FluidSystem::oilCompIdx, FluidSystem::gasCompIdx, solventCompIdx }; + const int phaseToComp[ 3 ] = { FluidSystem::waterCompIdx, FluidSystem::oilCompIdx, FluidSystem::gasCompIdx}; + if (phaseIdx > 2 ) + return phaseIdx; return phaseToComp[ phaseIdx ]; } @@ -469,24 +527,6 @@ namespace Opm { - template - int - StandardWellsDense:: - flowToEbosPvIdx( const int flowPv ) const - { - const int flowToEbos[ 4 ] = { - BlackoilIndices::pressureSwitchIdx, - BlackoilIndices::waterSaturationIdx, - BlackoilIndices::compositionSwitchIdx, - BlackoilIndices::solventSaturationIdx - }; - return flowToEbos[ flowPv ]; - } - - - - - template int StandardWellsDense:: @@ -497,10 +537,6 @@ namespace Opm { return flowToEbos[ phaseIdx ]; } - - - - template std::vector StandardWellsDense:: @@ -687,8 +723,8 @@ namespace Opm { { const int nw = wells().number_of_wells; // for two-phase numComp < numEq - const int numComp = numComponents(); - for (int eqIdx = 0; eqIdx < numComp; ++eqIdx) { + //const int numComp = numComponents(); + for (int eqIdx = 0; eqIdx < numWellEq; ++eqIdx) { for (int w = 0; w < nw; ++w) { const unsigned int idx = nw * eqIdx + w; assert( idx < wellVariables_.size() ); @@ -697,7 +733,7 @@ namespace Opm { eval = 0.0; eval.setValue( xw.wellSolutions()[ idx ] ); - eval.setDerivative(numWellEq + eqIdx, 1.0); + eval.setDerivative(numEq + eqIdx, 1.0); } } } @@ -765,7 +801,7 @@ namespace Opm { b_perfcells_dense[phase] = extendEval(fs.invB(ebosPhaseIdx)); } if (has_solvent_) { - b_perfcells_dense[solventCompIdx] = extendEval(intQuants.solventInverseFormationVolumeFactor()); + b_perfcells_dense[solventSaturationIdx] = extendEval(intQuants.solventInverseFormationVolumeFactor()); } // Pressure drawdown (also used to determine direction of flow) @@ -817,7 +853,7 @@ namespace Opm { } if (has_solvent_) { - volumeRatio += cmix_s[solventCompIdx] / b_perfcells_dense[solventCompIdx]; + volumeRatio += cmix_s[solventSaturationIdx] / b_perfcells_dense[solventSaturationIdx]; } if (active_[Oil] && active_[Gas]) { @@ -952,8 +988,9 @@ namespace Opm { const int nw = wells().number_of_wells; const int numComp = numComponents(); - std::vector res(numComp*nw); + std::vector res(numEq*nw, 0.0); for( int compIdx = 0; compIdx < numComp; ++compIdx) { + for (int wellIdx = 0; wellIdx < nw; ++wellIdx) { int idx = wellIdx + nw*compIdx; res[idx] = resWell_[ wellIdx ][ compIdx ]; @@ -1006,7 +1043,7 @@ namespace Opm { B += 1 / fs.invB(ebosPhaseIdx).value(); } if (has_solvent_) { - auto& B = B_avg[ solventCompIdx ]; + auto& B = B_avg[ solventSaturationIdx ]; B += 1 / intQuants.solventInverseFormationVolumeFactor().value(); } } @@ -1204,8 +1241,8 @@ namespace Opm { // 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(); + b_perf[numComp*perf + solventSaturationIdx] = intQuants.solventInverseFormationVolumeFactor().value(); + surf_dens_perf[numComp*perf + solventSaturationIdx] = intQuants.solventRefDensity(); } } } @@ -1235,20 +1272,20 @@ namespace Opm { // update the second and third well variable (The flux fractions) std::vector F(np,0.0); if (active_[ Water ]) { - const int sign2 = dwells[w][WFrac] > 0 ? 1: -1; - const double dx2_limited = sign2 * std::min(std::abs(dwells[w][WFrac]),dFLimit); + 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 (active_[ Gas ]) { - const int sign3 = dwells[w][GFrac] > 0 ? 1: -1; - const double dx3_limited = sign3 * std::min(std::abs(dwells[w][GFrac]),dFLimit); + 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; } if (has_solvent_) { - const int sign4 = dwells[w][SFrac] > 0 ? 1: -1; - const double dx4_limited = sign4 * std::min(std::abs(dwells[w][SFrac]),dFLimit); + 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; } @@ -1352,7 +1389,7 @@ namespace Opm { case THP: // The BHP and THP both uses the total rate as first well variable. case BHP: { - well_state.wellSolutions()[nw*XvarWell + w] = xvar_well_old[nw*XvarWell + w] - dwells[w][XvarWell]; + well_state.wellSolutions()[nw*XvarWell + w] = xvar_well_old[nw*XvarWell + w] - dwells[w][flowPhaseToEbosCompIdx(XvarWell)]; switch (wells().type[w]) { case INJECTOR: @@ -1420,8 +1457,8 @@ namespace Opm { case SURFACE_RATE: // Both rate controls use bhp as first well variable case RESERVOIR_RATE: { - const int sign1 = dwells[w][XvarWell] > 0 ? 1: -1; - const double dx1_limited = sign1 * std::min(std::abs(dwells[w][XvarWell]),std::abs(xvar_well_old[nw*XvarWell + w])*dBHPLimit); + 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]; @@ -1764,7 +1801,7 @@ namespace Opm { perfRates[perf*numComponent + phase] = xw.perfPhaseRates()[perf*np + phase]; } if(has_solvent_) { - perfRates[perf*numComponent + solventCompIdx] = xw.perfRateSolvent()[perf]; + perfRates[perf*numComponent + solventSaturationIdx] = xw.perfRateSolvent()[perf]; } } well_perforation_densities_ = @@ -2143,7 +2180,7 @@ namespace Opm { if (wells().type[wellIdx] == INJECTOR) { if (has_solvent_ ) { double comp_frac = 0.0; - if (compIdx == solventCompIdx) { // solvent + if (has_solvent_ && compIdx == solventSaturationIdx) { // 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)); @@ -2209,7 +2246,7 @@ namespace Opm { 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); + wellVolumeFractionScaledPhaseUnderControl += wellVolumeFractionScaled(wellIdx, solventSaturationIdx); } if (compIdx == phase_under_control) { @@ -2274,7 +2311,7 @@ namespace Opm { return wellVariables_[GFrac * nw + wellIdx]; } - if (compIdx == solventCompIdx) { + if (has_solvent_ && compIdx == solventSaturationIdx) { return wellVariables_[SFrac * nw + wellIdx]; } @@ -2305,7 +2342,7 @@ namespace Opm { const WellControls* wc = wells().ctrls[wellIdx]; if (well_controls_get_current_type(wc) == RESERVOIR_RATE) { - if (has_solvent_ && compIdx == solventCompIdx) { + if (has_solvent_ && compIdx == solventSaturationIdx) { return wellVolumeFraction(wellIdx, compIdx); } const double* distr = well_controls_get_current_distr(wc); @@ -2710,10 +2747,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] * (1.0 - wsolvent(well_index)) * xw.wellRates()[np*well_index + Gas] / tot_well_rate ; + xw.wellSolutions()[GFrac*nw + well_index] = g[Gas] * (xw.wellRates()[np*well_index + Gas] - xw.solventWellRate(well_index)) / 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 ; + xw.wellSolutions()[SFrac*nw + well_index] = g[Gas] * xw.solventWellRate(well_index) / tot_well_rate ; } } else { const WellType& well_type = wells().type[well_index]; @@ -3027,6 +3064,171 @@ namespace Opm { } + template + double + StandardWellsDense:: + wpolymer(const int well_index) const { + + if (!has_polymer_) { + 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_); + WellPolymerProperties polymer = well->getPolymerProperties(current_timeIdx_); + if (injection.injectorType == WellInjector::WATER) { + + double polymerFraction = polymer.m_polymerConcentration; + + // Look until we find the correct well + if (well->name() == wells().name[well_index]) { + return polymerFraction; + } + } + } + // we didn't find it return 0; + assert(false); + return 0.0; + } + + + template + void + StandardWellsDense:: + setupCompressedToCartesian(const int* global_cell, int number_of_cells, std::map& cartesian_to_compressed ) const + { + if (global_cell) { + for (int i = 0; i < number_of_cells; ++i) { + cartesian_to_compressed.insert(std::make_pair(global_cell[i], i)); + } + } + else { + for (int i = 0; i < number_of_cells; ++i) { + cartesian_to_compressed.insert(std::make_pair(i, i)); + } + } + + } + + template + void + StandardWellsDense:: + computeRepRadiusPerfLength(const auto& grid) + { + + // TODO, the function does not work for parallel running + // to be fixed later. + int number_of_cells = Opm::UgGridHelpers::numCells(grid); + const int* global_cell = Opm::UgGridHelpers::globalCell(grid); + const int* cart_dims = Opm::UgGridHelpers::cartDims(grid); + auto cell_to_faces = Opm::UgGridHelpers::cell2Faces(grid); + auto begin_face_centroids = Opm::UgGridHelpers::beginFaceCentroids(grid); + + if (wells_ecl_.size() == 0) { + OPM_MESSAGE("No wells specified in Schedule section, " + "initializing no wells"); + return; + } + + const int nw = wells().number_of_wells; + const int nperf = wells().well_connpos[nw]; + + const size_t timeStep = current_timeIdx_; + + wells_rep_radius_.clear(); + wells_perf_length_.clear(); + wells_bore_diameter_.clear(); + + wells_rep_radius_.reserve(nperf); + wells_perf_length_.reserve(nperf); + wells_bore_diameter_.reserve(nperf); + + std::map cartesian_to_compressed; + + setupCompressedToCartesian(global_cell, number_of_cells, + cartesian_to_compressed); + + int well_index = 0; + + for (auto wellIter= wells_ecl_.begin(); wellIter != wells_ecl_.end(); ++wellIter) { + const auto* well = (*wellIter); + + if (well->getStatus(timeStep) == WellCommon::SHUT) { + continue; + } + { // COMPDAT handling + const auto& completionSet = well->getCompletions(timeStep); + for (size_t c=0; c::const_iterator cgit = cartesian_to_compressed.find(cart_grid_indx); + if (cgit == cartesian_to_compressed.end()) { + OPM_THROW(std::runtime_error, "Cell with i,j,k indices " << i << ' ' << j << ' ' + << k << " not found in grid (well = " << well->name() << ')'); + } + int cell = cgit->second; + + { + double radius = 0.5*completion.getDiameter(); + if (radius <= 0.0) { + radius = 0.5*unit::feet; + OPM_MESSAGE("**** Warning: Well bore internal radius set to " << radius); + } + + const std::array cubical = + WellsManagerDetail::getCubeDim<3>(cell_to_faces, begin_face_centroids, cell); + + WellCompletion::DirectionEnum direction = completion.getDirection(); + + double re; // area equivalent radius of the grid block + double perf_length; // the length of the well perforation + + switch (direction) { + case Opm::WellCompletion::DirectionEnum::X: + re = std::sqrt(cubical[1] * cubical[2] / M_PI); + perf_length = cubical[0]; + break; + case Opm::WellCompletion::DirectionEnum::Y: + re = std::sqrt(cubical[0] * cubical[2] / M_PI); + perf_length = cubical[1]; + break; + case Opm::WellCompletion::DirectionEnum::Z: + re = std::sqrt(cubical[0] * cubical[1] / M_PI); + perf_length = cubical[2]; + break; + default: + OPM_THROW(std::runtime_error, " Dirtecion of well is not supported "); + } + + double repR = std::sqrt(re * radius); + wells_rep_radius_.push_back(repR); + wells_perf_length_.push_back(perf_length); + wells_bore_diameter_.push_back(2. * radius); + } + } else { + if (completion.getState() != WellCompletion::SHUT) { + OPM_THROW(std::runtime_error, "Completion state: " << WellCompletion::StateEnum2String( completion.getState() ) << " not handled"); + } + } + + } + } + well_index++; + } + } + + } // namespace Opm From 17ada607eb3afedff949955161dd8b0423160fd1 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Thu, 22 Jun 2017 16:02:59 +0200 Subject: [PATCH 2/6] Assume polymer and water is fully mixed in the well --- opm/autodiff/StandardWellsDense_impl.hpp | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/opm/autodiff/StandardWellsDense_impl.hpp b/opm/autodiff/StandardWellsDense_impl.hpp index 854aeb868..bd66466bc 100644 --- a/opm/autodiff/StandardWellsDense_impl.hpp +++ b/opm/autodiff/StandardWellsDense_impl.hpp @@ -205,6 +205,16 @@ namespace Opm { std::vector cq_s(numComp,0.0); std::vector mob(numComp, 0.0); getMobility(ebosSimulator, perf, cell_idx, mob); + + if (has_polymer_) { + // assume fully mixture for wells. + EvalWell polymerConcentration = extendEval(intQuants.polymerConcentration()); + + if (wells().type[w] == INJECTOR) { + const auto& viscosityMultiplier = PolymerModule::plyviscViscosityMultiplierTable(intQuants.pvtRegionIndex()); + mob[ Water ] /= (extendEval(intQuants.waterViscosityCorrection()) * viscosityMultiplier.eval(polymerConcentration, /*extrapolate=*/true) ); + } + } computeWellFlux(w, wells().WI[perf], intQuants, mob, bhp, wellPerforationPressureDiffs()[perf], allow_cf, cq_s); if (has_polymer_) { @@ -220,7 +230,6 @@ namespace Opm { // guard against zero porosity and no water const EvalWell denom = Opm::max( (area * poro * (Sw - Swcr)), 1e-12); EvalWell waterVelocity = cq_s[ Water ] / denom * extendEval(intQuants.fluidState().invB(flowPhaseToEbosPhaseIdx(Water))); - EvalWell polymerConcentration = extendEval(intQuants.polymerConcentration()); if (PolymerModule::hasShrate()) { // TODO Use the same conversion as for the reservoar equations. @@ -228,6 +237,7 @@ namespace Opm { // For now use the same formula as in legacy. waterVelocity *= PolymerModule::shrate( intQuants.pvtRegionIndex() ) / wells_bore_diameter_[perf]; } + EvalWell polymerConcentration = extendEval(intQuants.polymerConcentration()); EvalWell shearFactor = PolymerModule::computeShearFactor(polymerConcentration, intQuants.pvtRegionIndex(), waterVelocity); From 8088347c96232c5e313f9bf187a1d188f79b9127 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Fri, 23 Jun 2017 08:22:30 +0200 Subject: [PATCH 3/6] Move adjustment of water mobility caused by polymer to getMobility() --- opm/autodiff/StandardWellsDense.hpp | 1 + opm/autodiff/StandardWellsDense_impl.hpp | 92 ++++++++++++------------ 2 files changed, 49 insertions(+), 44 deletions(-) diff --git a/opm/autodiff/StandardWellsDense.hpp b/opm/autodiff/StandardWellsDense.hpp index 2557e5cde..114175433 100644 --- a/opm/autodiff/StandardWellsDense.hpp +++ b/opm/autodiff/StandardWellsDense.hpp @@ -151,6 +151,7 @@ enum WellVariablePositions { void getMobility(const Simulator& ebosSimulator, + const int w, const int perf, const int cell_idx, std::vector& mob) const; diff --git a/opm/autodiff/StandardWellsDense_impl.hpp b/opm/autodiff/StandardWellsDense_impl.hpp index bd66466bc..8d79a8a90 100644 --- a/opm/autodiff/StandardWellsDense_impl.hpp +++ b/opm/autodiff/StandardWellsDense_impl.hpp @@ -204,50 +204,9 @@ namespace Opm { const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); std::vector cq_s(numComp,0.0); std::vector mob(numComp, 0.0); - getMobility(ebosSimulator, perf, cell_idx, mob); - - if (has_polymer_) { - // assume fully mixture for wells. - EvalWell polymerConcentration = extendEval(intQuants.polymerConcentration()); - - if (wells().type[w] == INJECTOR) { - const auto& viscosityMultiplier = PolymerModule::plyviscViscosityMultiplierTable(intQuants.pvtRegionIndex()); - mob[ Water ] /= (extendEval(intQuants.waterViscosityCorrection()) * viscosityMultiplier.eval(polymerConcentration, /*extrapolate=*/true) ); - } - } + getMobility(ebosSimulator, w, perf, cell_idx, mob); computeWellFlux(w, wells().WI[perf], intQuants, mob, bhp, wellPerforationPressureDiffs()[perf], allow_cf, cq_s); - if (has_polymer_) { - if (PolymerModule::hasPlyshlog()) { - // compute the well water velocity based on the perforation rates. - double area = 2 * M_PI * wells_rep_radius_[perf] * wells_perf_length_[perf]; - const auto& materialLawManager = ebosSimulator.problem().materialLawManager(); - const auto& scaledDrainageInfo = - materialLawManager->oilWaterScaledEpsInfoDrainage(cell_idx); - const Scalar& Swcr = scaledDrainageInfo.Swcr; - const EvalWell poro = extendEval(intQuants.porosity()); - const EvalWell Sw = extendEval(intQuants.fluidState().saturation(flowPhaseToEbosPhaseIdx(Water))); - // guard against zero porosity and no water - const EvalWell denom = Opm::max( (area * poro * (Sw - Swcr)), 1e-12); - EvalWell waterVelocity = cq_s[ Water ] / denom * extendEval(intQuants.fluidState().invB(flowPhaseToEbosPhaseIdx(Water))); - - if (PolymerModule::hasShrate()) { - // TODO Use the same conversion as for the reservoar equations. - // Need the "permeability" of the well? - // For now use the same formula as in legacy. - waterVelocity *= PolymerModule::shrate( intQuants.pvtRegionIndex() ) / wells_bore_diameter_[perf]; - } - EvalWell polymerConcentration = extendEval(intQuants.polymerConcentration()); - EvalWell shearFactor = PolymerModule::computeShearFactor(polymerConcentration, - intQuants.pvtRegionIndex(), - waterVelocity); - - // modify the mobility with the shear factor and recompute the well fluxes. - mob[ Water ] /= shearFactor; - computeWellFlux(w, wells().WI[perf], intQuants, mob, bhp, wellPerforationPressureDiffs()[perf], allow_cf, cq_s); - } - } - for (int componentIdx = 0; componentIdx < numComp; ++componentIdx) { // the cq_s entering mass balance equations need to consider the efficiency factors. @@ -309,7 +268,7 @@ namespace Opm { template void StandardWellsDense:: - getMobility(const Simulator& ebosSimulator, const int perf, const int cell_idx, std::vector& mob) const + getMobility(const Simulator& ebosSimulator, const int w, const int perf, const int cell_idx, std::vector& mob) const { const int np = wells().number_of_phases; @@ -350,6 +309,51 @@ namespace Opm { OPM_THROW(std::runtime_error, "individual mobility for wells does not work in combination with solvent"); } } + + // modify the water mobility if polymer is present + if (has_polymer_) { + // assume fully mixture for wells. + EvalWell polymerConcentration = extendEval(intQuants.polymerConcentration()); + + if (wells().type[w] == INJECTOR) { + const auto& viscosityMultiplier = PolymerModule::plyviscViscosityMultiplierTable(intQuants.pvtRegionIndex()); + mob[ Water ] /= (extendEval(intQuants.waterViscosityCorrection()) * viscosityMultiplier.eval(polymerConcentration, /*extrapolate=*/true) ); + } + + if (PolymerModule::hasPlyshlog()) { + // compute the well water velocity with out shear effects. + const int numComp = numComponents(); + bool allow_cf = allow_cross_flow(w, ebosSimulator); + const EvalWell& bhp = getBhp(w); + std::vector cq_s(numComp,0.0); + computeWellFlux(w, wells().WI[perf], intQuants, mob, bhp, wellPerforationPressureDiffs()[perf], allow_cf, cq_s); + double area = 2 * M_PI * wells_rep_radius_[perf] * wells_perf_length_[perf]; + const auto& materialLawManager = ebosSimulator.problem().materialLawManager(); + const auto& scaledDrainageInfo = + materialLawManager->oilWaterScaledEpsInfoDrainage(cell_idx); + const Scalar& Swcr = scaledDrainageInfo.Swcr; + const EvalWell poro = extendEval(intQuants.porosity()); + const EvalWell Sw = extendEval(intQuants.fluidState().saturation(flowPhaseToEbosPhaseIdx(Water))); + // guard against zero porosity and no water + const EvalWell denom = Opm::max( (area * poro * (Sw - Swcr)), 1e-12); + EvalWell waterVelocity = cq_s[ Water ] / denom * extendEval(intQuants.fluidState().invB(flowPhaseToEbosPhaseIdx(Water))); + + if (PolymerModule::hasShrate()) { + // TODO Use the same conversion as for the reservoar equations. + // Need the "permeability" of the well? + // For now use the same formula as in legacy. + waterVelocity *= PolymerModule::shrate( intQuants.pvtRegionIndex() ) / wells_bore_diameter_[perf]; + } + EvalWell polymerConcentration = extendEval(intQuants.polymerConcentration()); + EvalWell shearFactor = PolymerModule::computeShearFactor(polymerConcentration, + intQuants.pvtRegionIndex(), + waterVelocity); + + // modify the mobility with the shear factor and recompute the well fluxes. + mob[ Water ] /= shearFactor; + } + } + } @@ -2844,7 +2848,7 @@ namespace Opm { // flux for each perforation std::vector cq_s(numComp, 0.0); std::vector mob(numComp, 0.0); - getMobility(ebosSimulator, perf, cell_index, mob); + getMobility(ebosSimulator, well_index, perf, cell_index, mob); computeWellFlux(well_index, wells().WI[perf], intQuants, mob, bhp, wellPerforationPressureDiffs()[perf], allow_cf, cq_s); From e9a1aa2a832ffb02a4106a4fccd0404bbd790c2c Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Fri, 23 Jun 2017 10:51:34 +0200 Subject: [PATCH 4/6] Clean up after rebase --- opm/autodiff/BlackoilModelEbos.hpp | 12 ++++- opm/autodiff/StandardWellsDense_impl.hpp | 57 ++++++++++++++++++------ 2 files changed, 54 insertions(+), 15 deletions(-) diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 4d48e28bc..740e116fa 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -887,7 +887,7 @@ namespace Opm { maxCoeff[ contiSolventEqIdx ] = std::max( maxCoeff[ contiSolventEqIdx ], std::abs( R2 ) / pvValue ); } if (has_polymer_ ) { - B_avg[ contiPolymerEqIdx ] += 1.0 / intQuants.solventInverseFormationVolumeFactor().value(); + B_avg[ contiPolymerEqIdx ] += 1.0 / fs.invB(FluidSystem::waterPhaseIdx).value(); const auto R2 = ebosResid[cell_idx][contiPolymerEqIdx]; R_sum[ contiPolymerEqIdx ] += R2; maxCoeff[ contiPolymerEqIdx ] = std::max( maxCoeff[ contiPolymerEqIdx ], std::abs( R2 ) / pvValue ); @@ -1280,6 +1280,11 @@ namespace Opm { } VectorType& ssol = has_solvent_ ? simData.getCellData( "SSOL" ) : zero; + if (has_polymer_) { + simData.registerCellData( "POLYMER", 1 ); + } + VectorType& cpolymer = has_polymer_ ? simData.getCellData( "POLYMER" ) : zero; + std::vector failed_cells_pb; std::vector failed_cells_pd; const auto& gridView = ebosSimulator().gridView(); @@ -1368,6 +1373,11 @@ namespace Opm { ssol[cellIdx] = intQuants.solventSaturation().value(); } + if (has_polymer_) + { + cpolymer[cellIdx] = intQuants.polymerConcentration().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 diff --git a/opm/autodiff/StandardWellsDense_impl.hpp b/opm/autodiff/StandardWellsDense_impl.hpp index 8d79a8a90..e651902de 100644 --- a/opm/autodiff/StandardWellsDense_impl.hpp +++ b/opm/autodiff/StandardWellsDense_impl.hpp @@ -219,16 +219,23 @@ namespace Opm { } // subtract sum of phase fluxes in the well equations. - resWell_[w][flowPhaseToEbosCompIdx(componentIdx)] -= cq_s[componentIdx].value(); + resWell_[w][componentIdx] -= cq_s[componentIdx].value(); // assemble the jacobians for (int pvIdx = 0; pvIdx < numWellEq; ++pvIdx) { if (!only_wells) { // also need to consider the efficiency factor when manipulating the jacobians. - ebosJac[cell_idx][cell_idx][flowPhaseToEbosCompIdx(componentIdx)][flowToEbosPvIdx(pvIdx)] -= cq_s_effective.derivative(pvIdx); - duneB_[w][cell_idx][flowToEbosPvIdx(pvIdx)][flowPhaseToEbosCompIdx(componentIdx)] -= cq_s_effective.derivative(pvIdx+numEq); // intput in transformed matrix + duneB_[w][cell_idx][pvIdx][flowPhaseToEbosCompIdx(componentIdx)] -= cq_s_effective.derivative(pvIdx+numEq); // intput in transformed matrix + } + invDuneD_[w][w][componentIdx][pvIdx] -= cq_s[componentIdx].derivative(pvIdx+numEq); + } + + for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) { + if (!only_wells) { + // also need to consider the efficiency factor when manipulating the jacobians. + ebosJac[cell_idx][cell_idx][flowPhaseToEbosCompIdx(componentIdx)][flowToEbosPvIdx(pvIdx)] -= cq_s_effective.derivative(pvIdx); + duneC_[w][cell_idx][componentIdx][flowToEbosPvIdx(pvIdx)] -= cq_s_effective.derivative(pvIdx); } - invDuneD_[w][w][flowPhaseToEbosCompIdx(componentIdx)][pvIdx] -= cq_s[componentIdx].derivative(pvIdx+numEq); } // add trivial equation for 2p cases (Only support water + oil) @@ -245,6 +252,21 @@ namespace Opm { } } + if (has_polymer_) { + EvalWell cq_s_poly = cq_s[Water]; + if (wells().type[w] == INJECTOR) { + cq_s_poly *= wpolymer(w); + } else { + cq_s_poly *= extendEval(intQuants.polymerConcentration()); + } + if (!only_wells) { + for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) { + ebosJac[cell_idx][cell_idx][contiPolymerEqIdx][flowToEbosPvIdx(pvIdx)] -= cq_s_poly.derivative(pvIdx); + } + ebosResid[cell_idx][contiPolymerEqIdx] -= cq_s_poly.value(); + } + } + // Store the perforation pressure for later usage. well_state.perfPress()[perf] = well_state.bhp()[w] + wellPerforationPressureDiffs()[perf]; } @@ -258,8 +280,15 @@ namespace Opm { } resWell_[w][componentIdx] += resWell_loc.value(); } + + // add trivial equation for polymer + if (has_polymer_) { + invDuneD_[w][w][contiPolymerEqIdx][polymerConcentrationIdx] = 1.0; // + } } + + // do the local inversion of D. localInvert( invDuneD_ ); } @@ -287,7 +316,7 @@ namespace Opm { mob[phase] = extendEval(intQuants.mobility(ebosPhaseIdx)); } if (has_solvent_) { - mob[solventCompIdx] = extendEval(intQuants.solventMobility()); + mob[solventSaturationIdx] = extendEval(intQuants.solventMobility()); } } else { @@ -1286,20 +1315,20 @@ namespace Opm { // update the second and third well variable (The flux fractions) std::vector F(np,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); + const int sign2 = dwells[w][WFrac] > 0 ? 1: -1; + const double dx2_limited = sign2 * std::min(std::abs(dwells[w][WFrac]),dFLimit); well_state.wellSolutions()[WFrac*nw + w] = xvar_well_old[WFrac*nw + w] - dx2_limited; } 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); + const int sign3 = dwells[w][GFrac] > 0 ? 1: -1; + const double dx3_limited = sign3 * std::min(std::abs(dwells[w][GFrac]),dFLimit); 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); + const int sign4 = dwells[w][SFrac] > 0 ? 1: -1; + const double dx4_limited = sign4 * std::min(std::abs(dwells[w][SFrac]),dFLimit); well_state.wellSolutions()[SFrac*nw + w] = xvar_well_old[SFrac*nw + w] - dx4_limited; } @@ -1403,7 +1432,7 @@ namespace Opm { case THP: // The BHP and THP both uses the total rate as first well variable. case BHP: { - well_state.wellSolutions()[nw*XvarWell + w] = xvar_well_old[nw*XvarWell + w] - dwells[w][flowPhaseToEbosCompIdx(XvarWell)]; + well_state.wellSolutions()[nw*XvarWell + w] = xvar_well_old[nw*XvarWell + w] - dwells[w][XvarWell]; switch (wells().type[w]) { case INJECTOR: @@ -1471,8 +1500,8 @@ namespace Opm { case SURFACE_RATE: // Both rate controls use bhp as first well variable case RESERVOIR_RATE: { - 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); + const int sign1 = dwells[w][XvarWell] > 0 ? 1: -1; + const double dx1_limited = sign1 * std::min(std::abs(dwells[w][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]; From 1bc2550541590fb73bb80aa80731b7dcac5363ff Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Fri, 23 Jun 2017 11:13:00 +0200 Subject: [PATCH 5/6] Fix polymer in well model for producers. --- opm/autodiff/BlackoilModelEbos.hpp | 1 - opm/autodiff/StandardWellsDense_impl.hpp | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 740e116fa..6c84f8d70 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -835,7 +835,6 @@ namespace Opm { const double tol_cnv = param_.tolerance_cnv_; const double tol_wells = param_.tolerance_wells_; - const int nc = Opm::AutoDiffGrid::numCells(grid_); const int np = numPhases(); const int numComp = numComponents(); diff --git a/opm/autodiff/StandardWellsDense_impl.hpp b/opm/autodiff/StandardWellsDense_impl.hpp index e651902de..c7f50afe8 100644 --- a/opm/autodiff/StandardWellsDense_impl.hpp +++ b/opm/autodiff/StandardWellsDense_impl.hpp @@ -257,7 +257,7 @@ namespace Opm { if (wells().type[w] == INJECTOR) { cq_s_poly *= wpolymer(w); } else { - cq_s_poly *= extendEval(intQuants.polymerConcentration()); + cq_s_poly *= extendEval(intQuants.polymerConcentration() * intQuants.polymerViscosityCorrection()); } if (!only_wells) { for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) { From 1c054022095b1e41079b9f4c8ae43ec8b7b1c47b Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Mon, 26 Jun 2017 08:40:30 +0200 Subject: [PATCH 6/6] Fix 2p case after rebase --- opm/autodiff/StandardWellsDense_impl.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/opm/autodiff/StandardWellsDense_impl.hpp b/opm/autodiff/StandardWellsDense_impl.hpp index c7f50afe8..87a954c4e 100644 --- a/opm/autodiff/StandardWellsDense_impl.hpp +++ b/opm/autodiff/StandardWellsDense_impl.hpp @@ -765,9 +765,9 @@ namespace Opm { setWellVariables(const WellState& xw) { const int nw = wells().number_of_wells; - // for two-phase numComp < numEq - //const int numComp = numComponents(); - for (int eqIdx = 0; eqIdx < numWellEq; ++eqIdx) { + // for two-phase numComp < numWellEq + const int numComp = numComponents(); + for (int eqIdx = 0; eqIdx < numComp; ++eqIdx) { for (int w = 0; w < nw; ++w) { const unsigned int idx = nw * eqIdx + w; assert( idx < wellVariables_.size() );