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 remote-tracking branch 'upstream/master'

Browse Source
This commit is contained in:
Markus Blatt 2013-06-18 10:48:28 +02:00
commit 914844055b
10 changed files with 350 additions and 201 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
opm/core/fluid/BlackoilPropertiesBasic.hpp
View File

@ -25,6 +25,7 @@
#include <opm/core/fluid/RockBasic.hpp>
#include <opm/core/fluid/PvtPropertiesBasic.hpp>
#include <opm/core/fluid/SaturationPropsBasic.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
namespace Opm
{

9
opm/core/simulator/WellState.hpp
View File

@ -42,12 +42,15 @@ namespace Opm
// to pressure in first perforation cell.
for (int w = 0; w < nw; ++w) {
const WellControls* ctrl = wells->ctrls[w];
if (ctrl->type[ctrl->current] == BHP) {
bhp_[w] = ctrl->target[ctrl->current];
} else {
if ((ctrl->current < 0) || // SHUT
(ctrl->type[ctrl->current] != BHP)) {
const int cell = wells->well_cells[wells->well_connpos[w]];
bhp_[w] = state.pressure()[cell];
}
else {
bhp_[w] = ctrl->target[ctrl->current];
}
}
perfrates_.resize(wells->well_connpos[nw], 0.0);
perfpress_.resize(wells->well_connpos[nw], -1e100);

35
opm/core/transport/reorder/DGBasis.cpp
View File

@ -25,12 +25,30 @@
namespace Opm
{
// ---------------- Methods for class DGBasisInterface ----------------
/// Virtual destructor.
DGBasisInterface::~DGBasisInterface()
{
}
/// Evaluate function f = sum_i c_i b_i at the point x.
/// Note that this function is not virtual, but implemented in
/// terms of the virtual functions of the class.
/// \param[in] cell Cell index
/// \param[in] coefficients Coefficients {c_i} for a single cell.
/// \param[in] x Point at which to compute f(x).
double DGBasisInterface::evalFunc(const int cell,
const double* coefficients,
const double* x) const
{
bvals_.resize(numBasisFunc());
eval(cell, x, &bvals_[0]);
return std::inner_product(bvals_.begin(), bvals_.end(), coefficients, 0.0);
}
// ---------------- Methods for class DGBasisBoundedTotalDegree ----------------
@ -154,6 +172,12 @@ namespace Opm
}
}
/// Compute the average of the function f = sum_i c_i b_i.
/// \param[in] coefficients Coefficients {c_i} for a single cell.
double DGBasisBoundedTotalDegree::functionAverage(const double* coefficients) const
{
return coefficients[0];
}
@ -300,11 +324,18 @@ namespace Opm
double* coefficients) const
{
const int nb = numBasisFunc();
const double average = std::accumulate(coefficients, coefficients + nb, 0.0)/double(nb);
const double aver = functionAverage(coefficients);
for (int ix = 0; ix < nb; ++ix) {
coefficients[ix] = factor*(coefficients[ix] - average) + average;
coefficients[ix] = factor*(coefficients[ix] - aver) + aver;
}
}
/// Compute the average of the function f = sum_i c_i b_i.
/// \param[in] coefficients Coefficients {c_i} for a single cell.
double DGBasisMultilin::functionAverage(const double* coefficients) const
{
const int nb = numBasisFunc();
return std::accumulate(coefficients, coefficients + nb, 0.0)/double(nb);
}
} // namespace Opm

27
opm/core/transport/reorder/DGBasis.hpp
View File

