Nektar::SolverUtils::DriverArpack Class Reference

Base class for the development of solvers. More...

#include <DriverArpack.h>

Inheritance diagram for Nektar::SolverUtils::DriverArpack:

Collaboration diagram for Nektar::SolverUtils::DriverArpack:

Static Public Member Functions
static DriverSharedPtr	create (const LibUtilities::SessionReaderSharedPtr &pSession)
	Creates an instance of this class. More...

Static Public Attributes
static std::string	className = GetDriverFactory().RegisterCreatorFunction("Arpack", DriverArpack::create)
	Name of the class. More...

Protected Member Functions
	DriverArpack (const LibUtilities::SessionReaderSharedPtr pSession)
	Constructor. More...
virtual	~DriverArpack ()
	Destructor. More...
virtual void	v_InitObject (ostream &out=cout)
	Virtual function for initialisation implementation. More...
virtual void	v_Execute (ostream &out=cout)
	Virtual function for solve implementation. More...
Protected Member Functions inherited from Nektar::SolverUtils::DriverArnoldi
	DriverArnoldi (const LibUtilities::SessionReaderSharedPtr pSession)
	Constructor. More...
virtual	~DriverArnoldi ()
	Destructor. More...
void	CopyArnoldiArrayToField (Array< OneD, NekDouble > &array)
	Copy Arnoldi storage to fields. More...
void	CopyFieldToArnoldiArray (Array< OneD, NekDouble > &array)
	Copy fields to Arnoldi storage. More...
void	CopyFwdToAdj ()
	Copy the forward field to the adjoint system in transient growth calculations. More...
void	WriteFld (std::string file, std::vector< Array< OneD, NekDouble > > coeffs)
void	WriteFld (std::string file, Array< OneD, NekDouble > coeffs)
void	WriteEvs (ostream &evlout, const int k, const NekDouble real, const NekDouble imag, NekDouble resid=NekConstants::kNekUnsetDouble)
virtual Array< OneD, NekDouble >	v_GetRealEvl (void)
virtual Array< OneD, NekDouble >	v_GetImagEvl (void)
Protected Member Functions inherited from Nektar::SolverUtils::Driver
	Driver (const LibUtilities::SessionReaderSharedPtr pSession)
	Initialises EquationSystem class members. More...

Protected Attributes
int	m_maxn
int	m_maxnev
int	m_maxncv
Protected Attributes inherited from Nektar::SolverUtils::DriverArnoldi
int	m_kdim
int	m_nvec
	Dimension of Krylov subspace. More...
int	m_nits
	Number of vectors to test. More...
NekDouble	m_evtol
	Maxmum number of iterations. More...
NekDouble	m_period
	Tolerance of iteratiosn. More...
bool	m_timeSteppingAlgorithm
	Period of time stepping algorithm. More...
int	m_infosteps
	underlying operator is time stepping More...
int	m_nfields
	interval to dump information if required. More...
NekDouble	m_realShift
NekDouble	m_imagShift
Array< OneD, NekDouble >	m_real_evl
Array< OneD, NekDouble >	m_imag_evl
Protected Attributes inherited from Nektar::SolverUtils::Driver
LibUtilities::CommSharedPtr	m_comm
	Communication object. More...
LibUtilities::SessionReaderSharedPtr	m_session
	Session reader object. More...
LibUtilities::SessionReaderSharedPtr	session_LinNS
	I the Coupling between SFD and arnoldi. More...
Array< OneD, EquationSystemSharedPtr >	m_equ
	Equation system to solve. More...
int	m_nequ
	number of equations More...
enum EvolutionOperatorType	m_EvolutionOperator
	Evolution Operator. More...

Static Protected Attributes
static std::string	arpackProblemTypeLookupIds []
static std::string	arpackProblemTypeDefault = LibUtilities::SessionReader::RegisterDefaultSolverInfo("ArpackProblemType","LargestMag")
static std::string	driverLookupId = LibUtilities::SessionReader::RegisterEnumValue("Driver","Arpack",0)
Static Protected Attributes inherited from Nektar::SolverUtils::Driver
static std::string	evolutionOperatorLookupIds []
static std::string	evolutionOperatorDef
static std::string	driverDefault

