From d873ae165d65ae66005decbc7b4861dc4ff157bd Mon Sep 17 00:00:00 2001 From: "Halvor M. Nilsen" Date: Wed, 9 Oct 2019 15:24:23 +0200 Subject: [PATCH 1/9] Enable single-phase runs. --- ebos/ecloutputblackoilmodule.hh | 11 +- ebos/eclproblem.hh | 17 +- opm/core/props/phaseUsageFromDeck.hpp | 7 +- opm/simulators/flow/BlackoilModelEbos.hpp | 43 +-- .../timestepping/AdaptiveSimulatorTimer.cpp | 6 +- .../timestepping/AdaptiveTimeSteppingEbos.hpp | 5 +- .../timestepping/TimeStepControl.cpp | 8 +- .../wells/BlackoilWellModel_impl.hpp | 20 +- opm/simulators/wells/RateConverter.hpp | 12 +- opm/simulators/wells/StandardWell.hpp | 3 + opm/simulators/wells/StandardWell_impl.hpp | 294 ++++++++++-------- opm/simulators/wells/WellInterface.hpp | 9 +- 12 files changed, 264 insertions(+), 171 deletions(-) diff --git a/ebos/ecloutputblackoilmodule.hh b/ebos/ecloutputblackoilmodule.hh index 46ec2840f..e963ee957 100644 --- a/ebos/ecloutputblackoilmodule.hh +++ b/ebos/ecloutputblackoilmodule.hh @@ -395,7 +395,16 @@ public: } if (oilPressure_.size() > 0) { - oilPressure_[globalDofIdx] = Opm::getValue(fs.pressure(oilPhaseIdx)); + if (FluidSystem::phaseIsActive(oilPhaseIdx)) { + oilPressure_[globalDofIdx] = Opm::getValue(fs.pressure(oilPhaseIdx)); + }else{ + // put pressure in oil pressure for output + if (FluidSystem::phaseIsActive(waterPhaseIdx)) { + oilPressure_[globalDofIdx] = Opm::getValue(fs.pressure(waterPhaseIdx)); + } else { + oilPressure_[globalDofIdx] = Opm::getValue(fs.pressure(gasPhaseIdx)); + } + } Opm::Valgrind::CheckDefined(oilPressure_[globalDofIdx]); } diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index 1adfc1c88..640668f30 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -665,17 +665,20 @@ public: minTimeStepSize_ = tuning.getTSMINZ(0); } + initFluidSystem_(); + // deal with DRSDT unsigned ntpvt = eclState.runspec().tabdims().getNumPVTTables(); - maxDRs_.resize(ntpvt, 1e30); - dRsDtOnlyFreeGas_.resize(ntpvt, false); size_t numDof = this->model().numGridDof(); - lastRs_.resize(numDof, 0.0); - maxDRv_.resize(ntpvt, 1e30); - lastRv_.resize(numDof, 0.0); - maxOilSaturation_.resize(numDof, 0.0); + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { + maxDRs_.resize(ntpvt, 1e30); + dRsDtOnlyFreeGas_.resize(ntpvt, false); + lastRs_.resize(numDof, 0.0); + maxDRv_.resize(ntpvt, 1e30); + lastRv_.resize(numDof, 0.0); + maxOilSaturation_.resize(numDof, 0.0); + } - initFluidSystem_(); updateElementDepths_(); readRockParameters_(); readMaterialParameters_(); diff --git a/opm/core/props/phaseUsageFromDeck.hpp b/opm/core/props/phaseUsageFromDeck.hpp index edab3133d..b96d2a8ad 100644 --- a/opm/core/props/phaseUsageFromDeck.hpp +++ b/opm/core/props/phaseUsageFromDeck.hpp @@ -54,14 +54,9 @@ namespace Opm } } - // Only 2 or 3 phase systems handled. - if (pu.num_phases < 2 || pu.num_phases > 3) { - OPM_THROW(std::runtime_error, "Cannot handle cases with " << pu.num_phases << " phases."); - } - // We need oil systems, since we do not support the keywords needed for // water-gas systems. - if (!pu.phase_used[BlackoilPhases::Liquid]) { + if (!pu.phase_used[BlackoilPhases::Liquid] && !(pu.num_phases == 1)) { OPM_THROW(std::runtime_error, "Cannot handle cases with no OIL, i.e. water-gas systems."); } diff --git a/opm/simulators/flow/BlackoilModelEbos.hpp b/opm/simulators/flow/BlackoilModelEbos.hpp index 45057db76..edcb40672 100644 --- a/opm/simulators/flow/BlackoilModelEbos.hpp +++ b/opm/simulators/flow/BlackoilModelEbos.hpp @@ -422,29 +422,34 @@ namespace Opm { Scalar saturationsOld[FluidSystem::numPhases] = { 0.0 }; Scalar oilSaturationOld = 1.0; - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - saturationsOld[FluidSystem::waterPhaseIdx] = priVarsOld[Indices::waterSaturationIdx]; - oilSaturationOld -= saturationsOld[FluidSystem::waterPhaseIdx]; - } - - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && priVarsOld.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) { - saturationsOld[FluidSystem::gasPhaseIdx] = priVarsOld[Indices::compositionSwitchIdx]; - oilSaturationOld -= saturationsOld[FluidSystem::gasPhaseIdx]; - } - - if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { - saturationsOld[FluidSystem::oilPhaseIdx] = oilSaturationOld; - } + if (FluidSystem::numActivePhases() > 1) { + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + saturationsOld[FluidSystem::waterPhaseIdx] = priVarsOld[Indices::waterSaturationIdx]; + oilSaturationOld -= saturationsOld[FluidSystem::waterPhaseIdx]; + } + + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && + priVarsOld.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) + { + saturationsOld[FluidSystem::gasPhaseIdx] = priVarsOld[Indices::compositionSwitchIdx]; + oilSaturationOld -= saturationsOld[FluidSystem::gasPhaseIdx]; + } + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { + saturationsOld[FluidSystem::oilPhaseIdx] = oilSaturationOld; + } + for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++ phaseIdx) { + Scalar tmp = saturationsNew[phaseIdx] - saturationsOld[phaseIdx]; + resultDelta += tmp*tmp; + resultDenom += saturationsNew[phaseIdx]*saturationsNew[phaseIdx]; + assert(std::isfinite(resultDelta)); + assert(std::isfinite(resultDenom)); + } + } + // NB fix me! adding pressures changes to satutation changes does not make sense Scalar tmp = pressureNew - pressureOld; resultDelta += tmp*tmp; resultDenom += pressureNew*pressureNew; - - for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++ phaseIdx) { - const Scalar tmpSat = saturationsNew[phaseIdx] - saturationsOld[phaseIdx]; - resultDelta += tmpSat * tmpSat; - resultDenom += saturationsNew[phaseIdx]*saturationsNew[phaseIdx]; - } } resultDelta = gridView.comm().sum(resultDelta); diff --git a/opm/simulators/timestepping/AdaptiveSimulatorTimer.cpp b/opm/simulators/timestepping/AdaptiveSimulatorTimer.cpp index 43486f91f..93c74c79f 100644 --- a/opm/simulators/timestepping/AdaptiveSimulatorTimer.cpp +++ b/opm/simulators/timestepping/AdaptiveSimulatorTimer.cpp @@ -60,6 +60,7 @@ AdaptiveSimulatorTimer& AdaptiveSimulatorTimer::operator++ () { ++current_step_; current_time_ += dt_; + assert(dt_ > 0); // store used time step sizes steps_.push_back( dt_ ); return *this; @@ -71,7 +72,7 @@ AdaptiveSimulatorTimer& AdaptiveSimulatorTimer::operator++ () double remaining = (total_time_ - current_time_); // apply max time step if it was set dt_ = std::min( dt_estimate, max_time_step_ ); - + assert(dt_ > 0); if( remaining > 0 ) { // set new time step (depending on remaining time) @@ -81,6 +82,7 @@ AdaptiveSimulatorTimer& AdaptiveSimulatorTimer::operator++ () if( dt_ > max_time_step_ ) { dt_ = 0.5 * remaining; } + assert(dt_ > 0); return; } @@ -89,6 +91,7 @@ AdaptiveSimulatorTimer& AdaptiveSimulatorTimer::operator++ () if( 1.5 * dt_ > remaining ) { dt_ = 0.5 * remaining; + assert(dt_ > 0); return; } } @@ -102,6 +105,7 @@ AdaptiveSimulatorTimer& AdaptiveSimulatorTimer::operator++ () double AdaptiveSimulatorTimer::currentStepLength () const { + assert(dt_ > 0); return dt_; } diff --git a/opm/simulators/timestepping/AdaptiveTimeSteppingEbos.hpp b/opm/simulators/timestepping/AdaptiveTimeSteppingEbos.hpp index f3587231a..632ede844 100644 --- a/opm/simulators/timestepping/AdaptiveTimeSteppingEbos.hpp +++ b/opm/simulators/timestepping/AdaptiveTimeSteppingEbos.hpp @@ -293,9 +293,10 @@ namespace Opm { double dtEstimate = timeStepControl_->computeTimeStepSize(dt, iterations, relativeChange, substepTimer.simulationTimeElapsed()); + assert(dtEstimate > 0); // limit the growth of the timestep size by the growth factor dtEstimate = std::min(dtEstimate, double(maxGrowth_ * dt)); - + assert(dtEstimate > 0); // further restrict time step size growth after convergence problems if (restarts > 0) { dtEstimate = std::min(growthFactor_ * dt, dtEstimate); @@ -309,7 +310,7 @@ namespace Opm { const double maxPredictionTHPTimestep = 16.0 * unit::day; dtEstimate = std::min(dtEstimate, maxPredictionTHPTimestep); } - + assert(dtEstimate > 0); if (timestepVerbose_) { std::ostringstream ss; substepReport.reportStep(ss); diff --git a/opm/simulators/timestepping/TimeStepControl.cpp b/opm/simulators/timestepping/TimeStepControl.cpp index ff3521a83..c900603e2 100644 --- a/opm/simulators/timestepping/TimeStepControl.cpp +++ b/opm/simulators/timestepping/TimeStepControl.cpp @@ -128,13 +128,17 @@ namespace Opm { // shift errors for( int i=0; i<2; ++i ) { - errors_[ i ] = errors_[i+1]; + errors_[ i ] = errors_[i+1]; } // store new error const double error = relChange.relativeChange(); errors_[ 2 ] = error; - + for( int i=0; i<2; ++i ) { + assert(std::isfinite(errors_[i])); + assert(errors_[i]>0); + } + if( error > tol_ ) { // adjust dt by given tolerance diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index 4e6fc5db5..c17fba4b2 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -251,9 +251,19 @@ namespace Opm { const unsigned cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0); const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0); const auto& fs = intQuants.fluidState(); - - const double p = fs.pressure(FluidSystem::oilPhaseIdx).value(); - cellPressures[cellIdx] = p; + // copy of get perfpressure in Standard well + // exept for value + double perf_pressure = 0.0; + if (Indices::oilEnabled) { + perf_pressure = fs.pressure(FluidSystem::oilPhaseIdx).value(); + } else { + if (Indices::waterEnabled) { + perf_pressure = fs.pressure(FluidSystem::waterPhaseIdx).value(); + } else { + perf_pressure = fs.pressure(FluidSystem::gasPhaseIdx).value(); + } + } + cellPressures[cellIdx] = perf_pressure; } well_state_.init(wells(), cellPressures, schedule(), wells_ecl_, timeStepIdx, &previous_well_state_, phase_usage_); @@ -1525,8 +1535,8 @@ namespace Opm { int BlackoilWellModel::numComponents() const { - if (numPhases() == 2) { - return 2; + if (numWells() > 0 && numPhases() < 3) { + return numPhases(); } int numComp = FluidSystem::numComponents; if (has_solvent_) { diff --git a/opm/simulators/wells/RateConverter.hpp b/opm/simulators/wells/RateConverter.hpp index 921fdac3f..1d58729c7 100644 --- a/opm/simulators/wells/RateConverter.hpp +++ b/opm/simulators/wells/RateConverter.hpp @@ -493,11 +493,13 @@ namespace Opm { // sum p, rs, rv, and T. double hydrocarbonPV = pv_cell*hydrocarbon; - pv += hydrocarbonPV; - p += fs.pressure(FluidSystem::oilPhaseIdx).value()*hydrocarbonPV; - rs += fs.Rs().value()*hydrocarbonPV; - rv += fs.Rv().value()*hydrocarbonPV; - T += fs.temperature(FluidSystem::oilPhaseIdx).value()*hydrocarbonPV; + if (hydrocarbonPV > 0) { + pv += hydrocarbonPV; + p += fs.pressure(FluidSystem::oilPhaseIdx).value()*hydrocarbonPV; + rs += fs.Rs().value()*hydrocarbonPV; + rv += fs.Rv().value()*hydrocarbonPV; + T += fs.temperature(FluidSystem::oilPhaseIdx).value()*hydrocarbonPV; + } } for (const auto& reg : rmap_.activeRegions()) { diff --git a/opm/simulators/wells/StandardWell.hpp b/opm/simulators/wells/StandardWell.hpp index 7b2979e79..d2c31ced8 100644 --- a/opm/simulators/wells/StandardWell.hpp +++ b/opm/simulators/wells/StandardWell.hpp @@ -59,6 +59,7 @@ namespace Opm using typename Base::Indices; using typename Base::RateConverterType; using typename Base::SparseMatrixAdapter; + using typename Base::FluidState; using Base::numEq; @@ -293,6 +294,8 @@ namespace Opm EvalWell extendEval(const Eval& in) const; + Eval getPerfPressure(const FluidState& fs) const; + // xw = inv(D)*(rw - C*x) void recoverSolutionWell(const BVector& x, BVectorWell& xw) const; diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index 533cedf8f..f8bfe2005 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -238,6 +238,10 @@ namespace Opm StandardWell:: wellVolumeFraction(const unsigned compIdx) const { + if (FluidSystem::numActivePhases() == 1) { + return EvalWell(numWellEq_ + numEq, 1.0); + } + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx)) { return primary_variables_evaluation_[WFrac]; } @@ -305,6 +309,24 @@ namespace Opm + template + typename StandardWell::Eval + StandardWell::getPerfPressure(const typename StandardWell::FluidState& fs) const + { + Eval pressure; + if (Indices::oilEnabled) { + pressure = fs.pressure(FluidSystem::oilPhaseIdx); + } else { + if (Indices::waterEnabled) { + pressure = fs.pressure(FluidSystem::waterPhaseIdx); + } else { + pressure = fs.pressure(FluidSystem::gasPhaseIdx); + } + } + return pressure; + } + + template void StandardWell:: @@ -321,7 +343,7 @@ namespace Opm { const auto& fs = intQuants.fluidState(); - const EvalWell pressure = extendEval(fs.pressure(FluidSystem::oilPhaseIdx)); + const EvalWell pressure = extendEval(getPerfPressure(fs)); const EvalWell rs = extendEval(fs.Rs()); const EvalWell rv = extendEval(fs.Rv()); std::vector b_perfcells_dense(num_components_, EvalWell{numWellEq_ + numEq, 0.0}); @@ -669,7 +691,9 @@ namespace Opm || (pu.phase_pos[Gas] == p && summaryConfig.hasSummaryKey("WPIG:" + name()))) { const unsigned int compIdx = flowPhaseToEbosCompIdx(p); - const double drawdown = well_state.perfPress()[first_perf_ + perf] - intQuants.fluidState().pressure(FluidSystem::oilPhaseIdx).value(); + const auto& fs = intQuants.fluidState(); + Eval perf_pressure = getPerfPressure(fs); + const double drawdown = well_state.perfPress()[first_perf_ + perf] - perf_pressure.value(); const bool new_well = schedule.hasWellEvent(name(), ScheduleEvents::NEW_WELL, current_step_); double productivity_index = cq_s[compIdx].value() / drawdown; scaleProductivityIndex(perf, productivity_index, new_well, deferred_logger); @@ -679,12 +703,15 @@ namespace Opm } - // add vol * dF/dt + Q to the well equations; for (int componentIdx = 0; componentIdx < numWellConservationEq; ++componentIdx) { // TODO: following the development in MSW, we need to convert the volume of the wellbore to be surface volume // since all the rates are under surface condition - EvalWell resWell_loc = (wellSurfaceVolumeFraction(componentIdx) - F0_[componentIdx]) * volume / dt; + EvalWell resWell_loc(numWellEq_ + numEq, 0.0); + if (FluidSystem::numActivePhases() > 1) { + assert(dt > 0); + resWell_loc += (wellSurfaceVolumeFraction(componentIdx) - F0_[componentIdx]) * volume / dt; + } resWell_loc -= getQs(componentIdx) * well_efficiency_factor_; for (int pvIdx = 0; pvIdx < numWellEq_; ++pvIdx) { invDuneD_[0][0][componentIdx][pvIdx] += resWell_loc.derivative(pvIdx+numEq); @@ -919,29 +946,30 @@ namespace Opm const double relaxation_factor_fractions = (well_type_ == PRODUCER) ? relaxationFactorFractionsProducer(old_primary_variables, dwells) : 1.0; + if (FluidSystem::numActivePhases() > 1) { + // update the second and third well variable (The flux fractions) + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + const int sign2 = dwells[0][WFrac] > 0 ? 1: -1; + const double dx2_limited = sign2 * std::min(std::abs(dwells[0][WFrac] * relaxation_factor_fractions), dFLimit); + // primary_variables_[WFrac] = old_primary_variables[WFrac] - dx2_limited; + primary_variables_[WFrac] = old_primary_variables[WFrac] - dx2_limited; + } - // update the second and third well variable (The flux fractions) - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - const int sign2 = dwells[0][WFrac] > 0 ? 1: -1; - const double dx2_limited = sign2 * std::min(std::abs(dwells[0][WFrac] * relaxation_factor_fractions), dFLimit); - // primary_variables_[WFrac] = old_primary_variables[WFrac] - dx2_limited; - primary_variables_[WFrac] = old_primary_variables[WFrac] - dx2_limited; + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const int sign3 = dwells[0][GFrac] > 0 ? 1: -1; + const double dx3_limited = sign3 * std::min(std::abs(dwells[0][GFrac] * relaxation_factor_fractions), dFLimit); + primary_variables_[GFrac] = old_primary_variables[GFrac] - dx3_limited; + } + + if (has_solvent) { + const int sign4 = dwells[0][SFrac] > 0 ? 1: -1; + const double dx4_limited = sign4 * std::min(std::abs(dwells[0][SFrac]) * relaxation_factor_fractions, dFLimit); + primary_variables_[SFrac] = old_primary_variables[SFrac] - dx4_limited; + } + + processFractions(); } - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - const int sign3 = dwells[0][GFrac] > 0 ? 1: -1; - const double dx3_limited = sign3 * std::min(std::abs(dwells[0][GFrac] * relaxation_factor_fractions), dFLimit); - primary_variables_[GFrac] = old_primary_variables[GFrac] - dx3_limited; - } - - if (has_solvent) { - const int sign4 = dwells[0][SFrac] > 0 ? 1: -1; - const double dx4_limited = sign4 * std::min(std::abs(dwells[0][SFrac]) * relaxation_factor_fractions, dFLimit); - primary_variables_[SFrac] = old_primary_variables[SFrac] - dx4_limited; - } - - processFractions(); - // updating the total rates Q_t const double relaxation_factor_rate = relaxationFactorRate(old_primary_variables, dwells); primary_variables_[WQTotal] = old_primary_variables[WQTotal] - dwells[0][WQTotal] * relaxation_factor_rate; @@ -956,6 +984,10 @@ namespace Opm } updateExtraPrimaryVariables(dwells); + + for (int i = 0; i < primary_variables_.size(); ++i) { + assert(Opm::isfinite(primary_variables_[i])); + } } @@ -1073,49 +1105,53 @@ namespace Opm updateWellStateFromPrimaryVariables(WellState& well_state, Opm::DeferredLogger& deferred_logger) const { const PhaseUsage& pu = phaseUsage(); - assert( FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) ); - const int oil_pos = pu.phase_pos[Oil]; std::vector F(number_of_phases_, 0.0); - F[oil_pos] = 1.0; + if (FluidSystem::numActivePhases() == 1) { + F[0] = 1.0; + } else { + assert( FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) ); + const int oil_pos = pu.phase_pos[Oil]; + F[oil_pos] = 1.0; - if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) { - const int water_pos = pu.phase_pos[Water]; - F[water_pos] = primary_variables_[WFrac]; - F[oil_pos] -= F[water_pos]; - } - - if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) { - const int gas_pos = pu.phase_pos[Gas]; - F[gas_pos] = primary_variables_[GFrac]; - F[oil_pos] -= F[gas_pos]; - } - - double F_solvent = 0.0; - if (has_solvent) { - F_solvent = primary_variables_[SFrac]; - F[oil_pos] -= F_solvent; - } - - // convert the fractions to be Q_p / G_total to calculate the phase rates - for (int p = 0; p < number_of_phases_; ++p) { - const double scal = scalingFactor(p); - // for injection wells, there should only one non-zero scaling factor - if (scal > 0) { - F[p] /= scal ; - } else { - // this should only happens to injection wells - F[p] = 0.; + if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) { + const int water_pos = pu.phase_pos[Water]; + F[water_pos] = primary_variables_[WFrac]; + F[oil_pos] -= F[water_pos]; } - } - // 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_solvent /= scalingFactor(contiSolventEqIdx); - F[pu.phase_pos[Gas]] += F_solvent; - } + if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) { + const int gas_pos = pu.phase_pos[Gas]; + F[gas_pos] = primary_variables_[GFrac]; + F[oil_pos] -= F[gas_pos]; + } + double F_solvent = 0.0; + if (has_solvent) { + F_solvent = primary_variables_[SFrac]; + F[oil_pos] -= F_solvent; + } + + // convert the fractions to be Q_p / G_total to calculate the phase rates + for (int p = 0; p < number_of_phases_; ++p) { + const double scal = scalingFactor(p); + // for injection wells, there should only one non-zero scaling factor + if (scal > 0) { + F[p] /= scal ; + } else { + // this should only happens to injection wells + F[p] = 0.; + } + } + + // 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_solvent /= scalingFactor(contiSolventEqIdx); + F[pu.phase_pos[Gas]] += F_solvent; + } + + } well_state.bhp()[index_of_well_] = primary_variables_[Bhp]; // calculate the phase rates based on the primary variables @@ -1292,6 +1328,10 @@ namespace Opm break; } // end of switch + + for (int i = 0; i < primary_variables_.size(); ++i) { + assert(Opm::isfinite(primary_variables_[i])); + } } @@ -1324,7 +1364,8 @@ namespace Opm const auto& int_quantities = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)); const auto& fs = int_quantities.fluidState(); // the pressure of the reservoir grid block the well connection is in - const double p_r = fs.pressure(FluidSystem::oilPhaseIdx).value(); + Eval perf_pressure = getPerfPressure(fs); + double p_r = perf_pressure.value(); // calculating the b for the connection std::vector b_perf(num_components_); @@ -1550,8 +1591,9 @@ namespace Opm const int cell_idx = well_cells_[perf]; const auto& intQuants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); const auto& fs = intQuants.fluidState(); + Eval perf_pressure = getPerfPressure(fs); + const double pressure = perf_pressure.value(); - const double pressure = (fs.pressure(FluidSystem::oilPhaseIdx)).value(); const double bhp = getBhp().value(); // Pressure drawdown (also used to determine direction of flow) @@ -1755,6 +1797,7 @@ namespace Opm // TODO: to check why should be perf - 1 const double p_above = perf == 0 ? well_state.bhp()[w] : well_state.perfPress()[first_perf_ + perf - 1]; const double p_avg = (well_state.perfPress()[first_perf_ + perf] + p_above)/2; + // strange choice const double temperature = fs.temperature(FluidSystem::oilPhaseIdx).value(); if (waterPresent) { @@ -2518,57 +2561,57 @@ namespace Opm } } - if (std::abs(total_well_rate) > 0.) { - const auto pu = phaseUsage(); - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - primary_variables_[WFrac] = scalingFactor(pu.phase_pos[Water]) * well_state.wellRates()[np*well_index + pu.phase_pos[Water]] / total_well_rate; - } - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - primary_variables_[GFrac] = scalingFactor(pu.phase_pos[Gas]) * (well_state.wellRates()[np*well_index + pu.phase_pos[Gas]] - well_state.solventWellRate(well_index)) / total_well_rate ; - } - if (has_solvent) { - primary_variables_[SFrac] = scalingFactor(pu.phase_pos[Gas]) * well_state.solventWellRate(well_index) / total_well_rate ; - } - } else { // total_well_rate == 0 - if (well_type_ == INJECTOR) { - auto phase = well_ecl_.getInjectionProperties().injectorType; - // only single phase injection handled + if (FluidSystem::numActivePhases() > 1) { + if (std::abs(total_well_rate) > 0.) { + const auto pu = phaseUsage(); if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - if (phase == Well2::InjectorType::WATER) { - primary_variables_[WFrac] = 1.0; - } else { - primary_variables_[WFrac] = 0.0; - } + primary_variables_[WFrac] = scalingFactor(pu.phase_pos[Water]) * well_state.wellRates()[np*well_index + pu.phase_pos[Water]] / total_well_rate; } - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - if (phase == Well2::InjectorType::GAS) { - primary_variables_[GFrac] = 1.0 - wsolvent(); - if (has_solvent) { - primary_variables_[SFrac] = wsolvent(); + primary_variables_[GFrac] = scalingFactor(pu.phase_pos[Gas]) * (well_state.wellRates()[np*well_index + pu.phase_pos[Gas]] - well_state.solventWellRate(well_index)) / total_well_rate ; + } + if (has_solvent) { + primary_variables_[SFrac] = scalingFactor(pu.phase_pos[Gas]) * well_state.solventWellRate(well_index) / total_well_rate ; + } + } else { // total_well_rate == 0 + if (well_type_ == INJECTOR) { + auto phase = well_ecl_.getInjectionProperties().injectorType; + // only single phase injection handled + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + if (phase == Well2::InjectorType::WATER) { + primary_variables_[WFrac] = 1.0; + } else { + primary_variables_[WFrac] = 0.0; + } + } + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + if (phase == Well2::InjectorType::GAS) { + primary_variables_[GFrac] = 1.0 - wsolvent(); + if (has_solvent) { + primary_variables_[SFrac] = wsolvent(); + } + } else { + primary_variables_[GFrac] = 0.0; } - } else { - primary_variables_[GFrac] = 0.0; } - } - // TODO: it is possible to leave injector as a oil well, - // when F_w and F_g both equals to zero, not sure under what kind of circumstance - // this will happen. - } else if (well_type_ == PRODUCER) { // producers - // TODO: the following are not addressed for the solvent case yet - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - primary_variables_[WFrac] = 1.0 / np; + // TODO: it is possible to leave injector as a oil well, + // when F_w and F_g both equals to zero, not sure under what kind of circumstance + // this will happen. + } else if (well_type_ == PRODUCER) { // producers + // TODO: the following are not addressed for the solvent case yet + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + primary_variables_[WFrac] = 1.0 / np; + } + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + primary_variables_[GFrac] = 1.0 / np; + } + } else { + OPM_DEFLOG_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well", deferred_logger); } - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - primary_variables_[GFrac] = 1.0 / np; - } - } else { - OPM_DEFLOG_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well", deferred_logger); } } - // BHP primary_variables_[Bhp] = well_state.bhp()[index_of_well_]; @@ -2579,6 +2622,10 @@ namespace Opm primary_variables_[Bhp + 1 + number_of_perforations_ + perf] = well_state.perfSkinPressure()[first_perf_ + perf]; } } + + for (int i = 0; i < primary_variables_.size(); ++i) { + assert(Opm::isfinite(primary_variables_[i])); + } } @@ -2824,28 +2871,31 @@ namespace Opm // 0.95 is a experimental value, which remains to be optimized double relaxation_factor = 1.0; - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - const double relaxation_factor_w = relaxationFactorFraction(primary_variables[WFrac], dwells[0][WFrac]); - relaxation_factor = std::min(relaxation_factor, relaxation_factor_w); - } + if (FluidSystem::numActivePhases() > 1) { + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + const double relaxation_factor_w = relaxationFactorFraction(primary_variables[WFrac], dwells[0][WFrac]); + relaxation_factor = std::min(relaxation_factor, relaxation_factor_w); + } - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - const double relaxation_factor_g = relaxationFactorFraction(primary_variables[GFrac], dwells[0][GFrac]); - relaxation_factor = std::min(relaxation_factor, relaxation_factor_g); - } + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const double relaxation_factor_g = relaxationFactorFraction(primary_variables[GFrac], dwells[0][GFrac]); + relaxation_factor = std::min(relaxation_factor, relaxation_factor_g); + } - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - // We need to make sure the even with the relaxation_factor, the sum of F_w and F_g is below one, so there will - // not be negative oil fraction later - const double original_sum = primary_variables[WFrac] + primary_variables[GFrac]; - const double relaxed_update = (dwells[0][WFrac] + dwells[0][GFrac]) * relaxation_factor; - const double possible_updated_sum = original_sum - relaxed_update; + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) + && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + // We need to make sure the even with the relaxation_factor, the sum of F_w and F_g is below one, so + // there will not be negative oil fraction later + const double original_sum = primary_variables[WFrac] + primary_variables[GFrac]; + const double relaxed_update = (dwells[0][WFrac] + dwells[0][GFrac]) * relaxation_factor; + const double possible_updated_sum = original_sum - relaxed_update; - if (possible_updated_sum > 1.0) { - assert(relaxed_update != 0.); + if (possible_updated_sum > 1.0) { + assert(relaxed_update != 0.); - const double further_relaxation_factor = std::abs((1. - original_sum) / relaxed_update) * 0.95; - relaxation_factor *= further_relaxation_factor; + const double further_relaxation_factor = std::abs((1. - original_sum) / relaxed_update) * 0.95; + relaxation_factor *= further_relaxation_factor; + } } } diff --git a/opm/simulators/wells/WellInterface.hpp b/opm/simulators/wells/WellInterface.hpp index 6af7ca93a..9e50fa3a8 100644 --- a/opm/simulators/wells/WellInterface.hpp +++ b/opm/simulators/wells/WellInterface.hpp @@ -94,6 +94,7 @@ namespace Opm static const bool has_solvent = GET_PROP_VALUE(TypeTag, EnableSolvent); static const bool has_polymer = GET_PROP_VALUE(TypeTag, EnablePolymer); static const bool has_energy = GET_PROP_VALUE(TypeTag, EnableEnergy); + static const bool has_temperature = GET_PROP_VALUE(TypeTag, EnableTemperature); // flag for polymer molecular weight related static const bool has_polymermw = GET_PROP_VALUE(TypeTag, EnablePolymerMW); static const bool has_foam = GET_PROP_VALUE(TypeTag, EnableFoam); @@ -106,7 +107,13 @@ namespace Opm // For the conversion between the surface volume rate and reservoir voidage rate using RateConverterType = RateConverter:: SurfaceToReservoirVoidage >; - + static const bool compositionSwitchEnabled = Indices::gasEnabled; + using FluidState = Opm::BlackOilFluidState; /// Constructor WellInterface(const Well2& well, const int time_step, const Wells* wells, const ModelParameters& param, From 2d9f41b1aa730565e51fc04443931e41fbb6bb3b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 10 Oct 2019 16:10:33 +0200 Subject: [PATCH 2/9] Rename getPerfPressure() -> getPerfCellPressure(). --- opm/simulators/wells/StandardWell.hpp | 2 +- opm/simulators/wells/StandardWell_impl.hpp | 10 +++++----- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/opm/simulators/wells/StandardWell.hpp b/opm/simulators/wells/StandardWell.hpp index d2c31ced8..0e72faecd 100644 --- a/opm/simulators/wells/StandardWell.hpp +++ b/opm/simulators/wells/StandardWell.hpp @@ -294,7 +294,7 @@ namespace Opm EvalWell extendEval(const Eval& in) const; - Eval getPerfPressure(const FluidState& fs) const; + Eval getPerfCellPressure(const FluidState& fs) const; // xw = inv(D)*(rw - C*x) void recoverSolutionWell(const BVector& x, BVectorWell& xw) const; diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index f8bfe2005..04bf02303 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -311,7 +311,7 @@ namespace Opm template typename StandardWell::Eval - StandardWell::getPerfPressure(const typename StandardWell::FluidState& fs) const + StandardWell::getPerfCellPressure(const typename StandardWell::FluidState& fs) const { Eval pressure; if (Indices::oilEnabled) { @@ -343,7 +343,7 @@ namespace Opm { const auto& fs = intQuants.fluidState(); - const EvalWell pressure = extendEval(getPerfPressure(fs)); + const EvalWell pressure = extendEval(getPerfCellPressure(fs)); const EvalWell rs = extendEval(fs.Rs()); const EvalWell rv = extendEval(fs.Rv()); std::vector b_perfcells_dense(num_components_, EvalWell{numWellEq_ + numEq, 0.0}); @@ -692,7 +692,7 @@ namespace Opm const unsigned int compIdx = flowPhaseToEbosCompIdx(p); const auto& fs = intQuants.fluidState(); - Eval perf_pressure = getPerfPressure(fs); + Eval perf_pressure = getPerfCellPressure(fs); const double drawdown = well_state.perfPress()[first_perf_ + perf] - perf_pressure.value(); const bool new_well = schedule.hasWellEvent(name(), ScheduleEvents::NEW_WELL, current_step_); double productivity_index = cq_s[compIdx].value() / drawdown; @@ -1364,7 +1364,7 @@ namespace Opm const auto& int_quantities = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)); const auto& fs = int_quantities.fluidState(); // the pressure of the reservoir grid block the well connection is in - Eval perf_pressure = getPerfPressure(fs); + Eval perf_pressure = getPerfCellPressure(fs); double p_r = perf_pressure.value(); // calculating the b for the connection @@ -1591,7 +1591,7 @@ namespace Opm const int cell_idx = well_cells_[perf]; const auto& intQuants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); const auto& fs = intQuants.fluidState(); - Eval perf_pressure = getPerfPressure(fs); + Eval perf_pressure = getPerfCellPressure(fs); const double pressure = perf_pressure.value(); const double bhp = getBhp().value(); From 3ddcdf2349d29c5096a611c51556bbcb32e8980f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 10 Oct 2019 16:11:09 +0200 Subject: [PATCH 3/9] Only run relperm diagnostics for > 1 phase cases. --- opm/simulators/flow/FlowMainEbos.hpp | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/opm/simulators/flow/FlowMainEbos.hpp b/opm/simulators/flow/FlowMainEbos.hpp index 9bb9017ad..fb9455323 100644 --- a/opm/simulators/flow/FlowMainEbos.hpp +++ b/opm/simulators/flow/FlowMainEbos.hpp @@ -432,9 +432,11 @@ namespace Opm return; } - // Run relperm diagnostics - RelpermDiagnostics diagnostic; - diagnostic.diagnosis(eclState(), deck(), this->grid()); + // Run relperm diagnostics if we have more than one phase. + if (FluidSystem::numActivePhases() > 1) { + RelpermDiagnostics diagnostic; + diagnostic.diagnosis(eclState(), deck(), this->grid()); + } } // Run the simulator. From 7270aa30135993bbc70270570b436b714b26d9ad Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 10 Oct 2019 16:11:56 +0200 Subject: [PATCH 4/9] Update flow_tag.hpp to match current flow.cpp. The flow_tag.hpp file is used for simple construction of simulators for particular models or other compile-time choices. This updates it to get identical behaviour as that of mainline (multi-model) Flow. --- flow/flow_blackoil_dunecpr.cpp | 4 +- flow/flow_tag.hpp | 235 +++++++++++++++++++++++++++++---- 2 files changed, 213 insertions(+), 26 deletions(-) diff --git a/flow/flow_blackoil_dunecpr.cpp b/flow/flow_blackoil_dunecpr.cpp index f6f08f472..ffa90f28f 100644 --- a/flow/flow_blackoil_dunecpr.cpp +++ b/flow/flow_blackoil_dunecpr.cpp @@ -99,7 +99,5 @@ namespace Opm { int main(int argc, char** argv) { typedef TTAG(EclFlowProblemSimple) TypeTag; - bool outputCout = true; - bool outputFiles = true; - return mainFlow(argc, argv, outputCout, outputFiles); + return mainFlow(argc, argv); } diff --git a/flow/flow_tag.hpp b/flow/flow_tag.hpp index ccc2ea16d..a568eaad6 100644 --- a/flow/flow_tag.hpp +++ b/flow/flow_tag.hpp @@ -29,10 +29,15 @@ #include #include +#include +#include +#include + #include #include #include #include +#include //#include //#include @@ -118,13 +123,162 @@ namespace detail throw std::invalid_argument( "Cannot find input case " + casename ); } + + + // This function is an extreme special case, if the program has been invoked + // *exactly* as: + // + // flow --version + // + // the call is intercepted by this function which will print "flow $version" + // on stdout and exit(0). + void handleVersionCmdLine(int argc, char** argv) { + for ( int i = 1; i < argc; ++i ) + { + if (std::strcmp(argv[i], "--version") == 0) { + std::cout << "flow " << Opm::moduleVersionName() << std::endl; + std::exit(EXIT_SUCCESS); + } + } + } + +} + + + +enum class FileOutputMode { + //! \brief No output to files. + OUTPUT_NONE = 0, + //! \brief Output only to log files, no eclipse output. + OUTPUT_LOG_ONLY = 1, + //! \brief Output to all files. + OUTPUT_ALL = 3 +}; + + + +void ensureOutputDirExists(const std::string& cmdline_output_dir) +{ + if (!boost::filesystem::is_directory(cmdline_output_dir)) { + try { + boost::filesystem::create_directories(cmdline_output_dir); + } + catch (...) { + throw std::runtime_error("Creation of output directory '" + cmdline_output_dir + "' failed\n"); + } + } +} + + +// Setup the OpmLog backends +FileOutputMode setupLogging(int mpi_rank_, const std::string& deck_filename, const std::string& cmdline_output_dir, const std::string& cmdline_output, bool output_cout_, const std::string& stdout_log_id) { + + if (!cmdline_output_dir.empty()) { + ensureOutputDirExists(cmdline_output_dir); + } + + // create logFile + using boost::filesystem::path; + path fpath(deck_filename); + std::string baseName; + std::ostringstream debugFileStream; + std::ostringstream logFileStream; + + // Strip extension "." or ".DATA" + std::string extension = boost::to_upper_copy(fpath.extension().string()); + if (extension == ".DATA" || extension == ".") { + baseName = boost::to_upper_copy(fpath.stem().string()); + } else { + baseName = boost::to_upper_copy(fpath.filename().string()); + } + + std::string output_dir = cmdline_output_dir; + if (output_dir.empty()) { + output_dir = absolute(path(baseName).parent_path()).string(); + } + + logFileStream << output_dir << "/" << baseName; + debugFileStream << output_dir << "/" << baseName; + + if (mpi_rank_ != 0) { + // Added rank to log file for non-zero ranks. + // This prevents message loss. + debugFileStream << "." << mpi_rank_; + // If the following file appears then there is a bug. + logFileStream << "." << mpi_rank_; + } + logFileStream << ".PRT"; + debugFileStream << ".DBG"; + + FileOutputMode output; + { + static std::map stringToOutputMode = + { {"none", FileOutputMode::OUTPUT_NONE }, + {"false", FileOutputMode::OUTPUT_LOG_ONLY }, + {"log", FileOutputMode::OUTPUT_LOG_ONLY }, + {"all" , FileOutputMode::OUTPUT_ALL }, + {"true" , FileOutputMode::OUTPUT_ALL }}; + auto outputModeIt = stringToOutputMode.find(cmdline_output); + if (outputModeIt != stringToOutputMode.end()) { + output = outputModeIt->second; + } + else { + output = FileOutputMode::OUTPUT_ALL; + std::cerr << "Value " << cmdline_output << + " is not a recognized output mode. Using \"all\" instead." + << std::endl; + } + } + + if (output > FileOutputMode::OUTPUT_NONE) { + std::shared_ptr prtLog = std::make_shared(logFileStream.str(), Opm::Log::NoDebugMessageTypes, false, output_cout_); + Opm::OpmLog::addBackend("ECLIPSEPRTLOG", prtLog); + prtLog->setMessageLimiter(std::make_shared()); + prtLog->setMessageFormatter(std::make_shared(false)); + } + + if (output >= FileOutputMode::OUTPUT_LOG_ONLY) { + std::string debugFile = debugFileStream.str(); + std::shared_ptr debugLog = std::make_shared(debugFileStream.str(), Opm::Log::DefaultMessageTypes, false, output_cout_); + Opm::OpmLog::addBackend("DEBUGLOG", debugLog); + } + + std::shared_ptr streamLog = std::make_shared(std::cout, Opm::Log::StdoutMessageTypes); + Opm::OpmLog::addBackend(stdout_log_id, streamLog); + streamLog->setMessageFormatter(std::make_shared(true)); + + return output; +} + + + +void setupMessageLimiter(const Opm::MessageLimits msgLimits, const std::string& stdout_log_id) { + std::shared_ptr stream_log = Opm::OpmLog::getBackend(stdout_log_id); + + const std::map limits = {{Opm::Log::MessageType::Note, + msgLimits.getCommentPrintLimit(0)}, + {Opm::Log::MessageType::Info, + msgLimits.getMessagePrintLimit(0)}, + {Opm::Log::MessageType::Warning, + msgLimits.getWarningPrintLimit(0)}, + {Opm::Log::MessageType::Error, + msgLimits.getErrorPrintLimit(0)}, + {Opm::Log::MessageType::Problem, + msgLimits.getProblemPrintLimit(0)}, + {Opm::Log::MessageType::Bug, + msgLimits.getBugPrintLimit(0)}}; + stream_log->setMessageLimiter(std::make_shared(10, limits)); } // ----------------- Main program ----------------- template -int mainFlow(int argc, char** argv, bool outputCout, bool outputFiles) +int mainFlow(int argc, char** argv) { + Dune::Timer externalSetupTimer; + externalSetupTimer.start(); + + detail::handleVersionCmdLine(argc, argv); // MPI setup. #if HAVE_DUNE_FEM Dune::Fem::MPIManager::initialize(argc, argv); @@ -151,15 +305,21 @@ int mainFlow(int argc, char** argv, bool outputCout, bool outputFiles) typedef TTAG(FlowEarlyBird) PreTypeTag; typedef GET_PROP_TYPE(PreTypeTag, Problem) PreProblem; - PreProblem::setBriefDescription("Simple Flow, an advanced reservoir simulator for ECL-decks provided by the Open Porous Media project."); + PreProblem::setBriefDescription("Flow, an advanced reservoir simulator for ECL-decks provided by the Open Porous Media project."); int status = Opm::FlowMainEbos::setupParameters_(argc, argv); - if (status != 0) + if (status != 0) { // if setupParameters_ returns a value smaller than 0, there was no error, but // the program should abort. This is the case e.g. for the --help and the // --print-properties parameters. +#if HAVE_MPI + MPI_Finalize(); +#endif return (status >= 0)?status:0; + } + FileOutputMode outputMode = FileOutputMode::OUTPUT_NONE; + bool outputCout = false; if (mpiRank == 0) outputCout = EWOMS_GET_PARAM(PreTypeTag, bool, EnableTerminalOutput); @@ -176,28 +336,57 @@ int mainFlow(int argc, char** argv, bool outputCout, bool outputFiles) // Create Deck and EclipseState. try { - Opm::Parser parser; - typedef std::pair ParseModePair; - typedef std::vector ParseModePairs; - ParseModePairs tmp; - tmp.push_back(ParseModePair(Opm::ParseContext::PARSE_RANDOM_SLASH, Opm::InputError::IGNORE)); - tmp.push_back(ParseModePair(Opm::ParseContext::PARSE_MISSING_DIMS_KEYWORD, Opm::InputError::WARN)); - tmp.push_back(ParseModePair(Opm::ParseContext::SUMMARY_UNKNOWN_WELL, Opm::InputError::WARN)); - tmp.push_back(ParseModePair(Opm::ParseContext::SUMMARY_UNKNOWN_GROUP, Opm::InputError::WARN)); - Opm::ParseContext parseContext(tmp); - Opm::ErrorGuard errorGuard; + if (outputCout) { + std::cout << "Reading deck file '" << deckFilename << "'\n"; + std::cout.flush(); + } + std::shared_ptr deck; + std::shared_ptr eclipseState; + std::shared_ptr schedule; + std::shared_ptr summaryConfig; + { + Opm::Parser parser; + Opm::ParseContext parseContext; + Opm::ErrorGuard errorGuard; + outputMode = setupLogging(mpiRank, + deckFilename, + EWOMS_GET_PARAM(PreTypeTag, std::string, OutputDir), + EWOMS_GET_PARAM(PreTypeTag, std::string, OutputMode), + outputCout, "STDOUT_LOGGER"); - std::shared_ptr deck = std::make_shared< Opm::Deck >( parser.parseFile(deckFilename , parseContext, errorGuard) ); - if ( outputCout ) { - Opm::checkDeck(*deck, parser, parseContext, errorGuard); - Opm::MissingFeatures::checkKeywords(*deck); - } - //Opm::Runspec runspec( *deck ); - //const auto& phases = runspec.phases(); + if (EWOMS_GET_PARAM(PreTypeTag, bool, EclStrictParsing)) + parseContext.update( Opm::InputError::DELAYED_EXIT1); + else { + parseContext.update(Opm::ParseContext::PARSE_RANDOM_SLASH, Opm::InputError::IGNORE); + parseContext.update(Opm::ParseContext::PARSE_MISSING_DIMS_KEYWORD, Opm::InputError::WARN); + parseContext.update(Opm::ParseContext::SUMMARY_UNKNOWN_WELL, Opm::InputError::WARN); + parseContext.update(Opm::ParseContext::SUMMARY_UNKNOWN_GROUP, Opm::InputError::WARN); + } - std::shared_ptr eclipseState = std::make_shared< Opm::EclipseState > ( *deck, parseContext, errorGuard ); - Opm::flowEbosSetDeck(*deck, *eclipseState); - return Opm::flowEbosMain(argc, argv, outputCout, outputFiles); + Opm::FlowMainEbos::printPRTHeader(outputCout); + + deck.reset( new Opm::Deck( parser.parseFile(deckFilename , parseContext, errorGuard))); + Opm::MissingFeatures::checkKeywords(*deck, parseContext, errorGuard); + if ( outputCout ) + Opm::checkDeck(*deck, parser, parseContext, errorGuard); + + eclipseState.reset( new Opm::EclipseState(*deck, parseContext, errorGuard )); + schedule.reset(new Opm::Schedule(*deck, *eclipseState, parseContext, errorGuard)); + summaryConfig.reset( new Opm::SummaryConfig(*deck, *schedule, eclipseState->getTableManager(), parseContext, errorGuard)); + setupMessageLimiter(schedule->getMessageLimits(), "STDOUT_LOGGER"); + + Opm::checkConsistentArrayDimensions(*eclipseState, *schedule, parseContext, errorGuard); + + if (errorGuard) { + errorGuard.dump(); + errorGuard.clear(); + + throw std::runtime_error("Unrecoverable errors were encountered while loading input."); + } + } + bool outputFiles = (outputMode != FileOutputMode::OUTPUT_NONE); + Opm::flowEbosSetDeck(*deck, *eclipseState); + return Opm::flowEbosMain(argc, argv, outputCout, outputFiles); } catch (const std::invalid_argument& e) { From 10e2200a25d27fad616d745126ab0ff3e59e228a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 10 Oct 2019 16:14:44 +0200 Subject: [PATCH 5/9] Add two new flow variants for one-phase (water, water+energy) simulation. --- CMakeLists.txt | 16 +++++++ flow/flow_onephase.cpp | 84 ++++++++++++++++++++++++++++++++++ flow/flow_onephase_energy.cpp | 86 +++++++++++++++++++++++++++++++++++ 3 files changed, 186 insertions(+) create mode 100644 flow/flow_onephase.cpp create mode 100644 flow/flow_onephase_energy.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 0d83dce1a..f93137085 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -163,6 +163,22 @@ opm_add_test(flow_blackoil_dunecpr DEPENDS "opmsimulators" LIBRARIES "opmsimulators") +opm_add_test(flow_onephase + ONLY_COMPILE + DEFAULT_ENABLE_IF ${FLOW_DEFAULT_ENABLE_IF} + SOURCES flow/flow_onephase.cpp + EXE_NAME flow_onephase + DEPENDS "opmsimulators" + LIBRARIES "opmsimulators") + +opm_add_test(flow_onephase_energy + ONLY_COMPILE + DEFAULT_ENABLE_IF ${FLOW_DEFAULT_ENABLE_IF} + SOURCES flow/flow_onephase_energy.cpp + EXE_NAME flow_onephase_energy + DEPENDS "opmsimulators" + LIBRARIES "opmsimulators") + if (BUILD_FLOW) diff --git a/flow/flow_onephase.cpp b/flow/flow_onephase.cpp new file mode 100644 index 000000000..4c3497d0b --- /dev/null +++ b/flow/flow_onephase.cpp @@ -0,0 +1,84 @@ +/* + Copyright 2013, 2014, 2015 SINTEF ICT, Applied Mathematics. + Copyright 2014 Dr. Blatt - HPC-Simulation-Software & Services + Copyright 2015, 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 . +*/ +#include "config.h" +#include "flow/flow_tag.hpp" +#include + + +BEGIN_PROPERTIES + NEW_TYPE_TAG(EclFlowProblemSimple, INHERITS_FROM(EclFlowProblem)); + NEW_PROP_TAG(FluidState); + NEW_PROP_TAG(FluidSystem); + //! The indices required by the model + SET_PROP(EclFlowProblemSimple, Indices) + { + private: + // it is unfortunately not possible to simply use 'TypeTag' here because this leads + // to cyclic definitions of some properties. if this happens the compiler error + // messages unfortunately are *really* confusing and not really helpful. + typedef TTAG(EclFlowProblem) BaseTypeTag; + typedef typename GET_PROP_TYPE(BaseTypeTag, FluidSystem) FluidSystem; + + public: + typedef Ewoms::BlackOilOnePhaseIndices type; + + }; + SET_PROP(EclFlowProblemSimple, FluidState) + { + private: + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; + enum { enableTemperature = GET_PROP_VALUE(TypeTag, EnableTemperature) }; + enum { enableSolvent = GET_PROP_VALUE(TypeTag, EnableSolvent) }; + enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) }; + enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) }; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; + static const bool compositionSwitchEnabled = Indices::gasEnabled; + + public: + //typedef Opm::BlackOilFluidSystemSimple type; + typedef Opm::BlackOilFluidState type; + }; + + // SET_PROP(EclFlowProblemSimple, FluidSystem) + // { + // private: + // //typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + // typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + // typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; + // typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; + + // public: + // typedef Opm::BlackOilFluidSystem type; + // }; +END_PROPERTIES + +int main(int argc, char** argv) +{ + typedef TTAG(EclFlowProblemSimple) TypeTag; + return mainFlow(argc, argv); +} diff --git a/flow/flow_onephase_energy.cpp b/flow/flow_onephase_energy.cpp new file mode 100644 index 000000000..7b91a9939 --- /dev/null +++ b/flow/flow_onephase_energy.cpp @@ -0,0 +1,86 @@ +/* + Copyright 2013, 2014, 2015 SINTEF ICT, Applied Mathematics. + Copyright 2014 Dr. Blatt - HPC-Simulation-Software & Services + Copyright 2015, 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 . +*/ +#include "config.h" +#include "flow/flow_tag.hpp" +#include + + +BEGIN_PROPERTIES + NEW_TYPE_TAG(EclFlowProblemSimple, INHERITS_FROM(EclFlowProblem)); + SET_BOOL_PROP(EclFlowProblemSimple, EnableEnergy, true); + NEW_PROP_TAG(FluidState); + NEW_PROP_TAG(FluidSystem); + //! The indices required by the model + SET_PROP(EclFlowProblemSimple, Indices) + { + private: + // it is unfortunately not possible to simply use 'TypeTag' here because this leads + // to cyclic definitions of some properties. if this happens the compiler error + // messages unfortunately are *really* confusing and not really helpful. + typedef TTAG(EclFlowProblem) BaseTypeTag; + typedef typename GET_PROP_TYPE(BaseTypeTag, FluidSystem) FluidSystem; + + public: + typedef Ewoms::BlackOilOnePhaseIndices type; + + }; + SET_PROP(EclFlowProblemSimple, FluidState) + { + private: + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; + enum { enableTemperature = GET_PROP_VALUE(TypeTag, EnableTemperature) }; + enum { enableSolvent = GET_PROP_VALUE(TypeTag, EnableSolvent) }; + enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) }; + enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) }; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; + static const bool compositionSwitchEnabled = Indices::gasEnabled; + + public: + //typedef Opm::BlackOilFluidSystemSimple type; + typedef Opm::BlackOilFluidState type; + }; + + // SET_PROP(EclFlowProblemSimple, FluidSystem) + // { + // private: + // //typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + // typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + // typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; + // typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; + + // public: + // typedef Opm::BlackOilFluidSystem type; + // }; +END_PROPERTIES + + +int main(int argc, char** argv) +{ + typedef TTAG(EclFlowProblemSimple) TypeTag; + return mainFlow(argc, argv); +} From 71c5129d3d90353c5ada7eb2d2d309ef476468ea Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 11 Oct 2019 08:46:17 +0200 Subject: [PATCH 6/9] Do not build the extra Flow variants by default. --- CMakeLists.txt | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index f93137085..3fbd46f2b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -19,6 +19,7 @@ cmake_minimum_required (VERSION 2.8) option(SIBLING_SEARCH "Search for other modules in sibling directories?" ON) set( USE_OPENMP_DEFAULT OFF ) # Use of OpenMP is considered experimental option(BUILD_FLOW "Build the production oriented flow simulator?" ON) +option(BUILD_FLOW_VARIANTS "Build the variants for flow by default?" OFF) option(BUILD_EBOS "Build the research oriented ebos simulator?" ON) option(BUILD_EBOS_EXTENSIONS "Build the variants for various extensions of ebos by default?" OFF) option(BUILD_EBOS_DEBUG_EXTENSIONS "Build the ebos variants which are purely for debugging by default?" OFF) @@ -155,9 +156,16 @@ opm_add_test(flow flow/flow_ebos_oilwater_polymer.cpp flow/flow_ebos_oilwater_polymer_injectivity.cpp) +if (NOT BUILD_FLOW_VARIANTS) + set(FLOW_VARIANTS_DEFAULT_ENABLE_IF "FALSE") +else() + set(FLOW_VARIANTS_DEFAULT_ENABLE_IF "TRUE") +endif() + +# Variant versions of Flow. opm_add_test(flow_blackoil_dunecpr ONLY_COMPILE - DEFAULT_ENABLE_IF ${FLOW_DEFAULT_ENABLE_IF} + DEFAULT_ENABLE_IF ${FLOW_VARIANTS_DEFAULT_ENABLE_IF} SOURCES flow/flow_blackoil_dunecpr.cpp EXE_NAME flow_blackoil_dunecpr DEPENDS "opmsimulators" @@ -165,7 +173,7 @@ opm_add_test(flow_blackoil_dunecpr opm_add_test(flow_onephase ONLY_COMPILE - DEFAULT_ENABLE_IF ${FLOW_DEFAULT_ENABLE_IF} + DEFAULT_ENABLE_IF ${FLOW_VARIANTS_DEFAULT_ENABLE_IF} SOURCES flow/flow_onephase.cpp EXE_NAME flow_onephase DEPENDS "opmsimulators" @@ -173,14 +181,12 @@ opm_add_test(flow_onephase opm_add_test(flow_onephase_energy ONLY_COMPILE - DEFAULT_ENABLE_IF ${FLOW_DEFAULT_ENABLE_IF} + DEFAULT_ENABLE_IF ${FLOW_VARIANTS_DEFAULT_ENABLE_IF} SOURCES flow/flow_onephase_energy.cpp EXE_NAME flow_onephase_energy DEPENDS "opmsimulators" LIBRARIES "opmsimulators") - - if (BUILD_FLOW) install(TARGETS flow DESTINATION bin) opm_add_bash_completion(flow) From 0b9b20695c8430650ad060ddceb733ef447ed343 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 11 Oct 2019 15:29:49 +0200 Subject: [PATCH 7/9] Change back relative change computation to old order. Changing this caused tiny timestepping differences. --- opm/simulators/flow/BlackoilModelEbos.hpp | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/opm/simulators/flow/BlackoilModelEbos.hpp b/opm/simulators/flow/BlackoilModelEbos.hpp index edcb40672..0b06b58c3 100644 --- a/opm/simulators/flow/BlackoilModelEbos.hpp +++ b/opm/simulators/flow/BlackoilModelEbos.hpp @@ -422,6 +422,12 @@ namespace Opm { Scalar saturationsOld[FluidSystem::numPhases] = { 0.0 }; Scalar oilSaturationOld = 1.0; + + // NB fix me! adding pressures changes to satutation changes does not make sense + Scalar tmp = pressureNew - pressureOld; + resultDelta += tmp*tmp; + resultDenom += pressureNew*pressureNew; + if (FluidSystem::numActivePhases() > 1) { if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { saturationsOld[FluidSystem::waterPhaseIdx] = priVarsOld[Indices::waterSaturationIdx]; @@ -446,10 +452,6 @@ namespace Opm { assert(std::isfinite(resultDenom)); } } - // NB fix me! adding pressures changes to satutation changes does not make sense - Scalar tmp = pressureNew - pressureOld; - resultDelta += tmp*tmp; - resultDenom += pressureNew*pressureNew; } resultDelta = gridView.comm().sum(resultDelta); From 5211217c946645b38aced3f5385f70e114ef43bb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 11 Oct 2019 15:57:51 +0200 Subject: [PATCH 8/9] Whitespace fixes (tabs->spaces, reformatted new files). --- flow/flow_blackoil_dunecpr.cpp | 8 +- flow/flow_onephase.cpp | 113 +++++++++-------- flow/flow_onephase_energy.cpp | 115 ++++++++++-------- flow/flow_tag.hpp | 14 +-- opm/simulators/flow/BlackoilModelEbos.hpp | 48 ++++---- .../timestepping/AdaptiveSimulatorTimer.cpp | 6 +- .../timestepping/AdaptiveTimeSteppingEbos.hpp | 6 +- .../timestepping/TimeStepControl.cpp | 10 +- .../wells/BlackoilWellModel_impl.hpp | 24 ++-- opm/simulators/wells/RateConverter.hpp | 16 +-- opm/simulators/wells/StandardWell_impl.hpp | 8 +- opm/simulators/wells/WellInterface.hpp | 10 +- 12 files changed, 196 insertions(+), 182 deletions(-) diff --git a/flow/flow_blackoil_dunecpr.cpp b/flow/flow_blackoil_dunecpr.cpp index ffa90f28f..f8104be48 100644 --- a/flow/flow_blackoil_dunecpr.cpp +++ b/flow/flow_blackoil_dunecpr.cpp @@ -43,7 +43,7 @@ SET_PROP(EclFlowProblemSimple, FluidState) typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; static const bool compositionSwitchEnabled = Indices::gasEnabled; - + public: //typedef Opm::BlackOilFluidSystemSimple type; typedef Opm::BlackOilFluidState type; @@ -89,7 +89,7 @@ namespace Opm { #endif SET_BOOL_PROP(EclFlowProblemSimple, EnableStorageCache, true); SET_BOOL_PROP(EclFlowProblemSimple, EnableIntensiveQuantityCache, true); - + //SET_INT_PROP(EclFlowProblemSimple, NumWellAdjoint, 1); //SET_BOOL_PROP(EclFlowProblem, EnableStorageCache, true); //SET_BOOL_PROP(EclFlowProblem, EnableIntensiveQuantityCache, true); @@ -98,6 +98,6 @@ namespace Opm { int main(int argc, char** argv) { - typedef TTAG(EclFlowProblemSimple) TypeTag; - return mainFlow(argc, argv); + typedef TTAG(EclFlowProblemSimple) TypeTag; + return mainFlow(argc, argv); } diff --git a/flow/flow_onephase.cpp b/flow/flow_onephase.cpp index 4c3497d0b..b0521a8f8 100644 --- a/flow/flow_onephase.cpp +++ b/flow/flow_onephase.cpp @@ -24,61 +24,68 @@ BEGIN_PROPERTIES - NEW_TYPE_TAG(EclFlowProblemSimple, INHERITS_FROM(EclFlowProblem)); - NEW_PROP_TAG(FluidState); - NEW_PROP_TAG(FluidSystem); - //! The indices required by the model - SET_PROP(EclFlowProblemSimple, Indices) - { - private: - // it is unfortunately not possible to simply use 'TypeTag' here because this leads - // to cyclic definitions of some properties. if this happens the compiler error - // messages unfortunately are *really* confusing and not really helpful. - typedef TTAG(EclFlowProblem) BaseTypeTag; - typedef typename GET_PROP_TYPE(BaseTypeTag, FluidSystem) FluidSystem; - - public: - typedef Ewoms::BlackOilOnePhaseIndices type; - - }; - SET_PROP(EclFlowProblemSimple, FluidState) - { - private: - typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; - typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; - enum { enableTemperature = GET_PROP_VALUE(TypeTag, EnableTemperature) }; - enum { enableSolvent = GET_PROP_VALUE(TypeTag, EnableSolvent) }; - enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) }; - enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) }; - typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; - typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; - static const bool compositionSwitchEnabled = Indices::gasEnabled; - - public: - //typedef Opm::BlackOilFluidSystemSimple type; - typedef Opm::BlackOilFluidState type; - }; +NEW_TYPE_TAG(EclFlowProblemSimple, INHERITS_FROM(EclFlowProblem)); +NEW_PROP_TAG(FluidState); +NEW_PROP_TAG(FluidSystem); +//! The indices required by the model +SET_PROP(EclFlowProblemSimple, Indices) +{ +private: + // it is unfortunately not possible to simply use 'TypeTag' here because this leads + // to cyclic definitions of some properties. if this happens the compiler error + // messages unfortunately are *really* confusing and not really helpful. + typedef TTAG(EclFlowProblem) BaseTypeTag; + typedef typename GET_PROP_TYPE(BaseTypeTag, FluidSystem) FluidSystem; - // SET_PROP(EclFlowProblemSimple, FluidSystem) - // { - // private: - // //typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; - // typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; - // typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; - // typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; - - // public: - // typedef Opm::BlackOilFluidSystem type; - // }; +public: + typedef Ewoms::BlackOilOnePhaseIndices + type; +}; +SET_PROP(EclFlowProblemSimple, FluidState) +{ +private: + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; + enum { enableTemperature = GET_PROP_VALUE(TypeTag, EnableTemperature) }; + enum { enableSolvent = GET_PROP_VALUE(TypeTag, EnableSolvent) }; + enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) }; + enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) }; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; + static const bool compositionSwitchEnabled = Indices::gasEnabled; + +public: + // typedef Opm::BlackOilFluidSystemSimple type; + typedef Opm::BlackOilFluidState + type; +}; + +// SET_PROP(EclFlowProblemSimple, FluidSystem) +// { +// private: +// //typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; +// typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; +// typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; +// typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; + +// public: +// typedef Opm::BlackOilFluidSystem type; +// }; END_PROPERTIES -int main(int argc, char** argv) +int +main(int argc, char** argv) { - typedef TTAG(EclFlowProblemSimple) TypeTag; - return mainFlow(argc, argv); + typedef TTAG(EclFlowProblemSimple) TypeTag; + return mainFlow(argc, argv); } diff --git a/flow/flow_onephase_energy.cpp b/flow/flow_onephase_energy.cpp index 7b91a9939..de8b8cbb0 100644 --- a/flow/flow_onephase_energy.cpp +++ b/flow/flow_onephase_energy.cpp @@ -24,63 +24,70 @@ BEGIN_PROPERTIES - NEW_TYPE_TAG(EclFlowProblemSimple, INHERITS_FROM(EclFlowProblem)); - SET_BOOL_PROP(EclFlowProblemSimple, EnableEnergy, true); - NEW_PROP_TAG(FluidState); - NEW_PROP_TAG(FluidSystem); - //! The indices required by the model - SET_PROP(EclFlowProblemSimple, Indices) - { - private: - // it is unfortunately not possible to simply use 'TypeTag' here because this leads - // to cyclic definitions of some properties. if this happens the compiler error - // messages unfortunately are *really* confusing and not really helpful. - typedef TTAG(EclFlowProblem) BaseTypeTag; - typedef typename GET_PROP_TYPE(BaseTypeTag, FluidSystem) FluidSystem; - - public: - typedef Ewoms::BlackOilOnePhaseIndices type; - - }; - SET_PROP(EclFlowProblemSimple, FluidState) - { - private: - typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; - typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; - enum { enableTemperature = GET_PROP_VALUE(TypeTag, EnableTemperature) }; - enum { enableSolvent = GET_PROP_VALUE(TypeTag, EnableSolvent) }; - enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) }; - enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) }; - typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; - typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; - static const bool compositionSwitchEnabled = Indices::gasEnabled; - - public: - //typedef Opm::BlackOilFluidSystemSimple type; - typedef Opm::BlackOilFluidState type; - }; +NEW_TYPE_TAG(EclFlowProblemSimple, INHERITS_FROM(EclFlowProblem)); +SET_BOOL_PROP(EclFlowProblemSimple, EnableEnergy, true); +NEW_PROP_TAG(FluidState); +NEW_PROP_TAG(FluidSystem); +//! The indices required by the model +SET_PROP(EclFlowProblemSimple, Indices) +{ +private: + // it is unfortunately not possible to simply use 'TypeTag' here because this leads + // to cyclic definitions of some properties. if this happens the compiler error + // messages unfortunately are *really* confusing and not really helpful. + typedef TTAG(EclFlowProblem) BaseTypeTag; + typedef typename GET_PROP_TYPE(BaseTypeTag, FluidSystem) FluidSystem; - // SET_PROP(EclFlowProblemSimple, FluidSystem) - // { - // private: - // //typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; - // typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; - // typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; - // typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; - - // public: - // typedef Opm::BlackOilFluidSystem type; - // }; +public: + typedef Ewoms::BlackOilOnePhaseIndices + type; +}; +SET_PROP(EclFlowProblemSimple, FluidState) +{ +private: + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; + enum { enableTemperature = GET_PROP_VALUE(TypeTag, EnableTemperature) }; + enum { enableSolvent = GET_PROP_VALUE(TypeTag, EnableSolvent) }; + enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) }; + enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) }; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; + static const bool compositionSwitchEnabled = Indices::gasEnabled; + +public: + // typedef Opm::BlackOilFluidSystemSimple type; + typedef Opm::BlackOilFluidState + type; +}; + +// SET_PROP(EclFlowProblemSimple, FluidSystem) +// { +// private: +// //typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; +// typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; +// typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; +// typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; + +// public: +// typedef Opm::BlackOilFluidSystem type; +// }; END_PROPERTIES -int main(int argc, char** argv) +int +main(int argc, char** argv) { - typedef TTAG(EclFlowProblemSimple) TypeTag; - return mainFlow(argc, argv); + typedef TTAG(EclFlowProblemSimple) TypeTag; + return mainFlow(argc, argv); } diff --git a/flow/flow_tag.hpp b/flow/flow_tag.hpp index a568eaad6..6b238a29a 100644 --- a/flow/flow_tag.hpp +++ b/flow/flow_tag.hpp @@ -40,7 +40,7 @@ #include //#include -//#include +//#include #include #include //#include @@ -71,7 +71,7 @@ namespace Opm { Vanguard::setExternalDeck(&deck); Vanguard::setExternalEclState(&eclState); } - + // ----------------- Main program ----------------- template int flowEbosMain(int argc, char** argv, bool outputCout, bool outputFiles) @@ -79,7 +79,7 @@ namespace Opm { // we always want to use the default locale, and thus spare us the trouble // with incorrect locale settings. Opm::resetLocale(); - + #if HAVE_DUNE_FEM Dune::Fem::MPIManager::initialize(argc, argv); #else @@ -306,7 +306,7 @@ int mainFlow(int argc, char** argv) typedef GET_PROP_TYPE(PreTypeTag, Problem) PreProblem; PreProblem::setBriefDescription("Flow, an advanced reservoir simulator for ECL-decks provided by the Open Porous Media project."); - + int status = Opm::FlowMainEbos::setupParameters_(argc, argv); if (status != 0) { // if setupParameters_ returns a value smaller than 0, there was no error, but @@ -391,12 +391,12 @@ int mainFlow(int argc, char** argv) catch (const std::invalid_argument& e) { if (outputCout) { - std::cerr << "Failed to create valid EclipseState object." << std::endl; - std::cerr << "Exception caught: " << e.what() << std::endl; + std::cerr << "Failed to create valid EclipseState object." << std::endl; + std::cerr << "Exception caught: " << e.what() << std::endl; } throw; } - + return EXIT_SUCCESS; } #endif diff --git a/opm/simulators/flow/BlackoilModelEbos.hpp b/opm/simulators/flow/BlackoilModelEbos.hpp index 0b06b58c3..0c77a3b91 100644 --- a/opm/simulators/flow/BlackoilModelEbos.hpp +++ b/opm/simulators/flow/BlackoilModelEbos.hpp @@ -423,35 +423,35 @@ namespace Opm { Scalar saturationsOld[FluidSystem::numPhases] = { 0.0 }; Scalar oilSaturationOld = 1.0; - // NB fix me! adding pressures changes to satutation changes does not make sense + // NB fix me! adding pressures changes to satutation changes does not make sense Scalar tmp = pressureNew - pressureOld; resultDelta += tmp*tmp; resultDenom += pressureNew*pressureNew; - if (FluidSystem::numActivePhases() > 1) { - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - saturationsOld[FluidSystem::waterPhaseIdx] = priVarsOld[Indices::waterSaturationIdx]; - oilSaturationOld -= saturationsOld[FluidSystem::waterPhaseIdx]; - } - - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && - priVarsOld.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) - { - saturationsOld[FluidSystem::gasPhaseIdx] = priVarsOld[Indices::compositionSwitchIdx]; - oilSaturationOld -= saturationsOld[FluidSystem::gasPhaseIdx]; - } + if (FluidSystem::numActivePhases() > 1) { + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + saturationsOld[FluidSystem::waterPhaseIdx] = priVarsOld[Indices::waterSaturationIdx]; + oilSaturationOld -= saturationsOld[FluidSystem::waterPhaseIdx]; + } - if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { - saturationsOld[FluidSystem::oilPhaseIdx] = oilSaturationOld; - } - for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++ phaseIdx) { - Scalar tmp = saturationsNew[phaseIdx] - saturationsOld[phaseIdx]; - resultDelta += tmp*tmp; - resultDenom += saturationsNew[phaseIdx]*saturationsNew[phaseIdx]; - assert(std::isfinite(resultDelta)); - assert(std::isfinite(resultDenom)); - } - } + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && + priVarsOld.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) + { + saturationsOld[FluidSystem::gasPhaseIdx] = priVarsOld[Indices::compositionSwitchIdx]; + oilSaturationOld -= saturationsOld[FluidSystem::gasPhaseIdx]; + } + + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { + saturationsOld[FluidSystem::oilPhaseIdx] = oilSaturationOld; + } + for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++ phaseIdx) { + Scalar tmp = saturationsNew[phaseIdx] - saturationsOld[phaseIdx]; + resultDelta += tmp*tmp; + resultDenom += saturationsNew[phaseIdx]*saturationsNew[phaseIdx]; + assert(std::isfinite(resultDelta)); + assert(std::isfinite(resultDenom)); + } + } } resultDelta = gridView.comm().sum(resultDelta); diff --git a/opm/simulators/timestepping/AdaptiveSimulatorTimer.cpp b/opm/simulators/timestepping/AdaptiveSimulatorTimer.cpp index 93c74c79f..26c0037c6 100644 --- a/opm/simulators/timestepping/AdaptiveSimulatorTimer.cpp +++ b/opm/simulators/timestepping/AdaptiveSimulatorTimer.cpp @@ -60,7 +60,7 @@ AdaptiveSimulatorTimer& AdaptiveSimulatorTimer::operator++ () { ++current_step_; current_time_ += dt_; - assert(dt_ > 0); + assert(dt_ > 0); // store used time step sizes steps_.push_back( dt_ ); return *this; @@ -82,7 +82,7 @@ AdaptiveSimulatorTimer& AdaptiveSimulatorTimer::operator++ () if( dt_ > max_time_step_ ) { dt_ = 0.5 * remaining; } - assert(dt_ > 0); + assert(dt_ > 0); return; } @@ -91,7 +91,7 @@ AdaptiveSimulatorTimer& AdaptiveSimulatorTimer::operator++ () if( 1.5 * dt_ > remaining ) { dt_ = 0.5 * remaining; - assert(dt_ > 0); + assert(dt_ > 0); return; } } diff --git a/opm/simulators/timestepping/AdaptiveTimeSteppingEbos.hpp b/opm/simulators/timestepping/AdaptiveTimeSteppingEbos.hpp index 632ede844..69ff7caf3 100644 --- a/opm/simulators/timestepping/AdaptiveTimeSteppingEbos.hpp +++ b/opm/simulators/timestepping/AdaptiveTimeSteppingEbos.hpp @@ -293,10 +293,10 @@ namespace Opm { double dtEstimate = timeStepControl_->computeTimeStepSize(dt, iterations, relativeChange, substepTimer.simulationTimeElapsed()); - assert(dtEstimate > 0); + assert(dtEstimate > 0); // limit the growth of the timestep size by the growth factor dtEstimate = std::min(dtEstimate, double(maxGrowth_ * dt)); - assert(dtEstimate > 0); + assert(dtEstimate > 0); // further restrict time step size growth after convergence problems if (restarts > 0) { dtEstimate = std::min(growthFactor_ * dt, dtEstimate); @@ -310,7 +310,7 @@ namespace Opm { const double maxPredictionTHPTimestep = 16.0 * unit::day; dtEstimate = std::min(dtEstimate, maxPredictionTHPTimestep); } - assert(dtEstimate > 0); + assert(dtEstimate > 0); if (timestepVerbose_) { std::ostringstream ss; substepReport.reportStep(ss); diff --git a/opm/simulators/timestepping/TimeStepControl.cpp b/opm/simulators/timestepping/TimeStepControl.cpp index c900603e2..9ff748c3b 100644 --- a/opm/simulators/timestepping/TimeStepControl.cpp +++ b/opm/simulators/timestepping/TimeStepControl.cpp @@ -134,11 +134,11 @@ namespace Opm // store new error const double error = relChange.relativeChange(); errors_[ 2 ] = error; - for( int i=0; i<2; ++i ) { - assert(std::isfinite(errors_[i])); - assert(errors_[i]>0); - } - + for( int i=0; i<2; ++i ) { + assert(std::isfinite(errors_[i])); + assert(errors_[i]>0); + } + if( error > tol_ ) { // adjust dt by given tolerance diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index c17fba4b2..90800c810 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -251,18 +251,18 @@ namespace Opm { const unsigned cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0); const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0); const auto& fs = intQuants.fluidState(); - // copy of get perfpressure in Standard well - // exept for value - double perf_pressure = 0.0; - if (Indices::oilEnabled) { - perf_pressure = fs.pressure(FluidSystem::oilPhaseIdx).value(); - } else { - if (Indices::waterEnabled) { - perf_pressure = fs.pressure(FluidSystem::waterPhaseIdx).value(); - } else { - perf_pressure = fs.pressure(FluidSystem::gasPhaseIdx).value(); - } - } + // copy of get perfpressure in Standard well + // exept for value + double perf_pressure = 0.0; + if (Indices::oilEnabled) { + perf_pressure = fs.pressure(FluidSystem::oilPhaseIdx).value(); + } else { + if (Indices::waterEnabled) { + perf_pressure = fs.pressure(FluidSystem::waterPhaseIdx).value(); + } else { + perf_pressure = fs.pressure(FluidSystem::gasPhaseIdx).value(); + } + } cellPressures[cellIdx] = perf_pressure; } well_state_.init(wells(), cellPressures, schedule(), wells_ecl_, timeStepIdx, &previous_well_state_, phase_usage_); diff --git a/opm/simulators/wells/RateConverter.hpp b/opm/simulators/wells/RateConverter.hpp index 1d58729c7..38595d0ad 100644 --- a/opm/simulators/wells/RateConverter.hpp +++ b/opm/simulators/wells/RateConverter.hpp @@ -493,13 +493,13 @@ namespace Opm { // sum p, rs, rv, and T. double hydrocarbonPV = pv_cell*hydrocarbon; - if (hydrocarbonPV > 0) { - pv += hydrocarbonPV; - p += fs.pressure(FluidSystem::oilPhaseIdx).value()*hydrocarbonPV; - rs += fs.Rs().value()*hydrocarbonPV; - rv += fs.Rv().value()*hydrocarbonPV; - T += fs.temperature(FluidSystem::oilPhaseIdx).value()*hydrocarbonPV; - } + if (hydrocarbonPV > 0) { + pv += hydrocarbonPV; + p += fs.pressure(FluidSystem::oilPhaseIdx).value()*hydrocarbonPV; + rs += fs.Rs().value()*hydrocarbonPV; + rv += fs.Rv().value()*hydrocarbonPV; + T += fs.temperature(FluidSystem::oilPhaseIdx).value()*hydrocarbonPV; + } } for (const auto& reg : rmap_.activeRegions()) { @@ -523,7 +523,7 @@ namespace Opm { } } - /** + /** * Region identifier. * * Integral type. diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index 04bf02303..0465e4727 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -238,10 +238,10 @@ namespace Opm StandardWell:: wellVolumeFraction(const unsigned compIdx) const { - if (FluidSystem::numActivePhases() == 1) { - return EvalWell(numWellEq_ + numEq, 1.0); - } - + if (FluidSystem::numActivePhases() == 1) { + return EvalWell(numWellEq_ + numEq, 1.0); + } + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx)) { return primary_variables_evaluation_[WFrac]; } diff --git a/opm/simulators/wells/WellInterface.hpp b/opm/simulators/wells/WellInterface.hpp index 9e50fa3a8..1dee6e106 100644 --- a/opm/simulators/wells/WellInterface.hpp +++ b/opm/simulators/wells/WellInterface.hpp @@ -109,11 +109,11 @@ namespace Opm SurfaceToReservoirVoidage >; static const bool compositionSwitchEnabled = Indices::gasEnabled; using FluidState = Opm::BlackOilFluidState; + FluidSystem, + has_temperature, + has_energy, + compositionSwitchEnabled, + Indices::numPhases >; /// Constructor WellInterface(const Well2& well, const int time_step, const Wells* wells, const ModelParameters& param, From 8ef1958c18da488deb02c0ac3ba2629258f31494 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 11 Oct 2019 20:18:53 +0200 Subject: [PATCH 9/9] Silence warnings. --- opm/simulators/flow/BlackoilModelEbos.hpp | 4 ++-- opm/simulators/wells/StandardWell_impl.hpp | 12 ++++++------ 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/opm/simulators/flow/BlackoilModelEbos.hpp b/opm/simulators/flow/BlackoilModelEbos.hpp index 0c77a3b91..08e227308 100644 --- a/opm/simulators/flow/BlackoilModelEbos.hpp +++ b/opm/simulators/flow/BlackoilModelEbos.hpp @@ -445,8 +445,8 @@ namespace Opm { saturationsOld[FluidSystem::oilPhaseIdx] = oilSaturationOld; } for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++ phaseIdx) { - Scalar tmp = saturationsNew[phaseIdx] - saturationsOld[phaseIdx]; - resultDelta += tmp*tmp; + Scalar tmpSat = saturationsNew[phaseIdx] - saturationsOld[phaseIdx]; + resultDelta += tmpSat*tmpSat; resultDenom += saturationsNew[phaseIdx]*saturationsNew[phaseIdx]; assert(std::isfinite(resultDelta)); assert(std::isfinite(resultDenom)); diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index 0465e4727..6c557b0d3 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -985,8 +985,8 @@ namespace Opm updateExtraPrimaryVariables(dwells); - for (int i = 0; i < primary_variables_.size(); ++i) { - assert(Opm::isfinite(primary_variables_[i])); + for (double v : primary_variables_) { + assert(Opm::isfinite(v)); } } @@ -1329,8 +1329,8 @@ namespace Opm break; } // end of switch - for (int i = 0; i < primary_variables_.size(); ++i) { - assert(Opm::isfinite(primary_variables_[i])); + for (double v : primary_variables_) { + assert(Opm::isfinite(v)); } } @@ -2623,8 +2623,8 @@ namespace Opm } } - for (int i = 0; i < primary_variables_.size(); ++i) { - assert(Opm::isfinite(primary_variables_[i])); + for (double v : primary_variables_) { + assert(Opm::isfinite(v)); } }