diff --git a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp index 289912d32..672a84077 100644 --- a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp +++ b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp @@ -834,10 +834,10 @@ namespace { e = residual_.mass_balance.end(); b != e; ++b) { - r = std::max(r, (*b).value().matrix().norm()); + r = std::max(r, (*b).value().matrix().lpNorm()); } - r = std::max(r, residual_.well_flux_eq.value().matrix().norm()); - r = std::max(r, residual_.well_eq.value().matrix().norm()); + r = std::max(r, residual_.well_flux_eq.value().matrix().lpNorm()); + r = std::max(r, residual_.well_eq.value().matrix().lpNorm()); return r; } @@ -848,8 +848,8 @@ namespace { ADB FullyImplicitCompressiblePolymerSolver::fluidViscosity(const int phase, - const ADB& p , - const std::vector& cells) const + const ADB& p , + const std::vector& cells) const { const ADB null = ADB::constant(V::Zero(grid_.number_of_cells, 1), p.blockPattern()); switch (phase) { @@ -869,8 +869,8 @@ namespace { ADB FullyImplicitCompressiblePolymerSolver::fluidReciprocFVF(const int phase, - const ADB& p , - const std::vector& cells) const + const ADB& p , + const std::vector& cells) const { const ADB null = ADB::constant(V::Zero(grid_.number_of_cells, 1), p.blockPattern()); switch (phase) { @@ -890,8 +890,8 @@ namespace { ADB FullyImplicitCompressiblePolymerSolver::fluidDensity(const int phase, - const ADB& p , - const std::vector& cells) const + const ADB& p , + const std::vector& cells) const { const double* rhos = fluid_.surfaceDensity(); ADB b = fluidReciprocFVF(phase, p, cells); diff --git a/opm/polymer/fullyimplicit/FullyImplicitTwoPhaseSolver.cpp b/opm/polymer/fullyimplicit/FullyImplicitTwoPhaseSolver.cpp deleted file mode 100644 index d258b4bc3..000000000 --- a/opm/polymer/fullyimplicit/FullyImplicitTwoPhaseSolver.cpp +++ /dev/null @@ -1,649 +0,0 @@ -/**/ - -#include - -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include -namespace Opm { - -typedef AutoDiffBlock ADB; -typedef ADB::V V; -typedef ADB::M M; -typedef Eigen::Array DataBlock; - -namespace { - - std::vector - buildAllCells(const int nc) - { - std::vector all_cells(nc); - for (int c = 0; c < nc; ++c) { all_cells[c] = c; } - - return all_cells; - } - struct Chop01 { - 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 - - - - - - - - - - FullyImplicitTwoPhaseSolver:: - FullyImplicitTwoPhaseSolver(const UnstructuredGrid& grid, - const IncompPropsAdInterface& fluid, - const Wells& wells, - const LinearSolverInterface& linsolver, - const double* gravity) - : grid_ (grid) - , fluid_(fluid) - , wells_(wells) - , linsolver_(linsolver) - , grav_(gravity) - , cells_ (buildAllCells(grid.number_of_cells)) - , ops_(grid) - , 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, - WellState& xw) - { - - V pvol(grid_.number_of_cells); - // Pore volume - const typename V::Index nc = grid_.number_of_cells; - V rho = V::Constant(pvol.size(), 1, *fluid_.porosity()); - std::transform(grid_.cell_volumes, grid_.cell_volumes + nc, - rho.data(), pvol.data(), - std::multiplies()); - - const V pvdt = pvol / dt; - - const SolutionState old_state = constantState(x, xw); - const double atol = 1.0e-12; - const double rtol = 5.0e-8; - const int maxit = 15; - - assemble(pvdt, old_state, x, xw, src); - - const double r0 = residualNorm(); - int it = 0; - std::cout << "\nIteration Residual\n" - << std::setw(9) << it << std::setprecision(9) - << std::setw(18) << r0 << std::endl; - bool resTooLarge = r0 > atol; - while (resTooLarge && (it < maxit)) { - const V dx = solveJacobianSystem(); - updateState(dx, x, xw); - - assemble(pvdt, old_state, x, xw, src); - - const double r = residualNorm(); - - resTooLarge = (r > atol) && (r > rtol*r0); - - it += 1; - std::cout << std::setw(9) << it << std::setprecision(9) - << std::setw(18) << r << std::endl; - } - - if (resTooLarge) { - std::cerr << "Failed to compute converged solution in " << it << " iterations. Ignoring!\n"; - // OPM_THROW(std::runtime_error, "Failed to compute converged solution in " << it << " iterations."); - } - } - - - - - - FullyImplicitTwoPhaseSolver::SolutionState::SolutionState(const int np) - : pressure ( ADB::null()) - , saturation (np, ADB::null()) - , bhp ( ADB::null()) - { - } - - - - - - FullyImplicitTwoPhaseSolver::SolutionState - FullyImplicitTwoPhaseSolver::constantState(const TwophaseState& x, - const WellState& xw) - { - const int nc = grid_.number_of_cells; - const int np = x.numPhases(); - std::vector 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, 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), 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; - } - - - - - - FullyImplicitTwoPhaseSolver::SolutionState - FullyImplicitTwoPhaseSolver::variableState(const TwophaseState& x, - const WellState& xw) - { - const int nc = grid_.number_of_cells; - const int np = x.numPhases(); - - std::vector vars0; - vars0.reserve(np); - - // Initial pressure. - assert (not x.pressure().empty()); - const V p = Eigen::Map(& x.pressure()[0], nc); - vars0.push_back(p); - - // Initial saturation. - assert (not x.saturation().empty()); - const DataBlock s_all = Eigen::Map(& x.saturation()[0], nc, np); - const V sw = s_all.col(0); - vars0.push_back(sw); - - // Initial Bottom-hole Pressure - assert (not xw.bhp().empty()); - const V bhp = Eigen::Map(& xw.bhp()[0], xw.bhp().size()); - vars0.push_back(bhp); - - std::vector vars = ADB::variables(vars0); - - SolutionState state(np); - - // Pressure. - int nextvar = 0; - state.pressure = vars[ nextvar++ ]; - - // Saturation. - const std::vector& bpat = vars[0].blockPattern(); - { - ADB so = ADB::constant(V::Ones(nc, 1), bpat); - ADB& sw = vars[ nextvar++ ]; - state.saturation[0] = sw; - so = so - sw; - state.saturation[1] = so; - } - // Bottom-hole pressure. - state.pressure = vars[ nextvar++ ]; - - assert(nextvar == int(vars.size())); - - return state; - } - - - - - - void - FullyImplicitTwoPhaseSolver:: - assemble(const V& pvdt, - const SolutionState& old_state, - const TwophaseState& x , - const WellState& xw, - const std::vector& src) - { - // Create the primary variables. - 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); - residual_.mass_balance[phase] = - pvdt*(state.saturation[phase] - old_state.saturation[phase]) - + ops_.div*mflux; - if (not src.empty()) { - ADB source = accumSource(phase, kr, src); - residual_.mass_balance[phase] = residual_.mass_balance[phase] - source; - } - } -#if 0 - // -------- 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; -#endif - } - - - - - - - ADB - FullyImplicitTwoPhaseSolver::accumSource(const int phase, - const std::vector& kr, - const std::vector& src) const - { - //extract the source to out and in source. - std::vector outsrc; - std::vector insrc; - std::vector::const_iterator it; - for (it = src.begin(); it != src.end(); ++it) { - if (*it < 0) { - outsrc.push_back(*it); - insrc.push_back(0.0); - } else if (*it > 0) { - insrc.push_back(*it); - outsrc.push_back(0.0); - } else { - outsrc.emplace_back(0); - insrc.emplace_back(0); - } - } - const V source = Eigen::Map(& src[0], grid_.number_of_cells); - const V outSrc = Eigen::Map(& outsrc[0], grid_.number_of_cells); - const V inSrc = Eigen::Map(& insrc[0], grid_.number_of_cells); - - // compute the out-fracflow. - ADB f_out = computeFracFlow(phase); - // compute the in-fracflow. - V f_in; - if (phase == 1) { - f_in = V::Zero(grid_.number_of_cells); - } else if (phase == 0) { - f_in = V::Ones(grid_.number_of_cells); - } - return f_out * outSrc + f_in * inSrc; - } - - - - - - ADB - FullyImplicitTwoPhaseSolver::computeFracFlow(const int phase) const - { - ADB total_mob = mob_[0] + mob_[1]; - ADB f = mob_[phase] / total_mob; - - return f; - } - - - - - - V - FullyImplicitTwoPhaseSolver::solveJacobianSystem() const - { - const int np = fluid_.numPhases(); - if (np != 2) { - OPM_THROW(std::logic_error, "Only two-phase ok in FullyImplicitTwoPhaseSolver."); - } - 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(), - total_res.value().data(), dx.data()); - if (!rep.converged) { - OPM_THROW(std::runtime_error, - "FullyImplicitBlackoilSolver::solveJacobianSystem(): " - "Linear solver convergence failure."); - } - return dx; - } - - - - - - void FullyImplicitTwoPhaseSolver::updateState(const V& dx, - TwophaseState& state, - WellState& well_state) const - { - const int np = fluid_.numPhases(); - const int nc = grid_.number_of_cells; - const int nw = wells_.number_of_wells; - const V null; - assert(null.size() == 0); - const V zero = V::Zero(nc); - const V one = V::Constant(nc, 1.0); - - // Extract parts of dx corresponding to each part. - const V dp = subset(dx, Span(nc)); - int varstart = nc; - const V dsw = subset(dx, Span(nc, 1, varstart)); - varstart += dsw.size(); - const V dbhp = subset(dx, Span(nc, 1, varstart)); - varstart += dbhp.size(); - - assert(varstart == dx.size()); - - - // Pressure update. - const V p_old = Eigen::Map(&state.pressure()[0], nc); - 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(& state.saturation()[0], nc, np); - V so = one; - const V sw_old = s_old.col(0); - const V dsw_limited = sign(dsw) * dsw.abs().min(dsmax); - const V sw = (sw_old - dsw_limited).unaryExpr(Chop01()); - so -= sw; - for (int c = 0; c < nc; ++c) { - state.saturation()[c*np] = sw[c]; - } - for (int c = 0; c < nc; ++c) { - state.saturation()[c*np + 1] = so[c]; - } - - // Bhp update. - const V bhp_old = Eigen::Map(&well_state.bhp()[0], nw, 1); - const V bhp = bhp_old - dbhp; - std::copy(&p[0], &p[0] + nc, well_state.bhp().begin()); - } - - - - - - std::vector - FullyImplicitTwoPhaseSolver::computeRelPerm(const SolutionState& state) const - { - - const ADB sw = state.saturation[0]; - const ADB so = state.saturation[1]; - - return fluid_.relperm(sw, so, cells_); - } - - - - - - - - - - - ADB - FullyImplicitTwoPhaseSolver::computeMassFlux(const int phase , - const V& trans , - const std::vector& kr , - const SolutionState& state ) - { -// const ADB tr_mult = transMult(state.pressure); - const double* mus = fluid_.viscosity(); - // 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_[phase]) * head; - } - - - - - - double - FullyImplicitTwoPhaseSolver::residualNorm() const - { - double r = 0; - for (std::vector::const_iterator - 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; - } - - - - - - V - FullyImplicitTwoPhaseSolver::transmissibility() const - { - const V::Index nc = grid_.number_of_cells; - V htrans(grid_.cell_facepos[nc]); - V trans(grid_.cell_facepos[nc]); - UnstructuredGrid* ug = const_cast(& grid_); - tpfa_htrans_compute(ug, fluid_.permeability(), htrans.data()); - tpfa_trans_compute (ug, htrans.data() , trans.data()); - - return trans; - } - - ADB - FullyImplicitTwoPhaseSolver::fluidDensity(const int phase, - const ADB& p) const - { - const double* rhos = fluid_.surfaceDensity(); - ADB rho = ADB::constant(V::Constant(grid_.number_of_cells, 1, rhos[phase]), p.blockPattern()); - - return rho; - } - - - - -}//namespace Opm diff --git a/opm/polymer/fullyimplicit/FullyImplicitTwoPhaseSolver.hpp b/opm/polymer/fullyimplicit/FullyImplicitTwoPhaseSolver.hpp deleted file mode 100644 index 16ad6b37c..000000000 --- a/opm/polymer/fullyimplicit/FullyImplicitTwoPhaseSolver.hpp +++ /dev/null @@ -1,127 +0,0 @@ -/**/ - -#ifndef OPM_FULLYIMPLICITTWOPHASESOLVER_HEADER_INCLUDED -#define OPM_FULLYIMPLICITTWOPHASESOLVER_HEADER_INCLUDED - -#include -#include -#include -#include -#include -#include - -struct UnstructuredGrid; -namespace Opm { -// struct HelperOps; - class LinearSolverInterface; - class TwophaseState; - - - class FullyImplicitTwoPhaseSolver - { - public: - FullyImplicitTwoPhaseSolver(const UnstructuredGrid& grid, - const IncompPropsAdInterface& fluid, - // const Wells& wells, - const LinearSolverInterface& linsolver); - // const double* gravity); - - void step(const double dt, - TwophaseState& state, - const std::vector& src); -// WellState& wstate); - private: - typedef AutoDiffBlock ADB; - typedef ADB::V V; - typedef ADB::M M; - typedef Eigen::Array DataBlock; - struct SolutionState { - 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_; - // const WellOps wops_; - - std::vector mob_; - - struct { - std::vector mass_balance; - ADB well_flux_eq; - ADB well_eq; - } residual_; - - - SolutionState - constantState(const TwophaseState& x, - const WellState& xw); - SolutionState - 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, - WellState& xw)const; - std::vector - computeRelPerm(const SolutionState& state) const; - V - transmissibility() const; - ADB - computeFracFlow(const int phase) const; - ADB - accumSource(const int phase, - const std::vector& kr, - const std::vector& src) const; - ADB - computeMassFlux(const int phase, - const V& trans, - const std::vector& kr, - const SolutionState& state); - double - residualNorm() const; - - ADB - rockPorosity(const ADB& p) const; - ADB - rockPermeability(const ADB& p) const; - const double - 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 diff --git a/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp b/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp index b9a9b475d..d869b3133 100644 --- a/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp +++ b/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp @@ -1,4 +1,3 @@ -/**/ #include @@ -760,11 +759,11 @@ namespace { e = residual_.mass_balance.end(); b != e; ++b) { - r = std::max(r, (*b).value().matrix().norm()); + r = std::max(r, (*b).value().matrix().lpNorm()); } - r = std::max(r, residual_.well_flux_eq.value().matrix().norm()); - r = std::max(r, residual_.well_eq.value().matrix().norm()); + r = std::max(r, residual_.well_flux_eq.value().matrix().lpNorm()); + r = std::max(r, residual_.well_eq.value().matrix().lpNorm()); return r; } diff --git a/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp b/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp index 5e646a269..865332ec1 100644 --- a/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp +++ b/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp @@ -1,5 +1,3 @@ -/**/ - #ifndef OPM_FULLYIMPLICITTWOPHASEPOLYMERSOLVER_HEADER_INCLUDED #define OPM_FULLYIMPLICITTWOPHASEPOLYMERSOLVER_HEADER_INCLUDED diff --git a/opm/polymer/fullyimplicit/SimulatorFullyImplicitTwophase.cpp b/opm/polymer/fullyimplicit/SimulatorFullyImplicitTwophase.cpp deleted file mode 100644 index be908b0a4..000000000 --- a/opm/polymer/fullyimplicit/SimulatorFullyImplicitTwophase.cpp +++ /dev/null @@ -1,470 +0,0 @@ -/* - Copyright 2013 SINTEF ICT, Applied Mathematics. - - This file is part of the Open Porous Media project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#include -#include -#include - -#include -#include - -#include -#include -#include - -#include -#include -#include -#include - - -#include - -#include -#include - -#include -#include -#include - -#include -#include -#include -#include -namespace Opm -{ - - - - class SimulatorFullyImplicitTwophase::Impl - { - public: - Impl(const parameter::ParameterGroup& param, - const UnstructuredGrid& grid, - const IncompPropsAdInterface& props, - WellsManager& wells_manager, - LinearSolverInterface& linsolver, - std::vector& src, - const double* gravity); - - SimulatorReport run(SimulatorTimer& timer, - TwophaseState& state, - std::vector& src, - WellState& well_state); - - private: - - // Parameters for output. - bool output_; - bool output_vtk_; - std::string output_dir_; - int output_interval_; - // Parameters for well control - bool check_well_controls_; - int max_well_control_iterations_; - // Observed objects. - const UnstructuredGrid& grid_; - const IncompPropsAdInterface& props_; - WellsManager& wells_manager_; - const Wells* wells_; - const std::vector& src_; - const double* gravity_; - // Solvers - FullyImplicitTwoPhaseSolver solver_; - // Misc. data - std::vector allcells_; - }; - - - - - SimulatorFullyImplicitTwophase::SimulatorFullyImplicitTwophase(const parameter::ParameterGroup& param, - const UnstructuredGrid& grid, - const IncompPropsAdInterface& props, - WellsManager& wells_manager, - LinearSolverInterface& linsolver, - std::vector& src, - const double* gravity) - { - pimpl_.reset(new Impl(param, grid, props, wells_manager, linsolver, src, gravity)); - } - - - - - - SimulatorReport SimulatorFullyImplicitTwophase::run(SimulatorTimer& timer, - TwophaseState& state, - std::vector& src, - WellState& well_state) - { - return pimpl_->run(timer, state, src, well_state); - } - - - - static void outputStateVtk(const UnstructuredGrid& grid, - const Opm::TwophaseState& state, - const int step, - const std::string& output_dir) - { - // Write data in VTK format. - std::ostringstream vtkfilename; - vtkfilename << output_dir << "/vtk_files"; - boost::filesystem::path fpath(vtkfilename.str()); - try { - create_directories(fpath); - } - catch (...) { - OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath); - } - vtkfilename << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu"; - std::ofstream vtkfile(vtkfilename.str().c_str()); - if (!vtkfile) { - OPM_THROW(std::runtime_error, "Failed to open " << vtkfilename.str()); - } - Opm::DataMap dm; - dm["saturation"] = &state.saturation(); - dm["pressure"] = &state.pressure(); - std::vector cell_velocity; - Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity); - dm["velocity"] = &cell_velocity; - Opm::writeVtkData(grid, dm, vtkfile); - } - - - static void outputStateMatlab(const UnstructuredGrid& grid, - const Opm::TwophaseState& state, - const int step, - const std::string& output_dir) - { - Opm::DataMap dm; - dm["saturation"] = &state.saturation(); - dm["pressure"] = &state.pressure(); - std::vector cell_velocity; - Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity); - dm["velocity"] = &cell_velocity; - - // Write data (not grid) in Matlab format - for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) { - std::ostringstream fname; - fname << output_dir << "/" << it->first; - boost::filesystem::path fpath = fname.str(); - try { - create_directories(fpath); - } - catch (...) { - OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath); - } - fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt"; - std::ofstream file(fname.str().c_str()); - if (!file) { - OPM_THROW(std::runtime_error, "Failed to open " << fname.str()); - } - file.precision(15); - const std::vector& d = *(it->second); - std::copy(d.begin(), d.end(), std::ostream_iterator(file, "\n")); - } - } - - static void outputWellStateMatlab(const Opm::WellState& well_state, - const int step, - const std::string& output_dir) - { - Opm::DataMap dm; - dm["bhp"] = &well_state.bhp(); - dm["wellrates"] = &well_state.wellRates(); - - // Write data (not grid) in Matlab format - for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) { - std::ostringstream fname; - fname << output_dir << "/" << it->first; - boost::filesystem::path fpath = fname.str(); - try { - create_directories(fpath); - } - catch (...) { - OPM_THROW(std::runtime_error,"Creating directories failed: " << fpath); - } - fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt"; - std::ofstream file(fname.str().c_str()); - if (!file) { - OPM_THROW(std::runtime_error,"Failed to open " << fname.str()); - } - file.precision(15); - const std::vector& d = *(it->second); - std::copy(d.begin(), d.end(), std::ostream_iterator(file, "\n")); - } - } - - - - -/* - static void outputWaterCut(const Opm::Watercut& watercut, - const std::string& output_dir) - { - // Write water cut curve. - std::string fname = output_dir + "/watercut.txt"; - std::ofstream os(fname.c_str()); - if (!os) { - OPM_THROW(std::runtime_error, "Failed to open " << fname); - } - watercut.write(os); - } - static void outputWellReport(const Opm::WellReport& wellreport, - const std::string& output_dir) - { - // Write well report. - std::string fname = output_dir + "/wellreport.txt"; - std::ofstream os(fname.c_str()); - if (!os) { - OPM_THROW(std::runtime_error, "Failed to open " << fname); - } - wellreport.write(os); - } - */ - - - - SimulatorFullyImplicitTwophase::Impl::Impl(const parameter::ParameterGroup& param, - const UnstructuredGrid& grid, - const IncompPropsAdInterface& props, - WellsManager& wells_manager, - LinearSolverInterface& linsolver, - std::vector& src, - const double* gravity) - : grid_(grid), - props_(props), - wells_manager_(wells_manager), - wells_ (wells_manager.c_wells()), - src_ (src), - gravity_(gravity), - solver_(grid_, props_, *wells_manager.c_wells(), linsolver, gravity_) - - { - // For output. - output_ = param.getDefault("output", true); - if (output_) { - output_vtk_ = param.getDefault("output_vtk", true); - output_dir_ = param.getDefault("output_dir", std::string("output")); - // Ensure that output dir exists - boost::filesystem::path fpath(output_dir_); - try { - create_directories(fpath); - } - catch (...) { - OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath); - } - output_interval_ = param.getDefault("output_interval", 1); - } - // Well control related init. - check_well_controls_ = param.getDefault("check_well_controls", false); - max_well_control_iterations_ = param.getDefault("max_well_control_iterations", 10); - - // Misc init. - const int num_cells = grid.number_of_cells; - allcells_.resize(num_cells); - for (int cell = 0; cell < num_cells; ++cell) { - allcells_[cell] = cell; - } - } - - SimulatorReport SimulatorFullyImplicitTwophase::Impl::run(SimulatorTimer& timer, - TwophaseState& state, - std::vector& src, - WellState& well_state) - { - - // Initialisation. - std::vector porevol; - Opm::computePorevolume(grid_, props_.porosity(), porevol); - - // const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0); - std::vector initial_porevol = porevol; - - // Main simulation loop. - Opm::time::StopWatch solver_timer; - double stime = 0.0; - Opm::time::StopWatch step_timer; - Opm::time::StopWatch total_timer; - total_timer.start(); -#if 0 - // These must be changed for three-phase. - double init_surfvol[2] = { 0.0 }; - double inplace_surfvol[2] = { 0.0 }; - double tot_injected[2] = { 0.0 }; - double tot_produced[2] = { 0.0 }; - Opm::computeSaturatedVol(porevol, state.surfacevol(), init_surfvol); - Opm::Watercut watercut; - watercut.push(0.0, 0.0, 0.0); -#endif - std::vector fractional_flows; - std::vector well_resflows_phase; - if (wells_) { - well_resflows_phase.resize((wells_->number_of_phases)*(wells_->number_of_wells), 0.0); - } - std::fstream tstep_os; - if (output_) { - std::string filename = output_dir_ + "/step_timing.param"; - tstep_os.open(filename.c_str(), std::fstream::out | std::fstream::app); - } - while (!timer.done()) { - // Report timestep and (optionally) write state to disk. - step_timer.start(); - timer.report(std::cout); - if (output_ && (timer.currentStepNum() % output_interval_ == 0)) { - if (output_vtk_) { - outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_); - } - outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_); - outputWellStateMatlab(well_state,timer.currentStepNum(), output_dir_); - - } - - SimulatorReport sreport; - - bool well_control_passed = !check_well_controls_; - int well_control_iteration = 0; - do { - // Run solver. - solver_timer.start(); - std::vector initial_pressure = state.pressure(); - solver_.step(timer.currentStepLength(), state, src, well_state); - - // Stop timer and report. - solver_timer.stop(); - const double st = solver_timer.secsSinceStart(); - std::cout << "Fully implicit solver took: " << st << " seconds." << std::endl; - - stime += st; - sreport.pressure_time = st; - - // Optionally, check if well controls are satisfied. - if (check_well_controls_) { - Opm::computePhaseFlowRatesPerWell(*wells_, - well_state.perfRates(), - fractional_flows, - well_resflows_phase); - std::cout << "Checking well conditions." << std::endl; - // For testing we set surface := reservoir - well_control_passed = wells_manager_.conditionsMet(well_state.bhp(), well_resflows_phase, well_resflows_phase); - ++well_control_iteration; - if (!well_control_passed && well_control_iteration > max_well_control_iterations_) { - OPM_THROW(std::runtime_error, "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); - - // Update pore volumes if rock is compressible. - initial_porevol = porevol; - - // The reports below are geared towards two phases only. -#if 0 - // Report mass balances. - double injected[2] = { 0.0 }; - double produced[2] = { 0.0 }; - Opm::computeInjectedProduced(props_, state, transport_src, stepsize, - injected, produced); - Opm::computeSaturatedVol(porevol, state.surfacevol(), inplace_surfvol); - tot_injected[0] += injected[0]; - tot_injected[1] += injected[1]; - tot_produced[0] += produced[0]; - tot_produced[1] += produced[1]; - std::cout.precision(5); - const int width = 18; - std::cout << "\nMass balance report.\n"; - std::cout << " Injected surface volumes: " - << std::setw(width) << injected[0] - << std::setw(width) << injected[1] << std::endl; - std::cout << " Produced surface volumes: " - << std::setw(width) << produced[0] - << std::setw(width) << produced[1] << std::endl; - std::cout << " Total inj surface volumes: " - << std::setw(width) << tot_injected[0] - << std::setw(width) << tot_injected[1] << std::endl; - std::cout << " Total prod surface volumes: " - << std::setw(width) << tot_produced[0] - << std::setw(width) << tot_produced[1] << std::endl; - const double balance[2] = { init_surfvol[0] - inplace_surfvol[0] - tot_produced[0] + tot_injected[0], - init_surfvol[1] - inplace_surfvol[1] - tot_produced[1] + tot_injected[1] }; - std::cout << " Initial - inplace + inj - prod: " - << std::setw(width) << balance[0] - << std::setw(width) << balance[1] - << std::endl; - std::cout << " Relative mass error: " - << std::setw(width) << balance[0]/(init_surfvol[0] + tot_injected[0]) - << std::setw(width) << balance[1]/(init_surfvol[1] + tot_injected[1]) - << std::endl; - std::cout.precision(8); - - // Make well reports. - watercut.push(timer.currentTime() + timer.currentStepLength(), - produced[0]/(produced[0] + produced[1]), - tot_produced[0]/tot_porevol_init); - if (wells_) { - wellreport.push(props_, *wells_, - state.pressure(), state.surfacevol(), state.saturation(), - timer.currentTime() + timer.currentStepLength(), - well_state.bhp(), well_state.perfRates()); - } -#endif - sreport.total_time = step_timer.secsSinceStart(); - if (output_) { - sreport.reportParam(tstep_os); - - if (output_vtk_) { - outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_); - } - outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_); - outputWellStateMatlab(well_state,timer.currentStepNum(), output_dir_); -#if 0 - outputWaterCut(watercut, output_dir_); - if (wells_) { - outputWellReport(wellreport, output_dir_); - } -#endif - tstep_os.close(); - } - - // advance to next timestep before reporting at this location - ++timer; - - // write an output file for later inspection - } - - total_timer.stop(); - - SimulatorReport report; - report.pressure_time = stime; - report.transport_time = 0.0; - report.total_time = total_timer.secsSinceStart(); - return report; - } - - - - - -} // namespace Opm diff --git a/opm/polymer/fullyimplicit/SimulatorFullyImplicitTwophase.hpp b/opm/polymer/fullyimplicit/SimulatorFullyImplicitTwophase.hpp deleted file mode 100644 index c70a65e77..000000000 --- a/opm/polymer/fullyimplicit/SimulatorFullyImplicitTwophase.hpp +++ /dev/null @@ -1,90 +0,0 @@ -/* - Copyright 2013 SINTEF ICT, Applied Mathematics. - - This file is part of the Open Porous Media project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#ifndef OPM_SIMULATORFULLYIMPLICITTWOPHASE_HEADER_INCLUDED -#define OPM_SIMULATORFULLYIMPLICITTWOPHASE_HEADER_INCLUDED - -#include -#include - -struct UnstructuredGrid; - -namespace Opm -{ - namespace parameter { class ParameterGroup; } - class IncompPropsAdInterface; - class LinearSolverInterface; - class SimulatorTimer; - class TwophaseState; - class WellsManager; - class WellState; - struct SimulatorReport; - - /// Class collecting all necessary components for a two-phase simulation. - class SimulatorFullyImplicitTwophase - { - public: - /// Initialise from parameters and objects to observe. - /// \param[in] param parameters, this class accepts the following: - /// parameter (default) effect - /// ----------------------------------------------------------- - /// output (true) write output to files? - /// output_dir ("output") output directoty - /// output_interval (1) output every nth step - /// nl_pressure_residual_tolerance (0.0) pressure solver residual tolerance (in Pascal) - /// nl_pressure_change_tolerance (1.0) pressure solver change tolerance (in Pascal) - /// nl_pressure_maxiter (10) max nonlinear iterations in pressure - /// nl_maxiter (30) max nonlinear iterations in transport - /// nl_tolerance (1e-9) transport solver absolute residual tolerance - /// num_transport_substeps (1) number of transport steps per pressure step - /// use_segregation_split (false) solve for gravity segregation (if false, - /// segregation is ignored). - /// - /// \param[in] grid grid data structure - /// \param[in] props fluid and rock properties - /// \param[in] linsolver linear solver - SimulatorFullyImplicitTwophase(const parameter::ParameterGroup& param, - const UnstructuredGrid& grid, - const IncompPropsAdInterface& props, - WellsManager& wells_manager, - LinearSolverInterface& linsolver, - std::vector& src, - const double* gravity); - - /// Run the simulation. - /// This will run succesive timesteps until timer.done() is true. It will - /// modify the reservoir and well states. - /// \param[in,out] timer governs the requested reporting timesteps - /// \param[in,out] state state of reservoir: pressure, fluxes - /// \param[in,out] well_state state of wells: bhp, perforation rates - /// \return simulation report, with timing data - SimulatorReport run(SimulatorTimer& timer, - TwophaseState& state, - std::vector& src, - WellState& well_state); - - private: - class Impl; - // Using shared_ptr instead of scoped_ptr since scoped_ptr requires complete type for Impl. - boost::shared_ptr pimpl_; - }; - -} // namespace Opm - -#endif // OPM_SIMULATORFULLYIMPLICITBLACKOIL_HEADER_INCLUDED