Static Private Attributes
static std::string	ArpackProblemTypeTrans []

Friends
class	MemoryManager< DriverArpack >

Additional Inherited Members
Public Member Functions inherited from Nektar::SolverUtils::DriverArnoldi
SOLVER_UTILS_EXPORT void	ArnoldiSummary (std::ostream &out)
Public Member Functions inherited from Nektar::SolverUtils::Driver
virtual	~Driver ()
	Destructor. More...
SOLVER_UTILS_EXPORT void	InitObject (ostream &out=cout)
	Initialise Object. More...
SOLVER_UTILS_EXPORT void	Execute (ostream &out=cout)
	Execute driver. More...
SOLVER_UTILS_EXPORT Array< OneD, EquationSystemSharedPtr >	GetEqu ()
SOLVER_UTILS_EXPORT Array< OneD, NekDouble >	GetRealEvl (void)
SOLVER_UTILS_EXPORT Array< OneD, NekDouble >	GetImagEvl (void)

Detailed Description

Base class for the development of solvers.

Definition at line 48 of file DriverArpack.h.

Constructor & Destructor Documentation

Nektar::SolverUtils::DriverArpack::DriverArpack ( const LibUtilities::SessionReaderSharedPtr  pSession )

protected

Constructor.

Definition at line 62 of file DriverArpack.cpp.

            : DriverArnoldi(pSession)
        {
        }

Nektar::SolverUtils::DriverArpack::~DriverArpack ( )

protectedvirtual

Destructor.

Definition at line 71 of file DriverArpack.cpp.

        {
        }

Member Function Documentation

static DriverSharedPtr Nektar::SolverUtils::DriverArpack::create ( const LibUtilities::SessionReaderSharedPtr &  pSession )

inlinestatic

Creates an instance of this class.

Definition at line 54 of file DriverArpack.h.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr().

                                                                                              {
                DriverSharedPtr p = MemoryManager<DriverArpack>::AllocateSharedPtr(pSession);
                p->InitObject();
                return p;
            }

void Nektar::SolverUtils::DriverArpack::v_Execute ( ostream &  out = cout )

protectedvirtual

Virtual function for solve implementation.

Implements Nektar::SolverUtils::Driver.

Definition at line 109 of file DriverArpack.cpp.

