OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

merge.

Browse Source
This commit is contained in:
Xavier Raynaud 2012-03-19 15:05:36 +01:00
commit 94b0c76050
5 changed files with 251 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
Makefile.am
View File

@ -29,6 +29,7 @@ opm/core/fluid/IncompPropertiesFromDeck.cpp \
opm/core/fluid/PvtPropertiesBasic.cpp \
opm/core/fluid/PvtPropertiesIncompFromDeck.cpp \
opm/core/fluid/RockBasic.cpp \
opm/core/fluid/RockCompressibility.cpp \
opm/core/fluid/RockFromDeck.cpp \
opm/core/fluid/SaturationPropsBasic.cpp \
opm/core/fluid/SaturationPropsFromDeck.cpp \
@ -107,6 +108,7 @@ opm/core/fluid/IncompPropertiesFromDeck.hpp \
opm/core/fluid/PvtPropertiesBasic.hpp \
opm/core/fluid/PvtPropertiesIncompFromDeck.hpp \
opm/core/fluid/RockBasic.hpp \
opm/core/fluid/RockCompressibility.hpp \
opm/core/fluid/RockFromDeck.hpp \
opm/core/fluid/SaturationPropsBasic.hpp \
opm/core/fluid/SaturationPropsFromDeck.hpp \

83
opm/core/fluid/RockCompressibility.cpp Normal file
View File

