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 <IncNavierStokesSolver/EquationSystems/IncNavierStokes.h>
36 #include <MultiRegions/ExpList.h>
37 #include <MultiRegions/GlobalLinSysKey.h>
38 #include <SolverUtils/Driver.h>
39 #include <VortexWaveInteraction/VortexWaveInteraction.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 =
65  LibUtilities::SessionReader::CreateInstance(argc, argv, 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  m_sessionVWI->LoadParameter("NeutralPointTolerance", m_neutralPointTol,
75  1e-4);
76  m_sessionVWI->LoadParameter("MaxOuterIterations", m_maxOuterIterations,
77  100);
78 
79  m_sessionVWI->LoadParameter("StartIteration", m_iterStart, 0);
80  m_sessionVWI->LoadParameter("EndIteration", m_iterEnd, 0);
81 
82  m_sessionVWI->LoadParameter("WaveForceMagStep", m_waveForceMagStep, 0.01);
83  m_sessionVWI->LoadParameter("DAlphaDWaveForceMag", m_dAlphaDWaveForceMag,
84  0.0);
85  m_sessionVWI->LoadParameter("MaxWaveForceMagIter", m_maxWaveForceMagIter,
86  1);
87  m_sessionVWI->LoadParameter("RollForceScale", m_rollForceScale, 1.0);
88 
89  if (m_sessionVWI->DefinesSolverInfo("DeltaFcnApprox"))
90  {
91  m_deltaFcnApprox = true;
92  m_sessionVWI->LoadParameter("DeltaFcnDecay", m_deltaFcnDecay,
93  1.0 / 500);
94  }
95  else
96  {
97  m_deltaFcnApprox = false;
98  m_deltaFcnDecay = 0.0;
99  }
100 
101  if (m_sessionVWI->DefinesSolverInfo("LinfPressureNorm"))
102  {
103  m_useLinfPressureNorm = true;
104  }
105  else
106  {
107  m_useLinfPressureNorm = false;
108  }
109 
110  if (m_sessionVWI->DefinesSolverInfo("INTERFACE"))
111  {
112  m_iterinterface = true;
113  }
114  else
115  {
116  m_iterinterface = false;
117  }
118 
119  if (m_sessionVWI->DefinesSolverInfo("MoveMeshToCriticalLayer"))
120  {
121  m_moveMeshToCriticalLayer = true;
122  }
123  else
124  {
125  m_moveMeshToCriticalLayer = false;
126  }
127 
128  m_alpha = Array<OneD, NekDouble>(storesize);
129  m_alpha[0] = m_sessionVWI->GetParameter("Alpha");
130  m_waveForceMag = Array<OneD, NekDouble>(storesize);
131  m_waveForceMag[0] = m_sessionVWI->GetParameter("WaveForceMag");
132 
133  m_leading_real_evl = Array<OneD, NekDouble>(storesize);
134  m_leading_imag_evl = Array<OneD, NekDouble>(storesize);
135 
136  if (m_sessionVWI->DefinesParameter("Relaxation"))
137  {
138  m_vwiRelaxation = m_sessionVWI->GetParameter("Relaxation");
139  // fix minimum number of iterations to be number of
140  // iterations required to make contribution of innitial
141  // forcing to 0.1
142  m_minInnerIterations = (int)(log(0.1) / log(m_vwiRelaxation));
143  }
144  else
145  {
146  m_vwiRelaxation = 0.0;
147  m_minInnerIterations = 1;
148  }
149 
150  // Initialise NS Roll solver
151  std::string IncCondFile(argv[argc - 1]);
152  IncCondFile += "_IncNSCond.xml";
153  std::vector<std::string> IncNSFilenames;
154  IncNSFilenames.push_back(meshfile);
155  IncNSFilenames.push_back(IncCondFile);
156 
157  // Create Incompressible NavierStokesSolver session reader.
158  m_sessionRoll = LibUtilities::SessionReader::CreateInstance(
159  argc, argv, IncNSFilenames, m_sessionVWI->GetComm());
160  m_graphRoll = SpatialDomains::MeshGraph::Read(m_sessionRoll);
161  std::string vEquation = m_sessionRoll->GetSolverInfo("SolverType");
162  m_solverRoll = GetEquationSystemFactory().CreateInstance(
163  vEquation, m_sessionRoll, m_graphRoll);
164  m_solverRoll->PrintSummary(cout);
165 
166  if (m_sessionRoll->DefinesSolverInfo("INTERFACE"))
167  {
168  m_sessionVWI->LoadParameter("EigenvalueRelativeTol", m_eigRelTol, 1e-2);
169  }
170 
171  int ncoeffs = m_solverRoll->UpdateFields()[0]->GetNcoeffs();
172  m_vwiForcing = Array<OneD, Array<OneD, NekDouble>>(4);
173  m_vwiForcing[0] = Array<OneD, NekDouble>(4 * ncoeffs);
174  for (int i = 1; i < 4; ++i)
175  {
176  m_vwiForcing[i] = m_vwiForcing[i - 1] + ncoeffs;
177  }
178 
179  // Fill forcing into m_vwiForcing
180  // Has forcing even been initialised yet?
181  // Vmath::Vcopy(ncoeffs,m_solverRoll->UpdateForces()[0]->GetCoeffs(),1,m_vwiForcing[0],1);
182  // Vmath::Vcopy(ncoeffs,m_solverRoll->UpdateForces()[1]->GetCoeffs(),1,m_vwiForcing[1],1);
183 
184  // Create AdvDiff Streak solver
185  std::string AdvDiffCondFile(argv[argc - 1]);
186  AdvDiffCondFile += "_AdvDiffCond.xml";
187  std::vector<std::string> AdvDiffFilenames;
188  AdvDiffFilenames.push_back(meshfile);
189  AdvDiffFilenames.push_back(AdvDiffCondFile);
190 
191  // Create AdvDiffusion session reader.
192  m_sessionStreak = LibUtilities::SessionReader::CreateInstance(
193  argc, argv, AdvDiffFilenames, m_sessionVWI->GetComm());
194  m_graphStreak = SpatialDomains::MeshGraph::Read(m_sessionStreak);
195 
196  // Initialise LinNS solver
197  std::string LinNSCondFile(argv[argc - 1]);
198  LinNSCondFile += "_LinNSCond.xml";
199  std::vector<std::string> LinNSFilenames;
200  LinNSFilenames.push_back(meshfile);
201  LinNSFilenames.push_back(LinNSCondFile);
202 
203  // Create Linearised NS stability session reader.
204  m_sessionWave = LibUtilities::SessionReader::CreateInstance(
205  argc, argv, LinNSFilenames, m_sessionVWI->GetComm());
206  m_graphWave = SpatialDomains::MeshGraph::Read(m_sessionWave);
207 
208  // Set the initial beta value in stability to be equal to VWI file
209  std::string LZstr("LZ");
210  NekDouble LZ = 2 * M_PI / m_alpha[0];
211  cout << "Setting LZ in Linearised solver to " << LZ << endl;
212  m_sessionWave->SetParameter(LZstr, LZ);
213 
214  // Check for iteration type
215  if (m_sessionVWI->DefinesSolverInfo("VWIIterationType"))
216  {
217  std::string IterationTypeStr =
218  m_sessionVWI->GetSolverInfo("VWIIterationType");
219  for (int i = 0; i < (int)eVWIIterationTypeSize; ++i)
220  {
221  if (boost::iequals(VWIIterationTypeMap[i], IterationTypeStr))
222  {
223  m_VWIIterationType = (VWIIterationType)i;
224  break;
225  }
226  }
227  }
228  else
229  {
230  m_VWIIterationType = eFixedAlphaWaveForcing;
231  }
232 
233  // Check for restart
234  bool restart;
235  m_sessionVWI->MatchSolverInfo("RestartIteration", "True", restart, false);
236  if (restart)
237  {
238  switch (m_VWIIterationType)
239  {
240  case eFixedAlpha:
241  break;
242  case eFixedWaveForcing:
243  {
244  FILE *fp;
245  // Check for OuterIter.his file to read
246  if ((fp = fopen("OuterIter.his", "r")))
247  {
248  char buf[BUFSIZ];
249  std::vector<NekDouble> Alpha, Growth, Phase;
250  NekDouble alpha, growth, phase;
251  while (fgets(buf, BUFSIZ, fp))
252  {
253  sscanf(buf, "%*d:%lf%lf%lf", &alpha, &growth, &phase);
254  Alpha.push_back(alpha);
255  Growth.push_back(growth);
256  Phase.push_back(phase);
257  }
258 
259  m_nOuterIterations = Alpha.size();
260 
261  int nvals =
262  std::min(m_nOuterIterations, (int)m_alpha.size());
263 
264  for (int i = 0; i < nvals; ++i)
265  {
266  m_alpha[nvals - 1 - i] =
267  Alpha[m_nOuterIterations - nvals + i];
268  m_leading_real_evl[nvals - 1 - i] =
269  Growth[m_nOuterIterations - nvals + i];
270  m_leading_imag_evl[nvals - 1 - i] =
271  Phase[m_nOuterIterations - nvals + i];
272  }
273 
274  UpdateAlpha(m_nOuterIterations++);
275  }
276  else
277  {
278  cout << " No File OuterIter.his to restart from" << endl;
279  }
280  }
281  break;
282  case eFixedAlphaWaveForcing:
283  {
284  string nstr = boost::lexical_cast<std::string>(m_iterStart);
285  cout << "Restarting from iteration " << m_iterStart << endl;
286  std::string rstfile = "cp -f Save/" + m_sessionName + ".rst." +
287  nstr + " " + m_sessionName + ".rst";
288  cout << " " << rstfile << endl;
289  if (system(rstfile.c_str()))
290  {
291  ASSERTL0(false, rstfile.c_str());
292  }
293  std::string vwifile = "cp -f Save/" + m_sessionName + ".vwi." +
294  nstr + " " + m_sessionName + ".vwi";
295  cout << " " << vwifile << endl;
296  if (system(vwifile.c_str()))
297  {
298  ASSERTL0(false, vwifile.c_str());
299  }
300  }
301  break;
302  default:
303  ASSERTL0(false, "Unknown VWIITerationType in restart");
304  }
305  }
306 
307  // Check for ConveredSoln to update DAlphaDWaveForce
308  {
309  FILE *fp;
310  if ((fp = fopen("ConvergedSolns", "r")))
311  {
312  char buf[BUFSIZ];
313  std::vector<NekDouble> WaveForce, Alpha;
314  NekDouble waveforce, alpha;
315  while (fgets(buf, BUFSIZ, fp))
316  {
317  sscanf(buf, "%*d:%lf%lf", &waveforce, &alpha);
318  WaveForce.push_back(waveforce);
319  Alpha.push_back(alpha);
320  }
321 
322  if (Alpha.size() > 1)
323  {
324  int min_i = 0;
325  NekDouble min_alph = fabs(m_alpha[0] - Alpha[min_i]);
326  // find nearest point
327  for (int i = 1; i < Alpha.size(); ++i)
328  {
329  if (fabs(m_alpha[0] - Alpha[min_i]) < min_alph)
330  {
331  min_i = i;
332  min_alph = fabs(m_alpha[0] - Alpha[min_i]);
333  }
334  }
335 
336  // find next nearest point
337  int min_j = (min_i == 0) ? 1 : 0;
338  min_alph = fabs(m_alpha[0] - Alpha[min_j]);
339  for (int i = 0; i < Alpha.size(); ++i)
340  {
341  if (i != min_i)
342  {
343  if (fabs(m_alpha[0] - Alpha[min_j]) < min_alph)
344  {
345  min_j = i;
346  min_alph = fabs(m_alpha[0] - Alpha[min_j]);
347  }
348  }
349  }
350 
351  if (fabs(Alpha[min_i] - Alpha[min_j]) > 1e-4)
352  {
353  m_dAlphaDWaveForceMag =
354  (Alpha[min_i] - Alpha[min_j]) /
355  (WaveForce[min_i] - WaveForce[min_j]);
356  }
357  }
358  }
359  }
360 }
361 
362 VortexWaveInteraction::~VortexWaveInteraction()
363 {
364  m_sessionVWI->Finalise();
365 }
366 
367 void VortexWaveInteraction::ExecuteRoll(void)
368 {
369  // set up the equation system to update the mesh
370  if (m_sessionRoll->DefinesSolverInfo("INTERFACE"))
371  {
372  string vEquation = m_sessionRoll->GetSolverInfo("solvertype");
373  EquationSystemSharedPtr solverRoll =
374  GetEquationSystemFactory().CreateInstance(vEquation, m_sessionRoll,
375  m_graphRoll);
376  // The forcing terms are inserted as N bcs
377  // Execute Roll
378  cout << "Executing Roll solver" << endl;
379  solverRoll->DoInitialise();
380  solverRoll->DoSolve();
381  solverRoll->Output();
382  m_rollField = solverRoll->UpdateFields();
383  for (int g = 0; g < solverRoll->GetNvariables(); ++g)
384  {
385  NekDouble vL2Error = solverRoll->L2Error(g, false);
386  NekDouble vLinfError = solverRoll->LinfError(g);
387  cout << "L 2 error (variable " << solverRoll->GetVariable(g)
388  << ") : " << vL2Error << endl;
389  cout << "L inf error (variable " << solverRoll->GetVariable(g)
390  << ") : " << vLinfError << endl;
391  }
392  }
393  else
394  {
395  if (m_moveMeshToCriticalLayer)
396  {
397  string vEquation = m_sessionRoll->GetSolverInfo("solvertype");
398  EquationSystemSharedPtr solverRoll =
399  GetEquationSystemFactory().CreateInstance(
400  vEquation, m_sessionRoll, m_graphRoll);
401  }
402  else
403  {
404  const int npoints = m_solverRoll->GetNpoints();
405  const int ncoeffs = m_solverRoll->GetNcoeffs();
406 
407  static int init = 1;
408  if (init)
409  {
410  m_solverRoll->DoInitialise();
411  m_vwiForcingObj =
412  std::dynamic_pointer_cast<SolverUtils::ForcingProgrammatic>(
413  GetForcingFactory().CreateInstance(
414  "Programmatic", m_sessionRoll, m_solverRoll,
415  m_solverRoll->UpdateFields(),
416  m_solverRoll->UpdateFields().size() - 1, 0));
417 
418  std::vector<std::string> vFieldNames =
419  m_sessionRoll->GetVariables();
420  vFieldNames.erase(vFieldNames.end() - 1);
421 
422  m_solverRoll->GetFunction("BodyForce")
423  ->Evaluate(vFieldNames, m_vwiForcingObj->UpdateForces());
424 
425  // Scale forcing
426  for (int i = 0; i < m_vwiForcingObj->UpdateForces().size(); ++i)
427  {
428  m_solverRoll->UpdateFields()[0]->FwdTrans(
429  m_vwiForcingObj->UpdateForces()[i],
430  m_vwiForcing[2 + i]);
431  Vmath::Smul(npoints, m_rollForceScale,
432  m_vwiForcingObj->UpdateForces()[i], 1,
433  m_vwiForcingObj->UpdateForces()[i], 1);
434  }
435 
436  IncNavierStokesSharedPtr ins =
437  m_solverRoll->as<IncNavierStokes>();
438  ins->AddForcing(m_vwiForcingObj);
439 
440  init = 0;
441  }
442  else // use internal definition of forcing in m_vwiForcing
443  {
444  // Scale forcing
445  for (int i = 0; i < m_vwiForcingObj->UpdateForces().size(); ++i)
446  {
447  m_solverRoll->UpdateFields()[i]->BwdTrans(
448  m_vwiForcing[i], m_vwiForcingObj->UpdateForces()[i]);
449  Vmath::Smul(npoints, m_rollForceScale,
450  m_vwiForcingObj->UpdateForces()[i], 1,
451  m_vwiForcingObj->UpdateForces()[i], 1);
452  }
453 
454  // Shift m_vwiForcing for new restart in case of relaxation
455  Vmath::Vcopy(ncoeffs, m_vwiForcing[0], 1, m_vwiForcing[2], 1);
456  Vmath::Vcopy(ncoeffs, m_vwiForcing[1], 1, m_vwiForcing[3], 1);
457  }
458  }
459 
460  // Execute Roll
461  cout << "Executing Roll solver" << endl;
462  m_solverRoll->DoSolve();
463  m_solverRoll->Output();
464  m_rollField = m_solverRoll->UpdateFields();
465  for (int g = 0; g < m_solverRoll->GetNvariables(); ++g)
466  {
467  NekDouble vL2Error = m_solverRoll->L2Error(g, false);
468  NekDouble vLinfError = m_solverRoll->LinfError(g);
469  cout << "L 2 error (variable " << m_solverRoll->GetVariable(g)
470  << ") : " << vL2Error << endl;
471  cout << "L inf error (variable " << m_solverRoll->GetVariable(g)
472  << ") : " << vLinfError << endl;
473  }
474  }
475 
476  // Copy .fld file to .rst and base.fld
477  cout << "Executing cp -f session.fld session.rst" << endl;
478  CopyFile(".fld", ".rst");
479 
480  // Write out data into base flow with variable Vx,Vy
481  cout << "Writing data to session-Base.fld" << endl;
482 
483  std::vector<std::string> variables(2);
484  variables[0] = "Vx";
485  variables[1] = "Vy";
486  std::vector<Array<OneD, NekDouble>> outfield(2);
487  outfield[0] = m_solverRoll->UpdateFields()[0]->UpdateCoeffs();
488  outfield[1] = m_solverRoll->UpdateFields()[1]->UpdateCoeffs();
489  std::string outname = m_sessionName + "-Base.fld";
490  m_solverRoll->WriteFld(outname, m_solverRoll->UpdateFields()[0], outfield,
491  variables);
492 }
493 
494 void VortexWaveInteraction::ExecuteStreak(void)
495 {
496  // Create driver
497 #if 1
498  std::string vDriverModule;
499  m_sessionStreak->LoadSolverInfo("Driver", vDriverModule, "Standard");
500 
501  DriverSharedPtr solverStreak = GetDriverFactory().CreateInstance(
502  vDriverModule, m_sessionStreak, m_graphStreak);
503  solverStreak->Execute();
504 
505  m_streakField = solverStreak->GetEqu()[0]->UpdateFields();
506 #else
507  // Setup and execute Advection Diffusion solver
508  string vEquation = m_sessionStreak->GetSolverInfo("EqType");
509  EquationSystemSharedPtr solverStreak =
510  GetEquationSystemFactory().CreateInstance(vEquation, m_sessionStreak,
511  m_graphStreak);
512 
513  cout << "Executing Streak Solver" << endl;
514  solverStreak->DoInitialise();
515  solverStreak->DoSolve();
516  solverStreak->Output();
517 
518  m_streakField = solverStreak->UpdateFields();
519 
520  if (m_sessionVWI->DefinesSolverInfo("INTERFACE"))
521  {
522  for (int g = 0; g < solverStreak->GetNvariables(); ++g)
523  {
524  NekDouble vL2Error = solverStreak->L2Error(g, false);
525  NekDouble vLinfError = solverStreak->LinfError(g);
526  cout << "L 2 error (variable " << solverStreak->GetVariable(g)
527  << ") : " << vL2Error << endl;
528  cout << "L inf error (variable " << solverStreak->GetVariable(g)
529  << ") : " << vLinfError << endl;
530  }
531  }
532 #endif
533 
534  cout << "Executing cp -f session.fld session_streak.fld" << endl;
535  CopyFile(".fld", "_streak.fld");
536 }
537 
538 void VortexWaveInteraction::ExecuteWave(void)
539 {
540 
541  // Set the initial beta value in stability to be equal to VWI file
542  std::string LZstr("LZ");
543  NekDouble LZ = 2 * M_PI / m_alpha[0];
544  cout << "Setting LZ in Linearised solver to " << LZ << endl;
545  m_sessionWave->SetParameter(LZstr, LZ);
546 
547  // Create driver
548  std::string vDriverModule;
549  m_sessionWave->LoadSolverInfo("Driver", vDriverModule, "ModifiedArnoldi");
550  cout << "Setting up linearised NS Solver" << endl;
551  DriverSharedPtr solverWave = GetDriverFactory().CreateInstance(
552  vDriverModule, m_sessionWave, m_graphWave);
553 
554  /// Do linearised NavierStokes Session with Modified Arnoldi
555  cout << "Executing wave solution " << endl;
556  solverWave->Execute();
557 
558  // Copy file to a rst location for next restart
559  cout << "Executing cp -f session_eig_0 session_eig_0.rst" << endl;
560  CopyFile("_eig_0", "_eig_0.rst");
561 
562  // Store data relevant to other operations
563  m_leading_real_evl[0] = solverWave->GetRealEvl()[0];
564  m_leading_imag_evl[0] = solverWave->GetImagEvl()[0];
565 
566  // note this will only be true for modified Arnoldi
567  NekDouble realShift = 0.0;
568  if (m_sessionWave->DefinesParameter("RealShift"))
569  {
570  bool defineshift;
571  // only use shift in modifiedArnoldi solver since
572  // implicitly handled in Arpack.
573  m_sessionWave->MatchSolverInfo("Driver", "ModifiedArnoldi", defineshift,
574  true);
575  if (defineshift)
576  {
577  realShift = m_sessionWave->GetParameter("RealShift");
578  }
579  }
580 
581  // Set up leading eigenvalue from inverse
582  NekDouble invmag = 1.0 / (m_leading_real_evl[0] * m_leading_real_evl[0] +
583  m_leading_imag_evl[0] * m_leading_imag_evl[0]);
584  m_leading_real_evl[0] *= -invmag;
585  m_leading_real_evl[0] += realShift;
586  m_leading_imag_evl[0] *= invmag;
587 
588  m_waveVelocities = solverWave->GetEqu()[0]->UpdateFields();
589  m_wavePressure = solverWave->GetEqu()[0]->GetPressure();
590 
591  if (m_sessionVWI->DefinesSolverInfo("INTERFACE"))
592  {
593  cout << "Growth =" << m_leading_real_evl[0] << endl;
594  cout << "Phase =" << m_leading_imag_evl[0] << endl;
595  }
596 }
597 
598 void VortexWaveInteraction::CalcNonLinearWaveForce(void)
599 {
600  if (m_sessionRoll->DefinesSolverInfo("INTERFACE"))
601  {
602  static int cnt = 0;
603  string wavefile = m_sessionName + ".fld";
604  string movedmesh = m_sessionName + "_advPost_moved.xml";
605  string filestreak = m_sessionName + "_streak.fld";
606  char c[16] = "";
607  sprintf(c, "%d", cnt);
608  char c_alpha[16] = "";
609  sprintf(c_alpha, "%f", m_alpha[0]);
610  string syscall;
611  if (m_sessionVWI->GetSolverInfo("INTERFACE") == "phase")
612  {
613  string filePost = m_sessionName + "_advPost.xml";
614  syscall = "../../utilities/PostProcessing/Extras/FldCalcBCs " +
615  filePost + " " +
616  "meshhalf_pos_Spen_stability_moved.fld "
617  "meshhalf_pos_Spen_advPost_moved.fld " +
618  c_alpha + " > data_alpha0";
619  cout << syscall.c_str() << endl;
620  if (system(syscall.c_str()))
621  {
622  ASSERTL0(false, syscall.c_str());
623  }
624 
625  syscall = "cp -f meshhalf_pos_Spen_stability_moved_u_5.bc " +
626  m_sessionName + "_u_5.bc";
627  cout << syscall.c_str() << endl;
628  if (system(syscall.c_str()))
629  {
630  ASSERTL0(false, syscall.c_str());
631  }
632  syscall = "cp -f meshhalf_pos_Spen_stability_moved_v_5.bc " +
633  m_sessionName + "_v_5.bc";
634  cout << syscall.c_str() << endl;
635  if (system(syscall.c_str()))
636  {
637  ASSERTL0(false, syscall.c_str());
638  }
639  }
640  else
641  {
642  syscall = "../../utilities/PostProcessing/Extras/FldCalcBCs " +
643  movedmesh + " " + wavefile + " " + filestreak + " " +
644  c_alpha + " > datasub_" + c;
645  cout << syscall.c_str() << endl;
646  if (system(syscall.c_str()))
647  {
648  ASSERTL0(false, syscall.c_str());
649  }
650  }
651 
652  // relaxation for different alpha values? does it make sense?
653 
654  // save the wave
655  string wave_subalp = m_sessionName + "_wave_subalp_" + c + ".fld";
656  syscall = "cp -f " + wavefile + " " + wave_subalp;
657  cout << syscall.c_str() << endl;
658  if (system(syscall.c_str()))
659  {
660  ASSERTL0(false, syscall.c_str());
661  }
662  // FileRelaxation(3);
663  cnt++;
664  }
665  else
666  {
667  int npts = m_waveVelocities[0]->GetPlane(0)->GetNpoints();
668  int ncoeffs = m_waveVelocities[0]->GetPlane(0)->GetNcoeffs();
669  Array<OneD, NekDouble> val(npts), der1(2 * npts);
670  Array<OneD, NekDouble> der2 = der1 + npts;
671  Array<OneD, NekDouble> streak;
672  static int projectfield = -1;
673 
674  if (m_deltaFcnApprox)
675  {
676  streak = Array<OneD, NekDouble>(npts);
677  m_streakField[0]->BwdTrans(m_streakField[0]->GetCoeffs(), streak);
678  }
679 
680  // Set project field to be first field that has a Neumann
681  // boundary since this not impose any condition on the vertical
682  // boundaries Othersise set to zero.
683  if (projectfield == -1)
684  {
685  Array<OneD, const SpatialDomains::BoundaryConditionShPtr> BndConds;
686 
687  for (int i = 0; i < m_waveVelocities.size(); ++i)
688  {
689  BndConds = m_waveVelocities[i]->GetBndConditions();
690  for (int j = 0; j < BndConds.size(); ++j)
691  {
692  if (BndConds[j]->GetBoundaryConditionType() ==
693  SpatialDomains::eNeumann)
694  {
695  projectfield = i;
696  break;
697  }
698  }
699  if (projectfield != -1)
700  {
701  break;
702  }
703  }
704  if (projectfield == -1)
705  {
706  cout << "using first field to project non-linear forcing which "
707  "imposes a Dirichlet condition"
708  << endl;
709  projectfield = 0;
710  }
711  }
712 
713  // Shift m_vwiForcing in case of relaxation
714  Vmath::Vcopy(ncoeffs, m_vwiForcing[0], 1, m_vwiForcing[2], 1);
715  Vmath::Vcopy(ncoeffs, m_vwiForcing[1], 1, m_vwiForcing[3], 1);
716 
717  // determine inverse of area normalised field.
718  m_wavePressure->GetPlane(0)->BwdTrans(
719  m_wavePressure->GetPlane(0)->GetCoeffs(),
720  m_wavePressure->GetPlane(0)->UpdatePhys());
721  m_wavePressure->GetPlane(1)->BwdTrans(
722  m_wavePressure->GetPlane(1)->GetCoeffs(),
723  m_wavePressure->GetPlane(1)->UpdatePhys());
724  NekDouble invnorm;
725 
726  if (m_useLinfPressureNorm)
727  {
728  Vmath::Vmul(npts, m_wavePressure->GetPlane(0)->UpdatePhys(), 1,
729  m_wavePressure->GetPlane(0)->UpdatePhys(), 1, der1, 1);
730  Vmath::Vvtvp(npts, m_wavePressure->GetPlane(1)->UpdatePhys(), 1,
731  m_wavePressure->GetPlane(1)->UpdatePhys(), 1, der1, 1,
732  der1, 1);
733  Vmath::Vsqrt(npts, der1, 1, der1, 1);
734 
735  NekDouble Linf = Vmath::Vmax(npts, der1, 1);
736 
737  invnorm = 1.0 / Linf;
738  }
739  else
740  {
741  // Determine normalisation of pressure so that |P|/A = 1
742  NekDouble l2;
743  l2 = m_wavePressure->GetPlane(0)->L2(
744  m_wavePressure->GetPlane(0)->GetPhys());
745  invnorm = l2 * l2;
746  l2 = m_wavePressure->GetPlane(1)->L2(
747  m_wavePressure->GetPlane(1)->GetPhys());
748  invnorm += l2 * l2;
749  Vmath::Fill(2 * npts, 1.0, der1, 1);
750  NekDouble area = m_waveVelocities[0]->GetPlane(0)->Integral(der1);
751  cout << "Area: " << area << endl;
752  invnorm = sqrt(area / invnorm);
753  }
754 
755  // Get hold of arrays.
756  m_waveVelocities[0]->GetPlane(0)->BwdTrans(
757  m_waveVelocities[0]->GetPlane(0)->GetCoeffs(),
758  m_waveVelocities[0]->GetPlane(0)->UpdatePhys());
759  Array<OneD, NekDouble> u_real =
760  m_waveVelocities[0]->GetPlane(0)->UpdatePhys();
761  Vmath::Smul(npts, invnorm, u_real, 1, u_real, 1);
762  m_waveVelocities[0]->GetPlane(1)->BwdTrans(
763  m_waveVelocities[0]->GetPlane(1)->GetCoeffs(),
764  m_waveVelocities[0]->GetPlane(1)->UpdatePhys());
765  Array<OneD, NekDouble> u_imag =
766  m_waveVelocities[0]->GetPlane(1)->UpdatePhys();
767  Vmath::Smul(npts, invnorm, u_imag, 1, u_imag, 1);
768  m_waveVelocities[1]->GetPlane(0)->BwdTrans(
769  m_waveVelocities[1]->GetPlane(0)->GetCoeffs(),
770  m_waveVelocities[1]->GetPlane(0)->UpdatePhys());
771  Array<OneD, NekDouble> v_real =
772  m_waveVelocities[1]->GetPlane(0)->UpdatePhys();
773  Vmath::Smul(npts, invnorm, v_real, 1, v_real, 1);
774  m_waveVelocities[1]->GetPlane(1)->BwdTrans(
775  m_waveVelocities[1]->GetPlane(1)->GetCoeffs(),
776  m_waveVelocities[1]->GetPlane(1)->UpdatePhys());
777  Array<OneD, NekDouble> v_imag =
778  m_waveVelocities[1]->GetPlane(1)->UpdatePhys();
779  Vmath::Smul(npts, invnorm, v_imag, 1, v_imag, 1);
780 
781  // normalise wave for output
782  Vmath::Smul(2 * ncoeffs, invnorm, m_waveVelocities[0]->UpdateCoeffs(),
783  1, m_waveVelocities[0]->UpdateCoeffs(), 1);
784  Vmath::Smul(2 * ncoeffs, invnorm, m_waveVelocities[1]->UpdateCoeffs(),
785  1, m_waveVelocities[1]->UpdateCoeffs(), 1);
786  Vmath::Smul(2 * ncoeffs, invnorm, m_waveVelocities[2]->UpdateCoeffs(),
787  1, m_waveVelocities[2]->UpdateCoeffs(), 1);
788 
789  // dump field
790  {
791  std::vector<std::string> variables(3);
792  std::vector<Array<OneD, NekDouble>> outfield(3);
793  variables[0] = "u_w";
794  variables[1] = "v_w";
795  variables[2] = "w_w";
796  outfield[0] = m_waveVelocities[0]->UpdateCoeffs();
797  outfield[1] = m_waveVelocities[1]->UpdateCoeffs();
798  outfield[2] = m_waveVelocities[2]->UpdateCoeffs();
799  std::string outname = m_sessionName + "_wave.fld";
800  m_solverRoll->WriteFld(outname, m_waveVelocities[0], outfield,
801  variables);
802  }
803 
804 #if 1
805  int ncoeffs_p = m_wavePressure->GetPlane(0)->GetNcoeffs();
806  Vmath::Smul(ncoeffs_p, invnorm,
807  m_wavePressure->GetPlane(0)->UpdateCoeffs(), 1,
808  m_wavePressure->GetPlane(0)->UpdateCoeffs(), 1);
809  Vmath::Smul(ncoeffs_p, invnorm,
810  m_wavePressure->GetPlane(1)->UpdateCoeffs(), 1,
811  m_wavePressure->GetPlane(1)->UpdateCoeffs(), 1);
812 #else
813  m_wavePressure->GetPlane(0)->BwdTrans(
814  m_wavePressure->GetPlane(0)->GetCoeffs(),
815  m_wavePressure->GetPlane(0)->UpdatePhys());
816  Vmath::Smul(npts, invnorm, m_wavePressure->GetPlane(0)->UpdatePhys(), 1,
817  m_wavePressure->GetPlane(0)->UpdatePhys(), 1);
818  m_wavePressure->GetPlane(0)->FwdTrans(
819  m_wavePressure->GetPlane(0)->UpdatePhys(),
820  m_wavePressure->GetPlane(0)->UpdateCoeffs());
821 
822  m_wavePressure->GetPlane(1)->BwdTrans(
823  m_wavePressure->GetPlane(1)->GetCoeffs(),
824  m_wavePressure->GetPlane(1)->UpdatePhys());
825  Vmath::Smul(npts, invnorm, m_wavePressure->GetPlane(1)->UpdatePhys(), 1,
826  m_wavePressure->GetPlane(1)->UpdatePhys(), 1);
827  m_wavePressure->GetPlane(1)->FwdTrans(
828  m_wavePressure->GetPlane(1)->UpdatePhys(),
829  m_wavePressure->GetPlane(1)->UpdateCoeffs());
830 #endif
831 
832  // Calculate non-linear terms for x and y directions
833  // d/dx(u u* + u* u)
834  Vmath::Vmul(npts, u_real, 1, u_real, 1, val, 1);
835  Vmath::Vvtvp(npts, u_imag, 1, u_imag, 1, val, 1, val, 1);
836  Vmath::Smul(npts, 2.0, val, 1, val, 1);
837  m_waveVelocities[0]->GetPlane(0)->PhysDeriv(0, val, der1);
838 
839  // d/dy(v u* + v* u)
840  Vmath::Vmul(npts, u_real, 1, v_real, 1, val, 1);
841  Vmath::Vvtvp(npts, u_imag, 1, v_imag, 1, val, 1, val, 1);
842  Vmath::Smul(npts, 2.0, val, 1, val, 1);
843  m_waveVelocities[0]->GetPlane(0)->PhysDeriv(1, val, der2);
844 
845  Vmath::Vadd(npts, der1, 1, der2, 1, der1, 1);
846 #if 1
847  m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(der1,
848  m_vwiForcing[0]);
849 #else
850  m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(
851  der1, m_vwiForcing[0]);
852 #endif
853  Vmath::Smul(ncoeffs, -m_waveForceMag[0], m_vwiForcing[0], 1,
854  m_vwiForcing[0], 1);
855 
856  // d/dx(u v* + u* v)
857  m_waveVelocities[0]->GetPlane(0)->PhysDeriv(0, val, der1);
858 
859  // d/dy(v v* + v* v)
860  Vmath::Vmul(npts, v_real, 1, v_real, 1, val, 1);
861  Vmath::Vvtvp(npts, v_imag, 1, v_imag, 1, val, 1, val, 1);
862  Vmath::Smul(npts, 2.0, val, 1, val, 1);
863  m_waveVelocities[0]->GetPlane(0)->PhysDeriv(1, val, der2);
864 
865  Vmath::Vadd(npts, der1, 1, der2, 1, der1, 1);
866 
867 #if 1
868  m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(der1,
869  m_vwiForcing[1]);
870 #else
871  m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(
872  der1, m_vwiForcing[1]);
873 #endif
874 
875  Vmath::Smul(ncoeffs, -m_waveForceMag[0], m_vwiForcing[1], 1,
876  m_vwiForcing[1], 1);
877 
878  // by default the symmetrization is on
879  bool symm = true;
880  m_sessionVWI->MatchSolverInfo("Symmetrization", "True", symm, true);
881 #if 0
882  if(symm== true )
883  {
884 
885  // Symmetrise forcing
886  //-> Get coordinates
887  Array<OneD, NekDouble> coord(2);
888  Array<OneD, NekDouble> coord_x(npts);
889  Array<OneD, NekDouble> coord_y(npts);
890 
891  //-> Impose symmetry (x -> -x + Lx/2, y-> -y) on coordinates
892  m_waveVelocities[0]->GetPlane(0)->GetCoords(coord_x,coord_y);
893  NekDouble xmax = Vmath::Vmax(npts,coord_x,1);
894  Vmath::Neg(npts,coord_x,1);
895  Vmath::Sadd(npts,xmax,coord_x,1,coord_x,1);
896  Vmath::Neg(npts,coord_y,1);
897 
898  int i, physoffset;
899 
900  //-> Obtain list of expansion element ids for each point.
901  Array<OneD, int> Eid(npts);
902  // This search may not be necessary every iteration
903  for(i = 0; i < npts; ++i)
904  {
905  coord[0] = coord_x[i];
906  coord[1] = coord_y[i];
907 
908  // Note this will not quite be symmetric.
909  Eid[i] = m_waveVelocities[0]->GetPlane(0)->GetExpIndex(coord,1e-6);
910  }
911 
912  // Interpolate field 0
913  m_waveVelocities[0]->GetPlane(0)->BwdTrans(m_vwiForcing[0],der1);
914  for(i = 0; i < npts; ++i)
915  {
916  physoffset = m_waveVelocities[0]->GetPlane(0)->GetPhys_Offset(Eid[i]);
917  coord[0] = coord_x[i];
918  coord[1] = coord_y[i];
919  der2 [i] = m_waveVelocities[0]->GetPlane(0)->GetExp(Eid[i])->PhysEvaluate(coord, der1 + physoffset);
920  }
921  //-> Average field 0
922  Vmath::Vsub(npts,der1,1,der2,1,der2,1);
923  Vmath::Smul(npts,0.5,der2,1,der2,1);
924 #if 1
925  m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(der2,m_vwiForcing[0]);
926 #else
927  m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(der2, m_vwiForcing[0]);
928 #endif
929 
930  //-> Interpoloate field 1
931  m_waveVelocities[0]->GetPlane(0)->BwdTrans(m_vwiForcing[1],der1);
932  for(i = 0; i < npts; ++i)
933  {
934  physoffset = m_waveVelocities[0]->GetPlane(0)->GetPhys_Offset(Eid[i]);
935  coord[0] = coord_x[i];
936  coord[1] = coord_y[i];
937  der2[i] = m_waveVelocities[0]->GetPlane(0)->GetExp(Eid[i])->PhysEvaluate(coord, der1 + physoffset);
938  }
939 
940  //-> Average field 1
941  Vmath::Vsub(npts,der1,1,der2,1,der2,1);
942  Vmath::Smul(npts,0.5,der2,1,der2,1);
943 #if 1
944  m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(der2,m_vwiForcing[1]);
945 #else
946  m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(der2, m_vwiForceing[1]);
947 #endif
948  }
949 #else
950  int i;
951  if (symm == true)
952  {
953  cout << "symmetrization is active" << endl;
954  static Array<OneD, int> index = GetReflectionIndex();
955 
956  m_waveVelocities[0]->GetPlane(0)->BwdTrans(m_vwiForcing[0], der1);
957 
958  for (i = 0; i < npts; ++i)
959  {
960  if (index[i] != -1)
961  {
962  val[i] = 0.5 * (der1[i] - der1[index[i]]);
963  }
964  }
965 #if 1
966  m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(
967  val, m_vwiForcing[0]);
968 #else
969  m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(
970  val, m_vwiForcing[0]);
971 #endif
972 
973  m_waveVelocities[0]->GetPlane(0)->BwdTrans(m_vwiForcing[1], der2);
974  for (i = 0; i < npts; ++i)
975  {
976  if (index[i] != -1)
977  {
978  val[i] = 0.5 * (der2[i] - der2[index[i]]);
979  }
980  }
981 #if 1
982  m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(
983  val, m_vwiForcing[1]);
984 #else
985  m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(
986  val, m_vwiForcing[1]);
987 #endif
988  }
989 
990  Vmath::Vmul(npts, der1, 1, der1, 1, val, 1);
991  Vmath::Vvtvp(npts, der2, 1, der2, 1, val, 1, val, 1);
992  Vmath::Vsqrt(npts, val, 1, val, 1);
993  cout << "F_Linf: " << Vmath::Vmax(npts, val, 1) << endl;
994 
995 #endif
996 
997  if (m_vwiRelaxation)
998  {
999  Vmath::Smul(ncoeffs, 1.0 - m_vwiRelaxation, m_vwiForcing[0], 1,
1000  m_vwiForcing[0], 1);
1001  Vmath::Svtvp(ncoeffs, m_vwiRelaxation, m_vwiForcing[2], 1,
1002  m_vwiForcing[0], 1, m_vwiForcing[0], 1);
1003 
1004  Vmath::Smul(ncoeffs, 1.0 - m_vwiRelaxation, m_vwiForcing[1], 1,
1005  m_vwiForcing[1], 1);
1006  Vmath::Svtvp(ncoeffs, m_vwiRelaxation, m_vwiForcing[3], 1,
1007  m_vwiForcing[1], 1, m_vwiForcing[1], 1);
1008  }
1009 
1010  if (m_deltaFcnApprox)
1011  {
1012  for (int j = 0; j < 2; ++j)
1013  {
1014 
1015  m_waveVelocities[projectfield]->GetPlane(0)->BwdTrans(
1016  m_vwiForcing[j], der1);
1017  for (int i = 0; i < npts; ++i)
1018  {
1019  der1[i] *= exp(-streak[i] * streak[i] / m_deltaFcnDecay);
1020  }
1021  m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(
1022  der1, m_vwiForcing[j]);
1023  }
1024  }
1025 
1026  // dump output
1027  std::vector<std::string> variables(4);
1028  std::vector<Array<OneD, NekDouble>> outfield(4);
1029  variables[0] = "u";
1030  variables[1] = "v";
1031  variables[2] = "pr";
1032  variables[3] = "pi";
1033  outfield[0] = m_vwiForcing[0];
1034  outfield[1] = m_vwiForcing[1];
1035  Array<OneD, NekDouble> soln(npts, 0.0);
1036  m_wavePressure->GetPlane(0)->BwdTrans(
1037  m_wavePressure->GetPlane(0)->GetCoeffs(),
1038  m_wavePressure->GetPlane(0)->UpdatePhys());
1039  outfield[2] = Array<OneD, NekDouble>(ncoeffs);
1040  m_waveVelocities[0]->GetPlane(0)->FwdTransLocalElmt(
1041  m_wavePressure->GetPlane(0)->GetPhys(), outfield[2]);
1042  m_wavePressure->GetPlane(1)->BwdTrans(
1043  m_wavePressure->GetPlane(1)->GetCoeffs(),
1044  m_wavePressure->GetPlane(1)->UpdatePhys());
1045 
1046  Vmath::Vmul(npts, m_wavePressure->GetPlane(0)->UpdatePhys(), 1,
1047  m_wavePressure->GetPlane(0)->UpdatePhys(), 1, val, 1);
1048  Vmath::Vvtvp(npts, m_wavePressure->GetPlane(1)->UpdatePhys(), 1,
1049  m_wavePressure->GetPlane(1)->UpdatePhys(), 1, val, 1, val,
1050  1);
1051  cout << "int P^2: " << m_wavePressure->GetPlane(0)->Integral(val)
1052  << endl;
1053  Vmath::Vsqrt(npts, val, 1, val, 1);
1054  cout << "PLinf: " << Vmath::Vmax(npts, val, 1) << endl;
1055 
1056  outfield[3] = Array<OneD, NekDouble>(ncoeffs);
1057  m_waveVelocities[1]->GetPlane(0)->FwdTransLocalElmt(
1058  m_wavePressure->GetPlane(1)->GetPhys(), outfield[3]);
1059 
1060  std::string outname = m_sessionName + ".vwi";
1061 
1062  m_solverRoll->WriteFld(outname, m_waveVelocities[0]->GetPlane(0),
1063  outfield, variables);
1064  }
1065 }
1066 
1067 void VortexWaveInteraction::CalcL2ToLinfPressure(void)
1068 {
1069 
1070  ExecuteWave();
1071 
1072  m_wavePressure->GetPlane(0)->BwdTrans(
1073  m_wavePressure->GetPlane(0)->GetCoeffs(),
1074  m_wavePressure->GetPlane(0)->UpdatePhys());
1075  m_wavePressure->GetPlane(1)->BwdTrans(
1076  m_wavePressure->GetPlane(1)->GetCoeffs(),
1077  m_wavePressure->GetPlane(1)->UpdatePhys());
1078 
1079  int npts = m_waveVelocities[0]->GetPlane(0)->GetNpoints();
1080  NekDouble Linf;
1081  Array<OneD, NekDouble> val(2 * npts, 0.0);
1082 
1083  Vmath::Vmul(npts, m_wavePressure->GetPlane(0)->UpdatePhys(), 1,
1084  m_wavePressure->GetPlane(0)->UpdatePhys(), 1, val, 1);
1085  Vmath::Vvtvp(npts, m_wavePressure->GetPlane(1)->UpdatePhys(), 1,
1086  m_wavePressure->GetPlane(1)->UpdatePhys(), 1, val, 1, val, 1);
1087  cout << "int P^2 " << m_wavePressure->GetPlane(0)->Integral(val) << endl;
1088  Vmath::Vsqrt(npts, val, 1, val, 1);
1089 
1090  Linf = Vmath::Vmax(npts, val, 1);
1091  cout << "Linf: " << Linf << endl;
1092 
1093  NekDouble l2, norm;
1094  l2 =
1095  m_wavePressure->GetPlane(0)->L2(m_wavePressure->GetPlane(0)->GetPhys());
1096  norm = l2 * l2;
1097  l2 =
1098  m_wavePressure->GetPlane(1)->L2(m_wavePressure->GetPlane(1)->GetPhys());
1099  norm += l2 * l2;
1100 
1101  Vmath::Fill(npts, 1.0, val, 1);
1102  NekDouble area = m_waveVelocities[0]->GetPlane(0)->Integral(val);
1103 
1104  l2 = sqrt(norm / area);
1105 
1106  cout << "L2: " << l2 << endl;
1107 
1108  cout << "Ratio Linf/L2: " << Linf / l2 << endl;
1109 }
1110 
1111 void VortexWaveInteraction::SaveFile(string file, string dir, int n)
1112 {
1113  static map<string, int> opendir;
1114 
1115  if (opendir.find(dir) == opendir.end())
1116  {
1117  // make directory and presume will fail if it already exists
1118  string mkdir = "mkdir " + dir;
1119  ASSERTL0(system(mkdir.c_str()) == 0,
1120  "Failed to make directory '" + dir + "'");
1121 
1122  opendir[dir] = 1;
1123  }
1124 
1125  string savefile =
1126  dir + "/" + file + "." + boost::lexical_cast<std::string>(n);
1127  string syscall = "cp -f " + file + " " + savefile;
1128 
1129  ASSERTL0(system(syscall.c_str()) == 0, syscall.c_str());
1130 }
1131 
1132 void VortexWaveInteraction::MoveFile(string file, string dir, int n)
1133 {
1134  static map<string, int> opendir;
1135 
1136  if (opendir.find(dir) == opendir.end())
1137  {
1138  // make directory and presume will fail if it already exists
1139  string mkdir = "mkdir " + dir;
1140  ASSERTL0(system(mkdir.c_str()) == 0,
1141  "Failed to make directory '" + dir + "'");
1142  opendir[dir] = 1;
1143  }
1144 
1145  string savefile =
1146  dir + "/" + file + "." + boost::lexical_cast<std::string>(n);
1147  string syscall = "mv -f " + file + " " + savefile;
1148 
1149  ASSERTL0(system(syscall.c_str()) == 0, syscall.c_str());
1150 }
1151 
1152 void VortexWaveInteraction::CopyFile(string file1end, string file2end)
1153 {
1154  string cpfile1 = m_sessionName + file1end;
1155  string cpfile2 = m_sessionName + file2end;
1156  string syscall = "cp -f " + cpfile1 + " " + cpfile2;
1157 
1158  if (system(syscall.c_str()))
1159  {
1160  ASSERTL0(false, syscall.c_str());
1161  }
1162 }
1163 
1164 void VortexWaveInteraction::AppendEvlToFile(string file, int n)
1165 {
1166  FILE *fp;
1167  fp = fopen(file.c_str(), "a");
1168  fprintf(fp, "%d: %lf %16.12le %16.12le\n", n, m_alpha[0],
1169  m_leading_real_evl[0], m_leading_imag_evl[0]);
1170  fclose(fp);
1171 }
1172 
1173 void VortexWaveInteraction::AppendEvlToFile(string file, NekDouble WaveForceMag)
1174 {
1175  FILE *fp;
1176  fp = fopen(file.c_str(), "a");
1177  fprintf(fp, "%lf %lf %16.12le %16.12le\n", WaveForceMag, m_alpha[0],
1178  m_leading_real_evl[0], m_leading_imag_evl[0]);
1179  fclose(fp);
1180 }
1181 
1182 void VortexWaveInteraction::SaveLoopDetails(std::string SaveDir, int i)
1183 
1184 {
1185  // Save NS restart file
1186  SaveFile(m_sessionName + ".rst", SaveDir, i + 1);
1187  // Save Streak Solution
1188  SaveFile(m_sessionName + "_streak.fld", SaveDir, i);
1189  // Save Wave solution output
1190  SaveFile(m_sessionName + ".evl", SaveDir, i);
1191  SaveFile(m_sessionName + "_eig_0", SaveDir, i);
1192  // Save new field file of eigenvalue
1193  SaveFile(m_sessionName + ".fld", SaveDir, i);
1194  if (!(m_sessionVWI->DefinesSolverInfo("INTERFACE")))
1195  {
1196  // Save new forcing file
1197  SaveFile(m_sessionName + ".vwi", SaveDir, i + 1);
1198  }
1199 }
1200 
1201 void VortexWaveInteraction::ExecuteLoop(bool CalcWaveForce)
1202 {
1203 
1204  // the global info has to be written in the
1205  // roll session file to use the interface loop
1206  if (m_sessionRoll->DefinesSolverInfo("INTERFACE"))
1207  {
1208  static int cnt = 0;
1209  bool skiprollstreak = false;
1210  if (cnt == 0 && m_sessionVWI->GetParameter("rollstreakfromit") == 1)
1211  {
1212  skiprollstreak = true;
1213  cout << "skip roll-streak at the first iteration" << endl;
1214  }
1215 
1216  if (skiprollstreak != true)
1217  {
1218 
1219  LibUtilities::NekManager<
1220  MultiRegions::GlobalLinSysKey,
1221  MultiRegions::GlobalLinSys>::EnableManagement("GlobalLinSys");
1222  ExecuteRoll();
1223  LibUtilities::NekManager<
1224  MultiRegions::GlobalLinSysKey,
1225  MultiRegions::GlobalLinSys>::DisableManagement("GlobalLinSys");
1226 #ifndef _WIN32
1227  sleep(3);
1228 #endif
1229  ExecuteStreak();
1230 #ifndef _WIN32
1231  sleep(3);
1232 #endif
1233  }
1234 
1235  string syscall;
1236  char c[16] = "";
1237  string movedmesh = m_sessionName + "_advPost_moved.xml";
1238  string movedinterpmesh = m_sessionName + "_interp_moved.xml";
1239  // rewrite the Rollsessionfile (we start from the waleffe forcing)
1240  // string meshbndjumps = m_sessionName +"_bndjumps.xml";
1241  // if(cnt==0)
1242  //{
1243  // take the conditions tag from meshbndjumps and copy into
1244  // the rolls session file
1245  //}
1246 
1247  sprintf(c, "%d", cnt);
1248  // save old roll solution
1249  string oldroll = m_sessionName + "_roll_" + c + ".fld";
1250  syscall = "cp -f " + m_sessionName + "-Base.fld" + " " + oldroll;
1251  cout << syscall.c_str() << endl;
1252  if (system(syscall.c_str()))
1253  {
1254  ASSERTL0(false, syscall.c_str());
1255  }
1256  // define file names
1257  string filePost = m_sessionName + "_advPost.xml";
1258  string filestreak = m_sessionName + "_streak.fld";
1259  string filewave = m_sessionName + "_wave.fld";
1260  string filewavepressure = m_sessionName + "_wave_p_split.fld";
1261  string fileinterp = m_sessionName + "_interp.xml";
1262  string interpstreak = m_sessionName + "_interpstreak_" + c + ".fld";
1263  string interwavepressure =
1264  m_sessionName + "_wave_p_split_interp_" + c + ".fld";
1265  char alpchar[16] = "";
1266  cout << "alpha = " << m_alpha[0] << endl;
1267  sprintf(alpchar, "%f", m_alpha[0]);
1268 
1269  if (m_sessionVWI->GetSolverInfo("INTERFACE") != "phase")
1270  {
1271  cout << "zerophase" << endl;
1272 
1273  syscall = "../../utilities/PostProcessing/Extras/MoveMesh " +
1274  filePost + " " + filestreak + " " + fileinterp + " " +
1275  alpchar;
1276 
1277  cout << syscall.c_str() << endl;
1278  if (system(syscall.c_str()))
1279  {
1280  ASSERTL0(false, syscall.c_str());
1281  }
1282 
1283  // move the advPost mesh (remark update alpha!!!)
1284  syscall = "../../utilities/PostProcessing/Extras/MoveMesh " +
1285  filePost + " " + filestreak + " " + filePost + " " +
1286  alpchar;
1287  cout << syscall.c_str() << endl;
1288  if (system(syscall.c_str()))
1289  {
1290  ASSERTL0(false, syscall.c_str());
1291  }
1292 
1293  // save oldstreak
1294  string oldstreak = m_sessionName + "_streak_" + c + ".fld";
1295  syscall = "cp -f " + filestreak + " " + oldstreak;
1296  cout << syscall.c_str() << endl;
1297  if (system(syscall.c_str()))
1298  {
1299  ASSERTL0(false, syscall.c_str());
1300  }
1301 
1302  // interpolate the streak field into the new mesh
1303  string movedmesh = m_sessionName + "_advPost_moved.xml";
1304  string movedinterpmesh = m_sessionName + "_interp_moved.xml";
1305 
1306  // create the interp streak
1307 
1308  syscall = "../../utilities/PostProcessing/Extras/FieldToField " +
1309  fileinterp + " " + filestreak + " " + movedinterpmesh +
1310  " " + interpstreak;
1311 
1312  cout << syscall.c_str() << endl;
1313  if (system(syscall.c_str()))
1314  {
1315  ASSERTL0(false, syscall.c_str());
1316  }
1317 
1318  // save the old mesh
1319  string meshfile = m_sessionName + ".xml";
1320  string meshold = m_sessionName + "_" + c + ".xml";
1321  syscall = "cp -f " + meshfile + " " + meshold;
1322  cout << syscall.c_str() << endl;
1323  if (system(syscall.c_str()))
1324  {
1325  ASSERTL0(false, syscall.c_str());
1326  }
1327 
1328  // overwriting the meshfile with the new mesh
1329  syscall = "cp -f " + movedmesh + " " + meshfile;
1330  cout << syscall.c_str() << endl;
1331  if (system(syscall.c_str()))
1332  {
1333  ASSERTL0(false, syscall.c_str());
1334  }
1335 
1336  // overwriting the streak file!!
1337  syscall = "cp -f " + interpstreak + " " + filestreak;
1338  cout << syscall.c_str() << endl;
1339  if (system(syscall.c_str()))
1340  {
1341  ASSERTL0(false, syscall.c_str());
1342  }
1343 
1344  // calculate the wave
1345  ExecuteWave();
1346 
1347  // save the wave field:
1348  string oldwave = m_sessionName + "_wave_" + c + ".fld";
1349  syscall = "cp -f " + m_sessionName + ".fld" + " " + oldwave;
1350  cout << syscall.c_str() << endl;
1351  if (system(syscall.c_str()))
1352  {
1353  ASSERTL0(false, syscall.c_str());
1354  }
1355 
1356  // save old jump conditions:
1357  string ujump = m_sessionName + "_u_5.bc";
1358  syscall = "cp -f " + ujump + " " + m_sessionName + "_u_5.bc_" + c;
1359  cout << syscall.c_str() << endl;
1360  if (system(syscall.c_str()))
1361  {
1362  ASSERTL0(false, syscall.c_str());
1363  }
1364 
1365  string vjump = m_sessionName + "_v_5.bc";
1366  syscall = "cp -f " + vjump + " " + m_sessionName + "_v_5.bc_" + c;
1367  cout << syscall.c_str() << endl;
1368  if (system(syscall.c_str()))
1369  {
1370  ASSERTL0(false, syscall.c_str());
1371  }
1372  cnt++;
1373 
1374  // use relaxation
1375  if (GetVWIIterationType() !=
1376  eFixedWaveForcingWithSubIterationOnAlpha)
1377  {
1378  // the critical layer should be the bnd region 3
1379  // int reg =3;
1380  // FileRelaxation(reg);
1381  }
1382  char c1[16] = "";
1383  sprintf(c1, "%d", cnt);
1384  // calculate the jump conditions
1385  string wavefile = m_sessionName + ".fld";
1386  syscall = "../../utilities/PostProcessing/Extras/FldCalcBCs " +
1387  movedmesh + " " + wavefile + " " + interpstreak +
1388  "> data" + c1;
1389  cout << syscall.c_str() << endl;
1390  if (system(syscall.c_str()))
1391  {
1392  ASSERTL0(false, syscall.c_str());
1393  }
1394 
1395  // move the new name_interp_moved.xml into name_interp.xml
1396  syscall = "cp -f " + movedinterpmesh + " " + fileinterp;
1397  cout << syscall.c_str() << endl;
1398  if (system(syscall.c_str()))
1399  {
1400  ASSERTL0(false, syscall.c_str());
1401  }
1402  // move the new name_advPost_moved.xml into name_advPost.xml
1403  syscall = "cp -f " + movedmesh + " " + filePost;
1404  cout << syscall.c_str() << endl;
1405  if (system(syscall.c_str()))
1406  {
1407  ASSERTL0(false, syscall.c_str());
1408  }
1409  }
1410  else if (m_sessionVWI->GetSolverInfo("INTERFACE") == "phase")
1411  {
1412  cout << "phase" << endl;
1413  // determine cr:
1414  NekDouble cr;
1415  string cr_str;
1416  stringstream st;
1417 
1418  // calculate the wave
1419  ExecuteWave();
1420 
1421  // save oldstreak
1422  string oldstreak = m_sessionName + "_streak_" + c + ".fld";
1423  syscall = "cp -f " + filestreak + " " + oldstreak;
1424  cout << syscall.c_str() << endl;
1425  if (system(syscall.c_str()))
1426  {
1427  ASSERTL0(false, syscall.c_str());
1428  }
1429 
1430  // save wave
1431  syscall = "cp -f " + m_sessionName + ".fld" + " " + filewave;
1432  cout << syscall.c_str() << endl;
1433  if (system(syscall.c_str()))
1434  {
1435  ASSERTL0(false, syscall.c_str());
1436  }
1437  // save the old mesh
1438  string meshfile = m_sessionName + ".xml";
1439  string meshold = m_sessionName + "_" + c + ".xml";
1440  syscall = "cp -f " + meshfile + " " + meshold;
1441  cout << syscall.c_str() << endl;
1442  if (system(syscall.c_str()))
1443  {
1444  ASSERTL0(false, syscall.c_str());
1445  }
1446 
1447  // save the oldwave field:
1448  string oldwave = m_sessionName + "_wave_" + c + ".fld";
1449  syscall = "cp -f " + m_sessionName + ".fld" + " " + oldwave;
1450  cout << syscall.c_str() << endl;
1451  if (system(syscall.c_str()))
1452  {
1453  ASSERTL0(false, syscall.c_str());
1454  }
1455 
1456  // save old jump conditions:
1457  string ujump = m_sessionName + "_u_5.bc";
1458  syscall = "cp -f " + ujump + " " + m_sessionName + "_u_5.bc_" + c;
1459  cout << syscall.c_str() << endl;
1460  if (system(syscall.c_str()))
1461  {
1462  ASSERTL0(false, syscall.c_str());
1463  }
1464 
1465  string vjump = m_sessionName + "_v_5.bc";
1466  syscall = "cp -f " + vjump + " " + m_sessionName + "_v_5.bc_" + c;
1467  cout << syscall.c_str() << endl;
1468  if (system(syscall.c_str()))
1469  {
1470  ASSERTL0(false, syscall.c_str());
1471  }
1472  cnt++;
1473 
1474  cr = m_leading_imag_evl[0] / m_alpha[0];
1475  st << cr;
1476  cr_str = st.str();
1477  cout << "cr=" << cr_str << endl;
1478  // NB -g or NOT!!!
1479  // move the mesh around the critical layer
1480  syscall = "../../utilities/PostProcessing/Extras/MoveMesh " +
1481  filePost + " " + filestreak + " " + fileinterp + " " +
1482  alpchar + " " + cr_str;
1483 
1484  cout << syscall.c_str() << endl;
1485  if (system(syscall.c_str()))
1486  {
1487  ASSERTL0(false, syscall.c_str());
1488  }
1489  // NB -g or NOT!!!
1490  // move the advPost mesh (remark update alpha!!!)
1491  syscall = "../../utilities/PostProcessing/Extras/MoveMesh " +
1492  filePost + " " + filestreak + " " + filePost + " " +
1493  alpchar + " " + cr_str;
1494  cout << syscall.c_str() << endl;
1495  if (system(syscall.c_str()))
1496  {
1497  ASSERTL0(false, syscall.c_str());
1498  }
1499 
1500  // interp streak into the new mesh
1501  syscall = "../../utilities/PostProcessing/Extras/FieldToField " +
1502  fileinterp + " " + filestreak + " " + movedinterpmesh +
1503  " " + interpstreak;
1504 
1505  cout << syscall.c_str() << endl;
1506  if (system(syscall.c_str()))
1507  {
1508  ASSERTL0(false, syscall.c_str());
1509  }
1510 
1511  // split wave sol
1512  syscall = "../../utilities/PostProcessing/Extras/SplitFld " +
1513  filePost + " " + filewave;
1514 
1515  cout << syscall.c_str() << endl;
1516  if (system(syscall.c_str()))
1517  {
1518  ASSERTL0(false, syscall.c_str());
1519  }
1520  // interp wave
1521  syscall = "../../utilities/PostProcessing/Extras/FieldToField " +
1522  filePost + " " + filewavepressure + " " + movedmesh +
1523  " " + interwavepressure;
1524 
1525  cout << syscall.c_str() << endl;
1526  if (system(syscall.c_str()))
1527  {
1528  ASSERTL0(false, syscall.c_str());
1529  }
1530 
1531  // use relaxation
1532  if (GetVWIIterationType() !=
1533  eFixedWaveForcingWithSubIterationOnAlpha)
1534  {
1535  // the critical layer should be the bnd region 3
1536  // int reg =3;
1537  // FileRelaxation(reg);
1538  }
1539  char c1[16] = "";
1540  sprintf(c1, "%d", cnt);
1541 
1542  // cp wavepressure to m_sessionName.fld(to get
1543  // the right bcs names using FldCalcBCs)
1544  syscall =
1545  "cp -f " + interwavepressure + " " + m_sessionName + ".fld";
1546  cout << syscall.c_str() << endl;
1547  if (system(syscall.c_str()))
1548  {
1549  ASSERTL0(false, syscall.c_str());
1550  }
1551 
1552  // calculate the jump conditions
1553  // NB -g or NOT!!!
1554  syscall = "../../utilities/PostProcessing/Extras/FldCalcBCs " +
1555  movedmesh + " " + m_sessionName + ".fld" + " " +
1556  interpstreak + "> data" + c1;
1557  cout << syscall.c_str() << endl;
1558  if (system(syscall.c_str()))
1559  {
1560  ASSERTL0(false, syscall.c_str());
1561  }
1562 
1563  // overwriting the meshfile with the new mesh
1564  syscall = "cp -f " + movedmesh + " " + meshfile;
1565  cout << syscall.c_str() << endl;
1566  if (system(syscall.c_str()))
1567  {
1568  ASSERTL0(false, syscall.c_str());
1569  }
1570 
1571  // overwriting the streak file!! (maybe is useless)
1572  syscall = "cp -f " + interpstreak + " " + filestreak;
1573  cout << syscall.c_str() << endl;
1574  if (system(syscall.c_str()))
1575  {
1576  ASSERTL0(false, syscall.c_str());
1577  }
1578  // move the new name_interp_moved.xml into name_interp.xml
1579  syscall = "cp -f " + movedinterpmesh + " " + fileinterp;
1580  cout << syscall.c_str() << endl;
1581  if (system(syscall.c_str()))
1582  {
1583  ASSERTL0(false, syscall.c_str());
1584  }
1585  // move the new name_advPost_moved.xml into name_advPost.xml
1586  syscall = "cp -f " + movedmesh + " " + filePost;
1587  cout << syscall.c_str() << endl;
1588  if (system(syscall.c_str()))
1589  {
1590  ASSERTL0(false, syscall.c_str());
1591  }
1592  }
1593  }
1594  else
1595  {
1596  LibUtilities::NekManager<
1597  MultiRegions::GlobalLinSysKey,
1598  MultiRegions::GlobalLinSys>::EnableManagement("GlobalLinSys");
1599  ExecuteRoll();
1600  LibUtilities::NekManager<
1601  MultiRegions::GlobalLinSysKey,
1602  MultiRegions::GlobalLinSys>::DisableManagement("GlobalLinSys");
1603 
1604 #ifndef _WIN32
1605  sleep(3);
1606 #endif
1607  ExecuteStreak();
1608 #ifndef _WIN32
1609  sleep(3);
1610 #endif
1611 
1612  if (m_moveMeshToCriticalLayer)
1613  {
1614  string syscall;
1615  char alpchar[16] = "";
1616  sprintf(alpchar, "%f", m_alpha[0]);
1617 
1618  string filePost = m_sessionName + "_advPost.xml";
1619  string filestreak = m_sessionName + "_streak.fld";
1620  string filewave = m_sessionName + "_wave.fld";
1621  string filewavepressure = m_sessionName + "_wave_p_split.fld";
1622  string fileinterp = m_sessionName + "_interp.xml";
1623  string interpstreak = m_sessionName + "_interpstreak.fld";
1624  string interwavepressure =
1625  m_sessionName + "_wave_p_split_interp.fld";
1626  syscall = "../../utilities/PostProcessing/Extras/MoveMesh " +
1627  filePost + " " + filestreak + " " + fileinterp + " " +
1628  alpchar;
1629 
1630  cout << syscall.c_str() << endl;
1631  if (system(syscall.c_str()))
1632  {
1633  ASSERTL0(false, syscall.c_str());
1634  }
1635 
1636  // move the advPost mesh (remark update alpha!!!)
1637  syscall = "../../utilities/PostProcessing/Extras/MoveMesh " +
1638  filePost + " " + filestreak + " " + filePost + " " +
1639  alpchar;
1640  cout << syscall.c_str() << endl;
1641  if (system(syscall.c_str()))
1642  {
1643  ASSERTL0(false, syscall.c_str());
1644  }
1645 
1646  // save oldstreak
1647  string oldstreak = m_sessionName + "_streak.fld";
1648  syscall = "cp -f " + filestreak + " " + oldstreak;
1649  cout << syscall.c_str() << endl;
1650  if (system(syscall.c_str()))
1651  {
1652  ASSERTL0(false, syscall.c_str());
1653  }
1654 
1655  // interpolate the streak field into the new mesh
1656  string movedmesh = m_sessionName + "_advPost_moved.xml";
1657  string movedinterpmesh = m_sessionName + "_interp_moved.xml";
1658 
1659  // create the interp streak
1660 
1661  syscall = "../../utilities/PostProcessing/Extras/FieldToField " +
1662  fileinterp + " " + filestreak + " " + movedinterpmesh +
1663  " " + interpstreak;
1664 
1665  cout << syscall.c_str() << endl;
1666  if (system(syscall.c_str()))
1667  {
1668  ASSERTL0(false, syscall.c_str());
1669  }
1670 
1671  // save the old mesh
1672  string meshfile = m_sessionName + ".xml";
1673  string meshold = m_sessionName + ".xml";
1674  syscall = "cp -f " + meshfile + " " + meshold;
1675  cout << syscall.c_str() << endl;
1676  if (system(syscall.c_str()))
1677  {
1678  ASSERTL0(false, syscall.c_str());
1679  }
1680 
1681  // overwriting the meshfile with the new mesh
1682  syscall = "cp -f " + movedmesh + " " + meshfile;
1683  cout << syscall.c_str() << endl;
1684  if (system(syscall.c_str()))
1685  {
1686  ASSERTL0(false, syscall.c_str());
1687  }
1688 
1689  // overwriting the streak file!!
1690  syscall = "cp -f " + interpstreak + " " + filestreak;
1691  cout << syscall.c_str() << endl;
1692  if (system(syscall.c_str()))
1693  {
1694  ASSERTL0(false, syscall.c_str());
1695  }
1696  }
1697 
1698  ExecuteWave();
1699 
1700  if (CalcWaveForce)
1701  {
1702  CalcNonLinearWaveForce();
1703  }
1704  }
1705 }
1706 
1707 bool VortexWaveInteraction::CheckEigIsStationary(bool reset)
1708 {
1709  static NekDouble previous_real_evl = -1.0;
1710  static NekDouble previous_imag_evl = -1.0;
1711  static int min_iter = 0;
1712 
1713  if (reset)
1714  {
1715  previous_real_evl = -1.0;
1716  min_iter = 0;
1717  }
1718 
1719  if (previous_real_evl == -1.0 || min_iter < m_minInnerIterations)
1720  {
1721  previous_real_evl = m_leading_real_evl[0];
1722  previous_imag_evl = m_leading_imag_evl[0];
1723  min_iter++;
1724  return false;
1725  }
1726 
1727  cout << "Growth tolerance: "
1728  << fabs((m_leading_real_evl[0] - previous_real_evl) /
1729  m_leading_real_evl[0])
1730  << endl;
1731  cout << "Phase tolerance: "
1732  << fabs((m_leading_imag_evl[0] - previous_imag_evl) /
1733  m_leading_imag_evl[0])
1734  << endl;
1735 
1736  // See if real and imaginary growth have converged to with m_eigRelTol
1737  if ((fabs((m_leading_real_evl[0] - previous_real_evl) /
1738  m_leading_real_evl[0]) < m_eigRelTol) ||
1739  (fabs(m_leading_real_evl[0]) < 1e-6))
1740  {
1741  previous_real_evl = m_leading_real_evl[0];
1742  previous_imag_evl = m_leading_imag_evl[0];
1743  if ((fabs((m_leading_imag_evl[0] - previous_imag_evl) /
1744  m_leading_imag_evl[0]) < m_eigRelTol) ||
1745  (fabs(m_leading_imag_evl[0]) < 1e-6))
1746  {
1747  return true;
1748  }
1749  }
1750  else
1751  {
1752  if (fabs(m_leading_imag_evl[0]) > 1e-2)
1753  {
1754  cout << "Warning: imaginary eigenvalue is greater than 1e-2"
1755  << endl;
1756  }
1757  previous_real_evl = m_leading_real_evl[0];
1758  previous_imag_evl = m_leading_imag_evl[0];
1759  return false;
1760  }
1761  return false;
1762 }
1763 
1764 // Check to see if leading eigenvalue is within tolerance defined
1765 // in m_neutralPointTol
1766 bool VortexWaveInteraction::CheckIfAtNeutralPoint(void)
1767 {
1768  bool returnval = false;
1769  if (m_sessionRoll->DefinesSolverInfo("INTERFACE"))
1770  {
1771  if (m_sessionVWI->GetSolverInfo("INTERFACE") == "phase")
1772  {
1773  if (abs(m_leading_real_evl[0]) < 1e-4)
1774  {
1775  returnval = true;
1776  }
1777  }
1778  else
1779  {
1780  if (abs(m_leading_real_evl[0]) < 1e-4 &&
1781  abs(m_leading_imag_evl[0]) < 2e-6)
1782  {
1783  returnval = true;
1784  }
1785  }
1786  }
1787  else
1788  {
1789  if ((m_leading_real_evl[0] * m_leading_real_evl[0] +
1790  m_leading_imag_evl[0] * m_leading_imag_evl[0]) <
1791  m_neutralPointTol * m_neutralPointTol)
1792  {
1793  returnval = true;
1794  }
1795  }
1796  if (fabs(m_leading_imag_evl[0]) > 1e-2)
1797  {
1798  cout << "Warning: imaginary eigenvalue is greater than 1e-2" << endl;
1799  }
1800 
1801  return returnval;
1802 }
1803 
1804 // Similar routine to UpdateAlpha
1805 
1806 void VortexWaveInteraction::UpdateWaveForceMag(int outeriter)
1807 {
1808  NekDouble wavef_new = 0.0;
1809 
1810  if (outeriter == 1)
1811  {
1812  m_waveForceMag[1] = m_waveForceMag[0];
1813  if (m_leading_real_evl[0] > 0.0)
1814  {
1815  wavef_new = m_waveForceMag[0] - m_waveForceMagStep;
1816  }
1817  else
1818  {
1819  wavef_new = m_waveForceMag[0] + m_waveForceMagStep;
1820  }
1821  }
1822  else
1823  {
1824  int i;
1825  int nstore = (m_waveForceMag.size() < outeriter) ? m_waveForceMag.size()
1826  : outeriter;
1827  Array<OneD, NekDouble> WaveForce(nstore);
1828  Array<OneD, NekDouble> Growth(nstore);
1829 
1830  Vmath::Vcopy(nstore, m_waveForceMag, 1, WaveForce, 1);
1831  Vmath::Vcopy(nstore, m_leading_real_evl, 1, Growth, 1);
1832 
1833  // Sort WaveForce Growth values;
1834  NekDouble store;
1835  int k;
1836  for (i = 0; i < nstore; ++i)
1837  {
1838  k = Vmath::Imin(nstore - i, &WaveForce[i], 1);
1839 
1840  store = WaveForce[i];
1841  WaveForce[i] = WaveForce[i + k];
1842  WaveForce[i + k] = store;
1843 
1844  store = Growth[i];
1845  Growth[i] = Growth[i + k];
1846  Growth[i + k] = store;
1847  }
1848 
1849  // See if we have any values that cross zero
1850  for (i = 0; i < nstore - 1; ++i)
1851  {
1852  if (Growth[i] * Growth[i + 1] < 0.0)
1853  {
1854  break;
1855  }
1856  }
1857 
1858  if (i != nstore - 1)
1859  {
1860  if (nstore == 2)
1861  {
1862  wavef_new =
1863  (WaveForce[0] * Growth[1] - WaveForce[1] * Growth[0]) /
1864  (Growth[1] - Growth[0]);
1865  }
1866  else
1867  {
1868  // use a quadratic fit and step through 10000 points
1869  // to find zero.
1870  int j;
1871  int nsteps = 10000;
1872  int idx = (i == 0) ? 1 : i;
1873  NekDouble da = WaveForce[idx + 1] - WaveForce[idx - 1];
1874  NekDouble gval_m1 = Growth[idx - 1], a, gval;
1875  NekDouble c1 = Growth[idx - 1] /
1876  (WaveForce[idx - 1] - WaveForce[idx]) /
1877  (WaveForce[idx - 1] - WaveForce[idx + 1]);
1878  NekDouble c2 = Growth[idx] /
1879  (WaveForce[idx] - WaveForce[idx - 1]) /
1880  (WaveForce[idx] - WaveForce[idx + 1]);
1881  NekDouble c3 = Growth[idx + 1] /
1882  (WaveForce[idx + 1] - WaveForce[idx - 1]) /
1883  (WaveForce[idx + 1] - WaveForce[idx]);
1884 
1885  for (j = 1; j < nsteps + 1; ++j)
1886  {
1887  a = WaveForce[i] + j * da / (NekDouble)nsteps;
1888  gval =
1889  c1 * (a - WaveForce[idx]) * (a - WaveForce[idx + 1]) +
1890  c2 * (a - WaveForce[idx - 1]) *
1891  (a - WaveForce[idx + 1]) +
1892  c3 * (a - WaveForce[idx - 1]) * (a - WaveForce[idx]);
1893 
1894  if (gval * gval_m1 < 0.0)
1895  {
1896  wavef_new = ((a + da / (NekDouble)nsteps) * gval -
1897  a * gval_m1) /
1898  (gval - gval_m1);
1899  break;
1900  }
1901  }
1902  }
1903  }
1904  else // step backward/forward
1905  {
1906  if (Growth[i] > 0.0)
1907  {
1908  wavef_new = m_waveForceMag[0] - m_waveForceMagStep;
1909  }
1910  else
1911  {
1912  wavef_new = m_waveForceMag[0] + m_waveForceMagStep;
1913  }
1914  }
1915  }
1916 
1917  for (int i = m_waveForceMag.size() - 1; i > 0; --i)
1918  {
1919  m_waveForceMag[i] = m_waveForceMag[i - 1];
1920  m_leading_real_evl[i] = m_leading_real_evl[i - 1];
1921  m_leading_imag_evl[i] = m_leading_imag_evl[i - 1];
1922  }
1923 
1924  m_waveForceMag[0] = wavef_new;
1925 }
1926 
1927 void VortexWaveInteraction::UpdateDAlphaDWaveForceMag(NekDouble alpha_init)
1928 {
1929  m_dAlphaDWaveForceMag = (m_alpha[0] - alpha_init) / m_waveForceMagStep;
1930 }
1931 
1932 void VortexWaveInteraction::UpdateAlpha(int outeriter)
1933 {
1934  NekDouble alp_new = 0.0;
1935 
1936  if (outeriter == 1)
1937  {
1938  m_alpha[1] = m_alpha[0];
1939  if (m_leading_real_evl[0] > 0.0)
1940  {
1941  alp_new = m_alpha[0] + m_alphaStep;
1942  }
1943  else
1944  {
1945  alp_new = m_alpha[0] - m_alphaStep;
1946  }
1947  }
1948  else
1949  {
1950  int i;
1951  int nstore = (m_alpha.size() < outeriter) ? m_alpha.size() : outeriter;
1952  Array<OneD, NekDouble> Alpha(nstore);
1953  Array<OneD, NekDouble> Growth(nstore);
1954 
1955  Vmath::Vcopy(nstore, m_alpha, 1, Alpha, 1);
1956  Vmath::Vcopy(nstore, m_leading_real_evl, 1, Growth, 1);
1957 
1958  // Sort Alpha Growth values;
1959  NekDouble store;
1960  int k;
1961  for (i = 0; i < nstore; ++i)
1962  {
1963  k = Vmath::Imin(nstore - i, &Alpha[i], 1);
1964 
1965  store = Alpha[i];
1966  Alpha[i] = Alpha[i + k];
1967  Alpha[i + k] = store;
1968 
1969  store = Growth[i];
1970  Growth[i] = Growth[i + k];
1971  Growth[i + k] = store;
1972  }
1973 
1974  // See if we have any values that cross zero
1975  for (i = 0; i < nstore - 1; ++i)
1976  {
1977  if (Growth[i] * Growth[i + 1] < 0.0)
1978  {
1979  break;
1980  }
1981  }
1982 
1983  if (i != nstore - 1)
1984  {
1985  if (nstore == 2)
1986  {
1987  alp_new = (Alpha[0] * Growth[1] - Alpha[1] * Growth[0]) /
1988  (Growth[1] - Growth[0]);
1989  }
1990  else
1991  {
1992  // use a quadratic fit and step through 10000 points
1993  // to find zero.
1994  int j;
1995  int nsteps = 10000;
1996  int idx = (i == 0) ? 1 : i;
1997  NekDouble da = Alpha[idx + 1] - Alpha[idx - 1];
1998  NekDouble gval_m1 = Growth[idx - 1], a, gval;
1999  NekDouble c1 = Growth[idx - 1] / (Alpha[idx - 1] - Alpha[idx]) /
2000  (Alpha[idx - 1] - Alpha[idx + 1]);
2001  NekDouble c2 = Growth[idx] / (Alpha[idx] - Alpha[idx - 1]) /
2002  (Alpha[idx] - Alpha[idx + 1]);
2003  NekDouble c3 = Growth[idx + 1] /
2004  (Alpha[idx + 1] - Alpha[idx - 1]) /
2005  (Alpha[idx + 1] - Alpha[idx]);
2006 
2007  for (j = 1; j < nsteps + 1; ++j)
2008  {
2009  a = Alpha[i] + j * da / (NekDouble)nsteps;
2010  gval = c1 * (a - Alpha[idx]) * (a - Alpha[idx + 1]) +
2011  c2 * (a - Alpha[idx - 1]) * (a - Alpha[idx + 1]) +
2012  c3 * (a - Alpha[idx - 1]) * (a - Alpha[idx]);
2013 
2014  if (gval * gval_m1 < 0.0)
2015  {
2016  alp_new = ((a + da / (NekDouble)nsteps) * gval -
2017  a * gval_m1) /
2018  (gval - gval_m1);
2019  break;
2020  }
2021  }
2022  }
2023  }
2024  else // step backward/forward
2025  {
2026  if (Growth[i] > 0.0)
2027  {
2028  alp_new = m_alpha[0] + m_alphaStep;
2029  }
2030  else
2031  {
2032  alp_new = m_alpha[0] - m_alphaStep;
2033  }
2034  }
2035  }
2036 
2037  for (int i = m_alpha.size() - 1; i > 0; --i)
2038  {
2039  m_alpha[i] = m_alpha[i - 1];
2040  m_leading_real_evl[i] = m_leading_real_evl[i - 1];
2041  m_leading_imag_evl[i] = m_leading_imag_evl[i - 1];
2042  }
2043 
2044  m_alpha[0] = alp_new;
2045 }
2046 
2047 Array<OneD, int> VortexWaveInteraction::GetReflectionIndex(void)
2048 {
2049  int i, j;
2050  int npts = m_waveVelocities[0]->GetPlane(0)->GetNpoints();
2051  int nel = m_waveVelocities[0]->GetNumElmts();
2052  Array<OneD, int> index(npts);
2053 
2054  Array<OneD, NekDouble> coord(2);
2055  Array<OneD, NekDouble> coord_x(npts);
2056  Array<OneD, NekDouble> coord_y(npts);
2057 
2058  //-> Dermine the point which is on coordinate (x -> -x + Lx/2, y-> -y)
2059  m_waveVelocities[0]->GetPlane(0)->GetCoords(coord_x, coord_y);
2060  NekDouble xmax = Vmath::Vmax(npts, coord_x, 1);
2061  // NekDouble tol =
2062  // NekConstants::kGeomFactorsTol*NekConstants::kGeomFactorsTol;
2063  NekDouble tol = 1e-5;
2064  NekDouble xnew, ynew;
2065 
2066  int start = npts - 1;
2067  int e_npts;
2068 
2069  bool useOnlyQuads = false;
2070  if (m_sessionVWI->DefinesSolverInfo("SymmetriseOnlyQuads"))
2071  {
2072  useOnlyQuads = true;
2073  }
2074 
2075  int cnt;
2076  for (int e = 0; e < nel; ++e)
2077  {
2078  e_npts = m_waveVelocities[0]->GetExp(e)->GetTotPoints();
2079  cnt = m_waveVelocities[0]->GetPhys_Offset(e);
2080 
2081  if (useOnlyQuads)
2082  {
2083  if (m_waveVelocities[0]->GetExp(e)->DetShapeType() ==
2084  LibUtilities::eTriangle)
2085  {
2086  for (i = 0; i < e_npts; ++i)
2087  {
2088  index[cnt + i] = -1;
2089  }
2090  continue;
2091  }
2092  }
2093 
2094  for (i = cnt; i < cnt + e_npts; ++i)
2095  {
2096  xnew = -coord_x[i] + xmax;
2097  ynew = -coord_y[i];
2098 
2099  for (j = start; j >= 0; --j)
2100  {
2101  if ((coord_x[j] - xnew) * (coord_x[j] - xnew) +
2102  (coord_y[j] - ynew) * (coord_y[j] - ynew) <
2103  tol)
2104  {
2105  index[i] = j;
2106  start = j;
2107  break;
2108  }
2109  }
2110 
2111  if (j == -1)
2112  {
2113 
2114  for (j = npts - 1; j > start; --j)
2115  {
2116 
2117  if ((coord_x[j] - xnew) * (coord_x[j] - xnew) +
2118  (coord_y[j] - ynew) * (coord_y[j] - ynew) <
2119  tol)
2120  {
2121  index[i] = j;
2122  break;
2123  }
2124  }
2125  ASSERTL0(j != start, "Failed to find matching point");
2126  }
2127  }
2128  }
2129  return index;
2130 }
2131 
2132 void VortexWaveInteraction::FileRelaxation(int reg)
2133 {
2134  cout << "relaxation..." << endl;
2135  static int cnt = 0;
2136  Array<OneD, MultiRegions::ExpListSharedPtr> Iexp =
2137  m_rollField[0]->GetBndCondExpansions();
2138  // cast to 1D explist (otherwise appenddata doesn't work)
2139  MultiRegions::ExpListSharedPtr Ilayer;
2140  Ilayer = MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(
2141  *std::static_pointer_cast<MultiRegions::ExpList>(Iexp[reg]));
2142  int nq = Ilayer->GetTotPoints();
2143  if (cnt == 0)
2144  {
2145  m_bcsForcing = Array<OneD, Array<OneD, NekDouble>>(4);
2146  m_bcsForcing[0] = Array<OneD, NekDouble>(4 * nq);
2147  for (int i = 1; i < 4; ++i)
2148  {
2149  m_bcsForcing[i] = m_bcsForcing[i - 1] + nq;
2150  }
2151  }
2152 
2153  // Read in mesh from input file
2154 
2155  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef_u;
2156  std::vector<std::vector<NekDouble>> FieldData_u;
2157  string file = m_sessionName;
2158 
2159  file += "_u_5.bc";
2160  LibUtilities::FieldIOSharedPtr fld =
2161  LibUtilities::FieldIO::CreateDefault(m_sessionVWI);
2162  fld->Import(file, FieldDef_u, FieldData_u);
2163  Ilayer->ExtractDataToCoeffs(FieldDef_u[0], FieldData_u[0],
2164  FieldDef_u[0]->m_fields[0],
2165  Ilayer->UpdateCoeffs());
2166  Ilayer->BwdTrans(Ilayer->GetCoeffs(), Ilayer->UpdatePhys());
2167 
2168  if (cnt == 0)
2169  {
2170  Vmath::Vcopy(nq, Ilayer->UpdatePhys(), 1, m_bcsForcing[2], 1);
2171  }
2172  Vmath::Vcopy(nq, Ilayer->UpdatePhys(), 1, m_bcsForcing[0], 1);
2173 
2174  // case relaxation==0 an additional array is needed
2175  Array<OneD, NekDouble> tmp_forcing(nq, 0.0);
2176 
2177  if (cnt != 0)
2178  {
2179  if (m_vwiRelaxation == 1.0)
2180  {
2181  Vmath::Vcopy(nq, m_bcsForcing[0], 1, tmp_forcing, 1);
2182  }
2183  Vmath::Smul(nq, 1.0 - m_vwiRelaxation, m_bcsForcing[0], 1,
2184  m_bcsForcing[0], 1);
2185  Vmath::Svtvp(nq, m_vwiRelaxation, m_bcsForcing[2], 1, m_bcsForcing[0],
2186  1, Ilayer->UpdatePhys(), 1);
2187  // generate again the bcs files:
2188 
2189  Array<OneD, Array<OneD, NekDouble>> fieldcoeffs(1);
2190  Ilayer->FwdTransLocalElmt(Ilayer->GetPhys(), Ilayer->UpdateCoeffs());
2191  fieldcoeffs[0] = Ilayer->UpdateCoeffs();
2192  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef1 =
2193  Ilayer->GetFieldDefinitions();
2194  std::vector<std::vector<NekDouble>> FieldData_1(FieldDef1.size());
2195  ;
2196  FieldDef1[0]->m_fields.push_back("u");
2197  Ilayer->AppendFieldData(FieldDef1[0], FieldData_1[0]);
2198  fld->Write(file, FieldDef1, FieldData_1);
2199  // save the bcs for the next iteration
2200  if (m_vwiRelaxation != 1.0)
2201  {
2202  Vmath::Smul(nq, 1. / (1.0 - m_vwiRelaxation), m_bcsForcing[0], 1,
2203  m_bcsForcing[0], 1);
2204  Vmath::Vcopy(nq, m_bcsForcing[0], 1, m_bcsForcing[2], 1);
2205  }
2206  else
2207  {
2208  Vmath::Vcopy(nq, tmp_forcing, 1, m_bcsForcing[2], 1);
2209  }
2210  }
2211 
2212  file = m_sessionName + "_v_5.bc";
2213 
2214  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef_v;
2215  std::vector<std::vector<NekDouble>> FieldData_v;
2216  fld->Import(file, FieldDef_v, FieldData_v);
2217  Ilayer->ExtractDataToCoeffs(FieldDef_v[0], FieldData_v[0],
2218  FieldDef_v[0]->m_fields[0],
2219  Ilayer->UpdateCoeffs());
2220  Ilayer->BwdTrans(Ilayer->GetCoeffs(), Ilayer->UpdatePhys());
2221  if (cnt == 0)
2222  {
2223  Vmath::Vcopy(nq, Ilayer->UpdatePhys(), 1, m_bcsForcing[3], 1);
2224  }
2225  Vmath::Vcopy(nq, Ilayer->UpdatePhys(), 1, m_bcsForcing[1], 1);
2226  if (cnt != 0)
2227  {
2228  if (m_vwiRelaxation == 1.0)
2229  {
2230  Vmath::Vcopy(nq, m_bcsForcing[1], 1, tmp_forcing, 1);
2231  }
2232  Vmath::Smul(nq, 1.0 - m_vwiRelaxation, m_bcsForcing[1], 1,
2233  m_bcsForcing[1], 1);
2234  Vmath::Svtvp(nq, m_vwiRelaxation, m_bcsForcing[3], 1, m_bcsForcing[1],
2235  1, Ilayer->UpdatePhys(), 1);
2236  // generate again the bcs files:
2237  Array<OneD, Array<OneD, NekDouble>> fieldcoeffs(1);
2238  Ilayer->FwdTransLocalElmt(Ilayer->GetPhys(), Ilayer->UpdateCoeffs());
2239  fieldcoeffs[0] = Ilayer->UpdateCoeffs();
2240  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef2 =
2241  Ilayer->GetFieldDefinitions();
2242  std::vector<std::vector<NekDouble>> FieldData_2(FieldDef2.size());
2243  ;
2244  FieldDef2[0]->m_fields.push_back("v");
2245  Ilayer->AppendFieldData(FieldDef2[0], FieldData_2[0]);
2246  fld->Write(file, FieldDef2, FieldData_2);
2247  // save the bcs for the next iteration
2248  if (m_vwiRelaxation != 1.0)
2249  {
2250  Vmath::Smul(nq, 1. / (1.0 - m_vwiRelaxation), m_bcsForcing[1], 1,
2251  m_bcsForcing[1], 1);
2252  Vmath::Vcopy(nq, m_bcsForcing[1], 1, m_bcsForcing[3], 1);
2253  }
2254  else
2255  {
2256  Vmath::Vcopy(nq, tmp_forcing, 1, m_bcsForcing[3], 1);
2257  }
2258  }
2259 
2260  cnt++;
2261 }
2262 } // namespace Nektar
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
Nektar::IncNavierStokes
This class is the base class for Navier Stokes problems.
Definition: IncNavierStokes.h:135
Nektar::IncNavierStokes::AddForcing
void AddForcing(const SolverUtils::ForcingSharedPtr &pForce)
Definition: IncNavierStokes.cpp:1132
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:198
Nektar::LibUtilities::NekFactory::CreateInstance
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:144
Nektar::LibUtilities::SessionReader::CreateInstance
static SessionReaderSharedPtr CreateInstance(int argc, char *argv[])
Creates an instance of the SessionReader class.
Definition: SessionReader.h:135
Nektar::MemoryManager::AllocateSharedPtr
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:168
Nektar::MultiRegions::GlobalLinSys
A global linear system.
Definition: GlobalLinSys.h:72
Nektar::SpatialDomains::MeshGraph::Read
static MeshGraphSharedPtr Read(const LibUtilities::SessionReaderSharedPtr pSession, LibUtilities::DomainRangeShPtr rng=LibUtilities::NullDomainRangeShPtr, bool fillGraph=true)
Definition: MeshGraph.cpp:110
Nektar::VortexWaveInteraction::m_streakField
Array< OneD, MultiRegions::ExpListSharedPtr > m_streakField
Definition: VortexWaveInteraction.h:254
Nektar::VortexWaveInteraction::m_vwiForcing
Array< OneD, Array< OneD, NekDouble > > m_vwiForcing
Definition: VortexWaveInteraction.h:249
Nektar::VortexWaveInteraction::CheckEigIsStationary
bool CheckEigIsStationary(bool reset=false)
Definition: VortexWaveInteraction.cpp:1707
Nektar::VortexWaveInteraction::m_sessionVWI
LibUtilities::SessionReaderSharedPtr m_sessionVWI
Definition: VortexWaveInteraction.h:259
Nektar::VortexWaveInteraction::m_sessionRoll
LibUtilities::SessionReaderSharedPtr m_sessionRoll
Definition: VortexWaveInteraction.h:261
Nektar::VortexWaveInteraction::m_alpha
Array< OneD, NekDouble > m_alpha
< Leading imaginary eigenvalue
Definition: VortexWaveInteraction.h:233
Nektar::VortexWaveInteraction::m_sessionStreak
LibUtilities::SessionReaderSharedPtr m_sessionStreak
Definition: VortexWaveInteraction.h:265
Nektar::VortexWaveInteraction::SaveFile
void SaveFile(std::string fileend, std::string dir, int n)
Definition: VortexWaveInteraction.cpp:1111
Nektar::VortexWaveInteraction::CopyFile
void CopyFile(std::string file1end, std::string file2end)
Definition: VortexWaveInteraction.cpp:1152
Nektar::VortexWaveInteraction::GetReflectionIndex
Array< OneD, int > GetReflectionIndex(void)
Definition: VortexWaveInteraction.cpp:2047
Nektar::VortexWaveInteraction::SaveLoopDetails
void SaveLoopDetails(std::string dir, int i)
Definition: VortexWaveInteraction.cpp:1182
Nektar::VortexWaveInteraction::m_vwiForcingObj
SolverUtils::ForcingProgrammaticSharedPtr m_vwiForcingObj
Definition: VortexWaveInteraction.h:250
Nektar::VortexWaveInteraction::m_waveVelocities
Array< OneD, MultiRegions::ExpListSharedPtr > m_waveVelocities
Definition: VortexWaveInteraction.h:246
Nektar::VortexWaveInteraction::ExecuteLoop
void ExecuteLoop(bool CalcWaveForce=true)
Definition: VortexWaveInteraction.cpp:1201
Nektar::VortexWaveInteraction::m_graphStreak
SpatialDomains::MeshGraphSharedPtr m_graphStreak
Definition: VortexWaveInteraction.h:266
Nektar::VortexWaveInteraction::m_bcsForcing
Array< OneD, Array< OneD, NekDouble > > m_bcsForcing
Definition: VortexWaveInteraction.h:252
Nektar::VortexWaveInteraction::m_graphRoll
SpatialDomains::MeshGraphSharedPtr m_graphRoll
Definition: VortexWaveInteraction.h:262
Nektar::VortexWaveInteraction::m_rollField
Array< OneD, MultiRegions::ExpListSharedPtr > m_rollField
Definition: VortexWaveInteraction.h:256
Nektar::VortexWaveInteraction::MoveFile
void MoveFile(std::string fileend, std::string dir, int n)
Definition: VortexWaveInteraction.cpp:1132
Nektar::VortexWaveInteraction::m_graphWave
SpatialDomains::MeshGraphSharedPtr m_graphWave
Definition: VortexWaveInteraction.h:268
Nektar::VortexWaveInteraction::m_waveForceMag
Array< OneD, NekDouble > m_waveForceMag
Definition: VortexWaveInteraction.h:224
Nektar::VortexWaveInteraction::AppendEvlToFile
void AppendEvlToFile(std::string file, int n)
Definition: VortexWaveInteraction.cpp:1164
Nektar::VortexWaveInteraction::m_sessionWave
LibUtilities::SessionReaderSharedPtr m_sessionWave
Definition: VortexWaveInteraction.h:267
Nektar::VortexWaveInteraction::UpdateDAlphaDWaveForceMag
void UpdateDAlphaDWaveForceMag(NekDouble alphainit)
Definition: VortexWaveInteraction.cpp:1927
Nektar::VortexWaveInteraction::m_wavePressure
MultiRegions::ExpListSharedPtr m_wavePressure
Definition: VortexWaveInteraction.h:247
Nektar::VortexWaveInteraction::m_leading_real_evl
Array< OneD, NekDouble > m_leading_real_evl
Definition: VortexWaveInteraction.h:229
Nektar::VortexWaveInteraction::m_leading_imag_evl
Array< OneD, NekDouble > m_leading_imag_evl
< Leading real eigenvalue
Definition: VortexWaveInteraction.h:231
Nektar::VortexWaveInteraction::GetVWIIterationType
VWIIterationType GetVWIIterationType(void)
Definition: VortexWaveInteraction.h:106
Nektar::VortexWaveInteraction::m_solverRoll
EquationSystemSharedPtr m_solverRoll
Definition: VortexWaveInteraction.h:263
CG_Iterations.savefile
string savefile
Definition: CG_Iterations.py:312
Nektar::LibUtilities::FieldIOSharedPtr
std::shared_ptr< FieldIO > FieldIOSharedPtr
Definition: FieldIO.h:301
Nektar::LibUtilities::SessionReaderSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: SessionReader.h:115
Nektar::MultiRegions::ExpListSharedPtr
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Definition: LocTraceToTraceMap.h:55
Nektar::SolverUtils::DriverSharedPtr
std::shared_ptr< Driver > DriverSharedPtr
A shared pointer to a Driver object.
Definition: Driver.h:51
Nektar::SolverUtils::GetDriverFactory
DriverFactory & GetDriverFactory()
Definition: Driver.cpp:64
Nektar::SolverUtils::EquationSystemSharedPtr
std::shared_ptr< EquationSystem > EquationSystemSharedPtr
A shared pointer to an EquationSystem object.
Definition: EquationSystem.h:63
Nektar::SolverUtils::GetEquationSystemFactory
EquationSystemFactory & GetEquationSystemFactory()
Definition: EquationSystem.cpp:85
Nektar::SolverUtils::GetForcingFactory
ForcingFactory & GetForcingFactory()
Declaration of the forcing factory singleton.
Definition: Forcing.cpp:44
Nektar
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
Nektar::IncNavierStokesSharedPtr
std::shared_ptr< IncNavierStokes > IncNavierStokesSharedPtr
Definition: IncNavierStokes.h:297
Nektar::VWIIterationTypeMap
const std::string VWIIterationTypeMap[]
Definition: VortexWaveInteraction.h:62
Nektar::eFixedWaveForcingWithSubIterationOnAlpha
@ eFixedWaveForcingWithSubIterationOnAlpha
Definition: VortexWaveInteraction.h:58
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Vmath::Vsqrt
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:534
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:209
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:622
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:518
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:574
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:359
Vmath::Smul
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:248
Vmath::Imin
int Imin(int n, const T *x, const int incx)
Return the index of the minimum element in x.
Definition: Vmath.cpp:1023
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
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:384
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:945
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1255
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:419
tinysimd::log
scalarT< T > log(scalarT< T > in)
Definition: scalar.hpp:300
tinysimd::abs
scalarT< T > abs(scalarT< T > in)
Definition: scalar.hpp:295
tinysimd::sqrt
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:291
Generated on Tue Aug 23 2022 07:27:15 for Nektar++ by doxygen 1.9.1