multistagetimestepper.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:
//
// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
// SPDX-License-Identifier: GPL-3.0-or-later
//
#ifndef DUMUX_TIMESTEPPING_MULTISTAGE_TIMESTEPPER_HH
#define DUMUX_TIMESTEPPING_MULTISTAGE_TIMESTEPPER_HH
 
#include <memory>
#include <vector>
#include <cmath>
#include <iostream>
 
#include <dumux/io/format.hh>
 
#include <dumux/common/variablesbackend.hh>
#include <dumux/common/timeloop.hh>
#include <dumux/experimental/timestepping/multistagemethods.hh>
 
namespace Dumux::Experimental {
 
template<class Scalar>
class MultiStageParams
{
    struct Params {
        Scalar alpha, betaDt, timeAtStage, dtFraction;
        bool skipTemporal, skipSpatial;
    };
public:
    MultiStageParams(const MultiStageMethod<Scalar>& m, std::size_t i, const Scalar t, const Scalar dt)
    : size_(i+1)
    {
        params_.resize(size_);
        for (std::size_t k = 0; k < size_; ++k)
        {
            auto& p = params_[k];
            p.alpha = m.temporalWeight(i, k);
            p.betaDt = m.spatialWeight(i, k)*dt;
            p.timeAtStage = t + m.timeStepWeight(k)*dt;
            p.dtFraction = m.timeStepWeight(k);
 
            using std::abs;
            p.skipTemporal = (abs(p.alpha) < 1e-6);
            p.skipSpatial = (abs(p.betaDt) < 1e-6);
        }
    }
 
    std::size_t size () const
    { return size_; }
 
    Scalar temporalWeight (std::size_t k) const
    { return params_[k].alpha; }
 
    Scalar spatialWeight (std::size_t k) const
    { return params_[k].betaDt; }
 
    Scalar timeAtStage (std::size_t k) const
    { return params_[k].timeAtStage; }
 
    Scalar timeStepFraction (std::size_t k) const
    { return params_[k].dtFraction; }
 
    Scalar skipTemporal (std::size_t k) const
    { return params_[k].skipTemporal; }
 
    Scalar skipSpatial (std::size_t k) const
    { return params_[k].skipSpatial; }
 
private:
    std::size_t size_;
    std::vector<Params> params_;
};
 
template<class PDESolver, class Scalar = double>
class MultiStageTimeStepper
{
    using Variables = typename PDESolver::Variables;
    using StageParams = MultiStageParams<Scalar>;
    using Backend = VariablesBackend<Variables>;
 
public:
 
    MultiStageTimeStepper(std::shared_ptr<PDESolver> pdeSolver,
                          std::shared_ptr<const MultiStageMethod<Scalar>> msMethod,
                          const std::string& paramGroup = "")
    : pdeSolver_(pdeSolver)
    , msMethod_(msMethod)
    {
        initParams_(paramGroup);
 
        if (pdeSolver_->verbosity() >= 1)
            std::cout << "Initialize time stepper with method " << msMethod_->name()
                    << Fmt::format(" ({} stage{})", msMethod_->numStages(), (msMethod_->numStages() > 1 ? "s" : ""))
                    << std::endl;
    }
 
    void step(Variables& vars, const Scalar t, const Scalar dt)
    {
        // make sure there are no traces of previous stages
        pdeSolver_->assembler().clearStages();
 
        // do time integration
        bool converged = step_(vars, t, dt);
 
        // clear traces of previously registered stages
        pdeSolver_->assembler().clearStages();
 
        // if the solver didn't converge we can't recover
        if (!converged)
            DUNE_THROW(NumericalProblem, "Solver did not converge!");
    }
 
    void step(Variables& vars, TimeLoopBase<Scalar>& timeLoop)
    {
        // make sure there are no traces of previous stages
        pdeSolver_->assembler().clearStages();
 
        // try solving the non-linear system
        for (std::size_t i = 0; i <= maxTimeStepDivisions_; ++i)
        {
            // try solving the non-linear system
            bool converged = step_(vars, timeLoop.time(), timeLoop.timeStepSize());
            if (converged)
            {
                // clear traces of previously registered stages
                pdeSolver_->assembler().clearStages();
                return;
            }
 
            else if (!converged && i < maxTimeStepDivisions_)
            {
                // set solution to previous solution & reset time step
                Backend::update(vars, pdeSolver_->assembler().prevSol());
                pdeSolver_->assembler().resetTimeStep(Backend::dofs(vars));
 
                if (pdeSolver_->verbosity() >= 1)
                {
                    const auto dt = timeLoop.timeStepSize();
                    std::cout << Fmt::format("Solver did not converge with dt = {} seconds. ", dt)
                            << Fmt::format("Retrying with time step of dt = {} seconds.\n", dt*retryTimeStepReductionFactor_);
                }
 
                // try again with dt = dt * retryTimeStepReductionFactor_
                timeLoop.setTimeStepSize(timeLoop.timeStepSize() * retryTimeStepReductionFactor_);
            }
 
            else
            {
                pdeSolver_->assembler().clearStages();
                DUNE_THROW(NumericalProblem,
                    Fmt::format("Solver didn't converge after {} time-step divisions; dt = {}.\n",
                                maxTimeStepDivisions_, timeLoop.timeStepSize()));
            }
        }
 
        DUNE_THROW(Dune::InvalidStateException, "Unreachable");
    }
 
private:
    bool step_(Variables& vars, const Scalar t, const Scalar dt)
    {
        for (auto stageIdx = 1UL; stageIdx <= msMethod_->numStages(); ++stageIdx)
        {
            // prepare the assembler for this stage
            pdeSolver_->assembler().prepareStage(
                vars,
                std::make_shared<StageParams>(*msMethod_, stageIdx, t, dt)
            );
 
            // assemble & solve
            bool converged = pdeSolver_->apply(vars);
            if (!converged)
                return false;
        }
 
        return true;
    }
 
    void initParams_(const std::string& group = "")
    {
        maxTimeStepDivisions_ = getParamFromGroup<std::size_t>(group, "TimeStepper.MaxTimeStepDivisions", 10);
        retryTimeStepReductionFactor_ = getParamFromGroup<Scalar>(group, "TimeStepper.RetryTimeStepReductionFactor", 0.5);
    }
 
    std::shared_ptr<PDESolver> pdeSolver_;
    std::shared_ptr<const MultiStageMethod<Scalar>> msMethod_;
 
    Scalar maxTimeStepDivisions_;
    Scalar retryTimeStepReductionFactor_;
};
 
} // end namespace Dumux::Experimental
 
#endif