From c5ad446a0a15b49e528b119e3408968772aef8c7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Fri, 24 May 2013 15:27:19 +0200 Subject: [PATCH 1/3] Include well support in FIBOSolver interface. Update callers accordingly. --- examples/test_implicit_ad.cpp | 46 +++++++++++++++++-- opm/autodiff/FullyImplicitBlackoilSolver.cpp | 25 ++++++---- opm/autodiff/FullyImplicitBlackoilSolver.hpp | 13 ++++-- .../SimulatorFullyImplicitBlackoil.cpp | 4 +- 4 files changed, 71 insertions(+), 17 deletions(-) diff --git a/examples/test_implicit_ad.cpp b/examples/test_implicit_ad.cpp index a36632429..09e0ca474 100644 --- a/examples/test_implicit_ad.cpp +++ b/examples/test_implicit_ad.cpp @@ -25,6 +25,8 @@ #include #include +#include + #include #include @@ -35,10 +37,43 @@ #include #include +#include #include #include +#include + +namespace { + boost::shared_ptr + createWellConfig() + { + boost::shared_ptr 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; + } +} + int main(int argc, char* argv[]) { @@ -50,18 +85,23 @@ main(int argc, char* argv[]) const Opm::BlackoilPropertiesBasic props0(param, 2, nc); const Opm::BlackoilPropsAd props(props0); + boost::shared_ptr 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); + Opm::WellState well_state; + well_state.init(wells.get(), state); + + solver.step(1.0, state, well_state); #if 0 Opm::WellState well_state; diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.cpp b/opm/autodiff/FullyImplicitBlackoilSolver.cpp index b24a46964..55a40a8fa 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.cpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.cpp @@ -26,6 +26,7 @@ #include #include +#include #include #include @@ -133,12 +134,15 @@ 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)) @@ -150,8 +154,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 +172,7 @@ namespace Opm { const int maxit = 15; #endif - assemble(dtpv, x); + assemble(dtpv, x, xw); #if 0 const double r0 = residualNorm(); @@ -175,7 +181,7 @@ namespace Opm { while (resTooLarge && (it < maxit)) { solveJacobianSystem(x); - assemble(dtpv, x); + assemble(dtpv, x, xw); const double r = residualNorm(); @@ -354,7 +360,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); diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.hpp b/opm/autodiff/FullyImplicitBlackoilSolver.hpp index ea737df7f..e1f539f0a 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.hpp @@ -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 @@ -86,6 +88,7 @@ namespace Opm { const UnstructuredGrid& grid_; const BlackoilPropsAdInterface& fluid_; const DerivedGeology& geo_; + const Wells& wells_; // For each canonical phase -> true if active const std::vector active_; // Size = # active faces. Maps active -> canonical phase indices. @@ -115,7 +118,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 computeRelPerm(const SolutionState& state); diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoil.cpp b/opm/autodiff/SimulatorFullyImplicitBlackoil.cpp index 01db9f573..e70ef665e 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoil.cpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoil.cpp @@ -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 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(); From 2dffbb3a70577a006fdec1c13f5c3ce0ec14e0d6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Fri, 24 May 2013 17:22:35 +0200 Subject: [PATCH 2/3] Add gather/scatter support for wells --- opm/autodiff/FullyImplicitBlackoilSolver.cpp | 29 ++++++++++++++++++++ opm/autodiff/FullyImplicitBlackoilSolver.hpp | 7 +++++ 2 files changed, 36 insertions(+) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.cpp b/opm/autodiff/FullyImplicitBlackoilSolver.cpp index 55a40a8fa..f9a477b2b 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.cpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.cpp @@ -147,6 +147,7 @@ namespace Opm { , canph_ (active2Canonical(fluid.phaseUsage())) , cells_ (buildAllCells(grid.number_of_cells)) , ops_ (grid) + , wops_ (wells) , grav_ (gravityOperator(grid_, ops_, geo_)) , rq_ (fluid.numPhases()) { @@ -215,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 Tri; + + std::vector 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() { diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.hpp b/opm/autodiff/FullyImplicitBlackoilSolver.hpp index e1f539f0a..7773e50bf 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.hpp @@ -80,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 }; @@ -95,6 +101,7 @@ namespace Opm { const std::vector canph_; const std::vector cells_; // All grid cells HelperOps ops_; + const WellOps wops_; const M grav_; std::vector rq_; From 26e965e858d2b7de564a3d4fd53f262bae2c733d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Fri, 24 May 2013 19:15:14 +0200 Subject: [PATCH 3/3] Add an operator<< for outputting vector results. Use it out print reservoir and well pressures at end of solution. --- examples/test_implicit_ad.cpp | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/examples/test_implicit_ad.cpp b/examples/test_implicit_ad.cpp index 09e0ca474..acdffd4a0 100644 --- a/examples/test_implicit_ad.cpp +++ b/examples/test_implicit_ad.cpp @@ -72,6 +72,15 @@ namespace { return wells; } + + template + Ostream& + operator<<(Ostream& os, const std::vector& v) + { + std::copy(v.begin(), v.end(), std::ostream_iterator(os, " ")); + + return os; + } } int @@ -103,17 +112,8 @@ main(int argc, char* argv[]) solver.step(1.0, state, well_state); -#if 0 - Opm::WellState well_state; - well_state.init(wells, state); - - ps.solve(1.0, state, well_state); - - std::copy(state.pressure().begin(), state.pressure().end(), std::ostream_iterator(std::cout, " ")); - std::cout << std::endl; - std::copy(well_state.bhp().begin(), well_state.bhp().end(), std::ostream_iterator(std::cout, " ")); - std::cout << std::endl; -#endif + std::cout << state.pressure() << '\n' + << well_state.bhp() << '\n'; return 0; }