diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index a061a4478..d796bf541 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -28,6 +28,7 @@ list (APPEND MAIN_SOURCE_FILES opm/autodiff/BlackoilPropsAd.cpp opm/autodiff/BlackoilPropsAdInterface.cpp + opm/autodiff/NewtonIterationBlackoilSimple.cpp opm/autodiff/GridHelpers.cpp opm/autodiff/ImpesTPFAAD.cpp opm/autodiff/SimulatorCompressibleAd.cpp @@ -103,6 +104,9 @@ list (APPEND PUBLIC_HEADER_FILES opm/autodiff/ImpesTPFAAD.hpp opm/autodiff/FullyImplicitBlackoilSolver.hpp opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp + opm/autodiff/NewtonIterationBlackoilInterface.hpp + opm/autodiff/NewtonIterationBlackoilSimple.hpp + opm/autodiff/LinearisedBlackoilResidual.hpp opm/autodiff/SimulatorCompressibleAd.hpp opm/autodiff/SimulatorFullyImplicitBlackoil.hpp opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp diff --git a/examples/sim_fibo_ad.cpp b/examples/sim_fibo_ad.cpp index 243926788..890ca5bee 100644 --- a/examples/sim_fibo_ad.cpp +++ b/examples/sim_fibo_ad.cpp @@ -37,6 +37,7 @@ #include #include +#include #include #include @@ -141,6 +142,7 @@ try // Linear solver. LinearSolverFactory linsolver(param); + NewtonIterationBlackoilSimple fis_solver(linsolver); // Write parameters used for later reference. bool output = param.getDefault("output", true); @@ -210,7 +212,7 @@ try *new_props, rock_comp->isActive() ? rock_comp.get() : 0, wells, - linsolver, + fis_solver, grav); SimulatorReport episodeReport = simulator.run(simtimer, state, well_state); diff --git a/examples/sim_fibo_ad_cp.cpp b/examples/sim_fibo_ad_cp.cpp index 63856026e..dfa96e3fb 100644 --- a/examples/sim_fibo_ad_cp.cpp +++ b/examples/sim_fibo_ad_cp.cpp @@ -56,6 +56,7 @@ #include #include +#include #include #include @@ -189,6 +190,7 @@ try // Linear solver. LinearSolverFactory linsolver(param); + NewtonIterationBlackoilSimple fis_solver(linsolver); // Write parameters used for later reference. bool output = param.getDefault("output", true); @@ -263,7 +265,7 @@ try *new_props, rock_comp->isActive() ? rock_comp.get() : 0, wells, - linsolver, + fis_solver, grav); SimulatorReport episodeReport = simulator.run(simtimer, state, well_state); diff --git a/examples/test_implicit_ad.cpp b/examples/test_implicit_ad.cpp index 650366366..563100adf 100644 --- a/examples/test_implicit_ad.cpp +++ b/examples/test_implicit_ad.cpp @@ -24,6 +24,7 @@ #include #include #include +#include #include #include @@ -103,8 +104,9 @@ try Opm::DerivedGeology geo(*g, props, grav); Opm::LinearSolverFactory linsolver(param); + Opm::NewtonIterationBlackoilSimple fis_solver(linsolver); - Opm::FullyImplicitBlackoilSolver solver(*g, props, geo, 0, *wells, linsolver); + Opm::FullyImplicitBlackoilSolver solver(*g, props, geo, 0, *wells, fis_solver); Opm::BlackoilState state; initStateBasic(*g, props0, param, 0.0, state); diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.hpp b/opm/autodiff/FullyImplicitBlackoilSolver.hpp index a21757cf4..1ded4c579 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.hpp @@ -23,6 +23,8 @@ #include #include #include +#include +#include struct UnstructuredGrid; struct Wells; @@ -31,7 +33,7 @@ namespace Opm { class DerivedGeology; class RockCompressibility; - class LinearSolverInterface; + class NewtonIterationBlackoilInterface; class BlackoilState; class WellStateFullyImplicitBlackoil; @@ -65,7 +67,7 @@ namespace Opm { const DerivedGeology& geo , const RockCompressibility* rock_comp_props, const Wells& wells, - const LinearSolverInterface& linsolver); + const NewtonIterationBlackoilInterface& linsolver); /// Take a single forward step, modifiying /// state.pressure() @@ -126,7 +128,7 @@ namespace Opm { const DerivedGeology& geo_; const RockCompressibility* rock_comp_props_; const Wells& wells_; - const LinearSolverInterface& linsolver_; + const NewtonIterationBlackoilInterface& linsolver_; // For each canonical phase -> true if active const std::vector active_; // Size = # active faces. Maps active -> canonical phase indices. @@ -140,14 +142,7 @@ namespace Opm { std::vector phaseCondition_; V well_perforation_pressure_diffs_; // Diff to bhp for each well perforation. - // 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 diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp index 110ffc744..4601e9adb 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp @@ -203,7 +203,7 @@ namespace { const DerivedGeology& geo , const RockCompressibility* rock_comp_props, const Wells& wells, - const LinearSolverInterface& linsolver) + const NewtonIterationBlackoilInterface& linsolver) : grid_ (grid) , fluid_ (fluid) , geo_ (geo) @@ -662,13 +662,13 @@ namespace { // std::cout << "===== rq_[" << phase << "].mflux = \n" << std::endl; // std::cout << rq_[phase].mflux; - residual_.mass_balance[ phaseIdx ] = + residual_.material_balance_eq[ phaseIdx ] = pvdt*(rq_[phaseIdx].accum[1] - rq_[phaseIdx].accum[0]) + ops_.div*rq_[phaseIdx].mflux; // DUMP(ops_.div*rq_[phase].mflux); - // DUMP(residual_.mass_balance[phase]); + // DUMP(residual_.material_balance_eq[phase]); } // -------- Extra (optional) rs and rv contributions to the mass balance equations -------- @@ -682,16 +682,16 @@ namespace { rq_[po].head.value()); const ADB rs_face = upwindOil.select(state.rs); - residual_.mass_balance[ Gas ] += ops_.div * (rs_face * rq_[po].mflux); + residual_.material_balance_eq[ Gas ] += ops_.div * (rs_face * rq_[po].mflux); const int pg = fluid_.phaseUsage().phase_pos[ Gas ]; const UpwindSelector upwindGas(grid_, ops_, rq_[pg].head.value()); const ADB rv_face = upwindGas.select(state.rv); - residual_.mass_balance[ Oil ] += ops_.div * (rv_face * rq_[pg].mflux); + residual_.material_balance_eq[ Oil ] += ops_.div * (rv_face * rq_[pg].mflux); - // DUMP(residual_.mass_balance[ Gas ]); + // DUMP(residual_.material_balance_eq[ Gas ]); } @@ -870,7 +870,7 @@ namespace { // Add well contributions to mass balance equations for (int phase = 0; phase < np; ++phase) { - residual_.mass_balance[phase] -= superset(cq_s[phase],well_cells,nc); + residual_.material_balance_eq[phase] -= superset(cq_s[phase],well_cells,nc); } // Add WELL EQUATIONS @@ -1154,7 +1154,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]; } if (active_[Gas] && active_[Oil]) { const int oilpos = pu.phase_pos[Oil]; @@ -1164,8 +1164,8 @@ namespace { well_rates_all += superset(wops_.p2w * (well_perf_rates[oilpos]*rs_perf), Span(nw, 1, gaspos*nw), nw*np); well_rates_all += superset(wops_.p2w * (well_perf_rates[gaspos]*rv_perf), Span(nw, 1, oilpos*nw), nw*np); // DUMP(well_contribs[gaspos] + well_contribs[oilpos]*state.rs); - residual_.mass_balance[gaspos] += well_contribs[oilpos]*state.rs; - residual_.mass_balance[oilpos] += well_contribs[gaspos]*state.rv; + residual_.material_balance_eq[gaspos] += well_contribs[oilpos]*state.rs; + residual_.material_balance_eq[oilpos] += well_contribs[gaspos]*state.rv; } // Set the well flux equation @@ -1180,40 +1180,7 @@ namespace { template V FullyImplicitBlackoilSolver::solveJacobianSystem() const { - const int np = fluid_.numPhases(); - ADB mass_res = residual_.mass_balance[0]; - for (int phase = 1; phase < np; ++phase) { - mass_res = vertcat(mass_res, residual_.mass_balance[phase]); - } - 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()); - /* - std::ofstream filestream("matrix.out"); - filestream << matr; - filestream.close(); - std::ofstream filestream2("sol.out"); - filestream2 << dx; - filestream2.close(); - std::ofstream filestream3("r.out"); - filestream3 << total_residual.value(); - filestream3.close(); */ - - - if (!rep.converged) { - OPM_THROW(std::runtime_error, - "FullyImplicitBlackoilSolver::solveJacobianSystem(): " - "Linear solver convergence failure."); - } - return dx; + return linsolver_.computeNewtonIncrement(residual_); } @@ -1630,10 +1597,9 @@ namespace { FullyImplicitBlackoilSolver::residualNorm() const { double globalNorm = 0; - std::vector::const_iterator quantityIt = residual_.mass_balance.begin(); - const std::vector::const_iterator endQuantityIt = residual_.mass_balance.end(); - for (; quantityIt != endQuantityIt; ++quantityIt) - { + std::vector::const_iterator quantityIt = residual_.material_balance_eq.begin(); + const std::vector::const_iterator endQuantityIt = residual_.material_balance_eq.end(); + for (; quantityIt != endQuantityIt; ++quantityIt) { const double quantityResid = (*quantityIt).value().matrix().norm(); if (!std::isfinite(quantityResid)) { OPM_THROW(Opm::NumericalProblem, diff --git a/opm/autodiff/LinearisedBlackoilResidual.hpp b/opm/autodiff/LinearisedBlackoilResidual.hpp new file mode 100644 index 000000000..f37983f85 --- /dev/null +++ b/opm/autodiff/LinearisedBlackoilResidual.hpp @@ -0,0 +1,70 @@ +/* + Copyright 2014 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_LINEARISEDBLACKOILRESIDUAL_HEADER_INCLUDED +#define OPM_LINEARISEDBLACKOILRESIDUAL_HEADER_INCLUDED + +#include + +namespace Opm +{ + + /// Residual structure of the fully implicit solver. + /// All equations are given as AD types, with multiple + /// jacobian blocks corresponding to the primary unknowns. The + /// primary unknowns are for a three-phase simulation, in order: + /// p (pressure) + /// sw (water saturation) + /// xvar (gas saturation, gas-oil ratio or oil-gas ratio) + /// qs (well outflows by well and phase) + /// bhp (bottom hole pressures) + /// In the above, the xvar variable will have a different + /// meaning from cell to cell, corresponding to the state in + /// that cell (saturated, undersaturated oil or undersaturated + /// gas). In a two-phase simulation, either sw or xvar is not + /// used, depending on which phase is missing. + /// + /// Note: this class is strongly coupled to the class + /// FullyImplicitBlackoilSolver, and is separated from that + /// class to facilitate the development of linear solver + /// strategies outside that class. + struct LinearisedBlackoilResidual { + /// A type alias for the automatic differentiation type. + typedef AutoDiffBlock ADB; + /// The material_balance_eq vector has one element for each + /// active phase, each of which has size equal to the number + /// of cells. Each material balance equation is given in terms + /// of surface volumes (in SI units, that is standard m^3). + std::vector material_balance_eq; + /// The well_flux_eq has size equal to the number of wells + /// times the number of phases. It contains the well flow + /// equations, relating the total well flows to + /// bottom-hole pressures and reservoir conditions. + ADB well_flux_eq; + /// The well_eq has size equal to the number of wells. It + /// contains the well control equations, that is for each + /// well either a rate specification or bottom hole + /// pressure specification. + ADB well_eq; + }; + +} // namespace Opm + + +#endif // OPM_LINEARISEDBLACKOILRESIDUAL_HEADER_INCLUDED diff --git a/opm/autodiff/NewtonIterationBlackoilInterface.hpp b/opm/autodiff/NewtonIterationBlackoilInterface.hpp new file mode 100644 index 000000000..277b36d9b --- /dev/null +++ b/opm/autodiff/NewtonIterationBlackoilInterface.hpp @@ -0,0 +1,46 @@ +/* + Copyright 2014 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_NEWTONITERATIONBLACKOILINTERFACE_HEADER_INCLUDED +#define OPM_NEWTONITERATIONBLACKOILINTERFACE_HEADER_INCLUDED + +#include + +namespace Opm +{ + + /// Interface class for (linear) solvers for the fully implicit black-oil system. + class NewtonIterationBlackoilInterface + { + public: + /// Return type for linearSolve(). A simple, non-ad vector type. + typedef LinearisedBlackoilResidual::ADB::V SolutionVector; + + /// Solve the linear system Ax = b, with A being the + /// combined derivative matrix of the residual and b + /// being the residual itself. + /// \param[in] residual residual object containing A and b. + /// \return the solution x + virtual SolutionVector computeNewtonIncrement(const LinearisedBlackoilResidual& residual) const = 0; + }; + +} // namespace Opm + + +#endif // OPM_NEWTONITERATIONBLACKOILINTERFACE_HEADER_INCLUDED diff --git a/opm/autodiff/NewtonIterationBlackoilSimple.cpp b/opm/autodiff/NewtonIterationBlackoilSimple.cpp new file mode 100644 index 000000000..92e885903 --- /dev/null +++ b/opm/autodiff/NewtonIterationBlackoilSimple.cpp @@ -0,0 +1,69 @@ +/* + Copyright 2014 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#include + +#include +#include +#include + +namespace Opm +{ + + /// Construct a system solver. + /// \param[in] linsolver linear solver to use + NewtonIterationBlackoilSimple::NewtonIterationBlackoilSimple(const LinearSolverInterface& linsolver) + : linsolver_(linsolver) + { + } + + /// Solve the linear system Ax = b, with A being the + /// combined derivative matrix of the residual and b + /// being the residual itself. + /// \param[in] residual residual object containing A and b. + /// \return the solution x + NewtonIterationBlackoilSimple::SolutionVector + NewtonIterationBlackoilSimple::computeNewtonIncrement(const LinearisedBlackoilResidual& residual) const + { + typedef LinearisedBlackoilResidual::ADB ADB; + const int np = residual.material_balance_eq.size(); + ADB mass_res = residual.material_balance_eq[0]; + for (int phase = 1; phase < np; ++phase) { + mass_res = vertcat(mass_res, residual.material_balance_eq[phase]); + } + const ADB well_res = vertcat(residual.well_flux_eq, residual.well_eq); + const ADB total_residual = collapseJacs(vertcat(mass_res, well_res)); + + const Eigen::SparseMatrix matr = total_residual.derivative()[0]; + + SolutionVector dx(SolutionVector::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, + "FullyImplicitBlackoilSolver::solveJacobianSystem(): " + "Linear solver convergence failure."); + } + return dx; + } + +} // namespace Opm + diff --git a/opm/autodiff/NewtonIterationBlackoilSimple.hpp b/opm/autodiff/NewtonIterationBlackoilSimple.hpp new file mode 100644 index 000000000..ce17ef8b1 --- /dev/null +++ b/opm/autodiff/NewtonIterationBlackoilSimple.hpp @@ -0,0 +1,55 @@ +/* + Copyright 2014 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_NEWTONITERATIONBLACKOILSIMPLE_HEADER_INCLUDED +#define OPM_NEWTONITERATIONBLACKOILSIMPLE_HEADER_INCLUDED + + +#include +#include + +namespace Opm +{ + + /// This class solves the fully implicit black-oil system by + /// simply concatenating the Jacobian matrices and passing the + /// resulting system to a linear solver. The linear solver used + /// can be passed in as a constructor argument. + class NewtonIterationBlackoilSimple : public NewtonIterationBlackoilInterface + { + public: + /// Construct a system solver. + /// \param[in] linsolver linear solver to use + NewtonIterationBlackoilSimple(const LinearSolverInterface& linsolver); + + /// Solve the system of linear equations Ax = b, with A being the + /// combined derivative matrix of the residual and b + /// being the residual itself. + /// \param[in] residual residual object containing A and b. + /// \return the solution x + virtual SolutionVector computeNewtonIncrement(const LinearisedBlackoilResidual& residual) const; + + private: + const LinearSolverInterface& linsolver_; + }; + +} // namespace Opm + + +#endif // OPM_NEWTONITERATIONBLACKOILSIMPLE_HEADER_INCLUDED diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoil.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoil.hpp index 5f17e969c..4551dfb8b 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoil.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoil.hpp @@ -33,7 +33,7 @@ namespace Opm class BlackoilPropsAdInterface; class RockCompressibility; class WellsManager; - class LinearSolverInterface; + class NewtonIterationBlackoilInterface; class SimulatorTimer; class BlackoilState; class WellStateFullyImplicitBlackoil; @@ -73,7 +73,7 @@ namespace Opm BlackoilPropsAdInterface& props, const RockCompressibility* rock_comp_props, WellsManager& wells_manager, - LinearSolverInterface& linsolver, + NewtonIterationBlackoilInterface& linsolver, const double* gravity); /// Run the simulation. diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp index 84c072508..b1862ee01 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp @@ -65,7 +65,7 @@ namespace Opm BlackoilPropsAdInterface& props, const RockCompressibility* rock_comp_props, WellsManager& wells_manager, - LinearSolverInterface& linsolver, + NewtonIterationBlackoilInterface& linsolver, const double* gravity); SimulatorReport run(SimulatorTimer& timer, @@ -106,7 +106,7 @@ namespace Opm BlackoilPropsAdInterface& props, const RockCompressibility* rock_comp_props, WellsManager& wells_manager, - LinearSolverInterface& linsolver, + NewtonIterationBlackoilInterface& linsolver, const double* gravity) { @@ -191,7 +191,7 @@ namespace Opm BlackoilPropsAdInterface& props, const RockCompressibility* rock_comp_props, WellsManager& wells_manager, - LinearSolverInterface& linsolver, + NewtonIterationBlackoilInterface& linsolver, const double* gravity) : grid_(grid), props_(props),