Nektar++
Public Member Functions | Private Attributes | List of all members
Nektar::VortexWaveInteraction Class Reference

#include <VortexWaveInteraction.h>

Public Member Functions

 VortexWaveInteraction (int argc, char *argv[])
 
 ~VortexWaveInteraction (void)
 
void ExecuteRoll (void)
 
void ExecuteStreak (void)
 
void ExecuteWave (void)
 
void ExecuteLoop (bool CalcWaveForce=true)
 
void SaveLoopDetails (std::string dir, int i)
 
void CalcNonLinearWaveForce (void)
 
void CalcL2ToLinfPressure (void)
 
void SaveFile (std::string fileend, std::string dir, int n)
 
void MoveFile (std::string fileend, std::string dir, int n)
 
void CopyFile (std::string file1end, std::string file2end)
 
bool CheckEigIsStationary (bool reset=false)
 
bool CheckIfAtNeutralPoint (void)
 
void UpdateAlpha (int n)
 
void UpdateWaveForceMag (int n)
 
void UpdateDAlphaDWaveForceMag (NekDouble alphainit)
 
void AppendEvlToFile (std::string file, int n)
 
void AppendEvlToFile (std::string file, NekDouble WaveForceMag)
 
int GetIterStart ()
 
int GetIterEnd ()
 
VWIIterationType GetVWIIterationType (void)
 
int GetNOuterIterations (void)
 
int GetMaxOuterIterations (void)
 
NekDouble GetAlpha (void)
 
NekDouble GetAlphaStep (void)
 
NekDouble GetWaveForceMag (void)
 
NekDouble GetWaveForceMagStep (void)
 
NekDouble GetDAlphaDWaveForceMag (void)
 
int GetMaxWaveForceMagIter (void)
 
NekDouble GetEigRelTol (void)
 
int GetMinInnerIterations (void)
 
NekDouble GetPrevAlpha (void)
 
void SetAlpha (NekDouble alpha)
 
void SetWaveForceMag (NekDouble mag)
 
void SetEigRelTol (NekDouble tol)
 
void SetAlphaStep (NekDouble step)
 
void SetMinInnerIterations (int niter)
 
void SetPrevAlpha (NekDouble alpha)
 
bool IfIterInterface (void)
 
Array< OneD, int > GetReflectionIndex (void)
 
void FileRelaxation (int reg)
 

Private Attributes

int m_iterStart
 
int m_iterEnd
 
int m_nOuterIterations
 
int m_maxOuterIterations
 
int m_minInnerIterations
 
int m_maxWaveForceMagIter
 
bool m_deltaFcnApprox
 
bool m_useLinfPressureNorm
 
bool m_moveMeshToCriticalLayer
 
NekDouble m_deltaFcnDecay
 
Array< OneD, NekDoublem_waveForceMag
 
NekDouble m_waveForceMagStep
 
NekDouble m_rollForceScale
 
Array< OneD, NekDoublem_leading_real_evl
 
Array< OneD, NekDoublem_leading_imag_evl
 < Leading real eigenvalue More...
 
Array< OneD, NekDoublem_alpha
 < Leading imaginary eigenvalue More...
 
NekDouble m_alphaStep
 
NekDouble m_neutralPointTol
 
NekDouble m_eigRelTol
 
NekDouble m_vwiRelaxation
 
NekDouble m_dAlphaDWaveForceMag
 
NekDouble m_prevAlpha
 
bool m_iterinterface
 
VWIIterationType m_VWIIterationType
 
Array< OneD, MultiRegions::ExpListSharedPtrm_waveVelocities
 
MultiRegions::ExpListSharedPtr m_wavePressure
 
Array< OneD, Array< OneD, NekDouble > > m_vwiForcing
 
SolverUtils::ForcingProgrammaticSharedPtr m_vwiForcingObj
 
Array< OneD, Array< OneD, NekDouble > > m_bcsForcing
 
Array< OneD, MultiRegions::ExpListSharedPtrm_streakField
 
Array< OneD, MultiRegions::ExpListSharedPtrm_rollField
 
std::string m_sessionName
 
LibUtilities::SessionReaderSharedPtr m_sessionVWI
 
LibUtilities::SessionReaderSharedPtr m_sessionRoll
 
SpatialDomains::MeshGraphSharedPtr m_graphRoll
 
EquationSystemSharedPtr m_solverRoll
 
LibUtilities::SessionReaderSharedPtr m_sessionStreak
 
SpatialDomains::MeshGraphSharedPtr m_graphStreak
 
LibUtilities::SessionReaderSharedPtr m_sessionWave
 
SpatialDomains::MeshGraphSharedPtr m_graphWave
 

Detailed Description

Definition at line 71 of file VortexWaveInteraction.h.

Constructor & Destructor Documentation

◆ VortexWaveInteraction()

Nektar::VortexWaveInteraction::VortexWaveInteraction ( int  argc,
char *  argv[] 
)

Definition at line 46 of file VortexWaveInteraction.cpp.

46  :
48  {
49 
50  int storesize;// number of previous iterations to store
51 
52  m_sessionName = argv[argc-1];
53  string meshfile = m_sessionName + ".xml";
54 
56 
57  std::string VWICondFile(argv[argc-1]);
58  VWICondFile += "_VWI.xml";
59  std::vector<std::string> VWIFilenames;
60  VWIFilenames.push_back(meshfile);
61  VWIFilenames.push_back(VWICondFile);
62 
63  // Create Incompressible NavierStokesSolver session reader.
65  VWIFilenames);
66  m_sessionVWI->InitSession();
67 
68  m_sessionVWI->LoadParameter("AlphaStep", m_alphaStep,0.05);
69  m_sessionVWI->LoadParameter("OuterIterationStoreSize",storesize,10);
70 
71  //set a low tol for interfaceVWI
72  m_sessionVWI->LoadParameter("EigenvalueRelativeTol", m_eigRelTol,1e-3);
73 
74 
75  m_sessionVWI->LoadParameter("NeutralPointTolerance", m_neutralPointTol,1e-4);
76  m_sessionVWI->LoadParameter("MaxOuterIterations", m_maxOuterIterations,100);
77 
78  m_sessionVWI->LoadParameter("StartIteration",m_iterStart, 0);
79  m_sessionVWI->LoadParameter("EndIteration", m_iterEnd, 0);
80 
81  m_sessionVWI->LoadParameter("WaveForceMagStep", m_waveForceMagStep,0.01);
82  m_sessionVWI->LoadParameter("DAlphaDWaveForceMag", m_dAlphaDWaveForceMag,0.0);
83  m_sessionVWI->LoadParameter("MaxWaveForceMagIter",m_maxWaveForceMagIter,1);
84  m_sessionVWI->LoadParameter("RollForceScale", m_rollForceScale,1.0);
85 
86  if(m_sessionVWI->DefinesSolverInfo("DeltaFcnApprox"))
87  {
88  m_deltaFcnApprox = true;
89  m_sessionVWI->LoadParameter("DeltaFcnDecay", m_deltaFcnDecay,1.0/500);
90  }
91  else
92  {
93  m_deltaFcnApprox = false;
94  m_deltaFcnDecay = 0.0;
95  }
96 
97  if(m_sessionVWI->DefinesSolverInfo("LinfPressureNorm"))
98  {
99  m_useLinfPressureNorm = true;
100  }
101  else
102  {
103  m_useLinfPressureNorm = false;
104  }
105 
106  if( m_sessionVWI->DefinesSolverInfo("INTERFACE") )
107  {
108  m_iterinterface = true;
109  }
110  else
111  {
112  m_iterinterface = false;
113  }
114 
115  if(m_sessionVWI->DefinesSolverInfo("MoveMeshToCriticalLayer"))
116  {
118  }
119  else
120  {
122  }
123 
124  m_alpha = Array<OneD, NekDouble> (storesize);
125  m_alpha[0] = m_sessionVWI->GetParameter("Alpha");
126  m_waveForceMag = Array<OneD, NekDouble> (storesize);
127  m_waveForceMag[0] = m_sessionVWI->GetParameter("WaveForceMag");
128 
129  m_leading_real_evl = Array<OneD, NekDouble> (storesize);
130  m_leading_imag_evl = Array<OneD, NekDouble> (storesize);
131 
132  if(m_sessionVWI->DefinesParameter("Relaxation"))
133  {
134  m_vwiRelaxation = m_sessionVWI->GetParameter("Relaxation");
135  // fix minimum number of iterations to be number of
136  // iterations required to make contribution of innitial
137  // forcing to 0.1
139  }
140  else
141  {
142  m_vwiRelaxation = 0.0;
144  }
145 
146  // Initialise NS Roll solver
147  std::string IncCondFile(argv[argc-1]);
148  IncCondFile += "_IncNSCond.xml";
149  std::vector<std::string> IncNSFilenames;
150  IncNSFilenames.push_back(meshfile);
151  IncNSFilenames.push_back(IncCondFile);
152 
153  // Create Incompressible NavierStokesSolver session reader.
154  m_sessionRoll = LibUtilities::SessionReader::CreateInstance(argc, argv, IncNSFilenames,
155  m_sessionVWI->GetComm());
157  std::string vEquation = m_sessionRoll->GetSolverInfo("SolverType");
159  m_solverRoll->PrintSummary(cout);
160 
161 
162  if(m_sessionRoll->DefinesSolverInfo("INTERFACE"))
163  {
164  m_sessionVWI->LoadParameter("EigenvalueRelativeTol", m_eigRelTol,1e-2);
165  }
166 
167  int ncoeffs = m_solverRoll->UpdateFields()[0]->GetNcoeffs();
168  m_vwiForcing = Array<OneD, Array<OneD, NekDouble> > (4);
169  m_vwiForcing[0] = Array<OneD, NekDouble> (4*ncoeffs);
170  for(int i = 1; i < 4; ++i)
171  {
172  m_vwiForcing[i] = m_vwiForcing[i-1] + ncoeffs;
173  }
174 
175  // Fill forcing into m_vwiForcing
176  // Has forcing even been initialised yet?
177 // Vmath::Vcopy(ncoeffs,m_solverRoll->UpdateForces()[0]->GetCoeffs(),1,m_vwiForcing[0],1);
178 // Vmath::Vcopy(ncoeffs,m_solverRoll->UpdateForces()[1]->GetCoeffs(),1,m_vwiForcing[1],1);
179 
180 
181  // Create AdvDiff Streak solver
182  std::string AdvDiffCondFile(argv[argc-1]);
183  AdvDiffCondFile += "_AdvDiffCond.xml";
184  std::vector<std::string> AdvDiffFilenames;
185  AdvDiffFilenames.push_back(meshfile);
186  AdvDiffFilenames.push_back(AdvDiffCondFile);
187 
188  // Create AdvDiffusion session reader.
189  m_sessionStreak = LibUtilities::SessionReader::CreateInstance(argc, argv, AdvDiffFilenames, m_sessionVWI->GetComm());
191 
192  // Initialise LinNS solver
193  std::string LinNSCondFile(argv[argc-1]);
194  LinNSCondFile += "_LinNSCond.xml";
195  std::vector<std::string> LinNSFilenames;
196  LinNSFilenames.push_back(meshfile);
197  LinNSFilenames.push_back(LinNSCondFile);
198 
199  // Create Linearised NS stability session reader.
200  m_sessionWave = LibUtilities::SessionReader::CreateInstance(argc, argv, LinNSFilenames, m_sessionVWI->GetComm());
202 
203  // Set the initial beta value in stability to be equal to VWI file
204  std::string LZstr("LZ");
205  NekDouble LZ = 2*M_PI/m_alpha[0];
206  cout << "Setting LZ in Linearised solver to " << LZ << endl;
207  m_sessionWave->SetParameter(LZstr,LZ);
208 
209  // Check for iteration type
210  if(m_sessionVWI->DefinesSolverInfo("VWIIterationType"))
211  {
212  std::string IterationTypeStr = m_sessionVWI->GetSolverInfo("VWIIterationType");
213  for(int i = 0; i < (int) eVWIIterationTypeSize; ++i)
214  {
215  if(boost::iequals(VWIIterationTypeMap[i],IterationTypeStr))
216  {
218  break;
219  }
220  }
221  }
222  else
223  {
225  }
226 
227 
228 
229  // Check for restart
230  bool restart;
231  m_sessionVWI->MatchSolverInfo("RestartIteration","True",restart,false);
232  if(restart)
233  {
234  switch(m_VWIIterationType)
235  {
236  case eFixedAlpha:
237  break;
238  case eFixedWaveForcing:
239  {
240  FILE *fp;
241  // Check for OuterIter.his file to read
242  if((fp = fopen("OuterIter.his","r")))
243  {
244  char buf[BUFSIZ];
245  std::vector<NekDouble> Alpha, Growth, Phase;
246  NekDouble alpha,growth,phase;
247  while(fgets(buf,BUFSIZ,fp))
248  {
249  sscanf(buf,"%*d:%lf%lf%lf",&alpha,&growth,&phase);
250  Alpha.push_back(alpha);
251  Growth.push_back(growth);
252  Phase.push_back(phase);
253  }
254 
255  m_nOuterIterations = Alpha.size();
256 
257  int nvals = std::min(m_nOuterIterations,(int)m_alpha.size());
258 
259  for(int i = 0; i < nvals; ++i)
260  {
261  m_alpha[nvals-1-i] = Alpha[m_nOuterIterations-nvals+i];
262  m_leading_real_evl[nvals-1-i] = Growth[m_nOuterIterations-nvals+i];
263  m_leading_imag_evl[nvals-1-i] = Phase [m_nOuterIterations-nvals+i];
264  }
265 
267  }
268  else
269  {
270  cout << " No File OuterIter.his to restart from" << endl;
271  }
272  }
273  break;
275  {
276  string nstr = boost::lexical_cast<std::string>(m_iterStart);
277  cout << "Restarting from iteration " << m_iterStart << endl;
278  std::string rstfile = "cp -f Save/" + m_sessionName + ".rst." + nstr + " " + m_sessionName + ".rst";
279  cout << " " << rstfile << endl;
280  if(system(rstfile.c_str()))
281  {
282  ASSERTL0(false,rstfile.c_str());
283  }
284  std::string vwifile = "cp -f Save/" + m_sessionName + ".vwi." + nstr + " " + m_sessionName + ".vwi";
285  cout << " " << vwifile << endl;
286  if(system(vwifile.c_str()))
287  {
288  ASSERTL0(false,vwifile.c_str());
289  }
290  }
291  break;
292  default:
293  ASSERTL0(false,"Unknown VWIITerationType in restart");
294  }
295  }
296 
297  // Check for ConveredSoln to update DAlphaDWaveForce
298  {
299  FILE *fp;
300  if((fp = fopen("ConvergedSolns","r")))
301  {
302  char buf[BUFSIZ];
303  std::vector<NekDouble> WaveForce, Alpha;
304  NekDouble waveforce,alpha;
305  while(fgets(buf,BUFSIZ,fp))
306  {
307  sscanf(buf,"%*d:%lf%lf",&waveforce,&alpha);
308  WaveForce.push_back(waveforce);
309  Alpha.push_back(alpha);
310  }
311 
312  if(Alpha.size() > 1)
313  {
314  int min_i = 0;
315  NekDouble min_alph = fabs(m_alpha[0]-Alpha[min_i]);
316  // find nearest point
317  for(int i = 1; i < Alpha.size(); ++i)
318  {
319  if(fabs(m_alpha[0]-Alpha[min_i]) < min_alph)
320  {
321  min_i = i;
322  min_alph = fabs(m_alpha[0]-Alpha[min_i]);
323  }
324  }
325 
326  // find next nearest point
327  int min_j = (min_i == 0)? 1:0;
328  min_alph = fabs(m_alpha[0]-Alpha[min_j]);
329  for(int i = 0; i < Alpha.size(); ++i)
330  {
331  if(i != min_i)
332  {
333  if(fabs(m_alpha[0]-Alpha[min_j]) < min_alph)
334  {
335  min_j = i;
336  min_alph = fabs(m_alpha[0]-Alpha[min_j]);
337  }
338  }
339  }
340 
341  if(fabs(Alpha[min_i] - Alpha[min_j]) > 1e-4)
342  {
343  m_dAlphaDWaveForceMag = (Alpha[min_i]-Alpha[min_j])/(WaveForce[min_i]-WaveForce[min_j]);
344  }
345  }
346  }
347  }
348  }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:145
static SessionReaderSharedPtr CreateInstance(int argc, char *argv[])
Creates an instance of the SessionReader class.
static MeshGraphSharedPtr Read(const LibUtilities::SessionReaderSharedPtr pSession, LibUtilities::DomainRangeShPtr rng=LibUtilities::NullDomainRangeShPtr, bool fillGraph=true)
Definition: MeshGraph.cpp:113
LibUtilities::SessionReaderSharedPtr m_sessionVWI
LibUtilities::SessionReaderSharedPtr m_sessionRoll
Array< OneD, NekDouble > m_alpha
< Leading imaginary eigenvalue
LibUtilities::SessionReaderSharedPtr m_sessionStreak
Array< OneD, Array< OneD, NekDouble > > m_vwiForcing
SpatialDomains::MeshGraphSharedPtr m_graphStreak
SpatialDomains::MeshGraphSharedPtr m_graphRoll
SpatialDomains::MeshGraphSharedPtr m_graphWave
Array< OneD, NekDouble > m_waveForceMag
LibUtilities::SessionReaderSharedPtr m_sessionWave
Array< OneD, NekDouble > m_leading_real_evl
Array< OneD, NekDouble > m_leading_imag_evl
< Leading real eigenvalue
EquationSystemSharedPtr m_solverRoll
std::shared_ptr< SessionReader > SessionReaderSharedPtr
EquationSystemFactory & GetEquationSystemFactory()
const std::string VWIIterationTypeMap[]
double NekDouble
scalarT< T > log(scalarT< T > in)
Definition: scalar.hpp:278

