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

Minor fixes for comments in PR

Browse Source
This commit is contained in:
babrodtk 2015-08-19 08:07:51 +02:00
parent b75baac8fc
commit 6b3356e74d
6 changed files with 57 additions and 42 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

28
opm/autodiff/BlackoilModelBase_impl.hpp
View File

@ -1244,7 +1244,6 @@ namespace detail {
const int np = wells().number_of_phases;
const int nw = wells().number_of_wells;
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
const std::vector<double> wr = xw.wellRates();
for (int w = 0; w < nw; ++w) {
const WellControls* wc = wells().ctrls[w];
// The current control in the well state overrides
@ -1300,13 +1299,13 @@ namespace detail {
double vapour = 0.0;
if (active_[ Water ]) {
aqua = wr[w*np + pu.phase_pos[ Water ] ];
aqua = xw.wellRates()[w*np + pu.phase_pos[ Water ] ];
}
if (active_[ Oil ]) {
liquid = wr[w*np + pu.phase_pos[ Oil ] ];
liquid = xw.wellRates()[w*np + pu.phase_pos[ Oil ] ];
}
if (active_[ Gas ]) {
vapour = wr[w*np + pu.phase_pos[ Gas ] ];
vapour = xw.wellRates()[w*np + pu.phase_pos[ Gas ] ];
}
const int vfp = well_controls_iget_vfp(wc, current);
@ -1459,9 +1458,19 @@ namespace detail {
const int np = wells().number_of_phases;
const int nw = wells().number_of_wells;
const ADB& aqua = subset(state.qs, Span(nw, 1, BlackoilPhases::Aqua*nw));
const ADB& liquid = subset(state.qs, Span(nw, 1, BlackoilPhases::Liquid*nw));
const ADB& vapour = subset(state.qs, Span(nw, 1, BlackoilPhases::Vapour*nw));
ADB aqua = ADB::constant(ADB::V::Zero(nw));
ADB liquid = ADB::constant(ADB::V::Zero(nw));
ADB vapour = ADB::constant(ADB::V::Zero(nw));
if (active_[Water]) {
aqua += subset(state.qs, Span(nw, 1, BlackoilPhases::Aqua*nw));
}
if (active_[Oil]) {
liquid += subset(state.qs, Span(nw, 1, BlackoilPhases::Liquid*nw));
}
if (active_[Gas]) {
vapour += subset(state.qs, Span(nw, 1, BlackoilPhases::Vapour*nw));
}
//THP calculation variables
std::vector<int> inj_table_id(nw, -1);
@ -2332,6 +2341,9 @@ namespace detail {
converged_Well = converged_Well && (residualWell < Opm::unit::barsa);
const bool converged = converged_MB && converged_CNV && converged_Well;
// Residual in Pascal can have high values and still be ok.
const double maxWellResidualAllowed = 1000.0 * maxResidualAllowed();
// if one of the residuals is NaN, throw exception, so that the solver can be restarted
if ( std::isnan(mass_balance_residual[Water]) || mass_balance_residual[Water] > maxResidualAllowed() ||
std::isnan(mass_balance_residual[Oil]) || mass_balance_residual[Oil] > maxResidualAllowed() ||
@ -2342,7 +2354,7 @@ namespace detail {
std::isnan(well_flux_residual[Water]) || well_flux_residual[Water] > maxResidualAllowed() ||
std::isnan(well_flux_residual[Oil]) || well_flux_residual[Oil] > maxResidualAllowed() ||
std::isnan(well_flux_residual[Gas]) || well_flux_residual[Gas] > maxResidualAllowed() ||
std::isnan(residualWell) || residualWell > 100*maxResidualAllowed() )
std::isnan(residualWell) || residualWell > maxWellResidualAllowed )
{
OPM_THROW(Opm::NumericalProblem,"One of the residuals is NaN or too large!");
}

7
opm/autodiff/SimulatorBase_impl.hpp
View File

@ -439,9 +439,12 @@ namespace Opm
well_controls_clear(ctrl);
well_controls_assert_number_of_phases(ctrl, int(np));
static const double invalid_alq = -std::numeric_limits<double>::max();
static const int invalid_vfp = -std::numeric_limits<int>::max();
const int ok_resv =
well_controls_add_new(RESERVOIR_RATE, target,
-std::numeric_limits<int>::max(), -std::numeric_limits<int>::max(),
invalid_alq, invalid_vfp,
& distr[0], ctrl);
// For WCONHIST/RESV the BHP limit is set to 1 atm.
@ -449,7 +452,7 @@ namespace Opm
// the WELTARG keyword
const int ok_bhp =
well_controls_add_new(BHP, unit::convert::from(1.0, unit::atm),
-std::numeric_limits<int>::max(), -std::numeric_limits<int>::max(),
invalid_alq, invalid_vfp,
NULL, ctrl);
if (ok_resv != 0 && ok_bhp != 0) {

38
opm/autodiff/VFPHelpers.hpp
View File

@ -256,8 +256,8 @@ inline InterpData findInterpData(const double& value, const std::vector<double>&
/**
* An "ADB-like" structure with a single value and a set of derivatives
*/
struct adb_like {
adb_like() : value(0.0), dthp(0.0), dwfr(0.0), dgfr(0.0), dalq(0.0), dflo(0.0) {};
struct VFPEvaluation {
VFPEvaluation() : value(0.0), dthp(0.0), dwfr(0.0), dgfr(0.0), dalq(0.0), dflo(0.0) {};
double value;
double dthp;
double dwfr;
@ -266,9 +266,9 @@ struct adb_like {
double dflo;
};
inline adb_like operator+(
adb_like lhs,
const adb_like& rhs) {
inline VFPEvaluation operator+(
VFPEvaluation lhs,
const VFPEvaluation& rhs) {
lhs.value += rhs.value;
lhs.dthp += rhs.dthp;
lhs.dwfr += rhs.dwfr;
@ -278,9 +278,9 @@ inline adb_like operator+(
return lhs;
}
inline adb_like operator-(
adb_like lhs,
const adb_like& rhs) {
inline VFPEvaluation operator-(
VFPEvaluation lhs,
const VFPEvaluation& rhs) {
lhs.value -= rhs.value;
lhs.dthp -= rhs.dthp;
lhs.dwfr -= rhs.dwfr;
@ -290,10 +290,10 @@ inline adb_like operator-(
return lhs;
}
inline adb_like operator*(
inline VFPEvaluation operator*(
double lhs,
const adb_like& rhs) {
adb_like retval;
const VFPEvaluation& rhs) {
VFPEvaluation retval;
retval.value = rhs.value * lhs;
retval.dthp = rhs.dthp * lhs;
retval.dwfr = rhs.dwfr * lhs;
@ -317,7 +317,7 @@ inline adb_like operator*(
#endif
inline adb_like interpolate(
inline VFPEvaluation interpolate(
const VFPProdTable::array_type& array,
const InterpData& flo_i,
const InterpData& thp_i,
@ -326,7 +326,7 @@ inline adb_like interpolate(
const InterpData& alq_i) {
//Values and derivatives in a 5D hypercube
adb_like nn[2][2][2][2][2];
VFPEvaluation nn[2][2][2][2][2];
//Pick out nearest neighbors (nn) to our evaluation point
@ -431,13 +431,13 @@ inline adb_like interpolate(
inline adb_like interpolate(
inline VFPEvaluation interpolate(
const VFPInjTable::array_type& array,
const InterpData& flo_i,
const InterpData& thp_i) {
//Values and derivatives in a 5D hypercube
adb_like nn[2][2];
VFPEvaluation nn[2][2];
//Pick out nearest neighbors (nn) to our evaluation point
@ -503,7 +503,7 @@ inline adb_like interpolate(
inline adb_like bhp(const VFPProdTable* table,
inline VFPEvaluation bhp(const VFPProdTable* table,
const double& aqua,
const double& liquid,
const double& vapour,
@ -521,7 +521,7 @@ inline adb_like bhp(const VFPProdTable* table,
auto gfr_i = detail::findInterpData(gfr, table->getGFRAxis());
auto alq_i = detail::findInterpData(alq, table->getALQAxis());
detail::adb_like retval = detail::interpolate(table->getTable(), flo_i, thp_i, wfr_i, gfr_i, alq_i);
detail::VFPEvaluation retval = detail::interpolate(table->getTable(), flo_i, thp_i, wfr_i, gfr_i, alq_i);
return retval;
}
@ -530,7 +530,7 @@ inline adb_like bhp(const VFPProdTable* table,
inline adb_like bhp(const VFPInjTable* table,
inline VFPEvaluation bhp(const VFPInjTable* table,
const double& aqua,
const double& liquid,
const double& vapour,
@ -543,7 +543,7 @@ inline adb_like bhp(const VFPInjTable* table,
auto thp_i = detail::findInterpData(thp, table->getTHPAxis());
//Then perform the interpolation itself
detail::adb_like retval = detail::interpolate(table->getTable(), flo_i, thp_i);
detail::VFPEvaluation retval = detail::interpolate(table->getTable(), flo_i, thp_i);
return retval;
}

4
opm/autodiff/VFPInjProperties.cpp
View File

@ -133,7 +133,7 @@ VFPInjProperties::ADB VFPInjProperties::bhp(const std::vector<int>& table_id,
auto flo_i = detail::findInterpData(flo.value()[i], table->getFloAxis());
auto thp_i = detail::findInterpData(thp.value()[i], table->getTHPAxis());
detail::adb_like bhp_val = detail::interpolate(table->getTable(), flo_i, thp_i);
detail::VFPEvaluation bhp_val = detail::interpolate(table->getTable(), flo_i, thp_i);
value[i] = bhp_val.value;
dthp[i] = bhp_val.dthp;
@ -179,7 +179,7 @@ double VFPInjProperties::bhp(int table_id,
const double& thp) const {
const VFPInjTable* table = detail::getTable(m_tables, table_id);
detail::adb_like retval = detail::bhp(table, aqua, liquid, vapour, thp);
detail::VFPEvaluation retval = detail::bhp(table, aqua, liquid, vapour, thp);
return retval.value;
}

6
opm/autodiff/VFPProdProperties.cpp
View File

@ -128,7 +128,7 @@ VFPProdProperties::ADB VFPProdProperties::bhp(const std::vector<int>& table_id,
auto gfr_i = detail::findInterpData(gfr.value()[i], table->getGFRAxis());
auto alq_i = detail::findInterpData(alq.value()[i], table->getALQAxis());
detail::adb_like bhp_val = detail::interpolate(table->getTable(), flo_i, thp_i, wfr_i, gfr_i, alq_i);
detail::VFPEvaluation bhp_val = detail::interpolate(table->getTable(), flo_i, thp_i, wfr_i, gfr_i, alq_i);
/*
static const int N=40;
@ -145,7 +145,7 @@ VFPProdProperties::ADB VFPProdProperties::bhp(const std::vector<int>& table_id,
const double dist = end - start;
double flo_d = start + (j/static_cast<double>(N-1)) * dist;
auto flo_i = detail::findInterpData(flo_d, table->getFloAxis());
detail::adb_like bhp_val = detail::interpolate(table->getTable(), flo_i, thp_i, wfr_i, gfr_i, alq_i);
detail::VFPEvaluation bhp_val = detail::interpolate(table->getTable(), flo_i, thp_i, wfr_i, gfr_i, alq_i);
std::cout << bhp_val.value << ",";
}
std::cout << "];" << std::endl;
@ -220,7 +220,7 @@ double VFPProdProperties::bhp(int table_id,
const double& alq) const {
const VFPProdTable* table = detail::getTable(m_tables, table_id);
detail::adb_like retval = detail::bhp(table, aqua, liquid, vapour, thp, alq);
detail::VFPEvaluation retval = detail::bhp(table, aqua, liquid, vapour, thp, alq);
return retval.value;
}

16
tests/test_vfpproperties.cpp
View File

@ -229,7 +229,7 @@ BOOST_AUTO_TEST_SUITE_END() // unit tests
*/
struct TrivialFixture {
typedef Opm::VFPProdProperties::ADB ADB;
typedef Opm::detail::adb_like adb_like;
typedef Opm::detail::VFPEvaluation VFPEvaluation;
TrivialFixture() : table_ids(1, 1),
thp_axis{0.0, 1.0},
@ -408,7 +408,7 @@ BOOST_AUTO_TEST_CASE(GetTable)
ADB qs_adb = ADB::constant(qs_adb_v);
//Check that our reference has not changed
Opm::detail::adb_like ref= Opm::detail::bhp(&table, aqua_d, liquid_d, vapour_d, thp_d, alq_d);
Opm::detail::VFPEvaluation ref= Opm::detail::bhp(&table, aqua_d, liquid_d, vapour_d, thp_d, alq_d);
BOOST_CHECK_CLOSE(ref.value, 1.0923565702101556, max_d_tol);
BOOST_CHECK_CLOSE(ref.dthp, 0.13174065498177251, max_d_tol);
BOOST_CHECK_CLOSE(ref.dwfr, -1.2298177745501071, max_d_tol);
@ -810,8 +810,8 @@ BOOST_AUTO_TEST_CASE(PartialDerivatives)
initProperties();
//Temps used to store reference and actual variables
adb_like sad;
adb_like max_d;
VFPEvaluation sad;
VFPEvaluation max_d;
//Check interpolation
for (int i=0; i<=n; ++i) {
@ -831,7 +831,7 @@ BOOST_AUTO_TEST_CASE(PartialDerivatives)
double gfr = Opm::detail::getGFR(aqua, liquid, vapour, table.getGFRType());
//Calculate reference
adb_like reference;
VFPEvaluation reference;
reference.value = thp + 2*wfr + 3*gfr+ 4*alq + 5*flo;
reference.dthp = 1;
reference.dwfr = 2;
@ -841,9 +841,9 @@ BOOST_AUTO_TEST_CASE(PartialDerivatives)
//Calculate actual
//Note order of arguments: id, aqua, liquid, vapour, thp, alq
adb_like actual = Opm::detail::bhp(&table, aqua, liquid, vapour, thp, alq);
VFPEvaluation actual = Opm::detail::bhp(&table, aqua, liquid, vapour, thp, alq);
adb_like abs_diff = actual - reference;
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);
@ -1010,7 +1010,7 @@ BOOST_AUTO_TEST_CASE(ParseInterpolateRealisticVFPPROD)
Opm::ParserPtr parser(new Opm::Parser());
Opm::ParseMode parse_mode;
boost::filesystem::path file("tests/VFPPROD2");
boost::filesystem::path file("VFPPROD2");
deck = parser->parseFile(file.string(), parse_mode);
Opm::checkDeck(deck);