umfpack solver added. TODO: seems not to work yet. Also the integration in the solver-map is ugly.

dune/amdis/CMakeLists.txt

@@ -49,6 +49,10 @@ if (MTL_FOUND)
   foreach (feature ${CXX_ELEVEN_FEATURE_LIST})
       target_compile_definitions("duneamdis" PUBLIC MTL_WITH_${feature})
   endforeach ()
+  
+  if (HAVE_UMFPACK OR ENABLE_SUITESPARSE OR SuiteSparse_FOUND)
+      target_compile_definitions("duneamdis" PUBLIC MTL_HAS_UMFPACK)
+  endif (HAVE_UMFPACK)
 endif (MTL_FOUND)
 
 

dune/amdis/linear_algebra/mtl/BITL_Solver.hpp

@@ -38,6 +38,10 @@ namespace AMDiS
         return std::make_unique<Solver<gcr_type>>(prefix);
 //       else if (name == "preonly")
 //         return std::make_unique<Solver<preonly_type>>(prefix);
+#ifdef HAVE_UMFPACK
+      else if (name == "umfpack" || name == "direct")
+        return std::make_unique< LinearSolver<Matrix, Vector, UmfpackRunner<Matrix, BaseVector<Vector>>> >(prefix);
+#endif
       else
         AMDIS_ERROR_EXIT("Unknown Solver-name!");
     }

dune/amdis/linear_algebra/mtl/BlockMTLMatrix.hpp

-/** \file BlockMTLMatrix.h */
+/** \file BlockMTLMatrix.hpp */
 
 #pragma once
 
@@ -24,6 +24,9 @@ namespace AMDiS
     
     /// The type of the elements of the MTLMatrix
     using value_type = typename MTLMatrix::value_type;
+    
+    /// The underlying mtl matrix type
+    using BaseMatrix = MTLMatrix;
 
   public:
     /// Return the (R,C)'th matrix block

dune/amdis/linear_algebra/mtl/Copy.hpp

+#pragma once
+
+#include <dune/amdis/linear_algebra/mtl/BlockMTLVector.hpp>
+#include <dune/amdis/linear_algebra/mtl/BlockMTLMatrix.hpp>
+
+#include <dune/amdis/linear_algebra/mtl/SystemVector.hpp>
+#include <dune/amdis/linear_algebra/mtl/SystemMatrix.hpp>
+
+#include <dune/amdis/linear_algebra/mtl/Mapper.hpp>
+
+namespace AMDiS 
+{
+  // requires MatrixType to have an inserter
+  template <class Value, size_t _N, size_t _M, class TargetType>
+  void copy(BlockMTLMatrix<Value, _N, _M> const& source, TargetType& target)
+  {
+    using Mapper = std::conditional_t<(_N == _M), BlockMapper, RectangularMapper>;
+    using value_type = typename TargetType::value_type;
+    using BaseInserter = mtl::mat::inserter<TargetType, mtl::operations::update_plus<value_type>>;
+    using MappedInserter = typename Mapper::template Inserter<BaseInserter>::type;
+      
+    size_t nnz = 0;
+    for (size_t rb = 0; rb < _N; ++rb)
+      for (size_t cb = 0; cb < _M; ++cb)
+        nnz += source[rb][cb].nnz();	  
+      
+    {
+      target.change_dim(num_rows(source), num_cols(source));
+      set_to_zero(target);
+      
+      Mapper mapper(source);
+      MappedInserter ins(target, mapper, int(1.2 * nnz / target.dim1()));
+      
+      for (size_t rb = 0; rb < _N; ++rb) {
+        mapper.setRow(rb);
+        for (size_t cb = 0; cb < _M; ++cb) {
+	  mapper.setCol(cb);
+          ins << source[rb][cb];
+        }
+      }
+    }
+  }
+  
+  template <class FeSpaces, class Value, class TargetType>
+  void copy(SystemMatrix<FeSpaces, Value> const& source, TargetType& target)
+  {
+    copy(source.getMatrix(), target);
+  }
+  
+} // end namespace AMDiS
\ No newline at end of file

dune/amdis/linear_algebra/mtl/ITL_Solver.hpp

 /** \file ITL_Solver.h */
 
 #pragma once
+#ifdef HAVE_CONFIG_H
+# include "config.h"
+#endif
 
 // MTL4 includes
 #include <boost/numeric/itl/itl.hpp>
