impet.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_IMPET_HH
#define DUMUX_IMPET_HH

 
#include "impetproperties.hh"
 
namespace Dumux
{
template<class TypeTag> class IMPET
{
    using Problem = GetPropType<TypeTag, Properties::Problem>;
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
 
    using SolutionTypes = GetProp<TypeTag, Properties::SolutionTypes>;
    using ElementMapper = typename SolutionTypes::ElementMapper;
    using TransportSolutionType = GetPropType<TypeTag, Properties::TransportSolutionType>;
 
    enum IterationType
        {
            noIter,
            iterToNumIter,
            iterToConverged,
        };
 
public:
    using SolutionType = typename SolutionTypes::ScalarSolution;

    void initialize()
    {
        //initial values of transported quantity
        problem_.transportModel().initialize();
        //call function with true to get a first initialisation of the pressure field
        problem_.pressureModel().initialize();
 
        //update constitutive functions
        problem_.pressureModel().updateMaterialLaws();
 
        Dune::dinfo << "IMPET: initialization done." << std::endl;
 
        return;
    }


    void update(const Scalar t, Scalar& dt, TransportSolutionType& updateVec)
    {
        if (iterFlag_ == noIter)
        {
            //update pressure
            problem_.pressureModel().update();
 
            //calculate defect of transported quantity
            problem_.transportModel().update(t, dt, updateVec, true);
 
            dt *= cFLFactor_;
        }
        else if (iterFlag_ == iterToNumIter || iterFlag_ == iterToConverged)
        {
            bool converg = false;
            int iter = 0;
            int iterTot = 0;
 
            // the method is valid for any transported quantity.
            TransportSolutionType transValueOldIter;
            problem_.transportModel().getTransportedQuantity(transValueOldIter);
            TransportSolutionType updateOldIter(transValueOldIter);
            updateOldIter = 0;
            TransportSolutionType transportedQuantity(transValueOldIter);
            TransportSolutionType updateHelp(transValueOldIter);
            TransportSolutionType updateDiff(transValueOldIter);
 
            while (!converg)
            {
                iter++;
                iterTot++;
 
                problem_.pressureModel().update();
 
                //calculate defect of transported quantity
                problem_.transportModel().update(t, dt, updateVec, true);
 
                updateHelp = updateVec;
                problem_.transportModel().getTransportedQuantity(transportedQuantity);
                transportedQuantity += (updateHelp *= (dt * cFLFactor_));
                transportedQuantity *= omega_;
                transValueOldIter *= (1 - omega_);
                transportedQuantity += transValueOldIter;
                updateDiff = updateVec;
                updateDiff -= updateOldIter;
                transValueOldIter = transportedQuantity;
                updateOldIter = updateVec;
                Dune::dinfo << " defect = " << dt * updateDiff.two_norm() / transportedQuantity.two_norm();
                // break criteria for iteration loop
                if (iterFlag_ == iterToConverged && dt * updateDiff.two_norm() / transportedQuantity.two_norm() <= maxDefect_)
                {
                    converg = true;
                }
                else if (iterFlag_ == iterToNumIter && iter > nIter_)
                {
                    converg = true;
                }
 
                if (iterFlag_ == iterToConverged && transportedQuantity.infinity_norm() > (1 + maxDefect_))
                {
                    converg = false;
                }
                if (!converg && iter > nIter_)
                {
                    converg = true;
                    std::cout << "Nonlinear loop in IMPET.update exceeded nIter = " << nIter_ << " iterations."
                              << std::endl;
                    std::cout << transportedQuantity.infinity_norm() << std::endl;
                }
            }
            // outputs
            if (iterFlag_ == iterToConverged)
                std::cout << "Iteration steps: " << iterTot << std::endl;
            std::cout.setf(std::ios::scientific, std::ios::floatfield);
 
            dt *= cFLFactor_;
        }
        else
        {
            DUNE_THROW(Dune::NotImplemented,"IMPET: Iterationtype not implemented!");
        }
 
    }
 
    void updateTransport(const Scalar t, Scalar& dt, TransportSolutionType& updateVec)
    {
        problem_.pressureModel().updateVelocity();
 
        problem_.transportModel().update(t, dt, updateVec, true);
    }


    template<class MultiWriter>
    void addOutputVtkFields(MultiWriter &writer)
    {
        problem_.pressureModel().addOutputVtkFields(writer);
        problem_.transportModel().addOutputVtkFields(writer);
 
        return;
    }

    const ElementMapper &dofMapper() const
    {
        return problem_.elementMapper();
    }


    IMPET(Problem& problem) :
        problem_(problem)
    {
        cFLFactor_ = getParam<Scalar>("Impet.CFLFactor");
        iterFlag_ = getParam<int>("Impet.IterationFlag", 0);
        nIter_ = getParam<int>("Impet.IterationNumber", 2);
        maxDefect_ = getParam<Scalar>("Impet.MaximumDefect", 1e-5);
        omega_ = getParam<Scalar>("Impet.RelaxationFactor", 1.0);
    }
 
private:
    Problem& problem_; 
    Scalar cFLFactor_;
    int iterFlag_; 
    int nIter_; 
    Scalar maxDefect_; 
    Scalar omega_; 
};
}
#endif