diff --git a/CMakeLists.txt b/CMakeLists.txt index ee87f128b..8af981a26 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -90,6 +90,7 @@ macro (sources_hook) if(DUNE_CORNERPOINT_FOUND OR dune-cornerpoint_FOUND) list (APPEND examples_SOURCES ${PROJECT_SOURCE_DIR}/examples/flow_mpi.cpp + ${PROJECT_SOURCE_DIR}/examples/flow_multisegment_mpi.cpp ) endif() endmacro (sources_hook) diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 0a48770a3..5f96fcc35 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -46,6 +46,7 @@ list (APPEND MAIN_SOURCE_FILES opm/autodiff/VFPProperties.cpp opm/autodiff/VFPProdProperties.cpp opm/autodiff/VFPInjProperties.cpp + opm/autodiff/WellMultiSegment.cpp ) @@ -77,6 +78,7 @@ list (APPEND TEST_DATA_FILES list (APPEND EXAMPLE_SOURCE_FILES examples/find_zero.cpp examples/flow.cpp + examples/flow_multisegment.cpp examples/flow_solvent.cpp examples/sim_2p_incomp_ad.cpp examples/sim_simple.cpp @@ -114,6 +116,8 @@ list (APPEND PUBLIC_HEADER_FILES opm/autodiff/BlackoilSolventModel.hpp opm/autodiff/BlackoilSolventModel_impl.hpp opm/autodiff/BlackoilSolventState.hpp + opm/autodiff/BlackoilMultiSegmentModel.hpp + opm/autodiff/BlackoilMultiSegmentModel_impl.hpp opm/autodiff/fastSparseProduct.hpp opm/autodiff/DuneMatrix.hpp opm/autodiff/ExtractParallelGridInformationToISTL.hpp @@ -139,6 +143,8 @@ list (APPEND PUBLIC_HEADER_FILES opm/autodiff/SimulatorFullyImplicitBlackoil.hpp opm/autodiff/SimulatorFullyImplicitBlackoilSolvent.hpp opm/autodiff/SimulatorFullyImplicitBlackoilSolvent_impl.hpp + opm/autodiff/SimulatorFullyImplicitBlackoilMultiSegment.hpp + opm/autodiff/SimulatorFullyImplicitBlackoilMultiSegment_impl.hpp opm/autodiff/SimulatorIncompTwophaseAd.hpp opm/autodiff/TransportSolverTwophaseAd.hpp opm/autodiff/WellDensitySegmented.hpp @@ -149,5 +155,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/autodiff/VFPHelpers.hpp opm/autodiff/VFPProdProperties.hpp opm/autodiff/VFPInjProperties.hpp + opm/autodiff/WellStateMultiSegment.hpp + opm/autodiff/WellMultiSegment.hpp ) diff --git a/examples/flow_multisegment.cpp b/examples/flow_multisegment.cpp new file mode 100644 index 000000000..30b9247af --- /dev/null +++ b/examples/flow_multisegment.cpp @@ -0,0 +1,463 @@ +/* + Copyright 2013 SINTEF ICT, Applied Mathematics. + Copyright 2014 Dr. Blatt - HPC-Simulation-Software & Services + Copyright 2015 IRIS AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + + +#if HAVE_CONFIG_H +#include "config.h" +#endif // HAVE_CONFIG_H + + +#include + +#include + +#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 3) +#include +#else +#include +#endif + +#if HAVE_DUNE_CORNERPOINT && WANT_DUNE_CORNERPOINTGRID +#define USE_DUNE_CORNERPOINTGRID 1 +#include +#include +#else +#undef USE_DUNE_CORNERPOINTGRID +#endif + +#include + +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include // Note: the GridHelpers must be included before this (to make overloads available). \TODO: Fix. + +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#ifdef _OPENMP +#include +#endif + +#include +#include +#include +#include +#include +#include +#include + +namespace +{ + void warnIfUnusedParams(const Opm::parameter::ParameterGroup& param) + { + if (param.anyUnused()) { + std::cout << "-------------------- Unused parameters: --------------------\n"; + param.displayUsage(); + std::cout << "----------------------------------------------------------------" << std::endl; + } + } +} // anon namespace + + + +// ----------------- Main program ----------------- +int +main(int argc, char** argv) +try +{ + using namespace Opm; +#if USE_DUNE_CORNERPOINTGRID + // Must ensure an instance of the helper is created to initialise MPI. + const Dune::MPIHelper& mpi_helper = Dune::MPIHelper::instance(argc, argv); + const int mpi_rank = mpi_helper.rank(); + const int mpi_size = mpi_helper.size(); +#else + // default values for serial run + const int mpi_rank = 0; + const int mpi_size = 1; +#endif + + // Write parameters used for later reference. (only if rank is zero) + const bool output_cout = ( mpi_rank == 0 ); + + if (output_cout) + { + std::string version = moduleVersionName(); + std::cout << "**********************************************************************\n"; + std::cout << "* *\n"; + std::cout << "* This is Flow (version " << version << ")" + << std::string(26 - version.size(), ' ') << "*\n"; + std::cout << "* *\n"; + std::cout << "* Flow is a simulator for fully implicit three-phase black-oil flow, *\n"; + std::cout << "* and is part of OPM. For more information see: *\n"; + std::cout << "* http://opm-project.org *\n"; + std::cout << "* *\n"; + std::cout << "**********************************************************************\n\n"; + } + +#ifdef _OPENMP + if (!getenv("OMP_NUM_THREADS")) { + //Default to at most 4 threads, regardless of + //number of cores (unless ENV(OMP_NUM_THREADS) is defined) + int num_cores = omp_get_num_procs(); + int num_threads = std::min(4, num_cores); + omp_set_num_threads(num_threads); + } +#pragma omp parallel + if (omp_get_thread_num() == 0){ + //opm_get_num_threads() only works as expected within a parallel region. + std::cout << "OpenMP using " << omp_get_num_threads() << " threads." << std::endl; + } +#endif + + // Read parameters, see if a deck was specified on the command line. + if ( output_cout ) + { + std::cout << "--------------- Reading parameters ---------------" << std::endl; + } + + parameter::ParameterGroup param(argc, argv, false, output_cout); + if( !output_cout ) + { + param.disableOutput(); + } + + if (!param.unhandledArguments().empty()) { + if (param.unhandledArguments().size() != 1) { + std::cerr << "You can only specify a single input deck on the command line.\n"; + return EXIT_FAILURE; + } else { + param.insertParameter("deck_filename", param.unhandledArguments()[0]); + } + } + + // We must have an input deck. Grid and props will be read from that. + if (!param.has("deck_filename")) { + std::cerr << "This program must be run with an input deck.\n" + "Specify the deck filename either\n" + " a) as a command line argument by itself\n" + " b) as a command line parameter with the syntax deck_filename=, or\n" + " c) as a parameter in a parameter file (.param or .xml) passed to the program.\n"; + return EXIT_FAILURE; + } + + // bool check_well_controls = false; + // int max_well_control_iterations = 0; + double gravity[3] = { 0.0 }; + std::string deck_filename = param.get("deck_filename"); + + // Write parameters used for later reference. (only if rank is zero) + bool output = ( mpi_rank == 0 ) && param.getDefault("output", true); + std::string output_dir; + if (output) { + // Create output directory if needed. + output_dir = + param.getDefault("output_dir", std::string("output")); + boost::filesystem::path fpath(output_dir); + try { + create_directories(fpath); + } + catch (...) { + std::cerr << "Creating directories failed: " << fpath << std::endl; + return EXIT_FAILURE; + } + // Write simulation parameters. + param.writeParam(output_dir + "/simulation.param"); + } + + std::string logFile = output_dir + "/LOGFILE.txt"; + Opm::ParserPtr parser(new Opm::Parser()); + { + std::shared_ptr streamLog = std::make_shared(logFile , Opm::Log::DefaultMessageTypes); + std::shared_ptr counterLog = std::make_shared(Opm::Log::DefaultMessageTypes); + + Opm::OpmLog::addBackend( "STREAM" , streamLog ); + Opm::OpmLog::addBackend( "COUNTER" , counterLog ); + } + + Opm::ParseMode parseMode({{ ParseMode::PARSE_RANDOM_SLASH , InputError::IGNORE }}); + Opm::DeckConstPtr deck; + std::shared_ptr eclipseState; + try { + deck = parser->parseFile(deck_filename, parseMode); + Opm::checkDeck(deck); + eclipseState.reset(new Opm::EclipseState(deck , parseMode)); + } + catch (const std::invalid_argument& e) { + std::cerr << "Failed to create valid ECLIPSESTATE object. See logfile: " << logFile << std::endl; + std::cerr << "Exception caught: " << e.what() << std::endl; + return EXIT_FAILURE; + } + + std::vector porv = eclipseState->getDoubleGridProperty("PORV")->getData(); +#if USE_DUNE_CORNERPOINTGRID + // Dune::CpGrid as grid manager + typedef Dune::CpGrid Grid; + // Grid init + Grid grid; + grid.processEclipseFormat(deck, false, false, false, porv); +#else + // UnstructuredGrid as grid manager + typedef UnstructuredGrid Grid; + GridManager gridManager( eclipseState->getEclipseGrid(), porv ); + const Grid& grid = *(gridManager.c_grid()); +#endif + + // Possibly override IOConfig setting (from deck) for how often RESTART files should get written to disk (every N report step) + if (param.has("output_interval")) { + int output_interval = param.get("output_interval"); + IOConfigPtr ioConfig = eclipseState->getIOConfig(); + ioConfig->overrideRestartWriteInterval((size_t)output_interval); + } + + const PhaseUsage pu = Opm::phaseUsageFromDeck(deck); + + std::vector compressedToCartesianIdx; + Opm::createGlobalCellArray(grid, compressedToCartesianIdx); + + typedef BlackoilPropsAdFromDeck::MaterialLawManager MaterialLawManager; + auto materialLawManager = std::make_shared(); + materialLawManager->initFromDeck(deck, eclipseState, compressedToCartesianIdx); + + // Rock and fluid init + BlackoilPropertiesFromDeck props( deck, eclipseState, materialLawManager, + Opm::UgGridHelpers::numCells(grid), + Opm::UgGridHelpers::globalCell(grid), + Opm::UgGridHelpers::cartDims(grid), + param); + + BlackoilPropsAdFromDeck new_props( deck, eclipseState, materialLawManager, grid ); + // check_well_controls = param.getDefault("check_well_controls", false); + // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10); + // Rock compressibility. + + RockCompressibility rock_comp(deck, eclipseState); + + // Gravity. + gravity[2] = deck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity; + + BlackoilState state; + // Init state variables (saturation and pressure). + if (param.has("init_saturation")) { + initStateBasic(Opm::UgGridHelpers::numCells(grid), + Opm::UgGridHelpers::globalCell(grid), + Opm::UgGridHelpers::cartDims(grid), + Opm::UgGridHelpers::numFaces(grid), + Opm::UgGridHelpers::faceCells(grid), + Opm::UgGridHelpers::beginFaceCentroids(grid), + Opm::UgGridHelpers::beginCellCentroids(grid), + Opm::UgGridHelpers::dimensions(grid), + props, param, gravity[2], state); + + initBlackoilSurfvol(Opm::UgGridHelpers::numCells(grid), props, state); + + enum { Oil = BlackoilPhases::Liquid, Gas = BlackoilPhases::Vapour }; + if (pu.phase_used[Oil] && pu.phase_used[Gas]) { + const int numPhases = props.numPhases(); + const int numCells = Opm::UgGridHelpers::numCells(grid); + for (int c = 0; c < numCells; ++c) { + state.gasoilratio()[c] = state.surfacevol()[c*numPhases + pu.phase_pos[Gas]] + / state.surfacevol()[c*numPhases + pu.phase_pos[Oil]]; + } + } + } else if (deck->hasKeyword("EQUIL") && props.numPhases() == 3) { + state.init(Opm::UgGridHelpers::numCells(grid), + Opm::UgGridHelpers::numFaces(grid), + props.numPhases()); + const double grav = param.getDefault("gravity", unit::gravity); + initStateEquil(grid, props, deck, eclipseState, grav, state); + state.faceflux().resize(Opm::UgGridHelpers::numFaces(grid), 0.0); + } else { + initBlackoilStateFromDeck(Opm::UgGridHelpers::numCells(grid), + Opm::UgGridHelpers::globalCell(grid), + Opm::UgGridHelpers::numFaces(grid), + Opm::UgGridHelpers::faceCells(grid), + Opm::UgGridHelpers::beginFaceCentroids(grid), + Opm::UgGridHelpers::beginCellCentroids(grid), + Opm::UgGridHelpers::dimensions(grid), + props, deck, gravity[2], state); + } + + + // The capillary pressure is scaled in new_props to match the scaled capillary pressure in props. + if (deck->hasKeyword("SWATINIT")) { + const int numCells = Opm::UgGridHelpers::numCells(grid); + std::vector cells(numCells); + for (int c = 0; c < numCells; ++c) { cells[c] = c; } + std::vector pc = state.saturation(); + props.capPress(numCells, state.saturation().data(), cells.data(), pc.data(),NULL); + new_props.setSwatInitScaling(state.saturation(),pc); + } + + bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0); + const double *grav = use_gravity ? &gravity[0] : 0; + + const bool use_local_perm = param.getDefault("use_local_perm", true); + + DerivedGeology geoprops(grid, new_props, eclipseState, use_local_perm, grav); + boost::any parallel_information; + + // At this point all properties and state variables are correctly initialized + // If there are more than one processors involved, we now repartition the grid + // and initilialize new properties and states for it. + if( mpi_size > 1 ) + { + Opm::distributeGridAndData( grid, deck, eclipseState, state, new_props, geoprops, materialLawManager, parallel_information, use_local_perm ); + } + + // create output writer after grid is distributed, otherwise the parallel output + // won't work correctly since we need to create a mapping from the distributed to + // the global view + Opm::BlackoilOutputWriter outputWriter(grid, param, eclipseState, pu, new_props.permeability() ); + + // Solver for Newton iterations. + std::unique_ptr fis_solver; + { + const std::string cprSolver = "cpr"; + const std::string interleavedSolver = "interleaved"; + const std::string directSolver = "direct"; + const std::string flowDefaultSolver = interleavedSolver; + + std::shared_ptr simCfg = eclipseState->getSimulationConfig(); + std::string solver_approach = flowDefaultSolver; + + if (param.has("solver_approach")) { + solver_approach = param.get("solver_approach"); + } else { + if (simCfg->useCPR()) { + solver_approach = cprSolver; + } + } + + if (solver_approach == cprSolver) { + fis_solver.reset(new NewtonIterationBlackoilCPR(param, parallel_information)); + } else if (solver_approach == interleavedSolver) { + fis_solver.reset(new NewtonIterationBlackoilInterleaved(param, parallel_information)); + } else if (solver_approach == directSolver) { + fis_solver.reset(new NewtonIterationBlackoilSimple(param, parallel_information)); + } else { + OPM_THROW( std::runtime_error , "Internal error - solver approach " << solver_approach << " not recognized."); + } + + } + + Opm::ScheduleConstPtr schedule = eclipseState->getSchedule(); + Opm::TimeMapConstPtr timeMap(schedule->getTimeMap()); + SimulatorTimer simtimer; + + // initialize variables + simtimer.init(timeMap); + + std::map, double> maxDp; + computeMaxDp(maxDp, deck, eclipseState, grid, state, props, gravity[2]); + std::vector threshold_pressures = thresholdPressures(deck, eclipseState, grid, maxDp); + + SimulatorFullyImplicitBlackoilMultiSegment< Grid > simulator(param, + grid, + geoprops, + new_props, + rock_comp.isActive() ? &rock_comp : 0, + *fis_solver, + grav, + deck->hasKeyword("DISGAS"), + deck->hasKeyword("VAPOIL"), + eclipseState, + outputWriter, + threshold_pressures); + + if (!schedule->initOnly()){ + if( output_cout ) + { + std::cout << "\n\n================ Starting main simulation loop ===============\n" + << std::flush; + } + + SimulatorReport fullReport = simulator.run(simtimer, state); + + if( output_cout ) + { + std::cout << "\n\n================ End of simulation ===============\n\n"; + fullReport.reportFullyImplicit(std::cout); + } + + if (output) { + std::string filename = output_dir + "/walltime.txt"; + std::fstream tot_os(filename.c_str(),std::fstream::trunc | std::fstream::out); + fullReport.reportParam(tot_os); + warnIfUnusedParams(param); + } + } else { + outputWriter.writeInit( simtimer ); + if ( output_cout ) + { + std::cout << "\n\n================ Simulation turned off ===============\n" << std::flush; + } + } +} +catch (const std::exception &e) { + std::cerr << "Program threw an exception: " << e.what() << "\n"; + return EXIT_FAILURE; +} + diff --git a/examples/flow_multisegment_mpi.cpp b/examples/flow_multisegment_mpi.cpp new file mode 100644 index 000000000..b54903ef6 --- /dev/null +++ b/examples/flow_multisegment_mpi.cpp @@ -0,0 +1,2 @@ +#define WANT_DUNE_CORNERPOINTGRID 1 +#include "flow_multisegment.cpp" diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp index 61d32942d..ec18441ac 100644 --- a/opm/autodiff/BlackoilModelBase.hpp +++ b/opm/autodiff/BlackoilModelBase.hpp @@ -273,6 +273,7 @@ namespace Opm { WellOps(const Wells* wells); Eigen::SparseMatrix w2p; // well -> perf (scatter) Eigen::SparseMatrix p2w; // perf -> well (gather) + std::vector well_cells; // the set of perforated cells }; // --------- Data members --------- @@ -344,6 +345,8 @@ namespace Opm { // return wells object const Wells& wells () const { assert( bool(wells_ != 0) ); return *wells_; } + int numWellVars() const; + void makeConstantState(SolutionState& state) const; @@ -388,6 +391,10 @@ namespace Opm { assembleMassBalanceEq(const SolutionState& state); void + extractWellPerfProperties(std::vector& mob_perfcells, + std::vector& b_perfcells) const; + + bool solveWellEq(const std::vector& mob_perfcells, const std::vector& b_perfcells, SolutionState& state, @@ -398,12 +405,12 @@ namespace Opm { const std::vector& mob_perfcells, const std::vector& b_perfcells, V& aliveWells, - std::vector& cq_s); + std::vector& cq_s) const; void updatePerfPhaseRatesAndPressures(const std::vector& cq_s, const SolutionState& state, - WellState& xw); + WellState& xw) const; void addWellFluxEq(const std::vector& cq_s, @@ -534,9 +541,8 @@ namespace Opm { /// maximum of tempV for the phase i. /// \param[out] B_avg An array of size MaxNumPhases where entry i contains the average /// of B for the phase i. - /// \param[out] maxNormWell The maximum of the well equations for each phase. + /// \param[out] maxNormWell The maximum of the well flux equations for each phase. /// \param[in] nc The number of cells of the local grid. - /// \param[in] nw The number of wells on the local grid. /// \return The total pore volume over all cells. double convergenceReduction(const Eigen::Array& B, @@ -546,8 +552,7 @@ namespace Opm { std::vector& maxCoeff, std::vector& B_avg, std::vector& maxNormWell, - int nc, - int nw) const; + int nc) const; double dpMaxRel() const { return param_.dp_max_rel_; } double dsMax() const { return param_.ds_max_; } diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index f53fe75cf..f1e82ff70 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -407,7 +407,8 @@ namespace detail { BlackoilModelBase:: WellOps::WellOps(const Wells* wells) : w2p(), - p2w() + p2w(), + well_cells() { if( wells ) { @@ -432,6 +433,8 @@ namespace detail { w2p.setFromTriplets(scatter.begin(), scatter.end()); p2w.setFromTriplets(gather .begin(), gather .end()); + + well_cells.assign(wells->well_cells, wells->well_cells + wells->well_connpos[wells->number_of_wells]); } } @@ -439,6 +442,19 @@ namespace detail { + template + int + BlackoilModelBase::numWellVars() const + { + // For each well, we have a bhp variable, and one flux per phase. + const int nw = localWellsActive() ? wells().number_of_wells : 0; + return (numPhases() + 1) * nw; + } + + + + + template void BlackoilModelBase::makeConstantState(SolutionState& state) const @@ -496,7 +512,7 @@ namespace detail { // and bhp and Q for the wells vars0.reserve(np + 1); variableReservoirStateInitials(x, vars0); - variableWellStateInitials(xw, vars0); + asImpl().variableWellStateInitials(xw, vars0); return vars0; } @@ -681,7 +697,7 @@ namespace detail { } } // wells - variableStateExtractWellsVars(indices, vars, state); + asImpl().variableStateExtractWellsVars(indices, vars, state); return state; } @@ -778,7 +794,6 @@ namespace detail { // taking std::vector arguments, and not Eigen objects. const int nperf = wells().well_connpos[wells().number_of_wells]; const int nw = wells().number_of_wells; - const std::vector well_cells(wells().well_cells, wells().well_cells + nperf); // Compute the average pressure in each well block const V perf_press = Eigen::Map(xw.perfPress().data(), nperf); @@ -791,6 +806,8 @@ namespace detail { } } + const std::vector& well_cells = wops_.well_cells; + // Use cell values for the temperature as the wells don't knows its temperature yet. const ADB perf_temp = subset(state.temperature, well_cells); @@ -882,12 +899,12 @@ namespace detail { SolutionState state = asImpl().variableState(reservoir_state, well_state); SolutionState state0 = state; asImpl().makeConstantState(state0); - computeWellConnectionPressures(state0, well_state); + asImpl().computeWellConnectionPressures(state0, well_state); } // Possibly switch well controls and updating well state to // get reasonable initial conditions for the wells - updateWellControls(well_state); + asImpl().updateWellControls(well_state); // Create the primary variables. SolutionState state = asImpl().variableState(reservoir_state, well_state); @@ -899,7 +916,7 @@ namespace detail { // Compute initial accumulation contributions // and well connection pressures. asImpl().computeAccum(state0, 0); - computeWellConnectionPressures(state0, well_state); + asImpl().computeWellConnectionPressures(state0, well_state); } // OPM_AD_DISKVAL(state.pressure); @@ -920,30 +937,20 @@ namespace detail { return; } - V aliveWells; - const int np = wells().number_of_phases; - std::vector cq_s(np, ADB::null()); - - const int nw = wells().number_of_wells; - const int nperf = wells().well_connpos[nw]; - const std::vector well_cells(wells().well_cells, wells().well_cells + nperf); - - std::vector mob_perfcells(np, ADB::null()); - std::vector b_perfcells(np, ADB::null()); - for (int phase = 0; phase < np; ++phase) { - mob_perfcells[phase] = subset(rq_[phase].mob, well_cells); - b_perfcells[phase] = subset(rq_[phase].b, well_cells); - } + std::vector mob_perfcells; + std::vector b_perfcells; + asImpl().extractWellPerfProperties(mob_perfcells, b_perfcells); if (param_.solve_welleq_initially_ && initial_assembly) { // solve the well equations as a pre-processing step - solveWellEq(mob_perfcells, b_perfcells, state, well_state); + asImpl().solveWellEq(mob_perfcells, b_perfcells, state, well_state); } - + V aliveWells; + std::vector cq_s; asImpl().computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s); asImpl().updatePerfPhaseRatesAndPressures(cq_s, state, well_state); asImpl().addWellFluxEq(cq_s, state); asImpl().addWellContributionToMassBalanceEq(cq_s, state, well_state); - addWellControlEq(state, well_state, aliveWells); + asImpl().addWellControlEq(state, well_state, aliveWells); } @@ -1056,25 +1063,52 @@ namespace detail { { // Add well contributions to mass balance equations const int nc = Opm::AutoDiffGrid::numCells(grid_); - const int nw = wells().number_of_wells; - const int nperf = wells().well_connpos[nw]; - const int np = wells().number_of_phases; - const std::vector well_cells(wells().well_cells, wells().well_cells + nperf); + const int np = asImpl().numPhases(); for (int phase = 0; phase < np; ++phase) { - residual_.material_balance_eq[phase] -= superset(cq_s[phase], well_cells, nc); + residual_.material_balance_eq[phase] -= superset(cq_s[phase], wops_.well_cells, nc); } } + + + template + void + BlackoilModelBase::extractWellPerfProperties(std::vector& mob_perfcells, + std::vector& b_perfcells) const + { + // If we have wells, extract the mobilities and b-factors for + // the well-perforated cells. + if (!asImpl().localWellsActive()) { + mob_perfcells.clear(); + b_perfcells.clear(); + return; + } else { + const int np = asImpl().numPhases(); + const std::vector& well_cells = wops_.well_cells; + mob_perfcells.resize(np, ADB::null()); + b_perfcells.resize(np, ADB::null()); + for (int phase = 0; phase < np; ++phase) { + mob_perfcells[phase] = subset(rq_[phase].mob, well_cells); + b_perfcells[phase] = subset(rq_[phase].b, well_cells); + } + } + } + + + + + + template void BlackoilModelBase::computeWellFlux(const SolutionState& state, const std::vector& mob_perfcells, const std::vector& b_perfcells, V& aliveWells, - std::vector& cq_s) + std::vector& cq_s) const { if( ! localWellsActive() ) return ; @@ -1083,7 +1117,7 @@ namespace detail { const int nperf = wells().well_connpos[nw]; const Opm::PhaseUsage& pu = fluid_.phaseUsage(); V Tw = Eigen::Map(wells().WI, nperf); - const std::vector well_cells(wells().well_cells, wells().well_cells + nperf); + const std::vector& well_cells = wops_.well_cells; // pressure diffs computed already (once per step, not changing per iteration) const V& cdp = well_perforation_pressure_diffs_; @@ -1201,6 +1235,7 @@ namespace detail { ADB cqt_is = cqt_i/volumeRatio; // connection phase volumerates at standard conditions + cq_s.resize(np, ADB::null()); for (int phase = 0; phase < np; ++phase) { cq_s[phase] = cq_ps[phase] + cmix_s[phase]*cqt_is; } @@ -1221,7 +1256,7 @@ namespace detail { template void BlackoilModelBase::updatePerfPhaseRatesAndPressures(const std::vector& cq_s, const SolutionState& state, - WellState& xw) + WellState& xw) const { // Update the perforation phase rates (used to calculate the pressure drop in the wellbore). const int np = wells().number_of_phases; @@ -1553,7 +1588,7 @@ namespace detail { template - void BlackoilModelBase::solveWellEq(const std::vector& mob_perfcells, + bool BlackoilModelBase::solveWellEq(const std::vector& mob_perfcells, const std::vector& b_perfcells, SolutionState& state, WellState& well_state) @@ -1563,6 +1598,7 @@ namespace detail { std::vector cq_s(np, ADB::null()); std::vector indices = variableWellStateIndices(); SolutionState state0 = state; + WellState well_state0 = well_state; asImpl().makeConstantState(state0); std::vector mob_perfcells_const(np, ADB::null()); @@ -1578,15 +1614,15 @@ namespace detail { // bhp and Q for the wells std::vector vars0; vars0.reserve(2); - variableWellStateInitials(well_state, vars0); + asImpl().variableWellStateInitials(well_state, vars0); std::vector vars = ADB::variables(vars0); SolutionState wellSolutionState = state0; - variableStateExtractWellsVars(indices, vars, wellSolutionState); + asImpl().variableStateExtractWellsVars(indices, vars, wellSolutionState); asImpl().computeWellFlux(wellSolutionState, mob_perfcells_const, b_perfcells_const, aliveWells, cq_s); asImpl().updatePerfPhaseRatesAndPressures(cq_s, wellSolutionState, well_state); asImpl().addWellFluxEq(cq_s, wellSolutionState); - addWellControlEq(wellSolutionState, well_state, aliveWells); + asImpl().addWellControlEq(wellSolutionState, well_state, aliveWells); converged = getWellConvergence(it); if (converged) { @@ -1608,9 +1644,9 @@ namespace detail { const Eigen::SparseLU< Sp > solver(Jn0); ADB::V total_residual_v = total_residual.value(); const Eigen::VectorXd& dx = solver.solve(total_residual_v.matrix()); - assert(dx.size() == (well_state.numWells() * (well_state.numPhases()+1))); - updateWellState(dx.array(), well_state); - updateWellControls(well_state); + assert(dx.size() == total_residual_v.size()); + asImpl().updateWellState(dx.array(), well_state); + asImpl().updateWellControls(well_state); } } while (it < 15); @@ -1640,6 +1676,11 @@ namespace detail { asImpl().computeWellConnectionPressures(state, well_state); } + if (!converged) { + well_state = well_state0; + } + + return converged; } @@ -1931,7 +1972,6 @@ namespace detail { using namespace Opm::AutoDiffGrid; const int np = fluid_.numPhases(); const int nc = numCells(grid_); - const int nw = localWellsActive() ? wells().number_of_wells : 0; const V null; assert(null.size() == 0); const V zero = V::Zero(nc); @@ -1946,7 +1986,7 @@ namespace detail { varstart += dxvar.size(); // Extract well parts np phase rates + bhp - const V dwells = subset(dx, Span((np+1)*nw, 1, varstart)); + const V dwells = subset(dx, Span(asImpl().numWellVars(), 1, varstart)); varstart += dwells.size(); assert(varstart == dx.size()); @@ -2136,7 +2176,7 @@ namespace detail { } - updateWellState(dwells,well_state); + asImpl().updateWellState(dwells,well_state); // Update phase conditions used for property calculations. updatePhaseCondFromPrimalVariable(); @@ -2488,11 +2528,12 @@ namespace detail { std::vector& maxCoeff, std::vector& B_avg, std::vector& maxNormWell, - int nc, - int nw) const + int nc) const { const int np = asImpl().numPhases(); const int nm = asImpl().numMaterials(); + const int nw = residual_.well_flux_eq.size() / np; + assert(nw * np == int(residual_.well_flux_eq.size())); // Do the global reductions #if HAVE_MPI @@ -2573,7 +2614,6 @@ namespace detail { const double tol_wells = param_.tolerance_wells_; const int nc = Opm::AutoDiffGrid::numCells(grid_); - const int nw = localWellsActive() ? wells().number_of_wells : 0; const int np = asImpl().numPhases(); const int nm = asImpl().numMaterials(); assert(int(rq_.size()) == nm); @@ -2598,7 +2638,7 @@ namespace detail { const double pvSum = convergenceReduction(B, tempV, R, R_sum, maxCoeff, B_avg, maxNormWell, - nc, nw); + nc); std::vector CNV(nm); std::vector mass_balance_residual(nm); @@ -2661,6 +2701,7 @@ namespace detail { for (int idx = 0; idx < np; ++idx) { std::cout << " W-FLUX(" << materialName(idx).substr(0, 1) << ")"; } + // std::cout << " WELL-CONT "; std::cout << '\n'; } const std::streamsize oprec = std::cout.precision(3); @@ -2675,6 +2716,7 @@ namespace detail { for (int idx = 0; idx < np; ++idx) { std::cout << std::setw(11) << well_flux_residual[idx]; } + // std::cout << std::setw(11) << residualWell; std::cout << std::endl; std::cout.precision(oprec); std::cout.flags(oflags); @@ -2693,7 +2735,6 @@ namespace detail { const double tol_wells = param_.tolerance_wells_; const int nc = Opm::AutoDiffGrid::numCells(grid_); - const int nw = localWellsActive() ? wells().number_of_wells : 0; const int np = asImpl().numPhases(); const int nm = asImpl().numMaterials(); @@ -2713,7 +2754,7 @@ namespace detail { tempV.col(idx) = R.col(idx).abs()/pv; } - convergenceReduction(B, tempV, R, R_sum, maxCoeff, B_avg, maxNormWell, nc, nw); + convergenceReduction(B, tempV, R, R_sum, maxCoeff, B_avg, maxNormWell, nc); std::vector well_flux_residual(np); bool converged_Well = true; diff --git a/opm/autodiff/BlackoilMultiSegmentModel.hpp b/opm/autodiff/BlackoilMultiSegmentModel.hpp new file mode 100644 index 000000000..535537e91 --- /dev/null +++ b/opm/autodiff/BlackoilMultiSegmentModel.hpp @@ -0,0 +1,312 @@ +/* + Copyright 2013, 2015 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_BLACKOILMULTISEGMENTMODEL_HEADER_INCLUDED +#define OPM_BLACKOILMULTISEGMENTMODEL_HEADER_INCLUDED + +#include +#include +#include +#include +#include + +namespace Opm { + + struct BlackoilMultiSegmentSolutionState : public DefaultBlackoilSolutionState + { + explicit BlackoilMultiSegmentSolutionState(const int np) + : DefaultBlackoilSolutionState(np) + , segp ( ADB::null()) + , segqs ( ADB::null()) + { + } + ADB segp; // the segment pressures + ADB segqs; // the segment phase rate in surface volume + }; + + /// A model implementation for three-phase black oil with support + /// for multi-segment wells. + /// + /// It uses automatic differentiation via the class AutoDiffBlock + /// to simplify assembly of the jacobian matrix. + /// \tparam Grid UnstructuredGrid or CpGrid. + /// \tparam Implementation Provides concrete state types. + template + class BlackoilMultiSegmentModel : public BlackoilModelBase> + { + public: + + typedef BlackoilModelBase > Base; // base class + typedef typename Base::ReservoirState ReservoirState; + typedef typename Base::WellState WellState; + typedef BlackoilMultiSegmentSolutionState SolutionState; + + friend Base; + + // --------- Public methods --------- + + /// Construct the model. It will retain references to the + /// arguments of this functions, and they are expected to + /// remain in scope for the lifetime of the solver. + /// \param[in] param parameters + /// \param[in] grid grid data structure + /// \param[in] fluid fluid properties + /// \param[in] geo rock properties + /// \param[in] rock_comp_props if non-null, rock compressibility properties + /// \param[in] wells well structure + /// \param[in] vfp_properties Vertical flow performance tables + /// \param[in] linsolver linear solver + /// \param[in] eclState eclipse state + /// \param[in] has_disgas turn on dissolved gas + /// \param[in] has_vapoil turn on vaporized oil feature + /// \param[in] terminal_output request output to cout/cerr + /// \param[in] wells_multisegment a vector of multisegment wells + BlackoilMultiSegmentModel(const typename Base::ModelParameters& param, + const Grid& grid , + const BlackoilPropsAdInterface& fluid, + const DerivedGeology& geo , + const RockCompressibility* rock_comp_props, + const Wells* wells, + const NewtonIterationBlackoilInterface& linsolver, + Opm::EclipseStateConstPtr eclState, + const bool has_disgas, + const bool has_vapoil, + const bool terminal_output, + const std::vector& wells_multisegment); + + /// Called once before each time step. + /// \param[in] dt time step size + /// \param[in, out] reservoir_state reservoir state variables + /// \param[in, out] well_state well state variables + void prepareStep(const double dt, + ReservoirState& reservoir_state, + WellState& well_state); + + + /// Assemble the residual and Jacobian of the nonlinear system. + /// \param[in] reservoir_state reservoir state variables + /// \param[in, out] well_state well state variables + /// \param[in] initial_assembly pass true if this is the first call to assemble() in this timestep + void assemble(const ReservoirState& reservoir_state, + WellState& well_state, + const bool initial_assembly); + + using Base::numPhases; + using Base::numMaterials; + using Base::materialName; + + protected: + // --------- Data members --------- + + // For non-segmented wells, it should be the density calculated with AVG or SEG way. + // while usually SEG way by default. + using Base::well_perforation_densities_; //Density of each well perforation + using Base::pvdt_; + using Base::geo_; + using Base::active_; + using Base::rq_; + using Base::fluid_; + using Base::terminal_output_; + using Base::grid_; + using Base::canph_; + using Base::residual_; + using Base::isSg_; + using Base::isRs_; + using Base::isRv_; + using Base::has_disgas_; + using Base::has_vapoil_; + using Base::primalVariable_; + using Base::cells_; + using Base::param_; + using Base::linsolver_; + + // Diff to bhp for each well perforation. only for usual wells. + // For segmented wells, they are zeros. + using Base::well_perforation_pressure_diffs_; + + // Pressure correction due to the different depth of the perforation + // and the cell center of the grid block + // For the non-segmented wells, since the perforation are forced to be + // at the center of the grid cell, it should be ZERO. + // It only applies to the mutli-segmented wells. + V well_perforation_cell_pressure_diffs_; + + // Pressure correction due to the depth differennce between segment depth and perforation depth. + ADB well_segment_perforation_pressure_diffs_; + + // The depth difference between segment nodes and perforations + V well_segment_perforation_depth_diffs_; + + // the average of the fluid densities in the grid block + // which is used to calculate the hydrostatic head correction due to the depth difference of the perforation + // and the cell center of the grid block + V well_perforation_cell_densities_; + + // the density of the fluid mixture in the segments + // which is calculated in an implicit way + ADB well_segment_densities_; + + // the hydrostatic pressure drop between segment nodes + // calculated with the above density of fluid mixtures + // for the top segment, they should always be zero for the moment. + ADB well_segment_pressures_delta_; + + // the surface volume of components in the segments + // the initial value at the beginning of the time step + std::vector segment_comp_surf_volume_initial_; + + // the value within the current iteration. + std::vector segment_comp_surf_volume_current_; + + // the mass flow rate in the segments + ADB segment_mass_flow_rates_; + + // the viscosity of the fluid mixture in the segments + // TODO: it is only used to calculate the Reynolds number as we know + // maybe it is not better just to store the Reynolds number? + ADB segment_viscosities_; + + const std::vector wells_multisegment_; + + std::vector top_well_segments_; + + // segment volume by dt (time step) + // to handle the volume effects of the segment + V segvdt_; + + // Well operations and data needed. + struct MultiSegmentWellOps { + explicit MultiSegmentWellOps(const std::vector& wells_ms); + Eigen::SparseMatrix w2p; // well -> perf (scatter) + Eigen::SparseMatrix p2w; // perf -> well (gather) + Eigen::SparseMatrix w2s; // well -> segment (scatter) + Eigen::SparseMatrix s2w; // segment -> well (gather) + Eigen::SparseMatrix s2p; // segment -> perf (scatter) + Eigen::SparseMatrix p2s; // perf -> segment (gather) + Eigen::SparseMatrix s2s_inlets; // segment -> its inlet segments + Eigen::SparseMatrix s2s_outlet; // segment -> its outlet segment + Eigen::SparseMatrix topseg2w; // top segment -> well + AutoDiffMatrix eliminate_topseg; // change the top segment related to be zero + std::vector well_cells; // the set of perforated cells + V conn_trans_factors; // connection transmissibility factors + }; + + MultiSegmentWellOps wops_ms_; + + + // return wells object + // TODO: remove this wells structure + using Base::wells; + using Base::updatePrimalVariableFromState; + using Base::wellsActive; + using Base::phaseCondition; + using Base::fluidRvSat; + using Base::fluidRsSat; + using Base::fluidDensity; + using Base::updatePhaseCondFromPrimalVariable; + using Base::computeGasPressure; + using Base::dpMaxRel; + using Base::dsMax; + using Base::drMaxRel; + using Base::convergenceReduction; + using Base::maxResidualAllowed; + using Base::variableState; + using Base::asImpl; + + const std::vector& wellsMultiSegment() const { return wells_multisegment_; } + + void updateWellControls(WellState& xw) const; + + + void updateWellState(const V& dwells, + WellState& well_state); + + void + variableWellStateInitials(const WellState& xw, + std::vector& vars0) const; + + void computeWellConnectionPressures(const SolutionState& state, + const WellState& xw); + + bool + solveWellEq(const std::vector& mob_perfcells, + const std::vector& b_perfcells, + SolutionState& state, + WellState& well_state); + + void + computeWellFlux(const SolutionState& state, + const std::vector& mob_perfcells, + const std::vector& b_perfcells, + V& aliveWells, + std::vector& cq_s) const; + + void + updatePerfPhaseRatesAndPressures(const std::vector& cq_s, + const SolutionState& state, + WellState& xw) const; + + void + addWellFluxEq(const std::vector& cq_s, + const SolutionState& state); + + void + addWellControlEq(const SolutionState& state, + const WellState& xw, + const V& aliveWells); + + int numWellVars() const; + + void + makeConstantState(SolutionState& state) const; + + void + variableStateExtractWellsVars(const std::vector& indices, + std::vector& vars, + SolutionState& state) const; + + // Calculate the density of the mixture in the segments + // And the surface volume of the components in the segments by dt + void + computeSegmentFluidProperties(const SolutionState& state); + + void + computeSegmentPressuresDelta(const SolutionState& state); + + + }; + + /// Providing types by template specialisation of ModelTraits for BlackoilMultiSegmentModel. + template + struct ModelTraits< BlackoilMultiSegmentModel > + { + typedef BlackoilState ReservoirState; + typedef WellStateMultiSegment WellState; + typedef BlackoilModelParameters ModelParameters; + typedef BlackoilMultiSegmentSolutionState SolutionState; + }; + + + + +} // namespace Opm + +#include "BlackoilMultiSegmentModel_impl.hpp" + +#endif // OPM_BLACKOILMULTISEGMENTMODEL_HEADER_INCLUDED diff --git a/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp b/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp new file mode 100644 index 000000000..b09a64fe1 --- /dev/null +++ b/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp @@ -0,0 +1,1505 @@ +/* + Copyright 2013, 2015 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_BLACKOIMULTISEGMENTLMODEL_IMPL_HEADER_INCLUDED +#define OPM_BLACKOIMULTISEGMENTLMODEL_IMPL_HEADER_INCLUDED + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +//#include + + +namespace Opm { + + + namespace detail + { + ADB onlyWellDerivs(const ADB& x) + { + V val = x.value(); + const int nb = x.numBlocks(); + if (nb < 2) { + OPM_THROW(std::logic_error, "Called onlyWellDerivs() with argument that has " << nb << " blocks."); + } + std::vector derivs = { x.derivative()[nb - 2], x.derivative()[nb - 1] }; + return ADB::function(std::move(val), std::move(derivs)); + } + } // namespace detail + + + + + template + BlackoilMultiSegmentModel:: + BlackoilMultiSegmentModel(const typename Base::ModelParameters& param, + const Grid& grid , + const BlackoilPropsAdInterface& fluid, + const DerivedGeology& geo , + const RockCompressibility* rock_comp_props, + const Wells* wells_arg, + const NewtonIterationBlackoilInterface& linsolver, + Opm::EclipseStateConstPtr eclState, + const bool has_disgas, + const bool has_vapoil, + const bool terminal_output, + const std::vector& wells_multisegment) + : Base(param, grid, fluid, geo, rock_comp_props, wells_arg, linsolver, + eclState, has_disgas, has_vapoil, terminal_output) + , well_segment_perforation_pressure_diffs_(ADB::null()) + , well_segment_densities_(ADB::null()) + , well_segment_pressures_delta_(ADB::null()) + , segment_comp_surf_volume_initial_(fluid.numPhases()) + , segment_comp_surf_volume_current_(fluid.numPhases(), ADB::null()) + , segment_mass_flow_rates_(ADB::null()) + , segment_viscosities_(ADB::null()) + , wells_multisegment_(wells_multisegment) + , wops_ms_(wells_multisegment) + { + } + + + + + + + template + BlackoilMultiSegmentModel:: + MultiSegmentWellOps::MultiSegmentWellOps(const std::vector& wells_ms) + { + if (wells_ms.empty()) { + return; + } + + // Count the total number of perforations and segments. + const int nw = wells_ms.size(); + int total_nperf = 0; + int total_nseg = 0; + for (int w = 0; w < nw; ++w) { + total_nperf += wells_ms[w]->numberOfPerforations(); + total_nseg += wells_ms[w]->numberOfSegments(); + } + + // Create well_cells and conn_trans_factors. + well_cells.reserve(total_nperf); + conn_trans_factors.resize(total_nperf); + int well_perf_start = 0; + for (int w = 0; w < nw; ++w) { + WellMultiSegmentConstPtr well = wells_ms[w]; + well_cells.insert(well_cells.end(), well->wellCells().begin(), well->wellCells().end()); + const std::vector& perf_trans = well->wellIndex(); + std::copy(perf_trans.begin(), perf_trans.end(), conn_trans_factors.data() + well_perf_start); + well_perf_start += well->numberOfPerforations(); + } + assert(well_perf_start == total_nperf); + assert(int(well_cells.size()) == total_nperf); + + // Create all the operator matrices, + // using the setFromTriplets() method. + s2s_inlets = Eigen::SparseMatrix(total_nseg, total_nseg); + s2s_outlet = Eigen::SparseMatrix(total_nseg, total_nseg); + s2w = Eigen::SparseMatrix(nw, total_nseg); + w2s = Eigen::SparseMatrix(total_nseg, nw); + topseg2w = Eigen::SparseMatrix(nw, total_nseg); + s2p = Eigen::SparseMatrix(total_nperf, total_nseg); + p2s = Eigen::SparseMatrix(total_nseg, total_nperf); + typedef Eigen::Triplet Tri; + std::vector s2s_inlets_vector; + std::vector s2s_outlet_vector; + std::vector s2w_vector; + std::vector w2s_vector; + std::vector topseg2w_vector; + std::vector s2p_vector; + std::vector p2s_vector; + V topseg_zero = V::Ones(total_nseg); + s2s_inlets_vector.reserve(total_nseg); + s2s_outlet_vector.reserve(total_nseg); + s2w_vector.reserve(total_nseg); + w2s_vector.reserve(total_nseg); + topseg2w_vector.reserve(nw); + s2p_vector.reserve(total_nperf); + p2s_vector.reserve(total_nperf); + int seg_start = 0; + int perf_start = 0; + for (int w = 0; w < nw; ++w) { + const int ns = wells_ms[w]->numberOfSegments(); + const int np = wells_ms[w]->numberOfPerforations(); + for (int seg = 0; seg < ns; ++seg) { + const int seg_ind = seg_start + seg; + w2s_vector.push_back(Tri(seg_ind, w, 1.0)); + s2w_vector.push_back(Tri(w, seg_ind, 1.0)); + if (seg == 0) { + topseg2w_vector.push_back(Tri(w, seg_ind, 1.0)); + topseg_zero(seg_ind) = 0.0; + } + int seg_outlet = wells_ms[w]->outletSegment()[seg]; + if (seg_outlet >= 0) { + const int outlet_ind = seg_start + seg_outlet; + s2s_inlets_vector.push_back(Tri(outlet_ind, seg_ind, 1.0)); + s2s_outlet_vector.push_back(Tri(seg_ind, outlet_ind, 1.0)); + } + + const auto& seg_perf = wells_ms[w]->segmentPerforations()[seg]; + // the number of perforations related to this segment + const int npseg = seg_perf.size(); + for (int perf = 0; perf < npseg; ++perf) { + const int perf_ind = perf_start + seg_perf[perf]; + s2p_vector.push_back(Tri(perf_ind, seg_ind, 1.0)); + p2s_vector.push_back(Tri(seg_ind, perf_ind, 1.0)); + } + } + seg_start += ns; + perf_start += np; + } + + s2s_inlets.setFromTriplets(s2s_inlets_vector.begin(), s2s_inlets_vector.end()); + s2s_outlet.setFromTriplets(s2s_outlet_vector.begin(), s2s_outlet_vector.end()); + w2s.setFromTriplets(w2s_vector.begin(), w2s_vector.end()); + s2w.setFromTriplets(s2w_vector.begin(), s2w_vector.end()); + topseg2w.setFromTriplets(topseg2w_vector.begin(), topseg2w_vector.end()); + s2p.setFromTriplets(s2p_vector.begin(), s2p_vector.end()); + p2s.setFromTriplets(p2s_vector.begin(), p2s_vector.end()); + + w2p = Eigen::SparseMatrix(total_nperf, nw); + p2w = Eigen::SparseMatrix(nw, total_nperf); + w2p = s2p * w2s; + p2w = s2w * p2s; + + eliminate_topseg = AutoDiffMatrix(topseg_zero.matrix().asDiagonal()); + } + + + + + + + template + void + BlackoilMultiSegmentModel:: + prepareStep(const double dt, + ReservoirState& reservoir_state, + WellState& well_state) + { + pvdt_ = geo_.poreVolume() / dt; + if (active_[Gas]) { + updatePrimalVariableFromState(reservoir_state); + } + + top_well_segments_ = well_state.topSegmentLoc(); + + const int nw = wellsMultiSegment().size(); + const int nseg_total = well_state.numSegments(); + std::vector segment_volume; + segment_volume.reserve(nseg_total); + for (int w = 0; w < nw; ++w) { + WellMultiSegmentConstPtr well = wellsMultiSegment()[w]; + const std::vector& segment_volume_well = well->segmentVolume(); + segment_volume.insert(segment_volume.end(), segment_volume_well.begin(), segment_volume_well.end()); + } + assert(int(segment_volume.size()) == nseg_total); + segvdt_ = Eigen::Map(segment_volume.data(), nseg_total) / dt; + } + + + + + + template + int + BlackoilMultiSegmentModel::numWellVars() const + { + // For each segment, we have a pressure variable, and one flux per phase. + const int nseg = wops_ms_.p2s.rows(); + return (numPhases() + 1) * nseg; + } + + + + + + template + void + BlackoilMultiSegmentModel::makeConstantState(SolutionState& state) const + { + Base::makeConstantState(state); + state.segp = ADB::constant(state.segp.value()); + state.segqs = ADB::constant(state.segqs.value()); + } + + + + + + template + void + BlackoilMultiSegmentModel::variableWellStateInitials(const WellState& xw, std::vector& vars0) const + { + // Initial well rates + if ( wellsMultiSegment().size() > 0 ) + { + // Need to reshuffle well segment rates, from phase running fastest + const int nseg = xw.numSegments(); + const int np = xw.numPhases(); + + // The transpose() below switches the ordering of the segment rates + const DataBlock segrates = Eigen::Map(& xw.segPhaseRates()[0], nseg, np).transpose(); + // segment phase rates in surface volume + const V segqs = Eigen::Map(segrates.data(), nseg * np); + vars0.push_back(segqs); + + // for the pressure of the segments + const V segp = Eigen::Map(& xw.segPress()[0], xw.segPress().size()); + vars0.push_back(segp); + } + else + { + // push null sates for segqs and segp + vars0.push_back(V()); + vars0.push_back(V()); + } + } + + + + + + template + void + BlackoilMultiSegmentModel::variableStateExtractWellsVars(const std::vector& indices, + std::vector& vars, + SolutionState& state) const + { + // TODO: using the original Qs for the segment rates for now, to be fixed eventually. + // TODO: using the original Bhp for the segment pressures for now, to be fixed eventually. + + // segment phase rates in surface volume + state.segqs = std::move(vars[indices[Qs]]); + + // segment pressures + state.segp = std::move(vars[indices[Bhp]]); + + // The qs and bhp are no longer primary variables, but could + // still be used in computations. They are identical to the + // pressures and flows of the top segments. + const int np = numPhases(); + const int ns = state.segp.size(); + const int nw = top_well_segments_.size(); + state.qs = ADB::constant(ADB::V::Zero(np*nw)); + for (int phase = 0; phase < np; ++phase) { + // Extract segment fluxes for this phase (ns consecutive elements). + ADB segqs_phase = subset(state.segqs, Span(ns, 1, ns*phase)); + // Extract top segment fluxes (= well fluxes) + ADB wellqs_phase = subset(segqs_phase, top_well_segments_); + // Expand to full size of qs (which contains all phases) and add. + state.qs += superset(wellqs_phase, Span(nw, 1, nw*phase), nw*np); + } + state.bhp = subset(state.segp, top_well_segments_); + } + + + + + // TODO: This is just a preliminary version, remains to be improved later when we decide a better way + // TODO: to intergrate the usual wells and multi-segment wells. + template + void BlackoilMultiSegmentModel::computeWellConnectionPressures(const SolutionState& state, + const WellState& xw) + { + if( ! wellsActive() ) return ; + + using namespace Opm::AutoDiffGrid; + // 1. Compute properties required by computeConnectionPressureDelta(). + // Note that some of the complexity of this part is due to the function + // taking std::vector arguments, and not Eigen objects. + const int nperf_total = xw.numPerforations(); + const int nw = xw.numWells(); + + std::vector& well_cells = wops_ms_.well_cells; + + well_perforation_densities_ = V::Zero(nperf_total); + + const V perf_press = Eigen::Map(xw.perfPress().data(), nperf_total); + + V avg_press = perf_press * 0.0; + + // for the non-segmented/regular wells, calculated the average pressures. + // If it is the top perforation, then average with the bhp(). + // If it is not the top perforation, then average with the perforation above it(). + int start_segment = 0; + for (int w = 0; w < nw; ++w) { + const int nseg = wellsMultiSegment()[w]->numberOfSegments(); + if (wellsMultiSegment()[w]->isMultiSegmented()) { + // maybe we should give some reasonable values to prevent the following calculations fail + start_segment += nseg; + continue; + } + + std::string well_name(wellsMultiSegment()[w]->name()); + typedef typename WellStateMultiSegment::SegmentedWellMapType::const_iterator const_iterator; + const_iterator it_well = xw.segmentedWellMap().find(well_name); + assert(it_well != xw.segmentedWellMap().end()); + + const int start_perforation = (*it_well).second.start_perforation; + const int end_perforation = start_perforation + (*it_well).second.number_of_perforations; + for (int perf = start_perforation; perf < end_perforation; ++perf) { + const double p_above = perf == start_perforation ? state.segp.value()[start_segment] : perf_press[perf - 1]; + const double p_avg = (perf_press[perf] + p_above)/2; + avg_press[perf] = p_avg; + } + start_segment += nseg; + } + assert(start_segment == xw.numSegments()); + + // Use cell values for the temperature as the wells don't knows its temperature yet. + const ADB perf_temp = subset(state.temperature, well_cells); + + // Compute b, rsmax, rvmax values for perforations. + // Evaluate the properties using average well block pressures + // and cell values for rs, rv, phase condition and temperature. + const ADB avg_press_ad = ADB::constant(avg_press); + std::vector perf_cond(nperf_total); + const std::vector& pc = phaseCondition(); + for (int perf = 0; perf < nperf_total; ++perf) { + perf_cond[perf] = pc[well_cells[perf]]; + } + const PhaseUsage& pu = fluid_.phaseUsage(); + DataBlock b(nperf_total, pu.num_phases); + std::vector rsmax_perf(nperf_total, 0.0); + std::vector rvmax_perf(nperf_total, 0.0); + if (pu.phase_used[BlackoilPhases::Aqua]) { + const V bw = fluid_.bWat(avg_press_ad, perf_temp, well_cells).value(); + b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw; + } + assert(active_[Oil]); + const V perf_so = subset(state.saturation[pu.phase_pos[Oil]].value(), well_cells); + if (pu.phase_used[BlackoilPhases::Liquid]) { + const ADB perf_rs = subset(state.rs, well_cells); + const V bo = fluid_.bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).value(); + b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo; + const V rssat = fluidRsSat(avg_press, perf_so, well_cells); + rsmax_perf.assign(rssat.data(), rssat.data() + nperf_total); + } + if (pu.phase_used[BlackoilPhases::Vapour]) { + const ADB perf_rv = subset(state.rv, well_cells); + const V bg = fluid_.bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value(); + b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg; + const V rvsat = fluidRvSat(avg_press, perf_so, well_cells); + rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf_total); + } + // b is row major, so can just copy data. + std::vector b_perf(b.data(), b.data() + nperf_total * pu.num_phases); + // Extract well connection depths. + const V depth = cellCentroidsZToEigen(grid_); + const V perfcelldepth = subset(depth, well_cells); + std::vector perf_cell_depth(perfcelldepth.data(), perfcelldepth.data() + nperf_total); + + // Surface density. + // The compute density segment wants the surface densities as + // an np * number of wells cells array + V rho = superset(fluid_.surfaceDensity(0 , well_cells), Span(nperf_total, pu.num_phases, 0), nperf_total * pu.num_phases); + for (int phase = 1; phase < pu.num_phases; ++phase) { + rho += superset(fluid_.surfaceDensity(phase , well_cells), Span(nperf_total, pu.num_phases, phase), nperf_total * pu.num_phases); + } + std::vector surf_dens_perf(rho.data(), rho.data() + nperf_total * pu.num_phases); + + // Gravity + double grav = detail::getGravity(geo_.gravity(), dimensions(grid_)); + + // 2. Compute densities + std::vector cd = + WellDensitySegmented::computeConnectionDensities( + wells(), xw, fluid_.phaseUsage(), + b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); + + // 3. Compute pressure deltas + std::vector cdp = + WellDensitySegmented::computeConnectionPressureDelta( + wells(), perf_cell_depth, cd, grav); + + // 4. Store the results + well_perforation_densities_ = Eigen::Map(cd.data(), nperf_total); // This one is not useful for segmented wells at all + well_perforation_pressure_diffs_ = Eigen::Map(cdp.data(), nperf_total); + + // compute the average of the fluid densites in the well blocks. + // the average is weighted according to the fluid relative permeabilities. + const std::vector kr_adb = Base::computeRelPerm(state); + size_t temp_size = kr_adb.size(); + std::vector perf_kr; + for(size_t i = 0; i < temp_size; ++i) { + // const ADB kr_phase_adb = subset(kr_adb[i], well_cells); + const V kr_phase = (subset(kr_adb[i], well_cells)).value(); + perf_kr.push_back(kr_phase); + } + + + // compute the averaged density for the well block + // TODO: for the non-segmented wells, they should be set to zero + // TODO: for the moment, they are still calculated, while not used later. + for (int i = 0; i < nperf_total; ++i) { + double sum_kr = 0.; + int np = perf_kr.size(); // make sure it is 3 + for (int p = 0; p < np; ++p) { + sum_kr += perf_kr[p][i]; + } + + for (int p = 0; p < np; ++p) { + perf_kr[p][i] /= sum_kr; + } + } + + V rho_avg_perf = V::Constant(nperf_total, 0.0); + // TODO: make sure the order of the density and the order of the kr are the same. + for (int phaseIdx = 0; phaseIdx < fluid_.numPhases(); ++phaseIdx) { + const int canonicalPhaseIdx = canph_[phaseIdx]; + const ADB fluid_density = fluidDensity(canonicalPhaseIdx, rq_[phaseIdx].b, state.rs, state.rv); + const V rho_perf = subset(fluid_density, well_cells).value(); + // TODO: phaseIdx or canonicalPhaseIdx ? + rho_avg_perf += rho_perf * perf_kr[phaseIdx]; + } + + well_perforation_cell_densities_ = Eigen::Map(rho_avg_perf.data(), nperf_total); + + // We should put this in a global class + std::vector perf_depth_vec; + perf_depth_vec.reserve(nperf_total); + for (int w = 0; w < nw; ++w) { + WellMultiSegmentConstPtr well = wellsMultiSegment()[w]; + const std::vector& perf_depth_well = well->perfDepth(); + perf_depth_vec.insert(perf_depth_vec.end(), perf_depth_well.begin(), perf_depth_well.end()); + } + assert(int(perf_depth_vec.size()) == nperf_total); + const V perf_depth = Eigen::Map(perf_depth_vec.data(), nperf_total); + + const V perf_cell_depth_diffs = perf_depth - perfcelldepth; + + well_perforation_cell_pressure_diffs_ = grav * well_perforation_cell_densities_ * perf_cell_depth_diffs; + + + // Calculating the depth difference between segment nodes and perforations. + // TODO: should be put somewhere else for better clarity later + well_segment_perforation_depth_diffs_ = V::Constant(nperf_total, -1e100); + + int start_perforation = 0; + for (int w = 0; w < nw; ++w) { + WellMultiSegmentConstPtr well = wellsMultiSegment()[w]; + const int nseg = well->numberOfSegments(); + const int nperf = well->numberOfPerforations(); + const std::vector>& segment_perforations = well->segmentPerforations(); + for (int s = 0; s < nseg; ++s) { + const int nperf_seg = segment_perforations[s].size(); + const double segment_depth = well->segmentDepth()[s]; + for (int perf = 0; perf < nperf_seg; ++perf) { + const int perf_number = segment_perforations[s][perf] + start_perforation; + well_segment_perforation_depth_diffs_[perf_number] = segment_depth - perf_depth[perf_number]; + } + } + start_perforation += nperf; + } + assert(start_perforation == nperf_total); + } + + + + + + + + template + void + BlackoilMultiSegmentModel:: + assemble(const ReservoirState& reservoir_state, + WellState& well_state, + const bool initial_assembly) + { + using namespace Opm::AutoDiffGrid; + + // TODO: include VFP effect. + // If we have VFP tables, we need the well connection + // pressures for the "simple" hydrostatic correction + // between well depth and vfp table depth. + // if (isVFPActive()) { + // SolutionState state = asImpl().variableState(reservoir_state, well_state); + // SolutionState state0 = state; + // asImpl().makeConstantState(state0); + // asImpl().computeWellConnectionPressures(state0, well_state); + // } + + // Possibly switch well controls and updating well state to + // get reasonable initial conditions for the wells + asImpl().updateWellControls(well_state); + + // Create the primary variables. + SolutionState state = asImpl().variableState(reservoir_state, well_state); + + if (initial_assembly) { + // Create the (constant, derivativeless) initial state. + SolutionState state0 = state; + asImpl().makeConstantState(state0); + // Compute initial accumulation contributions + // and well connection pressures. + asImpl().computeAccum(state0, 0); + asImpl().computeSegmentFluidProperties(state0); + const int np = numPhases(); + assert(np == int(segment_comp_surf_volume_initial_.size())); + for (int phase = 0; phase < np; ++phase) { + segment_comp_surf_volume_initial_[phase] = segment_comp_surf_volume_current_[phase].value(); + } + asImpl().computeWellConnectionPressures(state0, well_state); + } + + // OPM_AD_DISKVAL(state.pressure); + // OPM_AD_DISKVAL(state.saturation[0]); + // OPM_AD_DISKVAL(state.saturation[1]); + // OPM_AD_DISKVAL(state.saturation[2]); + // OPM_AD_DISKVAL(state.rs); + // OPM_AD_DISKVAL(state.rv); + // OPM_AD_DISKVAL(state.qs); + // OPM_AD_DISKVAL(state.bhp); + + // -------- Mass balance equations -------- + asImpl().assembleMassBalanceEq(state); + + // -------- Well equations ---------- + + if ( ! wellsActive() ) { + return; + } + + asImpl().computeSegmentFluidProperties(state); + asImpl().computeSegmentPressuresDelta(state); + + std::vector mob_perfcells; + std::vector b_perfcells; + asImpl().extractWellPerfProperties(mob_perfcells, b_perfcells); + if (param_.solve_welleq_initially_ && initial_assembly) { + // solve the well equations as a pre-processing step + asImpl().solveWellEq(mob_perfcells, b_perfcells, state, well_state); + } + + // the perforation flux here are different + // it is related to the segment location + V aliveWells; + std::vector cq_s; + asImpl().computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s); + asImpl().updatePerfPhaseRatesAndPressures(cq_s, state, well_state); + asImpl().addWellFluxEq(cq_s, state); + asImpl().addWellContributionToMassBalanceEq(cq_s, state, well_state); + asImpl().addWellControlEq(state, well_state, aliveWells); + } + + + + + + template + void + BlackoilMultiSegmentModel::computeWellFlux(const SolutionState& state, + const std::vector& mob_perfcells, + const std::vector& b_perfcells, + V& aliveWells, + std::vector& cq_s) const + { + // if( ! wellsActive() ) return ; + if (wellsMultiSegment().size() == 0) return; + + const int nw = wellsMultiSegment().size(); + const Opm::PhaseUsage& pu = fluid_.phaseUsage(); + + aliveWells = V::Constant(nw, 1.0); + + const int np = numPhases(); + const int nseg = wops_ms_.s2p.cols(); + const int nperf = wops_ms_.s2p.rows(); + + cq_s.resize(np, ADB::null()); + + { + const V& Tw = wops_ms_.conn_trans_factors; + const std::vector& well_cells = wops_ms_.well_cells; + + // determining in-flow (towards well-bore) or out-flow (towards reservoir) + // for mutli-segmented wells and non-segmented wells, the calculation of the drawdown are different. + const ADB& p_perfcells = subset(state.pressure, well_cells); + const ADB& rs_perfcells = subset(state.rs, well_cells); + const ADB& rv_perfcells = subset(state.rv, well_cells); + + const ADB& seg_pressures = state.segp; + + const ADB seg_pressures_perf = wops_ms_.s2p * seg_pressures; + + // Create selector for perforations of multi-segment vs. regular wells. + V is_multisegment_well(nw); + for (int w = 0; w < nw; ++w) { + is_multisegment_well[w] = double(wellsMultiSegment()[w]->isMultiSegmented()); + } + // Take one flag per well and expand to one flag per perforation. + V is_multisegment_perf = wops_ms_.w2p * is_multisegment_well.matrix(); + Selector msperf_selector(is_multisegment_perf, Selector::NotEqualZero); + + // Compute drawdown. + ADB h_nc = msperf_selector.select(well_segment_perforation_pressure_diffs_, + ADB::constant(well_perforation_pressure_diffs_)); + const V h_cj = msperf_selector.select(well_perforation_cell_pressure_diffs_, V::Zero(nperf)); + + // Special handling for when we are called from solveWellEq(). + // TODO: restructure to eliminate need for special treatmemt. + if ((h_nc.numBlocks() != 0) && (h_nc.numBlocks() != seg_pressures_perf.numBlocks())) { + assert(seg_pressures_perf.numBlocks() == 2); + assert(h_nc.numBlocks() > 2); + h_nc = detail::onlyWellDerivs(h_nc); + assert(h_nc.numBlocks() == 2); + } + + ADB drawdown = (p_perfcells + h_cj - seg_pressures_perf - h_nc); + + // selects injection perforations + V selectInjectingPerforations = V::Zero(nperf); + // selects producing perforations + V selectProducingPerforations = V::Zero(nperf); + for (int c = 0; c < nperf; ++c){ + if (drawdown.value()[c] < 0) + selectInjectingPerforations[c] = 1; + else + selectProducingPerforations[c] = 1; + } + + // handling flow into wellbore + // maybe there are something to do there make the procedure easier. + std::vector cq_ps(np, ADB::null()); + for (int phase = 0; phase < np; ++phase) { + const ADB cq_p = -(selectProducingPerforations * Tw) * (mob_perfcells[phase] * drawdown); + cq_ps[phase] = b_perfcells[phase] * cq_p; + } + + if (active_[Oil] && active_[Gas]) { + const int oilpos = pu.phase_pos[Oil]; + const int gaspos = pu.phase_pos[Gas]; + const ADB cq_psOil = cq_ps[oilpos]; + const ADB cq_psGas = cq_ps[gaspos]; + cq_ps[gaspos] += rs_perfcells * cq_psOil; + cq_ps[oilpos] += rv_perfcells * cq_psGas; + } + + // hadling flow out from wellbore + ADB total_mob = mob_perfcells[0]; + for (int phase = 1; phase < np; ++phase) { + total_mob += mob_perfcells[phase]; + } + + // injection perforations total volume rates + const ADB cqt_i = -(selectInjectingPerforations * Tw) * (total_mob * drawdown); + + // compute wellbore mixture for injecting perforations + // The wellbore mixture depends on the inflow from the reservoir + // and the well injection rates. + // TODO: should this based on the segments? + // TODO: for the usual wells, the well rates are the sum of the perforations. + // TODO: for multi-segmented wells, the segment rates are not the sum of the perforations. + + // TODO: two options here + // TODO: 1. for each segment, only the inflow from the perforations related to this segment are considered. + // TODO: 2. for each segment, the inflow from the perforrations related to this segment and also all the inflow + // TODO: from the upstreaming sgments and their perforations need to be considered. + // TODO: This way can be the more consistent way, while let us begin with the first option. The second option + // TODO: involves one operations that are not valid now. (i.e. how to transverse from the leaves to the root, + // TODO: although we can begin from the brutal force way) + + // TODO: stop using wells() here. + const DataBlock compi = Eigen::Map(wells().comp_frac, nw, np); + std::vector wbq(np, ADB::null()); + ADB wbqt = ADB::constant(V::Zero(nseg)); + + for (int phase = 0; phase < np; ++phase) { + const ADB& q_ps = wops_ms_.p2s * cq_ps[phase]; + const ADB& q_s = subset(state.segqs, Span(nseg, 1, phase * nseg)); + Selector injectingPhase_selector(q_s.value(), Selector::GreaterZero); + + const int pos = pu.phase_pos[phase]; + + // this is per segment + wbq[phase] = (wops_ms_.w2s * ADB::constant(compi.col(pos)) * injectingPhase_selector.select(q_s, ADB::constant(V::Zero(nseg)))) - q_ps; + + // TODO: it should be a single value for this certain well. + // TODO: it need to be changed later to handle things more consistently + // or there should be an earsier way to decide if the well is dead. + wbqt += wbq[phase]; + } + + // Set aliveWells. + // the first value of the wbqt is the one to decide if the well is dead + // or there should be some dead segments? + { + int topseg = 0; + for (int w = 0; w < nw; ++w) { + if (wbqt.value()[topseg] == 0.0) { // yes we really mean == here, no fuzzyness + aliveWells[w] = 0.0; + } + topseg += wells_multisegment_[w]->numberOfSegments(); + } + } + + // compute wellbore mixture at standard conditions. + // before, the determination of alive wells is based on wells. + // now, will there be any dead segment? I think no. + // TODO: it is not clear if the cmix_s should be based on segment or the well + std::vector cmix_s(np, ADB::null()); + Selector aliveWells_selector(aliveWells, Selector::NotEqualZero); + for (int phase = 0; phase < np; ++phase) { + const int pos = pu.phase_pos[phase]; + const ADB phase_fraction = wops_ms_.topseg2w * (wbq[phase] / wbqt); + cmix_s[phase] = wops_ms_.w2p * aliveWells_selector.select(phase_fraction, ADB::constant(compi.col(pos))); + } + + // compute volume ration between connection at standard conditions + ADB volumeRatio = ADB::constant(V::Zero(nperf)); + const ADB d = V::Constant(nperf,1.0) - rv_perfcells * rs_perfcells; + + for (int phase = 0; phase < np; ++phase) { + ADB tmp = cmix_s[phase]; + if (phase == Oil && active_[Gas]) { + const int gaspos = pu.phase_pos[Gas]; + tmp = tmp - rv_perfcells * cmix_s[gaspos] / d; + } + if (phase == Gas && active_[Oil]) { + const int oilpos = pu.phase_pos[Oil]; + tmp = tmp - rs_perfcells * cmix_s[oilpos] / d; + } + volumeRatio += tmp / b_perfcells[phase]; + } + + // injecting connections total volumerates at standard conditions + ADB cqt_is = cqt_i/volumeRatio; + + // connection phase volumerates at standard conditions + for (int phase = 0; phase < np; ++phase) { + cq_s[phase] = cq_ps[phase] + cmix_s[phase]*cqt_is; + } + } + } + + + + + + template + void BlackoilMultiSegmentModel::updatePerfPhaseRatesAndPressures(const std::vector& cq_s, + const SolutionState& state, + WellState& xw) const + { + // Update the perforation phase rates (used to calculate the pressure drop in the wellbore). + const int np = numPhases(); + const int nw = wellsMultiSegment().size(); + const int nperf_total = xw.perfPress().size(); + + V cq = superset(cq_s[0].value(), Span(nperf_total, np, 0), nperf_total * np); + for (int phase = 1; phase < np; ++phase) { + cq += superset(cq_s[phase].value(), Span(nperf_total, np, phase), nperf_total * np); + } + xw.perfPhaseRates().assign(cq.data(), cq.data() + nperf_total * np); + + // Update the perforation pressures for usual wells first to recover the resutls + // without mutlti segment wells. For segment wells, it has not been decided if + // we need th concept of preforation pressures + xw.perfPress().resize(nperf_total, -1.e100); + + const V& cdp = well_perforation_pressure_diffs_; + int start_segment = 0; + int start_perforation = 0; + for (int i = 0; i < nw; ++i) { + WellMultiSegmentConstPtr well = wellsMultiSegment()[i]; + const int nperf = well->numberOfPerforations(); + const int nseg = well->numberOfSegments(); + if (well->isMultiSegmented()) { + start_segment += nseg; + start_perforation += nperf; + continue; + } + const V cdp_well = subset(cdp, Span(nperf, 1, start_perforation)); + const ADB segp = subset(state.segp, Span(nseg, 1, start_segment)); + const V perfpressure = (well->wellOps().s2p * segp.value().matrix()).array() + cdp_well; + std::copy(perfpressure.data(), perfpressure.data() + nperf, &xw.perfPress()[start_perforation]); + + start_segment += nseg; + start_perforation += nperf; + } + } + + + + template + void BlackoilMultiSegmentModel::addWellFluxEq(const std::vector& cq_s, + const SolutionState& state) + { + // the well flux equations are for each segment and each phase. + // /delta m_p_n / dt - /sigma Q_pi - /sigma q_pj + Q_pn = 0 + // 1. It is the gain of the amount of the component p in the segment n during the + // current time step under stock-tank conditions. + // It is used to handle the volume storage effects of the wellbore. + // We need the information from the previous step and the crrent time step. + // 2. for the second term, it is flow into the segment from the inlet segments, + // which are unknown and treated implictly. + // 3. for the third term, it is the inflow through the perforations. + // 4. for the last term, it is the outlet rates and also the segment rates, + // which are the primary variable. + const int np = numPhases(); + const int nseg_total = state.segp.size(); + + ADB segqs = state.segqs; + + std::vector segment_volume_change_dt(np, ADB::null()); + for (int phase = 0; phase < np; ++phase) { + // Gain of the surface volume of each component in the segment by dt + segment_volume_change_dt[phase] = segment_comp_surf_volume_current_[phase] - + segment_comp_surf_volume_initial_[phase]; + + // Special handling for when we are called from solveWellEq(). + // TODO: restructure to eliminate need for special treatmemt. + if (segment_volume_change_dt[phase].numBlocks() != segqs.numBlocks()) { + assert(segment_volume_change_dt[phase].numBlocks() > 2); + assert(segqs.numBlocks() == 2); + segment_volume_change_dt[phase] = detail::onlyWellDerivs(segment_volume_change_dt[phase]); + assert(segment_volume_change_dt[phase].numBlocks() == 2); + } + + const ADB cq_s_seg = wops_ms_.p2s * cq_s[phase]; + const ADB segqs_phase = subset(segqs, Span(nseg_total, 1, phase * nseg_total)); + segqs -= superset(cq_s_seg + wops_ms_.s2s_inlets * segqs_phase + segment_volume_change_dt[phase], + Span(nseg_total, 1, phase * nseg_total), np * nseg_total); + } + + residual_.well_flux_eq = segqs; + } + + + + template + void BlackoilMultiSegmentModel::updateWellControls(WellState& xw) const + { + if( ! wellsActive() ) return ; + + std::string modestring[4] = { "BHP", "THP", "RESERVOIR_RATE", "SURFACE_RATE" }; + // Find, for each well, if any constraints are broken. If so, + // switch control to first broken constraint. + const int np = wellsMultiSegment()[0]->numberOfPhases(); + const int nw = wellsMultiSegment().size(); + for (int w = 0; w < nw; ++w) { + const WellControls* wc = wellsMultiSegment()[w]->wellControls(); + // The current control in the well state overrides + // the current control set in the Wells struct, which + // is instead treated as a default. + int current = xw.currentControls()[w]; + // Loop over all controls except the current one, and also + // skip any RESERVOIR_RATE controls, since we cannot + // handle those. + const int nwc = well_controls_get_num(wc); + int ctrl_index = 0; + for (; ctrl_index < nwc; ++ctrl_index) { + if (ctrl_index == current) { + // This is the currently used control, so it is + // used as an equation. So this is not used as an + // inequality constraint, and therefore skipped. + continue; + } + if (detail::constraintBroken( + xw.bhp(), xw.thp(), xw.wellRates(), + w, np, wellsMultiSegment()[w]->wellType(), wc, ctrl_index)) { + // ctrl_index will be the index of the broken constraint after the loop. + break; + } + } + + if (ctrl_index != nwc) { + // Constraint number ctrl_index was broken, switch to it. + if (terminal_output_) + { + std::cout << "Switching control mode for well " << wellsMultiSegment()[w]->name() + << " from " << modestring[well_controls_iget_type(wc, current)] + << " to " << modestring[well_controls_iget_type(wc, ctrl_index)] << std::endl; + } + xw.currentControls()[w] = ctrl_index; + current = xw.currentControls()[w]; + } + + // Get gravity for THP hydrostatic corrrection + // const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); + + // Updating well state and primary variables. + // Target values are used as initial conditions for BHP, THP, and SURFACE_RATE + const double target = well_controls_iget_target(wc, current); + const double* distr = well_controls_iget_distr(wc, current); + switch (well_controls_iget_type(wc, current)) { + case BHP: + xw.bhp()[w] = target; + xw.segPress()[xw.topSegmentLoc()[w]] = target; + break; + + case THP: { + OPM_THROW(std::runtime_error, "THP control is not implemented for multi-sgement wells yet!!"); + } + + case RESERVOIR_RATE: + // No direct change to any observable quantity at + // surface condition. In this case, use existing + // flow rates as initial conditions as reservoir + // rate acts only in aggregate. + break; + + case SURFACE_RATE: + for (int phase = 0; phase < np; ++phase) { + if (distr[phase] > 0.0) { + xw.wellRates()[np * w + phase] = target * distr[phase]; + // TODO: consider changing all (not just top) segment rates + // to make them consistent, it could possibly improve convergence. + xw.segPhaseRates()[np * xw.topSegmentLoc()[w] + phase] = target * distr[phase]; + } + } + break; + } + + } + } + + + + + + template + bool BlackoilMultiSegmentModel::solveWellEq(const std::vector& mob_perfcells, + const std::vector& b_perfcells, + SolutionState& state, + WellState& well_state) + { + const bool converged = Base::solveWellEq(mob_perfcells, b_perfcells, state, well_state); + + if (converged) { + // We must now update the state.segp and state.segqs members, + // that the base version does not know about. + const int np = numPhases(); + const int nseg_total =well_state.numSegments(); + { + // We will set the segp primary variable to the new ones, + // but we do not change the derivatives here. + ADB::V new_segp = Eigen::Map(well_state.segPress().data(), nseg_total); + // Avoiding the copy below would require a value setter method + // in AutoDiffBlock. + std::vector old_segp_derivs = state.segp.derivative(); + state.segp = ADB::function(std::move(new_segp), std::move(old_segp_derivs)); + } + { + // Need to reshuffle well rates, from phase running fastest + // to wells running fastest. + // The transpose() below switches the ordering. + const DataBlock segrates = Eigen::Map(well_state.segPhaseRates().data(), nseg_total, np).transpose(); + ADB::V new_segqs = Eigen::Map(segrates.data(), nseg_total * np); + std::vector old_segqs_derivs = state.segqs.derivative(); + state.segqs = ADB::function(std::move(new_segqs), std::move(old_segqs_derivs)); + } + + // This is also called by the base version, but since we have updated + // state.segp we must call it again. + asImpl().computeWellConnectionPressures(state, well_state); + } + + return converged; + } + + + + + + template + void BlackoilMultiSegmentModel::addWellControlEq(const SolutionState& state, + const WellState& xw, + const V& aliveWells) + { + // the name of the function is a a little misleading. + // Basically it is the function for the pressure equation. + // And also, it work as the control equation when it is the segment + if( wellsMultiSegment().empty() ) return; + + const int np = numPhases(); + const int nw = wellsMultiSegment().size(); + const int nseg_total = xw.numSegments(); + + ADB aqua = ADB::constant(ADB::V::Zero(nseg_total)); + ADB liquid = ADB::constant(ADB::V::Zero(nseg_total)); + ADB vapour = ADB::constant(ADB::V::Zero(nseg_total)); + + if (active_[Water]) { + aqua += subset(state.segqs, Span(nseg_total, 1, BlackoilPhases::Aqua * nseg_total)); + } + if (active_[Oil]) { + liquid += subset(state.segqs, Span(nseg_total, 1, BlackoilPhases::Liquid * nseg_total)); + } + if (active_[Gas]) { + vapour += subset(state.segqs, Span(nseg_total, 1, BlackoilPhases::Vapour * nseg_total)); + } + + // THP control is not implemented for the moment. + + // Hydrostatic correction variables + ADB::V rho_v = ADB::V::Zero(nw); + ADB::V vfp_ref_depth_v = ADB::V::Zero(nw); + + // Target vars + ADB::V bhp_targets = ADB::V::Zero(nw); + ADB::V rate_targets = ADB::V::Zero(nw); + Eigen::SparseMatrix rate_distr(nw, np*nw); + + // Selection variables + // well selectors + std::vector bhp_well_elems; + std::vector rate_well_elems; + // segment selectors + std::vector bhp_top_elems; + std::vector rate_top_elems; + std::vector rate_top_phase_elems; + std::vector others_elems; + + //Run through all wells to calculate BHP/RATE targets + //and gather info about current control + int start_segment = 0; + for (int w = 0; w < nw; ++w) { + const struct WellControls* wc = wellsMultiSegment()[w]->wellControls(); + + // The current control in the well state overrides + // the current control set in the Wells struct, which + // is instead treated as a default. + const int current = xw.currentControls()[w]; + + const int nseg = wellsMultiSegment()[w]->numberOfSegments(); + + switch (well_controls_iget_type(wc, current)) { + case BHP: + { + bhp_well_elems.push_back(w); + bhp_top_elems.push_back(start_segment); + bhp_targets(w) = well_controls_iget_target(wc, current); + rate_targets(w) = -1e100; + for (int p = 0; p < np; ++p) { + rate_top_phase_elems.push_back(np * start_segment + p); + } + } + break; + + case THP: + { + OPM_THROW(std::runtime_error, "THP control is not implemented for multi-sgement wells yet!!"); + } + break; + + case RESERVOIR_RATE: // Intentional fall-through + case SURFACE_RATE: + { + rate_well_elems.push_back(w); + rate_top_elems.push_back(start_segment); + for (int p = 0; p < np; ++p) { + rate_top_phase_elems.push_back(np * start_segment + p); + } + // RESERVOIR and SURFACE rates look the same, from a + // high-level point of view, in the system of + // simultaneous linear equations. + + const double* const distr = + well_controls_iget_distr(wc, current); + + for (int p = 0; p < np; ++p) { + rate_distr.insert(w, p*nw + w) = distr[p]; + } + + bhp_targets(w) = -1.0e100; + rate_targets(w) = well_controls_iget_target(wc, current); + } + break; + + } + + for (int i = 1; i < nseg; ++i) { + others_elems.push_back(i + start_segment); + } + start_segment += nseg; + } + + // for each segment: 1, if the segment is the top segment, then control equation + // 2, if the segment is not the top segment, then the pressure equation + const ADB bhp_residual = subset(state.segp, bhp_top_elems) - subset(bhp_targets, bhp_well_elems); + const ADB rate_residual = subset(rate_distr * subset(state.segqs, rate_top_phase_elems) - rate_targets, rate_well_elems); + + ADB others_residual = ADB::constant(V::Zero(nseg_total)); + + // Special handling for when we are called from solveWellEq(). + // TODO: restructure to eliminate need for special treatmemt. + ADB wspd = (state.segp.numBlocks() == 2) + ? detail::onlyWellDerivs(well_segment_pressures_delta_) + : well_segment_pressures_delta_; + + others_residual = wops_ms_.eliminate_topseg * (state.segp - wops_ms_.s2s_outlet * state.segp + wspd); + + // all the control equations + // TODO: can be optimized better + ADB well_eq_topsegment = subset(superset(bhp_residual, bhp_top_elems, nseg_total) + + superset(rate_residual, rate_top_elems, nseg_total), xw.topSegmentLoc()); + + // For wells that are dead (not flowing), and therefore not communicating + // with the reservoir, we set the equation to be equal to the well's total + // flow. This will be a solution only if the target rate is also zero. + Eigen::SparseMatrix rate_summer(nw, np*nw); + for (int w = 0; w < nw; ++w) { + for (int phase = 0; phase < np; ++phase) { + rate_summer.insert(w, phase*nw + w) = 1.0; + } + } + Selector alive_selector(aliveWells, Selector::NotEqualZero); + // TODO: Here only handles the wells, or the top segments + // should we also handle some non-alive non-top segments? + // should we introduce the cocept of non-alive segments? + // At the moment, we only handle the control equations + well_eq_topsegment = alive_selector.select(well_eq_topsegment, rate_summer * subset(state.segqs, rate_top_phase_elems)); + + /* residual_.well_eq = superset(bhp_residual, bhp_top_elems, nseg_total) + + superset(rate_residual, rate_top_elems, nseg_total) + + superset(others_residual, others_elems, nseg_total); */ + residual_.well_eq = superset(well_eq_topsegment, xw.topSegmentLoc(), nseg_total) + + others_residual; + + } + + + + + + template + void + BlackoilMultiSegmentModel::updateWellState(const V& dwells, + WellState& well_state) + { + + if (!wellsMultiSegment().empty()) + { + const int np = numPhases(); + const int nw = wellsMultiSegment().size(); + const int nseg_total = well_state.numSegments(); + + // Extract parts of dwells corresponding to each part. + int varstart = 0; + const V dsegqs = subset(dwells, Span(np * nseg_total, 1, varstart)); + varstart += dsegqs.size(); + const V dsegp = subset(dwells, Span(nseg_total, 1, varstart)); + varstart += dsegp.size(); + assert(varstart == dwells.size()); + const double dpmaxrel = dpMaxRel(); + + + // segment phase rates update + // in dwells, the phase rates are ordered by phase. + // while in WellStateMultiSegment, the phase rates are ordered by segments + const DataBlock wsr = Eigen::Map(dsegqs.data(), np, nseg_total).transpose(); + const V dwsr = Eigen::Map(wsr.data(), nseg_total * np); + const V wsr_old = Eigen::Map(&well_state.segPhaseRates()[0], nseg_total * np); + const V sr = wsr_old - dwsr; + std::copy(&sr[0], &sr[0] + sr.size(), well_state.segPhaseRates().begin()); + + + // segment pressure updates + const V segp_old = Eigen::Map(&well_state.segPress()[0], nseg_total, 1); + // TODO: applying the pressure change limiter to all the segments, not sure if it is the correct thing to do + const V dsegp_limited = sign(dsegp) * dsegp.abs().min(segp_old.abs() * dpmaxrel); + const V segp = segp_old - dsegp_limited; + std::copy(&segp[0], &segp[0] + segp.size(), well_state.segPress().begin()); + + // update the well rates and bhps, which are not anymore primary vabriables. + // they are updated directly from the updated segment phase rates and segment pressures. + + // Bhp update. + V bhp = V::Zero(nw); + V wr = V::Zero(nw * np); + // it is better to use subset + + int start_segment = 0; + for (int w = 0; w < nw; ++w) { + bhp[w] = well_state.segPress()[start_segment]; + // insert can be faster + for (int p = 0; p < np; ++p) { + wr[p + np * w] = well_state.segPhaseRates()[p + np * start_segment]; + } + + const int nseg = wellsMultiSegment()[w]->numberOfSegments(); + start_segment += nseg; + } + + assert(start_segment == nseg_total); + std::copy(&bhp[0], &bhp[0] + bhp.size(), well_state.bhp().begin()); + std::copy(&wr[0], &wr[0] + wr.size(), well_state.wellRates().begin()); + + // TODO: handling the THP control related. + } + } + + + + + template + void + BlackoilMultiSegmentModel::computeSegmentFluidProperties(const SolutionState& state) + { + const int nw = wellsMultiSegment().size(); + const int nseg_total = state.segp.size(); + const int np = numPhases(); + + // although we will calculate segment density for non-segmented wells at the same time, + // while under most of the cases, they will not be used, + // since for most of the cases, the density calculation for non-segment wells are + // set to be 'SEG' way, which is not a option for multi-segment wells. + // When the density calcuation for non-segmented wells are set to 'AVG', then + // the density calculation of the mixtures can be the same, while it remains to be verified. + + // The grid cells associated with segments. + // TODO: shoud be computed once and stored in WellState or global Wells structure or class. + std::vector segment_cells; + segment_cells.reserve(nseg_total); + for (int w = 0; w < nw; ++w) { + const std::vector& segment_cells_well = wellsMultiSegment()[w]->segmentCells(); + segment_cells.insert(segment_cells.end(), segment_cells_well.begin(), segment_cells_well.end()); + } + assert(int(segment_cells.size()) == nseg_total); + + const ADB segment_temp = subset(state.temperature, segment_cells); + // using the segment pressure or the average pressure + // using the segment pressure first + const ADB& segment_press = state.segp; + + // Compute PVT properties for segments. + std::vector segment_cond(nseg_total); + const std::vector& pc = phaseCondition(); + for (int s = 0; s < nseg_total; ++s) { + segment_cond[s] = pc[segment_cells[s]]; + } + std::vector b_seg(np, ADB::null()); + // Viscosities for different phases + std::vector mu_seg(np, ADB::null()); + ADB rsmax_seg = ADB::null(); + ADB rvmax_seg = ADB::null(); + const PhaseUsage& pu = fluid_.phaseUsage(); + if (pu.phase_used[Water]) { + b_seg[pu.phase_pos[Water]] = fluid_.bWat(segment_press, segment_temp, segment_cells); + mu_seg[pu.phase_pos[Water]] = fluid_.muWat(segment_press, segment_temp, segment_cells); + } + assert(active_[Oil]); + const ADB segment_so = subset(state.saturation[pu.phase_pos[Oil]], segment_cells); + if (pu.phase_used[Oil]) { + const ADB segment_rs = subset(state.rs, segment_cells); + b_seg[pu.phase_pos[Oil]] = fluid_.bOil(segment_press, segment_temp, segment_rs, + segment_cond, segment_cells); + rsmax_seg = fluidRsSat(segment_press, segment_so, segment_cells); + mu_seg[pu.phase_pos[Oil]] = fluid_.muOil(segment_press, segment_temp, segment_rs, + segment_cond, segment_cells); + } + assert(active_[Gas]); + if (pu.phase_used[Gas]) { + const ADB segment_rv = subset(state.rv, segment_cells); + b_seg[pu.phase_pos[Gas]] = fluid_.bGas(segment_press, segment_temp, segment_rv, + segment_cond, segment_cells); + rvmax_seg = fluidRvSat(segment_press, segment_so, segment_cells); + mu_seg[pu.phase_pos[Gas]] = fluid_.muGas(segment_press, segment_temp, segment_rv, + segment_cond, segment_cells); + } + + // Extract segment flow by phase (segqs) and compute total surface rate. + ADB tot_surface_rate = ADB::constant(V::Zero(nseg_total)); + std::vector segqs(np, ADB::null()); + for (int phase = 0; phase < np; ++phase) { + segqs[phase] = subset(state.segqs, Span(nseg_total, 1, phase * nseg_total)); + tot_surface_rate += segqs[phase]; + } + + // TODO: later this will be implmented as a global mapping + std::vector> comp_frac(np, std::vector(nseg_total, 0.0)); + int start_segment = 0; + for (int w = 0; w < nw; ++w) { + WellMultiSegmentConstPtr well = wellsMultiSegment()[w]; + const int nseg = well->numberOfSegments(); + const std::vector& comp_frac_well = well->compFrac(); + for (int phase = 0; phase < np; ++phase) { + for (int s = 0; s < nseg; ++s) { + comp_frac[phase][s + start_segment] = comp_frac_well[phase]; + } + } + start_segment += nseg; + } + assert(start_segment == nseg_total); + + // Compute mix. + // 'mix' contains the component fractions under surface conditions. + std::vector mix(np, ADB::null()); + for (int phase = 0; phase < np; ++phase) { + // initialize to be the compFrac for each well, + // then update only the one with non-zero total volume rate + mix[phase] = ADB::constant(Eigen::Map(comp_frac[phase].data(), nseg_total)); + } + // There should be a better way to do this. + Selector non_zero_tot_rate(tot_surface_rate.value(), Selector::NotEqualZero); + for (int phase = 0; phase < np; ++phase) { + mix[phase] = non_zero_tot_rate.select(segqs[phase] / tot_surface_rate, mix[phase]); + } + + // Calculate rs and rv. + ADB rs = ADB::constant(V::Zero(nseg_total)); + ADB rv = rs; + const int gaspos = pu.phase_pos[Gas]; + const int oilpos = pu.phase_pos[Oil]; + Selector non_zero_mix_oilpos(mix[oilpos].value(), Selector::GreaterZero); + Selector non_zero_mix_gaspos(mix[gaspos].value(), Selector::GreaterZero); + // What is the better way to do this? + // big values should not be necessary + ADB big_values = ADB::constant(V::Constant(nseg_total, 1.e100)); + ADB mix_gas_oil = non_zero_mix_oilpos.select(mix[gaspos] / mix[oilpos], big_values); + ADB mix_oil_gas = non_zero_mix_gaspos.select(mix[oilpos] / mix[gaspos], big_values); + if (active_[Oil]) { + V selectorUnderRsmax = V::Zero(nseg_total); + V selectorAboveRsmax = V::Zero(nseg_total); + for (int s = 0; s < nseg_total; ++s) { + if (mix_gas_oil.value()[s] > rsmax_seg.value()[s]) { + selectorAboveRsmax[s] = 1.0; + } else { + selectorUnderRsmax[s] = 1.0; + } + } + rs = non_zero_mix_oilpos.select(selectorAboveRsmax * rsmax_seg + selectorUnderRsmax * mix_gas_oil, rs); + } + if (active_[Gas]) { + V selectorUnderRvmax = V::Zero(nseg_total); + V selectorAboveRvmax = V::Zero(nseg_total); + for (int s = 0; s < nseg_total; ++s) { + if (mix_oil_gas.value()[s] > rvmax_seg.value()[s]) { + selectorAboveRvmax[s] = 1.0; + } else { + selectorUnderRvmax[s] = 1.0; + } + } + rv = non_zero_mix_gaspos.select(selectorAboveRvmax * rvmax_seg + selectorUnderRvmax * mix_oil_gas, rv); + } + + // Calculate the phase fraction under reservoir conditions. + std::vector x(np, ADB::null()); + for (int phase = 0; phase < np; ++phase) { + x[phase] = mix[phase]; + } + if (active_[Gas] && active_[Oil]) { + x[gaspos] = (mix[gaspos] - mix[oilpos] * rs) / (V::Ones(nseg_total) - rs * rv); + x[oilpos] = (mix[oilpos] - mix[gaspos] * rv) / (V::Ones(nseg_total) - rs * rv); + } + + // Compute total reservoir volume to surface volume ratio. + ADB volrat = ADB::constant(V::Zero(nseg_total)); + for (int phase = 0; phase < np; ++phase) { + volrat += x[phase] / b_seg[phase]; + } + + // Compute segment densities. + ADB dens = ADB::constant(V::Zero(nseg_total)); + for (int phase = 0; phase < np; ++phase) { + const V surface_density = fluid_.surfaceDensity(phase, segment_cells); + dens += surface_density * mix[phase]; + } + well_segment_densities_ = dens / volrat; + + // Calculating the surface volume of each component in the segment + assert(np == int(segment_comp_surf_volume_current_.size())); + const ADB segment_surface_volume = segvdt_ / volrat; + for (int phase = 0; phase < np; ++phase) { + segment_comp_surf_volume_current_[phase] = segment_surface_volume * mix[phase]; + } + + // Mass flow rate of the segments + segment_mass_flow_rates_ = ADB::constant(V::Zero(nseg_total)); + for (int phase = 0; phase < np; ++phase) { + // TODO: how to remove one repeated surfaceDensity() + const V surface_density = fluid_.surfaceDensity(phase, segment_cells); + segment_mass_flow_rates_ += surface_density * segqs[phase]; + } + + // Viscosity of the fluid mixture in the segments + segment_viscosities_ = ADB::constant(V::Zero(nseg_total)); + for (int phase = 0; phase < np; ++phase) { + segment_viscosities_ += x[phase] * mu_seg[phase]; + } + } + + + template + void + BlackoilMultiSegmentModel::computeSegmentPressuresDelta(const SolutionState& state) + { + const int nw = wellsMultiSegment().size(); + const int nseg_total = state.segp.size(); + + // calculate the depth difference of the segments + // TODO: we need to store the following values somewhere to avoid recomputation. + V segment_depth_delta = V::Zero(nseg_total); + int start_segment = 0; + for (int w = 0; w < nw; ++w) { + WellMultiSegmentConstPtr well = wellsMultiSegment()[w]; + const int nseg = well->numberOfSegments(); + for (int s = 1; s < nseg; ++s) { + const int s_outlet = well->outletSegment()[s]; + assert(s_outlet >= 0 && s_outlet < nseg); + segment_depth_delta[s + start_segment] = well->segmentDepth()[s_outlet] - well->segmentDepth()[s]; + } + start_segment += nseg; + } + assert(start_segment == nseg_total); + + const double grav = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); + const ADB grav_adb = ADB::constant(V::Constant(nseg_total, grav)); + well_segment_pressures_delta_ = segment_depth_delta * grav_adb * well_segment_densities_; + + ADB well_segment_perforation_densities = wops_ms_.s2p * well_segment_densities_; + well_segment_perforation_pressure_diffs_ = grav * well_segment_perforation_depth_diffs_ * well_segment_perforation_densities; + + } + +} // namespace Opm + +#endif // OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilMultiSegment.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilMultiSegment.hpp new file mode 100644 index 000000000..edac5d653 --- /dev/null +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilMultiSegment.hpp @@ -0,0 +1,104 @@ +/* + Copyright 2013 SINTEF ICT, Applied Mathematics. + Copyright 2015 Andreas Lauser + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_SIMULATORFULLYIMPLICITBLACKOILMULTISEGMENT_HEADER_INCLUDED +#define OPM_SIMULATORFULLYIMPLICITBLACKOILMULTISEGMENT_HEADER_INCLUDED + +#include + + +#include +#include +#include + +namespace Opm { + +template +class SimulatorFullyImplicitBlackoilMultiSegment; + +template +struct SimulatorTraits > +{ + typedef WellStateMultiSegment WellState; + typedef BlackoilState ReservoirState; + typedef BlackoilOutputWriter OutputWriter; + typedef GridT Grid; + typedef BlackoilMultiSegmentModel Model; + typedef NonlinearSolver Solver; +}; + +/// a simulator for the blackoil model +template +class SimulatorFullyImplicitBlackoilMultiSegment + : public SimulatorBase > +{ + typedef SimulatorBase > Base; + typedef SimulatorFullyImplicitBlackoilMultiSegment ThisType; + typedef SimulatorTraits Traits; + typedef typename Traits::ReservoirState ReservoirState; + typedef typename Traits::WellState WellState; + typedef typename Traits::Solver Solver; +public: + // forward the constructor to the base class + SimulatorFullyImplicitBlackoilMultiSegment(const parameter::ParameterGroup& param, + const GridT& grid, + DerivedGeology& geo, + BlackoilPropsAdInterface& props, + const RockCompressibility* rock_comp_props, + NewtonIterationBlackoilInterface& linsolver, + const double* gravity, + const bool disgas, + const bool vapoil, + std::shared_ptr eclipse_state, + BlackoilOutputWriter& output_writer, + const std::vector& threshold_pressures_by_face) + : Base(param, grid, geo, props, rock_comp_props, linsolver, gravity, disgas, vapoil, + eclipse_state, output_writer, threshold_pressures_by_face) + {} + + + SimulatorReport run(SimulatorTimer& timer, + ReservoirState& state); + +protected: + + std::unique_ptr createSolver(const Wells* wells, std::vector& wells_multisegment); + + using Base::output_writer_; + using Base::param_; + using Base::solver_; + using Base::terminal_output_; + using Base::eclipse_state_; + using Base::grid_; + using Base::props_; + using Base::is_parallel_run_; + using Base::allcells_; + using Base::model_param_; + using Base::geo_; + using Base::rock_comp_props_; + using Base::has_disgas_; + using Base::has_vapoil_; +}; + +} // namespace Opm + +#include "SimulatorFullyImplicitBlackoilMultiSegment_impl.hpp" + +#endif // OPM_SIMULATORFULLYIMPLICITBLACKOILMULTISEGMENT_HEADER_INCLUDED diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilMultiSegment_impl.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilMultiSegment_impl.hpp new file mode 100644 index 000000000..77ec6c5ea --- /dev/null +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilMultiSegment_impl.hpp @@ -0,0 +1,213 @@ +/* + Copyright 2013 SINTEF ICT, Applied Mathematics. + Copyright 2014 IRIS AS + Copyright 2015 Andreas Lauser + + 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 . +*/ + + +namespace Opm +{ + + template + auto SimulatorFullyImplicitBlackoilMultiSegment:: + createSolver(const Wells* wells, std::vector& wells_multisegment) + -> std::unique_ptr + { + typedef typename Traits::Model Model; + + auto model = std::unique_ptr(new Model(model_param_, + grid_, + props_, + geo_, + rock_comp_props_, + wells, + solver_, + eclipse_state_, + has_disgas_, + has_vapoil_, + terminal_output_, + wells_multisegment)); + + if (!Base::threshold_pressures_by_face_.empty()) { + model->setThresholdPressures(Base::threshold_pressures_by_face_); + } + + return std::unique_ptr(new Solver(Base::solver_param_, std::move(model))); + + } + + template + SimulatorReport SimulatorFullyImplicitBlackoilMultiSegment::run(SimulatorTimer& timer, + ReservoirState& state) + { + WellState prev_well_state; + + // Create timers and file for writing timing info. + Opm::time::StopWatch solver_timer; + double stime = 0.0; + Opm::time::StopWatch step_timer; + Opm::time::StopWatch total_timer; + total_timer.start(); + std::string tstep_filename = output_writer_.outputDirectory() + "/step_timing.txt"; + std::ofstream tstep_os(tstep_filename.c_str()); + + // adaptive time stepping + std::unique_ptr< AdaptiveTimeStepping > adaptiveTimeStepping; + if( param_.getDefault("timestep.adaptive", true ) ) + { + adaptiveTimeStepping.reset( new AdaptiveTimeStepping( param_, terminal_output_ ) ); + } + + // init output writer + output_writer_.writeInit( timer ); + + std::string restorefilename = param_.getDefault("restorefile", std::string("") ); + if( ! restorefilename.empty() ) + { + // -1 means that we'll take the last report step that was written + const int desiredRestoreStep = param_.getDefault("restorestep", int(-1) ); + output_writer_.restore( timer, state, prev_well_state, restorefilename, desiredRestoreStep ); + } + + unsigned int totalNonlinearIterations = 0; + unsigned int totalLinearIterations = 0; + + // Main simulation loop. + while (!timer.done()) { + // Report timestep. + step_timer.start(); + if ( terminal_output_ ) + { + timer.report(std::cout); + } + + // Create wells and well state. + WellsManager wells_manager(eclipse_state_, + timer.currentStepNum(), + Opm::UgGridHelpers::numCells(grid_), + Opm::UgGridHelpers::globalCell(grid_), + Opm::UgGridHelpers::cartDims(grid_), + Opm::UgGridHelpers::dimensions(grid_), + Opm::UgGridHelpers::cell2Faces(grid_), + Opm::UgGridHelpers::beginFaceCentroids(grid_), + props_.permeability(), + is_parallel_run_); + const Wells* wells = wells_manager.c_wells(); + WellState well_state; + // well_state.init(wells, state, prev_well_state); + + const std::vector& wells_ecl = eclipse_state_->getSchedule()->getWells(timer.currentStepNum()); + std::vector wells_multisegment; + wells_multisegment.reserve(wells_ecl.size()); + for (size_t i = 0; i < wells_ecl.size(); ++i) { + // not processing SHUT wells. + if (wells_ecl[i]->getStatus(timer.currentStepNum()) == WellCommon::SHUT) { + continue; + } + // checking if the well can be found in the wells + const std::string& well_name = wells_ecl[i]->name(); + // number of wells in wells + const int nw_wells = wells->number_of_wells; + int index_well; + for (index_well = 0; index_well < nw_wells; ++index_well) { + if (well_name == std::string(wells->name[index_well])) { + break; + } + } + + if (index_well != nw_wells) { // found in the wells + wells_multisegment.push_back(std::make_shared(wells_ecl[i], timer.currentStepNum(), wells)); + } + } + + well_state.init(wells_multisegment, state, prev_well_state); + + // give the polymer and surfactant simulators the chance to do their stuff + Base::asImpl().handleAdditionalWellInflow(timer, wells_manager, well_state, wells); + + // write simulation state at the report stage + output_writer_.writeTimeStep( timer, state, well_state ); + + // Max oil saturation (for VPPARS), hysteresis update. + props_.updateSatOilMax(state.saturation()); + props_.updateSatHyst(state.saturation(), allcells_); + + // Compute reservoir volumes for RESV controls. + Base::asImpl().computeRESV(timer.currentStepNum(), wells, state, well_state); + + // Run a multiple steps of the solver depending on the time step control. + solver_timer.start(); + + auto solver = createSolver(wells, wells_multisegment); + + // If sub stepping is enabled allow the solver to sub cycle + // in case the report steps are too large for the solver to converge + // + // \Note: The report steps are met in any case + // \Note: The sub stepping will require a copy of the state variables + if( adaptiveTimeStepping ) { + adaptiveTimeStepping->step( timer, *solver, state, well_state, output_writer_ ); + } + else { + // solve for complete report step + solver->step(timer.currentStepLength(), state, well_state); + } + + // take time that was used to solve system for this reportStep + solver_timer.stop(); + + // accumulate the number of nonlinear and linear Iterations + totalNonlinearIterations += solver->nonlinearIterations(); + totalLinearIterations += solver->linearIterations(); + + // Report timing. + const double st = solver_timer.secsSinceStart(); + + if ( terminal_output_ ) + { + std::cout << "Fully implicit solver took: " << st << " seconds." << std::endl; + } + + stime += st; + if ( output_writer_.output() ) { + SimulatorReport step_report; + step_report.pressure_time = st; + step_report.total_time = step_timer.secsSinceStart(); + step_report.reportParam(tstep_os); + } + + // Increment timer, remember well state. + ++timer; + prev_well_state = well_state; + } + + // Write final simulation state. + output_writer_.writeTimeStep( timer, state, prev_well_state ); + + // Stop timer and create timing report + total_timer.stop(); + SimulatorReport report; + report.pressure_time = stime; + report.transport_time = 0.0; + report.total_time = total_timer.secsSinceStart(); + report.total_newton_iterations = totalNonlinearIterations; + report.total_linear_iterations = totalLinearIterations; + return report; + } + +} // namespace Opm diff --git a/opm/autodiff/WellMultiSegment.cpp b/opm/autodiff/WellMultiSegment.cpp new file mode 100644 index 000000000..97cc29926 --- /dev/null +++ b/opm/autodiff/WellMultiSegment.cpp @@ -0,0 +1,400 @@ + +/* + Copyright 2015 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#include + + +namespace Opm +{ + + WellMultiSegment::WellMultiSegment(WellConstPtr well, size_t time_step, const Wells* wells) { + m_well_name_ = well->name(); + if (well->isMultiSegment(time_step)) { + initMultiSegmentWell(well, time_step, wells); + } else { + initNonMultiSegmentWell(well, time_step, wells); + } + updateWellOps(); + } + + void WellMultiSegment::initMultiSegmentWell(WellConstPtr well, size_t time_step, const Wells* wells) { + + CompletionSetConstPtr completion_set = well->getCompletions(time_step); + + m_is_multi_segment_ = true; + SegmentSetConstPtr segment_set = well->getSegmentSet(time_step); + m_number_of_segments_ = segment_set->numberSegment(); + m_number_of_perforations_ = completion_set->size(); + m_comp_pressure_drop_ = segment_set->compPressureDrop(); + m_multiphase_model_ = segment_set->multiPhaseModel(); + + m_outlet_segment_.resize(m_number_of_segments_); + m_inlet_segments_.resize(m_number_of_segments_); + m_segment_length_.resize(m_number_of_segments_); + m_segment_depth_.resize(m_number_of_segments_); + m_segment_internal_diameter_.resize(m_number_of_segments_); + m_segment_roughness_.resize(m_number_of_segments_); + m_segment_cross_area_.resize(m_number_of_segments_, 0.); + m_segment_volume_.resize(m_number_of_segments_); + m_segment_perforations_.resize(m_number_of_segments_); + + // we change the ID to location now for easier use later. + for (int i = 0; i < m_number_of_segments_; ++i) { + // The segment number for top segment is 0, the segment number of its outlet segment will be -1 + m_outlet_segment_[i] = segment_set->numberToLocation((*segment_set)[i]->outletSegment()); + m_segment_length_[i] = (*segment_set)[i]->totalLength(); + m_segment_depth_[i] = (*segment_set)[i]->depth(); + m_segment_internal_diameter_[i] = (*segment_set)[i]->internalDiameter(); + m_segment_roughness_[i] = (*segment_set)[i]->roughness(); + m_segment_cross_area_[i] = (*segment_set)[i]->crossArea(); + m_segment_volume_[i] = (*segment_set)[i]->volume(); + } + + // update the completion related information + // find the location of the well in wells + int index_well; + for (index_well = 0; index_well < wells->number_of_wells; ++index_well) { + if (m_well_name_ == std::string(wells->name[index_well])) { + break; + } + } + + std::vector temp_well_cell; + std::vector temp_well_index; + + if (index_well == wells->number_of_wells) { + throw std::runtime_error(" did not find the well " + m_well_name_ + "\n"); + } else { + + m_well_type_ = wells->type[index_well]; + m_well_controls_ = wells->ctrls[index_well]; + m_number_of_phases_ = wells->number_of_phases; + m_comp_frac_.resize(m_number_of_phases_); + std::copy(wells->comp_frac + index_well * m_number_of_phases_, + wells->comp_frac + (index_well + 1) * m_number_of_phases_, m_comp_frac_.begin()); + + int index_begin = wells->well_connpos[index_well]; + int index_end = wells->well_connpos[index_well + 1]; + + + for(int i = index_begin; i < index_end; ++i) { + // copy the WI and well_cell_ informatin to m_well_index_ and m_well_cell_ + // maybe also the depth of the perforations. + temp_well_cell.push_back(wells->well_cells[i]); + temp_well_index.push_back(wells->WI[i]); + } + } + + std::vector temp_perf_depth; + temp_perf_depth.resize(m_number_of_perforations_); + + for (int i = 0; i < (int)completion_set->size(); ++i) { + int i_segment = completion_set->get(i)->getSegmentNumber(); + // using the location of the segment in the array as the segment number/id. + // TODO: it can be helpful for output or postprocessing if we can keep the original number. + i_segment = segment_set->numberToLocation(i_segment); + m_segment_perforations_[i_segment].push_back(i); + temp_perf_depth[i] = completion_set->get(i)->getCenterDepth(); + } + + // reordering the perforation related informations + // so that the perforations belong to the same segment will be continuous + m_well_cell_.resize(m_number_of_perforations_); + m_well_index_.resize(m_number_of_perforations_); + m_perf_depth_.resize(m_number_of_perforations_); + m_segment_cell_.resize(m_number_of_segments_, -1); + + int perf_count = 0; + for (int is = 0; is < m_number_of_segments_; ++is) { + // TODO: the grid cell related to a segment should be calculated based on the location + // of the segment node. + // As the current temporary solution, the grid cell related to a segment determined by the + // first perforation cell related to the segment. + // when no perforation is related to the segment, use it outlet segment's cell. + const int nperf = m_segment_perforations_[is].size(); + if (nperf > 0) { + const int first_perf_number = m_segment_perforations_[is][0]; + m_segment_cell_[is] = temp_well_cell[first_perf_number]; + for (int iperf = 0; iperf < nperf; ++iperf) { + const int perf_number = m_segment_perforations_[is][iperf]; + m_well_cell_[perf_count] = temp_well_cell[perf_number]; + m_well_index_[perf_count] = temp_well_index[perf_number]; + m_perf_depth_[perf_count] = temp_perf_depth[perf_number]; + m_segment_perforations_[is][iperf] = perf_count; + ++perf_count; + } + } else { + // using the cell of its outlet segment + const int i_outlet_segment = m_outlet_segment_[is]; + if (i_outlet_segment < 0) { + assert(is == 0); // it must be the top segment + OPM_THROW(std::logic_error, "Top segment is not related to any perforation, its related cell must be calculated based the location of its segment node, which is not implemented yet \n"); + } else { + if (m_well_cell_[i_outlet_segment] < 0) { + OPM_THROW(std::logic_error, "The segment cell of its outlet segment is not determined yet, the current implementation does not support this \n"); + } else { + m_segment_cell_[is] = m_segment_cell_[i_outlet_segment]; + } + } + } + } + + assert(perf_count == m_number_of_perforations_); + + // update m_inlet_segments_ + // top segment does not have a outlet segment + for (int is = 1; is < m_number_of_segments_; ++is) { + const int index_outlet = m_outlet_segment_[is]; + m_inlet_segments_[index_outlet].push_back(is); + } + + } + + void WellMultiSegment::initNonMultiSegmentWell(WellConstPtr well, size_t time_step, const Wells* wells) { + + CompletionSetConstPtr completion_set = well->getCompletions(time_step); + + m_is_multi_segment_ = false; + m_number_of_segments_ = 1; + m_comp_pressure_drop_ = WellSegment::H__; + m_multiphase_model_ = WellSegment::HO; + + m_outlet_segment_.resize(m_number_of_segments_, -1); + m_segment_length_.resize(m_number_of_segments_, 0.); + m_segment_depth_.resize(m_number_of_segments_, 0.); + m_segment_internal_diameter_.resize(m_number_of_segments_, 0.); + m_segment_roughness_.resize(m_number_of_segments_, 0.); + m_segment_cross_area_.resize(m_number_of_segments_, 0.); + m_segment_volume_.resize(m_number_of_segments_, 0.); + m_segment_perforations_.resize(m_number_of_segments_); + + // update the completion related information + int index_well; + for (index_well = 0; index_well < wells->number_of_wells; ++index_well) { + if (m_well_name_ == std::string(wells->name[index_well])) { + break; + } + } + + if (index_well == wells->number_of_wells) { + throw std::runtime_error(" did not find the well " + m_well_name_ + "\n"); + } else { + m_well_type_ = wells->type[index_well]; + m_well_controls_ = wells->ctrls[index_well]; + m_number_of_phases_ = wells->number_of_phases; + // set the segment depth to be the bhp reference depth + m_segment_depth_[0] = wells->depth_ref[index_well]; + m_comp_frac_.resize(m_number_of_phases_); + std::copy(wells->comp_frac + index_well * m_number_of_phases_, + wells->comp_frac + (index_well + 1) * m_number_of_phases_, m_comp_frac_.begin()); + + int index_begin = wells->well_connpos[index_well]; + int index_end = wells->well_connpos[index_well + 1]; + m_number_of_perforations_ = index_end - index_begin; + + for(int i = index_begin; i < index_end; ++i) { + m_well_cell_.push_back(wells->well_cells[i]); + m_well_index_.push_back(wells->WI[i]); + } + m_segment_cell_.resize(1, -1); + m_segment_cell_[0] = m_well_cell_[0]; + } + + // TODO: not sure if we need the perf_depth_. + m_perf_depth_.resize(m_number_of_perforations_, 0.); + m_segment_perforations_[0].resize(m_number_of_perforations_); + + for (int i = 0; i < m_number_of_perforations_; ++i) { + m_segment_perforations_[0][i] = i; + m_perf_depth_[i] = completion_set->get(i)->getCenterDepth(); + } + + m_inlet_segments_.resize(m_number_of_segments_); + } + + void WellMultiSegment::updateWellOps() { + m_wops_.s2p = Matrix(m_number_of_perforations_, m_number_of_segments_); + m_wops_.p2s = Matrix(m_number_of_segments_, m_number_of_perforations_); + + typedef Eigen::Triplet Tri; + + std::vector s2p; + std::vector p2s; + + s2p.reserve(m_number_of_perforations_); + p2s.reserve(m_number_of_perforations_); + + for(int s = 0; s < (int)m_number_of_segments_; ++s) { + int temp_nperf = m_segment_perforations_[s].size(); + // some segment may not have any perforation + assert(temp_nperf >= 0); + for (int perf = 0; perf < temp_nperf; ++perf) { + const int index_perf = m_segment_perforations_[s][perf]; + s2p.push_back(Tri(index_perf, s, 1.0)); + p2s.push_back(Tri(s, index_perf, 1.0)); + } + } + m_wops_.s2p.setFromTriplets(s2p.begin(), s2p.end()); + m_wops_.p2s.setFromTriplets(p2s.begin(), p2s.end()); + + m_wops_.s2s_gather = Matrix(m_number_of_segments_, m_number_of_segments_); + std::vector s2s_gather; + + s2s_gather.reserve(m_number_of_segments_ * m_number_of_segments_); + + std::vector s2s_inlets; + s2s_inlets.reserve(m_number_of_segments_); + std::vector s2s_outlet; + s2s_outlet.reserve(m_number_of_segments_); + + for (int s = 0; s < (int)m_number_of_segments_; ++s) { + s2s_gather.push_back(Tri(s, s, 1.0)); + int s_outlet = m_outlet_segment_[s]; + if (s_outlet >=0) { + s2s_inlets.push_back(Tri(s_outlet, s, 1.0)); + s2s_outlet.push_back(Tri(s, s_outlet, 1.0)); + } + int temp_s = s; + while (m_outlet_segment_[temp_s] >=0) { + s2s_gather.push_back(Tri(m_outlet_segment_[temp_s], s, 1.0)); + temp_s = m_outlet_segment_[temp_s]; + } + } + + m_wops_.s2s_gather.setFromTriplets(s2s_gather.begin(), s2s_gather.end()); + + m_wops_.p2s_gather = Matrix(m_number_of_segments_, m_number_of_perforations_); + m_wops_.p2s_gather = m_wops_.s2s_gather * m_wops_.p2s; + + m_wops_.s2s_inlets = Matrix(m_number_of_segments_, m_number_of_segments_); + m_wops_.s2s_inlets.setFromTriplets(s2s_inlets.begin(), s2s_inlets.end()); + + m_wops_.s2s_outlet = Matrix(m_number_of_segments_, m_number_of_segments_); + m_wops_.s2s_outlet.setFromTriplets(s2s_outlet.begin(), s2s_outlet.end()); + + m_wops_.p2s_average = Matrix(m_number_of_segments_, m_number_of_perforations_); + std::vector p2s_average; + p2s_average.reserve(m_number_of_segments_); + + for (int s = 0; s < (int)m_number_of_segments_; ++s) { + const int nperf = m_segment_perforations_[s].size(); + if (nperf > 0) { + p2s_average.push_back(Tri(s, s, 1.0/nperf)); + } + } + + // constructing the diagonal matrix to do the averaging for p2s + Matrix temp_averaging_p2s = Matrix(m_number_of_segments_, m_number_of_segments_); + temp_averaging_p2s.setFromTriplets(p2s_average.begin(), p2s_average.end()); + m_wops_.p2s_average = temp_averaging_p2s * m_wops_.p2s; + } + + const std::string& WellMultiSegment::name() const { + return m_well_name_; + } + + bool WellMultiSegment::isMultiSegmented() const { + return m_is_multi_segment_; + } + + WellType WellMultiSegment::wellType() const { + return m_well_type_; + } + + const WellControls* WellMultiSegment::wellControls() const { + return m_well_controls_; + } + + int WellMultiSegment::numberOfPerforations() const { + return m_number_of_perforations_; + } + + int WellMultiSegment::numberOfSegments() const { + return m_number_of_segments_; + } + + std::string WellMultiSegment::compPressureDrop() const { + return WellSegment::CompPressureDropEnumToString(m_comp_pressure_drop_); + } + + const std::vector& WellMultiSegment::compFrac() const { + return m_comp_frac_; + } + + int WellMultiSegment::numberOfPhases() const { + return m_number_of_phases_; + } + + const std::vector& WellMultiSegment::wellIndex() const { + return m_well_index_; + } + + const std::vector& WellMultiSegment::perfDepth() const { + return m_perf_depth_; + } + + const std::vector& WellMultiSegment::wellCells() const { + return m_well_cell_; + } + + const std::vector& WellMultiSegment::segmentCells() const { + return m_segment_cell_; + } + + const std::vector& WellMultiSegment::outletSegment() const { + return m_outlet_segment_; + } + + const std::vector>& WellMultiSegment::inletSegments() const { + return m_inlet_segments_; + } + + const std::vector& WellMultiSegment::segmentLength() const { + return m_segment_length_; + } + + const std::vector& WellMultiSegment::segmentDepth() const { + return m_segment_depth_; + } + + const std::vector& WellMultiSegment::segmentDiameter() const { + return m_segment_internal_diameter_; + } + + const std::vector& WellMultiSegment::segmentCrossArea() const { + return m_segment_cross_area_; + } + + const std::vector& WellMultiSegment::segmentRoughness() const { + return m_segment_roughness_; + } + + const std::vector& WellMultiSegment::segmentVolume() const { + return m_segment_volume_; + } + + const std::vector>& WellMultiSegment::segmentPerforations() const { + return m_segment_perforations_; + } + + const WellMultiSegment::WellOps& WellMultiSegment::wellOps() const { + return m_wops_; + } +} diff --git a/opm/autodiff/WellMultiSegment.hpp b/opm/autodiff/WellMultiSegment.hpp new file mode 100644 index 000000000..428d2ea44 --- /dev/null +++ b/opm/autodiff/WellMultiSegment.hpp @@ -0,0 +1,227 @@ +/* + Copyright 2015 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_WELLMULTISEGMENT_HEADER_INCLUDED +#define OPM_WELLMULTISEGMENT_HEADER_INCLUDED + + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace Opm +{ + + class WellMultiSegment + { + public: + + typedef Eigen::SparseMatrix Matrix; + + /// Constructor of WellMultiSegment + /// \param[in] well information from EclipseState + /// \param[in] current time step + /// \param[in[ pointer to Wells structure, to be removed eventually + WellMultiSegment(WellConstPtr well, size_t time_step, const Wells* wells); + + /// Well name. + const std::string& name() const; + + /// Flag indicating if the well is a multi-segment well. + bool isMultiSegmented() const; + + /// Number of the perforations. + int numberOfPerforations() const; + + /// Number of the segments. + int numberOfSegments() const; + + /// Components of the pressure drop invloved. + /// HFA Hydrostatic + friction + acceleration + /// HF- Hydrostatic + friction + /// H-- Hydrostatic only. + std::string compPressureDrop() const; + + /// Well control. + const WellControls* wellControls() const; + + /// Component fractions for each well. + const std::vector& compFrac() const; + + /// Number of phases. + int numberOfPhases() const; + + /// Well type. + WellType wellType() const; + + /// Well productivity index. + const std::vector& wellIndex() const; + + /// Depth of the perforations. + const std::vector& perfDepth() const; + + /// Indices of the grid blocks that perforations are completed in. + const std::vector& wellCells() const; + + /// Indices of the gird blocks that segments locate at. + const std::vector& segmentCells() const; + + /// Outlet segments, a segment (except top segment) can only have one outlet segment. + /// For top segment, its outlet segments is -1 always, which means no outlet segment for top segment. + const std::vector& outletSegment() const; + + /// Inlet segments. a segment can have more than one inlet segments. + const std::vector>& inletSegments() const; + + /// The length of the segment nodes down the wellbore from the reference point. + const std::vector& segmentLength() const; + + /// The depth of the segment nodes. + const std::vector& segmentDepth() const; + + /// Tubing internal diameter. + const std::vector& segmentDiameter() const; + + /// Cross-sectional area. + const std::vector& segmentCrossArea() const; + + /// Effective absolute roughness of the tube. + const std::vector& segmentRoughness() const; + + /// Volume of segment. + const std::vector& segmentVolume() const; + + /// Perforations related to each segment. + const std::vector>& segmentPerforations() const; + + /// Struct for the well operator matrices. + /// All the operator matrics only apply to the one specifi well. + struct WellOps { + Matrix s2p; // segment -> perf (scatter) + Matrix p2s; // perf -> segment (gather) + Matrix p2s_average; // perf -> segment (avarage) + Matrix s2s_gather; // segment -> segment (in an accumlative way) + // means the outlet segments will gather all the contribution + // from all the inlet segments in a recurisive way + Matrix p2s_gather; // perforation -> segment (in an accumative way) + Matrix s2s_inlets; // segment -> its inlet segments + Matrix s2s_outlet; // segment -> its outlet segment + }; + + /// Well operator matrics + const WellOps& wellOps() const; + + private: + // for the moment, we use the information from wells. + // TODO: remove the dependency on wells from opm-core. + void initMultiSegmentWell(WellConstPtr well, size_t time_step, const Wells* wells); + void initNonMultiSegmentWell(WellConstPtr well, size_t time_step, const Wells* wells); + void updateWellOps(); + + private: + // name of the well + std::string m_well_name_; + // flag to indicate if this well is a + // multi-segmented well + bool m_is_multi_segment_; + // well type + // INJECTOR or PRODUCER + enum WellType m_well_type_; + // number of phases + int m_number_of_phases_; + // component fractions for each well + std::vector m_comp_frac_; + // controls for this well + // using pointer for temporary + // changing it when figuring out how to using it + struct WellControls *m_well_controls_; + // components of the pressure drop to be included + WellSegment::CompPressureDropEnum m_comp_pressure_drop_; + // multi-phase flow model + WellSegment::MultiPhaseModelEnum m_multiphase_model_; + // number of perforation for this well + int m_number_of_perforations_; + // number of segments for this well + int m_number_of_segments_; + // well index for each completion + std::vector m_well_index_; + // depth for each completion // form the keyword COMPSEGS? + std::vector m_perf_depth_; + // well cell for each completion + std::vector m_well_cell_; + // cell for each segment + std::vector m_segment_cell_; + // how to organize the segment structure here? + // indicate the outlet segment for each segment + // maybe here we can use the location in the vector + // at the moment, we still use the ID number + // then a mapping from the ID number to the actual location will be required + // The ID is already changed to the location now. + std::vector m_outlet_segment_; + // for convinience, we store the inlet segments for each segment + std::vector> m_inlet_segments_; + // this one is not necessary any more, since the segment number is its location. + // std::map m_number_to_location_; + // has not decided to use the absolute length from the well head + // or the length of this single segment + // using the absolute length first + std::vector m_segment_length_; + // the depth of the segmnet node + std::vector m_segment_depth_; + // the internal diameter of the segment + std::vector m_segment_internal_diameter_; + // the roughness of the segment + std::vector m_segment_roughness_; + // the cross area of the segment + std::vector m_segment_cross_area_; + // the volume of the segment + std::vector m_segment_volume_; + // the completions that is related to each segment + // the completions's ids are their location in the vector m_well_index_ + // m_well_cell_ + // This is also assuming the order of the completions in Well is the same with + // the order of the completions in wells. + std::vector> m_segment_perforations_; + + WellOps m_wops_; + }; + + typedef std::shared_ptr WellMultiSegmentPtr; + typedef std::shared_ptr WellMultiSegmentConstPtr; + +} // namespace Opm + + +#endif // OPM_WELLMULTISEGMENT_HEADER_INCLUDE diff --git a/opm/autodiff/WellStateMultiSegment.hpp b/opm/autodiff/WellStateMultiSegment.hpp new file mode 100644 index 000000000..a64111de8 --- /dev/null +++ b/opm/autodiff/WellStateMultiSegment.hpp @@ -0,0 +1,359 @@ +/* + Copyright 2015 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_WELLSTATEMULTISEGMENT_HEADER_INCLUDED +#define OPM_WELLSTATEMULTISEGMENT_HEADER_INCLUDED + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace Opm +{ + + /// The state of a set of multi-sgemnet wells + // Since we are avoiding to use the old Wells structure, + // it might be a good idea not to relate this State to the WellState much. + class WellStateMultiSegment + : public WellStateFullyImplicitBlackoil + { + public: + + typedef WellStateFullyImplicitBlackoil Base; + + typedef struct { + int well_number; + int start_segment; + int number_of_segments; + int start_perforation; + int number_of_perforations; + std::vector start_perforation_segment; // the starting position of perforation inside the segment + std::vector number_of_perforations_segment; // the numbers for perforations for the segments + } SegmentedMapentryType; + + typedef std::map SegmentedWellMapType; + + /// Allocate and initialize if wells is non-null. Also tries + /// to give useful initial values to the bhp(), wellRates() + /// and perfPhaseRates() fields, depending on controls + template + void init(const std::vector& wells, const ReservoirState& state, const PrevWellState& prevState) + { + const int nw = wells.size(); + nseg_ = 0; + nperf_ = 0; + if (nw == 0) { + perfPhaseRates().clear(); + perfPress().clear(); + segphaserates_.clear(); + segpress_.clear(); + return; + } + + const int np = wells[0]->numberOfPhases(); // number of the phases + + for (int iw = 0; iw < nw; ++iw) { + nperf_ += wells[iw]->numberOfPerforations(); + nseg_ += wells[iw]->numberOfSegments(); + } + + bhp().resize(nw); + thp().resize(nw); + top_segment_loc_.resize(nw); + temperature().resize(nw, 273.15 + 20); // standard temperature for now + + wellRates().resize(nw * np, 0.0); + + currentControls().resize(nw); + for(int iw = 0; iw < nw; ++iw) { + currentControls()[iw] = well_controls_get_current(wells[iw]->wellControls()); + } + + for (int iw = 0; iw < nw; ++iw) { + assert((wells[iw]->wellType() == INJECTOR) || (wells[iw]->wellType() == PRODUCER)); + } + + int start_segment = 0; + int start_perforation = 0; + + perfPhaseRates().resize(nperf_ * np, 0.0); + perfPress().resize(nperf_, -1.0e100); + segphaserates_.resize(nseg_ * np, 0.0); + segpress_.resize(nseg_, -1.0e100); + + + for (int w = 0; w < nw; ++w) { + assert((wells[w]->wellType() == INJECTOR) || (wells[w]->wellType() == PRODUCER)); + const WellControls* ctrl = wells[w]->wellControls(); + + std::string well_name(wells[w]->name()); + // Initialize the wellMap_ here + SegmentedMapentryType& wellMapEntry = segmentedWellMap_[well_name]; + wellMapEntry.well_number = w; + wellMapEntry.start_segment = start_segment; + wellMapEntry.number_of_segments = wells[w]->numberOfSegments(); + wellMapEntry.start_perforation = start_perforation; + wellMapEntry.number_of_perforations = wells[w]->numberOfPerforations(); + + top_segment_loc_[w] = start_segment; + + int start_perforation_segment = 0; + wellMapEntry.start_perforation_segment.resize(wellMapEntry.number_of_segments); + wellMapEntry.number_of_perforations_segment.resize(wellMapEntry.number_of_segments); + + for (int i = 0; i < wellMapEntry.number_of_segments; ++i) { + wellMapEntry.start_perforation_segment[i] = start_perforation_segment; + wellMapEntry.number_of_perforations_segment[i] = wells[w]->segmentPerforations()[i].size(); + start_perforation_segment += wellMapEntry.number_of_perforations_segment[i]; + } + assert(start_perforation_segment == wellMapEntry.number_of_perforations); + + + if (well_controls_well_is_stopped(ctrl)) { + // 1. WellRates: 0 + // 2. Bhp: assign bhp equal to bhp control, if applicable, otherwise + // assign equal to first perforation cell pressure. + if (well_controls_get_current_type(ctrl) == BHP) { + bhp()[w] = well_controls_get_current_target(ctrl); + } else { + const int first_cell = wells[0]->wellCells()[0]; + bhp()[w] = state.pressure()[first_cell]; + } + // 3. Thp: assign thp equal to thp control, if applicable, + // otherwise assign equal to bhp value. + if (well_controls_get_current_type(ctrl) == THP) { + thp()[w] = well_controls_get_current_target( ctrl ); + } else { + thp()[w] = bhp()[w]; + } + // 4. Perforation pressures and phase rates + // 5. Segment pressures and phase rates + + } else { + // Open Wells + // 1. Rates: initialize well rates to match controls if type is SURFACE_RATE. Otherwise, we + // cannot set the correct value here, so we aasign a small rate with the correct sign so that any + // logic depending on that sign will work as expected. + if (well_controls_get_current_type(ctrl) == SURFACE_RATE) { + const double rate_target = well_controls_get_current_target(ctrl); + const double * distr = well_controls_get_current_distr( ctrl ); + for (int p = 0; p < np; ++p) { + wellRates()[np * w + p] = rate_target * distr[p]; + } + } else { + const double small_rate = 1e-14; + const double sign = (wells[w]->wellType() == INJECTOR) ? 1.0 : -1.0; + for (int p = 0; p < np; ++p) { + wellRates()[np * w + p] = small_rate * sign; + } + } + + // 2. Bhp: + if (well_controls_get_current_type(ctrl) == BHP) { + bhp()[w] = well_controls_get_current_target(ctrl); + } else { + const int first_cell = wells[w]->wellCells()[0]; + const double safety_factor = (wells[w]->wellType() == INJECTOR) ? 1.01 : 0.99; + bhp()[w] = safety_factor* state.pressure()[first_cell]; + } + // 3. Thp: + if (well_controls_get_current_type(ctrl) == THP) { + thp()[w] = well_controls_get_current_target(ctrl); + } else { + thp()[w] = bhp()[w]; + } + + // 4. Perf rates and pressures + int number_of_perforations = wellMapEntry.number_of_perforations; + for (int i = 0; i < number_of_perforations; ++i) { + for (int p = 0; p < np; ++p) { + perfPhaseRates()[np * (i + start_perforation) + p] = wellRates()[np * w + p] / double(number_of_perforations); + } + if (wells[w]->isMultiSegmented()) { + const double safety_factor = (wells[w]->wellType() == INJECTOR) ? 1.01 : 0.99; + perfPress()[i + start_perforation] = safety_factor * state.pressure()[wells[w]->wellCells()[i]]; + } else { + perfPress()[i + start_perforation] = state.pressure()[wells[w]->wellCells()[i]]; + } + } + + // 5. Segment rates and pressures + int number_of_segments = wellMapEntry.number_of_segments; + // the seg_pressure is the same with the first perf_pressure. For the top segment, it is the same with bhp, + // when under bhp control. + // the seg_rates will related to the sum of the perforation rates, and also trying to keep consistent with the + // well rates. Most importantly, the segment rates of the top segment is the same with the well rates. + segpress_[start_segment] = bhp()[w]; + for (int i = 1; i < number_of_segments; ++i) { + int first_perforation_segment = start_perforation + wellMapEntry.start_perforation_segment[i]; + segpress_[i + start_segment] = perfPress()[first_perforation_segment]; + // the segmnent pressure of the top segment should be the bhp + } + + for (int p = 0; p < np; ++p) { + Eigen::VectorXd v_perf_rates(number_of_perforations); + for (int i = 0; i < number_of_perforations; ++i) { + v_perf_rates[i] = perfPhaseRates()[np * (i + start_perforation) + p]; + } + + Eigen::VectorXd v_segment_rates = wells[w]->wellOps().p2s_gather * v_perf_rates; + + for (int i = 0; i < number_of_segments; ++i) { + segphaserates_[np * (i + start_segment) + p] = v_segment_rates[i]; + } + } + } + + start_segment += wellMapEntry.number_of_segments; + start_perforation += wellMapEntry.number_of_perforations; + + } + + // Initialize current_controls_. + // The controls set in the Wells object are treated as defaults, + // and also used for initial values. + currentControls().resize(nw); + for (int w = 0; w < nw; ++w) { + currentControls()[w] = well_controls_get_current(wells[w]->wellControls()); + } + + // initialize wells that have been there before + // order can change so the mapping is based on the well names + if ( !(prevState.segmentedWellMap().empty()) ) + { + typedef typename SegmentedWellMapType::const_iterator const_iterator; + const_iterator end_old = prevState.segmentedWellMap().end(); + + for (int w = 0; w < nw; ++w) { + std::string well_name(wells[w]->name()); + const_iterator it_old = prevState.segmentedWellMap().find(well_name); + const_iterator it_this = segmentedWellMap().find(well_name); + + assert(it_this != segmentedWellMap().end()); // the current well must be present in the current well map + + if (it_old != end_old) { + const int oldIndex = (*it_old).second.well_number; + const int newIndex = w; + + // bhp + bhp()[newIndex] = prevState.bhp()[oldIndex]; + + // well rates + for( int i=0, idx=newIndex*np, oldidx=oldIndex*np; iwellControls())) { + // If the set of controls have changed, this may not be identical + // to the last control, but it must be a valid control. + currentControls()[ newIndex ] = old_control_index; + } + } + } + } + } + } + + + std::vector& segPress() { return segpress_; } + const std::vector& segPress() const { return segpress_; } + + std::vector& segPhaseRates() { return segphaserates_; } + const std::vector& segPhaseRates() const { return segphaserates_; } + + const std::vector& topSegmentLoc() const { return top_segment_loc_; }; + + const SegmentedWellMapType& segmentedWellMap() const { return segmentedWellMap_; } + SegmentedWellMapType& segmentedWellMap() { return segmentedWellMap_; } + + int numSegments() const { return nseg_; } + int numPerforations() const { return nperf_; } + + private: + // pressure for the segment nodes + std::vector segpress_; + // phase rates for the segments + std::vector segphaserates_; + + // the location of the top segments within the whole segment list + // it is better in the Wells class if we have a class instead of + // using a vector for all the wells + std::vector top_segment_loc_; + + SegmentedWellMapType segmentedWellMap_; + + int nseg_; + int nperf_; + }; + +} // namespace Opm + + +#endif // OPM_WELLSTATEMULTISEGMENT_HEADER_INCLUDE