References ASSERTL0, Nektar::LibUtilities::SessionReader::CreateInstance(), Nektar::LibUtilities::NekFactory< tKey, tBase, tParam >::CreateInstance(), Nektar::eFixedAlpha, Nektar::eFixedAlphaWaveForcing, Nektar::eFixedWaveForcing, Nektar::eVWIIterationTypeSize, Nektar::SolverUtils::GetEquationSystemFactory(), tinysimd::log(), m_alpha, m_alphaStep, m_dAlphaDWaveForceMag, m_deltaFcnApprox, m_deltaFcnDecay, m_eigRelTol, m_graphRoll, m_graphStreak, m_graphWave, m_iterEnd, m_iterinterface, m_iterStart, m_leading_imag_evl, m_leading_real_evl, m_maxOuterIterations, m_maxWaveForceMagIter, m_minInnerIterations, m_moveMeshToCriticalLayer, m_neutralPointTol, m_nOuterIterations, m_rollForceScale, m_sessionName, m_sessionRoll, m_sessionStreak, m_sessionVWI, m_sessionWave, m_solverRoll, m_useLinfPressureNorm, m_vwiForcing, m_VWIIterationType, m_vwiRelaxation, m_waveForceMag, m_waveForceMagStep, Nektar::SpatialDomains::MeshGraph::Read(), UpdateAlpha(), and Nektar::VWIIterationTypeMap.

◆ ~VortexWaveInteraction()

Nektar::VortexWaveInteraction::~VortexWaveInteraction ( void  )

Definition at line 351 of file VortexWaveInteraction.cpp.

352  {
353  m_sessionVWI->Finalise();
354  }

References m_sessionVWI.

Member Function Documentation

◆ AppendEvlToFile() [1/2]

void Nektar::VortexWaveInteraction::AppendEvlToFile ( std::string  file,
int  n 
)

Definition at line 1060 of file VortexWaveInteraction.cpp.

1061  {
1062  FILE *fp;
1063  fp = fopen(file.c_str(),"a");
1064  fprintf(fp, "%d: %lf %16.12le %16.12le\n",n, m_alpha[0], m_leading_real_evl[0],m_leading_imag_evl[0]);
1065  fclose(fp);
1066  }

References m_alpha, m_leading_imag_evl, and m_leading_real_evl.

Referenced by DoFixedForcingIteration(), and main().

◆ AppendEvlToFile() [2/2]

void Nektar::VortexWaveInteraction::AppendEvlToFile ( std::string  file,
NekDouble  WaveForceMag 
)

Definition at line 1068 of file VortexWaveInteraction.cpp.

1069  {
1070  FILE *fp;
1071  fp = fopen(file.c_str(),"a");
1072  fprintf(fp, "%lf %lf %16.12le %16.12le\n",WaveForceMag, m_alpha[0], m_leading_real_evl[0],m_leading_imag_evl[0]);
1073  fclose(fp);
1074  }

References m_alpha, m_leading_imag_evl, and m_leading_real_evl.

◆ CalcL2ToLinfPressure()

void Nektar::VortexWaveInteraction::CalcL2ToLinfPressure ( void  )

Definition at line 968 of file VortexWaveInteraction.cpp.

969  {
970 
971  ExecuteWave();
972 
973  m_wavePressure->GetPlane(0)->BwdTrans(m_wavePressure->GetPlane(0)->GetCoeffs(),
974  m_wavePressure->GetPlane(0)->UpdatePhys());
975  m_wavePressure->GetPlane(1)->BwdTrans(m_wavePressure->GetPlane(1)->GetCoeffs(),
976  m_wavePressure->GetPlane(1)->UpdatePhys());
977 
978  int npts = m_waveVelocities[0]->GetPlane(0)->GetNpoints();
979  NekDouble Linf;
980  Array<OneD, NekDouble> val(2*npts,0.0);
981 
982  Vmath::Vmul(npts,m_wavePressure->GetPlane(0)->UpdatePhys(),1,m_wavePressure->GetPlane(0)->UpdatePhys(),1,val,1);
983  Vmath::Vvtvp(npts,m_wavePressure->GetPlane(1)->UpdatePhys(),1,m_wavePressure->GetPlane(1)->UpdatePhys(),1,val,1,val,1);
984  cout << "int P^2 " << m_wavePressure->GetPlane(0)->Integral(val) << endl;
985  Vmath::Vsqrt(npts,val,1,val,1);
986 
987 
988  Linf = Vmath::Vmax(npts,val,1);
989  cout << "Linf: " << Linf << endl;
990 
991  NekDouble l2,norm;
992  l2 = m_wavePressure->GetPlane(0)->L2(m_wavePressure->GetPlane(0)->GetPhys());
993  norm = l2*l2;
994  l2 = m_wavePressure->GetPlane(1)->L2(m_wavePressure->GetPlane(1)->GetPhys());
995  norm += l2*l2;
996 
997 
998  Vmath::Fill(npts,1.0,val,1);
999  NekDouble area = m_waveVelocities[0]->GetPlane(0)->Integral(val);
1000 
1001  l2 = sqrt(norm/area);
1002 
1003  cout << "L2: " << l2 << endl;
1004 
1005  cout << "Ratio Linf/L2: "<< Linf/l2 << endl;
1006  }
Array< OneD, MultiRegions::ExpListSharedPtr > m_waveVelocities
MultiRegions::ExpListSharedPtr m_wavePressure
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:475
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.cpp:192
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:513
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:45
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
Definition: Vmath.cpp:892
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:267

References ExecuteWave(), Vmath::Fill(), m_wavePressure, m_waveVelocities, tinysimd::sqrt(), Vmath::Vmax(), Vmath::Vmul(), Vmath::Vsqrt(), and Vmath::Vvtvp().

Referenced by main().

◆ CalcNonLinearWaveForce()

void Nektar::VortexWaveInteraction::CalcNonLinearWaveForce ( void  )

Definition at line 563 of file VortexWaveInteraction.cpp.

