gsFlowLinSystSolver.h

Go to the documentation of this file.

 
#pragma once
 
#include <gsIncompressibleFlow/src/gsFlowUtils.h>
#include <gsIncompressibleFlow/src/gsFlowSolverParams.h>
#include <gsIncompressibleFlow/src/gsINSAssembler.h>
#include <gsSolver/gsGMRes.h>
 
namespace gismo
{
 
template<class T, int MatOrder>
class gsFlowLinSystSolver
{
 
protected: // *** Class members ***
 
    const gsFlowSolverParams<T>& m_paramsRef;
    gsStopwatch m_clock;
    T m_setupT, m_solveT;
 
 
public: // *** Constructor/destructor ***
 
    gsFlowLinSystSolver(const gsFlowSolverParams<T>& params):
    m_paramsRef(params)
    {
        m_setupT = 0;
        m_solveT = 0;
    }
 
    virtual ~gsFlowLinSystSolver()
    { }
 
 
protected: // *** Member functions ***
 
    T stopwatchStart();
 
    T stopwatchStop();
 
 
public: // *** Member functions ***
 
    virtual void setupSolver(const gsSparseMatrix<T, MatOrder>& mat)
    {GISMO_NO_IMPLEMENTATION}
 
    virtual void applySolver(const gsSparseMatrix<T, MatOrder>& mat, const gsMatrix<T>& rhs, gsMatrix<T>& solution)
    {GISMO_NO_IMPLEMENTATION}
 
    virtual void applySolver(const gsSparseMatrix<T, MatOrder>& mat, const gsMatrix<T>& rhs, gsMatrix<T>& solution, real_t alpha_u, real_t alpha_p, index_t usize, index_t pdofs);
 
 
public: // *** Getters ***
 
    virtual const T getSolverSetupTime() const { return m_setupT; }
 
    virtual const T getSolveTime() const { return m_solveT; }
 
 
}; // gsFlowLinSystSolver
 
// ===================================================================================================================
// ===================================================================================================================
 
template<class T, int MatOrder>
class gsFlowLinSystSolver_direct: public gsFlowLinSystSolver<T, MatOrder>
{
 
public:
    typedef gsFlowLinSystSolver<T, MatOrder> Base;
 
 
protected: // *** Class members ***
 
#ifdef GISMO_WITH_PARDISO
    typename gsSparseSolver<T>::PardisoLU m_solver;
#else
    typename gsSparseSolver<T>::LU m_solver;
#endif
 
 
protected: // *** Base class members ***
 
    using Base::m_paramsRef;
    using Base::m_setupT;
    using Base::m_solveT;
    using Base::stopwatchStart;
    using Base::stopwatchStop;
 
 
public: // *** Constructor/destructor ***
 
    gsFlowLinSystSolver_direct(const gsFlowSolverParams<T>& params):
    Base(params)
    {
        #ifdef GISMO_WITH_PARDISO
        pardisoSetup<T>(m_solver);
        #endif
    }
 
 
public: // *** Member functions ***
 
    virtual void setupSolver(const gsSparseMatrix<T, MatOrder>& mat);
 
    virtual void applySolver(const gsSparseMatrix<T, MatOrder>& mat, const gsMatrix<T>& rhs, gsMatrix<T>& solution);
 
 
}; // gsFlowLinSystSolver_direct
 
// ===================================================================================================================
 
template<class T, int MatOrder, class SolverType>
class gsFlowLinSystSolver_iter: public gsFlowLinSystSolver<T, MatOrder>
{
 
public:
    typedef gsFlowLinSystSolver<T, MatOrder> Base;
 
 
protected: // *** Class members ***
 
    typename gsLinearOperator<T>::Ptr m_precPtr;
    std::vector<index_t> m_linIterVector;
 
 
protected: // *** Base class members ***
 
    using Base::m_paramsRef;
    using Base::m_setupT;
    using Base::m_solveT;
    using Base::stopwatchStart;
    using Base::stopwatchStop;
 
 
public: // *** Constructor/destructor ***
 
    gsFlowLinSystSolver_iter(const gsFlowSolverParams<T>& params):
    Base(params)
    { }
 
 
public: // *** Member functions ***
 
