From febdc04e806b7aa5faa86efc8a649549f6374102 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Wed, 6 Dec 2017 12:07:33 +0100 Subject: [PATCH] remove tutorials which have been moved to opm-simulators --- CMakeLists.txt | 5 - CMakeLists_files.cmake | 6 +- tutorials/CMakeLists.txt | 17 -- tutorials/generate_doc_figures.py | 151 --------- tutorials/tutorial1.cpp | 134 -------- tutorials/tutorial2.cpp | 222 -------------- tutorials/tutorial3.cpp | 364 ---------------------- tutorials/tutorial4.cpp | 490 ------------------------------ 8 files changed, 1 insertion(+), 1388 deletions(-) delete mode 100644 tutorials/CMakeLists.txt delete mode 100755 tutorials/generate_doc_figures.py delete mode 100644 tutorials/tutorial1.cpp delete mode 100644 tutorials/tutorial2.cpp delete mode 100644 tutorials/tutorial3.cpp delete mode 100644 tutorials/tutorial4.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index d0363df5..1cef8aaf 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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) diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index c0282329..2779f49b 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -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: diff --git a/tutorials/CMakeLists.txt b/tutorials/CMakeLists.txt deleted file mode 100644 index ac3baa54..00000000 --- a/tutorials/CMakeLists.txt +++ /dev/null @@ -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() diff --git a/tutorials/generate_doc_figures.py b/tutorials/generate_doc_figures.py deleted file mode 100755 index 7afe8b51..00000000 --- a/tutorials/generate_doc_figures.py +++ /dev/null @@ -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) - diff --git a/tutorials/tutorial1.cpp b/tutorials/tutorial1.cpp deleted file mode 100644 index 43216d17..00000000 --- a/tutorials/tutorial1.cpp +++ /dev/null @@ -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 . -*/ -///\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 -#include -// 17.03.2016 Temporarily removed while moving functionality to opm-output -#ifdef DISABLE_OUTPUT -#include -#endif -#include -#include -#include -/// \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 nx, ny, nz, dx, dy, dz - /// 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 - diff --git a/tutorials/tutorial2.cpp b/tutorials/tutorial2.cpp deleted file mode 100644 index 84984f63..00000000 --- a/tutorials/tutorial2.cpp +++ /dev/null @@ -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 . -*/ -/// \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 -#include -// 17.03.2016 Temporarily removed while moving functionality to opm-output -#ifdef DISABLE_OUTPUT -#include -#endif -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -/// \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 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 viscosities(num_phases, 1.0*centi*Poise); - std::vector 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 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 - /// solve() method of class - /// Opm::IncompTPFA. - /// \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 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 diff --git a/tutorials/tutorial3.cpp b/tutorials/tutorial3.cpp deleted file mode 100644 index e8961584..00000000 --- a/tutorials/tutorial3.cpp +++ /dev/null @@ -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 . -*/ -/// \endcond -#ifdef HAVE_CONFIG_H -#include "config.h" -#endif - -#include -#include -#include -#include -#include -#include -#include -// 17.03.2016 Temporarily removed while moving functionality to opm-output -#ifdef DISABLE_OUTPUT -#include -#endif -#include -#include -#include -#include - -#include - -#include -#include -#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 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 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_w}{\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 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 density(num_phases, 1000.0); - std::vector 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 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 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 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 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. -/// -/// -/// -/// -/// -/// -/// -/// -/// -/// -/// -/// -/// -///
\image html tutorial3-000.png \image html tutorial3-005.png
\image html tutorial3-010.png \image html tutorial3-015.png
\image html tutorial3-019.png
- - -/// \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 diff --git a/tutorials/tutorial4.cpp b/tutorials/tutorial4.cpp deleted file mode 100644 index e37e86f1..00000000 --- a/tutorials/tutorial4.cpp +++ /dev/null @@ -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 . -*/ -/// \endcond -#ifdef HAVE_CONFIG_H -#include "config.h" -#endif - -#include -#include -#include -#include -#include -#include -#include -// 17.03.2016 Temporarily removed while moving functionality to opm-output -#ifdef DISABLE_OUTPUT -#include -#endif -#include -#include -#include -#include - -#include - -#include -#include -#include - -#include -#include -#include -#include - -/// \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 rho(2, 1000.); - /// \internal[density] - /// \endinternal - - /// \page tutorial4 - /// \details viscosity vector (one component per phase). - /// \snippet tutorial4.cpp viscosity - /// \internal[viscosity] - std::vector 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 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 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 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 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 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 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 well_resflowrates_phase; - std::vector well_surflowrates_phase; - std::vector 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. -/// -/// -/// -/// -/// -/// -/// -/// -/// -/// -/// -/// -/// -///
\image html tutorial4-000.png \image html tutorial4-005.png
\image html tutorial4-010.png \image html tutorial4-015.png
\image html tutorial4-019.png
- - -/// \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