diff --git a/examples/sim_2p_fincomp_ad.cpp b/examples/sim_2p_fincomp_ad.cpp index 17f75425b..91b026fad 100644 --- a/examples/sim_2p_fincomp_ad.cpp +++ b/examples/sim_2p_fincomp_ad.cpp @@ -6,6 +6,9 @@ #include #include #include + +#include +#include #include #include #include @@ -24,8 +27,17 @@ int main (int argc, char** argv) try { + using namespace Opm; parameter::ParameterGroup param(argc, argv, false); + + bool use_deck = param.has("deck_filename"); + if (!use_deck) { + OPM_THROW(std::runtime_error, "FullyImplicitTwoPhaseSolver cannot run without deckfile."); + } + double gravity[3] = { 0.0 }; + std::string deck_filename = param.get("deck_filename"); + EclipseGridParser deck = EclipseGridParser(deck_filename); int nx = param.getDefault("nx", 30); int ny = param.getDefault("ny", 1); int nz = 1; @@ -47,14 +59,20 @@ try SaturationPropsBasic::RelPermFunc rel_perm_func = SaturationPropsBasic::Linear; IncompPropsAdBasic props(num_phases, rel_perm_func, density, viscosity, porosity, permeability, grid.dimensions, num_cells); - std::vector omega; +/* std::vector src(num_cells, 0.0); src[0] = 1. / day; src[num_cells-1] = -1. / day; + */ FlowBCManager bcs; LinearSolverUmfpack linsolver; - FullyImplicitTwoPhaseSolver solver(grid, props, linsolver); + TwophaseState state; + state.init(grid, 2); + WellState well_state; + WellsManager wells(deck, grid, props.permeability()); + well_state.init(wells.c_wells(), state); + FullyImplicitTwoPhaseSolver solver(grid, props, *wells.c_wells(), linsolver, gravity); std::vector porevol; Opm::computePorevolume(grid, props.porosity(), porevol); const double dt = param.getDefault("dt", 10.) * day; @@ -63,10 +81,9 @@ try for (int cell = 0; cell < num_cells; ++cell) { allcells[cell] = cell; } - TwophaseState state; - state.init(grid, 2); - - //initial sat + std::vector src; // empty src term. + gravity[2] = param.getDefault("gravity", 0.0); + //initial sat for (int c = 0; c < num_cells; ++c) { state.saturation()[2*c] = 0.2; state.saturation()[2*c+1] = 0.8; @@ -74,8 +91,6 @@ try std::vector p(num_cells, 100*Opm::unit::barsa); state.pressure() = p; 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()); @@ -84,7 +99,7 @@ try dm["pressure"] = &state.pressure(); Opm::writeVtkData(grid, dm, vtkfile); for (int i = 0; i < num_time_steps; ++i) { - solver.step(dt, state, src); + solver.step(dt, state, src, well_state); vtkfilename.str(""); vtkfilename << "sim_2p_fincomp_" << std::setw(3) << std::setfill('0') << i + 1 << ".vtu"; std::ofstream vtkfile(vtkfilename.str().c_str()); @@ -92,7 +107,7 @@ try dm["saturation"] = &state.saturation(); dm["pressure"] = &state.pressure(); Opm::writeVtkData(grid, dm, vtkfile); - } + } } catch (const std::exception &e) { std::cerr << "Program threw an exception: " << e.what() << "\n"; diff --git a/opm/polymer/fullyimplicit/.FullyImplicitTwoPhaseSolver.cpp.swp b/opm/polymer/fullyimplicit/.FullyImplicitTwoPhaseSolver.cpp.swp new file mode 100644 index 000000000..776e79ab9 Binary files /dev/null and b/opm/polymer/fullyimplicit/.FullyImplicitTwoPhaseSolver.cpp.swp differ diff --git a/opm/polymer/fullyimplicit/FullyImplicitTwoPhaseSolver.cpp b/opm/polymer/fullyimplicit/FullyImplicitTwoPhaseSolver.cpp index cf03d56f3..b5211a366 100644 --- a/opm/polymer/fullyimplicit/FullyImplicitTwoPhaseSolver.cpp +++ b/opm/polymer/fullyimplicit/FullyImplicitTwoPhaseSolver.cpp @@ -22,6 +22,14 @@ #include namespace Opm { +typedef AutoDiffBlock ADB; +typedef ADB::V V; +typedef ADB::M M; +typedef Eigen::Array DataBlock; + namespace { std::vector @@ -36,20 +44,37 @@ namespace { double operator()(double x) const { return std::max(std::min(x, 1.0), 0.0); } }; + + V computePerfPress(const UnstructuredGrid& g, + const Wells& wells, + const V& rho, + const double grav) + { + const int nw = wells.number_of_wells; + const int nperf = wells.well_connpos[nw]; + const int dim = g.dimensions; + + V wdp = V::Zero(nperf, 1); + assert(wdp.size() == rho.size()); + + for (int w = 0; w < nw; ++w) { + const double ref_depth = wells.depth_ref[w]; + for (int j = wells.well_connpos[w]; j < wells.well_connpos[nw + 1]; ++j) { + const int cell = wells.well_cells[j]; + const double cell_depth = g.cell_centroids[dim * cell + dim - 1]; + wdp(j) = rho(j) * grav * (cell_depth - ref_depth); + } + } + + return wdp; + } + }//anonymous namespace -typedef AutoDiffBlock ADB; -typedef ADB::V V; -typedef ADB::M M; -typedef Eigen::Array DataBlock; - @@ -57,25 +82,61 @@ typedef Eigen::Array(fluid.numPhases(), ADB::null())) + , wops_(wells) + , mob_ (fluid.numPhases(), ADB::null()) + , residual_ ({ std::vector(fluid.numPhases(), ADB::null()), + ADB::null(), + ADB::null() }) { } + FullyImplicitTwoPhaseSolver:: + 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 FullyImplicitTwoPhaseSolver:: step(const double dt, TwophaseState& x, - const std::vector& src) + const std::vector& src, + WellState& xw) { V pvol(grid_.number_of_cells); @@ -88,12 +149,12 @@ typedef Eigen::Array atol; while (resTooLarge && (it < maxit)) { const V dx = solveJacobianSystem(); - updateState(dx, x); + updateState(dx, x, xw); - assemble(pvdt, old_state, x, src); + assemble(pvdt, old_state, x, xw, src); const double r = residualNorm(); @@ -129,6 +190,7 @@ typedef Eigen::Array bpat(np ,nc); + bpat.push_back(xw.bhp().size()); SolutionState state(np); // Pressure. assert (not x.pressure().empty()); const V p = Eigen::Map(& x.pressure()[0], nc); - state.pressure = ADB::constant(p); + state.pressure = ADB::constant(p, bpat); // Saturation. assert (not x.saturation().empty()); const DataBlock s_all = Eigen::Map(& x.saturation()[0], nc, np); for (int phase = 0; phase < np; ++phase) { - state.saturation[phase] = ADB::constant(s_all.col(phase)); + state.saturation[phase] = ADB::constant(s_all.col(phase), bpat); // state.saturation[1] = ADB::constant(s_all.col(1)); } + + // BHP + assert (not x.bhp().empty()); + const V bhp = Eigen::Map(& xw.bhp()[0], xw.bhp().size()); + state.bhp = ADB::constant(bhp, bpat); + return state; } @@ -164,7 +235,8 @@ typedef Eigen::Array(& xw.bhp()[0], xw.bhp().size()); + vars0.push_back(bhp); std::vector vars = ADB::variables(vars0); @@ -201,6 +277,8 @@ typedef Eigen::Array& src) { // Create the primary variables. - const SolutionState state = variableState(x); + const SolutionState state = variableState(x, xw); // -------- Mass balance equations -------- const V trans = subset(transmissibility(), ops_.internal_faces); const std::vector kr = computeRelPerm(state); for (int phase = 0; phase < fluid_.numPhases(); ++phase) { const ADB mflux = computeMassFlux(phase, trans, kr, state); - ADB source = accumSource(phase, kr, src); - residual_[phase] = + residual_.mass_balance[phase] = pvdt*(state.saturation[phase] - old_state.saturation[phase]) - + ops_.div*mflux - source; + + ops_.div*mflux; + if (not src.empty()) { + ADB source = accumSource(phase, kr, src); + residual_.mass_balance[phase] = residual_.mass_balance[phase] - source; + } } + // -------- Well equation, and well contributions to the mass balance equations -------- + + // Contribution to mass balance will have to wait. + const int nc = grid_.number_of_cells; + const int np = wells_.number_of_phases; + const int nw = wells_.number_of_wells; + const int nperf = wells_.well_connpos[nw]; + + const std::vector well_cells(wells_.well_cells, wells_.well_cells + nperf); + const V transw = Eigen::Map(wells_.WI, nperf); + + const ADB& bhp = state.bhp; + + const DataBlock well_s = wops_.w2p * Eigen::Map(wells_.comp_frac, nw, np).matrix(); + + // Extract variables for perforation cell pressures + // and corresponding perforation well pressures. + const ADB p_perfcell = subset(state.pressure, well_cells); + // Finally construct well perforation pressures and well flows. + + // Compute well pressure differentials. + // Construct pressure difference vector for wells. + const int dim = grid_.dimensions; + const double* g = gravity(); + if (g) { + // Guard against gravity in anything but last dimension. + for (int dd = 0; dd < dim - 1; ++dd) { + assert(g[dd] == 0.0); + } + } + ADB cell_rho_total = ADB::constant(V::Zero(nc), state.pressure.blockPattern()); + for (int phase = 0; phase < 2; ++phase) { + const ADB cell_rho = fluidDensity(phase, state.pressure); + cell_rho_total += state.saturation[phase] * cell_rho; + } + ADB inj_rho_total = ADB::constant(V::Zero(nperf), state.pressure.blockPattern()); + assert(np == wells_.number_of_phases); + const DataBlock compi = Eigen::Map(wells_.comp_frac, nw, np); + for (int phase = 0; phase < 2; ++phase) { + const ADB cell_rho = fluidDensity(phase, state.pressure); + const V fraction = compi.col(phase); + inj_rho_total += (wops_.w2p * fraction.matrix()).array() * subset(cell_rho, well_cells); + } + const V rho_perf_cell = subset(cell_rho_total, well_cells).value(); + const V rho_perf_well = inj_rho_total.value(); + V prodperfs = V::Constant(nperf, -1.0); + for (int w = 0; w < nw; ++w) { + if (wells_.type[w] == PRODUCER) { + std::fill(prodperfs.data() + wells_.well_connpos[w], + prodperfs.data() + wells_.well_connpos[w+1], 1.0); + } + } + const Selector producer(prodperfs); + const V rho_perf = producer.select(rho_perf_cell, rho_perf_well); + const V well_perf_dp = computePerfPress(grid_, wells_, rho_perf, g ? g[dim-1] : 0.0); + + const ADB p_perfwell = wops_.w2p * bhp + well_perf_dp; + const ADB nkgradp_well = transw * (p_perfcell - p_perfwell); + + const Selector cell_to_well_selector(nkgradp_well.value()); + ADB well_rates_all = ADB::constant(V::Zero(nw*np), state.bhp.blockPattern()); + ADB perf_total_mob = subset(mob_[0], well_cells) + + subset(mob_[1], well_cells); + std::vector well_contribs(np, ADB::null()); + std::vector well_perf_rates(np, ADB::null()); + for (int phase = 0; phase < np; ++phase) { + // const ADB& cell_b = rq_[phase].b; + // const ADB perf_b = subset(cell_b, well_cells); + const ADB& cell_mob = mob_[phase]; + const V well_fraction = compi.col(phase); + // Using total mobilities for all phases for injection. + const ADB perf_mob_injector = (wops_.w2p * well_fraction.matrix()).array() * perf_total_mob; + const ADB perf_mob = producer.select(subset(cell_mob, well_cells), + perf_mob_injector); + const ADB perf_flux = perf_mob * (nkgradp_well); // No gravity term for perforations. + well_contribs[phase] = superset(perf_flux, well_cells, nc); + residual_.mass_balance[phase] += well_contribs[phase]; + } + + + // Handling BHP and SURFACE_RATE wells. + V bhp_targets(nw); + for (int w = 0; w < nw; ++w) { + const WellControls* wc = wells_.ctrls[w]; + if (wc->type[wc->current] == BHP) { + bhp_targets[w] = wc->target[wc->current]; + } else { + OPM_THROW(std::runtime_error, "Can only handle BHP type controls."); + } + } + const ADB bhp_residual = bhp - bhp_targets; + // Choose bhp residual for positive bhp targets. + residual_.well_eq = bhp_residual; } @@ -265,7 +440,7 @@ typedef Eigen::Array(& insrc[0], grid_.number_of_cells); // compute the out-fracflow. - ADB f_out = computeFracFlow(phase, kr); + ADB f_out = computeFracFlow(phase); // compute the in-fracflow. V f_in; if (phase == 1) { @@ -281,15 +456,10 @@ typedef Eigen::Array& kr) const + FullyImplicitTwoPhaseSolver::computeFracFlow(const int phase) const { - const double* mus = fluid_.viscosity(); - ADB mob_phase = kr[phase] / V::Constant(kr[phase].size(), 1, mus[phase]); - ADB mob_wat = kr[0] / V::Constant(kr[0].size(), 1, mus[0]); - ADB mob_oil= kr[1] / V::Constant(kr[1].size(), 1, mus[1]); - ADB total_mob = mob_wat + mob_oil; - ADB f = mob_phase / total_mob; + ADB total_mob = mob_[0] + mob_[1]; + ADB f = mob_[phase] / total_mob; return f; } @@ -305,13 +475,15 @@ typedef Eigen::Array matr = mass_res.derivative()[0]; - V dx(V::Zero(mass_res.size())); + const ADB mass_res = vertcat(residual_.mass_balance[0], residual_.mass_balance[1]); + const ADB total_res = collapseJacs(vertcat(mass_res, residual_.well_eq)); + + const Eigen::SparseMatrix matr = total_res.derivative()[0]; + V dx(V::Zero(total_res.size())); Opm::LinearSolverInterface::LinearSolverReport rep = linsolver_.solve(matr.rows(), matr.nonZeros(), matr.outerIndexPtr(), matr.innerIndexPtr(), matr.valuePtr(), - mass_res.value().data(), dx.data()); + total_res.value().data(), dx.data()); if (!rep.converged) { OPM_THROW(std::runtime_error, "FullyImplicitBlackoilSolver::solveJacobianSystem(): " @@ -324,11 +496,13 @@ typedef Eigen::Array(&state.pressure()[0], nc); const V p = p_old - dp; @@ -362,6 +540,10 @@ typedef Eigen::Array(&well_state.bhp()[0], nw, 1); + const V bhp = bhp_old - dbhp; + std::copy(&p[0], &p[0] + nc, well_state.bhp().begin()); } @@ -389,21 +571,28 @@ typedef Eigen::Array& kr , - const SolutionState& state ) const + const SolutionState& state ) { // 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; + // ADB& mob = mob_[phase]; + mob_[phase] = kr[phase] / V::Constant(kr[phase].size(), 1, mus[phase]); + // ADB mob = kr[phase] / V::Constant(kr[phase].size(), 1, mus[phase]); + V z(grid_.number_of_cells); + for (int c = 0; c < grid_.number_of_cells; ++c) { + z(c) = grid_.cell_centroids[c * 3 + 2]; + } + const double* grav = gravity(); + const ADB rho = fluidDensity(phase, state.pressure); + const ADB rhoavg = ops_.caver * rho; + const ADB dp = ops_.ngrad * state.pressure - grav[2] * (rhoavg * (ops_.ngrad * z.matrix())); const ADB head = trans * dp; UpwindSelector upwind(grid_, ops_, head.value()); - return upwind.select(mob) * head; + return upwind.select(mob_[phase]) * head; } @@ -415,13 +604,14 @@ typedef Eigen::Array::const_iterator - b = residual_.begin(), - e = residual_.end(); + b = residual_.mass_balance.begin(), + e = residual_.mass_balance.end(); b != e; ++b) { r = std::max(r, (*b).value().matrix().norm()); } - + + r = std::max(r, residual_.well_eq.value().matrix().norm()); return r; } @@ -442,7 +632,15 @@ typedef Eigen::Array #include #include - +#include +#include struct UnstructuredGrid; namespace Opm { @@ -21,11 +22,14 @@ namespace Opm { public: FullyImplicitTwoPhaseSolver(const UnstructuredGrid& grid, const IncompPropsAdInterface& fluid, - const LinearSolverInterface& linsolver); + const Wells& wells, + const LinearSolverInterface& linsolver, + const double* gravity); void step(const double dt, TwophaseState& state, - const std::vector& src); + const std::vector& src, + WellState& wstate); private: typedef AutoDiffBlock ADB; typedef ADB::V V; @@ -38,34 +42,60 @@ namespace Opm { SolutionState(const int np); ADB pressure; std::vector saturation; + ADB bhp; }; + /* + struct Source { + Wells& wells; + std::vector src; + } source; + */ + struct WellOps { + WellOps(const Wells& wells); + M w2p; // well->perf + M p2w; // perf->well + }; + const UnstructuredGrid& grid_; const IncompPropsAdInterface& fluid_; + const Wells& wells_; const LinearSolverInterface& linsolver_; + const double* grav_; const std::vector cells_; HelperOps ops_; - std::vector residual_; + const WellOps wops_; + + std::vector mob_; + + struct { + std::vector mass_balance; + ADB well_flux_eq; + ADB well_eq; + } residual_; SolutionState - constantState(const TwophaseState& x); + constantState(const TwophaseState& x, + const WellState& xw); SolutionState - variableState(const TwophaseState& x); + variableState(const TwophaseState& x, + const WellState& xw); void assemble(const V& pvdt, const SolutionState& old_state, const TwophaseState& x, + const WellState& xw, const std::vector& src); V solveJacobianSystem() const; - void updateState(const V& dx, - TwophaseState& x) const; + void updateState(const V& dx, + TwophaseState& x, + WellState& xw)const; std::vector computeRelPerm(const SolutionState& state) const; V transmissibility() const; ADB - computeFracFlow(int phase, - const std::vector& kr) const; + computeFracFlow(const int phase) const; ADB accumSource(const int phase, const std::vector& kr, @@ -74,7 +104,7 @@ namespace Opm { computeMassFlux(const int phase, const V& trans, const std::vector& kr, - const SolutionState& state) const; + const SolutionState& state); double residualNorm() const; @@ -86,6 +116,11 @@ namespace Opm { fluidDensity(const int phase) const; ADB transMult(const ADB& p) const; + + ADB + fluidDensity(const int phase, + const ADB& p) const; + const double* gravity() const { return grav_; } }; } // namespace Opm #endif// OPM_FULLYIMPLICITTWOPHASESOLVER_HEADER_INCLUDED