Nektar++
Static Public Member Functions | Static Public Attributes | Protected Member Functions | Static Protected Attributes | Private Member Functions | Friends | List of all members
Nektar::SolverUtils::DriverModifiedArnoldi Class Reference

#include <DriverModifiedArnoldi.h>

Inheritance diagram for Nektar::SolverUtils::DriverModifiedArnoldi:
[legend]

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

 DriverModifiedArnoldi (const LibUtilities::SessionReaderSharedPtr pSession, const SpatialDomains::MeshGraphSharedPtr pGraph)
 Constructor. More...
 
 ~DriverModifiedArnoldi () override=default
 Destructor. More...
 
void v_InitObject (std::ostream &out=std::cout) override
 Virtual function for initialisation implementation. More...
 
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...
 
 ~DriverArnoldi () override=default
 Destructor. More...
 
void v_InitObject (std::ostream &out=std::cout) override
 Virtual function for initialisation implementation. More...
 
void v_Execute (std::ostream &out=std::cout) override
 Virtual function for solve implementation. 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 GetMaskInfo (std::vector< std::vector< LibUtilities::EquationSharedPtr > > &selectedDomains, std::set< int > &unselectedVariables)
 
SOLVER_UTILS_EXPORT void ArnoldiSummary (std::ostream &out)
 
- Protected Member Functions inherited from Nektar::SolverUtils::Driver
 Driver (const LibUtilities::SessionReaderSharedPtr pSession, const SpatialDomains::MeshGraphSharedPtr pGraph)
 Initialises EquationSystem class members. More...
 
virtual SOLVER_UTILS_EXPORT void v_InitObject (std::ostream &out=std::cout)
 Virtual function for initialisation implementation. More...
 
virtual SOLVER_UTILS_EXPORT void v_Execute (std::ostream &out=std::cout)=0
 Virtual function for solve implementation. More...
 

Static Protected Attributes

static std::string driverLookupId
 
- Static Protected Attributes inherited from Nektar::SolverUtils::Driver
static std::string evolutionOperatorLookupIds []
 
static std::string evolutionOperatorDef
 
static std::string driverDefault
 

Private Member Functions

void EV_update (Array< OneD, NekDouble > &src, Array< OneD, NekDouble > &tgt)
 Generates a new vector in the sequence by applying the linear operator. More...
 
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. More...
 
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. More...
 
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. More...
 
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)
 

Friends

class MemoryManager< DriverModifiedArnoldi >
 

Additional Inherited Members

- Public Member Functions inherited from Nektar::SolverUtils::DriverArnoldi
Array< OneD, NekDoubleGetRealEvl (void)
 
Array< OneD, NekDoubleGetImagEvl (void)
 
- Public Member Functions inherited from Nektar::SolverUtils::Driver
virtual ~Driver ()=default
 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, EquationSystemSharedPtrGetEqu ()
 
- 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 iterations. 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, NekDoublem_real_evl
 Operator in solve call is negated. More...
 
Array< OneD, NekDoublem_imag_evl
 
bool m_useMask
 
Array< OneD, NekDoublem_maskCoeffs
 
Array< OneD, NekDoublem_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
 Coupling between SFD and arnoldi. More...
 
SpatialDomains::MeshGraphSharedPtr m_graph
 MeshGraph object. More...
 
Array< OneD, EquationSystemSharedPtrm_equ
 Equation system to solve. More...
 
int m_nequ
 number of equations More...
 
enum EvolutionOperatorType m_EvolutionOperator
 Evolution Operator. More...
 

Detailed Description

Definition at line 44 of file DriverModifiedArnoldi.h.

Constructor & Destructor Documentation

◆ DriverModifiedArnoldi()

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

Constructor.

Definition at line 54 of file DriverModifiedArnoldi.cpp.

57 : DriverArnoldi(pSession, pGraph)
58{
59}
DriverArnoldi(const LibUtilities::SessionReaderSharedPtr pSession, const SpatialDomains::MeshGraphSharedPtr pGraph)
Constructor.

◆ ~DriverModifiedArnoldi()

Nektar::SolverUtils::DriverModifiedArnoldi::~DriverModifiedArnoldi ( )
overrideprotecteddefault

Destructor.

Member Function Documentation

◆ create()

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

Creates an instance of this class.

Definition at line 50 of file DriverModifiedArnoldi.h.

