OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge pull request #711 from atgeirr/eikonal-solver

Browse Source 
Solver for the anisotropic Eikonal equation using fast marching method
This commit is contained in:
Bård Skaflestad
2015-01-07 11:00:51 +01:00
parent a14b2fd434 14bc0aa7ee
commit 0820cab627
5 changed files with 793 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
CMakeLists_files.cmake
View File

@@ -108,6 +108,7 @@ list (APPEND MAIN_SOURCE_FILES
opm/core/simulator/SimulatorReport.cpp
opm/core/simulator/SimulatorState.cpp
opm/core/simulator/SimulatorTimer.cpp
opm/core/tof/AnisotropicEikonal.cpp
opm/core/tof/DGBasis.cpp
opm/core/tof/TofReorder.cpp
opm/core/tof/TofDiscGalReorder.cpp
@@ -189,6 +190,7 @@ list (APPEND TEST_SOURCE_FILES
tests/test_timer.cpp
tests/test_minpvprocessor.cpp
tests/test_gridutilities.cpp
tests/test_anisotropiceikonal.cpp
)
# originally generated with the command:
@@ -225,6 +227,7 @@ list (APPEND TEST_DATA_FILES
# originally generated with the command:
# find tutorials examples -name '*.c*' -printf '\t%p\n' | sort
list (APPEND EXAMPLE_SOURCE_FILES
examples/compute_eikonal_from_files.cpp
examples/compute_initial_state.cpp
examples/compute_tof.cpp
examples/compute_tof_from_files.cpp
@@ -371,6 +374,7 @@ list (APPEND PUBLIC_HEADER_FILES
opm/core/simulator/initState_impl.hpp
opm/core/simulator/initStateEquil.hpp
opm/core/simulator/initStateEquil_impl.hpp
opm/core/tof/AnisotropicEikonal.hpp
opm/core/tof/DGBasis.hpp
opm/core/tof/TofReorder.hpp
opm/core/tof/TofDiscGalReorder.hpp

130
examples/compute_eikonal_from_files.cpp Normal file
View File

@@ -0,0 +1,130 @@
/*
Copyright 2014 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/>.
*/
#if HAVE_CONFIG_H
#include "config.h"
#endif // HAVE_CONFIG_H
#include <opm/core/tof/AnisotropicEikonal.hpp>
#include <opm/core/grid.h>
#include <opm/core/grid/GridManager.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/StopWatch.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <boost/filesystem.hpp>
#include <memory>
#include <algorithm>
#include <iostream>
#include <vector>
#include <numeric>
#include <iterator>
namespace
{
void warnIfUnusedParams(const Opm::parameter::ParameterGroup& param)
{
if (param.anyUnused()) {
std::cout << "-------------------- Warning: unused parameters: --------------------\n";
param.displayUsage();
std::cout << "-------------------------------------------------------------------------" << std::endl;
}
}
} // anon namespace
// ----------------- Main program -----------------
int
main(int argc, char** argv)
try
{
using namespace Opm;
parameter::ParameterGroup param(argc, argv, false);
// Read grid.
GridManager grid_manager(param.get<std::string>("grid_filename"));
const UnstructuredGrid& grid = *grid_manager.c_grid();
// Read metric tensor.
std::vector<double> metric;
{
std::ifstream metric_stream(param.get<std::string>("metric_filename").c_str());
std::istream_iterator<double> beg(metric_stream);
std::istream_iterator<double> end;
metric.assign(beg, end);
if (int(metric.size()) != grid.number_of_cells*grid.dimensions*grid.dimensions) {
OPM_THROW(std::runtime_error, "Size of metric field differs from (dim^2 * number of cells).");
}
}
// Read starting cells.
std::vector<int> startcells;
{
std::ifstream start_stream(param.get<std::string>("startcells_filename").c_str());
std::istream_iterator<int> beg(start_stream);
std::istream_iterator<int> end;
startcells.assign(beg, end);
}
// Write parameters used for later reference.
bool output = param.getDefault("output", true);
std::string output_dir;
if (output) {
output_dir =
param.getDefault("output_dir", std::string("output"));
boost::filesystem::path fpath(output_dir);
try {
create_directories(fpath);
}
catch (...) {
OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
}
param.writeParam(output_dir + "/eikonal.param");
}
// Issue a warning if any parameters were unused.
warnIfUnusedParams(param);
// Solve eikonal equation.
Opm::time::StopWatch timer;
timer.start();
std::vector<double> solution;
AnisotropicEikonal2d ae(grid);
ae.solve(metric.data(), startcells, solution);
timer.stop();
double tt = timer.secsSinceStart();
std::cout << "Eikonal solver took: " << tt << " seconds." << std::endl;
// Output.
if (output) {
std::string filename = output_dir + "/solution.txt";
std::ofstream stream(filename.c_str());
stream.precision(16);
std::copy(solution.begin(), solution.end(), std::ostream_iterator<double>(stream, "\n"));
}
}
catch (const std::exception &e) {
std::cerr << "Program threw an exception: " << e.what() << "\n";
throw;
}

461
opm/core/tof/AnisotropicEikonal.cpp Normal file
View File

@@ -0,0 +1,461 @@
/*
Copyright 2014 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/>.
*/
#include <opm/core/tof/AnisotropicEikonal.hpp>
#include <opm/core/grid/GridUtilities.hpp>
#include <opm/core/grid.h>
#include <opm/core/utility/RootFinders.hpp>
#if BOOST_HEAP_AVAILABLE
namespace Opm
{
namespace
{
/// Euclidean (isotropic) distance.
double distanceIso(const double v1[2],
const double v2[2])
{
const double d[2] = { v2[0] - v1[0], v2[1] - v1[1] };
const double dist = std::sqrt(d[0]*d[0] + d[1]*d[1]);
return dist;
}
/// Anisotropic distance with respect to a metric g.
/// If d = v2 - v1, the distance is sqrt(d^T g d).
double distanceAniso(const double v1[2],
const double v2[2],
const double g[4])
{
const double d[2] = { v2[0] - v1[0], v2[1] - v1[1] };
const double dist = std::sqrt(+ g[0] * d[0] * d[0]
+ g[1] * d[0] * d[1]
+ g[2] * d[1] * d[0]
+ g[3] * d[1] * d[1]);
return dist;
}
} // anonymous namespace
/// Construct solver.
/// \param[in] grid A 2d grid.
AnisotropicEikonal2d::AnisotropicEikonal2d(const UnstructuredGrid& grid)
: grid_(grid),
safety_factor_(1.2)
{
if (grid.dimensions != 2) {
OPM_THROW(std::logic_error, "Grid for AnisotropicEikonal2d must be 2d.");
}
cell_neighbours_ = cellNeighboursAcrossVertices(grid);
orderCounterClockwise(grid, cell_neighbours_);
computeGridRadius();
}
/// Solve the eikonal equation.
/// \param[in] metric Array of metric tensors, M, for each cell.
/// \param[in] startcells Array of cells where u = 0 at the centroid.
/// \param[out] solution Array of solution to the eikonal equation.
void AnisotropicEikonal2d::solve(const double* metric,
const std::vector<int>& startcells,
std::vector<double>& solution)
{
// Compute anisotropy ratios to be used by isClose().
computeAnisoRatio(metric);
// The algorithm used is described in J.A. Sethian and A. Vladimirsky,
// "Ordered Upwind Methods for Static Hamilton-Jacobi Equations".
// Notation in comments is as used in that paper: U is the solution,
// and q is the boundary condition. One difference is that we talk about
// grid cells instead of mesh points.
//
// Algorithm summary:
// 1. Put all cells in Far. U_i = \inf.
// 2. Move the startcells to Accepted. U_i = q(x_i)
// 3. Move cells adjacent to startcells to Considered, evaluate
// U_i = min_{(x_j,x_k) \in NF(x_i)} G_{j,k}
// 4. Find the Considered cell with the smallest value: r.
// 5. Move cell r to Accepted. Update AcceptedFront.
// 6. Recompute the value for all Considered cells within
// distance h * F_2/F1 from x_r. Use min of previous and new.
// 7. Move cells adjacent to r from Far to Considered.
// 8. If Considered is not empty, go to step 4.
// 1. Put all cells in Far. U_i = \inf.
const int num_cells = grid_.number_of_cells;
const double inf = 1e100;
solution.clear();
solution.resize(num_cells, inf);
is_accepted_.clear();
is_accepted_.resize(num_cells, false);
accepted_front_.clear();
considered_.clear();
considered_handles_.clear();
is_considered_.clear();
is_considered_.resize(num_cells, false);
// 2. Move the startcells to Accepted. U_i = q(x_i)
const int num_startcells = startcells.size();
for (int ii = 0; ii < num_startcells; ++ii) {
is_accepted_[startcells[ii]] = true;
solution[startcells[ii]] = 0.0;
}
accepted_front_.insert(startcells.begin(), startcells.end());
// 3. Move cells adjacent to startcells to Considered, evaluate
// U_i = min_{(x_j,x_k) \in NF(x_i)} G_{j,k}
for (int ii = 0; ii < num_startcells; ++ii) {
const int scell = startcells[ii];
const int num_nb = cell_neighbours_[scell].size();
for (int nb = 0; nb < num_nb; ++nb) {
const int nb_cell = cell_neighbours_[scell][nb];
if (!is_accepted_[nb_cell] && !is_considered_[nb_cell]) {
const double value = computeValue(nb_cell, metric, solution.data());
pushConsidered(std::make_pair(value, nb_cell));
}
}
}
while (!considered_.empty()) {
// 4. Find the Considered cell with the smallest value: r.
const ValueAndCell r = topConsidered();
// std::cout << "Accepting cell " << r.second << std::endl;
// 5. Move cell r to Accepted. Update AcceptedFront.
const int rcell = r.second;
is_accepted_[rcell] = true;
solution[rcell] = r.first;
popConsidered();
accepted_front_.insert(rcell);
for (auto it = accepted_front_.begin(); it != accepted_front_.end();) {
// Note that loop increment happens in the body of this loop.
const int cell = *it;
bool on_front = false;
for (auto it2 = cell_neighbours_[cell].begin(); it2 != cell_neighbours_[cell].end(); ++it2) {
if (!is_accepted_[*it2]) {
on_front = true;
break;
}
}
if (!on_front) {
accepted_front_.erase(it++);
} else {
++it;
}
}
// 6. Recompute the value for all Considered cells within
// distance h * F_2/F1 from x_r. Use min of previous and new.
for (auto it = considered_.begin(); it != considered_.end(); ++it) {
const int ccell = it->second;
if (isClose(rcell, ccell)) {
const double value = computeValueUpdate(ccell, metric, solution.data(), rcell);
if (value < it->first) {
// Update value for considered cell.
// Note that as solution values decrease, their
// goodness w.r.t. the heap comparator increase,
// therefore we may safely call the increase()
// modificator below.
considered_.increase(considered_handles_[ccell], std::make_pair(value, ccell));
}
}
}
// 7. Move cells adjacent to r from Far to Considered.
for (auto it = cell_neighbours_[rcell].begin(); it != cell_neighbours_[rcell].end(); ++it) {
const int nb_cell = *it;
if (!is_accepted_[nb_cell] && !is_considered_[nb_cell]) {
assert(solution[nb_cell] == inf);
const double value = computeValue(nb_cell, metric, solution.data());
pushConsidered(std::make_pair(value, nb_cell));
}
}
// 8. If Considered is not empty, go to step 4.
}
}
bool AnisotropicEikonal2d::isClose(const int c1,
const int c2) const
{
const double* v[] = { grid_.cell_centroids + 2*c1,
grid_.cell_centroids + 2*c2 };
return distanceIso(v[0], v[1]) < safety_factor_ * aniso_ratio_[c1] * grid_radius_[c1];
}
double AnisotropicEikonal2d::computeValue(const int cell,
const double* metric,
const double* solution) const
{
// std::cout << "++++ computeValue(), cell = " << cell << std::endl;
const auto& nbs = cell_neighbours_[cell];
const int num_nbs = nbs.size();
const double inf = 1e100;
double val = inf;
for (int ii = 0; ii < num_nbs; ++ii) {
const int n[2] = { nbs[ii], nbs[(ii+1) % num_nbs] };
if (accepted_front_.count(n[0]) && accepted_front_.count(n[1])) {
const double cand_val = computeFromTri(cell, n[0], n[1], metric, solution);
val = std::min(val, cand_val);
}
}
if (val == inf) {
// Failed to find two accepted front nodes adjacent to this,
// so we go for a single-neighbour update.
for (int ii = 0; ii < num_nbs; ++ii) {
if (accepted_front_.count(nbs[ii])) {
const double cand_val = computeFromLine(cell, nbs[ii], metric, solution);
val = std::min(val, cand_val);
}
}
}
assert(val != inf);
// std::cout << "---> " << val << std::endl;
return val;
}
double AnisotropicEikonal2d::computeValueUpdate(const int cell,
const double* metric,
const double* solution,
const int new_cell) const
{
// std::cout << "++++ computeValueUpdate(), cell = " << cell << std::endl;
const auto& nbs = cell_neighbours_[cell];
const int num_nbs = nbs.size();
const double inf = 1e100;
double val = inf;
for (int ii = 0; ii < num_nbs; ++ii) {
const int n[2] = { nbs[ii], nbs[(ii+1) % num_nbs] };
if ((n[0] == new_cell || n[1] == new_cell)
&& accepted_front_.count(n[0]) && accepted_front_.count(n[1])) {
const double cand_val = computeFromTri(cell, n[0], n[1], metric, solution);
val = std::min(val, cand_val);
}
}
if (val == inf) {
// Failed to find two accepted front nodes adjacent to this,
// so we go for a single-neighbour update.
for (int ii = 0; ii < num_nbs; ++ii) {
if (nbs[ii] == new_cell && accepted_front_.count(nbs[ii])) {
const double cand_val = computeFromLine(cell, nbs[ii], metric, solution);
val = std::min(val, cand_val);
}
}
}
// std::cout << "---> " << val << std::endl;
return val;
}
double AnisotropicEikonal2d::computeFromLine(const int cell,
const int from,
const double* metric,
const double* solution) const
{
assert(!is_accepted_[cell]);
assert(is_accepted_[from]);
// Applying the first fundamental form to compute geodesic distance.
// Using the metric of 'cell', not 'from'.
const double dist = distanceAniso(grid_.cell_centroids + 2 * cell,
grid_.cell_centroids + 2 * from,
metric + 4 * cell);
return solution[from] + dist;
}
struct DistanceDerivative
{
const double* x1;
const double* x2;
const double* x;
double u1;
double u2;
const double* g;
double operator()(const double theta) const
{
const double xt[2] = { (1-theta)*x1[0] + theta*x2[0], (1-theta)*x1[1] + theta*x2[1] };
const double a[2] = { x[0] - xt[0], x[1] - xt[1] };
const double b[2] = { x1[0] - x2[0], x1[1] - x2[1] };
const double dQdtheta = 2*(a[0]*b[0]*g[0] + a[0]*b[1]*g[1] + a[1]*b[0]*g[2] + a[1]*b[1]*g[3]);
const double val = u2 - u1 + dQdtheta/(2*distanceAniso(x, xt, g));
// std::cout << theta << " " << val << std::endl;
return val;
}
};
double AnisotropicEikonal2d::computeFromTri(const int cell,
const int n0,
const int n1,
const double* metric,
const double* solution) const
{
// std::cout << "==== cell = " << cell << " n0 = " << n0 << " n1 = " << n1 << std::endl;
assert(!is_accepted_[cell]);
assert(is_accepted_[n0]);
assert(is_accepted_[n1]);
DistanceDerivative dd;
dd.x1 = grid_.cell_centroids + 2 * n0;
dd.x2 = grid_.cell_centroids + 2 * n1;
dd.x = grid_.cell_centroids + 2 * cell;
dd.u1 = solution[n0];
dd.u2 = solution[n1];
dd.g = metric + 4 * cell;
int iter = 0;
const double theta = RegulaFalsi<ContinueOnError>::solve(dd, 0.0, 1.0, 15, 1e-8, iter);
const double xt[2] = { (1-theta)*dd.x1[0] + theta*dd.x2[0],
(1-theta)*dd.x1[1] + theta*dd.x2[1] };
const double d1 = distanceAniso(dd.x1, dd.x, dd.g) + solution[n0];
const double d2 = distanceAniso(dd.x2, dd.x, dd.g) + solution[n1];
const double dt = distanceAniso(xt, dd.x, dd.g) + (1-theta)*solution[n0] + theta*solution[n1];
return std::min(d1, std::min(d2, dt));
}
const AnisotropicEikonal2d::ValueAndCell& AnisotropicEikonal2d::topConsidered() const
{
return considered_.top();
}
void AnisotropicEikonal2d::pushConsidered(const ValueAndCell& vc)
{
HeapHandle h = considered_.push(vc);
considered_handles_[vc.second] = h;
is_considered_[vc.second] = true;
}
void AnisotropicEikonal2d::popConsidered()
{
is_considered_[considered_.top().second] = false;
considered_handles_.erase(considered_.top().second);
considered_.pop();
}
void AnisotropicEikonal2d::computeGridRadius()
{
const int num_cells = cell_neighbours_.size();
grid_radius_.resize(num_cells);
for (int cell = 0; cell < num_cells; ++cell) {
double radius = 0.0;
const double* v1 = grid_.cell_centroids + 2*cell;
const auto& nb = cell_neighbours_[cell];
for (auto it = nb.begin(); it != nb.end(); ++it) {
const double* v2 = grid_.cell_centroids + 2*(*it);
radius = std::max(radius, distanceIso(v1, v2));
}
grid_radius_[cell] = radius;
}
}
void AnisotropicEikonal2d::computeAnisoRatio(const double* metric)
{
const int num_cells = cell_neighbours_.size();
aniso_ratio_.resize(num_cells);
for (int cell = 0; cell < num_cells; ++cell) {
const double* m = metric + 4*cell;
// Find the two eigenvalues from trace and determinant.
const double t = m[0] + m[3];
const double d = m[0]*m[3] - m[1]*m[2];
const double sd = std::sqrt(t*t/4.0 - d);
const double eig[2] = { t/2.0 - sd, t/2.0 + sd };
// Anisotropy ratio is the max ratio of the eigenvalues.
aniso_ratio_[cell] = std::max(eig[0]/eig[1], eig[1]/eig[0]);
}
}
} // namespace Opm
#else // BOOST_HEAP_AVAILABLE is false
namespace {
const char* AnisotropicEikonal2derrmsg =
"\n********************************************************************************\n"
"This library has not been compiled with support for the AnisotropicEikonal2d\n"
"class, due to too old version of the boost libraries (Boost.Heap from boost\n"
"version 1.49 or newer is required.\n"
"To use this class you must recompile opm-core on a system with sufficiently new\n"
"version of the boost libraries."
"\n********************************************************************************\n";
}
namespace Opm
{
AnisotropicEikonal2d::AnisotropicEikonal2d(const UnstructuredGrid&)
{
OPM_THROW(std::logic_error, AnisotropicEikonal2derrmsg);
}
void AnisotropicEikonal2d::solve(const double*,
const std::vector<int>&,
std::vector<double>&)
{
OPM_THROW(std::logic_error, AnisotropicEikonal2derrmsg);
}
}
#endif // BOOST_HEAP_AVAILABLE

101
opm/core/tof/AnisotropicEikonal.hpp Normal file
View File

@@ -0,0 +1,101 @@
/*
Copyright 2014 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_ANISOTROPICEIKONAL_HEADER_INCLUDED
#define OPM_ANISOTROPICEIKONAL_HEADER_INCLUDED
#include <opm/core/utility/SparseTable.hpp>
#include <vector>
#include <set>
#include <map>
#include <boost/version.hpp>
#define BOOST_HEAP_AVAILABLE ((BOOST_VERSION / 100 % 1000) >= 49)
#if BOOST_HEAP_AVAILABLE
#include <boost/heap/fibonacci_heap.hpp>
#endif
struct UnstructuredGrid;
namespace Opm
{
/// A solver for the anisotropic eikonal equation:
/// \f[ || \nabla u^T M^{-1}(x) \nabla u || = 1 \qquad x \in \Omega \f]
/// where M(x) is a symmetric positive definite matrix.
/// The boundary conditions are assumed to be
/// \f[ u(x) = 0 \qquad x \in \partial\Omega \f].
class AnisotropicEikonal2d
{
public:
/// Construct solver.
/// \param[in] grid A 2d grid.
explicit AnisotropicEikonal2d(const UnstructuredGrid& grid);
/// Solve the eikonal equation.
/// \param[in] metric Array of metric tensors, M, for each cell.
/// \param[in] startcells Array of cells where u = 0 at the centroid.
/// \param[out] solution Array of solution to the eikonal equation.
void solve(const double* metric,
const std::vector<int>& startcells,
std::vector<double>& solution);
private:
#if BOOST_HEAP_AVAILABLE
// Grid and topology.
const UnstructuredGrid& grid_;
SparseTable<int> cell_neighbours_;
// Keep track of accepted cells.
std::vector<char> is_accepted_;
std::set<int> accepted_front_;
// Quantities relating to anisotropy.
std::vector<double> grid_radius_;
std::vector<double> aniso_ratio_;
const double safety_factor_;
// Keep track of considered cells.
typedef std::pair<double, int> ValueAndCell;
typedef boost::heap::compare<std::greater<ValueAndCell>> Comparator;
typedef boost::heap::fibonacci_heap<ValueAndCell, Comparator> Heap;
Heap considered_;
typedef Heap::handle_type HeapHandle;
std::map<int, HeapHandle> considered_handles_;
std::vector<char> is_considered_;
bool isClose(const int c1, const int c2) const;
double computeValue(const int cell, const double* metric, const double* solution) const;
double computeValueUpdate(const int cell, const double* metric, const double* solution, const int new_cell) const;
double computeFromLine(const int cell, const int from, const double* metric, const double* solution) const;
double computeFromTri(const int cell, const int n0, const int n1, const double* metric, const double* solution) const;
const ValueAndCell& topConsidered() const;
void pushConsidered(const ValueAndCell& vc);
void popConsidered();
void computeGridRadius();
void computeAnisoRatio(const double* metric);
#endif // BOOST_HEAP_AVAILABLE
};
} // namespace Opm
#endif // OPM_ANISOTROPICEIKONAL_HEADER_INCLUDED

97
tests/test_anisotropiceikonal.cpp Normal file
View File

@@ -0,0 +1,97 @@
/*
Copyright 2014 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/>.
*/
#include <config.h>
#if HAVE_DYNAMIC_BOOST_TEST
#define BOOST_TEST_DYN_LINK
#endif
#define NVERBOSE // to suppress our messages when throwing
#define BOOST_TEST_MODULE AnisotropicEikonalTest
#include <boost/test/unit_test.hpp>
#include <opm/core/tof/AnisotropicEikonal.hpp>
#include <opm/core/grid/GridManager.hpp>
#include <opm/core/grid.h>
#include <cmath>
using namespace Opm;
#if BOOST_HEAP_AVAILABLE
BOOST_AUTO_TEST_CASE(cartesian_2d_a)
{
const GridManager gm(2, 2);
const UnstructuredGrid& grid = *gm.c_grid();
AnisotropicEikonal2d ae(grid);
const std::vector<double> metric = {
1, 0, 0, 1,
1, 0, 0, 1,
1, 0, 0, 1,
1, 0, 0, 1
};
BOOST_REQUIRE_EQUAL(metric.size(), grid.number_of_cells*grid.dimensions*grid.dimensions);
const std::vector<int> start = { 0 };
std::vector<double> sol;
ae.solve(metric.data(), start, sol);
BOOST_REQUIRE(!sol.empty());
BOOST_CHECK_EQUAL(sol.size(), grid.number_of_cells);
std::vector<double> truth = { 0, 1, 1, std::sqrt(2) };
BOOST_CHECK_EQUAL_COLLECTIONS(sol.begin(), sol.end(), truth.begin(), truth.end());
}
BOOST_AUTO_TEST_CASE(cartesian_2d_b)
{
const GridManager gm(3, 2, 1.0, 2.0);
const UnstructuredGrid& grid = *gm.c_grid();
AnisotropicEikonal2d ae(grid);
const std::vector<double> metric = {
1, 0, 0, 1,
1, 0, 0, 1,
1, 0, 0, 1,
1, 0, 0, 1,
1, 0, 0, 1,
1, 0, 0, 1
};
BOOST_REQUIRE_EQUAL(metric.size(), grid.number_of_cells*grid.dimensions*grid.dimensions);
const std::vector<int> start = { 0 };
std::vector<double> sol;
ae.solve(metric.data(), start, sol);
BOOST_REQUIRE(!sol.empty());
BOOST_CHECK_EQUAL(sol.size(), grid.number_of_cells);
// The test below works as a regression test, but does not test
// that cell 5 is close to the truth, which is sqrt(8).
std::vector<double> expected = { 0, 1, 2, 2, std::sqrt(5), 3.0222193552572132 };
for (int cell = 0; cell < grid.number_of_cells; ++cell) {
BOOST_CHECK_CLOSE(sol[cell], expected[cell], 1e-5);
}
}
#else // BOOST_HEAP_AVAILABLE is false
BOOST_AUTO_TEST_CASE(dummy)
{
BOOST_CHECK(true);
}
#endif // BOOST_HEAP_AVAILABLE