564  {
565  if(m_sessionRoll->DefinesSolverInfo("INTERFACE") )
566  {
567  static int cnt=0;
568  string wavefile = m_sessionName +".fld";
569  string movedmesh = m_sessionName + "_advPost_moved.xml";
570  string filestreak = m_sessionName + "_streak.fld";
571  char c[16]="";
572  sprintf(c,"%d",cnt);
573  char c_alpha[16]="";
574  sprintf(c_alpha,"%f",m_alpha[0]);
575  string syscall;
576  if( m_sessionVWI->GetSolverInfo("INTERFACE")=="phase" )
577  {
578  string filePost = m_sessionName + "_advPost.xml";
579  syscall = "../../utilities/PostProcessing/Extras/FldCalcBCs "
580  + filePost +" "+
581  "meshhalf_pos_Spen_stability_moved.fld meshhalf_pos_Spen_advPost_moved.fld "
582  +c_alpha +" > data_alpha0";
583  cout<<syscall.c_str()<<endl;
584  if(system(syscall.c_str()))
585  {
586  ASSERTL0(false,syscall.c_str());
587  }
588 
589  syscall = "cp -f meshhalf_pos_Spen_stability_moved_u_5.bc "+m_sessionName+"_u_5.bc";
590  cout<<syscall.c_str()<<endl;
591  if(system(syscall.c_str()))
592  {
593  ASSERTL0(false,syscall.c_str());
594  }
595  syscall = "cp -f meshhalf_pos_Spen_stability_moved_v_5.bc "+m_sessionName+"_v_5.bc";
596  cout<<syscall.c_str()<<endl;
597  if(system(syscall.c_str()))
598  {
599  ASSERTL0(false,syscall.c_str());
600  }
601  }
602  else
603  {
604  syscall = "../../utilities/PostProcessing/Extras/FldCalcBCs "
605  + movedmesh + " " + wavefile + " " + filestreak + " "+c_alpha +" > datasub_"+c;
606  cout<<syscall.c_str()<<endl;
607  if(system(syscall.c_str()))
608  {
609  ASSERTL0(false,syscall.c_str());
610  }
611  }
612 
613 
614 
615  //relaxation for different alpha values? does it make sense?
616 
617  //save the wave
618  string wave_subalp = m_sessionName + "_wave_subalp_"+c+".fld";
619  syscall = "cp -f " + wavefile + " " + wave_subalp;
620  cout<<syscall.c_str()<<endl;
621  if(system(syscall.c_str()))
622  {
623  ASSERTL0(false,syscall.c_str());
624  }
625  //FileRelaxation(3);
626  cnt++;
627 
628  }
629  else
630  {
631  int npts = m_waveVelocities[0]->GetPlane(0)->GetNpoints();
632  int ncoeffs = m_waveVelocities[0]->GetPlane(0)->GetNcoeffs();
633  Array<OneD, NekDouble> val(npts), der1(2*npts);
634  Array<OneD, NekDouble> der2 = der1 + npts;
635  Array<OneD, NekDouble> streak;
636  static int projectfield = -1;
637 
638  if(m_deltaFcnApprox)
639  {
640  streak = Array<OneD, NekDouble> (npts);
641  m_streakField[0]->BwdTrans(m_streakField[0]->GetCoeffs(), streak);
642  }
643 
644  // Set project field to be first field that has a Neumann
645  // boundary since this not impose any condition on the vertical boundaries
646  // Othersise set to zero.
647  if(projectfield == -1)
648  {
649  Array<OneD, const SpatialDomains::BoundaryConditionShPtr > BndConds;
650 
651  for(int i = 0; i < m_waveVelocities.size(); ++i)
652  {
653  BndConds = m_waveVelocities[i]->GetBndConditions();
654  for(int j = 0; j < BndConds.size(); ++j)
655  {
656  if(BndConds[j]->GetBoundaryConditionType() == SpatialDomains::eNeumann)
657  {
658  projectfield = i;
659  break;
660  }
661  }
662  if(projectfield != -1)
663  {
664  break;
665  }
666  }
667  if(projectfield == -1)
668  {
669  cout << "using first field to project non-linear forcing which imposes a Dirichlet condition" << endl;
670  projectfield = 0;
671  }
672  }
673 
674  // Shift m_vwiForcing in case of relaxation
675  Vmath::Vcopy(ncoeffs,m_vwiForcing[0],1,m_vwiForcing[2],1);
676  Vmath::Vcopy(ncoeffs,m_vwiForcing[1],1,m_vwiForcing[3],1);
677 
678  // determine inverse of area normalised field.
679  m_wavePressure->GetPlane(0)->BwdTrans(m_wavePressure->GetPlane(0)->GetCoeffs(), m_wavePressure->GetPlane(0)->UpdatePhys());
680  m_wavePressure->GetPlane(1)->BwdTrans(m_wavePressure->GetPlane(1)->GetCoeffs(), m_wavePressure->GetPlane(1)->UpdatePhys());
681  NekDouble invnorm;
682 
684  {
685  Vmath::Vmul(npts,m_wavePressure->GetPlane(0)->UpdatePhys(),1,m_wavePressure->GetPlane(0)->UpdatePhys(),1,der1,1);
686  Vmath::Vvtvp(npts,m_wavePressure->GetPlane(1)->UpdatePhys(),1,m_wavePressure->GetPlane(1)->UpdatePhys(),1,der1,1,der1,1);
687  Vmath::Vsqrt(npts,der1,1,der1,1);
688 
689  NekDouble Linf = Vmath::Vmax(npts,der1,1);
690 
691  invnorm = 1.0/Linf;
692  }
693  else
694  {
695  // Determine normalisation of pressure so that |P|/A = 1
696  NekDouble l2;
697  l2 = m_wavePressure->GetPlane(0)->L2(m_wavePressure->GetPlane(0)->GetPhys());
698  invnorm = l2*l2;
699  l2 = m_wavePressure->GetPlane(1)->L2(m_wavePressure->GetPlane(1)->GetPhys());
700  invnorm += l2*l2;
701  Vmath::Fill(2*npts,1.0,der1,1);
702  NekDouble area = m_waveVelocities[0]->GetPlane(0)->Integral(der1);
703  cout << "Area: " << area << endl;
704  invnorm = sqrt(area/invnorm);
705  }
706 
707  // Get hold of arrays.
708  m_waveVelocities[0]->GetPlane(0)->BwdTrans(m_waveVelocities[0]->GetPlane(0)->GetCoeffs(),m_waveVelocities[0]->GetPlane(0)->UpdatePhys());
709  Array<OneD, NekDouble> u_real = m_waveVelocities[0]->GetPlane(0)->UpdatePhys();
710  Vmath::Smul(npts,invnorm,u_real,1,u_real,1);
711  m_waveVelocities[0]->GetPlane(1)->BwdTrans(m_waveVelocities[0]->GetPlane(1)->GetCoeffs(),m_waveVelocities[0]->GetPlane(1)->UpdatePhys());
712  Array<OneD, NekDouble> u_imag = m_waveVelocities[0]->GetPlane(1)->UpdatePhys();
713  Vmath::Smul(npts,invnorm,u_imag,1,u_imag,1);
714  m_waveVelocities[1]->GetPlane(0)->BwdTrans(m_waveVelocities[1]->GetPlane(0)->GetCoeffs(),m_waveVelocities[1]->GetPlane(0)->UpdatePhys());
715  Array<OneD, NekDouble> v_real = m_waveVelocities[1]->GetPlane(0)->UpdatePhys();
716  Vmath::Smul(npts,invnorm,v_real,1,v_real,1);
717  m_waveVelocities[1]->GetPlane(1)->BwdTrans(m_waveVelocities[1]->GetPlane(1)->GetCoeffs(),m_waveVelocities[1]->GetPlane(1)->UpdatePhys());
718  Array<OneD, NekDouble> v_imag = m_waveVelocities[1]->GetPlane(1)->UpdatePhys();
719  Vmath::Smul(npts,invnorm,v_imag,1,v_imag,1);
720 
721  // normalise wave for output
722  Vmath::Smul(2*ncoeffs,invnorm,m_waveVelocities[0]->UpdateCoeffs(),1,m_waveVelocities[0]->UpdateCoeffs(),1);
723  Vmath::Smul(2*ncoeffs,invnorm,m_waveVelocities[1]->UpdateCoeffs(),1,m_waveVelocities[1]->UpdateCoeffs(),1);
724  Vmath::Smul(2*ncoeffs,invnorm,m_waveVelocities[2]->UpdateCoeffs(),1,m_waveVelocities[2]->UpdateCoeffs(),1);
725 
726  // dump field
727  {
728  std::vector<std::string> variables(3);
729  std::vector<Array<OneD, NekDouble> > outfield(3);
730  variables[0] = "u_w";
731  variables[1] = "v_w";
732  variables[2] = "w_w";
733  outfield[0] = m_waveVelocities[0]->UpdateCoeffs();
734  outfield[1] = m_waveVelocities[1]->UpdateCoeffs();
735  outfield[2] = m_waveVelocities[2]->UpdateCoeffs();
736  std::string outname = m_sessionName + "_wave.fld";
737  m_solverRoll->WriteFld(outname, m_waveVelocities[0], outfield, variables);
738  }
739 
740 #if 1
741  int ncoeffs_p = m_wavePressure->GetPlane(0)->GetNcoeffs();
742  Vmath::Smul(ncoeffs_p,invnorm,m_wavePressure->GetPlane(0)->UpdateCoeffs(),1,m_wavePressure->GetPlane(0)->UpdateCoeffs(),1);
743  Vmath::Smul(ncoeffs_p,invnorm,m_wavePressure->GetPlane(1)->UpdateCoeffs(),1,m_wavePressure->GetPlane(1)->UpdateCoeffs(),1);
744 #else
745  m_wavePressure->GetPlane(0)->BwdTrans(m_wavePressure->GetPlane(0)->GetCoeffs(),m_wavePressure->GetPlane(0)->UpdatePhys());
746  Vmath::Smul(npts,invnorm,m_wavePressure->GetPlane(0)->UpdatePhys(),1,m_wavePressure->GetPlane(0)->UpdatePhys(),1);
747  m_wavePressure->GetPlane(0)->FwdTrans(m_wavePressure->GetPlane(0)->UpdatePhys(),m_wavePressure->GetPlane(0)->UpdateCoeffs());
748 
749  m_wavePressure->GetPlane(1)->BwdTrans(m_wavePressure->GetPlane(1)->GetCoeffs(),m_wavePressure->GetPlane(1)->UpdatePhys());
750  Vmath::Smul(npts,invnorm,m_wavePressure->GetPlane(1)->UpdatePhys(),1,m_wavePressure->GetPlane(1)->UpdatePhys(),1);
751  m_wavePressure->GetPlane(1)->FwdTrans(m_wavePressure->GetPlane(1)->UpdatePhys(),m_wavePressure->GetPlane(1)->UpdateCoeffs());
752 #endif
753 
754  // Calculate non-linear terms for x and y directions
755  // d/dx(u u* + u* u)
756  Vmath::Vmul (npts,u_real,1,u_real,1,val,1);
757  Vmath::Vvtvp(npts,u_imag,1,u_imag,1,val,1,val,1);
758  Vmath::Smul (npts,2.0,val,1,val,1);
759  m_waveVelocities[0]->GetPlane(0)->PhysDeriv(0,val,der1);
760 
761  // d/dy(v u* + v* u)
762  Vmath::Vmul (npts,u_real,1,v_real,1,val,1);
763  Vmath::Vvtvp(npts,u_imag,1,v_imag,1,val,1,val,1);
764  Vmath::Smul (npts,2.0,val,1,val,1);
765  m_waveVelocities[0]->GetPlane(0)->PhysDeriv(1,val,der2);
766 
767  Vmath::Vadd(npts,der1,1,der2,1,der1,1);
768 #if 1
769  m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(der1,m_vwiForcing[0]);
770 #else
771  m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(der1,m_vwiForcing[0]);
772 #endif
773  Vmath::Smul(ncoeffs,-m_waveForceMag[0],m_vwiForcing[0],1,m_vwiForcing[0],1);
774 
775  // d/dx(u v* + u* v)
776  m_waveVelocities[0]->GetPlane(0)->PhysDeriv(0,val,der1);
777 
778  // d/dy(v v* + v* v)
779  Vmath::Vmul(npts,v_real,1,v_real,1,val,1);
780  Vmath::Vvtvp(npts,v_imag,1,v_imag,1,val,1,val,1);
781  Vmath::Smul (npts,2.0,val,1,val,1);
782  m_waveVelocities[0]->GetPlane(0)->PhysDeriv(1,val,der2);
783 
784  Vmath::Vadd(npts,der1,1,der2,1,der1,1);
785 
786 #if 1
787  m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(der1,m_vwiForcing[1]);
788 #else
789  m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(der1,m_vwiForcing[1]);
790 #endif
791 
792  Vmath::Smul(ncoeffs,-m_waveForceMag[0],m_vwiForcing[1],1,m_vwiForcing[1],1);
793 
794  //by default the symmetrization is on
795  bool symm=true;
796  m_sessionVWI->MatchSolverInfo("Symmetrization","True",symm,true);
797 #if 0
798  if(symm== true )
799  {
800 
801  // Symmetrise forcing
802  //-> Get coordinates
803  Array<OneD, NekDouble> coord(2);
804  Array<OneD, NekDouble> coord_x(npts);
805  Array<OneD, NekDouble> coord_y(npts);
806 
807  //-> Impose symmetry (x -> -x + Lx/2, y-> -y) on coordinates
808  m_waveVelocities[0]->GetPlane(0)->GetCoords(coord_x,coord_y);
809  NekDouble xmax = Vmath::Vmax(npts,coord_x,1);
810  Vmath::Neg(npts,coord_x,1);
811  Vmath::Sadd(npts,xmax,coord_x,1,coord_x,1);
812  Vmath::Neg(npts,coord_y,1);
813 
814  int i, physoffset;
815 
816  //-> Obtain list of expansion element ids for each point.
817  Array<OneD, int> Eid(npts);
818  // This search may not be necessary every iteration
819  for(i = 0; i < npts; ++i)
820  {
821  coord[0] = coord_x[i];
822  coord[1] = coord_y[i];
823 
824  // Note this will not quite be symmetric.
825  Eid[i] = m_waveVelocities[0]->GetPlane(0)->GetExpIndex(coord,1e-6);
826  }
827 
828  // Interpolate field 0
829  m_waveVelocities[0]->GetPlane(0)->BwdTrans_IterPerExp(m_vwiForcing[0],der1);
830  for(i = 0; i < npts; ++i)
831  {
832  physoffset = m_waveVelocities[0]->GetPlane(0)->GetPhys_Offset(Eid[i]);
833  coord[0] = coord_x[i];
834  coord[1] = coord_y[i];
835  der2 [i] = m_waveVelocities[0]->GetPlane(0)->GetExp(Eid[i])->PhysEvaluate(coord, der1 + physoffset);
836  }
837  //-> Average field 0
838  Vmath::Vsub(npts,der1,1,der2,1,der2,1);
839  Vmath::Smul(npts,0.5,der2,1,der2,1);
840 #if 1
841  m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(der2,m_vwiForcing[0]);
842 #else
843  m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(der2, m_vwiForcing[0]);
844 #endif
845 
846  //-> Interpoloate field 1
847  m_waveVelocities[0]->GetPlane(0)->BwdTrans_IterPerExp(m_vwiForcing[1],der1);
848  for(i = 0; i < npts; ++i)
849  {
850  physoffset = m_waveVelocities[0]->GetPlane(0)->GetPhys_Offset(Eid[i]);
851  coord[0] = coord_x[i];
852  coord[1] = coord_y[i];
853  der2[i] = m_waveVelocities[0]->GetPlane(0)->GetExp(Eid[i])->PhysEvaluate(coord, der1 + physoffset);
854  }
855 
856  //-> Average field 1
857  Vmath::Vsub(npts,der1,1,der2,1,der2,1);
858  Vmath::Smul(npts,0.5,der2,1,der2,1);
859 #if 1
860  m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(der2,m_vwiForcing[1]);
861 #else
862  m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(der2, m_vwiForceing[1]);
863 #endif
864  }
865 #else
866  int i;
867  if(symm== true )
868  {
869  cout<<"symmetrization is active"<<endl;
870  static Array<OneD, int> index = GetReflectionIndex();
871 
872  m_waveVelocities[0]->GetPlane(0)->BwdTrans_IterPerExp(m_vwiForcing[0],der1);
873  for(i = 0; i < npts; ++i)
874  {
875  if(index[i] != -1)
876  {
877  val[i] = 0.5*(der1[i] - der1[index[i]]);
878  }
879  }
880 #if 1
881  m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(val,m_vwiForcing[0]);
882 #else
883  m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(val, m_vwiForcing[0]);
884 #endif
885 
886  m_waveVelocities[0]->GetPlane(0)->BwdTrans_IterPerExp(m_vwiForcing[1],der2);
887  for(i = 0; i < npts; ++i)
888  {
889  if(index[i] != -1)
890  {
891  val[i] = 0.5*(der2[i] - der2[index[i]]);
892  }
893  }
894 #if 1
895  m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(val,m_vwiForcing[1]);
896 #else
897  m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(val, m_vwiForcing[1]);
898 #endif
899  }
900 
901 
902  Vmath::Vmul(npts,der1,1,der1,1,val,1);
903  Vmath::Vvtvp(npts,der2,1,der2,1,val,1,val,1);
904  Vmath::Vsqrt(npts,val,1,val,1);
905  cout << "F_Linf: " << Vmath::Vmax(npts,val,1) << endl;
906 
907 #endif
908 
909  if(m_vwiRelaxation)
910  {
911  Vmath::Smul(ncoeffs,1.0-m_vwiRelaxation,
912  m_vwiForcing[0],1,m_vwiForcing[0],1);
914  m_vwiForcing[0],1,m_vwiForcing[0],1);
915 
916  Vmath::Smul(ncoeffs,1.0-m_vwiRelaxation,
917  m_vwiForcing[1],1,m_vwiForcing[1],1);
919  m_vwiForcing[1],1,m_vwiForcing[1],1);
920  }
921 
922  if(m_deltaFcnApprox)
923  {
924  for(int j = 0; j < 2; ++j)
925  {
926 
927  m_waveVelocities[projectfield]->GetPlane(0)->BwdTrans(m_vwiForcing[j],der1);
928  for(int i = 0; i < npts; ++i)
929  {
930  der1[i] *= exp(-streak[i]*streak[i]/m_deltaFcnDecay);
931  }
932  m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(der1,m_vwiForcing[j]);
933  }
934  }
935 
936 
937  // dump output
938  std::vector<std::string> variables(4);
939  std::vector<Array<OneD, NekDouble> > outfield(4);
940  variables[0] = "u"; variables[1] = "v";
941  variables[2] = "pr"; variables[3] = "pi";
942  outfield[0] = m_vwiForcing[0];
943  outfield[1] = m_vwiForcing[1];
944  Array<OneD,NekDouble> soln(npts,0.0);
945  m_wavePressure->GetPlane(0)->BwdTrans(m_wavePressure->GetPlane(0)->GetCoeffs(),m_wavePressure->GetPlane(0)->UpdatePhys());
946  outfield[2] = Array<OneD,NekDouble>(ncoeffs);
947  m_waveVelocities[0]->GetPlane(0)->FwdTrans_IterPerExp(m_wavePressure->GetPlane(0)->GetPhys(),outfield[2]);
948  m_wavePressure->GetPlane(1)->BwdTrans(m_wavePressure->GetPlane(1)->GetCoeffs(),m_wavePressure->GetPlane(1)->UpdatePhys());
949 
950  Vmath::Vmul(npts,m_wavePressure->GetPlane(0)->UpdatePhys(),1,
951  m_wavePressure->GetPlane(0)->UpdatePhys(),1,val,1);
952  Vmath::Vvtvp(npts,m_wavePressure->GetPlane(1)->UpdatePhys(),1,
953  m_wavePressure->GetPlane(1)->UpdatePhys(),1,val,1,val,1);
954  cout << "int P^2: " << m_wavePressure->GetPlane(0)->Integral(val) << endl;
955  Vmath::Vsqrt(npts,val,1,val,1);
956  cout << "PLinf: " << Vmath::Vmax(npts,val,1) << endl;
957 
958  outfield[3] = Array<OneD,NekDouble>(ncoeffs);
959  m_waveVelocities[1]->GetPlane(0)->FwdTrans_IterPerExp(m_wavePressure->GetPlane(1)->GetPhys(),outfield[3]);
960 
961  std::string outname = m_sessionName + ".vwi";
962 
963  m_solverRoll->WriteFld(outname, m_waveVelocities[0]->GetPlane(0), outfield, variables);
964 
965  }
966  }
Array< OneD, MultiRegions::ExpListSharedPtr > m_streakField
Array< OneD, int > GetReflectionIndex(void)
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:565
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:461
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:322
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.cpp:225
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha - x.
Definition: Vmath.cpp:341
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1199
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:372

