BlackoilModelBase: added parameter singlePrecision and print residual to large at right

place.
This commit is contained in:
Robert Kloefkorn 2016-02-08 17:13:43 +01:00
parent 915177979c
commit 4111797822
4 changed files with 33 additions and 27 deletions

View File

@ -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;
}

View File

@ -65,6 +65,8 @@ namespace Opm
std::vector<double> matbalscale;
bool singlePrecision ;
/// The size of the non-linear system.
int sizeNonLinear() const;
};

View File

@ -214,7 +214,7 @@ namespace Opm
template <class Operator>
std::unique_ptr<SeqPreconditioner> constructPrecond(Operator& opA, const Dune::Amg::SequentialInformation&) const
{
const double relax = 1.0;
const double relax = 0.9;
std::unique_ptr<SeqPreconditioner> 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<ParPreconditioner> 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 );

View File

@ -189,7 +189,8 @@ namespace {
, residual_ ( { std::vector<ADB>(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<const V>(&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<V> vars0;
// Initial pressure.
assert (not x.pressure().empty());
const V p = Eigen::Map<const V>(& 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<const DataBlock>(& 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();