References ArpackProblemTypeTrans, ASSERTL0, Nektar::SolverUtils::DriverArnoldi::CopyArnoldiArrayToField(), Nektar::SolverUtils::DriverArnoldi::CopyFieldToArnoldiArray(), Nektar::SolverUtils::DriverArnoldi::CopyFwdToAdj(), Arpack::Dnaupd(), Arpack::Dneupd(), Nektar::SolverUtils::eTransientGrowth, Nektar::NekConstants::kNekZeroTol, Nektar::SolverUtils::Driver::m_comm, Nektar::SolverUtils::Driver::m_equ, Nektar::SolverUtils::Driver::m_EvolutionOperator, Nektar::SolverUtils::DriverArnoldi::m_evtol, Nektar::SolverUtils::DriverArnoldi::m_imag_evl, Nektar::SolverUtils::DriverArnoldi::m_imagShift, Nektar::SolverUtils::DriverArnoldi::m_infosteps, Nektar::SolverUtils::DriverArnoldi::m_kdim, m_maxn, Nektar::SolverUtils::DriverArnoldi::m_nfields, Nektar::SolverUtils::DriverArnoldi::m_nits, Nektar::SolverUtils::DriverArnoldi::m_nvec, Nektar::SolverUtils::DriverArnoldi::m_real_evl, Nektar::SolverUtils::DriverArnoldi::m_realShift, Nektar::SolverUtils::Driver::m_session, Nektar::SolverUtils::DriverArnoldi::m_timeSteppingAlgorithm, Nektar::SolverUtils::DriverArnoldi::WriteEvs(), and Nektar::SolverUtils::DriverArnoldi::WriteFld().

        {
            Array<OneD, NekDouble> tmpworkd;

            int  nq     = m_equ[0]->UpdateFields()[0]->GetNcoeffs(); // Number of points in the mesh
            int  n      = m_nfields*nq;    // Number of points in eigenvalue calculation
            int lworkl  = 3*m_kdim*(m_kdim+2); // Size of work array
            int       ido ;     //REVERSE COMMUNICATION parameter. At the first call must be initialised at 0
            int       info;     // do not set initial vector (info=0 random initial vector, info=1 read initial vector from session file)
        
            int iparam[11];
            int ipntr[14];
        
            Array<OneD, int> ritzSelect;
            Array<OneD, NekDouble> dr;
            Array<OneD, NekDouble> di;
            Array<OneD, NekDouble> workev;
            Array<OneD, NekDouble> z;
            NekDouble sigmar, sigmai;

            Array<OneD, NekDouble> resid(n);
            Array<OneD, NekDouble> v(n*m_kdim);
            Array<OneD, NekDouble> workl(lworkl);
            Array<OneD, NekDouble> workd(3*n, 0.0);

            ASSERTL0(n <= m_maxn,  "N is greater than   MAXN");

            if(m_session->DefinesFunction("InitialConditions"))
            {
            
                out << "\tInital vector       : input file  " << endl;
                info = 1;
                CopyFieldToArnoldiArray(resid);
            
            }
            else
            {
                out << "\tInital vector       : random  " << endl;
                info = 0;
            }
        
    
            iparam[0] = 1;      // strategy for shift-invert
            iparam[1] = 0;      // (deprecated)
            iparam[2] = m_nits; // maximum number of iterations allowed/taken
            iparam[3] = 1;      // blocksize to be used for recurrence
            iparam[4] = 0;      // number of converged ritz eigenvalues
            iparam[5] = 0;      // (deprecated)
            if((fabs(m_realShift) > NekConstants::kNekZeroTol)|| // use shift if m_realShift > 1e-12
               (fabs(m_imagShift) > NekConstants::kNekZeroTol))        
            {
                iparam[6] = 3;
            }
            else
            {
                iparam[6] = 1;      // computation mode 1=> matrix-vector prod
            }
            iparam[7] = 0;      // (for shift-invert)
            iparam[8] = 0;      // number of MV operations
            iparam[9] = 0;      // number of BV operations
            iparam[10]= 0;      // number of reorthogonalisation steps

            int cycle = 0;
            const char* problem = ArpackProblemTypeTrans[m_session->GetSolverInfoAsEnum<int>("ArpackProblemType")].c_str();
       
            std::string name = m_session->GetSessionName() + ".evl"; 
            ofstream pFile(name.c_str());
        
            ido     = 0;    //At the first call must be initialisedat 0
        
            while(ido != 99)//ido==-1 || ido==1 || ido==0)
            {
                //Routine for eigenvalue evaluation for non-symmetric operators
                Arpack::Dnaupd( ido, "I",       // B='I' for std eval problem
                                n, problem,  m_nvec,
                                m_evtol, &resid[0], m_kdim, 
                                &v[0], n, iparam, ipntr, &workd[0],
                                &workl[0], lworkl, info);
            
                //Plotting of real and imaginary part of the eigenvalues from workl
                out << "\rIteration " << cycle << ", output: " << info << ", ido=" << ido << " " << std::flush;

                if(!((cycle-1)%m_kdim)&&(cycle> m_kdim))
                {
                    pFile << "Krylov spectrum at iteration: " <<  cycle << endl;

                    if(m_timeSteppingAlgorithm)
                    {
                        pFile << "EV  Magnitude   Angle       Growth      Frequency   Residual"    << endl;
                    }
                    else
                    {
                        pFile << "EV  Real        Imaginary   inverse real  inverse imag  Residual"   << endl;
                    }

                    out << endl;
                    for(int k=0; k<=m_kdim-1; ++k)
                    {                    
                        // write m_nvec eigs to screen
                        if(m_kdim-1-k < m_nvec)
                        {
                            WriteEvs(out,k, workl[ipntr[5]-1+k],workl[ipntr[6]-1+k]);
                        }
                        // write m_kdim eigs to screen
                        WriteEvs(pFile,k, workl[ipntr[5]-1+k],workl[ipntr[6]-1+k]);
                    }
                }
            
                cycle++;
            
                if (ido == 99) break;
                        
                ASSERTL0(ido == 1, "Unexpected reverse communication request.");
   
                //workd[inptr[0]-1] copied into operator fields
                CopyArnoldiArrayToField(tmpworkd = workd + (ipntr[0]-1));

            m_equ[0]->TransCoeffToPhys();

                m_equ[0]->DoSolve();

                if(!(cycle%m_infosteps))
                {
                    out << endl;
                    m_equ[0]->Output();
                }

                if(m_EvolutionOperator == eTransientGrowth)
                {
                    //start Adjoint with latest fields of direct 
                    CopyFwdToAdj();
        
                    m_equ[1]->TransCoeffToPhys();
                    m_equ[1]->DoSolve();
                }
            
                // operated fields are copied into workd[inptr[1]-1] 
                CopyFieldToArnoldiArray(tmpworkd = workd + (ipntr[1]-1));
            
            }
        
            out<< endl << "Converged in " << iparam[8] << " iterations" << endl;
    
            ASSERTL0(info >= 0," Error with Dnaupd");
    
            ritzSelect = Array<OneD, int>       (m_kdim,0);
            dr         = Array<OneD, NekDouble> (m_nvec+1,0.0);
            di         = Array<OneD, NekDouble> (m_nvec+1,0.0);
            workev     = Array<OneD, NekDouble> (3*m_kdim);
            z          = Array<OneD, NekDouble> (n*(m_nvec+1));
        
            sigmar     = m_realShift; 
            sigmai     = m_imagShift;
    
            //Setting 'A', Ritz vectors are computed. 'S' for Shur vectors
            Arpack::Dneupd(1, "A", ritzSelect.get(), dr.get(), di.get(), z.get(), n, sigmar, sigmai, workev.get(), "I", n, problem, m_nvec, m_evtol, resid.get(), m_kdim, v.get(), n, iparam, ipntr, workd.get(), workl.get(),lworkl,info);
        
            ASSERTL0(info == 0, " Error with Dneupd");
            int nconv=iparam[4];
            Array<OneD, MultiRegions::ExpListSharedPtr>  fields = m_equ[0]->UpdateFields();
        
            out << "Converged Eigenvalues: " << nconv << endl;
            pFile << "Converged Eigenvalues:"<< nconv << endl;

            if(m_timeSteppingAlgorithm)
            {
                pFile << "EV  Magnitude   Angle       Growth      Frequency" << endl;
            }
            else
            {
                pFile << "EV  Real        Imaginary   inverse real  inverse imag" << endl;
            }

        
            for(int i= 0; i< nconv; ++i)
            {
                WriteEvs(out,i,dr[i],di[i]);
                WriteEvs(pFile,i,dr[i],di[i]);
        
                std::string file = m_session->GetSessionName() + "_eig_"
                                        + boost::lexical_cast<std::string>(i);
                WriteFld(file,z + i*nq);
            }
        
            m_real_evl = dr;
            m_imag_evl = di;
       
            pFile.close();
        
            for(int j = 0; j < m_equ[0]->GetNvariables(); ++j)
            {
                NekDouble vL2Error = m_equ[0]->L2Error(j,false);
                NekDouble vLinfError = m_equ[0]->LinfError(j);
                if (m_comm->GetRank() == 0)
                {
                    out << "L 2 error (variable " << m_equ[0]->GetVariable(j) << ") : " << vL2Error << endl;
                    out << "L inf error (variable " << m_equ[0]->GetVariable(j) << ") : " << vLinfError << endl;
                }   
            }
        }

