diff --git a/Doxyfile b/Doxyfile index 4b1c43ee..3f1fb256 100644 --- a/Doxyfile +++ b/Doxyfile @@ -610,8 +610,8 @@ WARN_LOGFILE = # directories like "/usr/src/myproject". Separate the files or directories # with spaces. -#INPUT = opm/core tutorials examples -INPUT = tutorials +INPUT = opm/core tutorials/tutorial1.cpp tutorials/tutorial2.cpp tutorials/tutorial3.cpp examples +#INPUT = tutorials # This tag can be used to specify the character encoding of the source files # that doxygen parses. Internally doxygen uses the UTF-8 encoding, which is diff --git a/Makefile.am b/Makefile.am index a445e380..3d79e38f 100644 --- a/Makefile.am +++ b/Makefile.am @@ -51,6 +51,8 @@ opm/core/fluid/RockFromDeck.cpp \ opm/core/fluid/SaturationPropsBasic.cpp \ opm/core/fluid/SaturationPropsFromDeck.cpp \ opm/core/grid.c \ +opm/core/newwells.c \ +opm/core/GridManager.cpp \ opm/core/grid/cart_grid.c \ opm/core/grid/cornerpoint_grid.c \ opm/core/grid/cpgpreprocess/geometry.c \ @@ -72,14 +74,12 @@ opm/core/utility/miscUtilities.cpp \ opm/core/utility/miscUtilitiesBlackoil.cpp \ opm/core/utility/SimulatorTimer.cpp \ opm/core/utility/StopWatch.cpp \ -opm/core/utility/newwells.c \ opm/core/utility/writeVtkData.cpp \ -opm/core/GridManager.cpp \ -opm/core/WellsManager.cpp \ -opm/core/WellsGroup.cpp \ -opm/core/WellCollection.cpp \ -opm/core/InjectionSpecification.cpp \ -opm/core/ProductionSpecification.cpp \ +opm/core/wells/WellsManager.cpp \ +opm/core/wells/WellsGroup.cpp \ +opm/core/wells/WellCollection.cpp \ +opm/core/wells/InjectionSpecification.cpp \ +opm/core/wells/ProductionSpecification.cpp \ opm/core/linalg/sparse_sys.c \ opm/core/linalg/LinearSolverInterface.cpp \ opm/core/linalg/LinearSolverFactory.cpp \ @@ -109,6 +109,7 @@ opm/core/pressure/ifsh.c \ opm/core/pressure/mimetic/mimetic.c \ opm/core/pressure/mimetic/hybsys_global.c \ opm/core/pressure/mimetic/hybsys.c \ +opm/core/simulator/SimulatorTwophase.cpp \ opm/core/transport/spu_explicit.c \ opm/core/transport/spu_implicit.c \ opm/core/transport/transport_source.c \ @@ -158,6 +159,7 @@ opm/core/grid/cpgpreprocess/geometry.h \ opm/core/grid/cpgpreprocess/facetopology.h \ opm/core/grid/cpgpreprocess/grdecl.h \ opm/core/utility/Average.hpp \ +opm/core/utility/ColumnExtract.hpp \ opm/core/utility/ErrorMacros.hpp \ opm/core/utility/Factory.hpp \ opm/core/utility/MonotCubicInterpolator.hpp \ @@ -193,17 +195,13 @@ opm/core/linalg/LinearSolverFactory.hpp \ opm/core/grid.h \ opm/core/well.h \ opm/core/newwells.h \ -opm/core/WellsGroup.hpp \ -opm/core/WellCollection.hpp \ -opm/core/InjectionSpecification.hpp \ -opm/core/ProductionSpecification.hpp \ -opm/core/BlackoilState.hpp \ -opm/core/ColumnExtract.hpp \ +opm/core/wells/WellsManager.hpp \ +opm/core/wells/WellsGroup.hpp \ +opm/core/wells/WellCollection.hpp \ +opm/core/wells/InjectionSpecification.hpp \ +opm/core/wells/ProductionSpecification.hpp \ opm/core/GridAdapter.hpp \ opm/core/GridManager.hpp \ -opm/core/TwophaseState.hpp \ -opm/core/WellsManager.hpp \ -opm/core/WellState.hpp \ opm/core/pressure/fsh.h \ opm/core/pressure/HybridPressureSolver.hpp \ opm/core/pressure/IncompTpfa.hpp \ @@ -230,6 +228,10 @@ opm/core/pressure/fsh_common_impl.h \ opm/core/pressure/mimetic/hybsys_global.h \ opm/core/pressure/mimetic/hybsys.h \ opm/core/pressure/mimetic/mimetic.h \ +opm/core/simulator/BlackoilState.hpp \ +opm/core/simulator/SimulatorTwophase.hpp \ +opm/core/simulator/TwophaseState.hpp \ +opm/core/simulator/WellState.hpp \ opm/core/transport/ImplicitTransport.hpp \ opm/core/transport/CSRMatrixUmfpackSolver.hpp \ opm/core/transport/CSRMatrixBlockAssembler.hpp \ @@ -271,3 +273,17 @@ opm/core/linalg/LinearSolverIstl.cpp nobase_include_HEADERS += \ opm/core/linalg/LinearSolverIstl.hpp endif + + +if BUILD_AGMG +libopmcore_la_SOURCES += \ +$(AGMG_SRCDIR)/dagmg.f90 \ +$(AGMG_SRCDIR)/dagmg_mumps.f90 \ +opm/core/linalg/LinearSolverAGMG.cpp + +nobase_include_HEADERS += \ +opm/core/linalg/LinearSolverAGMG.hpp + +libopmcore_la_LDFLAGS += \ +$(FCLIBS) +endif diff --git a/configure.ac b/configure.ac index 5b85b794..7f82ebcc 100644 --- a/configure.ac +++ b/configure.ac @@ -8,32 +8,33 @@ AM_INIT_AUTOMAKE([-Wall -Werror foreign subdir-objects]) m4_ifdef([AM_SILENT_RULES], [AM_SILENT_RULES([yes])]) -m4_ifdef([LT_INIT], [LT_INIT]) - AC_CONFIG_MACRO_DIR([m4]) AC_CONFIG_SRCDIR([opm/core/grid.h]) AC_CONFIG_HEADERS([config.h]) # Checks for programs. -AC_PROG_CXX AC_PROG_CC AM_PROG_CC_C_O -m4_ifdef([LT_INIT], [], [AC_PROG_LIBTOOL]) -AC_PROG_RANLIB -# AX_LAPACK requires a working F77 compiler or, rather, its runtime -# support libraries. -AC_PROG_F77 - -# F77 name mangling -AC_F77_WRAPPERS +m4_ifdef([LT_INIT], + [LT_INIT[]dnl + LT_LANG([C++])dnl + LT_LANG([Fortran 77])dnl + LT_LANG([Fortran])dnl + ],dnl + [AC_PROG_LIBTOOL[]dnl + AC_PROG_CXX[]dnl + AC_PROG_F77[]dnl + AC_PROG_FC[]dnl + ])[]dnl # Checks for libraries. # Bring in numerics support (standard library component) AC_SEARCH_LIBS([sqrt], [m]) -AX_LAPACK +OPM_LAPACK + AX_BOOST_BASE([1.37]) AX_BOOST_SYSTEM AX_BOOST_DATE_TIME @@ -41,6 +42,7 @@ AX_BOOST_FILESYSTEM AX_BOOST_UNIT_TEST_FRAMEWORK AX_DUNE_ISTL +OPM_AGMG # Checks for header files. AC_CHECK_HEADERS([float.h limits.h stddef.h stdlib.h string.h]) diff --git a/examples/Makefile.am b/examples/Makefile.am index a60d05ea..a6c29e80 100644 --- a/examples/Makefile.am +++ b/examples/Makefile.am @@ -1,31 +1,40 @@ +# Build-time flags needed to form example programs AM_CPPFLAGS = \ -I$(top_srcdir) \ $(BOOST_CPPFLAGS) +# All targets link to the library LDADD = $(top_builddir)/libopmcore.la +# ---------------------------------------------------------------------- +# Declare products (i.e., the example programs). +# +# Please keep the list sorted. + noinst_PROGRAMS = \ refine_wells \ scaneclipsedeck \ +sim_2p_incomp_reorder \ sim_wateroil \ wells_example -refine_wells_SOURCES = refine_wells.cpp +# ---------------------------------------------------------------------- +# Product constituents. Must be specified for every product that's +# built from more than a single ".c" file and/or that link to anything +# more than the OPM-Core library. +# +# Please maintain sort order from "noinst_PROGRAMS". -sim_wateroil_SOURCES = sim_wateroil.cpp -sim_wateroil_LDADD = \ -$(LDADD) $(LIBS) \ -$(BOOST_SYSTEM_LIB) $(BOOST_FILESYSTEM_LIB) \ -$(LAPACK_LIBS) $(LIBS) $(LIBS) +refine_wells_SOURCES = refine_wells.cpp +sim_2p_incomp_reorder_SOURCES = sim_2p_incomp_reorder.cpp +sim_wateroil_SOURCES = sim_wateroil.cpp +wells_example_SOURCES = wells_example.cpp -wells_example_SOURCES = wells_example.cpp +# ---------------------------------------------------------------------- +# Optional examples, or examples that use optional add-on components. if UMFPACK noinst_PROGRAMS += spu_2p -spu_2p_SOURCES = spu_2p.cpp -spu_2p_LDADD = \ -$(LDADD) $(LIBS) \ -$(BOOST_SYSTEM_LIB) $(BOOST_FILESYSTEM_LIB) \ -$(LAPACK_LIBS) $(LIBS) $(LIBS) +spu_2p_SOURCES = spu_2p.cpp endif diff --git a/examples/sim_2p_incomp_reorder.cpp b/examples/sim_2p_incomp_reorder.cpp new file mode 100644 index 00000000..0b1ed7fc --- /dev/null +++ b/examples/sim_2p_incomp_reorder.cpp @@ -0,0 +1,220 @@ +/* + 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 . +*/ + + +#if HAVE_CONFIG_H +#include "config.h" +#endif // HAVE_CONFIG_H + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + + + + +// ----------------- Main program ----------------- +int +main(int argc, char** argv) +{ + using namespace Opm; + + std::cout << "\n================ Test program for incompressible two-phase flow ===============\n\n"; + Opm::parameter::ParameterGroup param(argc, argv, false); + std::cout << "--------------- Reading parameters ---------------" << std::endl; + + // If we have a "deck_filename", grid and props will be read from that. + bool use_deck = param.has("deck_filename"); + boost::scoped_ptr grid; + boost::scoped_ptr props; + boost::scoped_ptr wells; + boost::scoped_ptr rock_comp; + Opm::SimulatorTimer simtimer; + Opm::TwophaseState state; + // bool check_well_controls = false; + // int max_well_control_iterations = 0; + double gravity[3] = { 0.0 }; + if (use_deck) { + std::string deck_filename = param.get("deck_filename"); + Opm::EclipseGridParser deck(deck_filename); + // Grid init + grid.reset(new Opm::GridManager(deck)); + // Rock and fluid init + const int* gc = grid->c_grid()->global_cell; + std::vector global_cell(gc, gc + grid->c_grid()->number_of_cells); + props.reset(new Opm::IncompPropertiesFromDeck(deck, global_cell)); + // Wells init. + wells.reset(new Opm::WellsManager(deck, *grid->c_grid(), props->permeability())); + // check_well_controls = param.getDefault("check_well_controls", false); + // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10); + // Timer init. + if (deck.hasField("TSTEP")) { + simtimer.init(deck); + } else { + simtimer.init(param); + } + // Rock compressibility. + rock_comp.reset(new Opm::RockCompressibility(deck)); + // Gravity. + gravity[2] = deck.hasField("NOGRAV") ? 0.0 : Opm::unit::gravity; + // Init state variables (saturation and pressure). + if (param.has("init_saturation")) { + initStateBasic(*grid->c_grid(), *props, param, gravity[2], state); + } else { + initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state); + } + } else { + // Grid init. + const int nx = param.getDefault("nx", 100); + const int ny = param.getDefault("ny", 100); + const int nz = param.getDefault("nz", 1); + const double dx = param.getDefault("dx", 1.0); + const double dy = param.getDefault("dy", 1.0); + const double dz = param.getDefault("dz", 1.0); + grid.reset(new Opm::GridManager(nx, ny, nz, dx, dy, dz)); + // Rock and fluid init. + props.reset(new Opm::IncompPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells)); + // Wells init. + wells.reset(new Opm::WellsManager()); + // Timer init. + simtimer.init(param); + // Rock compressibility. + rock_comp.reset(new Opm::RockCompressibility(param)); + // Gravity. + gravity[2] = param.getDefault("gravity", 0.0); + // Init state variables (saturation and pressure). + initStateBasic(*grid->c_grid(), *props, param, gravity[2], state); + } + + // Warn if gravity but no density difference. + bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0); + if (use_gravity) { + if (props->density()[0] == props->density()[1]) { + std::cout << "**** Warning: nonzero gravity, but zero density difference." << std::endl; + } + } + + // Source-related variables init. + int num_cells = grid->c_grid()->number_of_cells; + std::vector totmob; + std::vector omega; // Will remain empty if no gravity. + std::vector rc; // Will remain empty if no rock compressibility. + + // Extra rock init. + std::vector porevol; + if (rock_comp->isActive()) { + computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol); + } else { + computePorevolume(*grid->c_grid(), props->porosity(), porevol); + } + double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0); + + // Initialising src + std::vector src(num_cells, 0.0); + if (wells->c_wells()) { + // Do nothing, wells will be the driving force, not source terms. + // Opm::wellsToSrc(*wells->c_wells(), num_cells, src); + } else { + const double default_injection = use_gravity ? 0.0 : 0.1; + const double flow_per_sec = param.getDefault("injected_porevolumes_per_day", default_injection) + *tot_porevol_init/Opm::unit::day; + src[0] = flow_per_sec; + src[num_cells - 1] = -flow_per_sec; + } + + // Boundary conditions. + Opm::FlowBCManager bcs; + if (param.getDefault("use_pside", false)) { + int pside = param.get("pside"); + double pside_pressure = param.get("pside_pressure"); + bcs.pressureSide(*grid->c_grid(), Opm::FlowBCManager::Side(pside), pside_pressure); + } + + // Linear solver. + Opm::LinearSolverFactory linsolver(param); + + const double *grav = use_gravity ? &gravity[0] : 0; + + + Opm::SimulatorTwophase simulator(param, + *grid->c_grid(), + *props, + rock_comp->isActive() ? rock_comp.get() : 0, + wells->c_wells(), + src, + bcs.c_bcs(), + linsolver, + grav); + + // Warn if any parameters are unused. + if (param.anyUnused()) { + std::cout << "-------------------- Unused parameters: --------------------\n"; + param.displayUsage(); + std::cout << "----------------------------------------------------------------" << std::endl; + } + + // Write parameters used for later reference. + // if (output) { + // param.writeParam(output_dir + "/spu_2p.param"); + // } + + WellState well_state; + well_state.init(wells->c_wells(), state); + + simulator.run(simtimer, state, well_state); +} diff --git a/examples/sim_wateroil.cpp b/examples/sim_wateroil.cpp index 6c26c0fa..e2421bb7 100644 --- a/examples/sim_wateroil.cpp +++ b/examples/sim_wateroil.cpp @@ -27,7 +27,7 @@ #include #include #include -#include +#include #include #include #include @@ -44,9 +44,9 @@ #include -#include -#include -#include +#include +#include +#include #include #include diff --git a/examples/spu_2p.cpp b/examples/spu_2p.cpp index 07a7021b..702ec3b8 100644 --- a/examples/spu_2p.cpp +++ b/examples/spu_2p.cpp @@ -43,7 +43,7 @@ #include #include #include -#include +#include #include #include #include @@ -69,8 +69,8 @@ #include #include -#include -#include +#include +#include #include #include diff --git a/examples/wells_example.cpp b/examples/wells_example.cpp index 18c051aa..a8bdaae9 100644 --- a/examples/wells_example.cpp +++ b/examples/wells_example.cpp @@ -5,14 +5,14 @@ #include "opm/core/utility/initState.hpp" #include "opm/core/utility/SimulatorTimer.hpp" -#include +#include #include #include #include #include #include #include -#include +#include #include #include #include diff --git a/m4/agmg.m4 b/m4/agmg.m4 new file mode 100644 index 00000000..9e09157b --- /dev/null +++ b/m4/agmg.m4 @@ -0,0 +1,26 @@ +AC_DEFUN([OPM_AGMG],dnl +[AC_ARG_WITH([agmg],dnl + [AS_HELP_STRING([--with-agmg=SRCDIR],dnl + [Include DOUBLE PRECISION version Notay's of AGMG Algebraic + Multigrid solver from specified source location SRCDIR. Note: this + option requires a complete, working Fortran 90 environment.])], + [],dnl + [with_agmg=no])[]dnl + + AS_IF([test x"$with_agmg" != x"no" -a -d "$with_agmg"],dnl + [AS_IF([test -f "$with_agmg/dagmg.f90"],dnl + [AC_SUBST([AGMG_SRCDIR], [$with_agmg])[]dnl + AC_DEFINE([HAVE_AGMG], [1],dnl + [Define to `1' if Notay's AGMG solver is included])[]dnl + build_agmg="yes"],dnl + [build_agmg="no"])],dnl + [build_agmg="no"])[]dnl + + AS_IF([test x"$build_agmg" = x"yes"],dnl + [AC_PROG_FC_C_O[]dnl + AC_FC_WRAPPERS[]dnl + AC_FC_LIBRARY_LDFLAGS[]dnl + ], [:])[]dnl + + AM_CONDITIONAL([BUILD_AGMG], [test x"$build_agmg" = x"yes"]) +])[]dnl diff --git a/m4/opm_lapack.m4 b/m4/opm_lapack.m4 new file mode 100644 index 00000000..b5fa5ca5 --- /dev/null +++ b/m4/opm_lapack.m4 @@ -0,0 +1,4 @@ +AC_DEFUN([OPM_LAPACK], +[AC_REQUIRE([AC_F77_WRAPPERS])dnl + AC_REQUIRE([AX_LAPACK])dnl +])[]dnl diff --git a/opm/core/doxygen_main.hpp b/opm/core/doxygen_main.hpp new file mode 100644 index 00000000..e33af59a --- /dev/null +++ b/opm/core/doxygen_main.hpp @@ -0,0 +1,35 @@ +/* + 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 . +*/ + +#ifndef OPM_DOXYGEN_MAIN_HEADER_INCLUDED +#define OPM_DOXYGEN_MAIN_HEADER_INCLUDED + + +/** \mainpage Documentation for the opm-core library. + +The following are the main library features: +

