changed impes base clases to impet

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@4155 2fb0f335-1f38-0410-981e-8018bf24f1b0

dumux/decoupled/common/impet.hh

+// $Id$
+/*****************************************************************************
+ *   Copyright (C) 2007-2009 by Bernd Flemisch                               *
+ *   Copyright (C) 2008-2009 by Markus Wolff                                 *
+ *   Institute of Hydraulic Engineering                                      *
+ *   University of Stuttgart, Germany                                        *
+ *   email: <givenname>.<name>@iws.uni-stuttgart.de                          *
+ *                                                                           *
+ *   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 2 of the License, or       *
+ *   (at your option) any later version, as long as this copyright notice    *
+ *   is included in its original form.                                       *
+ *                                                                           *
+ *   This program is distributed WITHOUT ANY WARRANTY.                       *
+ *****************************************************************************/
+#ifndef DUMUX_IMPET_HH
+#define DUMUX_IMPET_HH
+
+/**
+ * @file
+ * @brief  IMPET scheme
+ * @author Bernd Flemisch, Markus Wolff
+ */
+
+#include "impetproperties.hh"
+
+namespace Dumux
+{
+/**
+ * \ingroup impet
+ * \brief IMplicit Pressure Explicit Transport (IMPET) scheme for the solution of weakly coupled diffusion/transport problems.
+ *
+ * The model implements the decoupled equations of two-phase flow.
+ * These equations can be derived from the two-phase flow equations shown for the two-phase box model (TwoPBoxModel).
+ * The first equation to solve is a pressure equation of elliptic character. The second one is a transport equation (e.g. for saturation, concentration,...),
+ * which can be hyperbolic or parabolic.
+ *
+ * As the equations are only weakly coupled they do not have to be solved simultaneously
+ * but can be solved sequentially. First the pressure equation is solved implicitly,
+ * second the transport equation can be solved explicitly. This solution procedure is called IMPES algorithm
+ * (IMplicit Pressure Explicit Saturation) for immiscible flow or IMPEC algorithm
+ * (IMplicit Pressure Explicit Concentration) for miscible flow.
+ */
+template<class TypeTag> class IMPET
+{
+    typedef IMPET<TypeTag> ThisType;
+    typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+
+    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+
+
+    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+
+    typedef typename GET_PROP(TypeTag, PTAG(SolutionTypes)) SolutionTypes;
+    typedef typename GET_PROP_TYPE(TypeTag, PTAG(TransportSolutionType)) TransportSolutionType;
+    enum
+    {
+        dim = GridView::dimension, dimWorld = GridView::dimensionworld
+    };
+
+    typedef typename SolutionTypes::ScalarSolution ScalarSolutionType;
+
+    typedef Dune::FieldVector<Scalar, dim> LocalPosition;
+    typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+public:
+    typedef typename SolutionTypes::ScalarSolution SolutionType;
+
+    //! Set initial solution and initialize parameters
+    virtual 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();
+        problem.pressureModel().calculateVelocity();
+
+        //update constitutive functions
+        problem.pressureModel().updateMaterialLaws();
+
+        Dune::dinfo << "IMPET: initialization done." << std::endl;
+
+        return;
+    }
+
+    //! Calculate the update.
+    /*!
+     *  \param  t         time
+     *  \param dt         time step size
+     *  \param updateVec  vector for the update values
+     *
+     *  Calculates the new pressure and velocity and determines the time step size and the update of the transported quantity for the explicit time step
+     */
+    void update(const Scalar t, Scalar& dt, TransportSolutionType& updateVec)
+    {
+    	// the method is valid for any transported quantity.
+        int transSize = problem.variables().gridSize();
+        TransportSolutionType transportedQuantity(problem.variables().transportedQuantity());
+        TransportSolutionType transValueOldIter(problem.variables().transportedQuantity());
+        TransportSolutionType transValueHelp(transSize);
+        TransportSolutionType transValueDiff(transSize);
+        TransportSolutionType updateOldIter(transSize);
+        TransportSolutionType updateHelp(transSize);
+        TransportSolutionType updateDiff(transSize);
+
+        bool converg = false;
+        int iter = 0;
+        int iterTot = 0;
+        updateOldIter = 0;
+        while (!converg)
+        {
+            iter++;
+            iterTot++;
+
+            problem.pressureModel().pressure(false);
+
+            //calculate velocities
+            problem.pressureModel().calculateVelocity();
+
+            //calculate defect of transported quantity
+            problem.transportModel().update(t, dt, updateVec, true);
+
+            if (iterFlag_)
+            { // only needed if iteration has to be done
+                updateHelp = updateVec;
+                transportedQuantity = problem.variables().transportedQuantity();
+                transportedQuantity += (updateHelp *= (dt * cFLFactor_));
+                transportedQuantity *= omega_;
+                transValueHelp = transValueOldIter;
+                transValueHelp *= (1 - omega_);
+                transportedQuantity += transValueHelp;
+                updateDiff = updateVec;
+                updateDiff -= updateOldIter;
+                transValueOldIter = transportedQuantity;
+                updateOldIter = updateVec;
+                //                problem.transportModel().updateMaterialLaws(transportedQuantity, true);
+            }
+            // break criteria for iteration loop
+            if (iterFlag_ == 2 && dt * updateDiff.two_norm() / transportedQuantity.two_norm() <= maxDefect_)
+            {
+                converg = true;
+            }
+            else if (iterFlag_ == 1 && iter > nIter_)
+            {
+                converg = true;
+            }
+            else if (iterFlag_ == 0)
+            {
+                converg = true;
+            }
+            if (iterFlag_ == 2 && transportedQuantity.infinity_norm() > (1 + maxDefect_))
+            {
+                converg = false;
+            }
+            if (!converg && iter > nIter_)
+            {
+                std::cout << "Nonlinear loop in IMPET.update exceeded nIter = " << nIter_ << " iterations." << std::endl;
+                std::cout << transportedQuantity.infinity_norm() << std::endl;
+            }
+        }
+        // outputs
+        if (iterFlag_ == 2)
+            std::cout << "Iteration steps: " << iterTot << std::endl;
+        std::cout.setf(std::ios::scientific, std::ios::floatfield);
+
+        dt *= cFLFactor_;
+    }
+
+    //! \brief Write data files
+    /*!
+     *  \param name file name
+     *  \param k format parameter
+     */
+    template<class MultiWriter>
+    void addOutputVtkFields(MultiWriter &writer)
+    {
+        problem.pressureModel().addOutputVtkFields(writer);
+        problem.transportModel().addOutputVtkFields(writer);
+
+        return;
+    }
+
+    // serialization methods
+    template<class Restarter>
+    void serialize(Restarter &res)
+    {
+        problem.variables().serialize<Restarter> (res);
+    }
+    template<class Restarter>
+    void deserialize(Restarter &res)
+    {
+        problem.variables().deserialize<Restarter> (res);
+    }
+
+    //! Constructs an IMPET object
+    /**
+     * \param prob Problem
+     */
+    IMPET(Problem& prob) :
+        problem(prob)
+    {
+        cFLFactor_ = GET_PROP_VALUE(TypeTag, PTAG(CFLFactor));
+        iterFlag_ = GET_PROP_VALUE(TypeTag, PTAG(IterationFlag));
+        nIter_ = GET_PROP_VALUE(TypeTag, PTAG(IterationNumber));
+        maxDefect_ = GET_PROP_VALUE(TypeTag, PTAG(MaximumDefect));
+        omega_ = GET_PROP_VALUE(TypeTag, PTAG(RelaxationFactor));
+    }
+
+protected:
+
+    Problem& problem; //!< object of type Problem including problem
+private:
+    Scalar cFLFactor_;
+    int iterFlag_;
+    int nIter_;
+    Scalar maxDefect_;
+    Scalar omega_;
+};
+}
+#endif