53 {
56 pGraph);
57 p->InitObject();
58 return p;
59 }
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
std::shared_ptr< Driver > DriverSharedPtr
A shared pointer to a Driver object.
Definition: Driver.h:52

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

◆ EV_big()

void Nektar::SolverUtils::DriverModifiedArnoldi::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 
)
private

Compute estimates of the eigenvalues/eigenvectors of the original system.

Definition at line 602 of file DriverModifiedArnoldi.cpp.

608{
609 NekDouble wgt, norm;
610 Array<OneD, Array<OneD, NekDouble>> btmp(nvec);
611 Array<OneD, Array<OneD, NekDouble>> etmp(nvec);
612 for (int i = 0; i < nvec; ++i)
613 {
614 if (m_useMask)
615 {
616 btmp[i] = Array<OneD, NekDouble>(ntot, 0.);
617 etmp[i] = Array<OneD, NekDouble>(ntot, 0.);
618 }
619 else
620 {
621 btmp[i] = evecs[i];
622 }
623 }
624
625 // Generate big eigenvectors
626 for (int j = 0; j < nvec; ++j)
627 {
628 Vmath::Zero(ntot, btmp[j], 1);
629 if (m_useMask)
630 {
631 Vmath::Zero(ntot, etmp[j], 1);
632 }
633 for (int i = 0; i < kdim; ++i)
634 {
635 wgt = zvec[i + j * kdim];
636 Vmath::Svtvp(ntot, wgt, bvecs[i], 1, btmp[j], 1, btmp[j], 1);
637 if (m_useMask)
638 {
639 Vmath::Svtvp(ntot, wgt, evecs[i], 1, etmp[j], 1, etmp[j], 1);
640 }
641 }
642 }
643
644 // Normalise the big eigenvectors
645 for (int i = 0; i < nvec; ++i)
646 {
647 if (wi[i] == 0.0) // Real mode
648 {
649 if (m_session->GetComm()->GetRank() == 0)
650 {
651 std::cout << "eigenvalue " << i << ": real mode" << std::endl;
652 }
653 norm = Blas::Ddot(ntot, &btmp[i][0], 1, &btmp[i][0], 1);
654 m_comm->AllReduce(norm, Nektar::LibUtilities::ReduceSum);
655 norm = std::sqrt(norm);
656 if (m_useMask)
657 {
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);
660 }
661 else
662 {
663 Vmath::Smul(ntot, 1.0 / norm, btmp[i], 1, evecs[i], 1);
664 }
665 }
666 else
667 {
668 if (m_session->GetComm()->GetRank() == 0)
669 {
670 std::cout << "eigenvalues " << i << ", " << i + 1
671 << ": complex modes" << std::endl;
672 }
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);
675 m_comm->AllReduce(norm, Nektar::LibUtilities::ReduceSum);
676 norm = std::sqrt(norm);
677
678 if (m_useMask)
679 {
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);
684 }
685 else
686 {
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);
689 }
690
691 i++;
692 }
693 }
694}
LibUtilities::SessionReaderSharedPtr m_session
Session reader object.
Definition: Driver.h:83
LibUtilities::CommSharedPtr m_comm
Communication object.
Definition: Driver.h:80
static double Ddot(const int &n, const double *x, const int &incx, const double *y, const int &incy)
BLAS level 1: output = .
Definition: Blas.hpp:163
double NekDouble
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.
Definition: Vmath.hpp:396
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.hpp:100
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.hpp:273
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:285

References Blas::Ddot(), Nektar::SolverUtils::Driver::m_comm, Nektar::SolverUtils::Driver::m_session, Nektar::SolverUtils::DriverArnoldi::m_useMask, Nektar::LibUtilities::ReduceSum, Vmath::Smul(), tinysimd::sqrt(), Vmath::Svtvp(), and Vmath::Zero().

Referenced by EV_post().

◆ EV_post()

void Nektar::SolverUtils::DriverModifiedArnoldi::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 
)
private

Post-process the Ritz eigenvalues/eigenvectors of the matrix H, to compute estimations of the leading eigenvalues and eigenvectors of the original matrix.

Definition at line 547 of file DriverModifiedArnoldi.cpp.

