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

#include <DriverSteadyState.h>

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

Public Member Functions

void PrintSummarySFD ()
 
void ConvergenceHistory (const Array< OneD, const Array< OneD, NekDouble > > &qBar1, const Array< OneD, const Array< OneD, NekDouble > > &q0, NekDouble &MaxNormDiff_q_qBar, NekDouble &MaxNormDiff_q1_q0)
 
void ComputeSFD (const int i, const Array< OneD, const Array< OneD, NekDouble > > &q0, const Array< OneD, const Array< OneD, NekDouble > > &qBar0, Array< OneD, Array< OneD, NekDouble > > &q1, Array< OneD, Array< OneD, NekDouble > > &qBar1)
 
void EvalEV_ScalarSFD (const NekDouble &X_input, const NekDouble &Delta_input, const std::complex< NekDouble > &alpha, NekDouble &MaxEV)
 
void GradientDescentMethod (const std::complex< NekDouble > &alpha, NekDouble &X_output, NekDouble &Delta_output)
 
void ReadEVfile (int &KrylovSubspaceDim, NekDouble &growthEV, NekDouble &frequencyEV)
 
void ComputeOptimization ()
 
void SetSFDOperator (const NekDouble X_input, const NekDouble Delta_input)
 
- Public Member Functions inherited from Nektar::SolverUtils::DriverArnoldi
SOLVER_UTILS_EXPORT void ArnoldiSummary (std::ostream &out)
 
- Public Member Functions inherited from Nektar::SolverUtils::Driver
virtual ~Driver ()
 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 ()
 
SOLVER_UTILS_EXPORT Array< OneD, NekDoubleGetRealEvl (void)
 
SOLVER_UTILS_EXPORT Array< OneD, NekDoubleGetImagEvl (void)
 

Static Public Member Functions

static DriverSharedPtr create (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
 Creates an instance of this class. More...
 
- Static Public Member Functions inherited from Nektar::SolverUtils::DriverModifiedArnoldi
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...
 
- Static Public Attributes inherited from Nektar::SolverUtils::DriverModifiedArnoldi
static std::string className
 Name of the class. More...
 

Protected Member Functions

SOLVER_UTILS_EXPORT DriverSteadyState (const LibUtilities::SessionReaderSharedPtr pSession, const SpatialDomains::MeshGraphSharedPtr pGraph)
 Constructor. More...
 
virtual SOLVER_UTILS_EXPORT ~DriverSteadyState ()
 Destructor. More...
 
virtual SOLVER_UTILS_EXPORT void v_InitObject (std::ostream &out=std::cout)
 Second-stage initialisation. More...
 
virtual SOLVER_UTILS_EXPORT void v_Execute (std::ostream &out=std::cout)
 Virtual function for solve implementation. More...
 
- Protected Member Functions inherited from Nektar::SolverUtils::DriverModifiedArnoldi
 DriverModifiedArnoldi (const LibUtilities::SessionReaderSharedPtr pSession, const SpatialDomains::MeshGraphSharedPtr pGraph)
 Constructor. More...
 
virtual ~DriverModifiedArnoldi ()
 Destructor. More...
 
- Protected Member Functions inherited from Nektar::SolverUtils::DriverArnoldi
 DriverArnoldi (const LibUtilities::SessionReaderSharedPtr pSession, const SpatialDomains::MeshGraphSharedPtr pGraph)
 Constructor. More...
 
virtual ~DriverArnoldi ()
 Destructor. 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)
 
virtual Array< OneD, NekDoublev_GetRealEvl (void)
 
virtual Array< OneD, NekDoublev_GetImagEvl (void)
 
- Protected Member Functions inherited from Nektar::SolverUtils::Driver
 Driver (const LibUtilities::SessionReaderSharedPtr pSession, const SpatialDomains::MeshGraphSharedPtr pGraph)
 Initialises EquationSystem class members. 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 Attributes

int m_stepCounter
 
int m_Check
 
int m_Check_BaseFlow
 
NekDouble TOL
 
int m_infosteps
 
int m_checksteps
 
int NumVar_SFD
 
LibUtilities::Timer timer
 
NekDouble cpuTime
 
NekDouble totalTime
 
NekDouble elapsed
 
std::ofstream m_file
 
NekDouble m_Delta
 Definition of the SFD parameters: More...
 
NekDouble m_X
 
NekDouble m_dt
 
NekDouble M11
 Definition of the SFD operator. More...
 
NekDouble M12
 
NekDouble M21
 
NekDouble M22
 
NekDouble Diff_q_qBar
 
NekDouble Diff_q1_q0
 
bool FlowPartiallyConverged
 For adaptive SFD method. More...
 
NekDouble AdaptiveTOL
 
NekDouble AdaptiveTime
 
int m_NonConvergingStepsCounter
 
NekDouble GrowthRateEV
 
NekDouble FrequencyEV
 

Friends

class MemoryManager< DriverSteadyState >
 

Additional Inherited Members

- 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 iteratiosn. 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
 
- 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
 I the 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 48 of file DriverSteadyState.h.

Constructor & Destructor Documentation

◆ DriverSteadyState()

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

Constructor.

Definition at line 59 of file DriverSteadyState.cpp.

