From 125cb33c368a95dfb07af69088817187f7dbb976 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 17 Jan 2012 10:25:49 +0100 Subject: [PATCH] Added SaturationPropsBasic class. --- opm/core/fluid/SaturationPropsBasic.cpp | 196 ++++++++++++++++++++++++ opm/core/fluid/SaturationPropsBasic.hpp | 86 +++++++++++ 2 files changed, 282 insertions(+) create mode 100644 opm/core/fluid/SaturationPropsBasic.cpp create mode 100644 opm/core/fluid/SaturationPropsBasic.hpp diff --git a/opm/core/fluid/SaturationPropsBasic.cpp b/opm/core/fluid/SaturationPropsBasic.cpp new file mode 100644 index 000000000..077576b04 --- /dev/null +++ b/opm/core/fluid/SaturationPropsBasic.cpp @@ -0,0 +1,196 @@ +/* + 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 + +namespace Opm +{ + + + // ---------- Helper functions ---------- + + namespace { + + struct KrFunConstant + { + double kr(double) + { + return 1.0; + } + double dkrds(double) + { + return 0.0; + } + }; + + struct KrFunLinear + { + double kr(double s) + { + return s; + } + double dkrds(double) + { + return 1.0; + } + }; + + struct KrFunQuadratic + { + double kr(double s) + { + return s*s; + } + double dkrds(double s) + { + return 2.0*s; + } + }; + + + template + static inline void evalAllKrDeriv(const int n, const int np, + const double* s, double* kr, double* dkrds, Fun fun) + { + if (dkrds == 0) { +#pragma omp parallel for + for (int i = 0; i < n*np; ++i) { + kr[i] = fun.kr(s[i]); + } + return; + } +#pragma omp parallel for + for (int i = 0; i < n; ++i) { + std::fill(dkrds + i*np*np, dkrds + (i+1)*np*np, 0.0); + for (int phase = 0; phase < np; ++phase) { + kr[i*np + phase] = fun.kr(s[i*np + phase]); + // Only diagonal elements in derivative. + dkrds[i*np*np + phase*np + phase] = fun.dkrds(s[i*np + phase]); + } + } + } + + + } // anon namespace + + + + // ---------- Class methods ---------- + + + + /// Default constructor. + SaturationPropsBasic::SaturationPropsBasic() + { + } + + + + + /// Initialize from parameters. + void SaturationPropsBasic::init(const Dune::parameter::ParameterGroup& param) + { + int num_phases = param.getDefault("num_phases", 2); + if (num_phases > 2 || num_phases < 1) { + THROW("SaturationPropsBasic::init() illegal num_phases: " << num_phases); + } + num_phases_ = num_phases; + std::string rpf = param.getDefault("relperm_func", std::string("Linear")); + if (rpf == "Constant") { + relperm_func_ = Constant; + } else if (rpf == "Linear") { + relperm_func_ = Linear; + } else if (rpf == "Quadratic") { + relperm_func_ = Quadratic; + } else { + THROW("SaturationPropsBasic::init() illegal relperm_func: " << rpf); + } + } + + + + + /// Relative permeability. + /// \param[in] n Number of data points. + /// \param[in] s Array of nP saturation values. + /// \param[out] kr Array of nP relperm values, array must be valid before calling. + /// \param[out] dkrds If non-null: array of nP^2 relperm derivative values, + /// array must be valid before calling. + /// The P^2 derivative matrix is + /// m_{ij} = \frac{dkr_i}{ds^j}, + /// and is output in Fortran order (m_00 m_10 m_20 m01 ...) + void SaturationPropsBasic::relperm(const int n, + const double* s, + double* kr, + double* dkrds) const + { + switch (relperm_func_) { + case Constant: + { + evalAllKrDeriv(n, num_phases_, s, kr, dkrds, KrFunConstant()); + break; + } + case Linear: + { + evalAllKrDeriv(n, num_phases_, s, kr, dkrds, KrFunLinear()); + break; + } + case Quadratic: + { + evalAllKrDeriv(n, num_phases_, s, kr, dkrds, KrFunQuadratic()); + break; + } + default: + THROW("SaturationPropsBasic::relperm() unhandled relperm func type: " << relperm_func_); + } + } + + + + + /// Capillary pressure. + /// \param[in] n Number of data points. + /// \param[in] s Array of nP saturation values. + /// \param[out] pc Array of nP capillary pressure values, array must be valid before calling. + /// \param[out] dpcds If non-null: array of nP^2 derivative values, + /// array must be valid before calling. + /// The P^2 derivative matrix is + /// m_{ij} = \frac{dpc_i}{ds^j}, + /// and is output in Fortran order (m_00 m_10 m_20 m01 ...) + void SaturationPropsBasic::capPress(const int n, + const double* /*s*/, + double* pc, + double* dpcds) const + { + std::fill(pc, pc + num_phases_*n, 0.0); + if (dpcds) { + std::fill(dpcds, dpcds + num_phases_*num_phases_*n, 0.0); + } + } + + + + + + +} // namespace Opm + + diff --git a/opm/core/fluid/SaturationPropsBasic.hpp b/opm/core/fluid/SaturationPropsBasic.hpp new file mode 100644 index 000000000..3a8e97fa4 --- /dev/null +++ b/opm/core/fluid/SaturationPropsBasic.hpp @@ -0,0 +1,86 @@ +/* + 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_SATURATIONPROPSBASIC_HEADER_INCLUDED +#define OPM_SATURATIONPROPSBASIC_HEADER_INCLUDED + +#include + +namespace Opm +{ + + + /// Class encapsulating basic saturation function behaviour, + /// by which we mean constant, linear or quadratic relative + /// permeability functions for a maximum of two phases, + /// and zero capillary pressure. + class SaturationPropsBasic + { + public: + /// Default constructor. + SaturationPropsBasic(); + + /// Initialize from parameters. + /// The following parameters are accepted (defaults): + /// num_phases (2) Must be 1 or 2. + /// relperm_func ("Linear") Must be "Constant", "Linear" or "Quadratic". + void init(const Dune::parameter::ParameterGroup& param); + + /// Relative permeability. + /// \param[in] n Number of data points. + /// \param[in] s Array of nP saturation values. + /// \param[out] kr Array of nP relperm values, array must be valid before calling. + /// \param[out] dkrds If non-null: array of nP^2 relperm derivative values, + /// array must be valid before calling. + /// The P^2 derivative matrix is + /// m_{ij} = \frac{dkr_i}{ds^j}, + /// and is output in Fortran order (m_00 m_10 m_20 m01 ...) + void relperm(const int n, + const double* s, + double* kr, + double* dkrds) const; + + /// Capillary pressure. + /// \param[in] n Number of data points. + /// \param[in] s Array of nP saturation values. + /// \param[out] pc Array of nP capillary pressure values, array must be valid before calling. + /// \param[out] dpcds If non-null: array of nP^2 derivative values, + /// array must be valid before calling. + /// The P^2 derivative matrix is + /// m_{ij} = \frac{dpc_i}{ds^j}, + /// and is output in Fortran order (m_00 m_10 m_20 m01 ...) + void capPress(const int n, + const double* s, + double* pc, + double* dpcds) const; + + private: + enum RelPermFunc { Constant, Linear, Quadratic }; + int num_phases_; + RelPermFunc relperm_func_; + }; + + + +} // namespace Opm + + + + +#endif // OPM_SATURATIONPROPSBASIC_HEADER_INCLUDED