From 928513b5619bbebbbdf70a4cda4bd14242d86e75 Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Tue, 28 Aug 2012 16:42:26 +0200 Subject: [PATCH] Added compressible polymer transport solver. Not yet done. --- Makefile.am | 2 + opm/polymer/CompressibleTpfaPolymer.cpp | 14 +- opm/polymer/CompressibleTpfaPolymer.hpp | 1 - .../TransportModelCompressiblePolymer.cpp | 1643 +++++++++++++++++ .../TransportModelCompressiblePolymer.hpp | 206 +++ 5 files changed, 1854 insertions(+), 12 deletions(-) create mode 100644 opm/polymer/TransportModelCompressiblePolymer.cpp create mode 100644 opm/polymer/TransportModelCompressiblePolymer.hpp diff --git a/Makefile.am b/Makefile.am index 5f92f99de..cd97ac5b6 100644 --- a/Makefile.am +++ b/Makefile.am @@ -10,6 +10,7 @@ libopmpolymer_la_SOURCES = \ opm/polymer/IncompTpfaPolymer.cpp \ opm/polymer/SimulatorPolymer.cpp \ opm/polymer/TransportModelPolymer.cpp \ +opm/polymer/TransportModelCompressiblePolymer.cpp \ opm/polymer/PolymerProperties.cpp \ opm/polymer/polymerUtilities.cpp @@ -23,4 +24,5 @@ opm/polymer/PolymerState.hpp \ opm/polymer/SinglePointUpwindTwoPhasePolymer.hpp \ opm/polymer/SimulatorPolymer.hpp \ opm/polymer/TransportModelPolymer.hpp \ +opm/polymer/TransportModelCompressiblePolymer.hpp \ opm/polymer/polymerUtilities.hpp diff --git a/opm/polymer/CompressibleTpfaPolymer.cpp b/opm/polymer/CompressibleTpfaPolymer.cpp index edda93435..4b15a0ede 100644 --- a/opm/polymer/CompressibleTpfaPolymer.cpp +++ b/opm/polymer/CompressibleTpfaPolymer.cpp @@ -102,14 +102,8 @@ namespace Opm const int nc = grid_.number_of_cells; const int np = props_.numPhases(); cell_relperm_.resize(nc*np); - cell_eff_relperm_.resize(nc*np); const double* cell_s = &state.saturation()[0]; props_.relperm(nc, cell_s, &allcells_[0], &cell_relperm_[0], 0); - std::copy(cell_relperm_.begin(), cell_relperm_.end(), cell_eff_relperm_.begin()); - for (int cell; cell < grid_.number_of_cells; ++cell) { - // only the water phase is modified by the presence og polymer. - poly_props_.effectiveRelperm((*c_)[cell], (*cmax_)[cell], &cell_relperm_[nc + 0], cell_eff_relperm_[nc + 0]); - } computeWellPotentials(state); if (rock_comp_props_ && rock_comp_props_->isActive()) { computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), initial_porevol_); @@ -143,14 +137,12 @@ namespace Opm cell_viscosity_.resize(nc*np); props_.viscosity(nc, cell_p, cell_z, &allcells_[0], &cell_viscosity_[0], 0); cell_phasemob_.resize(nc*np); + poly_props_.effective for (int cell; cell < nc; ++cell) { - // Only the water viscosity is modified by polymer. poly_props_.effectiveVisc((*c_)[cell], &cell_viscosity_[nc + 0], cell_eff_viscosity_[nc + 0]); + poly_props_.effectiveMobilities((*c_)[cell], (*cmax_)[cell], &cell_viscosity_[nc + 0], &cell_relperm_[nc + 0], &cell_phasemob_[nc + 0]); } - std::transform(cell_eff_relperm_.begin(), cell_eff_relperm_.end(), - cell_eff_viscosity_.begin(), - cell_phasemob_.begin(), - std::divides()); + // Volume discrepancy: we have that // z = Au, voldiscr = sum(u) - 1, // but I am not sure it is actually needed. diff --git a/opm/polymer/CompressibleTpfaPolymer.hpp b/opm/polymer/CompressibleTpfaPolymer.hpp index 09c2dd4f8..98346f150 100644 --- a/opm/polymer/CompressibleTpfaPolymer.hpp +++ b/opm/polymer/CompressibleTpfaPolymer.hpp @@ -106,7 +106,6 @@ namespace Opm const std::vector* cmax_; std::vector cell_eff_viscosity_; std::vector cell_relperm_; - std::vector cell_eff_relperm_; }; } // namespace Opm diff --git a/opm/polymer/TransportModelCompressiblePolymer.cpp b/opm/polymer/TransportModelCompressiblePolymer.cpp new file mode 100644 index 000000000..8556640ff --- /dev/null +++ b/opm/polymer/TransportModelCompressiblePolymer.cpp @@ -0,0 +1,1643 @@ +/* + Copyright 2012 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 +// Choose error policy for scalar solves here. +typedef Opm::RegulaFalsi RootFinder; + + +class Opm::TransportModelCompressiblePolymer::ResidualEquation +{ +public: + GradientMethod gradient_method; + int cell; + double s0; + double c0; + double cmax0; + double influx; // sum_j min(v_ij, 0)*f(s_j) + double influx_polymer; // sum_j min(v_ij, 0)*f(s_j)*mc(c_j) + double outflux; // sum_j max(v_ij, 0) + double porosity0; + double porosity; + double B_cell0; + double B_cell; + double dtpv; // dt/pv(i) + double dps; + double rhor; + double ads0; + + TransportModelCompressiblePolymer& tm; + + ResidualEquation(TransportModelCompressiblePolymer& tmodel, int cell_index); + void computeResidual(const double* x, double* res) const; + void computeResidual(const double* x, double* res, double& mc, double& ff) const; + double computeResidualS(const double* x) const; + double computeResidualC(const double* x) const; + void computeGradientResS(const double* x, double* res, double* gradient) const; + void computeGradientResC(const double* x, double* res, double* gradient) const; + void computeJacobiRes(const double* x, double* dres_s_dsdc, double* dres_c_dsdc) const; + +private: + void computeResAndJacobi(const double* x, const bool if_res_s, const bool if_res_c, + const bool if_dres_s_dsdc, const bool if_dres_c_dsdc, + double* res, double* dres_s_dsdc, + double* dres_c_dsdc, double& mc, double& ff) const; +}; + +class Opm::TransportModelCompressiblePolymer::ResidualCGrav { +public: + const TransportModelCompressiblePolymer& tm; + const int cell; + const double s0; + const double c0; + const double cmax0; + const double porosity; + const double dtpv; // dt/pv(i) + const double dps; + const double rhor; + double c_ads0; + double gf[2]; + int nbcell[2]; + + ResidualCGrav(const TransportModelCompressiblePolymer& tmodel, + const std::vector& cells, + const int pos, + const double* gravflux); + + double operator()(double c) const; + double computeGravResidualS(double s, double c) const; + double computeGravResidualC(double s, double c) const; +}; + +class Opm::TransportModelCompressiblePolymer::ResidualSGrav { +public: + const ResidualCGrav& res_c_eq_; + double c_; + + ResidualSGrav(const ResidualCGrav& res_c_eq, const double c = 0.0); + double sOfc(double c); + double operator()(double s) const; +}; + + +namespace +{ + bool check_interval(const double* xmin, const double* xmax, double* x); + + double norm(double* res) + { + return std::max(std::abs(res[0]), std::abs(res[1])); + } + + bool solveNewtonStepSC(const double* , const Opm::TransportModelCompressiblePolymer::ResidualEquation&, + const double*, double*); + bool solveNewtonStepC(const double* , const Opm::TransportModelCompressiblePolymer::ResidualEquation&, + const double*, double*); + + + + // Define a piecewise linear curve along which we will look for zero of the "s" or "r" residual. + // The curve starts at "x", goes along the direction "direction" until it hits the boundary of the box of + // admissible values for "s" and "x" (which is given by "[x_min[0], x_max[0]]x[x_min[1], x_max[1]]"). + // Then it joins in a straight line the point "end_point". + class CurveInSCPlane{ + public: + CurveInSCPlane(); + void setup(const double* x, const double* direction, + const double* end_point, const double* x_min, + const double* x_max, const double tol, + double& t_max_out, double& t_out_out); + void computeXOfT(double*, const double) const; + + private: + double direction_[2]; + double end_point_[2]; + double x_max_[2]; + double x_min_[2]; + double t_out_; + double t_max_; // t_max = t_out + 1 + double x_out_[2]; + double x_[2]; + }; + + + // Compute the "s" residual along the curve "curve" for a given residual equation "res_eq". + // The operator() is sent to a root solver. + class ResSOnCurve + { + public: + ResSOnCurve(const Opm::TransportModelCompressiblePolymer::ResidualEquation& res_eq); + double operator()(const double t) const; + CurveInSCPlane curve; + private: + const Opm::TransportModelCompressiblePolymer::ResidualEquation& res_eq_; + }; + + // Compute the "c" residual along the curve "curve" for a given residual equation "res_eq". + // The operator() is sent to a root solver. + class ResCOnCurve + { + public: + ResCOnCurve(const Opm::TransportModelCompressiblePolymer::ResidualEquation& res_eq); + double operator()(const double t) const; + CurveInSCPlane curve; + private: + const Opm::TransportModelCompressiblePolymer::ResidualEquation& res_eq_; + }; + +} + + +namespace Opm +{ + TransportModelCompressiblePolymer::TransportModelCompressiblePolymer(const UnstructuredGrid& grid, + const BlackoilPropertiesInterface& props, + const PolymerProperties& polyprops, + const RockCompressibility& rock_comp, + const SingleCellMethod method, + const double tol, + const int maxit) + : grid_(grid), + porosity_standard_(props.porosity()), + props_(props), + polyprops_(polyprops), + rock_comp_(rock_comp), + darcyflux_(0), + source_(0), + dt_(0.0), + inflow_c_(0.0), + tol_(tol), + maxit_(maxit), + method_(method), + adhoc_safety_(1.0), + concentration_(0), + cmax_(0), + fractionalflow_(grid.number_of_cells, -1.0), + mc_(grid.number_of_cells, -1.0) + { + if (props.numPhases() != 2) { + THROW("Property object must have 2 phases"); + } + int np = props.numPhases(); + int num_cells = props.numCells(); + visc_.resize(np*num_cells); + A_.resize(np*np*num_cells); + A0_.resize(np*np*num_cells); + smin_.resize(np*num_cells); + smax_.resize(np*num_cells); + allcells_.resize(num_cells); + for (int i = 0; i < num_cells; ++i) { + allcells_[i] = i; + } + props.satRange(props.numCells(), &allcells_[0], &smin_[0], &smax_[0]); + } + + + + + void TransportModelCompressiblePolymer::setPreferredMethod(SingleCellMethod method) + { + method_ = method; + } + + + + + void TransportModelCompressiblePolymer::solve(const double* darcyflux, + const std::vector& pressure0, + const std::vector& pressure, + const double* source, + const double dt, + const double inflow_c, + std::vector& saturation, + std::vector& concentration, + std::vector& cmax) + { + darcyflux_ = darcyflux; + source_ = source; + dt_ = dt; + inflow_c_ = inflow_c; + toWaterSat(saturation, saturation_); + concentration_ = &concentration[0]; + cmax_ = &cmax[0]; + +#if PROFILING + res_counts.clear(); +#endif + + props_.viscosity(props_.numCells(), &pressure[0], NULL, &allcells_[0], &visc_[0], NULL); + props_.matrix(props_.numCells(), &pressure0[0], NULL, &allcells_[0], &A0_[0], NULL); + props_.matrix(props_.numCells(), &pressure[0], NULL, &allcells_[0], &A_[0], NULL); + computePorosity(grid_, porosity_standard_, rock_comp_, pressure0, porosity0_); + computePorosity(grid_, porosity_standard_, rock_comp_, pressure, porosity_); + + // Check immiscibility requirement (only done for first cell). + if (A_[1] != 0.0 || A_[2] != 0.0) { + THROW("TransportCompressibleModelCompressibleTwophase requires a property object without miscibility."); + } + std::vector seq(grid_.number_of_cells); + std::vector comp(grid_.number_of_cells + 1); + int ncomp; + compute_sequence_graph(&grid_, darcyflux_, + &seq[0], &comp[0], &ncomp, + &ia_upw_[0], &ja_upw_[0]); + const int nf = grid_.number_of_faces; + std::vector neg_darcyflux(nf); + std::transform(darcyflux, darcyflux + nf, neg_darcyflux.begin(), std::negate()); + compute_sequence_graph(&grid_, &neg_darcyflux[0], + &seq[0], &comp[0], &ncomp, + &ia_downw_[0], &ja_downw_[0]); + reorderAndTransport(grid_, darcyflux); + toBothSat(saturation_, saturation); + + } + + + + + // Residual for saturation equation, single-cell implicit Euler transport + // + // r(s) = s - s0 + dt/pv*( influx + outflux*f(s) ) + // + // where influx is water influx, outflux is total outflux. + // Influxes are negative, outfluxes positive. + struct TransportModelCompressiblePolymer::ResidualS + { + TransportModelCompressiblePolymer::ResidualEquation& res_eq_; + const double c_; + explicit ResidualS(TransportModelCompressiblePolymer::ResidualEquation& res_eq, + const double c) + : res_eq_(res_eq), + c_(c) + { + } + + double operator()(double s) const + { + double x[2]; + x[0] = s; + x[1] = c_; + return res_eq_.computeResidualS(x); + } + }; + + // Residual for concentration equation, single-cell implicit Euler transport + // + // \TODO doc me + // where ... + // Influxes are negative, outfluxes positive. + struct TransportModelCompressiblePolymer::ResidualC + { + mutable double s; // Mutable in order to change it with every operator() call to be the last computed s value. + TransportModelCompressiblePolymer::ResidualEquation& res_eq_; + explicit ResidualC(TransportModelCompressiblePolymer::ResidualEquation& res_eq) + : res_eq_(res_eq) + {} + + void computeBothResiduals(const double s_arg, const double c_arg, double& res_s, double& res_c, double& mc, double& ff) const + { + double x[2]; + double res[2]; + x[0] = s_arg; + x[1] = c_arg; + res_eq_.computeResidual(x, res, mc, ff); + res_s = res[0]; + res_c = res[1]; + } + + double operator()(double c) const + { + ResidualS res_s(res_eq_, c); + int iters_used; + // Solve for s first. + // s = modifiedRegulaFalsi(res_s, std::max(tm.smin_[2*cell], dps), tm.smax_[2*cell], + // tm.maxit_, tm.tol_, iters_used); + s = RootFinder::solve(res_s, res_eq_.s0, 0.0, 1.0, + res_eq_.tm.maxit_, res_eq_.tm.tol_, iters_used); + double x[2]; + x[0] = s; + x[1] = c; + double res = res_eq_.computeResidualC(x); +#ifdef EXTRA_DEBUG_OUTPUT + std::cout << "c = " << c << " s = " << s << " c-residual = " << res << std::endl; +#endif + return res; + } + + double lastSaturation() const + { + return s; + } + }; + + + // ResidualEquation gathers parameters to construct the residual, computes its + // value and the values of its derivatives. + + TransportModelCompressiblePolymer::ResidualEquation::ResidualEquation(TransportModelCompressiblePolymer& tmodel, int cell_index) + : tm(tmodel) + { + gradient_method = Analytic; + cell = cell_index; + const int np = tm.props_.numPhases(); + s0 = tm.saturation_[cell]; + c0 = tm.concentration_[cell]; + cmax0 = tm.cmax_[cell]; + double src_flux = -tm.source_[cell]; + bool src_is_inflow = src_flux < 0.0; + // Not clear why we multiply by B_cell source terms. + influx = src_is_inflow ? B_cell*src_flux : 0.0; + outflux = !src_is_inflow ? B_cell*src_flux : 0.0; + porosity0 = tm.porosity0_[cell]; + porosity = tm.porosity_[cell]; + B_cell0 = 1.0/tm.A0_[np*np*cell + 0]; + B_cell = 1.0/tm.A_[np*np*cell + 0]; + dtpv = tm.dt_/(tm.grid_.cell_volumes[cell]*porosity) ; + dps = tm.polyprops_.deadPoreVol(); + rhor = tm.polyprops_.rockDensity(); + tm.polyprops_.adsorption(c0, cmax0, ads0); + double mc; + tm.computeMc(tm.inflow_c_, mc); + influx_polymer = src_is_inflow ? B_cell*src_flux*mc : 0.0; + for (int i = tm.grid_.cell_facepos[cell]; i < tm.grid_.cell_facepos[cell+1]; ++i) { + int f = tm.grid_.cell_faces[i]; + double flux; + int other; + // Compute cell flux + if (cell == tm.grid_.face_cells[2*f]) { + flux = tm.darcyflux_[f]; + other = tm.grid_.face_cells[2*f+1]; + } else { + flux =-tm.darcyflux_[f]; + other = tm.grid_.face_cells[2*f]; + } + // Add flux to influx or outflux, if interior. + if (other != -1) { + if (flux < 0.0) { + const double b_face =tm.A_[np*np*other+ 0]; + influx += B_cell*b_face*flux*tm.fractionalflow_[other]; + influx_polymer += flux*tm.fractionalflow_[other]*tm.mc_[other]; + } else { + outflux += flux; // Because B_cell*b_face = 1 for outflow faces + } + } + } + } + + + void TransportModelCompressiblePolymer::ResidualEquation::computeResidual(const double* x, double* res) const + { + double dres_s_dsdc[2]; + double dres_c_dsdc[2]; + double mc; + double ff; + computeResAndJacobi(x, true, true, false, false, res, dres_s_dsdc, dres_c_dsdc, mc, ff); + } + + void TransportModelCompressiblePolymer::ResidualEquation::computeResidual(const double* x, double* res, double& mc, double& ff) const + { + double dres_s_dsdc[2]; + double dres_c_dsdc[2]; + computeResAndJacobi(x, true, true, false, false, res, dres_s_dsdc, dres_c_dsdc, mc, ff); + } + + + double TransportModelCompressiblePolymer::ResidualEquation::computeResidualS(const double* x) const + { + double res[2]; + double dres_s_dsdc[2]; + double dres_c_dsdc[2]; + double mc; + double ff; + computeResAndJacobi(x, true, false, false, false, res, dres_s_dsdc, dres_c_dsdc, mc, ff); + return res[0]; + } + + double TransportModelCompressiblePolymer::ResidualEquation::computeResidualC(const double* x) const + { + double res[2]; + double dres_s_dsdc[2]; + double dres_c_dsdc[2]; + double mc; + double ff; + computeResAndJacobi(x, false, true, false, false, res, dres_s_dsdc, dres_c_dsdc, mc, ff); + return res[1]; + } + + void TransportModelCompressiblePolymer::ResidualEquation::computeGradientResS(const double* x, double* res, double* gradient) const + // If gradient_method == FinDif, use finite difference + // If gradient_method == Analytic, use analytic expresions + { + double dres_c_dsdc[2]; + double mc; + double ff; + computeResAndJacobi(x, true, true, true, false, res, gradient, dres_c_dsdc, mc, ff); + } + + void TransportModelCompressiblePolymer::ResidualEquation::computeGradientResC(const double* x, double* res, double* gradient) const + // If gradient_method == FinDif, use finite difference + // If gradient_method == Analytic, use analytic expresions + { + double dres_s_dsdc[2]; + double mc; + double ff; + computeResAndJacobi(x, true, true, false, true, res, dres_s_dsdc, gradient, mc, ff); + } + + // Compute the Jacobian of the residual equations. + void TransportModelCompressiblePolymer::ResidualEquation::computeJacobiRes(const double* x, double* dres_s_dsdc, double* dres_c_dsdc) const + { + double res[2]; + double mc; + double ff; + computeResAndJacobi(x, false, false, true, true, res, dres_s_dsdc, dres_c_dsdc, mc, ff); + } + + void TransportModelCompressiblePolymer::ResidualEquation::computeResAndJacobi(const double* x, const bool if_res_s, const bool if_res_c, + const bool if_dres_s_dsdc, const bool if_dres_c_dsdc, + double* res, double* dres_s_dsdc, + double* dres_c_dsdc, double& mc, double& ff) const + { + if ((if_dres_s_dsdc || if_dres_c_dsdc) && gradient_method == Analytic) { + double s = x[0]; + double c = x[1]; + double dff_dsdc[2]; + double mc_dc; + double ads_dc; + double ads; + tm.fracFlowWithDer(s, c, cmax0, cell, ff, dff_dsdc); + if (if_dres_c_dsdc) { + tm.polyprops_.adsorptionWithDer(c, cmax0, ads, ads_dc); + tm.computeMcWithDer(c, mc, mc_dc); + } else { + tm.polyprops_.adsorption(c, cmax0, ads); + tm.computeMc(c, mc); + } + if (if_res_s) { + res[0] = s - B_cell/B_cell0*porosity0/porosity*s0 + dtpv*(outflux*ff + influx); +#if PROFILING + tm.res_counts.push_back(Newton_Iter(true, cell, x[0], x[1])); +#endif + } + if (if_res_c) { + // Not clear if we the rockcompressibility should be + // considerede as a constant in the adsorption term. + res[1] = (1 - dps)*s*c - (1 - dps)*B_cell/B_cell0*porosity0/porosity*s0*c0 + + rhor*B_cell/porosity*((1.0 - porosity)*ads - (1.0 - porosity0)*ads0) + + dtpv*(outflux*ff*mc + influx_polymer); +#if PROFILING + tm.res_counts.push_back(Newton_Iter(false, cell, x[0], x[1])); +#endif + } + if (if_dres_s_dsdc) { + dres_s_dsdc[0] = 1 + dtpv*outflux*dff_dsdc[0]; + dres_s_dsdc[1] = dtpv*outflux*dff_dsdc[1]; + } + if (if_dres_c_dsdc) { + dres_c_dsdc[0] = (1.0 - dps)*c + dtpv*outflux*dff_dsdc[0]*mc; + dres_c_dsdc[1] = (1 - dps)*s + rhor*B_cell/porosity*(1.0 - porosity)*ads_dc + + dtpv*outflux*(dff_dsdc[1]*mc + ff*mc_dc); + } + + } else if (if_res_c || if_res_s) { + double s = x[0]; + double c = x[1]; + tm.fracFlow(s, c, cmax0, cell, ff); + if (if_res_s) { + res[0] = s - B_cell/B_cell0*porosity0/porosity*s0 + dtpv*(outflux*ff + influx); +#if PROFILING + tm.res_counts.push_back(Newton_Iter(true, cell, x[0], x[1])); +#endif + } + if (if_res_c) { + tm.computeMc(c, mc); + double ads; + tm.polyprops_.adsorption(c, cmax0, ads); + res[1] = (1 - dps)*s*c - (1 - dps)*B_cell/B_cell0*porosity0/porosity*s0*c0 + + rhor*B_cell/porosity*((1.0 - porosity)*ads - (1.0 - porosity0)*ads0) + + dtpv*(outflux*ff*mc + influx_polymer); +#if PROFILING + tm.res_counts.push_back(Newton_Iter(false, cell, x[0], x[1])); +#endif + } + } + + if ((if_dres_c_dsdc || if_dres_s_dsdc) && gradient_method == FinDif) { + double epsi = 1e-8; + double res_epsi[2]; + double res_0[2]; + double x_epsi[2]; + computeResidual(x, res_0); + if (if_dres_s_dsdc) { + x_epsi[0] = x[0] + epsi; + x_epsi[1] = x[1]; + computeResidual(x_epsi, res_epsi); + dres_s_dsdc[0] = (res_epsi[0] - res_0[0])/epsi; + x_epsi[0] = x[0]; + x_epsi[1] = x[1] + epsi; + computeResidual(x_epsi, res_epsi); + dres_s_dsdc[1] = (res_epsi[0] - res_0[0])/epsi; + } + if (if_dres_c_dsdc) { + x_epsi[0] = x[0] + epsi; + x_epsi[1] = x[1]; + computeResidual(x_epsi, res_epsi); + dres_c_dsdc[0] = (res_epsi[1] - res_0[1])/epsi; + x_epsi[0] = x[0]; + x_epsi[1] = x[1] + epsi; + computeResidual(x_epsi, res_epsi); + dres_c_dsdc[1] = (res_epsi[1] - res_0[1])/epsi; + } + } + } + + void TransportModelCompressiblePolymer::solveSingleCell(const int cell) + { + switch (method_) { + case Bracketing: + solveSingleCellBracketing(cell); + break; + case Newton: + solveSingleCellNewton(cell); + break; + case Gradient: + solveSingleCellGradient(cell); + break; + case NewtonSimpleSC: + solveSingleCellNewtonSimple(cell,true); + break; + case NewtonSimpleC: + solveSingleCellNewtonSimple(cell,false); + break; + default: + THROW("Unknown method " << method_); + } + } + + + void TransportModelCompressiblePolymer::solveSingleCellBracketing(int cell) + { + + ResidualEquation res_eq(*this, cell); + ResidualC res(res_eq); + const double a = 0.0; + const double b = polyprops_.cMax()*adhoc_safety_; // Add 10% to account for possible non-monotonicity of hyperbolic system. + int iters_used; + + // Check if current state is an acceptable solution. + double res_sc[2]; + double mc, ff; + res.computeBothResiduals(saturation_[cell], concentration_[cell], res_sc[0], res_sc[1], mc, ff); + if (norm(res_sc) < tol_) { + fractionalflow_[cell] = ff; + mc_[cell] = mc; + return; + } + + concentration_[cell] = RootFinder::solve(res, a, b, maxit_, tol_, iters_used); + cmax_[cell] = std::max(cmax_[cell], concentration_[cell]); + saturation_[cell] = res.lastSaturation(); + fracFlow(saturation_[cell], concentration_[cell], cmax_[cell], cell, + fractionalflow_[cell]); + computeMc(concentration_[cell], mc_[cell]); + } + + + + // Newton method, where we first try a Newton step. Then, if it does not work well, we look for + // the zero of either the residual in s or the residual in c along a specified piecewise linear + // curve. In these cases, we can use a robust 1d solver. + void TransportModelCompressiblePolymer::solveSingleCellGradient(int cell) + { + int iters_used_falsi = 0; + const int max_iters_split = maxit_; + int iters_used_split = 0; + + // Check if current state is an acceptable solution. + ResidualEquation res_eq(*this, cell); + double x[2] = {saturation_[cell], saturation_[cell]*concentration_[cell]}; + double res[2]; + double mc; + double ff; + double x_c[2]; + scToc(x, x_c); + res_eq.computeResidual(x_c, res, mc, ff); + if (norm(res) <= tol_) { + cmax_[cell] = std::max(cmax_[cell], concentration_[cell]); + fractionalflow_[cell] = ff; + mc_[cell] = mc; + return; + } + + double x_min[2] = { 0.0, 0.0 }; + double x_max[2] = { 1.0, polyprops_.cMax()*adhoc_safety_ }; + double x_min_res_s[2] = { x_min[0], x_min[1] }; + double x_max_res_s[2] = { x_max[0], x_min[0] }; + double x_min_res_sc[2] = { x_min[0], x_min[1] }; + double x_max_res_sc[2] = { x_max[0], x_max[1] }; + double t; + double t_max; + double t_out; + double direction[2]; + double end_point[2]; + double gradient[2]; + ResSOnCurve res_s_on_curve(res_eq); + ResCOnCurve res_c_on_curve(res_eq); + bool if_res_s; + + while ((norm(res) > tol_) && (iters_used_split < max_iters_split)) { + if (std::abs(res[0]) < std::abs(res[1])) { + if (res[0] < -tol_) { + direction[0] = x_max_res_s[0] - x[0]; + direction[1] = x_max_res_s[1] - x[1]; + if_res_s = true; + } else if (res[0] > tol_) { + direction[0] = x_min_res_s[0] - x[0]; + direction[1] = x_min_res_s[1] - x[1]; + if_res_s = true; + } else { + scToc(x, x_c); + res_eq.computeGradientResS(x_c, res, gradient); + // dResS/d(s_) = dResS/ds - c/s*dResS/ds + // dResS/d(sc_) = -1/s*dResS/dc + if (x[0] > 1e-2*tol_) { + // With s,c variables, we should have + // direction[0] = -gradient[1]; + // direction[1] = gradient[0]; + // With s, sc variables, we get: + scToc(x, x_c); + direction[0] = 1.0/x[0]*gradient[1]; + direction[1] = gradient[0] - x_c[1]/x[0]*gradient[1]; + } else { + // acceptable approximation for nonlinear relative permeability. + direction[0] = 0.0; + direction[1] = gradient[0]; + } + if_res_s = false; + } + } else { + if (res[1] < -tol_) { + direction[0] = x_max_res_sc[0] - x[0]; + direction[1] = x_max_res_sc[1] - x[1]; + if_res_s = false; + } else if (res[1] > tol_) { + direction[0] = x_min_res_sc[0] - x[0]; + direction[1] = x_min_res_sc[1] - x[1]; + if_res_s = false; + } else { + res_eq.computeGradientResC(x, res, gradient); + // dResC/d(s_) = dResC/ds - c/s*dResC/ds + // dResC/d(sc_) = -1/s*dResC/dc + if (x[0] > 1e-2*tol_) { + // With s,c variables, we should have + // direction[0] = -gradient[1]; + // direction[1] = gradient[0]; + // With s, sc variables, we get: + scToc(x, x_c); + direction[0] = 1.0/x[0]*gradient[1]; + direction[1] = gradient[0] - x_c[1]/x[0]*gradient[1]; + } else { + // We take 1.0/s*gradient[1]: wrong for linear permeability, + // acceptable for nonlinear relative permeability. + direction[0] = 1.0 - res_eq.dps; + direction[1] = gradient[0]; + } + if_res_s = true; + } + } + if (if_res_s) { + if (res[0] < 0) { + end_point[0] = x_max_res_s[0]; + end_point[1] = x_max_res_s[1]; + res_s_on_curve.curve.setup(x, direction, end_point, x_min, x_max, tol_, t_max, t_out); + if (res_s_on_curve(t_out) >= 0) { + t_max = t_out; + } + } else { + end_point[0] = x_min_res_s[0]; + end_point[1] = x_min_res_s[1]; + res_s_on_curve.curve.setup(x, direction, end_point, x_min, x_max, tol_, t_max, t_out); + if (res_s_on_curve(t_out) <= 0) { + t_max = t_out; + } + } + // Note: In some experiments modifiedRegularFalsi does not yield a result under the given tolerance. + t = RootFinder::solve(res_s_on_curve, 0., t_max, maxit_, tol_, iters_used_falsi); + res_s_on_curve.curve.computeXOfT(x, t); + } else { + if (res[1] < 0) { + end_point[0] = x_max_res_sc[0]; + end_point[1] = x_max_res_sc[1]; + res_c_on_curve.curve.setup(x, direction, end_point, x_min, x_max, tol_, t_max, t_out); + if (res_c_on_curve(t_out) >= 0) { + t_max = t_out; + } + } else { + end_point[0] = x_min_res_sc[0]; + end_point[1] = x_min_res_sc[1]; + res_c_on_curve.curve.setup(x, direction, end_point, x_min, x_max, tol_, t_max, t_out); + if (res_c_on_curve(t_out) <= 0) { + t_max = t_out; + } + } + t = RootFinder::solve(res_c_on_curve, 0., t_max, maxit_, tol_, iters_used_falsi); + res_c_on_curve.curve.computeXOfT(x, t); + + } + scToc(x, x_c); + res_eq.computeResidual(x_c, res, mc, ff); + iters_used_split += 1; + } + + + + if ((iters_used_split >= max_iters_split) && (norm(res) > tol_)) { + MESSAGE("Newton for single cell did not work in cell number " << cell); + solveSingleCellBracketing(cell); + } else { + scToc(x, x_c); + concentration_[cell] = x_c[1]; + cmax_[cell] = std::max(cmax_[cell], concentration_[cell]); + saturation_[cell] = x[0]; + fractionalflow_[cell] = ff; + mc_[cell] = mc; + } + } + + void TransportModelCompressiblePolymer::solveSingleCellNewton(int cell) + { + const int max_iters_split = maxit_; + int iters_used_split = 0; + + // Check if current state is an acceptable solution. + ResidualEquation res_eq(*this, cell); + double x[2] = {saturation_[cell], concentration_[cell]}; + double res[2]; + double mc; + double ff; + res_eq.computeResidual(x, res, mc, ff); + if (norm(res) <= tol_) { + cmax_[cell] = std::max(cmax_[cell], concentration_[cell]); + fractionalflow_[cell] = ff; + mc_[cell] = mc; + return; + } else if (0.99 < x[0]) { + // x[0] = 0.5; + // x[1] = polyprops_.cMax()/2.0; + // res_eq.computeResidual(x, res, mc, ff); + } + + const double x_min[2] = { 0.0, 0.0 }; + const double x_max[2] = { 1.0, polyprops_.cMax()*adhoc_safety_ }; + bool successfull_newton_step = true; + // initialize x_new to avoid warning + double x_new[2] = {0.0, 0.0}; + double res_new[2]; + ResSOnCurve res_s_on_curve(res_eq); + ResCOnCurve res_c_on_curve(res_eq); + + // We switch to s-sc variable + x[1] = x[0]*x[1]; + // x_c will contain the s-c variable + double x_c[2]; + + while ((norm(res) > tol_) && + (iters_used_split < max_iters_split) && + successfull_newton_step) { + double dres_s_dsdc[2]; + double dres_c_dsdc[2]; + scToc(x, x_c); + res_eq.computeJacobiRes(x_c, dres_s_dsdc, dres_c_dsdc); + double dFx_dx = (dres_s_dsdc[0]-x_c[1]*dres_s_dsdc[1]); + double dFx_dy; + if (x[0] < 1e-2*tol_) { + dFx_dy = 0.0; + } else { + dFx_dy = (dres_s_dsdc[1]/x[0]); + } + double dFy_dx = (dres_c_dsdc[0]-x_c[1]*dres_c_dsdc[1]); + double dFy_dy; + if (x[0] < 1e-2*tol_) { + dFy_dy = 1.0 - res_eq.dps; + } else { + dFy_dy = (dres_c_dsdc[1]/x[0]); + } + double det = dFx_dx*dFy_dy - dFy_dx*dFx_dy; + double alpha = 1.0; + int max_lin_it = 100; + int lin_it = 0; + res_new[0] = res[0]*2; + res_new[1] = res[1]*2; + while((norm(res_new)>norm(res)) && (lin_it=max_lin_it) { + successfull_newton_step = false; + } else { + scToc(x_new, x_c); + x[0] = x_c[0]; + x[1] = x_c[1]; + res[0] = res_new[0]; + res[1] = res_new[1]; + iters_used_split += 1; + successfull_newton_step = true;; + } + } + + if ((iters_used_split >= max_iters_split) && (norm(res) > tol_)) { + MESSAGE("Newton for single cell did not work in cell number " << cell); + solveSingleCellBracketing(cell); + } else { + concentration_[cell] = x[1]; + cmax_[cell] = std::max(cmax_[cell], concentration_[cell]); + saturation_[cell] = x[0]; + fractionalflow_[cell] = ff; + mc_[cell] = mc; + } + } + + void TransportModelCompressiblePolymer::solveSingleCellNewtonSimple(int cell,bool use_sc) + { + const int max_iters_split = maxit_; + int iters_used_split = 0; + + // Check if current state is an acceptable solution. + ResidualEquation res_eq(*this, cell); + double x[2] = {saturation_[cell], concentration_[cell]}; + double res[2]; + double mc; + double ff; + res_eq.computeResidual(x, res, mc, ff); + if (norm(res) <= tol_) { + cmax_[cell] = std::max(cmax_[cell], concentration_[cell]); + fractionalflow_[cell] = ff; + mc_[cell] = mc; + return; + }else{ + //* + x[0] = saturation_[cell]-res[0]; + if((x[0]>1) || (x[0]<0)){ + x[0] = 0.5; + x[1] = x[1]; + } + if(x[0]>0){ + x[1] = concentration_[cell]*saturation_[cell]-res[1]; + x[1] = x[1]/x[0]; + if(x[1]> polyprops_.cMax()){ + x[1]= polyprops_.cMax()/2.0; + } + if(x[1]<0){ + x[1]=0; + } + }else{ + x[1]=0; + } + //x[0]=0.5;x[1]=polyprops_.cMax()/2.0; + res_eq.computeResidual(x, res, mc, ff); + //*/ + } + + + // double x_min[2] = { std::max(polyprops_.deadPoreVol(), smin_[2*cell]), 0.0 }; + double x_min[2] = { 0.0, 0.0 }; + double x_max[2] = { 1.0, polyprops_.cMax()*adhoc_safety_ }; + bool successfull_newton_step = true; + double x_new[2]; + double res_new[2]; + ResSOnCurve res_s_on_curve(res_eq); + ResCOnCurve res_c_on_curve(res_eq); + + while ((norm(res) > tol_) && + (iters_used_split < max_iters_split) && + successfull_newton_step) { + // We first try a Newton step + double dres_s_dsdc[2]; + double dres_c_dsdc[2]; + double dx=tol_; + double tmp_x[2]; + if(!(x[0]>0)){ + tmp_x[0]=dx; + tmp_x[1]=0; + }else{ + tmp_x[0]=x[0]; + tmp_x[1]=x[1]; + } + res_eq.computeJacobiRes(tmp_x, dres_s_dsdc, dres_c_dsdc); + double dFx_dx,dFx_dy,dFy_dx,dFy_dy; + double det = dFx_dx*dFy_dy - dFy_dx*dFx_dy; + if(use_sc){ + dFx_dx=(dres_s_dsdc[0]-tmp_x[1]*dres_s_dsdc[1]/tmp_x[0]); + dFx_dy= (dres_s_dsdc[1]/tmp_x[0]); + dFy_dx=(dres_c_dsdc[0]-tmp_x[1]*dres_c_dsdc[1]/tmp_x[0]); + dFy_dy= (dres_c_dsdc[1]/tmp_x[0]); + }else{ + dFx_dx= dres_s_dsdc[0]; + dFx_dy= dres_s_dsdc[1]; + dFy_dx= dres_c_dsdc[0]; + dFy_dy= dres_c_dsdc[1]; + } + det = dFx_dx*dFy_dy - dFy_dx*dFx_dy; + if(det==0){ + THROW("det is zero"); + } + + double alpha=1; + int max_lin_it=10; + int lin_it=0; + /* + std::cout << "Nonlinear " << iters_used_split + << " " << norm(res_new) + << " " << norm(res) << std::endl;*/ + /* + std::cout << "x" << std::endl; + std::cout << " " << x[0] << " " << x[1] << std::endl; + std::cout << "dF" << std::endl; + std::cout << " " << dFx_dx << " " << dFx_dy << std::endl; + std::cout << " " << dFy_dx << " " << dFy_dy << std::endl; + std::cout << "dres" << std::endl; + std::cout << " " << dres_s_dsdc[0] << " " << dres_s_dsdc[1] << std::endl; + std::cout << " " << dres_c_dsdc[0] << " " << dres_c_dsdc[1] << std::endl; + std::cout << std::endl; + */ + res_new[0]=res[0]*2; + res_new[1]=res[1]*2; + double update[2]={(res[0]*dFy_dy - res[1]*dFx_dy)/det, + (res[1]*dFx_dx - res[0]*dFy_dx)/det}; + while((norm(res_new)>norm(res)) && (lin_it0) ? x_new[1]/x_new[0] : 0.0; + }else{ + x_new[1] = x[1] - alpha*update[1]; + } + check_interval(x_min, x_max, x_new); + res_eq.computeResidual(x_new, res_new, mc, ff); + // std::cout << " " << res_new[0] << " " << res_new[1] << std::endl; + alpha=alpha/2.0; + lin_it=lin_it+1; + // std::cout << "Linear iterations" << lin_it << " " << norm(res_new) << std::endl; + } + if (lin_it>=max_lin_it) { + successfull_newton_step = false; + } else { + x[0] = x_new[0]; + x[1] = x_new[1]; + res[0] = res_new[0]; + res[1] = res_new[1]; + iters_used_split += 1; + successfull_newton_step = true;; + } + // std::cout << "Nonlinear " << iters_used_split << " " << norm(res) << std::endl; + } + + if ((iters_used_split >= max_iters_split) || (norm(res) > tol_)) { + MESSAGE("NewtonSimple for single cell did not work in cell number " << cell); + solveSingleCellBracketing(cell); + } else { + concentration_[cell] = x[1]; + cmax_[cell] = std::max(cmax_[cell], concentration_[cell]); + saturation_[cell] = x[0]; + fractionalflow_[cell] = ff; + mc_[cell] = mc; + } + } + + + + void TransportModelCompressiblePolymer::solveMultiCell(const int num_cells, const int* cells) + { + double max_s_change = 0.0; + double max_c_change = 0.0; + int num_iters = 0; + // Must store state variables before we start. + std::vector s0(num_cells); + std::vector c0(num_cells); + std::vector cmax0(num_cells); + // Must set initial fractional flows etc. before we start. + for (int i = 0; i < num_cells; ++i) { + const int cell = cells[i]; + fracFlow(saturation_[cell], concentration_[cell], cmax_[cell], + cell, fractionalflow_[cell]); + computeMc(concentration_[cell], mc_[cell]); + s0[i] = saturation_[cell]; + c0[i] = concentration_[cell]; + cmax0[i] = cmax_[i]; + } + do { + // int max_s_change_cell = -1; + // int max_c_change_cell = -1; + max_s_change = 0.0; + max_c_change = 0.0; + for (int i = 0; i < num_cells; ++i) { + const int cell = cells[i]; + const double old_s = saturation_[cell]; + const double old_c = concentration_[cell]; + saturation_[cell] = s0[i]; + concentration_[cell] = c0[i]; + cmax_[cell] = cmax0[i]; + solveSingleCell(cell); + // std::cout << "cell = " << cell << " delta s = " << saturation_[cell] - old_s << std::endl; + // if (max_s_change < std::fabs(saturation_[cell] - old_s)) { + // max_s_change_cell = cell; + // } + // if (max_c_change < std::fabs(concentration_[cell] - old_c)) { + // max_c_change_cell = cell; + // } + max_s_change = std::max(max_s_change, std::fabs(saturation_[cell] - old_s)); + max_c_change = std::max(max_c_change, std::fabs(concentration_[cell] - old_c)); + } + // std::cout << "Iter = " << num_iters << " max_s_change = " << max_s_change + // << " in cell " << max_change_cell << std::endl; + } while (((max_s_change > tol_) || (max_c_change > tol_)) && ++num_iters < maxit_); + if (max_s_change > tol_) { + THROW("In solveMultiCell(), we did not converge after " + << num_iters << " iterations. Delta s = " << max_s_change); + } + if (max_c_change > tol_) { + THROW("In solveMultiCell(), we did not converge after " + << num_iters << " iterations. Delta c = " << max_c_change); + } + std::cout << "Solved " << num_cells << " cell multicell problem in " + << num_iters << " iterations." << std::endl; + } + + void TransportModelCompressiblePolymer::fracFlow(double s, double c, double cmax, + int cell, double& ff) const + { + double dummy[2]; + fracFlowBoth(s, c, cmax, cell, ff, dummy, false); + } + + void TransportModelCompressiblePolymer::fracFlowWithDer(double s, double c, double cmax, + int cell, double& ff, + double* dff_dsdc) const + { + fracFlowBoth(s, c, cmax, cell, ff, dff_dsdc, true); + } + + void TransportModelCompressiblePolymer::fracFlowBoth(double s, double c, double cmax, int cell, + double& ff, double* dff_dsdc, + bool if_with_der) const + { + double relperm[2]; + double drelperm_ds[4]; + double sat[2] = {s, 1 - s}; + if (if_with_der) { + props_.relperm(1, sat, &cell, relperm, drelperm_ds); + } else { + props_.relperm(1, sat, &cell, relperm, 0); + } + double mob[2]; + double dmob_ds[4]; + double dmob_dc[2]; + double dmobwat_dc; + polyprops_.effectiveMobilitiesBoth(c, cmax, &visc_[cell], relperm, drelperm_ds, + mob, dmob_ds, dmobwat_dc, if_with_der); + + ff = mob[0]/(mob[0] + mob[1]); + if (if_with_der) { + dmob_dc[0] = dmobwat_dc; + dmob_dc[1] = 0.; + //dff_dsdc[0] = (dmob_ds[0]*mob[1] + dmob_ds[3]*mob[0])/((mob[0] + mob[1])*(mob[0] + mob[1])); // derivative with respect to s + // at the moment the dmob_ds only have diagonal elements since the saturation is derivated out in effectiveMobilitiesBouth + dff_dsdc[0] = ((dmob_ds[0]-dmob_ds[2])*mob[1] - (dmob_ds[1]-dmob_ds[3])*mob[0])/((mob[0] + mob[1])*(mob[0] + mob[1])); // derivative with respect to s + dff_dsdc[1] = (dmob_dc[0]*mob[1] - dmob_dc[1]*mob[0])/((mob[0] + mob[1])*(mob[0] + mob[1])); // derivative with respect to c + } + } + + void TransportModelCompressiblePolymer::computeMc(double c, double& mc) const + { + polyprops_.computeMc(c, mc); + } + + void TransportModelCompressiblePolymer::computeMcWithDer(double c, double& mc, + double &dmc_dc) const + { + polyprops_.computeMcWithDer(c, mc, dmc_dc); + } + + + + TransportModelCompressiblePolymer::ResidualSGrav::ResidualSGrav(const ResidualCGrav& res_c_eq, + const double c) + : res_c_eq_(res_c_eq), + c_(c) + { + } + + double TransportModelCompressiblePolymer::ResidualSGrav::operator()(double s) const + { + return res_c_eq_.computeGravResidualS(s, c_); + } + + double TransportModelCompressiblePolymer::ResidualSGrav::sOfc(double c) // for a given c, returns the value of s for which residual in s vanishes. + { + c_ = c; + int iters_used; + return RootFinder::solve(*this, res_c_eq_.s0, 0.0, 1.0, + res_c_eq_.tm.maxit_, res_c_eq_.tm.tol_, + iters_used); + } + + + + // Residual for concentration equation for gravity segregation + // + // res_c = s*(1 - dps)*c - s0*( - dps)*c0 + // + dtpv*sum_{j adj i}( mc * gravmod_ij * gf_ij ). + // \TODO doc me + // where ... + // Influxes are negative, outfluxes positive. + + TransportModelCompressiblePolymer::ResidualCGrav::ResidualCGrav(const TransportModelCompressiblePolymer& tmodel, + const std::vector& cells, + const int pos, + const double* gravflux) // Always oriented towards next in column. Size = colsize - 1. + : tm(tmodel), + cell(cells[pos]), + s0(tm.saturation_[cell]), + c0(tm.concentration_[cell]), + cmax0(tm.cmax0_[cell]), + porosity(tm.porosity_[cell]), + dtpv(tm.dt_/tm.porevolume_[cell]), + dps(tm.polyprops_.deadPoreVol()), + rhor(tm.polyprops_.rockDensity()) + + { + nbcell[0] = -1; + gf[0] = 0.0; + if (pos > 0) { + nbcell[0] = cells[pos - 1]; + gf[0] = -gravflux[pos - 1]; + } + nbcell[1] = -1; + gf[1] = 0.0; + if (pos < int(cells.size() - 1)) { + nbcell[1] = cells[pos + 1]; + gf[1] = gravflux[pos]; + } + + tm.polyprops_.adsorption(c0, cmax0, c_ads0); + } + + double TransportModelCompressiblePolymer::ResidualCGrav::operator()(double c) const + { + + ResidualSGrav res_s(*this); + return computeGravResidualC(res_s.sOfc(c), c); + + } + + double TransportModelCompressiblePolymer::ResidualCGrav::computeGravResidualS(double s, double c) const + { + + double mobcell[2]; + tm.mobility(s, c, cell, mobcell); + + double res = s - s0; + + for (int nb = 0; nb < 2; ++nb) { + if (nbcell[nb] != -1) { + double m[2]; + if (gf[nb] < 0.0) { + m[0] = mobcell[0]; + m[1] = tm.mob_[2*nbcell[nb] + 1]; + } else { + m[0] = tm.mob_[2*nbcell[nb]]; + m[1] = mobcell[1]; + } + if (m[0] + m[1] > 0.0) { + res += -dtpv*gf[nb]*m[0]*m[1]/(m[0] + m[1]); + } + } + } + return res; + } + + double TransportModelCompressiblePolymer::ResidualCGrav::computeGravResidualC(double s, double c) const + { + + double mobcell[2]; + tm.mobility(s, c, cell, mobcell); + double c_ads; + tm.polyprops_.adsorption(c, cmax0, c_ads); + + double res = (1 - dps)*s*c - (1 - dps)*s0*c0 + + rhor*((1.0 - porosity)/porosity)*(c_ads - c_ads0); + for (int nb = 0; nb < 2; ++nb) { + if (nbcell[nb] != -1) { + double m[2]; + double mc; + if (gf[nb] < 0.0) { + m[0] = mobcell[0]; + tm.computeMc(c, mc); + m[1] = tm.mob_[2*nbcell[nb] + 1]; + } else { + m[0] = tm.mob_[2*nbcell[nb]]; + mc = tm.mc_[nbcell[nb]]; + m[1] = mobcell[1]; + } + if (m[0] + m[1] > 0.0) { + res += -dtpv*gf[nb]*mc*m[0]*m[1]/(m[0] + m[1]); + } + } + } + return res; + } + + + void TransportModelCompressiblePolymer::mobility(double s, double c, int cell, double* mob) const + { + double sat[2] = { s, 1.0 - s }; + double relperm[2]; + props_.relperm(1, sat, &cell, relperm, 0); + polyprops_.effectiveMobilities(c, cmax0_[cell], visc_, relperm, mob); + } + + + void TransportModelCompressiblePolymer::initGravity(const double* grav) + { + // Set up gravflux_ = T_ij g (rho_w - rho_o) (z_i - z_j) + std::vector htrans(grid_.cell_facepos[grid_.number_of_cells]); + const int nf = grid_.number_of_faces; + const int dim = grid_.dimensions; + gravflux_.resize(nf); + tpfa_htrans_compute(const_cast(&grid_), props_.permeability(), &htrans[0]); + tpfa_trans_compute(const_cast(&grid_), &htrans[0], &gravflux_[0]); + const double delta_rho = props_.density()[0] - props_.density()[1]; + for (int f = 0; f < nf; ++f) { + const int* c = &grid_.face_cells[2*f]; + double gdz = 0.0; + if (c[0] != -1 && c[1] != -1) { + for (int d = 0; d < dim; ++d) { + gdz += grav[d]*(grid_.cell_centroids[dim*c[0] + d] - grid_.cell_centroids[dim*c[1] + d]); + } + } + gravflux_[f] *= delta_rho*gdz; + } + } + + + void TransportModelCompressiblePolymer::solveSingleCellGravity(const std::vector& cells, + const int pos, + const double* gravflux) + { + const int cell = cells[pos]; + ResidualCGrav res_c(*this, cells, pos, gravflux); + + // Check if current state is an acceptable solution. + double res_sc[2]; + res_sc[0]=res_c.computeGravResidualS(saturation_[cell], concentration_[cell]); + res_sc[1]=res_c.computeGravResidualC(saturation_[cell], concentration_[cell]); + + if (norm(res_sc) < tol_) { + return; + } + + const double a = 0.0; + const double b = polyprops_.cMax()*adhoc_safety_; // Add 10% to account for possible non-monotonicity of hyperbolic system. + int iters_used; + concentration_[cell] = RootFinder::solve(res_c, a, b, maxit_, tol_, iters_used); + ResidualSGrav res_s(res_c); + saturation_[cell] = res_s.sOfc(concentration_[cell]); + cmax_[cell] = std::max(cmax0_[cell], concentration_[cell]); + saturation_[cell] = std::min(std::max(saturation_[cell], smin_[2*cell]), smax_[2*cell]); + computeMc(concentration_[cell], mc_[cell]); + mobility(saturation_[cell], concentration_[cell], cell, &mob_[2*cell]); + } + + int TransportModelCompressiblePolymer::solveGravityColumn(const std::vector& cells) + { + // Set up column gravflux. + const int nc = cells.size(); + std::vector col_gravflux(nc - 1); + for (int ci = 0; ci < nc - 1; ++ci) { + const int cell = cells[ci]; + const int next_cell = cells[ci + 1]; + for (int j = grid_.cell_facepos[cell]; j < grid_.cell_facepos[cell+1]; ++j) { + const int face = grid_.cell_faces[j]; + const int c1 = grid_.face_cells[2*face + 0]; + const int c2 = grid_.face_cells[2*face + 1]; + if (c1 == next_cell || c2 == next_cell) { + const double gf = gravflux_[face]; + col_gravflux[ci] = (c1 == cell) ? gf : -gf; + } + } + } + + // Store initial saturation s0 + s0_.resize(nc); + c0_.resize(nc); + for (int ci = 0; ci < nc; ++ci) { + s0_[ci] = saturation_[cells[ci]]; + c0_[ci] = concentration_[cells[ci]]; + } + + // Solve single cell problems, repeating if necessary. + double max_sc_change = 0.0; + int num_iters = 0; + do { + max_sc_change = 0.0; + for (int ci = 0; ci < nc; ++ci) { + const int ci2 = nc - ci - 1; + double old_s[2] = { saturation_[cells[ci]], + saturation_[cells[ci2]] }; + double old_c[2] = { concentration_[cells[ci]], + concentration_[cells[ci2]] }; + saturation_[cells[ci]] = s0_[ci]; + concentration_[cells[ci]] = c0_[ci]; + solveSingleCellGravity(cells, ci, &col_gravflux[0]); + saturation_[cells[ci2]] = s0_[ci2]; + concentration_[cells[ci2]] = c0_[ci2]; + solveSingleCellGravity(cells, ci2, &col_gravflux[0]); + max_sc_change = std::max(max_sc_change, 0.25*(std::fabs(saturation_[cells[ci]] - old_s[0]) + + std::fabs(concentration_[cells[ci]] - old_c[0]) + + std::fabs(saturation_[cells[ci2]] - old_s[1]) + + std::fabs(concentration_[cells[ci2]] - old_c[1]))); + } + // std::cout << "Iter = " << num_iters << " max_s_change = " << max_s_change << std::endl; + } while (max_sc_change > tol_ && ++num_iters < maxit_); + + if (max_sc_change > tol_) { + THROW("In solveGravityColumn(), we did not converge after " + << num_iters << " iterations. Delta s = " << max_sc_change); + } + return num_iters + 1; + } + + + void TransportModelCompressiblePolymer::solveGravity(const std::vector >& columns, + const double* porevolume, + const double dt, + std::vector& saturation, + std::vector& concentration, + std::vector& cmax) + { + // initialize variables. + porevolume_ = porevolume; + dt_ = dt; + toWaterSat(saturation, saturation_); + concentration_ = &concentration[0]; + cmax_ = &cmax[0]; + const int nc = grid_.number_of_cells; + cmax0_.resize(nc); + std::copy(cmax.begin(), cmax.end(), &cmax0_[0]); + + // Initialize mobilities. + mob_.resize(2*nc); + + for (int cell = 0; cell < nc; ++cell) { + mobility(saturation_[cell], concentration_[cell], cell, &mob_[2*cell]); + } + + + // Solve on all columns. + int num_iters = 0; + // std::cout << "Gauss-Seidel column solver # columns: " << columns.size() << std::endl; + for (std::vector >::size_type i = 0; i < columns.size(); i++) { + // std::cout << "==== new column" << std::endl; + num_iters += solveGravityColumn(columns[i]); + } + std::cout << "Gauss-Seidel column solver average iterations: " + << double(num_iters)/double(columns.size()) << std::endl; + + toBothSat(saturation_, saturation); + } + + void TransportModelCompressiblePolymer::scToc(const double* x, double* x_c) const { + x_c[0] = x[0]; + if (x[0] < 1e-2*tol_) { + x_c[1] = 0.5*polyprops_.cMax(); + } else { + x_c[1] = x[1]/x[0]; + } + } + +} // namespace Opm + + +namespace +{ + bool check_interval(const double* xmin, const double* xmax, double* x) { + bool test = false; + if (x[0] < xmin[0]) { + test = true; + x[0] = xmin[0]; + } else if (x[0] > xmax[0]) { + test = true; + x[0] = xmax[0]; + } + if (x[1] < xmin[1]) { + test = true; + x[1] = xmin[1]; + } else if (x[1] > xmax[1]) { + test = true; + x[1] = xmax[1]; + } + return test; + } + + + CurveInSCPlane::CurveInSCPlane() + { + } + + // Setup the curve (see comment above). + // The curve is parametrized by t in [0, t_max], t_out is equal to t when the curve hits the bounding + // rectangle. x_out=(s_out, c_out) denotes the values of s and c at that point. + void CurveInSCPlane::setup(const double* x, const double* direction, + const double* end_point, const double* x_min, + const double* x_max, const double tol, + double& t_max_out, double& t_out_out) + { + x_[0] = x[0]; + x_[1] = x[1]; + x_max_[0] = x_max[0]; + x_max_[1] = x_max[1]; + x_min_[0] = x_min[0]; + x_min_[1] = x_min[1]; + direction_[0] = direction[0]; + direction_[1] = direction[1]; + end_point_[0] = end_point[0]; + end_point_[1] = end_point[1]; + const double size_direction = std::abs(direction_[0]) + std::abs(direction_[1]); + if (size_direction < tol) { + direction_[0] = end_point_[0]-x_[0]; + direction_[1] = end_point_[1]-x_[1]; + } else if ((end_point_[0]-x_[0])*direction_[0] + (end_point_[1]-x_[1])*direction_[1] < 0) { + direction_[0] *= -1.0; + direction_[1] *= -1.0; + } + bool t0_exists = true; + double t0 = 0; // dummy default value (so that compiler does not complain). + if (direction_[0] > 0) { + t0 = (x_max_[0] - x_[0])/direction_[0]; + } else if (direction_[0] < 0) { + t0 = (x_min_[0] - x_[0])/direction_[0]; + } else { + t0_exists = false; + } + bool t1_exists = true; + double t1 = 0; // dummy default value. + if (direction_[1] > 0) { + t1 = (x_max_[1] - x_[1])/direction_[1]; + } else if (direction_[1] < 0) { + t1 = (x_min_[1] - x_[1])/direction_[1]; + } else { + t1_exists = false; + } + if (t0_exists) { + if (t1_exists) { + t_out_ = std::min(t0, t1); + } else { + t_out_ = t0; + } + } else if (t1_exists) { + t_out_ = t1; + } else { + THROW("Direction illegal: is a zero vector."); + } + x_out_[0] = x_[0] + t_out_*direction_[0]; + x_out_[1] = x_[1] + t_out_*direction_[1]; + t_max_ = t_out_ + 1; + t_max_out = t_max_; + t_out_out = t_out_; + } + + + // Compute x=(s,c) for a given t (t is the parameter for the piecewise linear curve) + void CurveInSCPlane::computeXOfT(double* x_of_t, const double t) const { + if (t <= t_out_) { + x_of_t[0] = x_[0] + t*direction_[0]; + x_of_t[1] = x_[1] + t*direction_[1]; + } else { + x_of_t[0] = 1/(t_max_-t_out_)*((t_max_ - t)*x_out_[0] + end_point_[0]*(t - t_out_)); + x_of_t[1] = 1/(t_max_-t_out_)*((t_max_ - t)*x_out_[1] + end_point_[1]*(t - t_out_)); + } + } + + + ResSOnCurve::ResSOnCurve(const Opm::TransportModelCompressiblePolymer::ResidualEquation& res_eq) + : res_eq_(res_eq) + { + } + + double ResSOnCurve::operator()(const double t) const + { + double x_of_t[2]; + double x_c[2]; + curve.computeXOfT(x_of_t, t); + res_eq_.tm.scToc(x_of_t, x_c); + return res_eq_.computeResidualS(x_c); + } + + ResCOnCurve::ResCOnCurve(const Opm::TransportModelCompressiblePolymer::ResidualEquation& res_eq) + : res_eq_(res_eq) + { + } + + double ResCOnCurve::operator()(const double t) const + { + double x_of_t[2]; + double x_c[2]; + curve.computeXOfT(x_of_t, t); + res_eq_.tm.scToc(x_of_t, x_c); + return res_eq_.computeResidualC(x_c); + } + + bool solveNewtonStepSC(const double* xx, const Opm::TransportModelCompressiblePolymer::ResidualEquation& res_eq, + const double* res, double* x_new) { + + double dres_s_dsdc[2]; + double dres_c_dsdc[2]; + double dx=1e-8; + double x[2]; + if(!(x[0]>0)){ + x[0] = dx; + x[1] = 0; + } else { + x[0] = xx[0]; + x[1] = xx[1]; + } + res_eq.computeJacobiRes(x, dres_s_dsdc, dres_c_dsdc); + + double dFx_dx=(dres_s_dsdc[0]-x[1]*dres_s_dsdc[1]); + double dFx_dy= (dres_s_dsdc[1]/x[0]); + double dFy_dx=(dres_c_dsdc[0]-x[1]*dres_c_dsdc[1]); + double dFy_dy= (dres_c_dsdc[1]/x[0]); + + double det = dFx_dx*dFy_dy - dFy_dx*dFx_dy; + if (std::abs(det) < 1e-8) { + return false; + } else { + x_new[0] = x[0] - (res[0]*dFy_dy - res[1]*dFx_dy)/det; + x_new[1] = x[1]*x[0] - (res[1]*dFx_dx - res[0]*dFy_dx)/det; + x_new[1] = (x_new[0]>0) ? x_new[1]/x_new[0] : 0.0; + return true; + } + } + bool solveNewtonStepC(const double* xx, const Opm::TransportModelCompressiblePolymer::ResidualEquation& res_eq, + const double* res, double* x_new) { + + double dres_s_dsdc[2]; + double dres_c_dsdc[2]; + double dx=1e-8; + double x[2]; + if(!(x[0]>0)){ + x[0] = dx; + x[1] = 0; + } else { + x[0] = xx[0]; + x[1] = xx[1]; + } + res_eq.computeJacobiRes(x, dres_s_dsdc, dres_c_dsdc); + double det = dres_s_dsdc[0]*dres_c_dsdc[1] - dres_c_dsdc[0]*dres_s_dsdc[1]; + if (std::abs(det) < 1e-8) { + return false; + } else { + x_new[0] = x[0] - (res[0]*dres_c_dsdc[1] - res[1]*dres_s_dsdc[1])/det; + x_new[1] = x[1] - (res[1]*dres_s_dsdc[0] - res[0]*dres_c_dsdc[0])/det; + return true; + } + } + + +} // Anonymous namespace + + +/* Local Variables: */ +/* c-basic-offset:4 */ +/* End: */ diff --git a/opm/polymer/TransportModelCompressiblePolymer.hpp b/opm/polymer/TransportModelCompressiblePolymer.hpp new file mode 100644 index 000000000..d88a46e2a --- /dev/null +++ b/opm/polymer/TransportModelCompressiblePolymer.hpp @@ -0,0 +1,206 @@ +/* + Copyright 2012 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_TRANSPORTMODELPOLYMER_HEADER_INCLUDED +#define OPM_TRANSPORTMODELPOLYMER_HEADER_INCLUDED + +#include +#include +#include +#include +#include +#include + +class UnstructuredGrid; + +namespace Opm +{ + + class BlackoilPropertiesInterface; + + /// Implements a reordering transport solver for incompressible two-phase flow + /// with polymer in the water phase. + /// \TODO Include permeability reduction effect. + class TransportModelCompressiblePolymer : public TransportModelInterface + { + public: + + enum SingleCellMethod { Bracketing, Newton, Gradient, NewtonSimpleSC, NewtonSimpleC}; + enum GradientMethod { Analytic, FinDif }; // Analytic is chosen (hard-coded) + + /// Construct solver. + /// \param[in] grid A 2d or 3d grid. + /// \param[in] props Rock and fluid properties. + /// \param[in] polyprops Polymer properties. + /// \param[in] rock_comp Rock compressibility properties + /// \param[in] method Bracketing: solve for c in outer loop, s in inner loop, + /// each solve being bracketed for robustness. + /// Newton: solve simultaneously for c and s with Newton's method. + /// (using gradient variant and bracketing as fallbacks). + /// \param[in] tol Tolerance used in the solver. + /// \param[in] maxit Maximum number of non-linear iterations used. + TransportModelCompressiblePolymer(const UnstructuredGrid& grid, + const BlackoilPropertiesInterface& props, + const PolymerProperties& polyprops, + const RockCompressibility& rock_comp, + const SingleCellMethod method, + const double tol, + const int maxit); + + /// Set the preferred method, Bracketing or Newton. + void setPreferredMethod(SingleCellMethod method); + + /// Solve for saturation, concentration and cmax at next timestep. + /// Using implicit Euler scheme, reordered. + /// \param[in] darcyflux Array of signed face fluxes. + /// \param[in] pressure0 Array with pressure at start of timestep. + /// \param[in] pressure Array with pressure. + /// \param[in] source Transport source term. + /// \param[in] dt Time step. + /// \param[in] inflow_c Inflow polymer. + /// \param[in, out] saturation Phase saturations. + /// \param[in, out] concentration Polymer concentration. + /// \param[in, out] cmax Highest concentration that has occured in a given cell. + void solve(const double* darcyflux, + const std::vector& pressure0, + const std::vector& pressure, + const double* source, + const double dt, + const double inflow_c, + std::vector& saturation, + std::vector& concentration, + std::vector& cmax); + + /// Solve for gravity segregation. + /// This uses a column-wise nonlinear Gauss-Seidel approach. + /// It assumes that the input columns contain cells in a single + /// vertical stack, that do not interact with other columns (for + /// gravity segregation. + /// \param[in] columns Vector of cell-columns. + /// \param[in] porevolume Array of pore volumes. + /// \param[in] dt Time step. + /// \param[in, out] saturation Phase saturations. + /// \param[in, out] concentration Polymer concentration. + /// \param[in, out] cmax Highest concentration that has occured in a given cell. + void solveGravity(const std::vector >& columns, + const double* porevolume, + const double dt, + std::vector& saturation, + std::vector& concentration, + std::vector& cmax); + + public: // But should be made private... + virtual void solveSingleCell(const int cell); + virtual void solveMultiCell(const int num_cells, const int* cells); + void solveSingleCellBracketing(int cell); + void solveSingleCellNewton(int cell); + void solveSingleCellGradient(int cell); + void solveSingleCellNewtonSimple(int cell,bool use_sc); + class ResidualEquation; + + void initGravity(const double* grav); + void solveSingleCellGravity(const std::vector& cells, + const int pos, + const double* gravflux); + int solveGravityColumn(const std::vector& cells); + void scToc(const double* x, double* x_c) const; + + #ifdef PROFILING + class Newton_Iter { + public: + bool res_s; + int cell; + double s; + double c; + + Newton_Iter(bool res_s_val, int cell_val, double s_val, double c_val) { + res_s = res_s_val; + cell = cell_val; + s = s_val; + c = c_val; + } + }; + + std::list res_counts; + #endif + + + private: + const UnstructuredGrid& grid_; + const double* porosity_standard_; + const BlackoilPropertiesInterface& props_; + const PolymerProperties& polyprops_; + const RockCompressibility& rock_comp_; + const double* darcyflux_; // one flux per grid face + const double* source_; // one source per cell + double dt_; + double inflow_c_; + double tol_; + double maxit_; + SingleCellMethod method_; + double adhoc_safety_; + + std::vector saturation_; // one per cell, only water saturation! + std::vector allcells_; + double* concentration_; + double* cmax_; + std::vector fractionalflow_; // one per cell + std::vector mc_; // one per cell + std::vector visc_; + std::vector A_; + std::vector A0_; + std::vector porosity0_; + std::vector porosity_; + std::vector smin_; + std::vector smax_; + + // For gravity segregation. + std::vector gravflux_; + std::vector mob_; + std::vector cmax0_; + // For gravity segregation, column variables + std::vector s0_; + std::vector c0_; + + // Storing the upwind and downwind graphs for experiments. + std::vector ia_upw_; + std::vector ja_upw_; + std::vector ia_downw_; + std::vector ja_downw_; + + struct ResidualC; + struct ResidualS; + + class ResidualCGrav; + class ResidualSGrav; + + + void fracFlow(double s, double c, double cmax, int cell, double& ff) const; + void fracFlowWithDer(double s, double c, double cmax, int cell, double& ff, + double* dff_dsdc) const; + void fracFlowBoth(double s, double c, double cmax, int cell, double& ff, + double* dff_dsdc, bool if_with_der) const; + void computeMc(double c, double& mc) const; + void computeMcWithDer(double c, double& mc, double& dmc_dc) const; + void mobility(double s, double c, int cell, double* mob) const; + }; + +} // namespace Opm + +#endif // OPM_TRANSPORTMODELPOLYMER_HEADER_INCLUDED