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

use Newton iteration solver interface from opm-autodiff, prepare for CPR supporting.

Browse Source
This commit is contained in:
Liu Ming 2014-09-26 14:54:57 +08:00
parent ca70d67d83
commit 979c2dc0f7
5 changed files with 32 additions and 46 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
examples/sim_poly_fi2p_comp_ad.cpp
View File

@ -41,6 +41,8 @@
#include <opm/core/props/rock/RockCompressibility.hpp>
#include <opm/core/linalg/LinearSolverFactory.hpp>
#include <opm/autodiff/NewtonIterationBlackoilSimple.hpp>
#include <opm/autodiff/NewtonIterationBlackoilCPR.hpp>
#include <opm/polymer/PolymerBlackoilState.hpp>
#include <opm/core/simulator/WellState.hpp>
@ -156,8 +158,13 @@ try
bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0);
const double *grav = use_gravity ? &gravity[0] : 0;
// Linear solver.
LinearSolverFactory linsolver(param);
// Solver for Newton iterations.
std::unique_ptr<NewtonIterationBlackoilInterface> fis_solver;
if (param.getDefault("use_cpr", true)) {
fis_solver.reset(new NewtonIterationBlackoilCPR(param));
} else {
fis_solver.reset(new NewtonIterationBlackoilSimple(param));
}
// Write parameters used for later reference.
bool output = param.getDefault("output", true);
@ -234,7 +241,7 @@ try
rock_comp->isActive() ? rock_comp.get() : 0,
wells,
*polymer_inflow,
linsolver,
*fis_solver,
grav);
if (reportStepIdx == 0) {
warnIfUnusedParams(param);

40
opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp
View File

@ -168,7 +168,7 @@ namespace {
const RockCompressibility* rock_comp_props,
const PolymerPropsAd& polymer_props_ad,
const Wells& wells,
const LinearSolverInterface& linsolver)
const NewtonIterationBlackoilInterface& linsolver)
: grid_ (grid)
, fluid_ (fluid)
, geo_ (geo)
@ -212,7 +212,7 @@ namespace {
assemble(dt, x, xw, polymer_inflow, src);
const double r0 = residualNorm();
const double r_polymer = residual_.mass_balance[2].value().matrix().lpNorm<Eigen::Infinity>();
const double r_polymer = residual_.material_balance_eq[2].value().matrix().lpNorm<Eigen::Infinity>();
int it = 0;
std::cout << "\nIteration Residual Polymer Res\n"
<< std::setw(9) << it << std::setprecision(9)
@ -227,7 +227,7 @@ namespace {
const double r = residualNorm();
const double rr_polymer = residual_.mass_balance[2].value().matrix().lpNorm<Eigen::Infinity>();
const double rr_polymer = residual_.material_balance_eq[2].value().matrix().lpNorm<Eigen::Infinity>();
resTooLarge = (r > atol) && (r > rtol*r0);
it += 1;
@ -521,11 +521,11 @@ namespace {
const ADB krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration, cmax, kr[0], state.saturation[0]);
const ADB mc = computeMc(state);
computeMassFlux(trans, mc, kr[1], krw_eff, state);
residual_.mass_balance[0] = pvdt*(rq_[0].accum[1] - rq_[0].accum[0])
residual_.material_balance_eq[0] = pvdt*(rq_[0].accum[1] - rq_[0].accum[0])
+ ops_.div*rq_[0].mflux;
residual_.mass_balance[1] = pvdt*(rq_[1].accum[1] - rq_[1].accum[0])
residual_.material_balance_eq[1] = pvdt*(rq_[1].accum[1] - rq_[1].accum[0])
+ ops_.div*rq_[1].mflux;
residual_.mass_balance[2] = pvdt*(rq_[2].accum[1] - rq_[2].accum[0]) //+ cell / dt * (rq_[2].ads[1] - rq_[2].ads[0])
residual_.material_balance_eq[2] = pvdt*(rq_[2].accum[1] - rq_[2].accum[0]) //+ cell / dt * (rq_[2].ads[1] - rq_[2].ads[0])
+ ops_.div*rq_[2].mflux;
// -------- Extra (optional) sg or rs equation, and rs contributions to the mass balance equations --------
@ -620,7 +620,7 @@ namespace {
// const ADB well_contrib = superset(perf_flux*perf_b, well_cells, nc);
well_contribs[phase] = superset(perf_flux*perf_b, well_cells, nc);
// DUMP(well_contribs[phase]);
residual_.mass_balance[phase] += well_contribs[phase];
residual_.material_balance_eq[phase] += well_contribs[phase];
for (int cell = 0; cell < nc; ++cell) {
src[cell] += well_contribs[phase].value()[cell];
}
@ -634,7 +634,7 @@ namespace {
const V poly_in_c = poly_in_perf;// * poly_mc_cell;
const V poly_mc = producer.select(poly_mc_cell, poly_in_c);
residual_.mass_balance[2] += superset(well_perf_rates[0] * poly_mc, well_cells, nc);
residual_.material_balance_eq[2] += superset(well_perf_rates[0] * poly_mc, well_cells, nc);
// Set the well flux equation
residual_.well_flux_eq = state.qs + well_rates_all;
// DUMP(residual_.well_flux_eq);
@ -675,25 +675,7 @@ namespace {
V FullyImplicitCompressiblePolymerSolver::solveJacobianSystem() const
{
ADB mass_res = vertcat(residual_.mass_balance[0], residual_.mass_balance[1]);
mass_res = vertcat(mass_res, residual_.mass_balance[2]);
const ADB well_res = vertcat(residual_.well_flux_eq, residual_.well_eq);
const ADB total_residual = collapseJacs(vertcat(mass_res, well_res));
// DUMP(total_residual);
const Eigen::SparseMatrix<double, Eigen::RowMajor> matr = total_residual.derivative()[0];
V dx(V::Zero(total_residual.size()));
Opm::LinearSolverInterface::LinearSolverReport rep
= linsolver_.solve(matr.rows(), matr.nonZeros(),
matr.outerIndexPtr(), matr.innerIndexPtr(), matr.valuePtr(),
total_residual.value().data(), dx.data());
if (!rep.converged) {
OPM_THROW(std::runtime_error,
"FullyImplicitCompressiblePolymerSolver::solveJacobianSystem(): "
"Linear solver convergence failure.");
}
return dx;
return linsolver_.computeNewtonIncrement(residual_);
}
@ -898,8 +880,8 @@ namespace {
{
double r = 0;
for (std::vector<ADB>::const_iterator
b = residual_.mass_balance.begin(),
e = residual_.mass_balance.end();
b = residual_.material_balance_eq.begin(),
e = residual_.material_balance_eq.end();
b != e; ++b)
{
r = std::max(r, (*b).value().matrix().lpNorm<Eigen::Infinity>());

15
opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp
View File

@ -24,6 +24,8 @@
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
#include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
#include <opm/autodiff/LinearisedBlackoilResidual.hpp>
#include <opm/polymer/PolymerProperties.hpp>
#include <opm/polymer/fullyimplicit/PolymerPropsAd.hpp>
@ -34,7 +36,7 @@ namespace Opm {
class DerivedGeology;
class RockCompressibility;
class LinearSolverInterface;
class NewtonIterationBlackoilInterface;
class PolymerBlackoilState;
class WellState;
@ -65,7 +67,7 @@ namespace Opm {
const RockCompressibility* rock_comp_props,
const PolymerPropsAd& polymer_props_ad,
const Wells& wells,
const LinearSolverInterface& linsolver);
const NewtonIterationBlackoilInterface& linsolver);
/// Take a single forward step, modifiying
/// state.pressure()
@ -129,7 +131,7 @@ namespace Opm {
const RockCompressibility* rock_comp_props_;
const PolymerPropsAd& polymer_props_ad_;
const Wells& wells_;
const LinearSolverInterface& linsolver_;
const NewtonIterationBlackoilInterface& linsolver_;
const std::vector<int> cells_; // All grid cells
HelperOps ops_;
const WellOps wops_;
@ -140,12 +142,7 @@ namespace Opm {
// The mass_balance vector has one element for each active phase,
// each of which has size equal to the number of cells.
// The well_eq has size equal to the number of wells.
struct {
std::vector<ADB> mass_balance;
ADB well_flux_eq;
ADB well_eq;
} residual_;
LinearisedBlackoilResidual residual_;
// Private methods.
SolutionState
constantState(const PolymerBlackoilState& x,

6
opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.cpp
View File

@ -92,7 +92,7 @@ namespace Opm
const RockCompressibility* rock_comp_props,
WellsManager& wells_manager,
PolymerInflowInterface& polymer_inflow,
LinearSolverInterface& linsolver,
NewtonIterationBlackoilInterface& linsolver,
const double* gravity);
SimulatorReport run(SimulatorTimer& timer,
@ -138,7 +138,7 @@ namespace Opm
const RockCompressibility* rock_comp_props,
WellsManager& wells_manager,
PolymerInflowInterface& polymer_inflow,
LinearSolverInterface& linsolver,
NewtonIterationBlackoilInterface& linsolver,
const double* gravity)
{
@ -169,7 +169,7 @@ namespace Opm
const RockCompressibility* rock_comp_props,
WellsManager& wells_manager,
PolymerInflowInterface& polymer_inflow,
LinearSolverInterface& linsolver,
NewtonIterationBlackoilInterface& linsolver,
const double* gravity)
: grid_(grid),
props_(props),

4
opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.hpp
View File

@ -34,7 +34,7 @@ namespace Opm
class RockCompressibility;
class DerivedGeology;
class WellsManager;
class LinearSolverInterface;
class NewtonIterationBlackoilInterface;
class SimulatorTimer;
class PolymerBlackoilState;
class WellState;
@ -78,7 +78,7 @@ namespace Opm
const RockCompressibility* rock_comp_props,
WellsManager& wells_manager,
PolymerInflowInterface& polymer_inflow,
LinearSolverInterface& linsolver,
NewtonIterationBlackoilInterface& linsolver,
const double* gravity);
/// Run the simulation.