OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Examples and tutorials follow change to IncompTpfa interface.

Browse Source
This commit is contained in:
Atgeirr Flø Rasmussen
2012-06-12 15:28:53 +02:00
parent e1d5e55f1b
commit 92f1ce9b83
5 changed files with 118 additions and 281 deletions
Download Patch File Download Diff File Expand all files Collapse all files

84
examples/spu_2p.cpp
View File

@@ -71,6 +71,7 @@
#include <opm/core/utility/ColumnExtract.hpp>
#include <opm/core/simulator/TwophaseState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/transport/GravityColumnSolver.hpp>
#include <opm/core/transport/reorder/TransportModelTwophase.hpp>
@@ -379,13 +380,15 @@ main(int argc, char** argv)
// Check that rock compressibility is not used with solvers that do not handle it.
int nl_pressure_maxiter = 0;
double nl_pressure_tolerance = 0.0;
double nl_pressure_residual_tolerance = 0.0;
double nl_pressure_change_tolerance = 0.0;
if (rock_comp->isActive()) {
if (!use_reorder) {
THROW("Cannot run implicit (non-reordering) transport solver with rock compressibility yet.");
}
nl_pressure_residual_tolerance = param.getDefault("nl_pressure_residual_tolerance", 1e-8);
nl_pressure_change_tolerance = param.getDefault("nl_pressure_change_tolerance", 1.0); // In Pascal.
nl_pressure_maxiter = param.getDefault("nl_pressure_maxiter", 10);
nl_pressure_tolerance = param.getDefault("nl_pressure_tolerance", 1.0); // in Pascal
}
// Source-related variables init.
@@ -442,7 +445,10 @@ main(int argc, char** argv)
Opm::LinearSolverFactory linsolver(param);
// Pressure solver.
const double *grav = use_gravity ? &gravity[0] : 0;
Opm::IncompTpfa psolver(*grid->c_grid(), props->permeability(), grav, linsolver, wells->c_wells());
Opm::IncompTpfa psolver(*grid->c_grid(), *props, rock_comp.get(), linsolver,
nl_pressure_residual_tolerance, nl_pressure_change_tolerance,
nl_pressure_maxiter,
grav, wells->c_wells(), src, bcs.c_bcs());
// Reordering solver.
const double nl_tolerance = param.getDefault("nl_tolerance", 1e-9);
const int nl_maxiter = param.getDefault("nl_maxiter", 30);
@@ -515,17 +521,15 @@ main(int argc, char** argv)
Opm::Watercut watercut;
watercut.push(0.0, 0.0, 0.0);
Opm::WellReport wellreport;
std::vector<double> well_bhp;
std::vector<double> well_perfrates;
Opm::WellState well_state;
well_state.init(wells->c_wells(), state);
std::vector<double> fractional_flows;
std::vector<double> well_resflows_phase;
int num_wells = 0;
if (wells->c_wells()) {
num_wells = wells->c_wells()->number_of_wells;
well_bhp.resize(num_wells, 0.0);
well_perfrates.resize(wells->c_wells()->well_connpos[num_wells], 0.0);
well_resflows_phase.resize((wells->c_wells()->number_of_phases)*(wells->c_wells()->number_of_wells), 0.0);
wellreport.push(*props, *wells->c_wells(), state.saturation(), 0.0, well_bhp, well_perfrates);
wellreport.push(*props, *wells->c_wells(), state.saturation(), 0.0, well_state.bhp(), well_state.perfRates());
}
for (; !simtimer.done(); ++simtimer) {
// Report timestep and (optionally) write state to disk.
@@ -540,10 +544,6 @@ main(int argc, char** argv)
} else {
computeTotalMobility(*props, allcells, state.saturation(), totmob);
}
std::vector<double> wdp;
if (wells->c_wells()) {
Opm::computeWDP(*wells->c_wells(), *grid->c_grid(), state.saturation(), props->density(), gravity[2], true, wdp);
}
if (check_well_controls) {
computeFractionalFlow(*props, allcells, state.saturation(), fractional_flows);
}
@@ -554,52 +554,9 @@ main(int argc, char** argv)
int well_control_iteration = 0;
do {
pressure_timer.start();
if (rock_comp->isActive()) {
rc.resize(num_cells);
std::vector<double> initial_pressure = state.pressure();
std::vector<double> initial_porevolume(num_cells);
computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, initial_pressure, initial_porevolume);
std::vector<double> pressure_increment(num_cells + num_wells);
std::vector<double> prev_pressure(num_cells + num_wells);
for (int iter = 0; iter < nl_pressure_maxiter; ++iter) {
for (int cell = 0; cell < num_cells; ++cell) {
rc[cell] = rock_comp->rockComp(state.pressure()[cell]);
}
computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol);
// prev_pressure = state.pressure();
std::copy(state.pressure().begin(), state.pressure().end(), prev_pressure.begin());
std::copy(well_bhp.begin(), well_bhp.end(), prev_pressure.begin() + num_cells);
// compute pressure increment
psolver.solveIncrement(totmob, omega, src, wdp, bcs.c_bcs(), porevol, rc,
prev_pressure, initial_porevolume, simtimer.currentStepLength(),
pressure_increment);
double max_change = 0.0;
for (int cell = 0; cell < num_cells; ++cell) {
state.pressure()[cell] += pressure_increment[cell];
max_change = std::max(max_change, std::fabs(pressure_increment[cell]));
}
for (int well = 0; well < num_wells; ++well) {
well_bhp[well] += pressure_increment[num_cells + well];
max_change = std::max(max_change, std::fabs(pressure_increment[num_cells + well]));
}
std::cout << "Pressure iter " << iter << " max change = " << max_change << std::endl;
if (max_change < nl_pressure_tolerance) {
break;
}
}
psolver.computeFaceFlux(totmob, omega, src, wdp, bcs.c_bcs(), state.pressure(), state.faceflux(),
well_bhp, well_perfrates);
} else {
std::vector<double> initial_pressure = state.pressure();
psolver.solve(totmob, omega, src, wdp, bcs.c_bcs(), state.pressure(), state.faceflux(),
well_bhp, well_perfrates);
std::vector<double> initial_pressure = state.pressure();
psolver.solve(simtimer.currentStepLength(), state, well_state);
if (!rock_comp->isActive()) {
// Compute average pressures of previous and last
// step, and total volume.
double av_prev_press = 0.;
@@ -610,14 +567,13 @@ main(int argc, char** argv)
av_press += state.pressure()[cell]*grid->c_grid()->cell_volumes[cell];
tot_vol += grid->c_grid()->cell_volumes[cell];
}
// Renormalization constant
const double ren_const = (av_prev_press - av_press)/tot_vol;
for (int cell = 0; cell < num_cells; ++cell) {
state.pressure()[cell] += ren_const;
}
for (int well = 0; well < num_wells; ++well) {
well_bhp[well] += ren_const;
well_state.bhp()[well] += ren_const;
}
}
pressure_timer.stop();
@@ -629,11 +585,11 @@ main(int argc, char** argv)
if (check_well_controls) {
Opm::computePhaseFlowRatesPerWell(*wells->c_wells(),
fractional_flows,
well_perfrates,
well_state.perfRates(),
well_resflows_phase);
std::cout << "Checking well conditions." << std::endl;
// For testing we set surface := reservoir
well_control_passed = wells->conditionsMet(well_bhp, well_resflows_phase, well_resflows_phase);
well_control_passed = wells->conditionsMet(well_state.bhp(), well_resflows_phase, well_resflows_phase);
++well_control_iteration;
if (!well_control_passed && well_control_iteration > max_well_control_iterations) {
THROW("Could not satisfy well conditions in " << max_well_control_iterations << " tries.");
@@ -648,7 +604,7 @@ main(int argc, char** argv)
// Process transport sources (to include bdy terms and well flows).
Opm::computeTransportSource(*grid->c_grid(), src, state.faceflux(), 1.0,
wells->c_wells(), well_perfrates, reorder_src);
wells->c_wells(), well_state.perfRates(), reorder_src);
if (!use_reorder) {
clear_transport_source(tsrc);
for (int cell = 0; cell < num_cells; ++cell) {
@@ -737,7 +693,7 @@ main(int argc, char** argv)
if (wells->c_wells()) {
wellreport.push(*props, *wells->c_wells(), state.saturation(),
simtimer.currentTime() + simtimer.currentStepLength(),
well_bhp, well_perfrates);
well_state.bhp(), well_state.perfRates());
}
}
total_timer.stop();

112
examples/wells_example.cpp
View File

@@ -13,6 +13,7 @@
#include <opm/core/grid.h>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/simulator/TwophaseState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/pressure/FlowBCManager.hpp>
#include <opm/core/linalg/LinearSolverFactory.hpp>
#include <opm/core/fluid/RockCompressibility.hpp>
@@ -45,10 +46,22 @@ int main(int argc, char** argv)
RockCompressibility rock_comp(parser);
Opm::LinearSolverFactory linsolver(parameters);
double nl_pressure_residual_tolerance = 1e-8;
double nl_pressure_change_tolerance = 0.0;
int nl_pressure_maxiter = 100;
if (rock_comp.isActive()) {
nl_pressure_residual_tolerance = parameters.getDefault("nl_pressure_residual_tolerance", 1e-8);
nl_pressure_change_tolerance = parameters.getDefault("nl_pressure_change_tolerance", 1.0); // in Pascal
nl_pressure_maxiter = parameters.getDefault("nl_pressure_maxiter", 10);
}
std::vector<double> src;
Opm::FlowBCManager bcs;
// EXPERIMENT_ISTL
IncompTpfa pressure_solver(*grid.c_grid(), incomp_properties.permeability(),
gravity, linsolver, wells.c_wells());
IncompTpfa pressure_solver(*grid.c_grid(), incomp_properties, &rock_comp, linsolver,
nl_pressure_residual_tolerance, nl_pressure_change_tolerance, nl_pressure_maxiter,
gravity, wells.c_wells(), src, bcs.c_bcs());
std::vector<int> all_cells;
@@ -60,67 +73,10 @@ int main(int argc, char** argv)
initStateFromDeck(*grid.c_grid(), incomp_properties, parser, gravity[2], state);
// Compute phase mobilities
std::vector<double> phase_mob;
computePhaseMobilities(incomp_properties, all_cells, state.saturation(), phase_mob);
// Compute total mobility and omega
std::vector<double> totmob;
std::vector<double> omega;
computeTotalMobilityOmega(incomp_properties, all_cells, state.saturation(), totmob, omega);
Opm::WellState well_state;
well_state.init(wells.c_wells(), state);
std::vector<double> wdp;
computeWDP(*wells.c_wells(), *grid.c_grid(), state.saturation(), incomp_properties.density(), gravity[2], true, wdp);
std::vector<double> src;
Opm::FlowBCManager bcs;
std::vector<double> pressure;
std::vector<double> face_flux;
std::vector<double> well_bhp;
std::vector<double> well_rate_per_cell;
std::vector<double> rc;
rc.resize(grid.c_grid()->number_of_cells);
int nl_pressure_maxiter = 100;
double nl_pressure_tolerance = 0.0;
if (rock_comp.isActive()) {
nl_pressure_maxiter = parameters.getDefault("nl_pressure_maxiter", 10);
nl_pressure_tolerance = parameters.getDefault("nl_pressure_tolerance", 1.0); // in Pascal
}
const int num_cells = grid.c_grid()->number_of_cells;
std::vector<double> porevol;
if (rock_comp.isActive()) {
computePorevolume(*grid.c_grid(), incomp_properties.porosity(), rock_comp, state.pressure(), porevol);
} else {
computePorevolume(*grid.c_grid(), incomp_properties.porosity(), porevol);
}
if (rock_comp.isActive()) {
std::vector<double> initial_pressure = state.pressure();
std::vector<double> prev_pressure;
for (int iter = 0; iter < nl_pressure_maxiter; ++iter) {
prev_pressure = state.pressure();
for (int cell = 0; cell < num_cells; ++cell) {
rc[cell] = rock_comp.rockComp(state.pressure()[cell]);
}
state.pressure() = initial_pressure;
pressure_solver.solve(totmob, omega, src, wdp, bcs.c_bcs(), porevol, rc, simtimer.currentStepLength(),
state.pressure(), state.faceflux(), well_bhp, well_rate_per_cell);
double max_change = 0.0;
for (int cell = 0; cell < num_cells; ++cell) {
max_change = std::max(max_change, std::fabs(state.pressure()[cell] - prev_pressure[cell]));
}
std::cout << "Pressure iter " << iter << " max change = " << max_change << std::endl;
if (max_change < nl_pressure_tolerance) {
break;
}
}
computePorevolume(*grid.c_grid(), incomp_properties.porosity(), rock_comp, state.pressure(), porevol);
} else {
pressure_solver.solve(totmob, omega, src, wdp, bcs.c_bcs(), state.pressure(), state.faceflux(),
well_bhp, well_rate_per_cell);
}
pressure_solver.solve(simtimer.currentStepLength(), state, well_state);
const int np = incomp_properties.numPhases();
std::vector<double> fractional_flows(grid.c_grid()->number_of_cells*np, 0.0);
@@ -128,38 +84,14 @@ int main(int argc, char** argv)
// This will be refactored into a separate function once done
std::vector<double> well_resflows(wells.c_wells()->number_of_wells*np, 0.0);
computePhaseFlowRatesPerWell(*wells.c_wells(), well_rate_per_cell, fractional_flows, well_resflows);
computePhaseFlowRatesPerWell(*wells.c_wells(), well_state.perfRates(), fractional_flows, well_resflows);
// We approximate (for _testing_ that resflows = surfaceflows)
for (int wc_iter = 0; wc_iter < 10 && !wells.conditionsMet(well_bhp, well_resflows, well_resflows); ++wc_iter) {
for (int wc_iter = 0; wc_iter < 10 && !wells.conditionsMet(well_state.bhp(), well_resflows, well_resflows); ++wc_iter) {
std::cout << "Conditions not met for well, trying again" << std::endl;
if (rock_comp.isActive()) {
std::vector<double> initial_pressure = state.pressure();
std::vector<double> prev_pressure;
for (int iter = 0; iter < nl_pressure_maxiter; ++iter) {
prev_pressure = state.pressure();
for (int cell = 0; cell < num_cells; ++cell) {
rc[cell] = rock_comp.rockComp(state.pressure()[cell]);
}
state.pressure() = initial_pressure;
pressure_solver.solve(totmob, omega, src, wdp, bcs.c_bcs(), porevol, rc, simtimer.currentStepLength(),
state.pressure(), state.faceflux(), well_bhp, well_rate_per_cell);
double max_change = 0.0;
for (int cell = 0; cell < num_cells; ++cell) {
max_change = std::max(max_change, std::fabs(state.pressure()[cell] - prev_pressure[cell]));
}
std::cout << "Pressure iter " << iter << " max change = " << max_change << std::endl;
if (max_change < nl_pressure_tolerance) {
break;
}
}
computePorevolume(*grid.c_grid(), incomp_properties.porosity(), rock_comp, state.pressure(), porevol);
} else {
pressure_solver.solve(totmob, omega, src, wdp, bcs.c_bcs(), state.pressure(), state.faceflux(),
well_bhp, well_rate_per_cell);
}
pressure_solver.solve(simtimer.currentStepLength(), state, well_state);
std::cout << "Solved" << std::endl;
computePhaseFlowRatesPerWell(*wells.c_wells(), well_rate_per_cell, fractional_flows, well_resflows);
computePhaseFlowRatesPerWell(*wells.c_wells(), well_state.perfRates(), fractional_flows, well_resflows);
}
#if 0

82
tutorials/tutorial2.cpp
View File

@@ -39,11 +39,14 @@
#include <iostream>
#include <fstream>
#include <vector>
#include <opm/core/fluid/IncompPropertiesBasic.hpp>
#include <opm/core/linalg/LinearSolverUmfpack.hpp>
#include <opm/core/pressure/IncompTpfa.hpp>
#include <opm/core/pressure/FlowBCManager.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/core/simulator/TwophaseState.hpp>
#include <opm/core/simulator/WellState.hpp>
/// \page tutorial2
/// \section commentedcode2 Program walkthrough.
@@ -68,16 +71,21 @@ int main()
int num_cells = grid.c_grid()->number_of_cells;
int num_faces = grid.c_grid()->number_of_faces;
/// \endcode
/// \page tutorial2
/// \details
/// We define a fluid viscosity equal to 1 cP.
/// We define a fluid viscosity equal to 1 cP and density equal
/// to 1000 kg/m^3.
/// The <opm/core/utility/Units.hpp> header contains support
/// for common units and prefixes, in the namespaces Opm::unit
/// and Opm::prefix.
/// \code
using namespace Opm::unit;
using namespace Opm::prefix;
double mu = 1.0*centi*Poise;
int num_phases = 1;
std::vector<double> mu(num_phases, 1.0*centi*Poise);
std::vector<double> rho(num_phases, 1000.0*kilogram/cubic(meter));
/// \endcode
/// \page tutorial2
/// \details
@@ -87,17 +95,13 @@ int main()
/// \endcode
/// \page tutorial2
/// \details
/// We set up a diagonal permeability tensor and compute the mobility for each cell.
/// The resulting permeability matrix is flattened in a vector.
/// \page tutorial2
/// \details
/// We set up a simple property object for a single-phase situation.
/// \code
std::vector<double> permeability(num_cells*dim*dim, 0.);
std::vector<double> mob(num_cells);
for (int cell = 0; cell < num_cells; ++cell) {
permeability[9*cell + 0] = k;
permeability[9*cell + 4] = k;
permeability[9*cell + 8] = k;
mob[cell] = 1/mu;
}
Opm::IncompPropertiesBasic props(1, Opm::SaturationPropsBasic::Constant, rho,
mu, 1.0, k, dim, num_cells);
/// \endcode
/// \page tutorial2
@@ -109,14 +113,6 @@ int main()
/// \endcode
/// \page tutorial2
/// We set up a pressure solver for the incompressible problem,
/// using the two-point flux approximation discretization. The
/// third argument which corresponds to gravity is set to a null
/// pointer (no gravity). The final argument would be a pointer to
/// a Wells data structure, again we use a null pointer to
/// indicate that we have no wells.
/// \code
Opm::IncompTpfa psolver(*grid.c_grid(), &permeability[0], 0, linsolver, 0);
/// \endcode
/// \page tutorial2
/// We define the source term.
@@ -132,39 +128,33 @@ int main()
Opm::FlowBCManager bcs;
/// \endcode
/// We set up a pressure solver for the incompressible problem,
/// using the two-point flux approximation discretization. The
/// null pointers correspond to arguments for gravity, wells and
/// boundary conditions, which are all defaulted (to zero gravity,
/// no wells, and no-flow boundaries).
/// \code
Opm::IncompTpfa psolver(*grid.c_grid(), props, linsolver, NULL, NULL, src, NULL);
/// \page tutorial2
/// We declare the solution vectors, i.e., the pressure and face
/// flux vectors we are going to compute. The well solution
/// vectors are needed for interface compatibility with the
/// We declare the state object, that will contain the pressure and face
/// flux vectors we are going to compute. The well state
/// object is needed for interface compatibility with the
/// <CODE>solve()</CODE> method of class
/// <CODE>Opm::IncompTPFA</CODE>.
/// \code
std::vector<double> pressure(num_cells);
std::vector<double> faceflux(num_faces);
std::vector<double> well_bhp;
std::vector<double> well_flux;
/// \endcode
/// \page tutorial2
/// \details
/// We declare the gravity term which is required by the pressure solver (see
/// Opm::IncompTpfa.solve()). In the absence of gravity, an empty vector is required.
/// \code
std::vector<double> omega;
/// \endcode
/// \page tutorial2
/// \details
/// We declare the wdp term which is required by the pressure solver (see
/// Opm::IncompTpfa.solve()). In the absence of wells, an empty vector is required.
/// \code
std::vector<double> wdp;
Opm::TwophaseState state;
state.pressure().resize(num_cells);
state.faceflux().resize(num_faces);
Opm::WellState well_state;
/// \endcode
/// \page tutorial2
/// We call the pressure solver.
/// The first (timestep) argument does not matter for this
/// incompressible case.
/// \code
psolver.solve(mob, omega, src, wdp, bcs.c_bcs(),
pressure, faceflux, well_bhp, well_flux);
psolver.solve(1.0*day, state, well_state);
/// \endcode
/// \page tutorial2
@@ -175,9 +165,9 @@ int main()
/// \code
std::ofstream vtkfile("tutorial2.vtu");
Opm::DataMap dm;
dm["pressure"] = &pressure;
dm["pressure"] = &state.pressure();
std::vector<double> cell_velocity;
Opm::estimateCellVelocity(*grid.c_grid(), faceflux, cell_velocity);
Opm::estimateCellVelocity(*grid.c_grid(), state.faceflux(), cell_velocity);
dm["velocity"] = &cell_velocity;
Opm::writeVtkData(*grid.c_grid(), dm, vtkfile);
}

59
tutorials/tutorial3.cpp
View File

@@ -38,6 +38,7 @@
#include <opm/core/transport/reorder/TransportModelTwophase.hpp>
#include <opm/core/simulator/TwophaseState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/Units.hpp>
@@ -161,16 +162,6 @@ int main ()
std::vector<double> omega;
/// \endcode
/// \page tutorial3
/// \details We may now set up the pressure solver. At this point,
/// unchanging parameters such as transmissibility are computed
/// and stored internally by the IncompTpfa class. The final (null pointer)
/// constructor argument is for wells, which are now used in this tutorial.
/// \code
LinearSolverUmfpack linsolver;
IncompTpfa psolver(grid, props.permeability(), grav, linsolver, 0);
/// \endcode
/// \page tutorial3
/// \details We set up the source term. Positive numbers indicate that the cell is a source,
/// while negative numbers indicate a sink.
@@ -181,14 +172,28 @@ int main ()
/// \endcode
/// \page tutorial3
/// \details We set up data vectors for the wells. Here, there are
/// no wells and we let them be empty dummies.
/// \details We set up the boundary conditions. Letting bcs be empty is equivalent
/// to no-flow boundary conditions.
/// \code
std::vector<double> empty_wdp;
std::vector<double> empty_well_bhp;
std::vector<double> empty_well_flux;
FlowBCManager bcs;
/// \endcode
/// \page tutorial3
/// \details We may now set up the pressure solver. At this point,
/// unchanging parameters such as transmissibility are computed
/// and stored internally by the IncompTpfa class. The null pointer
/// constructor argument is for wells, which are not used in this tutorial.
/// \code
LinearSolverUmfpack linsolver;
IncompTpfa psolver(grid, props, linsolver, grav, NULL, src, bcs.c_bcs());
/// \endcode
/// \page tutorial3
/// \details We set up a state object for the wells. Here, there are
/// no wells and we let it remain empty.
/// \code
WellState well_state;
/// \endcode
/// \page tutorial3
/// \details We compute the pore volume
@@ -223,14 +228,6 @@ int main ()
}
/// \endcode
/// \page tutorial3
/// \details We set up the boundary conditions. Letting bcs be empty is equivalent
/// to no-flow boundary conditions.
/// \code
FlowBCManager bcs;
/// \endcode
/// \page tutorial3
/// \details
/// We set up a two-phase state object, and
@@ -240,13 +237,6 @@ int main ()
state.init(grid, 2);
state.setFirstSat(allcells, props, TwophaseState::MinSat);
/// \endcode
/// \page tutorial3
/// \details We introduce a vector which contains the total mobility
/// on all cells.
/// \code
std::vector<double> totmob;
/// \endcode
/// \page tutorial3
/// \details This string stream will be used to construct a new
@@ -263,18 +253,11 @@ int main ()
/// \endcode
/// \page tutorial3
/// \details Compute the total mobility. It is needed by the
/// pressure solver and must be recomputed every time step
/// since it depends on the saturation.
/// \code
computeTotalMobility(props, allcells, state.saturation(), totmob);
/// \endcode
/// \page tutorial3
/// \details Solve the pressure equation
/// \code
psolver.solve(totmob, omega, src, empty_wdp, bcs.c_bcs(),
state.pressure(), state.faceflux(), empty_well_bhp,
empty_well_flux);
psolver.solve(dt, state, well_state);
/// \endcode
/// \page tutorial3
/// \details Solve the transport equation.