void Nektar::SolverUtils::DriverArpack::v_InitObject ( ostream &  out = cout )

protectedvirtual

Virtual function for initialisation implementation.

Reimplemented from Nektar::SolverUtils::DriverArnoldi.

Definition at line 78 of file DriverArpack.cpp.

References Nektar::SolverUtils::DriverArnoldi::ArnoldiSummary(), ArpackProblemTypeTrans, ASSERTL0, Nektar::SolverUtils::Driver::m_equ, Nektar::SolverUtils::DriverArnoldi::m_kdim, m_maxn, m_maxncv, m_maxnev, Nektar::SolverUtils::Driver::m_nequ, Nektar::SolverUtils::DriverArnoldi::m_nvec, Nektar::SolverUtils::Driver::m_session, and Nektar::SolverUtils::DriverArnoldi::v_InitObject().

        {
            DriverArnoldi::v_InitObject(out);
        
            //Initialisation of Arnoldi parameters
            m_maxn   = 1000000; // Maximum size of the problem
            m_maxnev = 200;      // maximum number of eigenvalues requested
            m_maxncv = 500;     // Largest number of basis vector used in Implicitly Restarted Arnoldi      
                    
            // Error alerts
            ASSERTL0(m_nvec <= m_maxnev,"NEV is greater than MAXNEV");
            ASSERTL0(m_kdim <= m_maxncv,"NEV is greater than MAXNEV");
            ASSERTL0(2      <= m_kdim-m_nvec,"NCV-NEV is less than 2");
    
            m_equ[0]->PrintSummary(out);
        
            // Print session parameters
            out << "\tArnoldi solver type    : Arpack" << endl;

            out << "\tArpack problem type    : ";
            out << ArpackProblemTypeTrans[m_session->GetSolverInfoAsEnum<int>("ArpackProblemType")] << endl;
            DriverArnoldi::ArnoldiSummary(out);
        
            m_equ[m_nequ - 1]->DoInitialise();
        
            //FwdTrans Initial conditions to be in Coefficient Space
            m_equ[m_nequ-1] ->TransPhysToCoeff();

        }

