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

Merge pull request #32 from atgeirr/minor_improvements

Browse Source 
A collection of minor improvements
This commit is contained in:
Bård Skaflestad
2013-09-24 01:31:15 -07:00
parent 9109dad0de cc58bc3cef
commit 243432142e
31 changed files with 588 additions and 430 deletions
Download Patch File Download Diff File Expand all files Collapse all files

8
CMakeLists_files.cmake
View File

@@ -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
)

26
doc/doxygen/Doxylocal Normal file
View File

@@ -0,0 +1,26 @@
# 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/

2
examples/find_zero.cpp
View File

@@ -80,7 +80,7 @@ public:
{
double x = initial_guess;
iterations_used = 0;
typedef AutoDiff::Forward<double> AD;
typedef Opm::AutoDiff<double> AD;
while (std::abs(f(x)) > tolerance && ++iterations_used < max_iter) {
AD xfad = AD::variable(x);
AD rfad = f(xfad);

37
examples/sim_2p_comp_ad.cpp
View File

@@ -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<EclipseGridParser> deck;
boost::scoped_ptr<GridManager> grid;
boost::scoped_ptr<BlackoilPropertiesInterface> 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<double> 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<double> porevol;
@@ -150,8 +152,21 @@ try
const double default_injection = use_gravity ? 0.0 : 0.1;
const double flow_per_sec = param.getDefault<double>("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; // 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);
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);

38
examples/sim_2p_incomp_adfi.cpp → examples/sim_2p_incomp_ad.cpp
View File

@@ -44,7 +44,7 @@
#include <opm/core/simulator/TwophaseState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/autodiff/SimulatorIncompTwophaseAdfi.hpp>
#include <opm/autodiff/SimulatorIncompTwophaseAd.hpp>
#include <boost/scoped_ptr.hpp>
#include <boost/filesystem.hpp>
@@ -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);
}

223
examples/sim_fibo_ad.cpp
View File

