OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

code is now included using \snippet. Apparently this looks better with the new Doxygen version. The HTML_EXTRA_STYLESHEET is now used rather then the HTML_STYLESHEET in order to include used-defined styles for the same reason

Browse Source
This commit is contained in:
Tor Harald Sandve
2013-03-06 10:17:27 +01:00
parent 7285f2bfd1
commit 20e7c63d69
5 changed files with 355 additions and 177 deletions
Download Patch File Download Diff File Expand all files Collapse all files

9
tutorials/generate_doc_figures.py
View File

@@ -38,7 +38,7 @@ collected_garbage_file = []
if not isdir(figure_path):
mkdir(figure_path)
## [tutorial1]
# tutorial 1
data_file_name = join(tutorial_data_path, "tutorial1.vtu")
# grid = servermanager.sources.XMLUnstructuredGridReader(FileName = data_file_name)
@@ -60,7 +60,9 @@ camera.SetFocalPoint(1.5, 1.5, 1)
Render()
WriteImage(join(figure_path, "tutorial1.png"))
Hide(grid)
## [tutorial1]
## [tutorial2]
# tutorial 2
data_file_name = join(tutorial_data_path, "tutorial2.vtu")
grid = XMLUnstructuredGridReader(FileName = data_file_name)
@@ -83,7 +85,9 @@ camera.SetFocalPoint(20, 20, 0.5)
Render()
WriteImage(join(figure_path, "tutorial2.png"))
Hide(grid)
## [tutorial2]
## [tutorial3]
# tutorial 3
for case in range(0,20):
data_file_name = join(tutorial_data_path, "tutorial3-"+"%(case)03d"%{"case": case}+".vtu")
@@ -110,7 +114,9 @@ for case in cases:
Render()
WriteImage(join(figure_path, "tutorial3-"+case+".png"))
Hide(grid)
## [tutorial3]
## [tutorial4]
# tutorial 4
for case in range(0,20):
data_file_name = join(tutorial_data_path, "tutorial4-"+"%(case)03d"%{"case": case}+".vtu")
@@ -137,6 +143,7 @@ for case in cases:
Render()
WriteImage(join(figure_path, "tutorial4-"+case+".png"))
Hide(grid)
## [tutorial4]
# remove temporary files
for f in collected_garbage_file:

62
tutorials/tutorial1.cpp
View File