@@ -27,6 +30,7 @@
 #include <dune/amdis/Initfile.hpp>
 #include <dune/amdis/linear_algebra/mtl/LinearSolver.hpp>
 #include <dune/amdis/linear_algebra/mtl/KrylovRunner.hpp>
+#include <dune/amdis/linear_algebra/mtl/UmfpackRunner.hpp>
 
 namespace AMDiS
 {
@@ -482,6 +486,10 @@ namespace AMDiS
         return std::make_unique<Solver<gcr_type>>(prefix);
       else if (name == "preonly")
         return std::make_unique<Solver<preonly_type>>(prefix);
+#ifdef HAVE_UMFPACK
+      else if (name == "umfpack" || name == "direct")
+        return std::make_unique< LinearSolver<Matrix, Vector, UmfpackRunner<Matrix, BaseVector<Vector>>> >(prefix);
+#endif
       else
         AMDIS_ERROR_EXIT("Unknown Solver-name!");
     }

dune/amdis/linear_algebra/mtl/MTLDenseVector.hpp

@@ -40,6 +40,7 @@ namespace AMDiS
     
   } // end namespace Impl
   
+  /// wrap the mtl dense_vector class to provide resize and size member-functions
   template <class Value>
   Impl::MTLWrapper<MTLDenseVector<Value>> wrapper(MTLDenseVector<Value>& vec)
   {

dune/amdis/linear_algebra/mtl/Mapper.hpp

+#pragma once
+
+#include <vector>
+
+#include <boost/numeric/mtl/matrix/mapped_inserter.hpp>
+
+#include <dune/amdis/linear_algebra/mtl/BlockMTLMatrix.hpp>
+#include <dune/amdis/linear_algebra/mtl/SystemMatrix.hpp>
+
+namespace AMDiS 
+{
+  
+  /**
+   * \brief BaseClass for all mapper types
+   * 
+   * A mapper assigned a local row/column index of a cell of
+   * a block matrix to a global matrix index.
+   * Call \ref setRow and \ref setCol first, to select a block
+   * in the block matrix. Then get with \ref row, respective \ref col,
+   * the global matrix index to the assigned local index.
+   **/  
+  template <class Derived>
+  struct MapperBase
+  {
+    using size_type = std::size_t;
+    
+    template <class BaseInserter>
+    struct Inserter {
+      using type = mtl::mat::mapped_inserter<BaseInserter, Derived>;
+    };
+    
+    /// set the current block row
+    void setRow(size_type rb) { self().setRow(rb); }
+    
+    /// set the current block columns
+    void setCol(size_type cb) { self().setCol(cb); }
+    
+    /// return global matrix row, for local row in the current matrix block
+    size_type row(size_type r) const { return self().row(r); }
+    
+    /// return global matrix column, for local column in the current matrix block
+    size_type col(size_type c) const { return self().col(c); }
+    
+    /// return overall number of rows
+    size_type getNumRows() const { return self().getNumRows(); }
+    
+    /// return overall number of columns
+    size_type getNumCols() const { return self().getNumCols(); }
+
+    size_type getNumRows(size_type comp) const { return self().getNumRows(comp); }
+    size_type getNumCols(size_type comp) const { return self().getNumCols(comp); }
+    
+    /// return number of components/blocks
+    size_type getNumComponents() const { return self().getNumComponents(); }
+    
+    size_type row(size_type rb, size_type cb, size_type r)
+    {
+      setRow(rb);
+      setCol(cb);
+      return row(r);
+    }
+        
+    size_type col(size_type rb, size_type cb, size_type c)
+    {
+      setRow(rb);
+      setCol(cb);
+      return col(c);
+    }
+    
+    Derived&       self()       { return static_cast<Derived&>(*this); }
+    Derived const& self() const { return static_cast<Derived const&>(*this); }
+  };
+
+  
+  /// Mapper implementation for non-parallel block matrices
+  struct BlockMapper : public MapperBase<BlockMapper>
+  {
+    /// Default constructor
+    BlockMapper() = default;
+    BlockMapper(BlockMapper const& other) = default;
+    BlockMapper& operator=(BlockMapper const& other) = default;
+    
+    /// Constructor for BlockMTLMatrix
+    template <class Value, size_t _N>
+    BlockMapper(BlockMTLMatrix<Value, _N, _N> const& mat)
+      : nComp(_N) 
+      , rowOffset(0), colOffset(0)
+      , nrow(num_rows(mat)), ncol(nrow)
+      , sizes(nComp)
+    {      
+      for (size_t i = 0; i < _N; ++i)
+        sizes[i] = num_rows(mat[i][i]);
+    }
+    
+    /// Constructor for SystemMatrix
+    template <class FeSpaces, class Value>
+    BlockMapper(SystemMatrix<FeSpaces, Value> const& mat)
+      : BlockMapper(mat.getMatrix())
+    {}
+    
+    /// Constructor for single DOFMatrix
+    template <class FeSpace, class Value>
+    BlockMapper(DOFMatrix<FeSpace, Value> const& mat)
+      : nComp(1) 
+      , rowOffset(0), colOffset(0)
+      , nrow(0), ncol(0)
+      , sizes(nComp)
+    {
+      sizes[0] = mat.N();
+      nrow += sizes[0];
+      ncol = nrow;
+    }
+
+    /// Constructor for system with equal components
+    BlockMapper(size_type nComp, size_type nDOFperComp)
+    : nComp(nComp), 
+      rowOffset(0), colOffset(0), 
+      nrow(nComp*nDOFperComp), ncol(nrow),
+      sizes(nComp)
+    {
+       for (size_type i = 0; i < nComp; ++i)
+         sizes[i] = nDOFperComp;
+    }
+      
+    /// calculate row offset for row component \param r
+    void setRow(size_type rb)
+    { 
+      AMDIS_TEST_EXIT_DBG(rb <= sizes.size(), "row nr out of range!");
+      rowOffset = sum(rb); 
+    }
+
+    /// calculate column offset for col component \param c
+    void setCol(size_type cb)
+    {
+      AMDIS_TEST_EXIT_DBG(rb <= sizes.size(), "column nr out of range!");
+      colOffset = sum(cb);
+    }
+
+    size_type row(size_type r) const { return r + rowOffset; }
+    size_type col(size_type c) const { return c + colOffset; }
+
+    size_type getNumRows() const { return nrow; }
+    size_type getNumCols() const { return ncol; }
+
+    size_type getNumRows(size_type comp) const { return sizes[comp]; }
+    size_type getNumCols(size_type comp) const { return sizes[comp]; }
+
+    size_type getNumComponents() const { return nComp; }
+
+  private: // methods
+
+    ///compute the sum of sizes from [0, end)
+    size_type sum(size_type end) const 
+    {
+      size_type ret = 0;
+      for (size_type i = 0; i < end; ++i)
+	ret += sizes[i];
+      return ret;
+    }
+    
+  private: // variables
+
+    size_type nComp = 0;
+    size_type rowOffset = 0;
+    size_type colOffset = 0;
+    size_type nrow = 0;
+    size_type ncol = 0;
+    
+    std::vector<size_type> sizes;
+  };
+
+  
+  /// Mapper implementation for non-parallel rectangular block matrices
+  struct RectangularMapper : public MapperBase<RectangularMapper>
+  {
+    /// Constructor for block-matrices
+    template <class Value, size_t _N, size_t _M>
+    RectangularMapper(BlockMTLMatrix<Value, _N, _M> const& mat)
+      : nRowComp(_N), nColComp(_M)
+      , rowOffset(0), colOffset(0)
+      , nrow(num_rows(mat)), ncol(num_cols(mat))
+      , sizes_rows(nRowComp), sizes_cols(nColComp)
+    {
+      for (size_type i = 0; i < _N; ++i)
+        sizes_rows[i] = num_rows(mat[i][0]);
+      for (size_type j = 0; j < _M; ++j)
+        sizes_rows[j] = num_rows(mat[0][j]);
+    }
+      
+    /// calculate row offset for row component \param r
+    void setRow(size_type rb)
+    { 
+      AMDIS_TEST_EXIT_DBG(rb <= sizes_rows.size(), "row nr out of range!");
+      rowOffset = sum(rb, sizes_rows); 
+    }
+
+    /// calculate column offset for col component \param c
+    void setCol(size_type cb)
+    {
+      AMDIS_TEST_EXIT_DBG(cb <= sizes_cols.size(), "column nr out of range!");
+      colOffset = sum(cb, sizes_cols);
+    }
+
+    size_type row(size_type r) const { return r + rowOffset; }
+    size_type col(size_type c) const { return c + colOffset; }
+
+    size_type getNumRows() const { return nrow; }
+    size_type getNumCols() const { return ncol; }
+
+    size_type getNumRows(size_type comp) const { return sizes_rows[comp]; }
+    size_type getNumCols(size_type comp) const { return sizes_cols[comp]; }
+
+    size_type getNumComponents() const { return nRowComp; }
+    size_type getNumRowComponents() const { return nRowComp; }
+    size_type getNumColComponents() const { return nColComp; }
+
+  private: // methods
+
+    ///compute the sum of sizes from [0, end)
+    size_type sum(size_type end, std::vector<size_type>& sizes) const 
+    {
+      size_type ret = 0;
+      for (size_type i = 0; i < end; ++i)
+	ret += sizes[i];
+      return ret;
+    }
+    
+  private: // variables
+
+    size_type nRowComp = 0;
+    size_type nColComp = 0;
+    size_type rowOffset = 0;
+    size_type colOffset = 0;
+    size_type nrow = 0;
+    size_type ncol = 0;
+    
+    std::vector<size_type> sizes_rows;
+    std::vector<size_type> sizes_cols;
+  };
+  
+} // end namespace AMDiS

dune/amdis/linear_algebra/mtl/UmfpackRunner.hpp

+#pragma once
+
+#ifdef HAVE_UMFPACK
+
+#include <iostream>
+
+#include <boost/numeric/mtl/operation/two_norm.hpp>
+#include <boost/numeric/mtl/interface/umfpack_solve.hpp>
+
+#include <dune/amdis/linear_algebra/RunnerInterface.hpp>
+#include <dune/amdis/linear_algebra/SolverInfo.hpp>
+
+#include <dune/amdis/linear_algebra/mtl/Copy.hpp>
+
+namespace AMDiS 
+{
+  using namespace mtl::mat::umfpack;
+  
+  /**
+   * \ingroup Solver
+   * \class AMDiS::UmfPackSolver
+   * \brief \implements LinearSolver
+   * Wrapper for the external UMFPACK solver:
+   *   http://www.cise.ufl.edu/research/sparse/umfpack/
+   *
+   * This is a direct solver for large sparse matrices.
+   */
+  template <class Matrix, class Vector>
+  class UmfpackRunner : public RunnerInterface<Matrix, Vector, Vector>
+  {
+    using Super = RunnerInterface<Matrix, Vector>;
+    using PreconBase = typename Super::PreconBase;
+    
+    // extract the matrix type used by UMFPACK
+    using FullMatrix = typename mtl::mat::umfpack::matrix_copy<
+        typename Matrix::BaseMatrix, 
+        typename Matrix::value_type, 
+        typename Matrix::BaseMatrix::parameters::orientation
+      >::matrix_type;
+    using SolverType = mtl::mat::umfpack::solver<FullMatrix>;
+    
+  public:
+    /// Constructor. Reads UMFPACK parameters from initfile
+    UmfpackRunner(std::string prefix)
+      : solver(NULL)
+      , store_symbolic(0)
+      , symmetric_strategy(0)
+      , alloc_init(0.7)
+    {
+      Parameters::get(prefix + "->store symbolic", store_symbolic); // ?
+      Parameters::get(prefix + "->symmetric strategy", symmetric_strategy);
+      Parameters::get(prefix + "->alloc init", alloc_init);
+    }
+
+    /// Destructor
+    ~UmfpackRunner() = default;
+
+    /// Implementation of \ref RunnerInterface::init()
+    virtual void init(Matrix const& A) override
+    {
+      AMDIS_MSG("copy matrix...");
+      copy(A, fullMatrix);
+      AMDIS_MSG("size(fullMatrix) = " << num_rows(fullMatrix) << " x " << num_cols(fullMatrix));
+      
+      {
+        mtl::io::matrix_market_ostream out("matrix.mtx");
+        out << fullMatrix;
+      }
+      
+      AMDIS_MSG("create factorization...");
+      try {
+	solver.reset(new SolverType(fullMatrix, symmetric_strategy, alloc_init));
+      } catch (mtl::mat::umfpack::error& e) {
+	AMDIS_ERROR_EXIT("UMFPACK_ERROR(factorize, " << e.code << ") = " << e.what());
+      }
+    }
+
+    /// Implementation of \ref RunnerBase::solve()
+    virtual int solve(Matrix const& A, Vector& x, Vector const& b, 
+                      SolverInfo& solverInfo) override
+    { 
+      AMDIS_FUNCNAME("Umfpack_Runner::solve()");
+      AMDIS_TEST_EXIT(solver, "The umfpack solver was not initialized\n");
+      
+      AMDIS_MSG("solve system...");
+      int code = 0;
+      try {
+	code = (*solver)(x, b);
+      } catch (mtl::mat::umfpack::error& e) {
+	AMDIS_ERROR_EXIT("UMFPACK_ERROR(solve, " << e.code << ") = " << e.what());
+      }
+      
+      auto r = Vector(b);
+      if (two_norm(x) != 0)
+        r -= A * x;
+      
+      solverInfo.setAbsResidual(two_norm(r));
+      solverInfo.setError(code);
+      
+      return code;
+    }
+
+    /// Implementation of \ref RunnerInterface::adjoint_solve()
+    virtual void exit() override {}
+
+  private:
+    FullMatrix fullMatrix;
+    std::shared_ptr<SolverType> solver;
+
+    int store_symbolic;
+    int symmetric_strategy;
+    double alloc_init;
+  };
+  
+}
+
+#endif // HAVE_UMFPACK

init/ellipt.json.3d

@@ -2,7 +2,7 @@
   "dimension of world":             3,
 
   "elliptMesh": {
-    "macro file name":    "./macro/macro.stand.3d",
+    "macro file name":    "./macro/BG1_withBoundary.3d",
     "global refinements": 0
   },
 
@@ -19,7 +19,7 @@
     },
 
     "output": {
-      "filename":         "output/ellipt.3d",
+      "filename":         "ellipt.3d",
       "output directory": "output"
     }
   }

