diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index aaa414c8d..3798e759d 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -29,6 +29,7 @@ list (APPEND MAIN_SOURCE_FILES opm/autodiff/BlackoilPropsAd.cpp opm/autodiff/BlackoilPropsAdInterface.cpp opm/autodiff/FullyImplicitBlackoilSolver.cpp + opm/autodiff/FullyImplicitSystemSolverSimple.cpp opm/autodiff/ImpesTPFAAD.cpp opm/autodiff/SimulatorCompressibleAd.cpp opm/autodiff/SimulatorFullyImplicitBlackoil.cpp @@ -100,7 +101,10 @@ list (APPEND PUBLIC_HEADER_FILES opm/autodiff/BlackoilPropsAdInterface.hpp opm/autodiff/GeoProps.hpp opm/autodiff/ImpesTPFAAD.hpp + opm/autodiff/FullyImplicitBlackoilResidual.hpp opm/autodiff/FullyImplicitBlackoilSolver.hpp + opm/autodiff/FullyImplicitBlackoilSystemSolverInterface.hpp + opm/autodiff/FullyImplicitBlackoilSystemSolverSimple.hpp opm/autodiff/SimulatorCompressibleAd.hpp opm/autodiff/SimulatorFullyImplicitBlackoil.hpp opm/autodiff/SimulatorIncompTwophaseAd.hpp diff --git a/examples/sim_fibo_ad.cpp b/examples/sim_fibo_ad.cpp index 348f983cb..a331c819f 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); + FullyImplicitSystemSolverSimple 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/test_implicit_ad.cpp b/examples/test_implicit_ad.cpp index 11558b396..48c586d40 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::FullyImplicitSystemSolverSimple 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/FullyImplicitBlackoilResidual.hpp b/opm/autodiff/FullyImplicitBlackoilResidual.hpp new file mode 100644 index 000000000..757270212 --- /dev/null +++ b/opm/autodiff/FullyImplicitBlackoilResidual.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_FULLYIMPLICITBLACKOILRESIDUAL_HEADER_INCLUDED +#define OPM_FULLYIMPLICITBLACKOILRESIDUAL_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 face 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 FullyImplicitBlackoilResidual { + /// 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. + 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_FULLYIMPLICITBLACKOILRESIDUAL_HEADER_INCLUDED diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.cpp b/opm/autodiff/FullyImplicitBlackoilSolver.cpp index d36d47590..6201b9be4 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.cpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.cpp @@ -195,7 +195,7 @@ namespace { const DerivedGeology& geo , const RockCompressibility* rock_comp_props, const Wells& wells, - const LinearSolverInterface& linsolver) + const FullyImplicitSystemSolverInterface& linsolver) : grid_ (grid) , fluid_ (fluid) , geo_ (geo) @@ -643,13 +643,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 -------- @@ -663,16 +663,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 ]); } @@ -850,7 +850,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 @@ -1133,7 +1133,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]; @@ -1143,8 +1143,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 @@ -1158,40 +1158,7 @@ namespace { 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_.linearSolve(residual_); } @@ -1599,10 +1566,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/FullyImplicitBlackoilSolver.hpp b/opm/autodiff/FullyImplicitBlackoilSolver.hpp index 15c27e35a..873308fbd 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 FullyImplicitSystemSolverInterface; class BlackoilState; class WellStateFullyImplicitBlackoil; @@ -62,7 +64,7 @@ namespace Opm { const DerivedGeology& geo , const RockCompressibility* rock_comp_props, const Wells& wells, - const LinearSolverInterface& linsolver); + const FullyImplicitSystemSolverInterface& linsolver); /// Take a single forward step, modifiying /// state.pressure() @@ -123,7 +125,7 @@ namespace Opm { const DerivedGeology& geo_; const RockCompressibility* rock_comp_props_; const Wells& wells_; - const LinearSolverInterface& linsolver_; + const FullyImplicitSystemSolverInterface& linsolver_; // For each canonical phase -> true if active const std::vector active_; // Size = # active faces. Maps active -> canonical phase indices. @@ -137,14 +139,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_; + FullyImplicitBlackoilResidual residual_; // Private methods. SolutionState diff --git a/opm/autodiff/FullyImplicitSystemSolverInterface.hpp b/opm/autodiff/FullyImplicitSystemSolverInterface.hpp new file mode 100644 index 000000000..d3b4b980d --- /dev/null +++ b/opm/autodiff/FullyImplicitSystemSolverInterface.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_FULLYIMPLICITSYSTEMSOLVERINTERFACE_HEADER_INCLUDED +#define OPM_FULLYIMPLICITSYSTEMSOLVERINTERFACE_HEADER_INCLUDED + +#include + +namespace Opm +{ + + /// Interface class for (linear) solvers for the fully implicit black-oil system. + class FullyImplicitSystemSolverInterface + { + public: + /// Return type for linearSolve(). A simple, non-ad vector type. + typedef FullyImplicitBlackoilResidual::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 linearSolve(const FullyImplicitBlackoilResidual& residual) const = 0; + }; + +} // namespace Opm + + +#endif // OPM_FULLYIMPLICITSYSTEMSOLVERINTERFACE_HEADER_INCLUDED diff --git a/opm/autodiff/FullyImplicitSystemSolverSimple.cpp b/opm/autodiff/FullyImplicitSystemSolverSimple.cpp new file mode 100644 index 000000000..ad5a5ec0b --- /dev/null +++ b/opm/autodiff/FullyImplicitSystemSolverSimple.cpp @@ -0,0 +1,68 @@ +/* + 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 + +namespace Opm +{ + + /// Construct a system solver. + /// \param[in] linsolver linear solver to use + FullyImplicitSystemSolverSimple::FullyImplicitSystemSolverSimple(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 + FullyImplicitSystemSolverSimple::SolutionVector + FullyImplicitSystemSolverSimple::linearSolve(const FullyImplicitBlackoilResidual& residual) const + { + typedef FullyImplicitBlackoilResidual::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/FullyImplicitSystemSolverSimple.hpp b/opm/autodiff/FullyImplicitSystemSolverSimple.hpp new file mode 100644 index 000000000..2c2fadd50 --- /dev/null +++ b/opm/autodiff/FullyImplicitSystemSolverSimple.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_FULLYIMPLICITSYSTEMSOLVERSIMPLE_HEADER_INCLUDED +#define OPM_FULLYIMPLICITSYSTEMSOLVERSIMPLE_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 FullyImplicitSystemSolverSimple : public FullyImplicitSystemSolverInterface + { + public: + /// Construct a system solver. + /// \param[in] linsolver linear solver to use + FullyImplicitSystemSolverSimple(const LinearSolverInterface& 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 + virtual SolutionVector linearSolve(const FullyImplicitBlackoilResidual& residual) const; + + private: + const LinearSolverInterface& linsolver_; + }; + +} // namespace Opm + + +#endif // OPM_FULLYIMPLICITSYSTEMSOLVERSIMPLE_HEADER_INCLUDED diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoil.cpp b/opm/autodiff/SimulatorFullyImplicitBlackoil.cpp index a36ae887e..cedfe9b71 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoil.cpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoil.cpp @@ -70,7 +70,7 @@ namespace Opm BlackoilPropsAdInterface& props, const RockCompressibility* rock_comp_props, WellsManager& wells_manager, - LinearSolverInterface& linsolver, + FullyImplicitSystemSolverInterface& linsolver, const double* gravity); SimulatorReport run(SimulatorTimer& timer, @@ -110,7 +110,7 @@ namespace Opm BlackoilPropsAdInterface& props, const RockCompressibility* rock_comp_props, WellsManager& wells_manager, - LinearSolverInterface& linsolver, + FullyImplicitSystemSolverInterface& linsolver, const double* gravity) { @@ -257,7 +257,7 @@ namespace Opm BlackoilPropsAdInterface& props, const RockCompressibility* rock_comp_props, WellsManager& wells_manager, - LinearSolverInterface& linsolver, + FullyImplicitSystemSolverInterface& linsolver, const double* gravity) : grid_(grid), props_(props), diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoil.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoil.hpp index 07e80fda2..0ce3aa8cf 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 FullyImplicitSystemSolverInterface; class SimulatorTimer; class BlackoilState; class WellStateFullyImplicitBlackoil; @@ -70,7 +70,7 @@ namespace Opm BlackoilPropsAdInterface& props, const RockCompressibility* rock_comp_props, WellsManager& wells_manager, - LinearSolverInterface& linsolver, + FullyImplicitSystemSolverInterface& linsolver, const double* gravity); /// Run the simulation.