opm-simulators/tests/test_vfpproperties.cpp

/*
  Copyright 2015 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>

#define BOOST_TEST_MODULE VFPTest

#include <algorithm>
#include <memory>
#include <map>
#include <sstream>
#include <limits>
#include <vector>

#include <opm/common/utility/platform_dependent/disable_warnings.h>
#include <boost/test/unit_test.hpp>
#include <boost/filesystem.hpp>
#include <opm/common/utility/platform_dependent/reenable_warnings.h>

#include <opm/core/wells.h>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/EclipseState/checkDeck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/VFPProdTable.hpp>
#include <opm/parser/eclipse/Units/UnitSystem.hpp>

#include <opm/simulators/wells/VFPHelpers.hpp>
#include <opm/simulators/wells/VFPProdProperties.hpp>



const double max_d_tol = 1.0e-10;
const double sad_tol = 1.0e-8;










BOOST_AUTO_TEST_SUITE( HelperTests )

BOOST_AUTO_TEST_CASE(findInterpData)
{
    std::vector<double> values = {1, 5, 7, 9, 11, 15};
    double exact = 9.0;
    double interpolate = 6.0;
    double extrapolate_left = -1.0;
    double extrapolate_right = 19;
    double first = 1;
    double last = 15;

    Opm::detail::InterpData eval0 = Opm::detail::findInterpData(exact, values);
    Opm::detail::InterpData eval1 = Opm::detail::findInterpData(interpolate, values);
    Opm::detail::InterpData eval2 = Opm::detail::findInterpData(extrapolate_left, values);
    Opm::detail::InterpData eval3 = Opm::detail::findInterpData(extrapolate_right, values);
    Opm::detail::InterpData eval4 = Opm::detail::findInterpData(first, values);
    Opm::detail::InterpData eval5 = Opm::detail::findInterpData(last, values);

    BOOST_CHECK_EQUAL(eval0.ind_[0], 2);
    BOOST_CHECK_EQUAL(eval0.ind_[1], 3);
    BOOST_CHECK_EQUAL(eval0.factor_, 1.0);

    BOOST_CHECK_EQUAL(eval1.ind_[0], 1);
    BOOST_CHECK_EQUAL(eval1.ind_[1], 2);
    BOOST_CHECK_EQUAL(eval1.factor_, 0.5);

    BOOST_CHECK_EQUAL(eval2.ind_[0], 0);
    BOOST_CHECK_EQUAL(eval2.ind_[1], 1);
    BOOST_CHECK_EQUAL(eval2.factor_, -0.25);

    BOOST_CHECK_EQUAL(eval3.ind_[0], 4);
    BOOST_CHECK_EQUAL(eval3.ind_[1], 5);
    BOOST_CHECK_EQUAL(eval3.factor_, 2.0);

    BOOST_CHECK_EQUAL(eval4.ind_[0], 0);
    BOOST_CHECK_EQUAL(eval4.ind_[1], 1);
    BOOST_CHECK_EQUAL(eval4.factor_, 0.0);

    BOOST_CHECK_EQUAL(eval5.ind_[0], 4);
    BOOST_CHECK_EQUAL(eval5.ind_[1], 5);
    BOOST_CHECK_EQUAL(eval5.factor_, 1.0);
}

BOOST_AUTO_TEST_SUITE_END() // HelperTests


/**
 * Test fixture to set up axis etc.
 * All of our axes go from 0 to 1, but with a varying number of
 * values data is given at
 */
struct TrivialFixture {
    typedef Opm::detail::VFPEvaluation VFPEvaluation;

    TrivialFixture() : table_ids(1, 1),
            thp_axis{0.0, 1.0},
            wfr_axis{0.0, 0.5, 1.0},
            gfr_axis{0.0, 0.25, 0.5, 0.75, 1},
            alq_axis{0.0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1},
            flo_axis{0.0, 0.0625, 0.125, 0.1875, 0.25, 0.3125, 0.375, 0.4375,
                0.5, 0.5625, 0.625, 0.6875, 0.75, 0.8125, 0.875, 0.9375, 1},
            nx(thp_axis.size()),
            ny(wfr_axis.size()),
            nz(gfr_axis.size()),
            nu(alq_axis.size()),
            nv(flo_axis.size()),
            size{{ nx, ny, nz, nu, nv }},
            data(size)
    {

    }

    ~TrivialFixture() {

    }