62  : DriverModifiedArnoldi(pSession, pGraph)
63 {
64 }
DriverModifiedArnoldi(const LibUtilities::SessionReaderSharedPtr pSession, const SpatialDomains::MeshGraphSharedPtr pGraph)
Constructor.

◆ ~DriverSteadyState()

Nektar::SolverUtils::DriverSteadyState::~DriverSteadyState ( )
protectedvirtual

Destructor.

Definition at line 70 of file DriverSteadyState.cpp.

71 {
72 }

Member Function Documentation

◆ ComputeOptimization()

void Nektar::SolverUtils::DriverSteadyState::ComputeOptimization ( )

Definition at line 361 of file DriverSteadyState.cpp.

References GradientDescentMethod(), m_Delta, Nektar::SolverUtils::DriverArnoldi::m_kdim, m_X, ReadEVfile(), and SetSFDOperator().

Referenced by v_Execute().

362 {
363  NekDouble growthEV(0.0);
364  NekDouble frequencyEV(0.0);
365 
366  // m_kdim is the dimension of Krylov subspace (defined in the xml file and
367  // used in DriverArnoldi.cpp)
368  ReadEVfile(m_kdim, growthEV, frequencyEV);
369 
370  cout << "\n\tgrowthEV = " << growthEV << endl;
371  cout << "\tfrequencyEV = " << frequencyEV << endl;
372 
373  complex<NekDouble> ApproxEV = polar(exp(growthEV), frequencyEV);
374 
375  NekDouble X_new = m_X;
376  NekDouble Delta_new = m_Delta;
377 
378  GradientDescentMethod(ApproxEV, X_new, Delta_new);
379 
380  m_X = X_new;
381  m_Delta = Delta_new;
382 
384 }
void ReadEVfile(int &KrylovSubspaceDim, NekDouble &growthEV, NekDouble &frequencyEV)
void GradientDescentMethod(const std::complex< NekDouble > &alpha, NekDouble &X_output, NekDouble &Delta_output)
double NekDouble
void SetSFDOperator(const NekDouble X_input, const NekDouble Delta_input)
NekDouble m_Delta
Definition of the SFD parameters:

◆ ComputeSFD()

void Nektar::SolverUtils::DriverSteadyState::ComputeSFD ( const int  i,
const Array< OneD, const Array< OneD, NekDouble > > &  q0,
const Array< OneD, const Array< OneD, NekDouble > > &  qBar0,
Array< OneD, Array< OneD, NekDouble > > &  q1,
Array< OneD, Array< OneD, NekDouble > > &  qBar1 
)

Encapsulated SFD method

Definition at line 341 of file DriverSteadyState.cpp.

References M11, M12, M21, M22, Nektar::SolverUtils::Driver::m_equ, Nektar::SolverUtils::Driver::m_nequ, and Vmath::Svtvp().

Referenced by v_Execute().

346 {
347  q1[i] = Array<OneD, NekDouble> (m_equ[m_nequ - 1]->GetTotPoints(),0.0);
348  qBar1[i] = Array<OneD, NekDouble> (m_equ[m_nequ - 1]->GetTotPoints(),0.0);
349 
350  ///Encapsulated SFD method
351  Vmath::Svtvp(q1[i].num_elements(), M11, q0[i], 1, q1[i], 1, q1[i], 1 );
352  Vmath::Svtvp(q1[i].num_elements(), M12, qBar0[i], 1, q1[i], 1, q1[i], 1 );
353 
354  Vmath::Svtvp(qBar1[i].num_elements(), M21, q0[i], 1, qBar1[i], 1,
355  qBar1[i], 1 );
356  Vmath::Svtvp(qBar1[i].num_elements(), M22, qBar0[i], 1, qBar1[i], 1,
357  qBar1[i], 1 );
358 }
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.cpp:488
Array< OneD, EquationSystemSharedPtr > m_equ
Equation system to solve.
Definition: Driver.h:98
NekDouble M11
Definition of the SFD operator.
int m_nequ
number of equations
Definition: Driver.h:101

◆ ConvergenceHistory()

void Nektar::SolverUtils::DriverSteadyState::ConvergenceHistory ( const Array< OneD, const Array< OneD, NekDouble > > &  qBar1,
const Array< OneD, const Array< OneD, NekDouble > > &  q0,
NekDouble MaxNormDiff_q_qBar,
NekDouble MaxNormDiff_q1_q0 
)

This routine evaluates |q-qBar|_inf (and |q1-q0|_inf) and writes the values in "ConvergenceHistory.txt"

Definition at line 585 of file DriverSteadyState.cpp.

References cpuTime, elapsed, Nektar::SolverUtils::Driver::m_comm, m_Delta, Nektar::SolverUtils::Driver::m_equ, m_file, Nektar::SolverUtils::Driver::m_nequ, m_stepCounter, m_X, NumVar_SFD, Nektar::LibUtilities::Timer::Start(), Nektar::LibUtilities::Timer::Stop(), Nektar::LibUtilities::Timer::TimePerTest(), timer, and totalTime.

Referenced by v_Execute().