@@ -23,31 +23,25 @@
#endif // HAVE_CONFIG_H
/// \page tutorial1 A simple cartesian grid
/// This tutorial explains how to construct a simple cartesian grid,
/// This tutorial explains how to construct a simple Cartesian grid,
/// and we will take a look at some output facilities.
/// \page tutorial1
/// \section commentedsource1 Program walkthrough.
/// \section commentedsource1 Program walk-through.
/// All headers from opm-core are found in the opm/core/ directory.
/// Some important headers are at the root, other headers are found
/// in subdirectories.
#include <opm/core/grid.h>
#include <opm/core/GridManager.hpp>
#include <opm/core/utility/writeVtkData.hpp>
#include <iostream>
#include <fstream>
#include <vector>
/// \snippet tutorial1.cpp including headers
/**
\code
/// \internal [including headers]
#include <opm/core/grid.h>
#include <opm/core/GridManager.hpp>
#include <opm/core/utility/writeVtkData.hpp>
#include <iostream>
#include <fstream>
#include <vector>
\endcode
*/
/// \internal [including headers]
/// \endinternal
// ----------------- Main program -----------------
@@ -55,18 +49,22 @@ int main()
{
/// \page tutorial1
/// We set the number of blocks in each direction.
/// \code
/// \snippet tutorial1.cpp num blocks
/// \internal [num blocks]
int nx = 4;
int ny = 3;
int nz = 2;
/// \endcode
/// \internal [num blocks]
/// \endinternal
/// The size of each block is 1m x 1m x 1m. The default units are always the
/// standard units (SI). But other units can easily be dealt with, see Opm::unit.
/// \code
/// \snippet tutorial1.cpp dim
/// \internal [dim]
double dx = 1.0;
double dy = 1.0;
double dz = 1.0;
/// \endcode
/// \internal [dim]
/// \endinternal
/// \page tutorial1
/// In opm-core, grid information is accessed via the UnstructuredGrid data structure.
/// This data structure has a pure C API, including helper functions to construct and
@@ -74,29 +72,38 @@ int main()
/// which is a C++ class that wraps the UnstructuredGrid and takes care of
/// object lifetime issues.
/// One of the constructors of the class Opm::GridManager takes <code>nx, ny, nz, dx, dy, dz</code>
/// and construct the corresponding cartesian grid.
/// \code
/// and construct the corresponding Cartesian grid.
/// \snippet tutorial1.cpp grid manager
/// \internal [grid manager]
Opm::GridManager grid(nx, ny, nz, dx, dy, dz);
/// \endcode
/// \internal [grid manager]
/// \endinternal
/// \page tutorial1
/// We open an output file stream for the output
/// \code
/// \snippet tutorial1.cpp output stream
/// \internal [output stream]
std::ofstream vtkfile("tutorial1.vtu");
/// \endcode
/// \internal [output stream]
/// \endinternal
/// \page tutorial1
/// The Opm::writeVtkData() function writes a grid together with
/// data to a stream. Here, we just want to visualize the grid. We
/// construct an empty Opm::DataMap object, which we send to
/// Opm::writeVtkData() together with the grid
/// \code
/// \snippet tutorial1.cpp data map
/// \internal [data map]
Opm::DataMap dm;
/// \endcode
/// \internal [data map]
/// \endinternal
/// \page tutorial1
/// Call Opm::writeVtkData() to write the output file.
/// \code
/// \snippet tutorial1.cpp write vtk
/// \internal [write vtk]
Opm::writeVtkData(*grid.c_grid(), dm, vtkfile);
/// \internal [write vtk]
/// \endinternal
}
/// \endcode
/// \page tutorial1
/// We read the vtu output file in \a Paraview and obtain the following grid.
/// \image html tutorial1.png
@@ -105,3 +112,8 @@ int main()
/// \section completecode1 Complete source code:
/// \include tutorial1.cpp
/// \page tutorial1
/// \details
/// \section pythonscript1 Python script to generate figures:
/// \snippet generate_doc_figures.py tutorial1

95
tutorials/tutorial2.cpp
View File