    /**
     * Fills our interpolation data with zeros
     */
    inline void fillData(double value) {
        for (int i=0; i<nx; ++i) {
            for (int j=0; j<ny; ++j) {
                for (int k=0; k<nz; ++k) {
                    for (int l=0; l<nu; ++l) {
                        for (int m=0; m<nv; ++m) {
                            data[i][j][k][l][m] = value;
                        }
                    }
                }
            }
        }
    }

    /**
     * Fills our interpolation data with an ND plane
     */
    inline void fillDataPlane() {
        for (int i=0; i<nx; ++i) {
            double x = i / static_cast<double>(nx-1);
            for (int j=0; j<ny; ++j) {
                double y = j / static_cast<double>(ny-1);
                for (int k=0; k<nz; ++k) {
                    double z = k / static_cast<double>(nz-1);
                    for (int l=0; l<nu; ++l) {
                        double u = l / static_cast<double>(nu-1);
                        for (int m=0; m<nv; ++m) {
                            double v = m / static_cast<double>(nv-1);
                            // table[thp_idx][wfr_idx][gfr_idx][alq_idx][flo_idx];
                            data[i][j][k][l][m] = x + 2*y + 3*z + 4*u + 5*v;
                        }
                    }
                }
            }
        }
    }



    /**
     * Fills our interpolation data with "random" values
     */
    inline void fillDataRandom() {
        unsigned long randx = 42;
        static double max_val = static_cast<double>(std::numeric_limits<unsigned long>::max());
        for (int i=0; i<nx; ++i) {
            for (int j=0; j<ny; ++j) {
                for (int k=0; k<nz; ++k) {
                    for (int l=0; l<nu; ++l) {
                        for (int m=0; m<nv; ++m) {
                            data[i][j][k][l][m] = randx / max_val;
                            randx = (randx*1103515245 + 12345);
                        }
                    }
                }
            }
        }
    }


    inline void initProperties() {
        table.reset(new Opm::VFPProdTable(1,
                                          1000.0,
                                          Opm::VFPProdTable::FLO_OIL,
                                          Opm::VFPProdTable::WFR_WOR,
                                          Opm::VFPProdTable::GFR_GOR,
                                          Opm::VFPProdTable::ALQ_UNDEF,
                                          flo_axis,
                                          thp_axis,
                                          wfr_axis,
                                          gfr_axis,
                                          alq_axis,
                                          data));

        properties.reset(new Opm::VFPProdProperties(table.get()));
    }



    std::shared_ptr<Opm::VFPProdProperties> properties;
    std::shared_ptr<Opm::VFPProdTable> table;
    std::vector<int> table_ids;

private:
    const std::vector<double> thp_axis;
    const std::vector<double> wfr_axis;
    const std::vector<double> gfr_axis;
    const std::vector<double> alq_axis;
    const std::vector<double> flo_axis;
    int nx;
    int ny;
    int nz;
    int nu;
    int nv;
    Opm::VFPProdTable::extents size;
    Opm::VFPProdTable::array_type data;
};





//Set F to be our test suite fixture for our "trivial" tests
BOOST_FIXTURE_TEST_SUITE( TrivialTests, TrivialFixture )



/**
 * Test that we can generate some dummy zero-data,
 * interpolate using doubles as input, and compare against the analytic solution
 */
BOOST_AUTO_TEST_CASE(InterpolateZero)
{
    fillData(0.0);
    initProperties();

    //Check interpolation
    double sum = 0.0;
    int n=5;
    for (int i=1; i<n; ++i) {
        const double x = i / static_cast<double>(n-1);
        for (int j=1; j<n; ++j) {
            const double y = j / static_cast<double>(n-1);
            for (int k=1; k<n; ++k) {
                const double z = k / static_cast<double>(n-1);
                for (int l=0; l<n; ++l) {
                    const double u = l / static_cast<double>(n-1);
                    for (int m=0; m<n; ++m) {
                        const double v = m / static_cast<double>(n-1);

                        //Note order of arguments!
                        sum += properties->bhp(1, v, x, y, z, u);
                    }
                }
            }
        }
    }

    BOOST_CHECK_EQUAL(sum, 0.0);
}


/**
 * Test that we can generate some dummy one-data,
 * interpolate using doubles as input, and compare against the analytic solution
 */