init/heat.json.3d

+{
+  "dimension of world": 3,
+
+  "heatMesh": {
+    "macro file name":    "./macro/macro.stand.3d",
+    "global refinements": 0
+  },
+
+  "heat": {
+    "mesh":                   "heatMesh",
+    "names": "u",
+ 
+    "solver" : {
+      "name":               "cg",
+      "max iteration":      1000,
+      "absolute tolerance": 1e-6,
+      "info":               1,
+      "left precon":        "diag"
+    },
+
+    "output": {
+      "filename":         "heat.3d",
+      "output directory": "output"
+    }
+  },
+  
+  "adapt": {
+    "timestep":    0.1,
+    "start time":  0.0,
+    "end time":    0.2
+  }
+}

init/pfc.json.2d

@@ -3,37 +3,34 @@
 
   "pfcMesh": {
     "macro file name":    "./macro/macro.stand.2d",
-    "global refinements": 6,
+    "global refinements": 0,
     "min corner": "0,0",
     "max corner": "1,1",
-    "num cells": "2,2",
+    "num cells":  "2,2",
     
-    "dimension": "10,10"
+    "dimension":  "10,10"
   },
 
   "pfc": {
     "mesh":                   "pfcMesh",
  
     "solver" : {
-      "name":           "bicgstab_ell",
+      "name":           "direct",
       "max iteration":  1000,
       "tolerance":      1e-8,
       "info":           10,
-      "ell": 3,
-      "left precon":    "ilu"
+      "ell":            3
     },
 
     "output": {
-      "filename":       "pfc.2d",
-      "output directory":       "output",
-      "ParaView format": 1,
-      "ParaView mode": 1
+      "filename":         "pfc.2d",
+      "output directory": "output"
     }
   },
   
   "adapt": {
     "timestep":    0.1,
     "start time":  0.0,
-    "end time":    10.0
+    "end time":    0.1
   }
 }

