fix bug when computing the maxconcentration that the cell experienced.

This commit is contained in:
Liu Ming 2014-01-10 17:20:08 +08:00
parent 139bfb7e62
commit 97f5c5ace5
4 changed files with 31 additions and 33 deletions

View File

@ -454,18 +454,18 @@ namespace {
ADB V
FullyImplicitCompressiblePolymerSolver:: FullyImplicitCompressiblePolymerSolver::
computeCmax(const ADB& c) const computeCmax(const PolymerBlackoilState& x) const
{ {
const int nc = c.value().size(); const int nc = c.value().size();
V cmax(nc); const V cmax = Eigen::Map<const V>(& x.maxconcentration()[0], nc);
const V c = Eigen::Map<const V>(& x.concentration()[0], nc);
for (int i = 0; i < nc; ++i) { for (int i = 0; i < nc; ++i) {
cmax(i) = (cmax(i) > c.value()(i)) ? cmax(i) : c.value()(i); cmax(i) = std::max(cmax(i), c(i));
} }
return ADB::constant(cmax, c.blockPattern()); return cmax;
} }
void void
@ -493,7 +493,8 @@ namespace {
const V trans = subset(geo_.transmissibility(), ops_.internal_faces); const V trans = subset(geo_.transmissibility(), ops_.internal_faces);
const std::vector<ADB> kr = computeRelPerm(state); const std::vector<ADB> kr = computeRelPerm(state);
const double dead_pore_vol = polymer_props_ad_.deadPoreVol(); const double dead_pore_vol = polymer_props_ad_.deadPoreVol();
const ADB cmax = computeCmax(state.concentration); const V cmax_v = computeCmax(x);
const ADB cmax = ADB::constant(cmax_v, state.concentration.blockPattern());
const ADB ads = polymer_props_ad_.adsorption(state.concentration, cmax); const ADB ads = polymer_props_ad_.adsorption(state.concentration, cmax);
const ADB krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration, cmax, kr[0], state.saturation[0]); const ADB krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration, cmax, kr[0], state.saturation[0]);
const ADB mc = computeMc(state); const ADB mc = computeMc(state);

View File

@ -191,8 +191,8 @@ namespace Opm {
computeFracFlow(const ADB& kro, computeFracFlow(const ADB& kro,
const ADB& krw_eff, const ADB& krw_eff,
const ADB& c) const; const ADB& c) const;
ADB V
computeCmax(const ADB& c) const; computeCmax(const PolymerBlackoilState& x) const;
ADB ADB
computeMc(const SolutionState& state) const; computeMc(const SolutionState& state) const;
ADB ADB

View File

@ -1,4 +1,3 @@
#include <opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp> #include <opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp>
#include <opm/core/pressure/tpfa/trans_tpfa.h> #include <opm/core/pressure/tpfa/trans_tpfa.h>
@ -335,6 +334,7 @@ namespace {
// Qs. // Qs.
state.qs = vars[ nextvar++ ]; state.qs = vars[ nextvar++ ];
// BHP. // BHP.
state.bhp = vars[ nextvar++ ]; state.bhp = vars[ nextvar++ ];
assert(nextvar == int(vars.size())); assert(nextvar == int(vars.size()));
@ -343,18 +343,18 @@ namespace {
} }
ADB V
FullyImplicitTwophasePolymerSolver:: FullyImplicitTwophasePolymerSolver::
computeCmax(const ADB& c) const computeCmax(const PolymerState& x) const
{ {
const int nc = c.value().size(); const int nc = c.value().size();
V cmax(nc); const V cmax = Eigen::Map<const V>(& x.maxconcentration()[0], nc);
const V c = Eigen::Map<const V>(& x.concentration()[0], nc);
for (int i = 0; i < nc; ++i) { for (int i = 0; i < nc; ++i) {
cmax(i) = (cmax(i) > c.value()(i)) ? cmax(i) : c.value()(i); cmax(i) = std::max(cmax(i), c(i));
} }
return ADB::constant(cmax, c.blockPattern()); return cmax;
} }
@ -374,7 +374,8 @@ namespace {
const V trans = subset(transmissibility(), ops_.internal_faces); const V trans = subset(transmissibility(), ops_.internal_faces);
const std::vector<ADB> kr = computeRelPerm(state); const std::vector<ADB> kr = computeRelPerm(state);
const ADB cmax = computeCmax(state.concentration); const V cmax_v = computeCmax(x);
const ADB cmax = ADB::constant(cmax, state.concentration.blockPattern());
const ADB ads = polymer_props_ad_.adsorption(state.concentration, cmax); const ADB ads = polymer_props_ad_.adsorption(state.concentration, cmax);
const ADB krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration, cmax, kr[0], state.saturation[0]); const ADB krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration, cmax, kr[0], state.saturation[0]);
const ADB mc = computeMc(state); const ADB mc = computeMc(state);
@ -594,7 +595,7 @@ namespace {
const V inSrc = Eigen::Map<const V>(& insrc[0], grid_.number_of_cells); const V inSrc = Eigen::Map<const V>(& insrc[0], grid_.number_of_cells);
const V polyin = Eigen::Map<const V>(& polymer_inflow_c[0], grid_.number_of_cells); const V polyin = Eigen::Map<const V>(& polymer_inflow_c[0], grid_.number_of_cells);
// compute the out-fracflow. // compute the out-fracflow.
const std::vector<ADB> f = computeFracFlow(kro, krw_eff, c); const std::vector<ADB> f = computeFracFlow();
// compute the in-fracflow. // compute the in-fracflow.
V zero = V::Zero(grid_.number_of_cells); V zero = V::Zero(grid_.number_of_cells);
V one = V::Ones(grid_.number_of_cells); V one = V::Ones(grid_.number_of_cells);
@ -614,9 +615,7 @@ namespace {
std::vector<ADB> std::vector<ADB>
FullyImplicitTwophasePolymerSolver::computeFracFlow(const ADB& kro, FullyImplicitTwophasePolymerSolver::computeFracFlow() const
const ADB& krw_eff,
const ADB& c) const
{ {
ADB total_mob = mob_[0] + mob_[1]; ADB total_mob = mob_[0] + mob_[1];

View File

@ -105,9 +105,7 @@ namespace Opm {
std::vector<ADB> std::vector<ADB>
computeFracFlow(const ADB& kro, computeFracFlow() const;
const ADB& krw_eff,
const ADB& c) const;
double double
residualNorm() const; residualNorm() const;
ADB ADB
@ -116,8 +114,8 @@ namespace Opm {
const std::vector<double>& polymer_inflow_c, const std::vector<double>& polymer_inflow_c,
const SolutionState& state) const; const SolutionState& state) const;
ADB V
computeCmax(const ADB& c) const; computeCmax(const PolymerState& x) const;
ADB ADB
computeMc(const SolutionState& state) const; computeMc(const SolutionState& state) const;
ADB ADB