From 080116c66ba780382846c9714b20a491e179f8e1 Mon Sep 17 00:00:00 2001 From: Liu Ming Date: Fri, 6 Dec 2013 23:35:13 +0800 Subject: [PATCH] Add fully implicit solver for incomp two phase with polymer and the polymer properties based on AD. --- .../FullyImplicitTwophasePolymerSolver.cpp | 485 ++++++++++++++++++ .../FullyImplicitTwophasePolymerSolver.hpp | 98 ++++ opm/polymer/fullyimplicit/PolymerPropsAd.cpp | 77 ++- opm/polymer/fullyimplicit/PolymerPropsAd.hpp | 9 +- 4 files changed, 655 insertions(+), 14 deletions(-) create mode 100644 opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp create mode 100644 opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp diff --git a/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp b/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp new file mode 100644 index 000000000..a2534ba1c --- /dev/null +++ b/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp @@ -0,0 +1,485 @@ +/**/ + +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +namespace Opm { + +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); } + }; + +}//anonymous namespace + + + + + +typedef AutoDiffBlock ADB; +typedef ADB::V V; +typedef ADB::M M; +typedef Eigen::Array DataBlock; + + + + + + FullyImplicitTwoPhaseSolver:: + FullyImplicitTwoPhaseSolver(const UnstructuredGrid& grid, + const IncompPropsAdInterface& fluid, + const PolymerProperties& polymer_props, + const PolymerPropsAd& polymer_props_ad, + const LinearSolverInterface& linsolver) + : grid_ (grid) + , fluid_(fluid) + , polymer_props_ (polymer_props) + , polymer_props_ad_ (polymer_props_ad_) + , linsolver_(linsolver) + , cells_ (buildAllCells(grid.number_of_cells)) + , ops_(grid) + , residual_(std::vector(fluid.numPhases() + 1, ADB::null())) + { + } + + + + + + void + FullyImplicitTwoPhaseSolver:: + step(const double dt, + PolymerState& x, + const std::vector& src) + { + + 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); + const double atol = 1.0e-12; + const double rtol = 5.0e-8; + const int maxit = 15; + + assemble(pvdt, old_state, x, 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); + + assemble(pvdt, old_state, x, 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()) + , concentration ( ADB::null()) + { + } + + + + + + FullyImplicitTwoPhaseSolver::SolutionState + FullyImplicitTwoPhaseSolver::constantState(const PolymerState& x) + { + const int nc = grid_.number_of_cells; + const int np = x.numPhases(); + + SolutionState state(np); + + // Pressure. + assert (not x.pressure().empty()); + const V p = Eigen::Map(& x.pressure()[0], nc); + state.pressure = ADB::constant(p); + + // Saturation. + assert (not x.saturation().empty()); + const DataBlock s_all = Eigen::Map(& x.saturation()[0], nc, np); + for (int phase = 0; phase < np; ++phase) { + state.saturation[phase] = ADB::constant(s_all.col(phase)); + // state.saturation[1] = ADB::constant(s_all.col(1)); + } + + // Concentration + + assert(not x.concentration().empty()); + const V c = Eigen::Map(&x.concentration()[0], nc); + + state.concentration = ADB::constant(c); + return state; + } + + + + + + FullyImplicitTwoPhaseSolver::SolutionState + FullyImplicitTwoPhaseSolver::variableState(const PolymerState& x) + { + 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 saturation. + assert (not x.concentration().empty()); + const V c = Eigen::Map(&x.concentration()[0], nc); + vars0.push_back(c); + + 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; + } + + // Concentration. + state.concentration = vars[nextvar++]; + + assert(nextvar == int(vars.size())); + + return state; + } + + + + + + void + FullyImplicitTwoPhaseSolver:: + assemble(const V& pvdt, + const SolutionState& old_state, + const PolymerState& x , + const std::vector& src) + { + // Create the primary variables. + const SolutionState state = variableState(x); + + // -------- Mass balance equations for water and oil -------- + const V trans = subset(transmissibility(), ops_.internal_faces); + const std::vector kr = computeRelPerm(state); + for (int phase = 0; phase < fluid_.numPhases(); ++phase) { + const ADB mflux = computeMassFlux(phase, trans, kr, state); + ADB source = accumSource(phase, kr, src); + residual_[phase] = + pvdt*(state.saturation[phase] - old_state.saturation[phase]) + + ops_.div*mflux - source; + } + // Mass balance equation for polymer + ADB polysrc = accumPolymerSource(kr, src); + ADB mc = computeMc(); + ADB mflux = computeMassFlux(0, trans, kr, state); + residual_[2] = pvdt * (state.saturation[0] * state.concentration - old_state.saturation[0] * old_state.concentration) + ops_.div * state.concentration * mc * mflux - polysrc; + + } + + + + + + + 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, kr); + // 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(int phase, + const std::vector& kr) const + { + const double* mus = fluid_.viscosity(); + ADB mob_phase = kr[phase] / V::Constant(kr[phase].size(), 1, mus[phase]); + ADB mob_wat = kr[0] / V::Constant(kr[0].size(), 1, mus[0]); + ADB mob_oil= kr[1] / V::Constant(kr[1].size(), 1, mus[1]); + ADB total_mob = mob_wat + mob_oil; + ADB f = mob_phase / total_mob; + + 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."); + } + ADB mass_phase_res = vertcat(residual_[0], residual_[1]); + ADB mass_res = collapseJacs(vertcat(mass_phase_res, residual_[2])); + + const Eigen::SparseMatrix matr = mass_res.derivative()[0]; + V dx(V::Zero(mass_res.size())); + Opm::LinearSolverInterface::LinearSolverReport rep + = linsolver_.solve(matr.rows(), matr.nonZeros(), + matr.outerIndexPtr(), matr.innerIndexPtr(), matr.valuePtr(), + mass_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, + PolymerState& state) const + { + const int np = fluid_.numPhases(); + const int nc = grid_.number_of_cells; + 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 dc = subset(dx, Span(nc ,1 varstart)); + varstart += dc.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]; + } + + // Concentration updates. + const V c_old = Eigen::Map(&state.concentration()[0], nc); + const V c = c_old - dc; + std::copy(&c[0], &c[0] + nc, state.concentration().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 + { +// const ADB tr_mult = transMult(state.pressure); + const double* mus = fluid_.viscosity(); + ADB mob = kr[phase] / V::Constant(kr[phase].size(), 1, mus[phase]); + + const ADB dp = ops_.ngrad * state.pressure; + const ADB head = trans * dp; + + UpwindSelector upwind(grid_, ops_, head.value()); + + return upwind.select(mob) * head; + } + + + + + + double + FullyImplicitTwoPhaseSolver::residualNorm() const + { + double r = 0; + for (std::vector::const_iterator + b = residual_.begin(), + e = residual_.end(); + b != e; ++b) + { + r = std::max(r, (*b).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; + } + + + + + + +}//namespace Opm diff --git a/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp b/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp new file mode 100644 index 000000000..8f5df13db --- /dev/null +++ b/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp @@ -0,0 +1,98 @@ +/**/ + +#ifndef OPM_FULLYIMPLICITTWOPHASEPOLYMERSOLVER_HEADER_INCLUDED +#define OPM_FULLYIMPLICITTWOPHASEPOLYMERSOLVER_HEADER_INCLUDED + +#include +#include +#include +#include +#include +#include + + +struct UnstructuredGrid; +namespace Opm { + class LinearSolverInterface; + class PolymerState; + + + class FullyImplicitTwoPhaseSolver + { + public: + FullyImplicitTwoPhaseSolver(const UnstructuredGrid& grid, + const IncompPropsAdInterface& fluid, + const PolymerProperties& polymer_props, + const PolymerPropsAd& polymer_props_ad, + const LinearSolverInterface& linsolver); + + void step(const double dt, + PolymerState& state, + const std::vector& src); + 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 concentration; +// ADB cmax; + }; + const UnstructuredGrid& grid_; + const IncompPropsAdInterface& fluid_; + const PolymerProperties& polymer_props_; + const PolymerPropsAd& polymer_props_ad_; + const LinearSolverInterface& linsolver_; + const std::vector cells_; + HelperOps ops_; + std::vector residual_; + + + SolutionState + constantState(const PolymerState& x); + SolutionState + variableState(const PolymerState& x); + void + assemble(const V& pvdt, + const SolutionState& old_state, + const PolymerState& x, + const std::vector& src); + V solveJacobianSystem() const; + void updateState(const V& dx, + PolymerState& x) const; + std::vector + computeRelPerm(const SolutionState& state) const; + V + transmissibility() const; + ADB + computeFracFlow(int phase, + const std::vector& kr) 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) const; + 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; + }; +} // namespace Opm +#endif// OPM_FULLYIMPLICITTWOPHASESOLVER_HEADER_INCLUDED diff --git a/opm/polymer/fullyimplicit/PolymerPropsAd.cpp b/opm/polymer/fullyimplicit/PolymerPropsAd.cpp index 333709a7c..15adf0149 100644 --- a/opm/polymer/fullyimplicit/PolymerPropsAd.cpp +++ b/opm/polymer/fullyimplicit/PolymerPropsAd.cpp @@ -163,13 +163,28 @@ namespace Opm { } */ + PolymerPropsAd::PolymerPropsAd(const PolymerProperties& polymer_props) + : polymer_props_ (polymer_props) + { + } + + + + + PolymerPropsAd::~PolymerPropsAd() + { + } + + + + V PolymerPropsAd::effectiveInvWaterVisc(const V& c, const double* visc) const { const int nc = c.size(); V inv_mu_w_eff(n); for (int i = 0; i < nc; ++i) { - double im; + double im = 0; polymer_props_.effectiveInvVisc(c(i), visc, im); inv_mu_w_eff(i) = im; } @@ -185,21 +200,21 @@ namespace Opm { ADB PolymerPropsAd::effectiveInvWaterVisc(const ADB& c, const double* visc) { - const int n = c.size(); + const int nc = c.size(); V inv_mu_w_eff(n); V dinv_mu_w_eff(n); - for (int i = 0; i < n; ++i) { - double im, dim; + for (int i = 0; i < nc; ++i) { + double im = 0, dim = 0; polymer_props_.effectiveInvViscWithDer(c.value()(i), visc, im, dim); inv_mu_w_eff(i) = im; dinv_mu_w_eff(i) = dim; } - ADB::M dim_diag = spdiag(dinv_mu_w_eff); - const int num_blocks = c.numBlocks(); - std::vector jacs(num_blocks); - for (int block = 0; block < num_blocks; ++block) { - jacs[block] = dim_diag * c.derivative()[block]; - } + ADB::M dim_diag = spdiag(dinv_mu_w_eff); + const int num_blocks = c.numBlocks(); + std::vector jacs(num_blocks); + for (int block = 0; block < num_blocks; ++block) { + jacs[block] = dim_diag * c.derivative()[block]; + } return ADB::function(inv_mu_w_eff, jacs); } @@ -209,8 +224,46 @@ namespace Opm { V PolymerPropsAd::polymerWaterVelocityRatio(const V& c) const { - const int nc + const int nc = c.size(); + V mc(n); + + for (int i = 0; i < nc; ++i) { + double m = 0; + polymer_props_.computeMc(c(i), m); + mc(i) = m; + } + return mc; } -} // namespace Opm + + + + + ADB PolymerPropsAd::polymerWaterVelocityRatio(const ADB& c) const + { + + const int nc = c.size(); + V mc(n); + V dmc(n); + + for (int i = 0; i < nc; ++i) { + double m = 0; + double dm = 0; + polymer_props_.computeMcWithDer(c(i), m, dm); + + mc(i) = m; + dmc(i) = dm; + } + + ADB::M dmc_diag = spdiag(dmc); + const int num_blocks = c.numBlocks(); + std::vector jacs(num_blocks); + for (int block = 0; block < num_blocks; ++block) { + jacs[block] = dmc_diag * c.derivative()[block]; + } + + return ADB::function(mc, jacs); + } + +}// namespace Opm diff --git a/opm/polymer/fullyimplicit/PolymerPropsAd.hpp b/opm/polymer/fullyimplicit/PolymerPropsAd.hpp index ad3f00471..2d77b228c 100644 --- a/opm/polymer/fullyimplicit/PolymerPropsAd.hpp +++ b/opm/polymer/fullyimplicit/PolymerPropsAd.hpp @@ -65,10 +65,15 @@ namespace Opm { */ typedef AutoDiffBlock ADB; typedef ADB::V V; - PolymerPropsAd(const PolymerPropperties& polyprops); + PolymerPropsAd(const PolymerPropperties& polymer_props); ~PolymerPropsAd(); + V PolymerPropsAd::effectiveInvWaterVisc(const V& c,const double* visc) const; + ADB PolymerPropsAd::effectiveInvWaterVisc(const ADB& c,const double* visc) const; + + V PolymerPropsAd::polymerWaterVelocityRatio(const V& c) const; + ADB PolymerPropsAd::polymerWaterVelocityRatio(const ADB& c) const; private: - const PolymerProperties polyprops_; + const PolymerProperties polymer_props_; };