@@ -24,7 +24,7 @@
/// \f${\bf u}\f$ denotes the velocity and \f$p\f$ the pressure. The permeability tensor is
/// given by \f$K\f$ and \f$\mu\f$ denotes the viscosity.
///
/// We solve the flow equations for a cartesian grid and we set the source term
/// We solve the flow equations for a Cartesian grid and we set the source term
/// \f$q\f$ be zero except at the left-lower and right-upper corner, where it is equal
/// with opposite sign (inflow equal to outflow).
@@ -49,28 +49,33 @@
#include <opm/core/simulator/WellState.hpp>
/// \page tutorial2
/// \section commentedcode2 Program walkthrough.
/// \code
/// \section commentedcode2 Program walk-through.
///
int main()
{
/// \endcode
/// \page tutorial2
/// We construct a cartesian grid
/// \code
/// We construct a Cartesian grid
/// \snippet tutorial2.cpp cartesian grid
/// \internal [cartesian grid]
int dim = 3;
int nx = 40;
int ny = 40;
int nz = 1;
Opm::GridManager grid(nx, ny, nz);
/// \endcode
/// \internal [cartesian grid]
/// \endinternal
/// \page tutorial2
/// \details We access the unstructured grid through
/// the pointer given by \c grid.c_grid(). For more details on the
/// UnstructuredGrid data structure, see grid.h.
/// \code
/// \snippet tutorial2.cpp access grid
/// \internal [access grid]
int num_cells = grid.c_grid()->number_of_cells;
int num_faces = grid.c_grid()->number_of_faces;
/// \endcode
/// \internal [access grid]
/// endinternal
/// \page tutorial2
@@ -80,90 +85,107 @@ int main()
/// The <opm/core/utility/Units.hpp> header contains support
/// for common units and prefixes, in the namespaces Opm::unit
/// and Opm::prefix.
/// \code
/// \snippet tutorial2.cpp fluid
/// \internal [fluid]
using namespace Opm::unit;
using namespace Opm::prefix;
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
/// \internal [fluid]
/// \endinternal
/// \page tutorial2
/// \details
/// We define a permeability equal to 100 mD.
/// \code
/// \snippet tutorial2.cpp perm
/// \internal [perm]
double k = 100.0*milli*darcy;
/// \endcode
/// \page tutorial2
/// \details
/// \internal [perm]
/// \endinternal
/// \page tutorial2
/// \details
/// We set up a simple property object for a single-phase situation.
/// \code
/// \snippet tutorial2.cpp single-phase property
/// \internal [single-phase property]
Opm::IncompPropertiesBasic props(1, Opm::SaturationPropsBasic::Constant, rho,
mu, 1.0, k, dim, num_cells);
/// \endcode
/// \internal [single-phase property]
/// /endinternal
/// \page tutorial2
/// \details
/// We take UMFPACK as the linear solver for the pressure solver
/// (this library has therefore to be installed).
/// \code
/// \snippet tutorial2.cpp linsolver
/// \internal [linsolver]
Opm::LinearSolverUmfpack linsolver;
/// \endcode
/// \page tutorial2
/// \internal [linsolver]
/// \endinternal
/// \endcode
/// \page tutorial2
/// We define the source term.
/// \code
/// \snippet tutorial2.cpp source
/// \internal [source]
std::vector<double> src(num_cells, 0.0);
src[0] = 100.;
src[num_cells-1] = -100.;
/// \endcode
/// \internal [source]
/// \endinternal
/// \page tutorial2
/// \details We set up the boundary conditions.
/// By default, we obtain no-flow boundary conditions.
/// \code
/// \snippet tutorial2.cpp boundary
/// \internal [boundary]
Opm::FlowBCManager bcs;
/// \endcode
/// \internal [boundary]
/// \endinternal
/// 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
/// \snippet tutorial2.cpp tpfa
/// \internal [tpfa]
Opm::IncompTpfa psolver(*grid.c_grid(), props, linsolver, NULL, NULL, src, NULL);
/// \internal [tpfa]
/// \endinternal
/// \page tutorial2
/// 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
/// \snippet tutorial2.cpp state
/// \internal [state]
Opm::TwophaseState state;
state.pressure().resize(num_cells, 0.0);
state.faceflux().resize(num_faces, 0.0);
state.saturation().resize(num_cells, 1.0);
Opm::WellState well_state;
/// \endcode
/// \internal [state]
/// \endinternal
/// \page tutorial2
/// We call the pressure solver.
/// The first (timestep) argument does not matter for this
/// incompressible case.
/// \code
/// \snippet tutorial2.cpp pressure solver
/// \internal [pressure solver]
psolver.solve(1.0*day, state, well_state);
/// \endcode
/// \internal [pressure solver]
/// \endinternal
/// \page tutorial2
/// We write the results to a file in VTK format.
/// The data vectors added to the Opm::DataMap must
/// contain cell data. They may be a scalar per cell
/// (pressure) or a vector per cell (cell_velocity).
/// \code
/// \snippet tutorial2.cpp write output
/// \internal [write output]
std::ofstream vtkfile("tutorial2.vtu");
Opm::DataMap dm;
dm["pressure"] = &state.pressure();
@@ -171,8 +193,10 @@ int main()
Opm::estimateCellVelocity(*grid.c_grid(), state.faceflux(), cell_velocity);
dm["velocity"] = &cell_velocity;
Opm::writeVtkData(*grid.c_grid(), dm, vtkfile);
/// \internal [write output]
/// \endinternal
}
/// \endcode
/// \page tutorial2
/// We read the vtu output file in \a Paraview and obtain the following pressure
/// distribution. \image html tutorial2.png
@@ -181,3 +205,8 @@ int main()
/// \page tutorial2
/// \section completecode2 Complete source code:
/// \include tutorial2.cpp
/// \page tutorial2
/// \details
/// \section pythonscript2 python script to generate figures:
/// \snippet generate_doc_figures.py tutorial2

136
tutorials/tutorial3.cpp
View File