src/CMakeLists.txt

@@ -2,19 +2,19 @@
 set(projects2d "heat" "pfc" "stokes" "navier_stokes")
 
 foreach(project ${projects2d})
-    add_executable(${project} ${project}.cc)
-    add_dune_alberta_flags(GRIDDIM 2 WORLDDIM 2 ${project})
-    target_link_dune_default_libraries(${project})
-    target_link_libraries(${project} "duneamdis")
-    target_compile_definitions(${project} PRIVATE AMDIS_DIM=2 AMDIS_DOW=2)
+    add_executable(${project}2d ${project}.cc)
+    add_dune_alberta_flags(GRIDDIM 2 WORLDDIM 2 ${project}2d)
+    target_link_dune_default_libraries(${project}2d)
+    target_link_libraries(${project}2d "duneamdis")
+    target_compile_definitions(${project}2d PRIVATE AMDIS_DIM=2 AMDIS_DOW=2)
 endforeach()
 
-set(projects3d "ellipt")
+set(projects3d "ellipt" "heat")
 
 foreach(project ${projects3d})
-    add_executable(${project} ${project}.cc)
-    add_dune_alberta_flags(GRIDDIM 3 WORLDDIM 3 ${project})
-    target_link_dune_default_libraries(${project})
-    target_link_libraries(${project} "duneamdis")
-    target_compile_definitions(${project} PRIVATE AMDIS_DIM=3 AMDIS_DOW=3)