590 {
591  Array<OneD, NekDouble > NormDiff_q_qBar(NumVar_SFD, 1.0);
592  Array<OneD, NekDouble > NormDiff_q1_q0(NumVar_SFD, 1.0);
593  MaxNormDiff_q_qBar=0.0;
594  MaxNormDiff_q1_q0=0.0;
595 
596  for(int i = 0; i < NumVar_SFD; ++i)
597  {
598  NormDiff_q_qBar[i] = m_equ[m_nequ - 1]->LinfError(i, qBar1[i]);
599  NormDiff_q1_q0[i] = m_equ[m_nequ - 1]->LinfError(i, q0[i]);
600 
601  if (MaxNormDiff_q_qBar < NormDiff_q_qBar[i])
602  {
603  MaxNormDiff_q_qBar = NormDiff_q_qBar[i];
604  }
605  if (MaxNormDiff_q1_q0 < NormDiff_q1_q0[i])
606  {
607  MaxNormDiff_q1_q0 = NormDiff_q1_q0[i];
608  }
609  }
610 
611  timer.Stop();
613  cpuTime += elapsed;
614  totalTime += elapsed;
615 
616  if (m_comm->GetRank() == 0)
617  {
618  cout << "SFD - Step: " << left << m_stepCounter+1
619  << ";\tTime: " << left << m_equ[m_nequ - 1]->GetFinalTime()
620  << ";\tCPU time = " << left << cpuTime << " s"
621  << ";\tTot time = " << left << totalTime << " s"
622  << ";\tX = " << left << m_X
623  << ";\tDelta = " << left << m_Delta
624  << ";\t|q-qBar|inf = " << left << MaxNormDiff_q_qBar << endl;
625  std::ofstream m_file( "ConvergenceHistory.txt", std::ios_base::app);
626  m_file << m_stepCounter+1 << "\t"
627  << m_equ[m_nequ - 1]->GetFinalTime() << "\t"
628  << totalTime << "\t"
629  << MaxNormDiff_q_qBar << "\t"
630  << MaxNormDiff_q1_q0 << endl;
631  m_file.close();
632  }
633 
634  cpuTime = 0.0;
635  timer.Start();
636 }
NekDouble TimePerTest(unsigned int n)
Returns amount of seconds per iteration in a test with n iterations.
Definition: Timer.cpp:57
Array< OneD, EquationSystemSharedPtr > m_equ
Equation system to solve.
Definition: Driver.h:98
LibUtilities::CommSharedPtr m_comm
Communication object.
Definition: Driver.h:86
NekDouble m_Delta
Definition of the SFD parameters:
int m_nequ
number of equations
Definition: Driver.h:101

◆ create()

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

Creates an instance of this class.

Definition at line 54 of file DriverSteadyState.h.

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

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

◆ EvalEV_ScalarSFD()

void Nektar::SolverUtils::DriverSteadyState::EvalEV_ScalarSFD ( const NekDouble X_input,
const NekDouble Delta_input,
const std::complex< NekDouble > &  alpha,
NekDouble MaxEV 
)

This routine evaluates the maximum eigenvalue of the SFD system when applied to the 1D model u(n+1) = alpha*u(n)

Definition at line 477 of file DriverSteadyState.cpp.

Referenced by GradientDescentMethod().

482 {
483  NekDouble A11 = ( 1.0 + X_input * Delta_input *
484  exp(-(X_input + 1.0/Delta_input)) )/(1.0 + X_input*Delta_input);
485  NekDouble A12 = ( X_input*Delta_input - X_input * Delta_input *
486  exp(-(X_input + 1.0/Delta_input)) )/(1.0 + X_input*Delta_input);
487  NekDouble A21 = ( 1.0 - 1.0 *
488  exp(-(X_input + 1.0/Delta_input)) )/(1 + X_input*Delta_input);
489  NekDouble A22 = ( X_input*Delta_input + 1.0 *
490  exp(-(X_input + 1.0/Delta_input)) )/(1.0 + X_input*Delta_input);
491 
492  complex<NekDouble> B11 = alpha;
493  NekDouble B12 = 0.0;
494  NekDouble B21 = 0.0;
495  NekDouble B22 = 1.0;
496 
497  complex<NekDouble> a = A11*B11 + A12*B21;
498  NekDouble b = A11*B12 + A12*B22;
499  complex<NekDouble> c = A21*B11 + A22*B21;
500  NekDouble d = A21*B12 + A22*B22;
501 
502  complex<NekDouble> delt = sqrt((a-d)*(a-d) + 4.0*b*c);
503  complex<NekDouble> lambda_1 = 0.5*(a+d + delt);
504  complex<NekDouble> lambda_2 = 0.5*(a+d - delt);
505 
506  NekDouble NORM_1 = abs(lambda_1);
507  NekDouble NORM_2 = abs(lambda_2);
508 
509  MaxEV = max(NORM_1, NORM_2);
510 }
double NekDouble

◆ GradientDescentMethod()

void Nektar::SolverUtils::DriverSteadyState::GradientDescentMethod ( const std::complex< NekDouble > &  alpha,
NekDouble X_output,
NekDouble Delta_output 
)

This routine implements a gradient descent method to find the parameters X end Delta which give the minimum eigenlavue of the SFD problem applied to the scalar case u(n+1) = *u(n).

Definition at line 391 of file DriverSteadyState.cpp.

References EvalEV_ScalarSFD().

Referenced by ComputeOptimization(), and v_Execute().