Grid handling

+

Well handling

+

Pressure solvers

+

Transport solvers

+

Various utilities

+ +*/ + +#endif // OPM_DOXYGEN_MAIN_HEADER_INCLUDED diff --git a/opm/core/eclipse/EclipseGridParser.cpp b/opm/core/eclipse/EclipseGridParser.cpp index 29d8bb05..a91629c2 100644 --- a/opm/core/eclipse/EclipseGridParser.cpp +++ b/opm/core/eclipse/EclipseGridParser.cpp @@ -41,6 +41,7 @@ #include #include #include +#include #include #include #include @@ -97,6 +98,7 @@ namespace EclipseKeywords string("SGOF"), string("SWOF"), string("ROCK"), string("ROCKTAB"), string("WELSPECS"), string("COMPDAT"), string("WCONINJE"), string("WCONPROD"), string("WELTARG"), + string("WELOPEN"), string("EQUIL"), string("PVCDO"), string("TSTEP"), string("PLYVISC"), string("PLYROCK"), string("PLYADS"), string("PLYMAX"), string("TLMIXPAR"), string("WPOLYMER"), @@ -147,6 +149,7 @@ namespace { enum FieldType { Integer, FloatingPoint, + Timestepping, SpecialField, IgnoreWithData, IgnoreNoData, @@ -161,6 +164,8 @@ namespace { return Integer; } else if (count(floating_fields, floating_fields + num_floating_fields, keyword)) { return FloatingPoint; + } else if (keyword == "TSTEP" || keyword == "DATES") { + return Timestepping; } else if (count(special_fields, special_fields + num_special_fields, keyword)) { return SpecialField; } else if (count(ignore_with_data, ignore_with_data + num_ignore_with_data, keyword)) { @@ -276,6 +281,10 @@ namespace { //--------------------------------------------------------------------------- EclipseGridParser::EclipseGridParser() //--------------------------------------------------------------------------- + : current_reading_mode_(Regular), + start_date_(boost::date_time::not_a_date_time), + current_time_days_(0.0), + current_epoch_(0) { } @@ -284,6 +293,10 @@ EclipseGridParser::EclipseGridParser() //--------------------------------------------------------------------------- EclipseGridParser::EclipseGridParser(const string& filename, bool convert_to_SI) //--------------------------------------------------------------------------- + : current_reading_mode_(Regular), + start_date_(boost::date_time::not_a_date_time), + current_time_days_(0.0), + current_epoch_(0) { // Store directory of filename boost::filesystem::path p(filename); @@ -304,13 +317,39 @@ void EclipseGridParser::read(istream& is, bool convert_to_SI) { integer_field_map_.clear(); floating_field_map_.clear(); - special_field_map_.clear(); + special_field_by_epoch_.clear(); + special_field_by_epoch_.push_back(SpecialMap()); readImpl(is); + + current_epoch_ = 0; + computeUnits(); if (convert_to_SI) { convertToSI(); } + +#define VERBOSE_LIST_FIELDS 1 +#if VERBOSE_LIST_FIELDS + std::cout << "\nInteger fields:\n"; + for (std::map >::iterator + i = integer_field_map_.begin(); i != integer_field_map_.end(); ++i) + std::cout << '\t' << i->first << '\n'; + + std::cout << "\nFloat fields:\n"; + for (std::map >::iterator + i = floating_field_map_.begin(); i != floating_field_map_.end(); ++i) + std::cout << '\t' << i->first << '\n'; + + std::cout << "\nSpecial fields:\n"; + for (int epoch = 0; epoch < numberOfEpochs(); ++epoch) { + std::cout << "Epoch " << epoch << '\n'; + const SpecialMap& sm = special_field_by_epoch_[epoch]; + for (SpecialMap::const_iterator i = sm.begin(); i != sm.end(); ++i) { + std::cout << '\t' << i->first << '\n'; + } + } +#endif } //--------------------------------------------------------------------------- @@ -329,7 +368,6 @@ void EclipseGridParser::readImpl(istream& is) // though (of course retaining the basic guarantee). map >& intmap = integer_field_map_; map >& floatmap = floating_field_map_; - map >& specialmap = special_field_map_; // Actually read the data std::string keyword; @@ -341,17 +379,77 @@ void EclipseGridParser::readImpl(istream& is) cout << "Keyword found: " << keyword << endl; //#endif FieldType type = classifyKeyword(keyword); + // std::cout << "Classification: " << type << std::endl; switch (type) { - case Integer: + case Integer: { readVectorData(is, intmap[keyword]); break; - case FloatingPoint: + } + case FloatingPoint: { readVectorData(is, floatmap[keyword]); break; + } + case Timestepping: { + SpecialMap& sm = special_field_by_epoch_[current_epoch_]; + if (start_date_.is_not_a_date()) { + // Set it to START date, or default if no START. + // This will only ever happen in the first epoch, + // upon first encountering a timestepping keyword. + SpecialMap::const_iterator it = sm.find("START"); + if (hasField("START")) { + start_date_ = getSTART().date; + } else { + start_date_ = boost::gregorian::date(1983, 1, 1); + } + } + if (current_reading_mode_ == Regular) { + current_reading_mode_ = Timesteps; + } + // Get current epoch's TSTEP, if it exists, create new if not. + SpecialMap::iterator it = sm.find("TSTEP"); + TSTEP* tstep = 0; + if (it != sm.end()) { + tstep = dynamic_cast(it->second.get()); + } else { + SpecialFieldPtr sb_ptr(new TSTEP()); + tstep = dynamic_cast(sb_ptr.get()); + sm["TSTEP"] = sb_ptr; + } + ASSERT(tstep != 0); + // Append new steps to current TSTEP object + if (keyword == "TSTEP") { + const int num_steps_old = tstep->tstep_.size(); + tstep->read(is); // This will append to the TSTEP object. + const double added_days + = std::accumulate(tstep->tstep_.begin() + num_steps_old, tstep->tstep_.end(), 0.0); + current_time_days_ += added_days; + } else if (keyword == "DATES") { + DATES dates; + dates.read(is); + for (std::size_t dix = 0; dix < dates.dates.size(); ++dix) { + boost::gregorian::date_duration since_start = dates.dates[dix] - start_date_; + double step = double(since_start.days()) - current_time_days_; + tstep->tstep_.push_back(step); + current_time_days_ = double(since_start.days()); + } + } else { + THROW("Keyword " << keyword << " cannot be handled here."); + } + break; + } case SpecialField: { - std::tr1::shared_ptr sb_ptr = createSpecialField(is, keyword); + if (current_reading_mode_ == Timesteps) { + // We have been reading timesteps, but have + // now encountered something else. + // That means we are in a new epoch. + current_reading_mode_ = Regular; + special_field_by_epoch_.push_back(SpecialMap()); + ++current_epoch_; + ASSERT(int(special_field_by_epoch_.size()) == current_epoch_ + 1); + } + SpecialFieldPtr sb_ptr = createSpecialField(is, keyword); if (sb_ptr) { - specialmap[keyword] = sb_ptr; + special_field_by_epoch_[current_epoch_][keyword] = sb_ptr; } else { THROW("Could not create field " << keyword); } @@ -398,24 +496,6 @@ void EclipseGridParser::readImpl(istream& is) is >> ignoreLine; } } - -#define VERBOSE_LIST_FIELDS 0 -#if VERBOSE_LIST_FIELDS - std::cout << "\nInteger fields:\n"; - for (std::map >::iterator - i = intmap.begin(); i != intmap.end(); ++i) - std::cout << '\t' << i->first << '\n'; - - std::cout << "\nFloat fields:\n"; - for (std::map >::iterator - i = floatmap.begin(); i != floatmap.end(); ++i) - std::cout << '\t' << i->first << '\n'; - - std::cout << "\nSpecial fields:\n"; - for (std::map >::iterator - i = specialmap.begin(); i != specialmap.end(); ++i) - std::cout << '\t' << i->first << '\n'; -#endif } @@ -425,9 +505,12 @@ void EclipseGridParser::convertToSI() //--------------------------------------------------------------------------- { // Convert all special fields. - typedef std::map >::iterator SpecialIt; - for (SpecialIt i = special_field_map_.begin(); i != special_field_map_.end(); ++i) { - i->second->convertToSI(units_); + typedef SpecialMap::iterator SpecialIt; + for (int epoch = 0; epoch < numberOfEpochs(); ++epoch) { + SpecialMap& sm = special_field_by_epoch_[epoch]; + for (SpecialIt i = sm.begin(); i != sm.end(); ++i) { + i->second->convertToSI(units_); + } } // Convert all floating point fields. @@ -479,7 +562,7 @@ bool EclipseGridParser::hasField(const string& keyword) const { string ukey = upcase(keyword); return integer_field_map_.count(ukey) || floating_field_map_.count(ukey) || - special_field_map_.count(ukey) || ignored_fields_.count(ukey); + special_field_by_epoch_[current_epoch_].count(ukey) || ignored_fields_.count(ukey); } @@ -505,7 +588,7 @@ vector EclipseGridParser::fieldNames() const vector names; names.reserve(integer_field_map_.size() + floating_field_map_.size() + - special_field_map_.size() + + special_field_by_epoch_[current_epoch_].size() + ignored_fields_.size()); { map >::const_iterator it = integer_field_map_.begin(); @@ -520,8 +603,8 @@ vector EclipseGridParser::fieldNames() const } } { - map >::const_iterator it = special_field_map_.begin(); - for (; it != special_field_map_.end(); ++it) { + SpecialMap::const_iterator it = special_field_by_epoch_[current_epoch_].begin(); + for (; it != special_field_by_epoch_[current_epoch_].end(); ++it) { names.push_back(it->first); } } @@ -534,6 +617,24 @@ vector EclipseGridParser::fieldNames() const return names; } + +//--------------------------------------------------------------------------- +int EclipseGridParser::numberOfEpochs() const +//--------------------------------------------------------------------------- +{ + return special_field_by_epoch_.size(); +} + + +//--------------------------------------------------------------------------- +void EclipseGridParser::setCurrentEpoch(int epoch) +//--------------------------------------------------------------------------- +{ + ASSERT(epoch >= 0 && epoch < numberOfEpochs()); + current_epoch_ = epoch; +} + + //--------------------------------------------------------------------------- const std::vector& EclipseGridParser::getIntegerValue(const std::string& keyword) const //--------------------------------------------------------------------------- @@ -574,8 +675,8 @@ const std::vector& EclipseGridParser::getFloatingPointValue(const std::s const std::tr1::shared_ptr EclipseGridParser::getSpecialValue(const std::string& keyword) const //--------------------------------------------------------------------------- { - map >::const_iterator it = special_field_map_.find(keyword); - if (it == special_field_map_.end()) { + SpecialMap::const_iterator it = special_field_by_epoch_[current_epoch_].find(keyword); + if (it == special_field_by_epoch_[current_epoch_].end()) { THROW("No such field: " << keyword); } else { return it->second; @@ -617,7 +718,7 @@ void EclipseGridParser::setSpecialField(const std::string& keyword, std::tr1::shared_ptr field) //--------------------------------------------------------------------------- { - special_field_map_[keyword] = field; + special_field_by_epoch_[current_epoch_][keyword] = field; } //--------------------------------------------------------------------------- diff --git a/opm/core/eclipse/EclipseGridParser.hpp b/opm/core/eclipse/EclipseGridParser.hpp index 5662ae87..659b3484 100644 --- a/opm/core/eclipse/EclipseGridParser.hpp +++ b/opm/core/eclipse/EclipseGridParser.hpp @@ -94,6 +94,14 @@ public: /// The keywords/fields found in the file. std::vector fieldNames() const; + /// Returns the number of distinct epochs found in the deck SCHEDULE. + int numberOfEpochs() const; + + /// Sets the current epoch. + /// Valid arguments are in [0, ..., numberOfEpochs() - 1]. + /// After reading, current epoch always starts at 0. + void setCurrentEpoch(int epoch); + /// Returns a reference to a vector containing the values /// corresponding to the given integer keyword. const std::vector& getIntegerValue(const std::string& keyword) const; @@ -102,10 +110,12 @@ public: /// corresponding to the given floating-point keyword. const std::vector& getFloatingPointValue(const std::string& keyword) const; + typedef std::tr1::shared_ptr SpecialFieldPtr; + /// Returns a reference to a vector containing pointers to the values /// corresponding to the given keyword when the values are not only integers /// or floats. - const std::tr1::shared_ptr getSpecialValue(const std::string& keyword) const; + const SpecialFieldPtr getSpecialValue(const std::string& keyword) const; // This macro implements support for a special field keyword. It requires that a subclass // of SpecialBase exists, that has the same name as the keyword. @@ -140,6 +150,7 @@ public: SPECIAL_FIELD(WCONINJE); SPECIAL_FIELD(WCONPROD); SPECIAL_FIELD(WELTARG); + SPECIAL_FIELD(WELOPEN); SPECIAL_FIELD(EQUIL); SPECIAL_FIELD(PVCDO); SPECIAL_FIELD(TSTEP); @@ -170,7 +181,7 @@ public: void setFloatingPointField(const std::string& keyword, const std::vector& field); /// Sets a special field to have a particular value. - void setSpecialField(const std::string& keyword, std::tr1::shared_ptr field); + void setSpecialField(const std::string& keyword, SpecialFieldPtr field); /// Compute the units used by the deck, depending on the presence /// of keywords such as METRIC, FIELD etc. It is an error to call @@ -181,17 +192,25 @@ public: const EclipseUnits& units() const; private: - std::tr1::shared_ptr createSpecialField(std::istream& is, const std::string& fieldname); + SpecialFieldPtr createSpecialField(std::istream& is, const std::string& fieldname); void readImpl(std::istream& is); std::string directory_; std::map > integer_field_map_; std::map > floating_field_map_; - std::map > special_field_map_; + // std::map special_field_map_; std::set ignored_fields_; - std::tr1::shared_ptr empty_special_field_; EclipseUnits units_; + + // For SCHEDULE handling. + enum ReadingMode { Regular, Timesteps }; + ReadingMode current_reading_mode_; + boost::gregorian::date start_date_; + double current_time_days_; + int current_epoch_; + typedef std::map SpecialMap; + std::vector special_field_by_epoch_; }; diff --git a/opm/core/eclipse/SpecialEclipseFields.hpp b/opm/core/eclipse/SpecialEclipseFields.hpp index d7e3548f..7e2fce9a 100644 --- a/opm/core/eclipse/SpecialEclipseFields.hpp +++ b/opm/core/eclipse/SpecialEclipseFields.hpp @@ -1658,6 +1658,67 @@ struct WELTARG : public SpecialBase }; + +struct WelopenLine +{ + std::string well_; // Well name or well name root + std::string openshutflag_; // What to do with the well. +}; + +/// Class for keyword WELOPEN +struct WELOPEN : public SpecialBase +{ + std::vector welopen; + + WELOPEN() + { + } + + virtual ~WELOPEN() + {} + + virtual std::string name() const {return std::string("WELOPEN");} + + virtual void read(std::istream& is) + { + while(is) { + std::string wellname = readString(is); + if (wellname[0] == '/') { + is >> ignoreLine; + break; + } + while (wellname.find("--") == 0) { + // This line is a comment + is >> ignoreLine; + wellname = readString(is); + } + WelopenLine welopen_line; + welopen_line.well_ = wellname; + welopen_line.openshutflag_ = readString(is); + ignoreSlashLine(is); + welopen.push_back(welopen_line); + } + } + + virtual void write(std::ostream& os) const + { + os << name() << std::endl; + for (int i=0; i<(int)welopen.size(); ++i) { + os << welopen[i].well_ << " " + << welopen[i].openshutflag_ + << std::endl; + } + os << std::endl; + } + + virtual void convertToSI(const EclipseUnits& /*units*/) + { + } +}; + + + + /// Class holding a data line of keyword EQUIL struct EquilLine { diff --git a/opm/core/linalg/LinearSolverAGMG.cpp b/opm/core/linalg/LinearSolverAGMG.cpp new file mode 100644 index 00000000..e028e2c5 --- /dev/null +++ b/opm/core/linalg/LinearSolverAGMG.cpp @@ -0,0 +1,124 @@ +/* + 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 . +*/ + + +#if HAVE_CONFIG_H +#include +#endif + +#include +#include + +#include +#include +#include +#include +#include + +#if HAVE_AGMG +// Manual prototype of main gateway routine to DOUBLE PRECISION, +// serial version of the AGMG software package. +// +// Note that both the matrix entries and column indices are writable. +// The solver may permute the matrix entries within each row during +// the setup phase. +#define DAGMG_ FC_FUNC(dagmg, DAGMG) + +extern "C" +void +DAGMG_(const int* N , // System size + double* sa , // Non-zero elements + int* ja , // Column indices + const int* ia , // Row pointers + const double* f , // Right-hand side + double* x , // Solution + int* ijob , // Main control parameter + int* iprint, // Message output unit + int* nrest , // Number of GCR restarts + int* iter , // Maximum (and actual) number of iterations + const double* tol ); // Residual reduction tolerance + +#endif // HAVE_AGMG + +namespace Opm +{ + LinearSolverAGMG::LinearSolverAGMG(const int max_it, + const double rtol , + const bool is_spd) + : max_it_(max_it), + rtol_ (rtol) , + is_spd_(is_spd) + { +#if !HAVE_AGMG + THROW("AGMG support is not enabled in this library"); +#endif // HAVE_AGMG + } + + LinearSolverAGMG::~LinearSolverAGMG() {} + + LinearSolverInterface::LinearSolverReport + LinearSolverAGMG::solve(const int size , + const int nonzeros, + const int* ia , + const int* ja , + const double* sa , + const double* rhs , + double* solution) const + { + const std::vector::size_type nnz = ia[size]; + + assert (nnz == std::vector::size_type(nonzeros)); + +#if defined(NDEBUG) + // Suppress warning about unused parameter. + static_cast(nonzeros); +#endif + + std::vector a(sa, sa + nnz); + + // Account for 1-based indexing. + std::vector i(ia, ia + std::vector::size_type(size + 1)); + std::transform(i.begin(), i.end(), i.begin(), + std::bind2nd(std::plus(), 1)); + + std::vector j(ja, ja + nnz); + std::transform(j.begin(), j.end(), j.begin(), + std::bind2nd(std::plus(), 1)); + + LinearSolverInterface::LinearSolverReport rpt = {}; + rpt.iterations = max_it_; + + int ijob = 0; // Setup + solution + cleanup, x0==0. + + int nrest; + if (is_spd_) { + nrest = 1; // Use CG algorithm + } + else { + nrest = 10; // Suggested default number of GCR restarts. + } + + int iprint = 0; // Suppress most output + DAGMG_(& size, & a[0], & j[0], & i[0], rhs, solution, + & ijob, & iprint, & nrest, & rpt.iterations, & rtol_); + + rpt.converged = rpt.iterations <= max_it_; + return rpt; + } +} // namespace Opm diff --git a/opm/core/linalg/LinearSolverAGMG.hpp b/opm/core/linalg/LinearSolverAGMG.hpp new file mode 100644 index 00000000..6056dc90 --- /dev/null +++ b/opm/core/linalg/LinearSolverAGMG.hpp @@ -0,0 +1,96 @@ +/* + 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 . +*/ + +#ifndef OPM_LINEARSOLVERAGMG_HEADER_INCLUDED +#define OPM_LINEARSOLVERAGMG_HEADER_INCLUDED + +/** + * \file + * Gateway to Notay's AGMG package implementing an aggregation-based + * algebraic multigrid method. + * + * References: + * * Y. Notay, + * An aggregation-based algebraic multigrid method, + * Electronic Transactions on Numerical Analysis, vol 37, + * pp. 123-146, 2010. + * + * * A. Napov and Y. Notay, + * An algebraic multigrid method with guaranteed convergence rate, + * Report GANMN 10-03, Université Libre de Bruxelles, Brussels, + * Belgium, 2010 (Revised 2011). + * + * * Y. Notay, + * Aggregation-based algebraic multigrid for convection-diffusion + * equations, Report GANMN 11-01, Université Libre de Bruxelles, + * Brussels, Belgium, 2011. + */ + +#include + +namespace Opm +{ + /// Concrete class encapsulating Notay's AGMG package + class LinearSolverAGMG : public LinearSolverInterface + { + public: + /** + * Constructor. + * \param[in] max_it Maximum number of solver iterations. + * \param[in] rtol Residual reduction tolerance. + * \param[in] is_spd Whether or not the matrix is SPD. SPD + * systems are solved using CG while other + * systems are solved using GCR. + */ + LinearSolverAGMG(const int max_it = 100 , + const double rtol = 1.0e-6, + const bool is_spd = false); + + /// Destructor. + virtual ~LinearSolverAGMG(); + + using LinearSolverInterface::solve; + + /// Solve a linear system, with a matrix given in compressed + /// sparse row format. + /// \param[in] size Number of rows (and colums). + /// \param[in] nonzeros Number of (structural) non-zeros. + /// \param[in] ia Row pointers. + /// \param[in] ja Column indices. + /// \param[in] sa (structurally) non-zero elements. + /// \param[in] rhs System right-hand side. + /// \param[inout] solution System solution. + /// \return Solver meta-data concerning most recent system solve. + virtual LinearSolverInterface::LinearSolverReport + solve(const int size, const int nonzeros, + const int* ia, const int* ja, const double* sa, + const double* rhs, double* solution) const; + + private: + int max_it_; + double rtol_ ; + bool is_spd_; + }; + + +} // namespace Opm + + + +#endif // OPM_LINEARSOLVERAGMG_HEADER_INCLUDED diff --git a/opm/core/utility/newwells.c b/opm/core/newwells.c similarity index 100% rename from opm/core/utility/newwells.c rename to opm/core/newwells.c diff --git a/opm/core/pressure/CompressibleTpfa.cpp b/opm/core/pressure/CompressibleTpfa.cpp index 0ee8bf8b..83d3b804 100644 --- a/opm/core/pressure/CompressibleTpfa.cpp +++ b/opm/core/pressure/CompressibleTpfa.cpp @@ -28,8 +28,8 @@ #include #include #include -#include -#include +#include +#include #include #include diff --git a/opm/core/BlackoilState.hpp b/opm/core/simulator/BlackoilState.hpp similarity index 100% rename from opm/core/BlackoilState.hpp rename to opm/core/simulator/BlackoilState.hpp diff --git a/opm/core/simulator/SimulatorTwophase.cpp b/opm/core/simulator/SimulatorTwophase.cpp new file mode 100644 index 00000000..1a617a4b --- /dev/null +++ b/opm/core/simulator/SimulatorTwophase.cpp @@ -0,0 +1,461 @@ +/* + 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 . +*/ + +#if HAVE_CONFIG_H +#include "config.h" +#endif // HAVE_CONFIG_H + +#include +#include +#include + +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include + +#include +#include + + +namespace Opm +{ + + class SimulatorTwophase::Impl + { + public: + Impl(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const IncompPropertiesInterface& props, + const RockCompressibility* rock_comp, + const Wells* wells, + const std::vector& src, + const FlowBoundaryConditions* bcs, + const LinearSolverInterface& linsolver, + const double* gravity); + + void run(SimulatorTimer& timer, + TwophaseState& state, + WellState& well_state); + + private: + // Data. + + // Parameters for output. + bool output_; + std::string output_dir_; + int output_interval_; + // Parameters for pressure solver. + int nl_pressure_maxiter_; + double nl_pressure_tolerance_; + // Parameters for transport solver. + int nl_maxiter_; + double nl_tolerance_; + int num_transport_substeps_; + bool use_segregation_split_; + // Observed objects. + const UnstructuredGrid& grid_; + const IncompPropertiesInterface& props_; + const RockCompressibility* rock_comp_; + const Wells* wells_; + const std::vector& src_; + const FlowBoundaryConditions* bcs_; + const LinearSolverInterface& linsolver_; + const double* gravity_; + // Solvers + IncompTpfa psolver_; + TransportModelTwophase tsolver_; + // Needed by column-based gravity segregation solver. + std::vector< std::vector > columns_; + // Misc. data + std::vector allcells_; + }; + + + + + SimulatorTwophase::SimulatorTwophase(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const IncompPropertiesInterface& props, + const RockCompressibility* rock_comp, + const Wells* wells, + const std::vector& src, + const FlowBoundaryConditions* bcs, + const LinearSolverInterface& linsolver, + const double* gravity) + { + pimpl_.reset(new Impl(param, grid, props, rock_comp, wells, src, bcs, linsolver, gravity)); + } + + + + + + void SimulatorTwophase::run(SimulatorTimer& timer, + TwophaseState& state, + WellState& well_state) + { + pimpl_->run(timer, state, well_state); + } + + + + static void outputState(const UnstructuredGrid& grid, + const Opm::TwophaseState& state, + const int step, + const std::string& output_dir) + { + // Write data in VTK format. + std::ostringstream vtkfilename; + vtkfilename << output_dir << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu"; + std::ofstream vtkfile(vtkfilename.str().c_str()); + if (!vtkfile) { + THROW("Failed to open " << vtkfilename.str()); + } + Opm::DataMap dm; + dm["saturation"] = &state.saturation(); + dm["pressure"] = &state.pressure(); + std::vector cell_velocity; + Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity); + dm["velocity"] = &cell_velocity; + Opm::writeVtkData(grid, dm, vtkfile); + + // Write data (not grid) in Matlab format + for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) { + std::ostringstream fname; + fname << output_dir << "/" << it->first << "-" << std::setw(3) << std::setfill('0') << step << ".dat"; + std::ofstream file(fname.str().c_str()); + if (!file) { + THROW("Failed to open " << fname.str()); + } + const std::vector& d = *(it->second); + std::copy(d.begin(), d.end(), std::ostream_iterator(file, "\n")); + } + } + + + static void outputWaterCut(const Opm::Watercut& watercut, + const std::string& output_dir) + { + // Write water cut curve. + std::string fname = output_dir + "/watercut.txt"; + std::ofstream os(fname.c_str()); + if (!os) { + THROW("Failed to open " << fname); + } + watercut.write(os); + } + + + static void outputWellReport(const Opm::WellReport& wellreport, + const std::string& output_dir) + { + // Write well report. + std::string fname = output_dir + "/wellreport.txt"; + std::ofstream os(fname.c_str()); + if (!os) { + THROW("Failed to open " << fname); + } + wellreport.write(os); + } + + + + + + SimulatorTwophase::Impl::Impl(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const IncompPropertiesInterface& props, + const RockCompressibility* rock_comp, + const Wells* wells, + const std::vector& src, + const FlowBoundaryConditions* bcs, + const LinearSolverInterface& linsolver, + const double* gravity) + : grid_(grid), + props_(props), + rock_comp_(rock_comp), + wells_(wells), + src_(src), + bcs_(bcs), + linsolver_(linsolver), + gravity_(gravity), + psolver_(grid, props.permeability(), gravity, linsolver), + tsolver_(grid, props, 1e-9, 30) + { + // For output. + output_ = param.getDefault("output", true); + if (output_) { + output_dir_ = param.getDefault("output_dir", std::string("output")); + // Ensure that output dir exists + boost::filesystem::path fpath(output_dir_); + try { + create_directories(fpath); + } + catch (...) { + THROW("Creating directories failed: " << fpath); + } + output_interval_ = param.getDefault("output_interval", 1); + } + + // For pressure solver + nl_pressure_maxiter_ = param.getDefault("nl_pressure_maxiter", 10); + nl_pressure_tolerance_ = param.getDefault("nl_pressure_tolerance", 1.0); // Pascal + + // For transport solver. + nl_maxiter_ = param.getDefault("nl_maxiter", 30); + nl_tolerance_ = param.getDefault("nl_tolerance", 1e-9); + num_transport_substeps_ = param.getDefault("num_transport_substeps", 1); + use_segregation_split_ = param.getDefault("use_segregation_split", false); + if (gravity != 0 && use_segregation_split_){ + tsolver_.initGravity(gravity); + extractColumn(grid_, columns_); + } + + // Misc init. + const int num_cells = grid.number_of_cells; + allcells_.resize(num_cells); + for (int cell = 0; cell < num_cells; ++cell) { + allcells_[cell] = cell; + } + } + + + + + void SimulatorTwophase::Impl::run(SimulatorTimer& timer, + TwophaseState& state, + WellState& well_state) + { + std::vector totmob; + std::vector omega; // Will remain empty if no gravity. + std::vector rc; // Will remain empty if no rock compressibility. + std::vector transport_src; + + // Initialisation. + std::vector porevol; + if (rock_comp_ && rock_comp_->isActive()) { + computePorevolume(grid_, props_.porosity(), *rock_comp_, state.pressure(), porevol); + } else { + computePorevolume(grid_, props_.porosity(), porevol); + } + const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0); + + + // Main simulation loop. + Opm::time::StopWatch pressure_timer; + double ptime = 0.0; + Opm::time::StopWatch transport_timer; + double ttime = 0.0; + Opm::time::StopWatch total_timer; + total_timer.start(); + std::cout << "\n\n================ Starting main simulation loop ===============" << std::endl; + double init_satvol[2] = { 0.0 }; + double satvol[2] = { 0.0 }; + double injected[2] = { 0.0 }; + double produced[2] = { 0.0 }; + double tot_injected[2] = { 0.0 }; + double tot_produced[2] = { 0.0 }; + Opm::computeSaturatedVol(porevol, state.saturation(), init_satvol); + std::cout << "\nInitial saturations are " << init_satvol[0]/tot_porevol_init + << " " << init_satvol[1]/tot_porevol_init << std::endl; + Opm::Watercut watercut; + watercut.push(0.0, 0.0, 0.0); + Opm::WellReport wellreport; + std::vector fractional_flows; + std::vector well_resflows_phase; + int num_wells = 0; + if (wells_) { + num_wells = wells_->number_of_wells; + well_resflows_phase.resize((wells_->number_of_phases)*(wells_->number_of_wells), 0.0); + wellreport.push(props_, *wells_, state.saturation(), 0.0, well_state.bhp(), well_state.perfRates()); + } + const int num_cells = grid_.number_of_cells; + for (; !timer.done(); ++timer) { + // Report timestep and (optionally) write state to disk. + timer.report(std::cout); + if (output_ && (timer.currentStepNum() % output_interval_ == 0)) { + outputState(grid_, state, timer.currentStepNum(), output_dir_); + } + + // Solve pressure. + if (gravity_) { + computeTotalMobilityOmega(props_, allcells_, state.saturation(), totmob, omega); + } else { + computeTotalMobility(props_, allcells_, state.saturation(), totmob); + } + std::vector wdp; + if (wells_) { + Opm::computeWDP(*wells_, grid_, state.saturation(), props_.density(), + gravity_ ? gravity_[2] : 0.0, true, wdp); + } + do { + pressure_timer.start(); + if (rock_comp_ && rock_comp_->isActive()) { + rc.resize(num_cells); + std::vector initial_pressure = state.pressure(); + std::vector initial_porevolume(num_cells); + computePorevolume(grid_, props_.porosity(), *rock_comp_, initial_pressure, initial_porevolume); + std::vector pressure_increment(num_cells + num_wells); + std::vector prev_pressure(num_cells + num_wells); + for (int iter = 0; iter < nl_pressure_maxiter_; ++iter) { + for (int cell = 0; cell < num_cells; ++cell) { + rc[cell] = rock_comp_->rockComp(state.pressure()[cell]); + } + computePorevolume(grid_, props_.porosity(), *rock_comp_, state.pressure(), porevol); + std::copy(state.pressure().begin(), state.pressure().end(), prev_pressure.begin()); + std::copy(well_state.bhp().begin(), well_state.bhp().end(), prev_pressure.begin() + num_cells); + // prev_pressure = state.pressure(); + + // compute pressure increment + psolver_.solveIncrement(totmob, omega, src_, wdp, bcs_, porevol, rc, + prev_pressure, initial_porevolume, timer.currentStepLength(), + pressure_increment); + + double max_change = 0.0; + for (int cell = 0; cell < num_cells; ++cell) { + state.pressure()[cell] += pressure_increment[cell]; + max_change = std::max(max_change, std::fabs(pressure_increment[cell])); + } + for (int well = 0; well < num_wells; ++well) { + well_state.bhp()[well] += pressure_increment[num_cells + well]; + max_change = std::max(max_change, std::fabs(pressure_increment[num_cells + well])); + } + + std::cout << "Pressure iter " << iter << " max change = " << max_change << std::endl; + if (max_change < nl_pressure_tolerance_) { + break; + } + } + psolver_.computeFaceFlux(totmob, omega, src_, wdp, bcs_, state.pressure(), state.faceflux(), + well_state.bhp(), well_state.perfRates()); + } else { + psolver_.solve(totmob, omega, src_, wdp, bcs_, state.pressure(), state.faceflux(), + well_state.bhp(), well_state.perfRates()); + } + pressure_timer.stop(); + double pt = pressure_timer.secsSinceStart(); + std::cout << "Pressure solver took: " << pt << " seconds." << std::endl; + ptime += pt; + } while (false); + + // Process transport sources (to include bdy terms and well flows). + Opm::computeTransportSource(grid_, src_, state.faceflux(), 1.0, + wells_, well_state.perfRates(), transport_src); + + // Solve transport. + transport_timer.start(); + double stepsize = timer.currentStepLength(); + if (num_transport_substeps_ != 1) { + stepsize /= double(num_transport_substeps_); + std::cout << "Making " << num_transport_substeps_ << " transport substeps." << std::endl; + } + for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) { + tsolver_.solve(&state.faceflux()[0], &porevol[0], &transport_src[0], + stepsize, state.saturation()); + Opm::computeInjectedProduced(props_, state.saturation(), transport_src, stepsize, injected, produced); + if (use_segregation_split_) { + tsolver_.solveGravity(columns_, &porevol[0], stepsize, state.saturation()); + } + } + transport_timer.stop(); + double tt = transport_timer.secsSinceStart(); + std::cout << "Transport solver took: " << tt << " seconds." << std::endl; + ttime += tt; + + // Report volume balances. + Opm::computeSaturatedVol(porevol, state.saturation(), satvol); + tot_injected[0] += injected[0]; + tot_injected[1] += injected[1]; + tot_produced[0] += produced[0]; + tot_produced[1] += produced[1]; + std::cout.precision(5); + const int width = 18; + std::cout << "\nVolume balance report (all numbers relative to total pore volume).\n"; + std::cout << " Saturated volumes: " + << std::setw(width) << satvol[0]/tot_porevol_init + << std::setw(width) << satvol[1]/tot_porevol_init << std::endl; + std::cout << " Injected volumes: " + << std::setw(width) << injected[0]/tot_porevol_init + << std::setw(width) << injected[1]/tot_porevol_init << std::endl; + std::cout << " Produced volumes: " + << std::setw(width) << produced[0]/tot_porevol_init + << std::setw(width) << produced[1]/tot_porevol_init << std::endl; + std::cout << " Total inj volumes: " + << std::setw(width) << tot_injected[0]/tot_porevol_init + << std::setw(width) << tot_injected[1]/tot_porevol_init << std::endl; + std::cout << " Total prod volumes: " + << std::setw(width) << tot_produced[0]/tot_porevol_init + << std::setw(width) << tot_produced[1]/tot_porevol_init << std::endl; + std::cout << " In-place + prod - inj: " + << std::setw(width) << (satvol[0] + tot_produced[0] - tot_injected[0])/tot_porevol_init + << std::setw(width) << (satvol[1] + tot_produced[1] - tot_injected[1])/tot_porevol_init << std::endl; + std::cout << " Init - now - pr + inj: " + << std::setw(width) << (init_satvol[0] - satvol[0] - tot_produced[0] + tot_injected[0])/tot_porevol_init + << std::setw(width) << (init_satvol[1] - satvol[1] - tot_produced[1] + tot_injected[1])/tot_porevol_init + << std::endl; + std::cout.precision(8); + + watercut.push(timer.currentTime() + timer.currentStepLength(), + produced[0]/(produced[0] + produced[1]), + tot_produced[0]/tot_porevol_init); + if (wells_) { + wellreport.push(props_, *wells_, state.saturation(), + timer.currentTime() + timer.currentStepLength(), + well_state.bhp(), well_state.perfRates()); + } + } + total_timer.stop(); + + std::cout << "\n\n================ End of simulation ===============\n" + << "Total time taken: " << total_timer.secsSinceStart() + << "\n Pressure time: " << ptime + << "\n Transport time: " << ttime << std::endl; + + if (output_) { + outputState(grid_, state, timer.currentStepNum(), output_dir_); + outputWaterCut(watercut, output_dir_); + if (wells_) { + outputWellReport(wellreport, output_dir_); + } + } + + + + } + +} // namespace Opm diff --git a/opm/core/simulator/SimulatorTwophase.hpp b/opm/core/simulator/SimulatorTwophase.hpp new file mode 100644 index 00000000..053a3196 --- /dev/null +++ b/opm/core/simulator/SimulatorTwophase.hpp @@ -0,0 +1,95 @@ +/* + 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 . +*/ + +#ifndef OPM_SIMULATORTWOPHASE_HEADER_INCLUDED +#define OPM_SIMULATORTWOPHASE_HEADER_INCLUDED + +#include +#include + +struct UnstructuredGrid; +struct Wells; +struct FlowBoundaryConditions; + +namespace Opm +{ + namespace parameter { class ParameterGroup; } + class IncompPropertiesInterface; + class RockCompressibility; + class LinearSolverInterface; + class SimulatorTimer; + class TwophaseState; + class WellState; + + /// Class collecting all necessary components for a two-phase simulation. + class SimulatorTwophase + { + public: + /// Initialise from parameters and objects to observe. + /// \param[in] param parameters, this class accepts the following: + /// parameter (default) effect + /// ----------------------------------------------------------- + /// output (true) write output to files? + /// output_dir ("output") output directoty + /// output_interval (1) output every nth step + /// nl_pressure_maxiter (10) max nonlinear iterations in pressure + /// nl_pressure_tolerance (1.0) pressure solver nonlinear tolerance (in Pascal) + /// nl_maxiter (30) max nonlinear iterations in transport + /// nl_tolerance (1e-9) transport solver absolute residual tolerance + /// num_transport_substeps (1) number of transport steps per pressure step + /// use_segregation_split (false) solve for gravity segregation (if false, + /// segregation is ignored). + /// + /// \param[in] grid grid data structure + /// \param[in] props fluid and rock properties + /// \param[in] rock_comp if non-null, rock compressibility properties + /// \param[in] wells if non-null, wells data structure + /// \param[in] src source terms + /// \param[in] bcs boundary conditions, treat as all noflow if null + /// \param[in] linsolver linear solver + /// \param[in] gravity if non-null, gravity vector + SimulatorTwophase(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const IncompPropertiesInterface& props, + const RockCompressibility* rock_comp, + const Wells* wells, + const std::vector& src, + const FlowBoundaryConditions* bcs, + const LinearSolverInterface& linsolver, + const double* gravity); + + /// Run the simulation. + /// This will run succesive timesteps until timer.done() is true. It will + /// modify the reservoir and well states. + /// \param[in,out] timer governs the requested reporting timesteps + /// \param[in,out] state state of reservoir: pressure, fluxes + /// \param[in,out] well_state state of wells: bhp, perforation rates + void run(SimulatorTimer& timer, + TwophaseState& state, + WellState& well_state); + + private: + class Impl; + // Using shared_ptr instead of scoped_ptr since scoped_ptr requires complete type for Impl. + boost::shared_ptr pimpl_; + }; + +} // namespace Opm + +#endif // OPM_SIMULATORTWOPHASE_HEADER_INCLUDED diff --git a/opm/core/TwophaseState.hpp b/opm/core/simulator/TwophaseState.hpp similarity index 100% rename from opm/core/TwophaseState.hpp rename to opm/core/simulator/TwophaseState.hpp diff --git a/opm/core/WellState.hpp b/opm/core/simulator/WellState.hpp similarity index 91% rename from opm/core/WellState.hpp rename to opm/core/simulator/WellState.hpp index 56ede664..59d9f060 100644 --- a/opm/core/WellState.hpp +++ b/opm/core/simulator/WellState.hpp @@ -26,9 +26,11 @@ namespace Opm { + /// The state of a set of wells. class WellState { public: + /// Allocate and initialize if wells is non-null. template void init(const Wells* wells, const State& state) { @@ -44,9 +46,11 @@ namespace Opm } } + /// One bhp pressure per well. std::vector& bhp() { return bhp_; } const std::vector& bhp() const { return bhp_; } + /// One rate per well connection. std::vector& perfRates() { return perfrates_; } const std::vector& perfRates() const { return perfrates_; } diff --git a/opm/core/ColumnExtract.hpp b/opm/core/utility/ColumnExtract.hpp similarity index 100% rename from opm/core/ColumnExtract.hpp rename to opm/core/utility/ColumnExtract.hpp diff --git a/opm/core/InjectionSpecification.cpp b/opm/core/wells/InjectionSpecification.cpp similarity index 88% rename from opm/core/InjectionSpecification.cpp rename to opm/core/wells/InjectionSpecification.cpp index dca7fc76..accae78c 100644 --- a/opm/core/InjectionSpecification.cpp +++ b/opm/core/wells/InjectionSpecification.cpp @@ -1,4 +1,5 @@ -#include +#include + namespace Opm { diff --git a/opm/core/InjectionSpecification.hpp b/opm/core/wells/InjectionSpecification.hpp similarity index 100% rename from opm/core/InjectionSpecification.hpp rename to opm/core/wells/InjectionSpecification.hpp diff --git a/opm/core/ProductionSpecification.cpp b/opm/core/wells/ProductionSpecification.cpp similarity index 88% rename from opm/core/ProductionSpecification.cpp rename to opm/core/wells/ProductionSpecification.cpp index eb0f7d83..4d4bdf4d 100644 --- a/opm/core/ProductionSpecification.cpp +++ b/opm/core/wells/ProductionSpecification.cpp @@ -1,4 +1,4 @@ -#include +#include namespace Opm { diff --git a/opm/core/ProductionSpecification.hpp b/opm/core/wells/ProductionSpecification.hpp similarity index 100% rename from opm/core/ProductionSpecification.hpp rename to opm/core/wells/ProductionSpecification.hpp diff --git a/opm/core/WellCollection.cpp b/opm/core/wells/WellCollection.cpp similarity index 99% rename from opm/core/WellCollection.cpp rename to opm/core/wells/WellCollection.cpp index 76b311e7..727a7242 100644 --- a/opm/core/WellCollection.cpp +++ b/opm/core/wells/WellCollection.cpp @@ -18,7 +18,7 @@ You should have received a copy of the GNU General Public License along with OpenRS. If not, see . */ -#include "WellCollection.hpp" +#include namespace Opm { diff --git a/opm/core/WellCollection.hpp b/opm/core/wells/WellCollection.hpp similarity index 99% rename from opm/core/WellCollection.hpp rename to opm/core/wells/WellCollection.hpp index 62abd4ee..1957576a 100644 --- a/opm/core/WellCollection.hpp +++ b/opm/core/wells/WellCollection.hpp @@ -23,7 +23,7 @@ along with OpenRS. If not, see . #define OPM_WELLCOLLECTION_HPP #include -#include +#include #include #include diff --git a/opm/core/WellsGroup.cpp b/opm/core/wells/WellsGroup.cpp similarity index 99% rename from opm/core/WellsGroup.cpp rename to opm/core/wells/WellsGroup.cpp index 903cc9b2..3514bb5d 100644 --- a/opm/core/WellsGroup.cpp +++ b/opm/core/wells/WellsGroup.cpp @@ -17,7 +17,7 @@ along with OPM. If not, see . */ -#include +#include #include #include #include diff --git a/opm/core/WellsGroup.hpp b/opm/core/wells/WellsGroup.hpp similarity index 99% rename from opm/core/WellsGroup.hpp rename to opm/core/wells/WellsGroup.hpp index 3c14950d..f1b1afbb 100644 --- a/opm/core/WellsGroup.hpp +++ b/opm/core/wells/WellsGroup.hpp @@ -20,8 +20,8 @@ #ifndef OPM_WELLSGROUP_HPP #define OPM_WELLSGROUP_HPP -#include -#include +#include +#include #include #include #include diff --git a/opm/core/WellsManager.cpp b/opm/core/wells/WellsManager.cpp similarity index 98% rename from opm/core/WellsManager.cpp rename to opm/core/wells/WellsManager.cpp index 6d0d83dd..d305ee9f 100644 --- a/opm/core/WellsManager.cpp +++ b/opm/core/wells/WellsManager.cpp @@ -18,13 +18,13 @@ */ -#include +#include #include #include #include #include #include -#include +#include #include #include @@ -435,13 +435,13 @@ namespace Opm // Add all controls that are present in well. int ok = 1; int control_pos[5] = { -1, -1, -1, -1, -1 }; - if (ok && wci_line.surface_flow_max_rate_ > 0.0) { + if (ok && wci_line.surface_flow_max_rate_ >= 0.0) { control_pos[InjectionControl::RATE] = w_->ctrls[wix]->num; const double distr[3] = { 1.0, 1.0, 1.0 }; ok = append_well_controls(SURFACE_RATE, wci_line.surface_flow_max_rate_, distr, wix, w_); } - if (ok && wci_line.reservoir_flow_max_rate_ > 0.0) { + if (ok && wci_line.reservoir_flow_max_rate_ >= 0.0) { control_pos[InjectionControl::RESV] = w_->ctrls[wix]->num; const double distr[3] = { 1.0, 1.0, 1.0 }; ok = append_well_controls(RESERVOIR_RATE, wci_line.reservoir_flow_max_rate_, @@ -515,7 +515,7 @@ namespace Opm // Add all controls that are present in well. int control_pos[9] = { -1, -1, -1, -1, -1, -1, -1, -1, -1 }; int ok = 1; - if (ok && wcp_line.oil_max_rate_ > 0.0) { + if (ok && wcp_line.oil_max_rate_ >= 0.0) { if (!pu.phase_used[BlackoilPhases::Liquid]) { THROW("Oil phase not active and ORAT control specified."); } @@ -525,7 +525,7 @@ namespace Opm ok = append_well_controls(SURFACE_RATE, -wcp_line.oil_max_rate_, distr, wix, w_); } - if (ok && wcp_line.water_max_rate_ > 0.0) { + if (ok && wcp_line.water_max_rate_ >= 0.0) { if (!pu.phase_used[BlackoilPhases::Aqua]) { THROW("Water phase not active and WRAT control specified."); } @@ -535,7 +535,7 @@ namespace Opm ok = append_well_controls(SURFACE_RATE, -wcp_line.water_max_rate_, distr, wix, w_); } - if (ok && wcp_line.gas_max_rate_ > 0.0) { + if (ok && wcp_line.gas_max_rate_ >= 0.0) { if (!pu.phase_used[BlackoilPhases::Vapour]) { THROW("Gas phase not active and GRAT control specified."); } @@ -545,7 +545,7 @@ namespace Opm ok = append_well_controls(SURFACE_RATE, -wcp_line.gas_max_rate_, distr, wix, w_); } - if (ok && wcp_line.liquid_max_rate_ > 0.0) { + if (ok && wcp_line.liquid_max_rate_ >= 0.0) { if (!pu.phase_used[BlackoilPhases::Aqua]) { THROW("Water phase not active and LRAT control specified."); } @@ -559,7 +559,7 @@ namespace Opm ok = append_well_controls(SURFACE_RATE, -wcp_line.liquid_max_rate_, distr, wix, w_); } - if (ok && wcp_line.reservoir_flow_max_rate_ > 0.0) { + if (ok && wcp_line.reservoir_flow_max_rate_ >= 0.0) { control_pos[ProductionControl::RESV] = w_->ctrls[wix]->num; double distr[3] = { 1.0, 1.0, 1.0 }; ok = append_well_controls(RESERVOIR_RATE, -wcp_line.reservoir_flow_max_rate_, diff --git a/opm/core/WellsManager.hpp b/opm/core/wells/WellsManager.hpp similarity index 98% rename from opm/core/WellsManager.hpp rename to opm/core/wells/WellsManager.hpp index 3fbdbd1a..6df65ec4 100644 --- a/opm/core/WellsManager.hpp +++ b/opm/core/wells/WellsManager.hpp @@ -21,8 +21,8 @@ #define OPM_WELLSMANAGER_HEADER_INCLUDED -#include -#include +#include +#include struct Wells; struct UnstructuredGrid; diff --git a/tests/Makefile.am b/tests/Makefile.am index a279e72d..064d385c 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -47,3 +47,9 @@ noinst_PROGRAMS += test_cfs_tpfa noinst_PROGRAMS += test_jacsys test_jacsys_SOURCES = test_jacsys.cpp endif + + +if BUILD_AGMG +noinst_PROGRAMS += test_agmg +test_agmg_SOURCES = test_agmg.cpp +endif diff --git a/tests/test_agmg.cpp b/tests/test_agmg.cpp new file mode 100644 index 00000000..51c4a9cf --- /dev/null +++ b/tests/test_agmg.cpp @@ -0,0 +1,126 @@ +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + Copyright 2012 Statoil ASA. + + 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 . +*/ + +#include +#include +#include +#include +#include +#include + +#include + +#include +#include + +namespace { + std::size_t + compute_nnz(const std::size_t m) + { + assert (m > 0); + + std::size_t nnz = m; // A(i,i) + + nnz += (m > 1) ? 2 : 0; // A(0,1), A(m-1,m-2) + nnz += (m > 2) ? 2 * (m - 2) : 0; // A(i,i-1), A(i,i+1) + + return nnz; + } + + + boost::shared_ptr + build_laplace_1d(const std::size_t m) + { + assert (m >= 2); + + const std::size_t nnz = compute_nnz(m); + + boost::shared_ptr + A(csrmatrix_new_known_nnz(m, nnz), csrmatrix_delete); + + A->ia[ 0 ] = 0; + + // First row + A->ia[ 0 + 1 ] = A->ia[ 0 ]; + A->ja[ A->ia[0 + 1] ] = 0 + 0; + A->sa[ A->ia[0 + 1] ++ ] = 2.0; + A->ja[ A->ia[0 + 1] ] = 0 + 1; + A->sa[ A->ia[0 + 1] ++ ] = - 1.0; + + // General rows + for (std::size_t i = 1; i < m - 1; ++i) { + A->ia[i + 1] = A->ia[i]; + + A->ja[ A->ia[i + 1] ] = int(i) - 1; + A->sa[ A->ia[i + 1] ++ ] = - 1.0; + + A->ja[ A->ia[i + 1] ] = int(i) ; + A->sa[ A->ia[i + 1] ++ ] = 2.0; + + A->ja[ A->ia[i + 1] ] = int(i) + 1; + A->sa[ A->ia[i + 1] ++ ] = - 1.0; + } + + // Last row + A->ia[ (m - 1) + 1 ] = A->ia[ m - 1 ]; + + A->ja[ A->ia[ (m - 1) + 1 ] ] = int(m - 1) - 1; + A->sa[ A->ia[ (m - 1) + 1 ] ++ ] = - 1.0; + + A->ja[ A->ia[ (m - 1) + 1 ] ] = int(m - 1) ; + A->sa[ A->ia[ (m - 1) + 1 ] ++ ] = 2.0; + + return A; + } +} + + +int main() +{ + const std::size_t m = 10; + + boost::shared_ptr A = build_laplace_1d(m); + + // Form right-hand side [1, 0, 0, ...., 0, 1] + std::vector b(m, 0.0); + b[0] = 1.0; b.back() = 1.0; + + // Allocate solution vector + std::vector x(m); + + // Create solver for SPD system. + Opm::LinearSolverAGMG linsolve(100, 1e-9, true); + + Opm::LinearSolverInterface::LinearSolverReport + rpt = linsolve.solve(A.get(), & b[0], & x[0]); + + double e = 0.0; + for (std::size_t i = 0; i < m; ++i) { + const double d = x[i] - 1.0; + e += d * d; + } + + std::cerr << "|| e ||_2 = " + << std::scientific + << std::setprecision(5) + << std::sqrt(e) / double(m) << '\n'; + + return 0; +} diff --git a/tests/test_column_extract.cpp b/tests/test_column_extract.cpp index 877f21a0..46b1d538 100644 --- a/tests/test_column_extract.cpp +++ b/tests/test_column_extract.cpp @@ -6,7 +6,7 @@ #define BOOST_TEST_MODULE ColumnExtractTest #include -#include +#include #include #include diff --git a/tutorials/tutorial3.cpp b/tutorials/tutorial3.cpp index 6420c199..fc15ab30 100644 --- a/tutorials/tutorial3.cpp +++ b/tutorials/tutorial3.cpp @@ -37,7 +37,7 @@ #include -#include +#include #include #include diff --git a/tutorials/tutorial4.cpp b/tutorials/tutorial4.cpp index b3996025..7bdbd1b9 100644 --- a/tutorials/tutorial4.cpp +++ b/tutorials/tutorial4.cpp @@ -44,12 +44,12 @@ #include #include -#include +#include #include #include #include -#include +#include /// \page tutorial4 Well controls ///