remove tutorials which have been moved to opm-simulators

This commit is contained in:
Arne Morten Kvarving 2017-12-06 12:07:33 +01:00
parent efc2915b74
commit febdc04e80
8 changed files with 1 additions and 1388 deletions

View File

@ -86,11 +86,6 @@ macro (sources_hook)
${PROJECT_SOURCE_DIR}/${opm-core_DIR}/core/linalg/call_umfpack.c
${PROJECT_SOURCE_DIR}/${opm-core_DIR}/core/linalg/LinearSolverUmfpack.cpp
)
list (REMOVE_ITEM examples_SOURCES
${PROJECT_SOURCE_DIR}/tutorials/tutorial2.cpp
${PROJECT_SOURCE_DIR}/tutorials/tutorial3.cpp
${PROJECT_SOURCE_DIR}/tutorials/tutorial4.cpp
)
endif (NOT SuiteSparse_FOUND)
if (NOT PETSC_FOUND)

View File

@ -156,17 +156,13 @@ list (APPEND TEST_DATA_FILES
)
# originally generated with the command:
# find tutorials examples -name '*.c*' -printf '\t%p\n' | sort
# find examples -name '*.c*' -printf '\t%p\n' | sort
list (APPEND EXAMPLE_SOURCE_FILES
examples/compute_eikonal_from_files.cpp
examples/compute_initial_state.cpp
examples/compute_tof.cpp
examples/compute_tof_from_files.cpp
examples/diagnose_relperm.cpp
tutorials/tutorial1.cpp
tutorials/tutorial2.cpp
tutorials/tutorial3.cpp
tutorials/tutorial4.cpp
)
# originally generated with the command:

View File

@ -1,17 +0,0 @@
# Build system for opm-core tutorials
# Note: This is independent of the main opm-core build system.
# It assumes an opm-core library is installed at the system level.
project(opm-core-tutorials)
cmake_minimum_required(VERSION 2.8)
find_package(opm-core REQUIRED)
include_directories(${opm-core_INCLUDE_DIRS})
add_definitions(${opm-core_DEFINITIONS})
foreach(tut 1 2 3 4)
add_executable(tutorial${tut} tutorial${tut}.cpp)
target_link_libraries(tutorial${tut} ${opm-core_LIBRARIES})
endforeach()

View File

