OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-11-25 10:40:21 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #996 from andlaus/no_eigen_in_VFP_and_FIP

Browse Source 
mostly eliminate Eigen in the FIP and VFP code
This commit is contained in:
Atgeirr Flø Rasmussen 2017-01-02 12:56:30 +01:00 committed by GitHub
commit edf883e747
14 changed files with 112 additions and 123 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
opm/autodiff/BlackoilModelBase.hpp
View File

@ -268,7 +268,7 @@ namespace Opm {
/// \param[in] ReservoirState
/// \param[in] FIPNUM for active cells not global cells.
/// \return fluid in place, number of fip regions, each region contains 5 values which are liquid, vapour, water, free gas and dissolved gas.
std::vector<V>
std::vector<std::vector<double> >
computeFluidInPlace(const ReservoirState& x,
const std::vector<int>& fipnum);

13
opm/autodiff/BlackoilModelBase_impl.hpp
View File

@ -2243,7 +2243,7 @@ typedef Eigen::Array<double,
template <class Grid, class WellModel, class Implementation>
std::vector<V>
std::vector<std::vector<double> >
BlackoilModelBase<Grid, WellModel, Implementation>::
computeFluidInPlace(const ReservoirState& x,
const std::vector<int>& fipnum)
@ -2282,7 +2282,10 @@ typedef Eigen::Array<double,
// For a parallel run this is just a local maximum and needs to be updated later
int dims = *std::max_element(fipnum.begin(), fipnum.end());
std::vector<V> values(dims, V::Zero(7));
std::vector<std::vector<double> > values(dims);
for (int i=0; i < dims; ++i) {
values[i].resize(7, 0.0);
}
const V hydrocarbon = saturation[Oil].value() + saturation[Gas].value();
V hcpv;
@ -2356,7 +2359,11 @@ typedef Eigen::Array<double,
auto comm = pinfo.communicator();
// Compute the global dims value and resize values accordingly.
dims = comm.max(dims);
values.resize(dims, V::Zero(7));
values.resize(dims);
for (int i=0; i < dims; ++i) {
values[i].resize(7);
std::fill(values[i].begin(), values[i].end(), 0.0);
}
//Accumulate phases for each region
for (int phase = 0; phase < maxnp; ++phase) {

2
opm/autodiff/BlackoilSequentialModel.hpp
View File

@ -263,7 +263,7 @@ namespace Opm {
/// \param[in] ReservoirState
/// \param[in] FIPNUM for active cells not global cells.
/// \return fluid in place, number of fip regions, each region contains 5 values which are liquid, vapour, water, free gas and dissolved gas.
std::vector<V>
std::vector<std::vector<double> >
computeFluidInPlace(const ReservoirState& x,
const std::vector<int>& fipnum) const
{

17
opm/autodiff/NonlinearSolver.hpp
View File

@ -21,11 +21,9 @@
#ifndef OPM_NONLINEARSOLVER_HEADER_INCLUDED
#define OPM_NONLINEARSOLVER_HEADER_INCLUDED
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/core/simulator/SimulatorReport.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/simulator/SimulatorTimerInterface.hpp>
#include <opm/autodiff/DuneMatrix.hpp>
#include <dune/common/fmatrix.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <memory>
@ -39,11 +37,6 @@ namespace Opm {
class NonlinearSolver
{
public:
// --------- Types and enums ---------
typedef AutoDiffBlock<double> ADB;
typedef ADB::V V;
typedef ADB::M M;
// Available relaxation scheme types.
enum RelaxType { DAMPEN, SOR };
@ -132,14 +125,13 @@ namespace Opm {
/// \param[in] ReservoirState
/// \param[in] FIPNUM for active cells not global cells.
/// \return fluid in place, number of fip regions, each region contains 5 values which are liquid, vapour, water, free gas and dissolved gas.
std::vector<V>
computeFluidInPlace(const ReservoirState& x,
const std::vector<int>& fipnum) const
std::vector<std::vector<double> >
computeFluidInPlace(const ReservoirState& x, const std::vector<int>& fipnum) const
{
return model_->computeFluidInPlace(x, fipnum);
}
std::vector<std::vector<double>>
std::vector<std::vector<double> >
computeFluidInPlace(const std::vector<int>& fipnum) const
{
return model_->computeFluidInPlace(fipnum);
@ -156,9 +148,6 @@ namespace Opm {
void detectOscillations(const std::vector<std::vector<double>>& residual_history,
const int it, bool& oscillate, bool& stagnate) const;
/// Apply a stabilization to dx, depending on dxOld and relaxation parameters.
void stabilizeNonlinearUpdate(V& dx, V& dxOld, const double omega) const;
/// Apply a stabilization to dx, depending on dxOld and relaxation parameters.
/// Implemention for Dune block vectors.
template <class BVector>

44
opm/autodiff/NonlinearSolver_impl.hpp
View File

@ -248,36 +248,6 @@ namespace Opm
}
template <class PhysicalModel>
void
NonlinearSolver<PhysicalModel>::stabilizeNonlinearUpdate(V& dx, V& dxOld, const double omega) const
{
// The dxOld is updated with dx.
// If omega is equal to 1., no relaxtion will be appiled.
const V tempDxOld = dxOld;
dxOld = dx;
switch (relaxType()) {
case DAMPEN:
if (omega == 1.) {
return;
}
dx = dx*omega;
return;
case SOR:
if (omega == 1.) {
return;
}
dx = dx*omega + (1.-omega)*tempDxOld;
return;
default:
OPM_THROW(std::runtime_error, "Can only handle DAMPEN and SOR relaxation type.");
}
return;
}
template <class PhysicalModel>
template <class BVector>
void
@ -294,15 +264,19 @@ namespace Opm
if (omega == 1.) {
return;
}
dx *= omega;
for (size_t i = 0; i < dx.size(); ++i) {
dx[i] *= omega;
}
return;
case SOR:
if (omega == 1.) {
return;
}
dx *= omega;
tempDxOld *= (1.-omega);
dx += tempDxOld;
for (size_t i = 0; i < dx.size(); ++i) {
dx[i] *= omega;
tempDxOld[i] *= (1.-omega);
dx[i] += tempDxOld[i];
}
return;
default:
OPM_THROW(std::runtime_error, "Can only handle DAMPEN and SOR relaxation type.");
@ -310,8 +284,6 @@ namespace Opm
return;
}
} // namespace Opm

16
opm/autodiff/SimulatorBase.hpp
View File

@ -21,6 +21,9 @@
#ifndef OPM_SIMULATORBASE_HEADER_INCLUDED
#define OPM_SIMULATORBASE_HEADER_INCLUDED
#include <opm/material/densead/Math.hpp>
#include <opm/autodiff/DuneMatrix.hpp>
#include <opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/common/ErrorMacros.hpp>
@ -159,17 +162,18 @@ namespace Opm
WellState& xw);
void
FIPUnitConvert(const UnitSystem& units, V& fip);
FIPUnitConvert(const UnitSystem& units,
std::vector<std::vector<double> >& fip);
void
FIPUnitConvert(const UnitSystem& units,
std::vector<V>& fip);
V
FIPTotals(const std::vector<V>& fip, const ReservoirState& state);
std::vector<double>& fip);
std::vector<double>
FIPTotals(const std::vector<std::vector<double> >& fip, const ReservoirState& state);
void
outputFluidInPlace(const V& oip, const V& cip, const UnitSystem& units, const int reg);
outputFluidInPlace(const std::vector<double>& oip, const std::vector<double>& cip, const UnitSystem& units, const int reg);
void computeWellPotentials(const Wells* wells,
const WellState& xw,

61
opm/autodiff/SimulatorBase_impl.hpp
View File

@ -149,7 +149,7 @@ namespace Opm
fipnum[c] = fipnum_global[AutoDiffGrid::globalCell(grid_)[c]];
}
}
std::vector<V> OOIP;
std::vector<std::vector<double> > OOIP;
// Main simulation loop.
while (!timer.done()) {
// Report timestep.
@ -275,10 +275,10 @@ namespace Opm
report.solver_time += solver_timer.secsSinceStart();
// Compute current FIP.
std::vector<V> COIP;
std::vector<std::vector<double> > COIP;
COIP = solver->computeFluidInPlace(state, fipnum);
V OOIP_totals = FIPTotals(OOIP, state);
V COIP_totals = FIPTotals(COIP, state);
std::vector<double> OOIP_totals = FIPTotals(OOIP, state);
std::vector<double> COIP_totals = FIPTotals(COIP, state);
//Convert to correct units
FIPUnitConvert(eclipse_state_->getUnits(), COIP);
@ -665,21 +665,19 @@ namespace Opm
}
}
template <class Implementation>
void
SimulatorBase<Implementation>::FIPUnitConvert(const UnitSystem& units,
std::vector<V>& fip)
std::vector<std::vector<double> >& fip)
{
for (size_t i = 0; i < fip.size(); ++i) {
FIPUnitConvert(units, fip[i]);
}
}
template <class Implementation>
void
SimulatorBase<Implementation>::FIPUnitConvert(const UnitSystem& units, V& fip)
SimulatorBase<Implementation>::FIPUnitConvert(const UnitSystem& units, std::vector<double>& fip)
{
if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_FIELD) {
fip[0] = unit::convert::to(fip[0], unit::stb);
@ -700,10 +698,10 @@ namespace Opm
template <class Implementation>
V
SimulatorBase<Implementation>::FIPTotals(const std::vector<V>& fip, const ReservoirState& state)
std::vector<double>
SimulatorBase<Implementation>::FIPTotals(const std::vector<std::vector<double> >& fip, const ReservoirState& state)
{
V totals(V::Zero(7));
std::vector<double> totals(7, 0.0);
for (int i = 0; i < 5; ++i) {
for (size_t reg = 0; reg < fip.size(); ++reg) {
totals[i] += fip[reg][i];
@ -713,15 +711,34 @@ namespace Opm
const int np = state.numPhases();
const PhaseUsage& pu = props_.phaseUsage();
const DataBlock s = Eigen::Map<const DataBlock>(& state.saturation()[0], nc, np);
const V so = pu.phase_used[BlackoilPhases::Liquid] ? V(s.col(BlackoilPhases::Liquid)) : V::Zero(nc);
const V sg = pu.phase_used[BlackoilPhases::Vapour] ? V(s.col(BlackoilPhases::Vapour)) : V::Zero(nc);
const V hydrocarbon = so + sg;
const V p = Eigen::Map<const V>(& state.pressure()[0], nc);
std::vector<double> so(nc);
std::vector<double> sg(nc);
std::vector<double> hydrocarbon(nc);
const auto& soCol = s.col(BlackoilPhases::Liquid);
const auto& sgCol = s.col(BlackoilPhases::Vapour);
for (unsigned c = 0; c < so.size(); ++ c) {
double mySo = 0.0;
if (pu.phase_used[BlackoilPhases::Liquid])
mySo = soCol[c];
double mySg = 0.0;
if (pu.phase_used[BlackoilPhases::Vapour])
mySg = sgCol[c];
so[c] = mySo;
sg[c] = mySg;
hydrocarbon[c] = mySo + mySg;
}
const std::vector<double> p = state.pressure();
if ( ! is_parallel_run_ )
{
double tmp = 0.0;
double tmp2 = 0.0;
for (unsigned i = 0; i < p.size(); ++i) {
tmp += p[i] * geo_.poreVolume()[i] * hydrocarbon[i];
tmp2 += geo_.poreVolume()[i] * hydrocarbon[i];
}
totals[5] = geo_.poreVolume().sum();
totals[6] = (p * geo_.poreVolume() * hydrocarbon).sum() / ((geo_.poreVolume() * hydrocarbon).sum());
}
totals[6] = unit::convert::to(tmp/tmp2, unit::barsa);
}
else
{
#if HAVE_MPI
@ -730,8 +747,12 @@ namespace Opm
auto operators = std::make_tuple(Opm::Reduction::makeGlobalSumFunctor<double>(),
Opm::Reduction::makeGlobalSumFunctor<double>(),
Opm::Reduction::makeGlobalSumFunctor<double>());
auto pav_nom = p * geo_.poreVolume() * hydrocarbon;
auto pav_denom = geo_.poreVolume() * hydrocarbon;
std::vector<double> pav_nom(p.size());
std::vector<double> pav_denom(pav_nom.size());
for (unsigned i = 0; i < p.size(); ++i) {
pav_nom[i] = p[i] * geo_.poreVolume()[i] * hydrocarbon[i];
pav_denom[i] = geo_.poreVolume()[i] * hydrocarbon[i];
}
// using ref cref to prevent copying
auto inputs = std::make_tuple(std::cref(geo_.poreVolume()),
@ -756,7 +777,7 @@ namespace Opm
template <class Implementation>
void
SimulatorBase<Implementation>::outputFluidInPlace(const V& oip, const V& cip, const UnitSystem& units, const int reg)
SimulatorBase<Implementation>::outputFluidInPlace(const std::vector<double>& oip, const std::vector<double>& cip, const UnitSystem& units, const int reg)
{
std::ostringstream ss;
if (!reg) {

22
opm/autodiff/VFPInjProperties.cpp
View File

@ -27,6 +27,7 @@
#include <opm/parser/eclipse/EclipseState/Tables/VFPProdTable.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/autodiff/VFPHelpers.hpp>
@ -67,14 +68,10 @@ VFPInjProperties::VFPInjProperties(const std::map<int, VFPInjTable>& tables) {
}
}
VFPInjProperties::ADB VFPInjProperties::bhp(const std::vector<int>& table_id,
const Wells& wells,
const ADB& qs,
const ADB& thp_val) const {
const Wells& wells,
const ADB& qs,
const ADB& thp_val) const {
const int nw = wells.number_of_wells;
//Short-hands for water / oil / gas phases
@ -87,14 +84,11 @@ VFPInjProperties::ADB VFPInjProperties::bhp(const std::vector<int>& table_id,
return bhp(table_id, w, o, g, thp_val);
}
VFPInjProperties::ADB VFPInjProperties::bhp(const std::vector<int>& table_id,
const ADB& aqua,
const ADB& liquid,
const ADB& vapour,
const ADB& thp_arg) const {
const ADB& aqua,
const ADB& liquid,
const ADB& vapour,
const ADB& thp_arg) const {
const int nw = thp_arg.size();
std::vector<int> block_pattern = detail::commonBlockPattern(aqua, liquid, vapour, thp_arg);

19
opm/autodiff/VFPInjProperties.hpp
View File

@ -22,7 +22,6 @@
#include <opm/parser/eclipse/EclipseState/Tables/VFPInjTable.hpp>
#include <opm/core/wells.h>
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/material/densead/Math.hpp>
#include <opm/material/densead/Evaluation.hpp>
#include <opm/autodiff/VFPHelpers.hpp>
@ -34,6 +33,8 @@
namespace Opm {
template <class Scalar>
class AutoDiffBlock;
class VFPInjProperties {
public:
@ -154,10 +155,10 @@ public:
* the above parameters from the values in the input table.
*/
double bhp(int table_id,
const double& aqua,
const double& liquid,
const double& vapour,
const double& thp) const;
const double& aqua,
const double& liquid,
const double& vapour,
const double& thp) const;
/**
@ -172,10 +173,10 @@ public:
* the above parameters from the values in the input table.
*/
double thp(int table_id,
const double& aqua,
const double& liquid,
const double& vapour,
const double& bhp) const;
const double& aqua,
const double& liquid,
const double& vapour,
const double& bhp) const;
/**
* Returns the table associated with the ID, or throws an exception if

22
opm/autodiff/VFPProdProperties.cpp
View File

@ -23,6 +23,7 @@
#include <opm/parser/eclipse/EclipseState/Tables/VFPProdTable.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/material/densead/Math.hpp>
#include <opm/material/densead/Evaluation.hpp>
@ -54,14 +55,11 @@ VFPProdProperties::VFPProdProperties(const std::map<int, VFPProdTable>& tables)
}
}
VFPProdProperties::ADB VFPProdProperties::bhp(const std::vector<int>& table_id,
const Wells& wells,
const ADB& qs,
const ADB& thp_arg,
const ADB& alq) const {
const Wells& wells,
const ADB& qs,
const ADB& thp_arg,
const ADB& alq) const {
const int nw = wells.number_of_wells;
//Short-hands for water / oil / gas phases
@ -76,11 +74,11 @@ VFPProdProperties::ADB VFPProdProperties::bhp(const std::vector<int>& table_id,
VFPProdProperties::ADB VFPProdProperties::bhp(const std::vector<int>& table_id,
const ADB& aqua,
const ADB& liquid,
const ADB& vapour,
const ADB& thp_arg,
const ADB& alq) const {
const ADB& aqua,
const ADB& liquid,
const ADB& vapour,
const ADB& thp_arg,
const ADB& alq) const {
const int nw = thp_arg.size();
std::vector<int> block_pattern = detail::commonBlockPattern(aqua, liquid, vapour, thp_arg, alq);

3
opm/autodiff/VFPProdProperties.hpp
View File

@ -22,7 +22,6 @@
#include <opm/parser/eclipse/EclipseState/Tables/VFPProdTable.hpp>
#include <opm/core/wells.h>
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/material/densead/Math.hpp>
#include <opm/material/densead/Evaluation.hpp>
#include <opm/autodiff/VFPHelpers.hpp>
@ -33,6 +32,8 @@
namespace Opm {
template <class Scalar>
class AutoDiffBlock;
/**
* Class which linearly interpolates BHP as a function of rate, tubing head pressure,

4
opm/autodiff/WellHelpers.hpp
View File

@ -23,7 +23,6 @@
#define OPM_WELLHELPERS_HEADER_INCLUDED
#include <opm/core/wells.h>
#include <opm/autodiff/AutoDiffBlock.hpp>
// #include <opm/autodiff/AutoDiffHelpers.hpp>
#include <vector>
@ -121,7 +120,6 @@ namespace Opm {
// --------- Types ---------
using Vector = AutoDiffBlock<double>::V;
/**
* Simple hydrostatic correction for VFP table
@ -149,7 +147,7 @@ namespace Opm {
return dp;
}
template <class Vector>
inline
Vector computeHydrostaticCorrection(const Wells& wells, const Vector vfp_ref_depth,
const Vector& well_perforation_densities, const double gravity) {

8
opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp
View File

@ -536,7 +536,7 @@ namespace {
std::vector<V>
std::vector<std::vector<double> >
FullyImplicitCompressiblePolymerSolver::computeFluidInPlace(const PolymerBlackoilState& x,
const std::vector<int>& fipnum)
{
@ -568,7 +568,11 @@ namespace {
const int dims = *std::max_element(fipnum.begin(), fipnum.end());
std::vector<V> values(dims, V::Zero(7));
std::vector<std::vector<double> > values(dims);
for (int i=0; i < dims; ++i) {
values[i].resize(7, 0.0);
}
V hcpv = V::Zero(nc);
V pres = V::Zero(nc);
for (int i = 0; i < 5; ++i) {

2
opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp
View File

@ -158,7 +158,7 @@ namespace Opm {
/// \param[in] WellState
/// \param[in] FIPNUM for active cells not global cells.
/// \return fluid in place, number of fip regions, each region contains 5 values which are liquid, vapour, water, free gas and dissolved gas.
std::vector<V>
std::vector<std::vector<double> >
computeFluidInPlace(const PolymerBlackoilState& x,
const std::vector<int>& fipnum);