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

Simulator for fullyimplicit two phase flow with polymer

Browse Source
This commit is contained in:
Liu Ming
2013-12-09 20:33:05 +08:00
parent 080116c66b
commit e226b5d95e
5 changed files with 111 additions and 65 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -7,6 +7,7 @@
#include <opm/polymer/fullyimplicit/AutoDiffHelpers.hpp> #include <opm/polymer/fullyimplicit/AutoDiffHelpers.hpp>
#include <opm/polymer/fullyimplicit/IncompPropsAdInterface.hpp> #include <opm/polymer/fullyimplicit/IncompPropsAdInterface.hpp>
#include <opm/polymer/PolymerProperties.hpp> #include <opm/polymer/PolymerProperties.hpp>
#include <opm/polymer/PolymerState.hpp>
#include <opm/polymer/fullyimplicit/PolymerPropsAd.hpp> #include <opm/polymer/fullyimplicit/PolymerPropsAd.hpp>
#include <opm/core/grid.h> #include <opm/core/grid.h>
@@ -56,15 +57,13 @@ typedef Eigen::Array<double,
FullyImplicitTwoPhaseSolver:: FullyImplicitTwophasePolymerSolver::
FullyImplicitTwoPhaseSolver(const UnstructuredGrid& grid, FullyImplicitTwophasePolymerSolver(const UnstructuredGrid& grid,
const IncompPropsAdInterface& fluid, const IncompPropsAdInterface& fluid,
const PolymerProperties& polymer_props,
const PolymerPropsAd& polymer_props_ad, const PolymerPropsAd& polymer_props_ad,
const LinearSolverInterface& linsolver) const LinearSolverInterface& linsolver)
: grid_ (grid) : grid_ (grid)
, fluid_(fluid) , fluid_(fluid)
, polymer_props_ (polymer_props)
, polymer_props_ad_ (polymer_props_ad_) , polymer_props_ad_ (polymer_props_ad_)
, linsolver_(linsolver) , linsolver_(linsolver)
, cells_ (buildAllCells(grid.number_of_cells)) , cells_ (buildAllCells(grid.number_of_cells))
@@ -78,10 +77,12 @@ typedef Eigen::Array<double,
void void
FullyImplicitTwoPhaseSolver:: FullyImplicitTwophasePolymerSolver::
step(const double dt, step(const double dt,
PolymerState& x, PolymerState& x,
const std::vector<double>& src) const std::vector<double>& src,
const std::vector<double>& polymer_inflow)
// const bool if_polymer_actived)
{ {
V pvol(grid_.number_of_cells); V pvol(grid_.number_of_cells);
@@ -99,7 +100,7 @@ typedef Eigen::Array<double,
const double rtol = 5.0e-8; const double rtol = 5.0e-8;
const int maxit = 15; const int maxit = 15;
assemble(pvdt, old_state, x, src); assemble(pvdt, old_state, x, src, polymer_inflow);//, if_polymer_actived);
const double r0 = residualNorm(); const double r0 = residualNorm();
int it = 0; int it = 0;
@@ -111,7 +112,7 @@ typedef Eigen::Array<double,
const V dx = solveJacobianSystem(); const V dx = solveJacobianSystem();
updateState(dx, x); updateState(dx, x);
assemble(pvdt, old_state, x, src); assemble(pvdt, old_state, x, src, polymer_inflow);//, if_polymer_actived);
const double r = residualNorm(); const double r = residualNorm();
@@ -132,7 +133,7 @@ typedef Eigen::Array<double,
FullyImplicitTwoPhaseSolver::SolutionState::SolutionState(const int np) FullyImplicitTwophasePolymerSolver::SolutionState::SolutionState(const int np)
: pressure ( ADB::null()) : pressure ( ADB::null())
, saturation (np, ADB::null()) , saturation (np, ADB::null())
, concentration ( ADB::null()) , concentration ( ADB::null())
@@ -143,8 +144,8 @@ typedef Eigen::Array<double,
FullyImplicitTwoPhaseSolver::SolutionState FullyImplicitTwophasePolymerSolver::SolutionState
FullyImplicitTwoPhaseSolver::constantState(const PolymerState& x) FullyImplicitTwophasePolymerSolver::constantState(const PolymerState& x)
{ {
const int nc = grid_.number_of_cells; const int nc = grid_.number_of_cells;
const int np = x.numPhases(); const int np = x.numPhases();
@@ -177,8 +178,8 @@ typedef Eigen::Array<double,
FullyImplicitTwoPhaseSolver::SolutionState FullyImplicitTwophasePolymerSolver::SolutionState
FullyImplicitTwoPhaseSolver::variableState(const PolymerState& x) FullyImplicitTwophasePolymerSolver::variableState(const PolymerState& x)
{ {
const int nc = grid_.number_of_cells; const int nc = grid_.number_of_cells;
const int np = x.numPhases(); const int np = x.numPhases();
@@ -197,7 +198,7 @@ typedef Eigen::Array<double,
const V sw = s_all.col(0); const V sw = s_all.col(0);
vars0.push_back(sw); vars0.push_back(sw);
// Initial saturation. // Initial concentration.
assert (not x.concentration().empty()); assert (not x.concentration().empty());
const V c = Eigen::Map<const V>(&x.concentration()[0], nc); const V c = Eigen::Map<const V>(&x.concentration()[0], nc);
vars0.push_back(c); vars0.push_back(c);
@@ -214,7 +215,7 @@ typedef Eigen::Array<double,
const std::vector<int>& bpat = vars[0].blockPattern(); const std::vector<int>& bpat = vars[0].blockPattern();
{ {
ADB so = ADB::constant(V::Ones(nc, 1), bpat); ADB so = ADB::constant(V::Ones(nc, 1), bpat);
ADB& sw = vars[ nextvar++ ]; ADB sw = vars[ nextvar++ ];
state.saturation[0] = sw; state.saturation[0] = sw;
so = so - sw; so = so - sw;
state.saturation[1] = so; state.saturation[1] = so;
@@ -233,11 +234,13 @@ typedef Eigen::Array<double,
void void
FullyImplicitTwoPhaseSolver:: FullyImplicitTwophasePolymerSolver::
assemble(const V& pvdt, assemble(const V& pvdt,
const SolutionState& old_state, const SolutionState& old_state,
const PolymerState& x , const PolymerState& x ,
const std::vector<double>& src) const std::vector<double>& src,
const std::vector<double>& polymer_inflow)
// const bool if_polymer_actived)
{ {
// Create the primary variables. // Create the primary variables.
const SolutionState state = variableState(x); const SolutionState state = variableState(x);
@@ -252,21 +255,42 @@ typedef Eigen::Array<double,
pvdt*(state.saturation[phase] - old_state.saturation[phase]) pvdt*(state.saturation[phase] - old_state.saturation[phase])
+ ops_.div*mflux - source; + ops_.div*mflux - source;
} }
// Mass balance equation for polymer
ADB polysrc = accumPolymerSource(kr, src);
ADB mc = computeMc();
ADB mflux = computeMassFlux(0, trans, kr, state);
residual_[2] = pvdt * (state.saturation[0] * state.concentration - old_state.saturation[0] * old_state.concentration) + ops_.div * state.concentration * mc * mflux - polysrc;
double amount = PolymerInjectedAmount(polymer_inflow);
bool use_polymer = (amount != 0);
std::cout << "\n Polymer Amount\n"
<< std::setprecision(9)
<< std::setw(18) << amount << std::endl;
const int nc = grid_.number_of_cells;
if (use_polymer) {
// Mass balance equation for polymer
const V src_polymer = Eigen::Map<const V>(&polymer_inflow[0], nc);
ADB mc = computeMc(state);
ADB mflux = computeMassFlux(0, trans, kr, state);
residual_[2] = pvdt * (state.saturation[0] * state.concentration
- old_state.saturation[0] * old_state.concentration)
+ ops_.div * state.concentration * mc * mflux - src_polymer;
} else {
residual_[2] = ADB::constant(V::Zero(nc));
}
}
double
FullyImplicitTwophasePolymerSolver::
PolymerInjectedAmount(const std::vector<double>& polymer_inflow) const
{
double amount = 0;
for (int i = 0; i < int(polymer_inflow.size()); ++i) {
amount += polymer_inflow[i];
}
return amount;
} }
ADB ADB
FullyImplicitTwoPhaseSolver::accumSource(const int phase, FullyImplicitTwophasePolymerSolver::accumSource(const int phase,
const std::vector<ADB>& kr, const std::vector<ADB>& kr,
const std::vector<double>& src) const const std::vector<double>& src) const
{ {
@@ -307,7 +331,7 @@ typedef Eigen::Array<double,
ADB ADB
FullyImplicitTwoPhaseSolver::computeFracFlow(int phase, FullyImplicitTwophasePolymerSolver::computeFracFlow(int phase,
const std::vector<ADB>& kr) const const std::vector<ADB>& kr) const
{ {
const double* mus = fluid_.viscosity(); const double* mus = fluid_.viscosity();
@@ -325,14 +349,14 @@ typedef Eigen::Array<double,
V V
FullyImplicitTwoPhaseSolver::solveJacobianSystem() const FullyImplicitTwophasePolymerSolver::solveJacobianSystem() const
{ {
const int np = fluid_.numPhases(); const int np = fluid_.numPhases();
if (np != 2) { if (np != 2) {
OPM_THROW(std::logic_error, "Only two-phase ok in FullyImplicitTwoPhaseSolver."); OPM_THROW(std::logic_error, "Only two-phase ok in FullyImplicitTwophasePolymerSolver.");
} }
ADB mass_phase_res = vertcat(residual_[0], residual_[1]); ADB mass_res = vertcat(residual_[0], residual_[1]);
ADB mass_res = collapseJacs(vertcat(mass_phase_res, residual_[2])); mass_res = collapseJacs(vertcat(mass_res, residual_[2]));
const Eigen::SparseMatrix<double, Eigen::RowMajor> matr = mass_res.derivative()[0]; const Eigen::SparseMatrix<double, Eigen::RowMajor> matr = mass_res.derivative()[0];
V dx(V::Zero(mass_res.size())); V dx(V::Zero(mass_res.size()));
@@ -352,7 +376,7 @@ typedef Eigen::Array<double,
void FullyImplicitTwoPhaseSolver::updateState(const V& dx, void FullyImplicitTwophasePolymerSolver::updateState(const V& dx,
PolymerState& state) const PolymerState& state) const
{ {
const int np = fluid_.numPhases(); const int np = fluid_.numPhases();
@@ -367,7 +391,8 @@ typedef Eigen::Array<double,
int varstart = nc; int varstart = nc;
const V dsw = subset(dx, Span(nc, 1, varstart)); const V dsw = subset(dx, Span(nc, 1, varstart));
varstart += dsw.size(); varstart += dsw.size();
const V dc = subset(dx, Span(nc ,1 varstart)); std::cout << dx.size() << " " << varstart << std::endl;
const V dc = subset(dx, Span(nc, 1, varstart));
varstart += dc.size(); varstart += dc.size();
assert(varstart == dx.size()); assert(varstart == dx.size());
@@ -405,7 +430,7 @@ typedef Eigen::Array<double,
std::vector<ADB> std::vector<ADB>
FullyImplicitTwoPhaseSolver::computeRelPerm(const SolutionState& state) const FullyImplicitTwophasePolymerSolver::computeRelPerm(const SolutionState& state) const
{ {
const ADB sw = state.saturation[0]; const ADB sw = state.saturation[0];
@@ -424,7 +449,7 @@ typedef Eigen::Array<double,
ADB ADB
FullyImplicitTwoPhaseSolver::computeMassFlux(const int phase , FullyImplicitTwophasePolymerSolver::computeMassFlux(const int phase ,
const V& trans, const V& trans,
const std::vector<ADB>& kr , const std::vector<ADB>& kr ,
const SolutionState& state ) const const SolutionState& state ) const
@@ -446,7 +471,7 @@ typedef Eigen::Array<double,
double double
FullyImplicitTwoPhaseSolver::residualNorm() const FullyImplicitTwophasePolymerSolver::residualNorm() const
{ {
double r = 0; double r = 0;
for (std::vector<ADB>::const_iterator for (std::vector<ADB>::const_iterator
@@ -465,7 +490,7 @@ typedef Eigen::Array<double,
V V
FullyImplicitTwoPhaseSolver::transmissibility() const FullyImplicitTwophasePolymerSolver::transmissibility() const
{ {
const V::Index nc = grid_.number_of_cells; const V::Index nc = grid_.number_of_cells;
V htrans(grid_.cell_facepos[nc]); V htrans(grid_.cell_facepos[nc]);
@@ -478,8 +503,12 @@ typedef Eigen::Array<double,
} }
ADB
FullyImplicitTwophasePolymerSolver::computeMc(const SolutionState& state) const
{
ADB c = state.concentration;
return polymer_props_ad_.polymerWaterVelocityRatio(c);
}
}//namespace Opm }//namespace Opm

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

@@ -17,18 +17,20 @@ namespace Opm {
class PolymerState; class PolymerState;
class FullyImplicitTwoPhaseSolver class FullyImplicitTwophasePolymerSolver
{ {
public: public:
FullyImplicitTwoPhaseSolver(const UnstructuredGrid& grid, FullyImplicitTwophasePolymerSolver(const UnstructuredGrid& grid,
const IncompPropsAdInterface& fluid, const IncompPropsAdInterface& fluid,
const PolymerProperties& polymer_props,
const PolymerPropsAd& polymer_props_ad, const PolymerPropsAd& polymer_props_ad,
const LinearSolverInterface& linsolver); const LinearSolverInterface& linsolver);
void step(const double dt, void step(const double dt,
PolymerState& state, PolymerState& state,
const std::vector<double>& src); const std::vector<double>& src,
const std::vector<double>& polymer_inflow
);
// const bool if_polymer_actived);
private: private:
typedef AutoDiffBlock<double> ADB; typedef AutoDiffBlock<double> ADB;
typedef ADB::V V; typedef ADB::V V;
@@ -46,7 +48,6 @@ namespace Opm {
}; };
const UnstructuredGrid& grid_; const UnstructuredGrid& grid_;
const IncompPropsAdInterface& fluid_; const IncompPropsAdInterface& fluid_;
const PolymerProperties& polymer_props_;
const PolymerPropsAd& polymer_props_ad_; const PolymerPropsAd& polymer_props_ad_;
const LinearSolverInterface& linsolver_; const LinearSolverInterface& linsolver_;
const std::vector<int> cells_; const std::vector<int> cells_;
@@ -62,7 +63,9 @@ namespace Opm {
assemble(const V& pvdt, assemble(const V& pvdt,
const SolutionState& old_state, const SolutionState& old_state,
const PolymerState& x, const PolymerState& x,
const std::vector<double>& src); const std::vector<double>& src,
const std::vector<double>& polymer_inflow);
// const bool if_polymer_actived);
V solveJacobianSystem() const; V solveJacobianSystem() const;
void updateState(const V& dx, void updateState(const V& dx,
PolymerState& x) const; PolymerState& x) const;
@@ -85,6 +88,9 @@ namespace Opm {
double double
residualNorm() const; residualNorm() const;
ADB
computeMc(const SolutionState& state) const;
ADB ADB
rockPorosity(const ADB& p) const; rockPorosity(const ADB& p) const;
ADB ADB
@@ -93,6 +99,8 @@ namespace Opm {
fluidDensity(const int phase) const; fluidDensity(const int phase) const;
ADB ADB
transMult(const ADB& p) const; transMult(const ADB& p) const;
double
PolymerInjectedAmount(const std::vector<double>& polymer_inflow) const;
}; };
} // namespace Opm } // namespace Opm
#endif// OPM_FULLYIMPLICITTWOPHASESOLVER_HEADER_INCLUDED #endif// OPM_FULLYIMPLICITTWOPHASESOLVER_HEADER_INCLUDED

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

@@ -2,7 +2,7 @@
#include <opm/polymer/fullyimplicit/IncompPropsAdBasic.hpp> #include <opm/polymer/fullyimplicit/IncompPropsAdBasic.hpp>
#include <opm/core/utility/Units.hpp> #include <opm/core/utility/Units.hpp>
#include <opm/polymer/fulluimplicit/AutoDiffBlock.hpp> #include <opm/polymer/fullyimplicit/AutoDiffBlock.hpp>
#include <opm/polymer/fullyimplicit/AutoDiffHelpers.hpp> #include <opm/polymer/fullyimplicit/AutoDiffHelpers.hpp>
#include <opm/core/utility/ErrorMacros.hpp> #include <opm/core/utility/ErrorMacros.hpp>
#include <iostream> #include <iostream>

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

@@ -182,7 +182,7 @@ namespace Opm {
const double* visc) const const double* visc) const
{ {
const int nc = c.size(); const int nc = c.size();
V inv_mu_w_eff(n); V inv_mu_w_eff(nc);
for (int i = 0; i < nc; ++i) { for (int i = 0; i < nc; ++i) {
double im = 0; double im = 0;
polymer_props_.effectiveInvVisc(c(i), visc, im); polymer_props_.effectiveInvVisc(c(i), visc, im);
@@ -198,11 +198,11 @@ namespace Opm {
ADB PolymerPropsAd::effectiveInvWaterVisc(const ADB& c, ADB PolymerPropsAd::effectiveInvWaterVisc(const ADB& c,
const double* visc) const double* visc) const
{ {
const int nc = c.size(); const int nc = c.size();
V inv_mu_w_eff(n); V inv_mu_w_eff(nc);
V dinv_mu_w_eff(n); V dinv_mu_w_eff(nc);
for (int i = 0; i < nc; ++i) { for (int i = 0; i < nc; ++i) {
double im = 0, dim = 0; double im = 0, dim = 0;
polymer_props_.effectiveInvViscWithDer(c.value()(i), visc, im, dim); polymer_props_.effectiveInvViscWithDer(c.value()(i), visc, im, dim);
@@ -221,11 +221,10 @@ namespace Opm {
V PolymerPropsAd::polymerWaterVelocityRatio(const V& c) const V PolymerPropsAd::polymerWaterVelocityRatio(const V& c) const
{ {
const int nc = c.size(); const int nc = c.size();
V mc(n); V mc(nc);
for (int i = 0; i < nc; ++i) { for (int i = 0; i < nc; ++i) {
double m = 0; double m = 0;
@@ -244,13 +243,13 @@ namespace Opm {
{ {
const int nc = c.size(); const int nc = c.size();
V mc(n); V mc(nc);
V dmc(n); V dmc(nc);
for (int i = 0; i < nc; ++i) { for (int i = 0; i < nc; ++i) {
double m = 0; double m = 0;
double dm = 0; double dm = 0;
polymer_props_.computeMcWithDer(c(i), m, dm); polymer_props_.computeMcWithDer(c.value()(i), m, dm);
mc(i) = m; mc(i) = m;
dmc(i) = dm; dmc(i) = dm;

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

@@ -6,7 +6,7 @@
#include <vector> #include <vector>
#include <opm/polymer/fullyimplicit/AutoDiffBlock.hpp> #include <opm/polymer/fullyimplicit/AutoDiffBlock.hpp>
#include <opm/polymer/fullyimplicit/AutoDiffHelpers.hpp> #include <opm/polymer/fullyimplicit/AutoDiffHelpers.hpp>
#include <opm/polymer/polymerProperties.hpp> #include <opm/polymer/PolymerProperties.hpp>
namespace Opm { namespace Opm {
@@ -65,15 +65,25 @@ namespace Opm {
*/ */
typedef AutoDiffBlock<double> ADB; typedef AutoDiffBlock<double> ADB;
typedef ADB::V V; typedef ADB::V V;
PolymerPropsAd(const PolymerPropperties& polymer_props);
~PolymerPropsAd();
V PolymerPropsAd::effectiveInvWaterVisc(const V& c,const double* visc) const;
ADB PolymerPropsAd::effectiveInvWaterVisc(const ADB& c,const double* visc) const;
V PolymerPropsAd::polymerWaterVelocityRatio(const V& c) const; PolymerPropsAd(const PolymerProperties& polymer_props);
ADB PolymerPropsAd::polymerWaterVelocityRatio(const ADB& c) const;
~PolymerPropsAd();
V
effectiveInvWaterVisc(const V& c,const double* visc) const;
ADB
effectiveInvWaterVisc(const ADB& c,const double* visc) const;
V
polymerWaterVelocityRatio(const V& c) const;
ADB
polymerWaterVelocityRatio(const ADB& c) const;
private: private:
const PolymerProperties polymer_props_; const PolymerProperties& polymer_props_;
}; };