fixpointcriterion.hh

Go to the documentation of this file.

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
  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 2 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 <http://www.gnu.org/licenses/>.

  Consult the COPYING file in the top-level source directory of this
  module for the precise wording of the license and the list of
  copyright holders.
*/
#ifndef EWOMS_ISTL_FIXPOINT_CRITERION_HH
#define EWOMS_ISTL_FIXPOINT_CRITERION_HH

#include "convergencecriterion.hh"

namespace Opm {
namespace Linear {

template <class Vector, class CollectiveCommunication>
class FixPointCriterion : public ConvergenceCriterion<Vector>
{
    using Scalar = typename Vector::field_type;
    using BlockType = typename Vector::block_type;

public:
    FixPointCriterion(const CollectiveCommunication& comm) : comm_(comm)
    {}

    FixPointCriterion(const CollectiveCommunication& comm,
                      const Vector& weightVec, Scalar reduction)
        : comm_(comm), weightVec_(weightVec), tolerance_(reduction)
    {}

    void setWeight(const Vector& weightVec)
    { weightVec_ = weightVec; }

    Scalar weight(int outerIdx, int innerIdx) const
    { return (weightVec_.size() == 0) ? 1.0 : weightVec_[outerIdx][innerIdx]; }

    void setTolerance(Scalar tol)
    { tolerance_ = tol; }

    Scalar tolerance() const
    { return tolerance_; }

    void setInitial(const Vector& curSol, const Vector&) override
    {
        lastSol_ = curSol;
        delta_ = 1000 * tolerance_;
    }

    void update(const Vector& curSol,
                const Vector&,
                const Vector&) override
    {
        assert(curSol.size() == lastSol_.size());

        delta_ = 0.0;
        for (size_t i = 0; i < curSol.size(); ++i) {
            for (size_t j = 0; j < BlockType::dimension; ++j) {
                delta_ =
                    std::max(delta_, weight(i, j)*std::abs(curSol[i][j] - lastSol_[i][j]));
            }
        }

        delta_ = comm_.max(delta_);
        lastSol_ = curSol;
    }

    bool converged() const override
    { return accuracy() < tolerance(); }

    Scalar accuracy() const
    { return delta_; }

private:
    const CollectiveCommunication& comm_;

    Vector lastSol_;   // solution of the last iteration
    Vector weightVec_; // solution of the last iteration
    Scalar delta_;     // the maximum of the absolute weighted difference of the
                       // last two iterations
    Scalar tolerance_; // the maximum allowed delta for the solution to be
                       // considered converged
};


}} // end namespace Linear, Opm

#endif