From 9fd882132b5e90ec53762e2be38bad728d2e6c7f Mon Sep 17 00:00:00 2001 From: Kjetil Olsen Lye Date: Wed, 16 May 2012 16:58:05 +0200 Subject: [PATCH] created first draft of well tutorial --- tutorials/tutorial4.cpp | 315 ++++++++++++++++++++++++++-------------- 1 file changed, 205 insertions(+), 110 deletions(-) diff --git a/tutorials/tutorial4.cpp b/tutorials/tutorial4.cpp index 0e9f4362a..0b09ce8a0 100644 --- a/tutorials/tutorial4.cpp +++ b/tutorials/tutorial4.cpp @@ -49,52 +49,13 @@ #include #include #include +#include -/// \page tutorial3 Multiphase flow -/// The Darcy law gives -/// \f[u_\alpha= -\frac1{\mu_\alpha} K_\alpha\nabla p_\alpha\f] -/// where \f$\mu_\alpha\f$ and \f$K_\alpha\f$ represent the viscosity -/// and the permeability tensor for each phase \f$\alpha\f$. In the two phase -/// case, we have either \f$\alpha=w\f$ or \f$\alpha=o\f$. -/// In this tutorial, we do not take into account capillary pressure so that -/// \f$p=p_w=p_o\f$ and gravity -/// effects. We denote by \f$K\f$ the absolute permeability tensor and each phase -/// permeability is defined through its relative permeability by the expression -/// \f[K_\alpha=k_{r\alpha}K.\f] -/// The phase mobility are defined as -/// \f[\lambda_\alpha=\frac{k_{r\alpha}}{\mu_\alpha}\f] -/// so that the Darcy law may be rewritten as -/// \f[u_\alpha= -\lambda_\alpha K\nabla p.\f] -/// The conservation of mass for each phase writes: -/// \f[\frac{\partial}{\partial t}(\phi\rho_\alpha s_\alpha)+\nabla\cdot(\rho_\alpha u_\alpha)=q_\alpha\f] -/// where \f$s_\alpha\f$ denotes the saturation of the phase \f$\alpha\f$ and \f$q_\alpha\f$ is a source term. Let -/// us consider a two phase flow with oil and water. We assume that the phases are incompressible. Since -/// \f$s_w+s_o=1\f$, we get -/// \f[ \nabla\cdot u=\frac{q_w}{\rho_w}+\frac{q_o}{\rho_o}.\f] -/// where we define -/// \f[u=u_w+u_o.\f] -/// Let the total mobility be equal to -/// \f[\lambda=\lambda_w+\lambda_o\f] -/// Then, we have -/// \f[u=-\lambda K\nabla p.\f] -/// The set of equations -/// \f[\nabla\cdot u=\frac{q_w}{\rho_w}+\frac{q_o}{\rho_o},\quad u=-\lambda K\nabla p.\f] -/// is referred to as the pressure equation. We introduce -/// the fractional flow \f$f_w\f$ -/// as -/// \f[f_w=\frac{\lambda_w}{\lambda_w+\lambda_o}\f] -/// and obtain -/// \f[\phi\frac{\partial s}{\partial t}+\nabla\cdot(f_w u)=\frac{q_w}{\rho_w}\f] -/// which is referred to as the transport equation. The pressure and -/// transport equation are coupled. In this tutorial, we implement a splitting scheme, -/// where, at each time step, we decouple the two equations. We solve first -/// the pressure equation and then update the water saturation by solving -/// the transport equation assuming that \f$u\f$ is constant in time in the time step -/// interval we are considering. - -/// \page tutorial3 -/// \section commentedcode3 Commented code: -/// \page tutorial3 +/// \page tutorial4 Well controls +/// +/// \page tutorial4 +/// \section commentedcode4 Commented code: +/// \page tutorial4 /// \details /// We define a class which computes mobility, capillary pressure and /// the minimum and maximum saturation value for each cell. @@ -103,20 +64,20 @@ class TwophaseFluid { public: /// \endcode - /// \page tutorial3 + /// \page tutorial4 /// \details Constructor operator. Takes in the fluid properties defined /// \c props /// \code TwophaseFluid(const Opm::IncompPropertiesInterface& props); /// \endcode - /// \page tutorial3 + /// \page tutorial4 /// \details Density for each phase. /// \code double density(int phase) const; /// \endcode - /// \page tutorial3 + /// \page tutorial4 /// \details Computes the mobility and its derivative with respect to saturation /// for a given cell \c c and saturation \c sat. The template classes \c Sat, /// \c Mob, \c DMob are typically arrays. By using templates, we do not have to @@ -127,7 +88,7 @@ public: void mobility(int c, const Sat& s, Mob& mob, DMob& dmob) const; /// \endcode - /// \page tutorial3 + /// \page tutorial4 /// \details Computes the capillary pressure and its derivative with respect /// to saturation /// for a given cell \c c and saturation \c sat. @@ -136,19 +97,19 @@ public: void pc(int c, const Sat& s, Pcap& pcap, DPcap& dpcap) const; /// \endcode - /// \page tutorial3 + /// \page tutorial4 /// \details Returns the minimum permitted saturation. /// \code double s_min(int c) const; /// \endcode - /// \page tutorial3 + /// \page tutorial4 /// \details Returns the maximum permitted saturation /// \code double s_max(int c) const; /// \endcode - /// \page tutorial3 + /// \page tutorial4 /// \details Private variables /// \code private: @@ -158,13 +119,13 @@ private: }; /// \endcode -/// \page tutorial3 +/// \page tutorial4 /// \details We set up the transport model. /// \code typedef Opm::SinglePointUpwindTwoPhase TransportModel; /// \endcode -/// \page tutorial3 +/// \page tutorial4 /// \details /// The transport equation is nonlinear. We use an implicit transport solver /// which implements a Newton-Raphson solver. @@ -185,7 +146,7 @@ public: }; /// \endcode -/// \page tutorial3 +/// \page tutorial4 /// \details /// We set up the solver. /// \code @@ -199,14 +160,14 @@ typedef Opm::ImplicitTransportnumber_of_cells; /// \endcode - /// \page tutorial3 + /// \page tutorial4 /// \details /// We define the properties of the fluid.\n /// Number of phases. @@ -231,24 +192,24 @@ int main () using namespace unit; using namespace prefix; /// \endcode - /// \page tutorial3 + /// \page tutorial4 /// \details density vector (one component per phase). /// \code std::vector rho(2, 1000.); /// \endcode - /// \page tutorial3 + /// \page tutorial4 /// \details viscosity vector (one component per phase). /// \code std::vector mu(2, 1.*centi*Poise); /// \endcode - /// \page tutorial3 + /// \page tutorial4 /// \details porosity and permeability of the rock. /// \code double porosity = 0.5; double k = 10*milli*darcy; /// \endcode - /// \page tutorial3 + /// \page tutorial4 /// \details We define the relative permeability function. We use a basic fluid /// description and set this function to be linear. For more realistic fluid, the /// saturation function is given by the data. @@ -257,7 +218,7 @@ int main () rel_perm_func = SaturationPropsBasic::Linear; /// \endcode - /// \page tutorial3 + /// \page tutorial4 /// \details We construct a basic fluid with the properties we have defined above. /// Each property is constant and hold for all cells. /// \code @@ -267,23 +228,17 @@ int main () /// \endcode - /// \page tutorial3 + /// \page tutorial4 /// \details Gravity parameters. Here, we set zero gravity. /// \code const double *grav = 0; std::vector omega; /// \endcode - /// \page tutorial3 - /// \details We set up the pressure solver. - /// \code - LinearSolverUmfpack linsolver; - IncompTpfa psolver(*grid.c_grid(), props.permeability(), grav, linsolver); - /// \endcode - + - /// \page tutorial3 + /// \page tutorial4 /// \details We set up the source term /// \code std::vector src(num_cells, 0.0); @@ -291,15 +246,19 @@ int main () src[num_cells-1] = -1.; /// \endcode - /// \page tutorial3 - /// \details We set up the wells. Here, there are no well and we let them empty. + /// \page tutorial4 + /// \details We set up necessary information for the wells /// \code - std::vector empty_wdp; - std::vector empty_well_bhp; - std::vector empty_well_flux; + std::vector wdp; + std::vector well_bhp; + std::vector well_flux; + std::vector well_resflowrates_phase; + std::vector well_surflowrates_phase; + std::vector fractional_flows; + /// \endcode - /// \page tutorial3 + /// \page tutorial4 /// \details We set up the source term for the transport solver. /// \code TransportSource* tsrc = create_transport_source(2, 2); @@ -315,14 +274,14 @@ int main () } /// \endcode - /// \page tutorial3 + /// \page tutorial4 /// \details We compute the pore volume /// \code std::vector porevol; Opm::computePorevolume(*grid.c_grid(), props.porosity(), porevol); /// \endcode - /// \page tutorial3 + /// \page tutorial4 /// \details We set up the transport solver. /// \code const bool guess_old_solution = true; @@ -330,14 +289,14 @@ int main () TransportSolver tsolver(model); /// \endcode - /// \page tutorial3 + /// \page tutorial4 /// \details Time integration parameters /// \code double dt = 0.1*day; int num_time_steps = 20; /// \endcode - /// \page tutorial3 + /// \page tutorial4 /// \details Control paramaters for the implicit solver. /// \code ImplicitTransportDetails::NRReport rpt; @@ -346,7 +305,7 @@ int main () - /// \page tutorial3 + /// \page tutorial4 /// \details We define a vector which contains all cell indexes. We use this /// vector to set up parameters on the whole domains. std::vector allcells(num_cells); @@ -355,14 +314,14 @@ int main () } - /// \page tutorial3 + /// \page tutorial4 /// \details We set up the boundary conditions. Letting bcs empty is equivalent /// to no flow boundary conditions. /// \code FlowBCManager bcs; /// \endcode - /// \page tutorial3 + /// \page tutorial4 /// \details /// Linear solver init. /// \code @@ -370,7 +329,7 @@ int main () CSRMatrixUmfpackSolver linsolve; /// \endcode - /// \page tutorial3 + /// \page tutorial4 /// \details /// We initialise water saturation to minimum everywhere. /// \code @@ -379,48 +338,184 @@ int main () state.setWaterSat(allcells, props, TwophaseState::MinSat); /// \endcode - /// \page tutorial3 + /// \page tutorial4 /// \details We introduce a vector which contains the total mobility /// on all cells. /// \code std::vector totmob; /// \endcode - /// \page tutorial3 + /// \page tutorial4 /// \details This string will contain the name of a VTK output vector. /// \code std::ostringstream vtkfilename; /// \endcode - /// \page tutorial3 + /// \page tutorial4 + /// To create wells we need an instance of the PhaseUsage-object + PhaseUsage phase_usage; + phase_usage.num_phases = num_phases; + phase_usage.phase_used[BlackoilPhases::Aqua] = 1; + phase_usage.phase_used[BlackoilPhases::Liquid] = 1; + phase_usage.phase_used[BlackoilPhases::Vapour] = 0; + + phase_usage.phase_pos[BlackoilPhases::Aqua] = 0; + phase_usage.phase_pos[BlackoilPhases::Liquid] = 1; + + /// \page tutorial4 + /// \details This will contain our well-specific information + /// \code + WellCollection well_collection; + /// \endcode + + /// \page tutorial4 + /// \details Create the production specification for our top well group. + /// We set a target limit for total reservoir rate, and set the controlling + /// mode of the group to be controlled by the reservoir rate. + /// \code + ProductionSpecification well_group_prod_spec; + well_group_prod_spec.reservoir_flow_max_rate_ = 0.1; + well_group_prod_spec.control_mode_ = ProductionSpecification::RESV; + /// \endcode + + /// \page tutorial4 + /// \details Create our well group. We hand it an empty injection specification, + /// as we don't want to control it's injection target. We use the shared_ptr-type because that's + /// what the interface expects. The first argument is the (unique) name of the group. + /// \code + std::tr1::shared_ptr well_group(new WellsGroup("group", well_group_prod_spec, InjectionSpecification(), + phase_usage)); + /// \endcode + + /// \page tutorial4 + /// \details We add our well_group to the well_collection + /// \code + well_collection.addChild(well_group); + /// \endcode + + + /// \page tutorial4 + /// \details Create the production specification and Well objects (total 4 wells). We set all our wells to be group controlled. + /// We pass in the string argument \C "group" to set the parent group. + /// \code + const int num_wells = 4; + for (int i = 0; i < num_wells; ++i) { + std::stringstream well_name; + well_name << "well" << i; + ProductionSpecification production_specification; + production_specification.control_mode_ = ProductionSpecification::GRUP; + std::tr1::shared_ptr well_leaf_node(new WellNode(well_name.str(), production_specification, InjectionSpecification(), + phase_usage)); + well_collection.addChild(well_leaf_node, "group"); + + } + /// \endcode + + /// \page tutorial4 + /// \details Now we create the C struct to hold our wells (this is to interface with the solver code). For now we + /// \code + Wells* wells = create_wells(num_phases, num_wells, num_wells /*number of perforations. We'll only have one perforation per well*/); + /// \endcode + + /// + + /// \page tutorial4 + /// \details We need to add each well to the C API. + /// To do this we need to specify the relevant cells the well will be located in (\C well_cells). + /// \code + for (int i = 0; i < num_wells; ++i) { + const int well_cells = i*nx; + const double well_index = 1; + add_well(PRODUCER, 0, 1, NULL, &well_cells, &well_index, wells); + } + /// \endcode + + /// \page tutorial4 + /// \details We need to make the well collection aware of our wells object + /// \code + well_collection.setWellsPointer(wells); + /// \endcode + /// \page tutorial4 + /// We're not using well controls, just group controls, so we need to apply them. + /// \code + well_collection.applyGroupControls(); + ///\endcode + + /// \page tutorial4 + /// \details We set up the pressure solver. We need to pass the wells pointer as the + /// last argument. + /// \code + LinearSolverUmfpack linsolver; + IncompTpfa psolver(*grid.c_grid(), props.permeability(), grav, linsolver, wells); + /// \endcode + + + /// \page tutorial4 /// \details Loop over the time steps. /// \code for (int i = 0; i < num_time_steps; ++i) { /// \endcode - /// \page tutorial3 + /// \page tutorial4 /// \details Compute the total mobility. It is needed by the pressure solver /// \code computeTotalMobility(props, allcells, state.saturation(), totmob); /// \endcode - /// \page tutorial3 - /// \details Solve the pressure equation + + /// \endcode + /// \page tutorial4 + /// \details In order to use the well controls, we need to generate the WDP for each well. /// \code - psolver.solve(totmob, omega, src, empty_wdp, bcs.c_bcs(), - state.pressure(), state.faceflux(), empty_well_bhp, - empty_well_flux); + Opm::computeWDP(*wells, *grid.c_grid(), state.saturation(), props.density(), gravity, true, wdp); /// \endcode - /// \page tutorial3 + + /// \page tutorial4 + /// \details We're solving the pressure until the well conditions are met + /// or until we reach the maximum number of iterations. + /// \code + const int max_well_iterations = 10; + int well_iter = 0; + bool well_conditions_met = false; + while (!well_conditions_met) { + /// \endcode + + /// \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); + + /// \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); + /// \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_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 /// \code tsolver.solve(*grid.c_grid(), tsrc, dt, ctrl, state, linsolve, rpt); /// \endcode - /// \page tutorial3 + /// \page tutorial4 /// \details Write the output to file. /// \code vtkfilename.str(""); - vtkfilename << "tutorial3-" << std::setw(3) << std::setfill('0') << i << ".vtu"; + vtkfilename << "tutorial4-" << std::setw(3) << std::setfill('0') << i << ".vtu"; std::ofstream vtkfile(vtkfilename.str().c_str()); Opm::DataMap dm; dm["saturation"] = &state.saturation(); @@ -430,7 +525,7 @@ int main () } /// \endcode -/// \page tutorial3 +/// \page tutorial4 /// \details Implementation of the TwophaseFluid class. /// \code TwophaseFluid::TwophaseFluid(const Opm::IncompPropertiesInterface& props) @@ -461,7 +556,7 @@ void TwophaseFluid::mobility(int c, const Sat& s, Mob& mob, DMob& dmob) const mob[0] /= mu[0]; mob[1] /= mu[1]; /// \endcode - /// \page tutorial3 + /// \page tutorial4 /// \details We /// recall that we use Fortran ordering for kr derivatives, /// therefore dmob[i*2 + j] is row j and column i of the @@ -496,26 +591,26 @@ double TwophaseFluid::s_max(int c) const /// \endcode -/// \page tutorial3 -/// \section results3 Results. +/// \page tutorial4 +/// \section results4 Results. /// /// -/// -/// +/// +/// /// /// -/// -/// +/// +/// /// /// -/// +/// /// /// ///
\image html tutorial3-000.png \image html tutorial3-005.png \image html tutorial4-000.png \image html tutorial4-005.png
\image html tutorial3-010.png \image html tutorial3-015.png \image html tutorial4-010.png \image html tutorial4-015.png
\image html tutorial3-019.png \image html tutorial4-019.png
-/// \page tutorial3 +/// \page tutorial4 /// \details -/// \section completecode3 Complete source code: -/// \include tutorial3.cpp +/// \section completecode4 Complete source code: +/// \include tutorial4.cpp /// \include generate_doc_figures.py