Nektar++
VortexWaveInteraction.cpp
Go to the documentation of this file.
1 //////////////////////////////////////////////////////////////////////////////
2 //
3 // File VortexWaveInteraction.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 //
20 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description: Vortex Wave Interaction class
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 #include <VortexWaveInteraction/VortexWaveInteraction.h>
36 #include <MultiRegions/GlobalLinSysKey.h>
37 #include <MultiRegions/ExpList1D.h>
38 #include <SolverUtils/Driver.h>
39 #include <IncNavierStokesSolver/EquationSystems/IncNavierStokes.h>
40 
41 using namespace std;
42 
43 namespace Nektar
44 {
45 
46  VortexWaveInteraction::VortexWaveInteraction(int argc, char * argv[]):
47  m_nOuterIterations(0)
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 
55  LibUtilities::SessionReaderSharedPtr session;
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.
64  m_sessionVWI = LibUtilities::SessionReader::CreateInstance(argc, argv,
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  {
117  m_moveMeshToCriticalLayer = true;
118  }
119  else
120  {
121  m_moveMeshToCriticalLayer = false;
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
138  m_minInnerIterations = (int) (log(0.1)/log(m_vwiRelaxation));
139  }
140  else
141  {
142  m_vwiRelaxation = 0.0;
143  m_minInnerIterations = 1;
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());
156  m_graphRoll = SpatialDomains::MeshGraph::Read(m_sessionRoll);
157  std::string vEquation = m_sessionRoll->GetSolverInfo("SolverType");
158  m_solverRoll = GetEquationSystemFactory().CreateInstance(vEquation, m_sessionRoll, m_graphRoll);
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());
190  m_graphStreak = SpatialDomains::MeshGraph::Read(m_sessionStreak);
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());
201  m_graphWave = SpatialDomains::MeshGraph::Read(m_sessionWave);
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  {
217  m_VWIIterationType = (VWIIterationType)i;
218  break;
219  }
220  }
221  }
222  else
223  {
224  m_VWIIterationType = eFixedAlphaWaveForcing;
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.num_elements());
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 
266  UpdateAlpha(m_nOuterIterations++);
267  }
268  else
269  {
270  cout << " No File OuterIter.his to restart from" << endl;
271  }
272  }
273  break;
274  case eFixedAlphaWaveForcing:
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  }
349 
350 
351  VortexWaveInteraction::~VortexWaveInteraction()
352  {
353  m_sessionVWI->Finalise();
354  }
355 
356  void VortexWaveInteraction::ExecuteRoll(void)
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");
362  EquationSystemSharedPtr solverRoll = GetEquationSystemFactory().CreateInstance(vEquation,m_sessionRoll,m_graphRoll);
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  {
380  if(m_moveMeshToCriticalLayer)
381  {
382  string vEquation = m_sessionRoll->GetSolverInfo("solvertype");
383  EquationSystemSharedPtr solverRoll = GetEquationSystemFactory().CreateInstance(vEquation,m_sessionRoll,m_graphRoll);
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().num_elements() - 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().num_elements(); ++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().num_elements(); ++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  }
460 
461 
462  void VortexWaveInteraction::ExecuteStreak(void)
463  {
464  // Create driver
465 #if 1
466  std::string vDriverModule;
467  m_sessionStreak->LoadSolverInfo("Driver", vDriverModule, "Standard");
468 
469  DriverSharedPtr solverStreak = GetDriverFactory().CreateInstance(vDriverModule, m_sessionStreak, m_graphStreak);
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");
476  EquationSystemSharedPtr solverStreak = GetEquationSystemFactory().CreateInstance(vEquation,m_sessionStreak, m_graphStreak);
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  }
500 
501  void VortexWaveInteraction::ExecuteWave(void)
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 sovler" << endl;
514  DriverSharedPtr solverWave = GetDriverFactory().CreateInstance(vDriverModule, m_sessionWave, m_graphWave);
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] +
545  m_leading_imag_evl[0]*m_leading_imag_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  }
562 
563  void VortexWaveInteraction::CalcNonLinearWaveForce(void)
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.num_elements(); ++i)
652  {
653  BndConds = m_waveVelocities[i]->GetBndConditions();
654  for(int j = 0; j < BndConds.num_elements(); ++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 
683  if(m_useLinfPressureNorm)
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)->PhysIntegral(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);
913  Vmath::Svtvp(ncoeffs,m_vwiRelaxation,m_vwiForcing[2],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);
918  Vmath::Svtvp(ncoeffs,m_vwiRelaxation,m_vwiForcing[3],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)->PhysIntegral(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  }
967 
968  void VortexWaveInteraction::CalcL2ToLinfPressure(void)
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)->PhysIntegral(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)->PhysIntegral(val);
1000 
1001  l2 = sqrt(norm/area);
1002 
1003  cout << "L2: " << l2 << endl;
1004 
1005  cout << "Ratio Linf/L2: "<< Linf/l2 << endl;
1006  }
1007 
1008  void VortexWaveInteraction::SaveFile(string file, string dir, int n)
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  }
1027 
1028 
1029  void VortexWaveInteraction::MoveFile(string file, string dir, int n)
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  }
1047 
1048  void VortexWaveInteraction::CopyFile(string file1end, string file2end)
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  }
1059 
1060  void VortexWaveInteraction::AppendEvlToFile(string file, int n)
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  }
1067 
1068  void VortexWaveInteraction::AppendEvlToFile(string file, NekDouble WaveForceMag)
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  }
1075 
1076  void VortexWaveInteraction::SaveLoopDetails(std::string SaveDir, int i)
1077 
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  }
1095 
1096  void VortexWaveInteraction::ExecuteLoop(bool CalcWaveForce)
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 
1116  LibUtilities::NekManager<MultiRegions::GlobalLinSysKey,MultiRegions::GlobalLinSys>::EnableManagement("GlobalLinSys");
1117  ExecuteRoll();
1118  LibUtilities::NekManager<MultiRegions::GlobalLinSysKey,MultiRegions::GlobalLinSys>::DisableManagement("GlobalLinSys");
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
1273  if(GetVWIIterationType()!=eFixedWaveForcingWithSubIterationOnAlpha
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
1435  if(GetVWIIterationType()!=eFixedWaveForcingWithSubIterationOnAlpha
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  {
1499  LibUtilities::NekManager<MultiRegions::GlobalLinSysKey,MultiRegions::GlobalLinSys>::EnableManagement("GlobalLinSys");
1500  ExecuteRoll();
1501  LibUtilities::NekManager<MultiRegions::GlobalLinSysKey,MultiRegions::GlobalLinSys>::DisableManagement("GlobalLinSys");
1502 
1503 #ifndef _WIN32
1504  sleep(3);
1505 #endif
1506  ExecuteStreak();
1507 #ifndef _WIN32
1508  sleep(3);
1509 #endif
1510 
1511  if(m_moveMeshToCriticalLayer)
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  {
1598  CalcNonLinearWaveForce();
1599  }
1600  }
1601  }
1602 
1603  bool VortexWaveInteraction::CheckEigIsStationary(bool reset)
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  }
1648 
1649  // Check to see if leading eigenvalue is within tolerance defined
1650  // in m_neutralPointTol
1651  bool VortexWaveInteraction::CheckIfAtNeutralPoint(void)
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  {
1674  if((m_leading_real_evl[0]*m_leading_real_evl[0] +
1675  m_leading_imag_evl[0]*m_leading_imag_evl[0]) <
1676  m_neutralPointTol*m_neutralPointTol)
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  }
1688 
1689  // Similar routine to UpdateAlpha
1690 
1691  void VortexWaveInteraction::UpdateWaveForceMag(int outeriter)
1692  {
1693  NekDouble wavef_new = 0.0;
1694 
1695 
1696  if(outeriter == 1)
1697  {
1698  m_waveForceMag[1] = m_waveForceMag[0];
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.num_elements() < outeriter)? m_waveForceMag.num_elements(): 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.num_elements()-1; i > 0; --i)
1795  {
1796  m_waveForceMag[i] = m_waveForceMag[i-1];
1797  m_leading_real_evl[i] = m_leading_real_evl[i-1];
1798  m_leading_imag_evl[i] = m_leading_imag_evl[i-1];
1799  }
1800 
1801  m_waveForceMag[0] = wavef_new;
1802 
1803  }
1804 
1805  void VortexWaveInteraction::UpdateDAlphaDWaveForceMag(NekDouble alpha_init)
1806  {
1807  m_dAlphaDWaveForceMag = (m_alpha[0]-alpha_init)/m_waveForceMagStep;
1808  }
1809 
1810  void VortexWaveInteraction::UpdateAlpha(int outeriter)
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.num_elements() < outeriter)? m_alpha.num_elements(): 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.num_elements()-1; i > 0; --i)
1914  {
1915  m_alpha[i] = m_alpha[i-1];
1916  m_leading_real_evl[i] = m_leading_real_evl[i-1];
1917  m_leading_imag_evl[i] = m_leading_imag_evl[i-1];
1918  }
1919 
1920  m_alpha[0] = alp_new;
1921  }
1922 
1923  Array<OneD, int> VortexWaveInteraction::GetReflectionIndex(void)
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  }
2001 
2002 
2003  void VortexWaveInteraction::FileRelaxation(int reg)
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)
2010  MultiRegions::ExpList1DSharedPtr Ilayer;
2011  Ilayer = MemoryManager<MultiRegions::ExpList1D>::
2012  AllocateSharedPtr(
2013  *std::static_pointer_cast<MultiRegions::ExpList1D>(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";
2033  LibUtilities::FieldIOSharedPtr fld = LibUtilities::FieldIO::CreateDefault(m_sessionVWI);
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  }
2053  Vmath::Smul(nq,1.0-m_vwiRelaxation,
2054  m_bcsForcing[0],1,m_bcsForcing[0],1);
2055  Vmath::Svtvp(nq,m_vwiRelaxation,m_bcsForcing[2],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  }
2100  Vmath::Smul(nq,1.0-m_vwiRelaxation,
2101  m_bcsForcing[1],1,m_bcsForcing[1],1);
2102  Vmath::Svtvp(nq,m_vwiRelaxation,m_bcsForcing[3],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  }
2131 }
Nektar::LibUtilities::NekManager::DisableManagement
static void DisableManagement(std::string whichPool="")
Definition: NekManager.hpp:300
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
Nektar::VortexWaveInteraction::m_sessionRoll
LibUtilities::SessionReaderSharedPtr m_sessionRoll
Definition: VortexWaveInteraction.h:271
Nektar::VortexWaveInteraction::m_wavePressure
MultiRegions::ExpListSharedPtr m_wavePressure
Definition: VortexWaveInteraction.h:257
Driver.h
Nektar::VortexWaveInteraction::m_graphWave
SpatialDomains::MeshGraphSharedPtr m_graphWave
Definition: VortexWaveInteraction.h:278
Vmath::Vsqrt
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:411
Nektar::IncNavierStokesSharedPtr
std::shared_ptr< IncNavierStokes > IncNavierStokesSharedPtr
Definition: IncNavierStokes.h:274
Vmath::Vmax
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:782
Nektar::VortexWaveInteraction::ExecuteLoop
void ExecuteLoop(bool CalcWaveForce=true)
Definition: VortexWaveInteraction.cpp:1096
Vmath::Fill
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:45
Nektar::VortexWaveInteraction::m_leading_imag_evl
Array< OneD, NekDouble > m_leading_imag_evl
< Leading real eigenvalue
Definition: VortexWaveInteraction.h:241
Vmath::Imin
int Imin(int n, const T *x, const int incx)
Return the index of the minimum element in x.
Definition: Vmath.cpp:850
Nektar::VortexWaveInteraction::m_rollField
Array< OneD, MultiRegions::ExpListSharedPtr > m_rollField
Definition: VortexWaveInteraction.h:266
Vmath::Svtvp
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
Vmath::Vvtvp
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:445
Nektar::SolverUtils::DriverSharedPtr
std::shared_ptr< Driver > DriverSharedPtr
A shared pointer to a Driver object.
Definition: Driver.h:51
Nektar::SolverUtils::GetForcingFactory
ForcingFactory & GetForcingFactory()
Declaration of the forcing factory singleton.
Definition: Forcing.cpp:44
std
STL namespace.
VortexWaveInteraction.h
Nektar::VortexWaveInteraction::CheckEigIsStationary
bool CheckEigIsStationary(bool reset=false)
Definition: VortexWaveInteraction.cpp:1603
Nektar::VortexWaveInteraction::GetVWIIterationType
VWIIterationType GetVWIIterationType(void)
Definition: VortexWaveInteraction.h:114
CG_Iterations.savefile
string savefile
Definition: CG_Iterations.py:312
Nektar::VortexWaveInteraction::SaveFile
void SaveFile(std::string fileend, std::string dir, int n)
Definition: VortexWaveInteraction.cpp:1008
Nektar::VortexWaveInteraction::m_waveVelocities
Array< OneD, MultiRegions::ExpListSharedPtr > m_waveVelocities
Definition: VortexWaveInteraction.h:256
Nektar::VortexWaveInteraction::UpdateDAlphaDWaveForceMag
void UpdateDAlphaDWaveForceMag(NekDouble alphainit)
Definition: VortexWaveInteraction.cpp:1805
Nektar::VortexWaveInteraction::MoveFile
void MoveFile(std::string fileend, std::string dir, int n)
Definition: VortexWaveInteraction.cpp:1029
Nektar::VortexWaveInteraction::m_solverRoll
EquationSystemSharedPtr m_solverRoll
Definition: VortexWaveInteraction.h:273
Nektar::VortexWaveInteraction::m_vwiForcingObj
SolverUtils::ForcingProgrammaticSharedPtr m_vwiForcingObj
Definition: VortexWaveInteraction.h:260
Nektar::LibUtilities::NekManager::EnableManagement
static void EnableManagement(std::string whichPool="")
Definition: NekManager.hpp:280
ExpList1D.h
Nektar::VortexWaveInteraction::m_leading_real_evl
Array< OneD, NekDouble > m_leading_real_evl
Definition: VortexWaveInteraction.h:240
Nektar::LibUtilities::SessionReader::CreateInstance
static SessionReaderSharedPtr CreateInstance(int argc, char *argv[])
Creates an instance of the SessionReader class.
Definition: SessionReader.h:132
Nektar::LibUtilities::NekFactory::CreateInstance
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:144
IncNavierStokes.h
Nektar::VortexWaveInteraction::m_bcsForcing
Array< OneD, Array< OneD, NekDouble > > m_bcsForcing
Definition: VortexWaveInteraction.h:262
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:216
Nektar::VortexWaveInteraction::m_sessionWave
LibUtilities::SessionReaderSharedPtr m_sessionWave
Definition: VortexWaveInteraction.h:277
Nektar::VWIIterationTypeMap
const std::string VWIIterationTypeMap[]
Definition: VortexWaveInteraction.h:63
Nektar::VortexWaveInteraction::m_waveForceMag
Array< OneD, NekDouble > m_waveForceMag
Definition: VortexWaveInteraction.h:235
GlobalLinSysKey.h
Nektar::MemoryManager::AllocateSharedPtr
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:161
Nektar::VortexWaveInteraction::SaveLoopDetails
void SaveLoopDetails(std::string dir, int i)
Definition: VortexWaveInteraction.cpp:1076
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:399
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Vmath::Sadd
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:318
Nektar::SolverUtils::GetEquationSystemFactory
EquationSystemFactory & GetEquationSystemFactory()
Definition: EquationSystem.cpp:90
Nektar::VortexWaveInteraction::m_vwiForcing
Array< OneD, Array< OneD, NekDouble > > m_vwiForcing
Definition: VortexWaveInteraction.h:259
Nektar::VortexWaveInteraction::m_sessionVWI
LibUtilities::SessionReaderSharedPtr m_sessionVWI
Definition: VortexWaveInteraction.h:269
Vmath::Vsub
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:346
Nektar::LibUtilities::FieldIO::CreateDefault
static std::shared_ptr< FieldIO > CreateDefault(const LibUtilities::SessionReaderSharedPtr session)
Returns an object for the default FieldIO method.
Definition: FieldIO.cpp:195
Nektar::IncNavierStokes::AddForcing
void AddForcing(const SolverUtils::ForcingSharedPtr &pForce)
Definition: IncNavierStokes.cpp:778
Nektar::VortexWaveInteraction::CopyFile
void CopyFile(std::string file1end, std::string file2end)
Definition: VortexWaveInteraction.cpp:1048
Nektar::IncNavierStokes
This class is the base class for Navier Stokes problems.
Definition: IncNavierStokes.h:136
Nektar::VortexWaveInteraction::m_graphRoll
SpatialDomains::MeshGraphSharedPtr m_graphRoll
Definition: VortexWaveInteraction.h:272
Nektar::SolverUtils::GetDriverFactory
DriverFactory & GetDriverFactory()
Definition: Driver.cpp:65
Nektar::VortexWaveInteraction::m_streakField
Array< OneD, MultiRegions::ExpListSharedPtr > m_streakField
Definition: VortexWaveInteraction.h:264
Nektar::VortexWaveInteraction::m_alpha
Array< OneD, NekDouble > m_alpha
< Leading imaginary eigenvalue
Definition: VortexWaveInteraction.h:243
Nektar::VortexWaveInteraction::GetReflectionIndex
Array< OneD, int > GetReflectionIndex(void)
Definition: VortexWaveInteraction.cpp:1923
Nektar::LibUtilities::FieldIOSharedPtr
std::shared_ptr< FieldIO > FieldIOSharedPtr
Definition: FieldIO.h:306
Nektar::VortexWaveInteraction::AppendEvlToFile
void AppendEvlToFile(std::string file, int n)
Definition: VortexWaveInteraction.cpp:1060
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1064
Nektar::MultiRegions::ExpList1DSharedPtr
std::shared_ptr< ExpList1D > ExpList1DSharedPtr
Shared pointer to an ExpList1D object.
Definition: ExpList1D.h:49
Nektar::LibUtilities::SessionReaderSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: SessionReader.h:112
Nektar::SpatialDomains::MeshGraph::Read
static MeshGraphSharedPtr Read(const LibUtilities::SessionReaderSharedPtr pSession, DomainRangeShPtr rng=NullDomainRangeShPtr, bool fillGraph=true)
Definition: MeshGraph.cpp:113
Vmath::Vadd
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:302
Vmath::Vmul
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:186
Nektar::VortexWaveInteraction::m_graphStreak
SpatialDomains::MeshGraphSharedPtr m_graphStreak
Definition: VortexWaveInteraction.h:276
Nektar::VortexWaveInteraction::m_sessionStreak
LibUtilities::SessionReaderSharedPtr m_sessionStreak
Definition: VortexWaveInteraction.h:275
Nektar::SolverUtils::EquationSystemSharedPtr
std::shared_ptr< EquationSystem > EquationSystemSharedPtr
A shared pointer to an EquationSystem object.
Definition: EquationSystem.h:59
Generated on Wed Jan 20 2021 03:33:43 for Nektar++ by   doxygen 1.8.13