@@ -45,7 +45,6 @@
#include <opm/core/simulator/WellState.hpp>
#include <opm/autodiff/SimulatorFullyImplicitBlackoil.hpp>
#include <opm/autodiff/BlackoilPropsAd.hpp>
#include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
#include <boost/scoped_ptr.hpp>
@@ -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<EclipseGridParser> deck;
boost::scoped_ptr<GridManager> grid;
boost::scoped_ptr<BlackoilPropertiesInterface> 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<std::string>("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<std::string>("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<double> 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<double> 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<double>("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";

8
examples/sim_simple.cpp
View File

@@ -111,7 +111,7 @@ fluxFunc(const std::vector<ADB>& m)
int main()
try
{
typedef AutoDiff::ForwardBlock<double> ADB;
typedef Opm::AutoDiffBlock<double> 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<double> ADB;
typedef Opm::AutoDiffBlock<double> 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<double> upwind(grid, ops, dflux);
const Opm::UpwindSelector<double> upwind(grid, ops, dflux);
const V pv = Eigen::Map<const V>(props.porosity(), nc, 1)
* Eigen::Map<const V>(grid.cell_volumes, nc, 1);
const double dt = 0.0005;

2
examples/test_impestpfa_ad.cpp
View File

@@ -58,8 +58,6 @@ try
const Opm::BlackoilPropertiesBasic oldprops(param, 2, nc);
const Opm::BlackoilPropsAd props(oldprops);
typedef AutoDiff::ForwardBlock<double> ADB;
Wells* wells = create_wells(2, 2, 5);
const double inj_frac[] = { 1.0, 0.0 };
const double prod_frac[] = { 0.0, 0.0 };

159
opm/autodiff/AutoDiff.hpp
View File

@@ -1,18 +1,3 @@
/*===========================================================================
//
// File: AutoDiff.hpp
//
// Created: 2013-04-29 10:51:23+0200
//
// Authors: Knut-Andreas Lie <Knut-Andreas.Lie@sintef.no>
// Halvor M. Nilsen <HalvorMoll.Nilsen@sintef.no>
// Atgeirr F. Rasmussen <atgeirr@sintef.no>
// Xavier Raynaud <Xavier.Raynaud@sintef.no>
// Bård Skaflestad <Bard.Skaflestad@sintef.no>
//
//==========================================================================*/
/*
Copyright 2013 SINTEF ICT, Applied Mathematics.
Copyright 2013 Statoil ASA.
@@ -33,33 +18,46 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <cmath>
#ifndef OPM_AUTODIFF_HPP_HEADER
#define OPM_AUTODIFF_HPP_HEADER
namespace AutoDiff {
#include <cmath>
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 <typename Scalar>
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 <class Ostream, typename Scalar>
Ostream&
operator<<(Ostream& os, const Forward<Scalar>& fw)
operator<<(Ostream& os, const AutoDiff<Scalar>& fw)
{
return fw.print(os);
}
template <typename Scalar>
Forward<Scalar>
operator +(const Forward<Scalar>& lhs,
const Forward<Scalar>& rhs)
AutoDiff<Scalar>
operator +(const AutoDiff<Scalar>& lhs,
const AutoDiff<Scalar>& rhs)
{
Forward<Scalar> ret = lhs;
AutoDiff<Scalar> ret = lhs;
ret += rhs;
return ret;
}
template <typename Scalar, typename T>
Forward<Scalar>
AutoDiff<Scalar>
operator +(const T lhs,
const Forward<Scalar>& rhs)
const AutoDiff<Scalar>& rhs)
{
Forward<Scalar> ret = rhs;
AutoDiff<Scalar> ret = rhs;
ret += Scalar(lhs);
return ret;
}
template <typename Scalar, typename T>
Forward<Scalar>
operator +(const Forward<Scalar>& lhs,
AutoDiff<Scalar>
operator +(const AutoDiff<Scalar>& lhs,
const T rhs)
{
Forward<Scalar> ret = lhs;
AutoDiff<Scalar> ret = lhs;
ret += Scalar(rhs);
return ret;
}
template <typename Scalar>
Forward<Scalar>
operator -(const Forward<Scalar>& lhs,
const Forward<Scalar>& rhs)
AutoDiff<Scalar>
operator -(const AutoDiff<Scalar>& lhs,
const AutoDiff<Scalar>& rhs)
{
Forward<Scalar> ret = lhs;
AutoDiff<Scalar> ret = lhs;
ret -= rhs;
return ret;
}
template <typename Scalar, typename T>
Forward<Scalar>
AutoDiff<Scalar>
operator -(const T lhs,
const Forward<Scalar>& rhs)
const AutoDiff<Scalar>& rhs)
{
return Forward<Scalar>::function(Scalar(lhs) - rhs.val(), -rhs.der());
return AutoDiff<Scalar>::function(Scalar(lhs) - rhs.val(), -rhs.der());
}
template <typename Scalar, typename T>
Forward<Scalar>
operator -(const Forward<Scalar>& lhs,
AutoDiff<Scalar>
operator -(const AutoDiff<Scalar>& lhs,
const T rhs)
{
Forward<Scalar> ret = lhs;
AutoDiff<Scalar> ret = lhs;
ret -= Scalar(rhs);
return ret;
}
template <typename Scalar>
Forward<Scalar>
operator *(const Forward<Scalar>& lhs,
const Forward<Scalar>& rhs)
AutoDiff<Scalar>
operator *(const AutoDiff<Scalar>& lhs,
const AutoDiff<Scalar>& rhs)
{
Forward<Scalar> ret = lhs;
AutoDiff<Scalar> ret = lhs;
ret *= rhs;
return ret;
}
template <typename Scalar, typename T>
Forward<Scalar>
AutoDiff<Scalar>
operator *(const T lhs,
const Forward<Scalar>& rhs)
const AutoDiff<Scalar>& rhs)
{
Forward<Scalar> ret = rhs;
AutoDiff<Scalar> ret = rhs;
ret *= Scalar(lhs);
return ret;
}
template <typename Scalar, typename T>
Forward<Scalar>
operator *(const Forward<Scalar>& lhs,
AutoDiff<Scalar>
operator *(const AutoDiff<Scalar>& lhs,
const T rhs)
{
Forward<Scalar> ret = lhs;
AutoDiff<Scalar> ret = lhs;
ret *= Scalar(rhs);
return ret;
}
template <typename Scalar>
Forward<Scalar>
operator /(const Forward<Scalar>& lhs,
const Forward<Scalar>& rhs)
AutoDiff<Scalar>
operator /(const AutoDiff<Scalar>& lhs,
const AutoDiff<Scalar>& rhs)
{
Forward<Scalar> ret = lhs;
AutoDiff<Scalar> ret = lhs;
ret /= rhs;
return ret;
}
template <typename Scalar, typename T>
Forward<Scalar>
AutoDiff<Scalar>
operator /(const T lhs,
const Forward<Scalar>& rhs)
const AutoDiff<Scalar>& rhs)
{
Scalar a = Scalar(lhs) / rhs.val();
Scalar b = -Scalar(lhs) / (rhs.val() * rhs.val());
return Forward<Scalar>::function(a, b);
return AutoDiff<Scalar>::function(a, b);
}
template <typename Scalar, typename T>
Forward<Scalar>
operator /(const Forward<Scalar>& lhs,
AutoDiff<Scalar>
operator /(const AutoDiff<Scalar>& lhs,
const T rhs)
{
Scalar a = lhs.val() / Scalar(rhs);
Scalar b = lhs.der() / Scalar(rhs);
return Forward<Scalar>::function(a, b);
return AutoDiff<Scalar>::function(a, b);
}
template <typename Scalar>
Forward<Scalar>
cos(const Forward<Scalar>& x)
AutoDiff<Scalar>
cos(const AutoDiff<Scalar>& x)
{
Scalar a = std::cos(x.val());
Scalar b = -std::sin(x.val()) * x.der();
return Forward<Scalar>::function(a, b);
return AutoDiff<Scalar>::function(a, b);
}
template <typename Scalar>
Forward<Scalar>
sqrt(const Forward<Scalar>& x)
AutoDiff<Scalar>
sqrt(const AutoDiff<Scalar>& x)
{
Scalar a = std::sqrt(x.val());
Scalar b = (Scalar(1.0) / (Scalar(2.0) * a)) * x.der();
return Forward<Scalar>::function(a, b);
return AutoDiff<Scalar>::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 */

210
opm/autodiff/AutoDiffBlock.hpp
View File

@@ -20,32 +20,86 @@
#ifndef OPM_AUTODIFFBLOCK_HEADER_INCLUDED
#define OPM_AUTODIFFBLOCK_HEADER_INCLUDED
#include <opm/autodiff/AutoDiff.hpp>
#include <Eigen/Eigen>
#include <Eigen/Sparse>
#include <vector>
#include <cassert>
namespace AutoDiff
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 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
/// jacobians. The overloaded operators are intended to behave in
/// 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.
///
/// 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:
/// - 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.
/// - 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 <typename Scalar>
class ForwardBlock
class AutoDiffBlock
{
public:
/// Underlying types for scalar vectors and jacobians.
/// Underlying type for values.
typedef Eigen::Array<Scalar, Eigen::Dynamic, 1> V;
/// Underlying type for jacobians.
typedef Eigen::SparseMatrix<Scalar> M;
/// Named constructor pattern used here.
static ForwardBlock null()
/// Construct an empty AutoDiffBlock.
static AutoDiffBlock null()
{
V val;
std::vector<M> jac;
return ForwardBlock(val, jac);
return AutoDiffBlock(val, jac);
}
static ForwardBlock constant(const V& val, const std::vector<int>& blocksizes)
/// Create an AutoDiffBlock representing a constant.
/// \param[in] val values
/// \param[in] blocksizes block pattern
static AutoDiffBlock constant(const V& val, const std::vector<int>& blocksizes)
{
std::vector<M> jac;
const int num_elem = val.size();
@@ -55,10 +109,17 @@ 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<int>& blocksizes)
/// 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<int>& blocksizes)
{
std::vector<M> jac;
const int num_elem = val.size();
@@ -74,32 +135,35 @@ 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<M>& 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<M>& jac)
{
return ForwardBlock(val, jac);
return AutoDiffBlock(val, jac);
}
/// Construct a set of primary variables,
/// each initialized to a given vector.
static std::vector<ForwardBlock> variables(const std::vector<V>& initial_values)
/// Construct a set of primary variables, each initialized to
/// a given vector.
static std::vector<AutoDiffBlock> variables(const std::vector<V>& initial_values)
{
const int num_vars = initial_values.size();
std::vector<int> bpat;
for (int v = 0; v < num_vars; ++v) {
bpat.push_back(initial_values[v].size());
}
std::vector<ForwardBlock> vars;
std::vector<AutoDiffBlock> vars;
for (int v = 0; v < num_vars; ++v) {
vars.emplace_back(variable(v, initial_values[v], bpat));
}
return vars;
}
/// Operator +=
ForwardBlock& operator+=(const ForwardBlock& rhs)
/// Elementwise operator +=
AutoDiffBlock& operator+=(const AutoDiffBlock& rhs)
{
assert (numBlocks() == rhs.numBlocks());
assert (value().size() == rhs.value().size());
@@ -116,8 +180,8 @@ namespace AutoDiff
return *this;
}
/// Operator +
ForwardBlock operator+(const ForwardBlock& rhs) const
/// Elementwise operator +
AutoDiffBlock operator+(const AutoDiffBlock& rhs) const
{
std::vector<M> jac = jac_;
assert(numBlocks() == rhs.numBlocks());
@@ -130,8 +194,8 @@ namespace AutoDiff
return function(val_ + rhs.val_, jac);
}
/// Operator -
ForwardBlock operator-(const ForwardBlock& rhs) const
/// Elementwise operator -
AutoDiffBlock operator-(const AutoDiffBlock& rhs) const
{
std::vector<M> jac = jac_;
assert(numBlocks() == rhs.numBlocks());
@@ -144,8 +208,8 @@ namespace AutoDiff
return function(val_ - rhs.val_, jac);
}
/// Operator *
ForwardBlock operator*(const ForwardBlock& rhs) const
/// Elementwise operator *
AutoDiffBlock operator*(const AutoDiffBlock& rhs) const
{
int num_blocks = numBlocks();
std::vector<M> jac(num_blocks);
@@ -161,8 +225,8 @@ namespace AutoDiff
return function(val_ * rhs.val_, jac);
}
/// Operator /
ForwardBlock operator/(const ForwardBlock& rhs) const
/// Elementwise operator /
AutoDiffBlock operator/(const AutoDiffBlock& rhs) const
{
int num_blocks = numBlocks();
std::vector<M> jac(num_blocks);
@@ -178,6 +242,7 @@ namespace AutoDiff
}
return function(val_ / rhs.val_, jac);
}
/// I/O.
template <class Ostream>
Ostream&
@@ -214,21 +279,21 @@ namespace AutoDiff
return bp;
}
/// Function value
/// Function value.
const V& value() const
{
return val_;
}
/// Function derivatives
/// Function derivatives.
const std::vector<M>& derivative() const
{
return jac_;
}
private:
ForwardBlock(const V& val,
const std::vector<M>& jac)
AutoDiffBlock(const V& val,
const std::vector<M>& jac)
: val_(val), jac_(jac)
{
#ifndef NDEBUG
@@ -242,12 +307,15 @@ namespace AutoDiff
V val_;
std::vector<M> jac_;
};
};
// --------- Free functions and operators for AutoDiffBlock ---------
/// Stream output.
template <class Ostream, typename Scalar>
Ostream&
operator<<(Ostream& os, const ForwardBlock<Scalar>& fw)
operator<<(Ostream& os, const AutoDiffBlock<Scalar>& fw)
{
return fw.print(os);
}
@@ -255,81 +323,89 @@ namespace AutoDiff
/// Multiply with sparse matrix from the left.
template <typename Scalar>
ForwardBlock<Scalar> operator*(const typename ForwardBlock<Scalar>::M& lhs,
const ForwardBlock<Scalar>& rhs)
AutoDiffBlock<Scalar> operator*(const typename AutoDiffBlock<Scalar>::M& lhs,
const AutoDiffBlock<Scalar>& rhs)
{
int num_blocks = rhs.numBlocks();
std::vector<typename ForwardBlock<Scalar>::M> jac(num_blocks);
std::vector<typename AutoDiffBlock<Scalar>::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<Scalar>::V val = lhs*rhs.value().matrix();
return ForwardBlock<Scalar>::function(val, jac);
typename AutoDiffBlock<Scalar>::V val = lhs*rhs.value().matrix();
return AutoDiffBlock<Scalar>::function(val, jac);
}
/// Elementwise multiplication with constant on the left.
template <typename Scalar>
ForwardBlock<Scalar> operator*(const typename ForwardBlock<Scalar>::V& lhs,
const ForwardBlock<Scalar>& rhs)
AutoDiffBlock<Scalar> operator*(const typename AutoDiffBlock<Scalar>::V& lhs,
const AutoDiffBlock<Scalar>& rhs)
{
return ForwardBlock<Scalar>::constant(lhs, rhs.blockPattern()) * rhs;
return AutoDiffBlock<Scalar>::constant(lhs, rhs.blockPattern()) * rhs;
}
/// Elementwise multiplication with constant on the right.
template <typename Scalar>
ForwardBlock<Scalar> operator*(const ForwardBlock<Scalar>& lhs,
const typename ForwardBlock<Scalar>::V& rhs)
AutoDiffBlock<Scalar> operator*(const AutoDiffBlock<Scalar>& lhs,
const typename AutoDiffBlock<Scalar>::V& rhs)
{
return rhs * lhs; // Commutative operation.
}
/// Elementwise addition with constant on the left.
template <typename Scalar>
ForwardBlock<Scalar> operator+(const typename ForwardBlock<Scalar>::V& lhs,
const ForwardBlock<Scalar>& rhs)
AutoDiffBlock<Scalar> operator+(const typename AutoDiffBlock<Scalar>::V& lhs,
const AutoDiffBlock<Scalar>& rhs)
{
return ForwardBlock<Scalar>::constant(lhs, rhs.blockPattern()) + rhs;
return AutoDiffBlock<Scalar>::constant(lhs, rhs.blockPattern()) + rhs;
}
/// Elementwise addition with constant on the right.
template <typename Scalar>
ForwardBlock<Scalar> operator+(const ForwardBlock<Scalar>& lhs,
const typename ForwardBlock<Scalar>::V& rhs)
AutoDiffBlock<Scalar> operator+(const AutoDiffBlock<Scalar>& lhs,
const typename AutoDiffBlock<Scalar>::V& rhs)
{
return rhs + lhs; // Commutative operation.
}
/// Elementwise subtraction with constant on the left.
template <typename Scalar>
ForwardBlock<Scalar> operator-(const typename ForwardBlock<Scalar>::V& lhs,
const ForwardBlock<Scalar>& rhs)
AutoDiffBlock<Scalar> operator-(const typename AutoDiffBlock<Scalar>::V& lhs,
const AutoDiffBlock<Scalar>& rhs)
{
return ForwardBlock<Scalar>::constant(lhs, rhs.blockPattern()) - rhs;
return AutoDiffBlock<Scalar>::constant(lhs, rhs.blockPattern()) - rhs;
}
/// Elementwise subtraction with constant on the right.
template <typename Scalar>
ForwardBlock<Scalar> operator-(const ForwardBlock<Scalar>& lhs,
const typename ForwardBlock<Scalar>::V& rhs)
AutoDiffBlock<Scalar> operator-(const AutoDiffBlock<Scalar>& lhs,
const typename AutoDiffBlock<Scalar>::V& rhs)
{
return lhs - ForwardBlock<Scalar>::constant(rhs, lhs.blockPattern());
return lhs - AutoDiffBlock<Scalar>::constant(rhs, lhs.blockPattern());
}
/// Elementwise division with constant on the left.
template <typename Scalar>
ForwardBlock<Scalar> operator/(const typename ForwardBlock<Scalar>::V& lhs,
const ForwardBlock<Scalar>& rhs)
AutoDiffBlock<Scalar> operator/(const typename AutoDiffBlock<Scalar>::V& lhs,
const AutoDiffBlock<Scalar>& rhs)
{
return ForwardBlock<Scalar>::constant(lhs, rhs.blockPattern()) / rhs;
return AutoDiffBlock<Scalar>::constant(lhs, rhs.blockPattern()) / rhs;
}
/// Elementwise division with constant on the right.
template <typename Scalar>
ForwardBlock<Scalar> operator/(const ForwardBlock<Scalar>& lhs,
const typename ForwardBlock<Scalar>::V& rhs)
AutoDiffBlock<Scalar> operator/(const AutoDiffBlock<Scalar>& lhs,
const typename AutoDiffBlock<Scalar>::V& rhs)
{
return lhs / ForwardBlock<Scalar>::constant(rhs, lhs.blockPattern());
return lhs / AutoDiffBlock<Scalar>::constant(rhs, lhs.blockPattern());
}
@@ -341,15 +417,15 @@ namespace AutoDiff
* @return The product
*/
template <typename Scalar>
ForwardBlock<Scalar> operator*(const ForwardBlock<Scalar>& lhs,
const Scalar& rhs)
AutoDiffBlock<Scalar> operator*(const AutoDiffBlock<Scalar>& lhs,
const Scalar& rhs)
{
std::vector< typename ForwardBlock<Scalar>::M > jac;
std::vector< typename AutoDiffBlock<Scalar>::M > jac;
jac.reserve( lhs.numBlocks() );
for (int block=0; block<lhs.numBlocks(); block++) {
jac.emplace_back( lhs.derivative()[block] * rhs );
}
return ForwardBlock<Scalar>::function( lhs.value() * rhs, jac );
return AutoDiffBlock<Scalar>::function( lhs.value() * rhs, jac );
}
@@ -361,14 +437,14 @@ namespace AutoDiff
* @return The product
*/
template <typename Scalar>
ForwardBlock<Scalar> operator*(const Scalar& lhs,
const ForwardBlock<Scalar>& rhs)
AutoDiffBlock<Scalar> operator*(const Scalar& lhs,
const AutoDiffBlock<Scalar>& rhs)
{
return rhs * lhs; // Commutative operation.
}
} // namespace Autodiff
} // namespace Opm

