diff --git a/opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver.hpp b/opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver.hpp
new file mode 100644
index 000000000..b938a3922
--- /dev/null
+++ b/opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver.hpp
@@ -0,0 +1,373 @@
+/*
+ Copyright 2013 SINTEF ICT, Applied Mathematics.
+ Copyright 2014 STATOIL ASA.
+
+
+ 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_FULLYIMPLICITBLACKOILPOLYMERSOLVER_HEADER_INCLUDED
+#define OPM_FULLYIMPLICITBLACKOILPOLYMERSOLVER_HEADER_INCLUDED
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+struct UnstructuredGrid;
+struct Wells;
+
+namespace Opm {
+
+ namespace parameter { class ParameterGroup; }
+ class DerivedGeology;
+ class RockCompressibility;
+ class NewtonIterationBlackoilInterface;
+ class PolymerBlackoilState;
+ class WellStateFullyImplicitBlackoil;
+
+
+ /// A fully implicit solver for the black-oil-polymer problem.
+ ///
+ /// The simulator is capable of handling three-phase problems
+ /// where gas can be dissolved in oil (but not vice versa). It
+ /// uses an industry-standard TPFA discretization with per-phase
+ /// upwind weighting of mobilities.
+ ///
+ /// It uses automatic differentiation via the class AutoDiffBlock
+ /// to simplify assembly of the jacobian matrix.
+ template
+ class FullyImplicitBlackoilPolymerSolver
+ {
+ public:
+ /// \brief The type of the grid that we use.
+ typedef T Grid;
+ /// Construct a solver. It will retain references to the
+ /// arguments of this functions, and they are expected to
+ /// remain in scope for the lifetime of the solver.
+ /// \param[in] param parameters
+ /// \param[in] grid grid data structure
+ /// \param[in] fluid fluid properties
+ /// \param[in] geo rock properties
+ /// \param[in] rock_comp_props if non-null, rock compressibility properties
+ /// \param[in] wells well structure
+ /// \param[in] linsolver linear solver
+ FullyImplicitBlackoilSolver(const parameter::ParameterGroup& param,
+ const Grid& grid ,
+ const BlackoilPropsAdInterface& fluid,
+ const DerivedGeology& geo ,
+ const RockCompressibility* rock_comp_props,
+ const PolymerPropsAd& polymer_props_ad,
+ const Wells& wells,
+ const NewtonIterationBlackoilInterface& linsolver,
+ const bool has_disgas,
+ const bool has_vapoil,
+ const bool has_polymer);
+
+ /// \brief Set threshold pressures that prevent or reduce flow.
+ /// This prevents flow across faces if the potential
+ /// difference is less than the threshold. If the potential
+ /// difference is greater, the threshold value is subtracted
+ /// before calculating flow. This is treated symmetrically, so
+ /// flow is prevented or reduced in both directions equally.
+ /// \param[in] threshold_pressures_by_face array of size equal to the number of faces
+ /// of the grid passed in the constructor.
+ void setThresholdPressures(const std::vector& threshold_pressures_by_face);
+
+ /// Take a single forward step, modifiying
+ /// state.pressure()
+ /// state.faceflux()
+ /// state.saturation()
+ /// state.gasoilratio()
+ /// wstate.bhp()
+ /// \param[in] dt time step size
+ /// \param[in] state reservoir state
+ /// \param[in] wstate well state
+ void
+ step(const double dt ,
+ PolymerBlackoilState& state ,
+ WellStateFullyImplicitBlackoil& wstate);
+
+ private:
+ // Types and enums
+ typedef AutoDiffBlock ADB;
+ typedef ADB::V V;
+ typedef ADB::M M;
+ typedef Eigen::Array DataBlock;
+
+ struct ReservoirResidualQuant {
+ ReservoirResidualQuant();
+ std::vector accum; // Accumulations
+ ADB mflux; // Mass flux (surface conditions)
+ ADB b; // Reciprocal FVF
+ ADB head; // Pressure drop across int. interfaces
+ ADB mob; // Phase mobility (per cell)
+ };
+
+ struct SolutionState {
+ SolutionState(const int np);
+ ADB pressure;
+ std::vector saturation;
+ ADB rs;
+ ADB rv;
+ ADB concentration;
+ ADB qs;
+ ADB bhp;
+ };
+
+ 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 };
+
+ // the Newton relaxation type
+ enum RelaxType { DAMPEN, SOR };
+ enum PrimalVariables { Sg = 0, RS = 1, RV = 2 };
+
+ // Member data
+ const Grid& grid_;
+ const BlackoilPropsAdInterface& fluid_;
+ const DerivedGeology& geo_;
+ const RockCompressibility* rock_comp_props_;
+ const PolymerPropsAd& polymer_props_ad_;
+ const Wells& wells_;
+ const NewtonIterationBlackoilInterface& linsolver_;
+ // For each canonical phase -> true if active
+ const std::vector active_;
+ // Size = # active phases. Maps active -> canonical phase indices.
+ const std::vector canph_;
+ const std::vector cells_; // All grid cells
+ HelperOps ops_;
+ const WellOps wops_;
+ V cmax_;
+ const bool has_disgas_;
+ const bool has_vapoil_;
+ const bool has_polymer_;
+ const int poly_pos_;
+ double dp_max_rel_;
+ double ds_max_;
+ double drs_max_rel_;
+ enum RelaxType relax_type_;
+ double relax_max_;
+ double relax_increment_;
+ double relax_rel_tol_;
+ int max_iter_;
+ bool use_threshold_pressure_;
+ V threshold_pressures_by_interior_face_;
+
+ std::vector rq_;
+ std::vector phaseCondition_;
+ V well_perforation_pressure_diffs_; // Diff to bhp for each well perforation.
+
+ LinearisedBlackoilResidual residual_;
+
+ std::vector primalVariable_;
+
+ // Private methods.
+ SolutionState
+ constantState(const PolymerBlackoilState& x,
+ const WellStateFullyImplicitBlackoil& xw);
+
+ SolutionState
+ variableState(const PolymerBlackoilState& x,
+ const WellStateFullyImplicitBlackoil& xw);
+
+ void
+ computeAccum(const SolutionState& state,
+ const int aix );
+
+ void computeWellConnectionPressures(const SolutionState& state,
+ const WellStateFullyImplicitBlackoil& xw);
+
+ void
+ addWellControlEq(const SolutionState& state,
+ const WellStateFullyImplicitBlackoil& xw,
+ const V& aliveWells);
+
+ void
+ addWellEq(const SolutionState& state,
+ WellStateFullyImplicitBlackoil& xw,
+ V& aliveWells);
+
+ void updateWellControls(ADB& bhp,
+ ADB& well_phase_flow_rate,
+ WellStateFullyImplicitBlackoil& xw) const;
+
+ void
+ assemble(const V& dtpv,
+ const PolymerBlackoilState& x,
+ WellStateFullyImplicitBlackoil& xw);
+
+ V solveJacobianSystem() const;
+
+ void updateState(const V& dx,
+ PolymerBlackoilState& state,
+ WellStateFullyImplicitBlackoil& well_state);
+
+ std::vector
+ computePressures(const SolutionState& state) const;
+
+ std::vector
+ computeRelPerm(const SolutionState& state) const;
+
+ std::vector
+ computeRelPermWells(const SolutionState& state,
+ const DataBlock& well_s,
+ const std::vector& well_cells) const;
+
+ void
+ computeMassFlux(const int actph ,
+ const V& transi,
+ const ADB& kr ,
+ const ADB& p ,
+ const SolutionState& state );
+
+ void
+ computeMassFlux(const V& trans,
+ const ADB& mc,
+ const ADB& kro,
+ const ADB& krw_eff,
+ const ADB& krg,
+ const SolutionState& state);
+
+ void
+ computeCmax(PolymerBlackoilState& state,
+ const ADB& c);
+
+ ADB
+ computeMc(const SolutionState& state) const;
+
+ ADB
+ rockPorosity(const ADB& p) const;
+
+ ADB
+ rockPermeability(const ADB& p) const;
+ void applyThresholdPressures(ADB& dp);
+
+ double
+ residualNorm() const;
+
+ std::vector residuals() const;
+
+ ADB
+ fluidViscosity(const int phase,
+ const ADB& p ,
+ const ADB& rs ,
+ const ADB& rv ,
+ const std::vector& cond,
+ const std::vector& cells) const;
+
+ ADB
+ fluidReciprocFVF(const int phase,
+ const ADB& p ,
+ const ADB& rs ,
+ const ADB& rv ,
+ const std::vector& cond,
+ const std::vector& cells) const;
+
+ ADB
+ fluidDensity(const int phase,
+ const ADB& p ,
+ const ADB& rs ,
+ const ADB& rv ,
+ const std::vector& cond,
+ const std::vector& cells) const;
+
+ V
+ fluidRsSat(const V& p,
+ const V& so,
+ const std::vector& cells) const;
+
+ ADB
+ fluidRsSat(const ADB& p,
+ const ADB& so,
+ const std::vector& cells) const;
+
+ V
+ fluidRvSat(const V& p,
+ const V& so,
+ const std::vector& cells) const;
+
+ ADB
+ fluidRvSat(const ADB& p,
+ const ADB& so,
+ const std::vector& cells) const;
+
+ ADB
+ computeMc(const SolutionState& state) const;
+
+ ADB
+ poroMult(const ADB& p) const;
+
+ ADB
+ transMult(const ADB& p) const;
+
+ void
+ classifyCondition(const SolutionState& state,
+ std::vector& cond ) const;
+
+ const std::vector
+ phaseCondition() const {return phaseCondition_;}
+
+ void
+ classifyCondition(const PolymerBlackoilState& state);
+
+
+ /// update the primal variable for Sg, Rv or Rs. The Gas phase must
+ /// be active to call this method.
+ void
+ updatePrimalVariableFromState(const PolymerBlackoilState& state);
+
+ /// Update the phaseCondition_ member based on the primalVariable_ member.
+ void
+ updatePhaseCondFromPrimalVariable();
+
+ /// Compute convergence based on total mass balance (tol_mb) and maximum
+ /// residual mass balance (tol_cnv).
+ bool getConvergence(const double dt);
+
+ void detectNewtonOscillations(const std::vector>& residual_history,
+ const int it, const double relaxRelTol,
+ bool& oscillate, bool& stagnate) const;
+
+ void stablizeNewton(V& dx, V& dxOld, const double omega, const RelaxType relax_type) const;
+
+ double dpMaxRel() const { return dp_max_rel_; }
+ double dsMax() const { return ds_max_; }
+ double drsMaxRel() const { return drs_max_rel_; }
+ enum RelaxType relaxType() const { return relax_type_; }
+ double relaxMax() const { return relax_max_; };
+ double relaxIncrement() const { return relax_increment_; };
+ double relaxRelTol() const { return relax_rel_tol_; };
+ double maxIter() const { return max_iter_; }
+
+ };
+} // namespace Opm
+
+#include "FullyImplicitBlackoilSolver_impl.hpp"
+
+#endif // OPM_FULLYIMPLICITBLACKOILSOLVER_HEADER_INCLUDED
diff --git a/opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver_impl.hpp b/opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver_impl.hpp
new file mode 100644
index 000000000..365c937bb
--- /dev/null
+++ b/opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver_impl.hpp
@@ -0,0 +1,2333 @@
+/*
+ Copyright 2013 SINTEF ICT, Applied Mathematics.
+ Copyright 2014 STATOIL ASA.
+
+ 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
+
+// A debugging utility.
+#define DUMP(foo) \
+ do { \
+ std::cout << "==========================================\n" \
+ << #foo ":\n" \
+ << collapseJacs(foo) << std::endl; \
+ } while (0)
+
+#define DUMPVAL(foo) \
+ do { \
+ std::cout << "==========================================\n" \
+ << #foo ":\n" \
+ << foo.value() << std::endl; \
+ } while (0)
+
+#define DISKVAL(foo) \
+ do { \
+ std::ofstream os(#foo); \
+ os.precision(16); \
+ os << foo.value() << std::endl; \
+ } while (0)
+
+
+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;
+ }
+
+
+
+ template
+ V computePerfPress(const Grid& grid, const Wells& wells, const V& rho, const double grav)
+ {
+ using namespace Opm::AutoDiffGrid;
+ const int nw = wells.number_of_wells;
+ const int nperf = wells.well_connpos[nw];
+ const int dim = dimensions(grid);
+ V wdp = V::Zero(nperf,1);
+ assert(wdp.size() == rho.size());
+
+ // Main loop, iterate over all perforations,
+ // using the following formula:
+ // wdp(perf) = g*(perf_z - well_ref_z)*rho(perf)
+ // where the total density rho(perf) is taken to be
+ // sum_p (rho_p*saturation_p) in the perforation cell.
+ // [although this is computed on the outside of this function].
+ 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[w + 1]; ++j) {
+ const int cell = wells.well_cells[j];
+ const double cell_depth = cellCentroid(grid, cell)[dim - 1];
+ wdp[j] = rho[j]*grav*(cell_depth - ref_depth);
+ }
+ }
+ return wdp;
+ }
+
+
+
+ template
+ std::vector
+ activePhases(const PU& pu)
+ {
+ const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
+ std::vector active(maxnp, false);
+
+ for (int p = 0; p < pu.MaxNumPhases; ++p) {
+ active[ p ] = pu.phase_used[ p ] != 0;
+ }
+
+ return active;
+ }
+
+
+
+ template
+ std::vector
+ active2Canonical(const PU& pu)
+ {
+ const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
+ std::vector act2can(maxnp, -1);
+
+ for (int phase = 0; phase < maxnp; ++phase) {
+ if (pu.phase_used[ phase ]) {
+ act2can[ pu.phase_pos[ phase ] ] = phase;
+ }
+ }
+
+ return act2can;
+ }
+
+
+ template
+ int polymerPos(const PU& pu)
+ {
+ const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
+ int pos = 0;
+ for (int p = 0; p < maxnp; ++p) {
+ pos += pu.phase_used[p];
+ }
+
+ return pos;
+ }
+} // Anonymous namespace
+
+
+
+ template
+ FullyImplicitBlackoilPolymerSolver::
+ FullyImplicitBlackoilPolymerSolver(const parameter::ParameterGroup& param,
+ const Grid& grid ,
+ const BlackoilPropsAdInterface& fluid,
+ const DerivedGeology& geo ,
+ const RockCompressibility* rock_comp_props,
+ const PolymerPropsAd& polymer_props_ad,
+ const Wells& wells,
+ const NewtonIterationBlackoilInterface& linsolver,
+ const bool has_disgas,
+ const bool has_vapoil,
+ const bool has_polymer)
+ : grid_ (grid)
+ , fluid_ (fluid)
+ , geo_ (geo)
+ , rock_comp_props_(rock_comp_props)
+ , polymer_props_ad_(polymer_props_ad)
+ , wells_ (wells)
+ , linsolver_ (linsolver)
+ , active_(activePhases(fluid.phaseUsage()))
+ , canph_ (active2Canonical(fluid.phaseUsage()))
+ , cells_ (buildAllCells(Opm::AutoDiffGrid::numCells(grid)))
+ , ops_ (grid)
+ , wops_ (wells)
+ , cmax_(V::Zero(numCells(grid)))
+ , has_disgas_(has_disgas)
+ , has_vapoil_(has_vapoil)
+ , has_polymer_(has_polymer)
+ , poly_pos_(polymerPos(fluid.phaseUsage()))
+ , dp_max_rel_ (1.0e9)
+ , ds_max_ (0.2)
+ , drs_max_rel_ (1.0e9)
+ , relax_type_ (DAMPEN)
+ , relax_max_ (0.5)
+ , relax_increment_ (0.1)
+ , relax_rel_tol_ (0.2)
+ , max_iter_ (15)
+ , use_threshold_pressure_(false)
+ , rq_ (fluid.numPhases() + 1)
+ , phaseCondition_(AutoDiffGrid::numCells(grid))
+ , residual_ ( { std::vector(fluid.numPhases() + 1, ADB::null()),
+ ADB::null(),
+ ADB::null() } )
+ {
+ if (has_polymer_) {
+ if (!active_[Water]) {
+ OPM_THROW(std::logic_error, "Polymer must solved in water!\n");
+ }
+ }
+ dp_max_rel_ = param.getDefault("dp_max_rel", dp_max_rel_);
+ ds_max_ = param.getDefault("ds_max", ds_max_);
+ drs_max_rel_ = param.getDefault("drs_max_rel", drs_max_rel_);
+ relax_max_ = param.getDefault("relax_max", relax_max_);
+ max_iter_ = param.getDefault("max_iter", max_iter_);
+
+ std::string relaxation_type = param.getDefault("relax_type", std::string("dampen"));
+ if (relaxation_type == "dampen") {
+ relax_type_ = DAMPEN;
+ } else if (relaxation_type == "sor") {
+ relax_type_ = SOR;
+ } else {
+ OPM_THROW(std::runtime_error, "Unknown Relaxtion Type " << relaxation_type);
+ }
+ }
+
+
+
+ template
+ void
+ FullyImplicitBlackoilPolymerSolver::
+ setThresholdPressures(const std::vector& threshold_pressures_by_face)
+ {
+ const int num_faces = AutoDiffGrid::numFaces(grid_);
+ if (int(threshold_pressures_by_face.size()) != num_faces) {
+ OPM_THROW(std::runtime_error, "Illegal size of threshold_pressures_by_face input, must be equal to number of faces.");
+ }
+ use_threshold_pressure_ = true;
+ // Map to interior faces.
+ const int num_ifaces = ops_.internal_faces.size();
+ threshold_pressures_by_interior_face_.resize(num_ifaces);
+ for (int ii = 0; ii < num_ifaces; ++ii) {
+ threshold_pressures_by_interior_face_[ii] = threshold_pressures_by_face[ops_.internal_faces[ii]];
+ }
+ }
+
+
+
+
+ template
+ void
+ FullyImplicitBlackoilPolymerSolver::
+ step(const double dt,
+ PolymerBlackoilState& x ,
+ WellStateFullyImplicitBlackoil& xw)
+ {
+ const V pvdt = geo_.poreVolume() / dt;
+
+ if (active_[Gas]) { updatePrimalVariableFromState(x); }
+
+ {
+ const SolutionState state = constantState(x, xw);
+ computeAccum(state, 0);
+ computeWellConnectionPressures(state, xw);
+ }
+
+
+ std::vector> residual_history;
+
+ assemble(pvdt, x, xw);
+
+
+ bool converged = false;
+ double omega = 1.;
+ const double r0 = residualNorm();
+
+ residual_history.push_back(residuals());
+
+ converged = getConvergence(dt);
+
+ int it = 0;
+ std::cout << "\nIteration Residual\n"
+ << std::setw(9) << it << std::setprecision(9)
+ << std::setw(18) << r0 << std::endl;
+
+ const int sizeNonLinear = residual_.sizeNonLinear();
+
+ V dxOld = V::Zero(sizeNonLinear);
+
+ bool isOscillate = false;
+ bool isStagnate = false;
+ const enum RelaxType relaxtype = relaxType();
+
+ while ((!converged) && (it < maxIter())) {
+ V dx = solveJacobianSystem();
+
+ detectNewtonOscillations(residual_history, it, relaxRelTol(), isOscillate, isStagnate);
+
+ if (isOscillate) {
+ omega -= relaxIncrement();
+ omega = std::max(omega, relaxMax());
+ std::cout << " Oscillating behavior detected: Relaxation set to " << omega << std::endl;
+ }
+
+ stablizeNewton(dx, dxOld, omega, relaxtype);
+
+ updateState(dx, x, xw);
+
+ assemble(pvdt, x, xw);
+
+ const double r = residualNorm();
+
+ residual_history.push_back(residuals());
+
+ converged = getConvergence(dt);
+
+ it += 1;
+ std::cout << std::setw(9) << it << std::setprecision(9)
+ << std::setw(18) << r << std::endl;
+ }
+
+ if (!converged) {
+ 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.");
+ }
+ }
+
+
+
+
+
+ template
+ FullyImplicitBlackoilPolymerSolver::ReservoirResidualQuant::ReservoirResidualQuant()
+ : accum(2, ADB::null())
+ , mflux( ADB::null())
+ , b ( ADB::null())
+ , head ( ADB::null())
+ , mob ( ADB::null())
+ {
+ }
+
+
+
+
+
+ template
+ FullyImplicitBlackoilPolymerSolver::SolutionState::SolutionState(const int np)
+ : pressure ( ADB::null())
+ , saturation(np, ADB::null())
+ , rs ( ADB::null())
+ , rv ( ADB::null())
+ , concentration( ADB::null())
+ , qs ( ADB::null())
+ , bhp ( ADB::null())
+ {
+ }
+
+
+
+
+
+ template
+ FullyImplicitBlackoilPolymerSolver::
+ 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());
+ }
+
+
+
+
+
+ template
+ typename FullyImplicitBlackoilPolymerSolver::SolutionState
+ FullyImplicitBlackoilPolymerSolver::constantState(const PolymerBlackoilState& x,
+ const WellStateFullyImplicitBlackoil& xw)
+ {
+ auto state = variableState(x, xw);
+
+ // HACK: throw away the derivatives. this may not be the most
+ // performant way to do things, but it will make the state
+ // automatically consistent with variableState() (and doing
+ // things automatically is all the rage in this module ;)
+ state.pressure = ADB::constant(state.pressure.value());
+ state.rs = ADB::constant(state.rs.value());
+ state.rv = ADB::constant(state.rv.value());
+ state.concentration = ADB::constant(state.concentration.value());
+ for (int phaseIdx= 0; phaseIdx < x.numPhases(); ++ phaseIdx)
+ state.saturation[phaseIdx] = ADB::constant(state.saturation[phaseIdx].value());
+ state.qs = ADB::constant(state.qs.value());
+ state.bhp = ADB::constant(state.bhp.value());
+
+ return state;
+ }
+
+
+
+
+
+ template
+ typename FullyImplicitBlackoilPolymerSolver::SolutionState
+ FullyImplicitBlackoilPolymerSolver::variableState(const PolymerBlackoilState& x,
+ const WellStateFullyImplicitBlackoil& xw)
+ {
+ using namespace Opm::AutoDiffGrid;
+ const int nc = numCells(grid_);
+ const int np = x.numPhases();
+
+ std::vector vars0;
+ // p, Sw and Rs, Rv or Sg, concentration are used as primary depending on solution conditions
+ vars0.reserve(np + 2);
+ // Initial pressure.
+ assert (not x.pressure().empty());
+ const V p = Eigen::Map(& x.pressure()[0], nc, 1);
+ vars0.push_back(p);
+
+ // Initial saturation.
+ assert (not x.saturation().empty());
+ const DataBlock s = Eigen::Map(& x.saturation()[0], nc, np);
+ const Opm::PhaseUsage pu = fluid_.phaseUsage();
+ // We do not handle a Water/Gas situation correctly, guard against it.
+ assert (active_[ Oil]);
+ if (active_[ Water ]) {
+ const V sw = s.col(pu.phase_pos[ Water ]);
+ vars0.push_back(sw);
+ }
+
+ // store cell status in vectors
+ V isRs = V::Zero(nc,1);
+ V isRv = V::Zero(nc,1);
+ V isSg = V::Zero(nc,1);
+
+ if (active_[ Gas ]){
+ for (int c = 0; c < nc ; c++ ) {
+ switch (primalVariable_[c]) {
+ case PrimalVariables::RS:
+ isRs[c] = 1;
+ break;
+
+ case PrimalVariables::RV:
+ isRv[c] = 1;
+ break;
+
+ default:
+ isSg[c] = 1;
+ break;
+ }
+ }
+
+
+ // define new primary variable xvar depending on solution condition
+ V xvar(nc);
+ const V sg = s.col(pu.phase_pos[ Gas ]);
+ const V rs = Eigen::Map(& x.gasoilratio()[0], x.gasoilratio().size());
+ const V rv = Eigen::Map(& x.rv()[0], x.rv().size());
+ xvar = isRs*rs + isRv*rv + isSg*sg;
+ vars0.push_back(xvar);
+ }
+
+ // Initial polymer concentration.
+ if (has_polymer_) {
+ assert (not x.concentration().empty());
+ const V c = Eigen::Map(& x.concentration()[0], nc, 1);
+ vars0.push_back(c);
+ }
+
+ // Initial well rates.
+ assert (not xw.wellRates().empty());
+ // Need to reshuffle well rates, from phase running fastest
+ // to wells running fastest.
+ const int nw = wells_.number_of_wells;
+ // The transpose() below switches the ordering.
+ const DataBlock wrates = Eigen::Map(& xw.wellRates()[0], nw, np).transpose();
+ const V qs = Eigen::Map(wrates.data(), nw*np);
+ vars0.push_back(qs);
+
+ // Initial well 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++ ];
+
+ // Saturations
+ const std::vector& bpat = vars[0].blockPattern();
+ {
+ ADB so = ADB::constant(V::Ones(nc, 1), bpat);
+
+ if (active_[ Water ]) {
+ ADB& sw = vars[ nextvar++ ];
+ state.saturation[pu.phase_pos[ Water ]] = sw;
+ so = so - sw;
+ }
+
+ if (active_[ Gas]) {
+ // Define Sg Rs and Rv in terms of xvar.
+ const ADB& xvar = vars[ nextvar++ ];
+ const ADB& sg = isSg*xvar + isRv* so;
+ state.saturation[ pu.phase_pos[ Gas ] ] = sg;
+ so = so - sg;
+ const ADB rsSat = fluidRsSat(state.pressure, so, cells_);
+ const ADB rvSat = fluidRvSat(state.pressure, so, cells_);
+
+ if (has_disgas_) {
+ state.rs = (1-isRs) * rsSat + isRs*xvar;
+ } else {
+ state.rs = rsSat;
+ }
+ if (has_vapoil_) {
+ state.rv = (1-isRv) * rvSat + isRv*xvar;
+ } else {
+ state.rv = rvSat;
+ }
+ }
+
+ if (active_[ Oil ]) {
+ // Note that so is never a primary variable.
+ state.saturation[ pu.phase_pos[ Oil ] ] = so;
+ }
+ }
+
+ // Concentration.
+ if (has_polymer_) {
+ state.concentration = vars[nextvar++];
+ }
+ // Qs.
+ state.qs = vars[ nextvar++ ];
+
+ // Bhp.
+ state.bhp = vars[ nextvar++ ];
+
+ assert(nextvar == int(vars.size()));
+
+ return state;
+ }
+
+
+
+
+
+ template
+ void
+ FullyImplicitBlackoilPolymerSolver::computeAccum(const SolutionState& state,
+ const int aix )
+ {
+ const Opm::PhaseUsage& pu = fluid_.phaseUsage();
+
+ const ADB& press = state.pressure;
+ const std::vector& sat = state.saturation;
+ const ADB& rs = state.rs;
+ const ADB& rv = state.rv;
+ const ADB& c = state.concentration;
+
+ const std::vector cond = phaseCondition();
+ std::vector pressure = computePressures(state);
+
+ const ADB pv_mult = poroMult(press);
+
+ const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
+ for (int phase = 0; phase < maxnp; ++phase) {
+ if (active_[ phase ]) {
+ const int pos = pu.phase_pos[ phase ];
+ rq_[pos].b = fluidReciprocFVF(phase, pressure, rs, rv, cond, cells_);
+ rq_[pos].accum[aix] = pv_mult * rq_[pos].b * sat[pos];
+ // DUMP(rq_[pos].b);
+ // DUMP(rq_[pos].accum[aix]);
+ }
+ }
+
+ if (active_[ Oil ] && active_[ Gas ]) {
+ // Account for gas dissolved in oil and vaporized oil
+ const int po = pu.phase_pos[ Oil ];
+ const int pg = pu.phase_pos[ Gas ];
+
+ rq_[pg].accum[aix] += state.rs * rq_[po].accum[aix];
+ rq_[po].accum[aix] += state.rv * rq_[pg].accum[aix];
+ //DUMP(rq_[pg].accum[aix]);
+ }
+ if (has_polymer_) {
+ // compute polymer properties.
+ const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern());
+ const ADB ads = polymer_props_ad_.adsorption(state.concentration, cmax);
+ const double rho_rock = polymer_props_ad_.rockDensity();
+ const V phi = Eigen::Map(& fluid_.porosity()[0], numCells(grid_), 1);
+ const double dead_pore_vol = polymer_props_ad_.deadPoreVol();
+ // compute total phases and determin polymer position.
+ rq_[poly_pos_].accum[aix] = pv_mult * rq_[pu.phase_pos[Water]].b * sat[pu.phase_pos[Water]] * c * (1. - dead_pore_vol) + pv_mult * rho_rock * (1. - phi) / phi * ads;
+ }
+
+ }
+
+
+
+
+
+ template
+ void FullyImplicitBlackoilPolymerSovler::computeCmax(PolymerBlackoilState& state,
+ const ADB& c)
+ {
+ const int nc = numCells(grid_);
+ for (int i = 0; i < nc; ++i) {
+ cmax_(i) = std::max(cmax_(i), c.value()(i));
+ }
+ std::copy(&cmax_[0], &cmax_[0] + nc, state.maxconcentration().begin());
+ }
+
+
+
+
+
+ template
+ void FullyImplicitBlackoilPolymerSolver::computeWellConnectionPressures(const SolutionState& state,
+ const WellStateFullyImplicitBlackoil& xw)
+ {
+ using namespace Opm::AutoDiffGrid;
+ // 1. Compute properties required by computeConnectionPressureDelta().
+ // Note that some of the complexity of this part is due to the function
+ // taking std::vector arguments, and not Eigen objects.
+ const int nperf = wells_.well_connpos[wells_.number_of_wells];
+ const std::vector well_cells(wells_.well_cells, wells_.well_cells + nperf);
+ // Compute b, rsmax, rvmax values for perforations.
+ const ADB perf_press = subset(state.pressure, well_cells);
+ std::vector perf_cond(nperf);
+ const std::vector& pc = phaseCondition();
+ for (int perf = 0; perf < nperf; ++perf) {
+ perf_cond[perf] = pc[well_cells[perf]];
+ }
+ const PhaseUsage& pu = fluid_.phaseUsage();
+ DataBlock b(nperf, pu.num_phases);
+ std::vector rssat_perf(nperf, 0.0);
+ std::vector rvsat_perf(nperf, 0.0);
+ if (pu.phase_used[BlackoilPhases::Aqua]) {
+ const ADB bw = fluid_.bWat(perf_press, well_cells);
+ b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw.value();
+ }
+ assert(active_[Oil]);
+ const ADB perf_so = subset(state.saturation[pu.phase_pos[Oil]], well_cells);
+ if (pu.phase_used[BlackoilPhases::Liquid]) {
+ const ADB perf_rs = subset(state.rs, well_cells);
+ const ADB bo = fluid_.bOil(perf_press, perf_rs, perf_cond, well_cells);
+ b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo.value();
+ const V rssat = fluidRsSat(perf_press.value(), perf_so.value(), well_cells);
+ rssat_perf.assign(rssat.data(), rssat.data() + nperf);
+ }
+ if (pu.phase_used[BlackoilPhases::Vapour]) {
+ const ADB perf_rv = subset(state.rv, well_cells);
+ const ADB bg = fluid_.bGas(perf_press, perf_rv, perf_cond, well_cells);
+ b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg.value();
+ const V rvsat = fluidRvSat(perf_press.value(), perf_so.value(), well_cells);
+ rvsat_perf.assign(rvsat.data(), rvsat.data() + nperf);
+ }
+ // b is row major, so can just copy data.
+ std::vector b_perf(b.data(), b.data() + nperf * pu.num_phases);
+ // Extract well connection depths.
+ const V depth = cellCentroidsZToEigen(grid_);
+ const V pdepth = subset(depth, well_cells);
+ std::vector perf_depth(pdepth.data(), pdepth.data() + nperf);
+ // Surface density.
+ std::vector surf_dens(fluid_.surfaceDensity(), fluid_.surfaceDensity() + pu.num_phases);
+ // Gravity
+ double grav = 0.0;
+ const double* g = geo_.gravity();
+ const int dim = dimensions(grid_);
+ if (g) {
+ // Guard against gravity in anything but last dimension.
+ for (int dd = 0; dd < dim - 1; ++dd) {
+ assert(g[dd] == 0.0);
+ }
+ grav = g[dim - 1];
+ }
+
+ // 2. Compute pressure deltas, and store the results.
+ std::vector cdp = WellDensitySegmented
+ ::computeConnectionPressureDelta(wells_, xw, fluid_.phaseUsage(),
+ b_perf, rssat_perf, rvsat_perf, perf_depth,
+ surf_dens, grav);
+ well_perforation_pressure_diffs_ = Eigen::Map(cdp.data(), nperf);
+ }
+
+
+
+
+
+ template
+ void
+ FullyImplicitBlackoilPolymerSolver::
+ assemble(const V& pvdt,
+ const BlackoilState& x ,
+ WellStateFullyImplicitBlackoil& xw )
+ {
+ using namespace Opm::AutoDiffGrid;
+ // Create the primary variables.
+ SolutionState state = variableState(x, xw);
+
+ // DISKVAL(state.pressure);
+ // DISKVAL(state.saturation[0]);
+ // DISKVAL(state.saturation[1]);
+ // DISKVAL(state.saturation[2]);
+ // DISKVAL(state.rs);
+ // DISKVAL(state.rv);
+ // DISKVAL(state.qs);
+ // DISKVAL(state.bhp);
+
+ // -------- Mass balance equations --------
+
+ // Compute b_p and the accumulation term b_p*s_p for each phase,
+ // except gas. For gas, we compute b_g*s_g + Rs*b_o*s_o.
+ // These quantities are stored in rq_[phase].accum[1].
+ // The corresponding accumulation terms from the start of
+ // the timestep (b^0_p*s^0_p etc.) were already computed
+ // in step() and stored in rq_[phase].accum[0].
+ computeAccum(state, 1);
+
+ // Set up the common parts of the mass balance equations
+ // for each active phase.
+ const V transi = subset(geo_.transmissibility(), ops_.internal_faces);
+ const std::vector kr = computeRelPerm(state);
+ const std::vector pressures = computePressures(state);
+ computeMassFlux(transi, kr, pressures, state);
+ for (int phaseIdx = 0; phaseIdx < fluid_.numPhases(); ++phaseIdx) {
+ // std::cout << "===== kr[" << phase << "] = \n" << std::endl;
+ // std::cout << kr[phase];
+ // std::cout << "===== rq_[" << phase << "].mflux = \n" << std::endl;
+ // std::cout << rq_[phase].mflux;
+ residual_.material_balance_eq[ phaseIdx ] =
+ pvdt*(rq_[phaseIdx].accum[1] - rq_[phaseIdx].accum[0])
+ + ops_.div*rq_[phaseIdx].mflux;
+
+
+ // DUMP(ops_.div*rq_[phase].mflux);
+ // DUMP(residual_.material_balance_eq[phase]);
+ }
+
+ // -------- Extra (optional) rs and rv contributions to the mass balance equations --------
+
+ // Add the extra (flux) terms to the mass balance equations
+ // From gas dissolved in the oil phase (rs) and oil vaporized in the gas phase (rv)
+ // The extra terms in the accumulation part of the equation are already handled.
+ if (active_[ Oil ] && active_[ Gas ]) {
+ const int po = fluid_.phaseUsage().phase_pos[ Oil ];
+ const int pg = fluid_.phaseUsage().phase_pos[ Gas ];
+
+ const UpwindSelector upwindOil(grid_, ops_,
+ rq_[po].head.value());
+ const ADB rs_face = upwindOil.select(state.rs);
+
+ const UpwindSelector upwindGas(grid_, ops_,
+ rq_[pg].head.value());
+ const ADB rv_face = upwindGas.select(state.rv);
+
+ residual_.material_balance_eq[ pg ] += ops_.div * (rs_face * rq_[po].mflux);
+ residual_.material_balance_eq[ po ] += ops_.div * (rv_face * rq_[pg].mflux);
+
+ // DUMP(residual_.material_balance_eq[ Gas ]);
+
+ }
+
+ const Opm::PhaseUsage& pu = fluid_.phaseUsage();
+ // Add polymer equation.
+ if (has_polymer) {
+ residula_.material_balance_eq[poly_pos_] = pvdt*(rq_[poly_pos_].accum[1] - rq_[poly_pos_].accum[0])
+ + ops_.div*rq_[poly_pos_].mflux;
+ }
+
+ // Note: updateWellControls() can change all its arguments if
+ // a well control is switched.
+ updateWellControls(state.bhp, state.qs, xw);
+ V aliveWells;
+ addWellEq(state, xw, aliveWells);
+ addWellControlEq(state, xw, aliveWells);
+ }
+
+
+
+
+
+ template
+ void FullyImplicitBlackoilPolymerSolver::addWellEq(const SolutionState& state,
+ WellStateFullyImplicitBlackoil& xw,
+ V& aliveWells)
+ {
+ const int nc = Opm::AutoDiffGrid::numCells(grid_);
+ const int np = wells_.number_of_phases;
+ const int nw = wells_.number_of_wells;
+ const int nperf = wells_.well_connpos[nw];
+ const Opm::PhaseUsage& pu = fluid_.phaseUsage();
+ V Tw = Eigen::Map(wells_.WI, nperf);
+ const std::vector well_cells(wells_.well_cells, wells_.well_cells + nperf);
+
+ // pressure diffs computed already (once per step, not changing per iteration)
+ const V& cdp = well_perforation_pressure_diffs_;
+
+ // Extract variables for perforation cell pressures
+ // and corresponding perforation well pressures.
+ const ADB p_perfcell = subset(state.pressure, well_cells);
+
+ // DUMPVAL(p_perfcell);
+ // DUMPVAL(state.bhp);
+ // DUMPVAL(ADB::constant(cdp));
+
+ // Pressure drawdown (also used to determine direction of flow)
+ const ADB drawdown = p_perfcell - (wops_.w2p * state.bhp + cdp);
+
+ // current injecting connections
+ auto connInjInx = drawdown.value() < 0;
+
+ // injector == 1, producer == 0
+ V isInj = V::Zero(nw);
+ for (int w = 0; w < nw; ++w) {
+ if (wells_.type[w] == INJECTOR) {
+ isInj[w] = 1;
+ }
+ }
+
+// // A cross-flow connection is defined as a connection which has opposite
+// // flow-direction to the well total flow
+// V isInjPerf = (wops_.w2p * isInj);
+// auto crossFlowConns = (connInjInx != isInjPerf);
+
+// bool allowCrossFlow = true;
+
+// if (not allowCrossFlow) {
+// auto closedConns = crossFlowConns;
+// for (int c = 0; c < nperf; ++c) {
+// if (closedConns[c]) {
+// Tw[c] = 0;
+// }
+// }
+// connInjInx = !closedConns;
+// }
+// TODO: not allow for crossflow
+
+
+ V isInjInx = V::Zero(nperf);
+ V isNotInjInx = V::Zero(nperf);
+ for (int c = 0; c < nperf; ++c){
+ if (connInjInx[c])
+ isInjInx[c] = 1;
+ else
+ isNotInjInx[c] = 1;
+ }
+
+
+ // HANDLE FLOW INTO WELLBORE
+
+ // compute phase volumerates standard conditions
+ std::vector cq_ps(np, ADB::null());
+ for (int phase = 0; phase < np; ++phase) {
+ const ADB& wellcell_mob = subset ( rq_[phase].mob, well_cells);
+ const ADB cq_p = -(isNotInjInx * Tw) * (wellcell_mob * drawdown);
+ cq_ps[phase] = subset(rq_[phase].b,well_cells) * cq_p;
+ }
+ if (active_[Oil] && active_[Gas]) {
+ const int oilpos = pu.phase_pos[Oil];
+ const int gaspos = pu.phase_pos[Gas];
+ ADB cq_psOil = cq_ps[oilpos];
+ ADB cq_psGas = cq_ps[gaspos];
+ cq_ps[gaspos] += subset(state.rs,well_cells) * cq_psOil;
+ cq_ps[oilpos] += subset(state.rv,well_cells) * cq_psGas;
+ }
+
+ // phase rates at std. condtions
+ std::vector q_ps(np, ADB::null());
+ for (int phase = 0; phase < np; ++phase) {
+ q_ps[phase] = wops_.p2w * cq_ps[phase];
+ }
+
+ // total rates at std
+ ADB qt_s = ADB::constant(V::Zero(nw), state.bhp.blockPattern());
+ for (int phase = 0; phase < np; ++phase) {
+ qt_s += subset(state.qs, Span(nw, 1, phase*nw));
+ }
+
+ // compute avg. and total wellbore phase volumetric rates at std. conds
+ const DataBlock compi = Eigen::Map(wells_.comp_frac, nw, np);
+ std::vector wbq(np, ADB::null());
+ ADB wbqt = ADB::constant(V::Zero(nw), state.pressure.blockPattern());
+ for (int phase = 0; phase < np; ++phase) {
+ const int pos = pu.phase_pos[phase];
+ wbq[phase] = (isInj * compi.col(pos)) * qt_s - q_ps[phase];
+ wbqt += wbq[phase];
+ }
+ // DUMPVAL(wbqt);
+
+ // check for dead wells
+ aliveWells = V::Constant(nw, 1.0);
+ for (int w = 0; w < nw; ++w) {
+ if (wbqt.value()[w] == 0) {
+ aliveWells[w] = 0.0;
+ }
+ }
+ // compute wellbore mixture at std conds
+ Selector notDeadWells_selector(wbqt.value(), Selector::Zero);
+ std::vector mix_s(np, ADB::null());
+ for (int phase = 0; phase < np; ++phase) {
+ const int pos = pu.phase_pos[phase];
+ mix_s[phase] = notDeadWells_selector.select(ADB::constant(compi.col(pos)), wbq[phase]/wbqt);
+ }
+
+
+ // HANDLE FLOW OUT FROM WELLBORE
+
+ // Total mobilities
+ ADB mt = subset(rq_[0].mob,well_cells);
+ for (int phase = 1; phase < np; ++phase) {
+ mt += subset(rq_[phase].mob,well_cells);
+ }
+
+ // DUMPVAL(ADB::constant(isInjInx));
+ // DUMPVAL(ADB::constant(Tw));
+ // DUMPVAL(mt);
+ // DUMPVAL(drawdown);
+
+ // injection connections total volumerates
+ ADB cqt_i = -(isInjInx * Tw) * (mt * drawdown);
+
+ // compute volume ratio between connection at standard conditions
+ ADB volRat = ADB::constant(V::Zero(nperf), state.pressure.blockPattern());
+ std::vector cmix_s(np, ADB::null());
+ for (int phase = 0; phase < np; ++phase) {
+ cmix_s[phase] = wops_.w2p * mix_s[phase];
+ }
+
+ ADB well_rv = subset(state.rv,well_cells);
+ ADB well_rs = subset(state.rs,well_cells);
+ ADB d = V::Constant(nperf,1.0) - well_rv * well_rs;
+
+ for (int phase = 0; phase < np; ++phase) {
+ ADB tmp = cmix_s[phase];
+
+ if (phase == Oil && active_[Gas]) {
+ const int gaspos = pu.phase_pos[Gas];
+ tmp = tmp - subset(state.rv,well_cells) * cmix_s[gaspos] / d;
+ }
+ if (phase == Gas && active_[Oil]) {
+ const int oilpos = pu.phase_pos[Oil];
+ tmp = tmp - subset(state.rs,well_cells) * cmix_s[oilpos] / d;
+ }
+ volRat += tmp / subset(rq_[phase].b,well_cells);
+ }
+
+ // DUMPVAL(cqt_i);
+ // DUMPVAL(volRat);
+
+ // injecting connections total volumerates at std cond
+ ADB cqt_is = cqt_i/volRat;
+
+ // connection phase volumerates at std cond
+ std::vector cq_s(np, ADB::null());
+ for (int phase = 0; phase < np; ++phase) {
+ cq_s[phase] = cq_ps[phase] + (wops_.w2p * mix_s[phase])*cqt_is;
+ }
+
+ // DUMPVAL(mix_s[2]);
+ // DUMPVAL(cq_ps[2]);
+
+ // Add well contributions to mass balance equations
+ for (int phase = 0; phase < np; ++phase) {
+ residual_.material_balance_eq[phase] -= superset(cq_s[phase],well_cells,nc);
+ }
+
+
+ // Add WELL EQUATIONS
+ ADB qs = state.qs;
+ for (int phase = 0; phase < np; ++phase) {
+ qs -= superset(wops_.p2w * cq_s[phase], Span(nw, 1, phase*nw), nw*np);
+
+ }
+
+
+ V cq = superset(cq_s[0].value(), Span(nperf, np, 0), nperf*np);
+ for (int phase = 1; phase < np; ++phase) {
+ cq += superset(cq_s[phase].value(), Span(nperf, np, phase), nperf*np);
+ }
+
+ std::vector cq_d(cq.data(), cq.data() + nperf*np);
+ xw.perfPhaseRates() = cq_d;
+
+ residual_.well_flux_eq = qs;
+ }
+
+
+
+
+
+ namespace
+ {
+ double rateToCompare(const ADB& well_phase_flow_rate,
+ const int well,
+ const int num_phases,
+ const double* distr)
+ {
+ const int num_wells = well_phase_flow_rate.size() / num_phases;
+ double rate = 0.0;
+ for (int phase = 0; phase < num_phases; ++phase) {
+ // Important: well_phase_flow_rate is ordered with all rates for first
+ // phase coming first, then all for second phase etc.
+ rate += well_phase_flow_rate.value()[well + phase*num_wells] * distr[phase];
+ }
+ return rate;
+ }
+
+ bool constraintBroken(const ADB& bhp,
+ const ADB& well_phase_flow_rate,
+ const int well,
+ const int num_phases,
+ const WellType& well_type,
+ const WellControls* wc,
+ const int ctrl_index)
+ {
+ const WellControlType ctrl_type = well_controls_iget_type(wc, ctrl_index);
+ const double target = well_controls_iget_target(wc, ctrl_index);
+ const double* distr = well_controls_iget_distr(wc, ctrl_index);
+
+ bool broken = false;
+
+ switch (well_type) {
+ case INJECTOR:
+ {
+ switch (ctrl_type) {
+ case BHP:
+ broken = bhp.value()[well] > target;
+ break;
+
+ case RESERVOIR_RATE: // Intentional fall-through
+ case SURFACE_RATE:
+ broken = rateToCompare(well_phase_flow_rate,
+ well, num_phases, distr) > target;
+ break;
+ }
+ }
+ break;
+
+ case PRODUCER:
+ {
+ switch (ctrl_type) {
+ case BHP:
+ broken = bhp.value()[well] < target;
+ break;
+
+ case RESERVOIR_RATE: // Intentional fall-through
+ case SURFACE_RATE:
+ // Note that the rates compared below are negative,
+ // so breaking the constraints means: too high flow rate
+ // (as for injection).
+ broken = rateToCompare(well_phase_flow_rate,
+ well, num_phases, distr) < target;
+ break;
+ }
+ }
+ break;
+
+ default:
+ OPM_THROW(std::logic_error, "Can only handle INJECTOR and PRODUCER wells.");
+ }
+
+ return broken;
+ }
+ } // anonymous namespace
+
+
+
+
+ template
+ void FullyImplicitBlackoilPolymerSolver::updateWellControls(ADB& bhp,
+ ADB& well_phase_flow_rate,
+ WellStateFullyImplicitBlackoil& xw) const
+ {
+ std::string modestring[3] = { "BHP", "RESERVOIR_RATE", "SURFACE_RATE" };
+ // Find, for each well, if any constraints are broken. If so,
+ // switch control to first broken constraint.
+ const int np = wells_.number_of_phases;
+ const int nw = wells_.number_of_wells;
+ bool bhp_changed = false;
+ bool rates_changed = false;
+ for (int w = 0; w < nw; ++w) {
+ const WellControls* wc = wells_.ctrls[w];
+ // The current control in the well state overrides
+ // the current control set in the Wells struct, which
+ // is instead treated as a default.
+ const int current = xw.currentControls()[w];
+ // Loop over all controls except the current one, and also
+ // skip any RESERVOIR_RATE controls, since we cannot
+ // handle those.
+ const int nwc = well_controls_get_num(wc);
+ int ctrl_index = 0;
+ for (; ctrl_index < nwc; ++ctrl_index) {
+ if (ctrl_index == current) {
+ // This is the currently used control, so it is
+ // used as an equation. So this is not used as an
+ // inequality constraint, and therefore skipped.
+ continue;
+ }
+ if (constraintBroken(bhp, well_phase_flow_rate, w, np, wells_.type[w], wc, ctrl_index)) {
+ // ctrl_index will be the index of the broken constraint after the loop.
+ break;
+ }
+ }
+ if (ctrl_index != nwc) {
+ // Constraint number ctrl_index was broken, switch to it.
+ std::cout << "Switching control mode for well " << wells_.name[w]
+ << " from " << modestring[well_controls_iget_type(wc, current)]
+ << " to " << modestring[well_controls_iget_type(wc, ctrl_index)] << std::endl;
+ xw.currentControls()[w] = ctrl_index;
+ // Also updating well state and primary variables.
+ // We can only be switching to BHP and SURFACE_RATE
+ // controls since we do not support RESERVOIR_RATE.
+ const double target = well_controls_iget_target(wc, ctrl_index);
+ const double* distr = well_controls_iget_distr(wc, ctrl_index);
+ switch (well_controls_iget_type(wc, ctrl_index)) {
+ case BHP:
+ xw.bhp()[w] = target;
+ bhp_changed = true;
+ break;
+
+ case RESERVOIR_RATE:
+ // No direct change to any observable quantity at
+ // surface condition. In this case, use existing
+ // flow rates as initial conditions as reservoir
+ // rate acts only in aggregate.
+ //
+ // Just record the fact that we need to recompute
+ // the 'well_phase_flow_rate'.
+ rates_changed = true;
+ break;
+
+ case SURFACE_RATE:
+ for (int phase = 0; phase < np; ++phase) {
+ if (distr[phase] > 0.0) {
+ xw.wellRates()[np*w + phase] = target * distr[phase];
+ }
+ }
+ rates_changed = true;
+ break;
+ }
+ }
+ }
+
+ // Update primary variables, if necessary.
+ if (bhp_changed) {
+ ADB::V new_bhp = Eigen::Map(xw.bhp().data(), nw);
+ bhp = ADB::function(new_bhp, bhp.derivative());
+ }
+ if (rates_changed) {
+ // Need to reshuffle well rates, from phase running fastest
+ // to wells running fastest.
+ // The transpose() below switches the ordering.
+ const DataBlock wrates = Eigen::Map(xw.wellRates().data(), nw, np).transpose();
+ const ADB::V new_qs = Eigen::Map(wrates.data(), nw*np);
+ well_phase_flow_rate = ADB::function(new_qs, well_phase_flow_rate.derivative());
+ }
+ }
+
+
+
+
+ template
+ void FullyImplicitBlackoilPolymerSolver::addWellControlEq(const SolutionState& state,
+ const WellStateFullyImplicitBlackoil& xw,
+ const V& aliveWells)
+ {
+ const int np = wells_.number_of_phases;
+ const int nw = wells_.number_of_wells;
+
+ V bhp_targets = V::Zero(nw);
+ V rate_targets = V::Zero(nw);
+ M rate_distr(nw, np*nw);
+ for (int w = 0; w < nw; ++w) {
+ const WellControls* wc = wells_.ctrls[w];
+ // The current control in the well state overrides
+ // the current control set in the Wells struct, which
+ // is instead treated as a default.
+ const int current = xw.currentControls()[w];
+
+ switch (well_controls_iget_type(wc, current)) {
+ case BHP:
+ {
+ bhp_targets (w) = well_controls_iget_target(wc, current);
+ rate_targets(w) = -1e100;
+ }
+ break;
+
+ case RESERVOIR_RATE: // Intentional fall-through
+ case SURFACE_RATE:
+ {
+ // RESERVOIR and SURFACE rates look the same, from a
+ // high-level point of view, in the system of
+ // simultaneous linear equations.
+
+ const double* const distr =
+ well_controls_iget_distr(wc, current);
+
+ for (int p = 0; p < np; ++p) {
+ rate_distr.insert(w, p*nw + w) = distr[p];
+ }
+
+ bhp_targets (w) = -1.0e100;
+ rate_targets(w) = well_controls_iget_target(wc, current);
+ }
+ break;
+ }
+ }
+ const ADB bhp_residual = state.bhp - bhp_targets;
+ const ADB rate_residual = rate_distr * state.qs - rate_targets;
+ // Choose bhp residual for positive bhp targets.
+ Selector bhp_selector(bhp_targets);
+ residual_.well_eq = bhp_selector.select(bhp_residual, rate_residual);
+ // For wells that are dead (not flowing), and therefore not communicating
+ // with the reservoir, we set the equation to be equal to the well's total
+ // flow. This will be a solution only if the target rate is also zero.
+ M rate_summer(nw, np*nw);
+ for (int w = 0; w < nw; ++w) {
+ for (int phase = 0; phase < np; ++phase) {
+ rate_summer.insert(w, phase*nw + w) = 1.0;
+ }
+ }
+ Selector alive_selector(aliveWells, Selector::NotEqualZero);
+ residual_.well_eq = alive_selector.select(residual_.well_eq, rate_summer * state.qs);
+ // DUMP(residual_.well_eq);
+ }
+
+
+
+
+
+ template
+ V FullyImplicitBlackoilPolymerSolver::solveJacobianSystem() const
+ {
+ return linsolver_.computeNewtonIncrement(residual_);
+ }
+
+
+
+
+
+ namespace {
+ struct Chop01 {
+ double operator()(double x) const { return std::max(std::min(x, 1.0), 0.0); }
+ };
+ }
+
+
+
+
+
+ template
+ void FullyImplicitBlackoilPolymerSolver::updateState(const V& dx,
+ PolymerBlackoilState& state,
+ WellStateFullyImplicitBlackoil& well_state)
+ {
+ using namespace Opm::AutoDiffGrid;
+ const int np = fluid_.numPhases();
+ const int nc = numCells(grid_);
+ 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);
+
+ // store cell status in vectors
+ V isRs = V::Zero(nc,1);
+ V isRv = V::Zero(nc,1);
+ V isSg = V::Zero(nc,1);
+ if (active_[Gas]) {
+ for (int c = 0; c < nc; ++c) {
+ switch (primalVariable_[c]) {
+ case PrimalVariables::RS:
+ isRs[c] = 1;
+ break;
+
+ case PrimalVariables::RV:
+ isRv[c] = 1;
+ break;
+
+ default:
+ isSg[c] = 1;
+ break;
+ }
+ }
+ }
+
+ // Extract parts of dx corresponding to each part.
+ const V dp = subset(dx, Span(nc));
+ int varstart = nc;
+ const V dsw = active_[Water] ? subset(dx, Span(nc, 1, varstart)) : null;
+ varstart += dsw.size();
+
+ const V dxvar = active_[Gas] ? subset(dx, Span(nc, 1, varstart)): null;
+ varstart += dxvar.size();
+
+ if (has_polymer_) {
+ const V dc = subset(dx, Span(nc, 1, varstart));
+ } else {
+ const V dc = null;
+ }
+ const V dqs = subset(dx, Span(np*nw, 1, varstart));
+ varstart += dqs.size();
+ const V dbhp = subset(dx, Span(nw, 1, varstart));
+ varstart += dbhp.size();
+ assert(varstart == dx.size());
+
+ // Pressure update.
+ const double dpmaxrel = dpMaxRel();
+ const V p_old = Eigen::Map(&state.pressure()[0], nc, 1);
+ const V absdpmax = dpmaxrel*p_old.abs();
+ const V dp_limited = sign(dp) * dp.abs().min(absdpmax);
+ const V p = (p_old - dp_limited).max(zero);
+ std::copy(&p[0], &p[0] + nc, state.pressure().begin());
+
+
+ // Saturation updates.
+ const Opm::PhaseUsage& pu = fluid_.phaseUsage();
+ const DataBlock s_old = Eigen::Map(& state.saturation()[0], nc, np);
+ const double dsmax = dsMax();
+ V so;
+ V sw;
+ V sg;
+
+ {
+ V maxVal = zero;
+ V dso = zero;
+ if (active_[Water]){
+ maxVal = dsw.abs().max(maxVal);
+ dso = dso - dsw;
+ }
+
+ V dsg;
+ if (active_[Gas]){
+ dsg = isSg * dxvar - isRv * dsw;
+ maxVal = dsg.abs().max(maxVal);
+ dso = dso - dsg;
+ }
+
+ maxVal = dso.abs().max(maxVal);
+
+ V step = dsmax/maxVal;
+ step = step.min(1.);
+
+ if (active_[Water]) {
+ const int pos = pu.phase_pos[ Water ];
+ const V sw_old = s_old.col(pos);
+ sw = sw_old - step * dsw;
+ }
+
+ if (active_[Gas]) {
+ const int pos = pu.phase_pos[ Gas ];
+ const V sg_old = s_old.col(pos);
+ sg = sg_old - step * dsg;
+ }
+
+ const int pos = pu.phase_pos[ Oil ];
+ const V so_old = s_old.col(pos);
+ so = so_old - step * dso;
+ }
+
+ // Appleyard chop process.
+ auto ixg = sg < 0;
+ for (int c = 0; c < nc; ++c) {
+ if (ixg[c]) {
+ sw[c] = sw[c] / (1-sg[c]);
+ so[c] = so[c] / (1-sg[c]);
+ sg[c] = 0;
+ }
+ }
+
+
+ auto ixo = so < 0;
+ for (int c = 0; c < nc; ++c) {
+ if (ixo[c]) {
+ sw[c] = sw[c] / (1-so[c]);
+ sg[c] = sg[c] / (1-so[c]);
+ so[c] = 0;
+ }
+ }
+
+ auto ixw = sw < 0;
+ for (int c = 0; c < nc; ++c) {
+ if (ixw[c]) {
+ so[c] = so[c] / (1-sw[c]);
+ sg[c] = sg[c] / (1-sw[c]);
+ sw[c] = 0;
+ }
+ }
+
+ const V sumSat = sw + so + sg;
+ sw = sw / sumSat;
+ so = so / sumSat;
+ sg = sg / sumSat;
+
+ // Update the state
+ for (int c = 0; c < nc; ++c) {
+ state.saturation()[c*np + pu.phase_pos[ Water ]] = sw[c];
+ }
+
+ for (int c = 0; c < nc; ++c) {
+ state.saturation()[c*np + pu.phase_pos[ Gas ]] = sg[c];
+ }
+
+ if (active_[ Oil ]) {
+ const int pos = pu.phase_pos[ Oil ];
+ for (int c = 0; c < nc; ++c) {
+ state.saturation()[c*np + pos] = so[c];
+ }
+ }
+
+ // Update rs and rv
+ const double drsmaxrel = drsMaxRel();
+ const double drvmax = 1e9;//% same as in Mrst
+ V rs;
+ if (has_disgas_) {
+ const V rs_old = Eigen::Map(&state.gasoilratio()[0], nc);
+ const V drs = isRs * dxvar;
+ const V drs_limited = sign(drs) * drs.abs().min(rs_old.abs()*drsmaxrel);
+ rs = rs_old - drs_limited;
+ }
+ V rv;
+ if (has_vapoil_) {
+ const V rv_old = Eigen::Map(&state.rv()[0], nc);
+ const V drv = isRv * dxvar;
+ const V drv_limited = sign(drv) * drv.abs().min(drvmax);
+ rv = rv_old - drv_limited;
+ }
+
+ // Update the state
+ if (has_disgas_)
+ std::copy(&rs[0], &rs[0] + nc, state.gasoilratio().begin());
+
+ if (has_vapoil_)
+ std::copy(&rv[0], &rv[0] + nc, state.rv().begin());
+
+ // Sg is used as primal variable for water only cells.
+ const double epsilon = std::sqrt(std::numeric_limits::epsilon());
+ auto watOnly = sw > (1 - epsilon);
+
+ // phase translation sg <-> rs
+ const V rsSat0 = fluidRsSat(p_old, s_old.col(pu.phase_pos[Oil]), cells_);
+ const V rsSat = fluidRsSat(p, so, cells_);
+
+ std::fill(primalVariable_.begin(), primalVariable_.end(), PrimalVariables::Sg);
+
+ if (has_disgas_) {
+ // The obvious case
+ auto hasGas = (sg > 0 && isRs == 0);
+
+ // keep oil saturated if previous sg is sufficient large:
+ const int pos = pu.phase_pos[ Gas ];
+ auto hadGas = (sg <= 0 && s_old.col(pos) > epsilon);
+ // Set oil saturated if previous rs is sufficiently large
+ const V rs_old = Eigen::Map(&state.gasoilratio()[0], nc);
+ auto gasVaporized = ( (rs > rsSat * (1+epsilon) && isRs == 1 ) && (rs_old > rsSat0 * (1-epsilon)) );
+
+ auto useSg = watOnly || hasGas || hadGas || gasVaporized;
+ for (int c = 0; c < nc; ++c) {
+ if (useSg[c]) { rs[c] = rsSat[c];}
+ else { primalVariable_[c] = PrimalVariables::RS; }
+
+ }
+ }
+
+ // phase transitions so <-> rv
+ const V rvSat0 = fluidRvSat(p_old, s_old.col(pu.phase_pos[Oil]), cells_);
+ const V rvSat = fluidRvSat(p, so, cells_);
+
+ if (has_vapoil_) {
+ // The obvious case
+ auto hasOil = (so > 0 && isRv == 0);
+
+ // keep oil saturated if previous so is sufficient large:
+ const int pos = pu.phase_pos[ Oil ];
+ auto hadOil = (so <= 0 && s_old.col(pos) > epsilon );
+ // Set oil saturated if previous rv is sufficiently large
+ const V rv_old = Eigen::Map(&state.rv()[0], nc);
+ auto oilCondensed = ( (rv > rvSat * (1+epsilon) && isRv == 1) && (rv_old > rvSat0 * (1-epsilon)) );
+ auto useSg = watOnly || hasOil || hadOil || oilCondensed;
+ for (int c = 0; c < nc; ++c) {
+ if (useSg[c]) { rv[c] = rvSat[c]; }
+ else {primalVariable_[c] = PrimalVariables::RV; }
+
+ }
+
+ }
+
+ //Polymer concentration updates.
+ if (has_polymer_) {
+ const V c_old = Eigen::Map(&state.concentration()[0], nc, 1);
+ const V c = (c_old - dc).max(zero);
+ std::copy(&c[0], &c[0] + nc, state.concentration().begin());
+ }
+
+ // Qs update.
+ // Since we need to update the wellrates, that are ordered by wells,
+ // from dqs which are ordered by phase, the simplest is to compute
+ // dwr, which is the data from dqs but ordered by wells.
+ const DataBlock wwr = Eigen::Map(dqs.data(), np, nw).transpose();
+ const V dwr = Eigen::Map(wwr.data(), nw*np);
+ const V wr_old = Eigen::Map(&well_state.wellRates()[0], nw*np);
+ const V wr = wr_old - dwr;
+ std::copy(&wr[0], &wr[0] + wr.size(), well_state.wellRates().begin());
+
+ // Bhp update.
+ const V bhp_old = Eigen::Map(&well_state.bhp()[0], nw, 1);
+ const V dbhp_limited = sign(dbhp) * dbhp.abs().min(bhp_old.abs()*dpmaxrel);
+ const V bhp = bhp_old - dbhp_limited;
+ std::copy(&bhp[0], &bhp[0] + bhp.size(), well_state.bhp().begin());
+
+ // Update phase conditions used for property calculations.
+ updatePhaseCondFromPrimalVariable();
+ }
+
+
+
+
+
+ template
+ std::vector
+ FullyImplicitBlackoilPolymerSolver::computeRelPerm(const SolutionState& state) const
+ {
+ using namespace Opm::AutoDiffGrid;
+ const int nc = numCells(grid_);
+ const std::vector& bpat = state.pressure.blockPattern();
+
+ const ADB null = ADB::constant(V::Zero(nc, 1), bpat);
+
+ const Opm::PhaseUsage& pu = fluid_.phaseUsage();
+ const ADB sw = (active_[ Water ]
+ ? state.saturation[ pu.phase_pos[ Water ] ]
+ : null);
+
+ const ADB so = (active_[ Oil ]
+ ? state.saturation[ pu.phase_pos[ Oil ] ]
+ : null);
+
+ const ADB sg = (active_[ Gas ]
+ ? state.saturation[ pu.phase_pos[ Gas ] ]
+ : null);
+
+ return fluid_.relperm(sw, so, sg, cells_);
+ }
+
+
+ template
+ std::vector
+ FullyImplicitBlackoilPolymerSolver::computePressures(const SolutionState& state) const
+ {
+ using namespace Opm::AutoDiffGrid;
+ const int nc = numCells(grid_);
+ const std::vector& bpat = state.pressure.blockPattern();
+
+ const ADB null = ADB::constant(V::Zero(nc, 1), bpat);
+
+ const Opm::PhaseUsage& pu = fluid_.phaseUsage();
+ const ADB sw = (active_[ Water ]
+ ? state.saturation[ pu.phase_pos[ Water ] ]
+ : null);
+
+ const ADB so = (active_[ Oil ]
+ ? state.saturation[ pu.phase_pos[ Oil ] ]
+ : null);
+
+ const ADB sg = (active_[ Gas ]
+ ? state.saturation[ pu.phase_pos[ Gas ] ]
+ : null);
+
+ // convert the pressure offsets to the capillary pressures
+ std::vector pressure = fluid_.capPress(sw, so, sg, cells_);
+ for (int phaseIdx = 0; phaseIdx < BlackoilPhases::MaxNumPhases; ++phaseIdx) {
+ // The reference pressure is always the liquid phase (oil) pressure.
+ if (phaseIdx == BlackoilPhases::Liquid)
+ continue;
+ pressure[phaseIdx] = pressure[phaseIdx] - pressure[BlackoilPhases::Liquid];
+ }
+
+ // Since pcow = po - pw, but pcog = pg - po,
+ // we have
+ // pw = po - pcow
+ // pg = po + pcgo
+ // This is an unfortunate inconsistency, but a convention we must handle.
+ for (int phaseIdx = 0; phaseIdx < BlackoilPhases::MaxNumPhases; ++phaseIdx) {
+ if (phaseIdx == BlackoilPhases::Aqua) {
+ pressure[phaseIdx] = state.pressure - pressure[phaseIdx];
+ } else {
+ pressure[phaseIdx] += state.pressure;
+ }
+ }
+
+ return pressure;
+ }
+
+
+
+
+ template
+ std::vector
+ FullyImplicitBlackoilPolymerSolver::computeRelPermWells(const SolutionState& state,
+ const DataBlock& well_s,
+ const std::vector& well_cells) const
+ {
+ const int nw = wells_.number_of_wells;
+ const int nperf = wells_.well_connpos[nw];
+ const std::vector& bpat = state.pressure.blockPattern();
+
+ const ADB null = ADB::constant(V::Zero(nperf), bpat);
+
+ const Opm::PhaseUsage& pu = fluid_.phaseUsage();
+ const ADB sw = (active_[ Water ]
+ ? ADB::constant(well_s.col(pu.phase_pos[ Water ]), bpat)
+ : null);
+
+ const ADB so = (active_[ Oil ]
+ ? ADB::constant(well_s.col(pu.phase_pos[ Oil ]), bpat)
+ : null);
+
+ const ADB sg = (active_[ Gas ]
+ ? ADB::constant(well_s.col(pu.phase_pos[ Gas ]), bpat)
+ : null);
+
+ return fluid_.relperm(sw, so, sg, well_cells);
+ }
+
+
+
+
+
+ template
+ void
+ FullyImplicitBlackoilPolymerSolver::computeMassFlux(const V& transi,
+ const std::vector& kr ,
+ const std::vector& phasePressure,
+ const SolutionState& state)
+ {
+ if (has_polymer_) {
+ const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern());
+ ADB krw_eff = ADB::null();
+ const ADB mc = computeMc(state);
+ ADB inv_wat_eff_vis = ADB::null();
+ }
+
+ for (int phase = 0; phase < fluid_.numPhases(); ++phase) {
+ const int canonicalPhaseIdx = canph_[phase];
+ const std::vector cond = phaseCondition();
+
+ const ADB tr_mult = transMult(state.pressure);
+ const ADB mu = fluidViscosity(canonicalPhaseIdx, phasePressure[canonicalPhaseIdx], state.rs, state.rv, cond, cells_);
+ rq_[canonicalPhaseIdx].mob = tr_mult * kr[canonicalPhaseIdx] / mu;
+ if (cononicalPhaseIdx == Water) {
+ if(has_polymer_) {
+ krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration,
+ cmax,
+ kr[cononicalPhaseIdx],
+ state.saturation[cononicalPhaseIdx]);
+ inv_wat_eff_vis = polymer_props_ad_.effectiveInvWaterVisc(state.concentration, mu.value().data());
+ rq_[canonicalPhaseIdx].mob = tr_mult * krw_eff * inv_wat_eff_visc;
+ }
+ }
+
+ const ADB rho = fluidDensity(canonicalPhaseIdx, phasePressure[canonicalPhaseIdx], state.rs, state.rv, cond, cells_);
+
+ ADB& head = rq_[canonicalPhaseIdx].head;
+
+ // compute gravity potensial using the face average as in eclipse and MRST
+ const ADB rhoavg = ops_.caver * rho;
+
+ ADB dp = ops_.ngrad * phasePressure[canonicalPhaseIdx] - geo_.gravity()[2] * (rhoavg * (ops_.ngrad * geo_.z().matrix()));
+
+ if (use_threshold_pressure_) {
+ applyThresholdPressures(dp);
+ }
+
+ head = transi*dp;
+ //head = transi*(ops_.ngrad * phasePressure) + gflux;
+
+ UpwindSelector upwind(grid_, ops_, head.value());
+
+ const ADB& b = rq_[canonicalPhaseIdx].b;
+ const ADB& mob = rq_[canonicalPhaseIdx].mob;
+ rq_[canonicalPhaseIdx].mflux = upwind.select(b * mob) * head;
+ }
+
+ if (has_polymer_) {
+ rq_[poly_pos_].mob = tr_mult * mc * krw_eff * inv_wat_eff_visc;
+ rq_[poly_pos_].b = rq_[canph_[Water]].b;
+ rq_[poly_pos_].head = rq_[canph_[Water]].head;
+ UpwindSelector upwind(grid_, ops_, rq_[poly_pos_].head.value());
+ rq_[poly_pos_].mflux = upwind.select(rq_[poly_pos_].b * rq_[poly_pos_].mob) * rq_[poly_pos_].head;
+ }
+ }
+
+
+
+ template
+ void
+ FullyImplicitBlackoilPolymerSolver::applyThresholdPressures(ADB& dp)
+ {
+ // We support reversible threshold pressures only.
+ // Method: if the potential difference is lower (in absolute
+ // value) than the threshold for any face, then the potential
+ // (and derivatives) is set to zero. If it is above the
+ // threshold, the threshold pressure is subtracted from the
+ // absolute potential (the potential is moved towards zero).
+
+ // Identify the set of faces where the potential is under the
+ // threshold, that shall have zero flow. Storing the bool
+ // Array as a V (a double Array) with 1 and 0 elements, a
+ // 1 where flow is allowed, a 0 where it is not.
+ const V high_potential = (dp.value().abs() >= threshold_pressures_by_interior_face_).template cast();
+
+ // Create a sparse vector that nullifies the low potential elements.
+ const M keep_high_potential = spdiag(high_potential);
+
+ // Find the current sign for the threshold modification
+ const V sign_dp = sign(dp.value());
+ const V threshold_modification = sign_dp * threshold_pressures_by_interior_face_;
+
+ // Modify potential and nullify where appropriate.
+ dp = keep_high_potential * (dp - threshold_modification);
+ }
+
+
+
+
+
+ template
+ double
+ FullyImplicitBlackoilPolymerSolver::residualNorm() const
+ {
+ double globalNorm = 0;
+ std::vector::const_iterator quantityIt = residual_.material_balance_eq.begin();
+ const std::vector::const_iterator endQuantityIt = residual_.material_balance_eq.end();
+ for (; quantityIt != endQuantityIt; ++quantityIt) {
+ const double quantityResid = (*quantityIt).value().matrix().norm();
+ if (!std::isfinite(quantityResid)) {
+ const int trouble_phase = quantityIt - residual_.material_balance_eq.begin();
+ OPM_THROW(Opm::NumericalProblem,
+ "Encountered a non-finite residual in material balance equation "
+ << trouble_phase);
+ }
+ globalNorm = std::max(globalNorm, quantityResid);
+ }
+ globalNorm = std::max(globalNorm, residual_.well_flux_eq.value().matrix().norm());
+ globalNorm = std::max(globalNorm, residual_.well_eq.value().matrix().norm());
+
+ return globalNorm;
+ }
+
+
+ template
+ std::vector
+ FullyImplicitBlackoilPolymerSolver::residuals() const
+ {
+ std::vector residual;
+
+ std::vector::const_iterator massBalanceIt = residual_.material_balance_eq.begin();
+ const std::vector::const_iterator endMassBalanceIt = residual_.material_balance_eq.end();
+
+ for (; massBalanceIt != endMassBalanceIt; ++massBalanceIt) {
+ const double massBalanceResid = (*massBalanceIt).value().matrix().template lpNorm();
+ if (!std::isfinite(massBalanceResid)) {
+ OPM_THROW(Opm::NumericalProblem,
+ "Encountered a non-finite residual");
+ }
+ residual.push_back(massBalanceResid);
+ }
+
+ // the following residuals are not used in the oscillation detection now
+ const double wellFluxResid = residual_.well_flux_eq.value().matrix().template lpNorm();
+ if (!std::isfinite(wellFluxResid)) {
+ OPM_THROW(Opm::NumericalProblem,
+ "Encountered a non-finite residual");
+ }
+ residual.push_back(wellFluxResid);
+
+ const double wellResid = residual_.well_eq.value().matrix().template lpNorm();
+ if (!std::isfinite(wellResid)) {
+ OPM_THROW(Opm::NumericalProblem,
+ "Encountered a non-finite residual");
+ }
+ residual.push_back(wellResid);
+
+ return residual;
+ }
+
+ template
+ void
+ FullyImplicitBlackoilPolymerSolver::detectNewtonOscillations(const std::vector>& residual_history,
+ const int it, const double relaxRelTol,
+ bool& oscillate, bool& stagnate) const
+ {
+ // The detection of oscillation in two primary variable results in the report of the detection
+ // of oscillation for the solver.
+ // Only the saturations are used for oscillation detection for the black oil model.
+ // Stagnate is not used for any treatment here.
+
+ if ( it < 2 ) {
+ oscillate = false;
+ stagnate = false;
+ return;
+ }
+
+ stagnate = true;
+ int oscillatePhase = 0;
+
+ for (int phaseIdx= 0; phaseIdx < fluid_.numPhases(); ++ phaseIdx){
+ if (active_[phaseIdx]) {
+ double relChange1 = std::fabs((residual_history[it][phaseIdx] - residual_history[it - 2][phaseIdx]) /
+ residual_history[it][phaseIdx]);
+ double relChange2 = std::fabs((residual_history[it][phaseIdx] - residual_history[it - 1][phaseIdx]) /
+ residual_history[it][phaseIdx]);
+ oscillatePhase += (relChange1 < relaxRelTol) && (relChange2 > relaxRelTol);
+
+ double relChange3 = std::fabs((residual_history[it - 1][phaseIdx] - residual_history[it - 2][phaseIdx]) /
+ residual_history[it - 2][phaseIdx]);
+ if (relChange3 > 1.e-3) {
+ stagnate = false;
+ }
+ }
+ }
+
+ oscillate = (oscillatePhase > 1);
+ }
+
+
+ template
+ void
+ FullyImplicitBlackoilPolymerSolver::stablizeNewton(V& dx, V& dxOld, const double omega,
+ const RelaxType relax_type) const
+ {
+ // The dxOld is updated with dx.
+ // If omega is equal to 1., no relaxtion will be appiled.
+
+ const V tempDxOld = dxOld;
+ dxOld = dx;
+
+ switch (relax_type) {
+ case DAMPEN:
+ if (omega == 1.) {
+ return;
+ }
+ dx = dx*omega;
+ return;
+ case SOR:
+ if (omega == 1.) {
+ return;
+ }
+ dx = dx*omega + (1.-omega)*tempDxOld;
+ return;
+ default:
+ OPM_THROW(std::runtime_error, "Can only handle DAMPEN and SOR relaxation type.");
+ }
+
+ return;
+ }
+
+ template
+ bool
+ FullyImplicitBlackoilPolymerSolver::getConvergence(const double dt)
+ {
+ const double tol_mb = 1.0e-7;
+ const double tol_cnv = 1.0e-3;
+
+ const int nc = Opm::AutoDiffGrid::numCells(grid_);
+ const Opm::PhaseUsage& pu = fluid_.phaseUsage();
+
+ const V pv = geo_.poreVolume();
+ const double pvSum = pv.sum();
+
+ const std::vector cond = phaseCondition();
+
+ double CNVW = 0.;
+ double CNVO = 0.;
+ double CNVG = 0.;
+
+ double RW_sum = 0.;
+ double RO_sum = 0.;
+ double RG_sum = 0.;
+
+ double BW_avg = 0.;
+ double BO_avg = 0.;
+ double BG_avg = 0.;
+
+ if (active_[Water]) {
+ const int pos = pu.phase_pos[Water];
+ const ADB& tempBW = rq_[pos].b;
+ V BW = 1./tempBW.value();
+ V RW = residual_.material_balance_eq[Water].value();
+ BW_avg = BW.sum()/nc;
+ const V tempV = RW.abs()/pv;
+
+ CNVW = BW_avg * dt * tempV.maxCoeff();
+ RW_sum = RW.sum();
+ }
+
+ if (active_[Oil]) {
+ // Omit the disgas here. We should add it.
+ const int pos = pu.phase_pos[Oil];
+ const ADB& tempBO = rq_[pos].b;
+ V BO = 1./tempBO.value();
+ V RO = residual_.material_balance_eq[Oil].value();
+ BO_avg = BO.sum()/nc;
+ const V tempV = RO.abs()/pv;
+
+ CNVO = BO_avg * dt * tempV.maxCoeff();
+ RO_sum = RO.sum();
+ }
+
+ if (active_[Gas]) {
+ // Omit the vapoil here. We should add it.
+ const int pos = pu.phase_pos[Gas];
+ const ADB& tempBG = rq_[pos].b;
+ V BG = 1./tempBG.value();
+ V RG = residual_.material_balance_eq[Gas].value();
+ BG_avg = BG.sum()/nc;
+ const V tempV = RG.abs()/pv;
+
+ CNVG = BG_avg * dt * tempV.maxCoeff();
+ RG_sum = RG.sum();
+ }
+
+ double tempValue = tol_mb * pvSum /dt;
+
+ bool converged_MB = (fabs(BW_avg*RW_sum) < tempValue)
+ && (fabs(BO_avg*RO_sum) < tempValue)
+ && (fabs(BG_avg*RG_sum) < tempValue);
+
+ bool converged_CNV = (CNVW < tol_cnv) && (CNVO < tol_cnv) && (CNVG < tol_cnv);
+
+ double residualWellFlux = residual_.well_flux_eq.value().matrix().template lpNorm();
+ double residualWell = residual_.well_eq.value().matrix().template lpNorm();
+
+ bool converged_Well = (residualWellFlux < 1./Opm::unit::day) && (residualWell < Opm::unit::barsa);
+
+ bool converged = converged_MB && converged_CNV && converged_Well;
+
+#ifdef OPM_VERBOSE
+ std::cout << " CNVW " << CNVW << " CNVO " << CNVO << " CNVG " << CNVG << std::endl;
+ std::cout << " converged_MB " << converged_MB << " converged_CNV " << converged_CNV
+ << " converged_Well " << converged_Well << " converged " << converged << std::endl;
+#endif
+ return converged;
+ }
+
+
+ template
+ ADB
+ FullyImplicitBlackoilPolymerSolver::fluidViscosity(const int phase,
+ const ADB& p ,
+ const ADB& rs ,
+ const ADB& rv ,
+ const std::vector& cond,
+ const std::vector& cells) const
+ {
+ switch (phase) {
+ case Water:
+ return fluid_.muWat(p, cells);
+ case Oil: {
+ return fluid_.muOil(p, rs, cond, cells);
+ }
+ case Gas:
+ return fluid_.muGas(p, rv, cond, cells);
+ default:
+ OPM_THROW(std::runtime_error, "Unknown phase index " << phase);
+ }
+ }
+
+
+
+
+
+ template
+ ADB
+ FullyImplicitBlackoilPolymerSolver::fluidReciprocFVF(const int phase,
+ const ADB& p ,
+ const ADB& rs ,
+ const ADB& rv ,
+ const std::vector& cond,
+ const std::vector& cells) const
+ {
+ switch (phase) {
+ case Water:
+ return fluid_.bWat(p, cells);
+ case Oil: {
+ return fluid_.bOil(p, rs, cond, cells);
+ }
+ case Gas:
+ return fluid_.bGas(p, rv, cond, cells);
+ default:
+ OPM_THROW(std::runtime_error, "Unknown phase index " << phase);
+ }
+ }
+
+
+
+
+
+ template
+ ADB
+ FullyImplicitBlackoilPolymerSolver::fluidDensity(const int phase,
+ const ADB& p ,
+ const ADB& rs ,
+ const ADB& rv ,
+ const std::vector& cond,
+ const std::vector& cells) const
+ {
+ const double* rhos = fluid_.surfaceDensity();
+ ADB b = fluidReciprocFVF(phase, p, rs, rv, cond, cells);
+ ADB rho = V::Constant(p.size(), 1, rhos[phase]) * b;
+ if (phase == Oil && active_[Gas]) {
+ // It is correct to index into rhos with canonical phase indices.
+ rho += V::Constant(p.size(), 1, rhos[Gas]) * rs * b;
+ }
+ if (phase == Gas && active_[Oil]) {
+ // It is correct to index into rhos with canonical phase indices.
+ rho += V::Constant(p.size(), 1, rhos[Oil]) * rv * b;
+ }
+ return rho;
+ }
+
+
+
+
+
+ template
+ V
+ FullyImplicitBlackoilPolymerSolver::fluidRsSat(const V& p,
+ const V& satOil,
+ const std::vector& cells) const
+ {
+ return fluid_.rsSat(p, satOil, cells);
+ }
+
+
+
+
+
+ template
+ ADB
+ FullyImplicitBlackoilPolymerSolver::fluidRsSat(const ADB& p,
+ const ADB& satOil,
+ const std::vector& cells) const
+ {
+ return fluid_.rsSat(p, satOil, cells);
+ }
+
+ template
+ V
+ FullyImplicitBlackoilPolymerSolver::fluidRvSat(const V& p,
+ const V& satOil,
+ const std::vector& cells) const
+ {
+ return fluid_.rvSat(p, satOil, cells);
+ }
+
+
+
+
+
+ template
+ ADB
+ FullyImplicitBlackoilPolymerSolver::fluidRvSat(const ADB& p,
+ const ADB& satOil,
+ const std::vector& cells) const
+ {
+ return fluid_.rvSat(p, satOil, cells);
+ }
+
+
+
+
+
+ template
+ ADB
+ FullyImplicitBlackoilPolymerSovler::computeMc(const SolutionState& state) const
+ {
+ return polymer_props_ad_.polymerWaterVelocityRatio(state.concentration);
+ }
+
+
+
+
+ template
+ ADB
+ FullyImplicitBlackoilPolymerSolver::poroMult(const ADB& p) const
+ {
+ const int n = p.size();
+ if (rock_comp_props_ && rock_comp_props_->isActive()) {
+ V pm(n);
+ V dpm(n);
+ for (int i = 0; i < n; ++i) {
+ pm[i] = rock_comp_props_->poroMult(p.value()[i]);
+ dpm[i] = rock_comp_props_->poroMultDeriv(p.value()[i]);
+ }
+ ADB::M dpm_diag = spdiag(dpm);
+ const int num_blocks = p.numBlocks();
+ std::vector jacs(num_blocks);
+ for (int block = 0; block < num_blocks; ++block) {
+ jacs[block] = dpm_diag * p.derivative()[block];
+ }
+ return ADB::function(pm, jacs);
+ else {
+ return ADB::constant(V::Constant(n, 1.0), p.blockPattern());
+ }
+ }
+
+
+
+
+
+ template
+ ADB
+ FullyImplicitBlackoilPolymerSolver::transMult(const ADB& p) const
+ {
+ const int n = p.size();
+ if (rock_comp_props_ && rock_comp_props_->isActive()) {
+ V tm(n);
+ V dtm(n);
+ for (int i = 0; i < n; ++i) {
+ tm[i] = rock_comp_props_->transMult(p.value()[i]);
+ dtm[i] = rock_comp_props_->transMultDeriv(p.value()[i]);
+ }
+ ADB::M dtm_diag = spdiag(dtm);
+ const int num_blocks = p.numBlocks();
+ std::vector jacs(num_blocks);
+ for (int block = 0; block < num_blocks; ++block) {
+ jacs[block] = dtm_diag * p.derivative()[block];
+ }
+ return ADB::function(tm, jacs);
+ } else {
+ return ADB::constant(V::Constant(n, 1.0), p.blockPattern());
+ }
+ }
+
+
+ /*
+ template
+ void
+ FullyImplicitBlackoilPolymerSolver::
+ classifyCondition(const SolutionState& state,
+ std::vector