45 std::string DriverArpack::arpackProblemTypeLookupIds[6] = {
46 LibUtilities::SessionReader::RegisterEnumValue(
"ArpackProblemType",
"LargestReal" ,0),
47 LibUtilities::SessionReader::RegisterEnumValue(
"ArpackProblemType",
"SmallestReal" ,1),
48 LibUtilities::SessionReader::RegisterEnumValue(
"ArpackProblemType",
"LargestImag" ,2),
49 LibUtilities::SessionReader::RegisterEnumValue(
"ArpackProblemType",
"SmallestImag" ,3),
50 LibUtilities::SessionReader::RegisterEnumValue(
"ArpackProblemType",
"LargestMag" ,4),
51 LibUtilities::SessionReader::RegisterEnumValue(
"ArpackProblemType",
"SmallestMag" ,5),
53 std::string DriverArpack::arpackProblemTypeDefault = LibUtilities::SessionReader::RegisterDefaultSolverInfo(
"ArpackProblemType",
"LargestMag");
54 std::string DriverArpack::driverLookupId = LibUtilities::SessionReader::RegisterEnumValue(
"Driver",
"Arpack",0);
58 std::string DriverArpack::ArpackProblemTypeTrans[6] =
59 {
"LR",
"SR",
"LI",
"SI",
"LM",
"SM" };
95 m_equ[0]->PrintSummary(out);
98 out <<
"\tArnoldi solver type : Arpack" << endl;
100 out <<
"\tArpack problem type : ";
104 for(
int i = 0; i <
m_nequ; ++i)
106 m_equ[i]->DoInitialise();
110 m_equ[m_nequ-1] ->TransPhysToCoeff();
120 int nq =
m_equ[0]->UpdateFields()[0]->GetNcoeffs();
143 if(
m_session->DefinesFunction(
"InitialConditions"))
145 out <<
"\tInital vector : input file " << endl;
151 out <<
"\tInital vector : random " << endl;
195 std::string name =
m_session->GetSessionName() +
".evl";
196 ofstream pFile(name.c_str());
206 &v[0], n, iparam, ipntr, &workd[0],
207 &workl[0], lworkl, info);
211 out <<
"\rIteration " << cycle <<
", output: " << info <<
", ido=" << ido <<
" " << std::flush;
213 if(!((cycle-1)%m_kdim)&&(cycle> m_kdim)&&(ido!=2))
215 pFile <<
"Krylov spectrum at iteration: " << cycle << endl;
219 pFile <<
"EV Magnitude Angle Growth Frequency Residual" << endl;
223 pFile <<
"EV Real Imaginary inverse real inverse imag Residual" << endl;
227 for(
int k = 0; k <
m_kdim; ++k)
230 WriteEvs(pFile,k, workl[ipntr[5]-1+k],workl[ipntr[6]-1+k]);
234 if (ido == 99)
break;
245 m_equ[0]->TransCoeffToPhys();
253 m_equ[1]->TransCoeffToPhys();
273 m_equ[0]->TransCoeffToPhys();
276 for (
int i = 0; i < fields.num_elements(); ++i)
278 fields[i]->IProductWRTBase(fields[i]->GetPhys(),
279 fields[i]->UpdateCoeffs());
287 ASSERTL0(
false,
"Unexpected reverse communication request.");
292 out<< endl <<
"Converged in " << iparam[8] <<
" iterations" << endl;
294 ASSERTL0(info >= 0,
" Error with Dnaupd");
318 z.get(), n, sigmar, sigmai, workev.get(), &B,
320 v.get(), n, iparam, ipntr, workd.get(),
321 workl.get(),lworkl,info);
323 ASSERTL0(info == 0,
" Error with Dneupd");
338 "Need to implement Ritz re-evaluation of"
339 "eigenvalue. Only one half of complex "
340 "value will be correct");
345 out <<
"Converged Eigenvalues: " << nconv << endl;
346 pFile <<
"Converged Eigenvalues: " << nconv << endl;
350 out <<
" Magnitude Angle Growth Frequency" << endl;
351 pFile <<
" Magnitude Angle Growth Frequency" << endl;
352 for(
int i= 0; i< nconv; ++i)
357 std::string file =
m_session->GetSessionName() +
"_eig_"
358 + boost::lexical_cast<std::string>(i)
365 out <<
" Real Imaginary " << endl;
366 pFile <<
" Real Imaginary " << endl;
367 for(
int i= 0; i< nconv; ++i)
374 std::string file =
m_session->GetSessionName() +
"_eig_"
375 + boost::lexical_cast<std::string>(i)
386 for(
int j = 0; j <
m_equ[0]->GetNvariables(); ++j)
390 if (
m_comm->GetRank() == 0)
392 out <<
"L 2 error (variable " <<
m_equ[0]->GetVariable(j) <<
") : " << vL2Error << endl;
393 out <<
"L inf error (variable " <<
m_equ[0]->GetVariable(j) <<
") : " << vLinfError << endl;
Array< OneD, NekDouble > m_imag_evl
static void Dneupd(const int &rvec, const char *howmny, const int *select, double *dr, double *di, double *z, const int &ldz, const double &sigmar, const double &sigmai, double *workev, const char *bmat, const int &n, const char *which, const int &nev, const double &tol, double *resid, const int &ncv, double *v, const int &ldv, int *iparam, int *ipntr, double *workd, double *workl, const int &lworkl, int &info)
Post-processing routine to computed eigenvector of computed eigenvalues in Dnaupd.
virtual ~DriverArpack()
Destructor.
#define ASSERTL0(condition, msg)
void CopyFwdToAdj()
Copy the forward field to the adjoint system in transient growth calculations.
NekDouble m_evtol
Maxmum number of iterations.
void CopyArnoldiArrayToField(Array< OneD, NekDouble > &array)
Copy Arnoldi storage to fields.
virtual void v_Execute(std::ostream &out=std::cout)
Virtual function for solve implementation.
SOLVER_UTILS_EXPORT void ArnoldiSummary(std::ostream &out)
void CopyFieldToArnoldiArray(Array< OneD, NekDouble > &array)
Copy fields to Arnoldi storage.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
void WriteFld(std::string file, std::vector< Array< OneD, NekDouble > > coeffs)
Write coefficients to file.
virtual void v_InitObject(std::ostream &out=std::cout)
Virtual function for initialisation implementation.
static const NekDouble kNekZeroTol
int m_nvec
Dimension of Krylov subspace.
Array< OneD, NekDouble > m_real_evl
Operator in solve call is negated.
void WriteEvs(std::ostream &evlout, const int k, const NekDouble real, const NekDouble imag, NekDouble resid=NekConstants::kNekUnsetDouble, bool DumpInverse=true)
#define WARNINGL0(condition, msg)
virtual void v_InitObject(std::ostream &out=std::cout)
static void Dnaupd(int &ido, const char *bmat, const int &n, const char *which, const int &nev, const double &tol, double *resid, const int &ncv, double *v, const int &ldv, int *iparam, int *ipntr, double *workd, double *workl, const int &lworkl, int &info)
Top level reverse communication interface to solve real double-precision non-symmetric problems...
enum EvolutionOperatorType m_EvolutionOperator
Evolution Operator.
static const NekDouble kNekUnsetDouble
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
int m_nfields
interval to dump information if required.
bool m_timeSteppingAlgorithm
Period of time stepping algorithm.
static std::string ArpackProblemTypeTrans[]
Array< OneD, EquationSystemSharedPtr > m_equ
Equation system to solve.
int m_infosteps
underlying operator is time stepping
DriverFactory & GetDriverFactory()
LibUtilities::CommSharedPtr m_comm
Communication object.
Base class for the development of solvers.
LibUtilities::SessionReaderSharedPtr m_session
Session reader object.
int m_nequ
number of equations
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
int m_nits
Number of vectors to test.