From 6e2cdfc33f65a04b8098b27a38497830d21deca1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 18 Sep 2013 14:32:09 +0200 Subject: [PATCH 01/14] Add helpers grad, fullngrad, fulldiv. --- opm/autodiff/AutoDiffHelpers.hpp | 25 +++++++++++++++++++++++-- 1 file changed, 23 insertions(+), 2 deletions(-) diff --git a/opm/autodiff/AutoDiffHelpers.hpp b/opm/autodiff/AutoDiffHelpers.hpp index 0e4c9f90b..30e4db08c 100644 --- a/opm/autodiff/AutoDiffHelpers.hpp +++ b/opm/autodiff/AutoDiffHelpers.hpp @@ -39,12 +39,19 @@ struct HelperOps typedef Eigen::Array IFaces; IFaces internal_faces; - /// Extract for each face the difference of its adjacent cells'values. + /// Extract for each internal face the difference of its adjacent cells'values (first - second). M ngrad; + /// Extract for each face the difference of its adjacent cells'values (second - first). + M grad; /// Extract for each face the average of its adjacent cells' values. M caver; - /// Extract for each cell the sum of its adjacent faces' (signed) values. + /// Extract for each cell the sum of its adjacent interior faces' (signed) values. M div; + /// Extract for each face the difference of its adjacent cells'values (first - second). + /// For boundary faces, one of the entries per row (corresponding to the outside) is zero. + M fullngrad; + /// Extract for each cell the sum of all its adjacent faces' (signed) values. + M fulldiv; /// Constructs all helper vectors and matrices. HelperOps(const UnstructuredGrid& grid) @@ -89,7 +96,21 @@ struct HelperOps } ngrad.setFromTriplets(ngrad_tri.begin(), ngrad_tri.end()); caver.setFromTriplets(caver_tri.begin(), caver_tri.end()); + grad = -ngrad; div = ngrad.transpose(); + std::vector fullngrad_tri; + fullngrad_tri.reserve(2*nf); + for (int i = 0; i < nf; ++i) { + if (nb(i,0) >= 0) { + fullngrad_tri.emplace_back(i, nb(i,0), 1.0); + } + if (nb(i,1) >= 0) { + fullngrad_tri.emplace_back(i, nb(i,1), -1.0); + } + } + fullngrad.resize(nf, nc); + fullngrad.setFromTriplets(fullngrad_tri.begin(), fullngrad_tri.end()); + fulldiv = fullngrad.transpose(); } }; From a33f7e964b8cae2d9ac5e3733bd80ce1fd94cd7a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 19 Sep 2013 10:09:53 +0200 Subject: [PATCH 02/14] Let sim_2p_comp_ad throw if not given an input deck. There is some code in place now to create wells for the no-deck case, but since it does not work correctly yet, the simulator intercepts this and throws. --- examples/sim_2p_comp_ad.cpp | 37 ++++++++++++++++-------- opm/autodiff/ImpesTPFAAD.cpp | 2 +- opm/autodiff/SimulatorCompressibleAd.cpp | 9 ++---- opm/autodiff/SimulatorCompressibleAd.hpp | 4 +-- 4 files changed, 29 insertions(+), 23 deletions(-) diff --git a/examples/sim_2p_comp_ad.cpp b/examples/sim_2p_comp_ad.cpp index 2bda013f6..a511d5d95 100644 --- a/examples/sim_2p_comp_ad.cpp +++ b/examples/sim_2p_comp_ad.cpp @@ -81,6 +81,11 @@ try // If we have a "deck_filename", grid and props will be read from that. bool use_deck = param.has("deck_filename"); + if (!use_deck) { + // This check should be removed when and if this simulator is verified and works without decks. + // The current code for the non-deck case fails for unknown reasons. + OPM_THROW(std::runtime_error, "This simulator cannot run without a deck with wells. Use deck_filename to specify deck."); + } boost::scoped_ptr deck; boost::scoped_ptr grid; boost::scoped_ptr props; @@ -132,12 +137,9 @@ try bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0); const double *grav = use_gravity ? &gravity[0] : 0; - // Initialising src - int num_cells = grid->c_grid()->number_of_cells; - std::vector src(num_cells, 0.0); - if (use_deck) { - // Do nothing, wells will be the driving force, not source terms. - } else { + // Initialising wells if not from deck. + Wells* simple_wells = 0; + if (!use_deck) { // Compute pore volumes, in order to enable specifying injection rate // terms of total pore volume. std::vector porevol; @@ -150,8 +152,21 @@ try 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/unit::day; - src[0] = flow_per_sec; - src[num_cells - 1] = -flow_per_sec; + simple_wells = create_wells(2, 2, 2); + const double inj_frac[2] = { 1.0, 0.0 }; + const int inj_cell = 0; + const double WI = 1e-8; + const double all_fluids[2] = { 1.0, 1.0 }; + int ok = add_well(INJECTOR, 0.0, 1, inj_frac, &inj_cell, &WI, "Injector", simple_wells); + ok = ok && append_well_controls(SURFACE_RATE, 0.01*flow_per_sec, all_fluids, 0, simple_wells); + const int prod_cell = grid->c_grid()->number_of_cells - 1; + ok = ok && add_well(PRODUCER, 0.0, 1, NULL, &prod_cell, &WI, "Producer", simple_wells); + ok = ok && append_well_controls(SURFACE_RATE, -0.01*flow_per_sec, all_fluids, 1, simple_wells); + if (!ok) { + OPM_THROW(std::runtime_error, "Simple well init failed."); + } + simple_wells->ctrls[0]->current = 0; + simple_wells->ctrls[1]->current = 0; } // Boundary conditions. @@ -196,13 +211,12 @@ try SimulatorReport rep; if (!use_deck) { // Simple simulation without a deck. - WellsManager wells; // no wells. + WellsManager wells(simple_wells); SimulatorCompressibleAd simulator(param, *grid->c_grid(), *props, rock_comp->isActive() ? rock_comp.get() : 0, wells, - src, bcs.c_bcs(), linsolver, grav); @@ -210,7 +224,7 @@ try simtimer.init(param); warnIfUnusedParams(param); WellState well_state; - well_state.init(0, state); + well_state.init(simple_wells, state); rep = simulator.run(simtimer, state, well_state); } else { // With a deck, we may have more epochs etc. @@ -257,7 +271,6 @@ try *props, rock_comp->isActive() ? rock_comp.get() : 0, wells, - src, bcs.c_bcs(), linsolver, grav); diff --git a/opm/autodiff/ImpesTPFAAD.cpp b/opm/autodiff/ImpesTPFAAD.cpp index 9fdf52685..4ffac1b34 100644 --- a/opm/autodiff/ImpesTPFAAD.cpp +++ b/opm/autodiff/ImpesTPFAAD.cpp @@ -183,7 +183,7 @@ namespace Opm { well_kr_ = fluid_.relperm(well_s.col(0), well_s.col(1), V::Zero(nperf,1), well_cells); const double atol = 1.0e-10; - const double rtol = 5.0e-8; + const double rtol = 5.0e-6; const int maxit = 15; assemble(dt, state, well_state); diff --git a/opm/autodiff/SimulatorCompressibleAd.cpp b/opm/autodiff/SimulatorCompressibleAd.cpp index fdd1977f2..c21ef8acf 100644 --- a/opm/autodiff/SimulatorCompressibleAd.cpp +++ b/opm/autodiff/SimulatorCompressibleAd.cpp @@ -70,7 +70,6 @@ namespace Opm const BlackoilPropertiesInterface& props, const RockCompressibility* rock_comp_props, WellsManager& wells_manager, - const std::vector& src, const FlowBoundaryConditions* bcs, LinearSolverInterface& linsolver, const double* gravity); @@ -99,7 +98,6 @@ namespace Opm const RockCompressibility* rock_comp_props_; WellsManager& wells_manager_; const Wells* wells_; - const std::vector& src_; const FlowBoundaryConditions* bcs_; const double* gravity_; // Solvers @@ -121,12 +119,11 @@ namespace Opm const BlackoilPropertiesInterface& props, const RockCompressibility* rock_comp_props, WellsManager& wells_manager, - const std::vector& src, const FlowBoundaryConditions* bcs, LinearSolverInterface& linsolver, const double* gravity) { - pimpl_.reset(new Impl(param, grid, props, rock_comp_props, wells_manager, src, bcs, linsolver, gravity)); + pimpl_.reset(new Impl(param, grid, props, rock_comp_props, wells_manager, bcs, linsolver, gravity)); } @@ -235,13 +232,12 @@ namespace Opm - // \TODO: make CompressibleTpfa take src and bcs. + // \TODO: make CompressibleTpfa take bcs. SimulatorCompressibleAd::Impl::Impl(const parameter::ParameterGroup& param, const UnstructuredGrid& grid, const BlackoilPropertiesInterface& props, const RockCompressibility* rock_comp_props, WellsManager& wells_manager, - const std::vector& src, const FlowBoundaryConditions* bcs, LinearSolverInterface& linsolver, const double* gravity) @@ -250,7 +246,6 @@ namespace Opm rock_comp_props_(rock_comp_props), wells_manager_(wells_manager), wells_(wells_manager.c_wells()), - src_(src), bcs_(bcs), gravity_(gravity), fluid_(props_), diff --git a/opm/autodiff/SimulatorCompressibleAd.hpp b/opm/autodiff/SimulatorCompressibleAd.hpp index 98f5470f9..0936ad346 100644 --- a/opm/autodiff/SimulatorCompressibleAd.hpp +++ b/opm/autodiff/SimulatorCompressibleAd.hpp @@ -62,8 +62,7 @@ namespace Opm /// \param[in] grid grid data structure /// \param[in] props fluid and rock properties /// \param[in] rock_comp_props if non-null, rock compressibility properties - /// \param[in] well_manager well manager, may manage no (null) wells - /// \param[in] src source terms + /// \param[in] well_manager well manager /// \param[in] bcs boundary conditions, treat as all noflow if null /// \param[in] linsolver linear solver /// \param[in] gravity if non-null, gravity vector @@ -72,7 +71,6 @@ namespace Opm const BlackoilPropertiesInterface& props, const RockCompressibility* rock_comp_props, WellsManager& wells_manager, - const std::vector& src, const FlowBoundaryConditions* bcs, LinearSolverInterface& linsolver, const double* gravity); From 86e9e04d2f8d608866d6fe2eca7ca9652eb9ac74 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 19 Sep 2013 11:32:29 +0200 Subject: [PATCH 03/14] Remove unnecessary include statement. --- opm/autodiff/AutoDiffBlock.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/opm/autodiff/AutoDiffBlock.hpp b/opm/autodiff/AutoDiffBlock.hpp index b0f03cbc2..f1cd9d86b 100644 --- a/opm/autodiff/AutoDiffBlock.hpp +++ b/opm/autodiff/AutoDiffBlock.hpp @@ -20,7 +20,6 @@ #ifndef OPM_AUTODIFFBLOCK_HEADER_INCLUDED #define OPM_AUTODIFFBLOCK_HEADER_INCLUDED -#include #include #include #include From e9b933bf4f249c5f53119d7cdf5298d376efa3b2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 19 Sep 2013 11:32:47 +0200 Subject: [PATCH 04/14] Rename AutoDiff::Forward -> Opm::AutoDiff. --- examples/find_zero.cpp | 2 +- opm/autodiff/AutoDiff.hpp | 159 +++++++++++++++++++------------------- tests/test_syntax.cpp | 16 ++-- 3 files changed, 88 insertions(+), 89 deletions(-) diff --git a/examples/find_zero.cpp b/examples/find_zero.cpp index 45f4a9920..bae50e2b8 100644 --- a/examples/find_zero.cpp +++ b/examples/find_zero.cpp @@ -80,7 +80,7 @@ public: { double x = initial_guess; iterations_used = 0; - typedef AutoDiff::Forward AD; + typedef Opm::AutoDiff AD; while (std::abs(f(x)) > tolerance && ++iterations_used < max_iter) { AD xfad = AD::variable(x); AD rfad = f(xfad); diff --git a/opm/autodiff/AutoDiff.hpp b/opm/autodiff/AutoDiff.hpp index 8a758c453..a542abe89 100644 --- a/opm/autodiff/AutoDiff.hpp +++ b/opm/autodiff/AutoDiff.hpp @@ -1,18 +1,3 @@ -/*=========================================================================== -// -// File: AutoDiff.hpp -// -// Created: 2013-04-29 10:51:23+0200 -// -// Authors: Knut-Andreas Lie -// Halvor M. Nilsen -// Atgeirr F. Rasmussen -// Xavier Raynaud -// Bård Skaflestad -// -//==========================================================================*/ - - /* Copyright 2013 SINTEF ICT, Applied Mathematics. Copyright 2013 Statoil ASA. @@ -33,33 +18,46 @@ along with OPM. If not, see . */ -#include - #ifndef OPM_AUTODIFF_HPP_HEADER #define OPM_AUTODIFF_HPP_HEADER -namespace AutoDiff { +#include + +namespace Opm +{ + + /// A simple class for forward-mode automatic differentiation. + /// + /// The class represents a single value and a single derivative. + /// Only basic arithmetic operators and a few functions are + /// implemented for it, it is mostly intended for simple + /// experimentation. template - class Forward { + class AutoDiff + { public: - static Forward + /// Create an AutoDiff object representing a constant, that + /// is, its derivative is zero. + static AutoDiff constant(const Scalar x) { - // Constant is function with zero derivative. return function(x, Scalar(0)); } - static Forward + /// Create an AutoDiff object representing a primary variable, + /// that is, its derivative is one. + static AutoDiff variable(const Scalar x) { - // Variable is function with unit derivative (wrt. itself). return function(x, Scalar(1)); } - static Forward + /// Create an AutoDiff object representing a function value + /// and its derivative. + static AutoDiff function(const Scalar x, const Scalar dx) { - return Forward(x, dx); + return AutoDiff(x, dx); } void @@ -69,7 +67,7 @@ namespace AutoDiff { } void - operator +=(const Forward& rhs) + operator +=(const AutoDiff& rhs) { val_ += rhs.val_; der_ += rhs.der_; @@ -82,7 +80,7 @@ namespace AutoDiff { } void - operator -=(const Forward& rhs) + operator -=(const AutoDiff& rhs) { val_ -= rhs.val_; der_ -= rhs.der_; @@ -96,7 +94,7 @@ namespace AutoDiff { } void - operator *=(const Forward& rhs) + operator *=(const AutoDiff& rhs) { der_ = der_*rhs.val_ + val_*rhs.der_; val_ *= rhs.val_; @@ -110,7 +108,7 @@ namespace AutoDiff { } void - operator /=(const Forward& rhs) + operator /=(const AutoDiff& rhs) { der_ = (der_*rhs.val_ - val_*rhs.der_) / (rhs.val_ * rhs.val_); val_ /= rhs.val_; @@ -129,7 +127,7 @@ namespace AutoDiff { const Scalar der() const { return der_; } private: - Forward(const Scalar x, const Scalar dx) + AutoDiff(const Scalar x, const Scalar dx) : val_(x), der_(dx) {} @@ -140,164 +138,165 @@ namespace AutoDiff { template Ostream& - operator<<(Ostream& os, const Forward& fw) + operator<<(Ostream& os, const AutoDiff& fw) { return fw.print(os); } template - Forward - operator +(const Forward& lhs, - const Forward& rhs) + AutoDiff + operator +(const AutoDiff& lhs, + const AutoDiff& rhs) { - Forward ret = lhs; + AutoDiff ret = lhs; ret += rhs; return ret; } template - Forward + AutoDiff operator +(const T lhs, - const Forward& rhs) + const AutoDiff& rhs) { - Forward ret = rhs; + AutoDiff ret = rhs; ret += Scalar(lhs); return ret; } template - Forward - operator +(const Forward& lhs, + AutoDiff + operator +(const AutoDiff& lhs, const T rhs) { - Forward ret = lhs; + AutoDiff ret = lhs; ret += Scalar(rhs); return ret; } template - Forward - operator -(const Forward& lhs, - const Forward& rhs) + AutoDiff + operator -(const AutoDiff& lhs, + const AutoDiff& rhs) { - Forward ret = lhs; + AutoDiff ret = lhs; ret -= rhs; return ret; } template - Forward + AutoDiff operator -(const T lhs, - const Forward& rhs) + const AutoDiff& rhs) { - return Forward::function(Scalar(lhs) - rhs.val(), -rhs.der()); + return AutoDiff::function(Scalar(lhs) - rhs.val(), -rhs.der()); } template - Forward - operator -(const Forward& lhs, + AutoDiff + operator -(const AutoDiff& lhs, const T rhs) { - Forward ret = lhs; + AutoDiff ret = lhs; ret -= Scalar(rhs); return ret; } template - Forward - operator *(const Forward& lhs, - const Forward& rhs) + AutoDiff + operator *(const AutoDiff& lhs, + const AutoDiff& rhs) { - Forward ret = lhs; + AutoDiff ret = lhs; ret *= rhs; return ret; } template - Forward + AutoDiff operator *(const T lhs, - const Forward& rhs) + const AutoDiff& rhs) { - Forward ret = rhs; + AutoDiff ret = rhs; ret *= Scalar(lhs); return ret; } template - Forward - operator *(const Forward& lhs, + AutoDiff + operator *(const AutoDiff& lhs, const T rhs) { - Forward ret = lhs; + AutoDiff ret = lhs; ret *= Scalar(rhs); return ret; } template - Forward - operator /(const Forward& lhs, - const Forward& rhs) + AutoDiff + operator /(const AutoDiff& lhs, + const AutoDiff& rhs) { - Forward ret = lhs; + AutoDiff ret = lhs; ret /= rhs; return ret; } template - Forward + AutoDiff operator /(const T lhs, - const Forward& rhs) + const AutoDiff& rhs) { Scalar a = Scalar(lhs) / rhs.val(); Scalar b = -Scalar(lhs) / (rhs.val() * rhs.val()); - return Forward::function(a, b); + return AutoDiff::function(a, b); } template - Forward - operator /(const Forward& lhs, + AutoDiff + operator /(const AutoDiff& lhs, const T rhs) { Scalar a = lhs.val() / Scalar(rhs); Scalar b = lhs.der() / Scalar(rhs); - return Forward::function(a, b); + return AutoDiff::function(a, b); } template - Forward - cos(const Forward& x) + AutoDiff + cos(const AutoDiff& x) { Scalar a = std::cos(x.val()); Scalar b = -std::sin(x.val()) * x.der(); - return Forward::function(a, b); + return AutoDiff::function(a, b); } template - Forward - sqrt(const Forward& x) + AutoDiff + sqrt(const AutoDiff& x) { Scalar a = std::sqrt(x.val()); Scalar b = (Scalar(1.0) / (Scalar(2.0) * a)) * x.der(); - return Forward::function(a, b); + return AutoDiff::function(a, b); } -} // namespace AutoDiff + +} // namespace Opm namespace std { - using AutoDiff::cos; - using AutoDiff::sqrt; + using Opm::cos; + using Opm::sqrt; } #endif /* OPM_AUTODIFF_HPP_HEADER */ diff --git a/tests/test_syntax.cpp b/tests/test_syntax.cpp index b3c23ee4f..e8d5bc044 100644 --- a/tests/test_syntax.cpp +++ b/tests/test_syntax.cpp @@ -15,7 +15,7 @@ BOOST_AUTO_TEST_CASE(Initialisation) { - typedef AutoDiff::Forward AdFW; + typedef Opm::AutoDiff AdFW; const double atol = 1.0e-14; @@ -36,7 +36,7 @@ BOOST_AUTO_TEST_CASE(Initialisation) BOOST_AUTO_TEST_CASE(Addition) { - typedef AutoDiff::Forward AdFW; + typedef Opm::AutoDiff AdFW; const double atol = 1.0e-14; @@ -71,7 +71,7 @@ BOOST_AUTO_TEST_CASE(Addition) BOOST_AUTO_TEST_CASE(Subtraction) { - typedef AutoDiff::Forward AdFW; + typedef Opm::AutoDiff AdFW; const double atol = 1.0e-14; @@ -110,7 +110,7 @@ BOOST_AUTO_TEST_CASE(Subtraction) BOOST_AUTO_TEST_CASE(Multiplication) { - typedef AutoDiff::Forward AdFW; + typedef Opm::AutoDiff AdFW; const double atol = 1.0e-14; @@ -149,7 +149,7 @@ BOOST_AUTO_TEST_CASE(Multiplication) BOOST_AUTO_TEST_CASE(Division) { - typedef AutoDiff::Forward AdFW; + typedef Opm::AutoDiff AdFW; const double atol = 1.0e-14; @@ -186,7 +186,7 @@ BOOST_AUTO_TEST_CASE(Division) BOOST_AUTO_TEST_CASE(Polynomial) { - typedef AutoDiff::Forward AdFW; + typedef Opm::AutoDiff AdFW; const double atol = 1.0e-14; @@ -204,7 +204,7 @@ BOOST_AUTO_TEST_CASE(Polynomial) BOOST_AUTO_TEST_CASE(Cosine) { - typedef AutoDiff::Forward AdFW; + typedef Opm::AutoDiff AdFW; const double atol = 1.0e-14; @@ -223,7 +223,7 @@ BOOST_AUTO_TEST_CASE(Cosine) BOOST_AUTO_TEST_CASE(SquareRoot) { - typedef AutoDiff::Forward AdFW; + typedef Opm::AutoDiff AdFW; const double atol = 1.0e-14; From 85f79c0e846041621ca4df28ec6d7d127fc2d134 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 19 Sep 2013 12:53:28 +0200 Subject: [PATCH 05/14] Rename AutoDiff::ForwardBlock -> Opm::AutoDiffBlock. Also moved AutoDiffHelpers.hpp content to Opm namespace, and modified other files as required by these two changes. --- examples/sim_simple.cpp | 8 +- examples/test_impestpfa_ad.cpp | 2 - opm/autodiff/AutoDiffBlock.hpp | 110 +++++++++---------- opm/autodiff/AutoDiffHelpers.hpp | 48 ++++---- opm/autodiff/BlackoilPropsAd.hpp | 2 +- opm/autodiff/BlackoilPropsAdFromDeck.hpp | 2 +- opm/autodiff/BlackoilPropsAdInterface.hpp | 2 +- opm/autodiff/FullyImplicitBlackoilSolver.cpp | 15 ++- opm/autodiff/FullyImplicitBlackoilSolver.hpp | 2 +- opm/autodiff/ImpesTPFAAD.cpp | 10 +- opm/autodiff/ImpesTPFAAD.hpp | 2 +- opm/autodiff/TransportSolverTwophaseAd.hpp | 2 +- tests/test_block.cpp | 12 +- tests/test_scalar_mult.cpp | 4 +- tests/test_span.cpp | 2 + 15 files changed, 118 insertions(+), 105 deletions(-) diff --git a/examples/sim_simple.cpp b/examples/sim_simple.cpp index 9ddc29832..b8403b7ef 100644 --- a/examples/sim_simple.cpp +++ b/examples/sim_simple.cpp @@ -111,7 +111,7 @@ fluxFunc(const std::vector& m) int main() try { - typedef AutoDiff::ForwardBlock ADB; + typedef Opm::AutoDiffBlock ADB; typedef ADB::V V; typedef ADB::M M; @@ -148,11 +148,11 @@ try std::cerr << "Opm core " << clock.secsSinceLast() << std::endl; // Define neighbourhood-derived operator matrices. - const HelperOps ops(grid); + const Opm::HelperOps ops(grid); const int num_internal = ops.internal_faces.size(); std::cerr << "Topology matrices " << clock.secsSinceLast() << std::endl; - typedef AutoDiff::ForwardBlock ADB; + typedef Opm::AutoDiffBlock ADB; typedef ADB::V V; // q @@ -248,7 +248,7 @@ try V sw1 = 0.5*V::Ones(nc,1); const V ndp = (ops.ngrad * p1.matrix()).array(); const V dflux = mobtransf * ndp; - const UpwindSelector upwind(grid, ops, dflux); + const Opm::UpwindSelector upwind(grid, ops, dflux); const V pv = Eigen::Map(props.porosity(), nc, 1) * Eigen::Map(grid.cell_volumes, nc, 1); const double dt = 0.0005; diff --git a/examples/test_impestpfa_ad.cpp b/examples/test_impestpfa_ad.cpp index 4179182a2..d0252694d 100644 --- a/examples/test_impestpfa_ad.cpp +++ b/examples/test_impestpfa_ad.cpp @@ -58,8 +58,6 @@ try const Opm::BlackoilPropertiesBasic oldprops(param, 2, nc); const Opm::BlackoilPropsAd props(oldprops); - typedef AutoDiff::ForwardBlock ADB; - Wells* wells = create_wells(2, 2, 5); const double inj_frac[] = { 1.0, 0.0 }; const double prod_frac[] = { 0.0, 0.0 }; diff --git a/opm/autodiff/AutoDiffBlock.hpp b/opm/autodiff/AutoDiffBlock.hpp index f1cd9d86b..ef73bd5a2 100644 --- a/opm/autodiff/AutoDiffBlock.hpp +++ b/opm/autodiff/AutoDiffBlock.hpp @@ -25,11 +25,11 @@ #include #include -namespace AutoDiff +namespace Opm { template - class ForwardBlock + class AutoDiffBlock { public: /// Underlying types for scalar vectors and jacobians. @@ -37,14 +37,14 @@ namespace AutoDiff typedef Eigen::SparseMatrix M; /// Named constructor pattern used here. - static ForwardBlock null() + static AutoDiffBlock null() { V val; std::vector jac; - return ForwardBlock(val, jac); + return AutoDiffBlock(val, jac); } - static ForwardBlock constant(const V& val, const std::vector& blocksizes) + static AutoDiffBlock constant(const V& val, const std::vector& blocksizes) { std::vector jac; const int num_elem = val.size(); @@ -54,10 +54,10 @@ namespace AutoDiff for (int i = 0; i < num_blocks; ++i) { jac[i] = M(num_elem, blocksizes[i]); } - return ForwardBlock(val, jac); + return AutoDiffBlock(val, jac); } - static ForwardBlock variable(const int index, const V& val, const std::vector& blocksizes) + static AutoDiffBlock variable(const int index, const V& val, const std::vector& blocksizes) { std::vector jac; const int num_elem = val.size(); @@ -73,24 +73,24 @@ namespace AutoDiff for (typename M::Index row = 0; row < val.size(); ++row) { jac[index].insert(row, row) = Scalar(1.0); } - return ForwardBlock(val, jac); + return AutoDiffBlock(val, jac); } - static ForwardBlock function(const V& val, const std::vector& jac) + static AutoDiffBlock function(const V& val, const std::vector& jac) { - return ForwardBlock(val, jac); + return AutoDiffBlock(val, jac); } /// Construct a set of primary variables, /// each initialized to a given vector. - static std::vector variables(const std::vector& initial_values) + static std::vector variables(const std::vector& initial_values) { const int num_vars = initial_values.size(); std::vector bpat; for (int v = 0; v < num_vars; ++v) { bpat.push_back(initial_values[v].size()); } - std::vector vars; + std::vector vars; for (int v = 0; v < num_vars; ++v) { vars.emplace_back(variable(v, initial_values[v], bpat)); } @@ -98,7 +98,7 @@ namespace AutoDiff } /// Operator += - ForwardBlock& operator+=(const ForwardBlock& rhs) + AutoDiffBlock& operator+=(const AutoDiffBlock& rhs) { assert (numBlocks() == rhs.numBlocks()); assert (value().size() == rhs.value().size()); @@ -116,7 +116,7 @@ namespace AutoDiff } /// Operator + - ForwardBlock operator+(const ForwardBlock& rhs) const + AutoDiffBlock operator+(const AutoDiffBlock& rhs) const { std::vector jac = jac_; assert(numBlocks() == rhs.numBlocks()); @@ -130,7 +130,7 @@ namespace AutoDiff } /// Operator - - ForwardBlock operator-(const ForwardBlock& rhs) const + AutoDiffBlock operator-(const AutoDiffBlock& rhs) const { std::vector jac = jac_; assert(numBlocks() == rhs.numBlocks()); @@ -144,7 +144,7 @@ namespace AutoDiff } /// Operator * - ForwardBlock operator*(const ForwardBlock& rhs) const + AutoDiffBlock operator*(const AutoDiffBlock& rhs) const { int num_blocks = numBlocks(); std::vector jac(num_blocks); @@ -161,7 +161,7 @@ namespace AutoDiff } /// Operator / - ForwardBlock operator/(const ForwardBlock& rhs) const + AutoDiffBlock operator/(const AutoDiffBlock& rhs) const { int num_blocks = numBlocks(); std::vector jac(num_blocks); @@ -226,8 +226,8 @@ namespace AutoDiff } private: - ForwardBlock(const V& val, - const std::vector& jac) + AutoDiffBlock(const V& val, + const std::vector& jac) : val_(val), jac_(jac) { #ifndef NDEBUG @@ -241,12 +241,12 @@ namespace AutoDiff V val_; std::vector jac_; - }; + }; template Ostream& - operator<<(Ostream& os, const ForwardBlock& fw) + operator<<(Ostream& os, const AutoDiffBlock& fw) { return fw.print(os); } @@ -254,81 +254,81 @@ namespace AutoDiff /// Multiply with sparse matrix from the left. template - ForwardBlock operator*(const typename ForwardBlock::M& lhs, - const ForwardBlock& rhs) + AutoDiffBlock operator*(const typename AutoDiffBlock::M& lhs, + const AutoDiffBlock& rhs) { int num_blocks = rhs.numBlocks(); - std::vector::M> jac(num_blocks); + std::vector::M> jac(num_blocks); assert(lhs.cols() == rhs.value().rows()); for (int block = 0; block < num_blocks; ++block) { jac[block] = lhs*rhs.derivative()[block]; } - typename ForwardBlock::V val = lhs*rhs.value().matrix(); - return ForwardBlock::function(val, jac); + typename AutoDiffBlock::V val = lhs*rhs.value().matrix(); + return AutoDiffBlock::function(val, jac); } template - ForwardBlock operator*(const typename ForwardBlock::V& lhs, - const ForwardBlock& rhs) + AutoDiffBlock operator*(const typename AutoDiffBlock::V& lhs, + const AutoDiffBlock& rhs) { - return ForwardBlock::constant(lhs, rhs.blockPattern()) * rhs; + return AutoDiffBlock::constant(lhs, rhs.blockPattern()) * rhs; } template - ForwardBlock operator*(const ForwardBlock& lhs, - const typename ForwardBlock::V& rhs) + AutoDiffBlock operator*(const AutoDiffBlock& lhs, + const typename AutoDiffBlock::V& rhs) { return rhs * lhs; // Commutative operation. } template - ForwardBlock operator+(const typename ForwardBlock::V& lhs, - const ForwardBlock& rhs) + AutoDiffBlock operator+(const typename AutoDiffBlock::V& lhs, + const AutoDiffBlock& rhs) { - return ForwardBlock::constant(lhs, rhs.blockPattern()) + rhs; + return AutoDiffBlock::constant(lhs, rhs.blockPattern()) + rhs; } template - ForwardBlock operator+(const ForwardBlock& lhs, - const typename ForwardBlock::V& rhs) + AutoDiffBlock operator+(const AutoDiffBlock& lhs, + const typename AutoDiffBlock::V& rhs) { return rhs + lhs; // Commutative operation. } template - ForwardBlock operator-(const typename ForwardBlock::V& lhs, - const ForwardBlock& rhs) + AutoDiffBlock operator-(const typename AutoDiffBlock::V& lhs, + const AutoDiffBlock& rhs) { - return ForwardBlock::constant(lhs, rhs.blockPattern()) - rhs; + return AutoDiffBlock::constant(lhs, rhs.blockPattern()) - rhs; } template - ForwardBlock operator-(const ForwardBlock& lhs, - const typename ForwardBlock::V& rhs) + AutoDiffBlock operator-(const AutoDiffBlock& lhs, + const typename AutoDiffBlock::V& rhs) { - return lhs - ForwardBlock::constant(rhs, lhs.blockPattern()); + return lhs - AutoDiffBlock::constant(rhs, lhs.blockPattern()); } template - ForwardBlock operator/(const typename ForwardBlock::V& lhs, - const ForwardBlock& rhs) + AutoDiffBlock operator/(const typename AutoDiffBlock::V& lhs, + const AutoDiffBlock& rhs) { - return ForwardBlock::constant(lhs, rhs.blockPattern()) / rhs; + return AutoDiffBlock::constant(lhs, rhs.blockPattern()) / rhs; } template - ForwardBlock operator/(const ForwardBlock& lhs, - const typename ForwardBlock::V& rhs) + AutoDiffBlock operator/(const AutoDiffBlock& lhs, + const typename AutoDiffBlock::V& rhs) { - return lhs / ForwardBlock::constant(rhs, lhs.blockPattern()); + return lhs / AutoDiffBlock::constant(rhs, lhs.blockPattern()); } @@ -340,15 +340,15 @@ namespace AutoDiff * @return The product */ template - ForwardBlock operator*(const ForwardBlock& lhs, - const Scalar& rhs) + AutoDiffBlock operator*(const AutoDiffBlock& lhs, + const Scalar& rhs) { - std::vector< typename ForwardBlock::M > jac; + std::vector< typename AutoDiffBlock::M > jac; jac.reserve( lhs.numBlocks() ); for (int block=0; block::function( lhs.value() * rhs, jac ); + return AutoDiffBlock::function( lhs.value() * rhs, jac ); } @@ -360,14 +360,14 @@ namespace AutoDiff * @return The product */ template - ForwardBlock operator*(const Scalar& lhs, - const ForwardBlock& rhs) + AutoDiffBlock operator*(const Scalar& lhs, + const AutoDiffBlock& rhs) { return rhs * lhs; // Commutative operation. } -} // namespace Autodiff +} // namespace Opm diff --git a/opm/autodiff/AutoDiffHelpers.hpp b/opm/autodiff/AutoDiffHelpers.hpp index 30e4db08c..8959c2c79 100644 --- a/opm/autodiff/AutoDiffHelpers.hpp +++ b/opm/autodiff/AutoDiffHelpers.hpp @@ -26,14 +26,17 @@ #include +namespace Opm +{ + // -------------------- class HelperOps -------------------- /// Contains vectors and sparse matrices that represent subsets or /// operations on (AD or regular) vectors of data. struct HelperOps { - typedef AutoDiff::ForwardBlock::M M; - typedef AutoDiff::ForwardBlock::V V; + typedef AutoDiffBlock::M M; + typedef AutoDiffBlock::V V; /// A list of internal faces. typedef Eigen::Array IFaces; @@ -202,7 +205,7 @@ namespace { template class UpwindSelector { public: - typedef AutoDiff::ForwardBlock ADB; + typedef AutoDiffBlock ADB; UpwindSelector(const UnstructuredGrid& g, const HelperOps& h, @@ -297,11 +300,11 @@ namespace { /// Returns x(indices). template -AutoDiff::ForwardBlock -subset(const AutoDiff::ForwardBlock& x, +AutoDiffBlock +subset(const AutoDiffBlock& x, const IntVec& indices) { - return ::constructSubsetSparseMatrix(x.value().size(), indices) * x; + return constructSubsetSparseMatrix(x.value().size(), indices) * x; } @@ -312,19 +315,19 @@ Eigen::Array subset(const Eigen::Array& x, const IntVec& indices) { - return (::constructSubsetSparseMatrix(x.size(), indices) * x.matrix()).array(); + return (constructSubsetSparseMatrix(x.size(), indices) * x.matrix()).array(); } /// Returns v where v(indices) == x, v(!indices) == 0 and v.size() == n. template -AutoDiff::ForwardBlock -superset(const AutoDiff::ForwardBlock& x, +AutoDiffBlock +superset(const AutoDiffBlock& x, const IntVec& indices, const int n) { - return ::constructSupersetSparseMatrix(n, indices) * x; + return constructSupersetSparseMatrix(n, indices) * x; } @@ -336,7 +339,7 @@ superset(const Eigen::Array& x, const IntVec& indices, const int n) { - return ::constructSupersetSparseMatrix(n, indices) * x.matrix(); + return constructSupersetSparseMatrix(n, indices) * x.matrix(); } @@ -345,10 +348,10 @@ superset(const Eigen::Array& x, /// elements of d on the diagonal. /// Need to mark this as inline since it is defined in a header and not a template. inline -AutoDiff::ForwardBlock::M -spdiag(const AutoDiff::ForwardBlock::V& d) +AutoDiffBlock::M +spdiag(const AutoDiffBlock::V& d) { - typedef AutoDiff::ForwardBlock::M M; + typedef AutoDiffBlock::M M; const int n = d.size(); M mat(n, n); @@ -367,7 +370,7 @@ spdiag(const AutoDiff::ForwardBlock::V& d) template class Selector { public: - typedef AutoDiff::ForwardBlock ADB; + typedef AutoDiffBlock ADB; Selector(const typename ADB::V& selection_basis) { @@ -421,10 +424,10 @@ spdiag(const AutoDiff::ForwardBlock::V& d) /// Returns the input expression, but with all Jacobians collapsed to one. inline -AutoDiff::ForwardBlock -collapseJacs(const AutoDiff::ForwardBlock& x) +AutoDiffBlock +collapseJacs(const AutoDiffBlock& x) { - typedef AutoDiff::ForwardBlock ADB; + typedef AutoDiffBlock ADB; const int nb = x.numBlocks(); typedef Eigen::Triplet Tri; int nnz = 0; @@ -457,9 +460,9 @@ collapseJacs(const AutoDiff::ForwardBlock& x) /// Returns the vertical concatenation [ x; y ] of the inputs. inline -AutoDiff::ForwardBlock -vertcat(const AutoDiff::ForwardBlock& x, - const AutoDiff::ForwardBlock& y) +AutoDiffBlock +vertcat(const AutoDiffBlock& x, + const AutoDiffBlock& y) { const int nx = x.size(); const int ny = y.size(); @@ -585,4 +588,7 @@ inline Eigen::ArrayXd sign (const Eigen::ArrayXd& x) return retval; } + +} // namespace Opm + #endif // OPM_AUTODIFFHELPERS_HEADER_INCLUDED diff --git a/opm/autodiff/BlackoilPropsAd.hpp b/opm/autodiff/BlackoilPropsAd.hpp index 28696961a..23aaf2159 100644 --- a/opm/autodiff/BlackoilPropsAd.hpp +++ b/opm/autodiff/BlackoilPropsAd.hpp @@ -67,7 +67,7 @@ namespace Opm // Fluid interface // //////////////////////////// - typedef AutoDiff::ForwardBlock ADB; + typedef AutoDiffBlock ADB; typedef ADB::V V; typedef std::vector Cells; diff --git a/opm/autodiff/BlackoilPropsAdFromDeck.hpp b/opm/autodiff/BlackoilPropsAdFromDeck.hpp index c14e45788..09d07459e 100644 --- a/opm/autodiff/BlackoilPropsAdFromDeck.hpp +++ b/opm/autodiff/BlackoilPropsAdFromDeck.hpp @@ -75,7 +75,7 @@ namespace Opm // Fluid interface // //////////////////////////// - typedef AutoDiff::ForwardBlock ADB; + typedef AutoDiffBlock ADB; typedef ADB::V V; typedef std::vector Cells; diff --git a/opm/autodiff/BlackoilPropsAdInterface.hpp b/opm/autodiff/BlackoilPropsAdInterface.hpp index fb376d1a9..98e63fa77 100644 --- a/opm/autodiff/BlackoilPropsAdInterface.hpp +++ b/opm/autodiff/BlackoilPropsAdInterface.hpp @@ -64,7 +64,7 @@ namespace Opm // Fluid interface // //////////////////////////// - typedef AutoDiff::ForwardBlock ADB; + typedef AutoDiffBlock ADB; typedef ADB::V V; typedef ADB::M M; typedef std::vector Cells; diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.cpp b/opm/autodiff/FullyImplicitBlackoilSolver.cpp index 974657e90..42bdcdd77 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.cpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.cpp @@ -37,6 +37,7 @@ #include #include +// A debugging utility. #define DUMP(foo) \ do { \ std::cout << "==========================================\n" \ @@ -44,7 +45,11 @@ << collapseJacs(foo) << std::endl; \ } while (0) -typedef AutoDiff::ForwardBlock ADB; + + +namespace Opm { + +typedef AutoDiffBlock ADB; typedef ADB::V V; typedef ADB::M M; typedef Eigen::Array - AutoDiff::ForwardBlock::M + AutoDiffBlock::M gravityOperator(const UnstructuredGrid& grid, const HelperOps& ops , const GeoProps& geo ) @@ -86,8 +91,8 @@ namespace { } } - typedef AutoDiff::ForwardBlock::V V; - typedef AutoDiff::ForwardBlock::M M; + typedef AutoDiffBlock::V V; + typedef AutoDiffBlock::M M; const V& gpot = geo.gravityPotential(); const V& trans = geo.transmissibility(); @@ -182,8 +187,6 @@ namespace { -namespace Opm { - FullyImplicitBlackoilSolver:: FullyImplicitBlackoilSolver(const UnstructuredGrid& grid , diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.hpp b/opm/autodiff/FullyImplicitBlackoilSolver.hpp index a3859f472..a909745cf 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.hpp @@ -60,7 +60,7 @@ namespace Opm { private: // Types and enums - typedef AutoDiff::ForwardBlock ADB; + typedef AutoDiffBlock ADB; typedef ADB::V V; typedef ADB::M M; typedef Eigen::Array #include +namespace Opm { // Repeated from inside ImpesTPFAAD for convenience. -typedef AutoDiff::ForwardBlock ADB; +typedef AutoDiffBlock ADB; typedef ADB::V V; typedef ADB::M M; @@ -48,7 +49,7 @@ namespace { } template - AutoDiff::ForwardBlock::M + AutoDiffBlock::M gravityOperator(const UnstructuredGrid& grid, const HelperOps& ops , const GeoProps& geo ) @@ -65,8 +66,8 @@ namespace { } } - typedef AutoDiff::ForwardBlock::V V; - typedef AutoDiff::ForwardBlock::M M; + typedef AutoDiffBlock::V V; + typedef AutoDiffBlock::M M; const V& gpot = geo.gravityPotential(); const V& trans = geo.transmissibility(); @@ -122,7 +123,6 @@ namespace { } // anonymous namespace -namespace Opm { diff --git a/opm/autodiff/ImpesTPFAAD.hpp b/opm/autodiff/ImpesTPFAAD.hpp index 298f0543a..edb2480e3 100644 --- a/opm/autodiff/ImpesTPFAAD.hpp +++ b/opm/autodiff/ImpesTPFAAD.hpp @@ -63,7 +63,7 @@ namespace Opm { ImpesTPFAAD& operator=(const ImpesTPFAAD& rhs); // Types - typedef AutoDiff::ForwardBlock ADB; + typedef AutoDiffBlock ADB; typedef ADB::V V; typedef ADB::M M; typedef Eigen::Array ADB; + typedef AutoDiffBlock ADB; typedef ADB::V V; typedef ADB::M M; diff --git a/tests/test_block.cpp b/tests/test_block.cpp index c6cb04154..6a8936e3c 100644 --- a/tests/test_block.cpp +++ b/tests/test_block.cpp @@ -32,6 +32,8 @@ #include #include +using namespace Opm; + namespace { template bool @@ -73,7 +75,7 @@ namespace { BOOST_AUTO_TEST_CASE(ConstantInitialisation) { - typedef AutoDiff::ForwardBlock ADB; + typedef AutoDiffBlock ADB; std::vector blocksizes = { 3, 1, 2 }; @@ -92,7 +94,7 @@ BOOST_AUTO_TEST_CASE(ConstantInitialisation) BOOST_AUTO_TEST_CASE(VariableInitialisation) { - typedef AutoDiff::ForwardBlock ADB; + typedef AutoDiffBlock ADB; std::vector blocksizes = { 3, 1, 2 }; @@ -119,7 +121,7 @@ BOOST_AUTO_TEST_CASE(VariableInitialisation) BOOST_AUTO_TEST_CASE(FunctionInitialisation) { - typedef AutoDiff::ForwardBlock ADB; + typedef AutoDiffBlock ADB; std::vector blocksizes = { 3, 1, 2 }; std::vector::size_type num_blocks = blocksizes.size(); @@ -150,7 +152,7 @@ BOOST_AUTO_TEST_CASE(FunctionInitialisation) BOOST_AUTO_TEST_CASE(Addition) { - typedef AutoDiff::ForwardBlock ADB; + typedef AutoDiffBlock ADB; std::vector blocksizes = { 3, 1, 2 }; ADB::V va(3); @@ -195,7 +197,7 @@ BOOST_AUTO_TEST_CASE(Addition) int main() try { - typedef AutoDiff::ForwardBlock ADB; + typedef AutoDiffBlock ADB; std::vector blocksizes = { 3, 1, 2 }; int num_blocks = blocksizes.size(); ADB::V v1(3); diff --git a/tests/test_scalar_mult.cpp b/tests/test_scalar_mult.cpp index ac2077dce..aa2741aac 100644 --- a/tests/test_scalar_mult.cpp +++ b/tests/test_scalar_mult.cpp @@ -34,9 +34,11 @@ #include +using namespace Opm; + BOOST_AUTO_TEST_CASE(ScalarMultiplication) { - typedef AutoDiff::ForwardBlock ADB; + typedef AutoDiffBlock ADB; std::vector blocksizes = { 3, 1, 2 }; ADB::V vx(3); diff --git a/tests/test_span.cpp b/tests/test_span.cpp index f4d0d8e18..5ddee5695 100644 --- a/tests/test_span.cpp +++ b/tests/test_span.cpp @@ -32,6 +32,8 @@ #include +using namespace Opm; + BOOST_AUTO_TEST_CASE(OneArgConstr) { const int num = 4; From 9a20c1ee02df62001e57165cfc92a077796b36d7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 19 Sep 2013 14:07:05 +0200 Subject: [PATCH 06/14] Documented AutoDiffBlock. --- opm/autodiff/AutoDiffBlock.hpp | 98 ++++++++++++++++++++++++++++++---- 1 file changed, 87 insertions(+), 11 deletions(-) diff --git a/opm/autodiff/AutoDiffBlock.hpp b/opm/autodiff/AutoDiffBlock.hpp index ef73bd5a2..161fa16ca 100644 --- a/opm/autodiff/AutoDiffBlock.hpp +++ b/opm/autodiff/AutoDiffBlock.hpp @@ -28,15 +28,66 @@ namespace Opm { + /// A class for forward-mode automatic differentiation with vector + /// values and sparse jacobian matrices. + /// + /// The class contains a (column) vector of values and multiple + /// sparse matrices representing its derivatives. Each such matrix + /// has a number of rows equal to the number of rows in the value + /// vector, and a number of columns equal to the number of + /// variables we want to compute the derivatives with respect + /// to. The reason to have multiple such jacobians instead of just + /// one is to allow simpler grouping of variables, making it + /// easier to implement various preconditioning schemes. Only + /// basic arithmetic operators are implemented for this class, + /// reflecting our needs so far. + /// + /// The class is built on the Eigen library, using an Eigen array + /// type to contain the values and Eigen sparse matrices for the + /// jacobians. The overloaded operators are intended to behave in + /// a similar way to Eigen arrays, meaning that the '*' operator + /// is elementwise multiplication. The only exception is + /// multiplication with a sparse matrix from the left, which is + /// treated as an Eigen matrix operation. + /// + /// There are no public constructors, instead we use the Named + /// Constructor pattern. In general, one needs to know which + /// variables one wants to compute the derivatives with respect to + /// before constructing an AutoDiffBlock. Some of the constructors + /// require you to pass a block pattern. This should be a vector + /// containing the number of columns you want for each jacobian + /// matrix. + /// + /// For example: you want the derivatives with respect to three + /// different variables p, r and s. Assuming that there are 10 + /// elements in p, and 20 in each of r and s, the block pattern is + /// { 10, 20, 20 }. When creating the variables p, r and s in your + /// program you have two options: + /// a) Use the variable() constructor three times, passing the + /// index (0 for p, 1 for r and 2 for s), initial value of + /// each variable and the block pattern. + /// b) Use the variables() constructor passing only the initial + /// values of each variable. The block pattern will be + /// inferred from the size of the initial value vectors. + /// This is usually the simplest option if you have multiple + /// variables. Note that this constructor returns a vector + /// of all three variables, so you need to use index access + /// (operator[]) to get the individual variables (that is p, + /// r and s). + /// After this, the r variable for example will have a size() of + /// 20 and three jacobian matrices. The first is a 20 by 10 zero + /// matrix, the second is a 20 by 20 identity matrix, and the + /// third is a 20 by 20 zero matrix. template class AutoDiffBlock { public: - /// Underlying types for scalar vectors and jacobians. + /// Underlying type for values. typedef Eigen::Array V; + /// Underlying type for jacobians. typedef Eigen::SparseMatrix M; - /// Named constructor pattern used here. + /// Construct an empty AutoDiffBlock. static AutoDiffBlock null() { V val; @@ -44,6 +95,9 @@ namespace Opm return AutoDiffBlock(val, jac); } + /// Create an AutoDiffBlock representing a constant. + /// \param[in] val values + /// \param[in] blocksizes block pattern static AutoDiffBlock constant(const V& val, const std::vector& blocksizes) { std::vector jac; @@ -57,6 +111,13 @@ namespace Opm return AutoDiffBlock(val, jac); } + /// Create an AutoDiffBlock representing a single variable block. + /// \param[in] index index of the variable you are constructing + /// \param[in] val values + /// \param[in] blocksizes block pattern + /// The resulting object will have size() equal to block_pattern[index]. + /// Its jacobians will all be zero, except for derivative()[index], which + /// will be an identity matrix. static AutoDiffBlock variable(const int index, const V& val, const std::vector& blocksizes) { std::vector jac; @@ -76,13 +137,16 @@ namespace Opm return AutoDiffBlock(val, jac); } + /// Create an AutoDiffBlock by directly specifying values and jacobians. + /// \param[in] val values + /// \param[in] jac vector of jacobians static AutoDiffBlock function(const V& val, const std::vector& jac) { return AutoDiffBlock(val, jac); } - /// Construct a set of primary variables, - /// each initialized to a given vector. + /// Construct a set of primary variables, each initialized to + /// a given vector. static std::vector variables(const std::vector& initial_values) { const int num_vars = initial_values.size(); @@ -97,7 +161,7 @@ namespace Opm return vars; } - /// Operator += + /// Elementwise operator += AutoDiffBlock& operator+=(const AutoDiffBlock& rhs) { assert (numBlocks() == rhs.numBlocks()); @@ -115,7 +179,7 @@ namespace Opm return *this; } - /// Operator + + /// Elementwise operator + AutoDiffBlock operator+(const AutoDiffBlock& rhs) const { std::vector jac = jac_; @@ -129,7 +193,7 @@ namespace Opm return function(val_ + rhs.val_, jac); } - /// Operator - + /// Elementwise operator - AutoDiffBlock operator-(const AutoDiffBlock& rhs) const { std::vector jac = jac_; @@ -143,7 +207,7 @@ namespace Opm return function(val_ - rhs.val_, jac); } - /// Operator * + /// Elementwise operator * AutoDiffBlock operator*(const AutoDiffBlock& rhs) const { int num_blocks = numBlocks(); @@ -160,7 +224,7 @@ namespace Opm return function(val_ * rhs.val_, jac); } - /// Operator / + /// Elementwise operator / AutoDiffBlock operator/(const AutoDiffBlock& rhs) const { int num_blocks = numBlocks(); @@ -177,6 +241,7 @@ namespace Opm } return function(val_ / rhs.val_, jac); } + /// I/O. template Ostream& @@ -213,13 +278,13 @@ namespace Opm return bp; } - /// Function value + /// Function value. const V& value() const { return val_; } - /// Function derivatives + /// Function derivatives. const std::vector& derivative() const { return jac_; @@ -244,6 +309,9 @@ namespace Opm }; + // --------- Free functions and operators for AutoDiffBlock --------- + + /// Stream output. template Ostream& operator<<(Ostream& os, const AutoDiffBlock& fw) @@ -268,6 +336,7 @@ namespace Opm } + /// Elementwise multiplication with constant on the left. template AutoDiffBlock operator*(const typename AutoDiffBlock::V& lhs, const AutoDiffBlock& rhs) @@ -276,6 +345,7 @@ namespace Opm } + /// Elementwise multiplication with constant on the right. template AutoDiffBlock operator*(const AutoDiffBlock& lhs, const typename AutoDiffBlock::V& rhs) @@ -284,6 +354,7 @@ namespace Opm } + /// Elementwise addition with constant on the left. template AutoDiffBlock operator+(const typename AutoDiffBlock::V& lhs, const AutoDiffBlock& rhs) @@ -292,6 +363,7 @@ namespace Opm } + /// Elementwise addition with constant on the right. template AutoDiffBlock operator+(const AutoDiffBlock& lhs, const typename AutoDiffBlock::V& rhs) @@ -300,6 +372,7 @@ namespace Opm } + /// Elementwise subtraction with constant on the left. template AutoDiffBlock operator-(const typename AutoDiffBlock::V& lhs, const AutoDiffBlock& rhs) @@ -308,6 +381,7 @@ namespace Opm } + /// Elementwise subtraction with constant on the right. template AutoDiffBlock operator-(const AutoDiffBlock& lhs, const typename AutoDiffBlock::V& rhs) @@ -316,6 +390,7 @@ namespace Opm } + /// Elementwise division with constant on the left. template AutoDiffBlock operator/(const typename AutoDiffBlock::V& lhs, const AutoDiffBlock& rhs) @@ -324,6 +399,7 @@ namespace Opm } + /// Elementwise division with constant on the right. template AutoDiffBlock operator/(const AutoDiffBlock& lhs, const typename AutoDiffBlock::V& rhs) From 0aa96af32921d0cbc343c40443f5a4dda33f95c3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 19 Sep 2013 14:45:40 +0200 Subject: [PATCH 07/14] Documentation refinement. Added main doc file. --- doc/doxygen/Doxylocal | 25 +++++++++++ opm/autodiff/AutoDiffBlock.hpp | 9 ++-- opm/autodiff/opm-autodiff_doxygen_main.hpp | 49 ++++++++++++++++++++++ 3 files changed, 79 insertions(+), 4 deletions(-) create mode 100644 doc/doxygen/Doxylocal create mode 100644 opm/autodiff/opm-autodiff_doxygen_main.hpp diff --git a/doc/doxygen/Doxylocal b/doc/doxygen/Doxylocal new file mode 100644 index 000000000..db859f51a --- /dev/null +++ b/doc/doxygen/Doxylocal @@ -0,0 +1,25 @@ +# This file contains local changes to the doxygen configuration +# please us '+=' to add file/directories to the lists + +# The INPUT tag can be used to specify the files and/or directories that contain +# documented source files. You may enter file names like "myfile.cpp" or +# directories like "/usr/src/myproject". Separate the files or directories +# with spaces. + +INPUT += @abs_top_srcdir@/opm/ + +# see e.g. dune-grid for the examples of mainpage and modules +# INPUT += @abs_top_srcdir@/doc/doxygen/mainpage # @srcdir@/modules + +# The EXCLUDE tag can be used to specify files and/or directories that should +# excluded from the INPUT source files. This way you can easily exclude a +# subdirectory from a directory tree whose root is specified with the INPUT tag. + +# EXCLUDE += @abs_top_srcdir@/tests/ + +# The EXAMPLE_PATH tag can be used to specify one or more files or +# directories that contain example code fragments that are included (see +# the \include command). + +# EXAMPLE_PATH += @abs_top_srcdir@/examples/ + diff --git a/opm/autodiff/AutoDiffBlock.hpp b/opm/autodiff/AutoDiffBlock.hpp index 161fa16ca..56fbfd451 100644 --- a/opm/autodiff/AutoDiffBlock.hpp +++ b/opm/autodiff/AutoDiffBlock.hpp @@ -45,8 +45,8 @@ namespace Opm /// The class is built on the Eigen library, using an Eigen array /// type to contain the values and Eigen sparse matrices for the /// jacobians. The overloaded operators are intended to behave in - /// a similar way to Eigen arrays, meaning that the '*' operator - /// is elementwise multiplication. The only exception is + /// a similar way to Eigen arrays, meaning for example that the * + /// operator is elementwise multiplication. The only exception is /// multiplication with a sparse matrix from the left, which is /// treated as an Eigen matrix operation. /// @@ -63,10 +63,10 @@ namespace Opm /// elements in p, and 20 in each of r and s, the block pattern is /// { 10, 20, 20 }. When creating the variables p, r and s in your /// program you have two options: - /// a) Use the variable() constructor three times, passing the + /// - Use the variable() constructor three times, passing the /// index (0 for p, 1 for r and 2 for s), initial value of /// each variable and the block pattern. - /// b) Use the variables() constructor passing only the initial + /// - Use the variables() constructor passing only the initial /// values of each variable. The block pattern will be /// inferred from the size of the initial value vectors. /// This is usually the simplest option if you have multiple @@ -74,6 +74,7 @@ namespace Opm /// of all three variables, so you need to use index access /// (operator[]) to get the individual variables (that is p, /// r and s). + /// /// After this, the r variable for example will have a size() of /// 20 and three jacobian matrices. The first is a 20 by 10 zero /// matrix, the second is a 20 by 20 identity matrix, and the diff --git a/opm/autodiff/opm-autodiff_doxygen_main.hpp b/opm/autodiff/opm-autodiff_doxygen_main.hpp new file mode 100644 index 000000000..80391d6c5 --- /dev/null +++ b/opm/autodiff/opm-autodiff_doxygen_main.hpp @@ -0,0 +1,49 @@ +/* + Copyright 2013 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_OPM-AUTODIFF_DOXYGEN_MAIN_HEADER_INCLUDED +#define OPM_OPM-AUTODIFF_DOXYGEN_MAIN_HEADER_INCLUDED + + +/** \mainpage Documentation for the opm-autodiff library. + +