    virtual void setupSolver(const gsSparseMatrix<T, MatOrder>& mat)
    { }
 
    virtual void setupPreconditioner(const gsSparseMatrix<T, MatOrder>& mat);
 
    virtual void applySolver(const gsSparseMatrix<T, MatOrder>& mat, const gsMatrix<T>& rhs, gsMatrix<T>& solution);
 
 
public: // *** Getters/setters ***
 
    std::vector<index_t> getLinIterVector() const { return m_linIterVector; }
 
    T getAvgLinIterations() const
    {
        index_t linIterSum = 0;
 
        for (size_t i = 0; i < m_linIterVector.size(); i++)
            linIterSum += m_linIterVector[i];
 
        return (T)linIterSum / m_linIterVector.size();
    }
 
}; // gsFlowLinSystSolver_iter
 
// ===================================================================================================================
 
template<class T, int MatOrder, class SolverType>
class gsFlowLinSystSolver_iterSP: public gsFlowLinSystSolver_iter<T, MatOrder, SolverType>
{
 
public:
    typedef gsFlowLinSystSolver_iter<T, MatOrder, SolverType> Base;
 
 
protected: // *** Class members ***
 
    std::string m_precType;
    gsOptionList m_precOpt;
    const gsINSAssembler<T, MatOrder>* m_assemblerPtr;
    std::map<std::string, gsSparseMatrix<T, MatOrder> > m_matrices;
    std::vector<index_t> m_linIterVector;
 
 
protected: // *** Base class members ***
 
    using Base::m_precPtr;
    using Base::m_paramsRef;
    using Base::m_setupT;
    using Base::m_solveT;
    using Base::stopwatchStart;
    using Base::stopwatchStop;
 
 
public: // *** Constructor/destructor ***
 
    gsFlowLinSystSolver_iterSP(const gsFlowSolverParams<T>& params, const gsINSAssembler<T, MatOrder>* assemblerPtr):
    Base(params), m_assemblerPtr(assemblerPtr)
    {
        m_precType = m_paramsRef.options().getString("lin.precType");
        m_precOpt = m_paramsRef.precOptions();
        m_precOpt.addInt("dim", "Problem dimension", m_assemblerPtr->getTarDim());
        m_precOpt.addReal("visc", "Viscosity", m_paramsRef.getPde().viscosity());
        m_precOpt.addInt("udofs", "Number of velocity dofs", m_assemblerPtr->getUdofs());
        m_precOpt.addInt("pdofs", "Number of pressure dofs", m_assemblerPtr->getPdofs());
    }
 
 
public: // *** Member functions ***
 
    virtual void setupPreconditioner(const gsSparseMatrix<T, MatOrder>& mat);
 
}; // gsFlowLinSystSolver_iterSP
 
// ===================================================================================================================
// ===================================================================================================================
 
template<class T, int MatOrder, class SolverType = gsGMRes<T> >
gsFlowLinSystSolver<T, MatOrder>* createLinSolver(const gsFlowSolverParams<T>& params, const gsFlowAssemblerBase<T, MatOrder>* assemblerPtr = NULL)
{
    const gsINSAssembler<T, MatOrder>* INSassembPtr = dynamic_cast< const gsINSAssembler<T, MatOrder>* >(assemblerPtr);
 
    std::string type = params.options().getString("lin.solver");
 
    if (type == "direct")
        return new gsFlowLinSystSolver_direct<T, MatOrder>(params);
    else if (type == "iter" && INSassembPtr == NULL) // assembler is not an INS assembler
        return new gsFlowLinSystSolver_iter<T, MatOrder, SolverType>(params);
    else if (type == "iter" && INSassembPtr != NULL) // assembler is an INS assembler
        return new gsFlowLinSystSolver_iterSP<T, MatOrder, SolverType>(params, INSassembPtr);
    // else if (type == _"petsc")
    //     return new gsFlowLinSystSolver_PETSc<T, MatOrder>(params);
    else
    {
        gsInfo << "Invalid linear system solver type, using direct.\n";
        return new gsFlowLinSystSolver_direct<T, MatOrder>(params);
    }
 
}
 
} // namespace gismo
 
#ifndef GISMO_BUILD_LIB
#include GISMO_HPP_HEADER(gsFlowLinSystSolver.hpp)
#endif