@ -20,6 +20,8 @@
#ifndef OPM_DGBASIS_HEADER_INCLUDED
#define OPM_DGBASIS_HEADER_INCLUDED
#include <vector>
struct UnstructuredGrid;
namespace Opm
@ -75,6 +77,23 @@ namespace Opm
/// \param[out] coefficients Coefficients {c_i} for a single cell.
virtual void multiplyGradient(const double factor,
double* coefficients) const = 0;
/// Evaluate function f = sum_i c_i b_i at the point x.
/// Note that this function is not virtual, but implemented in
/// terms of the virtual functions of the class.
/// \param[in] cell Cell index
/// \param[in] coefficients Coefficients {c_i} for a single cell.
/// \param[in] x Point at which to compute f(x).
double evalFunc(const int cell,
const double* coefficients,
const double* x) const;
/// Compute the average of the function f = sum_i c_i b_i.
/// \param[in] coefficients Coefficients {c_i} for a single cell.
virtual double functionAverage(const double* coefficients) const = 0;
private:
mutable std::vector<double> bvals_; // For evalFunc().
};
@ -144,6 +163,10 @@ namespace Opm
virtual void multiplyGradient(const double factor,
double* coefficients) const;
/// Compute the average of the function f = sum_i c_i b_i.
/// \param[in] coefficients Coefficients {c_i} for a single cell.
virtual double functionAverage(const double* coefficients) const;
private:
const UnstructuredGrid& grid_;
const int degree_;
@ -217,6 +240,10 @@ namespace Opm
virtual void multiplyGradient(const double factor,
double* coefficients) const;
/// Compute the average of the function f = sum_i c_i b_i.
/// \param[in] coefficients Coefficients {c_i} for a single cell.
virtual double functionAverage(const double* coefficients) const;
private:
const UnstructuredGrid& grid_;
const int degree_;

231
opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp
View File

