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
//
// ----------------------------------------------------------------------------

#ifndef _FCST_LINEAR_SOLVERS_H_
#define _FCST_LINEAR_SOLVERS_H_

//-- dealII
#include <deal.II/lac/vector.h>
#include <deal.II/lac/block_vector.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/block_sparse_matrix.h>
#include <deal.II/lac/block_matrix_array.h>
#include <deal.II/lac/vector_memory.h>
#include <deal.II/lac/sparse_direct.h>
#include <deal.II/lac/solver_gmres.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/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() { FcstUtilities::log << "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() { FcstUtilities::log << "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);
            
        }
        
        template<typename MATRIX, typename VECTOR, typename PRECONDITIONER>
        void solve(SolverControl& solver_control,
                   const MATRIX&         matrix,
                   VECTOR&               solution,
                   const VECTOR&         right_hand_side,
                   const PRECONDITIONER& preconditioner) const
        {
            SolverGMRES< VECTOR > solver(solver_control);
            
            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()),
        preconditioner(matrix.n_block_rows())
        {
            FcstUtilities::log << "ILUPreconditioner is used" << std::endl;
                        
            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