References ASSERTL0, Nektar::SpatialDomains::eNeumann, Vmath::Fill(), GetReflectionIndex(), m_alpha, m_deltaFcnApprox, m_deltaFcnDecay, m_sessionName, m_sessionRoll, m_sessionVWI, m_solverRoll, m_streakField, m_useLinfPressureNorm, m_vwiForcing, m_vwiRelaxation, m_waveForceMag, m_wavePressure, m_waveVelocities, Vmath::Neg(), Vmath::Sadd(), Vmath::Smul(), tinysimd::sqrt(), Vmath::Svtvp(), Vmath::Vadd(), Vmath::Vcopy(), Vmath::Vmax(), Vmath::Vmul(), Vmath::Vsqrt(), Vmath::Vsub(), and Vmath::Vvtvp().

Referenced by DoFixedForcingIteration(), ExecuteLoop(), and main().

◆ CheckEigIsStationary()

bool Nektar::VortexWaveInteraction::CheckEigIsStationary ( bool  reset = false)

Definition at line 1603 of file VortexWaveInteraction.cpp.

1604  {
1605  static NekDouble previous_real_evl = -1.0;
1606  static NekDouble previous_imag_evl = -1.0;
1607  static int min_iter = 0;
1608 
1609  if(reset)
1610  {
1611  previous_real_evl = -1.0;
1612  min_iter = 0;
1613  }
1614 
1615  if(previous_real_evl == -1.0 || min_iter < m_minInnerIterations)
1616  {
1617  previous_real_evl = m_leading_real_evl[0];
1618  previous_imag_evl = m_leading_imag_evl[0];
1619  min_iter++;
1620  return false;
1621  }
1622 
1623  cout << "Growth tolerance: " << fabs((m_leading_real_evl[0] - previous_real_evl)/m_leading_real_evl[0]) << endl;
1624  cout << "Phase tolerance: " << fabs((m_leading_imag_evl[0] - previous_imag_evl)/m_leading_imag_evl[0]) << endl;
1625 
1626  // See if real and imaginary growth have converged to with m_eigRelTol
1627  if((fabs((m_leading_real_evl[0] - previous_real_evl)/m_leading_real_evl[0]) < m_eigRelTol)||(fabs(m_leading_real_evl[0]) < 1e-6))
1628  {
1629  previous_real_evl = m_leading_real_evl[0];
1630  previous_imag_evl = m_leading_imag_evl[0];
1631  if((fabs((m_leading_imag_evl[0] - previous_imag_evl)/m_leading_imag_evl[0]) < m_eigRelTol)||(fabs(m_leading_imag_evl[0]) < 1e-6))
1632  {
1633  return true;
1634  }
1635  }
1636  else
1637  {
1638  if(fabs(m_leading_imag_evl[0]) > 1e-2)
1639  {
1640  cout << "Warning: imaginary eigenvalue is greater than 1e-2" << endl;
1641  }
1642  previous_real_evl = m_leading_real_evl[0];
1643  previous_imag_evl = m_leading_imag_evl[0];
1644  return false;
1645  }
1646  return false;
1647  }

References m_eigRelTol, m_leading_imag_evl, m_leading_real_evl, and m_minInnerIterations.

Referenced by DoFixedForcingIteration().

◆ CheckIfAtNeutralPoint()

bool Nektar::VortexWaveInteraction::CheckIfAtNeutralPoint ( void  )

Definition at line 1651 of file VortexWaveInteraction.cpp.

1652  {
1653  bool returnval = false;
1654  if(m_sessionRoll->DefinesSolverInfo("INTERFACE"))
1655  {
1656  if(m_sessionVWI->GetSolverInfo("INTERFACE")=="phase")
1657  {
1658  if( abs(m_leading_real_evl[0]) < 1e-4 )
1659  {
1660  returnval = true;
1661  }
1662  }
1663  else
1664  {
1665  if( abs(m_leading_real_evl[0]) < 1e-4 && abs(m_leading_imag_evl[0]) <2e-6 )
1666  {
1667  returnval = true;
1668  }
1669  }
1670 
1671  }
1672  else
1673  {
1677  {
1678  returnval = true;
1679  }
1680  }
1681  if(fabs(m_leading_imag_evl[0]) > 1e-2)
1682  {
1683  cout << "Warning: imaginary eigenvalue is greater than 1e-2" << endl;
1684  }
1685 
1686  return returnval;
1687  }
scalarT< T > abs(scalarT< T > in)
Definition: scalar.hpp:272

References tinysimd::abs(), m_leading_imag_evl, m_leading_real_evl, m_neutralPointTol, m_sessionRoll, and m_sessionVWI.

Referenced by DoFixedForcingIteration().

◆ CopyFile()

void Nektar::VortexWaveInteraction::CopyFile ( std::string  file1end,
std::string  file2end 
)

Definition at line 1048 of file VortexWaveInteraction.cpp.

1049  {
1050  string cpfile1 = m_sessionName + file1end;
1051  string cpfile2 = m_sessionName + file2end;
1052  string syscall = "cp -f " + cpfile1 + " " + cpfile2;
1053 
1054  if(system(syscall.c_str()))
1055  {
1056  ASSERTL0(false,syscall.c_str());
1057  }
1058  }

References ASSERTL0, and m_sessionName.

Referenced by ExecuteRoll(), ExecuteStreak(), and ExecuteWave().

◆ ExecuteLoop()

void Nektar::VortexWaveInteraction::ExecuteLoop ( bool  CalcWaveForce = true)

Definition at line 1096 of file VortexWaveInteraction.cpp.

