diff --git a/Makefile.am b/Makefile.am index cb7a5633..845c5fa9 100644 --- a/Makefile.am +++ b/Makefile.am @@ -50,6 +50,7 @@ opm/core/fluid/RockCompressibility.cpp \ opm/core/fluid/RockFromDeck.cpp \ opm/core/fluid/SaturationPropsBasic.cpp \ opm/core/fluid/SaturationPropsFromDeck.cpp \ +opm/core/fluid/SatFuncGwseg.cpp \ opm/core/fluid/SatFuncStone2.cpp \ opm/core/fluid/SatFuncSimple.cpp \ opm/core/fluid/blackoil/BlackoilPvtProperties.cpp \ @@ -148,6 +149,7 @@ opm/core/fluid/PvtPropertiesIncompFromDeck.hpp \ opm/core/fluid/RockBasic.hpp \ opm/core/fluid/RockCompressibility.hpp \ opm/core/fluid/RockFromDeck.hpp \ +opm/core/fluid/SatFuncGwseg.hpp \ opm/core/fluid/SatFuncStone2.hpp \ opm/core/fluid/SatFuncSimple.hpp \ opm/core/fluid/SaturationPropsBasic.hpp \ diff --git a/opm/core/fluid/SatFuncGwseg.cpp b/opm/core/fluid/SatFuncGwseg.cpp new file mode 100644 index 00000000..cd97e9af --- /dev/null +++ b/opm/core/fluid/SatFuncGwseg.cpp @@ -0,0 +1,440 @@ +/* + 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 + +namespace Opm +{ + + + + + void SatFuncGwsegUniform::init(const EclipseGridParser& deck, + const int table_num, + const PhaseUsage phase_usg, + const int samples) + { + phase_usage = phase_usg; + double swco = 0.0; + double swmax = 1.0; + if (phase_usage.phase_used[Aqua]) { + const SWOF::table_t& swof_table = deck.getSWOF().swof_; + const std::vector& sw = swof_table[table_num][0]; + const std::vector& krw = swof_table[table_num][1]; + const std::vector& krow = swof_table[table_num][2]; + const std::vector& pcow = swof_table[table_num][3]; + buildUniformMonotoneTable(sw, krw, samples, krw_); + buildUniformMonotoneTable(sw, krow, samples, krow_); + buildUniformMonotoneTable(sw, pcow, samples, pcow_); + krocw_ = krow[0]; // At connate water -> ecl. SWOF + swco = sw[0]; + smin_[phase_usage.phase_pos[Aqua]] = sw[0]; + swmax = sw.back(); + smax_[phase_usage.phase_pos[Aqua]] = sw.back(); + } + if (phase_usage.phase_used[Vapour]) { + const SGOF::table_t& sgof_table = deck.getSGOF().sgof_; + const std::vector& sg = sgof_table[table_num][0]; + const std::vector& krg = sgof_table[table_num][1]; + const std::vector& krog = sgof_table[table_num][2]; + const std::vector& pcog = sgof_table[table_num][3]; + buildUniformMonotoneTable(sg, krg, samples, krg_); + buildUniformMonotoneTable(sg, krog, samples, krog_); + buildUniformMonotoneTable(sg, pcog, samples, pcog_); + smin_[phase_usage.phase_pos[Vapour]] = sg[0]; + if (std::fabs(sg.back() + swco - 1.0) > 1e-3) { + THROW("Gas maximum saturation in SGOF table = " << sg.back() << + ", should equal (1.0 - connate water sat) = " << (1.0 - swco)); + } + smax_[phase_usage.phase_pos[Vapour]] = sg.back(); + } + // These only consider water min/max sats. Consider gas sats? + smin_[phase_usage.phase_pos[Liquid]] = 1.0 - swmax; + smax_[phase_usage.phase_pos[Liquid]] = 1.0 - swco; + } + + + void SatFuncGwsegUniform::evalKr(const double* s, double* kr) const + { + if (phase_usage.num_phases == 3) { + // Relative permeability model based on segregation of water + // and gas, with oil present in both water and gas zones. + const double swco = smin_[phase_usage.phase_pos[Aqua]]; + const double sw = std::max(s[Aqua], swco); + const double sg = s[Vapour]; + // xw and xg are the fractions occupied by water and gas zones. + const double eps = 1e-6; + const double xw = (sw - swco) / std::max(sg + sw - swco, eps); + const double xg = 1 - xw; + const double ssw = sg + sw; + const double ssg = sw - swco + sg; + const double krw = krw_(ssw); + const double krg = krg_(ssg); + const double krow = krow_(ssw); + const double krog = krog_(ssg); + kr[Aqua] = xw*krw; + kr[Vapour] = xg*krg; + kr[Liquid] = xw*krow + xg*krog; + return; + } + // We have a two-phase situation. We know that oil is active. + if (phase_usage.phase_used[Aqua]) { + int wpos = phase_usage.phase_pos[Aqua]; + int opos = phase_usage.phase_pos[Liquid]; + double sw = s[wpos]; + double krw = krw_(sw); + double krow = krow_(sw); + kr[wpos] = krw; + kr[opos] = krow; + } else { + ASSERT(phase_usage.phase_used[Vapour]); + int gpos = phase_usage.phase_pos[Vapour]; + int opos = phase_usage.phase_pos[Liquid]; + double sg = s[gpos]; + double krg = krg_(sg); + double krog = krog_(sg); + kr[gpos] = krg; + kr[opos] = krog; + } + } + + + void SatFuncGwsegUniform::evalKrDeriv(const double* s, double* kr, double* dkrds) const + { + const int np = phase_usage.num_phases; + std::fill(dkrds, dkrds + np*np, 0.0); + + if (np == 3) { + // Relative permeability model based on segregation of water + // and gas, with oil present in both water and gas zones. + const double swco = smin_[phase_usage.phase_pos[Aqua]]; + const double sw = std::max(s[Aqua], swco); + const double sg = s[Vapour]; + // xw and xg are the fractions occupied by water and gas zones. + const double eps = 1e-6; + const double xw = (sw - swco) / std::max(sg + sw - swco, eps); + const double xg = 1 - xw; + const double ssw = sg + sw; + const double ssg = sw - swco + sg; + const double krw = krw_(ssw); + const double krg = krg_(ssg); + const double krow = krow_(ssw); + const double krog = krog_(ssg); + kr[Aqua] = xw*krw; + kr[Vapour] = xg*krg; + kr[Liquid] = xw*krow + xg*krog; + + // Derivatives. + const double dkrww = krw_.derivative(ssw); + const double dkrgg = krg_.derivative(ssg); + const double dkrow = krow_.derivative(ssw); + const double dkrog = krog_.derivative(ssg); + const double d = ssg; // = sw - swco + sg (using 'd' for consistency with mrst docs). + dkrds[Aqua + Aqua*np] = (xg/d)*krw + xw*dkrww; + dkrds[Aqua + Vapour*np] = -(xw/d)*krw + xw*dkrww; + dkrds[Liquid + Aqua*np] = (xg/d)*krow + xw*dkrow - (xg/d)*krog + xg*dkrog; + dkrds[Liquid + Vapour*np] = -(xw/d)*krow + xw*dkrow + (xw/d)*krog + xg*dkrog; + dkrds[Vapour + Aqua*np] = -(xg/d)*krg + xg*dkrgg; + dkrds[Vapour + Vapour*np] = (xw/d)*krg + xg*dkrgg; + return; + } + // We have a two-phase situation. We know that oil is active. + if (phase_usage.phase_used[Aqua]) { + int wpos = phase_usage.phase_pos[Aqua]; + int opos = phase_usage.phase_pos[Liquid]; + double sw = s[wpos]; + double krw = krw_(sw); + double dkrww = krw_.derivative(sw); + double krow = krow_(sw); + double dkrow = krow_.derivative(sw); + kr[wpos] = krw; + kr[opos] = krow; + dkrds[wpos + wpos*np] = dkrww; + dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order. + } else { + ASSERT(phase_usage.phase_used[Vapour]); + int gpos = phase_usage.phase_pos[Vapour]; + int opos = phase_usage.phase_pos[Liquid]; + double sg = s[gpos]; + double krg = krg_(sg); + double dkrgg = krg_.derivative(sg); + double krog = krog_(sg); + double dkrog = krog_.derivative(sg); + kr[gpos] = krg; + kr[opos] = krog; + dkrds[gpos + gpos*np] = dkrgg; + dkrds[opos + gpos*np] = dkrog; + } + + } + + + void SatFuncGwsegUniform::evalPc(const double* s, double* pc) const + { + pc[phase_usage.phase_pos[Liquid]] = 0.0; + if (phase_usage.phase_used[Aqua]) { + int pos = phase_usage.phase_pos[Aqua]; + pc[pos] = pcow_(s[pos]); + } + if (phase_usage.phase_used[Vapour]) { + int pos = phase_usage.phase_pos[Vapour]; + pc[pos] = pcog_(s[pos]); + } + } + + void SatFuncGwsegUniform::evalPcDeriv(const double* s, double* pc, double* dpcds) const + { + // The problem of determining three-phase capillary pressures + // is very hard experimentally, usually one extends two-phase + // data (as for relative permeability). + // In our approach the derivative matrix is quite sparse, only + // the diagonal elements corresponding to non-oil phases are + // (potentially) nonzero. + const int np = phase_usage.num_phases; + std::fill(dpcds, dpcds + np*np, 0.0); + pc[phase_usage.phase_pos[Liquid]] = 0.0; + if (phase_usage.phase_used[Aqua]) { + int pos = phase_usage.phase_pos[Aqua]; + pc[pos] = pcow_(s[pos]); + dpcds[np*pos + pos] = pcow_.derivative(s[pos]); + } + if (phase_usage.phase_used[Vapour]) { + int pos = phase_usage.phase_pos[Vapour]; + pc[pos] = pcog_(s[pos]); + dpcds[np*pos + pos] = pcog_.derivative(s[pos]); + } + } + + + + + // ====== Methods for SatFuncGwsegNonuniform ====== + + + + + + void SatFuncGwsegNonuniform::init(const EclipseGridParser& deck, + const int table_num, + const PhaseUsage phase_usg, + const int /*samples*/) + { + phase_usage = phase_usg; + double swco = 0.0; + double swmax = 1.0; + if (phase_usage.phase_used[Aqua]) { + const SWOF::table_t& swof_table = deck.getSWOF().swof_; + const std::vector& sw = swof_table[table_num][0]; + const std::vector& krw = swof_table[table_num][1]; + const std::vector& krow = swof_table[table_num][2]; + const std::vector& pcow = swof_table[table_num][3]; + krw_ = NonuniformTableLinear(sw, krw); + krow_ = NonuniformTableLinear(sw, krow); + pcow_ = NonuniformTableLinear(sw, pcow); + krocw_ = krow[0]; // At connate water -> ecl. SWOF + swco = sw[0]; + smin_[phase_usage.phase_pos[Aqua]] = sw[0]; + swmax = sw.back(); + smax_[phase_usage.phase_pos[Aqua]] = sw.back(); + } + if (phase_usage.phase_used[Vapour]) { + const SGOF::table_t& sgof_table = deck.getSGOF().sgof_; + const std::vector& sg = sgof_table[table_num][0]; + const std::vector& krg = sgof_table[table_num][1]; + const std::vector& krog = sgof_table[table_num][2]; + const std::vector& pcog = sgof_table[table_num][3]; + krg_ = NonuniformTableLinear(sg, krg); + krog_ = NonuniformTableLinear(sg, krog); + pcog_ = NonuniformTableLinear(sg, pcog); + smin_[phase_usage.phase_pos[Vapour]] = sg[0]; + if (std::fabs(sg.back() + swco - 1.0) > 1e-3) { + THROW("Gas maximum saturation in SGOF table = " << sg.back() << + ", should equal (1.0 - connate water sat) = " << (1.0 - swco)); + } + smax_[phase_usage.phase_pos[Vapour]] = sg.back(); + } + // These only consider water min/max sats. Consider gas sats? + smin_[phase_usage.phase_pos[Liquid]] = 1.0 - swmax; + smax_[phase_usage.phase_pos[Liquid]] = 1.0 - swco; + } + + + void SatFuncGwsegNonuniform::evalKr(const double* s, double* kr) const + { + if (phase_usage.num_phases == 3) { + // Relative permeability model based on segregation of water + // and gas, with oil present in both water and gas zones. + const double swco = smin_[phase_usage.phase_pos[Aqua]]; + const double sw = std::max(s[Aqua], swco); + const double sg = s[Vapour]; + // xw and xg are the fractions occupied by water and gas zones. + const double eps = 1e-6; + const double xw = (sw - swco) / std::max(sg + sw - swco, eps); + const double xg = 1 - xw; + const double ssw = sg + sw; + const double ssg = sw - swco + sg; + const double krw = krw_(ssw); + const double krg = krg_(ssg); + const double krow = krow_(ssw); + const double krog = krog_(ssg); + kr[Aqua] = xw*krw; + kr[Vapour] = xg*krg; + kr[Liquid] = xw*krow + xg*krog; + return; + } + // We have a two-phase situation. We know that oil is active. + if (phase_usage.phase_used[Aqua]) { + int wpos = phase_usage.phase_pos[Aqua]; + int opos = phase_usage.phase_pos[Liquid]; + double sw = s[wpos]; + double krw = krw_(sw); + double krow = krow_(sw); + kr[wpos] = krw; + kr[opos] = krow; + } else { + ASSERT(phase_usage.phase_used[Vapour]); + int gpos = phase_usage.phase_pos[Vapour]; + int opos = phase_usage.phase_pos[Liquid]; + double sg = s[gpos]; + double krg = krg_(sg); + double krog = krog_(sg); + kr[gpos] = krg; + kr[opos] = krog; + } + } + + + void SatFuncGwsegNonuniform::evalKrDeriv(const double* s, double* kr, double* dkrds) const + { + const int np = phase_usage.num_phases; + std::fill(dkrds, dkrds + np*np, 0.0); + + if (np == 3) { + // Relative permeability model based on segregation of water + // and gas, with oil present in both water and gas zones. + const double swco = smin_[phase_usage.phase_pos[Aqua]]; + const double sw = std::max(s[Aqua], swco); + const double sg = s[Vapour]; + // xw and xg are the fractions occupied by water and gas zones. + const double eps = 1e-6; + const double xw = (sw - swco) / std::max(sg + sw - swco, eps); + const double xg = 1 - xw; + const double ssw = sg + sw; + const double ssg = sw - swco + sg; + const double krw = krw_(ssw); + const double krg = krg_(ssg); + const double krow = krow_(ssw); + const double krog = krog_(ssg); + kr[Aqua] = xw*krw; + kr[Vapour] = xg*krg; + kr[Liquid] = xw*krow + xg*krog; + + // Derivatives. + const double dkrww = krw_.derivative(ssw); + const double dkrgg = krg_.derivative(ssg); + const double dkrow = krow_.derivative(ssw); + const double dkrog = krog_.derivative(ssg); + const double d = ssg; // = sw - swco + sg (using 'd' for consistency with mrst docs). + dkrds[Aqua + Aqua*np] = (xg/d)*krw + xw*dkrww; + dkrds[Aqua + Vapour*np] = -(xw/d)*krw + xw*dkrww; + dkrds[Liquid + Aqua*np] = (xg/d)*krow + xw*dkrow - (xg/d)*krog + xg*dkrog; + dkrds[Liquid + Vapour*np] = -(xw/d)*krow + xw*dkrow + (xw/d)*krog + xg*dkrog; + dkrds[Vapour + Aqua*np] = -(xg/d)*krg + xg*dkrgg; + dkrds[Vapour + Vapour*np] = (xw/d)*krg + xg*dkrgg; + return; + } + // We have a two-phase situation. We know that oil is active. + if (phase_usage.phase_used[Aqua]) { + int wpos = phase_usage.phase_pos[Aqua]; + int opos = phase_usage.phase_pos[Liquid]; + double sw = s[wpos]; + double krw = krw_(sw); + double dkrww = krw_.derivative(sw); + double krow = krow_(sw); + double dkrow = krow_.derivative(sw); + kr[wpos] = krw; + kr[opos] = krow; + dkrds[wpos + wpos*np] = dkrww; + dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order. + } else { + ASSERT(phase_usage.phase_used[Vapour]); + int gpos = phase_usage.phase_pos[Vapour]; + int opos = phase_usage.phase_pos[Liquid]; + double sg = s[gpos]; + double krg = krg_(sg); + double dkrgg = krg_.derivative(sg); + double krog = krog_(sg); + double dkrog = krog_.derivative(sg); + kr[gpos] = krg; + kr[opos] = krog; + dkrds[gpos + gpos*np] = dkrgg; + dkrds[opos + gpos*np] = dkrog; + } + + } + + + void SatFuncGwsegNonuniform::evalPc(const double* s, double* pc) const + { + pc[phase_usage.phase_pos[Liquid]] = 0.0; + if (phase_usage.phase_used[Aqua]) { + int pos = phase_usage.phase_pos[Aqua]; + pc[pos] = pcow_(s[pos]); + } + if (phase_usage.phase_used[Vapour]) { + int pos = phase_usage.phase_pos[Vapour]; + pc[pos] = pcog_(s[pos]); + } + } + + void SatFuncGwsegNonuniform::evalPcDeriv(const double* s, double* pc, double* dpcds) const + { + // The problem of determining three-phase capillary pressures + // is very hard experimentally, usually one extends two-phase + // data (as for relative permeability). + // In our approach the derivative matrix is quite sparse, only + // the diagonal elements corresponding to non-oil phases are + // (potentially) nonzero. + const int np = phase_usage.num_phases; + std::fill(dpcds, dpcds + np*np, 0.0); + pc[phase_usage.phase_pos[Liquid]] = 0.0; + if (phase_usage.phase_used[Aqua]) { + int pos = phase_usage.phase_pos[Aqua]; + pc[pos] = pcow_(s[pos]); + dpcds[np*pos + pos] = pcow_.derivative(s[pos]); + } + if (phase_usage.phase_used[Vapour]) { + int pos = phase_usage.phase_pos[Vapour]; + pc[pos] = pcog_(s[pos]); + dpcds[np*pos + pos] = pcog_.derivative(s[pos]); + } + } + + + + + +} // namespace Opm diff --git a/opm/core/fluid/SatFuncGwseg.hpp b/opm/core/fluid/SatFuncGwseg.hpp new file mode 100644 index 00000000..1932e7f0 --- /dev/null +++ b/opm/core/fluid/SatFuncGwseg.hpp @@ -0,0 +1,80 @@ +/* + 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 SATFUNCGWSEG_HPP +#define SATFUNCGWSEG_HPP + +#include +#include +#include +#include +#include + +namespace Opm +{ + class SatFuncGwsegUniform : public BlackoilPhases + { + public: + void init(const EclipseGridParser& deck, + const int table_num, + const PhaseUsage phase_usg, + const int samples); + void evalKr(const double* s, double* kr) const; + void evalKrDeriv(const double* s, double* kr, double* dkrds) const; + void evalPc(const double* s, double* pc) const; + void evalPcDeriv(const double* s, double* pc, double* dpcds) const; + double smin_[PhaseUsage::MaxNumPhases]; + double smax_[PhaseUsage::MaxNumPhases]; + private: + PhaseUsage phase_usage; // A copy of the outer class' phase_usage_. + UniformTableLinear krw_; + UniformTableLinear krow_; + UniformTableLinear pcow_; + UniformTableLinear krg_; + UniformTableLinear krog_; + UniformTableLinear pcog_; + double krocw_; // = krow_(s_wc) + }; + + + class SatFuncGwsegNonuniform : public BlackoilPhases + { + public: + void init(const EclipseGridParser& deck, + const int table_num, + const PhaseUsage phase_usg, + const int samples); + void evalKr(const double* s, double* kr) const; + void evalKrDeriv(const double* s, double* kr, double* dkrds) const; + void evalPc(const double* s, double* pc) const; + void evalPcDeriv(const double* s, double* pc, double* dpcds) const; + double smin_[PhaseUsage::MaxNumPhases]; + double smax_[PhaseUsage::MaxNumPhases]; + private: + PhaseUsage phase_usage; // A copy of the outer class' phase_usage_. + NonuniformTableLinear krw_; + NonuniformTableLinear krow_; + NonuniformTableLinear pcow_; + NonuniformTableLinear krg_; + NonuniformTableLinear krog_; + NonuniformTableLinear pcog_; + double krocw_; // = krow_(s_wc) + }; + +} // namespace Opm +#endif // SATFUNCGWSEG_HPP