Added experimental Gauss-Seidel segregation solver. Not yet functioning.
This commit is contained in:
Atgeirr Flø Rasmussen
parent
fc3357e98a
commit
dd0cd85bf6
@ -439,10 +439,14 @@ main(int argc, char** argv)
|
||||
}
|
||||
bool use_segregation_split = false;
|
||||
bool use_column_solver = false;
|
||||
bool use_gauss_seidel_gravity = false;
|
||||
if (use_gravity && use_reorder) {
|
||||
use_segregation_split = param.getDefault("use_segregation_split", use_segregation_split);
|
||||
if (use_segregation_split) {
|
||||
use_column_solver = param.getDefault("use_column_solver", use_column_solver);
|
||||
if (use_column_solver) {
|
||||
use_gauss_seidel_gravity = param.getDefault("use_gauss_seidel_gravity", use_gauss_seidel_gravity);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@ -459,6 +463,9 @@ main(int argc, char** argv)
|
||||
const double nltol = param.getDefault("nl_tolerance", 1e-9);
|
||||
const int maxit = param.getDefault("nl_maxiter", 30);
|
||||
Opm::TransportModelTwophase reorder_model(*grid->c_grid(), &porevol[0], *props, nltol, maxit);
|
||||
if (use_gauss_seidel_gravity) {
|
||||
reorder_model.initGravity(grav);
|
||||
}
|
||||
// Non-reordering solver.
|
||||
TransportModel model (fluid, *grid->c_grid(), porevol, grav, guess_old_solution);
|
||||
if (use_gravity) {
|
||||
@ -701,7 +708,12 @@ main(int argc, char** argv)
|
||||
Opm::computeInjectedProduced(*props, state.saturation(), src, simtimer.currentStepLength(), injected, produced);
|
||||
if (use_segregation_split) {
|
||||
if (use_column_solver) {
|
||||
colsolver.solve(columns, simtimer.currentStepLength(), state.saturation());
|
||||
if (use_gauss_seidel_gravity) {
|
||||
reorder_model.solveGravity(columns, simtimer.currentStepLength(), &reorder_sat[0]);
|
||||
Opm::toBothSat(reorder_sat, state.saturation());
|
||||
} else {
|
||||
colsolver.solve(columns, simtimer.currentStepLength(), state.saturation());
|
||||
}
|
||||
} else {
|
||||
std::vector<double> fluxes = state.faceflux();
|
||||
std::fill(state.faceflux().begin(), state.faceflux().end(), 0.0);
|
||||
|
||||
@ -22,6 +22,7 @@
|
||||
#include <opm/core/grid.h>
|
||||
#include <opm/core/transport/reorder/reordersequence.h>
|
||||
#include <opm/core/utility/RootFinders.hpp>
|
||||
#include <opm/core/pressure/tpfa/trans_tpfa.h>
|
||||
|
||||
#include <fstream>
|
||||
#include <iterator>
|
||||
@ -49,7 +50,8 @@ namespace Opm
|
||||
source_(0),
|
||||
dt_(0.0),
|
||||
saturation_(0),
|
||||
fractionalflow_(grid.number_of_cells, -1.0)
|
||||
fractionalflow_(grid.number_of_cells, -1.0),
|
||||
mob_(2*grid.number_of_cells, -1.0)
|
||||
#ifdef EXPERIMENT_GAUSS_SEIDEL
|
||||
, ia_upw_(grid.number_of_cells + 1, -1),
|
||||
ja_upw_(grid.number_of_faces, -1),
|
||||
@ -430,6 +432,179 @@ namespace Opm
|
||||
|
||||
|
||||
|
||||
|
||||
// Residual function r(s) for a single-cell implicit Euler gravity segregation
|
||||
//
|
||||
// r(s) = s - s0 + dt/pv*sum_{j adj i}( gravmod_ij * gf_ij ).
|
||||
//
|
||||
struct TransportModelTwophase::GravityResidual
|
||||
{
|
||||
int cell;
|
||||
int nbcell[2];
|
||||
double s0;
|
||||
double dtpv; // dt/pv(i)
|
||||
double gf[2];
|
||||
const TransportModelTwophase& tm;
|
||||
explicit GravityResidual(const TransportModelTwophase& tmodel,
|
||||
const std::vector<int>& cells,
|
||||
const int pos,
|
||||
const double* gravflux) // Always oriented towards next in column. Size = colsize - 1.
|
||||
: tm(tmodel)
|
||||
{
|
||||
cell = cells[pos];
|
||||
nbcell[0] = -1;
|
||||
gf[0] = 0.0;
|
||||
if (pos > 0) {
|
||||
nbcell[0] = cells[pos - 1];
|
||||
gf[0] = -gravflux[pos - 1];
|
||||
}
|
||||
nbcell[1] = -1;
|
||||
gf[1] = 0.0;
|
||||
if (pos < int(cells.size() - 1)) {
|
||||
nbcell[1] = cells[pos + 1];
|
||||
gf[1] = gravflux[pos];
|
||||
}
|
||||
s0 = tm.saturation_[cell];
|
||||
dtpv = tm.dt_/tm.porevolume_[cell];
|
||||
|
||||
}
|
||||
double operator()(double s) const
|
||||
{
|
||||
double res = s - s0;
|
||||
|
||||
for (int nb = 0; nb < 2; ++nb) {
|
||||
if (nbcell[nb] != -1) {
|
||||
if (gf[nb] < 0.0) {
|
||||
double m[2] = { tm.mob_[2*cell], tm.mob_[2*nbcell[nb] + 1] };
|
||||
res += -dtpv*gf[nb]*m[0]*m[1]/(m[0] + m[1]);
|
||||
} else {
|
||||
double m[2] = { tm.mob_[2*nbcell[nb]], tm.mob_[2*cell + 1] };
|
||||
res += -dtpv*gf[nb]*m[0]*m[1]/(m[0] + m[1]);
|
||||
}
|
||||
}
|
||||
}
|
||||
return res;
|
||||
}
|
||||
};
|
||||
|
||||
void TransportModelTwophase::mobility(double s, int cell, double* mob) const
|
||||
{
|
||||
double sat[2] = { s, 1.0 - s };
|
||||
props_.relperm(1, sat, &cell, mob, 0);
|
||||
mob[0] /= visc_[0];
|
||||
mob[1] /= visc_[1];
|
||||
}
|
||||
|
||||
|
||||
|
||||
void TransportModelTwophase::initGravity(const double* grav)
|
||||
{
|
||||
// Set up gravflux_ = T_ij g (rho_w - rho_o) (z_i - z_j)
|
||||
std::vector<double> htrans(grid_.cell_facepos[grid_.number_of_cells]);
|
||||
const int nf = grid_.number_of_faces;
|
||||
const int dim = grid_.dimensions;
|
||||
gravflux_.resize(nf);
|
||||
tpfa_htrans_compute(const_cast<UnstructuredGrid*>(&grid_), props_.permeability(), &htrans[0]);
|
||||
tpfa_trans_compute(const_cast<UnstructuredGrid*>(&grid_), &htrans[0], &gravflux_[0]);
|
||||
const double delta_rho = props_.density()[0] - props_.density()[1];
|
||||
for (int f = 0; f < nf; ++f) {
|
||||
const int* c = &grid_.face_cells[2*f];
|
||||
double gdz = 0.0;
|
||||
for (int d = 0; d < dim; ++d) {
|
||||
gdz += grav[d]*(grid_.cell_centroids[dim*c[0] + d] - grid_.cell_centroids[dim*c[1] + d]);
|
||||
}
|
||||
gravflux_[f] *= delta_rho*gdz;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
void TransportModelTwophase::solveSingleCellGravity(const std::vector<int>& cells,
|
||||
const int pos,
|
||||
const double* gravflux)
|
||||
{
|
||||
const int cell = cells[pos];
|
||||
GravityResidual res(*this, cells, pos, gravflux);
|
||||
int iters_used;
|
||||
saturation_[cell] = modifiedRegulaFalsi(res, smin_[2*cell], smax_[2*cell], maxit_, tol_, iters_used);
|
||||
mobility(saturation_[cell], cell, &mob_[2*cell]);
|
||||
}
|
||||
|
||||
|
||||
|
||||
void TransportModelTwophase::solveGravityColumn(const std::vector<int>& cells)
|
||||
{
|
||||
// Set up column gravflux.
|
||||
const int nc = cells.size();
|
||||
std::vector<double> col_gravflux(nc - 1);
|
||||
for (int ci = 0; ci < nc - 1; ++ci) {
|
||||
const int cell = cells[ci];
|
||||
const int next_cell = cells[ci + 1];
|
||||
for (int j = grid_.cell_facepos[cell]; j < grid_.cell_facepos[cell+1]; ++j) {
|
||||
const int face = grid_.cell_faces[j];
|
||||
const int c1 = grid_.face_cells[2*face + 0];
|
||||
const int c2 = grid_.face_cells[2*face + 1];
|
||||
if (c1 == next_cell || c2 == next_cell) {
|
||||
const double gf = gravflux_[face];
|
||||
col_gravflux[ci] = (c1 == cell) ? gf : -gf;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Solve single cell problems, repeating if necessary.
|
||||
double max_s_change = 0.0;
|
||||
int num_iters = 0;
|
||||
do {
|
||||
max_s_change = 0.0;
|
||||
for (int ci = 0; ci < nc; ++ci) {
|
||||
double old_s[2] = { saturation_[cells[ci]],
|
||||
saturation_[cells[nc-ci]] };
|
||||
solveSingleCellGravity(cells, ci, &col_gravflux[0]);
|
||||
solveSingleCellGravity(cells, nc - ci, &col_gravflux[0]);
|
||||
max_s_change = std::max(max_s_change, std::max(std::fabs(saturation_[cells[ci]] - old_s[0]),
|
||||
std::fabs(saturation_[cells[nc-ci]] - old_s[1])));
|
||||
}
|
||||
std::cout << "Iter = " << num_iters << " max_s_change = " << max_s_change << std::endl;
|
||||
} while (max_s_change > tol_ && ++num_iters < maxit_);
|
||||
|
||||
if (max_s_change > tol_) {
|
||||
THROW("In solveGravityColumn(), we did not converge after "
|
||||
<< num_iters << " iterations. Delta s = " << max_s_change);
|
||||
}
|
||||
// Repeat if necessary.
|
||||
}
|
||||
|
||||
|
||||
|
||||
void TransportModelTwophase::solveGravity(const std::map<int, std::vector<int> >& columns,
|
||||
const double dt,
|
||||
double* saturation)
|
||||
{
|
||||
// Initialize mobilities.
|
||||
const int nc = grid_.number_of_cells;
|
||||
std::vector<int> cells(nc);
|
||||
for (int c = 0; c < nc; ++c) {
|
||||
cells[c] = c;
|
||||
}
|
||||
mob_.resize(2*nc);
|
||||
props_.relperm(cells.size(), saturation, &cells[0], &mob_[0], 0);
|
||||
const double* mu = props_.viscosity();
|
||||
for (int c = 0; c < nc; ++c) {
|
||||
mob_[2*c] /= mu[0];
|
||||
mob_[2*c + 1] /= mu[1];
|
||||
}
|
||||
|
||||
// Set up other variables.
|
||||
dt_ = dt;
|
||||
saturation_ = saturation;
|
||||
|
||||
// Solve on all columns.
|
||||
std::map<int, std::vector<int> >::const_iterator it;
|
||||
for (it = columns.begin(); it != columns.end(); ++it) {
|
||||
solveGravityColumn(it->second);
|
||||
}
|
||||
}
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
|
||||
|
||||
@ -22,6 +22,7 @@
|
||||
|
||||
#include <opm/core/transport/reorder/TransportModelInterface.hpp>
|
||||
#include <vector>
|
||||
#include <map>
|
||||
|
||||
struct UnstructuredGrid;
|
||||
|
||||
@ -47,6 +48,15 @@ namespace Opm
|
||||
virtual void solveSingleCell(const int cell);
|
||||
virtual void solveMultiCell(const int num_cells, const int* cells);
|
||||
|
||||
void initGravity(const double* grav);
|
||||
void solveSingleCellGravity(const std::vector<int>& cells,
|
||||
const int pos,
|
||||
const double* gravflux);
|
||||
void solveGravityColumn(const std::vector<int>& cells);
|
||||
void solveGravity(const std::map<int, std::vector<int> >& columns,
|
||||
const double dt,
|
||||
double* saturation);
|
||||
|
||||
private:
|
||||
const UnstructuredGrid& grid_;
|
||||
const double* porevolume_; // one volume per cell
|
||||
@ -61,7 +71,10 @@ namespace Opm
|
||||
const double* source_; // one source per cell
|
||||
double dt_;
|
||||
double* saturation_; // one per cell
|
||||
std::vector<double> fractionalflow_; // one per cell
|
||||
std::vector<double> fractionalflow_; // = m[0]/(m[0] + m[1]) per cell
|
||||
// For gravity segregation.
|
||||
std::vector<double> gravflux_;
|
||||
std::vector<double> mob_;
|
||||
|
||||
// Storing the upwind and downwind graphs for experiments.
|
||||
std::vector<int> ia_upw_;
|
||||
@ -71,6 +84,9 @@ namespace Opm
|
||||
|
||||
struct Residual;
|
||||
double fracFlow(double s, int cell) const;
|
||||
|
||||
struct GravityResidual;
|
||||
void mobility(double s, int cell, double* mob) const;
|
||||
};
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
Loading…
Reference in New Issue
Block a user