62
tutorials/tutorial4.cpp
View File

@@ -38,6 +38,7 @@
#include <opm/core/transport/reorder/TransportModelTwophase.hpp>
#include <opm/core/simulator/TwophaseState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/Units.hpp>
@@ -120,20 +121,6 @@ int main ()
std::vector<double> omega;
/// \endcode
/// \page tutorial4
/// \details We set up necessary information for the wells
/// \code
std::vector<double> wdp;
std::vector<double> well_bhp;
std::vector<double> well_flux;
std::vector<double> well_resflowrates_phase;
std::vector<double> well_surflowrates_phase;
std::vector<double> fractional_flows;
/// \endcode
/// \page tutorial4
/// \details We set up the source term. Positive numbers indicate that the cell is a source,
/// while negative numbers indicate a sink.
@@ -192,13 +179,6 @@ int main ()
state.setFirstSat(allcells, props, TwophaseState::MinSat);
/// \endcode
/// \page tutorial4
/// \details We introduce a vector which contains the total mobility
/// on all cells.
/// \code
std::vector<double> totmob;
/// \endcode
/// \page tutorial4
/// \details This string will contain the name of a VTK output vector.
/// \code
@@ -299,11 +279,21 @@ int main ()
///\endcode
/// \page tutorial4
/// \details We set up the pressure solver. We need to pass the wells pointer as the
/// last argument.
/// \details We set up necessary information for the wells
/// \code
WellState well_state;
well_state.init(wells, state);
std::vector<double> well_resflowrates_phase;
std::vector<double> well_surflowrates_phase;
std::vector<double> fractional_flows;
/// \endcode
/// \page tutorial4
/// \details We set up the pressure solver.
/// \code
LinearSolverUmfpack linsolver;
IncompTpfa psolver(grid, props.permeability(), grav, linsolver, wells);
IncompTpfa psolver(grid, props, linsolver,
grav, wells, src, bcs.c_bcs());
/// \endcode
@@ -312,18 +302,6 @@ int main ()
/// \code
for (int i = 0; i < num_time_steps; ++i) {
/// \endcode
/// \page tutorial4
/// \details Compute the total mobility. It is needed by the pressure solver
/// \code
computeTotalMobility(props, allcells, state.saturation(), totmob);
/// \endcode
/// \endcode
/// \page tutorial4
/// \details In order to use the well controls, we need to generate the WDP for each well.
/// \code
Opm::computeWDP(*wells, grid, state.saturation(), props.density(), gravity, true, wdp);
/// \endcode
/// \page tutorial4
/// \details We're solving the pressure until the well conditions are met
@@ -338,31 +316,29 @@ int main ()
/// \page tutorial4
/// \details Solve the pressure equation
/// \code
psolver.solve(totmob, omega, src, wdp, bcs.c_bcs(),
state.pressure(), state.faceflux(), well_bhp,
well_flux);
psolver.solve(dt, state, well_state);
/// \endcode
/// \page tutorial4
/// \details We compute the new well rates. Notice that we approximate (wrongly) surfflowsrates := resflowsrate
Opm::computeFractionalFlow(props, allcells, state.saturation(), fractional_flows);
Opm::computePhaseFlowRatesPerWell(*wells, well_flux, fractional_flows, well_resflowrates_phase);
Opm::computePhaseFlowRatesPerWell(*wells, well_flux, fractional_flows, well_surflowrates_phase);
Opm::computePhaseFlowRatesPerWell(*wells, well_state.perfRates(), fractional_flows, well_resflowrates_phase);
Opm::computePhaseFlowRatesPerWell(*wells, well_state.perfRates(), fractional_flows, well_surflowrates_phase);
/// \endcode
/// \page tutorial4
/// \details We check if the well conditions are met.
well_conditions_met = well_collection.conditionsMet(well_bhp, well_resflowrates_phase, well_surflowrates_phase);
well_conditions_met = well_collection.conditionsMet(well_state.bhp(), well_resflowrates_phase, well_surflowrates_phase);
++well_iter;
if (!well_conditions_met && well_iter == max_well_iterations) {
THROW("Conditions not met within " << max_well_iterations<< " iterations.");
}
}
/// \endcode
/// \page tutorial4
/// \details Transport solver
/// \TODO We must call computeTransportSource() here, since we have wells.
/// \code
transport_solver.solve(&state.faceflux()[0], &porevol[0], &src[0], dt, state.saturation());
/// \endcode