@ -411,10 +411,10 @@ namespace Opm
{
switch (limiter_method_) {
case MinUpwindFace:
applyMinUpwindFaceLimiter(cell, tof);
applyMinUpwindLimiter(cell, true, tof);
break;
case MinUpwindAverage:
applyMinUpwindAverageLimiter(cell, tof);
applyMinUpwindLimiter(cell, false, tof);
break;
default:
THROW("Limiter type not implemented: " << limiter_method_);
@ -424,48 +424,32 @@ namespace Opm
void TransportModelTracerTofDiscGal::applyMinUpwindFaceLimiter(const int cell, double* tof)
void TransportModelTracerTofDiscGal::applyMinUpwindLimiter(const int cell, const bool face_min, double* tof)
{
if (basis_func_->degree() != 1) {
THROW("This limiter only makes sense for our DG1 implementation.");
}
// Limiter principles:
// 1. Let M be the minimum TOF value on the upstream faces,
// evaluated in the upstream cells. Then the value at all
// points in this cell shall be at least M.
// Upstream faces whose flux does not exceed the relative
// flux threshold are not considered for this minimum.
// 1. Let M be either:
// - the minimum TOF value of all upstream faces,
// evaluated in the upstream cells
// (chosen if face_min is true).
// or:
// - the minimum average TOF value of all upstream cells
// (chosen if face_min is false).
// Then the value at all points in this cell shall be at
// least M. Upstream faces whose flux does not exceed the
// relative flux threshold are not considered for this
// minimum.
// 2. The TOF shall not be below zero in any point.
// Find total upstream/downstream fluxes.
double upstream_flux = 0.0;
double downstream_flux = 0.0;
for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) {
const int face = grid_.cell_faces[hface];
double flux = 0.0;
if (cell == grid_.face_cells[2*face]) {
flux = darcyflux_[face];
} else {
flux = -darcyflux_[face];
}
if (flux < 0.0) {
upstream_flux += flux;
} else {
downstream_flux += flux;
}
}
// In the presence of sources, significant fluxes may be missing from the computed fluxes,
// setting the total flux to the (positive) maximum avoids this: since source is either
// inflow or outflow, not both, either upstream_flux or downstream_flux must be correct.
const double total_flux = std::max(-upstream_flux, downstream_flux);
// Find minimum tof on upstream faces and for this cell.
const int dim = grid_.dimensions;
// Find minimum tof on upstream faces/cells and for this cell.
const int num_basis = basis_func_->numBasisFunc();
double min_upstream_tof = 1e100;
double min_here_tof = 1e100;
int num_upstream_faces = 0;
const double total_flux = totalFlux(cell);
for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) {
const int face = grid_.cell_faces[hface];
double flux = 0.0;
@ -478,30 +462,22 @@ namespace Opm
upstream_cell = grid_.face_cells[2*face];
}
const bool upstream = (flux < -total_flux*limiter_relative_flux_threshold_);
const bool interior = (upstream_cell >= 0);
// Find minimum tof in this cell and upstream.
// The meaning of minimum upstream tof depends on method.
min_here_tof = std::min(min_here_tof, minCornerVal(cell, face));
if (upstream) {
++num_upstream_faces;
}
bool interior = (upstream_cell >= 0);
// Evaluate the solution in all corners.
for (int fnode = grid_.face_nodepos[face]; fnode < grid_.face_nodepos[face+1]; ++fnode) {
const double* nc = grid_.node_coordinates + dim*grid_.face_nodes[fnode];
basis_func_->eval(cell, nc, &basis_[0]);
const double tof_here = std::inner_product(basis_.begin(), basis_.end(),
tof_coeff_ + num_basis*cell, 0.0);
min_here_tof = std::min(min_here_tof, tof_here);
if (upstream) {
if (interior) {
basis_func_->eval(upstream_cell, nc, &basis_nb_[0]);
const double tof_upstream
= std::inner_product(basis_nb_.begin(), basis_nb_.end(),
tof_coeff_ + num_basis*upstream_cell, 0.0);
min_upstream_tof = std::min(min_upstream_tof, tof_upstream);
double upstream_tof = 0.0;
if (interior) {
if (face_min) {
upstream_tof = minCornerVal(upstream_cell, face);
} else {
// Allow tof down to 0 on inflow boundaries.
min_upstream_tof = std::min(min_upstream_tof, 0.0);
upstream_tof = basis_func_->functionAverage(tof_coeff_ + num_basis*upstream_cell);
}
}
min_upstream_tof = std::min(min_upstream_tof, upstream_tof);
}
}
@ -513,7 +489,7 @@ namespace Opm
if (min_upstream_tof < 0.0) {
min_upstream_tof = 0.0;
}
const double tof_c = tof_coeff_[num_basis*cell];
const double tof_c = basis_func_->functionAverage(tof_coeff_ + num_basis*cell);
double limiter = (tof_c - min_upstream_tof)/(tof_c - min_here_tof);
if (tof_c < min_upstream_tof) {
// Handle by setting a flat solution.
@ -535,110 +511,6 @@ namespace Opm
void TransportModelTracerTofDiscGal::applyMinUpwindAverageLimiter(const int cell, double* tof)
{
if (basis_func_->degree() != 1) {
THROW("This limiter only makes sense for our DG1 implementation.");
}
// Limiter principles:
// 1. Let M be the average TOF value of the upstream cells.
/// Then the value at all points in this cell shall be at least M.
// Upstream faces whose flux does not exceed the relative
// flux threshold are not considered for this minimum.
// 2. The TOF shall not be below zero in any point.
// Find total upstream/downstream fluxes.
double upstream_flux = 0.0;
double downstream_flux = 0.0;
for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) {
const int face = grid_.cell_faces[hface];
double flux = 0.0;
if (cell == grid_.face_cells[2*face]) {
flux = darcyflux_[face];
} else {
flux = -darcyflux_[face];
}
if (flux < 0.0) {
upstream_flux += flux;
} else {
downstream_flux += flux;
}
}
// In the presence of sources, significant fluxes may be missing from the computed fluxes,
// setting the total flux to the (positive) maximum avoids this: since source is either
// inflow or outflow, not both, either upstream_flux or downstream_flux must be correct.
const double total_flux = std::max(-upstream_flux, downstream_flux);
// Find minimum tof on upstream faces and for this cell.
const int dim = grid_.dimensions;
const int num_basis = basis_func_->numBasisFunc();
double min_upstream_tof = 1e100;
double min_here_tof = 1e100;
int num_upstream_faces = 0;
for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) {
const int face = grid_.cell_faces[hface];
double flux = 0.0;
int upstream_cell = -1;
if (cell == grid_.face_cells[2*face]) {
flux = darcyflux_[face];
upstream_cell = grid_.face_cells[2*face+1];
} else {
flux = -darcyflux_[face];
upstream_cell = grid_.face_cells[2*face];
}
const bool upstream = (flux < -total_flux*limiter_relative_flux_threshold_);
if (upstream) {
++num_upstream_faces;
}
bool interior = (upstream_cell >= 0);
// Evaluate the solution in all corners.
for (int fnode = grid_.face_nodepos[face]; fnode < grid_.face_nodepos[face+1]; ++fnode) {
const double* nc = grid_.node_coordinates + dim*grid_.face_nodes[fnode];
basis_func_->eval(cell, nc, &basis_[0]);
const double tof_here = std::inner_product(basis_.begin(), basis_.end(),
tof_coeff_ + num_basis*cell, 0.0);
min_here_tof = std::min(min_here_tof, tof_here);
if (upstream) {
if (interior) {
min_upstream_tof = std::min(min_upstream_tof, tof_coeff_[num_basis*upstream_cell]);
} else {
// Allow tof down to 0 on inflow boundaries.
min_upstream_tof = std::min(min_upstream_tof, 0.0);
}
}
}
}
// Compute slope multiplier (limiter).
if (num_upstream_faces == 0) {
min_upstream_tof = 0.0;
min_here_tof = 0.0;
}
if (min_upstream_tof < 0.0) {
min_upstream_tof = 0.0;
}
const double tof_c = tof_coeff_[num_basis*cell];
double limiter = (tof_c - min_upstream_tof)/(tof_c - min_here_tof);
if (tof_c < min_upstream_tof) {
// Handle by setting a flat solution.
std::cout << "Trouble in cell " << cell << std::endl;
limiter = 0.0;
basis_func_->addConstant(min_upstream_tof - tof_c, tof + num_basis*cell);
}
ASSERT(limiter >= 0.0);
// Actually do the limiting (if applicable).
if (limiter < 1.0) {
// std::cout << "Applying limiter in cell " << cell << ", limiter = " << limiter << std::endl;
basis_func_->multiplyGradient(limiter, tof + num_basis*cell);
} else {
// std::cout << "Not applying limiter in cell " << cell << "!" << std::endl;
}
}
void TransportModelTracerTofDiscGal::applyLimiterAsPostProcess()
@ -675,4 +547,51 @@ namespace Opm
double TransportModelTracerTofDiscGal::totalFlux(const int cell) const
{
// Find total upstream/downstream fluxes.
double upstream_flux = 0.0;
double downstream_flux = 0.0;
for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) {
const int face = grid_.cell_faces[hface];
double flux = 0.0;
if (cell == grid_.face_cells[2*face]) {
flux = darcyflux_[face];
} else {
flux = -darcyflux_[face];
}
if (flux < 0.0) {
upstream_flux += flux;
} else {
downstream_flux += flux;
}
}
// In the presence of sources, significant fluxes may be missing from the computed fluxes,
// setting the total flux to the (positive) maximum avoids this: since source is either
// inflow or outflow, not both, either upstream_flux or downstream_flux must be correct.
return std::max(-upstream_flux, downstream_flux);
}
double TransportModelTracerTofDiscGal::minCornerVal(const int cell, const int face) const
{
// Evaluate the solution in all corners.
const int dim = grid_.dimensions;
const int num_basis = basis_func_->numBasisFunc();
double min_cornerval = 1e100;
for (int fnode = grid_.face_nodepos[face]; fnode < grid_.face_nodepos[face+1]; ++fnode) {
const double* nc = grid_.node_coordinates + dim*grid_.face_nodes[fnode];
basis_func_->eval(cell, nc, &basis_[0]);
const double tof_corner = std::inner_product(basis_.begin(), basis_.end(),
tof_coeff_ + num_basis*cell, 0.0);
min_cornerval = std::min(min_cornerval, tof_corner);
}
return min_cornerval;
}
} // namespace Opm