BOOST_AUTO_TEST_CASE(InterpolateOne)
{
    fillData(1.0);
    initProperties();

    //Check interpolation
    double sum = 0.0;
    int n=5;
    for (int i=1; i<n; ++i) {
        const double x = i / static_cast<double>(n-1);
        for (int j=1; j<n; ++j) {
            const double y = j / static_cast<double>(n-1);
            for (int k=1; k<n; ++k) {
                const double z = k / static_cast<double>(n-1);
                for (int l=0; l<n; ++l) {
                    const double u = l / static_cast<double>(n-1);
                    for (int m=0; m<n; ++m) {
                        const double v = m / static_cast<double>(n-1);

                        //Note order of arguments!
                        const double value = properties->bhp(1, v, x, y, z, u);

                        sum += value;
                    }
                }
            }
        }
    }

    double reference = (n-1)*(n-1)*(n-1)*n*n;
    BOOST_CHECK_EQUAL(sum, reference);
}


/**
 * Test that we can generate some dummy data representing an ND plane,
 * interpolate using doubles as input, and compare against the analytic solution
 */
BOOST_AUTO_TEST_CASE(InterpolatePlane)
{
    const int n=5;

    fillDataPlane();
    initProperties();

    //Temps used to store reference and actual variables
    double sad = 0.0;
    double max_d = 0.0;

    //Check interpolation
    for (int i=0; i<=n; ++i) {
        const double thp = i / static_cast<double>(n);
        for (int j=1; j<=n; ++j) {
            const double aqua = -j / static_cast<double>(n);
            for (int k=1; k<=n; ++k) {
                const double vapour = -k / static_cast<double>(n);
                for (int l=0; l<=n; ++l) {
                    const double alq = l / static_cast<double>(n);
                    for (int m=1; m<=n; ++m) {
                        const double liquid = -m / static_cast<double>(n);

                        //Find values that should be in table
                        double flo = Opm::detail::getFlo(aqua, liquid, vapour, table->getFloType());
                        double wfr = Opm::detail::getWFR(aqua, liquid, vapour, table->getWFRType());
                        double gfr = Opm::detail::getGFR(aqua, liquid, vapour, table->getGFRType());

                        //Calculate reference
                        double reference = thp + 2*wfr + 3*gfr+ 4*alq - 5*flo;

                        //Calculate actual
                        //Note order of arguments: id, aqua, liquid, vapour, thp, alq
                        double actual = properties->bhp(1, aqua, liquid, vapour, thp, alq);


                        double abs_diff = std::abs(actual - reference);
                        max_d = std::max(max_d, abs_diff);
                        sad = sad + abs_diff;
                    }
                }
            }
        }
    }

    BOOST_CHECK_SMALL(max_d, max_d_tol);
    BOOST_CHECK_SMALL(sad, sad_tol);
}



/**
 * Test that we can generate some dummy data representing an ND plane,
 * interpolate using doubles as input, and compare against the analytic solution
 */
BOOST_AUTO_TEST_CASE(ExtrapolatePlane)
{
    fillDataPlane();
    initProperties();

    //Check linear extrapolation (i.e., using values of x, y, etc. outside our interpolant domain)
    double sum = 0.0;
    double reference_sum = 0.0;
    double sad = 0.0; // Sum absolute difference
    double max_d = 0.0; // Maximum difference
    int n=1;
    int o=5;
    for (int i=0; i<=n+o; ++i) {
        const double x = i / static_cast<double>(n);
        for (int j=1; j<=n+o; ++j) {
            const double aqua = -j / static_cast<double>(n);
            for (int k=1; k<=n+o; ++k) {
                const double vapour = -k / static_cast<double>(n);
                for (int l=0; l<=n+o; ++l) {
                    const double u = l / static_cast<double>(n);
                    for (int m=1; m<=n+o; ++m) {
                        const double liquid = -m / static_cast<double>(n);

                        //Find values that should be in table
                        double v = Opm::detail::getFlo(aqua, liquid, vapour, table->getFloType());
                        double y = Opm::detail::getWFR(aqua, liquid, vapour, table->getWFRType());
                        double z = Opm::detail::getGFR(aqua, liquid, vapour, table->getGFRType());

                        double reference = x + 2*y + 3*z+ 4*u - 5*v;
                        reference_sum += reference;

                        //Note order of arguments! id, aqua, liquid, vapour, thp , alq
                        double value = properties->bhp(1, aqua, liquid, vapour, x, u);
                        sum += value;

                        double abs_diff = std::abs(value - reference);

                        sad += std::abs(abs_diff);
                        max_d = std::max(max_d, abs_diff);
                    }
                }
            }
        }
    }

    BOOST_CHECK_CLOSE(sum, reference_sum, 0.0001);
    BOOST_CHECK_SMALL(max_d, max_d_tol);
    BOOST_CHECK_SMALL(sad, sad_tol);
}





