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

#include <DriverArpack.h>

Inheritance diagram for Nektar::SolverUtils::DriverArpack:

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

Static Public Attributes
static std::string	className
	Name of the class. More...

Protected Member Functions
	DriverArpack (const LibUtilities::SessionReaderSharedPtr pSession, const SpatialDomains::MeshGraphSharedPtr pGraph)
	Constructor. More...

virtual	~DriverArpack ()
	Destructor. More...

virtual void	v_InitObject (std::ostream &out=std::cout) override
	Virtual function for initialisation implementation. More...

virtual void	v_Execute (std::ostream &out=std::cout) override
	Virtual function for solve implementation. More...

Protected Member Functions inherited from Nektar::SolverUtils::DriverArnoldi
	DriverArnoldi (const LibUtilities::SessionReaderSharedPtr pSession, const SpatialDomains::MeshGraphSharedPtr pGraph)
	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)
	Write coefficients to file. More...

void	WriteFld (std::string file, Array< OneD, NekDouble > coeffs)

void	WriteEvs (std::ostream &evlout, const int k, const NekDouble real, const NekDouble imag, NekDouble resid=NekConstants::kNekUnsetDouble, bool DumpInverse=true)

void	MaskInit ()
	init mask More...

void	GetUnmaskFunction (std::vector< std::vector< LibUtilities::EquationSharedPtr >> &unmaskfun)

virtual Array< OneD, NekDouble >	v_GetRealEvl (void) override

virtual Array< OneD, NekDouble >	v_GetImagEvl (void) override

Protected Member Functions inherited from Nektar::SolverUtils::Driver
	Driver (const LibUtilities::SessionReaderSharedPtr pSession, const SpatialDomains::MeshGraphSharedPtr pGraph)
	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

int	m_negatedOp

Array< OneD, NekDouble >	m_real_evl
	Operator in solve call is negated. More...

Array< OneD, NekDouble >	m_imag_evl

bool	m_useMask

Array< OneD, NekDouble >	m_maskCoeffs

Array< OneD, NekDouble >	m_maskPhys

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

SpatialDomains::MeshGraphSharedPtr	m_graph
	MeshGraph object. 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

static std::string	driverLookupId

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)

SOLVER_UTILS_EXPORT const Array< OneD, const NekDouble > &	GetMaskCoeff () const

SOLVER_UTILS_EXPORT const Array< OneD, const NekDouble > &	GetMaskPhys () const

Public Member Functions inherited from Nektar::SolverUtils::Driver
virtual	~Driver ()
	Destructor. More...

SOLVER_UTILS_EXPORT void	InitObject (std::ostream &out=std::cout)
	Initialise Object. More...

SOLVER_UTILS_EXPORT void	Execute (std::ostream &out=std::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

◆ DriverArpack()

Nektar::SolverUtils::DriverArpack::DriverArpack	(	const LibUtilities::SessionReaderSharedPtr	pSession,
		const SpatialDomains::MeshGraphSharedPtr	pGraph
	)

protected

Constructor.

Definition at line 73 of file DriverArpack.cpp.

     : DriverArnoldi(pSession, pGraph)
 {
 }

◆ ~DriverArpack()

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

protectedvirtual

Destructor.

Definition at line 82 of file DriverArpack.cpp.

83 {

84 }

Member Function Documentation

◆ create()

static DriverSharedPtr Nektar::SolverUtils::DriverArpack::create	(	const LibUtilities::SessionReaderSharedPtr &	pSession,
		const SpatialDomains::MeshGraphSharedPtr &	pGraph
	)

inlinestatic

Creates an instance of this class.

Definition at line 54 of file DriverArpack.h.

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

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), and CellMLToNektar.cellml_metadata::p.

◆ v_Execute()

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

overrideprotectedvirtual

Virtual function for solve implementation.

Implements Nektar::SolverUtils::Driver.