1097  {
1098 
1099 
1100  //the global info has to be written in the
1101  //roll session file to use the interface loop
1102  if(m_sessionRoll->DefinesSolverInfo("INTERFACE"))
1103  {
1104  static int cnt=0;
1105  bool skiprollstreak=false;
1106  if(cnt==0 && m_sessionVWI->GetParameter("rollstreakfromit")==1)
1107  {
1108  skiprollstreak =true;
1109  cout<<"skip roll-streak at the first iteration"<<endl;
1110  }
1111 
1112 
1113  if(skiprollstreak != true)
1114  {
1115 
1117  ExecuteRoll();
1119 #ifndef _WIN32
1120  sleep(3);
1121 #endif
1122  ExecuteStreak();
1123 #ifndef _WIN32
1124  sleep(3);
1125 #endif
1126  }
1127 
1128  string syscall;
1129  char c[16]="";
1130  string movedmesh = m_sessionName + "_advPost_moved.xml";
1131  string movedinterpmesh = m_sessionName + "_interp_moved.xml";
1132  //rewrite the Rollsessionfile (we start from the waleffe forcing)
1133  //string meshbndjumps = m_sessionName +"_bndjumps.xml";
1134  //if(cnt==0)
1135  //{
1136  //take the conditions tag from meshbndjumps and copy into
1137  // the rolls session file
1138  //}
1139 
1140 
1141  sprintf(c,"%d",cnt);
1142  //save old roll solution
1143  string oldroll = m_sessionName +"_roll_"+c +".fld";
1144  syscall = "cp -f " + m_sessionName+"-Base.fld" + " " + oldroll;
1145  cout<<syscall.c_str()<<endl;
1146  if(system(syscall.c_str()))
1147  {
1148  ASSERTL0(false,syscall.c_str());
1149  }
1150  //define file names
1151  string filePost = m_sessionName + "_advPost.xml";
1152  string filestreak = m_sessionName + "_streak.fld";
1153  string filewave = m_sessionName + "_wave.fld";
1154  string filewavepressure = m_sessionName + "_wave_p_split.fld";
1155  string fileinterp = m_sessionName + "_interp.xml";
1156  string interpstreak = m_sessionName +"_interpstreak_"+ c +".fld";
1157  string interwavepressure = m_sessionName +"_wave_p_split_interp_"+ c +".fld";
1158  char alpchar[16]="";
1159 cout<<"alpha = "<<m_alpha[0]<<endl;
1160  sprintf(alpchar, "%f", m_alpha[0]);
1161 
1162 
1163  if( m_sessionVWI->GetSolverInfo("INTERFACE")!="phase" )
1164  {
1165  cout<<"zerophase"<<endl;
1166 
1167  syscall = "../../utilities/PostProcessing/Extras/MoveMesh "
1168  + filePost +" "+ filestreak +" "+ fileinterp + " "+ alpchar;
1169 
1170  cout<<syscall.c_str()<<endl;
1171  if(system(syscall.c_str()))
1172  {
1173  ASSERTL0(false,syscall.c_str());
1174  }
1175 
1176  //move the advPost mesh (remark update alpha!!!)
1177  syscall = "../../utilities/PostProcessing/Extras/MoveMesh "
1178  + filePost + " " + filestreak + " " + filePost + " "+ alpchar;
1179  cout<<syscall.c_str()<<endl;
1180  if(system(syscall.c_str()))
1181  {
1182  ASSERTL0(false,syscall.c_str());
1183  }
1184 
1185 
1186 
1187  //save oldstreak
1188  string oldstreak = m_sessionName +"_streak_"+ c +".fld";
1189  syscall = "cp -f " + filestreak + " " + oldstreak;
1190  cout<<syscall.c_str()<<endl;
1191  if(system(syscall.c_str()))
1192  {
1193  ASSERTL0(false,syscall.c_str());
1194  }
1195 
1196  //interpolate the streak field into the new mesh
1197  string movedmesh = m_sessionName + "_advPost_moved.xml";
1198  string movedinterpmesh = m_sessionName + "_interp_moved.xml";
1199 
1200  //create the interp streak
1201 
1202  syscall = "../../utilities/PostProcessing/Extras/FieldToField "
1203  + fileinterp + " " + filestreak + " " + movedinterpmesh + " "
1204  + interpstreak;
1205 
1206  cout<<syscall.c_str()<<endl;
1207  if(system(syscall.c_str()))
1208  {
1209  ASSERTL0(false,syscall.c_str());
1210  }
1211 
1212 
1213  //save the old mesh
1214  string meshfile = m_sessionName + ".xml";
1215  string meshold = m_sessionName +"_"+ c +".xml";
1216  syscall = "cp -f " + meshfile + " " + meshold;
1217  cout<<syscall.c_str()<<endl;
1218  if(system(syscall.c_str()))
1219  {
1220  ASSERTL0(false,syscall.c_str());
1221  }
1222 
1223  //overwriting the meshfile with the new mesh
1224  syscall = "cp -f " + movedmesh + " " + meshfile;
1225  cout<<syscall.c_str()<<endl;
1226  if(system(syscall.c_str()))
1227  {
1228  ASSERTL0(false,syscall.c_str());
1229  }
1230 
1231  //overwriting the streak file!!
1232  syscall = "cp -f " + interpstreak + " " + filestreak;
1233  cout<<syscall.c_str()<<endl;
1234  if(system(syscall.c_str()))
1235  {
1236  ASSERTL0(false,syscall.c_str());
1237  }
1238 
1239  //calculate the wave
1240  ExecuteWave();
1241 
1242  //save the wave field:
1243  string oldwave = m_sessionName +"_wave_"+c +".fld";
1244  syscall = "cp -f " + m_sessionName+".fld" + " " + oldwave;
1245  cout<<syscall.c_str()<<endl;
1246  if(system(syscall.c_str()))
1247  {
1248  ASSERTL0(false,syscall.c_str());
1249  }
1250 
1251  //save old jump conditions:
1252  string ujump = m_sessionName+"_u_5.bc";
1253  syscall = "cp -f " + ujump + " " + m_sessionName+"_u_5.bc_"+c;
1254  cout<<syscall.c_str()<<endl;
1255  if(system(syscall.c_str()))
1256  {
1257  ASSERTL0(false,syscall.c_str());
1258  }
1259 
1260  string vjump = m_sessionName+"_v_5.bc";
1261  syscall = "cp -f " + vjump + " " + m_sessionName+"_v_5.bc_"+c;
1262  cout<<syscall.c_str()<<endl;
1263  if(system(syscall.c_str()))
1264  {
1265  ASSERTL0(false,syscall.c_str());
1266  }
1267  cnt++;
1268 
1269 
1270 
1271 
1272  //use relaxation
1274  )
1275  {
1276  // the critical layer should be the bnd region 3
1277  //int reg =3;
1278  //FileRelaxation(reg);
1279  }
1280  char c1[16]="";
1281  sprintf(c1,"%d",cnt);
1282  //calculate the jump conditions
1283  string wavefile = m_sessionName +".fld";
1284  syscall = "../../utilities/PostProcessing/Extras/FldCalcBCs "
1285  + movedmesh + " " + wavefile + " " + interpstreak + "> data"+c1;
1286  cout<<syscall.c_str()<<endl;
1287  if(system(syscall.c_str()))
1288  {
1289  ASSERTL0(false,syscall.c_str());
1290  }
1291 
1292  //move the new name_interp_moved.xml into name_interp.xml
1293  syscall = "cp -f " + movedinterpmesh + " " + fileinterp;
1294  cout<<syscall.c_str()<<endl;
1295  if(system(syscall.c_str()))
1296  {
1297  ASSERTL0(false,syscall.c_str());
1298  }
1299  //move the new name_advPost_moved.xml into name_advPost.xml
1300  syscall = "cp -f " + movedmesh + " " + filePost;
1301  cout<<syscall.c_str()<<endl;
1302  if(system(syscall.c_str()))
1303  {
1304  ASSERTL0(false,syscall.c_str());
1305  }
1306 
1307 
1308  }
1309  else if( m_sessionVWI->GetSolverInfo("INTERFACE")=="phase" )
1310  {
1311 cout<<"phase"<<endl;
1312  //determine cr:
1313  NekDouble cr;
1314  string cr_str;
1315  stringstream st;
1316 
1317  //calculate the wave
1318  ExecuteWave();
1319 
1320  //save oldstreak
1321  string oldstreak = m_sessionName +"_streak_"+ c +".fld";
1322  syscall = "cp -f " + filestreak + " " + oldstreak;
1323  cout<<syscall.c_str()<<endl;
1324  if(system(syscall.c_str()))
1325  {
1326  ASSERTL0(false,syscall.c_str());
1327  }
1328 
1329  //save wave
1330  syscall = "cp -f " + m_sessionName+".fld" + " " + filewave;
1331  cout<<syscall.c_str()<<endl;
1332  if(system(syscall.c_str()))
1333  {
1334  ASSERTL0(false,syscall.c_str());
1335  }
1336  //save the old mesh
1337  string meshfile = m_sessionName + ".xml";
1338  string meshold = m_sessionName +"_"+ c +".xml";
1339  syscall = "cp -f " + meshfile + " " + meshold;
1340  cout<<syscall.c_str()<<endl;
1341  if(system(syscall.c_str()))
1342  {
1343  ASSERTL0(false,syscall.c_str());
1344  }
1345 
1346  //save the oldwave field:
1347  string oldwave = m_sessionName +"_wave_"+c +".fld";
1348  syscall = "cp -f " + m_sessionName+".fld" + " " + oldwave;
1349  cout<<syscall.c_str()<<endl;
1350  if(system(syscall.c_str()))
1351  {
1352  ASSERTL0(false,syscall.c_str());
1353  }
1354 
1355  //save old jump conditions:
1356  string ujump = m_sessionName+"_u_5.bc";
1357  syscall = "cp -f " + ujump + " " + m_sessionName+"_u_5.bc_"+c;
1358  cout<<syscall.c_str()<<endl;
1359  if(system(syscall.c_str()))
1360  {
1361  ASSERTL0(false,syscall.c_str());
1362  }
1363 
1364  string vjump = m_sessionName+"_v_5.bc";
1365  syscall = "cp -f " + vjump + " " + m_sessionName+"_v_5.bc_"+c;
1366  cout<<syscall.c_str()<<endl;
1367  if(system(syscall.c_str()))
1368  {
1369  ASSERTL0(false,syscall.c_str());
1370  }
1371  cnt++;
1372 
1373 
1374  cr = m_leading_imag_evl[0]/m_alpha[0];
1375  st << cr;
1376  cr_str = st.str();
1377 cout<<"cr="<<cr_str<<endl;
1378  //NB -g or NOT!!!
1379  //move the mesh around the critical layer
1380  syscall = "../../utilities/PostProcessing/Extras/MoveMesh "
1381  + filePost +" "+ filestreak +" "+ fileinterp + " "+ alpchar
1382  +" "+cr_str;
1383 
1384  cout<<syscall.c_str()<<endl;
1385  if(system(syscall.c_str()))
1386  {
1387  ASSERTL0(false,syscall.c_str());
1388  }
1389  //NB -g or NOT!!!
1390  //move the advPost mesh (remark update alpha!!!)
1391  syscall = "../../utilities/PostProcessing/Extras/MoveMesh "
1392  + filePost + " " + filestreak + " " + filePost + " "+ alpchar
1393  +" "+cr_str;
1394  cout<<syscall.c_str()<<endl;
1395  if(system(syscall.c_str()))
1396  {
1397  ASSERTL0(false,syscall.c_str());
1398  }
1399 
1400  //interp streak into the new mesh
1401  syscall = "../../utilities/PostProcessing/Extras/FieldToField "
1402  + fileinterp + " " + filestreak + " " + movedinterpmesh + " "
1403  + interpstreak;
1404 
1405  cout<<syscall.c_str()<<endl;
1406  if(system(syscall.c_str()))
1407  {
1408  ASSERTL0(false,syscall.c_str());
1409  }
1410 
1411  //split wave sol
1412  syscall = "../../utilities/PostProcessing/Extras/SplitFld "
1413  + filePost + " " + filewave;
1414 
1415  cout<<syscall.c_str()<<endl;
1416  if(system(syscall.c_str()))
1417  {
1418  ASSERTL0(false,syscall.c_str());
1419  }
1420  //interp wave
1421  syscall = "../../utilities/PostProcessing/Extras/FieldToField "
1422  + filePost + " " + filewavepressure + " " + movedmesh + " "
1423  + interwavepressure;
1424 
1425  cout<<syscall.c_str()<<endl;
1426  if(system(syscall.c_str()))
1427  {
1428  ASSERTL0(false,syscall.c_str());
1429  }
1430 
1431 
1432 
1433 
1434  //use relaxation
1436  )
1437  {
1438  // the critical layer should be the bnd region 3
1439  //int reg =3;
1440  //FileRelaxation(reg);
1441  }
1442  char c1[16]="";
1443  sprintf(c1,"%d",cnt);
1444 
1445  //cp wavepressure to m_sessionName.fld(to get
1446  // the right bcs names using FldCalcBCs)
1447  syscall = "cp -f "+ interwavepressure +" "+m_sessionName+".fld";
1448  cout<<syscall.c_str()<<endl;
1449  if(system(syscall.c_str()))
1450  {
1451  ASSERTL0(false,syscall.c_str());
1452  }
1453 
1454  //calculate the jump conditions
1455  //NB -g or NOT!!!
1456  syscall = "../../utilities/PostProcessing/Extras/FldCalcBCs "
1457  + movedmesh + " " +m_sessionName+".fld" + " "
1458  + interpstreak + "> data"+c1;
1459  cout<<syscall.c_str()<<endl;
1460  if(system(syscall.c_str()))
1461  {
1462  ASSERTL0(false,syscall.c_str());
1463  }
1464 
1465 
1466  //overwriting the meshfile with the new mesh
1467  syscall = "cp -f " + movedmesh + " " + meshfile;
1468  cout<<syscall.c_str()<<endl;
1469  if(system(syscall.c_str()))
1470  {
1471  ASSERTL0(false,syscall.c_str());
1472  }
1473 
1474  //overwriting the streak file!! (maybe is useless)
1475  syscall = "cp -f " + interpstreak + " " + filestreak;
1476  cout<<syscall.c_str()<<endl;
1477  if(system(syscall.c_str()))
1478  {
1479  ASSERTL0(false,syscall.c_str());
1480  }
1481  //move the new name_interp_moved.xml into name_interp.xml
1482  syscall = "cp -f " + movedinterpmesh + " " + fileinterp;
1483  cout<<syscall.c_str()<<endl;
1484  if(system(syscall.c_str()))
1485  {
1486  ASSERTL0(false,syscall.c_str());
1487  }
1488  //move the new name_advPost_moved.xml into name_advPost.xml
1489  syscall = "cp -f " + movedmesh + " " + filePost;
1490  cout<<syscall.c_str()<<endl;
1491  if(system(syscall.c_str()))
1492  {
1493  ASSERTL0(false,syscall.c_str());
1494  }
1495  }
1496  }
1497  else
1498  {
1500  ExecuteRoll();
1502 
1503 #ifndef _WIN32
1504  sleep(3);
1505 #endif
1506  ExecuteStreak();
1507 #ifndef _WIN32
1508  sleep(3);
1509 #endif
1510 
1512  {
1513  string syscall;
1514  char alpchar[16]="";
1515  sprintf(alpchar, "%f", m_alpha[0]);
1516 
1517  string filePost = m_sessionName + "_advPost.xml";
1518  string filestreak = m_sessionName + "_streak.fld";
1519  string filewave = m_sessionName + "_wave.fld";
1520  string filewavepressure = m_sessionName + "_wave_p_split.fld";
1521  string fileinterp = m_sessionName + "_interp.xml";
1522  string interpstreak = m_sessionName +"_interpstreak.fld";
1523  string interwavepressure = m_sessionName +"_wave_p_split_interp.fld";
1524  syscall = "../../utilities/PostProcessing/Extras/MoveMesh "
1525  + filePost +" "+ filestreak +" "+ fileinterp + " "+ alpchar;
1526 
1527  cout<<syscall.c_str()<<endl;
1528  if(system(syscall.c_str()))
1529  {
1530  ASSERTL0(false,syscall.c_str());
1531  }
1532 
1533  //move the advPost mesh (remark update alpha!!!)
1534  syscall = "../../utilities/PostProcessing/Extras/MoveMesh "
1535  + filePost + " " + filestreak + " " + filePost + " "+ alpchar;
1536  cout<<syscall.c_str()<<endl;
1537  if(system(syscall.c_str()))
1538  {
1539  ASSERTL0(false,syscall.c_str());
1540  }
1541 
1542  //save oldstreak
1543  string oldstreak = m_sessionName +"_streak.fld";
1544  syscall = "cp -f " + filestreak + " " + oldstreak;
1545  cout<<syscall.c_str()<<endl;
1546  if(system(syscall.c_str()))
1547  {
1548  ASSERTL0(false,syscall.c_str());
1549  }
1550 
1551  //interpolate the streak field into the new mesh
1552  string movedmesh = m_sessionName + "_advPost_moved.xml";
1553  string movedinterpmesh = m_sessionName + "_interp_moved.xml";
1554 
1555  //create the interp streak
1556 
1557  syscall = "../../utilities/PostProcessing/Extras/FieldToField "
1558  + fileinterp + " " + filestreak + " " + movedinterpmesh
1559  + " " + interpstreak;
1560 
1561  cout<<syscall.c_str()<<endl;
1562  if(system(syscall.c_str()))
1563  {
1564  ASSERTL0(false,syscall.c_str());
1565  }
1566 
1567  //save the old mesh
1568  string meshfile = m_sessionName + ".xml";
1569  string meshold = m_sessionName + ".xml";
1570  syscall = "cp -f " + meshfile + " " + meshold;
1571  cout<<syscall.c_str()<<endl;
1572  if(system(syscall.c_str()))
1573  {
1574  ASSERTL0(false,syscall.c_str());
1575  }
1576 
1577  //overwriting the meshfile with the new mesh
1578  syscall = "cp -f " + movedmesh + " " + meshfile;
1579  cout<<syscall.c_str()<<endl;
1580  if(system(syscall.c_str()))
1581  {
1582  ASSERTL0(false,syscall.c_str());
1583  }
1584 
1585  //overwriting the streak file!!
1586  syscall = "cp -f " + interpstreak + " " + filestreak;
1587  cout<<syscall.c_str()<<endl;
1588  if(system(syscall.c_str()))
1589  {
1590  ASSERTL0(false,syscall.c_str());
1591  }
1592  }
1593 
1594  ExecuteWave();
1595 
1596  if(CalcWaveForce)
1597  {
1599  }
1600  }
1601  }
static void DisableManagement(std::string whichPool="")
Definition: NekManager.hpp:300
static void EnableManagement(std::string whichPool="")
Definition: NekManager.hpp:280
VWIIterationType GetVWIIterationType(void)
@ eFixedWaveForcingWithSubIterationOnAlpha

References ASSERTL0, CalcNonLinearWaveForce(), Nektar::LibUtilities::NekManager< KeyType, ValueT, opLessCreator >::DisableManagement(), Nektar::eFixedWaveForcingWithSubIterationOnAlpha, Nektar::LibUtilities::NekManager< KeyType, ValueT, opLessCreator >::EnableManagement(), ExecuteRoll(), ExecuteStreak(), ExecuteWave(), GetVWIIterationType(), m_alpha, m_leading_imag_evl, m_moveMeshToCriticalLayer, m_sessionName, m_sessionRoll, and m_sessionVWI.

Referenced by DoFixedForcingIteration(), and main().

◆ ExecuteRoll()

void Nektar::VortexWaveInteraction::ExecuteRoll ( void  )

Definition at line 356 of file VortexWaveInteraction.cpp.

357  {
358  //set up the equation system to update the mesh
359  if(m_sessionRoll->DefinesSolverInfo("INTERFACE"))
360  {
361  string vEquation = m_sessionRoll->GetSolverInfo("solvertype");
363  // The forcing terms are inserted as N bcs
364  // Execute Roll
365  cout << "Executing Roll solver" << endl;
366  solverRoll->DoInitialise();
367  solverRoll->DoSolve();
368  solverRoll->Output();
369  m_rollField = solverRoll->UpdateFields();
370  for(int g=0; g< solverRoll->GetNvariables(); ++g)
371  {
372  NekDouble vL2Error = solverRoll->L2Error(g,false);
373  NekDouble vLinfError = solverRoll->LinfError(g);
374  cout << "L 2 error (variable " << solverRoll->GetVariable(g) << ") : " << vL2Error << endl;
375  cout << "L inf error (variable " << solverRoll->GetVariable(g) << ") : " << vLinfError << endl;
376  }
377  }
378  else
379  {
381  {
382  string vEquation = m_sessionRoll->GetSolverInfo("solvertype");
384  }
385  else
386  {
387  const int npoints = m_solverRoll->GetNpoints();
388  const int ncoeffs = m_solverRoll->GetNcoeffs();
389 
390  static int init = 1;
391  if(init)
392  {
393  m_solverRoll->DoInitialise();
394  m_vwiForcingObj = std::dynamic_pointer_cast<SolverUtils::ForcingProgrammatic>(GetForcingFactory().CreateInstance("Programmatic", m_sessionRoll, m_solverRoll, m_solverRoll->UpdateFields(), m_solverRoll->UpdateFields().size() - 1, 0));
395 
396  std::vector<std::string> vFieldNames = m_sessionRoll->GetVariables();
397  vFieldNames.erase(vFieldNames.end()-1);
398 
399  m_solverRoll->GetFunction("BodyForce")->Evaluate(vFieldNames, m_vwiForcingObj->UpdateForces());
400 
401  // Scale forcing
402  for(int i = 0; i < m_vwiForcingObj->UpdateForces().size(); ++i)
403  {
404  m_solverRoll->UpdateFields()[0]->FwdTrans(m_vwiForcingObj->UpdateForces()[i], m_vwiForcing[2+i]);
405  Vmath::Smul(npoints,m_rollForceScale,m_vwiForcingObj->UpdateForces()[i],1,m_vwiForcingObj->UpdateForces()[i],1);
406  }
407 
408  IncNavierStokesSharedPtr ins = m_solverRoll->as<IncNavierStokes>();
409  ins->AddForcing(m_vwiForcingObj);
410 
411  init = 0;
412  }
413  else // use internal definition of forcing in m_vwiForcing
414  {
415  // Scale forcing
416  for(int i = 0; i < m_vwiForcingObj->UpdateForces().size(); ++i)
417  {
418  m_solverRoll->UpdateFields()[i]->BwdTrans(m_vwiForcing[i],m_vwiForcingObj->UpdateForces()[i]);
419  Vmath::Smul(npoints,m_rollForceScale,m_vwiForcingObj->UpdateForces()[i],1,m_vwiForcingObj->UpdateForces()[i],1);
420  }
421 
422  // Shift m_vwiForcing for new restart in case of relaxation
423  Vmath::Vcopy(ncoeffs,m_vwiForcing[0],1,m_vwiForcing[2],1);
424  Vmath::Vcopy(ncoeffs,m_vwiForcing[1],1,m_vwiForcing[3],1);
425  }
426  }
427 
428  // Execute Roll
429  cout << "Executing Roll solver" << endl;
430  m_solverRoll->DoSolve();
431  m_solverRoll->Output();
432  m_rollField = m_solverRoll->UpdateFields();
433  for(int g=0; g< m_solverRoll->GetNvariables(); ++g)
434  {
435  NekDouble vL2Error = m_solverRoll->L2Error(g,false);
436  NekDouble vLinfError = m_solverRoll->LinfError(g);
437  cout << "L 2 error (variable " << m_solverRoll->GetVariable(g) << ") : " << vL2Error << endl;
438  cout << "L inf error (variable " << m_solverRoll->GetVariable(g) << ") : " << vLinfError << endl;
439  }
440 
441 
442  }
443 
444  // Copy .fld file to .rst and base.fld
445  cout << "Executing cp -f session.fld session.rst" << endl;
446  CopyFile(".fld",".rst");
447 
448  // Write out data into base flow with variable Vx,Vy
449  cout << "Writing data to session-Base.fld" << endl;
450 
451  std::vector<std::string> variables(2);
452  variables[0] = "Vx"; variables[1] = "Vy";
453  std::vector<Array<OneD, NekDouble> > outfield(2);
454  outfield[0] = m_solverRoll->UpdateFields()[0]->UpdateCoeffs();
455  outfield[1] = m_solverRoll->UpdateFields()[1]->UpdateCoeffs();
456  std::string outname = m_sessionName + "-Base.fld";
457  m_solverRoll->WriteFld(outname, m_solverRoll->UpdateFields()[0],
458  outfield, variables);
459  }
void CopyFile(std::string file1end, std::string file2end)
SolverUtils::ForcingProgrammaticSharedPtr m_vwiForcingObj
Array< OneD, MultiRegions::ExpListSharedPtr > m_rollField
std::shared_ptr< EquationSystem > EquationSystemSharedPtr
A shared pointer to an EquationSystem object.
ForcingFactory & GetForcingFactory()
Declaration of the forcing factory singleton.
Definition: Forcing.cpp:44
std::shared_ptr< IncNavierStokes > IncNavierStokesSharedPtr

