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

Merge remote-tracking branch 'bska/fully-implicit' into fully-implicit

Browse Source
This commit is contained in:
Atgeirr Flø Rasmussen
2013-05-24 21:46:11 +02:00
parent dc83ecb18e 26e965e858
commit a3be33a312
4 changed files with 116 additions and 26 deletions
Download Patch File Download Diff File Expand all files Collapse all files

64
examples/test_implicit_ad.cpp
View File

@@ -25,6 +25,8 @@
#include <opm/autodiff/BlackoilPropsAd.hpp>
#include <opm/core/grid.h>
#include <opm/core/wells.h>
#include <opm/core/grid/GridManager.hpp>
#include <opm/core/linalg/LinearSolverFactory.hpp>
@@ -35,10 +37,52 @@
#include <opm/core/utility/Units.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/simulator/initState.hpp>
#include <algorithm>
#include <boost/shared_ptr.hpp>
namespace {
boost::shared_ptr<Wells>
createWellConfig()
{
boost::shared_ptr<Wells> wells(create_wells(2, 2, 2),
destroy_wells);
const double inj_frac[] = { 1.0, 0.0 };
const double prod_frac[] = { 0.0, 0.0 };
const int num_inj = 1;
const int inj_cells[num_inj] = { 0 };
const int num_prod = 1;
const int prod_cells[num_prod] = { 19 };
const double WI[3] = { 1e-12, 1e-12, 1e-12 };
bool ok = add_well(INJECTOR, 0.0, num_inj, inj_frac, inj_cells, WI, "Inj", wells.get());
ok = ok && add_well(PRODUCER, 0.0, num_prod, prod_frac, prod_cells, WI, "Prod", wells.get());
ok = ok && append_well_controls(BHP, 500.0*Opm::unit::barsa, 0, 0, wells.get());
// ok = ok && append_well_controls(BHP, 200.0*Opm::unit::barsa, 0, 1, wells);
double oildistr[2] = { 0.0, 1.0 };
ok = ok && append_well_controls(SURFACE_RATE, 1e-3, oildistr, 1, wells.get());
if (!ok) {
THROW("Something went wrong with well init.");
}
set_current_control(0, 0, wells.get());
set_current_control(1, 0, wells.get());
return wells;
}
template <class Ostream, typename T, class A>
Ostream&
operator<<(Ostream& os, const std::vector<T,A>& v)
{
std::copy(v.begin(), v.end(), std::ostream_iterator<T>(os, " "));
return os;
}
}
int
main(int argc, char* argv[])
{
@@ -50,30 +94,26 @@ main(int argc, char* argv[])
const Opm::BlackoilPropertiesBasic props0(param, 2, nc);
const Opm::BlackoilPropsAd props(props0);
boost::shared_ptr<Wells> wells = createWellConfig();
typedef Opm::FullyImplicitBlackoilSolver BOSolver;
double grav[] = { 1.0, 0.0 };
double grav[] = { 0.0, 0.0 };
Opm::DerivedGeology geo(*g, props, grav);
BOSolver solver(*g, props, geo);
BOSolver solver(*g, props, geo, *wells);
Opm::BlackoilState state;
initStateBasic(*g, props0, param, 0.0, state);
initBlackoilSurfvol(*g, props0, state);
solver.step(1.0, state);
#if 0
Opm::WellState well_state;
well_state.init(wells, state);
well_state.init(wells.get(), state);
ps.solve(1.0, state, well_state);
solver.step(1.0, state, well_state);
std::copy(state.pressure().begin(), state.pressure().end(), std::ostream_iterator<double>(std::cout, " "));
std::cout << std::endl;
std::copy(well_state.bhp().begin(), well_state.bhp().end(), std::ostream_iterator<double>(std::cout, " "));
std::cout << std::endl;
#endif
std::cout << state.pressure() << '\n'
<< well_state.bhp() << '\n';
return 0;
}

54
opm/autodiff/FullyImplicitBlackoilSolver.cpp
View File