554{
555 if (icon == 0)
556 {
557 // Not converged, write final Krylov vector
558 ASSERTL0(false, "Convergence was not achieved within the "
559 "prescribed number of iterations.");
560 }
561 else if (icon < 0)
562 {
563 // Minimum residual reached
564 ASSERTL0(false, "Minimum residual reached.");
565 }
566 else if (icon == nvec)
567 {
568 // Converged, write out eigenvectors
569 EV_big(Tseq, Kseq, ntot, kdim, icon, zvec, wr, wi);
570 Array<OneD, MultiRegions::ExpListSharedPtr> fields =
571 m_equ[0]->UpdateFields();
572
573 for (int j = 0; j < icon; ++j)
574 {
575 std::string file = m_session->GetSessionName() + "_eig_" +
576 std::to_string(j) + ".fld";
577
578 if (m_comm->GetRank() == 0)
579 {
580 WriteEvs(std::cout, j, wr[j], wi[j]);
581 }
582 WriteFld(file, Kseq[j]);
583 if (m_useMask)
584 {
585 std::string fileunmask = m_session->GetSessionName() +
586 "_eig_masked_" + std::to_string(j) +
587 ".fld";
588 WriteFld(fileunmask, Tseq[j]);
589 }
590 }
591 }
592 else
593 {
594 // Not recognised
595 ASSERTL0(false, "Unrecognised value.");
596 }
597}
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
void WriteFld(std::string file, std::vector< Array< OneD, NekDouble > > coeffs)
Write coefficients to file.
void WriteEvs(std::ostream &evlout, const int k, const NekDouble real, const NekDouble imag, NekDouble resid=NekConstants::kNekUnsetDouble, bool DumpInverse=true)
Array< OneD, EquationSystemSharedPtr > m_equ
Equation system to solve.
Definition: Driver.h:92
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)

References ASSERTL0, EV_big(), Nektar::SolverUtils::Driver::m_comm, Nektar::SolverUtils::Driver::m_equ, Nektar::SolverUtils::Driver::m_session, Nektar::SolverUtils::DriverArnoldi::m_useMask, Nektar::SolverUtils::DriverArnoldi::WriteEvs(), and Nektar::SolverUtils::DriverArnoldi::WriteFld().

Referenced by v_Execute().

◆ EV_small()

void Nektar::SolverUtils::DriverModifiedArnoldi::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 
)
private

Generates the upper Hessenberg matrix H and computes its eigenvalues.

Computes the Ritz eigenvalues and eigenvectors of the Hessenberg matrix constructed from the Krylov sequence.

Definition at line 381 of file DriverModifiedArnoldi.cpp.

387{
388 int kdimp = kdim + 1;
389 int lwork = 10 * kdim;
390 int ier;
391 Array<OneD, NekDouble> R(kdimp * kdimp, 0.0);
392 Array<OneD, NekDouble> H(kdimp * kdim, 0.0);
393 Array<OneD, NekDouble> rwork(lwork, 0.0);
394
395 // Modified G-S orthonormalisation
396 for (int i = 0; i < kdimp; ++i)
397 {
398 NekDouble gsc = Blas::Ddot(ntot, &Kseq[i][0], 1, &Kseq[i][0], 1);
400 gsc = std::sqrt(gsc);
401 ASSERTL0(gsc != 0.0, "Vectors are linearly independent.");
402
403 R[i * kdimp + i] = gsc;
404 Vmath::Smul(ntot, 1.0 / gsc, Kseq[i], 1, Kseq[i], 1);
405 if (m_useMask)
406 {
407 Vmath::Smul(ntot, 1.0 / gsc, Kseqcopy[i], 1, Kseqcopy[i], 1);
408 }
409
410 for (int j = i + 1; j < kdimp; ++j)
411 {
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);
415 if (m_useMask)
416 {
417 Vmath::Svtvp(ntot, -gsc, Kseqcopy[i], 1, Kseqcopy[j], 1,
418 Kseqcopy[j], 1);
419 }
420 R[j * kdimp + i] = gsc;
421 }
422 }
423
424 // Compute H matrix
425 for (int i = 0; i < kdim; ++i)
426 {
427 for (int j = 0; j < kdim; ++j)
428 {
429 H[j * kdim + i] =
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];
433 }
434 }
435
436 H[(kdim - 1) * kdim + kdim] =
437 alpha[kdim] *
438 std::fabs(R[kdim * kdimp + kdim] / R[(kdim - 1) * kdimp + kdim - 1]);
439
440 Lapack::dgeev_('N', 'V', kdim, &H[0], kdim, &wr[0], &wi[0], nullptr, 1,
441 &zvec[0], kdim, &rwork[0], lwork, ier);
442
443 ASSERTL0(!ier, "Error with dgeev");
444
445 resnorm = H[(kdim - 1) * kdim + kdim];
446}
T Dot(int n, const T *w, const T *x)
dot product
Definition: Vmath.hpp:761

