48 string DriverModifiedArnoldi::className =
50 DriverModifiedArnoldi::create);
51 string DriverModifiedArnoldi::driverLookupId =
52 LibUtilities::SessionReader::RegisterEnumValue(
"Driver",
58 DriverModifiedArnoldi::DriverModifiedArnoldi(
80 m_equ[0]->PrintSummary(out);
83 if (
m_comm->GetRank() == 0)
85 out <<
"\tArnoldi solver type : Modified Arnoldi" << endl;
90 for(
int i = 0; i <
m_nequ; ++i)
92 m_equ[i]->DoInitialise();
96 m_equ[m_nequ-1] ->TransPhysToCoeff();
108 int nq =
m_equ[0]->UpdateFields()[0]->GetNcoeffs();
113 std::string evlFile =
m_session->GetSessionName() +
".evl";
115 if (
m_comm->GetRank() == 0)
117 evlout.open(evlFile.c_str());
130 for (i = 0; i < m_kdim + 1; ++i)
137 if(
m_session->DefinesFunction(
"InitialConditions"))
139 if (
m_comm->GetRank() == 0)
141 out <<
"\tInital vector : specified in input file " << endl;
143 m_equ[0]->SetInitialConditions(0.0,
false);
149 if (
m_comm->GetRank() == 0)
151 out <<
"\tInital vector : random " << endl;
161 if (
m_comm->GetRank() == 0)
163 out <<
"Iteration: " << 0 << endl;
167 alpha[0] =
Blas::Ddot(ntot, &Kseq[0][0], 1, &Kseq[0][0], 1);
169 alpha[0] = std::sqrt(alpha[0]);
170 Vmath::Smul(ntot, 1.0/alpha[0], Kseq[0], 1, Kseq[0], 1);
174 for (i = 1; !converged && i <=
m_kdim; ++i)
180 alpha[i] =
Blas::Ddot(ntot, &Kseq[i][0], 1, &Kseq[i][0], 1);
182 alpha[i] = std::sqrt(alpha[i]);
185 Vmath::Smul(ntot, 1.0/alpha[i], Kseq[i], 1, Kseq[i], 1);
188 for (
int k = 0; k < i + 1; ++k)
194 EV_small(Tseq, ntot, alpha, i, zvec, wr, wi, resnorm);
197 converged =
EV_test(i, i, zvec, wr, wi, resnorm,
198 std::min(i,
m_nvec), evlout, resid0);
201 converged = max (converged, 0);
208 if (
m_comm->GetRank() == 0)
210 out <<
"Iteration: " << i <<
" (residual : " << resid0
218 for (i = m_kdim + 1; !converged && i <=
m_nits; ++i)
222 for (
int j = 1; j <=
m_kdim; ++j)
224 alpha[j-1] = alpha[j];
229 EV_update(Kseq[m_kdim - 1], Kseq[m_kdim]);
233 &Kseq[m_kdim][0], 1);
235 alpha[
m_kdim] = std::sqrt(alpha[m_kdim]);
236 Vmath::Smul(ntot, 1.0/alpha[m_kdim], Kseq[m_kdim], 1,
240 for (
int k = 0; k < m_kdim + 1; ++k)
246 EV_small(Tseq, ntot, alpha, m_kdim, zvec, wr, wi, resnorm);
249 converged =
EV_test(i, m_kdim, zvec, wr, wi, resnorm,
252 if (
m_comm->GetRank() == 0)
254 out <<
"Iteration: " << i <<
" (residual : "
255 << resid0 <<
")" <<endl;
266 for(j = 0; j <
m_equ[0]->GetNvariables(); ++j)
270 if (
m_comm->GetRank() == 0)
272 out <<
"L 2 error (variable " <<
m_equ[0]->GetVariable(j)
273 <<
") : " << vL2Error << endl;
274 out <<
"L inf error (variable " <<
m_equ[0]->GetVariable(j)
275 <<
") : " << vLinfError << endl;
280 EV_post(Tseq, Kseq, ntot, min(--i, m_kdim),
m_nvec, zvec, wr, wi,
284 "Need to implement Ritz re-evaluation of"
285 "eigenvalue. Only one half of complex "
286 "value will be correct");
293 if (
m_comm->GetRank() == 0)
309 m_equ[0]->TransCoeffToPhys();
316 fields =
m_equ[0]->UpdateFields();
320 m_equ[1]->TransCoeffToPhys();
344 int kdimp = kdim + 1;
352 for (
int i = 0; i < kdimp; ++i)
356 gsc = std::sqrt(gsc);
357 ASSERTL0(gsc != 0.0,
"Vectors are linearly independent.");
360 Vmath::Smul(ntot, 1.0/gsc, Kseq[i], 1, Kseq[i], 1);
362 for (
int j = i + 1; j < kdimp; ++j)
364 gsc =
Blas::Ddot(ntot, &Kseq[i][0], 1, &Kseq[j][0], 1);
366 Vmath::Svtvp(ntot, -gsc, Kseq[i], 1, Kseq[j], 1, Kseq[j], 1);
372 for (
int i = 0; i < kdim; ++i)
374 for (
int j = 0; j < kdim; ++j)
376 H[j*kdim+i] = alpha[j+1] * R[(j+1)*kdimp+i]
377 -
Vmath::Dot(j, &H[0] + i, kdim, &R[0] + j*kdimp, 1);
378 H[j*kdim+i] /= R[j*kdimp+j];
382 H[(kdim-1)*kdim+kdim] = alpha[kdim]
383 * std::fabs(R[kdim*kdimp+kdim] / R[(kdim-1)*kdimp + kdim-1]);
385 Lapack::dgeev_(
'N',
'V', kdim, &H[0], kdim, &wr[0], &wi[0], 0, 1,
386 &zvec[0], kdim, &rwork[0], lwork, ier);
390 resnorm = H[(kdim-1)*kdim + kdim];
412 for (
int i = 0; i < kdim; ++i)
415 &zvec[0] + i*kdim, 1));
416 resid[i] = resnorm * std::fabs(zvec[kdim - 1 + i*kdim]) / tmp;
419 resid[i-1] = resid[i] = hypot(resid[i-1], resid[i]);
423 EV_sort(zvec, wr, wi, resid, kdim);
430 if (
m_comm->GetRank() == 0)
432 evlout <<
"-- Iteration = " << itrn <<
", H(k+1, k) = "
435 evlout.setf(ios::scientific, ios::floatfield);
438 evlout <<
" Magnitude Angle Growth "
439 <<
"Frequency Residual" << endl;
443 evlout <<
" Real Imaginary inverse real "
444 <<
"inverse imag Residual" << endl;
447 for (
int i = 0; i < kdim; i++)
449 WriteEvs(evlout,i,wr[i],wi[i],resid[i]);
453 resid0 = resid[nvec-1];
470 for (
int j = 1; j < dim; ++j)
477 while (i >= 0 && test[i] > te_tmp)
482 Vmath::Vcopy(dim, &evec[0] + i*dim, 1, &evec[0] + (i+1)*dim, 1);
488 Vmath::Vcopy(dim, &z_tmp[0], 1, &evec[0] + (i+1)*dim, 1);
512 ASSERTL0(
false,
"Convergence was not achieved within the "
513 "prescribed number of iterations.");
518 ASSERTL0(
false,
"Minimum residual reached.");
520 else if (icon == nvec)
523 EV_big(Tseq, Kseq, ntot, kdim, icon, zvec, wr, wi);
525 =
m_equ[0]->UpdateFields();
527 for (
int j = 0; j < icon; ++j)
529 std::string file =
m_session->GetSessionName() +
"_eig_"
530 + boost::lexical_cast<std::string>(j)
540 ASSERTL0(
false,
"Unrecognised value.");
561 for (
int j = 0; j < nvec; ++j)
564 for (
int i = 0; i < kdim; ++i)
566 wgt = zvec[i + j*kdim];
567 Vmath::Svtvp(ntot, wgt, bvecs[i], 1, evecs[j], 1, evecs[j], 1);
572 for (
int i = 0; i < nvec; ++i)
576 norm =
Blas::Ddot(ntot, &evecs[i][0], 1, &evecs[i][0], 1);
578 norm = std::sqrt(norm);
579 Vmath::Smul(ntot, 1.0/norm, evecs[i], 1, evecs[i], 1);
583 norm =
Blas::Ddot(ntot, &evecs[i][0], 1, &evecs[i][0], 1);
584 norm +=
Blas::Ddot(ntot, &evecs[i+1][0], 1, &evecs[i+1][0], 1);
586 norm = std::sqrt(norm);
588 Vmath::Smul(ntot, 1.0/norm, evecs[i], 1, evecs[i], 1);
589 Vmath::Smul(ntot, 1.0/norm, evecs[i+1], 1, evecs[i+1], 1);
Array< OneD, NekDouble > m_imag_evl
#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.
int EV_test(const int itrn, const int kdim, Array< OneD, NekDouble > &zvec, Array< OneD, NekDouble > &wr, Array< OneD, NekDouble > &wi, const NekDouble resnorm, const int nvec, std::ofstream &evlout, NekDouble &resid0)
Tests for convergence of eigenvalues of H.
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
SOLVER_UTILS_EXPORT void ArnoldiSummary(std::ostream &out)
void EV_small(Array< OneD, Array< OneD, NekDouble > > &Kseq, const int ntot, const Array< OneD, NekDouble > &alpha, const int kdim, Array< OneD, NekDouble > &zvec, Array< OneD, NekDouble > &wr, Array< OneD, NekDouble > &wi, NekDouble &resnorm)
Generates the upper Hessenberg matrix H and computes its eigenvalues.
void CopyFieldToArnoldiArray(Array< OneD, NekDouble > &array)
Copy fields to Arnoldi storage.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
virtual void v_InitObject(std::ostream &out=std::cout)
Virtual function for initialisation implementation.
void WriteFld(std::string file, std::vector< Array< OneD, NekDouble > > coeffs)
Write coefficients to file.
void EV_sort(Array< OneD, NekDouble > &evec, Array< OneD, NekDouble > &wr, Array< OneD, NekDouble > &wi, Array< OneD, NekDouble > &test, const int dim)
Sorts a sequence of eigenvectors/eigenvalues by magnitude.
virtual ~DriverModifiedArnoldi()
Destructor.
void EV_update(Array< OneD, NekDouble > &src, Array< OneD, NekDouble > &tgt)
Generates a new vector in the sequence by applying the linear operator.
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
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)
enum EvolutionOperatorType m_EvolutionOperator
Evolution Operator.
T Ddot(int n, const Array< OneD, const T > &w, const int incw, const Array< OneD, const T > &x, const int incx, const Array< OneD, const int > &y, const int incy)
int m_nfields
interval to dump information if required.
bool m_timeSteppingAlgorithm
Period of time stepping algorithm.
T Dot(int n, const T *w, const T *x)
vvtvp (vector times vector times vector): z = w*x*y
Array< OneD, EquationSystemSharedPtr > m_equ
Equation system to solve.
void EV_post(Array< OneD, Array< OneD, NekDouble > > &Tseq, Array< OneD, Array< OneD, NekDouble > > &Kseq, const int ntot, const int kdim, const int nvec, Array< OneD, NekDouble > &zvec, Array< OneD, NekDouble > &wr, Array< OneD, NekDouble > &wi, const int icon)
void EV_big(Array< OneD, Array< OneD, NekDouble > > &bvecs, Array< OneD, Array< OneD, NekDouble > > &evecs, const int ntot, const int kdim, const int nvec, Array< OneD, NekDouble > &zvec, Array< OneD, NekDouble > &wr, Array< OneD, NekDouble > &wi)
void FillWhiteNoise(int n, const T eps, T *x, const int incx, int outseed)
Fills a vector with white noise.
DriverFactory & GetDriverFactory()
LibUtilities::CommSharedPtr m_comm
Communication object.
virtual void v_Execute(std::ostream &out=std::cout)
Virtual function for solve implementation.
Base class for the development of solvers.
void Zero(int n, T *x, const int incx)
Zero vector.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
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.