395 {
396  cout << "\n\tWe enter the Gradient Descent Method [...]" << endl;
397  bool OptParmFound = false;
398  bool Descending = true;
399  NekDouble X_input = X_output;
400  NekDouble Delta_input = Delta_output;
401 
402  NekDouble X_init = X_output;
403  NekDouble Delta_init = Delta_output;
404  int stepCounter(0);
405 
406  NekDouble F0(0.0);
407  NekDouble F0x(0.0);
408  NekDouble F0y(0.0);
409  NekDouble F1(0.0);
410  NekDouble dx = 0.00000001;
411  NekDouble dirx(0.0);
412  NekDouble diry(0.0);
413  NekDouble s(0.0);
414  NekDouble CurrentMin = 1.0;
415 
416  while (OptParmFound == false)
417  {
418  Descending = true;
419  EvalEV_ScalarSFD(X_input, Delta_input, alpha, F0);
420  EvalEV_ScalarSFD(X_input + dx, Delta_input, alpha, F0x);
421  EvalEV_ScalarSFD(X_input, Delta_input + dx, alpha, F0y);
422  dirx = (F0 - F0x);
423  diry = (F0 - F0y);
424  s = abs(0.000001/dirx);
425  X_output = X_input + s*dirx;
426  Delta_output = Delta_input + s*diry;
427  F1 = F0;
428 
429  while (Descending == true)
430  {
431  CurrentMin = F1;
432  X_input = X_output;
433  Delta_input = Delta_output;
434  EvalEV_ScalarSFD(X_output, Delta_output, alpha, F1);
435 
436  if (F1 > CurrentMin)
437  {
438  Descending = false;
439  }
440  else
441  {
442  s = s + s*0.01;
443  X_output = X_input + s*dirx;
444  Delta_output = Delta_input + s*diry;
445  }
446 
447  if(stepCounter > 9999999)
448  {
449  //We are stuck in this loop..
450  //Then we restart it with different initail conditions
451  Descending = false;
452  X_input = X_init;
453  Delta_init = Delta_init + Delta_init*0.1;
454  Delta_input = Delta_init;
455  stepCounter = 0;
456  }
457  stepCounter++;
458  }
459 
460  if (abs(F0-F1) < dx)
461  {
462  cout << "\tThe Gradient Descent Method has converged!" << endl;
463  EvalEV_ScalarSFD(X_output, Delta_output, alpha, F1);
464  cout << "\n\tThe updated parameters are: X_tilde = " << X_output
465  << " and Delta_tilde = " << Delta_output << endl;
466  OptParmFound = true;
467  }
468 
469  }
470 }
double NekDouble
void EvalEV_ScalarSFD(const NekDouble &X_input, const NekDouble &Delta_input, const std::complex< NekDouble > &alpha, NekDouble &MaxEV)

◆ PrintSummarySFD()

void Nektar::SolverUtils::DriverSteadyState::PrintSummarySFD ( )

Definition at line 639 of file DriverSteadyState.cpp.

References AdaptiveTime, AdaptiveTOL, Nektar::SolverUtils::eAdaptiveSFD, m_Delta, Nektar::SolverUtils::Driver::m_EvolutionOperator, m_X, and TOL.

Referenced by v_Execute().

640 {
641  cout << "\n====================================="
642  "=================================="
643  << endl;
644  cout << "Parameters for the SFD method:" << endl;
645  cout << "\tControl Coefficient: X = " << m_X << endl;
646  cout << "\tFilter Width: Delta = " << m_Delta << endl;
647  cout << "The simulation is stopped when |q-qBar|inf < " << TOL << endl;
649  {
650  cout << "\nWe use the adaptive SFD method:" << endl;
651  cout << " The parameters are updated every " << AdaptiveTime
652  << " time units;" << endl;
653  cout << " until |q-qBar|inf becomes smaller than " << AdaptiveTOL
654  << endl;
655  }
656  cout << "====================================="
657  "==================================\n" << endl;
658 }
enum EvolutionOperatorType m_EvolutionOperator
Evolution Operator.
Definition: Driver.h:104
NekDouble m_Delta
Definition of the SFD parameters:

◆ ReadEVfile()

void Nektar::SolverUtils::DriverSteadyState::ReadEVfile ( int &  KrylovSubspaceDim,
NekDouble growthEV,
NekDouble frequencyEV 
)

Definition at line 513 of file DriverSteadyState.cpp.

References ASSERTL0, and Nektar::SolverUtils::Driver::m_session.

Referenced by ComputeOptimization().

