68 m_equ[0]->PrintSummary(out);
71 if (
m_comm->GetRank() == 0)
73 out <<
"\tArnoldi solver type : Modified Arnoldi" << std::endl;
78 for (
int i = 0; i <
m_nequ; ++i)
80 m_equ[i]->DoInitialise();
94 int nm =
m_equ[0]->UpdateFields()[0]->GetNcoeffs();
99 std::string evlFile =
m_session->GetSessionName() +
".evl";
101 if (
m_comm->GetRank() == 0)
103 evlout.open(evlFile.c_str());
118 for (i = 0; i <
m_kdim + 1; ++i)
127 Kseqcopy[i] = Kseq[i];
133 if (
m_session->DefinesFunction(
"InitialConditions"))
135 if (
m_comm->GetRank() == 0)
137 out <<
"\tInital vector : specified in input file "
140 m_equ[0]->SetInitialConditions(0.0,
false);
146 if (
m_comm->GetRank() == 0)
148 out <<
"\tInital vector : random " << std::endl;
154 auto contfield = std::dynamic_pointer_cast<MultiRegions::ContField>(
155 m_equ[0]->UpdateFields()[0]);
159 contfield->GetLocalToGlobalMap()->GetNumGlobalCoeffs();
163 contfield->LocalToGlobal(Kseq[1] + n * nm, global);
164 contfield->GlobalToLocal(global, tmp = Kseq[1] + n * nm);
176 if (
m_comm->GetRank() == 0)
178 out <<
"Iteration: " << 0 << std::endl;
186 alpha[0] =
Blas::Ddot(ntot, &Kseqcopy[0][0], 1, &Kseqcopy[0][0], 1);
189 Vmath::Smul(ntot, 1.0 / alpha[0], Kseq[0], 1, Kseq[0], 1);
193 for (i = 1; !converged && i <=
m_kdim; ++i)
203 alpha[i] =
Blas::Ddot(ntot, &Kseqcopy[i][0], 1, &Kseqcopy[i][0], 1);
208 Vmath::Smul(ntot, 1.0 / alpha[i], Kseq[i], 1, Kseq[i], 1);
211 for (
int k = 0; k < i + 1; ++k)
226 EV_small(Tseq, Kseqcopy, ntot, alpha, i, zvec, wr, wi, resnorm);
229 converged =
EV_test(i, i, zvec, wr, wi, resnorm, std::min(i,
m_nvec),
233 converged = std::max(converged, 0);
240 if (
m_comm->GetRank() == 0)
242 out <<
"Iteration: " << i <<
" (residual : " << resid0 <<
")"
254 for (
int j = 1; j <=
m_kdim; ++j)
256 alpha[j - 1] = alpha[j];
277 for (
int k = 0; k <
m_kdim + 1; ++k)
299 if (
m_comm->GetRank() == 0)
301 out <<
"Iteration: " << i <<
" (residual : " << resid0 <<
")"
313 for (j = 0; j <
m_equ[0]->GetNvariables(); ++j)
317 if (
m_comm->GetRank() == 0)
319 out <<
"L 2 error (variable " <<
m_equ[0]->GetVariable(j)
320 <<
") : " << vL2Error << std::endl;
321 out <<
"L inf error (variable " <<
m_equ[0]->GetVariable(j)
322 <<
") : " << vLinfError << std::endl;
332 "Need to implement Ritz re-evaluation of"
333 "eigenvalue. Only one half of complex "
334 "value will be correct");
341 if (
m_comm->GetRank() == 0)
355 m_equ[0]->TransCoeffToPhys();
357 m_equ[0]->SetTime(0.);
363 fields =
m_equ[0]->UpdateFields();
367 m_equ[1]->TransCoeffToPhys();
369 m_equ[1]->SetTime(0.);
388 int kdimp = kdim + 1;
389 int lwork = 10 * kdim;
396 for (
int i = 0; i < kdimp; ++i)
401 ASSERTL0(gsc != 0.0,
"Vectors are linearly independent.");
403 R[i * kdimp + i] = gsc;
404 Vmath::Smul(ntot, 1.0 / gsc, Kseq[i], 1, Kseq[i], 1);
407 Vmath::Smul(ntot, 1.0 / gsc, Kseqcopy[i], 1, Kseqcopy[i], 1);
410 for (
int j = i + 1; j < kdimp; ++j)
412 gsc =
Blas::Ddot(ntot, &Kseq[i][0], 1, &Kseq[j][0], 1);
414 Vmath::Svtvp(ntot, -gsc, Kseq[i], 1, Kseq[j], 1, Kseq[j], 1);
417 Vmath::Svtvp(ntot, -gsc, Kseqcopy[i], 1, Kseqcopy[j], 1,
420 R[j * kdimp + i] = gsc;
425 for (
int i = 0; i < kdim; ++i)
427 for (
int j = 0; j < kdim; ++j)
430 alpha[j + 1] * R[(j + 1) * kdimp + i] -
431 Vmath::Dot(j, &H[0] + i, kdim, &R[0] + j * kdimp, 1);
432 H[j * kdim + i] /= R[j * kdimp + j];
436 H[(kdim - 1) * kdim + kdim] =
438 std::fabs(R[kdim * kdimp + kdim] / R[(kdim - 1) * kdimp + kdim - 1]);
440 Lapack::dgeev_(
'N',
'V', kdim, &H[0], kdim, &wr[0], &wi[0],
nullptr, 1,
441 &zvec[0], kdim, &rwork[0], lwork, ier);
445 resnorm = H[(kdim - 1) * kdim + kdim];
456 std::ofstream &evlout,
NekDouble &resid0)
462 for (
int i = 0; i < kdim; ++i)
465 Vmath::Dot(kdim, &zvec[0] + i * kdim, 1, &zvec[0] + i * kdim, 1));
466 resid[i] = resnorm * std::fabs(zvec[kdim - 1 + i * kdim]) / tmp;
469 resid[i - 1] = resid[i] = hypot(resid[i - 1], resid[i]);
473 EV_sort(zvec, wr, wi, resid, kdim);
475 while (nvec <= kdim && resid[nvec - 1] <
m_evtol)
482 if (
m_comm->GetRank() == 0)
484 evlout <<
"-- Iteration = " << itrn <<
", H(k+1, k) = " << resnorm
487 evlout.setf(std::ios::scientific, std::ios::floatfield);
490 evlout <<
" Magnitude Angle Growth "
491 <<
"Frequency Residual" << std::endl;
495 evlout <<
" Real Imaginary inverse real "
496 <<
"inverse imag Residual" << std::endl;
499 for (
int i = 0; i < kdim; i++)
501 WriteEvs(evlout, i, wr[i], wi[i], resid[i]);
505 resid0 = resid[nvec - 1];
519 for (
int j = 1; j < dim; ++j)
526 while (i >= 0 &&
test[i] > te_tmp)
531 Vmath::Vcopy(dim, &evec[0] + i * dim, 1, &evec[0] + (i + 1) * dim,
537 test[i + 1] = te_tmp;
538 Vmath::Vcopy(dim, &z_tmp[0], 1, &evec[0] + (i + 1) * dim, 1);
549 const int ntot,
const int kdim,
558 ASSERTL0(
false,
"Convergence was not achieved within the "
559 "prescribed number of iterations.");
564 ASSERTL0(
false,
"Minimum residual reached.");
566 else if (icon == nvec)
569 EV_big(Tseq, Kseq, ntot, kdim, icon, zvec, wr, wi);
571 m_equ[0]->UpdateFields();
573 for (
int j = 0; j < icon; ++j)
575 std::string file =
m_session->GetSessionName() +
"_eig_" +
576 std::to_string(j) +
".fld";
578 if (
m_comm->GetRank() == 0)
580 WriteEvs(std::cout, j, wr[j], wi[j]);
585 std::string fileunmask =
m_session->GetSessionName() +
586 "_eig_masked_" + std::to_string(j) +
595 ASSERTL0(
false,
"Unrecognised value.");
604 const int ntot,
const int kdim,
612 for (
int i = 0; i < nvec; ++i)
626 for (
int j = 0; j < nvec; ++j)
633 for (
int i = 0; i < kdim; ++i)
635 wgt = zvec[i + j * kdim];
636 Vmath::Svtvp(ntot, wgt, bvecs[i], 1, btmp[j], 1, btmp[j], 1);
639 Vmath::Svtvp(ntot, wgt, evecs[i], 1, etmp[j], 1, etmp[j], 1);
645 for (
int i = 0; i < nvec; ++i)
649 if (
m_session->GetComm()->GetRank() == 0)
651 std::cout <<
"eigenvalue " << i <<
": real mode" << std::endl;
653 norm =
Blas::Ddot(ntot, &btmp[i][0], 1, &btmp[i][0], 1);
658 Vmath::Smul(ntot, 1.0 / norm, btmp[i], 1, bvecs[i], 1);
659 Vmath::Smul(ntot, 1.0 / norm, etmp[i], 1, evecs[i], 1);
663 Vmath::Smul(ntot, 1.0 / norm, btmp[i], 1, evecs[i], 1);
668 if (
m_session->GetComm()->GetRank() == 0)
670 std::cout <<
"eigenvalues " << i <<
", " << i + 1
671 <<
": complex modes" << std::endl;
673 norm =
Blas::Ddot(ntot, &btmp[i][0], 1, &btmp[i][0], 1);
674 norm +=
Blas::Ddot(ntot, &btmp[i + 1][0], 1, &btmp[i + 1][0], 1);
680 Vmath::Smul(ntot, 1.0 / norm, btmp[i], 1, bvecs[i], 1);
681 Vmath::Smul(ntot, 1.0 / norm, btmp[i + 1], 1, bvecs[i + 1], 1);
682 Vmath::Smul(ntot, 1.0 / norm, etmp[i], 1, evecs[i], 1);
683 Vmath::Smul(ntot, 1.0 / norm, etmp[i + 1], 1, evecs[i + 1], 1);
687 Vmath::Smul(ntot, 1.0 / norm, btmp[i], 1, evecs[i], 1);
688 Vmath::Smul(ntot, 1.0 / norm, btmp[i + 1], 1, evecs[i + 1], 1);
#define ASSERTL0(condition, msg)
#define WARNINGL0(condition, msg)
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
static std::string RegisterEnumValue(std::string pEnum, std::string pString, int pEnumValue)
Registers an enumeration value.
Base class for the development of solvers.
void CopyFwdToAdj()
Copy the forward field to the adjoint system in transient growth calculations.
void v_InitObject(std::ostream &out=std::cout) override
Virtual function for initialisation implementation.
void WriteFld(std::string file, std::vector< Array< OneD, NekDouble > > coeffs)
Write coefficients to file.
Array< OneD, NekDouble > m_maskCoeffs
void CopyFieldToArnoldiArray(Array< OneD, NekDouble > &array)
Copy fields to Arnoldi storage.
int m_nvec
Dimension of Krylov subspace.
bool m_timeSteppingAlgorithm
Period of time stepping algorithm.
int m_nits
Number of vectors to test.
Array< OneD, NekDouble > m_imag_evl
void CopyArnoldiArrayToField(Array< OneD, NekDouble > &array)
Copy Arnoldi storage to fields.
NekDouble m_evtol
Maxmum number of iterations.
int m_nfields
interval to dump information if required.
SOLVER_UTILS_EXPORT void ArnoldiSummary(std::ostream &out)
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)
LibUtilities::SessionReaderSharedPtr m_session
Session reader object.
LibUtilities::CommSharedPtr m_comm
Communication object.
enum EvolutionOperatorType m_EvolutionOperator
Evolution Operator.
Array< OneD, EquationSystemSharedPtr > m_equ
Equation system to solve.
int m_nequ
number of equations
void EV_update(Array< OneD, NekDouble > &src, Array< OneD, NekDouble > &tgt)
Generates a new vector in the sequence by applying the linear operator.
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)
int EV_test(const int itrn, const int kdim, Array< OneD, NekDouble > &zvec, Array< OneD, NekDouble > &wr, Array< OneD, NekDouble > &wi, const NekDouble resnorm, int nvec, std::ofstream &evlout, NekDouble &resid0)
Tests for convergence of eigenvalues of H.
void EV_small(Array< OneD, Array< OneD, NekDouble > > &Kseq, Array< OneD, Array< OneD, NekDouble > > &Kseqcopy, 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 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.
static DriverSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
DriverModifiedArnoldi(const LibUtilities::SessionReaderSharedPtr pSession, const SpatialDomains::MeshGraphSharedPtr pGraph)
Constructor.
static std::string className
Name of the class.
static std::string driverLookupId
void v_InitObject(std::ostream &out=std::cout) override
Virtual function for initialisation implementation.
void v_Execute(std::ostream &out=std::cout) override
Virtual function for solve implementation.
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)
static double Ddot(const int &n, const double *x, const int &incx, const double *y, const int &incy)
BLAS level 1: output = .
std::shared_ptr< SessionReader > SessionReaderSharedPtr
DriverFactory & GetDriverFactory()
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
MultiRegions::ContFieldSharedPtr contfield
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
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.
T Dot(int n, const T *w, const T *x)
dot product
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
void Zero(int n, T *x, const int incx)
Zero vector.
void FillWhiteNoise(int n, const T eps, T *x, const int incx, int outseed)
Fills a vector with white noise.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
scalarT< T > sqrt(scalarT< T > in)