@@ -89,15 +89,19 @@
/// \page tutorial3
/// \section commentedsource1 Program walk-through.
/// \details
/// Main function
/// \code
/// \snippet tutorial3.cpp main
/// \internal [main]
int main ()
{
/// \endcode
/// \internal [main]
/// \endinternal
/// \page tutorial3
/// \details
/// We define the grid. A cartesian grid with 400 cells,
/// We define the grid. A Cartesian grid with 400 cells,
/// each being 10m along each side. Note that we treat the
/// grid as 3-dimensional, but have a thickness of only one
/// layer in the Z direction.
@@ -105,7 +109,8 @@ int main ()
/// The Opm::GridManager is responsible for creating and destroying the grid,
/// the UnstructuredGrid data structure contains the actual grid topology
/// and geometry.
/// \code
/// \snippet tutorial3.cpp grid
/// \internal [grid]
int nx = 20;
int ny = 20;
int nz = 1;
@@ -116,7 +121,8 @@ int main ()
GridManager grid_manager(nx, ny, nz, dx, dy, dz);
const UnstructuredGrid& grid = *grid_manager.c_grid();
int num_cells = grid.number_of_cells;
/// \endcode
/// \internal [grid]
/// \endinternal
/// \page tutorial3
/// \details
@@ -128,7 +134,8 @@ int main ()
/// available for use, however. They are stored as constants in
/// the Opm::unit namespace, while prefixes are in the Opm::prefix
/// namespace. See Units.hpp for more.
/// \code
/// \snippet tutorial3.cpp set properties
/// \internal [set properties]
int num_phases = 2;
using namespace Opm::unit;
using namespace Opm::prefix;
@@ -136,139 +143,173 @@ int main ()
std::vector<double> viscosity(num_phases, 1.0*centi*Poise);
double porosity = 0.5;
double permeability = 10.0*milli*darcy;
/// \endcode
/// \internal [set properties]
/// \endinternal
/// \page tutorial3
/// \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 may be interpolated from experimental data.
/// \code
/// \snippet tutorial3.cpp relperm
/// \internal [relperm]
SaturationPropsBasic::RelPermFunc rel_perm_func = SaturationPropsBasic::Linear;
/// \endcode
/// \internal [relperm]
/// \endinternal
/// \page tutorial3
/// \details We construct a basic fluid and rock property object
/// with the properties we have defined above. Each property is
/// constant and hold for all cells.
/// \code
/// \snippet tutorial3.cpp properties
/// \internal [properties]
IncompPropertiesBasic props(num_phases, rel_perm_func, density, viscosity,
porosity, permeability, grid.dimensions, num_cells);
/// \endcode
/// \internal [properties]
/// \endinternal
/// \page tutorial3
/// \details Gravity parameters. Here, we set zero gravity.
/// \code
/// \snippet tutorial3.cpp gravity
/// \internal [gravity]
const double *grav = 0;
std::vector<double> omega;
/// \endcode
/// \internal [gravity]
/// \endinternal
/// \page tutorial3
/// \details We set up the source term. Positive numbers indicate that the cell is a source,
/// while negative numbers indicate a sink.
/// \code
/// \snippet tutorial3.cpp source
/// \internal [source]
std::vector<double> src(num_cells, 0.0);
src[0] = 1.;
src[num_cells-1] = -1.;
/// \endcode
/// \internal [source]
/// \endinternal
/// \page tutorial3
/// \details We set up the boundary conditions. Letting bcs be empty is equivalent
/// to no-flow boundary conditions.
/// \code
/// \snippet tutorial3.cpp boundary
/// \internal [boundary]
FlowBCManager bcs;
/// \endcode
/// \internal [boundary]
/// \endinternal
/// \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
/// \snippet tutorial3.cpp pressure solver
/// \internal [pressure solver]
LinearSolverUmfpack linsolver;
IncompTpfa psolver(grid, props, linsolver, grav, NULL, src, bcs.c_bcs());
/// \endcode
/// \internal [pressure solver]
/// \endinternal
/// \page tutorial3
/// \details We set up a state object for the wells. Here, there are
/// no wells and we let it remain empty.
/// \code
/// \snippet tutorial3.cpp well
/// \internal [well]
WellState well_state;
/// \endcode
/// \internal [well]
/// \endinternal
/// \page tutorial3
/// \details We compute the pore volume
/// \code
/// \snippet tutorial3.cpp pore volume
/// \internal [pore volume]
std::vector<double> porevol;
Opm::computePorevolume(grid, props.porosity(), porevol);
/// \endcode
/// \internal [pore volume]
/// \endinternal
/// \page tutorial3
/// \details Set up the transport solver. This is a reordering implicit Euler transport solver.
/// \code
/// \snippet tutorial3.cpp transport solver
/// \internal [transport solver]
const double tolerance = 1e-9;
const int max_iterations = 30;
Opm::TransportModelTwophase transport_solver(grid, props, tolerance, max_iterations);
/// \endcode
/// \internal [transport solver]
/// \endinternal
/// \page tutorial3
/// \details Time integration parameters
/// \code
/// \snippet tutorial3.cpp time parameters
/// \internal [time parameters]
const double dt = 0.1*day;
const int num_time_steps = 20;
/// \endcode
/// \internal [time parameters]
/// \endinternal
/// \page tutorial3
/// \details We define a vector which contains all cell indexes. We use this
/// vector to set up parameters on the whole domain.
/// \code
/// \snippet tutorial3.cpp cell indexes
/// \internal [cell indexes]
std::vector<int> allcells(num_cells);
for (int cell = 0; cell < num_cells; ++cell) {
allcells[cell] = cell;
}
/// \endcode
/// \internal [cell indexes]
/// \endinternal
/// \page tutorial3
/// \details
/// We set up a two-phase state object, and
/// initialise water saturation to minimum everywhere.
/// \code
/// initialize water saturation to minimum everywhere.
/// \snippet tutorial3.cpp two-phase state
/// \internal [two-phase state]
TwophaseState state;
state.init(grid, 2);
state.setFirstSat(allcells, props, TwophaseState::MinSat);
/// \endcode
/// \internal [two-phase state]
/// \endinternal
/// \page tutorial3
/// \details This string stream will be used to construct a new
/// output filename at each timestep.
/// \code
/// \snippet tutorial3.cpp output stream
/// \internal [output stream]
std::ostringstream vtkfilename;
/// \endcode
/// \internal [output stream]
/// \endinternal
/// \page tutorial3
/// \details Loop over the time steps.
/// \code
/// \snippet tutorial3.cpp time loop
/// \internal [time loop]
for (int i = 0; i < num_time_steps; ++i) {
/// \endcode
/// \page tutorial3
/// \internal [time loop]
/// \endinternal
/// \endcode
/// \page tutorial3
/// \details Solve the pressure equation
/// \code
/// \snippet tutorial3.cpp solve pressure
/// \internal [solve pressure]
psolver.solve(dt, state, well_state);
/// \endcode
/// \internal [solve pressure]
/// \endinternal
/// \page tutorial3
/// \details Solve the transport equation.
/// \code
/// \snippet tutorial3.cpp transport solve
/// \internal [transport solve]
transport_solver.solve(&state.faceflux()[0], &porevol[0], &src[0],
dt, state.saturation());
/// \endcode
/// \internal [transport solve]
/// \endinternal
/// \page tutorial3
/// \details Write the output to file.
/// \code
/// \snippet tutorial3.cpp write output
/// \internal [write output]
vtkfilename.str("");
vtkfilename << "tutorial3-" << std::setw(3) << std::setfill('0') << i << ".vtu";
std::ofstream vtkfile(vtkfilename.str().c_str());
@@ -278,7 +319,8 @@ int main ()
Opm::writeVtkData(grid, dm, vtkfile);
}
}
/// \endcode
/// \internal [write output]
/// \endinternal
@@ -304,4 +346,8 @@ int main ()
/// \details
/// \section completecode3 Complete source code:
/// \include tutorial3.cpp
/// \include generate_doc_figures.py
/// \page tutorial3
/// \details
/// \section pythonscript3 python script to generate figures:
/// \snippet generate_doc_figures.py tutorial3