@ -1,151 +0,0 @@
# This script generate the illustration pictures for the documentation.
#
# To run this script, you have to install paraview, see:
#
# http://www.paraview.org/paraview/resources/software.php
#
# Eventually, set up the paths (figure_path, tutorial_data_path) according to your own installation.
# (The default values should be ok.)
#
# Make sure that pvpython is in your path of executables.
#
# After all the tutorial programs have been executed, run the following
# command in the same directory:
#
# pvpython generate_doc_figures.py Documentation/Figure
#
from paraview.simple import *
from os import remove, mkdir, curdir
from os.path import join, isdir
from sys import argv, exit
# we need at least the output directory
if len(argv) <= 1:
exit('Synopsis: pvpython generate_doc_figures.py dest-dir [src-dir]')
figure_path = argv[1]
# default for the input directory is the current one
if len(argv) <= 2:
tutorial_data_path = curdir
else:
tutorial_data_path = argv[2]
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)
grid = XMLUnstructuredGridReader(FileName = data_file_name)
collected_garbage_file.append(data_file_name)
grid.UpdatePipeline()
Show(grid)
dp = GetDisplayProperties(grid)
dp.Representation = 'Wireframe'
dp.LineWidth = 5
dp.AmbientColor = [1, 0, 0]
view = GetActiveView()
view.Background = [1, 1, 1]
camera = GetActiveCamera()
camera.SetPosition(4, -6, 5)
camera.SetViewUp(-0.19, 0.4, 0.9)
camera.SetViewAngle(30)
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)
collected_garbage_file.append(data_file_name)
grid.UpdatePipeline()
Show(grid)
dp = GetDisplayProperties(grid)
dp.Representation = 'Surface'
dp.ColorArrayName = 'pressure'
pres = grid.CellData.GetArray(0)
pres_lookuptable = GetLookupTableForArray( "pressure", 1, RGBPoints=[pres.GetRange()[0], 1, 0, 0, pres.GetRange()[1], 0, 0, 1] )
dp.LookupTable = pres_lookuptable
view = GetActiveView()
view.Background = [1, 1, 1]
camera = GetActiveCamera()
camera.SetPosition(20, 20, 110)
camera.SetViewUp(0, 1, 0)
camera.SetViewAngle(30)
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")
collected_garbage_file.append(data_file_name)
cases = ["000", "005", "010", "015", "019"]
for case in cases:
data_file_name = join(tutorial_data_path, "tutorial3-"+case+".vtu")
grid = XMLUnstructuredGridReader(FileName = data_file_name)
grid.UpdatePipeline()
Show(grid)
dp = GetDisplayProperties(grid)
dp.Representation = 'Surface'
dp.ColorArrayName = 'saturation'
sat = grid.CellData.GetArray(1)
sat_lookuptable = GetLookupTableForArray( "saturation", 1, RGBPoints=[0, 1, 0, 0, 1, 0, 0, 1])
dp.LookupTable = sat_lookuptable
view.Background = [1, 1, 1]
camera = GetActiveCamera()
camera.SetPosition(100, 100, 550)
camera.SetViewUp(0, 1, 0)
camera.SetViewAngle(30)
camera.SetFocalPoint(100, 100, 5)
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")
collected_garbage_file.append(data_file_name)
cases = ["000", "005", "010", "015", "019"]
for case in cases:
data_file_name = join(tutorial_data_path, "tutorial4-"+case+".vtu")
grid = XMLUnstructuredGridReader(FileName = data_file_name)
grid.UpdatePipeline()
Show(grid)
dp = GetDisplayProperties(grid)
dp.Representation = 'Surface'
dp.ColorArrayName = 'saturation'
sat = grid.CellData.GetArray(1)
sat_lookuptable = GetLookupTableForArray( "saturation", 1, RGBPoints=[0, 1, 0, 0, 1, 0, 0, 1])
dp.LookupTable = sat_lookuptable
view.Background = [1, 1, 1]
camera = GetActiveCamera()
camera.SetPosition(100, 100, 550)
camera.SetViewUp(0, 1, 0)
camera.SetViewAngle(30)
camera.SetFocalPoint(100, 100, 5)
Render()
WriteImage(join(figure_path, "tutorial4-"+case+".png"))
Hide(grid)
## [tutorial4]
# remove temporary files
for f in collected_garbage_file:
remove(f)

View File

@ -1,134 +0,0 @@
/// \cond SKIP
/*!
Copyright 2012 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/>.
*/
///\endcond
/// \page tutorial1 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 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.
/// \snippet tutorial1.cpp including headers
/// \internal [including headers]
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <opm/core/grid.h>
#include <opm/core/grid/GridManager.hpp>
// 17.03.2016 Temporarily removed while moving functionality to opm-output
#ifdef DISABLE_OUTPUT
#include <opm/core/io/vtk/writeVtkData.hpp>
#endif
#include <iostream>
#include <fstream>
#include <vector>
/// \internal [including headers]
/// \endinternal
// ----------------- Main program -----------------
int main()
try
{
/// \page tutorial1
/// We set the number of blocks in each direction.
/// \snippet tutorial1.cpp num blocks
/// \internal [num blocks]
int nx = 4;
int ny = 3;
int nz = 2;
/// \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.
/// \snippet tutorial1.cpp dim
/// \internal [dim]
double dx = 1.0;
double dy = 1.0;
double dz = 1.0;
/// \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
/// destroy the data structure. In this tutorial however, we will use Opm::GridManager,
/// 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.
/// \snippet tutorial1.cpp grid manager
/// \internal [grid manager]
Opm::GridManager grid(nx, ny, nz, dx, dy, dz);
/// \internal [grid manager]
/// \endinternal
/// \page tutorial1
/// We open an output file stream for the output
/// \snippet tutorial1.cpp output stream
/// \internal [output stream]
std::ofstream vtkfile("tutorial1.vtu");
/// \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
/// \snippet tutorial1.cpp data map
/// \internal [data map]
// 17.03.2016 Temporarily removed while moving functionality to opm-output
#ifdef DISABLE_OUTPUT
Opm::DataMap dm;
#endif
/// \internal [data map]
/// \endinternal
/// \page tutorial1
/// Call Opm::writeVtkData() to write the output file.
/// \snippet tutorial1.cpp write vtk
/// \internal [write vtk]
// 17.03.2016 Temporarily removed while moving functionality to opm-output
#ifdef DISABLE_OUTPUT
Opm::writeVtkData(*grid.c_grid(), dm, vtkfile);
#endif
/// \internal [write vtk]
/// \endinternal
}
catch (const std::exception &e) {
std::cerr << "Program threw an exception: " << e.what() << "\n";
throw;
}
/// \page tutorial1
/// We read the vtu output file in \a Paraview and obtain the following grid.
/// \image html tutorial1.png
/// \page tutorial1
/// \section completecode1 Complete source code:
/// \include tutorial1.cpp
/// \page tutorial1
/// \details
/// \section pythonscript1 Python script to generate figures:
/// \snippet generate_doc_figures.py tutorial1