517 {
518  // This routine reads the .evl file written by the Arnoldi algorithm
519  // (written in September 2014)
520  std::string line;
521  int NumLinesInFile = 0;
522  std::string EVfileName = m_session->GetSessionName() + ".evl";
523  std::ifstream EVfile(EVfileName.c_str());
524  ASSERTL0(EVfile.good(), "Cannot open .evl file.");
525 
526  if(EVfile)
527  {
528  // This block counts the total number of lines of the .evl file
529  // We keep going util we reach the end of the file
530  while(getline(EVfile, line))
531  {
532  ++NumLinesInFile;
533  }
534 
535  // It may happen that the Stability method that have produced the .elv
536  // file converges in less than m_kdim iterations. In this case,
537  // KrylovSubspaceDim has to be changed here
538  if(NumLinesInFile < KrylovSubspaceDim*2.0
539  + KrylovSubspaceDim*(KrylovSubspaceDim+1.0)/2.0)
540  {
541  for(int i = 1; i <= KrylovSubspaceDim; ++i)
542  {
543  if(NumLinesInFile == i*2.0 + i*(i+1.0)/2.0)
544  {
545  KrylovSubspaceDim = i;
546  }
547  }
548  }
549 
550  // go back to the beginning of the file
551  EVfile.clear();
552  EVfile.seekg(0, ios::beg);
553 
554  // We now want to go to the line where the most unstable eigenlavue was
555  // written
556  for(int i = 0; i < (NumLinesInFile - KrylovSubspaceDim); ++i)
557  {
558  std::getline(EVfile, line);
559  cout << "Discard line: " << line << endl;
560  }
561 
562  std::vector<std::string> tokens;
563  std::getline(EVfile, line);
564  boost::algorithm::split(tokens, line,
565  boost::is_any_of("\t "),boost::token_compress_on);
566 
567  ASSERTL0(tokens.size() >= 5,
568  "Unexpected formatting of .evl file while reading line:\n"
569  + line);
570  growthEV = boost::lexical_cast<NekDouble>(tokens[4]);
571  frequencyEV = boost::lexical_cast<NekDouble>(tokens[5]);
572  }
573  else
574  {
575  cout << "An error occurred when opening the .evl file" << endl;
576  }
577  EVfile.close();
578 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
double NekDouble
LibUtilities::SessionReaderSharedPtr m_session
Session reader object.
Definition: Driver.h:89

◆ SetSFDOperator()

void Nektar::SolverUtils::DriverSteadyState::SetSFDOperator ( const NekDouble  X_input,
const NekDouble  Delta_input 
)

This routine defines the encapsulated SFD operator with first-order splitting and exact resolution of the second subproblem. (See http://scitation.aip.org/content/aip/journal/pof2/26/3/10.1063/1.4867482 for details)

Definition at line 328 of file DriverSteadyState.cpp.

References M11, M12, M21, M22, and m_dt.

Referenced by ComputeOptimization(), and v_Execute().

330 {
331  NekDouble X = X_input*m_dt;
332  NekDouble Delta = Delta_input/m_dt;
333  NekDouble coeff = 1.0/(1.0 + X*Delta);
334  M11 = coeff*(1.0 + X*Delta*exp(-(X + 1.0/Delta)));
335  M12 = coeff*(X*Delta*(1.0 - exp(-(X + 1.0/Delta))));
336  M21 = coeff*(1.0 - exp(-(X + 1.0/Delta)));
337  M22 = coeff*(X*Delta + exp(-(X + 1.0/Delta)));
338 }
double NekDouble
NekDouble M11
Definition of the SFD operator.

◆ v_Execute()

void Nektar::SolverUtils::DriverSteadyState::v_Execute ( std::ostream &  out = std::cout)
protectedvirtual

Virtual function for solve implementation.

Definition of variables used in this algorithm

Call the Navier-Stokes solver for one time step

Copy the current flow field into q0

Apply the linear operator to the outcome of the solver

Update qBar

Copy the output of the SFD method into the current flow field

Loop for the adaptive SFD method

We save the final solution into a .fld file

Reimplemented from Nektar::SolverUtils::DriverModifiedArnoldi.

Definition at line 84 of file DriverSteadyState.cpp.

References AdaptiveTime, AdaptiveTOL, ComputeOptimization(), ComputeSFD(), ConvergenceHistory(), cpuTime, Diff_q1_q0, Diff_q_qBar, Nektar::SolverUtils::eAdaptiveSFD, elapsed, FlowPartiallyConverged, FrequencyEV, GradientDescentMethod(), GrowthRateEV, m_Check, m_Check_BaseFlow, m_checksteps, Nektar::SolverUtils::Driver::m_comm, m_Delta, m_dt, Nektar::SolverUtils::Driver::m_equ, Nektar::SolverUtils::Driver::m_EvolutionOperator, m_file, m_infosteps, Nektar::SolverUtils::Driver::m_nequ, m_NonConvergingStepsCounter, Nektar::SolverUtils::Driver::m_session, m_stepCounter, m_X, NumVar_SFD, PrintSummarySFD(), Nektar::LibUtilities::ReduceSum, SetSFDOperator(), Nektar::LibUtilities::Timer::Start(), timer, TOL, totalTime, and Nektar::SolverUtils::DriverModifiedArnoldi::v_Execute().

86 {
87  // With a loop over "DoSolve", this Driver implements the
88  // "encaplulated" Selective Frequency Damping method
89  // to find the steady state of a flow above the critical Reynolds
90  // number.
91  m_equ[m_nequ - 1]->PrintSummary(out);
92 
93  m_session->LoadParameter("IO_InfoSteps", m_infosteps, 1000);
94  m_session->LoadParameter("IO_CheckSteps", m_checksteps, 100000);
95  m_session->LoadParameter("ControlCoeff",m_X, 1);
96  m_session->LoadParameter("FilterWidth", m_Delta, 2);
97 
98  // To evaluate optimum SFD parameters if growth rate provided in the
99  // xml file
100  m_session->LoadParameter("GrowthRateEV", GrowthRateEV, 0.0);
101 
102  // To evaluate optimum SFD parameters if frequency provided in the xml
103  // file
104  m_session->LoadParameter("FrequencyEV", FrequencyEV, 0.0);
105 
106  // Determine when SFD method is converged
107  m_session->LoadParameter("TOL", TOL, 1.0e-08);
108 
109  // Used only for the Adaptive SFD method
110  m_session->LoadParameter("AdaptiveTOL", AdaptiveTOL, 1.0e-02);
111 
112  // Used only for the Adaptive SFD method
113  m_session->LoadParameter("AdaptiveTime", AdaptiveTime, 50.0*m_Delta);
114 
115  if (m_comm->GetRank() == 0)
116  {
117  PrintSummarySFD();
118  }
119 
120  timer.Start();
121 
122  // Definition of shared pointer (used only for the Adaptive SFD method)
123  AdvectionSystemSharedPtr A = m_equ[0]->as<AdvectionSystem>();
124 
125  // Condition necessary to run SFD for the compressible case
126  NumVar_SFD = m_equ[m_nequ - 1]->UpdateFields()[0]->GetCoordim(0);
127  if (m_session->GetSolverInfo("EqType") == "EulerCFE" ||
128  m_session->GetSolverInfo("EqType") == "NavierStokesCFE")
129  {
130  // Number of variables for the compressible equations
131  NumVar_SFD += 2;
132  }
133  if(m_session->DefinesSolverInfo("HOMOGENEOUS"))
134  {
135  if (m_session->GetSolverInfo("HOMOGENEOUS") == "1D")
136  {
137  NumVar_SFD += 1;
138  }
139  }
140 
141  // We store the time step
142  m_dt = m_equ[m_nequ - 1]->GetTimeStep();
143 
144  // Evaluate optimum SFD parameters if dominent EV given by xml file
145  if (GrowthRateEV != 0.0 && FrequencyEV != 0.0)
146  {
147  cout << "Besed on the dominant EV given in the xml file,"
148  << "a 1D model is used to evaluate the optumum parameters"
149  << "of the SFD method:" << endl;
150  complex<NekDouble> EV = polar(exp(GrowthRateEV), FrequencyEV);
152  }
153 
154 
155  // We set up the elements of the operator of the encapsulated
156  // formulation of the selective frequencive damping method
158 
159  // m_steps is set to 1. Then "m_equ[m_nequ - 1]->DoSolve()" will run
160  // for only one time step
161  m_equ[m_nequ - 1]->SetSteps(1);
162  ofstream m_file("ConvergenceHistory.txt", ios::out | ios::trunc);
163 
164  Array<OneD, Array<OneD, NekDouble> > q0(NumVar_SFD);
165  Array<OneD, Array<OneD, NekDouble> > q1(NumVar_SFD);
166  Array<OneD, Array<OneD, NekDouble> > qBar0(NumVar_SFD);
167  Array<OneD, Array<OneD, NekDouble> > qBar1(NumVar_SFD);
168 
169  for(int i = 0; i < NumVar_SFD; ++i)
170  {
171  q0[i] = Array<OneD, NekDouble> (m_equ[m_nequ-1]->GetTotPoints(),
172  0.0); //q0 is initialised
173  qBar0[i] = Array<OneD, NekDouble> (m_equ[m_nequ-1]->GetTotPoints(),
174  0.0);
175  m_equ[m_nequ - 1]->CopyFromPhysField(i, qBar0[i]);
176  }
177 
178  ///Definition of variables used in this algorithm
179  m_stepCounter = 0;
180  m_Check = 0;
181  m_Check_BaseFlow = 1;
183  Diff_q_qBar = 1.0;
184  Diff_q1_q0 = 1.0;
185  cpuTime = 0.0;
186  elapsed = 0.0;
187  totalTime = 0.0;
188  FlowPartiallyConverged = false;
189 
190  while (max(Diff_q_qBar, Diff_q1_q0) > TOL)
191  {
192  ///Call the Navier-Stokes solver for one time step
193  m_equ[m_nequ - 1]->DoSolve();
194 
195  for(int i = 0; i < NumVar_SFD; ++i)
196  {
197  ///Copy the current flow field into q0
198  m_equ[m_nequ - 1]->CopyFromPhysField(i, q0[i]);
199 
200  ///Apply the linear operator to the outcome of the solver
201  ComputeSFD(i, q0, qBar0, q1, qBar1);
202 
203  ///Update qBar
204  qBar0[i] = qBar1[i];
205 
206  ///Copy the output of the SFD method into the current flow field
207  m_equ[m_nequ - 1]->CopyToPhysField(i, q1[i]);
208  }
209 
211  {
213 
214  ///Loop for the adaptive SFD method
216  FlowPartiallyConverged == false)
217  {
218 
220 
221  if (Diff_q_qBar < AdaptiveTOL)
222  {
223  if (m_comm->GetRank() == 0)
224  {
225  cout << "\n\t The SFD method is converging: we compute "
226  << "stability analysis using the 'partially "
227  << "converged' steady state as base flow:\n" << endl;
228  }
229 
230  m_equ[m_nequ - 1]->Checkpoint_BaseFlow(m_Check_BaseFlow);
232 
233  A->GetAdvObject()->SetBaseFlow(q0,m_equ[0]->UpdateFields());
235 
236  if (m_comm->GetRank() == 0)
237  {
238  // Compute the update of the SFD parameters only on
239  // one processor
241  }
242  else
243  {
244  // On all the other processors, the parameters are set
245  // to 0
246  m_X = 0;
247  m_Delta = 0;
248  }
249  // The we give to all the processors the value of X and
250  // Delta of the first processor
253 
254  FlowPartiallyConverged = true;
255  }
257  >= AdaptiveTime)
258  {
259  if (m_comm->GetRank() == 0)
260  {
261  cout << "\n\t We compute stability analysis using"
262  << " the current flow field as base flow:\n"
263  << endl;
264  }
265 
266  m_equ[m_nequ - 1]->Checkpoint_BaseFlow(m_Check_BaseFlow);
268 
269  A->GetAdvObject()->SetBaseFlow(q0,m_equ[0]->UpdateFields());
271 
272  if (m_comm->GetRank() == 0)
273  {
274  // Compute the update of the SFD parameters only on
275  // one processor
277  }
278  else
279  {
280  // On all the other processors, the parameters are set
281  // to 0
282  m_X = 0;
283  m_Delta = 0;
284  }
285  // The we give to all the processors the value of X and
286  // Delta of the first processor
289 
291  }
292  }
293  }
294 
296  {
297  m_Check++;
298  m_equ[m_nequ - 1]->Checkpoint_Output(m_Check);
299  }
300  m_stepCounter++;
301  }
302 
303  m_file.close();
304 
305  ///We save the final solution into a .fld file
306  m_equ[m_nequ - 1]->Output();
307 
308  for(int j = 0; j < m_equ[m_nequ - 1]->GetNvariables(); ++j)
309  {
310  NekDouble vL2Error = m_equ[m_nequ - 1]->L2Error(j,false);
311  NekDouble vLinfError = m_equ[m_nequ - 1]->LinfError(j);
312  if (m_comm->GetRank() == 0)
313  {
314  out << "L 2 error (variable " << m_equ[m_nequ - 1]->GetVariable(j)
315  << ") : " << vL2Error << endl;
316  out << "L inf error (variable " << m_equ[m_nequ - 1]->GetVariable(j)
317  << ") : " << vLinfError << endl;
318  }
319  }
320 }
void ConvergenceHistory(const Array< OneD, const Array< OneD, NekDouble > > &qBar1, const Array< OneD, const Array< OneD, NekDouble > > &q0, NekDouble &MaxNormDiff_q_qBar, NekDouble &MaxNormDiff_q1_q0)
void ComputeSFD(const int i, const Array< OneD, const Array< OneD, NekDouble > > &q0, const Array< OneD, const Array< OneD, NekDouble > > &qBar0, Array< OneD, Array< OneD, NekDouble > > &q1, Array< OneD, Array< OneD, NekDouble > > &qBar1)
std::shared_ptr< AdvectionSystem > AdvectionSystemSharedPtr
Shared pointer to an AdvectionSystem class.
enum EvolutionOperatorType m_EvolutionOperator
Evolution Operator.
Definition: Driver.h:104
void GradientDescentMethod(const std::complex< NekDouble > &alpha, NekDouble &X_output, NekDouble &Delta_output)
double NekDouble
void SetSFDOperator(const NekDouble X_input, const NekDouble Delta_input)
Array< OneD, EquationSystemSharedPtr > m_equ
Equation system to solve.
Definition: Driver.h:98
LibUtilities::CommSharedPtr m_comm
Communication object.
Definition: Driver.h:86
NekDouble m_Delta
Definition of the SFD parameters:
virtual void v_Execute(std::ostream &out=std::cout)
Virtual function for solve implementation.
LibUtilities::SessionReaderSharedPtr m_session
Session reader object.
Definition: Driver.h:89
int m_nequ
number of equations
Definition: Driver.h:101
bool FlowPartiallyConverged
For adaptive SFD method.

