OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-12-28 18:21:00 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #541 from GitPaean/support_segment_well_rebased

Browse Source 
Support multi-segment wells
This commit is contained in:
dr-robertk 2015-12-02 11:18:33 -07:00
commit bd8586b477
13 changed files with 3694 additions and 54 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
CMakeLists.txt
View File

@ -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)

8
CMakeLists_files.cmake
View File

@ -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
)

463
examples/flow_multisegment.cpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#if HAVE_CONFIG_H
#include "config.h"
#endif // HAVE_CONFIG_H
#include <dune/common/version.hh>
#include <opm/common/utility/platform_dependent/disable_warnings.h>
#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 3)
#include <dune/common/parallel/mpihelper.hh>
#else
#include <dune/common/mpihelper.hh>
#endif
#if HAVE_DUNE_CORNERPOINT && WANT_DUNE_CORNERPOINTGRID
#define USE_DUNE_CORNERPOINTGRID 1
#include <dune/grid/CpGrid.hpp>
#include <dune/grid/common/GridAdapter.hpp>
#else
#undef USE_DUNE_CORNERPOINTGRID
#endif
#include <opm/common/utility/platform_dependent/reenable_warnings.h>
#include <opm/core/pressure/FlowBCManager.hpp>
#include <opm/core/grid.h>
#include <opm/core/grid/cornerpoint_grid.h>
#include <opm/core/grid/GridManager.hpp>
#include <opm/autodiff/GridHelpers.hpp>
#include <opm/autodiff/createGlobalCellArray.hpp>
#include <opm/core/wells.h>
#include <opm/core/wells/WellsManager.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/core/simulator/initState.hpp>
#include <opm/core/simulator/initStateEquil.hpp>
#include <opm/core/simulator/SimulatorReport.hpp>
#include <opm/core/simulator/SimulatorTimer.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/utility/thresholdPressures.hpp> // Note: the GridHelpers must be included before this (to make overloads available). \TODO: Fix.
#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
#include <opm/core/props/BlackoilPropertiesBasic.hpp>
#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
#include <opm/core/props/rock/RockCompressibility.hpp>
#include <opm/core/linalg/LinearSolverFactory.hpp>
#include <opm/autodiff/NewtonIterationBlackoilSimple.hpp>
#include <opm/autodiff/NewtonIterationBlackoilCPR.hpp>
#include <opm/autodiff/NewtonIterationBlackoilInterleaved.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/autodiff/SimulatorFullyImplicitBlackoilMultiSegment.hpp>
#include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
#include <opm/autodiff/RedistributeDataHandles.hpp>
#include <opm/autodiff/moduleVersion.hpp>
#include <opm/core/utility/share_obj.hpp>
#include <opm/parser/eclipse/OpmLog/OpmLog.hpp>
#include <opm/parser/eclipse/OpmLog/StreamLog.hpp>
#include <opm/parser/eclipse/OpmLog/CounterLog.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/Parser/ParseMode.hpp>
#include <opm/parser/eclipse/EclipseState/checkDeck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <boost/filesystem.hpp>
#include <boost/algorithm/string.hpp>
#ifdef _OPENMP
#include <omp.h>
#endif
#include <memory>
#include <algorithm>
#include <cstdlib>
#include <iostream>
#include <vector>
#include <numeric>
#include <cstdlib>
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=<path to your deck>, 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<std::string>("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<Opm::StreamLog> streamLog = std::make_shared<Opm::StreamLog>(logFile , Opm::Log::DefaultMessageTypes);
std::shared_ptr<Opm::CounterLog> counterLog = std::make_shared<Opm::CounterLog>(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> 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<double> 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<int>("output_interval");
IOConfigPtr ioConfig = eclipseState->getIOConfig();
ioConfig->overrideRestartWriteInterval((size_t)output_interval);
}
const PhaseUsage pu = Opm::phaseUsageFromDeck(deck);
std::vector<int> compressedToCartesianIdx;
Opm::createGlobalCellArray(grid, compressedToCartesianIdx);
typedef BlackoilPropsAdFromDeck::MaterialLawManager MaterialLawManager;
auto materialLawManager = std::make_shared<MaterialLawManager>();
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<int> cells(numCells);
for (int c = 0; c < numCells; ++c) { cells[c] = c; }
std::vector<double> 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<NewtonIterationBlackoilInterface> 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<const Opm::SimulationConfig> simCfg = eclipseState->getSimulationConfig();
std::string solver_approach = flowDefaultSolver;
if (param.has("solver_approach")) {
solver_approach = param.get<std::string>("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<std::pair<int, int>, double> maxDp;
computeMaxDp(maxDp, deck, eclipseState, grid, state, props, gravity[2]);
std::vector<double> 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;
}

2
examples/flow_multisegment_mpi.cpp Normal file
View File

@ -0,0 +1,2 @@
#define WANT_DUNE_CORNERPOINTGRID 1
#include "flow_multisegment.cpp"

17
opm/autodiff/BlackoilModelBase.hpp
View File

@ -273,6 +273,7 @@ namespace Opm {
WellOps(const Wells* wells);
Eigen::SparseMatrix<double> w2p; // well -> perf (scatter)
Eigen::SparseMatrix<double> p2w; // perf -> well (gather)
std::vector<int> 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<ADB>& mob_perfcells,
std::vector<ADB>& b_perfcells) const;
bool
solveWellEq(const std::vector<ADB>& mob_perfcells,
const std::vector<ADB>& b_perfcells,
SolutionState& state,
@ -398,12 +405,12 @@ namespace Opm {
const std::vector<ADB>& mob_perfcells,
const std::vector<ADB>& b_perfcells,
V& aliveWells,
std::vector<ADB>& cq_s);
std::vector<ADB>& cq_s) const;
void
updatePerfPhaseRatesAndPressures(const std::vector<ADB>& cq_s,
const SolutionState& state,
WellState& xw);
WellState& xw) const;
void
addWellFluxEq(const std::vector<ADB>& 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<double, Eigen::Dynamic, Eigen::Dynamic>& B,
@ -546,8 +552,7 @@ namespace Opm {
std::vector<double>& maxCoeff,
std::vector<double>& B_avg,
std::vector<double>& maxNormWell,
int nc,
int nw) const;
int nc) const;
double dpMaxRel() const { return param_.dp_max_rel_; }
double dsMax() const { return param_.ds_max_; }

137
opm/autodiff/BlackoilModelBase_impl.hpp
View File

@ -407,7 +407,8 @@ namespace detail {
BlackoilModelBase<Grid, Implementation>::
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 <class Grid, class Implementation>
int
BlackoilModelBase<Grid, Implementation>::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 <class Grid, class Implementation>
void
BlackoilModelBase<Grid, Implementation>::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<double> 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<int> well_cells(wells().well_cells, wells().well_cells + nperf);
// Compute the average pressure in each well block
const V perf_press = Eigen::Map<const V>(xw.perfPress().data(), nperf);
@ -791,6 +806,8 @@ namespace detail {
}
}
const std::vector<int>& 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<ADB> cq_s(np, ADB::null());
const int nw = wells().number_of_wells;
const int nperf = wells().well_connpos[nw];
const std::vector<int> well_cells(wells().well_cells, wells().well_cells + nperf);
std::vector<ADB> mob_perfcells(np, ADB::null());
std::vector<ADB> 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<ADB> mob_perfcells;
std::vector<ADB> 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<ADB> 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<int> 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 <class Grid, class Implementation>
void
BlackoilModelBase<Grid, Implementation>::extractWellPerfProperties(std::vector<ADB>& mob_perfcells,
std::vector<ADB>& 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<int>& 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 <class Grid, class Implementation>
void
BlackoilModelBase<Grid, Implementation>::computeWellFlux(const SolutionState& state,
const std::vector<ADB>& mob_perfcells,
const std::vector<ADB>& b_perfcells,
V& aliveWells,
std::vector<ADB>& cq_s)
std::vector<ADB>& 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<const V>(wells().WI, nperf);
const std::vector<int> well_cells(wells().well_cells, wells().well_cells + nperf);
const std::vector<int>& 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 <class Grid, class Implementation>
void BlackoilModelBase<Grid, Implementation>::updatePerfPhaseRatesAndPressures(const std::vector<ADB>& 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 <class Grid, class Implementation>
void BlackoilModelBase<Grid, Implementation>::solveWellEq(const std::vector<ADB>& mob_perfcells,
bool BlackoilModelBase<Grid, Implementation>::solveWellEq(const std::vector<ADB>& mob_perfcells,
const std::vector<ADB>& b_perfcells,
SolutionState& state,
WellState& well_state)
@ -1563,6 +1598,7 @@ namespace detail {
std::vector<ADB> cq_s(np, ADB::null());
std::vector<int> indices = variableWellStateIndices();
SolutionState state0 = state;
WellState well_state0 = well_state;
asImpl().makeConstantState(state0);
std::vector<ADB> mob_perfcells_const(np, ADB::null());
@ -1578,15 +1614,15 @@ namespace detail {
// bhp and Q for the wells
std::vector<V> vars0;
vars0.reserve(2);
variableWellStateInitials(well_state, vars0);
asImpl().variableWellStateInitials(well_state, vars0);
std::vector<ADB> 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<double>& maxCoeff,
std::vector<double>& B_avg,
std::vector<double>& 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<double> CNV(nm);
std::vector<double> 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<double> well_flux_residual(np);
bool converged_Well = true;

312
opm/autodiff/BlackoilMultiSegmentModel.hpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_BLACKOILMULTISEGMENTMODEL_HEADER_INCLUDED
#define OPM_BLACKOILMULTISEGMENTMODEL_HEADER_INCLUDED
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/autodiff/BlackoilModelBase.hpp>
#include <opm/autodiff/BlackoilModelParameters.hpp>
#include <opm/autodiff/WellStateMultiSegment.hpp>
#include <opm/autodiff/WellMultiSegment.hpp>
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 Grid>
class BlackoilMultiSegmentModel : public BlackoilModelBase<Grid, BlackoilMultiSegmentModel<Grid>>
{
public:
typedef BlackoilModelBase<Grid, BlackoilMultiSegmentModel<Grid> > 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<WellMultiSegmentConstPtr>& 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<V> segment_comp_surf_volume_initial_;
// the value within the current iteration.
std::vector<ADB> 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<WellMultiSegmentConstPtr> wells_multisegment_;
std::vector<int> 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<WellMultiSegmentConstPtr>& wells_ms);
Eigen::SparseMatrix<double> w2p; // well -> perf (scatter)
Eigen::SparseMatrix<double> p2w; // perf -> well (gather)
Eigen::SparseMatrix<double> w2s; // well -> segment (scatter)
Eigen::SparseMatrix<double> s2w; // segment -> well (gather)
Eigen::SparseMatrix<double> s2p; // segment -> perf (scatter)
Eigen::SparseMatrix<double> p2s; // perf -> segment (gather)
Eigen::SparseMatrix<double> s2s_inlets; // segment -> its inlet segments
Eigen::SparseMatrix<double> s2s_outlet; // segment -> its outlet segment
Eigen::SparseMatrix<double> topseg2w; // top segment -> well
AutoDiffMatrix eliminate_topseg; // change the top segment related to be zero
std::vector<int> 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<WellMultiSegmentConstPtr>& 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<V>& vars0) const;
void computeWellConnectionPressures(const SolutionState& state,
const WellState& xw);
bool
solveWellEq(const std::vector<ADB>& mob_perfcells,
const std::vector<ADB>& b_perfcells,
SolutionState& state,
WellState& well_state);
void
computeWellFlux(const SolutionState& state,
const std::vector<ADB>& mob_perfcells,
const std::vector<ADB>& b_perfcells,
V& aliveWells,
std::vector<ADB>& cq_s) const;
void
updatePerfPhaseRatesAndPressures(const std::vector<ADB>& cq_s,
const SolutionState& state,
WellState& xw) const;
void
addWellFluxEq(const std::vector<ADB>& 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<int>& indices,
std::vector<ADB>& 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 <class GridT>
struct ModelTraits< BlackoilMultiSegmentModel<GridT> >
{
typedef BlackoilState ReservoirState;
typedef WellStateMultiSegment WellState;
typedef BlackoilModelParameters ModelParameters;
typedef BlackoilMultiSegmentSolutionState SolutionState;
};
} // namespace Opm
#include "BlackoilMultiSegmentModel_impl.hpp"
#endif // OPM_BLACKOILMULTISEGMENTMODEL_HEADER_INCLUDED

1505
opm/autodiff/BlackoilMultiSegmentModel_impl.hpp Normal file
View File

File diff suppressed because it is too large Load Diff

104
opm/autodiff/SimulatorFullyImplicitBlackoilMultiSegment.hpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_SIMULATORFULLYIMPLICITBLACKOILMULTISEGMENT_HEADER_INCLUDED
#define OPM_SIMULATORFULLYIMPLICITBLACKOILMULTISEGMENT_HEADER_INCLUDED
#include <opm/autodiff/SimulatorBase.hpp>
#include <opm/autodiff/NonlinearSolver.hpp>
#include <opm/autodiff/BlackoilMultiSegmentModel.hpp>
#include <opm/autodiff/WellStateMultiSegment.hpp>
namespace Opm {
template <class GridT>
class SimulatorFullyImplicitBlackoilMultiSegment;
template <class GridT>
struct SimulatorTraits<SimulatorFullyImplicitBlackoilMultiSegment<GridT> >
{
typedef WellStateMultiSegment WellState;
typedef BlackoilState ReservoirState;
typedef BlackoilOutputWriter OutputWriter;
typedef GridT Grid;
typedef BlackoilMultiSegmentModel<Grid> Model;
typedef NonlinearSolver<Model> Solver;
};
/// a simulator for the blackoil model
template <class GridT>
class SimulatorFullyImplicitBlackoilMultiSegment
: public SimulatorBase<SimulatorFullyImplicitBlackoilMultiSegment<GridT> >
{
typedef SimulatorBase<SimulatorFullyImplicitBlackoilMultiSegment<GridT> > Base;
typedef SimulatorFullyImplicitBlackoilMultiSegment<GridT> ThisType;
typedef SimulatorTraits<ThisType> 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<EclipseState> eclipse_state,
BlackoilOutputWriter& output_writer,
const std::vector<double>& 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<Solver> createSolver(const Wells* wells, std::vector<WellMultiSegmentConstPtr>& 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

213
opm/autodiff/SimulatorFullyImplicitBlackoilMultiSegment_impl.hpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
namespace Opm
{
template <class GridT>
auto SimulatorFullyImplicitBlackoilMultiSegment<GridT>::
createSolver(const Wells* wells, std::vector<WellMultiSegmentConstPtr>& wells_multisegment)
-> std::unique_ptr<Solver>
{
typedef typename Traits::Model Model;
auto model = std::unique_ptr<Model>(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<ThisType::Solver>(new Solver(Base::solver_param_, std::move(model)));
}
template <class GridT>
SimulatorReport SimulatorFullyImplicitBlackoilMultiSegment<GridT>::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<WellConstPtr>& wells_ecl = eclipse_state_->getSchedule()->getWells(timer.currentStepNum());
std::vector<WellMultiSegmentConstPtr> 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<WellMultiSegment>(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

400
opm/autodiff/WellMultiSegment.cpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#include <opm/autodiff/WellMultiSegment.hpp>
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<int> temp_well_cell;
std::vector<double> 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<double> 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<double> Tri;
std::vector<Tri> s2p;
std::vector<Tri> 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<Tri> s2s_gather;
s2s_gather.reserve(m_number_of_segments_ * m_number_of_segments_);
std::vector<Tri> s2s_inlets;
s2s_inlets.reserve(m_number_of_segments_);
std::vector<Tri> 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<Tri> 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<double>& WellMultiSegment::compFrac() const {
return m_comp_frac_;
}
int WellMultiSegment::numberOfPhases() const {
return m_number_of_phases_;
}
const std::vector<double>& WellMultiSegment::wellIndex() const {
return m_well_index_;
}
const std::vector<double>& WellMultiSegment::perfDepth() const {
return m_perf_depth_;
}
const std::vector<int>& WellMultiSegment::wellCells() const {
return m_well_cell_;
}
const std::vector<int>& WellMultiSegment::segmentCells() const {
return m_segment_cell_;
}
const std::vector<int>& WellMultiSegment::outletSegment() const {
return m_outlet_segment_;
}
const std::vector<std::vector<int>>& WellMultiSegment::inletSegments() const {
return m_inlet_segments_;
}
const std::vector<double>& WellMultiSegment::segmentLength() const {
return m_segment_length_;
}
const std::vector<double>& WellMultiSegment::segmentDepth() const {
return m_segment_depth_;
}
const std::vector<double>& WellMultiSegment::segmentDiameter() const {
return m_segment_internal_diameter_;
}
const std::vector<double>& WellMultiSegment::segmentCrossArea() const {
return m_segment_cross_area_;
}
const std::vector<double>& WellMultiSegment::segmentRoughness() const {
return m_segment_roughness_;
}
const std::vector<double>& WellMultiSegment::segmentVolume() const {
return m_segment_volume_;
}
const std::vector<std::vector<int>>& WellMultiSegment::segmentPerforations() const {
return m_segment_perforations_;
}
const WellMultiSegment::WellOps& WellMultiSegment::wellOps() const {
return m_wops_;
}
}

227
opm/autodiff/WellMultiSegment.hpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_WELLMULTISEGMENT_HEADER_INCLUDED
#define OPM_WELLMULTISEGMENT_HEADER_INCLUDED
#include <opm/common/utility/platform_dependent/disable_warnings.h>
#include <Eigen/Eigen>
#include <Eigen/Sparse>
#include <opm/common/utility/platform_dependent/reenable_warnings.h>
#include <opm/core/wells.h>
#include <opm/core/well_controls.h>
#include <opm/core/simulator/WellState.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/ScheduleEnums.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/MSW/SegmentSet.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
#include <vector>
#include <cassert>
#include <string>
#include <utility>
#include <map>
#include <algorithm>
#include <array>
namespace Opm
{
class WellMultiSegment
{
public:
typedef Eigen::SparseMatrix<double> 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<double>& compFrac() const;
/// Number of phases.
int numberOfPhases() const;
/// Well type.
WellType wellType() const;
/// Well productivity index.
const std::vector<double>& wellIndex() const;
/// Depth of the perforations.
const std::vector<double>& perfDepth() const;
/// Indices of the grid blocks that perforations are completed in.
const std::vector<int>& wellCells() const;
/// Indices of the gird blocks that segments locate at.
const std::vector<int>& 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<int>& outletSegment() const;
/// Inlet segments. a segment can have more than one inlet segments.
const std::vector<std::vector<int>>& inletSegments() const;
/// The length of the segment nodes down the wellbore from the reference point.
const std::vector<double>& segmentLength() const;
/// The depth of the segment nodes.
const std::vector<double>& segmentDepth() const;
/// Tubing internal diameter.
const std::vector<double>& segmentDiameter() const;
/// Cross-sectional area.
const std::vector<double>& segmentCrossArea() const;
/// Effective absolute roughness of the tube.
const std::vector<double>& segmentRoughness() const;
/// Volume of segment.
const std::vector<double>& segmentVolume() const;
/// Perforations related to each segment.
const std::vector<std::vector<int>>& 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<double> 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<double> m_well_index_;
// depth for each completion // form the keyword COMPSEGS?
std::vector<double> m_perf_depth_;
// well cell for each completion
std::vector<int> m_well_cell_;
// cell for each segment
std::vector<int> 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<int> m_outlet_segment_;
// for convinience, we store the inlet segments for each segment
std::vector<std::vector<int>> m_inlet_segments_;
// this one is not necessary any more, since the segment number is its location.
// std::map<int, int> 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<double> m_segment_length_;
// the depth of the segmnet node
std::vector<double> m_segment_depth_;
// the internal diameter of the segment
std::vector<double> m_segment_internal_diameter_;
// the roughness of the segment
std::vector<double> m_segment_roughness_;
// the cross area of the segment
std::vector<double> m_segment_cross_area_;
// the volume of the segment
std::vector<double> 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<std::vector<int>> m_segment_perforations_;
WellOps m_wops_;
};
typedef std::shared_ptr<WellMultiSegment> WellMultiSegmentPtr;
typedef std::shared_ptr<const WellMultiSegment> WellMultiSegmentConstPtr;
} // namespace Opm
#endif // OPM_WELLMULTISEGMENT_HEADER_INCLUDE

359
opm/autodiff/WellStateMultiSegment.hpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_WELLSTATEMULTISEGMENT_HEADER_INCLUDED
#define OPM_WELLSTATEMULTISEGMENT_HEADER_INCLUDED
#include <opm/core/wells.h>
#include <opm/core/well_controls.h>
#include <opm/common/ErrorMacros.hpp>
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/autodiff/WellMultiSegment.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <vector>
#include <cassert>
#include <string>
#include <utility>
#include <map>
#include <algorithm>
#include <array>
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<int> start_perforation_segment; // the starting position of perforation inside the segment
std::vector<int> number_of_perforations_segment; // the numbers for perforations for the segments
} SegmentedMapentryType;
typedef std::map<std::string, SegmentedMapentryType> 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 <class ReservoirState, class PrevWellState>
void init(const std::vector<WellMultiSegmentConstPtr>& 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; i<np; ++i, ++idx, ++oldidx )
{
wellRates()[ idx ] = prevState.wellRates()[ oldidx ];
}
const int num_seg_old_well = (*it_old).second.number_of_segments;
const int num_perf_old_well = (*it_old).second.number_of_perforations;
const int num_seg_this_well = (*it_this).second.number_of_segments;
const int num_perf_this_well = (*it_this).second.number_of_perforations;
// determing if the structure of the wells has been changed by comparing the number of segments and perforations
// may not be very safe.
// The strategy HAS to be changed later with experiments and analysis
if ((num_perf_old_well == num_perf_this_well) && (num_seg_old_well == num_seg_this_well)) {
const int old_start_perforation = (*it_old).second.start_perforation;
const int old_start_segment = (*it_old).second.start_segment;
const int this_start_perforation = (*it_this).second.start_perforation;
const int this_start_segment = (*it_this).second.start_segment;
// this is not good when the the well rates changed dramatically
for (int i = 0; i < num_seg_this_well * np; ++i) {
segphaserates_[this_start_segment * np + i] = prevState.segPhaseRates()[old_start_segment * np + i];
}
for (int i = 0; i < num_perf_this_well * np; ++i) {
perfPhaseRates()[this_start_perforation * np + i] = prevState.perfPhaseRates()[old_start_perforation * np + i];
}
// perf_pressures_
for (int i = 0; i < num_perf_this_well; ++i) {
// p
perfPress()[this_start_perforation + i] = prevState.perfPress()[old_start_perforation + i];
}
// segpress_
for (int i = 0; i < num_seg_this_well; ++i) {
// p
segpress_[this_start_segment + i] = prevState.segPress()[old_start_segment + i];
}
// current controls
const int old_control_index = prevState.currentControls()[ oldIndex ];
if (old_control_index < well_controls_get_num(wells[w]->wellControls())) {
// 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<double>& segPress() { return segpress_; }
const std::vector<double>& segPress() const { return segpress_; }
std::vector<double>& segPhaseRates() { return segphaserates_; }
const std::vector<double>& segPhaseRates() const { return segphaserates_; }
const std::vector<int>& 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<double> segpress_;
// phase rates for the segments
std::vector<double> 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<int> top_segment_loc_;
SegmentedWellMapType segmentedWellMap_;
int nseg_;
int nperf_;
};
} // namespace Opm
#endif // OPM_WELLSTATEMULTISEGMENT_HEADER_INCLUDE