porousmediumflow/richards/newtonsolver.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:
/*****************************************************************************
 *   See the file COPYING for full copying permissions.                      *
 *                                                                           *
 *   This program 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.                                     *
 *                                                                           *
 *   This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.   *
 *****************************************************************************/
#ifndef DUMUX_RICHARDS_NEWTON_SOLVER_HH
#define DUMUX_RICHARDS_NEWTON_SOLVER_HH
 
#include <algorithm>
 
#include <dumux/common/properties.hh>
#include <dumux/nonlinear/newtonsolver.hh>
#include <dumux/discretization/elementsolution.hh>
 
namespace Dumux {
template <class Assembler, class LinearSolver>
class RichardsNewtonSolver : public NewtonSolver<Assembler, LinearSolver>
{
    using Scalar = typename Assembler::Scalar;
    using ParentType = NewtonSolver<Assembler, LinearSolver>;
    using Indices = typename Assembler::GridVariables::VolumeVariables::Indices;
    enum { pressureIdx = Indices::pressureIdx };
 
    using typename ParentType::Backend;
    using typename ParentType::SolutionVector;
 
public:
    using ParentType::ParentType;
    using typename ParentType::Variables;
 
private:
 
    void choppedUpdate_(Variables &varsCurrentIter,
                        const SolutionVector &uLastIter,
                        const SolutionVector &deltaU) final
    {
        auto uCurrentIter = uLastIter;
        uCurrentIter -= deltaU;
 
        // do not clamp anything after 5 iterations
        if (this->numSteps_ <= 4)
        {
            // clamp saturation change to at most 20% per iteration
            const auto& gridGeometry = this->assembler().gridGeometry();
            auto fvGeometry = localView(gridGeometry);
            for (const auto& element : elements(gridGeometry.gridView()))
            {
                fvGeometry.bindElement(element);
 
                for (auto&& scv : scvs(fvGeometry))
                {
                    auto dofIdxGlobal = scv.dofIndex();
 
                    // calculate the old wetting phase saturation
                    const auto& spatialParams = this->assembler().problem().spatialParams();
                    const auto elemSol = elementSolution(element, uCurrentIter, gridGeometry);
 
                    const auto fluidMatrixInteraction = spatialParams.fluidMatrixInteraction(element, scv, elemSol);
                    const Scalar pcMin = fluidMatrixInteraction.pc(1.0);
                    const Scalar pw = uLastIter[dofIdxGlobal][pressureIdx];
                    using std::max;
                    const Scalar pn = max(this->assembler().problem().nonwettingReferencePressure(), pw + pcMin);
                    const Scalar pcOld = pn - pw;
                    const Scalar SwOld = max(0.0, fluidMatrixInteraction.sw(pcOld));
 
                    // convert into minimum and maximum wetting phase pressures
                    const Scalar pwMin = pn - fluidMatrixInteraction.pc(SwOld - 0.2);
                    const Scalar pwMax = pn - fluidMatrixInteraction.pc(SwOld + 0.2);
 
                    // clamp the result
                    using std::clamp;
                    uCurrentIter[dofIdxGlobal][pressureIdx] = clamp(uCurrentIter[dofIdxGlobal][pressureIdx], pwMin, pwMax);
                }
            }
        }
 
        // update the variables
        this->solutionChanged_(varsCurrentIter, uCurrentIter);
 
        if (this->enableResidualCriterion())
            this->computeResidualReduction_(varsCurrentIter);
    }
};
 
} // end namespace Dumux
 
#endif