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;
105 int nq =
m_equ[0]->UpdateFields()[0]->GetNcoeffs();
110 std::string evlFile =
m_session->GetSessionName() +
".evl";
112 if (
m_comm->GetRank() == 0)
114 evlout.open(evlFile.c_str());
127 for (i = 0; i < m_kdim + 1; ++i)
134 if(
m_session->DefinesFunction(
"InitialConditions"))
136 if (
m_comm->GetRank() == 0)
138 out <<
"\tInital vector : specified in input file " << endl;
140 m_equ[0]->SetInitialConditions(0.0,
false);
146 if (
m_comm->GetRank() == 0)
148 out <<
"\tInital vector : random " << endl;
158 if (
m_comm->GetRank() == 0)
160 out <<
"Iteration: " << 0 << endl;
164 alpha[0] =
Blas::Ddot(ntot, &Kseq[0][0], 1, &Kseq[0][0], 1);
166 alpha[0] = std::sqrt(alpha[0]);
167 Vmath::Smul(ntot, 1.0/alpha[0], Kseq[0], 1, Kseq[0], 1);
171 for (i = 1; !converged && i <=
m_kdim; ++i)
177 alpha[i] =
Blas::Ddot(ntot, &Kseq[i][0], 1, &Kseq[i][0], 1);
179 alpha[i] = std::sqrt(alpha[i]);
182 Vmath::Smul(ntot, 1.0/alpha[i], Kseq[i], 1, Kseq[i], 1);
185 for (
int k = 0; k < i + 1; ++k)
191 EV_small(Tseq, ntot, alpha, i, zvec, wr, wi, resnorm);
194 converged =
EV_test(i, i, zvec, wr, wi, resnorm,
195 std::min(i,
m_nvec), evlout, resid0);
198 converged = max (converged, 0);
205 if (
m_comm->GetRank() == 0)
207 out <<
"Iteration: " << i <<
" (residual : " << resid0
215 for (i = m_kdim + 1; !converged && i <=
m_nits; ++i)
219 for (
int j = 1; j <=
m_kdim; ++j)
221 alpha[j-1] = alpha[j];
226 EV_update(Kseq[m_kdim - 1], Kseq[m_kdim]);
230 &Kseq[m_kdim][0], 1);
232 alpha[
m_kdim] = std::sqrt(alpha[m_kdim]);
233 Vmath::Smul(ntot, 1.0/alpha[m_kdim], Kseq[m_kdim], 1,
237 for (
int k = 0; k < m_kdim + 1; ++k)
243 EV_small(Tseq, ntot, alpha, m_kdim, zvec, wr, wi, resnorm);
246 converged =
EV_test(i, m_kdim, zvec, wr, wi, resnorm,
249 if (
m_comm->GetRank() == 0)
251 out <<
"Iteration: " << i <<
" (residual : "
252 << resid0 <<
")" <<endl;
263 for(j = 0; j <
m_equ[0]->GetNvariables(); ++j)
267 if (
m_comm->GetRank() == 0)
269 out <<
"L 2 error (variable " <<
m_equ[0]->GetVariable(j)
270 <<
") : " << vL2Error << endl;
271 out <<
"L inf error (variable " <<
m_equ[0]->GetVariable(j)
272 <<
") : " << vLinfError << endl;
277 EV_post(Tseq, Kseq, ntot, min(--i, m_kdim),
m_nvec, zvec, wr, wi,
281 "Need to implement Ritz re-evaluation of"
282 "eigenvalue. Only one half of complex "
283 "value will be correct");
290 if (
m_comm->GetRank() == 0)
306 m_equ[0]->TransCoeffToPhys();
313 fields =
m_equ[0]->UpdateFields();
317 m_equ[1]->TransCoeffToPhys();
341 int kdimp = kdim + 1;
349 for (
int i = 0; i < kdimp; ++i)
353 gsc = std::sqrt(gsc);
354 ASSERTL0(gsc != 0.0,
"Vectors are linearly independent.");
357 Vmath::Smul(ntot, 1.0/gsc, Kseq[i], 1, Kseq[i], 1);
359 for (
int j = i + 1; j < kdimp; ++j)
361 gsc =
Blas::Ddot(ntot, &Kseq[i][0], 1, &Kseq[j][0], 1);
363 Vmath::Svtvp(ntot, -gsc, Kseq[i], 1, Kseq[j], 1, Kseq[j], 1);
369 for (
int i = 0; i < kdim; ++i)
371 for (
int j = 0; j < kdim; ++j)
373 H[j*kdim+i] = alpha[j+1] * R[(j+1)*kdimp+i]
374 -
Vmath::Dot(j, &H[0] + i, kdim, &R[0] + j*kdimp, 1);
375 H[j*kdim+i] /= R[j*kdimp+j];
379 H[(kdim-1)*kdim+kdim] = alpha[kdim]
380 * std::fabs(R[kdim*kdimp+kdim] / R[(kdim-1)*kdimp + kdim-1]);
382 Lapack::dgeev_(
'N',
'V', kdim, &H[0], kdim, &wr[0], &wi[0], 0, 1,
383 &zvec[0], kdim, &rwork[0], lwork, ier);
387 resnorm = H[(kdim-1)*kdim + kdim];
409 for (
int i = 0; i < kdim; ++i)
412 &zvec[0] + i*kdim, 1));
413 resid[i] = resnorm * std::fabs(zvec[kdim - 1 + i*kdim]) / tmp;
416 resid[i-1] = resid[i] = hypot(resid[i-1], resid[i]);
420 EV_sort(zvec, wr, wi, resid, kdim);
427 if (
m_comm->GetRank() == 0)
429 evlout <<
"-- Iteration = " << itrn <<
", H(k+1, k) = "
432 evlout.setf(ios::scientific, ios::floatfield);
435 evlout <<
" Magnitude Angle Growth "
436 <<
"Frequency Residual" << endl;
440 evlout <<
" Real Imaginary inverse real "
441 <<
"inverse imag Residual" << endl;
444 for (
int i = 0; i < kdim; i++)
446 WriteEvs(evlout,i,wr[i],wi[i],resid[i]);
450 resid0 = resid[nvec-1];
467 for (
int j = 1; j < dim; ++j)
474 while (i >= 0 && test[i] > te_tmp)
479 Vmath::Vcopy(dim, &evec[0] + i*dim, 1, &evec[0] + (i+1)*dim, 1);
485 Vmath::Vcopy(dim, &z_tmp[0], 1, &evec[0] + (i+1)*dim, 1);
509 ASSERTL0(
false,
"Convergence was not achieved within the "
510 "prescribed number of iterations.");
515 ASSERTL0(
false,
"Minimum residual reached.");
517 else if (icon == nvec)
520 EV_big(Tseq, Kseq, ntot, kdim, icon, zvec, wr, wi);
522 =
m_equ[0]->UpdateFields();
524 for (
int j = 0; j < icon; ++j)
526 std::string file =
m_session->GetSessionName() +
"_eig_"
527 + boost::lexical_cast<std::string>(j)
537 ASSERTL0(
false,
"Unrecognised value.");
558 for (
int j = 0; j < nvec; ++j)
561 for (
int i = 0; i < kdim; ++i)
563 wgt = zvec[i + j*kdim];
564 Vmath::Svtvp(ntot, wgt, bvecs[i], 1, evecs[j], 1, evecs[j], 1);
569 for (
int i = 0; i < nvec; ++i)
573 norm =
Blas::Ddot(ntot, &evecs[i][0], 1, &evecs[i][0], 1);
575 norm = std::sqrt(norm);
576 Vmath::Smul(ntot, 1.0/norm, evecs[i], 1, evecs[i], 1);
580 norm =
Blas::Ddot(ntot, &evecs[i][0], 1, &evecs[i][0], 1);
581 norm +=
Blas::Ddot(ntot, &evecs[i+1][0], 1, &evecs[i+1][0], 1);
583 norm = std::sqrt(norm);
585 Vmath::Smul(ntot, 1.0/norm, evecs[i], 1, evecs[i], 1);
586 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.