LinearAlgebra backend for istl and mtl

CMakeLists.txt

@@ -7,6 +7,8 @@ set(CMAKE_PREFIX_PATH /opt/software/dune/lib/cmake)
 set(UG_DIR /opt/software/dune/lib/cmake/ug)
 set(Vc_DIR /opt/software/dune/lib/cmake/Vc)
 
+set(MTL_DIR /opt/sources/amdis2/lib/mtl4)
+
 if(NOT (dune-common_DIR OR dune-common_ROOT OR
       "${CMAKE_PREFIX_PATH}" MATCHES ".*dune-common.*"))
     string(REPLACE  ${CMAKE_PROJECT_NAME} dune-common dune-common_DIR

dune/amdis/AdaptBase.cpp

+#include "AdaptBase.hpp"
+
+namespace AMDiS
+{
+  int AdaptBase::info = 10;
+
+} // end namespace AMDiS

dune/amdis/AdaptBase.hpp

+#pragma once
+
+// std c++ headers
+#include <string>
+
+namespace AMDiS
+{
+  // forward declarations
+  class AdaptInfo;
+  class ProblemIterationInterface;
+  class ProblemTimeInterface;
+  
+  /// Interface for adaption loops.
+  class AdaptBase
+  {
+  public:
+    /// Constructor
+    AdaptBase(std::string sname,
+              ProblemIterationInterface* problemIteration,
+              AdaptInfo& adapt,
+              ProblemTimeInterface* problemTime = NULL,
+              AdaptInfo* initialAdaptInfo = NULL)
+      : name(sname)
+      , problemIteration(problemIteration)
+      , adaptInfo(adapt)
+      , problemTime(problemTime)
+      , initialAdaptInfo(initialAdaptInfo)
+    {}
+
+    /// Destructor
+    virtual ~AdaptBase() {}
+
+    /** \brief
+     * Pure virtual method. Must be overloaded by sub classes to perform
+     * a concrete adaption loop.
+     */
+    virtual int adapt() = 0;
+
+    /// Returns \ref name
+    std::string getName() const
+    {
+      return name;
+    }
+
+    /// Returns \ref problemIteration_
+    ProblemIterationInterface* getProblemIteration() const
+    {
+      return problemIteration;
+    }
+
+    ///
+    void setProblemIteration(ProblemIterationInterface* pii)
+    {
+      problemIteration = pii;
+    }
+
+    /// Returns \ref adaptInfo
+    AdaptInfo& getAdaptInfo() const
+    {
+      return adaptInfo;
+    }
+
+    /// Returns \ref problemTime_
+    ProblemTimeInterface* getProblemTime() const
+    {
+      return problemTime;
+    }
+
+    ///
+    void setProblemTime(ProblemTimeInterface* pti)
+    {
+      problemTime = pti;
+    }
+
+    /// Returns \ref initialAdaptInfo_
+    AdaptInfo& getInitialAdaptInfo() const
+    {
+      return *initialAdaptInfo;
+    }
+
+  protected:
+    /// Name of the adaption loop
+    std::string name;
+
+    /// Problem iteration interface
+    ProblemIterationInterface* problemIteration;
+
+    /// Main adapt info
+    AdaptInfo& adaptInfo;
+
+    /// problem time interface
+    ProblemTimeInterface* problemTime;
+
+    /** \brief
+     * Adapt info for initial adapt. Will be given to
+     * problemTime_->solveInitialProblem().
+     */
+    AdaptInfo* initialAdaptInfo;
+
+    /// Info level
+    static int info;
+  };
+
+} // end namespace AMDiS

dune/amdis/AdaptInstationary.cpp

+#include "AdaptInstationary.hpp"
+
+// AMDiS includes
+#include "AdaptInfo.hpp"
+#include "Flag.hpp"
+#include "Initfile.hpp"
+#include "Log.hpp"
+#include "ProblemIterationInterface.hpp"
+#include "ProblemTimeInterface.hpp"
+
+namespace AMDiS
+{
+  AdaptInstationary::AdaptInstationary(std::string name,
+                                       ProblemIterationInterface& problemStat,
+                                       AdaptInfo& adaptInfo,
+                                       ProblemTimeInterface& problemInstat,
+                                       AdaptInfo& initialInfo)
+    : AdaptBase(name, &problemStat, adaptInfo, &problemInstat, &initialInfo),
+      breakWhenStable(0)
+  {
+    strategy = 0;
+    timeDelta1 = 0.7071;
+    timeDelta2 = 1.4142;
+
+    Parameters::get(name + "->strategy", strategy);
+    Parameters::get(name + "->time delta 1", timeDelta1);
+    Parameters::get(name + "->time delta 2", timeDelta2);
+    Parameters::get(name + "->info", info);
+    Parameters::get(name + "->break when stable", breakWhenStable);
+
+    fixedTimestep = (adaptInfo.getMinTimestep() == adaptInfo.getMaxTimestep());
+  }
+
+
+  void AdaptInstationary::explicitTimeStrategy()
+  {
+    AMDIS_FUNCNAME("AdaptInstationary::explicitTimeStrategy()");
+
+    // estimate before first adaption
+    if (adaptInfo.getTime() <= adaptInfo.getStartTime())
+      problemIteration->oneIteration(adaptInfo, ESTIMATE);
+
+
+    // increment time
+    adaptInfo.setTime(adaptInfo.getTime() + adaptInfo.getTimestep());
+
+    problemTime->setTime(adaptInfo);
+
+    AMDIS_MSG("time = " << adaptInfo.getTime() << ", timestep = "
+                        << adaptInfo.getTimestep());
+
+    adaptInfo.setSpaceIteration(0);
+
+    // do the iteration
+    problemIteration->beginIteration(adaptInfo);
+    problemIteration->oneIteration(adaptInfo, FULL_ITERATION);
+    problemIteration->endIteration(adaptInfo);
+    adaptInfo.setLastProcessedTimestep(adaptInfo.getTimestep());
+  }
+
+
+  void AdaptInstationary::implicitTimeStrategy()
+  {
+    AMDIS_FUNCNAME("AdaptInstationary::implicitTimeStrategy()");
+
+    do
+    {
+      adaptInfo.setTime(adaptInfo.getTime() + adaptInfo.getTimestep());
+      problemTime->setTime(adaptInfo);
+
+      AMDIS_MSG("time = " << adaptInfo.getTime() << ", timestep = "
+                          << adaptInfo.getTimestep());
+
+      problemIteration->oneIteration(adaptInfo, NO_ADAPTION);
+
+      adaptInfo.incTimestepIteration();
+
+      if (!fixedTimestep &&
+          !adaptInfo.timeToleranceReached() &&
+          adaptInfo.getTimestepIteration() <= adaptInfo.getMaxTimestepIteration() &&
+          !(adaptInfo.getTimestep() <= adaptInfo.getMinTimestep()))
+      {
+        adaptInfo.setTime(adaptInfo.getTime() - adaptInfo.getTimestep());
+        adaptInfo.setTimestep(adaptInfo.getTimestep() * timeDelta1);
+        continue;
+      }
+
+
+      adaptInfo.setSpaceIteration(0);
+
+
+      // === Do space iterations only if the maximum is higher than 0. ===
+
+      if (adaptInfo.getMaxSpaceIteration() > 0)
+      {
+
+        // === Space iterations. ===
+        do
+        {
+          problemIteration->beginIteration(adaptInfo);
+
+          Flag adapted = problemIteration->oneIteration(adaptInfo, FULL_ITERATION);
+          if (adapted == Flag{0})
+          {
+            if (!fixedTimestep &&
+                !adaptInfo.timeToleranceReached() &&
+                !(adaptInfo.getTimestep() <= adaptInfo.getMinTimestep()))
+            {
+              adaptInfo.setTime(adaptInfo.getTime() - adaptInfo.getTimestep());
+              adaptInfo.setTimestep(adaptInfo.getTimestep() * timeDelta2);
+              problemIteration->endIteration(adaptInfo);
+              adaptInfo.incSpaceIteration();
+              break;
+            }
+          }
+
+          adaptInfo.incSpaceIteration();
+          problemIteration->endIteration(adaptInfo);
+
+        }
+        while (!adaptInfo.spaceToleranceReached() &&
+               adaptInfo.getSpaceIteration() <= adaptInfo.getMaxSpaceIteration());
+
+      }
+      else
+      {
+        problemIteration->endIteration(adaptInfo);
+      }
+
+
+    }
+    while(!adaptInfo.timeToleranceReached() &&
+          !(adaptInfo.getTimestep() <= adaptInfo.getMinTimestep()) &&
+          adaptInfo.getTimestepIteration() <= adaptInfo.getMaxTimestepIteration());
+
+    adaptInfo.setLastProcessedTimestep(adaptInfo.getTimestep());
+
+    // After successful iteration/timestep the timestep will be changed according
+    // adaption rules for next timestep.
+    // First, check for increase of timestep
+    if (!fixedTimestep && adaptInfo.timeErrorLow())
+      adaptInfo.setTimestep(adaptInfo.getTimestep() * timeDelta2);
+
+    // Second, check for decrease of timestep
+    if (!fixedTimestep &&
+        !adaptInfo.timeToleranceReached() &&
+        !(adaptInfo.getTimestep() <= adaptInfo.getMinTimestep()))
+      adaptInfo.setTimestep(adaptInfo.getTimestep() * timeDelta1);
+  }
+
+
+  void AdaptInstationary::simpleAdaptiveTimeStrategy()
+  {
+    AMDIS_FUNCNAME("AdaptInstationary::simpleAdaptiveTimeStrategy()");
+
+    // estimate before first adaption
+    if (adaptInfo.getTime() <= adaptInfo.getStartTime())
+      problemIteration->oneIteration(adaptInfo, ESTIMATE);
+
+    adaptInfo.setTime(adaptInfo.getTime() + adaptInfo.getTimestep());
+    problemTime->setTime(adaptInfo);
+
+    AMDIS_MSG("time = " << adaptInfo.getTime() << ", timestep = "
+                        << adaptInfo.getTimestep());
+
+    problemIteration->oneIteration(adaptInfo, FULL_ITERATION);
+
+    adaptInfo.setLastProcessedTimestep(adaptInfo.getTimestep());
+
+    // First, check for increase of timestep
+    if (!fixedTimestep && adaptInfo.timeErrorLow())
+      adaptInfo.setTimestep(adaptInfo.getTimestep() * timeDelta2);
+
+    // Second, check for decrease of timestep
+    if (!fixedTimestep &&
+        !adaptInfo.timeToleranceReached() &&
+        !(adaptInfo.getTimestep() <= adaptInfo.getMinTimestep()))
+      adaptInfo.setTimestep(adaptInfo.getTimestep() * timeDelta1);
+  }
+
+
+  void AdaptInstationary::oneTimestep()
+  {
+    AMDIS_FUNCNAME("AdaptInstationary::oneTimestep()");
+
+    adaptInfo.setTimestepIteration(0);
+
+    switch (strategy)
+    {
+    case 0:
+      explicitTimeStrategy();
+      break;
+    case 1:
+      implicitTimeStrategy();
+      break;
+    case 2:
+      simpleAdaptiveTimeStrategy();
+      break;
+    default:
+      AMDIS_ERROR_EXIT("Unknown strategy = " << strategy);
+    }
+
+    adaptInfo.incTimestepNumber();
+  }
+
+
+  int AdaptInstationary::adapt()
+  {
+    AMDIS_FUNCNAME("AdaptInstationary::adapt()");
+    int errorCode = 0;
+
+    AMDIS_TEST_EXIT(adaptInfo.getTimestep() >= adaptInfo.getMinTimestep(),
+      "timestep < min timestep");
+    AMDIS_TEST_EXIT(adaptInfo.getTimestep() <= adaptInfo.getMaxTimestep(),
+      "timestep > max timestep");
+
+    AMDIS_TEST_EXIT(adaptInfo.getTimestep() > 0, "timestep <= 0!");
+
+    if (adaptInfo.getTimestepNumber() == 0)
+    {
+      adaptInfo.setTime(adaptInfo.getStartTime());
+      initialAdaptInfo->setStartTime(adaptInfo.getStartTime());
+      initialAdaptInfo->setTime(adaptInfo.getStartTime());
+
+      problemTime->setTime(adaptInfo);
+
+      // initial adaption
+      problemTime->solveInitialProblem(*initialAdaptInfo);
+      problemTime->transferInitialSolution(adaptInfo);
+    }
+
+    while (!adaptInfo.reachedEndTime())
+    {
+      problemTime->initTimestep(adaptInfo);
+      oneTimestep();
+      problemTime->closeTimestep(adaptInfo);
+
+      if (breakWhenStable && (adaptInfo.getSolverIterations() == 0))
+        break;
+    }
+
+    return errorCode;
+  }
+
+} // end namespace AMDiS

dune/amdis/AdaptInstationary.hpp

+#pragma once
+
+// std c++ headers
+#include <string>
+
+// AMDiS includes
+#include "AdaptBase.hpp"
+
+namespace AMDiS
+{
+  // forward declarations
+  class AdaptInfo;
+  class ProblemInterationInterface;
+  class ProblemTimeInterface;
+
+  /** \ingroup Adaption
+   * \brief
+   * AdaptInstationary implements the adaptive procdure for time dependent
+   * problems (see ProblemInstat). It contains a pointer to a ProblemInstat
+   * object.
+   */
+  class AdaptInstationary : public AdaptBase
+  {
+  public:
+    /// Creates a AdaptInstationary object with the given name for the time
+    /// dependent problem problemInstat.
+    AdaptInstationary(std::string name,
+                      ProblemIterationInterface& problemStat,
+                      AdaptInfo& info,
+                      ProblemTimeInterface& problemInstat,
+                      AdaptInfo& initialInfo);
+
+    /// Sets \ref strategy to aStrategy
+    void setStrategy(int aStrategy)
+    {
+      strategy = aStrategy;
+    }
+
+    /// Returns \ref strategy
+    int getStrategy() const
+    {
+      return strategy;
+    }
+
+    /// Implementation of AdaptBase::adapt()
+    virtual int adapt() override;
+
+  protected:
+    /** \brief
+     * Implements one (maybe adaptive) timestep. Both the explicit and the
+     * implicit time strategy are implemented. The semi-implicit strategy
+     * is only a special case of the implicit strategy with a limited number of
+     * iterations (exactly one).
+     * The routine uses the parameter \ref strategy to select the strategy:
+     * strategy 0: Explicit strategy,
+     * strategy 1: Implicit strategy.
+     */
+    virtual void oneTimestep();
+
+    /// Initialisation of this AdaptInstationary object
+    void initialize(std::string aName);
+
+    /// Implements the explit time strategy. Used by \ref oneTimestep().
+    virtual void explicitTimeStrategy();
+
+    /// Implements the implicit time strategy. Used by \ref oneTimestep().
+    virtual void implicitTimeStrategy();
+
+    /** \brief
+     * This iteration strategy allows the timestep and the mesh to be adapted
+     * after each timestep solution. There are no inner loops for mesh adaption and
+     * no refused timesteps.
+     */
+    void simpleAdaptiveTimeStrategy();
+
+  protected:
+    /// Strategy for choosing one timestep
+    int strategy;
+
+    /// Parameter \f$ \delta_1 \f$ used in time step reduction
+    double timeDelta1;
+
+    /// Parameter \f$ \delta_2 \f$ used in time step enlargement
+    double timeDelta2;
+
+    /// If this parameter is 1 and the instationary problem is stable, hence the number
+    /// of solver iterations to solve the problem is zero, the adaption loop will stop.
+    int breakWhenStable;
+
+    ///
+    bool fixedTimestep;
+  };
+
+} // end namespace AMDiS

dune/amdis/AdaptStationary.cpp

+#include "AdaptStationary.hpp"
+
+// AMDiS includes
+#include "AdaptInfo.hpp"
+#include "Flag.hpp"
+#include "Initfile.hpp"
+#include "ProblemIterationInterface.hpp"
+
+namespace AMDiS
+{
+
+  AdaptStationary::AdaptStationary(std::string name,
+                                   ProblemIterationInterface& prob,
+                                   AdaptInfo& adaptInfo)
+    : AdaptBase(name, &prob, adaptInfo)
+  {
+    Parameters::get(name + "->info", info);
+  }
+
+
+  int AdaptStationary::adapt()
+  {
+    // initial iteration
+    if (adaptInfo.getSpaceIteration() == -1)
+    {
+      problemIteration->beginIteration(adaptInfo);
+      problemIteration->oneIteration(adaptInfo, NO_ADAPTION);
+      problemIteration->endIteration(adaptInfo);
+      adaptInfo.incSpaceIteration();
+    }
+
+    // adaption loop
+    while (!adaptInfo.spaceToleranceReached() &&
+           (adaptInfo.getSpaceIteration() < adaptInfo.getMaxSpaceIteration() ||
+            adaptInfo.getMaxSpaceIteration() < 0) )
+    {
+
+      problemIteration->beginIteration(adaptInfo);
+      Flag adapted = problemIteration->oneIteration(adaptInfo, FULL_ITERATION);
+      problemIteration->endIteration(adaptInfo);
+
+      if (adapted == Flag{0})
+        break;
+
+      adaptInfo.incSpaceIteration();
+    }
+
+    return 0;
+  }
+
+} // end namespace AMDiS

dune/amdis/AdaptStationary.hpp

+/** \defgroup Adaption Adaption module
+ * @{ <img src="adaption.png"> @}
+ *
+ * \brief
+ * Contains all classes needed for adaption.
+ */
+
+#pragma once
+
+// std c++ headers
+#include <string>
+
+// AMDiS includes
+#include "AdaptBase.hpp"
+
+namespace AMDiS
+{
+  // forward declarations
+  class AdaptInfo;
+  class ProblemIterationInterface;
+
+  /** \ingroup Adaption
+   * \brief
+   * AdaptStationary contains information about the adaptive procedure and the
+   * adapt procedure itself
+   */
+  class AdaptStationary : public AdaptBase
+  {
+  public:
+    /// Creates a AdaptStationary object with given name.
+    AdaptStationary(std::string name,
+                    ProblemIterationInterface& prob,
+                    AdaptInfo& info);
+
+    /// Implementation of AdaptBase::adapt()
+    virtual int adapt() override;
+  };
+
+} // end namespace AMDiS

dune/amdis/Basic.hpp

 #pragma once
 
 #include <cassert>
-#include <tuple>
-#include <map>
 #include <list>
+#include <map>
+#include <memory>
+#include <tuple>
+#include <type_traits>
+#include <vector>
 
 #include "IndexSeq.hpp"
 

dune/amdis/CMakeLists.txt

 #install headers
 
 dune_add_library("duneamdis" NO_EXPORT
+    AdaptBase.cpp
     AdaptInfo.cpp
+    AdaptInstationary.cpp
+    AdaptStationary.cpp
     AMDiS.cpp
     Initfile.cpp
+    ProblemInstatBase.cpp
+    ProblemInstat.cpp
     ProblemStat.cpp
-    SystemVector.cpp
+    StandardProblemIteration.cpp
+    linear_algebra/istl/SystemMatrix.cpp
+    linear_algebra/istl/SystemVector.cpp
+    linear_algebra/mtl/SystemMatrix.cpp
+    linear_algebra/mtl/SystemVector.cpp
     )
 add_dune_alberta_flags("duneamdis" OBJECT USE_GENERIC)
+target_compile_definitions("duneamdis" PUBLIC AMDIS_BACKEND_MTL=1)
 
 set(BOOST_VERSION "1.54")
 set(BOOST_LIBS_REQUIRED system program_options)
@@ -17,25 +27,40 @@ endif (NOT BUILD_SHARED_LIBS)
 set(Boost_NO_SYSTEM_PATHS ON)
 find_package(Boost ${BOOST_VERSION} REQUIRED ${BOOST_LIBS_REQUIRED})
 if (Boost_FOUND)
-    add_library(boost INTERFACE)
-    target_include_directories(boost INTERFACE ${Boost_INCLUDE_DIR})
-    target_link_libraries(boost INTERFACE  ${Boost_LIBRARIES})
-    target_link_libraries("duneamdis" INTERFACE boost)
-
-    if (MSVC_SHARED_LIBS)
-        link_directories(${Boost_LIBRARY_DIRS})
-        target_compile_definitions("duneamdis" INTERFACE ${Boost_LIB_DIAGNOSTIC_DEFINITIONS})
-    endif (MSVC_SHARED_LIBS)
+  add_library(boost INTERFACE)
+  target_include_directories(boost INTERFACE ${Boost_INCLUDE_DIR})
+  target_link_libraries(boost INTERFACE  ${Boost_LIBRARIES})
+  target_link_libraries("duneamdis" INTERFACE boost)
+
+  if (MSVC_SHARED_LIBS)
+    link_directories(${Boost_LIBRARY_DIRS})
+    target_compile_definitions("duneamdis" INTERFACE ${Boost_LIB_DIAGNOSTIC_DEFINITIONS})
+  endif (MSVC_SHARED_LIBS)
 endif (Boost_FOUND)
 
 
+find_package(MTL REQUIRED)
+if (MTL_FOUND)
+  target_include_directories("duneamdis" PUBLIC ${MTL_INCLUDE_DIRS})
+#   target_link_libraries("duneamdis" PUBLIC ${MTL_LIBRARIES})
+#   target_compile_options("duneamdis" PUBLIC ${MTL_CXX_DEFINITIONS})