fwdsolver.h

// -*-C++-*-
// ==========================================================================
// Forward model: finite element method
// ==========================================================================

#ifndef __FWDSOLVER_H
#define __FWDSOLVER_H

//#include "supermatrix.h"
#include "toasttype.h"
#include "mathlib.h"

class Solution;
class MWsolution;


enum LSOLVER{
    LSOLVER_UNDEFINED,      
    LSOLVER_DIRECT,         
    LSOLVER_ITERATIVE       
};


enum DataScale {
    DATA_DEFAULT,           
    DATA_LIN,               
    DATA_LOG                
};

// =========================================================================
// Nonmember declarations

template<class T> class TFwdSolver;

template<class T>
STOASTLIB TVector<T> ProjectSingle (const QMMesh *mesh, int q,
    const TCompRowMatrix<T> &mvec, const TVector<T> &phi,
    DataScale dscale = DATA_LIN);

template<class T>
STOASTLIB TVector<T> ProjectAll (const QMMesh *mesh,
    const TCompRowMatrix<T> &mvec, const TVector<T> *phi,
    DataScale dscale = DATA_LIN);

// =========================================================================

template<class T> class TFwdSolver {
public:
    TFwdSolver (const QMMesh *mesh, LSOLVER linsolver, double tol = 1e-10);

    TFwdSolver (const QMMesh *mesh, const char *solver, double tol = 1e-10);

    TFwdSolver (const QMMesh *mesh, ParamParser &pp);

    ~TFwdSolver ();

    void SetDataScaling (DataScale scl);

    DataScale GetDataScaling () const;

    void SetPrecon (PreconType type);

    void SetPhaseUnwrap (bool unwrap)
    { unwrap_phase = unwrap; }

    void ReadParams (ParamParser &pp);

    void WriteParams (ParamParser &pp);

    void SetLinSolver (const char *solver, double tol = 1e-10);

    inline LSOLVER LinSolver() const { return solvertp; }

    inline double GetLinSolverTol() const { return iterative_tol; }

    inline void SetLinSolverMaxit (int maxit) { iterative_maxit = maxit; }

    inline const QMMesh *MeshPtr() const { return meshptr; }

    void Allocate ();

    void AssembleSystemMatrix (const Solution &sol, double omega = 0,
        bool elbasis = false);

    void AssembleSystemMatrixComponent (const RVector &prm, int type);

    void AssembleMassMatrix (const Mesh *mesh = 0);

    void Reset (const Solution &sol, double omega = 0, bool elbasis = false);

    void CalcField (const TVector<T> &qvec, TVector<T> &phi,
        IterativeSolverResult *res = 0) const;

    void CalcFields (const TCompRowMatrix<T> &qvec, TVector<T> *phi,
        IterativeSolverResult *res = 0) const;

    void CalcFields (int q0, int q1, const TCompRowMatrix<T> &qvec,
        TVector<T> *phi, IterativeSolverResult *res = 0) const;

    TVector<T> ProjectSingle (int q, const TCompRowMatrix<T> &mvec,
        const TVector<T> &phi, DataScale scl = DATA_DEFAULT) const;

    TVector<T> ProjectAll (const TCompRowMatrix<T> &mvec,
        const TVector<T> *dphi, DataScale scl = DATA_DEFAULT);

    TVector<T> ProjectAll (const TCompRowMatrix<T> &qvec,
        const TCompRowMatrix<T> &mvec, const Solution &sol, double omega,
        DataScale scl = DATA_DEFAULT);

    RVector ProjectAll_real (const TCompRowMatrix<T> &mvec,
        const TVector<T> *dphi, DataScale scl = DATA_DEFAULT);
    FVector ProjectAll_singlereal (const TCompRowMatrix<T> &mvec,
        const TVector<T> *dphi, DataScale scl = DATA_DEFAULT);

    RVector ProjectAll_real (const TCompRowMatrix<T> &qvec,
        const TCompRowMatrix<T> &mvec, const Solution &sol, double omega,
        DataScale scl = DATA_DEFAULT);
    FVector ProjectAll_singlereal (const TCompRowMatrix<T> &qvec,
        const TCompRowMatrix<T> &mvec, const Solution &sol, double omega,
        DataScale scl = DATA_DEFAULT);

    const QMMesh *meshptr;  
    LSOLVER solvertp;       
    IterativeMethod method; 
#ifdef MPI_FWDSOLVER
    TCompRowMatrixMPI<T> *F;  
#else
    TCompRowMatrix<T> *F;   
#endif
    TCompRowMatrix<T> *FL;  
    TVector<T> *Fd;         
    TPreconditioner<T> *precon; 
    TCompRowMatrix<T> *B;   
    //mutable SuperLU_data<T> lu_data; ///< parameters for LU solver
    void *SuperLU;          
    double iterative_tol;   
    int iterative_maxit;    
    PreconType precontp;    

protected:
    void Setup ();
    void SetupType ();
    void DeleteType ();

    RVector UnfoldComplex (const TVector<T> &vec) const;
    FVector UnfoldSComplex (const TVector<T> &vec) const;

#ifdef UNDEF
    void UnfoldComplex (const TVector<T> &vec, RVector &res) const;
    void UnfoldComplex (const TVector<T> &vec, FVector &res) const;
#endif

    DataScale dscale;       
    TVector<T> *pphi;       
    bool unwrap_phase;      

    // ===============================================================
    // MPI-specific functions and data members

#ifdef MPI_FWDSOLVER

public:

    void Setup_MPI();

    void Cleanup_MPI();

    void DistributeSources_MPI (int nq) const;

    void CalcFields_proc (const TCompRowMatrix<T> &qvec,
        TVector<T> *phi) const;

    void ProjectAll_proc (const TCompRowMatrix<T> &mvec,
        const TVector<T> *phi, DataScale scl,  TVector<T> &proj);

protected:

    MPI_Datatype mpitp;   
    int rnk;              
    int sze;              
    mutable int nQ;       
    mutable int *Q0, *Q1; 

#endif // MPI_FWDSOLVER

#if THREAD_LEVEL==2
    int nthread;
#endif
};


// OBSOLETE
STOASTLIB void Project_cplx (const QMMesh &mesh, int q, const CVector &phi,
    CVector &proj);

// ==========================================================================
// Some macros

#ifdef MPI_FWDSOLVER
#define SETUP_MPI() Setup_MPI()
#define CLEANUP_MPI() Cleanup_MPI()
#else
#define SETUP_MPI()
#define CLEANUP_MPI()
#endif

// ==========================================================================
// template typedefs

typedef TFwdSolver<float>          FFwdSolver;
typedef TFwdSolver<double>         RFwdSolver;
typedef TFwdSolver<std::complex<float> > SCFwdSolver;
typedef TFwdSolver<std::complex<double> > CFwdSolver;

// ==========================================================================
// extern declarations of FwdSolver (only required for VS)

#ifndef __FWDSOLVER_CC
extern template class STOASTLIB TFwdSolver<float>;
extern template class STOASTLIB TFwdSolver<double>;
extern template class STOASTLIB TFwdSolver<std::complex<float> >;
extern template class STOASTLIB TFwdSolver<std::complex<double> >;
#endif // !__FWDSOLVER_CC

#endif // !__FWDSOLVER_H