From eb502faae79b4db670c6ff9bb1c094d143fce18b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 12 Apr 2012 14:10:47 +0200 Subject: [PATCH] Moved AdHocProps to new Inc.Props.DefaultPolymer class. Also moved some functions in polymer_reorder.cpp to be more similar to spu_2p.cpp. --- Makefile.am | 15 +- examples/polymer_reorder.cpp | 172 +++++------------- .../IncompPropertiesDefaultPolymer.hpp | 119 ++++++++++++ 3 files changed, 170 insertions(+), 136 deletions(-) create mode 100644 opm/polymer/IncompPropertiesDefaultPolymer.hpp diff --git a/Makefile.am b/Makefile.am index 5aff607a3..241ae66b6 100644 --- a/Makefile.am +++ b/Makefile.am @@ -10,11 +10,12 @@ libopmpolymer_la_SOURCES = \ opm/polymer/TransportModelPolymer.cpp \ opm/polymer/polymerUtilities.cpp -nobase_include_HEADERS = \ -opm/polymer/GravityColumnSolverPolymer.hpp \ -opm/polymer/GravityColumnSolverPolymer_impl.hpp \ -opm/polymer/PolymerProperties.hpp \ -opm/polymer/PolymerState.hpp \ -opm/polymer/SinglePointUpwindTwoPhasePolymer.hpp \ -opm/polymer/TransportModelPolymer.hpp \ +nobase_include_HEADERS = \ +opm/polymer/GravityColumnSolverPolymer.hpp \ +opm/polymer/GravityColumnSolverPolymer_impl.hpp \ +opm/polymer/IncompPropertiesDefaultPolymer.hpp \ +opm/polymer/PolymerProperties.hpp \ +opm/polymer/PolymerState.hpp \ +opm/polymer/SinglePointUpwindTwoPhasePolymer.hpp \ +opm/polymer/TransportModelPolymer.hpp \ opm/polymer/polymerUtilities.hpp diff --git a/examples/polymer_reorder.cpp b/examples/polymer_reorder.cpp index f45ee8aa6..e6ad646af 100644 --- a/examples/polymer_reorder.cpp +++ b/examples/polymer_reorder.cpp @@ -38,7 +38,7 @@ #include #include -#include +#include #include #include @@ -84,90 +84,54 @@ -class AdHocProps : public Opm::IncompPropertiesBasic + +static void outputState(const UnstructuredGrid& grid, + const Opm::PolymerState& state, + const int step, + const std::string& output_dir) { -public: - AdHocProps(const Opm::parameter::ParameterGroup& param, int dim, int num_cells) - : Opm::IncompPropertiesBasic(param, dim, num_cells) - { - ASSERT(numPhases() == 2); - sw_.resize(3); - sw_[0] = 0.2; - sw_[1] = 0.7; - sw_[2] = 1.0; - krw_.resize(3); - krw_[0] = 0.0; - krw_[1] = 0.7; - krw_[2] = 1.0; - so_.resize(2); - so_[0] = 0.3; - so_[1] = 0.8; - kro_.resize(2); - kro_[0] = 0.0; - kro_[1] = 1.0; + // Write data in VTK format. + std::ostringstream vtkfilename; + vtkfilename << output_dir << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu"; + std::ofstream vtkfile(vtkfilename.str().c_str()); + if (!vtkfile) { + THROW("Failed to open " << vtkfilename.str()); } + Opm::DataMap dm; + dm["saturation"] = &state.saturation(); + dm["pressure"] = &state.pressure(); + dm["concentration"] = &state.concentration(); + dm["cmax"] = &state.maxconcentration(); + std::vector cell_velocity; + Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity); + dm["velocity"] = &cell_velocity; + Opm::writeVtkData(grid, dm, vtkfile); - virtual void relperm(const int n, - const double* s, - const int* /*cells*/, - double* kr, - double* dkrds) const - { - // ASSERT(dkrds == 0); - // We assume two phases flow - for (int i = 0; i < n; ++i) { - kr[2*i] = krw(s[2*i]); - kr[2*i+1] = kro(s[2*i+1]); - if (dkrds != 0) { - dkrds[2*i] = krw_dsw(s[2*i]); - dkrds[2*i+3] = kro_dso(s[2*i+1]); - dkrds[2*i+1] = -dkrds[2*i+3]; - dkrds[2*i+2] = -dkrds[2*i]; - } + // Write data (not grid) in Matlab format + for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) { + std::ostringstream fname; + fname << output_dir << "/" << it->first << "-" << std::setw(3) << std::setfill('0') << step << ".dat"; + std::ofstream file(fname.str().c_str()); + if (!file) { + THROW("Failed to open " << fname.str()); } + const std::vector& d = *(it->second); + std::copy(d.begin(), d.end(), std::ostream_iterator(file, "\n")); } +} - virtual void satRange(const int n, - const int* /*cells*/, - double* smin, - double* smax) const - { - const int np = 2; - for (int i = 0; i < n; ++i) { - smin[np*i + 0] = sw_[0]; - smax[np*i + 0] = sw_.back(); - smin[np*i + 1] = 1.0 - sw_[0]; - smax[np*i + 1] = 1.0 - sw_.back(); - } + +static void outputWaterCut(const Opm::Watercut& watercut, + const std::string& output_dir) +{ + // Write water cut curve. + std::string fname = output_dir + "/watercut.txt"; + std::ofstream os(fname.c_str()); + if (!os) { + THROW("Failed to open " << fname); } - -private: - double krw(double s) const - { - return Opm::linearInterpolation(sw_, krw_, s); - } - - double krw_dsw(double s) const - { - return Opm::linearInterpolationDerivative(sw_, krw_, s); - } - - - double kro(double s) const - { - return Opm::linearInterpolation(so_, kro_, s); - } - - double kro_dso(double s) const - { - return Opm::linearInterpolationDerivative(so_, kro_, s); - } - - std::vector sw_; - std::vector krw_; - std::vector so_; - std::vector kro_; -}; + watercut.write(os); +} // --------------- Types needed to define transport solver --------------- @@ -301,42 +265,6 @@ typedef Opm::ImplicitTransport cell_velocity; - Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity); - dm["velocity"] = &cell_velocity; - Opm::writeVtkData(grid, dm, vtkfile); - - // Write data (not grid) in Matlab format - for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) { - std::ostringstream fname; - fname << output_dir << "/" << it->first << "-" << std::setw(3) << std::setfill('0') << step << ".dat"; - std::ofstream file(fname.str().c_str()); - if (!file) { - THROW("Failed to open " << fname.str()); - } - const std::vector& d = *(it->second); - std::copy(d.begin(), d.end(), std::ostream_iterator(file, "\n")); - } -} - - class PolymerInflow @@ -364,19 +292,6 @@ private: -static void outputWaterCut(const Opm::Watercut& watercut, - const std::string& output_dir) -{ - // Write water cut curve. - std::string fname = output_dir + "/watercut.txt"; - std::ofstream os(fname.c_str()); - if (!os) { - THROW("Failed to open " << fname); - } - watercut.write(os); -} - - // ----------------- Main program ----------------- int main(int argc, char** argv) @@ -418,7 +333,6 @@ main(int argc, char** argv) const int* gc = grid->c_grid()->global_cell; std::vector global_cell(gc, gc + grid->c_grid()->number_of_cells); props.reset(new Opm::IncompPropertiesFromDeck(deck, global_cell)); - // props.reset(new AdHocProps(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells)); // Wells init. wells.reset(new Opm::WellsManager(deck, *grid->c_grid(), props->permeability())); // Timer init. @@ -444,7 +358,7 @@ main(int argc, char** argv) grid.reset(new Opm::GridManager(nx, ny, nz, dx, dy, dz)); // Rock and fluid init. // props.reset(new Opm::IncompPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells)); - props.reset(new AdHocProps(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells)); + props.reset(new Opm::IncompPropertiesDefaultPolymer(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells)); // Wells init. wells.reset(new Opm::WellsManager()); // Timer init. diff --git a/opm/polymer/IncompPropertiesDefaultPolymer.hpp b/opm/polymer/IncompPropertiesDefaultPolymer.hpp new file mode 100644 index 000000000..411eb20c6 --- /dev/null +++ b/opm/polymer/IncompPropertiesDefaultPolymer.hpp @@ -0,0 +1,119 @@ +/* + 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_INCOMPPROPERTIESDEFAULTPOLYMER_HEADER_INCLUDED +#define OPM_INCOMPPROPERTIESDEFAULTPOLYMER_HEADER_INCLUDED + +#include +#include +#include +#include +#include + +namespace Opm +{ + + class IncompPropertiesDefaultPolymer : public Opm::IncompPropertiesBasic + { + public: + IncompPropertiesDefaultPolymer(const Opm::parameter::ParameterGroup& param, int dim, int num_cells) + : Opm::IncompPropertiesBasic(param, dim, num_cells) + { + ASSERT(numPhases() == 2); + sw_.resize(3); + sw_[0] = 0.2; + sw_[1] = 0.7; + sw_[2] = 1.0; + krw_.resize(3); + krw_[0] = 0.0; + krw_[1] = 0.7; + krw_[2] = 1.0; + so_.resize(2); + so_[0] = 0.3; + so_[1] = 0.8; + kro_.resize(2); + kro_[0] = 0.0; + kro_[1] = 1.0; + } + + virtual void relperm(const int n, + const double* s, + const int* /*cells*/, + double* kr, + double* dkrds) const + { + // ASSERT(dkrds == 0); + // We assume two phases flow + for (int i = 0; i < n; ++i) { + kr[2*i] = krw(s[2*i]); + kr[2*i+1] = kro(s[2*i+1]); + if (dkrds != 0) { + dkrds[2*i] = krw_dsw(s[2*i]); + dkrds[2*i+3] = kro_dso(s[2*i+1]); + dkrds[2*i+1] = -dkrds[2*i+3]; + dkrds[2*i+2] = -dkrds[2*i]; + } + } + } + + virtual void satRange(const int n, + const int* /*cells*/, + double* smin, + double* smax) const + { + const int np = 2; + for (int i = 0; i < n; ++i) { + smin[np*i + 0] = sw_[0]; + smax[np*i + 0] = sw_.back(); + smin[np*i + 1] = 1.0 - sw_[0]; + smax[np*i + 1] = 1.0 - sw_.back(); + } + } + + private: + double krw(double s) const + { + return Opm::linearInterpolation(sw_, krw_, s); + } + + double krw_dsw(double s) const + { + return Opm::linearInterpolationDerivative(sw_, krw_, s); + } + + + double kro(double s) const + { + return Opm::linearInterpolation(so_, kro_, s); + } + + double kro_dso(double s) const + { + return Opm::linearInterpolationDerivative(so_, kro_, s); + } + + std::vector sw_; + std::vector krw_; + std::vector so_; + std::vector kro_; + }; + +} // namespace Opm + +#endif // OPM_INCOMPPROPERTIESDEFAULTPOLYMER_HEADER_INCLUDED