◆ v_InitObject()

void Nektar::SolverUtils::DriverSteadyState::v_InitObject ( std::ostream &  out = std::cout)
protectedvirtual

Second-stage initialisation.

Reimplemented from Nektar::SolverUtils::DriverModifiedArnoldi.

Definition at line 78 of file DriverSteadyState.cpp.

References Nektar::SolverUtils::DriverModifiedArnoldi::v_InitObject().

79 {
81 }
virtual void v_InitObject(std::ostream &out=std::cout)
Virtual function for initialisation implementation.

Friends And Related Function Documentation

◆ MemoryManager< DriverSteadyState >

friend class MemoryManager< DriverSteadyState >
friend

Definition at line 51 of file DriverSteadyState.h.

Member Data Documentation

◆ AdaptiveTime

NekDouble Nektar::SolverUtils::DriverSteadyState::AdaptiveTime
private

Definition at line 154 of file DriverSteadyState.h.

Referenced by PrintSummarySFD(), and v_Execute().

◆ AdaptiveTOL

NekDouble Nektar::SolverUtils::DriverSteadyState::AdaptiveTOL
private

Definition at line 153 of file DriverSteadyState.h.

Referenced by PrintSummarySFD(), and v_Execute().

◆ className

string Nektar::SolverUtils::DriverSteadyState::className
static
Initial value:

Name of the class.

Definition at line 65 of file DriverSteadyState.h.

◆ cpuTime

NekDouble Nektar::SolverUtils::DriverSteadyState::cpuTime
private

Definition at line 132 of file DriverSteadyState.h.

Referenced by ConvergenceHistory(), and v_Execute().

◆ Diff_q1_q0

NekDouble Nektar::SolverUtils::DriverSteadyState::Diff_q1_q0
private

Definition at line 149 of file DriverSteadyState.h.

Referenced by v_Execute().

◆ Diff_q_qBar

NekDouble Nektar::SolverUtils::DriverSteadyState::Diff_q_qBar
private

Definition at line 148 of file DriverSteadyState.h.

Referenced by v_Execute().

◆ driverLookupId

string Nektar::SolverUtils::DriverSteadyState::driverLookupId
staticprotected
Initial value:

Definition at line 120 of file DriverSteadyState.h.

◆ elapsed

NekDouble Nektar::SolverUtils::DriverSteadyState::elapsed
private

Definition at line 134 of file DriverSteadyState.h.

Referenced by ConvergenceHistory(), and v_Execute().

◆ FlowPartiallyConverged

bool Nektar::SolverUtils::DriverSteadyState::FlowPartiallyConverged
private

For adaptive SFD method.

Definition at line 152 of file DriverSteadyState.h.

