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

remove some output commits.

Browse Source 
get values from command line for exapmes.
This commit is contained in:
Liu Ming 2013-12-11 18:27:05 +08:00
parent a358da7afa
commit 11ea882ed8
4 changed files with 54 additions and 36 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

33
examples/sim_2p_fincomp_ad.cpp
View File

@ -24,14 +24,14 @@
int main (int argc, char** argv)
try
{
int nx = 30;
int ny = 30;
int nz = 1;
double dx = 2.0;
double dy = 2.0;
double dz = 0.5;
using namespace Opm;
parameter::ParameterGroup param(argc, argv, false);
int nx = param.getDefault("nx", 30);
int ny = param.getDefault("ny", 1);
int nz = 1;
double dx = 10./nx;
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;
@ -49,15 +49,15 @@ try
porosity, permeability, grid.dimensions, num_cells);
std::vector<double> omega;
std::vector<double> src(num_cells, 0.0);
src[0] = 10. / day;
src[num_cells-1] = -10. / day;
src[0] = 1. / day;
src[num_cells-1] = -1. / 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", 10) * day;
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) {
@ -68,18 +68,25 @@ try
//initial sat
for (int c = 0; c < num_cells; ++c) {
state.saturation()[2*c] = 0;
state.saturation()[2*c+1] = 1;
state.saturation()[2*c] = 0.2;
state.saturation()[2*c+1] = 0.8;
}
std::vector<double> p(num_cells, 100*Opm::unit::barsa);
state.pressure() = p;
// state.setFirstSat(allcells, props, TwophaseState::MinSat);
std::ostringstream vtkfilename;
// Write the initial state.
vtkfilename.str("");
vtkfilename << "sim_2p_fincomp_" << std::setw(3) << std::setfill('0') << 0 << ".vtu";
std::ofstream vtkfile(vtkfilename.str().c_str());
Opm::DataMap dm;
dm["saturation"] = &state.saturation();
dm["pressure"] = &state.pressure();
Opm::writeVtkData(grid, dm, vtkfile);
for (int i = 0; i < num_time_steps; ++i) {
solver.step(dt, state, src);
vtkfilename.str("");
vtkfilename << "sim_2p_fincomp-" << std::setw(3) << std::setfill('0') << i << ".vtu";
vtkfilename << "sim_2p_fincomp_" << std::setw(3) << std::setfill('0') << i + 1 << ".vtu";
std::ofstream vtkfile(vtkfilename.str().c_str());
Opm::DataMap dm;
dm["saturation"] = &state.saturation();

20
examples/sim_poly2p_fincomp_ad.cpp
View File

@ -39,7 +39,7 @@ try
int nx = param.getDefault("nx", 30);
int ny = param.getDefault("ny", 30);
int nz = 1;
double dx = 10.0;
double dx = 10./nx;
double dy = 1.0;
double dz = 1.0;
GridManager grid_manager(nx, ny, nz, dx, dy, dz);
@ -60,7 +60,7 @@ try
// Init polymer properties.
// Setting defaults to provide a simple example case.
PolymerProperties polymer_props(deck);
PolymerProperties polymer_props(deck);
#if 0
if (use_poly_deck) {
} else {
@ -96,8 +96,8 @@ try
std::vector<double> omega;
std::vector<double> src(num_cells, 0.0);
std::vector<double> src_polymer(num_cells);
src[0] = 10. / day;
src[num_cells-1] = -10. / day;
src[0] = param.getDefault("insrc", 1.) / day;
src[num_cells-1] = -param.getDefault("insrc", 1.) / day;
PolymerInflowBasic polymer_inflow(param.getDefault("poly_start_days", 300.0)*Opm::unit::day,
param.getDefault("poly_end_days", 800.0)*Opm::unit::day,
@ -127,12 +127,22 @@ try
state.concentration() = c;
std::ostringstream vtkfilename;
double currentime = 0;
// Write the initial state.
vtkfilename.str("");
vtkfilename << "sim_poly2p_fincomp_ad_" << std::setw(3) << std::setfill('0') << 0<< ".vtu";
std::ofstream vtkfile(vtkfilename.str().c_str());
Opm::DataMap dm;
dm["saturation"] = &state.saturation();
dm["pressure"] = &state.pressure();
dm["concentration"] = &state.concentration();
Opm::writeVtkData(grid, dm, vtkfile);
for (int i = 0; i < num_time_steps; ++i) {
currentime += dt;
polymer_inflow.getInflowValues(currentime, currentime+dt, src_polymer);
solver.step(dt, state, src, src_polymer);
vtkfilename.str("");
vtkfilename << "sim_poly2p_fincomp_ad_" << std::setw(3) << std::setfill('0') << i << ".vtu";
vtkfilename << "sim_poly2p_fincomp_ad_" << std::setw(3) << std::setfill('0') << i + 1<< ".vtu";
std::ofstream vtkfile(vtkfilename.str().c_str());
Opm::DataMap dm;
dm["saturation"] = &state.saturation();

3
opm/polymer/fullyimplicit/FullyImplicitTwoPhaseSolver.cpp
View File

@ -396,7 +396,8 @@ typedef Eigen::Array<double,
// const ADB tr_mult = transMult(state.pressure);
const double* mus = fluid_.viscosity();
ADB mob = kr[phase] / V::Constant(kr[phase].size(), 1, mus[phase]);
if (phase ==0)
std::cout << "watetr mob\n" << mob.value() << std::endl;
const ADB dp = ops_.ngrad * state.pressure;
const ADB head = trans * dp;

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

@ -165,8 +165,8 @@ typedef Eigen::Array<double,
// Concentration
assert(not x.concentration().empty());
const V c = Eigen::Map<const V>(&x.concentration()[0], nc);
state.concentration = ADB::constant(c);
return state;
}
@ -246,25 +246,16 @@ typedef Eigen::Array<double,
for (int phase = 0; phase < fluid_.numPhases(); ++phase) {
const ADB mflux = computeMassFlux(phase, trans, kr, state);
ADB source = accumSource(phase, kr, src);
// std::cout << "phase-"<<phase<<"-source\n"<< source.value() << std::endl;
residual_[phase] =
pvdt*(state.saturation[phase] - old_state.saturation[phase])
+ ops_.div*mflux - source;
}
// 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
// std::cout << "polymer source \n" << src_polymer.value() << std::endl;
residual_[2] = pvdt * (state.saturation[0] * state.concentration
- old_state.saturation[0] * old_state.concentration)
+ ops_.div * poly_mflux - src_polymer;
@ -418,7 +409,6 @@ typedef Eigen::Array<double,
const V p = p_old - dp;
std::copy(&p[0], &p[0] + nc, state.pressure().begin());
// Saturation updates.
const double dsmax = 0.3;
const DataBlock s_old = Eigen::Map<const DataBlock>(& state.saturation()[0], nc, np);
@ -472,8 +462,18 @@ typedef Eigen::Array<double,
{
// const ADB tr_mult = transMult(state.pressure);
const double* mus = fluid_.viscosity();
ADB mob = kr[phase] / V::Constant(kr[phase].size(), 1, mus[phase]);
ADB mob = ADB::null();
if (phase == 0) {
ADB inv_wat_eff_vis = polymer_props_ad_.effectiveInvWaterVisc(state.concentration, mus);
// for (int i = 0; i < grid_.number_of_cells; ++i) {
// std::cout << state.concentration.value()(i) << " " << 1./inv_wat_eff_vis.value()(i) << std::endl;
// std::cout << "water effective vis\n"<<1./inv_wat_eff_v1./inv_wat_eff_vis.value()is.value()<<std::endl;
// }
mob = kr[0] * inv_wat_eff_vis;
// std::cout << "watetr mob\n" << mob.value() << std::endl;
} else if (phase == 1) {
mob = kr[phase] / V::Constant(kr[phase].size(), 1, mus[phase]);
}
const ADB dp = ops_.ngrad * state.pressure;
const ADB head = trans * dp;
@ -490,8 +490,8 @@ typedef Eigen::Array<double,
{
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;
ADB inv_wat_eff_vis = polymer_props_ad_.effectiveInvWaterVisc(state.concentration, mus);
ADB poly_mob = mc * kr[0] * inv_wat_eff_vis;
const ADB dp = ops_.ngrad * state.pressure;
const ADB head = trans * dp;