@@ -26,6 +26,7 @@
#include <opm/autodiff/GeoProps.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/grid.h>
#include <opm/core/utility/ErrorMacros.hpp>
@@ -133,16 +134,20 @@ namespace {
namespace Opm {
FullyImplicitBlackoilSolver::FullyImplicitBlackoilSolver(const UnstructuredGrid& grid ,
const BlackoilPropsAdInterface& fluid,
const DerivedGeology& geo )
FullyImplicitBlackoilSolver::
FullyImplicitBlackoilSolver(const UnstructuredGrid& grid ,
const BlackoilPropsAdInterface& fluid,
const DerivedGeology& geo ,
const Wells& wells)
: grid_ (grid)
, fluid_ (fluid)
, geo_ (geo)
, wells_ (wells)
, active_(activePhases(fluid.phaseUsage()))
, canph_ (active2Canonical(fluid.phaseUsage()))
, cells_ (buildAllCells(grid.number_of_cells))
, ops_ (grid)
, wops_ (wells)
, grav_ (gravityOperator(grid_, ops_, geo_))
, rq_ (fluid.numPhases())
{
@@ -150,8 +155,10 @@ namespace Opm {
}
void
FullyImplicitBlackoilSolver::step(const double dt,
BlackoilState& x)
FullyImplicitBlackoilSolver::
step(const double dt,
BlackoilState& x ,
WellState& xw)
{
const V dtpv = geo_.poreVolume() / dt;
@@ -166,7 +173,7 @@ namespace Opm {
const int maxit = 15;
#endif
assemble(dtpv, x);
assemble(dtpv, x, xw);
#if 0
const double r0 = residualNorm();
@@ -175,7 +182,7 @@ namespace Opm {
while (resTooLarge && (it < maxit)) {
solveJacobianSystem(x);
assemble(dtpv, x);
assemble(dtpv, x, xw);
const double r = residualNorm();
@@ -209,6 +216,34 @@ namespace Opm {
}
FullyImplicitBlackoilSolver::
WellOps::WellOps(const Wells& wells)
: w2p(wells.well_connpos[ wells.number_of_wells ],
wells.number_of_wells)
, p2w(wells.number_of_wells,
wells.well_connpos[ wells.number_of_wells ])
{
const int nw = wells.number_of_wells;
const int* const wpos = wells.well_connpos;
typedef Eigen::Triplet<double> Tri;
std::vector<Tri> scatter, gather;
scatter.reserve(wpos[nw]);
gather .reserve(wpos[nw]);
for (int w = 0, i = 0; w < nw; ++w) {
for (; i < wpos[ w + 1 ]; ++i) {
scatter.push_back(Tri(i, w, 1.0));
gather .push_back(Tri(w, i, 1.0));
}
}
w2p.setFromTriplets(scatter.begin(), scatter.end());
p2w.setFromTriplets(gather .begin(), gather .end());
}
void
FullyImplicitBlackoilSolver::allocateResidual()
{
@@ -354,7 +389,10 @@ namespace Opm {
}
void
FullyImplicitBlackoilSolver::assemble(const V& dtpv, const BlackoilState& x)
FullyImplicitBlackoilSolver::
assemble(const V& dtpv,
const BlackoilState& x ,
const WellState& xw )
{
const V transi = subset(geo_.transmissibility(),
ops_.internal_faces);

20
opm/autodiff/FullyImplicitBlackoilSolver.hpp
View File

@@ -41,7 +41,8 @@ namespace Opm {
public:
FullyImplicitBlackoilSolver(const UnstructuredGrid& grid ,
const BlackoilPropsAdInterface& fluid,
const DerivedGeology& geo );
const DerivedGeology& geo ,
const Wells& wells);
/// Take a single forward step, modifiying
/// state.pressure()
@@ -49,8 +50,9 @@ namespace Opm {
/// state.saturation()
/// state.surfacevol()
void
step(const double dt,
BlackoilState& state);
step(const double dt ,
BlackoilState& state ,
WellState& wstate);
private:
// Types and enums
@@ -78,6 +80,12 @@ namespace Opm {
ADB Rs;
};
struct WellOps {
WellOps(const Wells& wells);
M w2p; // well -> perf (scatter)
M p2w; // perf -> well (gather)
};
enum { Water = BlackoilPropsAdInterface::Water,
Oil = BlackoilPropsAdInterface::Oil ,
Gas = BlackoilPropsAdInterface::Gas };
@@ -86,12 +94,14 @@ namespace Opm {
const UnstructuredGrid& grid_;
const BlackoilPropsAdInterface& fluid_;
const DerivedGeology& geo_;
const Wells& wells_;
// For each canonical phase -> true if active
const std::vector<bool> active_;
// Size = # active faces. Maps active -> canonical phase indices.
const std::vector<int> canph_;
const std::vector<int> cells_; // All grid cells
HelperOps ops_;
const WellOps wops_;
const M grav_;
std::vector<ReservoirResidualQuant> rq_;
@@ -115,7 +125,9 @@ namespace Opm {
const int aix );
void
assemble(const V& dtpv, const BlackoilState& x);
assemble(const V& dtpv,
const BlackoilState& x ,
const WellState& xw );
std::vector<ADB>
computeRelPerm(const SolutionState& state);

4
opm/autodiff/SimulatorFullyImplicitBlackoil.cpp
View File

@@ -252,7 +252,7 @@ namespace Opm
gravity_(gravity),
fluid_(props_),
geo_(grid_, fluid_, gravity_),
solver_(grid_, fluid_, geo_/*, *wells_manager.c_wells(), linsolver*/)
solver_(grid_, fluid_, geo_, *wells_manager.c_wells()/*, linsolver*/)
/* param.getDefault("nl_pressure_residual_tolerance", 0.0),
param.getDefault("nl_pressure_change_tolerance", 1.0),
param.getDefault("nl_pressure_maxiter", 10),
@@ -361,7 +361,7 @@ namespace Opm
// Run solver.
solver_timer.start();
std::vector<double> initial_pressure = state.pressure();
solver_.step(timer.currentStepLength(), state/*, well_state*/);
solver_.step(timer.currentStepLength(), state, well_state);
// Stop timer and report.
solver_timer.stop();