230
tutorials/tutorial4.cpp
View File

@@ -46,18 +46,22 @@
#include <opm/core/wells/WellCollection.hpp>
/// \page tutorial4 Well controls
/// This tutorial explains how to construct an example with wells
/// \page tutorial4
/// \section commentedsource1 Program walk-through.
/// \details
/// Main function
/// \code
/// \snippet tutorial4.cpp main
/// \internal[main]
int main ()
{
/// \endcode
/// \internal[main]
/// \endinternal
/// \page tutorial4
/// \details
/// We define the grid. A cartesian grid with 1200 cells.
/// \code
/// We define the grid. A Cartesian grid with 1200 cells.
/// \snippet tutorial4.cpp cartesian grid
/// \internal[cartesian grid]
int dim = 3;
int nx = 20;
int ny = 20;
@@ -69,126 +73,162 @@ int main ()
GridManager grid_manager(nx, ny, nz, dx, dy, dz);
const UnstructuredGrid& grid = *grid_manager.c_grid();
int num_cells = grid.number_of_cells;
/// \endcode
/// \internal[cartesian grid]
/// \endinternal
/// \page tutorial4
/// \details
/// We define the properties of the fluid.\n
/// Number of phases.
/// \code
/// \snippet tutorial4.cpp Number of phases
/// \internal[Number of phases]
int num_phases = 2;
using namespace unit;
using namespace prefix;
/// \endcode
/// \internal[Number of phases]
/// \endinternal
/// \page tutorial4
/// \details density vector (one component per phase).
/// \code
/// \snippet tutorial4.cpp density
/// \internal[density]
std::vector<double> rho(2, 1000.);
/// \endcode
/// \internal[density]
/// \endinternal
/// \page tutorial4
/// \details viscosity vector (one component per phase).
/// \code
/// \snippet tutorial4.cpp viscosity
/// \internal[viscosity]
std::vector<double> mu(2, 1.*centi*Poise);
/// \endcode
/// \internal[viscosity]
/// \endinternal
/// \page tutorial4
/// \details porosity and permeability of the rock.
/// \code
/// \snippet tutorial4.cpp rock
/// \internal[rock]
double porosity = 0.5;
double k = 10*milli*darcy;
/// \endcode
/// \internal[rock]
/// \endinternal
/// \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.
/// \code
/// \snippet tutorial4.cpp relative permeability
/// \internal[relative permeability]
SaturationPropsBasic::RelPermFunc rel_perm_func = SaturationPropsBasic::Linear;
/// \endcode
/// \internal[relative permeability]
/// \endinternal
/// \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
/// \snippet tutorial4.cpp fluid properties
/// \internal[fluid properties]
IncompPropertiesBasic props(num_phases, rel_perm_func, rho, mu,
porosity, k, dim, num_cells);
/// \endcode
/// \internal[fluid properties]
/// \endinternal
/// \page tutorial4
/// \details Gravity parameters. Here, we set zero gravity.
/// \code
/// \snippet tutorial4.cpp Gravity
/// \internal[Gravity]
const double *grav = 0;
std::vector<double> omega;
/// \endcode
/// \internal[Gravity]
/// \endinternal
/// \page tutorial4
/// \details We set up the source term. Positive numbers indicate that the cell is a source,
/// while negative numbers indicate a sink.
/// \code
/// \snippet tutorial4.cpp source
/// \internal[source]
std::vector<double> src(num_cells, 0.0);
src[0] = 1.;
src[num_cells-1] = -1.;
/// \endcode
/// \internal[source]
/// \endinternal
/// \page tutorial4
/// \details We compute the pore volume
/// \code
/// \snippet tutorial4.cpp pore volume
/// \internal[pore volume]
std::vector<double> porevol;
Opm::computePorevolume(grid, props.porosity(), porevol);
/// \endcode
/// \internal[pore volume]
/// \endinternal
/// \page tutorial4
/// \details Set up the transport solver. This is a reordering implicit Euler transport solver.
/// \code
/// \snippet tutorial4.cpp transport solver
/// \internal[transport solver]
const double tolerance = 1e-9;
const int max_iterations = 30;
Opm::TransportModelTwophase transport_solver(grid, props, tolerance, max_iterations);
/// \endcode
/// \internal[transport solver]
/// \endinternal
/// \page tutorial4
/// \details Time integration parameters
/// \code
/// \snippet tutorial4.cpp Time integration
/// \internal[Time integration]
double dt = 0.1*day;
int num_time_steps = 20;
/// \endcode
/// \internal[Time integration]
/// \endinternal
/// \page tutorial4
/// \details We define a vector which contains all cell indexes. We use this
/// vector to set up parameters on the whole domains.
/// \snippet tutorial4.cpp cell indexes
/// \internal[cell indexes]
std::vector<int> allcells(num_cells);
for (int cell = 0; cell < num_cells; ++cell) {
allcells[cell] = cell;
}
/// \internal[cell indexes]
/// \endinternal
/// \page tutorial4
/// \details We set up the boundary conditions. Letting bcs empty is equivalent
/// to no flow boundary conditions.
/// \code
/// \snippet tutorial4.cpp boundary
/// \internal[boundary]
FlowBCManager bcs;
/// \endcode
/// \internal[boundary]
/// \endinternal
/// \page tutorial4
/// \details
/// We set up a two-phase state object, and
/// initialise water saturation to minimum everywhere.
/// \code
/// \snippet tutorial4.cpp two-phase state
/// \internal[two-phase state]
TwophaseState state;
state.init(grid, 2);
state.setFirstSat(allcells, props, TwophaseState::MinSat);
/// \endcode
/// \internal[two-phase state]
/// \endinternal
/// \page tutorial4
/// \details This string will contain the name of a VTK output vector.
/// \code
/// \snippet tutorial4.cpp VTK output
/// \internal[VTK output]
std::ostringstream vtkfilename;
/// \endcode
/// \internal[VTK output]
/// \endinternal
/// \page tutorial4
/// To create wells we need an instance of the PhaseUsage-object
/// \code
/// \snippet tutorial4.cpp PhaseUsage-object
/// \internal[PhaseUsage-object]
PhaseUsage phase_usage;
phase_usage.num_phases = num_phases;
phase_usage.phase_used[BlackoilPhases::Aqua] = 1;
@@ -197,44 +237,54 @@ int main ()
phase_usage.phase_pos[BlackoilPhases::Aqua] = 0;
phase_usage.phase_pos[BlackoilPhases::Liquid] = 1;
/// \endcode
/// \internal[PhaseUsage-object]
/// \endinternal
/// \page tutorial4
/// \details This will contain our well-specific information
/// \code
/// \snippet tutorial4.cpp well_collection
/// \internal[well_collection]
WellCollection well_collection;
/// \endcode
/// \internal[well_collection]
/// \endinternal
/// \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
/// \snippet tutorial4.cpp production specification
/// \internal[production specification]
ProductionSpecification well_group_prod_spec;
well_group_prod_spec.reservoir_flow_max_rate_ = 0.1;
well_group_prod_spec.control_mode_ = ProductionSpecification::RESV;
/// \endcode
/// \internal[production specification]
/// \endinternal
/// \page tutorial4
/// \details Create our well group. We hand it an empty injection specification,
/// as we don't want to control its 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
/// \snippet tutorial4.cpp injection specification
/// \internal[injection specification]
std::tr1::shared_ptr<WellsGroupInterface> well_group(new WellsGroup("group", well_group_prod_spec, InjectionSpecification(),
phase_usage));
/// \endcode
/// \internal[injection specification]
/// \endinternal
/// \page tutorial4
/// \details We add our well_group to the well_collection
/// \code
/// \snippet tutorial4.cpp well_collection
/// \internal[well_collection]
well_collection.addChild(well_group);
/// \endcode
/// \internal[well_collection]
/// \endinternal
/// \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
/// \snippet tutorial4.cpp create well objects
/// \internal[create well objects]
const int num_wells = 4;
for (int i = 0; i < num_wells; ++i) {
std::stringstream well_name;
@@ -246,20 +296,24 @@ int main ()
well_collection.addChild(well_leaf_node, "group");
}
/// \endcode
/// \internal[create well objects]
/// \endinternal
/// \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
/// \snippet tutorial4.cpp well struct
/// \internal[well struct]
Wells* wells = create_wells(num_phases, num_wells, num_wells /*number of perforations. We'll only have one perforation per well*/);
/// \endcode
/// \internal[well struct]
/// \endinternal
///
/// \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
/// \snippet tutorial4.cpp well cells
/// \internal[well cells]
for (int i = 0; i < num_wells; ++i) {
const int well_cells = i*nx;
const double well_index = 1;
@@ -268,87 +322,112 @@ int main ()
add_well(PRODUCER, 0, 1, NULL, &well_cells, &well_index,
well_name.str().c_str(), wells);
}
/// \endcode
/// \internal[well cells]
/// \endinternal
/// \page tutorial4
/// \details We need to make the well collection aware of our wells object
/// \code
/// \snippet tutorial4.cpp set well pointer
/// \internal[set well pointer]
well_collection.setWellsPointer(wells);
/// \endcode
/// \internal[set well pointer]
/// \endinternal
/// \page tutorial4
/// We're not using well controls, just group controls, so we need to apply them.
/// \code
/// \snippet tutorial4.cpp apply group controls
/// \internal[apply group controls]
well_collection.applyGroupControls();
///\endcode
/// \internal[apply group controls]
/// \endinternal
/// \page tutorial4
/// \details We set up necessary information for the wells
/// \code
/// \snippet tutorial4.cpp init wells
/// \internal[init wells]
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
/// \internal[init wells]
/// \endinternal
/// \page tutorial4
/// \details We set up the pressure solver.
/// \code
/// \snippet tutorial4.cpp pressure solver
/// \internal[pressure solver]
LinearSolverUmfpack linsolver;
IncompTpfa psolver(grid, props, linsolver,
grav, wells, src, bcs.c_bcs());
/// \endcode
/// \internal[pressure solver]
/// \endinternal
/// \page tutorial4
/// \details Loop over the time steps.
/// \code
/// \snippet tutorial4.cpp time loop
/// \internal[time loop]
for (int i = 0; i < num_time_steps; ++i) {
/// \endcode
/// \internal[time loop]
/// \endinternal
/// \page tutorial4
/// \details We're solving the pressure until the well conditions are met
/// or until we reach the maximum number of iterations.
/// \code
/// \snippet tutorial4.cpp well iterations
/// \internal[well iterations]
const int max_well_iterations = 10;
int well_iter = 0;
bool well_conditions_met = false;
while (!well_conditions_met) {
/// \endcode
/// \internal[well iterations]
/// \endinternal
/// \page tutorial4
/// \details Solve the pressure equation
/// \code
/// \snippet tutorial4.cpp pressure solve
/// \internal[pressure solve]
psolver.solve(dt, state, well_state);
/// \endcode
/// \internal[pressure solve]
/// \endinternal
/// \page tutorial4
/// \details We compute the new well rates. Notice that we approximate (wrongly) surfflowsrates := resflowsrate
/// \snippet tutorial4.cpp compute well rates
/// \internal[compute well rates]
Opm::computeFractionalFlow(props, allcells, state.saturation(), fractional_flows);
Opm::computePhaseFlowRatesPerWell(*wells, well_state.perfRates(), fractional_flows, well_resflowrates_phase);
Opm::computePhaseFlowRatesPerWell(*wells, well_state.perfRates(), fractional_flows, well_surflowrates_phase);
/// \endcode
/// \internal[compute well rates]
/// \endinternal
/// \page tutorial4
/// \details We check if the well conditions are met.
/// \snippet tutorial4.cpp check well conditions
/// \internal[check well conditions]
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
/// \internal[check well conditions]
/// \endinternal
/// \page tutorial4
/// \details Transport solver
/// \TODO We must call computeTransportSource() here, since we have wells.
/// \code
/// \snippet tutorial4.cpp tranport solver
/// \internal[tranport solver]
transport_solver.solve(&state.faceflux()[0], &porevol[0], &src[0], dt, state.saturation());
/// \endcode
/// \internal[tranport solver]
/// \endinternal
/// \page tutorial4
/// \details Write the output to file.
/// \code
/// \snippet tutorial4.cpp write output
/// \internal[write output]
vtkfilename.str("");
vtkfilename << "tutorial4-" << std::setw(3) << std::setfill('0') << i << ".vtu";
std::ofstream vtkfile(vtkfilename.str().c_str());
@@ -360,7 +439,8 @@ int main ()
destroy_wells(wells);
}
/// \endcode
/// \internal[write output]
/// \endinternal
@@ -386,4 +466,8 @@ int main ()
/// \details
/// \section completecode4 Complete source code:
/// \include tutorial4.cpp
/// \include generate_doc_figures.py
/// \page tutorial4
/// \details
/// \section pythonscript4 python script to generate figures:
/// \snippet generate_doc_figures.py tutorial4