Automatic differentiation

+ +This library implements automatic differentiation for vector data with +multiple blocks of sparse jacobians. This is contained in the class +Opm::AutoDiffBlock. Also available is Opm::AutoDiff, a much simpler +single-value single-derivative AD class. + +There are also some helper classes and functions that are intended to +aid in the development of solvers and simulators with AD, these +include Opm::HelperOps, Opm::UpwindSelector, Opm::subset, +Opm::superset, Opm::Selector, Opm::collapseJacs, Opm::vertcat, +Opm::Span and Opm::sign. + +

Solvers and simulators

+ +There are some solvers and simulators in opm-autodiff. They should all +be considered experimental prototypes at this point. Notable simulator +prototypes include +- examples/sim_fibo_ad.cpp, a fully implicit black-oil simulator. +- examples/sim_2p_incomp_ad.cpp, a sequential incompressible 2-phase simulator. + +*/ + +#endif // OPM_OPM-AUTODIFF_DOXYGEN_MAIN_HEADER_INCLUDED From 9337a6227b2e8859f8df7e117cfcafc4b15d348a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 19 Sep 2013 14:52:39 +0200 Subject: [PATCH 08/14] Renamed misleading adfi -> fi. Some classes and a program were renamed since they were not fully implicit solvers despite their names indicating it. --- CMakeLists_files.cmake | 8 ++-- ...p_incomp_adfi.cpp => sim_2p_incomp_ad.cpp} | 38 +++++++++---------- ...Adfi.cpp => SimulatorIncompTwophaseAd.cpp} | 12 +++--- ...Adfi.hpp => SimulatorIncompTwophaseAd.hpp} | 24 ++++++------ 4 files changed, 41 insertions(+), 41 deletions(-) rename examples/{sim_2p_incomp_adfi.cpp => sim_2p_incomp_ad.cpp} (90%) rename opm/autodiff/{SimulatorIncompTwophaseAdfi.cpp => SimulatorIncompTwophaseAd.cpp} (98%) rename opm/autodiff/{SimulatorIncompTwophaseAdfi.hpp => SimulatorIncompTwophaseAd.hpp} (83%) diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index bb86b3c49..a7b7a1b1a 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -32,7 +32,7 @@ list (APPEND MAIN_SOURCE_FILES opm/autodiff/ImpesTPFAAD.cpp opm/autodiff/SimulatorCompressibleAd.cpp opm/autodiff/SimulatorFullyImplicitBlackoil.cpp - opm/autodiff/SimulatorIncompTwophaseAdfi.cpp + opm/autodiff/SimulatorIncompTwophaseAd.cpp opm/autodiff/TransportSolverTwophaseAd.cpp opm/autodiff/BlackoilPropsAdFromDeck.cpp ) @@ -57,7 +57,7 @@ list (APPEND EXAMPLE_SOURCE_FILES examples/find_zero.cpp examples/sim_fibo_ad.cpp examples/sim_2p_comp_ad.cpp - examples/sim_2p_incomp_adfi.cpp + examples/sim_2p_incomp_ad.cpp examples/sim_simple.cpp examples/test_impestpfa_ad.cpp examples/test_implicit_ad.cpp @@ -66,7 +66,7 @@ list (APPEND EXAMPLE_SOURCE_FILES # programs listed here will not only be compiled, but also marked for # installation list (APPEND PROGRAM_SOURCE_FILES - examples/sim_2p_incomp_adfi.cpp + examples/sim_2p_incomp_ad.cpp examples/sim_fibo_ad.cpp ) @@ -84,6 +84,6 @@ list (APPEND PUBLIC_HEADER_FILES opm/autodiff/FullyImplicitBlackoilSolver.hpp opm/autodiff/SimulatorCompressibleAd.hpp opm/autodiff/SimulatorFullyImplicitBlackoil.hpp - opm/autodiff/SimulatorIncompTwophaseAdfi.hpp + opm/autodiff/SimulatorIncompTwophaseAd.hpp opm/autodiff/TransportSolverTwophaseAd.hpp ) diff --git a/examples/sim_2p_incomp_adfi.cpp b/examples/sim_2p_incomp_ad.cpp similarity index 90% rename from examples/sim_2p_incomp_adfi.cpp rename to examples/sim_2p_incomp_ad.cpp index 05e6f8bf9..7469bdca0 100644 --- a/examples/sim_2p_incomp_adfi.cpp +++ b/examples/sim_2p_incomp_ad.cpp @@ -44,7 +44,7 @@ #include #include -#include +#include #include #include @@ -216,15 +216,15 @@ try if (!use_deck) { // Simple simulation without a deck. WellsManager wells; // no wells. - SimulatorIncompTwophaseAdfi simulator(param, - *grid->c_grid(), - *props, - rock_comp->isActive() ? rock_comp.get() : 0, - wells, - src, - bcs.c_bcs(), - linsolver, - grav); + SimulatorIncompTwophaseAd simulator(param, + *grid->c_grid(), + *props, + rock_comp->isActive() ? rock_comp.get() : 0, + wells, + src, + bcs.c_bcs(), + linsolver, + grav); SimulatorTimer simtimer; simtimer.init(param); warnIfUnusedParams(param); @@ -271,15 +271,15 @@ try } // Create and run simulator. - SimulatorIncompTwophaseAdfi simulator(param, - *grid->c_grid(), - *props, - rock_comp->isActive() ? rock_comp.get() : 0, - wells, - src, - bcs.c_bcs(), - linsolver, - grav); + SimulatorIncompTwophaseAd simulator(param, + *grid->c_grid(), + *props, + rock_comp->isActive() ? rock_comp.get() : 0, + wells, + src, + bcs.c_bcs(), + linsolver, + grav); if (epoch == 0) { warnIfUnusedParams(param); } diff --git a/opm/autodiff/SimulatorIncompTwophaseAdfi.cpp b/opm/autodiff/SimulatorIncompTwophaseAd.cpp similarity index 98% rename from opm/autodiff/SimulatorIncompTwophaseAdfi.cpp rename to opm/autodiff/SimulatorIncompTwophaseAd.cpp index 5ad933fd7..d009927f6 100644 --- a/opm/autodiff/SimulatorIncompTwophaseAdfi.cpp +++ b/opm/autodiff/SimulatorIncompTwophaseAd.cpp @@ -22,7 +22,7 @@ #include "config.h" #endif // HAVE_CONFIG_H -#include +#include #include #include @@ -60,7 +60,7 @@ namespace Opm { - class SimulatorIncompTwophaseAdfi::Impl + class SimulatorIncompTwophaseAd::Impl { public: Impl(const parameter::ParameterGroup& param, @@ -109,7 +109,7 @@ namespace Opm - SimulatorIncompTwophaseAdfi::SimulatorIncompTwophaseAdfi(const parameter::ParameterGroup& param, + SimulatorIncompTwophaseAd::SimulatorIncompTwophaseAd(const parameter::ParameterGroup& param, const UnstructuredGrid& grid, const IncompPropertiesInterface& props, const RockCompressibility* rock_comp_props, @@ -126,7 +126,7 @@ namespace Opm - SimulatorReport SimulatorIncompTwophaseAdfi::run(SimulatorTimer& timer, + SimulatorReport SimulatorIncompTwophaseAd::run(SimulatorTimer& timer, TwophaseState& state, WellState& well_state) { @@ -310,7 +310,7 @@ namespace Opm - SimulatorIncompTwophaseAdfi::Impl::Impl(const parameter::ParameterGroup& param, + SimulatorIncompTwophaseAd::Impl::Impl(const parameter::ParameterGroup& param, const UnstructuredGrid& grid, const IncompPropertiesInterface& props, const RockCompressibility* rock_comp_props, @@ -409,7 +409,7 @@ namespace Opm - SimulatorReport SimulatorIncompTwophaseAdfi::Impl::run(SimulatorTimer& timer, + SimulatorReport SimulatorIncompTwophaseAd::Impl::run(SimulatorTimer& timer, TwophaseState& state, WellState& well_state) { diff --git a/opm/autodiff/SimulatorIncompTwophaseAdfi.hpp b/opm/autodiff/SimulatorIncompTwophaseAd.hpp similarity index 83% rename from opm/autodiff/SimulatorIncompTwophaseAdfi.hpp rename to opm/autodiff/SimulatorIncompTwophaseAd.hpp index e035a0e6a..70545e77e 100644 --- a/opm/autodiff/SimulatorIncompTwophaseAdfi.hpp +++ b/opm/autodiff/SimulatorIncompTwophaseAd.hpp @@ -17,8 +17,8 @@ along with OPM. If not, see . */ -#ifndef OPM_SIMULATORINCOMPTWOPHASEADFI_HEADER_INCLUDED -#define OPM_SIMULATORINCOMPTWOPHASEADFI_HEADER_INCLUDED +#ifndef OPM_SIMULATORINCOMPTWOPHASEAD_HEADER_INCLUDED +#define OPM_SIMULATORINCOMPTWOPHASEAD_HEADER_INCLUDED #include #include @@ -40,7 +40,7 @@ namespace Opm struct SimulatorReport; /// Class collecting all necessary components for a two-phase simulation. - class SimulatorIncompTwophaseAdfi + class SimulatorIncompTwophaseAd { public: /// Initialise from parameters and objects to observe. @@ -67,14 +67,14 @@ namespace Opm /// \param[in] bcs boundary conditions, treat as all noflow if null /// \param[in] linsolver linear solver /// \param[in] gravity if non-null, gravity vector - SimulatorIncompTwophaseAdfi(const parameter::ParameterGroup& param, - const UnstructuredGrid& grid, - const IncompPropertiesInterface& props, - const RockCompressibility* rock_comp_props, - WellsManager& wells_manager, - const std::vector& src, - const FlowBoundaryConditions* bcs, - LinearSolverInterface& linsolver, + SimulatorIncompTwophaseAd(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const IncompPropertiesInterface& props, + const RockCompressibility* rock_comp_props, + WellsManager& wells_manager, + const std::vector& src, + const FlowBoundaryConditions* bcs, + LinearSolverInterface& linsolver, const double* gravity); /// Run the simulation. @@ -96,4 +96,4 @@ namespace Opm } // namespace Opm -#endif // OPM_SIMULATORINCOMPTWOPHASEADFI_HEADER_INCLUDED +#endif // OPM_SIMULATORINCOMPTWOPHASEAD_HEADER_INCLUDED From ad78cb87139263be256a8af166fc5d710f10c07c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 20 Sep 2013 14:31:57 +0200 Subject: [PATCH 09/14] Document class BlackoilPropsAd. --- opm/autodiff/BlackoilPropsAd.hpp | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/opm/autodiff/BlackoilPropsAd.hpp b/opm/autodiff/BlackoilPropsAd.hpp index 23aaf2159..b00a00f29 100644 --- a/opm/autodiff/BlackoilPropsAd.hpp +++ b/opm/autodiff/BlackoilPropsAd.hpp @@ -29,9 +29,17 @@ namespace Opm class BlackoilPropertiesInterface; - /// This class is intended to present a fluid interface for - /// three-phase black-oil that is easy to use with the AD-using - /// simulators. + /// This class implements the AD-adapted fluid interface for + /// three-phase black-oil. + /// + /// It is implemented by wrapping a BlackoilPropertiesInterface + /// object (the interface class defined in opm-core) and calling + /// its methods. This approach works well for most methods, but + /// the rsMax() method cannot be implemented by such a wrapping, + /// without access to the underlying pvt objects. Therefore we + /// cannot use this class with any case that involves + /// miscibility. A rethinking of fluid interfaces is probably + /// necessary. /// /// Most methods are available in two overloaded versions, one /// taking a constant vector and returning the same, and one From a10eb61b665597da062589f6e06bc3cdb69879bb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 20 Sep 2013 14:32:29 +0200 Subject: [PATCH 10/14] Require deck and remove simple source term in fully implicit sim. --- examples/sim_fibo_ad.cpp | 223 +++++++----------- .../SimulatorFullyImplicitBlackoil.cpp | 13 +- .../SimulatorFullyImplicitBlackoil.hpp | 2 - 3 files changed, 87 insertions(+), 151 deletions(-) diff --git a/examples/sim_fibo_ad.cpp b/examples/sim_fibo_ad.cpp index b5f42ef25..33818da6d 100644 --- a/examples/sim_fibo_ad.cpp +++ b/examples/sim_fibo_ad.cpp @@ -45,7 +45,6 @@ #include #include -#include #include #include @@ -78,12 +77,16 @@ try { using namespace Opm; - std::cout << "\n================ Test program for weakly compressible two-phase flow ===============\n\n"; + std::cout << "\n================ Test program for fully implicit three-phase black-oil flow ===============\n\n"; 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"); + if (!use_deck) { + OPM_THROW(std::runtime_error, "This program must be run with an input deck. " + "Specify the deck with deck_filename=deckname.data (for example)."); + } boost::scoped_ptr deck; boost::scoped_ptr grid; boost::scoped_ptr props; @@ -93,83 +96,40 @@ try // 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"); - deck.reset(new EclipseGridParser(deck_filename)); - // Grid init - grid.reset(new GridManager(*deck)); - // Rock and fluid init - props.reset(new BlackoilPropertiesFromDeck(*deck, *grid->c_grid(), param)); - new_props.reset(new BlackoilPropsAdFromDeck(*deck, *grid->c_grid())); - // check_well_controls = param.getDefault("check_well_controls", false); - // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10); - // Rock compressibility. - rock_comp.reset(new RockCompressibility(*deck)); - // Gravity. - gravity[2] = deck->hasField("NOGRAV") ? 0.0 : unit::gravity; - // Init state variables (saturation and pressure). - if (param.has("init_saturation")) { - initStateBasic(*grid->c_grid(), *props, param, gravity[2], state); - initBlackoilSurfvol(*grid->c_grid(), *props, state); - enum { Oil = BlackoilPhases::Liquid, Gas = BlackoilPhases::Vapour }; - const PhaseUsage pu = props->phaseUsage(); - if (pu.phase_used[Oil] && pu.phase_used[Gas]) { - const int np = props->numPhases(); - const int nc = grid->c_grid()->number_of_cells; - for (int c = 0; c < nc; ++c) { - state.gasoilratio()[c] = state.surfacevol()[c*np + pu.phase_pos[Gas]] - / state.surfacevol()[c*np + pu.phase_pos[Oil]]; - } - } - } else { - initBlackoilStateFromDeck(*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 GridManager(nx, ny, nz, dx, dy, dz)); - // Rock and fluid init. - props.reset(new BlackoilPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells)); - new_props.reset(new BlackoilPropsAd(*props)); - // Rock compressibility. - rock_comp.reset(new RockCompressibility(param)); - // Gravity. - gravity[2] = param.getDefault("gravity", 0.0); - // Init state variables (saturation and pressure). + std::string deck_filename = param.get("deck_filename"); + deck.reset(new EclipseGridParser(deck_filename)); + // Grid init + grid.reset(new GridManager(*deck)); + // Rock and fluid init + props.reset(new BlackoilPropertiesFromDeck(*deck, *grid->c_grid(), param)); + new_props.reset(new BlackoilPropsAdFromDeck(*deck, *grid->c_grid())); + // check_well_controls = param.getDefault("check_well_controls", false); + // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10); + // Rock compressibility. + rock_comp.reset(new RockCompressibility(*deck)); + // Gravity. + gravity[2] = deck->hasField("NOGRAV") ? 0.0 : unit::gravity; + // Init state variables (saturation and pressure). + if (param.has("init_saturation")) { initStateBasic(*grid->c_grid(), *props, param, gravity[2], state); initBlackoilSurfvol(*grid->c_grid(), *props, state); + enum { Oil = BlackoilPhases::Liquid, Gas = BlackoilPhases::Vapour }; + const PhaseUsage pu = props->phaseUsage(); + if (pu.phase_used[Oil] && pu.phase_used[Gas]) { + const int np = props->numPhases(); + const int nc = grid->c_grid()->number_of_cells; + for (int c = 0; c < nc; ++c) { + state.gasoilratio()[c] = state.surfacevol()[c*np + pu.phase_pos[Gas]] + / state.surfacevol()[c*np + pu.phase_pos[Oil]]; + } + } + } else { + initBlackoilStateFromDeck(*grid->c_grid(), *props, *deck, gravity[2], state); } bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0); const double *grav = use_gravity ? &gravity[0] : 0; - // Initialising src - int num_cells = grid->c_grid()->number_of_cells; - std::vector src(num_cells, 0.0); - if (use_deck) { - // Do nothing, wells will be the driving force, not source terms. - } else { - // Compute pore volumes, in order to enable specifying injection rate - // terms of total pore volume. - 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); - } - const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0); - 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/unit::day; - src[0] = flow_per_sec; - src[num_cells - 1] = -flow_per_sec; - } - // Boundary conditions. FlowBCManager bcs; if (param.getDefault("use_pside", false)) { @@ -207,87 +167,66 @@ try std::cout << "\n\n================ Starting main simulation loop ===============\n" << " (number of epochs: " - << (use_deck ? deck->numberOfEpochs() : 1) << ")\n\n" << std::flush; + << (deck->numberOfEpochs()) << ")\n\n" << std::flush; SimulatorReport rep; - if (!use_deck) { - // Simple simulation without a deck. - WellsManager wells; // no wells. + // With a deck, we may have more epochs etc. + WellState well_state; + int step = 0; + SimulatorTimer simtimer; + // Use timer for last epoch to obtain total time. + deck->setCurrentEpoch(deck->numberOfEpochs() - 1); + simtimer.init(*deck); + const double total_time = simtimer.totalTime(); + for (int epoch = 0; epoch < deck->numberOfEpochs(); ++epoch) { + // Set epoch index. + deck->setCurrentEpoch(epoch); + + // Update the timer. + if (deck->hasField("TSTEP")) { + simtimer.init(*deck); + } else { + if (epoch != 0) { + OPM_THROW(std::runtime_error, "No TSTEP in deck for epoch " << epoch); + } + simtimer.init(param); + } + simtimer.setCurrentStepNum(step); + simtimer.setTotalTime(total_time); + + // Report on start of epoch. + std::cout << "\n\n-------------- Starting epoch " << epoch << " --------------" + << "\n (number of steps: " + << simtimer.numSteps() - step << ")\n\n" << std::flush; + + // Create new wells, well_state + WellsManager wells(*deck, *grid->c_grid(), props->permeability()); + // @@@ HACK: we should really make a new well state and + // properly transfer old well state to it every epoch, + // since number of wells may change etc. + if (epoch == 0) { + well_state.init(wells.c_wells(), state); + } + + // Create and run simulator. SimulatorFullyImplicitBlackoil simulator(param, *grid->c_grid(), *new_props, rock_comp->isActive() ? rock_comp.get() : 0, wells, - src, bcs.c_bcs(), linsolver, grav); - SimulatorTimer simtimer; - simtimer.init(param); - warnIfUnusedParams(param); - WellState well_state; - well_state.init(0, state); - rep = simulator.run(simtimer, state, well_state); - } else { - // With a deck, we may have more epochs etc. - WellState well_state; - int step = 0; - SimulatorTimer simtimer; - // Use timer for last epoch to obtain total time. - deck->setCurrentEpoch(deck->numberOfEpochs() - 1); - simtimer.init(*deck); - const double total_time = simtimer.totalTime(); - for (int epoch = 0; epoch < deck->numberOfEpochs(); ++epoch) { - // Set epoch index. - deck->setCurrentEpoch(epoch); - - // Update the timer. - if (deck->hasField("TSTEP")) { - simtimer.init(*deck); - } else { - if (epoch != 0) { - OPM_THROW(std::runtime_error, "No TSTEP in deck for epoch " << epoch); - } - simtimer.init(param); - } - simtimer.setCurrentStepNum(step); - simtimer.setTotalTime(total_time); - - // Report on start of epoch. - std::cout << "\n\n-------------- Starting epoch " << epoch << " --------------" - << "\n (number of steps: " - << simtimer.numSteps() - step << ")\n\n" << std::flush; - - // Create new wells, well_state - WellsManager wells(*deck, *grid->c_grid(), props->permeability()); - // @@@ HACK: we should really make a new well state and - // properly transfer old well state to it every epoch, - // since number of wells may change etc. - if (epoch == 0) { - well_state.init(wells.c_wells(), state); - } - - // Create and run simulator. - SimulatorFullyImplicitBlackoil simulator(param, - *grid->c_grid(), - *new_props, - rock_comp->isActive() ? rock_comp.get() : 0, - wells, - src, - bcs.c_bcs(), - linsolver, - grav); - if (epoch == 0) { - warnIfUnusedParams(param); - } - SimulatorReport epoch_rep = simulator.run(simtimer, state, well_state); - if (output) { - epoch_rep.reportParam(epoch_os); - } - // Update total timing report and remember step number. - rep += epoch_rep; - step = simtimer.currentStepNum(); + if (epoch == 0) { + warnIfUnusedParams(param); } + SimulatorReport epoch_rep = simulator.run(simtimer, state, well_state); + if (output) { + epoch_rep.reportParam(epoch_os); + } + // Update total timing report and remember step number. + rep += epoch_rep; + step = simtimer.currentStepNum(); } std::cout << "\n\n================ End of simulation ===============\n\n"; diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoil.cpp b/opm/autodiff/SimulatorFullyImplicitBlackoil.cpp index 6b1d6f931..5b1b6718f 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoil.cpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoil.cpp @@ -70,7 +70,6 @@ namespace Opm const BlackoilPropsAdInterface& props, const RockCompressibility* rock_comp_props, WellsManager& wells_manager, - const std::vector& src, const FlowBoundaryConditions* bcs, LinearSolverInterface& linsolver, const double* gravity); @@ -99,7 +98,6 @@ namespace Opm const RockCompressibility* rock_comp_props_; WellsManager& wells_manager_; const Wells* wells_; - const std::vector& src_; const FlowBoundaryConditions* bcs_; const double* gravity_; // Solvers @@ -117,12 +115,11 @@ namespace Opm const BlackoilPropsAdInterface& props, const RockCompressibility* rock_comp_props, WellsManager& wells_manager, - const std::vector& src, const FlowBoundaryConditions* bcs, LinearSolverInterface& linsolver, const double* gravity) { - pimpl_.reset(new Impl(param, grid, props, rock_comp_props, wells_manager, src, bcs, linsolver, gravity)); + pimpl_.reset(new Impl(param, grid, props, rock_comp_props, wells_manager, bcs, linsolver, gravity)); } @@ -231,13 +228,12 @@ namespace Opm #endif - // \TODO: make CompressibleTpfa take src and bcs. + // \TODO: Treat bcs properly. SimulatorFullyImplicitBlackoil::Impl::Impl(const parameter::ParameterGroup& param, const UnstructuredGrid& grid, const BlackoilPropsAdInterface& props, const RockCompressibility* rock_comp_props, WellsManager& wells_manager, - const std::vector& src, const FlowBoundaryConditions* bcs, LinearSolverInterface& linsolver, const double* gravity) @@ -246,7 +242,6 @@ namespace Opm rock_comp_props_(rock_comp_props), wells_manager_(wells_manager), wells_(wells_manager.c_wells()), - src_(src), bcs_(bcs), gravity_(gravity), geo_(grid_, props_, gravity_), @@ -256,6 +251,10 @@ namespace Opm param.getDefault("nl_pressure_maxiter", 10), gravity, */ { + // Intercept usage of bcs, since we do not handle it. + if (bcs) { + OPM_THROW(std::runtime_error, "SimulatorFullyImplicitBlackoil cannot handle boundary conditions other than no-flow. Not implemented yet."); + } // For output. output_ = param.getDefault("output", true); if (output_) { diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoil.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoil.hpp index dd1fab6cb..71619ceca 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoil.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoil.hpp @@ -63,7 +63,6 @@ namespace Opm /// \param[in] props fluid and rock properties /// \param[in] rock_comp_props if non-null, rock compressibility properties /// \param[in] well_manager well manager, may manage no (null) wells - /// \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 @@ -72,7 +71,6 @@ namespace Opm const BlackoilPropsAdInterface& props, const RockCompressibility* rock_comp_props, WellsManager& wells_manager, - const std::vector& src, const FlowBoundaryConditions* bcs, LinearSolverInterface& linsolver, const double* gravity); From 56e50f02fda180af4fb1c9a27b248757e03979f1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 20 Sep 2013 14:42:12 +0200 Subject: [PATCH 11/14] Documented BlackoilPropsAdFromDeck. --- opm/autodiff/BlackoilPropsAdFromDeck.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/opm/autodiff/BlackoilPropsAdFromDeck.hpp b/opm/autodiff/BlackoilPropsAdFromDeck.hpp index 09d07459e..f8cf5f7ad 100644 --- a/opm/autodiff/BlackoilPropsAdFromDeck.hpp +++ b/opm/autodiff/BlackoilPropsAdFromDeck.hpp @@ -35,9 +35,9 @@ namespace Opm class SinglePvtInterface; - /// This class is intended to present a fluid interface for - /// three-phase black-oil that is easy to use with the AD-using - /// simulators. + /// This class implements the AD-adapted fluid interface for + /// three-phase black-oil. It requires an input deck from which it + /// reads all relevant property data. /// /// Most methods are available in two overloaded versions, one /// taking a constant vector and returning the same, and one From cb892d6a18125fdcd1ec4ad20fb6f5c2f7c00be8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 20 Sep 2013 14:55:24 +0200 Subject: [PATCH 12/14] Documented class DerivedGeology. --- opm/autodiff/GeoProps.hpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/opm/autodiff/GeoProps.hpp b/opm/autodiff/GeoProps.hpp index a10dc5fc2..2b321c1c5 100644 --- a/opm/autodiff/GeoProps.hpp +++ b/opm/autodiff/GeoProps.hpp @@ -27,6 +27,11 @@ namespace Opm { + /// Class containing static geological properties that can be + /// derived from grid and petrophysical properties: + /// - pore volume + /// - transmissibilities + /// - gravity potentials class DerivedGeology { public: From edd7e1487bf3b13eea5d4337146babf37a93ed63 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 20 Sep 2013 14:55:43 +0200 Subject: [PATCH 13/14] Documented class FullyImplicitBlackoilSolver. --- opm/autodiff/FullyImplicitBlackoilSolver.hpp | 22 +++++++++++++++++++- 1 file changed, 21 insertions(+), 1 deletion(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.hpp b/opm/autodiff/FullyImplicitBlackoilSolver.hpp index a909745cf..399fb4051 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.hpp @@ -36,10 +36,27 @@ namespace Opm { class WellState; - /// A fully implicit TPFA-based solver for the black-oil problem. + /// A fully implicit solver for the black-oil problem. + /// + /// The simulator is capable of handling three-phase problems + /// where gas can be dissolved in oil (but not vice versa). It + /// uses an industry-standard TPFA discretization with per-phase + /// upwind weighting of mobilities. + /// + /// It uses automatic differentiation via the class AutoDiffBlock + /// to simplify assembly of the jacobian matrix. class FullyImplicitBlackoilSolver { public: + /// Construct a solver. It will retain references to the + /// arguments of this functions, and they are expected to + /// remain in scope for the lifetime of the solver. + /// \param[in] grid grid data structure + /// \param[in] fluid fluid properties + /// \param[in] geo rock properties + /// \param[in] rock_comp_props if non-null, rock compressibility properties + /// \param[in] wells well structure + /// \param[in] linsolver linear solver FullyImplicitBlackoilSolver(const UnstructuredGrid& grid , const BlackoilPropsAdInterface& fluid, const DerivedGeology& geo , @@ -53,6 +70,9 @@ namespace Opm { /// state.saturation() /// state.gasoilratio() /// wstate.bhp() + /// \param[in] dt time step size + /// \param[in] state reservoir state + /// \param[in] wstate well state void step(const double dt , BlackoilState& state , From cc58bc3cef2120b45649ceab48cb68a8948618ad Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 23 Sep 2013 13:02:56 +0200 Subject: [PATCH 14/14] Fix minor issues pointed out by bska. --- doc/doxygen/Doxylocal | 3 ++- examples/sim_2p_comp_ad.cpp | 2 +- opm/autodiff/AutoDiffBlock.hpp | 18 +++++++++--------- opm/autodiff/AutoDiffHelpers.hpp | 6 +++--- opm/autodiff/BlackoilPropsAd.hpp | 11 +++++------ opm/autodiff/GeoProps.hpp | 2 +- 6 files changed, 21 insertions(+), 21 deletions(-) diff --git a/doc/doxygen/Doxylocal b/doc/doxygen/Doxylocal index db859f51a..da261f76f 100644 --- a/doc/doxygen/Doxylocal +++ b/doc/doxygen/Doxylocal @@ -9,7 +9,8 @@ INPUT += @abs_top_srcdir@/opm/ # see e.g. dune-grid for the examples of mainpage and modules -# INPUT += @abs_top_srcdir@/doc/doxygen/mainpage # @srcdir@/modules +# INPUT += @abs_top_srcdir@/doc/doxygen/mainpage +# @srcdir@/modules # The EXCLUDE tag can be used to specify files and/or directories that should # excluded from the INPUT source files. This way you can easily exclude a diff --git a/examples/sim_2p_comp_ad.cpp b/examples/sim_2p_comp_ad.cpp index a511d5d95..8a06df3ca 100644 --- a/examples/sim_2p_comp_ad.cpp +++ b/examples/sim_2p_comp_ad.cpp @@ -155,7 +155,7 @@ try simple_wells = create_wells(2, 2, 2); const double inj_frac[2] = { 1.0, 0.0 }; const int inj_cell = 0; - const double WI = 1e-8; + const double WI = 1e-8; // This is a completely made-up number. const double all_fluids[2] = { 1.0, 1.0 }; int ok = add_well(INJECTOR, 0.0, 1, inj_frac, &inj_cell, &WI, "Injector", simple_wells); ok = ok && append_well_controls(SURFACE_RATE, 0.01*flow_per_sec, all_fluids, 0, simple_wells); diff --git a/opm/autodiff/AutoDiffBlock.hpp b/opm/autodiff/AutoDiffBlock.hpp index 56fbfd451..b10cbfdc6 100644 --- a/opm/autodiff/AutoDiffBlock.hpp +++ b/opm/autodiff/AutoDiffBlock.hpp @@ -32,15 +32,15 @@ namespace Opm /// values and sparse jacobian matrices. /// /// The class contains a (column) vector of values and multiple - /// sparse matrices representing its derivatives. Each such matrix - /// has a number of rows equal to the number of rows in the value - /// vector, and a number of columns equal to the number of - /// variables we want to compute the derivatives with respect - /// to. The reason to have multiple such jacobians instead of just - /// one is to allow simpler grouping of variables, making it - /// easier to implement various preconditioning schemes. Only - /// basic arithmetic operators are implemented for this class, - /// reflecting our needs so far. + /// sparse matrices representing its partial derivatives. Each + /// such matrix has a number of rows equal to the number of rows + /// in the value vector, and a number of columns equal to the + /// number of discrete variables we want to compute the + /// derivatives with respect to. The reason to have multiple such + /// jacobians instead of just one is to allow simpler grouping of + /// variables, making it easier to implement various + /// preconditioning schemes. Only basic arithmetic operators are + /// implemented for this class, reflecting our needs so far. /// /// The class is built on the Eigen library, using an Eigen array /// type to contain the values and Eigen sparse matrices for the diff --git a/opm/autodiff/AutoDiffHelpers.hpp b/opm/autodiff/AutoDiffHelpers.hpp index 8959c2c79..66472b83b 100644 --- a/opm/autodiff/AutoDiffHelpers.hpp +++ b/opm/autodiff/AutoDiffHelpers.hpp @@ -42,15 +42,15 @@ struct HelperOps typedef Eigen::Array IFaces; IFaces internal_faces; - /// Extract for each internal face the difference of its adjacent cells'values (first - second). + /// Extract for each internal face the difference of its adjacent cells' values (first - second). M ngrad; - /// Extract for each face the difference of its adjacent cells'values (second - first). + /// Extract for each face the difference of its adjacent cells' values (second - first). M grad; /// Extract for each face the average of its adjacent cells' values. M caver; /// Extract for each cell the sum of its adjacent interior faces' (signed) values. M div; - /// Extract for each face the difference of its adjacent cells'values (first - second). + /// Extract for each face the difference of its adjacent cells' values (first - second). /// For boundary faces, one of the entries per row (corresponding to the outside) is zero. M fullngrad; /// Extract for each cell the sum of all its adjacent faces' (signed) values. diff --git a/opm/autodiff/BlackoilPropsAd.hpp b/opm/autodiff/BlackoilPropsAd.hpp index b00a00f29..70545d4f6 100644 --- a/opm/autodiff/BlackoilPropsAd.hpp +++ b/opm/autodiff/BlackoilPropsAd.hpp @@ -34,12 +34,11 @@ namespace Opm /// /// It is implemented by wrapping a BlackoilPropertiesInterface /// object (the interface class defined in opm-core) and calling - /// its methods. This approach works well for most methods, but - /// the rsMax() method cannot be implemented by such a wrapping, - /// without access to the underlying pvt objects. Therefore we - /// cannot use this class with any case that involves - /// miscibility. A rethinking of fluid interfaces is probably - /// necessary. + /// its methods. This class does not implement rsMax() because the + /// required information is not available when wrapping a + /// BlackoilPropertiesInterface. Consequently, class + /// BlackoilPropsAd cannot be used to simulate problems involving + /// miscibility. /// /// Most methods are available in two overloaded versions, one /// taking a constant vector and returning the same, and one diff --git a/opm/autodiff/GeoProps.hpp b/opm/autodiff/GeoProps.hpp index 2b321c1c5..fbf646295 100644 --- a/opm/autodiff/GeoProps.hpp +++ b/opm/autodiff/GeoProps.hpp @@ -27,7 +27,7 @@ namespace Opm { - /// Class containing static geological properties that can be + /// Class containing static geological properties that are /// derived from grid and petrophysical properties: /// - pore volume /// - transmissibilities