dumux/decoupled/common/impetproblem.hh

+// $Id$
+/*****************************************************************************
+ *   Copyright (C) 2010 by Markus Wolff, Andreas Lauser                      *
+ *   Institute of Hydraulic Engineering                                      *
+ *   University of Stuttgart, Germany                                        *
+ *   email: <givenname>.<name>@iws.uni-stuttgart.de                          *
+ *                                                                           *
+ *   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 2 of the License, or       *
+ *   (at your option) any later version, as long as this copyright notice    *
+ *   is included in its original form.                                       *
+ *                                                                           *
+ *   This program is distributed WITHOUT ANY WARRANTY.                       *
+ *****************************************************************************/
+#ifndef DUMUX_IMPETPROBLEM_HH
+#define DUMUX_IMPETPROBLEM_HH
+
+#include "impetproperties.hh"
+#include <dumux/io/vtkmultiwriter.hh>
+#include <dumux/io/restart.hh>
+
+#include <dumux/common/timemanager.hh>
+
+/**
+ * @file
+ * @brief  Base class for defining an instance of the diffusion problem
+ * @author Bernd Flemisch
+ */
+
+namespace Dumux
+{
+
+/*! \ingroup impet
+ * @brief base class for problems using a sequential implicit-explicit strategy
+ *
+ *  Template parameters are:
+ *
+ *  - TypeTag      problem TypeTag
+ *  - Implementation        problem implementation
+ */
+template<class TypeTag, class Implementation>
+class IMPETProblem
+{
+private:
+    typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+    typedef Dumux::TimeManager<TypeTag>  TimeManager;
+
+    typedef Dumux::VtkMultiWriter<GridView>  VtkMultiWriter;
+
+    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Variables)) Variables;
+
+    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) IMPETModel;
+    typedef typename GET_PROP_TYPE(TypeTag, PTAG(TransportSolutionType)) TransportSolutionType;
+    typedef typename GET_PROP_TYPE(TypeTag, PTAG(PressureModel)) PressureModel;
+    typedef typename GET_PROP_TYPE(TypeTag, PTAG(TransportModel)) TransportModel;
+
+    typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+
+    enum
+    {
+        dim = GridView::dimension,
+        dimWorld = GridView::dimensionworld,
+    };
+    enum
+    {
+        wetting = 0, nonwetting = 1
+    };
+
+    typedef Dune::FieldVector<Scalar,dim> LocalPosition;
+    typedef Dune::FieldVector<Scalar,dimWorld> GlobalPosition;
+    typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
+
+public:
+
+    //! Constructs an object of type IMPETProblemProblem
+    /** @param gridView gridview to the grid.
+     *  @param verbose verbosity.
+     */
+    IMPETProblem(const GridView &gridView, bool verbose = true)
+        : gridView_(gridView),
+          bboxMin_(std::numeric_limits<double>::max()),
+          bboxMax_(-std::numeric_limits<double>::max()),
+          timeManager_(verbose),
+          variables_(gridView),
+          dt_(0),
+          resultWriter_(asImp_().name())
+    {
+        // calculate the bounding box of the grid view
+        VertexIterator vIt = gridView.template begin<dim>();
+        const VertexIterator vEndIt = gridView.template end<dim>();
+        for (; vIt!=vEndIt; ++vIt) {
+            for (int i=0; i<dim; i++) {
+                bboxMin_[i] = std::min(bboxMin_[i], vIt->geometry().corner(0)[i]);
+                bboxMax_[i] = std::max(bboxMax_[i], vIt->geometry().corner(0)[i]);
+            }
+        }
+
+        pressModel_ = new PressureModel(asImp_());
+        satModel_ = new TransportModel(asImp_());
+        model_ = new IMPETModel(asImp_()) ;
+    }
+
+    //! destructor
+    virtual ~IMPETProblem ()
+    {
+        delete pressModel_;
+        delete satModel_;
+        delete model_;
+    }
+
+    /*!
+     * \name Simulation steering
+     */
+    // \{
+
+    /*!
+     * \brief Called by the Dumux::TimeManager in order to
+     *        initialize the problem.
+     */
+    void init()
+    {
+        // set the initial condition of the model
+        model().initialize();
+    }
+
+    /*!
+     * \brief Called by the time manager before the time integration.
+     */
+    void preTimeStep()
+    {}
+
+    /*!
+     * \brief Called by Dumux::TimeManager in order to do a time
+     *        integration on the model.
+     *
+     * \note \a timeStepSize and \a nextStepSize are references and may
+     *       be modified by the timeIntegration(). On exit of this
+     *       function \a timeStepSize must contain the step size
+     *       actually used by the time integration for the current
+     *       steo, and \a nextStepSize must contain a suggestion for the
+     *       next time step size.
+     */
+    void timeIntegration()
+    {
+        // allocate temporary vectors for the updates
+        typedef TransportSolutionType Solution;
+        Solution k1 = asImp_().variables().transportedQuantity();
+
+        Scalar t = timeManager().time();
+
+        // obtain the first update and the time step size
+        model().update(t, dt_, k1);
+
+        //make sure t_old + dt is not larger than tend
+        dt_ = std::min(dt_, timeManager().episodeMaxTimeStepSize());
+
+        // check if we are in first TS and an initialDt was assigned
+        if (t==0. && timeManager().timeStepSize()!=0.)
+        {
+            // check if assigned initialDt is in accordance with dt_ from first transport step
+            if (timeManager().timeStepSize() > dt_)
+                Dune::dwarn << "initial timestep of size " << timeManager().timeStepSize()
+                            << "is larger then dt_= "<<dt_<<" from transport" << std::endl;
+            // internally assign next tiestep size
+            dt_ = std::min(dt_, timeManager().timeStepSize());
+        }
+
+        //assign next tiestep size
+        timeManager().setTimeStepSize(dt_);
+
+        // explicit Euler: Sat <- Sat + dt*N(Sat)
+        asImp_().variables().transportedQuantity() += (k1 *= timeManager().timeStepSize());
+    }
+
+    /*!
+     * \brief Called by Dumux::TimeManager whenever a solution for a
+     *        timestep has been computed and the simulation time has
+     *        been updated.
+     *
+     * This is used to do some janitorial tasks like writing the
+     * current solution to disk.
+     */
+    void postTimeStep()
+    {
+        asImp_().pressureModel().updateMaterialLaws();
+    };
+
+    /*!
+     * \brief Returns the current time step size [seconds].
+     */
+    Scalar timeStepSize() const
+    { return timeManager_.timeStepSize(); }
+
+    /*!
+     * \brief Sets the current time step size [seconds].
+     */
+    void setTimeStepSize(Scalar dt)
+    { return timeManager_.setTimeStepSize(dt); }
+
+    /*!
+     * \brief Called by Dumux::TimeManager whenever a solution for a
+     *        timestep has been computed and the simulation time has
+     *        been updated.
+     */
+    Scalar nextTimeStepSize()
+    { return dt_;}
+
+    /*!
+     * \brief Returns true if a restart file should be written to
+     *        disk.
+     *
+     * The default behaviour is to write one restart file every 5 time
+     * steps. This file is intented to be overwritten by the
+     * implementation.
+     */
+    bool shouldWriteRestartFile() const
+    {
+        return
+            timeManager().timeStepIndex() > 0 &&
+            (timeManager().timeStepIndex() % 5 == 0);
+    }
+
+    /*!
+     * \brief Returns true if the current solution should be written to
+     *        disk (i.e. as a VTK file)
+     *
+     * The default behaviour is to write out every the solution for
+     * very time step. This file is intented to be overwritten by the
+     * implementation.
+     */
+    bool shouldWriteOutput() const
+    { return true; }
+
+    /*!
+     * \brief Called when the end of an simulation episode is reached.
+     */
+    void episodeEnd()
+    {
+        std::cerr << "The end of an episode is reached, but the problem "
+                  << "does not override the episodeEnd() method. "
+                  << "Doing nothing!\n";
+    };
+
+    // \}
+
+    /*!
+     * \brief The problem name.
+     *
+     * This is used as a prefix for files generated by the simulation.
+     * It could be either overwritten by the problem files, or simply
+     * declared over the setName() function in the application file.
+     */
+    const char *name() const
+    {
+        return simname_.c_str();
+    }
+
+    /*!
+     * \brief Set the problem name.
+     *
+     * This function sets the simulation name, which should be called before
+     * the application porblem is declared! If not, the default name "sim"
+     * will be used.
+     */
+    static void setName(const char *newName)
+    {
+        simname_ = newName;
+    }
+
+    /*!
+     * \brief The GridView which used by the problem.
+     */
+    const GridView &gridView() const
+    { return gridView_; }
+
+    /*!
+     * \brief The coordinate of the corner of the GridView's bounding
+     *        box with the smallest values.
+     */
+    const GlobalPosition &bboxMin() const
+    { return bboxMin_; }
+
+    /*!
+     * \brief The coordinate of the corner of the GridView's bounding
+     *        box with the largest values.
+     */
+    const GlobalPosition &bboxMax() const
+    { return bboxMax_; }
+
+    /*!
+     * \brief Returns TimeManager object used by the simulation
+     */
+    TimeManager &timeManager()
+    { return timeManager_; }
+
+    /*!
+     * \copydoc timeManager()
+     */
+    const TimeManager &timeManager() const
+    { return timeManager_; }
+
+    Variables& variables ()
+    { return variables_; }
+
+    const Variables& variables () const
+    { return variables_; }
+
+    /*!
+     * \brief Returns numerical model used for the problem.
+     */
+    IMPETModel &model()
+    { return *model_; }
+
+    /*!
+     * \copydoc model()
+     */
+    const IMPETModel &model() const
+    { return *model_; }
+    // \}
+
+    /*!
+     * \brief Returns numerical model used for the problem.
+     */
+    PressureModel &pressureModel()
+    { return *pressModel_; }
+
+    /*!
+     * \copydoc model()
+     */
+    const PressureModel &pressureModel() const
+    { return *pressModel_; }
+    // \}
+
+    /*!
+     * \brief Returns numerical model used for the problem.
+     */
+    TransportModel &transportModel()
+    { return *satModel_; }
+
+    /*!
+     * \copydoc model()
+     */
+    const TransportModel &transportModel() const
+    { return *satModel_; }
+    // \}
+
+    /*!
+     * \name Restart mechanism
+     */
+    // \{
+
+    /*!
+     * \brief This method writes the complete state of the problem
+     *        to the harddisk.
+     *
+     * The file will start with the prefix returned by the name()
+     * method, has the current time of the simulation clock in it's
+     * name and uses the extension <tt>.drs</tt>. (Dumux ReStart
+     * file.)  See Dumux::Restart for details.
+     */
+    void serialize()
+    {
+        typedef Dumux::Restart Restarter;
+
+        Restarter res;
+        res.serializeBegin(asImp_());
+        std::cerr << "Serialize to file " << res.fileName() << "\n";
+
+        timeManager_.serialize(res);
+        resultWriter_.serialize(res);
+        model().serialize(res);
+
+        res.serializeEnd();
+    }
+
+    /*!
+     * \brief This method restores the complete state of the problem
+     *        from disk.
+     *
+     * It is the inverse of the serialize() method.
+     */
+    void deserialize(double t)
+    {
+        typedef Dumux::Restart Restarter;
+
+        Restarter res;
+        res.deserializeBegin(asImp_(), t);
+        std::cerr << "Deserialize from file " << res.fileName() << "\n";
+
+        timeManager_.deserialize(res);
+        resultWriter_.deserialize(res);
+        model().deserialize(res);
+
+        res.deserializeEnd();
+    };