References Nektar::IncNavierStokes::AddForcing(), CopyFile(), Nektar::LibUtilities::NekFactory< tKey, tBase, tParam >::CreateInstance(), Nektar::SolverUtils::GetEquationSystemFactory(), Nektar::SolverUtils::GetForcingFactory(), m_graphRoll, m_moveMeshToCriticalLayer, m_rollField, m_rollForceScale, m_sessionName, m_sessionRoll, m_solverRoll, m_vwiForcing, m_vwiForcingObj, Vmath::Smul(), and Vmath::Vcopy().

Referenced by ExecuteLoop().

◆ ExecuteStreak()

void Nektar::VortexWaveInteraction::ExecuteStreak ( void  )

Definition at line 462 of file VortexWaveInteraction.cpp.

463  {
464  // Create driver
465 #if 1
466  std::string vDriverModule;
467  m_sessionStreak->LoadSolverInfo("Driver", vDriverModule, "Standard");
468 
470  solverStreak->Execute();
471 
472  m_streakField = solverStreak->GetEqu()[0]->UpdateFields();
473 #else
474  // Setup and execute Advection Diffusion solver
475  string vEquation = m_sessionStreak->GetSolverInfo("EqType");
477 
478  cout << "Executing Streak Solver" << endl;
479  solverStreak->DoInitialise();
480  solverStreak->DoSolve();
481  solverStreak->Output();
482 
483  m_streakField = solverStreak->UpdateFields();
484 
485  if(m_sessionVWI->DefinesSolverInfo("INTERFACE"))
486  {
487  for(int g=0; g< solverStreak->GetNvariables(); ++g)
488  {
489  NekDouble vL2Error = solverStreak->L2Error(g,false);
490  NekDouble vLinfError = solverStreak->LinfError(g);
491  cout << "L 2 error (variable " << solverStreak->GetVariable(g) << ") : " << vL2Error << endl;
492  cout << "L inf error (variable " << solverStreak->GetVariable(g) << ") : " << vLinfError << endl;
493  }
494  }
495 #endif
496 
497  cout << "Executing cp -f session.fld session_streak.fld" << endl;
498  CopyFile(".fld","_streak.fld");
499  }
std::shared_ptr< Driver > DriverSharedPtr
A shared pointer to a Driver object.
Definition: Driver.h:51
DriverFactory & GetDriverFactory()
Definition: Driver.cpp:65

References CopyFile(), Nektar::LibUtilities::NekFactory< tKey, tBase, tParam >::CreateInstance(), Nektar::SolverUtils::GetDriverFactory(), Nektar::SolverUtils::GetEquationSystemFactory(), m_graphStreak, m_sessionStreak, m_sessionVWI, and m_streakField.

Referenced by ExecuteLoop().

◆ ExecuteWave()

void Nektar::VortexWaveInteraction::ExecuteWave ( void  )

Do linearised NavierStokes Session with Modified Arnoldi

Definition at line 501 of file VortexWaveInteraction.cpp.

502  {
503 
504  // Set the initial beta value in stability to be equal to VWI file
505  std::string LZstr("LZ");
506  NekDouble LZ = 2*M_PI/m_alpha[0];
507  cout << "Setting LZ in Linearised solver to " << LZ << endl;
508  m_sessionWave->SetParameter(LZstr,LZ);
509 
510  // Create driver
511  std::string vDriverModule;
512  m_sessionWave->LoadSolverInfo("Driver", vDriverModule, "ModifiedArnoldi");
513  cout << "Setting up linearised NS Solver" << endl;
515 
516  /// Do linearised NavierStokes Session with Modified Arnoldi
517  cout << "Executing wave solution " << endl;
518  solverWave->Execute();
519 
520  // Copy file to a rst location for next restart
521  cout << "Executing cp -f session_eig_0 session_eig_0.rst" << endl;
522  CopyFile("_eig_0","_eig_0.rst");
523 
524  // Store data relevant to other operations
525  m_leading_real_evl[0] = solverWave->GetRealEvl()[0];
526  m_leading_imag_evl[0] = solverWave->GetImagEvl()[0];
527 
528  // note this will only be true for modified Arnoldi
529  NekDouble realShift = 0.0;
530  if(m_sessionWave->DefinesParameter("RealShift"))
531  {
532  bool defineshift;
533  // only use shift in modifiedArnoldi solver since
534  // implicitly handled in Arpack.
535  m_sessionWave->MatchSolverInfo("Driver","ModifiedArnoldi",defineshift,true);
536  if(defineshift)
537  {
538  realShift = m_sessionWave->GetParameter("RealShift");
539  }
540 
541  }
542 
543  // Set up leading eigenvalue from inverse
544  NekDouble invmag = 1.0/(m_leading_real_evl[0]*m_leading_real_evl[0] +
546  m_leading_real_evl[0] *= -invmag;
547  m_leading_real_evl[0] += realShift;
548  m_leading_imag_evl[0] *= invmag;
549 
550 
551  m_waveVelocities = solverWave->GetEqu()[0]->UpdateFields();
552  m_wavePressure = solverWave->GetEqu()[0]->GetPressure();
553 
554 
555  if( m_sessionVWI->DefinesSolverInfo("INTERFACE") )
556  {
557  cout << "Growth =" <<m_leading_real_evl[0]<<endl;
558  cout << "Phase =" <<m_leading_imag_evl[0]<<endl;
559  }
560 
561  }

References CopyFile(), Nektar::LibUtilities::NekFactory< tKey, tBase, tParam >::CreateInstance(), Nektar::SolverUtils::GetDriverFactory(), m_alpha, m_graphWave, m_leading_imag_evl, m_leading_real_evl, m_sessionVWI, m_sessionWave, m_wavePressure, and m_waveVelocities.

Referenced by CalcL2ToLinfPressure(), DoFixedForcingIteration(), ExecuteLoop(), and main().

◆ FileRelaxation()

void Nektar::VortexWaveInteraction::FileRelaxation ( int  reg)

Definition at line 2003 of file VortexWaveInteraction.cpp.

2004  {
2005  cout<<"relaxation..."<<endl;
2006  static int cnt=0;
2007  Array<OneD, MultiRegions::ExpListSharedPtr> Iexp
2008  =m_rollField[0]->GetBndCondExpansions();
2009  //cast to 1D explist (otherwise appenddata doesn't work)
2013  *std::static_pointer_cast<MultiRegions::ExpList>(Iexp[reg]));
2014  int nq = Ilayer->GetTotPoints();
2015  if( cnt==0)
2016  {
2017  m_bcsForcing = Array<OneD, Array<OneD, NekDouble> > (4);
2018  m_bcsForcing[0] = Array<OneD, NekDouble> (4*nq);
2019  for(int i = 1; i < 4; ++i)
2020  {
2021  m_bcsForcing[i] = m_bcsForcing[i-1] + nq;
2022  }
2023  }
2024 
2025  // Read in mesh from input file
2026 
2027  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef_u;
2028  std::vector<std::vector<NekDouble> > FieldData_u;
2029  string file = m_sessionName;
2030 
2031 
2032  file += "_u_5.bc";
2034  fld->Import(file,FieldDef_u, FieldData_u);
2035  Ilayer->ExtractDataToCoeffs(FieldDef_u[0], FieldData_u[0], FieldDef_u[0]->m_fields[0],Ilayer->UpdateCoeffs());
2036  Ilayer->BwdTrans_IterPerExp(Ilayer->GetCoeffs(), Ilayer->UpdatePhys());
2037 
2038  if(cnt==0)
2039  {
2040  Vmath::Vcopy(nq,Ilayer->UpdatePhys(),1,m_bcsForcing[2],1);
2041  }
2042  Vmath::Vcopy(nq,Ilayer->UpdatePhys(),1,m_bcsForcing[0],1);
2043 
2044  //case relaxation==0 an additional array is needed
2045  Array<OneD, NekDouble> tmp_forcing(nq, 0.0);
2046 
2047  if(cnt!=0)
2048  {
2049  if(m_vwiRelaxation==1.0)
2050  {
2051  Vmath::Vcopy(nq, m_bcsForcing[0],1, tmp_forcing,1);
2052  }
2054  m_bcsForcing[0],1,m_bcsForcing[0],1);
2056  m_bcsForcing[0],1,Ilayer->UpdatePhys(),1);
2057  //generate again the bcs files:
2058 
2059  Array<OneD, Array<OneD, NekDouble> > fieldcoeffs(1);
2060  Ilayer->FwdTrans_IterPerExp(Ilayer->GetPhys(),Ilayer->UpdateCoeffs());
2061  fieldcoeffs[0] = Ilayer->UpdateCoeffs();
2062  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef1 = Ilayer->GetFieldDefinitions();
2063  std::vector<std::vector<NekDouble> > FieldData_1(FieldDef1.size());;
2064  FieldDef1[0]->m_fields.push_back("u");
2065  Ilayer->AppendFieldData(FieldDef1[0], FieldData_1[0]);
2066  fld->Write(file,FieldDef1,FieldData_1);
2067  //save the bcs for the next iteration
2068  if(m_vwiRelaxation!=1.0)
2069  {
2070  Vmath::Smul(nq,1./(1.0-m_vwiRelaxation),
2071  m_bcsForcing[0],1,m_bcsForcing[0],1);
2072  Vmath::Vcopy(nq,m_bcsForcing[0],1,m_bcsForcing[2],1);
2073  }
2074  else
2075  {
2076  Vmath::Vcopy(nq, tmp_forcing,1, m_bcsForcing[2],1);
2077  }
2078  }
2079 
2080 
2081 
2082  file = m_sessionName+ "_v_5.bc";
2083 
2084  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef_v;
2085  std::vector<std::vector<NekDouble> > FieldData_v;
2086  fld->Import(file,FieldDef_v, FieldData_v);
2087  Ilayer->ExtractDataToCoeffs(FieldDef_v[0], FieldData_v[0], FieldDef_v[0]->m_fields[0],Ilayer->UpdateCoeffs());
2088  Ilayer->BwdTrans_IterPerExp(Ilayer->GetCoeffs(), Ilayer->UpdatePhys());
2089  if(cnt==0)
2090  {
2091  Vmath::Vcopy(nq,Ilayer->UpdatePhys(),1,m_bcsForcing[3],1);
2092  }
2093  Vmath::Vcopy(nq,Ilayer->UpdatePhys(),1,m_bcsForcing[1],1);
2094  if(cnt!=0)
2095  {
2096  if(m_vwiRelaxation==1.0)
2097  {
2098  Vmath::Vcopy(nq, m_bcsForcing[1],1, tmp_forcing,1);
2099  }
2101  m_bcsForcing[1],1,m_bcsForcing[1],1);
2103  m_bcsForcing[1],1,Ilayer->UpdatePhys(),1);
2104  //generate again the bcs files:
2105  Array<OneD, Array<OneD, NekDouble> > fieldcoeffs(1);
2106  Ilayer->FwdTrans_IterPerExp(Ilayer->GetPhys(),Ilayer->UpdateCoeffs());
2107  fieldcoeffs[0] = Ilayer->UpdateCoeffs();
2108  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef2 = Ilayer->GetFieldDefinitions();
2109  std::vector<std::vector<NekDouble> > FieldData_2(FieldDef2.size());;
2110  FieldDef2[0]->m_fields.push_back("v");
2111  Ilayer->AppendFieldData(FieldDef2[0], FieldData_2[0]);
2112  fld->Write(file,FieldDef2,FieldData_2);
2113  //save the bcs for the next iteration
2114  if(m_vwiRelaxation!=1.0)
2115  {
2116  Vmath::Smul(nq,1./(1.0-m_vwiRelaxation),
2117  m_bcsForcing[1],1,m_bcsForcing[1],1);
2118  Vmath::Vcopy(nq,m_bcsForcing[1],1,m_bcsForcing[3],1);
2119  }
2120  else
2121  {
2122  Vmath::Vcopy(nq, tmp_forcing,1, m_bcsForcing[3],1);
2123  }
2124 
2125 
2126  }
2127 
2128 
2129  cnt++;
2130  }
static std::shared_ptr< FieldIO > CreateDefault(const LibUtilities::SessionReaderSharedPtr session)
Returns an object for the default FieldIO method.
Definition: FieldIO.cpp:195
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Array< OneD, Array< OneD, NekDouble > > m_bcsForcing
std::shared_ptr< FieldIO > FieldIOSharedPtr
Definition: FieldIO.h:306
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::LibUtilities::FieldIO::CreateDefault(), m_bcsForcing, m_rollField, m_sessionName, m_sessionVWI, m_vwiRelaxation, Vmath::Smul(), Vmath::Svtvp(), and Vmath::Vcopy().

◆ GetAlpha()

NekDouble Nektar::VortexWaveInteraction::GetAlpha ( void  )
inline

