OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-15 02:43:25 -06:00
Code Issues Packages Projects Releases Wiki Activity

Add license and documentation.

Browse Source
This commit is contained in:
Liu Ming 2014-02-25 09:52:10 +08:00
parent e63e318d40
commit 979d503b7d
17 changed files with 668 additions and 445 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -58,9 +58,11 @@ typedef Eigen::Array<double,
Eigen::RowMajor> DataBlock;
namespace {
namespace {
std::vector<int>
buildAllCells(const int nc)
{
@ -73,6 +75,8 @@ namespace {
template <class GeoProps>
AutoDiffBlock<double>::M
gravityOperator(const UnstructuredGrid& grid,
@ -123,6 +127,8 @@ namespace {
V computePerfPress(const UnstructuredGrid& grid, const Wells& wells, const V& rho, const double grav)
{
const int nw = wells.number_of_wells;
@ -148,23 +154,20 @@ namespace {
return wdp;
}
} // Anonymous namespace
FullyImplicitCompressiblePolymerSolver::
FullyImplicitCompressiblePolymerSolver(const UnstructuredGrid& grid ,
const BlackoilPropsAdInterface& fluid,
const DerivedGeology& geo ,
const RockCompressibility* rock_comp_props,
const PolymerPropsAd& polymer_props_ad,
const Wells& wells,
const LinearSolverInterface& linsolver
)
FullyImplicitCompressiblePolymerSolver(const UnstructuredGrid& grid,
const BlackoilPropsAdInterface& fluid,
const DerivedGeology& geo ,
const RockCompressibility* rock_comp_props,
const PolymerPropsAd& polymer_props_ad,
const Wells& wells,
const LinearSolverInterface& linsolver)
: grid_ (grid)
, fluid_ (fluid)
, geo_ (geo)
@ -464,6 +467,7 @@ namespace {
void
FullyImplicitCompressiblePolymerSolver::
computeCmax(PolymerBlackoilState& state,
@ -478,6 +482,10 @@ namespace {
}
void
FullyImplicitCompressiblePolymerSolver::
assemble(const double dt,
@ -515,7 +523,6 @@ namespace {
residual_.mass_balance[2] = pvdt*(rq_[2].accum[1] - rq_[2].accum[0]) //+ cell / dt * (rq_[2].ads[1] - rq_[2].ads[0])
+ ops_.div*rq_[2].mflux;
// -------- Extra (optional) sg or rs equation, and rs contributions to the mass balance equations --------
// Add the extra (flux) terms to the gas mass balance equations
@ -615,9 +622,7 @@ namespace {
// well rates contribs to polymer mass balance eqn.
// for injection wells.
const V polyin = Eigen::Map<const V>(& polymer_inflow[0], nc);
// std::cout<< "Polymer in flow:" << polyin << std::endl;
const V poly_in_perf = subset(polyin, well_cells);
//const V poly_c_cell = subset(state.concentration, well_cells).value();
const V poly_mc_cell = subset(mc, well_cells).value();
const V poly_in_c = poly_in_perf;// * poly_mc_cell;
const V poly_mc = producer.select(poly_mc_cell, poly_in_c);
@ -654,11 +659,6 @@ namespace {
// Choose bhp residual for positive bhp targets.
Selector<double> bhp_selector(bhp_targets);
residual_.well_eq = bhp_selector.select(bhp_residual, rate_residual);
// for (int i = 0; i < nc; ++i) {
// std::cout << src[i] << " ";
// if ((i+1) % 10 == 0)
// std::cout<<std::endl;
// }
// DUMP(residual_.well_eq);
}
@ -703,9 +703,10 @@ namespace {
void FullyImplicitCompressiblePolymerSolver::updateState(const V& dx,
PolymerBlackoilState& state,
WellState& well_state) const
void FullyImplicitCompressiblePolymerSolver::
updateState(const V& dx,
PolymerBlackoilState& state,
WellState& well_state) const
{
const int np = fluid_.numPhases();
const int nc = grid_.number_of_cells;
@ -740,7 +741,6 @@ namespace {
const V sw_old = s_old.col(0);
const V dsw_limited = sign(dsw) * dsw.abs().min(dsmax);
const V sw = (sw_old - dsw_limited).unaryExpr(Chop01());
// const V sw = (sw_old - dsw);
so -= sw;
for (int c = 0; c < nc; ++c) {
state.saturation()[c*np] = sw[c];
@ -814,6 +814,8 @@ namespace {
std::vector<ADB>
FullyImplicitCompressiblePolymerSolver::
computePressures(const SolutionState& state) const
@ -843,6 +845,8 @@ namespace {
void
FullyImplicitCompressiblePolymerSolver::computeMassFlux(
const V& transi,
@ -957,6 +961,9 @@ namespace {
}
// here mc means m(c) * c.
ADB
FullyImplicitCompressiblePolymerSolver::computeMc(const SolutionState& state) const
@ -967,6 +974,8 @@ namespace {
ADB
FullyImplicitCompressiblePolymerSolver::poroMult(const ADB& p) const
{
@ -1018,4 +1027,4 @@ namespace {
}
} // namespace Opm
} //namespace Opm

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

@ -1,5 +1,6 @@
/*
Copyright 2013 SINTEF ICT, Applied Mathematics.
Copyright 2014 SINTEF ICT, Applied Mathematics.
Copyright 2014 STATOIL.
This file is part of the Open Porous Media project (OPM).
@ -37,11 +38,10 @@ namespace Opm {
class PolymerBlackoilState;
class WellState;
/// A fully implicit solver for the black-oil problem.
/// A fully implicit solver for the oil-water with polymer problem.
///
/// The simulator is capable of handling three-phase problems
/// where gas can be dissolved in oil (but not vice versa). It
/// uses an industry-standard TPFA discretization with per-phase
/// The simulator is capable of handling oil-water-polymer problems
/// It uses an industry-standard TPFA discretization with per-phase
/// upwind weighting of mobilities.
///
/// It uses automatic differentiation via the class AutoDiffBlock
@ -56,34 +56,36 @@ namespace Opm {
/// \param[in] fluid fluid properties
/// \param[in] geo rock properties
/// \param[in] rock_comp_props if non-null, rock compressibility properties
/// \param[in] polymer_props_ad polymer properties
/// \param[in] wells well structure
/// \param[in] linsolver linear solver
FullyImplicitCompressiblePolymerSolver(const UnstructuredGrid& grid ,
const BlackoilPropsAdInterface& fluid,
const DerivedGeology& geo ,
const RockCompressibility* rock_comp_props,
const PolymerPropsAd& polymer_props_ad,
const Wells& wells,
const LinearSolverInterface& linsolver);
const BlackoilPropsAdInterface& fluid,
const DerivedGeology& geo ,
const RockCompressibility* rock_comp_props,
const PolymerPropsAd& polymer_props_ad,
const Wells& wells,
const LinearSolverInterface& linsolver);
/// Take a single forward step, modifiying
/// state.pressure()
/// state.faceflux()
/// state.saturation()
/// state.gasoilratio()
/// state.concentration()
/// wstate.bhp()
/// \param[in] dt time step size
/// \param[in] state reservoir state
/// \param[in] wstate well state
/// \param[in] polymer_inflow polymer influx
/// \param[in] src to caculate wc
void
step(const double dt ,
PolymerBlackoilState& state ,
WellState& wstate,
step(const double dt,
PolymerBlackoilState& state ,
WellState& wstate,
const std::vector<double>& polymer_inflow,
std::vector<double>& src);
std::vector<double>& src);
private:
// Types and enums
typedef AutoDiffBlock<double> ADB;
typedef ADB::V V;
typedef ADB::M M;
@ -99,7 +101,7 @@ namespace Opm {
ADB b; // Reciprocal FVF
ADB head; // Pressure drop across int. interfaces
ADB mob; // Phase mobility (per cell)
std::vector<ADB> ads; //
std::vector<ADB> ads; // Adsorption term.
};
struct SolutionState {
@ -158,11 +160,11 @@ namespace Opm {
const int aix );
void
assemble(const double dt,
assemble(const double dt,
const PolymerBlackoilState& x,
const WellState& xw,
const std::vector<double>& polymer_inflow,
std::vector<double>& src);
const WellState& xw,
const std::vector<double>& polymer_inflow,
std::vector<double>& src);
V solveJacobianSystem() const;
@ -185,6 +187,7 @@ namespace Opm {
const V& transi,
const std::vector<ADB>& kr ,
const SolutionState& state );
void
computeMassFlux(const V& trans,
const ADB& mc,
@ -196,15 +199,20 @@ namespace Opm {
computeFracFlow(const ADB& kro,
const ADB& krw_eff,
const ADB& c) const;
void
computeCmax(PolymerBlackoilState& state,
const ADB& c);
ADB
computeMc(const SolutionState& state) const;
ADB
rockPorosity(const ADB& p) const;
rockPorosity(const ADB& p) const;
ADB
rockPermeability(const ADB& p) const;
rockPermeability(const ADB& p) const;
double
residualNorm() const;

86
opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp
View File

@ -1,3 +1,23 @@
/*
Copyright 2014 SINTEF ICT, Applied Mathematics.
Copyright 2014 STATOIL.
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/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp>
#include <opm/core/pressure/tpfa/trans_tpfa.h>
@ -15,17 +35,20 @@
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/well_controls.h>
#include <cassert>
#include <cmath>
#include <iostream>
#include <iomanip>
#include <Eigen/Eigen>
#include <algorithm>
namespace Opm {
typedef AutoDiffBlock<double> ADB;
typedef ADB::V V;
typedef ADB::M M;
@ -37,6 +60,7 @@ typedef Eigen::Array<double,
namespace {
std::vector<int>
@ -80,8 +104,7 @@ namespace {
return wdp;
}
}//anonymous namespace
} //anonymous namespace
@ -201,6 +224,8 @@ namespace {
FullyImplicitTwophasePolymerSolver::ReservoirResidualQuant::ReservoirResidualQuant()
: accum(2, ADB::null())
, mflux( ADB::null())
@ -212,6 +237,8 @@ namespace {
FullyImplicitTwophasePolymerSolver::SolutionState::SolutionState(const int np)
: pressure ( ADB::null())
, saturation (np, ADB::null())
@ -243,7 +270,6 @@ namespace {
bpat.push_back(xw.bhp().size() * np);
bpat.push_back(xw.bhp().size());
SolutionState state(np);
// Pressure.
@ -357,6 +383,9 @@ namespace {
}
void
FullyImplicitTwophasePolymerSolver::
computeCmax(PolymerState& state,
@ -370,6 +399,10 @@ namespace {
}
void
FullyImplicitTwophasePolymerSolver::
computeAccum(const SolutionState& state,
@ -381,16 +414,20 @@ namespace {
rq_[0].accum[aix] = sat[0];
rq_[1].accum[aix] = sat[1];
const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern());
const ADB ads = polymer_props_ad_.adsorption(state.concentration, cmax);
const double rho_rock = polymer_props_ad_.rockDensity();
const V phi = Eigen::Map<const V>(&fluid_.porosity()[0], grid_.number_of_cells, 1);
const double dead_pore_vol = polymer_props_ad_.deadPoreVol();
rq_[2].accum[aix] = sat[0] * c * (1. - dead_pore_vol) + rho_rock * (1. - phi) / phi * ads;
}
void
FullyImplicitTwophasePolymerSolver::
assemble(const V& pvdt,
@ -547,6 +584,10 @@ namespace {
// residual_.well_eq = bhp_residual;
}
std::vector<ADB>
FullyImplicitTwophasePolymerSolver::
@ -567,6 +608,10 @@ namespace {
return pressure;
}
void
FullyImplicitTwophasePolymerSolver::computeMassFlux(const V& trans,
const ADB& mc,
@ -604,6 +649,9 @@ namespace {
}
std::vector<ADB>
FullyImplicitTwophasePolymerSolver::accumSource(const ADB& kro,
const ADB& krw_eff,
@ -650,6 +698,7 @@ namespace {
std::vector<ADB>
FullyImplicitTwophasePolymerSolver::computeFracFlow() const
{
@ -744,8 +793,6 @@ namespace {
const V c = (c_old - dc).max(zero);
std::copy(&c[0], &c[0] + nc, state.concentration().begin());
// Qs update.
// Since we need to update the wellrates, that are ordered by wells,
// from dqs which are ordered by phase, the simplest is to compute
@ -756,12 +803,10 @@ namespace {
const V wr = wr_old - dwr;
std::copy(&wr[0], &wr[0] + wr.size(), well_state.wellRates().begin());
// Bhp update.
const V bhp_old = Eigen::Map<const V>(&well_state.bhp()[0], nw, 1);
const V bhp = bhp_old - dbhp;
std::copy(&bhp[0], &bhp[0] + bhp.size(), well_state.bhp().begin());
}
@ -782,13 +827,6 @@ namespace {
double
FullyImplicitTwophasePolymerSolver::residualNorm() const
{
@ -803,15 +841,17 @@ namespace {
r = std::max(r, residual_.well_flux_eq.value().matrix().lpNorm<Eigen::Infinity>());
r = std::max(r, residual_.well_eq.value().matrix().lpNorm<Eigen::Infinity>());
return r;
}
ADB
FullyImplicitTwophasePolymerSolver::fluidDensity(const int phase,
const ADB p) const
const ADB p) const
{
const double* rhos = fluid_.surfaceDensity();
ADB rho = ADB::constant(V::Constant(grid_.number_of_cells, 1, rhos[phase]),
@ -819,7 +859,10 @@ namespace {
return rho;
}
V
FullyImplicitTwophasePolymerSolver::transmissibility() const
@ -834,6 +877,10 @@ namespace {
return trans;
}
// here mc means m(c) * c.
ADB
FullyImplicitTwophasePolymerSolver::computeMc(const SolutionState& state) const
@ -843,4 +890,7 @@ namespace {
}
}//namespace Opm
} //namespace Opm

70
opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp
View File

@ -1,3 +1,23 @@
/*
Copyright 2014 SINTEF ICT, Applied Mathematics.
Copyright 2014 STATOIL.
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_FULLYIMPLICITTWOPHASEPOLYMERSOLVER_HEADER_INCLUDED
#define OPM_FULLYIMPLICITTWOPHASEPOLYMERSOLVER_HEADER_INCLUDED
@ -8,24 +28,53 @@
#include <opm/polymer/fullyimplicit/PolymerPropsAd.hpp>
#include <opm/core/pressure/tpfa/trans_tpfa.h>
struct UnstructuredGrid;
struct Wells;
namespace Opm {
class LinearSolverInterface;
class PolymerState;
class WellState;
/// A fully implicit solver for the incompressible oil-water flow wtih polymer problem.
///
/// The simulator is capable of handling oil-water with polymer problems
/// It uses an industry-standard TPFA discretization with per-phase
/// upwind weighting of mobilities.
///
/// It uses automatic differentiation via the class AutoDiffBlock
/// to simplify assembly of the jacobian matrix.
class FullyImplicitTwophasePolymerSolver
{
public:
/// Construct a solver. It will retain references to the
/// arguments of this functions, and they are expected to
/// remain in scope for the lifetime of the solver.
/// \param[in] grid grid data structure
/// \param[in] fluid fluid properties
/// \param[in] polymer_props_ad polymer properties
/// \param[in] wells well structure
/// \param[in] linsolver linear solver
/// \param[in] gravity gravity
FullyImplicitTwophasePolymerSolver(const UnstructuredGrid& grid,
const IncompPropsAdInterface& fluid,
const PolymerPropsAd& polymer_props_ad,
const LinearSolverInterface& linsolver,
const Wells& wells,
const double* gravity);
const LinearSolverInterface& linsolver,
const Wells& wells,
const double* gravity);
/// Take a single forward step, modifiying
/// state.pressure()
/// state.faceflux()
/// state.saturation()
/// state.concentration()
/// wstate.bhp()
/// \param[in] dt time step size
/// \param[in] state reservoir state with polymer
/// \param[in] wstate well state
/// \param[in] polymer_inflow polymer influx
/// \param[in] src to caculate wc
void step(const double dt,
PolymerState& state,
WellState& well_state,
@ -75,6 +124,7 @@ namespace Opm {
const WellOps wops_;
V cmax_;
std::vector<ReservoirResidualQuant> rq_;
struct {
std::vector<ADB> mass_balance;
ADB well_eq;
@ -93,17 +143,21 @@ namespace Opm {
const WellState& xw,
const std::vector<double>& polymer_inflow,
std::vector<double>& src);
V solveJacobianSystem() const;
void updateState(const V& dx,
PolymerState& x,
WellState& xw) const;
std::vector<ADB>
computeRelPerm(const SolutionState& state) const;
V
transmissibility() const;
std::vector<ADB>
computePressures(const SolutionState& state) const;
void
computeMassFlux(const V& trans,
const ADB& mc,
@ -121,8 +175,10 @@ namespace Opm {
std::vector<ADB>
computeFracFlow() const;
double
residualNorm() const;
ADB
polymerSource(const std::vector<ADB>& kr,
const std::vector<double>& src,
@ -137,19 +193,25 @@ namespace Opm {
const int aix );
ADB
computeMc(const SolutionState& state) const;
ADB
rockPorosity(const ADB& p) const;
ADB
rockPermeability(const ADB& p) const;
const double
fluidDensity(const int phase) const;
ADB
fluidDensity(const int phase,
const ADB p) const;
ADB
transMult(const ADB& p) const;
};
} // namespace Opm
#endif// OPM_FULLYIMPLICITTWOPHASESOLVER_HEADER_INCLUDED

5
opm/polymer/fullyimplicit/IncompPropsAd.cpp
View File

@ -1,5 +0,0 @@
/*
Author : Liu Ming
Date : 2013-11-28 in Oslo
Email : miliu@statoil.com
*/

59
opm/polymer/fullyimplicit/IncompPropsAd.hpp
View File

@ -1,59 +0,0 @@
/*
Author : Liu Ming
Data : 2013-11-28 in Oslo
Email : miliu@statoil.com
Properties for incompressible immiscible two-phase flow
*/
#ifndef OPM_INCOMPROPSAD_HEADER_INCLUDED
#define OPM_INCOMPROPSAD_HEADER_INCLUDED
#include <opm/core/props/IncompPropertiesBasic.hpp>
#include <opm/core/props/rock/RockBasic.hpp>
#include <opm/polymer/fullyimplicit/AutoDiffBlock.hpp>
#include <opm/core/props/pvt/PvtPropertiesBasic.hpp>
#include <opm/core/props/satfunc/SaturationPropsBasic.hpp>
namespace Opm
{
// class BlackoilPhases;
class IncompropsAd : public IncompPropertiesBasic
{
public:
IncomPropsAd(const parameter::ParameterGroup& param,
const int dim,
const int num_cells);
IncompPropsAd(const int num_phases,
const SaturationPropsBasic::RelPermFunc& relpermfunc,
const std::vector<double>& rho,
const std::vector<double>& mu,
const double porosity,
const double permeability,
const int dim,
const int num_cells);
~IncompPropsAd();
typedef AutoDiffBlock<double> ADB;
typedef ADB::V V;
V relperm(const V& sw,
const V& so,
const std::vector<int>& cells);
ADB relperm(const ADB& sw,
const ADB& so,
const std::vector<int>& cells);
V capPress(const V& sw,
const V& so,
const std::vector<int>& cells);
ADB capPress(const ADB& sw,
const ADB& so,
const std::vector<int>& cells);
}
}
#endif// OPM_INCOMPROPSAD_HEADER_INCLUDED

57
opm/polymer/fullyimplicit/IncompPropsAdBasic.cpp
View File

@ -1,4 +1,22 @@
/**/
/*
Copyright 2014 SINTEF ICT, Applied Mathematics.
Copyright 2014 STATOIL
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/polymer/fullyimplicit/IncompPropsAdBasic.hpp>
#include <opm/core/utility/Units.hpp>
@ -9,6 +27,7 @@
namespace Opm
{
/// Constructor.
IncompPropsAdBasic::IncompPropsAdBasic(const parameter::ParameterGroup& param,
const int dim,
const int num_cells)
@ -28,7 +47,7 @@ namespace Opm
pvt_.mu(1, 0, 0, &viscosity_[0]);
}
/// Constructor.
IncompPropsAdBasic::IncompPropsAdBasic(const int num_phases,
const SaturationPropsBasic::RelPermFunc& relpermfunc,
const std::vector<double>& rho,
@ -48,10 +67,16 @@ namespace Opm
viscosity_.resize(pvt_.numPhases());
pvt_.mu(1, 0, 0, &viscosity_[0]);
}
/// Destructor.
IncompPropsAdBasic::~IncompPropsAdBasic()
{
}
////////////////////////////
// Rock interface //
////////////////////////////
/// \return D, the number of spatial dimensions.
int IncompPropsAdBasic::numDimensions() const
{
@ -79,7 +104,9 @@ namespace Opm
}
// ---- Fluid interface ----
////////////////////////////
// Fluid interface //
////////////////////////////
/// \return P, the number of phases (also the number of components).
int IncompPropsAdBasic::numPhases() const
@ -92,21 +119,34 @@ namespace Opm
{
return &viscosity_[0];
}
/// \return Array of P density values.
/// Densities of fluid phases at reservoir conditions.
/// \return Array of P density values.
const double* IncompPropsAdBasic::density() const
{
return pvt_.surfaceDensities();
}
/// Densities of fluid phases at surface conditions.
/// \return Array of P density values.
const double* IncompPropsAdBasic::surfaceDensity() const
{
return pvt_.surfaceDensities();
}
typedef IncompPropsAdBasic::ADB ADB;
typedef IncompPropsAdBasic::V V;
typedef Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> Block;
typedef std::vector<int> Cells;
// ------ Relative permeability ------
/// Relative permeabilities for all phases.
/// \param[in] sw Array of n water saturation values.
/// \param[in] so Array of n oil saturation values.
/// \param[in] cells Array of n cell indices to be associated with the saturation values.
/// \return An std::vector with 2 elements, each an array of n relperm values,
/// containing krw, kro.
std::vector<V>
IncompPropsAdBasic::relperm(const V& sw,
const V& so,
@ -130,6 +170,12 @@ namespace Opm
return relperms;
}
/// Relative permeabilities for all phases.
/// \param[in] sw Array of n water saturation values.
/// \param[in] so Array of n oil saturation values.
/// \param[in] cells Array of n cell indices to be associated with the saturation values.
/// \return An std::vector with 2 elements, each an array of n relperm values,
/// containing krw, kro.
std::vector<ADB>
IncompPropsAdBasic::relperm(const ADB& sw,
const ADB& so,
@ -169,4 +215,5 @@ namespace Opm
}
return relperms;
}
}
} //namespace Opm

127
opm/polymer/fullyimplicit/IncompPropsAdBasic.hpp
View File

@ -1,4 +1,22 @@
/**/
/*
Copyright 2014 SINTEF ICT, Applied Mathematics.
Copyright 2014 STATOIL
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_INCOMPPROPSADBASIC_HEADER_INCLUDED
#define OPM_INCOMPPROPSADBASIC_HEADER_INCLUDED
@ -8,47 +26,104 @@
#include <opm/core/props/pvt/PvtPropertiesBasic.hpp>
#include <opm/core/props/satfunc/SaturationPropsBasic.hpp>
#include <opm/polymer/fullyimplicit/AutoDiffBlock.hpp>
namespace Opm
{
/// This class implements the AD-adapted fluid interface for
/// two-phase oil-water. It requires an input deck from which it
/// reads all relevant property data.
///
/// Most methods are available in two overloaded versions, one
/// taking a constant vector and returning the same, and one
/// taking an AD type and returning the same. Derivatives are not
/// returned separately by any method, only implicitly with the AD
/// version of the methods.
class IncompPropsAdBasic : public IncompPropsAdInterface
{
public:
/// Constructor.
IncompPropsAdBasic(const parameter::ParameterGroup& param,
const int dim,
const int num_cells);
/// Constructor.
IncompPropsAdBasic(const int num_phases,
const SaturationPropsBasic::RelPermFunc& relpermfunc,
const std::vector<double>& rho,
const std::vector<double>& mu,
const double porosity,
const double permeability,
const int dim,
const int num_cells);
const SaturationPropsBasic::RelPermFunc& relpermfunc,
const std::vector<double>& rho,
const std::vector<double>& mu,
const double porosity,
const double permeability,
const int dim,
const int num_cells);
/// Destructor.
~IncompPropsAdBasic();
~IncompPropsAdBasic();
int numDimensions() const;
int numCells() const;
const double* porosity() const;
const double* permeability() const;
////////////////////////////
// Rock interface //
////////////////////////////
/// \return D, the number of spatial dimensions.
int numDimensions() const;
/// \return N, the number of cells.
int numCells() const;
/// \return Array of N porosity values.
const double* porosity() const;
/// \return Array of ND^2 permeability values.
/// The D^2 permeability values for a cell are organized as a matrix,
/// which is symmetric (so ordering does not matter).
const double* permeability() const;
int numPhases() const;
const double* viscosity() const;
const double* density() const;
const double* surfaceDensity() const;
////////////////////////////
// Fluid interface //
////////////////////////////
typedef AutoDiffBlock<double> ADB;
typedef ADB::V V;
std::vector<V> relperm(const V& sw,
const V& so,
const std::vector<int>& cells) const;
std::vector<ADB> relperm(const ADB& sw,
const ADB& so,
const std::vector<int>& cells) const;
typedef AutoDiffBlock<double> ADB;
typedef ADB::V V;
/// \return P, Number of active phases (also the number of components).
int numPhases() const;
/// \return Array of P viscosity values.
const double* viscosity() const;
/// Densities of fluid phases at reservoir conditions.
/// \return Array of P density values.
const double* density() const;
/// Densities of fluid phases at surface conditions.
/// \return Array of P density values.
const double* surfaceDensity() const;
// ------ Relative permeability ------
/// Relative permeabilities for all phases.
/// \param[in] sw Array of n water saturation values.
/// \param[in] so Array of n oil saturation values.
/// \param[in] cells Array of n cell indices to be associated with the saturation values.
/// \return An std::vector with 2 elements, each an array of n relperm values,
/// containing krw, kro.
std::vector<V> relperm(const V& sw,
const V& so,
const std::vector<int>& cells) const;
/// Relative permeabilities for all phases.
/// \param[in] sw Array of n water saturation values.
/// \param[in] so Array of n oil saturation values.
/// \param[in] cells Array of n cell indices to be associated with the saturation values.
/// \return An std::vector with 2 elements, each an array of n relperm values,
/// containing krw, kro.
std::vector<ADB> relperm(const ADB& sw,
const ADB& so,
const std::vector<int>& cells) const;
private:
RockBasic rock_;
PvtPropertiesBasic pvt_;
SaturationPropsBasic satprops_;
std::vector<double> viscosity_;
};
}
} //namespace Opm
#endif // OPM_INCOMPPROPSADBASIC_HEADER_INCLUDED

74
opm/polymer/fullyimplicit/IncompPropsAdFromDeck.cpp
View File

@ -1,4 +1,23 @@
/**/
/*
Copyright 2014 SINTEF ICT, Applied Mathematics.
Copyright 2014 STATOIL.
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/polymer/fullyimplicit/IncompPropsAdFromDeck.hpp>
#include <opm/polymer/fullyimplicit/AutoDiffHelpers.hpp>
#include <opm/core/utility/Units.hpp>
@ -6,7 +25,8 @@
#include <iostream>
namespace Opm
{
{
/// Constructor wrapping an opm-core two-phase interface.
IncompPropsAdFromDeck::IncompPropsAdFromDeck(const EclipseGridParser& deck,
const UnstructuredGrid& grid)
{
@ -22,44 +42,62 @@ namespace Opm
IncompPropsAdFromDeck::~IncompPropsAdFromDeck()
{
}
// rock interface
////////////////////////////
// Rock interface //
////////////////////////////
/// \return D, the number of spatial dimensions.
int IncompPropsAdFromDeck::numDimensions() const
{
return rock_.numDimensions();
}
/// \return N, the number of cells.
int IncompPropsAdFromDeck::numCells() const
{
return rock_.numCells();
}
/// \return Array of N porosity values.
const double* IncompPropsAdFromDeck::porosity() const
{
return rock_.porosity();
}
/// \return Array of ND^2 permeability values.
/// The D^2 permeability values for a cell are organized as a matrix,
/// which is symmetric (so ordering does not matter).
const double* IncompPropsAdFromDeck::permeability() const
{
return rock_.permeability();
}
// fluid interface
////////////////////////////
// Fluid interface //
////////////////////////////
/// \return P, Number of active phases (also the number of components).
int IncompPropsAdFromDeck::numPhases() const
{
return pvt_.numPhases();
}
/// \return Array of P viscosity values.
const double* IncompPropsAdFromDeck::viscosity() const
{
return pvt_.viscosity();
}
///Densities of fluid phases at reservoir conditions.
/// \return Array of P density values.
const double* IncompPropsAdFromDeck::density() const
{
return pvt_.reservoirDensities();
}
/// Densities of fluid phases at surface conditions.
/// \return Array of P density values.
const double* IncompPropsAdFromDeck::surfaceDensity() const
{
return pvt_.surfaceDensities();
@ -71,6 +109,14 @@ namespace Opm
typedef Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> Block;
// ------ Relative permeability ------
/// Relative permeabilities for all phases.
/// \param[in] sw Array of n water saturation values.
/// \param[in] so Array of n oil saturation values.
/// \param[in] cells Array of n cell indices to be associated with the saturation values.
/// \return An std::vector with 2 elements, each an array of n relperm values,
/// containing krw, kro.
std::vector<V>
IncompPropsAdFromDeck::relperm(const V& sw,
const V& so,
@ -94,6 +140,12 @@ namespace Opm
/// Relative permeabilities for all phases.
/// \param[in] sw Array of n water saturation values.
/// \param[in] so Array of n oil saturation values.
/// \param[in] cells Array of n cell indices to be associated with the saturation values.
/// \return An std::vector with 2 elements, each an array of n relperm values,
/// containing krw, kro.
std::vector<ADB>
IncompPropsAdFromDeck::relperm(const ADB& sw,
const ADB& so,
@ -132,6 +184,14 @@ namespace Opm
}
return relperms;
}
/// Capillary pressure for all phases.
/// \param[in] sw Array of n water saturation values.
/// \param[in] so Array of n oil saturation values.
/// \param[in] cells Array of n cell indices to be associated with the saturation values.
/// \return An std::vector with 2 elements, each an array of n capillary pressure values,
/// containing the offsets for each p_o, p_w. The capillary pressure between
/// two arbitrary phases alpha and beta is then given as p_alpha - p_beta.
std::vector<ADB>
IncompPropsAdFromDeck::capPress(const ADB& sw,
const ADB& so,

90
opm/polymer/fullyimplicit/IncompPropsAdFromDeck.hpp
View File

@ -1,6 +1,26 @@
/**/
/*
Copyright 2014 SINTEF ICT, Applied Mathematics.
Copyright 2014 STATOIL
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_INCOMPPROPSADFROMDECK_HEADER_INCLUDED
#define OPM_INCOMPPROPSADFROMDECK_HEADER_INCLUDED
#include <opm/polymer/fullyimplicit/IncompPropsAdInterface.hpp>
#include <opm/polymer/fullyimplicit/AutoDiffBlock.hpp>
#include <opm/core/props/rock/RockFromDeck.hpp>
@ -14,34 +34,84 @@
struct UnstructuredGrid;
namespace Opm
{
/// This class implements the AD-adapted fluid interface for
/// two-phase oil-water. It requires an input deck from which it
/// reads all relevant property data.
///
/// Most methods are available in two overloaded versions, one
/// taking a constant vector and returning the same, and one
/// taking an AD type and returning the same. Derivatives are not
/// returned separately by any method, only implicitly with the AD
/// version of the methods.
class IncompPropsAdFromDeck : public IncompPropsAdInterface
{
public:
/// Constructor wrapping an opm-core two-phase interface.
IncompPropsAdFromDeck(const EclipseGridParser& deck,
const UnstructuredGrid& grid);
const UnstructuredGrid& grid);
/// Destructor.
~IncompPropsAdFromDeck();
//--Rock interface--
////////////////////////////
// Rock interface //
////////////////////////////
/// \return D, the number of spatial dimensions.
int numDimensions() const;
/// \return N, the number of cells.
int numCells() const;
/// \return Array of N porosity values.
const double* porosity() const;
/// \return Array of ND^2 permeability values.
/// The D^2 permeability values for a cell are organized as a matrix,
/// which is symmetric (so ordering does not matter).
const double* permeability() const;
// --Fluid interface--
int numPhases() const;
const double* viscosity() const;
const double* density() const;
const double* surfaceDensity() const;
////////////////////////////
// Fluid interface //
////////////////////////////
typedef AutoDiffBlock<double> ADB;
typedef ADB::V V;
typedef std::vector<int> Cells;
/// \return P, Number of active phases (also the number of components).
int numPhases() const;
/// \return Array of P viscosity values.
const double* viscosity() const;
/// Densities of fluid phases at reservoir conditions.
/// \return Array of P density values.
const double* density() const;
/// Densities of fluid phases at surface conditions.
/// \return Array of P density values.
const double* surfaceDensity() const;
// ------ Relative permeability ------
/// Relative permeabilities for all phases.
/// \param[in] sw Array of n water saturation values.
/// \param[in] so Array of n oil saturation values.
/// \param[in] cells Array of n cell indices to be associated with the saturation values.
/// \return An std::vector with 2 elements, each an array of n relperm values,
/// containing krw, kro.
std::vector<V> relperm(const V& sw,
const V& so,
const Cells& cells) const;
/// Relative permeabilities for all phases.
/// \param[in] sw Array of n water saturation values.
/// \param[in] so Array of n oil saturation values.
/// \param[in] cells Array of n cell indices to be associated with the saturation values.
/// \return An std::vector with 2 elements, each an array of n relperm values,
/// containing krw, kro.
std::vector<ADB> relperm(const ADB& sw,
const ADB& so,
const Cells& cells) const;
@ -50,8 +120,8 @@ namespace Opm
/// \param[in] sw Array of n water saturation values.
/// \param[in] so Array of n oil saturation values.
/// \param[in] cells Array of n cell indices to be associated with the saturation values.
/// \return An std::vector with 3 elements, each an array of n capillary pressure values,
/// containing the offsets for each p_g, p_o, p_w. The capillary pressure between
/// \return An std::vector with 2 elements, each an array of n capillary pressure values,
/// containing the offsets for each p_o, p_w. The capillary pressure between
/// two arbitrary phases alpha and beta is then given as p_alpha - p_beta.
std::vector<ADB> capPress(const ADB& sw,
const ADB& so,

19
opm/polymer/fullyimplicit/IncompPropsAdInterface.cpp
View File

@ -1,3 +1,22 @@
/*
Copyright 2014 SINTEF ICT, Applied Mathematics.
Copyright 2014 STATOIL.
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/polymer/fullyimplicit/IncompPropsAdInterface.hpp>

40
opm/polymer/fullyimplicit/IncompPropsAdInterface.hpp
View File

@ -1,9 +1,26 @@
/*
Copyright 2014 SINTEF ICT, Applied Mathematics.
Copyright 2014 STATOIL.
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_INCOMPPROPSADINTERFACE_HEADER_INCLUDED
#define OPM_INCOMPPROPSADINTERFACE_HEADER_INCLUDED
#include <opm/polymer/fullyimplicit/AutoDiffBlock.hpp>
namespace Opm
@ -11,6 +28,7 @@ namespace Opm
class IncompPropsAdInterface
{
public:
/// Virtual destructor for inheritance.
virtual ~IncompPropsAdInterface();
////////////////////////////
@ -31,7 +49,6 @@ namespace Opm
/// which is symmetric (so ordering does not matter).
virtual const double* permeability() const = 0;
////////////////////////////
// Fluid interface //
////////////////////////////
@ -41,24 +58,21 @@ namespace Opm
typedef ADB::M M;
typedef std::vector<int> Cells;
/// \return Number of active phases (also the number of components).
/// \return P, Number of active phases (also the number of components).
virtual int numPhases() const = 0;
// ------ Density ------
/// Densities of stock components at surface conditions.
/// \return Array of 2 density values.
/// \return Array of P density values.
virtual const double* surfaceDensity() const = 0;
// ------ Viscosity ------
/// Viscosity of stock components at surface conditions.
/// \return Array of 2 viscosity values.
/// \return Array of P viscosity values.
virtual const double* viscosity() const = 0;
// ------ Relative permeability ------
/// Relative permeabilities for all phases.
@ -66,7 +80,7 @@ namespace Opm
/// \param[in] so Array of n oil saturation values.
/// \param[in] cells Array of n cell indices to be associated with the saturation values.
/// \return An std::vector with 2 elements, each an array of n relperm values,
/// containing krw, kro. Use PhaseIndex for indexing into the result.
/// containing krw, kro.
virtual
std::vector<V> relperm(const V& sw,
const V& so,
@ -77,24 +91,24 @@ namespace Opm
/// \param[in] so Array of n oil saturation values.
/// \param[in] cells Array of n cell indices to be associated with the saturation values.
/// \return An std::vector with 2 elements, each an array of n relperm values,
/// containing krw, kro. Use PhaseIndex for indexing into the result.
/// containing krw, kro.
virtual
std::vector<ADB> relperm(const ADB& sw,
const ADB& so,
const Cells& cells) const = 0;
/// Capillary pressure for all phases.
/// \param[in] sw Array of n water saturation values.
/// \param[in] so Array of n oil saturation values.
/// \param[in] cells Array of n cell indices to be associated with the saturation values.
/// \return An std::vector with 3 elements, each an array of n capillary pressure values,
/// containing the offsets for each p_g, p_o, p_w. The capillary pressure between
/// \return An std::vector with 2 elements, each an array of n capillary pressure values,
/// containing the offsets for each p_o, p_w. The capillary pressure between
/// two arbitrary phases alpha and beta is then given as p_alpha - p_beta.
virtual
std::vector<ADB> capPress(const ADB& sw,
const ADB& so,
const Cells& cells) const = 0;
};
} // namespace Opm

182
opm/polymer/fullyimplicit/PolymerPropsAd.cpp
View File

@ -7,161 +7,16 @@
namespace Opm {
typedef PolymerPropsAd::ADB ADB;
typedef PolymerPropsAd::V V;
/*
PolymerPropsAd::PolymerPropsAd()
{
}
PolymerPropsAd::PolymerPropsAd(const int num_cells,
const double c_max,
const double mix_param,
const double std::vector<double>& c_vals_visc,
const double std::vector<double>& visc_mult_vals)
: num_cells_ (num_cells)
, c_max_ (c_max)
, mix_param_(mix_param)
, c_vals_visc_ (c_vals_visc)
, visc_mult_vals_ (visc_mult_vals)
{
}
double PolymerPropsAd::num_cells() const
{
return num__cells_;
}
double PolymerPropsAd::cMax() const
{
return c_max_;
}
double PolymerPropsAd::mixParam() const
{
return mix_param_;
}
V PolymerPropsAd::muM(const V& c,
const double* visc) const
{
const int nc = num_cells();
assert(nc == c.size());
std::vector<double> mu(nc);
for (int i = 0; i < nc; ++i) {
mu[i] = Opm::linearInterpolation(c_vals_visc_, visc_mult_vals_, c(i));
}
const V muM = Eigen::Map<const V>(&mu[0], nc);
const double mu_w = visc[0];
return muM * mu_w;
}
ADB PolymerPropsAd::muM(const ADB& c,
const double* visc) const
{
const int nc = num_cells();
assert(nc == c.size());
V mu_m = muM(c.value());
std::vector<double> dmu_dc(nc);
for (int i = 0; i < nc; ++i) {
dmu_dc[i] = Opm::linearInterpolationDerivative(c_vals_visc_, visc_mult_vals_, c.value()(i));
}
const V dmu = Eigen::Map<const V>(&dmu_dc[0], nc);
ADB::M dmu_diag = spdiag(dmu);
const int num_blocks = c.numBlocks();
std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) {
jacs[block] = dmu_diag * c.derivative()[block];
}
const double mu_w = visc[0]
return ADB::function(mu_m, jacs) * mu_w;
}
V PolymerPropsAd::ToddLongstaff(const double mix_param,
const V& muM,
const V& mu) const
{
const int nc = num_cells();
assert(nc == muM.size());
std::vector<double> mueff(nc);
const double omega = mix_param;
for (int i = 0; i < nc; ++i) {
mueff[i] = std::pow(muM(i),omega) * std::pow(mu(i), 1. - omega);
}
const V muEff = Eigen::Map<const V>(&mueff[0], nc);
return muEff;
}
ADB PolymerPropsAd::ToddLongstaff(const double mix_param,
const ADB& muM,
const ADB& mu) const
{
const int nc = num_cells();
}
V PolymerPropsAd::muPolyEff(const double mix_param,
const V& muM,
const V& muPoly) const
{
return ToddLongstaff(mix_param, muM, muPoly);
}
V PolymerPropsAd::muWatEff(const double mix_param,
const std::vecotr<double>& c_max,
const V& c,
const V& muM,
const V& muWat,
const V& muPolyEff) const
{
const int nc = num_cells();
assert(nc == c.size());
V muWate = ToddLongstaff(mix_param, muM, muWat);
// V cmax = V::Constant(nc, 1, c_max);
const V cmax = Eigen::Map<const V>(&c_max[0], nc);
const V one = V::Ones(nc, 1);
V inv_muWatEff = (one - c / cmax) / muWate + c / cmax / muPolyEff;
V muWatEff = one / inv_muWatEff;
return muWatEff;
}
*/
double
PolymerPropsAd::rockDensity() const
{
@ -171,17 +26,25 @@ namespace Opm {
double
PolymerPropsAd::deadPoreVol() const
{
return polymer_props_.deadPoreVol();
}
double
PolymerPropsAd::cMax() const
{
return polymer_props_.cMax();
}
PolymerPropsAd::PolymerPropsAd(const PolymerProperties& polymer_props)
@ -192,6 +55,7 @@ namespace Opm {
PolymerPropsAd::~PolymerPropsAd()
{
}
@ -199,6 +63,7 @@ namespace Opm {
V PolymerPropsAd::effectiveInvWaterVisc(const V& c,
const double* visc) const
{
@ -217,7 +82,6 @@ namespace Opm {
ADB PolymerPropsAd::effectiveInvWaterVisc(const ADB& c,
const double* visc) const
{
@ -242,6 +106,7 @@ namespace Opm {
V PolymerPropsAd::polymerWaterVelocityRatio(const V& c) const
{
const int nc = c.size();
@ -288,6 +153,8 @@ namespace Opm {
V PolymerPropsAd::adsorption(const V& c, const V& cmax_cells) const
{
const int nc = c.size();
@ -303,6 +170,10 @@ namespace Opm {
return ads;
}
ADB PolymerPropsAd::adsorption(const ADB& c, const ADB& cmax_cells) const
{
const int nc = c.value().size();
@ -329,6 +200,9 @@ namespace Opm {
}
V
PolymerPropsAd::effectiveRelPerm(const V& c,
const V& cmax_cells,
@ -348,6 +222,9 @@ namespace Opm {
}
ADB
PolymerPropsAd::effectiveRelPerm(const ADB& c,
const ADB& cmax_cells,
@ -355,9 +232,7 @@ namespace Opm {
const ADB& sw) const
{
const int nc = c.value().size();
V one = V::Ones(nc);
ADB ads = adsorption(c, cmax_cells);
V krw_eff = effectiveRelPerm(c.value(), cmax_cells.value(), krw.value());
@ -365,15 +240,14 @@ namespace Opm {
double res_factor = polymer_props_.resFactor();
double factor = (res_factor - 1.) / max_ads;
ADB rk = one + ads * factor;
// ADB dkrw_ds = krw / rk.value();
ADB dkrw_ds = krw / rk;
// ADB dkrw_dc = -krw.value() / (rk.value() * rk.value()) * ads * factor;
ADB dkrw_dc = -factor * krw / (rk * rk) * ads ;
const int num_blocks = c.numBlocks();
std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) {
jacs[block] = dkrw_ds.derivative()[block] * sw.derivative()[block] + dkrw_dc.derivative()[block] * c.derivative()[block];
jacs[block] = dkrw_ds.derivative()[block] * sw.derivative()[block]
+ dkrw_dc.derivative()[block] * c.derivative()[block];
}
return ADB::function(krw_eff, jacs);

130
opm/polymer/fullyimplicit/PolymerPropsAd.hpp
View File

@ -1,120 +1,112 @@
/*
Copyright 2014 SINTEF ICT, Applied Mathematics.
Copyright 2014 STATOIL.
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_POLYMERPROPSAD_HEADED_INLCUDED
#define OPM_POLYMERPROPSAD_HEADED_INLCUDED
#include <cmath>
#include <vector>
#include <opm/polymer/fullyimplicit/AutoDiffBlock.hpp>
#include <opm/polymer/fullyimplicit/AutoDiffHelpers.hpp>
#include <opm/polymer/PolymerProperties.hpp>
namespace Opm {
namespace Opm {
class PolymerPropsAd
{
public:
/* PolymerPropsAd(const int num_cells;
const double mix_param,
const std::vector<double>& c_max,
const std::vector<double>& c_vals_visc,
const std::vector<double>& visc_mult_vals);
double num_cells() const;
double cMax() const;
double mixParam() const;
typedef AutoDiffBlock<double> ADB;
typedef ADB::V V;
V muM(const V& c,
const double* visc) const;
ADB muM(const ADB& c
const double* visc) const;
V ToddLongstaff(const double mix_param,
const V& muM,
const V& mu) const;
ADB ToddLongstaff(const double mix_param,
const ADB& muM,
const ADB& mu) const;
V muPolyEff(const double mix_param,
const V& muM,
const V& muPoly) const;
ADB muPolyEff(const double mix_param,
const ADB& muM,
const ADB& muPoly) const;
V muWatEff(const double mix_param,
const std::vector<double>& c_max,
const V& c,
const V& muM,
const V& muWat,
const V& muPolyEff) const;
ADB muWatEff(const double mix_param,
const std::vector<double>& c_max,
const ADB& c,
const ADB& muM,
const ADB& muWat,
const ADB& muPolyEff) const;
*/
/// \return Reference rock density.
double rockDensity() const;
/// \return The value of dead pore volume.
double deadPoreVol() const;
/// \return The max concentration injected.
double cMax() const;
typedef AutoDiffBlock<double> ADB;
typedef ADB::V V;
/// Constructor wrapping a polymer props.
PolymerPropsAd(const PolymerProperties& polymer_props);
/// Destructor.
~PolymerPropsAd();
/// \param[in] c Array of n polymer concentraion values.
/// \param[in] visc Array of 2 viscosity value.
/// \return value of inverse effective water viscosity.
V
effectiveInvWaterVisc(const V& c,const double* visc) const;
/// \param[in] c Array of n polymer concentraion values.
/// \param[in] visc Array of 2 viscosity value
/// \return value of inverse effective water viscosity.
ADB
effectiveInvWaterVisc(const ADB& c,const double* visc) const;
/// \param[in] c Array of n polymer concentraion values.
/// \return Array of n mc values, here mc means m(c) * c.
V
polymerWaterVelocityRatio(const V& c) const;
/// \param[in] c Array of n polymer concentraion values.
/// \return Array of n mc values, here mc means m(c) * c.
ADB
polymerWaterVelocityRatio(const ADB& c) const;
/// \param[in] c Array of n polymer concentraion values.
/// \param[in] cmax_cells Array of n polymer concentraion values
/// that the cell experienced.
/// \return Array of n adsorption values.
V
adsorption(const V& c, const V& cmax_cells) const;
/// \param[in] c Array of n polymer concentraion values.
/// \param[in] cmax_cells Array of n polymer concentraion values
/// that the cell experienced.
/// \return Array of n adsorption values.
ADB
adsorption(const ADB& c, const ADB& cmax_cells) const;
/// \param[in] c Array of n polymer concentraion values.
/// \param[in] cmax_cells Array of n polymer concentraion values
/// that the cell experienced.
/// \param[in] relperm Array of n relative water relperm values.
/// \return Array of n adsorption values.
V
effectiveRelPerm(const V& c, const V& cmax_cells, const V& relperm) const;
/// \param[in] c Array of n polymer concentraion values.
/// \param[in] cmax_cells Array of n polymer concentraion values
/// that the cell experienced.
/// \param[in] relperm Array of n relative water relperm values.
/// \return Array of n adsorption values.
ADB
effectiveRelPerm(const ADB& c, const ADB& cmax_cells, const ADB& krw, const ADB& sw) const;
private:
const PolymerProperties& polymer_props_;
};
}
} //namespace Opm
#endif// OPM_POLYMERPROPSAD_HEADED_INLCUDED

18
opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.cpp
View File

@ -242,7 +242,7 @@ namespace Opm
}
#endif
static void outputWaterCut(const Opm::Watercut& watercut,
const std::string& output_dir)
const std::string& output_dir)
{
// Write water cut curve.
std::string fname = output_dir + "/watercut.txt";
@ -270,14 +270,14 @@ namespace Opm
// \TODO: Treat bcs.
SimulatorFullyImplicitCompressiblePolymer::Impl::Impl(const parameter::ParameterGroup& param,
const UnstructuredGrid& grid,
const BlackoilPropsAdInterface& props,
const PolymerPropsAd& polymer_props,
const RockCompressibility* rock_comp_props,
WellsManager& wells_manager,
PolymerInflowInterface& polymer_inflow,
LinearSolverInterface& linsolver,
const double* gravity)
const UnstructuredGrid& grid,
const BlackoilPropsAdInterface& props,
const PolymerPropsAd& polymer_props,
const RockCompressibility* rock_comp_props,
WellsManager& wells_manager,
PolymerInflowInterface& polymer_inflow,
LinearSolverInterface& linsolver,
const double* gravity)
: grid_(grid),
props_(props),
polymer_props_(polymer_props),

23
opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.hpp
View File

@ -1,5 +1,6 @@
/*
Copyright 2013 SINTEF ICT, Applied Mathematics.
Copyright 2014 SINTEF ICT, Applied Mathematics.
Copyright 2014 STATOIL.
This file is part of the Open Porous Media project (OPM).
@ -62,19 +63,21 @@ namespace Opm
///
/// \param[in] grid grid data structure
/// \param[in] props fluid and rock properties
/// \param[in] rock_comp_props if non-null, rock compressibility properties
/// \param[in] polymer_props polymer properties
/// \param[in] rock_comp_props if non-null, rock compressibility properties
/// \param[in] well_manager well manager, may manage no (null) wells
/// \param[in] polymer_inflow polymer influx.
/// \param[in] linsolver linear solver
/// \param[in] gravity if non-null, gravity vector
SimulatorFullyImplicitCompressiblePolymer(const parameter::ParameterGroup& param,
const UnstructuredGrid& grid,
const BlackoilPropsAdInterface& props,
const PolymerPropsAd& polymer_props,
const RockCompressibility* rock_comp_props,
WellsManager& wells_manager,
PolymerInflowInterface& polymer_inflow,
LinearSolverInterface& linsolver,
const double* gravity);
const UnstructuredGrid& grid,
const BlackoilPropsAdInterface& props,
const PolymerPropsAd& polymer_props,
const RockCompressibility* rock_comp_props,
WellsManager& wells_manager,
PolymerInflowInterface& polymer_inflow,
LinearSolverInterface& linsolver,
const double* gravity);
/// Run the simulation.
/// This will run succesive timesteps until timer.done() is true. It will

18
opm/polymer/fullyimplicit/SimulatorFullyImplicitTwophasePolymer.hpp
View File

@ -60,15 +60,19 @@ namespace Opm
///
/// \param[in] grid grid data structure
/// \param[in] props fluid and rock properties
/// \param[in] polymer_props polymer properties
/// \param[in] linsolver linear solver
/// \param[in] well_manager well manager, may manage no (null) wells
/// \param[in] polymer_inflow polymer influx.
/// \param[in] gravity if non-null, gravity vector
SimulatorFullyImplicitTwophasePolymer(const parameter::ParameterGroup& param,
const UnstructuredGrid& grid,
const IncompPropsAdInterface& props,
const PolymerPropsAd& polymer_props,
LinearSolverInterface& linsolver,
WellsManager& wells_manager,
PolymerInflowInterface& polymer_inflow,
const double* gravity);
const UnstructuredGrid& grid,
const IncompPropsAdInterface& props,
const PolymerPropsAd& polymer_props,
LinearSolverInterface& linsolver,
WellsManager& wells_manager,
PolymerInflowInterface& polymer_inflow,
const double* gravity);
/// Run the simulation.
/// This will run succesive timesteps until timer.done() is true. It will