linear_solvers.h

Go to the documentation of this file.

// ----------------------------------------------------------------------------
//
// FCST: Fuel Cell Simulation Toolbox
//
// Copyright (C) 2006-2013 by Energy Systems Design Laboratory, University of Alberta
//
// This software is distributed under the MIT License
// For more information, see the README file in /doc/LICENSE
//
// - Class: linear_solvers.h
// - Description: This namespace contains various linear solvers and preconditioners
// - Developers: Valentin N. Zingan, University of Alberta
//               Marc Secanell Gallart, University of Alberta
// - Id: $Id: linear_solvers.h 1354 2013-08-17 00:01:22Z secanell $
//
// ----------------------------------------------------------------------------

#ifndef _FCST_LINEAR_SOLVERS_H_
#define _FCST_LINEAR_SOLVERS_H_

#include <lac/vector.h>
#include <lac/block_vector.h>

#include <lac/sparse_matrix.h>
#include <lac/block_sparse_matrix.h>
#include <lac/filtered_matrix.h>

#include <lac/block_matrix_array.h>
#include <lac/vector_memory.h>

#include <lac/sparse_direct.h>
#include <lac/solver_gmres.h>

#include <lac/precondition.h>
#include <lac/sparse_ilu.h>

using namespace dealii;

namespace LinearSolvers
{

static GrowingVectorMemory< Vector<double> > vector_pool;
static GrowingVectorMemory< BlockVector<double> > block_vector_pool;

class SparseDirectUMFPACKSolver
{
public:



  SparseDirectUMFPACKSolver() { deallog << "Solving the linear system using UMFPACK" << std::endl; }

 ~SparseDirectUMFPACKSolver() { }




  template<typename MATRIX, typename VECTOR>
  void solve(const MATRIX& matrix,
             VECTOR&       solution,
             const VECTOR& right_hand_side) const
  {
         Vector<double> tmp;
         tmp = right_hand_side;

         SparseDirectUMFPACK solver;
         solver.initialize(matrix);
         solver.solve(tmp);

         solution = tmp;
  }


};

class GMRESSolver
{
public:



  GMRESSolver() { deallog << "Solving the linear system using GMRES" << std::endl; }

 ~GMRESSolver() { }




  template<typename MATRIX, typename VECTOR, typename PRECONDITIONER>
  void solve(const MATRIX&         matrix,
             VECTOR&               solution,
             const VECTOR&         right_hand_side,
             const unsigned int&   n_iter,
             const double&         tolerance,
             const PRECONDITIONER& preconditioner) const
  {
         SolverControl solver_control(n_iter,
                                      tolerance);

         typename SolverGMRES< VECTOR >::AdditionalData additional_data(100,
                                                                        false,
                                                                        true);

         SolverGMRES< VECTOR > solver(solver_control,
                                      additional_data);

         solver.solve(matrix,
                      solution,
                      right_hand_side,
                      preconditioner);
  }


};

class ILUPreconditioner
{
public:



  ILUPreconditioner(const BlockSparseMatrix<double>& matrix)
  :
  preconditioner_matrix(matrix.n_block_rows()),
  preconditioner_pointer(matrix.n_block_rows())
  {
         deallog << "ILUPreconditioner is used" << std::endl;

         preconditioner.initialize(matrix.n_block_rows(),
                                   vector_pool);

         for(unsigned int i = 0; i < matrix.n_block_rows(); ++i)
         {
                preconditioner_matrix[i].initialize(matrix.block(i,i),
                                                    SparseILU<double>::AdditionalData());
                preconditioner_pointer[i].set_memory(&vector_pool);
                preconditioner_pointer[i] = &preconditioner_matrix[i];
         }

         for(unsigned int i = 0; i < matrix.n_block_rows(); ++i)
                for(unsigned int j = 0; j < matrix.n_block_cols(); ++j)
                {
                       if( i == j )
                              preconditioner.enter(preconditioner_pointer[i], i, i);
                       else
                              preconditioner.enter(matrix.block(i,j), i, j);
                }
  }

 ~ILUPreconditioner() { }




  BlockTrianglePrecondition<double> preconditioner;

  std::vector< SparseILU<double> > preconditioner_matrix;

  std::vector< PointerMatrixAux< SparseILU<double>, Vector<double> > > preconditioner_pointer;


};

} // LinearSolvers

#endif