/**
 * Test that the partial derivatives are reasonable
 */
BOOST_AUTO_TEST_CASE(PartialDerivatives)
{
    const int n=5;

    fillDataPlane();
    initProperties();

    //Temps used to store reference and actual variables
    VFPEvaluation sad;
    VFPEvaluation max_d;

    //Check interpolation
    for (int i=0; i<=n; ++i) {
        const double thp = i / static_cast<double>(n);
        for (int j=1; j<=n; ++j) {
            const double aqua = -j / static_cast<double>(n);
            for (int k=1; k<=n; ++k) {
                const double vapour = -k / static_cast<double>(n);
                for (int l=0; l<=n; ++l) {
                    const double alq = l / static_cast<double>(n);
                    for (int m=1; m<=n; ++m) {
                        const double liquid = -m / static_cast<double>(n);

                        //Find values that should be in table
                        double flo = Opm::detail::getFlo(aqua, liquid, vapour, table->getFloType());
                        double wfr = Opm::detail::getWFR(aqua, liquid, vapour, table->getWFRType());
                        double gfr = Opm::detail::getGFR(aqua, liquid, vapour, table->getGFRType());

                        //Calculate reference
                        VFPEvaluation reference;
                        reference.value = thp + 2*wfr + 3*gfr+ 4*alq - 5*flo;
                        reference.dthp = 1;
                        reference.dwfr = 2;
                        reference.dgfr = 3;
                        reference.dalq = 4;
                        reference.dflo = 5;

                        //Calculate actual
                        //Note order of arguments: id, aqua, liquid, vapour, thp, alq
                        VFPEvaluation actual = Opm::detail::bhp(table.get(), aqua, liquid, vapour, thp, alq);

                        VFPEvaluation abs_diff = actual - reference;
                        abs_diff.value = std::abs(abs_diff.value);
                        abs_diff.dthp = std::abs(abs_diff.dthp);
                        abs_diff.dwfr = std::abs(abs_diff.dwfr);
                        abs_diff.dgfr = std::abs(abs_diff.dgfr);
                        abs_diff.dalq = std::abs(abs_diff.dalq);
                        abs_diff.dflo = std::abs(abs_diff.dflo);

                        max_d.value = std::max(max_d.value, abs_diff.value);
                        max_d.dthp = std::max(max_d.dthp, abs_diff.dthp);
                        max_d.dwfr = std::max(max_d.dwfr, abs_diff.dwfr);
                        max_d.dgfr = std::max(max_d.dgfr, abs_diff.dgfr);
                        max_d.dalq = std::max(max_d.dalq, abs_diff.dalq);
                        max_d.dflo = std::max(max_d.dflo, abs_diff.dflo);

                        sad = sad + abs_diff;
                    }
                }
            }
        }
    }

    BOOST_CHECK_SMALL(max_d.value, max_d_tol);
    BOOST_CHECK_SMALL(max_d.dthp, max_d_tol);
    BOOST_CHECK_SMALL(max_d.dwfr, max_d_tol);
    BOOST_CHECK_SMALL(max_d.dgfr, max_d_tol);
    BOOST_CHECK_SMALL(max_d.dalq, max_d_tol);
    BOOST_CHECK_SMALL(max_d.dflo, max_d_tol);

    BOOST_CHECK_SMALL(sad.value, sad_tol);
    BOOST_CHECK_SMALL(sad.dthp, sad_tol);
    BOOST_CHECK_SMALL(sad.dwfr, sad_tol);
    BOOST_CHECK_SMALL(sad.dgfr, sad_tol);
    BOOST_CHECK_SMALL(sad.dalq, sad_tol);
    BOOST_CHECK_SMALL(sad.dflo, sad_tol);
}




BOOST_AUTO_TEST_CASE(THPToBHPAndBackPlane)
{
    fillDataPlane();
    initProperties();

    double aqua = -0.5;
    double liquid = -0.9;
    double vapour = -0.1;
    double thp = 50.0;
    double alq = 32.9;

    double bhp_val = properties->bhp(1, aqua, liquid, vapour, thp, alq);
    double thp_val = properties->thp(1, aqua, liquid, vapour, bhp_val, alq);

    BOOST_CHECK_CLOSE(thp_val, thp, max_d_tol);
}