Definition at line 129 of file VortexWaveInteraction.h.

130  {
131  return m_alpha[0];
132  }

Referenced by DoFixedForcingIteration(), and main().

◆ GetAlphaStep()

NekDouble Nektar::VortexWaveInteraction::GetAlphaStep ( void  )
inline

Definition at line 135 of file VortexWaveInteraction.h.

136  {
137  return m_alphaStep;
138  }

Referenced by main().

◆ GetDAlphaDWaveForceMag()

NekDouble Nektar::VortexWaveInteraction::GetDAlphaDWaveForceMag ( void  )
inline

Definition at line 150 of file VortexWaveInteraction.h.

151  {
152  return m_dAlphaDWaveForceMag;
153  }

◆ GetEigRelTol()

NekDouble Nektar::VortexWaveInteraction::GetEigRelTol ( void  )
inline

Definition at line 161 of file VortexWaveInteraction.h.

163  {
164  return m_eigRelTol;
165  }

Referenced by DoFixedForcingIteration().

◆ GetIterEnd()

int Nektar::VortexWaveInteraction::GetIterEnd ( )
inline

Definition at line 109 of file VortexWaveInteraction.h.

110  {
111  return m_iterEnd;
112  }

Referenced by DoFixedForcingIteration().

◆ GetIterStart()

int Nektar::VortexWaveInteraction::GetIterStart ( )
inline

Definition at line 104 of file VortexWaveInteraction.h.

105  {
106  return m_iterStart;
107  }

Referenced by DoFixedForcingIteration().

◆ GetMaxOuterIterations()

int Nektar::VortexWaveInteraction::GetMaxOuterIterations ( void  )
inline

Definition at line 124 of file VortexWaveInteraction.h.

125  {
126  return m_maxOuterIterations;
127  }

Referenced by DoFixedForcingIteration().

◆ GetMaxWaveForceMagIter()

int Nektar::VortexWaveInteraction::GetMaxWaveForceMagIter ( void  )
inline

Definition at line 155 of file VortexWaveInteraction.h.

156  {
157  return m_maxWaveForceMagIter;
158  }

Referenced by main().

◆ GetMinInnerIterations()

int Nektar::VortexWaveInteraction::GetMinInnerIterations ( void  )
inline

Definition at line 167 of file VortexWaveInteraction.h.

168  {
169  return m_minInnerIterations;
170  }

Referenced by DoFixedForcingIteration().

◆ GetNOuterIterations()

int Nektar::VortexWaveInteraction::GetNOuterIterations ( void  )
inline

Definition at line 119 of file VortexWaveInteraction.h.

120  {
121  return m_nOuterIterations;
122  }

Referenced by DoFixedForcingIteration().

◆ GetPrevAlpha()

NekDouble Nektar::VortexWaveInteraction::GetPrevAlpha ( void  )
inline

Definition at line 172 of file VortexWaveInteraction.h.

173  {
174  return m_prevAlpha;
175  }

Referenced by DoFixedForcingIteration().

◆ GetReflectionIndex()

Array< OneD, int > Nektar::VortexWaveInteraction::GetReflectionIndex ( void  )

Definition at line 1923 of file VortexWaveInteraction.cpp.

1924  {
1925  int i,j;
1926  int npts = m_waveVelocities[0]->GetPlane(0)->GetNpoints();
1927  int nel = m_waveVelocities[0]->GetNumElmts();
1928  Array<OneD, int> index(npts);
1929 
1930  Array<OneD, NekDouble> coord(2);
1931  Array<OneD, NekDouble> coord_x(npts);
1932  Array<OneD, NekDouble> coord_y(npts);
1933 
1934  //-> Dermine the point which is on coordinate (x -> -x + Lx/2, y-> -y)
1935  m_waveVelocities[0]->GetPlane(0)->GetCoords(coord_x,coord_y);
1936  NekDouble xmax = Vmath::Vmax(npts,coord_x,1);
1937  //NekDouble tol = NekConstants::kGeomFactorsTol*NekConstants::kGeomFactorsTol;
1938  NekDouble tol = 1e-5;
1939  NekDouble xnew,ynew;
1940 
1941  int start = npts-1;
1942  int e_npts;
1943 
1944  bool useOnlyQuads = false;
1945  if(m_sessionVWI->DefinesSolverInfo("SymmetriseOnlyQuads"))
1946  {
1947  useOnlyQuads = true;
1948  }
1949 
1950  int cnt;
1951  for(int e = 0; e < nel; ++e)
1952  {
1953  e_npts = m_waveVelocities[0]->GetExp(e)->GetTotPoints();
1954  cnt = m_waveVelocities[0]->GetPhys_Offset(e);
1955 
1956  if(useOnlyQuads)
1957  {
1958  if(m_waveVelocities[0]->GetExp(e)->DetShapeType() == LibUtilities::eTriangle)
1959  {
1960  for(i = 0; i < e_npts; ++i)
1961  {
1962  index[cnt+i] = -1;
1963  }
1964  continue;
1965  }
1966  }
1967 
1968  for(i = cnt; i < cnt+e_npts; ++i)
1969  {
1970  xnew = - coord_x[i] + xmax;
1971  ynew = - coord_y[i];
1972 
1973  for(j = start; j >=0 ; --j)
1974  {
1975  if((coord_x[j]-xnew)*(coord_x[j]-xnew) + (coord_y[j]-ynew)*(coord_y[j]-ynew) < tol)
1976  {
1977  index[i] = j;
1978  start = j;
1979  break;
1980  }
1981  }
1982 
1983  if(j == -1)
1984  {
1985 
1986  for(j = npts-1; j > start; --j)
1987  {
1988 
1989  if((coord_x[j]-xnew)*(coord_x[j]-xnew) + (coord_y[j]-ynew)*(coord_y[j]-ynew) < tol)
1990  {
1991  index[i] = j;
1992  break;
1993  }
1994  }
1995  ASSERTL0(j != start,"Failed to find matching point");
1996  }
1997  }
1998  }
1999  return index;
2000  }

References ASSERTL0, Nektar::LibUtilities::eTriangle, m_sessionVWI, m_waveVelocities, and Vmath::Vmax().

Referenced by CalcNonLinearWaveForce().

◆ GetVWIIterationType()

VWIIterationType Nektar::VortexWaveInteraction::GetVWIIterationType ( void  )
inline

Definition at line 114 of file VortexWaveInteraction.h.

115  {
116  return m_VWIIterationType;
117  }

Referenced by DoFixedForcingIteration(), ExecuteLoop(), and main().

◆ GetWaveForceMag()

NekDouble Nektar::VortexWaveInteraction::GetWaveForceMag ( void  )
inline

Definition at line 140 of file VortexWaveInteraction.h.

141  {
142  return m_waveForceMag[0];
143  }

Referenced by main().

◆ GetWaveForceMagStep()

NekDouble Nektar::VortexWaveInteraction::GetWaveForceMagStep ( void  )
inline

Definition at line 145 of file VortexWaveInteraction.h.

146  {
147  return m_waveForceMagStep;
148  }

Referenced by main().

◆ IfIterInterface()

bool Nektar::VortexWaveInteraction::IfIterInterface ( void  )
inline

Definition at line 208 of file VortexWaveInteraction.h.

209  {
210  return m_iterinterface;
211  }

Referenced by DoFixedForcingIteration().

◆ MoveFile()

void Nektar::VortexWaveInteraction::MoveFile ( std::string  fileend,
std::string  dir,
int  n 
)

Definition at line 1029 of file VortexWaveInteraction.cpp.

1030  {
1031  static map<string,int> opendir;
1032 
1033  if(opendir.find(dir) == opendir.end())
1034  {
1035  // make directory and presume will fail if it already exists
1036  string mkdir = "mkdir " + dir;
1037  ASSERTL0(system(mkdir.c_str()) == 0,
1038  "Failed to make directory '" + dir + "'");
1039  opendir[dir] = 1;
1040  }
1041 
1042  string savefile = dir + "/" + file + "." + boost::lexical_cast<std::string>(n);
1043  string syscall = "mv -f " + file + " " + savefile;
1044 
1045  ASSERTL0(system(syscall.c_str()) == 0, syscall.c_str());
1046  }

References ASSERTL0, and CG_Iterations::savefile.

Referenced by DoFixedForcingIteration().

◆ SaveFile()

void Nektar::VortexWaveInteraction::SaveFile ( std::string  fileend,
std::string  dir,
int  n 
)

Definition at line 1008 of file VortexWaveInteraction.cpp.

1009  {
1010  static map<string,int> opendir;
1011 
1012  if(opendir.find(dir) == opendir.end())
1013  {
1014  // make directory and presume will fail if it already exists
1015  string mkdir = "mkdir " + dir;
1016  ASSERTL0(system(mkdir.c_str()) == 0,
1017  "Failed to make directory '" + dir + "'");
1018 
1019  opendir[dir] = 1;
1020  }
1021 
1022  string savefile = dir + "/" + file + "." + boost::lexical_cast<std::string>(n);
1023  string syscall = "cp -f " + file + " " + savefile;
1024 
1025  ASSERTL0(system(syscall.c_str()) == 0, syscall.c_str());
1026  }

References ASSERTL0, and CG_Iterations::savefile.

Referenced by SaveLoopDetails().

◆ SaveLoopDetails()

void Nektar::VortexWaveInteraction::SaveLoopDetails ( std::string  dir,
int  i 
)

Definition at line 1076 of file VortexWaveInteraction.cpp.

1078  {
1079  // Save NS restart file
1080  SaveFile(m_sessionName + ".rst",SaveDir,i+1);
1081  // Save Streak Solution
1082  SaveFile(m_sessionName + "_streak.fld",SaveDir,i);
1083  // Save Wave solution output
1084  SaveFile(m_sessionName + ".evl",SaveDir,i);
1085  SaveFile(m_sessionName + "_eig_0",SaveDir,i);
1086  // Save new field file of eigenvalue
1087  SaveFile(m_sessionName + ".fld",SaveDir,i);
1088  if(!(m_sessionVWI->DefinesSolverInfo("INTERFACE")))
1089  {
1090  // Save new forcing file
1091  SaveFile(m_sessionName + ".vwi",SaveDir,i+1);
1092  }
1093 
1094  }
void SaveFile(std::string fileend, std::string dir, int n)

References m_sessionName, m_sessionVWI, and SaveFile().

Referenced by DoFixedForcingIteration().

◆ SetAlpha()

void Nektar::VortexWaveInteraction::SetAlpha ( NekDouble  alpha)
inline

Definition at line 177 of file VortexWaveInteraction.h.

178  {
179  m_alpha[0] = alpha;
180  }

Referenced by main().

◆ SetAlphaStep()

void Nektar::VortexWaveInteraction::SetAlphaStep ( NekDouble  step)
inline

Definition at line 193 of file VortexWaveInteraction.h.

194  {
195  m_alphaStep = step;
196  }

◆ SetEigRelTol()

void Nektar::VortexWaveInteraction::SetEigRelTol ( NekDouble  tol)
inline

Definition at line 188 of file VortexWaveInteraction.h.

189  {
190  m_eigRelTol = tol;
191  }

Referenced by DoFixedForcingIteration().

◆ SetMinInnerIterations()

void Nektar::VortexWaveInteraction::SetMinInnerIterations ( int  niter)
inline

Definition at line 198 of file VortexWaveInteraction.h.

199  {
200  m_minInnerIterations = niter;
201  }

Referenced by DoFixedForcingIteration().

◆ SetPrevAlpha()

void Nektar::VortexWaveInteraction::SetPrevAlpha ( NekDouble  alpha)
inline

Definition at line 203 of file VortexWaveInteraction.h.

204  {
205  m_prevAlpha = alpha;
206  }

Referenced by main().

◆ SetWaveForceMag()

void Nektar::VortexWaveInteraction::SetWaveForceMag ( NekDouble  mag)
inline

Definition at line 183 of file VortexWaveInteraction.h.

184  {
185  m_waveForceMag[0] = mag;
186  }

Referenced by main().

◆ UpdateAlpha()

void Nektar::VortexWaveInteraction::UpdateAlpha ( int  n)

Definition at line 1810 of file VortexWaveInteraction.cpp.

1811  {
1812  NekDouble alp_new = 0.0;
1813 
1814 
1815  if(outeriter == 1)
1816  {
1817  m_alpha[1] = m_alpha[0];
1818  if(m_leading_real_evl[0] > 0.0)
1819  {
1820  alp_new = m_alpha[0] + m_alphaStep;
1821  }
1822  else
1823  {
1824  alp_new = m_alpha[0] - m_alphaStep;
1825  }
1826  }
1827  else
1828  {
1829  int i;
1830  int nstore = (m_alpha.size() < outeriter)? m_alpha.size(): outeriter;
1831  Array<OneD, NekDouble> Alpha(nstore);
1832  Array<OneD, NekDouble> Growth(nstore);
1833 
1834  Vmath::Vcopy(nstore,m_alpha,1,Alpha,1);
1835  Vmath::Vcopy(nstore,m_leading_real_evl,1,Growth,1);
1836 
1837  // Sort Alpha Growth values;
1838  NekDouble store;
1839  int k;
1840  for(i = 0; i < nstore; ++i)
1841  {
1842  k = Vmath::Imin(nstore-i,&Alpha[i],1);
1843 
1844  store = Alpha[i];
1845  Alpha[i] = Alpha[i+k];
1846  Alpha[i+k] = store;
1847 
1848  store = Growth[i];
1849  Growth[i] = Growth[i+k];
1850  Growth[i+k] = store;
1851  }
1852 
1853  // See if we have any values that cross zero
1854  for(i = 0; i < nstore-1; ++i)
1855  {
1856  if(Growth[i]*Growth[i+1] < 0.0)
1857  {
1858  break;
1859  }
1860  }
1861 
1862  if(i != nstore-1)
1863  {
1864  if(nstore == 2)
1865  {
1866  alp_new = (Alpha[0]*Growth[1] - Alpha[1]*Growth[0])/(Growth[1]-Growth[0]);
1867  }
1868  else
1869  {
1870  // use a quadratic fit and step through 10000 points
1871  // to find zero.
1872  int j;
1873  int nsteps = 10000;
1874  int idx = (i == 0)?1:i;
1875  NekDouble da = Alpha[idx+1] - Alpha[idx-1];
1876  NekDouble gval_m1 = Growth[idx-1],a,gval;
1877  NekDouble c1 = Growth[idx-1]/(Alpha[idx-1]-Alpha[idx])/
1878  (Alpha[idx-1]-Alpha[idx+1]);
1879  NekDouble c2 = Growth[idx]/(Alpha[idx]-Alpha[idx-1])/
1880  (Alpha[idx]-Alpha[idx+1]);
1881  NekDouble c3 = Growth[idx+1]/(Alpha[idx+1]-Alpha[idx-1])/
1882  (Alpha[idx+1]-Alpha[idx]);
1883 
1884  for(j = 1; j < nsteps+1; ++j)
1885  {
1886  a = Alpha[i] + j*da/(NekDouble) nsteps;
1887  gval = c1*(a - Alpha[idx ])*(a - Alpha[idx+1])
1888  + c2*(a - Alpha[idx-1])*(a - Alpha[idx+1])
1889  + c3*(a - Alpha[idx-1])*(a - Alpha[idx]);
1890 
1891  if(gval*gval_m1 < 0.0)
1892  {
1893  alp_new = ((a+da/(NekDouble)nsteps)*gval - a*gval_m1)/
1894  (gval - gval_m1);
1895  break;
1896  }
1897  }
1898  }
1899  }
1900  else // step backward/forward
1901  {
1902  if(Growth[i] > 0.0)
1903  {
1904  alp_new = m_alpha[0] + m_alphaStep;
1905  }
1906  else
1907  {
1908  alp_new = m_alpha[0] - m_alphaStep;
1909  }
1910  }
1911  }
1912 
1913  for(int i = m_alpha.size()-1; i > 0; --i)
1914  {
1915  m_alpha[i] = m_alpha[i-1];
1918  }
1919 
1920  m_alpha[0] = alp_new;
1921  }
int Imin(int n, const T *x, const int incx)
Return the index of the minimum element in x.
Definition: Vmath.cpp:966