Referenced by v_Execute().

◆ FrequencyEV

NekDouble Nektar::SolverUtils::DriverSteadyState::FrequencyEV
private

Definition at line 157 of file DriverSteadyState.h.

Referenced by v_Execute().

◆ GrowthRateEV

NekDouble Nektar::SolverUtils::DriverSteadyState::GrowthRateEV
private

Definition at line 156 of file DriverSteadyState.h.

Referenced by v_Execute().

◆ M11

NekDouble Nektar::SolverUtils::DriverSteadyState::M11
private

Definition of the SFD operator.

Definition at line 143 of file DriverSteadyState.h.

Referenced by ComputeSFD(), and SetSFDOperator().

◆ M12

NekDouble Nektar::SolverUtils::DriverSteadyState::M12
private

Definition at line 144 of file DriverSteadyState.h.

Referenced by ComputeSFD(), and SetSFDOperator().

◆ M21

NekDouble Nektar::SolverUtils::DriverSteadyState::M21
private

Definition at line 145 of file DriverSteadyState.h.

Referenced by ComputeSFD(), and SetSFDOperator().

◆ M22

NekDouble Nektar::SolverUtils::DriverSteadyState::M22
private

Definition at line 146 of file DriverSteadyState.h.

Referenced by ComputeSFD(), and SetSFDOperator().

◆ m_Check

int Nektar::SolverUtils::DriverSteadyState::m_Check
private

Definition at line 124 of file DriverSteadyState.h.

Referenced by v_Execute().

◆ m_Check_BaseFlow

int Nektar::SolverUtils::DriverSteadyState::m_Check_BaseFlow
private

Definition at line 125 of file DriverSteadyState.h.

Referenced by v_Execute().

◆ m_checksteps

int Nektar::SolverUtils::DriverSteadyState::m_checksteps
private

Definition at line 128 of file DriverSteadyState.h.

Referenced by v_Execute().

◆ m_Delta

NekDouble Nektar::SolverUtils::DriverSteadyState::m_Delta
private

Definition of the SFD parameters:

Definition at line 138 of file DriverSteadyState.h.

Referenced by ComputeOptimization(), ConvergenceHistory(), PrintSummarySFD(), and v_Execute().

◆ m_dt

NekDouble Nektar::SolverUtils::DriverSteadyState::m_dt
private

Definition at line 140 of file DriverSteadyState.h.

Referenced by SetSFDOperator(), and v_Execute().

◆ m_file

std::ofstream Nektar::SolverUtils::DriverSteadyState::m_file
private

Definition at line 135 of file DriverSteadyState.h.

Referenced by ConvergenceHistory(), and v_Execute().

◆ m_infosteps

int Nektar::SolverUtils::DriverSteadyState::m_infosteps
private

Definition at line 127 of file DriverSteadyState.h.

Referenced by v_Execute().

◆ m_NonConvergingStepsCounter

int Nektar::SolverUtils::DriverSteadyState::m_NonConvergingStepsCounter
private

Definition at line 155 of file DriverSteadyState.h.

Referenced by v_Execute().

◆ m_stepCounter

int Nektar::SolverUtils::DriverSteadyState::m_stepCounter
private

Definition at line 123 of file DriverSteadyState.h.

Referenced by ConvergenceHistory(), and v_Execute().

◆ m_X

NekDouble Nektar::SolverUtils::DriverSteadyState::m_X
private

◆ NumVar_SFD

int Nektar::SolverUtils::DriverSteadyState::NumVar_SFD
private

Definition at line 129 of file DriverSteadyState.h.

Referenced by ConvergenceHistory(), and v_Execute().

◆ timer

LibUtilities::Timer Nektar::SolverUtils::DriverSteadyState::timer
private

Definition at line 131 of file DriverSteadyState.h.

Referenced by ConvergenceHistory(), and v_Execute().

◆ TOL

NekDouble Nektar::SolverUtils::DriverSteadyState::TOL
private

Definition at line 126 of file DriverSteadyState.h.

Referenced by PrintSummarySFD(), and v_Execute().

◆ totalTime

NekDouble Nektar::SolverUtils::DriverSteadyState::totalTime
private

Definition at line 133 of file DriverSteadyState.h.

Referenced by ConvergenceHistory(), and v_Execute().