BOOST_AUTO_TEST_CASE(THPToBHPAndBackNonTrivial)
{
    fillDataRandom();
    initProperties();

    double aqua = -0.5;
    double liquid = -0.9;
    double vapour = -0.1;
    double thp = 50.0;
    double alq = 32.9;

    double bhp_val = properties->bhp(1, aqua, liquid, vapour, thp, alq);
    double thp_val = properties->thp(1, aqua, liquid, vapour, bhp_val, alq);

    BOOST_CHECK_CLOSE(thp_val, thp, max_d_tol);
}




BOOST_AUTO_TEST_SUITE_END() // Trivial tests


BOOST_AUTO_TEST_SUITE(IntegrationTests)

extern const double reference[];


/**
 * Uses a VFP table that should be linear in THP vs BHP,
 * so that
 *   bhp(aqua, liquid, vapour, thp, alq) == 0.0 for thp == 0.0,
 * and
 *   bhp(aqua, liquid, vapour, thp, alq) == 1.0 for thp == 1.0,
 */
BOOST_AUTO_TEST_CASE(ParseInterpolateLine)
{
    std::string table_str = "\
-- VFP table that is basically the identity: BHP == THP   \n\
-- In our case, we simply degenerate all axes except the THP axis \n\
-- and set that bhp(flo, thp, wfr, gfr, alq) = 0 for thp = 0, and 1 for thp = 1. \n\
-- The value of flo, wfr, gfr, alq, should all be irrelevant. \n\
VFPPROD \n\
-- table_num, datum_depth, flo, wfr, gfr, pressure, alq, unit, table_vals \n\
42 7.0E+03 LIQ WCT GOR THP ' ' FIELD BHP / \n\
1.0 / flo axis \n\
0.0 1.0 / THP axis \n\
0.0 / WFR axis \n\
0.0 / GFR axis \n\
0.0 / ALQ axis \n\
-- Table itself: thp_idx wfr_idx gfr_idx alq_idx <vals> \n\
1 1 1 1 0.0 / \n\
2 1 1 1 1.0 / \n\
";

    auto units = Opm::UnitSystem::newFIELD();
    Opm::Parser parser;
    auto deck = parser.parseString(table_str);

    BOOST_REQUIRE(deck.hasKeyword("VFPPROD"));
    BOOST_CHECK_EQUAL(deck.count("VFPPROD"), 1);

    Opm::VFPProdTable table(deck.getKeyword("VFPPROD", 0), units);

    Opm::VFPProdProperties properties(&table);

    const int n = 5; //Number of points to check per axis
    double bhp_sad = 0.0; //Sum of absolute difference
    double bhp_max_d = 0.0; //Maximum difference
    double thp_sad = 0.0;
    double thp_max_d = 0.0;
    for (int w=0; w<n; ++w) { //water
        for (int o=0; o<n; ++o) { //oil
            for (int g=0; g<n; ++g) { //gas
                for (int t=0; t<n; ++t) { //thp
                    for (int a=0; a<n; ++a) { //alq
                        double aqua = w * 52.3;
                        double liquid = o * 9.9;
                        double vapour = g * 0.1;
                        double thp = t * 456.78;
                        double alq = a * 42.24;

                        double bhp_interp = properties.bhp(42, aqua, liquid, vapour, thp, alq);
                        double bhp_ref = thp;
                        double thp_interp = properties.thp(42, aqua, liquid, vapour, bhp_ref, alq);
                        double thp_ref = thp;

                        double bhp_diff = std::abs(bhp_interp - bhp_ref);
                        bhp_sad += bhp_diff;
                        bhp_max_d = std::max(bhp_diff, bhp_max_d);

                        double thp_diff = std::abs(thp_interp - thp_ref);
                        thp_sad += thp_diff;
                        thp_max_d = std::max(thp_diff, thp_max_d);
                    }
                }
            }
        }
    }

    BOOST_CHECK_SMALL(bhp_max_d, max_d_tol);
    BOOST_CHECK_SMALL(bhp_sad, sad_tol);

    BOOST_CHECK_SMALL(thp_max_d, max_d_tol);
    BOOST_CHECK_SMALL(thp_sad, sad_tol);
}

