From 979c2dc0f772866065e85ce0886ad5b3e25c800d Mon Sep 17 00:00:00 2001 From: Liu Ming Date: Fri, 26 Sep 2014 14:54:57 +0800 Subject: [PATCH] use Newton iteration solver interface from opm-autodiff, prepare for CPR supporting. --- examples/sim_poly_fi2p_comp_ad.cpp | 13 ++++-- ...FullyImplicitCompressiblePolymerSolver.cpp | 40 +++++-------------- ...FullyImplicitCompressiblePolymerSolver.hpp | 15 +++---- ...ulatorFullyImplicitCompressiblePolymer.cpp | 6 +-- ...ulatorFullyImplicitCompressiblePolymer.hpp | 4 +- 5 files changed, 32 insertions(+), 46 deletions(-) diff --git a/examples/sim_poly_fi2p_comp_ad.cpp b/examples/sim_poly_fi2p_comp_ad.cpp index bea3f57b4..4f5198cf7 100644 --- a/examples/sim_poly_fi2p_comp_ad.cpp +++ b/examples/sim_poly_fi2p_comp_ad.cpp @@ -41,6 +41,8 @@ #include #include +#include +#include #include #include @@ -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 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); diff --git a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp index 41721d932..e5a5dfd22 100644 --- a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp +++ b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp @@ -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(); + const double r_polymer = residual_.material_balance_eq[2].value().matrix().lpNorm(); 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(); + const double rr_polymer = residual_.material_balance_eq[2].value().matrix().lpNorm(); 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 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::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()); diff --git a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp index b93fe206e..393eabbbb 100644 --- a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp +++ b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp @@ -24,6 +24,8 @@ #include #include #include +#include +#include #include #include @@ -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 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 mass_balance; - ADB well_flux_eq; - ADB well_eq; - } residual_; - + LinearisedBlackoilResidual residual_; // Private methods. SolutionState constantState(const PolymerBlackoilState& x, diff --git a/opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.cpp b/opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.cpp index 1c58031c4..6152474af 100644 --- a/opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.cpp +++ b/opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.cpp @@ -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), diff --git a/opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.hpp b/opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.hpp index 515268844..11016a262 100644 --- a/opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.hpp +++ b/opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.hpp @@ -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.