References ASSERTL0, Blas::Ddot(), Vmath::Dot(), Nektar::SolverUtils::Driver::m_comm, Nektar::SolverUtils::DriverArnoldi::m_useMask, Nektar::LibUtilities::ReduceSum, Vmath::Smul(), tinysimd::sqrt(), and Vmath::Svtvp().

Referenced by v_Execute().

◆ EV_sort()

void Nektar::SolverUtils::DriverModifiedArnoldi::EV_sort ( Array< OneD, NekDouble > &  evec,
Array< OneD, NekDouble > &  wr,
Array< OneD, NekDouble > &  wi,
Array< OneD, NekDouble > &  test,
const int  dim 
)
private

Sorts a sequence of eigenvectors/eigenvalues by magnitude.

Sorts the computed eigenvalues by smallest residual first.

Definition at line 512 of file DriverModifiedArnoldi.cpp.

516{
517 Array<OneD, NekDouble> z_tmp(dim, 0.0);
518 NekDouble wr_tmp, wi_tmp, te_tmp;
519 for (int j = 1; j < dim; ++j)
520 {
521 wr_tmp = wr[j];
522 wi_tmp = wi[j];
523 te_tmp = test[j];
524 Vmath::Vcopy(dim, &evec[0] + j * dim, 1, &z_tmp[0], 1);
525 int i = j - 1;
526 while (i >= 0 && test[i] > te_tmp)
527 {
528 wr[i + 1] = wr[i];
529 wi[i + 1] = wi[i];
530 test[i + 1] = test[i];
531 Vmath::Vcopy(dim, &evec[0] + i * dim, 1, &evec[0] + (i + 1) * dim,
532 1);
533 i--;
534 }
535 wr[i + 1] = wr_tmp;
536 wi[i + 1] = wi_tmp;
537 test[i + 1] = te_tmp;
538 Vmath::Vcopy(dim, &z_tmp[0], 1, &evec[0] + (i + 1) * dim, 1);
539 }
540}
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825

References Nektar::UnitTests::test(), and Vmath::Vcopy().

Referenced by EV_test().

◆ EV_test()

int Nektar::SolverUtils::DriverModifiedArnoldi::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 
)
private

Tests for convergence of eigenvalues of H.

Check for convergence of the residuals of the eigenvalues of H.

Definition at line 451 of file DriverModifiedArnoldi.cpp.