View File

@ -1,222 +0,0 @@
/// \cond SKIP
/*!
Copyright 2012 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/>.
*/
/// \endcond
/// \page tutorial2 Flow Solver for a single phase
/// \details The flow equations consist of the mass conservation equation
/// \f[\nabla\cdot {\bf u}=q\f] and the Darcy law \f[{\bf u} =- \frac{1}{\mu}K\nabla p.\f] Here,
/// \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
/// \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).
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <opm/core/grid.h>
#include <opm/core/grid/GridManager.hpp>
// 17.03.2016 Temporarily removed while moving functionality to opm-output
#ifdef DISABLE_OUTPUT
#include <opm/core/io/vtk/writeVtkData.hpp>
#endif
#include <iostream>
#include <fstream>
#include <vector>
#include <opm/core/props/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/parser/eclipse/Units/Units.hpp>
#include <opm/core/simulator/TwophaseState.hpp>
#include <opm/core/simulator/WellState.hpp>
/// \page tutorial2
/// \section commentedcode2 Program walk-through.
///
int main()
try
{
/// \page tutorial2
/// 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);
/// \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.
/// \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;
/// \internal [access grid]
/// endinternal
/// \page tutorial2
/// \details
/// We define a fluid viscosity equal to 1 cP and density equal
/// to 1000 kg/m^3.
/// The <opm/parser/eclipse/Units/Units.hpp> header contains support
/// for common units and prefixes, in the namespaces Opm::unit
/// and Opm::prefix.
/// \snippet tutorial2.cpp fluid
/// \internal [fluid]
using namespace Opm::unit;
using namespace Opm::prefix;
int num_phases = 1;
std::vector<double> viscosities(num_phases, 1.0*centi*Poise);
std::vector<double> densities(num_phases, 1000.0*kilogram/cubic(meter));
/// \internal [fluid]
/// \endinternal
/// \page tutorial2
/// \details
/// We define a permeability equal to 100 mD.
/// \snippet tutorial2.cpp perm
/// \internal [perm]
double permeability = 100.0*milli*darcy;
/// \internal [perm]
/// \endinternal
/// \page tutorial2
/// \details
/// We set up a simple property object for a single-phase situation.
/// \snippet tutorial2.cpp single-phase property
/// \internal [single-phase property]
const double porosity = 1.;
Opm::IncompPropertiesBasic props(1, Opm::SaturationPropsBasic::Constant,
densities, viscosities, porosity,
permeability, dim, num_cells);
/// \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).
/// \snippet tutorial2.cpp linsolver
/// \internal [linsolver]
Opm::LinearSolverUmfpack linsolver;
/// \internal [linsolver]
/// \endinternal
/// \page tutorial2
/// We define the source term.
/// \snippet tutorial2.cpp source
/// \internal [source]
std::vector<double> src(num_cells, 0.0);
src[0] = 150.*cubic(meter)/day;
src[num_cells-1] = -src[0];
/// \internal [source]
/// \endinternal
/// \page tutorial2
/// \details We set up the boundary conditions.
/// By default, we obtain no-flow boundary conditions.
/// \snippet tutorial2.cpp boundary
/// \internal [boundary]
Opm::FlowBCManager bcs;
/// \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).
/// \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>.
/// \snippet tutorial2.cpp state
/// \internal [state]
Opm::TwophaseState state( num_cells , num_faces );
Opm::WellState well_state;
/// \internal [state]
/// \endinternal
/// \page tutorial2
/// We call the pressure solver.
/// The first (timestep) argument does not matter for this
/// incompressible case.
/// \snippet tutorial2.cpp pressure solver
/// \internal [pressure solver]
psolver.solve(1.0*day, state, well_state);
/// \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).
/// \snippet tutorial2.cpp write output
/// \internal [write output]
// 17.03.2016 Temporarily removed while moving functionality to opm-output
#ifdef DISABLE_OUTPUT
std::ofstream vtkfile("tutorial2.vtu");
Opm::DataMap dm;
dm["pressure"] = &state.pressure();
std::vector<double> cell_velocity;
Opm::estimateCellVelocity(*grid.c_grid(), state.faceflux(), cell_velocity);
dm["velocity"] = &cell_velocity;
Opm::writeVtkData(*grid.c_grid(), dm, vtkfile);
#endif
/// \internal [write output]
/// \endinternal
}
catch (const std::exception &e) {
std::cerr << "Program threw an exception: " << e.what() << "\n";
throw;
}
/// \page tutorial2
/// We read the vtu output file in \a Paraview and obtain the following pressure
/// distribution. \image html tutorial2.png
/// \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

