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

Added bound checks to Newton column solver.

Browse Source
This commit is contained in:
Xavier Raynaud
2012-05-23 16:22:48 +02:00
parent 263a41b84b
commit 4b805c191e
3 changed files with 119 additions and 37 deletions
Download Patch File Download Diff File Expand all files Collapse all files

32
examples/polymer_reorder.cpp
View File

@@ -237,6 +237,10 @@ public:
return smax_[2*c + 0];
}
double cMax() const
{
return polyprops_.cMax();
}
template <class PolyC,
class Mc,
@@ -255,7 +259,7 @@ private:
};
typedef PolymerFluid2pWrappingProps TwophaseFluidPolymer;
typedef Opm::SinglePointUpwindTwoPhasePolymer<TwophaseFluidPolymer> NewtonPolymerTransportModel;
typedef Opm::SinglePointUpwindTwoPhasePolymer<TwophaseFluidPolymer> FluxModel;
using namespace Opm::ImplicitTransportDefault;
@@ -271,7 +275,7 @@ public:
}
};
typedef Opm::ImplicitTransport<NewtonPolymerTransportModel,
typedef Opm::ImplicitTransport<FluxModel,
JacSys ,
MaxNorm ,
VectorNegater ,
@@ -403,12 +407,16 @@ main(int argc, char** argv)
if (param.has("poly_init")) {
double poly_init = param.getDefault("poly_init", 0.0);
for (int cell = 0; cell < grid->c_grid()->number_of_cells; ++cell) {
if (state.saturation()[2*cell] > 0) {
double smin[2], smax[2];
props->satRange(1, &cell, smin, smax);
if (state.saturation()[2*cell] > 0.5*(smin[0] + smax[0])) {
state.concentration()[cell] = poly_init;
state.maxconcentration()[cell] = poly_init;
} else {
state.saturation()[2*cell + 0] = 0.;
state.saturation()[2*cell + 1] = 1.;
state.concentration()[cell] = 0.;
state.maxconcentration()[cell] = poly_init;
state.maxconcentration()[cell] = 0.;
}
}
}
@@ -435,10 +443,10 @@ main(int argc, char** argv)
std::vector<double> ads_vals(3, -1e100);
ads_vals[0] = 0.0;
// polyprop.ads_vals[1] = param.getDefault("c_max_ads", 0.0025);
// ads_vals[1] = 0.0015;
// ads_vals[2] = 0.0025;
ads_vals[1] = 0.0;
ads_vals[2] = 0.0;
ads_vals[1] = 0.0015;
ads_vals[2] = 0.0025;
// ads_vals[1] = 0.0;
// ads_vals[2] = 0.0;
polyprop.set(c_max, mix_param, rock_density, dead_pore_vol, res_factor, c_max_ads,
static_cast<Opm::PolymerProperties::AdsorptionBehaviour>(ads_index),
c_vals_visc, visc_mult_vals, c_vals_ads, ads_vals);
@@ -552,17 +560,17 @@ main(int argc, char** argv)
reorder_model.initGravity(grav);
}
// Non-reordering solver.
NewtonPolymerTransportModel model(fluid, *grid->c_grid(), porevol, grav, guess_old_solution);
FluxModel fmodel(fluid, *grid->c_grid(), porevol, grav, guess_old_solution);
if (use_gravity) {
model.initGravityTrans(*grid->c_grid(), psolver.getHalfTrans());
fmodel.initGravityTrans(*grid->c_grid(), psolver.getHalfTrans());
}
TransportSolver tsolver(model);
TransportSolver tsolver(fmodel);
// Column-based gravity segregation solver.
std::vector<std::vector<int> > columns;
if (use_column_solver) {
Opm::extractColumn(*grid->c_grid(), columns);
}
Opm::GravityColumnSolverPolymer<NewtonPolymerTransportModel> colsolver(model, *grid->c_grid(), nl_tolerance, nl_maxiter);
Opm::GravityColumnSolverPolymer<FluxModel, TwophaseFluidPolymer> colsolver(fmodel, fluid, *grid->c_grid(), nl_tolerance, nl_maxiter);
// // // Not implemented for polymer.
// // Control init.

