From 352dccd5e957b987fb71f0e9ea93eacdfd281f9c Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Thu, 31 Aug 2017 07:50:28 +0200 Subject: [PATCH 1/5] Make 2p gas oil compile and runs -- use mapping from canonicalToActiveCompIdx from Ebos -- add guards againts non-existing components --- examples/flow.cpp | 29 ++++- examples/flow_ebos_2p.cpp | 2 +- opm/autodiff/BlackoilModelEbos.hpp | 61 ++++------ .../SimulatorFullyImplicitBlackoilEbos.hpp | 8 +- opm/autodiff/StandardWell.hpp | 17 +-- opm/autodiff/StandardWell_impl.hpp | 104 +++++++++++------- opm/autodiff/StandardWellsDense_impl.hpp | 12 +- opm/autodiff/WellInterface_impl.hpp | 51 ++++----- 8 files changed, 160 insertions(+), 124 deletions(-) diff --git a/examples/flow.cpp b/examples/flow.cpp index 2227410d9..f0f28a3a9 100644 --- a/examples/flow.cpp +++ b/examples/flow.cpp @@ -52,12 +52,17 @@ namespace Properties { // Twophase case /////////////////////////////////// - NEW_TYPE_TAG(EclFlowTwoPhaseProblem, INHERITS_FROM(EclFlowProblem)); + NEW_TYPE_TAG(EclFlowOilWaterProblem, INHERITS_FROM(EclFlowProblem)); //! The indices required by the model - SET_TYPE_PROP(EclFlowTwoPhaseProblem, Indices, - Ewoms::BlackOilTwoPhaseIndices); + SET_TYPE_PROP(EclFlowOilWaterProblem, Indices, + Ewoms::BlackOilTwoPhaseIndices); + NEW_TYPE_TAG(EclFlowGasOilProblem, INHERITS_FROM(EclFlowProblem)); + //! The indices required by the model + SET_TYPE_PROP(EclFlowGasOilProblem, Indices, + Ewoms::BlackOilTwoPhaseIndices); + /////////////////////////////////// // Polymer case /////////////////////////////////// @@ -168,8 +173,22 @@ int main(int argc, char** argv) // Twophase case if( phases.size() == 2 ) { - Opm::FlowMainEbos mainfunc; - return mainfunc.execute(argc, argv, deck, eclipseState ); + // oil-gas + if (phases.active( Opm::Phase::GAS )) + { + Opm::FlowMainEbos mainfunc; + return mainfunc.execute(argc, argv, deck, eclipseState ); + } + // oil-water + else if ( phases.active( Opm::Phase::WATER ) ) + { + Opm::FlowMainEbos mainfunc; + return mainfunc.execute(argc, argv, deck, eclipseState ); + } + else { + std::cerr << "No suitable configuration found, valid are Twophase (oilwater and oilgas), polymer, solvent, or blackoil" << std::endl; + return EXIT_FAILURE; + } } // Polymer case else if ( phases.active( Opm::Phase::POLYMER ) ) { diff --git a/examples/flow_ebos_2p.cpp b/examples/flow_ebos_2p.cpp index 0ef174ca0..16431b067 100644 --- a/examples/flow_ebos_2p.cpp +++ b/examples/flow_ebos_2p.cpp @@ -38,7 +38,7 @@ namespace Properties { NEW_TYPE_TAG(EclFlowTwoPhaseProblem, INHERITS_FROM(EclFlowProblem)); //! The indices required by the model SET_TYPE_PROP(EclFlowTwoPhaseProblem, Indices, - Ewoms::BlackOilTwoPhaseIndices); + Ewoms::BlackOilTwoPhaseIndices); }} // ----------------- Main program ----------------- diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 82797ea17..8a4061873 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -635,7 +635,7 @@ namespace Opm { PrimaryVariables& priVars = solution[ cell_idx ]; - const double& dp = dx[cell_idx][flowPhaseToEbosCompIdx(0)]; + const double& dp = dx[cell_idx][Indices::pressureSwitchIdx]; double& p = priVars[Indices::pressureSwitchIdx]; const double& dp_rel_max = dpMaxRel(); const int sign_dp = dp > 0 ? 1: -1; @@ -643,9 +643,9 @@ namespace Opm { p = std::max(p, 0.0); // Saturation updates. - const double dsw = active_[Water] ? dx[cell_idx][flowPhaseToEbosCompIdx(1)] : 0.0; + const double dsw = active_[Water] ? dx[cell_idx][Indices::waterSaturationIdx] : 0.0; const int xvar_ind = active_[Water] ? 2 : 1; - const double dxvar = active_[Gas] ? dx[cell_idx][flowPhaseToEbosCompIdx(xvar_ind)] : 0.0; + const double dxvar = active_[Gas] ? dx[cell_idx][Indices::compositionSwitchIdx] : 0.0; double dso = 0.0; double dsg = 0.0; @@ -1533,51 +1533,36 @@ namespace Opm { public: - int ebosCompToFlowPhaseIdx( const int compIdx ) const - { - assert(compIdx < 3); - const int compToPhase[ 3 ] = { Oil, Water, Gas }; - return compToPhase[ compIdx ]; - } - - int flowToEbosPvIdx( const int flowPv ) const - { - const int flowToEbos[] = { - Indices::pressureSwitchIdx, - Indices::waterSaturationIdx, - Indices::compositionSwitchIdx, - Indices::solventSaturationIdx - }; - - if (flowPv > 2 ) - return flowPv; - - return flowToEbos[ flowPv ]; - } - int flowPhaseToEbosCompIdx( const int phaseIdx ) const { - const int phaseToComp[] = { - FluidSystem::waterCompIdx, - FluidSystem::oilCompIdx, - FluidSystem::gasCompIdx - }; + const auto& pu = phaseUsage_; + if (active_[Water] && pu.phase_pos[Water] == phaseIdx) + return Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); + if (active_[Oil] && pu.phase_pos[Oil] == phaseIdx) + return Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); + if (active_[Gas] && pu.phase_pos[Gas] == phaseIdx) + return Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); - if (phaseIdx > 2 ) - return phaseIdx; - - return phaseToComp[ phaseIdx ]; + // for other phases return the index + return phaseIdx; } - private: - int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const { + const auto& pu = phaseUsage_; + if (active_[Water] && pu.phase_pos[Water] == phaseIdx) + return FluidSystem::waterPhaseIdx; + if (active_[Oil] && pu.phase_pos[Oil] == phaseIdx) + return FluidSystem::oilPhaseIdx; + if (active_[Gas] && pu.phase_pos[Gas] == phaseIdx) + return FluidSystem::gasPhaseIdx; + assert(phaseIdx < 3); - const int flowToEbos[ 3 ] = { FluidSystem::waterPhaseIdx, FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx}; - return flowToEbos[ phaseIdx ]; + // for other phases return the index + return phaseIdx; } + private: void updateRateConverter() { diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp index 5085427fa..e7200b2b6 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp @@ -900,7 +900,9 @@ protected: // set non-switching primary variables PrimaryVariables& cellPv = solution[ cellIdx ]; // set water saturation - cellPv[BlackoilIndices::waterSaturationIdx] = saturations[cellIdx*numPhases + pu.phase_pos[Water]]; + if ( active[Water] ) { + cellPv[BlackoilIndices::waterSaturationIdx] = saturations[cellIdx*numPhases + pu.phase_pos[Water]]; + } if (has_solvent) { cellPv[BlackoilIndices::solventSaturationIdx] = reservoirState.getCellData( reservoirState.SSOL )[cellIdx]; @@ -944,7 +946,9 @@ protected: /*storeViscosity=*/false, /*storeEnthalpy=*/false> SatOnlyFluidState; SatOnlyFluidState fluidState; - fluidState.setSaturation(FluidSystem::waterPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Water]]); + if ( active[Water] ) { + fluidState.setSaturation(FluidSystem::waterPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Water]]); + } fluidState.setSaturation(FluidSystem::oilPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Oil]]); fluidState.setSaturation(FluidSystem::gasPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Gas]]); diff --git a/opm/autodiff/StandardWell.hpp b/opm/autodiff/StandardWell.hpp index 036ba51c1..7c72823e5 100644 --- a/opm/autodiff/StandardWell.hpp +++ b/opm/autodiff/StandardWell.hpp @@ -47,20 +47,21 @@ namespace Opm using typename Base::BlackoilIndices; using typename Base::PolymerModule; + using Base::numEq; + // the positions of the primary variables for StandardWell // there are three primary variables, the second and the third ones are F_w and F_g // the first one can be total rate (G_t) or bhp, based on the control - enum WellVariablePositions { - XvarWell = 0, - WFrac = 1, - GFrac = 2, - SFrac = 3 - }; + + static const bool gasoil = numEq == 2 && (BlackoilIndices::compositionSwitchIdx >= 0); + static const int XvarWell = 0; + static const int WFrac = gasoil? -1000: 1; + static const int GFrac = gasoil? 1: 2; + static const int SFrac = 3; using typename Base::Scalar; using typename Base::ConvergenceReport; - using Base::numEq; using Base::has_solvent; using Base::has_polymer; @@ -296,6 +297,8 @@ namespace Opm void getMobility(const Simulator& ebosSimulator, const int perf, std::vector& mob) const; + + double scalingFactor(const int comp_idx) const; }; } diff --git a/opm/autodiff/StandardWell_impl.hpp b/opm/autodiff/StandardWell_impl.hpp index 53c7e24cc..88f816bf1 100644 --- a/opm/autodiff/StandardWell_impl.hpp +++ b/opm/autodiff/StandardWell_impl.hpp @@ -318,8 +318,7 @@ namespace Opm } } - std::vector g = {1, 1, 0.01, 0.01}; - return (wellVolumeFraction(compIdx) / g[compIdx]); + return (wellVolumeFraction(compIdx) / scalingFactor(compIdx) ); } @@ -331,11 +330,12 @@ namespace Opm StandardWell:: wellVolumeFraction(const int compIdx) const { - if (compIdx == Water) { + const auto pu = phaseUsage(); + if (active()[Water] && compIdx == pu.phase_pos[Water]) { return primary_variables_evaluation_[WFrac]; } - if (compIdx == Gas) { + if (active()[Gas] && compIdx == pu.phase_pos[Gas]) { return primary_variables_evaluation_[GFrac]; } @@ -390,7 +390,7 @@ namespace Opm EvalWell out = 0.0; out.setValue(in.value()); for(int eqIdx = 0; eqIdx < numEq;++eqIdx) { - out.setDerivative(eqIdx, in.derivative(flowToEbosPvIdx(eqIdx))); + out.setDerivative(eqIdx, in.derivative(eqIdx)); } return out; } @@ -585,8 +585,8 @@ namespace Opm 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); - duneB_[0][cell_idx][componentIdx][flowToEbosPvIdx(pvIdx)] -= cq_s_effective.derivative(pvIdx); + ebosJac[cell_idx][cell_idx][flowPhaseToEbosCompIdx(componentIdx)][pvIdx] -= cq_s_effective.derivative(pvIdx); + duneB_[0][cell_idx][componentIdx][pvIdx] -= cq_s_effective.derivative(pvIdx); } } @@ -613,7 +613,7 @@ namespace Opm } if (!only_wells) { for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) { - ebosJac[cell_idx][cell_idx][contiPolymerEqIdx][flowToEbosPvIdx(pvIdx)] -= cq_s_poly.derivative(pvIdx); + ebosJac[cell_idx][cell_idx][contiPolymerEqIdx][pvIdx] -= cq_s_poly.derivative(pvIdx); } ebosResid[cell_idx][contiPolymerEqIdx] -= cq_s_poly.value(); } @@ -787,6 +787,7 @@ namespace Opm const int np = number_of_phases_; const double dBHPLimit = param.dbhp_max_rel_; const double dFLimit = param.dwell_fraction_max_; + const auto pu = phaseUsage(); const std::vector xvar_well_old = primary_variables_; @@ -811,68 +812,68 @@ namespace Opm } assert(active()[ Oil ]); - F[Oil] = 1.0; + F[pu.phase_pos[Oil]] = 1.0; if (active()[ Water ]) { - F[Water] = primary_variables_[WFrac]; - F[Oil] -= F[Water]; + F[pu.phase_pos[Water]] = primary_variables_[WFrac]; + F[pu.phase_pos[Oil]] -= F[pu.phase_pos[Water]]; } if (active()[ Gas ]) { - F[Gas] = primary_variables_[GFrac]; - F[Oil] -= F[Gas]; + F[pu.phase_pos[Gas]] = primary_variables_[GFrac]; + F[pu.phase_pos[Oil]] -= F[pu.phase_pos[Gas]]; } double F_solvent = 0.0; if (has_solvent) { F_solvent = primary_variables_[SFrac]; - F[Oil] -= F_solvent; + F[pu.phase_pos[Oil]] -= F_solvent; } if (active()[ Water ]) { if (F[Water] < 0.0) { if (active()[ Gas ]) { - F[Gas] /= (1.0 - F[Water]); + F[pu.phase_pos[Gas]] /= (1.0 - F[pu.phase_pos[Water]]); } if (has_solvent) { - F_solvent /= (1.0 - F[Water]); + F_solvent /= (1.0 - F[pu.phase_pos[Water]]); } - F[Oil] /= (1.0 - F[Water]); - F[Water] = 0.0; + F[pu.phase_pos[Oil]] /= (1.0 - F[pu.phase_pos[Water]]); + F[pu.phase_pos[Water]] = 0.0; } } if (active()[ Gas ]) { - if (F[Gas] < 0.0) { + if (F[pu.phase_pos[Gas]] < 0.0) { if (active()[ Water ]) { - F[Water] /= (1.0 - F[Gas]); + F[pu.phase_pos[Water]] /= (1.0 - F[pu.phase_pos[Gas]]); } if (has_solvent) { - F_solvent /= (1.0 - F[Gas]); + F_solvent /= (1.0 - F[pu.phase_pos[Gas]]); } - F[Oil] /= (1.0 - F[Gas]); - F[Gas] = 0.0; + F[pu.phase_pos[Oil]] /= (1.0 - F[pu.phase_pos[Gas]]); + F[pu.phase_pos[Gas]] = 0.0; } } - if (F[Oil] < 0.0) { + if (F[pu.phase_pos[Oil]] < 0.0) { if (active()[ Water ]) { - F[Water] /= (1.0 - F[Oil]); + F[pu.phase_pos[Water]] /= (1.0 - F[pu.phase_pos[Oil]]); } if (active()[ Gas ]) { - F[Gas] /= (1.0 - F[Oil]); + F[pu.phase_pos[Gas]] /= (1.0 - F[pu.phase_pos[Oil]]); } if (has_solvent) { - F_solvent /= (1.0 - F[Oil]); + F_solvent /= (1.0 - F[pu.phase_pos[Oil]]); } - F[Oil] = 0.0; + F[pu.phase_pos[Oil]] = 0.0; } if (active()[ Water ]) { - primary_variables_[WFrac] = F[Water]; + primary_variables_[WFrac] = F[pu.phase_pos[Water]]; } if (active()[ Gas ]) { - primary_variables_[GFrac] = F[Gas]; + primary_variables_[GFrac] = F[pu.phase_pos[Gas]]; } if(has_solvent) { primary_variables_[SFrac] = F_solvent; @@ -881,7 +882,7 @@ namespace Opm // 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; + F[pu.phase_pos[Gas]] += F_solvent; } // The interpretation of the first well variable depends on the well control @@ -892,7 +893,6 @@ namespace Opm const int current = well_state.currentControls()[index_of_well_]; const double target_rate = well_controls_iget_target(wc, current); - std::vector g = {1,1,0.01}; if (well_controls_iget_type(wc, current) == RESERVOIR_RATE) { const double* distr = well_controls_iget_distr(wc, current); for (int p = 0; p < np; ++p) { @@ -904,7 +904,7 @@ namespace Opm } } else { for (int p = 0; p < np; ++p) { - F[p] /= g[p]; + F[p] /= scalingFactor(p); } } @@ -1225,7 +1225,7 @@ namespace Opm const int w = index_of_well_; //rs and rv are only used if both oil and gas is present - if (pu.phase_used[BlackoilPhases::Vapour] && pu.phase_pos[BlackoilPhases::Liquid]) { + if (pu.phase_used[BlackoilPhases::Vapour] && pu.phase_used[BlackoilPhases::Liquid]) { rsmax_perf.resize(nperf); rvmax_perf.resize(nperf); } @@ -1874,12 +1874,17 @@ namespace Opm const int well_index = index_of_well_; const WellControls* wc = well_controls_; const double* distr = well_controls_get_current_distr(wc); + const auto pu = phaseUsage(); - std::vector g = {1.0, 1.0, 0.01}; + std::vector g(np); if (well_controls_get_current_type(wc) == RESERVOIR_RATE) { for (int phase = 0; phase < np; ++phase) { g[phase] = distr[phase]; } + } else { + for (int phase = 0; phase < np; ++phase) { + g[phase] = scalingFactor(phase); + } } switch (well_controls_get_current_type(wc)) { @@ -1909,13 +1914,13 @@ namespace Opm } if(std::abs(tot_well_rate) > 0) { if (active()[ Water ]) { - primary_variables_[WFrac] = g[Water] * well_state.wellRates()[np*well_index + Water] / tot_well_rate; + primary_variables_[WFrac] = g[pu.phase_pos[Water]] * well_state.wellRates()[np*well_index + pu.phase_pos[Water]] / tot_well_rate; } if (active()[ Gas ]) { - primary_variables_[GFrac] = g[Gas] * (well_state.wellRates()[np*well_index + Gas] - well_state.solventWellRate(well_index)) / tot_well_rate ; + primary_variables_[GFrac] = g[ pu.phase_pos[Gas]] * (well_state.wellRates()[np*well_index + pu.phase_pos[Gas]] - well_state.solventWellRate(well_index)) / tot_well_rate ; } if (has_solvent) { - primary_variables_[SFrac] = g[Gas] * well_state.solventWellRate(well_index) / tot_well_rate ; + primary_variables_[SFrac] = g[ pu.phase_pos[Gas]] * well_state.solventWellRate(well_index) / tot_well_rate ; } } else { // tot_well_rate == 0 if (well_type_ == INJECTOR) { @@ -1929,7 +1934,7 @@ namespace Opm } if (active()[Gas]) { - if (distr[Gas] > 0.0) { + if (distr[pu.phase_pos[Gas]] > 0.0) { primary_variables_[GFrac] = 1.0 - wsolvent(); if (has_solvent) { primary_variables_[SFrac] = wsolvent(); @@ -2053,4 +2058,25 @@ namespace Opm return thp; } + template + double + StandardWell::scalingFactor(const int phaseIdx) const + { + //std::vector g = {1,1,0.01,0.01}; + const auto& pu = phaseUsage(); + if (active()[Water] && pu.phase_pos[Water] == phaseIdx) + return 1.0; + if (active()[Oil] && pu.phase_pos[Oil] == phaseIdx) + return 1.0; + if (active()[Gas] && pu.phase_pos[Gas] == phaseIdx) + return 0.01; + if (has_solvent && phaseIdx == contiSolventEqIdx ) + return 0.01; + + // we should come this fare + assert(false); + return 1.0; + } + + } diff --git a/opm/autodiff/StandardWellsDense_impl.hpp b/opm/autodiff/StandardWellsDense_impl.hpp index 84d47d00d..2514c2bce 100644 --- a/opm/autodiff/StandardWellsDense_impl.hpp +++ b/opm/autodiff/StandardWellsDense_impl.hpp @@ -302,9 +302,17 @@ namespace Opm { StandardWellsDense:: flowPhaseToEbosPhaseIdx( const int phaseIdx ) const { + const auto& pu = phase_usage_; + if (active_[Water] && pu.phase_pos[Water] == phaseIdx) + return FluidSystem::waterPhaseIdx; + if (active_[Oil] && pu.phase_pos[Oil] == phaseIdx) + return FluidSystem::oilPhaseIdx; + if (active_[Gas] && pu.phase_pos[Gas] == phaseIdx) + return FluidSystem::gasPhaseIdx; + assert(phaseIdx < 3); - const int flowToEbos[ 3 ] = { FluidSystem::waterPhaseIdx, FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx }; - return flowToEbos[ phaseIdx ]; + // for other phases return the index + return phaseIdx; } diff --git a/opm/autodiff/WellInterface_impl.hpp b/opm/autodiff/WellInterface_impl.hpp index 6133291e0..e2c958f1e 100644 --- a/opm/autodiff/WellInterface_impl.hpp +++ b/opm/autodiff/WellInterface_impl.hpp @@ -224,47 +224,38 @@ namespace Opm WellInterface:: flowPhaseToEbosCompIdx( const int phaseIdx ) const { - const int phaseToComp[ 3 ] = { FluidSystem::waterCompIdx, FluidSystem::oilCompIdx, FluidSystem::gasCompIdx}; - if (phaseIdx > 2 ) - return phaseIdx; - return phaseToComp[ phaseIdx ]; + const auto& pu = phaseUsage(); + if (active()[Water] && pu.phase_pos[Water] == phaseIdx) + return BlackoilIndices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); + if (active()[Oil] && pu.phase_pos[Oil] == phaseIdx) + return BlackoilIndices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); + if (active()[Gas] && pu.phase_pos[Gas] == phaseIdx) + return BlackoilIndices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + + // for other phases return the index + return phaseIdx; } - - template - int - WellInterface:: - 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 WellInterface:: flowPhaseToEbosPhaseIdx( const int phaseIdx ) const { - assert(phaseIdx < 3); - const int flowToEbos[ 3 ] = { FluidSystem::waterPhaseIdx, FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx }; - return flowToEbos[ phaseIdx ]; - } + const auto& pu = phaseUsage(); + if (active()[Water] && pu.phase_pos[Water] == phaseIdx) + return FluidSystem::waterPhaseIdx; + if (active()[Oil] && pu.phase_pos[Oil] == phaseIdx) + return FluidSystem::oilPhaseIdx; + if (active()[Gas] && pu.phase_pos[Gas] == phaseIdx) + return FluidSystem::gasPhaseIdx; + assert(phaseIdx < 3); + // for other phases return the index + return phaseIdx; + } From 69c608829f9cee6c2df205c7df9d794d92a829ec Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Tue, 5 Sep 2017 09:25:34 +0200 Subject: [PATCH 2/5] Include RESV in the scaling factor - solvent + RESV is not correct. Make the simulator throw until this is sorted out. - remove unused parameter --- examples/flow.cpp | 4 +- examples/flow_ebos_2p.cpp | 2 +- opm/autodiff/BlackoilModelEbos.hpp | 1 - opm/autodiff/StandardWell_impl.hpp | 72 +++++++++++------------------- 4 files changed, 28 insertions(+), 51 deletions(-) diff --git a/examples/flow.cpp b/examples/flow.cpp index f0f28a3a9..e5885177f 100644 --- a/examples/flow.cpp +++ b/examples/flow.cpp @@ -55,13 +55,13 @@ namespace Properties { NEW_TYPE_TAG(EclFlowOilWaterProblem, INHERITS_FROM(EclFlowProblem)); //! The indices required by the model SET_TYPE_PROP(EclFlowOilWaterProblem, Indices, - Ewoms::BlackOilTwoPhaseIndices); + Ewoms::BlackOilTwoPhaseIndices); NEW_TYPE_TAG(EclFlowGasOilProblem, INHERITS_FROM(EclFlowProblem)); //! The indices required by the model SET_TYPE_PROP(EclFlowGasOilProblem, Indices, - Ewoms::BlackOilTwoPhaseIndices); + Ewoms::BlackOilTwoPhaseIndices); /////////////////////////////////// // Polymer case diff --git a/examples/flow_ebos_2p.cpp b/examples/flow_ebos_2p.cpp index 16431b067..5e802331e 100644 --- a/examples/flow_ebos_2p.cpp +++ b/examples/flow_ebos_2p.cpp @@ -38,7 +38,7 @@ namespace Properties { NEW_TYPE_TAG(EclFlowTwoPhaseProblem, INHERITS_FROM(EclFlowProblem)); //! The indices required by the model SET_TYPE_PROP(EclFlowTwoPhaseProblem, Indices, - Ewoms::BlackOilTwoPhaseIndices); + Ewoms::BlackOilTwoPhaseIndices); }} // ----------------- Main program ----------------- diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 8a4061873..b3c5a862e 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -644,7 +644,6 @@ namespace Opm { // Saturation updates. const double dsw = active_[Water] ? dx[cell_idx][Indices::waterSaturationIdx] : 0.0; - const int xvar_ind = active_[Water] ? 2 : 1; const double dxvar = active_[Gas] ? dx[cell_idx][Indices::compositionSwitchIdx] : 0.0; double dso = 0.0; diff --git a/opm/autodiff/StandardWell_impl.hpp b/opm/autodiff/StandardWell_impl.hpp index 88f816bf1..529e2ba82 100644 --- a/opm/autodiff/StandardWell_impl.hpp +++ b/opm/autodiff/StandardWell_impl.hpp @@ -299,26 +299,13 @@ namespace Opm StandardWell:: wellVolumeFractionScaled(const int compIdx) const { - // TODO: we should be able to set the g for the well based on the control type - // instead of using explicit code for g all the times - const WellControls* wc = well_controls_; - if (well_controls_get_current_type(wc) == RESERVOIR_RATE) { - if (has_solvent && compIdx == contiSolventEqIdx) { - return wellVolumeFraction(compIdx); - } - const double* distr = well_controls_get_current_distr(wc); - assert(compIdx < 3); - if (distr[compIdx] > 0.) { - return wellVolumeFraction(compIdx) / distr[compIdx]; - } else { - // TODO: not sure why return EvalWell(0.) causing problem here - // Probably due to the wrong Jacobians. - return wellVolumeFraction(compIdx); - } - } + const double scal = scalingFactor(compIdx); + if (scal > 0) + return wellVolumeFraction(compIdx) / scal; - return (wellVolumeFraction(compIdx) / scalingFactor(compIdx) ); + // the scaling factor may be zero for RESV controlled wells. + return wellVolumeFraction(compIdx); } @@ -893,18 +880,12 @@ namespace Opm const int current = well_state.currentControls()[index_of_well_]; const double target_rate = well_controls_iget_target(wc, current); - if (well_controls_iget_type(wc, current) == RESERVOIR_RATE) { - const double* distr = well_controls_iget_distr(wc, current); - for (int p = 0; p < np; ++p) { - if (distr[p] > 0.) { // For injection wells, there only one non-zero distr value - F[p] /= distr[p]; - } else { - F[p] = 0.; - } - } - } else { - for (int p = 0; p < np; ++p) { - F[p] /= scalingFactor(p); + for (int p = 0; p < np; ++p) { + const double scal = scalingFactor(p); + if (scal > 0) { + F[p] /= scal ; + } else { + F[p] = 0.; } } @@ -1876,17 +1857,6 @@ namespace Opm const double* distr = well_controls_get_current_distr(wc); const auto pu = phaseUsage(); - std::vector g(np); - if (well_controls_get_current_type(wc) == RESERVOIR_RATE) { - for (int phase = 0; phase < np; ++phase) { - g[phase] = distr[phase]; - } - } else { - for (int phase = 0; phase < np; ++phase) { - g[phase] = scalingFactor(phase); - } - } - switch (well_controls_get_current_type(wc)) { case THP: case BHP: { @@ -1897,7 +1867,7 @@ namespace Opm } } else { for (int p = 0; p < np; ++p) { - primary_variables_[XvarWell] += g[p] * well_state.wellRates()[np*well_index + p]; + primary_variables_[XvarWell] += scalingFactor(p) * well_state.wellRates()[np*well_index + p]; } } break; @@ -1910,17 +1880,17 @@ namespace Opm double tot_well_rate = 0.0; for (int p = 0; p < np; ++p) { - tot_well_rate += g[p] * well_state.wellRates()[np*well_index + p]; + tot_well_rate += scalingFactor(p) * well_state.wellRates()[np*well_index + p]; } if(std::abs(tot_well_rate) > 0) { if (active()[ Water ]) { - primary_variables_[WFrac] = g[pu.phase_pos[Water]] * well_state.wellRates()[np*well_index + pu.phase_pos[Water]] / tot_well_rate; + primary_variables_[WFrac] = scalingFactor(pu.phase_pos[Water]) * well_state.wellRates()[np*well_index + pu.phase_pos[Water]] / tot_well_rate; } if (active()[ Gas ]) { - primary_variables_[GFrac] = g[ pu.phase_pos[Gas]] * (well_state.wellRates()[np*well_index + pu.phase_pos[Gas]] - well_state.solventWellRate(well_index)) / tot_well_rate ; + primary_variables_[GFrac] = scalingFactor(pu.phase_pos[Gas]) * (well_state.wellRates()[np*well_index + pu.phase_pos[Gas]] - well_state.solventWellRate(well_index)) / tot_well_rate ; } if (has_solvent) { - primary_variables_[SFrac] = g[ pu.phase_pos[Gas]] * well_state.solventWellRate(well_index) / tot_well_rate ; + primary_variables_[SFrac] = scalingFactor(pu.phase_pos[Gas]) * well_state.solventWellRate(well_index) / tot_well_rate ; } } else { // tot_well_rate == 0 if (well_type_ == INJECTOR) { @@ -2062,7 +2032,15 @@ namespace Opm double StandardWell::scalingFactor(const int phaseIdx) const { - //std::vector g = {1,1,0.01,0.01}; + const WellControls* wc = well_controls_; + const double* distr = well_controls_get_current_distr(wc); + + if (well_controls_get_current_type(wc) == RESERVOIR_RATE) { + if (has_solvent && phaseIdx == contiSolventEqIdx ) + OPM_THROW(std::runtime_error, "RESERVOIR_RATE control in combination with solvent is not implemented"); + + return distr[phaseIdx]; + } const auto& pu = phaseUsage(); if (active()[Water] && pu.phase_pos[Water] == phaseIdx) return 1.0; From b452d16f46b1c41ba060c2711da96432ca9bd1de Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Tue, 19 Sep 2017 11:16:34 +0200 Subject: [PATCH 3/5] Fix parallell run --- 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 2514c2bce..ba9d59204 100644 --- a/opm/autodiff/StandardWellsDense_impl.hpp +++ b/opm/autodiff/StandardWellsDense_impl.hpp @@ -45,13 +45,13 @@ namespace Opm { // has to be set always for the convergence check! global_nc_ = global_nc; + phase_usage_ = phase_usage_arg; + active_ = active_arg; + if ( ! localWellsActive() ) { return; } - phase_usage_ = phase_usage_arg; - active_ = active_arg; - calculateEfficiencyFactors(); #ifndef NDEBUG From 992ab84435a29f37640281e98603bd5a6a6895e7 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Wed, 20 Sep 2017 08:11:44 +0200 Subject: [PATCH 4/5] Remove support and tests for pseudo 2p model in flow ebos --- compareECLFiles.cmake | 1 - opm/autodiff/StandardWell_impl.hpp | 6 ------ 2 files changed, 7 deletions(-) diff --git a/compareECLFiles.cmake b/compareECLFiles.cmake index 392a7e27d..78e2903d6 100644 --- a/compareECLFiles.cmake +++ b/compareECLFiles.cmake @@ -106,7 +106,6 @@ set(rel_tol 1e-5) add_test_compareECLFiles(spe1 SPE1CASE2 flow ${abs_tol} ${rel_tol} compareECLFiles "") add_test_compareECLFiles(spe1 SPE1CASE2 flow_ebos ${abs_tol} ${rel_tol} compareECLFiles "") add_test_compareECLFiles(spe1_2p SPE1CASE2_2P flow ${abs_tol} ${rel_tol} compareECLFiles "" spe1) -add_test_compareECLFiles(spe1_2p SPE1CASE2_2P flow_ebos ${abs_tol} ${rel_tol} compareECLFiles "" spe1) add_test_compareECLFiles(spe1_2p SPE1CASE2_2P flow_ebos_2p ${abs_tol} ${rel_tol} compareECLFiles "" spe1) add_test_compareECLFiles(spe1 SPE1CASE2 flow_legacy ${abs_tol} ${rel_tol} compareECLFiles "") add_test_compareECLFiles(spe1_2p SPE1CASE2_2P flow_legacy ${abs_tol} ${rel_tol} compareECLFiles "" spe1) diff --git a/opm/autodiff/StandardWell_impl.hpp b/opm/autodiff/StandardWell_impl.hpp index 529e2ba82..5acdbcc5a 100644 --- a/opm/autodiff/StandardWell_impl.hpp +++ b/opm/autodiff/StandardWell_impl.hpp @@ -577,12 +577,6 @@ namespace Opm } } - // add a trivial equation for the dummy phase for 2p cases (Only support water + oil) - if ( numComp < numWellEq ) { - assert(!active()[ Gas ]); - invDuneD_[0][0][Gas][Gas] = 1.0; - } - // Store the perforation phase flux for later usage. if (has_solvent && componentIdx == contiSolventEqIdx) {// if (flowPhaseToEbosCompIdx(componentIdx) == Solvent) well_state.perfRateSolvent()[first_perf_ + perf] = cq_s[componentIdx].value(); From 5839962e31b1bade3f034de61de9f09538486263 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Thu, 21 Sep 2017 12:16:46 +0200 Subject: [PATCH 5/5] Some trivial cleaning --- opm/autodiff/StandardWell.hpp | 1 - opm/autodiff/StandardWell_impl.hpp | 2 +- opm/autodiff/WellInterface.hpp | 2 -- opm/autodiff/WellInterface_impl.hpp | 9 ++++++--- 4 files changed, 7 insertions(+), 7 deletions(-) diff --git a/opm/autodiff/StandardWell.hpp b/opm/autodiff/StandardWell.hpp index 7c72823e5..b32b25344 100644 --- a/opm/autodiff/StandardWell.hpp +++ b/opm/autodiff/StandardWell.hpp @@ -159,7 +159,6 @@ namespace Opm using Base::getAllowCrossFlow; using Base::phaseUsage; using Base::active; - using Base::flowToEbosPvIdx; using Base::flowPhaseToEbosPhaseIdx; using Base::flowPhaseToEbosCompIdx; using Base::numComponents; diff --git a/opm/autodiff/StandardWell_impl.hpp b/opm/autodiff/StandardWell_impl.hpp index 5acdbcc5a..c91f74c39 100644 --- a/opm/autodiff/StandardWell_impl.hpp +++ b/opm/autodiff/StandardWell_impl.hpp @@ -2045,7 +2045,7 @@ namespace Opm if (has_solvent && phaseIdx == contiSolventEqIdx ) return 0.01; - // we should come this fare + // we should not come this far assert(false); return 1.0; } diff --git a/opm/autodiff/WellInterface.hpp b/opm/autodiff/WellInterface.hpp index 950fca584..5d096edfc 100644 --- a/opm/autodiff/WellInterface.hpp +++ b/opm/autodiff/WellInterface.hpp @@ -268,8 +268,6 @@ namespace Opm int flowPhaseToEbosCompIdx( const int phaseIdx ) const; - int flowToEbosPvIdx( const int flowPv ) const; - int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const; // TODO: it is dumplicated with StandardWellsDense diff --git a/opm/autodiff/WellInterface_impl.hpp b/opm/autodiff/WellInterface_impl.hpp index e2c958f1e..def98243d 100644 --- a/opm/autodiff/WellInterface_impl.hpp +++ b/opm/autodiff/WellInterface_impl.hpp @@ -245,12 +245,15 @@ namespace Opm flowPhaseToEbosPhaseIdx( const int phaseIdx ) const { const auto& pu = phaseUsage(); - if (active()[Water] && pu.phase_pos[Water] == phaseIdx) + if (active()[Water] && pu.phase_pos[Water] == phaseIdx) { return FluidSystem::waterPhaseIdx; - if (active()[Oil] && pu.phase_pos[Oil] == phaseIdx) + } + if (active()[Oil] && pu.phase_pos[Oil] == phaseIdx) { return FluidSystem::oilPhaseIdx; - if (active()[Gas] && pu.phase_pos[Gas] == phaseIdx) + } + if (active()[Gas] && pu.phase_pos[Gas] == phaseIdx) { return FluidSystem::gasPhaseIdx; + } assert(phaseIdx < 3); // for other phases return the index