Definition at line 129 of file DriverArpack.cpp.

 {
     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, 0.0);
     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;
     }
  
     char B;
  
     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)
  
     // Use generalized B matrix for coupled solver.
     if (m_timeSteppingAlgorithm)
     {
         iparam[6] = 1; // computation mode 1=> matrix-vector prod
         B         = 'I';
     }
     else
     {
         iparam[6] = 3; // computation mode 1=> matrix-vector prod
         B         = 'G';
     }
 #if 0
     if((fabs(m_realShift) > NekConstants::kNekZeroTol)|| // use shift if m_realShift > 1e-12
        (fabs(m_imagShift) > NekConstants::kNekZeroTol))
     {
         iparam[6] = 3; // This was 3 need to know what to set it to
         B = 'G';
     }
     else
     {
         iparam[6] = 1;      // computation mode 1=> matrix-vector prod
         B = 'I';
     }
 #endif
     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, &B, // 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) && (ido != 2))
         {
             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; ++k)
             {
                 // write m_kdim eigs to screen
                 WriteEvs(pFile, k, workl[ipntr[5] - 1 + k],
                          workl[ipntr[6] - 1 + k]);
             }
         }
  
         if (ido == 99)
             break;
  
         switch (ido)
         {
             case -1:
             case 1: // Note that ido=1 we are using input x
                     // (workd[inptr[0]-1]) rather than Mx as
                     // recommended in manual since it is not
                     // possible to impose forcing directly.
                 CopyArnoldiArrayToField(tmpworkd = workd + (ipntr[0] - 1));
  
                 m_equ[0]->TransCoeffToPhys();
  
                 m_equ[0]->DoSolve();
                 if (m_EvolutionOperator == eTransientGrowth)
                 {
                     // start Adjoint with latest fields of direct
                     CopyFwdToAdj();
  
                     m_equ[1]->TransCoeffToPhys();
                     m_equ[1]->DoSolve();
                 }
  
                 if (!(cycle % m_infosteps))
                 {
                     out << endl;
                     m_equ[0]->Output();
                 }
  
                 // operated fields are copied into workd[inptr[1]-1]
                 CopyFieldToArnoldiArray(tmpworkd = workd + (ipntr[1] - 1));
  
                 cycle++;
                 break;
             case 2: // provide y = M x (bwd trans and iproduct);
             {
                 // workd[inptr[0]-1] copied into operator fields
                 CopyArnoldiArrayToField(tmpworkd = workd + (ipntr[0] - 1));
  
                 m_equ[0]->TransCoeffToPhys();
  
                 Array<OneD, MultiRegions::ExpListSharedPtr> fields =
                     m_equ[0]->UpdateFields();
                 for (int i = 0; i < fields.size(); ++i)
                 {
                     fields[i]->IProductWRTBase(fields[i]->GetPhys(),
                                                fields[i]->UpdateCoeffs());
                 }
  
                 // operated fields are copied into workd[inptr[1]-1]
                 CopyFieldToArnoldiArray(tmpworkd = workd + (ipntr[1] - 1));
                 break;
             }
             default:
                 ASSERTL0(false, "Unexpected reverse communication request.");
         }
     }
  
     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));
  
     if (m_negatedOp)
     {
         sigmar = -m_realShift;
     }
     else
     {
         sigmar = m_realShift;
     }
  
     // Do not pass imaginary shift to Arpack since we have not
     // used a Fortran complex number format and so processing
     // is mucked up. Need to do some processing afterwards.
     sigmai = 0;
  
     // 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(), &B, 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];
  
     // Subtract off complex shift if it exists
     if (m_negatedOp)
     {
         Vmath::Sadd(nconv, m_imagShift, di, 1, di, 1);
     }
     else
     {
         Vmath::Sadd(nconv, -m_imagShift, di, 1, di, 1);
     }
  
     WARNINGL0(m_imagShift == 0, "Complex Shift applied. "
                                 "Need to implement Ritz re-evaluation of"
                                 "eigenvalue. Only one half of complex "
                                 "value will be correct");
  
     Array<OneD, MultiRegions::ExpListSharedPtr> fields =
         m_equ[0]->UpdateFields();
  
     out << "Converged Eigenvalues: " << nconv << endl;
     pFile << "Converged Eigenvalues: " << nconv << endl;
  
     if (m_timeSteppingAlgorithm)
     {
         out << "         Magnitude   Angle       Growth      Frequency" << endl;
         pFile << "         Magnitude   Angle       Growth      Frequency"
               << 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) + ".fld";
             WriteFld(file, z + i * n);
         }
     }
     else
     {
         out << "         Real        Imaginary " << endl;
         pFile << "         Real        Imaginary " << endl;
         for (int i = 0; i < nconv; ++i)
         {
             WriteEvs(out, i, dr[i], di[i], NekConstants::kNekUnsetDouble,
                      false);
             WriteEvs(pFile, i, dr[i], di[i], NekConstants::kNekUnsetDouble,
                      false);
  
             std::string file = m_session->GetSessionName() + "_eig_" +
                                boost::lexical_cast<std::string>(i) + ".fld";
             WriteFld(file, z + i * n);
         }
     }
  
     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;
         }
     }
 }