9
opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp
View File

@ -119,8 +119,8 @@ namespace Opm
std::vector<double> orig_jac_; // single-cell jacobian (copy)
// Below: storage for quantities needed by solveSingleCell().
std::vector<double> coord_;
std::vector<double> basis_;
std::vector<double> basis_nb_;
mutable std::vector<double> basis_;
mutable std::vector<double> basis_nb_;
std::vector<double> grad_basis_;
std::vector<double> velocity_;
@ -130,10 +130,11 @@ namespace Opm
// (will read data from tof_coeff_, it is ok to call
// with tof_coeff as tof argument.
void applyLimiter(const int cell, double* tof);
void applyMinUpwindFaceLimiter(const int cell, double* tof);
void applyMinUpwindAverageLimiter(const int cell, double* tof);
void applyMinUpwindLimiter(const int cell, const bool face_min, double* tof);
void applyLimiterAsPostProcess();
void applyLimiterAsSimultaneousPostProcess();
double totalFlux(const int cell) const;
double minCornerVal(const int cell, const int face) const;
};
} // namespace Opm

26
opm/core/transport/reorder/tarjan.c
View File

@ -20,8 +20,8 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
*/
#include <string.h>
#include <assert.h>
#include <stddef.h>
#ifdef MATLAB_MEX_FILE
#include "tarjan.h"
@ -30,7 +30,13 @@ SOFTWARE.
#endif
static void
clear_vector(size_t n, int *v)
{
size_t i;
for (i = 0; i < n; i++) { v[i] = 0; }
}
static int min(int a, int b){ return a < b? a : b;}
@ -73,19 +79,21 @@ tarjan (int nv, const int *ia, const int *ja, int *vert, int *comp,
int *stack = comp + nv;
int *bottom = stack;
int *cstack = vert + nv-1;
#if !defined(NDEBUG)
int *cbottom = cstack;
#endif
int t = 0;
int pos = 0;
int *time = work;
int *link = (int *) time + nv;
int *status = (int *) link + nv; /* dual usage... */
int *link = time + nv;
int *status = link + nv; /* dual usage... */
(void) cbottom;
memset(work, 0, 3*nv * sizeof *work);
memset(vert, 0, nv * sizeof *vert );
memset(comp, 0, (nv+1) * sizeof *comp );
clear_vector(3 * ((size_t) nv), work);
clear_vector(1 * ((size_t) nv), vert);
clear_vector(1 + ((size_t) nv), comp);
/* Init status all vertices */
for (i=0; i<nv; ++i)
@ -188,7 +196,7 @@ tarjan (int nv, const int *ia, const int *ja, int *vert, int *comp,
}
else
{
assert(status[child] = DONE);
assert(status[child] == DONE);
}
}
}