457{
458 int idone = 0;
459 // NekDouble period = 0.1;
460
461 Array<OneD, NekDouble> resid(kdim);
462 for (int i = 0; i < kdim; ++i)
463 {
464 NekDouble tmp = std::sqrt(
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;
467 if (wi[i] < 0.0)
468 {
469 resid[i - 1] = resid[i] = hypot(resid[i - 1], resid[i]);
470 }
471 }
472
473 EV_sort(zvec, wr, wi, resid, kdim);
474
475 while (nvec <= kdim && resid[nvec - 1] < m_evtol)
476 {
477 idone = nvec;
478 ++nvec;
479 }
480 nvec -= (idone > 0);
481
482 if (m_comm->GetRank() == 0)
483 {
484 evlout << "-- Iteration = " << itrn << ", H(k+1, k) = " << resnorm
485 << std::endl;
486 evlout.precision(4);
487 evlout.setf(std::ios::scientific, std::ios::floatfield);
489 {
490 evlout << " Magnitude Angle Growth "
491 << "Frequency Residual" << std::endl;
492 }
493 else
494 {
495 evlout << " Real Imaginary inverse real "
496 << "inverse imag Residual" << std::endl;
497 }
498
499 for (int i = 0; i < kdim; i++)
500 {
501 WriteEvs(evlout, i, wr[i], wi[i], resid[i]);
502 }
503 }
504
505 resid0 = resid[nvec - 1];
506 return idone;
507}
bool m_timeSteppingAlgorithm
Period of time stepping algorithm.
Definition: DriverArnoldi.h:65
NekDouble m_evtol
Maxmum number of iterations.
Definition: DriverArnoldi.h:63
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.

References Vmath::Dot(), EV_sort(), Nektar::SolverUtils::Driver::m_comm, Nektar::SolverUtils::DriverArnoldi::m_evtol, Nektar::SolverUtils::DriverArnoldi::m_timeSteppingAlgorithm, tinysimd::sqrt(), and Nektar::SolverUtils::DriverArnoldi::WriteEvs().

Referenced by v_Execute().

◆ EV_update()

void Nektar::SolverUtils::DriverModifiedArnoldi::EV_update ( Array< OneD, NekDouble > &  src,
Array< OneD, NekDouble > &  tgt 
)
private

Generates a new vector in the sequence by applying the linear operator.

Definition at line 350 of file DriverModifiedArnoldi.cpp.

352{
353 // Copy starting vector into first sequence element.
355 m_equ[0]->TransCoeffToPhys();
356
357 m_equ[0]->SetTime(0.);
358 m_equ[0]->DoSolve();
359
361 {
362 Array<OneD, MultiRegions::ExpListSharedPtr> fields;
363 fields = m_equ[0]->UpdateFields();
364
365 // start Adjoint with latest fields of direct
366 CopyFwdToAdj();
367 m_equ[1]->TransCoeffToPhys();
368
369 m_equ[1]->SetTime(0.);
370 m_equ[1]->DoSolve();
371 }
372
373 // Copy starting vector into first sequence element.
375}
void CopyFwdToAdj()
Copy the forward field to the adjoint system in transient growth calculations.
void CopyFieldToArnoldiArray(Array< OneD, NekDouble > &array)
Copy fields to Arnoldi storage.
void CopyArnoldiArrayToField(Array< OneD, NekDouble > &array)
Copy Arnoldi storage to fields.
enum EvolutionOperatorType m_EvolutionOperator
Evolution Operator.
Definition: Driver.h:98

References Nektar::SolverUtils::DriverArnoldi::CopyArnoldiArrayToField(), Nektar::SolverUtils::DriverArnoldi::CopyFieldToArnoldiArray(), Nektar::SolverUtils::DriverArnoldi::CopyFwdToAdj(), Nektar::SolverUtils::eTransientGrowth, Nektar::SolverUtils::Driver::m_equ, and Nektar::SolverUtils::Driver::m_EvolutionOperator.

Referenced by v_Execute().

◆ v_Execute()

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

Virtual function for solve implementation.

Reimplemented from Nektar::SolverUtils::DriverArnoldi.

Reimplemented in Nektar::SolverUtils::DriverSteadyState.

Definition at line 90 of file DriverModifiedArnoldi.cpp.

91{
92 int i = 0;
93 int j = 0;
94 int nm = m_equ[0]->UpdateFields()[0]->GetNcoeffs();
95 int ntot = m_nfields * nm;
96 int converged = 0;
97 NekDouble resnorm = 0.0;
98 std::ofstream evlout;
99 std::string evlFile = m_session->GetSessionName() + ".evl";
100
101 if (m_comm->GetRank() == 0)
102 {
103 evlout.open(evlFile.c_str());
104 }
105
106 // Allocate memory
107 Array<OneD, NekDouble> alpha = Array<OneD, NekDouble>(m_kdim + 1, 0.0);
108 Array<OneD, NekDouble> wr = Array<OneD, NekDouble>(m_kdim, 0.0);
109 Array<OneD, NekDouble> wi = Array<OneD, NekDouble>(m_kdim, 0.0);
110 Array<OneD, NekDouble> zvec = Array<OneD, NekDouble>(m_kdim * m_kdim, 0.0);
111
112 Array<OneD, Array<OneD, NekDouble>> Kseq =
113 Array<OneD, Array<OneD, NekDouble>>(m_kdim + 1);
114 Array<OneD, Array<OneD, NekDouble>> Kseqcopy =
115 Array<OneD, Array<OneD, NekDouble>>(m_kdim + 1);
116 Array<OneD, Array<OneD, NekDouble>> Tseq =
117 Array<OneD, Array<OneD, NekDouble>>(m_kdim + 1);
118 for (i = 0; i < m_kdim + 1; ++i)
119 {
120 Kseq[i] = Array<OneD, NekDouble>(ntot, 0.0);
121 if (m_useMask)
122 {
123 Kseqcopy[i] = Array<OneD, NekDouble>(ntot, 0.0);
124 }
125 else
126 {
127 Kseqcopy[i] = Kseq[i];
128 }
129 Tseq[i] = Array<OneD, NekDouble>(ntot, 0.0);
130 }
131
132 // Copy starting vector into second sequence element (temporary).
133 if (m_session->DefinesFunction("InitialConditions"))
134 {
135 if (m_comm->GetRank() == 0)
136 {
137 out << "\tInital vector : specified in input file "
138 << std::endl;
139 }
140 m_equ[0]->SetInitialConditions(0.0, false);
141
143 }
144 else
145 {
146 if (m_comm->GetRank() == 0)
147 {
148 out << "\tInital vector : random " << std::endl;
149 }
150
151 NekDouble eps = 0.0001;
152 Vmath::FillWhiteNoise(ntot, eps, &Kseq[1][0], 1);
153
154 auto contfield = std::dynamic_pointer_cast<MultiRegions::ContField>(
155 m_equ[0]->UpdateFields()[0]);
156 if (contfield)
157 {
158 int nGlobal =
159 contfield->GetLocalToGlobalMap()->GetNumGlobalCoeffs();
160 Array<OneD, NekDouble> global(nGlobal), tmp;
161 for (int n = 0; n < m_nfields; ++n)
162 {
163 contfield->LocalToGlobal(Kseq[1] + n * nm, global);
164 contfield->GlobalToLocal(global, tmp = Kseq[1] + n * nm);
165 }
166 }
167 if (m_useMask)
168 {
169 Vmath::Vmul(ntot, Kseq[1], 1, m_maskCoeffs, 1, Kseq[1], 1);
170 }
171 }
172
173 // Perform one iteration to enforce boundary conditions.
174 // Set this as the initial value in the sequence.
175 EV_update(Kseq[1], Kseq[0]);
176 if (m_comm->GetRank() == 0)
177 {
178 out << "Iteration: " << 0 << std::endl;
179 }
180
181 // Normalise first vector in sequence
182 if (m_useMask)
183 {
184 Vmath::Vmul(ntot, Kseq[0], 1, m_maskCoeffs, 1, Kseqcopy[0], 1);
185 }
186 alpha[0] = Blas::Ddot(ntot, &Kseqcopy[0][0], 1, &Kseqcopy[0][0], 1);
187 m_comm->AllReduce(alpha[0], Nektar::LibUtilities::ReduceSum);
188 alpha[0] = std::sqrt(alpha[0]);
189 Vmath::Smul(ntot, 1.0 / alpha[0], Kseq[0], 1, Kseq[0], 1);
190
191 // Fill initial krylov sequence
192 NekDouble resid0;
193 for (i = 1; !converged && i <= m_kdim; ++i)
194 {
195 // Compute next vector
196 EV_update(Kseq[i - 1], Kseq[i]);
197
198 // Normalise
199 if (m_useMask)
200 {
201 Vmath::Vmul(ntot, Kseq[i], 1, m_maskCoeffs, 1, Kseqcopy[i], 1);
202 }
203 alpha[i] = Blas::Ddot(ntot, &Kseqcopy[i][0], 1, &Kseqcopy[i][0], 1);
204 m_comm->AllReduce(alpha[i], Nektar::LibUtilities::ReduceSum);
205 alpha[i] = std::sqrt(alpha[i]);
206
207 // alpha[i] = std::sqrt(alpha[i]);
208 Vmath::Smul(ntot, 1.0 / alpha[i], Kseq[i], 1, Kseq[i], 1);
209
210 // Copy Krylov sequence into temporary storage
211 for (int k = 0; k < i + 1; ++k)
212 {
213 if (m_useMask)
214 {
215 Vmath::Vmul(ntot, Kseq[k], 1, m_maskCoeffs, 1, Tseq[k], 1);
216 Vmath::Vcopy(ntot, Kseq[k], 1, Kseqcopy[k], 1);
217 }
218 else
219 {
220 Vmath::Vcopy(ntot, Kseq[k], 1, Tseq[k], 1);
221 ;
222 }
223 }
224
225 // Generate Hessenberg matrix and compute eigenvalues of it.
226 EV_small(Tseq, Kseqcopy, ntot, alpha, i, zvec, wr, wi, resnorm);
227
228 // Test for convergence.
229 converged = EV_test(i, i, zvec, wr, wi, resnorm, std::min(i, m_nvec),
230 evlout, resid0);
231 if (i >= m_nvec)
232 {
233 converged = std::max(converged, 0);
234 }
235 else
236 {
237 converged = 0;
238 }
239
240 if (m_comm->GetRank() == 0)
241 {
242 out << "Iteration: " << i << " (residual : " << resid0 << ")"
243 << std::endl;
244 }
245 }
246
247 // Continue with full sequence
248 if (!converged)
249 {
250 for (i = m_kdim + 1; !converged && i <= m_nits; ++i)
251 {
252 // Shift all the vectors in the sequence.
253 // First vector is removed.
254 for (int j = 1; j <= m_kdim; ++j)
255 {
256 alpha[j - 1] = alpha[j];
257 Vmath::Vcopy(ntot, Kseq[j], 1, Kseq[j - 1], 1);
258 }
259
260 // Compute next vector
261 EV_update(Kseq[m_kdim - 1], Kseq[m_kdim]);
262
263 // Compute new scale factor
264 if (m_useMask)
265 {
266 Vmath::Vmul(ntot, Kseq[m_kdim], 1, m_maskCoeffs, 1,
267 Kseqcopy[m_kdim], 1);
268 }
269 alpha[m_kdim] = Blas::Ddot(ntot, &Kseqcopy[m_kdim][0], 1,
270 &Kseqcopy[m_kdim][0], 1);
272 alpha[m_kdim] = std::sqrt(alpha[m_kdim]);
273 Vmath::Smul(ntot, 1.0 / alpha[m_kdim], Kseq[m_kdim], 1,
274 Kseq[m_kdim], 1);
275
276 // Copy Krylov sequence into temporary storage
277 for (int k = 0; k < m_kdim + 1; ++k)
278 {
279 if (m_useMask)
280 {
281 Vmath::Vmul(ntot, Kseq[k], 1, m_maskCoeffs, 1, Tseq[k], 1);
282 Vmath::Vcopy(ntot, Kseq[k], 1, Kseqcopy[k], 1);
283 }
284 else
285 {
286 Vmath::Vcopy(ntot, Kseq[k], 1, Tseq[k], 1);
287 ;
288 }
289 }
290
291 // Generate Hessenberg matrix and compute eigenvalues of it
292 EV_small(Tseq, Kseqcopy, ntot, alpha, m_kdim, zvec, wr, wi,
293 resnorm);
294
295 // Test for convergence.
296 converged = EV_test(i, m_kdim, zvec, wr, wi, resnorm, m_nvec,
297 evlout, resid0);
298
299 if (m_comm->GetRank() == 0)
300 {
301 out << "Iteration: " << i << " (residual : " << resid0 << ")"
302 << std::endl;
303 }
304 }
305 }
306
307 m_equ[0]->Output();
308
309 // Evaluate and output computation time and solution accuracy.
310 // The specific format of the error output is essential for the
311 // regression tests to work.
312 // Evaluate L2 Error
313 for (j = 0; j < m_equ[0]->GetNvariables(); ++j)
314 {
315 NekDouble vL2Error = m_equ[0]->L2Error(j, false);
316 NekDouble vLinfError = m_equ[0]->LinfError(j);
317 if (m_comm->GetRank() == 0)
318 {
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;
323 }
324 }
325
326 // Process eigenvectors and write out.
327 m_nvec = converged;
328 EV_post(Tseq, Kseqcopy, ntot, std::min(--i, m_kdim), m_nvec, zvec, wr, wi,
329 converged);
330
331 WARNINGL0(m_imagShift == 0, "Complex Shift applied. "
332 "Need to implement Ritz re-evaluation of"
333 "eigenvalue. Only one half of complex "
334 "value will be correct");
335
336 // store eigenvalues so they can be accessed from driver class
337 m_real_evl = wr;
338 m_imag_evl = wi;
339
340 // Close the runtime info file.
341 if (m_comm->GetRank() == 0)
342 {
343 evlout.close();
344 }
345}
#define WARNINGL0(condition, msg)
Definition: ErrorUtil.hpp:215
Array< OneD, NekDouble > m_maskCoeffs
Definition: DriverArnoldi.h:78
int m_nvec
Dimension of Krylov subspace.
Definition: DriverArnoldi.h:61
int m_nits
Number of vectors to test.
Definition: DriverArnoldi.h:62
Array< OneD, NekDouble > m_imag_evl
Definition: DriverArnoldi.h:75
int m_nfields
interval to dump information if required.
Definition: DriverArnoldi.h:69
Array< OneD, NekDouble > m_real_evl
Operator in solve call is negated.
Definition: DriverArnoldi.h:74
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.
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.
Definition: Vmath.hpp:72
void FillWhiteNoise(int n, const T eps, T *x, const int incx, int outseed)
Fills a vector with white noise.
Definition: Vmath.cpp:154

References Nektar::VarcoeffHashingTest::contfield, Nektar::SolverUtils::DriverArnoldi::CopyFieldToArnoldiArray(), Blas::Ddot(), EV_post(), EV_small(), EV_test(), EV_update(), Vmath::FillWhiteNoise(), Nektar::SolverUtils::Driver::m_comm, Nektar::SolverUtils::Driver::m_equ, Nektar::SolverUtils::DriverArnoldi::m_imag_evl, Nektar::SolverUtils::DriverArnoldi::m_imagShift, Nektar::SolverUtils::DriverArnoldi::m_kdim, Nektar::SolverUtils::DriverArnoldi::m_maskCoeffs, Nektar::SolverUtils::DriverArnoldi::m_nfields, Nektar::SolverUtils::DriverArnoldi::m_nits, Nektar::SolverUtils::DriverArnoldi::m_nvec, Nektar::SolverUtils::DriverArnoldi::m_real_evl, Nektar::SolverUtils::Driver::m_session, Nektar::SolverUtils::DriverArnoldi::m_useMask, Nektar::LibUtilities::ReduceSum, Vmath::Smul(), tinysimd::sqrt(), Vmath::Vcopy(), Vmath::Vmul(), and WARNINGL0.

Referenced by Nektar::SolverUtils::DriverSteadyState::v_Execute().

◆ v_InitObject()

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

Virtual function for initialisation implementation.

Reimplemented from Nektar::SolverUtils::DriverArnoldi.

Reimplemented in Nektar::SolverUtils::DriverSteadyState.

Definition at line 64 of file DriverModifiedArnoldi.cpp.

65{
67
68 m_equ[0]->PrintSummary(out);
69
70 // Print session parameters
71 if (m_comm->GetRank() == 0)
72 {
73 out << "\tArnoldi solver type : Modified Arnoldi" << std::endl;
74 }
75
77
78 for (int i = 0; i < m_nequ; ++i)
79 {
80 m_equ[i]->DoInitialise();
81 }
82
83 // FwdTrans Initial conditions to be in Coefficient Space
84 m_equ[m_nequ - 1]->TransPhysToCoeff();
85}
void v_InitObject(std::ostream &out=std::cout) override
Virtual function for initialisation implementation.
SOLVER_UTILS_EXPORT void ArnoldiSummary(std::ostream &out)
int m_nequ
number of equations
Definition: Driver.h:95

References Nektar::SolverUtils::DriverArnoldi::ArnoldiSummary(), Nektar::SolverUtils::Driver::m_comm, Nektar::SolverUtils::Driver::m_equ, Nektar::SolverUtils::Driver::m_nequ, and Nektar::SolverUtils::DriverArnoldi::v_InitObject().

Referenced by Nektar::SolverUtils::DriverSteadyState::v_InitObject().

Friends And Related Function Documentation

◆ MemoryManager< DriverModifiedArnoldi >

friend class MemoryManager< DriverModifiedArnoldi >
friend

Definition at line 1 of file DriverModifiedArnoldi.h.

Member Data Documentation

◆ className

std::string Nektar::SolverUtils::DriverModifiedArnoldi::className
static
Initial value:
=
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
static DriverSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
DriverFactory & GetDriverFactory()
Definition: Driver.cpp:64

Name of the class.

Definition at line 62 of file DriverModifiedArnoldi.h.

◆ driverLookupId

std::string Nektar::SolverUtils::DriverModifiedArnoldi::driverLookupId
staticprotected
Initial value:
=
0)
static std::string RegisterEnumValue(std::string pEnum, std::string pString, int pEnumValue)
Registers an enumeration value.

Definition at line 78 of file DriverModifiedArnoldi.h.