◆ v_InitObject()

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

overrideprotectedvirtual

Virtual function for initialisation implementation.

Reimplemented from Nektar::SolverUtils::DriverArnoldi.

Definition at line 89 of file DriverArpack.cpp.

 {
     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");
  
     ASSERTL0(m_comm->GetSize() == 1,
              "..ARPACK is not currently set-up for parallel execution...\n");
     m_equ[0]->PrintSummary(out);
  
     ASSERTL0(m_comm->GetSize() == 1,
              "ARPACK Arnoldi solver does not support execution in parallel.");
  
     // Print session parameters
     out << "\tArnoldi solver type    : Arpack" << endl;
  
     out << "\tArpack problem type    : ";
     out << ArpackProblemTypeTrans[m_session->GetSolverInfoAsEnum<int>(
                "ArpackProblemType")]
         << endl;
     DriverArnoldi::ArnoldiSummary(out);
  
     for (int i = 0; i < m_nequ; ++i)
     {
         m_equ[i]->DoInitialise();
     }
  
     // FwdTrans Initial conditions to be in Coefficient Space
     m_equ[m_nequ - 1]->TransPhysToCoeff();
 }

References Nektar::SolverUtils::DriverArnoldi::ArnoldiSummary(), ArpackProblemTypeTrans, ASSERTL0, Nektar::SolverUtils::Driver::m_comm, 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().

Friends And Related Function Documentation

◆ MemoryManager< DriverArpack >

friend class MemoryManager< DriverArpack >

friend

Definition at line 1 of file DriverArpack.h.

Member Data Documentation

◆ arpackProblemTypeDefault

std::string Nektar::SolverUtils::DriverArpack::arpackProblemTypeDefault

staticprotected

Initial value:

=
    LibUtilities::SessionReader::RegisterDefaultSolverInfo("ArpackProblemType",
                                                           "LargestMag")

Definition at line 87 of file DriverArpack.h.

◆ arpackProblemTypeLookupIds

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 86 of file DriverArpack.h.

◆ ArpackProblemTypeTrans

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

staticprivate

Initial value:

= {"LR", "SR", "LI",

"SI", "LM", "SM"}

Definition at line 91 of file DriverArpack.h.

Referenced by v_Execute(), and v_InitObject().

◆ className

std::string Nektar::SolverUtils::DriverArpack::className

static

Initial value:

=

GetDriverFactory().RegisterCreatorFunction("Arpack", DriverArpack::create)

Nektar::LibUtilities::NekFactory::RegisterCreatorFunction

tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")

Register a class with the factory.

Definition: NekFactory.hpp:198

Nektar::SolverUtils::DriverArpack::create

static DriverSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)

Creates an instance of this class.

Definition: DriverArpack.h:54

Nektar::SolverUtils::GetDriverFactory

DriverFactory & GetDriverFactory()

Definition: Driver.cpp:64

Name of the class.

Definition at line 65 of file DriverArpack.h.

◆ driverLookupId

std::string Nektar::SolverUtils::DriverArpack::driverLookupId

staticprotected

Initial value:

=

LibUtilities::SessionReader::RegisterEnumValue("Driver", "Arpack", 0)

Definition at line 88 of file DriverArpack.h.

◆ m_maxn

int Nektar::SolverUtils::DriverArpack::m_maxn

protected

Definition at line 68 of file DriverArpack.h.

Referenced by v_Execute(), and v_InitObject().

◆ m_maxncv

int Nektar::SolverUtils::DriverArpack::m_maxncv

protected

Definition at line 70 of file DriverArpack.h.

Referenced by v_InitObject().

◆ m_maxnev

int Nektar::SolverUtils::DriverArpack::m_maxnev

protected

Definition at line 69 of file DriverArpack.h.

Referenced by v_InitObject().

Static Public Member Functions

Static Public Attributes

Protected Member Functions

Protected Attributes

Static Protected Attributes

Static Private Attributes

Friends

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ DriverArpack()

◆ ~DriverArpack()

Member Function Documentation

◆ create()

◆ v_Execute()

◆ v_InitObject()

Friends And Related Function Documentation

◆ MemoryManager< DriverArpack >

Member Data Documentation

◆ arpackProblemTypeDefault

◆ arpackProblemTypeLookupIds

◆ ArpackProblemTypeTrans

◆ className

◆ driverLookupId

◆ m_maxn

◆ m_maxncv

◆ m_maxnev