View File

@ -1,364 +0,0 @@
/// \cond SKIP
/*!
Copyright 2012 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/>.
*/
/// \endcond
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <iostream>
#include <iomanip>
#include <fstream>
#include <vector>
#include <cassert>
#include <opm/core/grid.h>
#include <opm/core/grid/GridManager.hpp>
// 17.03.2016 Temporarily removed while moving functionality to opm-output
#ifdef DISABLE_OUTPUT
#include <opm/core/io/vtk/writeVtkData.hpp>
#endif
#include <opm/core/linalg/LinearSolverUmfpack.hpp>
#include <opm/core/pressure/IncompTpfa.hpp>
#include <opm/core/pressure/FlowBCManager.hpp>
#include <opm/core/props/IncompPropertiesBasic.hpp>
#include <opm/core/transport/reorder/TransportSolverTwophaseReorder.hpp>
#include <opm/core/simulator/initState.hpp>
#include <opm/core/simulator/TwophaseState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/parser/eclipse/Units/Units.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
/// \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 rock and both fluid phases are incompressible. Since
/// \f$s_w+s_o=1\f$, we may add the conservation equations to 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 <strong>pressure equation</strong>. 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_w}{\partial t}+\nabla\cdot(f_w u)=\frac{q_w}{\rho_w}\f]
/// which is referred to as the <strong>transport equation</strong>. 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 commentedsource1 Program walk-through.
/// \details
/// Main function
/// \snippet tutorial3.cpp main
/// \internal [main]
int main ()
try
{
/// \internal [main]
/// \endinternal
/// \page tutorial3
/// \details
/// 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.
///
/// The Opm::GridManager is responsible for creating and destroying the grid,
/// the UnstructuredGrid data structure contains the actual grid topology
/// and geometry.
/// \snippet tutorial3.cpp grid
/// \internal [grid]
int nx = 20;
int ny = 20;
int nz = 1;
double dx = 10.0;
double dy = 10.0;
double dz = 10.0;
using namespace Opm;
GridManager grid_manager(nx, ny, nz, dx, dy, dz);
const UnstructuredGrid& grid = *grid_manager.c_grid();
int num_cells = grid.number_of_cells;
/// \internal [grid]
/// \endinternal
/// \page tutorial3
/// \details
/// We define the properties of the fluid.\n
/// Number of phases, phase densities, phase viscosities,
/// rock porosity and permeability.
///
/// We always use SI units in the simulator. Many units are
/// 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.
/// \snippet tutorial3.cpp set properties
/// \internal [set properties]
int num_phases = 2;
using namespace Opm::unit;
using namespace Opm::prefix;
std::vector<double> density(num_phases, 1000.0);
std::vector<double> viscosity(num_phases, 1.0*centi*Poise);
double porosity = 0.5;
double permeability = 10.0*milli*darcy;
/// \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.
/// \snippet tutorial3.cpp relperm
/// \internal [relperm]
SaturationPropsBasic::RelPermFunc rel_perm_func = SaturationPropsBasic::Linear;
/// \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.
/// \snippet tutorial3.cpp properties
/// \internal [properties]
IncompPropertiesBasic props(num_phases, rel_perm_func, density, viscosity,
porosity, permeability, grid.dimensions, num_cells);
/// \internal [properties]
/// \endinternal
/// \page tutorial3
/// \details Gravity parameters. Here, we set zero gravity.
/// \snippet tutorial3.cpp gravity
/// \internal [gravity]
const double *grav = 0;
std::vector<double> omega;
/// \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.
/// \snippet tutorial3.cpp source
/// \internal [source]
std::vector<double> src(num_cells, 0.0);
src[0] = 1.;
src[num_cells-1] = -1.;
/// \internal [source]
/// \endinternal
/// \page tutorial3
/// \details We set up the boundary conditions. Letting bcs be empty is equivalent
/// to no-flow boundary conditions.
/// \snippet tutorial3.cpp boundary
/// \internal [boundary]
FlowBCManager bcs;
/// \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.
/// \snippet tutorial3.cpp pressure solver
/// \internal [pressure solver]
LinearSolverUmfpack linsolver;
IncompTpfa psolver(grid, props, linsolver, grav, NULL, src, bcs.c_bcs());
/// \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.
/// \snippet tutorial3.cpp well
/// \internal [well]
WellState well_state;
/// \internal [well]
/// \endinternal
/// \page tutorial3
/// \details We compute the pore volume
/// \snippet tutorial3.cpp pore volume
/// \internal [pore volume]
std::vector<double> porevol;
Opm::computePorevolume(grid, props.porosity(), porevol);
/// \internal [pore volume]
/// \endinternal
/// \page tutorial3
/// \details Set up the transport solver. This is a reordering implicit Euler transport solver.
/// \snippet tutorial3.cpp transport solver
/// \internal [transport solver]
const double tolerance = 1e-9;
const int max_iterations = 30;
Opm::TransportSolverTwophaseReorder transport_solver(grid, props, NULL, tolerance, max_iterations);
/// \internal [transport solver]
/// \endinternal
/// \page tutorial3
/// \details Time integration parameters
/// \snippet tutorial3.cpp time parameters
/// \internal [time parameters]
const double dt = 0.1*day;
const int num_time_steps = 20;
/// \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.
/// \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;
}
/// \internal [cell indexes]
/// \endinternal
/// \page tutorial3
/// \details
/// We set up a two-phase state object, and
/// initialize water saturation to minimum everywhere.
/// \snippet tutorial3.cpp two-phase state
/// \internal [two-phase state]
TwophaseState state( grid.number_of_cells , grid.number_of_faces );
initSaturation( allcells , props , state , MinSat );
/// \internal [two-phase state]
/// \endinternal
/// \page tutorial3
/// \details This string stream will be used to construct a new
/// output filename at each timestep.
/// \snippet tutorial3.cpp output stream
/// \internal [output stream]
std::ostringstream vtkfilename;
/// \internal [output stream]
/// \endinternal
/// \page tutorial3
/// \details Loop over the time steps.
/// \snippet tutorial3.cpp time loop
/// \internal [time loop]
for (int i = 0; i < num_time_steps; ++i) {
/// \internal [time loop]
/// \endinternal
/// \page tutorial3
/// \details Solve the pressure equation
/// \snippet tutorial3.cpp solve pressure
/// \internal [solve pressure]
psolver.solve(dt, state, well_state);
/// \internal [solve pressure]
/// \endinternal
/// \page tutorial3
/// \details Solve the transport equation.
/// \snippet tutorial3.cpp transport solve
/// \internal [transport solve]
transport_solver.solve(&porevol[0], &src[0], dt, state);
/// \internal [transport solve]
/// \endinternal
/// \page tutorial3
/// \details Write the output to file.
/// \snippet tutorial3.cpp write output
/// \internal [write output]
vtkfilename.str("");
vtkfilename << "tutorial3-" << std::setw(3) << std::setfill('0') << i << ".vtu";
// 17.03.2016 Temporarily removed while moving functionality to opm-output
#ifdef DISABLE_OUTPUT
std::ofstream vtkfile(vtkfilename.str().c_str());
Opm::DataMap dm;
dm["saturation"] = &state.saturation();
dm["pressure"] = &state.pressure();
Opm::writeVtkData(grid, dm, vtkfile);
#endif
}
}
catch (const std::exception &e) {
std::cerr << "Program threw an exception: " << e.what() << "\n";
throw;
}
/// \internal [write output]
/// \endinternal
/// \page tutorial3
/// \section results3 Results.
/// <TABLE>
/// <TR>
/// <TD> \image html tutorial3-000.png </TD>
/// <TD> \image html tutorial3-005.png </TD>
/// </TR>
/// <TR>
/// <TD> \image html tutorial3-010.png </TD>
/// <TD> \image html tutorial3-015.png </TD>
/// </TR>
/// <TR>
/// <TD> \image html tutorial3-019.png </TD>
/// <TD> </TD>
/// </TR>
/// </TABLE>
/// \page tutorial3
/// \details
/// \section completecode3 Complete source code:
/// \include tutorial3.cpp
/// \page tutorial3
/// \details
/// \section pythonscript3 python script to generate figures:
/// \snippet generate_doc_figures.py tutorial3

