OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

removed bugs. can compile. Not tested.

Browse Source
This commit is contained in:
Xavier Raynaud 2012-03-19 16:33:32 +01:00
parent e7e9c03a67
commit 542ed94b84
3 changed files with 49 additions and 42 deletions

7
opm/polymer/GravityColumnSolverPolymer.hpp
View File

@ -54,9 +54,8 @@ namespace Opm
const double dt,
std::vector<double>& s,
std::vector<double>& c,
std::vector<double>& sol_s_vec
std::vector<double>& sol_c_vec
);
std::vector<double>& sol_vec
);
Model& model_;
const UnstructuredGrid& grid_;
const double tol_;
@ -65,6 +64,6 @@ namespace Opm
} // namespace Opm
#include <opm/core/transport/GravityColumnSolverPolymer_impl.hpp>
#include <opm/polymer/GravityColumnSolverPolymer_impl.hpp>
#endif // OPM_GRAVITYCOLUMNSOLVERPOLYMER_HEADER_INCLUDED

20
opm/polymer/GravityColumnSolverPolymer_impl.hpp
View File

@ -93,15 +93,15 @@ namespace Opm
// compute the position where to store the coefficient of a matrix A_{i,j} (i,j=0,...,N-1)
// in a array which is sent to the band matrix solver of LAPACK.
int operator ()(int i, int j) {
int operator ()(int i, int j) const {
return kl_ + ku_ + i - j + j*nrow_;
}
const int N_;
const int ku_;
const int kl_;
const int nrow_;
const int N_;
}
};
} // anon namespace
@ -158,19 +158,21 @@ namespace Opm
/// \param[in] column_cells the cells on which to solve the segregation
/// problem. Must be in a single vertical column,
/// and ordered (direction doesn't matter).
template <class Model>
void GravityColumnSolver<Model>::solveSingleColumn(const std::vector<int>& column_cells,
const double dt,
std::vector<double>& s,
std::vector<double>& sol_vec)
void GravityColumnSolverPolymer<Model>::solveSingleColumn(const std::vector<int>& column_cells,
const double dt,
std::vector<double>& s,
std::vector<double>& c,
std::vector<double>& sol_vec)
{
// This is written only to work with SinglePointUpwindTwoPhase,
// not with arbitrary problem models.
const int col_size = column_cells.size();
StateWithZeroFlux state(s); // This holds s by reference.
StateWithZeroFlux state(s, c); // This holds s by reference.
// Assemble.
const int kl = 3;
@ -235,8 +237,8 @@ namespace Opm
// Solve.
const int num_rhs = 1;
int info = 0;
int ipiv;
// Solution will be written to rhs.
dgtsv_(&col_size, &num_rhs, DL, D, DU, &rhs[0], &col_size, &info);
dgbsv_(&col_size, &kl, &ku, &num_rhs, &hm[0], &nrow, &ipiv, &rhs[0], &col_size, &info);
if (info != 0) {
THROW("Lapack reported error in dgtsv: " << info);

64
opm/polymer/SinglePointUpwindTwoPhasePolymer.hpp
View File

@ -44,7 +44,7 @@
#include <iostream>
namespace Opm {
namespace spu_2p {
namespace polymer_reorder {
class ModelParameterStorage {
public:
ModelParameterStorage(int nc, int totconn)
@ -73,9 +73,9 @@ namespace Opm {
mob_ = &data_[0];
dmobds_ = mob_ + (2 * nc );
dmobwatdc_ = dmobwatds_ + (2 * nc );
mc_ = dmobdc_ + (1 * nc );
dmcdc_ = mc_ + (1 * nc );
dmobwatdc_ = dmobds_ + (2 * nc );
mc_ = dmobwatdc_ + (1 * nc );
dmcdc_ = mc_ + (1 * nc );
porevol_ = dmcdc_ + (1 * nc );
deadporespace_ = porevol_ + (1 * nc );
dg_ = deadporespace_ + (1 * nc );
@ -97,8 +97,8 @@ namespace Opm {
double* dmobds (int cell) { return dmobds_ + (2*cell + 0); }
const double* dmobds (int cell) const { return dmobds_ + (2*cell + 0); }
double& dmobwatdc (int cell) { return dmobdc_[cell]; }
double dmobwatdc (int cell) const { return dmobdc_[cell]; }
double& dmobwatdc (int cell) { return dmobwatdc_[cell]; }
double dmobwatdc (int cell) const { return dmobwatdc_[cell]; }
double& mc (int cell) { return mc_[cell]; }
double mc (int cell) const { return mc_[cell]; }
@ -140,13 +140,16 @@ namespace Opm {
double drho_ ;
double *mob_ ;
double *dmobds_ ;
double *dmobdc_ ;
double *dmobwatdc_ ;
double *mc_ ;
double *dmcdc_ ;
double *porevol_ ;
double *deadporespace_ ;
double *dg_ ;
double *sw_ ;
double *c_ ;
double *ds_ ;
double *dsc_ ;
double *dsc_ ;
double *pc_ ;
double *dpc_ ;
double *trans_ ;
@ -250,7 +253,7 @@ namespace Opm {
const int f ,
double* F , // F[0] = s-residual, F[1] = c-residual
double* dFd1 , //Jacobi matrix for residual with respect to variables in cell
double* dFd2 , //Jacobi matrix for residual with respect to variables in OTHER cell
double* dFd2 //Jacobi matrix for residual with respect to variables in OTHER cell
//dFd1[0]= d(F[0])/d(s1), dFd1[1]= d(F[0])/d(c1), dFd1[2]= d(F[1])/d(s1), dFd1[3]= d(F[1])/d(c1),
//dFd2[0]= d(F[0])/d(s2), dFd2[1]= d(F[0])/d(c2), dFd2[2]= d(F[1])/d(s2), dFd2[3]= d(F[1])/d(c2).
@ -296,22 +299,22 @@ namespace Opm {
dFcds[0] = &dFd1[2]; dFcds[1] = &dFd2[2];
dFcdc[0] = &dFd1[3]; dFcdc[1] = &dFd2[3];
// sign is positive
*dFd1[0] += sgn*dt * f1 * dpcflux[0] * m[1];
*dFd2[0] += sgn*dt * f1 * dpcflux[1] * m[1];
*dFd1[2] += sgn*dt * f1 * mc * dpcflux[0] * m[1];
*dFd2[2] += sgn*dt * f1 * mc * dpcflux[1] * m[1];
dFd1[0] += sgn*dt * f1 * dpcflux[0] * m[1];
dFd2[0] += sgn*dt * f1 * dpcflux[1] * m[1];
dFd1[2] += sgn*dt * f1 * mc * dpcflux[0] * m[1];
dFd2[2] += sgn*dt * f1 * mc * dpcflux[1] * m[1];
// We assume that the capillary pressure is independent of the polymer concentration.
// Hence, no more contributions.
} else {
dFsds[0] = dFsds2; dFsds[1] = dFsds1;
dFsdc[0] = dFsdc2; dFsdc[1] = dFsdc1;
dFcds[0] = dFcds2; dFcds[1] = dFcds1;
dFcdc[0] = dFcdc2; dFcdc[1] = dFcdc1;
dFsds[0] = &dFd2[0]; dFsds[1] = &dFd1[0];
dFsdc[0] = &dFd2[1]; dFsdc[1] = &dFd1[1];
dFcds[0] = &dFd2[2]; dFcds[1] = &dFd1[2];
dFcdc[0] = &dFd2[3]; dFcdc[1] = &dFd1[3];
// sign is negative
*dFd1[0] += sgn*dt * f1 * dpcflux[1] * m[1];
*dFd2[0] += sgn*dt * f1 * dpcflux[0] * m[1];
*dFd1[2] += sgn*dt * f1 * mc * dpcflux[1] * m[1];
*dFd2[2] += sgn*dt * f1 * mc * dpcflux[0] * m[1];
dFd1[0] += sgn*dt * f1 * dpcflux[1] * m[1];
dFd2[0] += sgn*dt * f1 * dpcflux[0] * m[1];
dFd1[2] += sgn*dt * f1 * mc * dpcflux[1] * m[1];
dFd2[2] += sgn*dt * f1 * mc * dpcflux[0] * m[1];
// We assume that the capillary pressure is independent of the polymer concentration.
// Hence, no more contributions.
}
@ -374,9 +377,9 @@ namespace Opm {
*F += dt * dflux * src->saturation[2*i + 0];
} else {
// cell -> src
const int c = src->cell[i];
const double* m = store_.mob (c);
const double* dm = store_.dmob(c);
const int cell = src->cell[i];
const double* m = store_.mob (cell);
const double* dm = store_.dmobds(cell);
const double mt = m[0] + m[1];
@ -460,7 +463,10 @@ namespace Opm {
const Grid& g ,
JacobianSystem& sys) {
double s[2], mob[2], dmobds[2 * 2], dmobdc[2 * 2];
std::vector<double> s(2, 0.);
std::vector<double> mob(2, 0.);
std::vector<double> dmobds(4, 0.);
double dmobwatdc;
double c;
double mc, dmcdc;
double pc, dpc;
@ -497,8 +503,8 @@ namespace Opm {
s[0] = std::min(s_max, s[0]);
s[1] = 1 - s[0];
fluid_.mobility(cell, s, mob, dmobds, dmobwatdc);
fluid_.mc(cell, s, mc, dmcdc);
fluid_.mobility(cell, s, c, mob, dmobds, dmobwatdc);
fluid_.mc(c, mc, dmcdc);
fluid_.pc(cell, s, pc, dpc);
store_.mob (cell)[0] = mob [0];
@ -654,10 +660,10 @@ namespace Opm {
dpcflux[1] = store_.trans(f)*store_.dpc(i2);
}
TwophaseFluid fluid_ ;
TwophaseFluidPolymer fluid_ ;
const double* gravity_;
std::vector<int> f2hf_ ;
spu_2p::ModelParameterStorage store_ ;
polymer_reorder::ModelParameterStorage store_ ;
bool init_step_use_previous_sol_;
double sat_tol_;
};