Friends And Related Function Documentation

friend class MemoryManager< DriverArpack >

friend

Definition at line 51 of file DriverArpack.h.

Member Data Documentation

std::string Nektar::SolverUtils::DriverArpack::arpackProblemTypeDefault = LibUtilities::SessionReader::RegisterDefaultSolverInfo("ArpackProblemType","LargestMag")

staticprotected

Definition at line 85 of file DriverArpack.h.

std::string Nektar::SolverUtils::DriverArpack::arpackProblemTypeLookupIds

staticprotected

Initial value:

= {
            LibUtilities::SessionReader::RegisterEnumValue("ArpackProblemType","LargestReal"    ,0),
            LibUtilities::SessionReader::RegisterEnumValue("ArpackProblemType","SmallestReal"   ,1),
            LibUtilities::SessionReader::RegisterEnumValue("ArpackProblemType","LargestImag"    ,2),
            LibUtilities::SessionReader::RegisterEnumValue("ArpackProblemType","SmallestImag"   ,3),
            LibUtilities::SessionReader::RegisterEnumValue("ArpackProblemType","LargestMag"     ,4),
            LibUtilities::SessionReader::RegisterEnumValue("ArpackProblemType","SmallestMag"    ,5),
        }

Definition at line 84 of file DriverArpack.h.

std::string Nektar::SolverUtils::DriverArpack::ArpackProblemTypeTrans

staticprivate

Initial value:

=
        { "LR", "SR", "LI", "SI", "LM", "SM" }

Definition at line 89 of file DriverArpack.h.

Referenced by v_Execute(), and v_InitObject().

std::string Nektar::SolverUtils::DriverArpack::className = GetDriverFactory().RegisterCreatorFunction("Arpack", DriverArpack::create)

static

Name of the class.

Definition at line 61 of file DriverArpack.h.

std::string Nektar::SolverUtils::DriverArpack::driverLookupId = LibUtilities::SessionReader::RegisterEnumValue("Driver","Arpack",0)

staticprotected

Definition at line 86 of file DriverArpack.h.

int Nektar::SolverUtils::DriverArpack::m_maxn

protected

Definition at line 66 of file DriverArpack.h.

Referenced by v_Execute(), and v_InitObject().

int Nektar::SolverUtils::DriverArpack::m_maxncv

protected

Definition at line 68 of file DriverArpack.h.

Referenced by v_InitObject().

int Nektar::SolverUtils::DriverArpack::m_maxnev

protected

Definition at line 67 of file DriverArpack.h.

Referenced by v_InitObject().