73
opm/autodiff/AutoDiffHelpers.hpp
View File

@@ -26,25 +26,35 @@
#include <iostream>
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<double>::M M;
typedef AutoDiff::ForwardBlock<double>::V V;
typedef AutoDiffBlock<double>::M M;
typedef AutoDiffBlock<double>::V V;
/// A list of internal faces.
typedef Eigen::Array<int, Eigen::Dynamic, 1> 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 +99,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<Tri> 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();
}
};
@@ -181,7 +205,7 @@ namespace {
template <typename Scalar>
class UpwindSelector {
public:
typedef AutoDiff::ForwardBlock<Scalar> ADB;
typedef AutoDiffBlock<Scalar> ADB;
UpwindSelector(const UnstructuredGrid& g,
const HelperOps& h,
@@ -276,11 +300,11 @@ namespace {
/// Returns x(indices).
template <typename Scalar, class IntVec>
AutoDiff::ForwardBlock<Scalar>
subset(const AutoDiff::ForwardBlock<Scalar>& x,
AutoDiffBlock<Scalar>
subset(const AutoDiffBlock<Scalar>& x,
const IntVec& indices)
{
return ::constructSubsetSparseMatrix<Scalar>(x.value().size(), indices) * x;
return constructSubsetSparseMatrix<Scalar>(x.value().size(), indices) * x;
}
@@ -291,19 +315,19 @@ Eigen::Array<Scalar, Eigen::Dynamic, 1>
subset(const Eigen::Array<Scalar, Eigen::Dynamic, 1>& x,
const IntVec& indices)
{
return (::constructSubsetSparseMatrix<Scalar>(x.size(), indices) * x.matrix()).array();
return (constructSubsetSparseMatrix<Scalar>(x.size(), indices) * x.matrix()).array();
}
/// Returns v where v(indices) == x, v(!indices) == 0 and v.size() == n.
template <typename Scalar, class IntVec>
AutoDiff::ForwardBlock<Scalar>
superset(const AutoDiff::ForwardBlock<Scalar>& x,
AutoDiffBlock<Scalar>
superset(const AutoDiffBlock<Scalar>& x,
const IntVec& indices,
const int n)
{
return ::constructSupersetSparseMatrix<Scalar>(n, indices) * x;
return constructSupersetSparseMatrix<Scalar>(n, indices) * x;
}
@@ -315,7 +339,7 @@ superset(const Eigen::Array<Scalar, Eigen::Dynamic, 1>& x,
const IntVec& indices,
const int n)
{
return ::constructSupersetSparseMatrix<Scalar>(n, indices) * x.matrix();
return constructSupersetSparseMatrix<Scalar>(n, indices) * x.matrix();
}
@@ -324,10 +348,10 @@ superset(const Eigen::Array<Scalar, Eigen::Dynamic, 1>& 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<double>::M
spdiag(const AutoDiff::ForwardBlock<double>::V& d)
AutoDiffBlock<double>::M
spdiag(const AutoDiffBlock<double>::V& d)
{
typedef AutoDiff::ForwardBlock<double>::M M;
typedef AutoDiffBlock<double>::M M;
const int n = d.size();
M mat(n, n);
@@ -346,7 +370,7 @@ spdiag(const AutoDiff::ForwardBlock<double>::V& d)
template <typename Scalar>
class Selector {
public:
typedef AutoDiff::ForwardBlock<Scalar> ADB;
typedef AutoDiffBlock<Scalar> ADB;
Selector(const typename ADB::V& selection_basis)
{
@@ -400,10 +424,10 @@ spdiag(const AutoDiff::ForwardBlock<double>::V& d)
/// Returns the input expression, but with all Jacobians collapsed to one.
inline
AutoDiff::ForwardBlock<double>
collapseJacs(const AutoDiff::ForwardBlock<double>& x)
AutoDiffBlock<double>
collapseJacs(const AutoDiffBlock<double>& x)
{
typedef AutoDiff::ForwardBlock<double> ADB;
typedef AutoDiffBlock<double> ADB;
const int nb = x.numBlocks();
typedef Eigen::Triplet<double> Tri;
int nnz = 0;
@@ -436,9 +460,9 @@ collapseJacs(const AutoDiff::ForwardBlock<double>& x)
/// Returns the vertical concatenation [ x; y ] of the inputs.
inline
AutoDiff::ForwardBlock<double>
vertcat(const AutoDiff::ForwardBlock<double>& x,
const AutoDiff::ForwardBlock<double>& y)
AutoDiffBlock<double>
vertcat(const AutoDiffBlock<double>& x,
const AutoDiffBlock<double>& y)
{
const int nx = x.size();
const int ny = y.size();
@@ -564,4 +588,7 @@ inline Eigen::ArrayXd sign (const Eigen::ArrayXd& x)
return retval;
}
} // namespace Opm
#endif // OPM_AUTODIFFHELPERS_HEADER_INCLUDED

15
opm/autodiff/BlackoilPropsAd.hpp
View File

@@ -29,9 +29,16 @@ 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 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
@@ -67,7 +74,7 @@ namespace Opm
// Fluid interface //
////////////////////////////
typedef AutoDiff::ForwardBlock<double> ADB;
typedef AutoDiffBlock<double> ADB;
typedef ADB::V V;
typedef std::vector<int> Cells;

8
opm/autodiff/BlackoilPropsAdFromDeck.hpp
View File

@@ -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
@@ -75,7 +75,7 @@ namespace Opm
// Fluid interface //
////////////////////////////
typedef AutoDiff::ForwardBlock<double> ADB;
typedef AutoDiffBlock<double> ADB;
typedef ADB::V V;
typedef std::vector<int> Cells;

2
opm/autodiff/BlackoilPropsAdInterface.hpp
View File

@@ -64,7 +64,7 @@ namespace Opm
// Fluid interface //
////////////////////////////
typedef AutoDiff::ForwardBlock<double> ADB;
typedef AutoDiffBlock<double> ADB;
typedef ADB::V V;
typedef ADB::M M;
typedef std::vector<int> Cells;

15
opm/autodiff/FullyImplicitBlackoilSolver.cpp
View File

@@ -37,6 +37,7 @@
#include <iostream>
#include <iomanip>
// A debugging utility.
#define DUMP(foo) \
do { \
std::cout << "==========================================\n" \
@@ -44,7 +45,11 @@
<< collapseJacs(foo) << std::endl; \
} while (0)
typedef AutoDiff::ForwardBlock<double> ADB;
namespace Opm {
typedef AutoDiffBlock<double> ADB;
typedef ADB::V V;
typedef ADB::M M;
typedef Eigen::Array<double,
@@ -69,7 +74,7 @@ namespace {
template <class GeoProps>
AutoDiff::ForwardBlock<double>::M
AutoDiffBlock<double>::M
gravityOperator(const UnstructuredGrid& grid,
const HelperOps& ops ,
const GeoProps& geo )
@@ -86,8 +91,8 @@ namespace {
}
}
typedef AutoDiff::ForwardBlock<double>::V V;
typedef AutoDiff::ForwardBlock<double>::M M;
typedef AutoDiffBlock<double>::V V;
typedef AutoDiffBlock<double>::M M;
const V& gpot = geo.gravityPotential();
const V& trans = geo.transmissibility();
@@ -182,8 +187,6 @@ namespace {
namespace Opm {
FullyImplicitBlackoilSolver::
FullyImplicitBlackoilSolver(const UnstructuredGrid& grid ,

24
opm/autodiff/FullyImplicitBlackoilSolver.hpp
View File

@@ -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 ,
@@ -60,7 +80,7 @@ namespace Opm {
private:
// Types and enums
typedef AutoDiff::ForwardBlock<double> ADB;
typedef AutoDiffBlock<double> ADB;
typedef ADB::V V;
typedef ADB::M M;
typedef Eigen::Array<double,

5
opm/autodiff/GeoProps.hpp
View File

@@ -27,6 +27,11 @@
namespace Opm
{
/// Class containing static geological properties that are
/// derived from grid and petrophysical properties:
/// - pore volume
/// - transmissibilities
/// - gravity potentials
class DerivedGeology
{
public:

12
opm/autodiff/ImpesTPFAAD.cpp
View File

@@ -29,9 +29,10 @@
#include <iostream>
#include <iomanip>
namespace Opm {
// Repeated from inside ImpesTPFAAD for convenience.
typedef AutoDiff::ForwardBlock<double> ADB;
typedef AutoDiffBlock<double> ADB;
typedef ADB::V V;
typedef ADB::M M;
@@ -48,7 +49,7 @@ namespace {
}
template <class GeoProps>
AutoDiff::ForwardBlock<double>::M
AutoDiffBlock<double>::M
gravityOperator(const UnstructuredGrid& grid,
const HelperOps& ops ,
const GeoProps& geo )
@@ -65,8 +66,8 @@ namespace {
}
}
typedef AutoDiff::ForwardBlock<double>::V V;
typedef AutoDiff::ForwardBlock<double>::M M;
typedef AutoDiffBlock<double>::V V;
typedef AutoDiffBlock<double>::M M;
const V& gpot = geo.gravityPotential();
const V& trans = geo.transmissibility();
@@ -122,7 +123,6 @@ namespace {
} // anonymous namespace
namespace Opm {
@@ -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);

2
opm/autodiff/ImpesTPFAAD.hpp
View File

@@ -63,7 +63,7 @@ namespace Opm {
ImpesTPFAAD& operator=(const ImpesTPFAAD& rhs);
// Types
typedef AutoDiff::ForwardBlock<double> ADB;
typedef AutoDiffBlock<double> ADB;
typedef ADB::V V;
typedef ADB::M M;
typedef Eigen::Array<double,

9
opm/autodiff/SimulatorCompressibleAd.cpp
View File

@@ -70,7 +70,6 @@ namespace Opm
const BlackoilPropertiesInterface& props,
const RockCompressibility* rock_comp_props,
WellsManager& wells_manager,
const std::vector<double>& 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<double>& 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<double>& 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<double>& 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_),

4
opm/autodiff/SimulatorCompressibleAd.hpp
View File

@@ -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<double>& src,
const FlowBoundaryConditions* bcs,
LinearSolverInterface& linsolver,
const double* gravity);

13
opm/autodiff/SimulatorFullyImplicitBlackoil.cpp
View File

@@ -70,7 +70,6 @@ namespace Opm
const BlackoilPropsAdInterface& props,
const RockCompressibility* rock_comp_props,
WellsManager& wells_manager,
const std::vector<double>& 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<double>& 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<double>& 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<double>& 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_) {

2
opm/autodiff/SimulatorFullyImplicitBlackoil.hpp
View File

@@ -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<double>& src,
const FlowBoundaryConditions* bcs,
LinearSolverInterface& linsolver,
const double* gravity);

12
opm/autodiff/SimulatorIncompTwophaseAdfi.cpp → opm/autodiff/SimulatorIncompTwophaseAd.cpp
View File

@@ -22,7 +22,7 @@
#include "config.h"
#endif // HAVE_CONFIG_H
#include <opm/autodiff/SimulatorIncompTwophaseAdfi.hpp>
#include <opm/autodiff/SimulatorIncompTwophaseAd.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
@@ -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)
{

24
opm/autodiff/SimulatorIncompTwophaseAdfi.hpp → opm/autodiff/SimulatorIncompTwophaseAd.hpp
View File

@@ -17,8 +17,8 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_SIMULATORINCOMPTWOPHASEADFI_HEADER_INCLUDED
#define OPM_SIMULATORINCOMPTWOPHASEADFI_HEADER_INCLUDED
#ifndef OPM_SIMULATORINCOMPTWOPHASEAD_HEADER_INCLUDED
#define OPM_SIMULATORINCOMPTWOPHASEAD_HEADER_INCLUDED
#include <boost/shared_ptr.hpp>
#include <vector>
@@ -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<double>& 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<double>& 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

2
opm/autodiff/TransportSolverTwophaseAd.hpp
View File

@@ -68,7 +68,7 @@ namespace Opm
TwophaseState& state);
private:
typedef AutoDiff::ForwardBlock<double> ADB;
typedef AutoDiffBlock<double> ADB;
typedef ADB::V V;
typedef ADB::M M;

49
opm/autodiff/opm-autodiff_doxygen_main.hpp Normal file
View File

@@ -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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_OPM-AUTODIFF_DOXYGEN_MAIN_HEADER_INCLUDED
#define OPM_OPM-AUTODIFF_DOXYGEN_MAIN_HEADER_INCLUDED
/** \mainpage Documentation for the opm-autodiff library.
<h3>Automatic differentiation</h3>
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.
<h3>Solvers and simulators</h3>
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

12
tests/test_block.cpp
View File

@@ -32,6 +32,8 @@
#include <Eigen/Eigen>
#include <Eigen/Sparse>
using namespace Opm;
namespace {
template <typename Scalar>
bool
@@ -73,7 +75,7 @@ namespace {
BOOST_AUTO_TEST_CASE(ConstantInitialisation)
{
typedef AutoDiff::ForwardBlock<double> ADB;
typedef AutoDiffBlock<double> ADB;
std::vector<int> blocksizes = { 3, 1, 2 };
@@ -92,7 +94,7 @@ BOOST_AUTO_TEST_CASE(ConstantInitialisation)
BOOST_AUTO_TEST_CASE(VariableInitialisation)
{
typedef AutoDiff::ForwardBlock<double> ADB;
typedef AutoDiffBlock<double> ADB;
std::vector<int> blocksizes = { 3, 1, 2 };
@@ -119,7 +121,7 @@ BOOST_AUTO_TEST_CASE(VariableInitialisation)
BOOST_AUTO_TEST_CASE(FunctionInitialisation)
{
typedef AutoDiff::ForwardBlock<double> ADB;
typedef AutoDiffBlock<double> ADB;
std::vector<int> blocksizes = { 3, 1, 2 };
std::vector<int>::size_type num_blocks = blocksizes.size();
@@ -150,7 +152,7 @@ BOOST_AUTO_TEST_CASE(FunctionInitialisation)
BOOST_AUTO_TEST_CASE(Addition)
{
typedef AutoDiff::ForwardBlock<double> ADB;
typedef AutoDiffBlock<double> ADB;
std::vector<int> blocksizes = { 3, 1, 2 };
ADB::V va(3);
@@ -195,7 +197,7 @@ BOOST_AUTO_TEST_CASE(Addition)
int main()
try
{
typedef AutoDiff::ForwardBlock<double> ADB;
typedef AutoDiffBlock<double> ADB;
std::vector<int> blocksizes = { 3, 1, 2 };
int num_blocks = blocksizes.size();
ADB::V v1(3);

4
tests/test_scalar_mult.cpp
View File

@@ -34,9 +34,11 @@
#include <iostream>
using namespace Opm;
BOOST_AUTO_TEST_CASE(ScalarMultiplication)
{
typedef AutoDiff::ForwardBlock<double> ADB;
typedef AutoDiffBlock<double> ADB;
std::vector<int> blocksizes = { 3, 1, 2 };
ADB::V vx(3);

2
tests/test_span.cpp
View File

@@ -32,6 +32,8 @@
#include <iostream>
using namespace Opm;
BOOST_AUTO_TEST_CASE(OneArgConstr)
{
const int num = 4;

16
tests/test_syntax.cpp
View File

@@ -15,7 +15,7 @@
BOOST_AUTO_TEST_CASE(Initialisation)
{
typedef AutoDiff::Forward<double> AdFW;
typedef Opm::AutoDiff<double> AdFW;
const double atol = 1.0e-14;
@@ -36,7 +36,7 @@ BOOST_AUTO_TEST_CASE(Initialisation)
BOOST_AUTO_TEST_CASE(Addition)
{
typedef AutoDiff::Forward<double> AdFW;
typedef Opm::AutoDiff<double> AdFW;
const double atol = 1.0e-14;
@@ -71,7 +71,7 @@ BOOST_AUTO_TEST_CASE(Addition)
BOOST_AUTO_TEST_CASE(Subtraction)
{
typedef AutoDiff::Forward<double> AdFW;
typedef Opm::AutoDiff<double> AdFW;
const double atol = 1.0e-14;
@@ -110,7 +110,7 @@ BOOST_AUTO_TEST_CASE(Subtraction)
BOOST_AUTO_TEST_CASE(Multiplication)
{
typedef AutoDiff::Forward<double> AdFW;
typedef Opm::AutoDiff<double> AdFW;
const double atol = 1.0e-14;
@@ -149,7 +149,7 @@ BOOST_AUTO_TEST_CASE(Multiplication)
BOOST_AUTO_TEST_CASE(Division)
{
typedef AutoDiff::Forward<double> AdFW;
typedef Opm::AutoDiff<double> AdFW;
const double atol = 1.0e-14;
@@ -186,7 +186,7 @@ BOOST_AUTO_TEST_CASE(Division)
BOOST_AUTO_TEST_CASE(Polynomial)
{
typedef AutoDiff::Forward<double> AdFW;
typedef Opm::AutoDiff<double> AdFW;
const double atol = 1.0e-14;
@@ -204,7 +204,7 @@ BOOST_AUTO_TEST_CASE(Polynomial)
BOOST_AUTO_TEST_CASE(Cosine)
{
typedef AutoDiff::Forward<double> AdFW;
typedef Opm::AutoDiff<double> AdFW;
const double atol = 1.0e-14;
@@ -223,7 +223,7 @@ BOOST_AUTO_TEST_CASE(Cosine)
BOOST_AUTO_TEST_CASE(SquareRoot)
{
typedef AutoDiff::Forward<double> AdFW;
typedef Opm::AutoDiff<double> AdFW;
const double atol = 1.0e-14;