References Vmath::Imin(), m_alpha, m_alphaStep, m_leading_imag_evl, m_leading_real_evl, and Vmath::Vcopy().

Referenced by DoFixedForcingIteration(), and VortexWaveInteraction().

◆ UpdateDAlphaDWaveForceMag()

void Nektar::VortexWaveInteraction::UpdateDAlphaDWaveForceMag ( NekDouble  alphainit)

Definition at line 1805 of file VortexWaveInteraction.cpp.

1806  {
1808  }

References m_alpha, m_dAlphaDWaveForceMag, and m_waveForceMagStep.

Referenced by DoFixedForcingIteration().

◆ UpdateWaveForceMag()

void Nektar::VortexWaveInteraction::UpdateWaveForceMag ( int  n)

Definition at line 1691 of file VortexWaveInteraction.cpp.

1692  {
1693  NekDouble wavef_new = 0.0;
1694 
1695 
1696  if(outeriter == 1)
1697  {
1699  if(m_leading_real_evl[0] > 0.0)
1700  {
1701  wavef_new = m_waveForceMag[0] - m_waveForceMagStep;
1702  }
1703  else
1704  {
1705  wavef_new = m_waveForceMag[0] + m_waveForceMagStep;
1706  }
1707  }
1708  else
1709  {
1710  int i;
1711  int nstore = (m_waveForceMag.size() < outeriter)? m_waveForceMag.size(): outeriter;
1712  Array<OneD, NekDouble> WaveForce(nstore);
1713  Array<OneD, NekDouble> Growth(nstore);
1714 
1715  Vmath::Vcopy(nstore,m_waveForceMag,1,WaveForce,1);
1716  Vmath::Vcopy(nstore,m_leading_real_evl,1,Growth,1);
1717 
1718  // Sort WaveForce Growth values;
1719  NekDouble store;
1720  int k;
1721  for(i = 0; i < nstore; ++i)
1722  {
1723  k = Vmath::Imin(nstore-i,&WaveForce[i],1);
1724 
1725  store = WaveForce[i];
1726  WaveForce[i] = WaveForce[i+k];
1727  WaveForce[i+k] = store;
1728 
1729  store = Growth[i];
1730  Growth[i] = Growth[i+k];
1731  Growth[i+k] = store;
1732  }
1733 
1734  // See if we have any values that cross zero
1735  for(i = 0; i < nstore-1; ++i)
1736  {
1737  if(Growth[i]*Growth[i+1] < 0.0)
1738  {
1739  break;
1740  }
1741  }
1742 
1743  if(i != nstore-1)
1744  {
1745  if(nstore == 2)
1746  {
1747  wavef_new = (WaveForce[0]*Growth[1] - WaveForce[1]*Growth[0])/(Growth[1]-Growth[0]);
1748  }
1749  else
1750  {
1751  // use a quadratic fit and step through 10000 points
1752  // to find zero.
1753  int j;
1754  int nsteps = 10000;
1755  int idx = (i == 0)?1:i;
1756  NekDouble da = WaveForce[idx+1] - WaveForce[idx-1];
1757  NekDouble gval_m1 = Growth[idx-1],a,gval;
1758  NekDouble c1 = Growth[idx-1]/(WaveForce[idx-1]-WaveForce[idx])/
1759  (WaveForce[idx-1]-WaveForce[idx+1]);
1760  NekDouble c2 = Growth[idx]/(WaveForce[idx]-WaveForce[idx-1])/
1761  (WaveForce[idx]-WaveForce[idx+1]);
1762  NekDouble c3 = Growth[idx+1]/(WaveForce[idx+1]-WaveForce[idx-1])/
1763  (WaveForce[idx+1]-WaveForce[idx]);
1764 
1765  for(j = 1; j < nsteps+1; ++j)
1766  {
1767  a = WaveForce[i] + j*da/(NekDouble) nsteps;
1768  gval = c1*(a - WaveForce[idx ])*(a - WaveForce[idx+1])
1769  + c2*(a - WaveForce[idx-1])*(a - WaveForce[idx+1])
1770  + c3*(a - WaveForce[idx-1])*(a - WaveForce[idx]);
1771 
1772  if(gval*gval_m1 < 0.0)
1773  {
1774  wavef_new = ((a+da/(NekDouble)nsteps)*gval - a*gval_m1)/
1775  (gval - gval_m1);
1776  break;
1777  }
1778  }
1779  }
1780  }
1781  else // step backward/forward
1782  {
1783  if(Growth[i] > 0.0)
1784  {
1785  wavef_new = m_waveForceMag[0] - m_waveForceMagStep;
1786  }
1787  else
1788  {
1789  wavef_new = m_waveForceMag[0] + m_waveForceMagStep;
1790  }
1791  }
1792  }
1793 
1794  for(int i = m_waveForceMag.size()-1; i > 0; --i)
1795  {
1796  m_waveForceMag[i] = m_waveForceMag[i-1];
1799  }
1800 
1801  m_waveForceMag[0] = wavef_new;
1802 
1803  }

References Vmath::Imin(), m_leading_imag_evl, m_leading_real_evl, m_waveForceMag, m_waveForceMagStep, and Vmath::Vcopy().

Referenced by DoFixedForcingIteration().

Member Data Documentation

◆ m_alpha

Array<OneD, NekDouble> Nektar::VortexWaveInteraction::m_alpha
private

◆ m_alphaStep

NekDouble Nektar::VortexWaveInteraction::m_alphaStep
private

Definition at line 245 of file VortexWaveInteraction.h.

Referenced by UpdateAlpha(), and VortexWaveInteraction().

◆ m_bcsForcing

Array<OneD, Array<OneD, NekDouble > > Nektar::VortexWaveInteraction::m_bcsForcing
private

Definition at line 262 of file VortexWaveInteraction.h.

Referenced by FileRelaxation().

◆ m_dAlphaDWaveForceMag

NekDouble Nektar::VortexWaveInteraction::m_dAlphaDWaveForceMag
private

Definition at line 249 of file VortexWaveInteraction.h.

Referenced by UpdateDAlphaDWaveForceMag(), and VortexWaveInteraction().

◆ m_deltaFcnApprox

bool Nektar::VortexWaveInteraction::m_deltaFcnApprox
private

Definition at line 228 of file VortexWaveInteraction.h.

Referenced by CalcNonLinearWaveForce(), and VortexWaveInteraction().

◆ m_deltaFcnDecay

NekDouble Nektar::VortexWaveInteraction::m_deltaFcnDecay
private

Definition at line 233 of file VortexWaveInteraction.h.

Referenced by CalcNonLinearWaveForce(), and VortexWaveInteraction().

◆ m_eigRelTol

NekDouble Nektar::VortexWaveInteraction::m_eigRelTol
private

Definition at line 247 of file VortexWaveInteraction.h.

Referenced by CheckEigIsStationary(), and VortexWaveInteraction().

◆ m_graphRoll

SpatialDomains::MeshGraphSharedPtr Nektar::VortexWaveInteraction::m_graphRoll
private

Definition at line 272 of file VortexWaveInteraction.h.

Referenced by ExecuteRoll(), and VortexWaveInteraction().

◆ m_graphStreak

SpatialDomains::MeshGraphSharedPtr Nektar::VortexWaveInteraction::m_graphStreak
private

Definition at line 276 of file VortexWaveInteraction.h.

Referenced by ExecuteStreak(), and VortexWaveInteraction().

◆ m_graphWave

SpatialDomains::MeshGraphSharedPtr Nektar::VortexWaveInteraction::m_graphWave
private

Definition at line 278 of file VortexWaveInteraction.h.

Referenced by ExecuteWave(), and VortexWaveInteraction().

◆ m_iterEnd

int Nektar::VortexWaveInteraction::m_iterEnd
private

Definition at line 221 of file VortexWaveInteraction.h.

Referenced by VortexWaveInteraction().

◆ m_iterinterface

bool Nektar::VortexWaveInteraction::m_iterinterface
private

Definition at line 252 of file VortexWaveInteraction.h.

Referenced by VortexWaveInteraction().

◆ m_iterStart

int Nektar::VortexWaveInteraction::m_iterStart
private

Definition at line 220 of file VortexWaveInteraction.h.

Referenced by VortexWaveInteraction().

◆ m_leading_imag_evl

Array<OneD, NekDouble> Nektar::VortexWaveInteraction::m_leading_imag_evl
private

◆ m_leading_real_evl

Array<OneD, NekDouble> Nektar::VortexWaveInteraction::m_leading_real_evl
private

◆ m_maxOuterIterations

int Nektar::VortexWaveInteraction::m_maxOuterIterations
private

Definition at line 224 of file VortexWaveInteraction.h.

Referenced by VortexWaveInteraction().

◆ m_maxWaveForceMagIter

int Nektar::VortexWaveInteraction::m_maxWaveForceMagIter
private

Definition at line 226 of file VortexWaveInteraction.h.

Referenced by VortexWaveInteraction().

◆ m_minInnerIterations

int Nektar::VortexWaveInteraction::m_minInnerIterations
private

Definition at line 225 of file VortexWaveInteraction.h.

Referenced by CheckEigIsStationary(), and VortexWaveInteraction().

◆ m_moveMeshToCriticalLayer

bool Nektar::VortexWaveInteraction::m_moveMeshToCriticalLayer
private

Definition at line 231 of file VortexWaveInteraction.h.

Referenced by ExecuteLoop(), ExecuteRoll(), and VortexWaveInteraction().

◆ m_neutralPointTol

NekDouble Nektar::VortexWaveInteraction::m_neutralPointTol
private

Definition at line 246 of file VortexWaveInteraction.h.

Referenced by CheckIfAtNeutralPoint(), and VortexWaveInteraction().

◆ m_nOuterIterations

int Nektar::VortexWaveInteraction::m_nOuterIterations
private

Definition at line 223 of file VortexWaveInteraction.h.

Referenced by VortexWaveInteraction().

◆ m_prevAlpha

NekDouble Nektar::VortexWaveInteraction::m_prevAlpha
private

Definition at line 250 of file VortexWaveInteraction.h.

◆ m_rollField

Array<OneD, MultiRegions::ExpListSharedPtr> Nektar::VortexWaveInteraction::m_rollField
private

Definition at line 266 of file VortexWaveInteraction.h.

Referenced by ExecuteRoll(), and FileRelaxation().

◆ m_rollForceScale

NekDouble Nektar::VortexWaveInteraction::m_rollForceScale
private

Definition at line 238 of file VortexWaveInteraction.h.

Referenced by ExecuteRoll(), and VortexWaveInteraction().

◆ m_sessionName

std::string Nektar::VortexWaveInteraction::m_sessionName
private

◆ m_sessionRoll

LibUtilities::SessionReaderSharedPtr Nektar::VortexWaveInteraction::m_sessionRoll
private

◆ m_sessionStreak

LibUtilities::SessionReaderSharedPtr Nektar::VortexWaveInteraction::m_sessionStreak
private

Definition at line 275 of file VortexWaveInteraction.h.

Referenced by ExecuteStreak(), and VortexWaveInteraction().

◆ m_sessionVWI

LibUtilities::SessionReaderSharedPtr Nektar::VortexWaveInteraction::m_sessionVWI
private

◆ m_sessionWave

LibUtilities::SessionReaderSharedPtr Nektar::VortexWaveInteraction::m_sessionWave
private

Definition at line 277 of file VortexWaveInteraction.h.

Referenced by ExecuteWave(), and VortexWaveInteraction().

◆ m_solverRoll

EquationSystemSharedPtr Nektar::VortexWaveInteraction::m_solverRoll
private

◆ m_streakField

Array<OneD, MultiRegions::ExpListSharedPtr> Nektar::VortexWaveInteraction::m_streakField
private

Definition at line 264 of file VortexWaveInteraction.h.

Referenced by CalcNonLinearWaveForce(), and ExecuteStreak().

◆ m_useLinfPressureNorm

bool Nektar::VortexWaveInteraction::m_useLinfPressureNorm
private

Definition at line 229 of file VortexWaveInteraction.h.

Referenced by CalcNonLinearWaveForce(), and VortexWaveInteraction().

◆ m_vwiForcing

Array<OneD, Array<OneD, NekDouble > > Nektar::VortexWaveInteraction::m_vwiForcing
private

◆ m_vwiForcingObj

SolverUtils::ForcingProgrammaticSharedPtr Nektar::VortexWaveInteraction::m_vwiForcingObj
private

Definition at line 260 of file VortexWaveInteraction.h.

Referenced by ExecuteRoll().

◆ m_VWIIterationType

VWIIterationType Nektar::VortexWaveInteraction::m_VWIIterationType
private

Definition at line 254 of file VortexWaveInteraction.h.

Referenced by VortexWaveInteraction().

◆ m_vwiRelaxation

NekDouble Nektar::VortexWaveInteraction::m_vwiRelaxation
private

◆ m_waveForceMag

Array<OneD, NekDouble> Nektar::VortexWaveInteraction::m_waveForceMag
private

◆ m_waveForceMagStep

NekDouble Nektar::VortexWaveInteraction::m_waveForceMagStep
private

◆ m_wavePressure

MultiRegions::ExpListSharedPtr Nektar::VortexWaveInteraction::m_wavePressure
private

◆ m_waveVelocities

Array<OneD, MultiRegions::ExpListSharedPtr> Nektar::VortexWaveInteraction::m_waveVelocities
private