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

Implemented solve() with rock compressibility (untested). Struct init warning suppression.

Browse Source
This commit is contained in:
Atgeirr Flø Rasmussen 2012-03-19 11:52:03 +01:00
parent 9f540e9cf7
commit bc04c0fa3e
2 changed files with 112 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -23,6 +23,7 @@
#include <opm/core/pressure/mimetic/mimetic.h> #include <opm/core/pressure/mimetic/mimetic.h>
#include <opm/core/pressure/flow_bc.h> #include <opm/core/pressure/flow_bc.h>
#include <opm/core/linalg/LinearSolverInterface.hpp> #include <opm/core/linalg/LinearSolverInterface.hpp>
#include <opm/core/linalg/sparse_sys.h>
#include <opm/core/utility/ErrorMacros.hpp> #include <opm/core/utility/ErrorMacros.hpp>
namespace Opm 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]; } if (! src.empty()) { F.src = &src[0]; }
F.bc = bcs; F.bc = bcs;
@ -120,7 +121,88 @@ namespace Opm
pressure.resize(grid_.number_of_cells); pressure.resize(grid_.number_of_cells);
faceflux.resize(grid_.number_of_faces); 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
int j = h_->A->ia[c];
for (; j < h_->A->ia[c+1]; ++j) {
if (h_->A->ja[j] == c) {
break;
}
}
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.cell_press = &pressure[0];
soln.face_flux = &faceflux[0]; soln.face_flux = &faceflux[0];
@ -129,5 +211,4 @@ namespace Opm
} // namespace Opm } // namespace Opm

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

@ -55,7 +55,7 @@ namespace Opm
/// Destructor. /// Destructor.
~IncompTpfa(); ~IncompTpfa();
/// Assemble and solve pressure system. /// Assemble and solve incompressible pressure system.
/// \param[in] totmob Must contain N total mobility values (one per cell). /// \param[in] totmob Must contain N total mobility values (one per cell).
/// totmob = \sum_{p} kr_p/mu_p. /// totmob = \sum_{p} kr_p/mu_p.
/// \param[in] omega Must be empty if constructor gravity argument was null. /// \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>& pressure,
std::vector<double>& faceflux); 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. /// Expose read-only reference to internal half-transmissibility.
const ::std::vector<double>& getHalfTrans() const { return htrans_; } const ::std::vector<double>& getHalfTrans() const { return htrans_; }