Using miscUtilities.hpp from opm-core instead of keeping own versions.

This commit is contained in:
Atgeirr Flø Rasmussen
2012-02-26 21:18:47 +01:00
parent 753ebedc0d
commit e9b17dee02
4 changed files with 17 additions and 268 deletions

View File

@@ -12,6 +12,4 @@ noinst_PROGRAMS = \
polymer_reorder
polymer_reorder_SOURCES = \
polymer_reorder.cpp \
Utilities.hpp \
Utilities.cpp
polymer_reorder.cpp

View File

@@ -1,164 +0,0 @@
/*
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 <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include "Utilities.hpp"
#include <opm/core/pressure/tpfa/ifs_tpfa.h>
#include <opm/core/pressure/tpfa/trans_tpfa.h>
#include <opm/core/utility/cart_grid.h>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/core/utility/cpgpreprocess/cgridinterface.h>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/fluid/SimpleFluid2p.hpp>
#include <opm/core/fluid/IncompPropertiesBasic.hpp>
#include <opm/core/fluid/IncompPropertiesFromDeck.hpp>
#include <opm/core/transport/CSRMatrixUmfpackSolver.hpp>
#include <boost/filesystem/convenience.hpp>
#include <boost/scoped_ptr.hpp>
#include <boost/lexical_cast.hpp>
#include <cassert>
#include <cstddef>
#include <algorithm>
#include <tr1/array>
#include <functional>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <iterator>
#include <vector>
#include <numeric>
namespace Opm
{
void
compute_porevolume(const UnstructuredGrid* g,
const Opm::IncompPropertiesInterface& props,
std::vector<double>& porevol)
{
int num_cells = g->number_of_cells;
porevol.resize(num_cells);
const double* poro = props.porosity();
::std::transform(poro, poro + num_cells,
g->cell_volumes,
porevol.begin(),
::std::multiplies<double>());
}
void
compute_totmob(const Opm::IncompPropertiesInterface& props,
const std::vector<double>& s,
std::vector<double>& totmob)
{
int num_cells = props.numCells();
int num_phases = props.numPhases();
totmob.resize(num_cells);
ASSERT(int(s.size()) == num_cells*num_phases);
std::vector<int> cells(num_cells);
for (int cell = 0; cell < num_cells; ++cell) {
cells[cell] = cell;
}
std::vector<double> kr(num_cells*num_phases);
props.relperm(num_cells, &s[0], &cells[0], &kr[0], 0);
const double* mu = props.viscosity();
for (int cell = 0; cell < num_cells; ++cell) {
totmob[cell] = 0;
for (int phase = 0; phase < num_phases; ++phase) {
totmob[cell] += kr[2*cell + phase]/mu[phase];
}
}
}
void
compute_totmob_omega(const Opm::IncompPropertiesInterface& props,
const std::vector<double>& s,
std::vector<double>& totmob,
std::vector<double>& omega)
{
int num_cells = props.numCells();
int num_phases = props.numPhases();
totmob.resize(num_cells);
omega.resize(num_cells);
ASSERT(int(s.size()) == num_cells*num_phases);
std::vector<int> cells(num_cells);
for (int cell = 0; cell < num_cells; ++cell) {
cells[cell] = cell;
}
std::vector<double> kr(num_cells*num_phases);
props.relperm(num_cells, &s[0], &cells[0], &kr[0], 0);
const double* mu = props.viscosity();
for (int cell = 0; cell < num_cells; ++cell) {
totmob[cell] = 0.0;
for (int phase = 0; phase < num_phases; ++phase) {
totmob[cell] += kr[2*cell + phase]/mu[phase];
}
}
const double* rho = props.density();
for (int cell = 0; cell < num_cells; ++cell) {
omega[cell] = 0.0;
for (int phase = 0; phase < num_phases; ++phase) {
omega[cell] += rho[phase]*(kr[2*cell + phase]/mu[phase])/totmob[cell];
}
}
}
void toWaterSat(const std::vector<double>& sboth, std::vector<double>& sw)
{
int num = sboth.size()/2;
sw.resize(num);
for (int i = 0; i < num; ++i) {
sw[i] = sboth[2*i];
}
}
void toBothSat(const std::vector<double>& sw, std::vector<double>& sboth)
{
int num = sw.size();
sboth.resize(2*num);
for (int i = 0; i < num; ++i) {
sboth[2*i] = sw[i];
sboth[2*i + 1] = 1.0 - sw[i];
}
}
} // namespace Opm

View File

@@ -1,92 +0,0 @@
/*
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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_UTILITIES_HEADER_INCLUDED
#define OPM_UTILITIES_HEADER_INCLUDED
#include <opm/core/pressure/tpfa/ifs_tpfa.h>
#include <opm/core/pressure/tpfa/trans_tpfa.h>
#include <opm/core/utility/cart_grid.h>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/core/utility/cpgpreprocess/cgridinterface.h>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/fluid/SimpleFluid2p.hpp>
#include <opm/core/fluid/IncompPropertiesBasic.hpp>
#include <opm/core/fluid/IncompPropertiesFromDeck.hpp>
#include <opm/core/transport/CSRMatrixUmfpackSolver.hpp>
#include <boost/filesystem/convenience.hpp>
#include <boost/scoped_ptr.hpp>
#include <boost/lexical_cast.hpp>
#include <cassert>
#include <cstddef>
#include <algorithm>
#include <tr1/array>
#include <functional>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <iterator>
#include <vector>
#include <numeric>
namespace Opm
{
void
compute_porevolume(const UnstructuredGrid* g,
const Opm::IncompPropertiesInterface& props,
std::vector<double>& porevol);
void
compute_totmob(const Opm::IncompPropertiesInterface& props,
const std::vector<double>& s,
std::vector<double>& totmob);
void
compute_totmob_omega(const Opm::IncompPropertiesInterface& props,
const std::vector<double>& s,
std::vector<double>& totmob,
std::vector<double>& omega);
void toWaterSat(const std::vector<double>& sboth, std::vector<double>& sw);
void toBothSat(const std::vector<double>& sw, std::vector<double>& sboth);
} // namespace Opm
#endif // OPM_UTILITIES_HEADER_INCLUDED

View File

@@ -20,14 +20,13 @@
#include "config.h"
#include "Utilities.hpp"
#include <opm/core/pressure/IncompTpfa.hpp>
#include <opm/core/grid.h>
#include <opm/core/GridManager.hpp>
#include <opm/core/utility/writeVtkData.hpp>
#include <opm/core/utility/linearInterpolation.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/StopWatch.hpp>
#include <opm/core/utility/Units.hpp>
@@ -319,7 +318,7 @@ main(int argc, char** argv)
// Extra rock init.
std::vector<double> porevol;
compute_porevolume(grid->c_grid(), *props, porevol);
computePorevolume(*grid->c_grid(), *props, porevol);
double tot_porevol = std::accumulate(porevol.begin(), porevol.end(), 0.0);
// Gravity init.
@@ -348,26 +347,34 @@ main(int argc, char** argv)
method, nltol, maxit);
// State-related and source-related variables init.
int num_cells = grid->c_grid()->number_of_cells;
std::vector<double> totmob;
std::vector<double> omega; // Empty dummy unless/until we include gravity here.
double init_sat = param.getDefault("init_sat", 0.0);
ReservoirState state(grid->c_grid(), props->numPhases(), init_sat);
// We need a separate reorder_sat, because the reorder
// code expects a scalar sw, not both sw and so.
std::vector<double> reorder_sat(grid->c_grid()->number_of_cells);
std::vector<double> reorder_sat(num_cells);
double flow_per_sec = 0.1*tot_porevol/Opm::unit::day;
if (param.has("injection_rate_per_day")) {
flow_per_sec = param.get<double>("injection_rate_per_day")/Opm::unit::day;
}
std::vector<double> src (grid->c_grid()->number_of_cells, 0.0);
src[0] = flow_per_sec;
src[grid->c_grid()->number_of_cells - 1] = -flow_per_sec;
std::vector<double> src(num_cells, 0.0);
src[0] = flow_per_sec;
src[num_cells - 1] = -flow_per_sec;
std::vector<double> reorder_src = src;
// Control init.
double current_time = 0.0;
double total_time = stepsize*num_psteps;
// The allcells vector is used in calls to computeTotalMobility()
// and computeTotalMobilityOmega().
std::vector<int> allcells(num_cells);
for (int cell = 0; cell < num_cells; ++cell) {
allcells[cell] = cell;
}
// Warn if any parameters are unused.
if (param.anyUnused()) {
std::cout << "-------------------- Unused parameters: --------------------\n";
@@ -401,9 +408,9 @@ main(int argc, char** argv)
}
if (use_gravity) {
compute_totmob_omega(*props, state.saturation(), totmob, omega);
computeTotalMobilityOmega(*props, allcells, state.saturation(), totmob, omega);
} else {
compute_totmob(*props, state.saturation(), totmob);
computeTotalMobility(*props, allcells, state.saturation(), totmob);
}
pressure_timer.start();
psolver.solve(totmob, omega, src, state.pressure(), state.faceflux());