OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

normalized tabification in spu_2p.

Browse Source
This commit is contained in:
Xavier Raynaud 2012-05-10 13:50:22 +02:00
parent 51c1da1df8
commit fd9a3318b0
1 changed files with 85 additions and 84 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

169
examples/spu_2p.cpp
View File

@ -162,9 +162,9 @@ class SimpleFluid2pWrappingProps
{ {
public: public:
SimpleFluid2pWrappingProps(const Opm::IncompPropertiesInterface& props) SimpleFluid2pWrappingProps(const Opm::IncompPropertiesInterface& props)
: props_(props), : props_(props),
smin_(props.numCells()*props.numPhases()), smin_(props.numCells()*props.numPhases()),
smax_(props.numCells()*props.numPhases()) smax_(props.numCells()*props.numPhases())
{ {
if (props.numPhases() != 2) { if (props.numPhases() != 2) {
THROW("SimpleFluid2pWrapper requires 2 phases."); THROW("SimpleFluid2pWrapper requires 2 phases.");
@ -183,8 +183,8 @@ public:
} }
template <class Sat, template <class Sat,
class Mob, class Mob,
class DMob> class DMob>
void mobility(int c, const Sat& s, Mob& mob, DMob& dmob) const void mobility(int c, const Sat& s, Mob& mob, DMob& dmob) const
{ {
props_.relperm(1, &s[0], &c, &mob[0], &dmob[0]); props_.relperm(1, &s[0], &c, &mob[0], &dmob[0]);
@ -203,8 +203,8 @@ public:
} }
template <class Sat, template <class Sat,
class Pcap, class Pcap,
class DPcap> class DPcap>
void pc(int c, const Sat& s, Pcap& pcap, DPcap& dpcap) const void pc(int c, const Sat& s, Pcap& pcap, DPcap& dpcap) const
{ {
double pcow[2]; double pcow[2];
@ -252,12 +252,12 @@ public:
}; };
typedef Opm::ImplicitTransport<TransportModel, typedef Opm::ImplicitTransport<TransportModel,
JacSys , JacSys ,
MaxNorm , MaxNorm ,
VectorNegater , VectorNegater ,
VectorZero , VectorZero ,
MatrixZero , MatrixZero ,
VectorAssign > TransportSolver; VectorAssign > TransportSolver;
@ -324,11 +324,11 @@ main(int argc, char** argv)
// Gravity. // Gravity.
gravity[2] = deck.hasField("NOGRAV") ? 0.0 : Opm::unit::gravity; gravity[2] = deck.hasField("NOGRAV") ? 0.0 : Opm::unit::gravity;
// Init state variables (saturation and pressure). // Init state variables (saturation and pressure).
if (param.has("init_saturation")) { if (param.has("init_saturation")) {
initStateTwophaseBasic(*grid->c_grid(), *props, param, gravity[2], state); initStateTwophaseBasic(*grid->c_grid(), *props, param, gravity[2], state);
} else { } else {
initStateTwophaseFromDeck(*grid->c_grid(), *props, deck, gravity[2], state); initStateTwophaseFromDeck(*grid->c_grid(), *props, deck, gravity[2], state);
} }
} else { } else {
// Grid init. // Grid init.
const int nx = param.getDefault("nx", 100); const int nx = param.getDefault("nx", 100);
@ -560,75 +560,76 @@ main(int argc, char** argv)
if (rock_comp->isActive()) { if (rock_comp->isActive()) {
rc.resize(num_cells); rc.resize(num_cells);
std::vector<double> initial_pressure = state.pressure(); std::vector<double> initial_pressure = state.pressure();
std::vector<double> initial_porevolume(num_cells); std::vector<double> initial_porevolume(num_cells);
computePorevolume(*grid->c_grid(), *props, *rock_comp, initial_pressure, initial_porevolume); computePorevolume(*grid->c_grid(), *props, *rock_comp, initial_pressure, initial_porevolume);
std::vector<double> pressure_increment(num_cells + num_wells); std::vector<double> pressure_increment(num_cells + num_wells);
std::vector<double> prev_pressure(num_cells + num_wells); std::vector<double> prev_pressure(num_cells + num_wells);
for (int iter = 0; iter < nl_pressure_maxiter; ++iter) { for (int iter = 0; iter < nl_pressure_maxiter; ++iter) {
for (int cell = 0; cell < num_cells; ++cell) {
rc[cell] = rock_comp->rockComp(state.pressure()[cell]);
}
computePorevolume(*grid->c_grid(), *props, *rock_comp, state.pressure(), porevol);
std::copy(state.pressure().begin(), state.pressure().end(), prev_pressure.begin());
std::copy(well_bhp.begin(), well_bhp.end(), prev_pressure.begin() + num_cells);
// prev_pressure = state.pressure();
// compute pressure increment
psolver.solveIncrement(totmob, omega, src, wdp, bcs.c_bcs(), porevol, rc,
prev_pressure, initial_porevolume, simtimer.currentStepLength(),
pressure_increment);
double max_change = 0.0; for (int cell = 0; cell < num_cells; ++cell) {
for (int cell = 0; cell < num_cells; ++cell) { rc[cell] = rock_comp->rockComp(state.pressure()[cell]);
state.pressure()[cell] += pressure_increment[cell]; }
max_change = std::max(max_change, std::fabs(pressure_increment[cell])); computePorevolume(*grid->c_grid(), *props, *rock_comp, state.pressure(), porevol);
} std::copy(state.pressure().begin(), state.pressure().end(), prev_pressure.begin());
for (int well = 0; well < num_wells; ++well) { std::copy(well_bhp.begin(), well_bhp.end(), prev_pressure.begin() + num_cells);
well_bhp[well] += pressure_increment[num_cells + well]; // prev_pressure = state.pressure();
max_change = std::max(max_change, std::fabs(pressure_increment[num_cells + well]));
}
std::cout << "Pressure iter " << iter << " max change = " << max_change << std::endl;
if (max_change < nl_pressure_tolerance) {
break;
}
}
psolver.computeFaceFlux(totmob, omega, src, wdp, bcs.c_bcs(), state.pressure(), state.faceflux(),
well_bhp, well_perfrates);
} else {
psolver.solve(totmob, omega, src, wdp, bcs.c_bcs(), state.pressure(), state.faceflux(),
well_bhp, well_perfrates);
}
pressure_timer.stop();
double pt = pressure_timer.secsSinceStart();
std::cout << "Pressure solver took: " << pt << " seconds." << std::endl;
ptime += pt;
if (check_well_controls) { // compute pressure increment
Opm::computePhaseFlowRatesPerWell(*wells->c_wells(), psolver.solveIncrement(totmob, omega, src, wdp, bcs.c_bcs(), porevol, rc,
fractional_flows, prev_pressure, initial_porevolume, simtimer.currentStepLength(),
well_perfrates, pressure_increment);
well_resflows_phase);
std::cout << "Checking well conditions." << std::endl; double max_change = 0.0;
// For testing we set surface := reservoir for (int cell = 0; cell < num_cells; ++cell) {
well_control_passed = wells->conditionsMet(well_bhp, well_resflows_phase, well_resflows_phase); state.pressure()[cell] += pressure_increment[cell];
++well_control_iteration; max_change = std::max(max_change, std::fabs(pressure_increment[cell]));
if (!well_control_passed && well_control_iteration > max_well_control_iterations) { }
THROW("Could not satisfy well conditions in " << max_well_control_iterations << " tries."); for (int well = 0; well < num_wells; ++well) {
} well_bhp[well] += pressure_increment[num_cells + well];
if (!well_control_passed) { max_change = std::max(max_change, std::fabs(pressure_increment[num_cells + well]));
std::cout << "Well controls not passed, solving again." << std::endl; }
} else {
std::cout << "Well conditions met." << std::endl; std::cout << "Pressure iter " << iter << " max change = " << max_change << std::endl;
} if (max_change < nl_pressure_tolerance) {
} break;
} while (!well_control_passed); }
}
// Process transport sources (to include bdy terms and well flows). psolver.computeFaceFlux(totmob, omega, src, wdp, bcs.c_bcs(), state.pressure(), state.faceflux(),
Opm::computeTransportSource(*grid->c_grid(), src, state.faceflux(), 1.0, well_bhp, well_perfrates);
wells->c_wells(), well_perfrates, reorder_src); } else {
if (!use_reorder) { psolver.solve(totmob, omega, src, wdp, bcs.c_bcs(), state.pressure(), state.faceflux(),
well_bhp, well_perfrates);
}
pressure_timer.stop();
double pt = pressure_timer.secsSinceStart();
std::cout << "Pressure solver took: " << pt << " seconds." << std::endl;
ptime += pt;
if (check_well_controls) {
Opm::computePhaseFlowRatesPerWell(*wells->c_wells(),
fractional_flows,
well_perfrates,
well_resflows_phase);
std::cout << "Checking well conditions." << std::endl;
// For testing we set surface := reservoir
well_control_passed = wells->conditionsMet(well_bhp, well_resflows_phase, well_resflows_phase);
++well_control_iteration;
if (!well_control_passed && well_control_iteration > max_well_control_iterations) {
THROW("Could not satisfy well conditions in " << max_well_control_iterations << " tries.");
}
if (!well_control_passed) {
std::cout << "Well controls not passed, solving again." << std::endl;
} else {
std::cout << "Well conditions met." << std::endl;
}
}
} while (!well_control_passed);
// Process transport sources (to include bdy terms and well flows).
Opm::computeTransportSource(*grid->c_grid(), src, state.faceflux(), 1.0,
wells->c_wells(), well_perfrates, reorder_src);
if (!use_reorder) {
clear_transport_source(tsrc); clear_transport_source(tsrc);
for (int cell = 0; cell < num_cells; ++cell) { for (int cell = 0; cell < num_cells; ++cell) {
if (reorder_src[cell] > 0.0) { if (reorder_src[cell] > 0.0) {