[next] introduce time-dependent problem

dumux/common/basicproperties.hh

@@ -97,7 +97,11 @@ NEW_PROP_TAG(PointSourceHelper);
  *
  * The default is to not limit the step size.
  */
-NEW_PROP_TAG(TimeManagerMaxTimeStepSize);
+NEW_PROP_TAG(TimeLoopMaxTimeStepSize);
+
+//! the maximum allowed number of timestep divisions for the
+//! Newton solver
+NEW_PROP_TAG(TimeLoopMaxTimeStepDivisions);
 
 //! Property to define the output level
 NEW_PROP_TAG(VtkOutputLevel);
@@ -117,7 +121,10 @@ SET_TYPE_PROP(NumericModel, NumEqVector, Dune::FieldVector<typename GET_PROP_TYP
 SET_TYPE_PROP(NumericModel, PrimaryVariables, typename GET_PROP_TYPE(TypeTag, NumEqVector));
 
 //! use an unlimited time step size by default
-SET_SCALAR_PROP(NumericModel, TimeManagerMaxTimeStepSize, std::numeric_limits<typename GET_PROP_TYPE(TypeTag,Scalar)>::max());
+SET_SCALAR_PROP(NumericModel, TimeLoopMaxTimeStepSize, std::numeric_limits<typename GET_PROP_TYPE(TypeTag,Scalar)>::max());
+
+//! set number of maximum timestep divisions to 10
+SET_INT_PROP(NumericModel, TimeLoopMaxTimeStepDivisions, 10);
 
 //! Set the ParameterTree property
 SET_PROP(NumericModel, ParameterTree)

dumux/common/timeloop.hh

+// -*- 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 2 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/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Manages the handling of time dependent problems
+ */
+#ifndef DUMUX_TIME_LOOP_HH
+#define DUMUX_TIME_LOOP_HH
+
+#include <algorithm>
+
+#include <dune/common/float_cmp.hh>
+#include <dune/common/timer.hh>
+#include <dune/common/parallel/mpihelper.hh>
+
+#include "propertysystem.hh"
+#include "parameters.hh"
+
+namespace Dumux
+{
+/*!
+ * \ingroup SimControl
+ * \brief Manages the handling of time dependent problems.
+ *
+ * This class facilitates the time management of the simulation.
+ * It doesn't manage any user data, but keeps track of what the
+ * current time, time step size and "episode" of the
+ * simulation is. It triggers the initialization of the problem and
+ * is responsible for the time control of a simulation run.
+ *
+ * The time manager allows to specify a sequence of "episodes" which
+ * determine the boundary conditions of a problem. This approach
+ * is handy if the problem is not static, i.e. the boundary
+ * conditions change over time.
+ *
+ * An episode is a span of simulated time in which
+ * the problem behaves in a specific way. It is characterized by
+ * the (simulation) time it starts, its length and a consecutive
+ * index starting at 0.
+ */
+template <class Scalar>
+class TimeLoop
+{
+public:
+    TimeLoop(Scalar startTime, Scalar dt, Scalar tEnd, bool verbose = true) : timer_(false)
+    {
+        verbose_ =
+            verbose &&
+            Dune::MPIHelper::getCollectiveCommunication().rank() == 0;
+
+        time_ = startTime;
+        endTime_ = tEnd;
+
+        timeStepSize_ = dt;
+        previousTimeStepSize_ = timeStepSize_;
+        maxTimeStepSize_ = std::numeric_limits<Scalar>::max();
+        timeStepIdx_ = 0;
+        finished_ = false;
+    }
+
+    /*!
+     *  \name Simulated time and time step management
+     * @{
+     */
+
+    /*!
+     * \brief Tells the time loop to start tracking the time.
+     */
+    void start()
+    {
+        timer_.start();
+        cpuTime_ = 0.0;
+    }
+
+    /*!
+     * \brief State info on cpu time.
+     */
+    void reportTimeStep()
+    {
+        auto timeStepCpuTime = timer_.elapsed();
+        cpuTime_ += timeStepCpuTime;
+
+        if (verbose_)
+        {
+                std::cout << "Time step " << timeStepIndex() << " done in "
+                          << timeStepCpuTime << " seconds. "
+                          << "Wall time: " << cpuTime_
+                          << ", time: " << time()
+                          << ", time step size: " << timeStepSize()
+                          << std::endl;
+        }
+
+        timer_.reset();
+    }
+
+    /*!
+     * \brief Advance time step.
+     */
+    void advanceTimeStep()
+    {
+        timeStepIdx_++;
+        time_ += timeStepSize_;
+    }
+
+    /*!
+     * \brief Print final status and stops tracking the time.
+     */
+    void finalize()
+    {
+        timer_.stop();
+        cpuTime_ += timer_.elapsed();
+
+        if (verbose_)
+        {
+            const auto& comm = Dune::MPIHelper::getCollectiveCommunication();
+            std::cout << "Simulation took " << timer_.elapsed() << " seconds on "
+                      << comm.size() << " processes.\n"
+                      << "The cumulative CPU time was " << comm.sum(cpuTime_) << " seconds.\n";
+        }
+    }
+
+    /*!
+     * \brief Set the current simulated time, don't change the current
+     *        time step index.
+     *
+     * \param t The time \f$\mathrm{[s]}\f$ which should be jumped to
+     */
+    void setTime(Scalar t)
+    { time_ = t; }
+
+    /*!
+     * \brief Set the current simulated time and the time step index.
+     *
+     * \param t The time \f$\mathrm{[s]}\f$ which should be jumped to
+     * \param stepIdx The new time step index
+     */
+    void setTime(Scalar t, int stepIdx)
+    { time_ = t; timeStepIdx_ = stepIdx; }
+
+    /*!
+     * \brief Return the time \f$\mathrm{[s]}\f$ before the time integration.
+     * To get the time after the time integration you have to add timeStepSize() to
+     * time().
+     */
+    Scalar time() const
+    { return time_; }
+
+    /*!
+     * \brief Returns the number of (simulated) seconds which the simulation runs.
+     */
+    Scalar endTime() const
+    { return endTime_; }
+
+    /*!
+     * \brief Set the time of simulated seconds at which the simulation runs.
+     *
+     * \param t The time \f$\mathrm{[s]}\f$ at which the simulation is finished
+     */
+    void setEndTime(Scalar t)
+    { endTime_ = t; }
+
+    /*!
+     * \brief Returns the current wall time (cpu time).
+     */
+    double wallTime() const
+    {  return cpuTime_; }
+
+    /*!
+     * \brief Set the current time step size to a given value.
+     *
+     * If the step size would exceed the length of the current
+     * episode, the timeStep() method will take care that the step
+     * size won't exceed the episode or the end of the simulation,
+     * though.
+     *
+     * \param dt The new value for the time step size \f$\mathrm{[s]}\f$
+     */
+    void setTimeStepSize(Scalar dt)
+    {
+        using std::min;
+        timeStepSize_ = min(dt, maxTimeStepSize());
+    }
+
+    /*!
+     * \brief Set the maximum time step size to a given value.
+     *
+     * \param dt The new value for the maximum time step size \f$\mathrm{[s]}\f$
+     */
+    void setMaxTimeStepSize(Scalar maxDt)
+    { maxTimeStepSize_ = maxDt; }
+
+    /*!
+     * \brief Returns the suggested time step length \f$\mathrm{[s]}\f$ so that we
+     *        don't miss the beginning of the next episode or cross
+     *        the end of the simulation.
+     */
+    Scalar timeStepSize() const
+    { return timeStepSize_; }
+
+    /*!
+     * \brief Returns the size of the previous time step \f$\mathrm{[s]}\f$.
+     */
+    Scalar previousTimeStepSize() const
+    { return previousTimeStepSize_; }
+
+    /*!
+     * \brief Returns number of time steps which have been
+     *        executed since the beginning of the simulation.
+     */
+    int timeStepIndex() const
+    { return timeStepIdx_; }
+
+    /*!
+     * \brief Specify whether the simulation is finished
+     *
+     * \param finished If true the simulation is considered finished
+     *                 before the end time is reached, else it is only
+     *                 considered finished if the end time is reached.
+     */
+    void setFinished(bool finished = true)
+    { finished_ = finished; }
+
+    /*!
+     * \brief Returns true if the simulation is finished.
+     *
+     * This is the case if either setFinished(true) has been called or
+     * if the end time is reached.
+     */
+    bool finished() const
+    { return finished_ || time() >= endTime(); }
+
+    /*!
+     * \brief Returns true if the simulation is finished after the
+     *        time level is incremented by the current time step size.
+     */
+    bool willBeFinished() const
+    { return finished_ || time() + timeStepSize() >= endTime(); }
+
+    /*!
+     * \brief Aligns dt to the episode boundary or the end time of the
+     *        simulation.
+     */
+    Scalar maxTimeStepSize() const
+    {
+        if (finished())
+            return 0.0;
+
+        using std::max;
+        using std::min;
+        // return min(min(episodeMaxTimeStepSize(),
+        //                problem_->maxTimeStepSize()),
+        //            max<Scalar>(0.0, endTime() - time()));
+        return min(maxTimeStepSize_,
+                   max<Scalar>(0.0, endTime() - time()));
+    }
+
+    /*
+     * @}
+     */
+
+private:
+    Dune::Timer timer_;
+    Scalar time_;
+    Scalar endTime_;
+    double cpuTime_;
+
+    Scalar timeStepSize_;
+    Scalar previousTimeStepSize_;
+    Scalar maxTimeStepSize_;
+    int timeStepIdx_;
+    bool finished_;
+    bool verbose_;
+};
+}
+
+#endif

dumux/implicit/cellcentered/assembler.hh

@@ -26,6 +26,7 @@
 
 #include <dune/istl/matrixindexset.hh>
 
+#include <dumux/common/timeloop.hh>
 #include <dumux/implicit/properties.hh>
 #include <dumux/discretization/methods.hh>
 
@@ -79,7 +80,7 @@ class CCImplicitAssembler
 public:
     using ResidualType = SolutionVector;
 
-    //! The constructor
+    //! The constructor for stationary problems
     CCImplicitAssembler(std::shared_ptr<const Problem> problem,
                         std::shared_ptr<const FVGridGeometry> fvGridGeometry,
                         std::shared_ptr<GridVariables> gridVariables)
@@ -88,12 +89,22 @@ public:
     , gridVariables_(gridVariables)
     {}
 
+    //! The constructor for instationary problems
+    CCImplicitAssembler(std::shared_ptr<const Problem> problem,
+                        std::shared_ptr<const FVGridGeometry> fvGridGeometry,
+                        std::shared_ptr<GridVariables> gridVariables,
+                        std::shared_ptr<TimeLoop<Scalar>> timeLoop)
+    : problem_(problem)
+    , fvGridGeometry_(fvGridGeometry)
+    , gridVariables_(gridVariables)
+    , localResidual_(timeLoop)
+    {}
+
     /*!
      * \brief Assembles the global Jacobian of the residual
      *        and the residual for the current solution.
      */
-    void assembleJacobianAndResidual(const SolutionVector& curSol,
-                                     const SolutionVector& prevSol)
+    void assembleJacobianAndResidual(const SolutionVector& curSol)
     {
         if(!matrix_)
             DUNE_THROW(Dune::InvalidStateException,
@@ -113,7 +124,7 @@ public:
         {
             // let the local assembler add the element contributions
             for (const auto element : elements(gridView()))
-                LocalAssembler::assemble(*this, residual(), element, curSol, prevSol);
+                LocalAssembler::assemble(*this, residual(), element, curSol);
 
             // if we get here, everything worked well
             succeeded = true;
@@ -137,8 +148,7 @@ public:
     /*!
      * \brief Assembles only the global Jacobian of the residual.
      */
-    void assembleJacobian(const SolutionVector& curSol,
-                          const SolutionVector& prevSol)
+    void assembleJacobian(const SolutionVector& curSol)
     {
         resetMatrix_();
 
@@ -148,7 +158,7 @@ public:
         {
             // let the local assembler add the element contributions
             for (const auto& element : elements(gridView()))
-                LocalAssembler::assemble(*this, element, curSol, prevSol);
+                LocalAssembler::assemble(*this, element, curSol);
 
             // if we get here, everything worked well
             succeeded = true;
@@ -170,11 +180,11 @@ public:
     }
 
     //! compute the residuals
-    void assembleResidual(const SolutionVector& curSol, const SolutionVector& prevSol) const
-    { assembleResidual(residual(), curSol, prevSol); }
+    void assembleResidual(const SolutionVector& curSol) const
+    { assembleResidual(residual(), curSol); }
 
     //! compute the residuals
-    void assembleResidual(ResidualType& r, const SolutionVector& curSol, const SolutionVector& prevSol) const
+    void assembleResidual(ResidualType& r, const SolutionVector& curSol) const
     {
         r = 0.0;
         for (const auto& element : elements(gridView()))
@@ -186,17 +196,16 @@ public:
                                                 element,
                                                 fvGridGeometry(),
                                                 gridVariables(),
-                                                curSol,
-                                                prevSol)[0];
+                                                curSol)[0];
             }
         }
     }
 
     //! computes the global residual
-    Scalar globalResidual(const SolutionVector& curSol, const SolutionVector& prevSol) const
+    Scalar globalResidual(const SolutionVector& curSol) const
     {
         ResidualType residual(numDofs());
-        assembleResidual(residual, curSol, prevSol);
+        assembleResidual(residual, curSol);
 
         // calculate the square norm of the residual
         Scalar result2 = residual.two_norm2();
@@ -212,8 +221,8 @@ public:
      *        This also resizes the containers to the required sizes and sets the
      *        sparsity pattern of the matrix.
      */
-    void setLinearSystem(std::shared_ptr<JacobianMatrix>& A,
-                         std::shared_ptr<SolutionVector>& r)
+    void setLinearSystem(std::shared_ptr<JacobianMatrix> A,
+                         std::shared_ptr<SolutionVector> r)
     {
         matrix_ = A;
         residual_ = r;
@@ -241,6 +250,19 @@ public:
         allocateLinearSystem();
     }
 
+    /*!
+     * \brief Sets the solution from which to start the time integration. Has to be
+     *        called prior to assembly for time-dependent problems.
+     */
+    void setPreviousSolution(const SolutionVector& u)
+    { localResidual_.setPreviousSolution(u); }
+
+    /*!
+     * \brief Return the solution that has been set as the previous one.
+     */
+    const SolutionVector& prevSol() const
+    { return localResidual_.prevSol(); }
+
     /*!
      * \brief Resizes the matrix and right hand side and sets the matrix' sparsity pattern.
      */

dumux/implicit/cellcentered/localassembler.hh

@@ -70,10 +70,10 @@ public:
      */
     template<class Assembler>
     static void assemble(Assembler& assembler, typename Assembler::ResidualType& res, const Element& element,
-                         const SolutionVector& curSol, const SolutionVector& prevSol)
+                         const SolutionVector& curSol)
     {
         const auto globalI = assembler.fvGridGeometry().elementMapper().index(element);
-        res[globalI] = Implementation::assemble_(assembler, element, curSol, prevSol);
+        res[globalI] = Implementation::assemble_(assembler, element, curSol);
     }
 
     /*!
@@ -82,9 +82,9 @@ public:
      */
     template<class Assembler>
     static void assemble(Assembler& assembler, const Element& element,
-                         const SolutionVector& curSol, const SolutionVector& prevSol)
+                         const SolutionVector& curSol)
     {
-        Implementation::assemble_(assembler, element, curSol, prevSol);
+        Implementation::assemble_(assembler, element, curSol);
     }
 
     /*!
@@ -119,8 +119,7 @@ private:
      * \return The element residual at the current solution.
      */
     template<class Assembler>