View File

@ -1,490 +0,0 @@
/// \cond SKIP
/*!
Copyright 2012 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/>.
*/
/// \endcond
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <iostream>
#include <iomanip>
#include <fstream>
#include <vector>
#include <cassert>
#include <opm/core/grid.h>
#include <opm/core/grid/GridManager.hpp>
// 17.03.2016 Temporarily removed while moving functionality to opm-output
#ifdef DISABLE_OUTPUT
#include <opm/core/io/vtk/writeVtkData.hpp>
#endif
#include <opm/core/linalg/LinearSolverUmfpack.hpp>
#include <opm/core/pressure/IncompTpfa.hpp>
#include <opm/core/pressure/FlowBCManager.hpp>
#include <opm/core/props/IncompPropertiesBasic.hpp>
#include <opm/core/transport/reorder/TransportSolverTwophaseReorder.hpp>
#include <opm/core/simulator/initState.hpp>
#include <opm/core/simulator/TwophaseState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/parser/eclipse/Units/Units.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#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
/// \snippet tutorial4.cpp main
/// \internal[main]
int main ()
try
{
/// \internal[main]
/// \endinternal
/// \page tutorial4
/// \details
/// 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;
int nz = 1;
double dx = 10.;
double dy = 10.;
double dz = 10.;
using namespace Opm;
GridManager grid_manager(nx, ny, nz, dx, dy, dz);
const UnstructuredGrid& grid = *grid_manager.c_grid();
int num_cells = grid.number_of_cells;
/// \internal[cartesian grid]
/// \endinternal
/// \page tutorial4
/// \details
/// We define the properties of the fluid.\n
/// Number of phases.
/// \snippet tutorial4.cpp Number of phases
/// \internal[Number of phases]
int num_phases = 2;
using namespace unit;
using namespace prefix;
/// \internal[Number of phases]
/// \endinternal
/// \page tutorial4
/// \details density vector (one component per phase).
/// \snippet tutorial4.cpp density
/// \internal[density]
std::vector<double> rho(2, 1000.);
/// \internal[density]
/// \endinternal
/// \page tutorial4
/// \details viscosity vector (one component per phase).
/// \snippet tutorial4.cpp viscosity
/// \internal[viscosity]
std::vector<double> mu(2, 1.*centi*Poise);
/// \internal[viscosity]
/// \endinternal
/// \page tutorial4
/// \details porosity and permeability of the rock.
/// \snippet tutorial4.cpp rock
/// \internal[rock]
double porosity = 0.5;
double k = 10*milli*darcy;
/// \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.
/// \snippet tutorial4.cpp relative permeability
/// \internal[relative permeability]
SaturationPropsBasic::RelPermFunc rel_perm_func = SaturationPropsBasic::Linear;
/// \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.
/// \snippet tutorial4.cpp fluid properties
/// \internal[fluid properties]
IncompPropertiesBasic props(num_phases, rel_perm_func, rho, mu,
porosity, k, dim, num_cells);
/// \internal[fluid properties]
/// \endinternal
/// \page tutorial4
/// \details Gravity parameters. Here, we set zero gravity.
/// \snippet tutorial4.cpp Gravity
/// \internal[Gravity]
const double *grav = 0;
std::vector<double> omega;
/// \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.
/// \snippet tutorial4.cpp source
/// \internal[source]
std::vector<double> src(num_cells, 0.0);
src[0] = 1.;
src[num_cells-1] = -1.;
/// \internal[source]
/// \endinternal
/// \page tutorial4
/// \details We compute the pore volume
/// \snippet tutorial4.cpp pore volume
/// \internal[pore volume]
std::vector<double> porevol;
Opm::computePorevolume(grid, props.porosity(), porevol);
/// \internal[pore volume]
/// \endinternal
/// \page tutorial4
/// \details Set up the transport solver. This is a reordering implicit Euler transport solver.
/// \snippet tutorial4.cpp transport solver
/// \internal[transport solver]
const double tolerance = 1e-9;
const int max_iterations = 30;
Opm::TransportSolverTwophaseReorder transport_solver(grid, props, NULL, tolerance, max_iterations);
/// \internal[transport solver]
/// \endinternal
/// \page tutorial4
/// \details Time integration parameters
/// \snippet tutorial4.cpp Time integration
/// \internal[Time integration]
double dt = 0.1*day;
int num_time_steps = 20;
/// \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.
/// \snippet tutorial4.cpp boundary
/// \internal[boundary]
FlowBCManager bcs;
/// \internal[boundary]
/// \endinternal
/// \page tutorial4
/// \details
/// We set up a two-phase state object, and
/// initialise water saturation to minimum everywhere.
/// \snippet tutorial4.cpp two-phase state
/// \internal[two-phase state]
TwophaseState state( grid.number_of_cells , grid.number_of_faces );
initSaturation( allcells , props , state , MinSat );
/// \internal[two-phase state]
/// \endinternal
/// \page tutorial4
/// \details This string will contain the name of a VTK output vector.
/// \snippet tutorial4.cpp VTK output
/// \internal[VTK output]
std::ostringstream vtkfilename;
/// \internal[VTK output]
/// \endinternal
/// \page tutorial4
/// To create wells we need an instance of the PhaseUsage-object
/// \snippet tutorial4.cpp PhaseUsage-object
/// \internal[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;
/// \internal[PhaseUsage-object]
/// \endinternal
/// \page tutorial4
/// \details This will contain our well-specific information
/// \snippet tutorial4.cpp well_collection
/// \internal[well_collection]
WellCollection well_collection;
/// \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.
/// \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;
/// \internal[production specification]
/// \endinternal
const double group_efficiency_factor = 0.9;
/// \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.
/// \snippet tutorial4.cpp injection specification
/// \internal[injection specification]
std::shared_ptr<WellsGroupInterface> well_group(new WellsGroup("group", group_efficiency_factor, well_group_prod_spec,
InjectionSpecification(), phase_usage));
/// \internal[injection specification]
/// \endinternal
/// \page tutorial4
/// \details We add our well_group to the well_collection
/// \snippet tutorial4.cpp well_collection
/// \internal[well_collection]
well_collection.addChild(well_group);
/// \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.
/// \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;
well_name << "well" << i;
ProductionSpecification production_specification;
production_specification.control_mode_ = ProductionSpecification::GRUP;
const double efficiency_factor = 0.8;
std::shared_ptr<WellsGroupInterface> well_leaf_node(new WellNode(well_name.str(), efficiency_factor, production_specification,
InjectionSpecification(), phase_usage));
well_collection.addChild(well_leaf_node, "group");
}
/// \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).
/// \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*/);
/// \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).
/// \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;
const int sat_table_id = -1;
std::stringstream well_name;
well_name << "well" << i;
bool allowCrossFlow = true;
add_well(PRODUCER, 0, 1, NULL, &well_cells, &well_index, &sat_table_id,
well_name.str().c_str(), allowCrossFlow, wells);
}
/// \internal[well cells]
/// \endinternal
/// \page tutorial4
/// \details We need to make the well collection aware of our wells object
/// \snippet tutorial4.cpp set well pointer
/// \internal[set well pointer]
well_collection.setWellsPointer(wells);
/// \internal[set well pointer]
/// \endinternal
/// \page tutorial4
/// We're not using well controls, just group controls, so we need to apply them.
/// \snippet tutorial4.cpp apply group controls
/// \internal[apply group controls]
well_collection.applyGroupControls();
/// \internal[apply group controls]
/// \endinternal
/// \page tutorial4
/// \details We set up necessary information for the wells
/// \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;
/// \internal[init wells]
/// \endinternal
/// \page tutorial4
/// \details We set up the pressure solver.
/// \snippet tutorial4.cpp pressure solver
/// \internal[pressure solver]
LinearSolverUmfpack linsolver;
IncompTpfa psolver(grid, props, linsolver,
grav, wells, src, bcs.c_bcs());
/// \internal[pressure solver]
/// \endinternal
/// \page tutorial4
/// \details Loop over the time steps.
/// \snippet tutorial4.cpp time loop
/// \internal[time loop]
for (int i = 0; i < num_time_steps; ++i) {
/// \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.
/// \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) {
/// \internal[well iterations]
/// \endinternal
/// \page tutorial4
/// \details Solve the pressure equation
/// \snippet tutorial4.cpp pressure solve
/// \internal[pressure solve]
psolver.solve(dt, state, well_state);
/// \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);
/// \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) {
OPM_THROW(std::runtime_error, "Conditions not met within " << max_well_iterations<< " iterations.");
}
}
/// \internal[check well conditions]
/// \endinternal
/// \page tutorial4
/// \details Transport solver
/// \TODO We must call computeTransportSource() here, since we have wells.
/// \snippet tutorial4.cpp tranport solver
/// \internal[tranport solver]
transport_solver.solve(&porevol[0], &src[0], dt, state);
/// \internal[tranport solver]
/// \endinternal
/// \page tutorial4
/// \details Write the output to file.
/// \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());
// 17.03.2016 Temporarily removed while moving functionality to opm-output
#ifdef DISABLE_OUTPUT
Opm::DataMap dm;
dm["saturation"] = &state.saturation();
dm["pressure"] = &state.pressure();
Opm::writeVtkData(grid, dm, vtkfile);
#endif
}
destroy_wells(wells);
}
catch (const std::exception &e) {
std::cerr << "Program threw an exception: " << e.what() << "\n";
throw;
}
/// \internal[write output]
/// \endinternal
/// \page tutorial4
/// \section results4 Results.
/// <TABLE>
/// <TR>
/// <TD> \image html tutorial4-000.png </TD>
/// <TD> \image html tutorial4-005.png </TD>
/// </TR>
/// <TR>
/// <TD> \image html tutorial4-010.png </TD>
/// <TD> \image html tutorial4-015.png </TD>
/// </TR>
/// <TR>
/// <TD> \image html tutorial4-019.png </TD>
/// <TD> </TD>
/// </TR>
/// </TABLE>
/// \page tutorial4
/// \details
/// \section completecode4 Complete source code:
/// \include tutorial4.cpp
/// \page tutorial4
/// \details
/// \section pythonscript4 python script to generate figures:
/// \snippet generate_doc_figures.py tutorial4