14
opm/polymer/GravityColumnSolverPolymer.hpp
View File

@@ -29,16 +29,17 @@ namespace Opm
/// Class for doing gravity segregation (only),
/// on a vertical column of cells.
template <class Model>
template <class FluxModel, class Model>
class GravityColumnSolverPolymer
{
public:
/// Note: the model will be changed since it stores computed
/// quantities in itself, such as mobilities.
GravityColumnSolverPolymer(Model& model,
const UnstructuredGrid& grid,
const double tol,
const int maxit);
GravityColumnSolverPolymer(FluxModel& fmodel,
const Model& model,
const UnstructuredGrid& grid,
const double tol,
const int maxit);
/// \param[in] columns for each column (with logical cartesian indices as key),
/// contains the cells on which to solve the segregation
@@ -59,7 +60,8 @@ namespace Opm
std::vector<double>& cmax,
std::vector<double>& sol_vec
);
Model& model_;
FluxModel& fmodel_;
const Model& model_;
const UnstructuredGrid& grid_;
const double tol_;
const int maxit_;

110
opm/polymer/GravityColumnSolverPolymer_impl.hpp
View File

@@ -41,15 +41,24 @@
#include <opm/core/utility/ErrorMacros.hpp>
#include <iterator>
// namespace {
// int sgn(double val) {
// return (0.0 < val) - (val < 0);
// }
// }
namespace Opm
{
template <class Model>
GravityColumnSolverPolymer<Model>::GravityColumnSolverPolymer(Model& model,
const UnstructuredGrid& grid,
const double tol,
const int maxit)
: model_(model), grid_(grid), tol_(tol), maxit_(maxit)
template <class FluxModel, class Model>
GravityColumnSolverPolymer<FluxModel, Model>::GravityColumnSolverPolymer(FluxModel& fmodel,
const Model& model,
const UnstructuredGrid& grid,
const double tol,
const int maxit)
: fmodel_(fmodel), model_(model), grid_(grid), tol_(tol), maxit_(maxit)
{
}
@@ -117,8 +126,8 @@ namespace Opm
/// problem. For each column, its cells must be in a single
/// vertical column, connected and ordered
/// (direction doesn't matter).
template <class Model>
void GravityColumnSolverPolymer<Model>::solve(const std::vector<std::vector<int> >& columns,
template <class FluxModel, class Model>
void GravityColumnSolverPolymer<FluxModel, Model>::solve(const std::vector<std::vector<int> >& columns,
const double dt,
std::vector<double>& s,
std::vector<double>& c,
@@ -129,20 +138,83 @@ namespace Opm
StateWithZeroFlux state(s, c, cmax); // This holds s, c and cmax by reference.
JacSys sys(2*grid_.number_of_cells);
std::vector<double> increment(2*grid_.number_of_cells, 0.0);
model_.initStep(state, grid_, sys);
fmodel_.initStep(state, grid_, sys);
int iter = 0;
double max_delta = 1e100;
const double cmax_cell = 2.0*model_.cMax();
const double tol_c_cell = 1e-2*cmax_cell;
while (iter < maxit_) {
model_.initIteration(state, grid_, sys);
fmodel_.initIteration(state, grid_, sys);
int size = columns.size();
for(int i = 0; i < size; ++i) {
solveSingleColumn(columns[i], dt, s, c, cmax, increment);
}
for (int cell = 0; cell < grid_.number_of_cells; ++cell) {
sys.vector().writableSolution()[2*cell + 0] += increment[2*cell + 0];
sys.vector().writableSolution()[2*cell + 1] += increment[2*cell + 1];
}
double& s_cell = sys.vector().writableSolution()[2*cell + 0];
double& c_cell = sys.vector().writableSolution()[2*cell + 1];
s_cell += increment[2*cell + 0];
c_cell += increment[2*cell + 1];
if (s_cell < 0.) {
double& incr = increment[2*cell + 0];
s_cell -= incr;
if (std::fabs(incr) < 1e-2) {
incr = -s_cell;
s_cell = 0.;
} else {
incr = -s_cell/2.0;
s_cell = s_cell/2.0;
}
}
if (s_cell > 1.) {
double& incr = increment[2*cell + 0];
s_cell -= incr;
if (std::fabs(incr) < 1e-2) {
incr = 1. - s_cell;
s_cell = 1.;
} else {
incr = (1 - s_cell)/2.0;
s_cell = (1 + s_cell)/2.0;
}
}
if (c_cell < 0.) {
double& incr = increment[2*cell + 1];
c_cell -= incr;
if (std::fabs(incr) < tol_c_cell) {
incr = -c_cell;
c_cell = 0.;
} else {
incr = -c_cell/2.0;
c_cell = c_cell/2.0;
}
}
if (c_cell > cmax_cell) {
double& incr = increment[2*cell + 1];
c_cell -= incr;
if (std::fabs(incr) < tol_c_cell) {
incr = cmax_cell - c_cell;
c_cell = cmax_cell;
} else {
incr = (cmax_cell - c_cell)/2.0;
c_cell = (cmax_cell + c_cell)/2.0;
}
}
// if (s_cell < 0.) {
// increment[2*cell + 0] = increment[2*cell + 0] - s_cell;
// s_cell = 0.;
// } else if (s_cell > 1.) {
// increment[2*cell + 0] = increment[2*cell + 0] - s_cell + 1.;
// s_cell = 1.;
// }
// if (c_cell < 0) {
// increment[2*cell + 1] = increment[2*cell + 1] - c_cell;
// c_cell = 0.;
// } else if (c_cell > cmax_cell) {
// increment[2*cell + 1] = increment[2*cell + 1] - c_cell + cmax_cell;
// c_cell = cmax_cell;
// }
}
const double maxelem = *std::max_element(increment.begin(), increment.end());
const double minelem = *std::min_element(increment.begin(), increment.end());
max_delta = std::max(maxelem, -minelem);
@@ -156,10 +228,10 @@ namespace Opm
THROW("Failed to converge!");
}
// Finalize.
// model_.finishIteration(); //
// fmodel_.finishIteration(); //
// finishStep() writes to state, which holds s by reference.
// This will update the entire grid's state... cmax is updated here.
model_.finishStep(grid_, sys.vector().solution(), state);
fmodel_.finishStep(grid_, sys.vector().solution(), state);
}
@@ -168,8 +240,8 @@ 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 GravityColumnSolverPolymer<Model>::solveSingleColumn(const std::vector<int>& column_cells,
template <class FluxModel, class Model>
void GravityColumnSolverPolymer<FluxModel, Model>::solveSingleColumn(const std::vector<int>& column_cells,
const double dt,
std::vector<double>& s,
std::vector<double>& c,
@@ -215,7 +287,7 @@ namespace Opm
F.assign(2, 0.);
dFd1.assign(4, 0.);
dFd2.assign(4, 0.);
model_.fluxConnection(state, grid_, dt, cell, face, &F[0], &dFd1[0], &dFd2[0]);
fmodel_.fluxConnection(state, grid_, dt, cell, face, &F[0], &dFd1[0], &dFd2[0]);
if (c1 == prev_cell || c2 == prev_cell) {
hm[bmc(2*ci + 0, 2*(ci - 1) + 0)] += dFd2[0];
hm[bmc(2*ci + 0, 2*(ci - 1) + 1)] += dFd2[1];
@@ -239,7 +311,7 @@ namespace Opm
}
F.assign(2, 0.);
dF.assign(4, 0.);
model_.accumulation(grid_, cell, &F[0], &dF[0]);
fmodel_.accumulation(grid_, cell, &F[0], &dF[0]);
hm[bmc(2*ci + 0, 2*ci + 0)] += dF[0];
hm[bmc(2*ci + 0, 2*ci + 1)] += dF[1];
hm[bmc(2*ci + 1, 2*ci + 0)] += dF[2];