mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
fix the polymer source bug,
warnning: water initial saturation should not bt zero when running this simulator.
This commit is contained in:
Liu Ming
parent
6fc24236df
commit
a358da7afa
Binary file not shown.
@ -24,12 +24,12 @@
|
||||
int main (int argc, char** argv)
|
||||
try
|
||||
{
|
||||
int nx = 20;
|
||||
int ny = 20;
|
||||
int nx = 30;
|
||||
int ny = 30;
|
||||
int nz = 1;
|
||||
double dx = 10.0;
|
||||
double dy = 10.0;
|
||||
double dz = 10.0;
|
||||
double dx = 2.0;
|
||||
double dy = 2.0;
|
||||
double dz = 0.5;
|
||||
using namespace Opm;
|
||||
parameter::ParameterGroup param(argc, argv, false);
|
||||
GridManager grid_manager(nx, ny, nz, dx, dy, dz);
|
||||
@ -40,23 +40,25 @@ try
|
||||
using namespace Opm::prefix;
|
||||
std::vector<double> density(num_phases, 1000.0);
|
||||
std::vector<double> viscosity(num_phases, 1.0*centi*Poise);
|
||||
double porosity = 0.5;
|
||||
viscosity[0] = 0.5 * centi * Poise;
|
||||
viscosity[1] = 5 * centi * Poise;
|
||||
double porosity = 0.35;
|
||||
double permeability = 10.0*milli*darcy;
|
||||
SaturationPropsBasic::RelPermFunc rel_perm_func = SaturationPropsBasic::Linear;
|
||||
IncompPropsAdBasic props(num_phases, rel_perm_func, density, viscosity,
|
||||
porosity, permeability, grid.dimensions, num_cells);
|
||||
std::vector<double> omega;
|
||||
std::vector<double> src(num_cells, 0.0);
|
||||
src[0] = 1.;
|
||||
src[num_cells-1] = -1.;
|
||||
src[0] = 10. / day;
|
||||
src[num_cells-1] = -10. / day;
|
||||
|
||||
FlowBCManager bcs;
|
||||
LinearSolverUmfpack linsolver;
|
||||
FullyImplicitTwoPhaseSolver solver(grid, props, linsolver);
|
||||
std::vector<double> porevol;
|
||||
Opm::computePorevolume(grid, props.porosity(), porevol);
|
||||
const double dt = param.getDefault("dt", 0.1) * day;
|
||||
const int num_time_steps = param.getDefault("nsteps", 20);
|
||||
const double dt = param.getDefault("dt", 10) * day;
|
||||
const int num_time_steps = param.getDefault("nsteps", 10);
|
||||
std::vector<int> allcells(num_cells);
|
||||
for (int cell = 0; cell < num_cells; ++cell) {
|
||||
allcells[cell] = cell;
|
||||
@ -66,10 +68,10 @@ try
|
||||
|
||||
//initial sat
|
||||
for (int c = 0; c < num_cells; ++c) {
|
||||
state.saturation()[2*c] = 0.2;
|
||||
state.saturation()[2*c+1] = 0.8;
|
||||
state.saturation()[2*c] = 0;
|
||||
state.saturation()[2*c+1] = 1;
|
||||
}
|
||||
std::vector<double> p(num_cells, 200*Opm::unit::barsa);
|
||||
std::vector<double> p(num_cells, 100*Opm::unit::barsa);
|
||||
state.pressure() = p;
|
||||
// state.setFirstSat(allcells, props, TwophaseState::MinSat);
|
||||
std::ostringstream vtkfilename;
|
||||
|
||||
@ -36,12 +36,12 @@ try
|
||||
}
|
||||
std::string deck_filename = param.get<std::string>("deck_filename");
|
||||
EclipseGridParser deck = EclipseGridParser(deck_filename);
|
||||
int nx = param.getDefault("nx", 20);
|
||||
int ny = param.getDefault("ny", 20);
|
||||
int nx = param.getDefault("nx", 30);
|
||||
int ny = param.getDefault("ny", 30);
|
||||
int nz = 1;
|
||||
double dx = 2.0;
|
||||
double dy = 2.0;
|
||||
double dz = 0.5;
|
||||
double dx = 10.0;
|
||||
double dy = 1.0;
|
||||
double dz = 1.0;
|
||||
GridManager grid_manager(nx, ny, nz, dx, dy, dz);
|
||||
const UnstructuredGrid& grid = *grid_manager.c_grid();
|
||||
int num_cells = grid.number_of_cells;
|
||||
@ -51,7 +51,7 @@ try
|
||||
std::vector<double> density(num_phases, 1000.0);
|
||||
std::vector<double> viscosity(num_phases, 1.0*centi*Poise);
|
||||
viscosity[0] = 0.5 * centi * Poise;
|
||||
viscosity[0] = 5 * centi * Poise;
|
||||
viscosity[1] = 5 * centi * Poise;
|
||||
double porosity = 0.35;
|
||||
double permeability = 10.0*milli*darcy;
|
||||
SaturationPropsBasic::RelPermFunc rel_perm_func = SaturationPropsBasic::Linear;
|
||||
@ -97,10 +97,10 @@ try
|
||||
std::vector<double> src(num_cells, 0.0);
|
||||
std::vector<double> src_polymer(num_cells);
|
||||
src[0] = 10. / day;
|
||||
src[num_cells-1] = -1. / day;
|
||||
src[num_cells-1] = -10. / day;
|
||||
|
||||
PolymerInflowBasic polymer_inflow(param.getDefault("poly_start_days", 30.0)*Opm::unit::day,
|
||||
param.getDefault("poly_end_days", 80.0)*Opm::unit::day,
|
||||
PolymerInflowBasic polymer_inflow(param.getDefault("poly_start_days", 300.0)*Opm::unit::day,
|
||||
param.getDefault("poly_end_days", 800.0)*Opm::unit::day,
|
||||
param.getDefault("poly_amount", polymer_props.cMax()));
|
||||
FlowBCManager bcs;
|
||||
LinearSolverUmfpack linsolver;
|
||||
@ -120,7 +120,7 @@ try
|
||||
state.saturation()[2*c] = 0.2;
|
||||
state.saturation()[2*c+1] = 0.8;
|
||||
}
|
||||
std::vector<double> p(num_cells, 200*Opm::unit::barsa);
|
||||
std::vector<double> p(num_cells, 100*Opm::unit::barsa);
|
||||
state.pressure() = p;
|
||||
|
||||
std::vector<double> c(num_cells, 0.0);
|
||||
@ -137,6 +137,7 @@ try
|
||||
Opm::DataMap dm;
|
||||
dm["saturation"] = &state.saturation();
|
||||
dm["pressure"] = &state.pressure();
|
||||
dm["concentration"] = &state.concentration();
|
||||
Opm::writeVtkData(grid, dm, vtkfile);
|
||||
}
|
||||
}
|
||||
|
||||
@ -97,13 +97,8 @@ typedef Eigen::Array<double,
|
||||
const SolutionState old_state = constantState(x);
|
||||
const double atol = 1.0e-12;
|
||||
const double rtol = 5.0e-8;
|
||||
const int maxit = 15;
|
||||
|
||||
std::cout << "Primary variables:\n";
|
||||
std::cout << "pressure\n"<<old_state.pressure<< std::endl;
|
||||
std::cout << "saturation\n"<<old_state.saturation[0] << std::endl;
|
||||
std::cout << "concentration\n"<<old_state.concentration << std::endl;
|
||||
assemble(pvdt, old_state, x, src, polymer_inflow);//, if_polymer_actived);
|
||||
const int maxit = 40;
|
||||
assemble(pvdt, old_state, x, src, polymer_inflow);
|
||||
|
||||
const double r0 = residualNorm();
|
||||
int it = 0;
|
||||
@ -115,7 +110,7 @@ typedef Eigen::Array<double,
|
||||
const V dx = solveJacobianSystem();
|
||||
updateState(dx, x);
|
||||
|
||||
assemble(pvdt, old_state, x, src, polymer_inflow);//, if_polymer_actived);
|
||||
assemble(pvdt, old_state, x, src, polymer_inflow);
|
||||
|
||||
const double r = residualNorm();
|
||||
|
||||
@ -238,10 +233,9 @@ typedef Eigen::Array<double,
|
||||
FullyImplicitTwophasePolymerSolver::
|
||||
assemble(const V& pvdt,
|
||||
const SolutionState& old_state,
|
||||
const PolymerState& x ,
|
||||
const PolymerState& x,
|
||||
const std::vector<double>& src,
|
||||
const std::vector<double>& polymer_inflow)
|
||||
// const bool if_polymer_actived)
|
||||
{
|
||||
// Create the primary variables.
|
||||
const SolutionState state = variableState(x);
|
||||
@ -256,44 +250,33 @@ typedef Eigen::Array<double,
|
||||
pvdt*(state.saturation[phase] - old_state.saturation[phase])
|
||||
+ ops_.div*mflux - source;
|
||||
}
|
||||
|
||||
// bool use_polymer = 1;//(amount != 0);
|
||||
// const int nc = grid_.number_of_cells;
|
||||
// if (use_polymer) {
|
||||
// Mass balance equation for polymer
|
||||
const ADB src_polymer = polymerSource(kr ,src, polymer_inflow, state);
|
||||
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));
|
||||
// }
|
||||
for (int i = 0; i < 3; ++i)
|
||||
std::cout<<"residual_["<<i<<"]\n"<<residual_[i] << std::endl;
|
||||
// Mass balance equation for polymer
|
||||
const ADB src_polymer = polymerSource(kr, src, polymer_inflow, state);
|
||||
// std::cout << "polymer src\n";
|
||||
// std::cout << src_polymer << std::endl;
|
||||
// const ADB src_polymer =ADB::function(V::Zero(grid_.number_of_cells), state.concentration.derivative());
|
||||
ADB mc = computeMc(state);
|
||||
|
||||
ADB poly_mflux = computePolymerMassFlux(trans, mc, kr, state);
|
||||
#if 0
|
||||
std::cout << "polymer mass\n";
|
||||
std::cout << pvdt * (state.saturation[0] * state.concentration
|
||||
- old_state.saturation[0] * old_state.concentration) <<std::endl;
|
||||
std::cout << "polymer mass flux\n";
|
||||
std::cout << poly_mflux << std::endl;
|
||||
#endif
|
||||
residual_[2] = pvdt * (state.saturation[0] * state.concentration
|
||||
- old_state.saturation[0] * old_state.concentration)
|
||||
+ ops_.div * poly_mflux - src_polymer;
|
||||
|
||||
}
|
||||
|
||||
|
||||
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
|
||||
FullyImplicitTwophasePolymerSolver::accumSource(const int phase,
|
||||
const std::vector<ADB>& kr,
|
||||
const std::vector<double>& src
|
||||
) const
|
||||
const std::vector<double>& src) const
|
||||
{
|
||||
//extract the source to out and in source.
|
||||
std::vector<double> outsrc;
|
||||
@ -329,8 +312,7 @@ typedef Eigen::Array<double,
|
||||
|
||||
ADB
|
||||
FullyImplicitTwophasePolymerSolver::
|
||||
polymerSource(
|
||||
const std::vector<ADB>& kr,
|
||||
polymerSource(const std::vector<ADB>& kr,
|
||||
const std::vector<double>& src,
|
||||
const std::vector<double>& polymer_inflow_c,
|
||||
const SolutionState& state) const
|
||||
@ -359,13 +341,15 @@ typedef Eigen::Array<double,
|
||||
ADB f_out = computeFracFlow(0, kr);
|
||||
// compute the in-fracflow.
|
||||
V f_in = V::Ones(grid_.number_of_cells);
|
||||
return f_out * outSrc * state.concentration + f_in * inSrc * polyin ;
|
||||
|
||||
// ADB polymer_insrc = ADB::function(f_in * inSrc * polyin, state.concentration.derivative());
|
||||
return f_out * outSrc * state.concentration + f_in * inSrc * polyin;
|
||||
}
|
||||
|
||||
|
||||
ADB
|
||||
FullyImplicitTwophasePolymerSolver::computeFracFlow(int phase,
|
||||
const std::vector<ADB>& kr) const
|
||||
FullyImplicitTwophasePolymerSolver::computeFracFlow(int phase,
|
||||
const std::vector<ADB>& kr) const
|
||||
{
|
||||
const double* mus = fluid_.viscosity();
|
||||
ADB mob_phase = kr[phase] / V::Constant(kr[phase].size(), 1, mus[phase]);
|
||||
@ -498,8 +482,25 @@ typedef Eigen::Array<double,
|
||||
return upwind.select(mob) * head;
|
||||
}
|
||||
|
||||
ADB
|
||||
FullyImplicitTwophasePolymerSolver::computePolymerMassFlux(const V& trans,
|
||||
const ADB& mc,
|
||||
const std::vector<ADB>& kr,
|
||||
const SolutionState& state) const
|
||||
|
||||
{
|
||||
const double* mus = fluid_.viscosity();
|
||||
ADB water_mob = kr[0] / V::Constant(kr[0].size(), 1, mus[0]);
|
||||
ADB poly_mob = state.concentration * mc * water_mob;
|
||||
|
||||
const ADB dp = ops_.ngrad * state.pressure;
|
||||
const ADB head = trans * dp;
|
||||
|
||||
UpwindSelector<double> upwind(grid_, ops_, head.value());
|
||||
|
||||
return upwind.select(poly_mob) * head;
|
||||
|
||||
}
|
||||
|
||||
|
||||
double
|
||||
|
||||
@ -30,7 +30,6 @@ namespace Opm {
|
||||
const std::vector<double>& src,
|
||||
const std::vector<double>& polymer_inflow
|
||||
);
|
||||
// const bool if_polymer_actived);
|
||||
private:
|
||||
typedef AutoDiffBlock<double> ADB;
|
||||
typedef ADB::V V;
|
||||
@ -44,7 +43,6 @@ namespace Opm {
|
||||
ADB pressure;
|
||||
std::vector<ADB> saturation;
|
||||
ADB concentration;
|
||||
// ADB cmax;
|
||||
};
|
||||
const UnstructuredGrid& grid_;
|
||||
const IncompPropsAdInterface& fluid_;
|
||||
@ -65,7 +63,6 @@ namespace Opm {
|
||||
const PolymerState& x,
|
||||
const std::vector<double>& src,
|
||||
const std::vector<double>& polymer_inflow);
|
||||
// const bool if_polymer_actived);
|
||||
V solveJacobianSystem() const;
|
||||
void updateState(const V& dx,
|
||||
PolymerState& x) const;
|
||||
@ -85,6 +82,11 @@ namespace Opm {
|
||||
const V& trans,
|
||||
const std::vector<ADB>& kr,
|
||||
const SolutionState& state) const;
|
||||
ADB
|
||||
computePolymerMassFlux(const V& trans,
|
||||
const ADB& mc,
|
||||
const std::vector<ADB>& kr,
|
||||
const SolutionState& state) const;
|
||||
double
|
||||
residualNorm() const;
|
||||
ADB
|
||||
@ -103,8 +105,8 @@ namespace Opm {
|
||||
fluidDensity(const int phase) const;
|
||||
ADB
|
||||
transMult(const ADB& p) const;
|
||||
double
|
||||
polymerInjectedAmount(const std::vector<double>& polymer_inflow) const;
|
||||
};
|
||||
|
||||
|
||||
} // namespace Opm
|
||||
#endif// OPM_FULLYIMPLICITTWOPHASESOLVER_HEADER_INCLUDED
|
||||
|
||||
Loading…
Reference in New Issue
Block a user