27
opm/core/wells/WellCollection.cpp
View File

@ -1,22 +1,21 @@
/*
Copyright 2011 SINTEF ICT, Applied Mathematics.
Copyright 2011 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
This file is part of The Open Reservoir Simulator Project (OpenRS).
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.
OpenRS 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.
OpenRS 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 OpenRS. If not, see <http://www.gnu.org/licenses/>.
*/
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/wells/WellCollection.hpp>

26
opm/core/wells/WellCollection.hpp
View File

@ -1,22 +1,22 @@
/*
Copyright 2011 SINTEF ICT, Applied Mathematics.
Copyright 2011 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
This file is part of The Open Reservoir Simulator Project (OpenRS).
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.
OpenRS 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.
OpenRS 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/>.
*/
You should have received a copy of the GNU General Public License
along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_WELLCOLLECTION_HPP

160
tests/test_dgbasis.cpp Normal file
View File

@ -0,0 +1,160 @@
/*
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/>.
*/
#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 DGBasisTest
#include <boost/test/unit_test.hpp>
#include <opm/core/transport/reorder/DGBasis.hpp>
#include <opm/core/GridManager.hpp>
#include <opm/core/grid.h>
#include <cmath>
using namespace Opm;
namespace
{
bool aequal(double a, double b)
{
const double eps = 1e-15;
return std::fabs(a - b) < eps;
}
} // anonymous namespace
namespace cart2d
{
static void test()
{
// Set up 2d 1-cell cartesian case.
GridManager g(1, 1);
const UnstructuredGrid& grid = *g.c_grid();
// Test DGBasisBoundedTotalDegree, degree 0.
{
DGBasisBoundedTotalDegree b(grid, 0);
BOOST_CHECK_EQUAL(b.numBasisFunc(), 1);
std::vector<double> bx(b.numBasisFunc(), 0.0);
b.eval(0, grid.cell_centroids, &bx[0]);
BOOST_CHECK(aequal(bx[0], 1.0));
double x[2] = { 0.123, 0.456 };
b.eval(0, x, &bx[0]);
BOOST_CHECK(aequal(bx[0], 1.0));
std::vector<double> c(b.numBasisFunc(), 0.0);
b.addConstant(0.789, &c[0]);
BOOST_CHECK(aequal(c[0], 0.789));
b.multiplyGradient(1.234, &c[0]);
BOOST_CHECK(aequal(c[0], 0.789));
}
// Test DGBasisBoundedTotalDegree, degree 1.
{
DGBasisBoundedTotalDegree b(grid, 1);
BOOST_CHECK_EQUAL(b.numBasisFunc(), 3);
std::vector<double> bx(b.numBasisFunc(), 0.0);
b.eval(0, grid.cell_centroids, &bx[0]);
BOOST_CHECK(aequal(bx[0], 1.0));
BOOST_CHECK(aequal(bx[1], 0.0));
BOOST_CHECK(aequal(bx[2], 0.0));
double x[2] = { 0.123, 0.456 };
b.eval(0, x, &bx[0]);
BOOST_CHECK(aequal(bx[0], 1.0));
BOOST_CHECK(aequal(bx[1], 0.123 - 0.5));
BOOST_CHECK(aequal(bx[2], 0.456 - 0.5));
std::vector<double> c(b.numBasisFunc(), 0.0);
c[0] = 1.0; c[1] = 2.0; c[2] = 3.0;
b.addConstant(0.789, &c[0]);
BOOST_CHECK(aequal(c[0], 1.789));
BOOST_CHECK(aequal(c[1], 2.0));
BOOST_CHECK(aequal(c[2], 3.0));
const double fx = c[0]*bx[0] + c[1]*bx[1] + c[2]*bx[2];
b.multiplyGradient(1.234, &c[0]);
BOOST_CHECK(aequal(c[0], 1.789));
BOOST_CHECK(aequal(c[1], 2.0*1.234));
BOOST_CHECK(aequal(c[2], 3.0*1.234));
const double fx2 = c[0]*bx[0] + c[1]*bx[1] + c[2]*bx[2];
BOOST_CHECK(aequal(fx2 - c[0], 1.234*(fx - c[0])));
}
// Test DGBasisMultilin, degree 0.
{
DGBasisMultilin b(grid, 0);
BOOST_CHECK_EQUAL(b.numBasisFunc(), 1);
std::vector<double> bx(b.numBasisFunc(), 0.0);
b.eval(0, grid.cell_centroids, &bx[0]);
BOOST_CHECK(aequal(bx[0], 1.0));
double x[2] = { 0.123, 0.456 };
b.eval(0, x, &bx[0]);
BOOST_CHECK(aequal(bx[0], 1.0));
std::vector<double> c(b.numBasisFunc(), 0.0);
b.addConstant(0.789, &c[0]);
BOOST_CHECK(aequal(c[0], 0.789));
b.multiplyGradient(1.234, &c[0]);
BOOST_CHECK(aequal(c[0], 0.789));
}
// Test DGBasisMultilin, degree 1.
{
DGBasisMultilin b(grid, 1);
BOOST_CHECK_EQUAL(b.numBasisFunc(), 4);
std::vector<double> bx(b.numBasisFunc(), 0.0);
b.eval(0, grid.cell_centroids, &bx[0]);
BOOST_CHECK(aequal(bx[0], 0.25));
BOOST_CHECK(aequal(bx[1], 0.25));
BOOST_CHECK(aequal(bx[2], 0.25));
BOOST_CHECK(aequal(bx[3], 0.25));
double x[2] = { 0.123, 0.456 };
b.eval(0, x, &bx[0]);
const double xm[2] = { 1.0 - x[0], x[0] };
const double ym[2] = { 1.0 - x[1], x[1] };
BOOST_CHECK(aequal(bx[0], xm[0]*ym[0]));
BOOST_CHECK(aequal(bx[1], xm[0]*ym[1]));
BOOST_CHECK(aequal(bx[2], xm[1]*ym[0]));
BOOST_CHECK(aequal(bx[3], xm[1]*ym[1]));
std::vector<double> c(b.numBasisFunc(), 0.0);
c[0] = -1.567; c[1] = 1.42; c[2] = 0.59; c[3] = 3.225;
std::vector<double> corig = c;
b.addConstant(0.789, &c[0]);
BOOST_CHECK(aequal(c[0], corig[0] + 0.25*0.789));
BOOST_CHECK(aequal(c[1], corig[1] + 0.25*0.789));
BOOST_CHECK(aequal(c[2], corig[2] + 0.25*0.789));
BOOST_CHECK(aequal(c[3], corig[3] + 0.25*0.789));
const double fx = c[0]*bx[0] + c[1]*bx[1] + c[2]*bx[2] + c[3]*bx[3];
const double fc = 0.25*(c[0] + c[1] + c[2] + c[3]);
b.multiplyGradient(1.234, &c[0]);
const double fx2 = c[0]*bx[0] + c[1]*bx[1] + c[2]*bx[2] + c[3]*bx[3];
BOOST_CHECK(aequal(fx2 - fc, 1.234*(fx - fc)));
}
}
} // namespace cart2d
BOOST_AUTO_TEST_CASE(test_dgbasis)
{
cart2d::test();
}