@ -0,0 +1,83 @@
/*
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 <opm/core/fluid/RockCompressibility.hpp>
#include <opm/core/eclipse/EclipseGridParser.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/linearInterpolation.hpp>
namespace Opm
{
RockCompressibility::RockCompressibility(const parameter::ParameterGroup& param)
: pref_(0.0),
rock_comp_(0.0)
{
pref_ = param.getDefault("rock_compressibility_pref", 100.0)*unit::barsa;
rock_comp_ = param.getDefault("rock_compressibility", 0.0)/unit::barsa;
}
RockCompressibility::RockCompressibility(const EclipseGridParser& deck)
: pref_(0.0),
rock_comp_(0.0)
{
if (deck.hasField("ROCKTAB")) {
const table_t& rt = deck.getROCKTAB().rocktab_;
int n = rt[0].size();
p_.resize(n);
poromult_.resize(n);
for (int i = 0; i < n; ++i) {
p_[i] = rt[0][i][0];
poromult_[i] = rt[0][i][1];
}
} else if (deck.hasField("ROCK")) {
const ROCK& r = deck.getROCK();
pref_ = r.rock_compressibilities_[0][0];
rock_comp_ = r.rock_compressibilities_[0][1];
} else {
std::cout << "**** warning: no rock compressibility data found in deck (ROCK or ROCKTAB)." << std::endl;
}
}
double RockCompressibility::poroMult(double pressure)
{
if (p_.empty()) {
// Approximating with a quadratic curve.
const double cpnorm = rock_comp_*(pressure - pref_);
return (1.0 + cpnorm + 0.5*cpnorm*cpnorm);
} else {
return Opm::linearInterpolation(p_, poromult_, pressure);
}
}
double RockCompressibility::rockComp(double pressure)
{
if (p_.empty()) {
return rock_comp_;
} else {
const double poromult = Opm::linearInterpolation(p_, poromult_, pressure);
const double dporomultdp = Opm::linearInterpolationDerivative(p_, poromult_, pressure);
return dporomultdp/poromult;
}
}
} // namespace Opm

59
opm/core/fluid/RockCompressibility.hpp Normal file
View File

@ -0,0 +1,59 @@
/*
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_ROCKCOMPRESSIBILITY_HEADER_INCLUDED
#define OPM_ROCKCOMPRESSIBILITY_HEADER_INCLUDED
#include <vector>
namespace Opm
{
class EclipseGridParser;
namespace parameter { class ParameterGroup; }
class RockCompressibility
{
public:
/// Construct from input deck.
RockCompressibility(const EclipseGridParser& deck);
/// Construct from parameters.
/// Accepts the following parameters (with defaults).
/// rock_compressibility_pref (100.0) [given in bar]
/// rock_compressibility (0.0) [given in bar^{-1}]
RockCompressibility(const parameter::ParameterGroup& param);
/// Porosity multiplier.
double poroMult(double pressure);
/// Rock compressibility = (d poro / d p)*(1 / poro).
double rockComp(double pressure);
private:
std::vector<double> p_;
std::vector<double> poromult_;
double pref_;
double rock_comp_;
};
} // namespace Opm
#endif // OPM_ROCKCOMPRESSIBILITY_HEADER_INCLUDED

82
opm/core/pressure/IncompTpfa.cpp
View File

@ -23,6 +23,7 @@
#include <opm/core/pressure/mimetic/mimetic.h>
#include <opm/core/pressure/flow_bc.h>
#include <opm/core/linalg/LinearSolverInterface.hpp>
#include <opm/core/linalg/sparse_sys.h>
#include <opm/core/utility/ErrorMacros.hpp>
namespace Opm
@ -109,7 +110,7 @@ namespace Opm
}
}
ifs_tpfa_forces F = { 0 };
ifs_tpfa_forces F = { NULL, NULL };
if (! src.empty()) { F.src = &src[0]; }
F.bc = bcs;
@ -120,7 +121,83 @@ namespace Opm
pressure.resize(grid_.number_of_cells);
faceflux.resize(grid_.number_of_faces);
ifs_tpfa_solution soln = { 0 };
ifs_tpfa_solution soln = { NULL, NULL };
soln.cell_press = &pressure[0];
soln.face_flux = &faceflux[0];
ifs_tpfa_press_flux(gg, &F, &trans_[0], h_, &soln);
}
/// Assemble and solve pressure system with rock compressibility (assumed constant per cell).
/// \param[in] totmob Must contain N total mobility values (one per cell).
/// totmob = \sum_{p} kr_p/mu_p.
/// \param[in] omega Must be empty if constructor gravity argument was null.
/// Otherwise must contain N mobility-weighted density values (one per cell).
/// omega = \frac{\sum_{p} mob_p rho_p}{\sum_p rho_p}.
/// \param[in] src Must contain N source rates (one per cell).
/// Positive values represent total inflow rates,
/// negative values represent total outflow rates.
/// \param[in] bcs If non-null, specifies boundary conditions.
/// If null, noflow conditions are assumed.
/// \param[in] porevol Must contain N pore volumes.
/// \param[in] rock_comp Must contain N rock compressibilities.
/// rock_comp = (d poro / d p)*(1/poro).
/// \param[in] dt Timestep.
/// \param[out] pressure Will contain N cell-pressure values.
/// \param[out] faceflux Will contain F signed face flux values.
void IncompTpfa::solve(const std::vector<double>& totmob,
const std::vector<double>& omega,
const std::vector<double>& src,
const FlowBoundaryConditions* bcs,
const std::vector<double>& porevol,
const std::vector<double>& rock_comp,
const double dt,
std::vector<double>& pressure,
std::vector<double>& faceflux)
{
UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
tpfa_eff_trans_compute(gg, &totmob[0], &htrans_[0], &trans_[0]);
if (!omega.empty()) {
if (gpress_.empty()) {
THROW("Nozero omega argument given, but gravity was null in constructor.");
}
mim_ip_density_update(gg->number_of_cells, gg->cell_facepos,
&omega[0],
&gpress_[0], &gpress_omegaweighted_[0]);
} else {
if (!gpress_.empty()) {
THROW("Empty omega argument given, but gravity was non-null in constructor.");
}
}
ifs_tpfa_forces F = { NULL, NULL };
if (! src.empty()) { F.src = &src[0]; }
F.bc = bcs;
ifs_tpfa_assemble(gg, &F, &trans_[0], &gpress_omegaweighted_[0], h_);
if (!rock_comp.empty()) {
// The extra term of the equation is
//
// porevol*rock_comp*(p - p0)/dt.
//
// The p part goes on the diagonal, the p0 on the rhs.
for (int c = 0; c < gg->number_of_cells; ++c) {
// Find diagonal
size_t j = csrmatrix_elm_index(c, c, h_->A);
h_->A->sa[j] += porevol[c]*rock_comp[c]/dt;
h_->b[c] += porevol[c]*rock_comp[c]*pressure[c]/dt;
}
}
linsolver_.solve(h_->A, h_->b, h_->x);
pressure.resize(grid_.number_of_cells);
faceflux.resize(grid_.number_of_faces);
ifs_tpfa_solution soln = { NULL, NULL };
soln.cell_press = &pressure[0];
soln.face_flux = &faceflux[0];
@ -129,5 +206,4 @@ namespace Opm
} // namespace Opm

29
opm/core/pressure/IncompTpfa.hpp
View File

@ -55,7 +55,7 @@ namespace Opm
/// Destructor.
~IncompTpfa();
/// Assemble and solve pressure system.
/// Assemble and solve incompressible pressure system.
/// \param[in] totmob Must contain N total mobility values (one per cell).
/// totmob = \sum_{p} kr_p/mu_p.
/// \param[in] omega Must be empty if constructor gravity argument was null.
@ -75,6 +75,33 @@ namespace Opm
std::vector<double>& pressure,
std::vector<double>& faceflux);
/// Assemble and solve pressure system with rock compressibility (assumed constant per cell).
/// \param[in] totmob Must contain N total mobility values (one per cell).
/// totmob = \sum_{p} kr_p/mu_p.
/// \param[in] omega Must be empty if constructor gravity argument was null.
/// Otherwise must contain N mobility-weighted density values (one per cell).
/// omega = \frac{\sum_{p} mob_p rho_p}{\sum_p rho_p}.
/// \param[in] src Must contain N source rates (one per cell).
/// Positive values represent total inflow rates,
/// negative values represent total outflow rates.
/// \param[in] bcs If non-null, specifies boundary conditions.
/// If null, noflow conditions are assumed.
/// \param[in] porevol Must contain N pore volumes.
/// \param[in] rock_comp Must contain N rock compressibilities.
/// rock_comp = (d poro / d p)*(1/poro).
/// \param[in] dt Timestep.
/// \param[out] pressure Will contain N cell-pressure values.
/// \param[out] faceflux Will contain F signed face flux values.
void solve(const std::vector<double>& totmob,
const std::vector<double>& omega,
const std::vector<double>& src,
const FlowBoundaryConditions* bcs,
const std::vector<double>& porevol,
const std::vector<double>& rock_comp,
const double dt,
std::vector<double>& pressure,
std::vector<double>& faceflux);
/// Expose read-only reference to internal half-transmissibility.
const ::std::vector<double>& getHalfTrans() const { return htrans_; }