/**
 * Tests that we can actually parse some input data, and interpolate within that space
 */
BOOST_AUTO_TEST_CASE(ParseInterpolateRealisticVFPPROD)
{
    auto units = Opm::UnitSystem::newMETRIC();

    Opm::ParseContext parseContext;
    Opm::ErrorGuard errorGuard;
    Opm::Parser parser;
    boost::filesystem::path file("VFPPROD2");

    auto deck = parser.parseFile(file.string());
    Opm::checkDeck(deck, parser, parseContext, errorGuard);

    BOOST_REQUIRE(deck.hasKeyword("VFPPROD"));
    BOOST_CHECK_EQUAL(deck.count("VFPPROD"), 1);

    Opm::VFPProdTable table(deck.getKeyword("VFPPROD", 0), units);

    Opm::VFPProdProperties properties(&table);

    //Do some rudimentary testing
    //Get the BHP as a function of rate, thp, wfr, gfr, alq
    double liq[] = {100, 2942.8571428571427, 5785.7142857142853, 8628.5714285714294, 11471.428571428571, 14314.285714285714, 17157.142857142859, 20000};
    double gor[] = {90, 1505.7142857142858, 2921.4285714285716, 4337.1428571428569, 5752.8571428571431, 7168.5714285714284, 8584.2857142857138, 10000};
    double wct[] = {0, 0.14285714285714285, 0.2857142857142857, 0.42857142857142855, 0.5714285714285714, 0.7142857142857143, 0.8571428571428571, 1};
    double thp[] = {16.010000000000002, 22.438571428571429, 28.867142857142859, 35.295714285714283, 41.724285714285713, 48.152857142857144, 54.581428571428575, 61.009999999999998};
    int n = sizeof(liq) / sizeof(liq[0]);

    int i = 0;
    double sad = 0.0; //Sum of absolute difference
    double max_d = 0.0; //Maximum difference
    for (int t=0; t<n; ++t) {
        for (int w=0; w<n; ++w) {
            for (int g=0; g<n; ++g) {
                //for (unsigned int a=0; a<n; ++a) { //n==1, skip this loop
                    for (int f=0; f<n; ++f) {
                        //Liq given as SM3/day => convert to SM3/second
                        double f_i = -liq[f]*1.1574074074074073e-05;

                        //THP given as BARSA => convert to Pascal
                        double t_i = thp[t]*100000.0;

                        //WCT given as fraction, SM3/SM3
                        double w_i = wct[w];

                        //GOR given as SM3 / SM3
                        double g_i = gor[g];

                        //ALQ unit not relevant in this case
                        double a_i = 0.0;

                        //Now reconstruct possible aqua, liquid,
                        //and vapour phase compositions that satisfy the above

                        //liq[f] = aqua + liquid
                        //wct[w] = aqua / (aqua + liquid)
                        //gor[g] = vapour / liquid
                        //
                        // aqua = wct[w] * liq[f]
                        // liquid = liq[f] - aqua
                        // vapour = gor[g] * liquid

                        double aqua = w_i * f_i;
                        double liquid = f_i - aqua;
                        double vapour = g_i * liquid;

                        if ((aqua + liquid) == 0.0 || liquid == 0.0) {
                            //FIXME: This skips some corner cases, when
                            //getWFR(...) and getGFR(...) are infinite
                        }
                        else {
                            //Value given as pascal, convert to barsa for comparison with reference
                            double value_i = properties.bhp(32, aqua, liquid, vapour, t_i, a_i) * 10.0e-6;

                            double abs_diff = std::abs(value_i - reference[i]);
                            sad += abs_diff;
                            max_d = std::max(max_d, abs_diff);
                        }
                        ++i;
                    }
                //}
            }
        }
    }

    BOOST_CHECK_SMALL(max_d, max_d_tol);
    BOOST_CHECK_SMALL(sad, sad_tol);
}

/**
 * Reference computed using MATLAB with the input above.
 */
