From 4111797822da79002cd3dfd53c737c0fadd260a3 Mon Sep 17 00:00:00 2001 From: Robert Kloefkorn Date: Mon, 8 Feb 2016 17:13:43 +0100 Subject: [PATCH] BlackoilModelBase: added parameter singlePrecision and print residual to large at right place. --- opm/autodiff/BlackoilModelBase_impl.hpp | 35 ++++++++++--------- opm/autodiff/LinearisedBlackoilResidual.hpp | 2 ++ .../NewtonIterationBlackoilInterleaved.cpp | 6 ++-- ...FullyImplicitCompressiblePolymerSolver.cpp | 17 ++++----- 4 files changed, 33 insertions(+), 27 deletions(-) diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index d290b807a..6eb4ece58 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -250,6 +250,8 @@ namespace detail { const bool converged = asImpl().getConvergence(dt, iteration); const bool must_solve = (iteration < nonlinear_solver.minIter()) || (!converged); if (must_solve) { + residual_.singlePrecision = false ; + // Compute the nonlinear update. V dx = asImpl().solveJacobianSystem(); @@ -2705,22 +2707,6 @@ namespace detail { // Residual in Pascal can have high values and still be ok. const double maxWellResidualAllowed = 1000.0 * maxResidualAllowed(); - for (int idx = 0; idx < nm; ++idx) { - if (std::isnan(mass_balance_residual[idx]) - || std::isnan(CNV[idx]) - || (idx < np && std::isnan(well_flux_residual[idx]))) { - OPM_THROW(Opm::NumericalProblem, "NaN residual for phase " << materialName(idx)); - } - if (mass_balance_residual[idx] > maxResidualAllowed() - || CNV[idx] > maxResidualAllowed() - || (idx < np && well_flux_residual[idx] > maxResidualAllowed())) { - OPM_THROW(Opm::NumericalProblem, "Too large residual for phase " << materialName(idx)); - } - } - if (std::isnan(residualWell) || residualWell > maxWellResidualAllowed) { - OPM_THROW(Opm::NumericalProblem, "NaN or too large residual for well control equation"); - } - if ( terminal_output_ ) { // Only rank 0 does print to std::cout @@ -2755,6 +2741,23 @@ namespace detail { std::cout.precision(oprec); std::cout.flags(oflags); } + + for (int idx = 0; idx < nm; ++idx) { + if (std::isnan(mass_balance_residual[idx]) + || std::isnan(CNV[idx]) + || (idx < np && std::isnan(well_flux_residual[idx]))) { + OPM_THROW(Opm::NumericalProblem, "NaN residual for phase " << materialName(idx)); + } + if (mass_balance_residual[idx] > maxResidualAllowed() + || CNV[idx] > maxResidualAllowed() + || (idx < np && well_flux_residual[idx] > maxResidualAllowed())) { + OPM_THROW(Opm::NumericalProblem, "Too large residual for phase " << materialName(idx)); + } + } + if (std::isnan(residualWell) || residualWell > maxWellResidualAllowed) { + OPM_THROW(Opm::NumericalProblem, "NaN or too large residual for well control equation"); + } + return converged; } diff --git a/opm/autodiff/LinearisedBlackoilResidual.hpp b/opm/autodiff/LinearisedBlackoilResidual.hpp index 91fd8478b..e740750fb 100644 --- a/opm/autodiff/LinearisedBlackoilResidual.hpp +++ b/opm/autodiff/LinearisedBlackoilResidual.hpp @@ -65,6 +65,8 @@ namespace Opm std::vector matbalscale; + bool singlePrecision ; + /// The size of the non-linear system. int sizeNonLinear() const; }; diff --git a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp index f77290cee..c816eb46c 100644 --- a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp +++ b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp @@ -214,7 +214,7 @@ namespace Opm template std::unique_ptr constructPrecond(Operator& opA, const Dune::Amg::SequentialInformation&) const { - const double relax = 1.0; + const double relax = 0.9; std::unique_ptr precond(new SeqPreconditioner(opA.getmat(), relax)); return precond; } @@ -227,7 +227,7 @@ namespace Opm constructPrecond(Operator& opA, const Comm& comm) const { typedef std::unique_ptr Pointer; - const double relax = 1.0; + const double relax = 0.9; return Pointer(new ParPreconditioner(opA.getmat(), comm, relax)); } #endif @@ -542,7 +542,7 @@ namespace Opm // get np and call appropriate template method const int np = residual.material_balance_eq.size(); #if ! HAVE_UMFPACK - const bool singlePrecision = false ; // residual.singlePrecision ; + const bool singlePrecision = residual.singlePrecision ; const NewtonIterationBlackoilInterface& newtonIncrement = singlePrecision ? detail::NewtonIncrement< maxNumberEquations_, float > :: get( newtonIncrementSinglePrecision_, parameters_, parallelInformation_, np ) : detail::NewtonIncrement< maxNumberEquations_, double > :: get( newtonIncrementDoublePrecision_, parameters_, parallelInformation_, np ); diff --git a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp index 84d4fdf0f..f82905eda 100644 --- a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp +++ b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp @@ -189,7 +189,8 @@ namespace { , residual_ ( { std::vector(fluid.numPhases() + 1, ADB::null()), ADB::null(), ADB::null(), - { 1.1169, 1.0031, 0.0031, 1.0 }} ) // default scaling + { 1.1169, 1.0031, 0.0031, 1.0 }, + false } ) // default scaling { } @@ -369,7 +370,7 @@ namespace { assert(not x.concentration().empty()); const V c = Eigen::Map(&x.concentration()[0], nc); state.concentration = ADB::constant(c); - + // Well rates. assert (not xw.wellRates().empty()); // Need to reshuffle well rates, from ordered by wells, then phase, @@ -400,7 +401,7 @@ namespace { const int np = x.numPhases(); std::vector vars0; - + // Initial pressure. assert (not x.pressure().empty()); const V p = Eigen::Map(& x.pressure()[0], nc, 1); @@ -496,12 +497,12 @@ namespace { const double dead_pore_vol = polymer_props_ad_.deadPoreVol(); rq_[2].accum[aix] = pv_mult * rq_[0].b * sat[0] * c * (1. - dead_pore_vol) + pv_mult * rho_rock * (1. - phi) / phi * ads; } - - void + + void FullyImplicitCompressiblePolymerSolver:: computeCmax(PolymerBlackoilState& state) { @@ -654,7 +655,7 @@ namespace { const V poly_mc_cell = subset(mc, well_cells).value(); const V poly_in_c = poly_in_perf;// * poly_mc_cell; const V poly_mc = producer.select(poly_mc_cell, poly_in_c); - + 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; @@ -960,7 +961,7 @@ namespace { - // here mc means m(c) * c. + // here mc means m(c) * c. ADB FullyImplicitCompressiblePolymerSolver::computeMc(const SolutionState& state) const { @@ -1028,7 +1029,7 @@ namespace { { const int nc = grid_.number_of_cells; const DataBlock s = Eigen::Map(& state.saturation()[0], nc, 2); - + const V so = s.col(1); for (V::Index c = 0, e = so.size(); c != e; ++c) { phaseCondition_[c].setFreeWater();