-    static NumEqVector assemble_(Assembler& assembler, const Element& element,
-                                 const SolutionVector& curSol, const SolutionVector& prevSol)
+    static NumEqVector assemble_(Assembler& assembler, const Element& element, const SolutionVector& curSol)
     {
         // get some references for convenience
         const auto& problem = assembler.problem();
@@ -137,12 +136,14 @@ private:
         auto curElemVolVars = localView(gridVariables.curGridVolVars());
         curElemVolVars.bind(element, fvGeometry, curSol);
 
-        auto prevElemVolVars = localView(gridVariables.prevGridVolVars());
-        prevElemVolVars.bindElement(element, fvGeometry, prevSol);
-
         auto elemFluxVarsCache = localView(gridVariables.gridFluxVarsCache());
         elemFluxVarsCache.bind(element, fvGeometry, curElemVolVars);
 
+        const bool isStationary = localResidual.isStationary();
+        auto prevElemVolVars = localView(gridVariables.prevGridVolVars());
+        if (!isStationary)
+            prevElemVolVars.bindElement(element, fvGeometry, localResidual.prevSol());
+
         // the global dof of the actual element
         const auto globalI = fvGridGeometry.elementMapper().index(element);
 
@@ -157,13 +158,28 @@ private:
         // the actual element's current residual
         NumEqVector residual(0.0);
         if (!isGhost)
-            residual = localResidual.eval(problem,
-                                          element,
-                                          fvGeometry,
-                                          prevElemVolVars,
-                                          curElemVolVars,
-                                          elemBcTypes,
-                                          elemFluxVarsCache)[0];
+        {
+            if (isStationary)
+            {
+                residual = localResidual.eval(problem,
+                                              element,
+                                              fvGeometry,
+                                              curElemVolVars,
+                                              elemBcTypes,
+                                              elemFluxVarsCache)[0];
+            }
+            else
+            {
+                residual = localResidual.eval(problem,
+                                              element,
+                                              fvGeometry,
+                                              prevElemVolVars,
+                                              curElemVolVars,
+                                              elemBcTypes,
+                                              elemFluxVarsCache)[0];
+            }
+        }
+
 
         // TODO Do we really need this??????????
         // this->model_().updatePVWeights(fvGeometry);
@@ -248,13 +264,27 @@ private:
 
                 // calculate the residual with the deflected primary variables
                 if (!isGhost)
-                    partialDeriv = localResidual.eval(problem,
+                {
+                    if (isStationary)
+                    {
+                        partialDeriv = localResidual.eval(problem,
                                                       element,
                                                       fvGeometry,
-                                                      prevElemVolVars,
                                                       curElemVolVars,
                                                       elemBcTypes,
                                                       elemFluxVarsCache)[0];
+                    }
+                    else
+                    {
+                        partialDeriv = localResidual.eval(problem,
+                                                          element,
+                                                          fvGeometry,
+                                                          prevElemVolVars,
+                                                          curElemVolVars,
+                                                          elemBcTypes,
+                                                          elemFluxVarsCache)[0];
+                    }
+                }
 
                 // calculate the fluxes in the neighbors with the deflected primary variables
                 for (std::size_t k = 0; k < numNeighbors; ++k)
@@ -293,13 +323,27 @@ private:
 
                 // calculate the residual with the deflected primary variables and subtract it
                 if (!isGhost)
-                    partialDeriv -= localResidual.eval(problem,
-                                                       element,
-                                                       fvGeometry,
-                                                       prevElemVolVars,
-                                                       curElemVolVars,
-                                                       elemBcTypes,
-                                                       elemFluxVarsCache)[0];
+                {
+                    if (isStationary)
+                    {
+                        partialDeriv -= localResidual.eval(problem,
+                                                           element,
+                                                           fvGeometry,
+                                                           curElemVolVars,
+                                                           elemBcTypes,
+                                                           elemFluxVarsCache)[0];
+                    }
+                    else
+                    {
+                        partialDeriv -= localResidual.eval(problem,
+                                                           element,
+                                                           fvGeometry,
+                                                           prevElemVolVars,
+                                                           curElemVolVars,
+                                                           elemBcTypes,
+                                                           elemFluxVarsCache)[0];
+                    }
+                }
 
                 // calculate the fluxes into element with the deflected primary variables
                 for (std::size_t k = 0; k < numNeighbors; ++k)

dumux/implicit/cellcentered/localresidual.hh

@@ -55,6 +55,7 @@ class CCLocalResidual : public ImplicitLocalResidual<TypeTag>
     using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
 
 public:
+    using ParentType::ParentType;
 
     void evalFlux(ElementResidualVector& residual,
                   const Problem& problem,
@@ -82,7 +83,7 @@ public:
         // inner faces
         if (!scvf.boundary())
         {
-            flux += this->asImp_().computeFlux(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
+            flux += this->asImp().computeFlux(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);