const double reference[] = {
        44.850000000000001,  30.36904761904762, 40.238928571428566, 53.221714285714292, 67.885571428571438, 83.100857142857151, 99.021000000000015,             115.14, 27.986149999999999, 76.394680612244898, 146.58456352040815, 227.73038755102044, 323.29702921768705, 433.02474183673473, 519.09348404081641, 596.70497857142868, 25.483874999999998, 135.96136479591837, 276.35715114795914, 452.68607295918366, 676.76830807823114, 945.81213010204078,  1108.470853367347, 1234.3462321428572,            22.9816, 195.52804897959183, 406.12973877551008, 677.64175836734682, 1030.2395869387753, 1458.5995183673467, 1697.8482226938775, 1871.9874857142856,           21.93648, 260.13428359183678, 566.16254540816328, 972.26681151020421, 1344.8295676190476, 1774.9372779591836, 2013.0210626204084, 2185.0791771428571, 22.174319999999998,  329.1776176326531, 752.83812551020401,  1328.232636244898, 1625.1862990476188, 1918.3098995918363, 2086.7690417469385, 2212.4194514285709, 22.412159999999997, 398.22095167346936, 939.51370561224473,  1684.198460979592, 1905.5430304761903, 2061.6825212244894, 2160.5170208734694, 2239.7597257142857, 22.649999999999999, 467.26428571428573, 1126.1892857142857, 2040.1642857142861, 2185.8997619047618, 2205.0551428571425, 2234.2649999999999, 2267.0999999999999, 45.935714285714283, 31.165795918367351, 39.937755102040811,  52.72722448979593, 66.603455782312921, 81.156040816326524, 96.557926530612264, 112.18714285714286, 29.541501020408166,  68.10885813411079, 127.61792838921284, 195.14191513119539, 273.51727814382895, 362.26484685131197, 464.80925879300298, 571.05809489795934, 26.418334183673469, 118.58690218658893, 235.65489568148689, 377.46362485422748, 554.83594890670554, 764.83203053935858, 1015.1273887900876, 1276.0984719387757, 23.295167346938779, 169.06494623906704, 343.69186297376086,  559.7853345772595, 836.15461966958196, 1167.3992142274051, 1565.4455187871722, 1981.1388489795918, 21.823392653061227, 222.83265130417885, 472.45379858600586, 795.44837871720142, 1174.8842335704569, 1543.4152790437317, 1945.2266863183679, 2358.5925306122449, 21.805595102040819, 279.49675619825075, 619.46314463556848,  1078.076089212828, 1564.1616228338191, 1896.0542676618074, 2174.8576629877552, 2447.6207346938772, 21.787797551020411, 336.16086109232265, 766.47249068513111, 1360.7037997084549, 1953.4390120971816, 2248.6932562798829, 2404.4886396571433,   2536.64893877551, 21.770000000000003,  392.8249659863946, 913.48183673469384,  1643.331510204082, 2342.7164013605443, 2601.3322448979588, 2634.1196163265308, 2625.6771428571428,  47.03857142857143, 32.294666666666664, 39.655204081632654, 52.333755102040818, 66.056870748299303, 79.655795918367346, 94.904048979591849, 110.69571428571429,   31.4618112244898, 60.027269698736646, 109.31101308309036, 164.18323034985428, 226.53348555879492, 296.57316871720116, 377.61696073469392,  461.6498928571429, 27.707640306122446, 101.52058415937805, 197.02186251822158, 308.25257120991256, 443.23278785228376, 602.26389059766768,  794.2130148979594, 995.03794642857144, 23.953469387755103,  143.0138986200194, 284.73271195335275, 452.32191206997084, 659.93209014577235, 907.95461247813398, 1210.8090690612244, 1528.4259999999999, 22.096767346938773, 186.56175919922259,  385.1060212390671, 631.63031250145787,  943.5575367035957, 1305.3664132711372, 1664.3889660991258, 2019.3338204081635, 21.910702040816325, 231.91855601943632,  496.6280617784256, 841.96513350437328, 1286.1084689135077, 1783.5345340174924, 2150.5314876034986, 2472.8396897959183, 21.724636734693878, 277.27535283965011, 608.15010231778422, 1052.2999545072887, 1628.6594011234206, 2261.7026547638479, 2636.6740091078714, 2926.3455591836728,  21.53857142857143,  322.6321496598639, 719.67214285714283, 1262.6347755102042,  1971.210333333333, 2739.8707755102037, 3122.8165306122451, 3379.8514285714286, 48.172142857142859,   34.1152925170068, 39.812755102040818,   51.5600612244898, 65.615465986394554, 79.399040816326519, 94.464034693877551, 109.96857142857144, 33.996096428571427,  52.22108403790088, 92.038355357142862,  135.3683000145773, 183.628
};

BOOST_AUTO_TEST_SUITE_END() // Integration tests//