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  string syscall = "";
607  string c = std::to_string(cnt);
608  string c_alpha = std::to_string(m_alpha[0]);
609  if (m_sessionVWI->GetSolverInfo("INTERFACE") == "phase")
610  {
611  string filePost = m_sessionName + "_advPost.xml";
612  syscall = "../../utilities/PostProcessing/Extras/FldCalcBCs " +
613  filePost + " " +
614  "meshhalf_pos_Spen_stability_moved.fld "
615  "meshhalf_pos_Spen_advPost_moved.fld " +
616  c_alpha + " > data_alpha0";
617  cout << syscall.c_str() << endl;
618  if (system(syscall.c_str()))
619  {
620  ASSERTL0(false, syscall.c_str());
621  }
622 
623  syscall = "cp -f meshhalf_pos_Spen_stability_moved_u_5.bc " +
624  m_sessionName + "_u_5.bc";
625  cout << syscall.c_str() << endl;
626  if (system(syscall.c_str()))
627  {
628  ASSERTL0(false, syscall.c_str());
629  }
630  syscall = "cp -f meshhalf_pos_Spen_stability_moved_v_5.bc " +
631  m_sessionName + "_v_5.bc";
632  cout << syscall.c_str() << endl;
633  if (system(syscall.c_str()))
634  {
635  ASSERTL0(false, syscall.c_str());
636  }
637  }
638  else
639  {
640  syscall = "../../utilities/PostProcessing/Extras/FldCalcBCs " +
641  movedmesh + " " + wavefile + " " + filestreak + " " +
642  c_alpha + " > datasub_" + c;
643  cout << syscall.c_str() << endl;
644  if (system(syscall.c_str()))
645  {
646  ASSERTL0(false, syscall.c_str());
647  }
648  }
649 
650  // relaxation for different alpha values? does it make sense?
651 
652  // save the wave
653  string wave_subalp = m_sessionName + "_wave_subalp_" + c + ".fld";
654  syscall = "cp -f " + wavefile + " " + wave_subalp;
655  cout << syscall.c_str() << endl;
656  if (system(syscall.c_str()))
657  {
658  ASSERTL0(false, syscall.c_str());
659  }
660  // FileRelaxation(3);
661  cnt++;
662  }
663  else
664  {
665  int npts = m_waveVelocities[0]->GetPlane(0)->GetNpoints();
666  int ncoeffs = m_waveVelocities[0]->GetPlane(0)->GetNcoeffs();
667  Array<OneD, NekDouble> val(npts), der1(2 * npts);
668  Array<OneD, NekDouble> der2 = der1 + npts;
669  Array<OneD, NekDouble> streak;
670  static int projectfield = -1;
671 
672  if (m_deltaFcnApprox)
673  {
674  streak = Array<OneD, NekDouble>(npts);
675  m_streakField[0]->BwdTrans(m_streakField[0]->GetCoeffs(), streak);
676  }
677 
678  // Set project field to be first field that has a Neumann
679  // boundary since this not impose any condition on the vertical
680  // boundaries Othersise set to zero.
681  if (projectfield == -1)
682  {
683  Array<OneD, const SpatialDomains::BoundaryConditionShPtr> BndConds;
684 
685  for (int i = 0; i < m_waveVelocities.size(); ++i)
686  {
687  BndConds = m_waveVelocities[i]->GetBndConditions();
688  for (int j = 0; j < BndConds.size(); ++j)
689  {
690  if (BndConds[j]->GetBoundaryConditionType() ==
691  SpatialDomains::eNeumann)
692  {
693  projectfield = i;
694  break;
695  }
696  }
697  if (projectfield != -1)
698  {
699  break;
700  }
701  }
702  if (projectfield == -1)
703  {
704  cout << "using first field to project non-linear forcing which "
705  "imposes a Dirichlet condition"
706  << endl;
707  projectfield = 0;
708  }
709  }
710 
711  // Shift m_vwiForcing in case of relaxation
712  Vmath::Vcopy(ncoeffs, m_vwiForcing[0], 1, m_vwiForcing[2], 1);
713  Vmath::Vcopy(ncoeffs, m_vwiForcing[1], 1, m_vwiForcing[3], 1);
714 
715  // determine inverse of area normalised field.
716  m_wavePressure->GetPlane(0)->BwdTrans(
717  m_wavePressure->GetPlane(0)->GetCoeffs(),
718  m_wavePressure->GetPlane(0)->UpdatePhys());
719  m_wavePressure->GetPlane(1)->BwdTrans(
720  m_wavePressure->GetPlane(1)->GetCoeffs(),
721  m_wavePressure->GetPlane(1)->UpdatePhys());
722  NekDouble invnorm;
723 
724  if (m_useLinfPressureNorm)
725  {
726  Vmath::Vmul(npts, m_wavePressure->GetPlane(0)->UpdatePhys(), 1,
727  m_wavePressure->GetPlane(0)->UpdatePhys(), 1, der1, 1);
728  Vmath::Vvtvp(npts, m_wavePressure->GetPlane(1)->UpdatePhys(), 1,
729  m_wavePressure->GetPlane(1)->UpdatePhys(), 1, der1, 1,
730  der1, 1);
731  Vmath::Vsqrt(npts, der1, 1, der1, 1);
732 
733  NekDouble Linf = Vmath::Vmax(npts, der1, 1);
734 
735  invnorm = 1.0 / Linf;
736  }
737  else
738  {
739  // Determine normalisation of pressure so that |P|/A = 1
740  NekDouble l2;
741  l2 = m_wavePressure->GetPlane(0)->L2(
742  m_wavePressure->GetPlane(0)->GetPhys());
743  invnorm = l2 * l2;
744  l2 = m_wavePressure->GetPlane(1)->L2(
745  m_wavePressure->GetPlane(1)->GetPhys());
746  invnorm += l2 * l2;
747  Vmath::Fill(2 * npts, 1.0, der1, 1);
748  NekDouble area = m_waveVelocities[0]->GetPlane(0)->Integral(der1);
749  cout << "Area: " << area << endl;
750  invnorm = sqrt(area / invnorm);
751  }
752 
753  // Get hold of arrays.
754  m_waveVelocities[0]->GetPlane(0)->BwdTrans(
755  m_waveVelocities[0]->GetPlane(0)->GetCoeffs(),
756  m_waveVelocities[0]->GetPlane(0)->UpdatePhys());
757  Array<OneD, NekDouble> u_real =
758  m_waveVelocities[0]->GetPlane(0)->UpdatePhys();
759  Vmath::Smul(npts, invnorm, u_real, 1, u_real, 1);
760  m_waveVelocities[0]->GetPlane(1)->BwdTrans(
761  m_waveVelocities[0]->GetPlane(1)->GetCoeffs(),
762  m_waveVelocities[0]->GetPlane(1)->UpdatePhys());
763  Array<OneD, NekDouble> u_imag =
764  m_waveVelocities[0]->GetPlane(1)->UpdatePhys();
765  Vmath::Smul(npts, invnorm, u_imag, 1, u_imag, 1);
766  m_waveVelocities[1]->GetPlane(0)->BwdTrans(
767  m_waveVelocities[1]->GetPlane(0)->GetCoeffs(),
768  m_waveVelocities[1]->GetPlane(0)->UpdatePhys());
769  Array<OneD, NekDouble> v_real =
770  m_waveVelocities[1]->GetPlane(0)->UpdatePhys();
771  Vmath::Smul(npts, invnorm, v_real, 1, v_real, 1);
772  m_waveVelocities[1]->GetPlane(1)->BwdTrans(
773  m_waveVelocities[1]->GetPlane(1)->GetCoeffs(),
774  m_waveVelocities[1]->GetPlane(1)->UpdatePhys());
775  Array<OneD, NekDouble> v_imag =
776  m_waveVelocities[1]->GetPlane(1)->UpdatePhys();
777  Vmath::Smul(npts, invnorm, v_imag, 1, v_imag, 1);
778 
779  // normalise wave for output
780  Vmath::Smul(2 * ncoeffs, invnorm, m_waveVelocities[0]->UpdateCoeffs(),
781  1, m_waveVelocities[0]->UpdateCoeffs(), 1);
782  Vmath::Smul(2 * ncoeffs, invnorm, m_waveVelocities[1]->UpdateCoeffs(),
783  1, m_waveVelocities[1]->UpdateCoeffs(), 1);
784  Vmath::Smul(2 * ncoeffs, invnorm, m_waveVelocities[2]->UpdateCoeffs(),
785  1, m_waveVelocities[2]->UpdateCoeffs(), 1);
786 
787  // dump field
788  {
789  std::vector<std::string> variables(3);
790  std::vector<Array<OneD, NekDouble>> outfield(3);
791  variables[0] = "u_w";
792  variables[1] = "v_w";
793  variables[2] = "w_w";
794  outfield[0] = m_waveVelocities[0]->UpdateCoeffs();
795  outfield[1] = m_waveVelocities[1]->UpdateCoeffs();
796  outfield[2] = m_waveVelocities[2]->UpdateCoeffs();
797  std::string outname = m_sessionName + "_wave.fld";
798  m_solverRoll->WriteFld(outname, m_waveVelocities[0], outfield,
799  variables);
800  }
801 
802 #if 1
803  int ncoeffs_p = m_wavePressure->GetPlane(0)->GetNcoeffs();
804  Vmath::Smul(ncoeffs_p, invnorm,
805  m_wavePressure->GetPlane(0)->UpdateCoeffs(), 1,
806  m_wavePressure->GetPlane(0)->UpdateCoeffs(), 1);
807  Vmath::Smul(ncoeffs_p, invnorm,
808  m_wavePressure->GetPlane(1)->UpdateCoeffs(), 1,
809  m_wavePressure->GetPlane(1)->UpdateCoeffs(), 1);
810 #else
811  m_wavePressure->GetPlane(0)->BwdTrans(
812  m_wavePressure->GetPlane(0)->GetCoeffs(),
813  m_wavePressure->GetPlane(0)->UpdatePhys());
814  Vmath::Smul(npts, invnorm, m_wavePressure->GetPlane(0)->UpdatePhys(), 1,
815  m_wavePressure->GetPlane(0)->UpdatePhys(), 1);
816  m_wavePressure->GetPlane(0)->FwdTrans(
817  m_wavePressure->GetPlane(0)->UpdatePhys(),
818  m_wavePressure->GetPlane(0)->UpdateCoeffs());
819 
820  m_wavePressure->GetPlane(1)->BwdTrans(
821  m_wavePressure->GetPlane(1)->GetCoeffs(),
822  m_wavePressure->GetPlane(1)->UpdatePhys());
823  Vmath::Smul(npts, invnorm, m_wavePressure->GetPlane(1)->UpdatePhys(), 1,
824  m_wavePressure->GetPlane(1)->UpdatePhys(), 1);
825  m_wavePressure->GetPlane(1)->FwdTrans(
826  m_wavePressure->GetPlane(1)->UpdatePhys(),
827  m_wavePressure->GetPlane(1)->UpdateCoeffs());
828 #endif
829 
830  // Calculate non-linear terms for x and y directions
831  // d/dx(u u* + u* u)
832  Vmath::Vmul(npts, u_real, 1, u_real, 1, val, 1);
833  Vmath::Vvtvp(npts, u_imag, 1, u_imag, 1, val, 1, val, 1);
834  Vmath::Smul(npts, 2.0, val, 1, val, 1);
835  m_waveVelocities[0]->GetPlane(0)->PhysDeriv(0, val, der1);
836 
837  // d/dy(v u* + v* u)
838  Vmath::Vmul(npts, u_real, 1, v_real, 1, val, 1);
839  Vmath::Vvtvp(npts, u_imag, 1, v_imag, 1, val, 1, val, 1);
840  Vmath::Smul(npts, 2.0, val, 1, val, 1);
841  m_waveVelocities[0]->GetPlane(0)->PhysDeriv(1, val, der2);
842 
843  Vmath::Vadd(npts, der1, 1, der2, 1, der1, 1);
844 #if 1
845  m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(der1,
846  m_vwiForcing[0]);
847 #else
848  m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(
849  der1, m_vwiForcing[0]);
850 #endif
851  Vmath::Smul(ncoeffs, -m_waveForceMag[0], m_vwiForcing[0], 1,
852  m_vwiForcing[0], 1);
853 
854  // d/dx(u v* + u* v)
855  m_waveVelocities[0]->GetPlane(0)->PhysDeriv(0, val, der1);
856 
857  // d/dy(v v* + v* v)
858  Vmath::Vmul(npts, v_real, 1, v_real, 1, val, 1);
859  Vmath::Vvtvp(npts, v_imag, 1, v_imag, 1, val, 1, val, 1);
860  Vmath::Smul(npts, 2.0, val, 1, val, 1);
861  m_waveVelocities[0]->GetPlane(0)->PhysDeriv(1, val, der2);
862 
863  Vmath::Vadd(npts, der1, 1, der2, 1, der1, 1);
864 
865 #if 1
866  m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(der1,
867  m_vwiForcing[1]);
868 #else
869  m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(
870  der1, m_vwiForcing[1]);
871 #endif
872 
873  Vmath::Smul(ncoeffs, -m_waveForceMag[0], m_vwiForcing[1], 1,
874  m_vwiForcing[1], 1);
875 
876  // by default the symmetrization is on
877  bool symm = true;
878  m_sessionVWI->MatchSolverInfo("Symmetrization", "True", symm, true);
879 #if 0
880  if(symm== true )
881  {
882 
883  // Symmetrise forcing
884  //-> Get coordinates
885  Array<OneD, NekDouble> coord(2);
886  Array<OneD, NekDouble> coord_x(npts);
887  Array<OneD, NekDouble> coord_y(npts);
888 
889  //-> Impose symmetry (x -> -x + Lx/2, y-> -y) on coordinates
890  m_waveVelocities[0]->GetPlane(0)->GetCoords(coord_x,coord_y);
891  NekDouble xmax = Vmath::Vmax(npts,coord_x,1);
892  Vmath::Neg(npts,coord_x,1);
893  Vmath::Sadd(npts,xmax,coord_x,1,coord_x,1);
894  Vmath::Neg(npts,coord_y,1);
895 
896  int i, physoffset;
897 
898  //-> Obtain list of expansion element ids for each point.
899  Array<OneD, int> Eid(npts);
900  // This search may not be necessary every iteration
901  for(i = 0; i < npts; ++i)
902  {
903  coord[0] = coord_x[i];
904  coord[1] = coord_y[i];
905 
906  // Note this will not quite be symmetric.
907  Eid[i] = m_waveVelocities[0]->GetPlane(0)->GetExpIndex(coord,1e-6);
908  }
909 
910  // Interpolate field 0
911  m_waveVelocities[0]->GetPlane(0)->BwdTrans(m_vwiForcing[0],der1);
912  for(i = 0; i < npts; ++i)
913  {
914  physoffset = m_waveVelocities[0]->GetPlane(0)->GetPhys_Offset(Eid[i]);
915  coord[0] = coord_x[i];
916  coord[1] = coord_y[i];
917  der2 [i] = m_waveVelocities[0]->GetPlane(0)->GetExp(Eid[i])->PhysEvaluate(coord, der1 + physoffset);
918  }
919  //-> Average field 0
920  Vmath::Vsub(npts,der1,1,der2,1,der2,1);
921  Vmath::Smul(npts,0.5,der2,1,der2,1);
922 #if 1
923  m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(der2,m_vwiForcing[0]);
924 #else
925  m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(der2, m_vwiForcing[0]);
926 #endif
927 
928  //-> Interpoloate field 1
929  m_waveVelocities[0]->GetPlane(0)->BwdTrans(m_vwiForcing[1],der1);
930  for(i = 0; i < npts; ++i)
931  {
932  physoffset = m_waveVelocities[0]->GetPlane(0)->GetPhys_Offset(Eid[i]);
933  coord[0] = coord_x[i];
934  coord[1] = coord_y[i];
935  der2[i] = m_waveVelocities[0]->GetPlane(0)->GetExp(Eid[i])->PhysEvaluate(coord, der1 + physoffset);
936  }
937 
938  //-> Average field 1
939  Vmath::Vsub(npts,der1,1,der2,1,der2,1);
940  Vmath::Smul(npts,0.5,der2,1,der2,1);
941 #if 1
942  m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(der2,m_vwiForcing[1]);
943 #else
944  m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(der2, m_vwiForceing[1]);
945 #endif
946  }
947 #else
948  int i;
949  if (symm == true)
950  {
951  cout << "symmetrization is active" << endl;
952  static Array<OneD, int> index = GetReflectionIndex();
953 
954  m_waveVelocities[0]->GetPlane(0)->BwdTrans(m_vwiForcing[0], der1);
955 
956  for (i = 0; i < npts; ++i)
957  {
958  if (index[i] != -1)
959  {
960  val[i] = 0.5 * (der1[i] - der1[index[i]]);
961  }
962  }
963 #if 1
964  m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(
965  val, m_vwiForcing[0]);
966 #else
967  m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(
968  val, m_vwiForcing[0]);
969 #endif
970 
971  m_waveVelocities[0]->GetPlane(0)->BwdTrans(m_vwiForcing[1], der2);
972  for (i = 0; i < npts; ++i)
973  {
974  if (index[i] != -1)
975  {
976  val[i] = 0.5 * (der2[i] - der2[index[i]]);
977  }
978  }
979 #if 1
980  m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(
981  val, m_vwiForcing[1]);
982 #else
983  m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(
984  val, m_vwiForcing[1]);
985 #endif
986  }
987 
988  Vmath::Vmul(npts, der1, 1, der1, 1, val, 1);
989  Vmath::Vvtvp(npts, der2, 1, der2, 1, val, 1, val, 1);
990  Vmath::Vsqrt(npts, val, 1, val, 1);
991  cout << "F_Linf: " << Vmath::Vmax(npts, val, 1) << endl;
992 
993 #endif
994 
995  if (m_vwiRelaxation)
996  {
997  Vmath::Smul(ncoeffs, 1.0 - m_vwiRelaxation, m_vwiForcing[0], 1,
998  m_vwiForcing[0], 1);
999  Vmath::Svtvp(ncoeffs, m_vwiRelaxation, m_vwiForcing[2], 1,
1000  m_vwiForcing[0], 1, m_vwiForcing[0], 1);
1001 
1002  Vmath::Smul(ncoeffs, 1.0 - m_vwiRelaxation, m_vwiForcing[1], 1,
1003  m_vwiForcing[1], 1);
1004  Vmath::Svtvp(ncoeffs, m_vwiRelaxation, m_vwiForcing[3], 1,
1005  m_vwiForcing[1], 1, m_vwiForcing[1], 1);
1006  }
1007 
1008  if (m_deltaFcnApprox)
1009  {
1010  for (int j = 0; j < 2; ++j)
1011  {
1012 
1013  m_waveVelocities[projectfield]->GetPlane(0)->BwdTrans(
1014  m_vwiForcing[j], der1);
1015  for (int i = 0; i < npts; ++i)
1016  {
1017  der1[i] *= exp(-streak[i] * streak[i] / m_deltaFcnDecay);
1018  }
1019  m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(
1020  der1, m_vwiForcing[j]);
1021  }
1022  }
1023 
1024  // dump output
1025  std::vector<std::string> variables(4);
1026  std::vector<Array<OneD, NekDouble>> outfield(4);
1027  variables[0] = "u";
1028  variables[1] = "v";
1029  variables[2] = "pr";
1030  variables[3] = "pi";
1031  outfield[0] = m_vwiForcing[0];
1032  outfield[1] = m_vwiForcing[1];
1033  Array<OneD, NekDouble> soln(npts, 0.0);
1034  m_wavePressure->GetPlane(0)->BwdTrans(
1035  m_wavePressure->GetPlane(0)->GetCoeffs(),
1036  m_wavePressure->GetPlane(0)->UpdatePhys());
1037  outfield[2] = Array<OneD, NekDouble>(ncoeffs);
1038  m_waveVelocities[0]->GetPlane(0)->FwdTransLocalElmt(
1039  m_wavePressure->GetPlane(0)->GetPhys(), outfield[2]);
1040  m_wavePressure->GetPlane(1)->BwdTrans(
1041  m_wavePressure->GetPlane(1)->GetCoeffs(),
1042  m_wavePressure->GetPlane(1)->UpdatePhys());
1043 
1044  Vmath::Vmul(npts, m_wavePressure->GetPlane(0)->UpdatePhys(), 1,
1045  m_wavePressure->GetPlane(0)->UpdatePhys(), 1, val, 1);
1046  Vmath::Vvtvp(npts, m_wavePressure->GetPlane(1)->UpdatePhys(), 1,
1047  m_wavePressure->GetPlane(1)->UpdatePhys(), 1, val, 1, val,
1048  1);
1049  cout << "int P^2: " << m_wavePressure->GetPlane(0)->Integral(val)
1050  << endl;
1051  Vmath::Vsqrt(npts, val, 1, val, 1);
1052  cout << "PLinf: " << Vmath::Vmax(npts, val, 1) << endl;
1053 
1054  outfield[3] = Array<OneD, NekDouble>(ncoeffs);
1055  m_waveVelocities[1]->GetPlane(0)->FwdTransLocalElmt(
1056  m_wavePressure->GetPlane(1)->GetPhys(), outfield[3]);
1057 
1058  std::string outname = m_sessionName + ".vwi";
1059 
1060  m_solverRoll->WriteFld(outname, m_waveVelocities[0]->GetPlane(0),
1061  outfield, variables);
1062  }
1063 }
1064 
1065 void VortexWaveInteraction::CalcL2ToLinfPressure(void)
1066 {
1067 
1068  ExecuteWave();
1069 
1070  m_wavePressure->GetPlane(0)->BwdTrans(
1071  m_wavePressure->GetPlane(0)->GetCoeffs(),
1072  m_wavePressure->GetPlane(0)->UpdatePhys());
1073  m_wavePressure->GetPlane(1)->BwdTrans(
1074  m_wavePressure->GetPlane(1)->GetCoeffs(),
1075  m_wavePressure->GetPlane(1)->UpdatePhys());
1076 
1077  int npts = m_waveVelocities[0]->GetPlane(0)->GetNpoints();
1078  NekDouble Linf;
1079  Array<OneD, NekDouble> val(2 * npts, 0.0);
1080 
1081  Vmath::Vmul(npts, m_wavePressure->GetPlane(0)->UpdatePhys(), 1,
1082  m_wavePressure->GetPlane(0)->UpdatePhys(), 1, val, 1);
1083  Vmath::Vvtvp(npts, m_wavePressure->GetPlane(1)->UpdatePhys(), 1,
1084  m_wavePressure->GetPlane(1)->UpdatePhys(), 1, val, 1, val, 1);
1085  cout << "int P^2 " << m_wavePressure->GetPlane(0)->Integral(val) << endl;
1086  Vmath::Vsqrt(npts, val, 1, val, 1);
1087 
1088  Linf = Vmath::Vmax(npts, val, 1);
1089  cout << "Linf: " << Linf << endl;
1090 
1091  NekDouble l2, norm;
1092  l2 =
1093  m_wavePressure->GetPlane(0)->L2(m_wavePressure->GetPlane(0)->GetPhys());
1094  norm = l2 * l2;
1095  l2 =
1096  m_wavePressure->GetPlane(1)->L2(m_wavePressure->GetPlane(1)->GetPhys());
1097  norm += l2 * l2;
1098 
1099  Vmath::Fill(npts, 1.0, val, 1);
1100  NekDouble area = m_waveVelocities[0]->GetPlane(0)->Integral(val);
1101 
1102  l2 = sqrt(norm / area);
1103 
1104  cout << "L2: " << l2 << endl;
1105 
1106  cout << "Ratio Linf/L2: " << Linf / l2 << endl;
1107 }
1108 
1109 void VortexWaveInteraction::SaveFile(string file, string dir, int n)
1110 {
1111  static map<string, int> opendir;
1112 
1113  if (opendir.find(dir) == opendir.end())
1114  {
1115  // make directory and presume will fail if it already exists
1116  string mkdir = "mkdir " + dir;
1117  ASSERTL0(system(mkdir.c_str()) == 0,
1118  "Failed to make directory '" + dir + "'");
1119 
1120  opendir[dir] = 1;
1121  }
1122 
1123  string savefile =
1124  dir + "/" + file + "." + boost::lexical_cast<std::string>(n);
1125  string syscall = "cp -f " + file + " " + savefile;
1126 
1127  ASSERTL0(system(syscall.c_str()) == 0, syscall.c_str());
1128 }
1129 
1130 void VortexWaveInteraction::MoveFile(string file, string dir, int n)
1131 {
1132  static map<string, int> opendir;
1133 
1134  if (opendir.find(dir) == opendir.end())
1135  {
1136  // make directory and presume will fail if it already exists
1137  string mkdir = "mkdir " + dir;
1138  ASSERTL0(system(mkdir.c_str()) == 0,
1139  "Failed to make directory '" + dir + "'");
1140  opendir[dir] = 1;
1141  }
1142 
1143  string savefile =
1144  dir + "/" + file + "." + boost::lexical_cast<std::string>(n);
1145  string syscall = "mv -f " + file + " " + savefile;
1146 
1147  ASSERTL0(system(syscall.c_str()) == 0, syscall.c_str());
1148 }
1149 
1150 void VortexWaveInteraction::CopyFile(string file1end, string file2end)
1151 {
1152  string cpfile1 = m_sessionName + file1end;
1153  string cpfile2 = m_sessionName + file2end;
1154  string syscall = "cp -f " + cpfile1 + " " + cpfile2;
1155 
1156  if (system(syscall.c_str()))
1157  {
1158  ASSERTL0(false, syscall.c_str());
1159  }
1160 }
1161 
1162 void VortexWaveInteraction::AppendEvlToFile(string file, int n)
1163 {
1164  FILE *fp;
1165  fp = fopen(file.c_str(), "a");
1166  fprintf(fp, "%d: %lf %16.12le %16.12le\n", n, m_alpha[0],
1167  m_leading_real_evl[0], m_leading_imag_evl[0]);
1168  fclose(fp);
1169 }
1170 
1171 void VortexWaveInteraction::AppendEvlToFile(string file, NekDouble WaveForceMag)
1172 {
1173  FILE *fp;
1174  fp = fopen(file.c_str(), "a");
1175  fprintf(fp, "%lf %lf %16.12le %16.12le\n", WaveForceMag, m_alpha[0],
1176  m_leading_real_evl[0], m_leading_imag_evl[0]);
1177  fclose(fp);
1178 }
1179 
1180 void VortexWaveInteraction::SaveLoopDetails(std::string SaveDir, int i)
1181 
1182 {
1183  // Save NS restart file
1184  SaveFile(m_sessionName + ".rst", SaveDir, i + 1);
1185  // Save Streak Solution
1186  SaveFile(m_sessionName + "_streak.fld", SaveDir, i);
1187  // Save Wave solution output
1188  SaveFile(m_sessionName + ".evl", SaveDir, i);
1189  SaveFile(m_sessionName + "_eig_0", SaveDir, i);
1190  // Save new field file of eigenvalue
1191  SaveFile(m_sessionName + ".fld", SaveDir, i);
1192  if (!(m_sessionVWI->DefinesSolverInfo("INTERFACE")))
1193  {
1194  // Save new forcing file
1195  SaveFile(m_sessionName + ".vwi", SaveDir, i + 1);
1196  }
1197 }
1198 
1199 void VortexWaveInteraction::ExecuteLoop(bool CalcWaveForce)
1200 {
1201 
1202  // the global info has to be written in the
1203  // roll session file to use the interface loop
1204  if (m_sessionRoll->DefinesSolverInfo("INTERFACE"))
1205  {
1206  static int cnt = 0;
1207  bool skiprollstreak = false;
1208  if (cnt == 0 && m_sessionVWI->GetParameter("rollstreakfromit") == 1)
1209  {
1210  skiprollstreak = true;
1211  cout << "skip roll-streak at the first iteration" << endl;
1212  }
1213 
1214  if (skiprollstreak != true)
1215  {
1216 
1217  LibUtilities::NekManager<
1218  MultiRegions::GlobalLinSysKey,
1219  MultiRegions::GlobalLinSys>::EnableManagement("GlobalLinSys");
1220  ExecuteRoll();
1221  LibUtilities::NekManager<
1222  MultiRegions::GlobalLinSysKey,
1223  MultiRegions::GlobalLinSys>::DisableManagement("GlobalLinSys");
1224 #ifndef _WIN32
1225  sleep(3);
1226 #endif
1227  ExecuteStreak();
1228 #ifndef _WIN32
1229  sleep(3);
1230 #endif
1231  }
1232 
1233  string syscall;
1234  string movedmesh = m_sessionName + "_advPost_moved.xml";
1235  string movedinterpmesh = m_sessionName + "_interp_moved.xml";
1236  // rewrite the Rollsessionfile (we start from the waleffe forcing)
1237  // string meshbndjumps = m_sessionName +"_bndjumps.xml";
1238  // if(cnt==0)
1239  //{
1240  // take the conditions tag from meshbndjumps and copy into
1241  // the rolls session file
1242  //}
1243 
1244  string c = std::to_string(cnt);
1245 
1246  // save old roll solution
1247  string oldroll = m_sessionName + "_roll_" + c + ".fld";
1248  syscall = "cp -f " + m_sessionName + "-Base.fld" + " " + oldroll;
1249  cout << syscall.c_str() << endl;
1250  if (system(syscall.c_str()))
1251  {
1252  ASSERTL0(false, syscall.c_str());
1253  }
1254  // define file names
1255  string filePost = m_sessionName + "_advPost.xml";
1256  string filestreak = m_sessionName + "_streak.fld";
1257  string filewave = m_sessionName + "_wave.fld";
1258  string filewavepressure = m_sessionName + "_wave_p_split.fld";
1259  string fileinterp = m_sessionName + "_interp.xml";
1260  string interpstreak = m_sessionName + "_interpstreak_" + c + ".fld";
1261  string interwavepressure =
1262  m_sessionName + "_wave_p_split_interp_" + c + ".fld";
1263  string c_alpha = std::to_string(m_alpha[0]);
1264  cout << "alpha = " << m_alpha[0] << endl;
1265 
1266  if (m_sessionVWI->GetSolverInfo("INTERFACE") != "phase")
1267  {
1268  cout << "zerophase" << endl;
1269 
1270  syscall = "../../utilities/PostProcessing/Extras/MoveMesh " +
1271  filePost + " " + filestreak + " " + fileinterp + " " +
1272  c_alpha;
1273 
1274  cout << syscall.c_str() << endl;
1275  if (system(syscall.c_str()))
1276  {
1277  ASSERTL0(false, syscall.c_str());
1278  }
1279 
1280  // move the advPost mesh (remark update alpha!!!)
1281  syscall = "../../utilities/PostProcessing/Extras/MoveMesh " +
1282  filePost + " " + filestreak + " " + filePost + " " +
1283  c_alpha;
1284  cout << syscall.c_str() << endl;
1285  if (system(syscall.c_str()))
1286  {
1287  ASSERTL0(false, syscall.c_str());
1288  }
1289 
1290  // save oldstreak
1291  string oldstreak = m_sessionName + "_streak_" + c + ".fld";
1292  syscall = "cp -f " + filestreak + " " + oldstreak;
1293  cout << syscall.c_str() << endl;
1294  if (system(syscall.c_str()))
1295  {
1296  ASSERTL0(false, syscall.c_str());
1297  }
1298 
1299  // interpolate the streak field into the new mesh
1300  string movedmesh = m_sessionName + "_advPost_moved.xml";
1301  string movedinterpmesh = m_sessionName + "_interp_moved.xml";
1302 
1303  // create the interp streak
1304 
1305  syscall = "../../utilities/PostProcessing/Extras/FieldToField " +
1306  fileinterp + " " + filestreak + " " + movedinterpmesh +
1307  " " + interpstreak;
1308 
1309  cout << syscall.c_str() << endl;
1310  if (system(syscall.c_str()))
1311  {
1312  ASSERTL0(false, syscall.c_str());
1313  }
1314 
1315  // save the old mesh
1316  string meshfile = m_sessionName + ".xml";
1317  string meshold = m_sessionName + "_" + c + ".xml";
1318  syscall = "cp -f " + meshfile + " " + meshold;
1319  cout << syscall.c_str() << endl;
1320  if (system(syscall.c_str()))
1321  {
1322  ASSERTL0(false, syscall.c_str());
1323  }
1324 
1325  // overwriting the meshfile with the new mesh
1326  syscall = "cp -f " + movedmesh + " " + meshfile;
1327  cout << syscall.c_str() << endl;
1328  if (system(syscall.c_str()))
1329  {
1330  ASSERTL0(false, syscall.c_str());
1331  }
1332 
1333  // overwriting the streak file!!
1334  syscall = "cp -f " + interpstreak + " " + filestreak;
1335  cout << syscall.c_str() << endl;
1336  if (system(syscall.c_str()))
1337  {
1338  ASSERTL0(false, syscall.c_str());
1339  }
1340 
1341  // calculate the wave
1342  ExecuteWave();
1343 
1344  // save the wave field:
1345  string oldwave = m_sessionName + "_wave_" + c + ".fld";
1346  syscall = "cp -f " + m_sessionName + ".fld" + " " + oldwave;
1347  cout << syscall.c_str() << endl;
1348  if (system(syscall.c_str()))
1349  {
1350  ASSERTL0(false, syscall.c_str());
1351  }
1352 
1353  // save old jump conditions:
1354  string ujump = m_sessionName + "_u_5.bc";
1355  syscall = "cp -f " + ujump + " " + m_sessionName + "_u_5.bc_" + c;
1356  cout << syscall.c_str() << endl;
1357  if (system(syscall.c_str()))
1358  {
1359  ASSERTL0(false, syscall.c_str());
1360  }
1361 
1362  string vjump = m_sessionName + "_v_5.bc";
1363  syscall = "cp -f " + vjump + " " + m_sessionName + "_v_5.bc_" + c;
1364  cout << syscall.c_str() << endl;
1365  if (system(syscall.c_str()))
1366  {
1367  ASSERTL0(false, syscall.c_str());
1368  }
1369  cnt++;
1370  c = std::to_string(cnt);
1371 
1372  // use relaxation
1373  if (GetVWIIterationType() !=
1374  eFixedWaveForcingWithSubIterationOnAlpha)
1375  {
1376  // the critical layer should be the bnd region 3
1377  // int reg =3;
1378  // FileRelaxation(reg);
1379  }
1380  // calculate the jump conditions
1381  string wavefile = m_sessionName + ".fld";
1382  syscall = "../../utilities/PostProcessing/Extras/FldCalcBCs " +
1383  movedmesh + " " + wavefile + " " + interpstreak +
1384  "> data" + c;
1385  cout << syscall.c_str() << endl;
1386  if (system(syscall.c_str()))
1387  {
1388  ASSERTL0(false, syscall.c_str());
1389  }
1390 
1391  // move the new name_interp_moved.xml into name_interp.xml
1392  syscall = "cp -f " + movedinterpmesh + " " + fileinterp;
1393  cout << syscall.c_str() << endl;
1394  if (system(syscall.c_str()))
1395  {
1396  ASSERTL0(false, syscall.c_str());
1397  }
1398  // move the new name_advPost_moved.xml into name_advPost.xml
1399  syscall = "cp -f " + movedmesh + " " + filePost;
1400  cout << syscall.c_str() << endl;
1401  if (system(syscall.c_str()))
1402  {
1403  ASSERTL0(false, syscall.c_str());
1404  }
1405  }
1406  else if (m_sessionVWI->GetSolverInfo("INTERFACE") == "phase")
1407  {
1408  cout << "phase" << endl;
1409  // determine cr:
1410  NekDouble cr;
1411  string cr_str;
1412  stringstream st;
1413 
1414  // calculate the wave
1415  ExecuteWave();
1416 
1417  // save oldstreak
1418  string oldstreak = m_sessionName + "_streak_" + c + ".fld";
1419  syscall = "cp -f " + filestreak + " " + oldstreak;
1420  cout << syscall.c_str() << endl;
1421  if (system(syscall.c_str()))
1422  {
1423  ASSERTL0(false, syscall.c_str());
1424  }
1425 
1426  // save wave
1427  syscall = "cp -f " + m_sessionName + ".fld" + " " + filewave;
1428  cout << syscall.c_str() << endl;
1429  if (system(syscall.c_str()))
1430  {
1431  ASSERTL0(false, syscall.c_str());
1432  }
1433  // save the old mesh
1434  string meshfile = m_sessionName + ".xml";
1435  string meshold = m_sessionName + "_" + c + ".xml";
1436  syscall = "cp -f " + meshfile + " " + meshold;
1437  cout << syscall.c_str() << endl;
1438  if (system(syscall.c_str()))
1439  {
1440  ASSERTL0(false, syscall.c_str());
1441  }
1442 
1443  // save the oldwave field:
1444  string oldwave = m_sessionName + "_wave_" + c + ".fld";
1445  syscall = "cp -f " + m_sessionName + ".fld" + " " + oldwave;
1446  cout << syscall.c_str() << endl;
1447  if (system(syscall.c_str()))
1448  {
1449  ASSERTL0(false, syscall.c_str());
1450  }
1451 
1452  // save old jump conditions:
1453  string ujump = m_sessionName + "_u_5.bc";
1454  syscall = "cp -f " + ujump + " " + m_sessionName + "_u_5.bc_" + c;
1455  cout << syscall.c_str() << endl;
1456  if (system(syscall.c_str()))
1457  {
1458  ASSERTL0(false, syscall.c_str());
1459  }
1460 
1461  string vjump = m_sessionName + "_v_5.bc";
1462  syscall = "cp -f " + vjump + " " + m_sessionName + "_v_5.bc_" + c;
1463  cout << syscall.c_str() << endl;
1464  if (system(syscall.c_str()))
1465  {
1466  ASSERTL0(false, syscall.c_str());
1467  }
1468  cnt++;
1469 
1470  cr = m_leading_imag_evl[0] / m_alpha[0];
1471  st << cr;
1472  cr_str = st.str();
1473  cout << "cr=" << cr_str << endl;
1474  // NB -g or NOT!!!
1475  // move the mesh around the critical layer
1476  syscall = "../../utilities/PostProcessing/Extras/MoveMesh " +
1477  filePost + " " + filestreak + " " + fileinterp + " " +
1478  c_alpha + " " + cr_str;
1479 
1480  cout << syscall.c_str() << endl;
1481  if (system(syscall.c_str()))
1482  {
1483  ASSERTL0(false, syscall.c_str());
1484  }
1485  // NB -g or NOT!!!
1486  // move the advPost mesh (remark update alpha!!!)
1487  syscall = "../../utilities/PostProcessing/Extras/MoveMesh " +
1488  filePost + " " + filestreak + " " + filePost + " " +
1489  c_alpha + " " + cr_str;
1490  cout << syscall.c_str() << endl;
1491  if (system(syscall.c_str()))
1492  {
1493  ASSERTL0(false, syscall.c_str());
1494  }
1495 
1496  // interp streak into the new mesh
1497  syscall = "../../utilities/PostProcessing/Extras/FieldToField " +
1498  fileinterp + " " + filestreak + " " + movedinterpmesh +
1499  " " + interpstreak;
1500 
1501  cout << syscall.c_str() << endl;
1502  if (system(syscall.c_str()))
1503  {
1504  ASSERTL0(false, syscall.c_str());
1505  }
1506 
1507  // split wave sol
1508  syscall = "../../utilities/PostProcessing/Extras/SplitFld " +
1509  filePost + " " + filewave;
1510 
1511  cout << syscall.c_str() << endl;
1512  if (system(syscall.c_str()))
1513  {
1514  ASSERTL0(false, syscall.c_str());
1515  }
1516  // interp wave
1517  syscall = "../../utilities/PostProcessing/Extras/FieldToField " +
1518  filePost + " " + filewavepressure + " " + movedmesh +
1519  " " + interwavepressure;
1520 
1521  cout << syscall.c_str() << endl;
1522  if (system(syscall.c_str()))
1523  {
1524  ASSERTL0(false, syscall.c_str());
1525  }
1526 
1527  // use relaxation
1528  if (GetVWIIterationType() !=
1529  eFixedWaveForcingWithSubIterationOnAlpha)
1530  {
1531  // the critical layer should be the bnd region 3
1532  // int reg =3;
1533  // FileRelaxation(reg);
1534  }
1535 
1536  // cp wavepressure to m_sessionName.fld(to get
1537  // the right bcs names using FldCalcBCs)
1538  syscall =
1539  "cp -f " + interwavepressure + " " + m_sessionName + ".fld";
1540  cout << syscall.c_str() << endl;
1541  if (system(syscall.c_str()))
1542  {
1543  ASSERTL0(false, syscall.c_str());
1544  }
1545 
1546  // calculate the jump conditions
1547  // NB -g or NOT!!!
1548  syscall = "../../utilities/PostProcessing/Extras/FldCalcBCs " +
1549  movedmesh + " " + m_sessionName + ".fld" + " " +
1550  interpstreak + "> data" + c;
1551  cout << syscall.c_str() << endl;
1552  if (system(syscall.c_str()))
1553  {
1554  ASSERTL0(false, syscall.c_str());
1555  }
1556 
1557  // overwriting the meshfile with the new mesh
1558  syscall = "cp -f " + movedmesh + " " + meshfile;
1559  cout << syscall.c_str() << endl;
1560  if (system(syscall.c_str()))
1561  {
1562  ASSERTL0(false, syscall.c_str());
1563  }
1564 
1565  // overwriting the streak file!! (maybe is useless)
1566  syscall = "cp -f " + interpstreak + " " + filestreak;
1567  cout << syscall.c_str() << endl;
1568  if (system(syscall.c_str()))
1569  {
1570  ASSERTL0(false, syscall.c_str());
1571  }
1572  // move the new name_interp_moved.xml into name_interp.xml
1573  syscall = "cp -f " + movedinterpmesh + " " + fileinterp;
1574  cout << syscall.c_str() << endl;
1575  if (system(syscall.c_str()))
1576  {
1577  ASSERTL0(false, syscall.c_str());
1578  }
1579  // move the new name_advPost_moved.xml into name_advPost.xml
1580  syscall = "cp -f " + movedmesh + " " + filePost;
1581  cout << syscall.c_str() << endl;
1582  if (system(syscall.c_str()))
1583  {
1584  ASSERTL0(false, syscall.c_str());
1585  }
1586  }
1587  }
1588  else
1589  {
1590  LibUtilities::NekManager<
1591  MultiRegions::GlobalLinSysKey,
1592  MultiRegions::GlobalLinSys>::EnableManagement("GlobalLinSys");
1593  ExecuteRoll();
1594  LibUtilities::NekManager<
1595  MultiRegions::GlobalLinSysKey,
1596  MultiRegions::GlobalLinSys>::DisableManagement("GlobalLinSys");
1597 
1598 #ifndef _WIN32
1599  sleep(3);
1600 #endif
1601  ExecuteStreak();
1602 #ifndef _WIN32
1603  sleep(3);
1604 #endif
1605 
1606  if (m_moveMeshToCriticalLayer)
1607  {
1608  string syscall;
1609  char alpchar[16] = "";
1610  snprintf(alpchar, 16, "%f", m_alpha[0]);
1611 
1612  string filePost = m_sessionName + "_advPost.xml";
1613  string filestreak = m_sessionName + "_streak.fld";
1614  string filewave = m_sessionName + "_wave.fld";
1615  string filewavepressure = m_sessionName + "_wave_p_split.fld";
1616  string fileinterp = m_sessionName + "_interp.xml";
1617  string interpstreak = m_sessionName + "_interpstreak.fld";
1618  string interwavepressure =
1619  m_sessionName + "_wave_p_split_interp.fld";
1620  syscall = "../../utilities/PostProcessing/Extras/MoveMesh " +
1621  filePost + " " + filestreak + " " + fileinterp + " " +
1622  alpchar;
1623 
1624  cout << syscall.c_str() << endl;
1625  if (system(syscall.c_str()))
1626  {
1627  ASSERTL0(false, syscall.c_str());
1628  }
1629 
1630  // move the advPost mesh (remark update alpha!!!)
1631  syscall = "../../utilities/PostProcessing/Extras/MoveMesh " +
1632  filePost + " " + filestreak + " " + filePost + " " +
1633  alpchar;
1634  cout << syscall.c_str() << endl;
1635  if (system(syscall.c_str()))
1636  {
1637  ASSERTL0(false, syscall.c_str());
1638  }
1639 
1640  // save oldstreak
1641  string oldstreak = m_sessionName + "_streak.fld";
1642  syscall = "cp -f " + filestreak + " " + oldstreak;
1643  cout << syscall.c_str() << endl;
1644  if (system(syscall.c_str()))
1645  {
1646  ASSERTL0(false, syscall.c_str());
1647  }
1648 
1649  // interpolate the streak field into the new mesh
1650  string movedmesh = m_sessionName + "_advPost_moved.xml";
1651  string movedinterpmesh = m_sessionName + "_interp_moved.xml";
1652 
1653  // create the interp streak
1654 
1655  syscall = "../../utilities/PostProcessing/Extras/FieldToField " +
1656  fileinterp + " " + filestreak + " " + movedinterpmesh +
1657  " " + interpstreak;
1658 
1659  cout << syscall.c_str() << endl;
1660  if (system(syscall.c_str()))
1661  {
1662  ASSERTL0(false, syscall.c_str());
1663  }
1664 
1665  // save the old mesh
1666  string meshfile = m_sessionName + ".xml";
1667  string meshold = m_sessionName + ".xml";
1668  syscall = "cp -f " + meshfile + " " + meshold;
1669  cout << syscall.c_str() << endl;
1670  if (system(syscall.c_str()))
1671  {
1672  ASSERTL0(false, syscall.c_str());
1673  }
1674 
1675  // overwriting the meshfile with the new mesh
1676  syscall = "cp -f " + movedmesh + " " + meshfile;
1677  cout << syscall.c_str() << endl;
1678  if (system(syscall.c_str()))
1679  {
1680  ASSERTL0(false, syscall.c_str());
1681  }
1682 
1683  // overwriting the streak file!!
1684  syscall = "cp -f " + interpstreak + " " + filestreak;
1685  cout << syscall.c_str() << endl;
1686  if (system(syscall.c_str()))
1687  {
1688  ASSERTL0(false, syscall.c_str());
1689  }
1690  }
1691 
1692  ExecuteWave();
1693 
1694  if (CalcWaveForce)
1695  {
1696  CalcNonLinearWaveForce();
1697  }
1698  }
1699 }
1700 
1701 bool VortexWaveInteraction::CheckEigIsStationary(bool reset)
1702 {
1703  static NekDouble previous_real_evl = -1.0;
1704  static NekDouble previous_imag_evl = -1.0;
1705  static int min_iter = 0;
1706 
1707  if (reset)
1708  {
1709  previous_real_evl = -1.0;
1710  min_iter = 0;
1711  }
1712 
1713  if (previous_real_evl == -1.0 || min_iter < m_minInnerIterations)
1714  {
1715  previous_real_evl = m_leading_real_evl[0];
1716  previous_imag_evl = m_leading_imag_evl[0];
1717  min_iter++;
1718  return false;
1719  }
1720 
1721  cout << "Growth tolerance: "
1722  << fabs((m_leading_real_evl[0] - previous_real_evl) /
1723  m_leading_real_evl[0])
1724  << endl;
1725  cout << "Phase tolerance: "
1726  << fabs((m_leading_imag_evl[0] - previous_imag_evl) /
1727  m_leading_imag_evl[0])
1728  << endl;
1729 
1730  // See if real and imaginary growth have converged to with m_eigRelTol
1731  if ((fabs((m_leading_real_evl[0] - previous_real_evl) /
1732  m_leading_real_evl[0]) < m_eigRelTol) ||
1733  (fabs(m_leading_real_evl[0]) < 1e-6))
1734  {
1735  previous_real_evl = m_leading_real_evl[0];
1736  previous_imag_evl = m_leading_imag_evl[0];
1737  if ((fabs((m_leading_imag_evl[0] - previous_imag_evl) /
1738  m_leading_imag_evl[0]) < m_eigRelTol) ||
1739  (fabs(m_leading_imag_evl[0]) < 1e-6))
1740  {
1741  return true;
1742  }
1743  }
1744  else
1745  {
1746  if (fabs(m_leading_imag_evl[0]) > 1e-2)
1747  {
1748  cout << "Warning: imaginary eigenvalue is greater than 1e-2"
1749  << endl;
1750  }
1751  previous_real_evl = m_leading_real_evl[0];
1752  previous_imag_evl = m_leading_imag_evl[0];
1753  return false;
1754  }
1755  return false;
1756 }
1757 
1758 // Check to see if leading eigenvalue is within tolerance defined
1759 // in m_neutralPointTol
1760 bool VortexWaveInteraction::CheckIfAtNeutralPoint(void)
1761 {
1762  bool returnval = false;
1763  if (m_sessionRoll->DefinesSolverInfo("INTERFACE"))
1764  {
1765  if (m_sessionVWI->GetSolverInfo("INTERFACE") == "phase")
1766  {
1767  if (abs(m_leading_real_evl[0]) < 1e-4)
1768  {
1769  returnval = true;
1770  }
1771  }
1772  else
1773  {
1774  if (abs(m_leading_real_evl[0]) < 1e-4 &&
1775  abs(m_leading_imag_evl[0]) < 2e-6)
1776  {
1777  returnval = true;
1778  }
1779  }
1780  }
1781  else
1782  {
1783  if ((m_leading_real_evl[0] * m_leading_real_evl[0] +
1784  m_leading_imag_evl[0] * m_leading_imag_evl[0]) <
1785  m_neutralPointTol * m_neutralPointTol)
1786  {
1787  returnval = true;
1788  }
1789  }
1790  if (fabs(m_leading_imag_evl[0]) > 1e-2)
1791  {
1792  cout << "Warning: imaginary eigenvalue is greater than 1e-2" << endl;
1793  }
1794 
1795  return returnval;
1796 }
1797 
1798 // Similar routine to UpdateAlpha
1799 
1800 void VortexWaveInteraction::UpdateWaveForceMag(int outeriter)
1801 {
1802  NekDouble wavef_new = 0.0;
1803 
1804  if (outeriter == 1)
1805  {
1806  m_waveForceMag[1] = m_waveForceMag[0];
1807  if (m_leading_real_evl[0] > 0.0)
1808  {
1809  wavef_new = m_waveForceMag[0] - m_waveForceMagStep;
1810  }
1811  else
1812  {
1813  wavef_new = m_waveForceMag[0] + m_waveForceMagStep;
1814  }
1815  }
1816  else
1817  {
1818  int i;
1819  int nstore = (m_waveForceMag.size() < outeriter) ? m_waveForceMag.size()
1820  : outeriter;
1821  Array<OneD, NekDouble> WaveForce(nstore);
1822  Array<OneD, NekDouble> Growth(nstore);
1823 
1824  Vmath::Vcopy(nstore, m_waveForceMag, 1, WaveForce, 1);
1825  Vmath::Vcopy(nstore, m_leading_real_evl, 1, Growth, 1);
1826 
1827  // Sort WaveForce Growth values;
1828  NekDouble store;
1829  int k;
1830  for (i = 0; i < nstore; ++i)
1831  {
1832  k = Vmath::Imin(nstore - i, &WaveForce[i], 1);
1833 
1834  store = WaveForce[i];
1835  WaveForce[i] = WaveForce[i + k];
1836  WaveForce[i + k] = store;
1837 
1838  store = Growth[i];
1839  Growth[i] = Growth[i + k];
1840  Growth[i + k] = store;
1841  }
1842 
1843  // See if we have any values that cross zero
1844  for (i = 0; i < nstore - 1; ++i)
1845  {
1846  if (Growth[i] * Growth[i + 1] < 0.0)
1847  {
1848  break;
1849  }
1850  }
1851 
1852  if (i != nstore - 1)
1853  {
1854  if (nstore == 2)
1855  {
1856  wavef_new =
1857  (WaveForce[0] * Growth[1] - WaveForce[1] * Growth[0]) /
1858  (Growth[1] - Growth[0]);
1859  }
1860  else
1861  {
1862  // use a quadratic fit and step through 10000 points
1863  // to find zero.
1864  int j;
1865  int nsteps = 10000;
1866  int idx = (i == 0) ? 1 : i;
1867  NekDouble da = WaveForce[idx + 1] - WaveForce[idx - 1];
1868  NekDouble gval_m1 = Growth[idx - 1], a, gval;
1869  NekDouble c1 = Growth[idx - 1] /
1870  (WaveForce[idx - 1] - WaveForce[idx]) /
1871  (WaveForce[idx - 1] - WaveForce[idx + 1]);
1872  NekDouble c2 = Growth[idx] /
1873  (WaveForce[idx] - WaveForce[idx - 1]) /
1874  (WaveForce[idx] - WaveForce[idx + 1]);
1875  NekDouble c3 = Growth[idx + 1] /
1876  (WaveForce[idx + 1] - WaveForce[idx - 1]) /
1877  (WaveForce[idx + 1] - WaveForce[idx]);
1878 
1879  for (j = 1; j < nsteps + 1; ++j)
1880  {
1881  a = WaveForce[i] + j * da / (NekDouble)nsteps;
1882  gval =
1883  c1 * (a - WaveForce[idx]) * (a - WaveForce[idx + 1]) +
1884  c2 * (a - WaveForce[idx - 1]) *
1885  (a - WaveForce[idx + 1]) +
1886  c3 * (a - WaveForce[idx - 1]) * (a - WaveForce[idx]);
1887 
1888  if (gval * gval_m1 < 0.0)
1889  {
1890  wavef_new = ((a + da / (NekDouble)nsteps) * gval -
1891  a * gval_m1) /
1892  (gval - gval_m1);
1893  break;
1894  }
1895  }
1896  }
1897  }
1898  else // step backward/forward
1899  {
1900  if (Growth[i] > 0.0)
1901  {
1902  wavef_new = m_waveForceMag[0] - m_waveForceMagStep;
1903  }
1904  else
1905  {
1906  wavef_new = m_waveForceMag[0] + m_waveForceMagStep;
1907  }
1908  }
1909  }
1910 
1911  for (int i = m_waveForceMag.size() - 1; i > 0; --i)
1912  {
1913  m_waveForceMag[i] = m_waveForceMag[i - 1];
1914  m_leading_real_evl[i] = m_leading_real_evl[i - 1];
1915  m_leading_imag_evl[i] = m_leading_imag_evl[i - 1];
1916  }
1917 
1918  m_waveForceMag[0] = wavef_new;
1919 }
1920 
1921 void VortexWaveInteraction::UpdateDAlphaDWaveForceMag(NekDouble alpha_init)
1922 {
1923  m_dAlphaDWaveForceMag = (m_alpha[0] - alpha_init) / m_waveForceMagStep;
1924 }
1925 
1926 void VortexWaveInteraction::UpdateAlpha(int outeriter)
1927 {
1928  NekDouble alp_new = 0.0;
1929 
1930  if (outeriter == 1)
1931  {
1932  m_alpha[1] = m_alpha[0];
1933  if (m_leading_real_evl[0] > 0.0)
1934  {
1935  alp_new = m_alpha[0] + m_alphaStep;
1936  }
1937  else
1938  {
1939  alp_new = m_alpha[0] - m_alphaStep;
1940  }
1941  }
1942  else
1943  {
1944  int i;
1945  int nstore = (m_alpha.size() < outeriter) ? m_alpha.size() : outeriter;
1946  Array<OneD, NekDouble> Alpha(nstore);
1947  Array<OneD, NekDouble> Growth(nstore);
1948 
1949  Vmath::Vcopy(nstore, m_alpha, 1, Alpha, 1);
1950  Vmath::Vcopy(nstore, m_leading_real_evl, 1, Growth, 1);
1951 
1952  // Sort Alpha Growth values;
1953  NekDouble store;
1954  int k;
1955  for (i = 0; i < nstore; ++i)
1956  {
1957  k = Vmath::Imin(nstore - i, &Alpha[i], 1);
1958 
1959  store = Alpha[i];
1960  Alpha[i] = Alpha[i + k];
1961  Alpha[i + k] = store;
1962 
1963  store = Growth[i];
1964  Growth[i] = Growth[i + k];
1965  Growth[i + k] = store;
1966  }
1967 
1968  // See if we have any values that cross zero
1969  for (i = 0; i < nstore - 1; ++i)
1970  {
1971  if (Growth[i] * Growth[i + 1] < 0.0)
1972  {
1973  break;
1974  }
1975  }
1976 
1977  if (i != nstore - 1)
1978  {
1979  if (nstore == 2)
1980  {
1981  alp_new = (Alpha[0] * Growth[1] - Alpha[1] * Growth[0]) /
1982  (Growth[1] - Growth[0]);
1983  }
1984  else
1985  {
1986  // use a quadratic fit and step through 10000 points
1987  // to find zero.
1988  int j;
1989  int nsteps = 10000;
1990  int idx = (i == 0) ? 1 : i;
1991  NekDouble da = Alpha[idx + 1] - Alpha[idx - 1];
1992  NekDouble gval_m1 = Growth[idx - 1], a, gval;
1993  NekDouble c1 = Growth[idx - 1] / (Alpha[idx - 1] - Alpha[idx]) /
1994  (Alpha[idx - 1] - Alpha[idx + 1]);
1995  NekDouble c2 = Growth[idx] / (Alpha[idx] - Alpha[idx - 1]) /
1996  (Alpha[idx] - Alpha[idx + 1]);
1997  NekDouble c3 = Growth[idx + 1] /
1998  (Alpha[idx + 1] - Alpha[idx - 1]) /
1999  (Alpha[idx + 1] - Alpha[idx]);
2000 
2001  for (j = 1; j < nsteps + 1; ++j)
2002  {
2003  a = Alpha[i] + j * da / (NekDouble)nsteps;
2004  gval = c1 * (a - Alpha[idx]) * (a - Alpha[idx + 1]) +
2005  c2 * (a - Alpha[idx - 1]) * (a - Alpha[idx + 1]) +
2006  c3 * (a - Alpha[idx - 1]) * (a - Alpha[idx]);
2007 
2008  if (gval * gval_m1 < 0.0)
2009  {
2010  alp_new = ((a + da / (NekDouble)nsteps) * gval -
2011  a * gval_m1) /
2012  (gval - gval_m1);
2013  break;
2014  }
2015  }
2016  }
2017  }
2018  else // step backward/forward
2019  {
2020  if (Growth[i] > 0.0)
2021  {
2022  alp_new = m_alpha[0] + m_alphaStep;
2023  }
2024  else
2025  {
2026  alp_new = m_alpha[0] - m_alphaStep;
2027  }
2028  }
2029  }
2030 
2031  for (int i = m_alpha.size() - 1; i > 0; --i)
2032  {
2033  m_alpha[i] = m_alpha[i - 1];
2034  m_leading_real_evl[i] = m_leading_real_evl[i - 1];
2035  m_leading_imag_evl[i] = m_leading_imag_evl[i - 1];
2036  }
2037 
2038  m_alpha[0] = alp_new;
2039 }
2040 
2041 Array<OneD, int> VortexWaveInteraction::GetReflectionIndex(void)
2042 {
2043  int i, j;
2044  int npts = m_waveVelocities[0]->GetPlane(0)->GetNpoints();
2045  int nel = m_waveVelocities[0]->GetNumElmts();
2046  Array<OneD, int> index(npts);
2047 
2048  Array<OneD, NekDouble> coord(2);
2049  Array<OneD, NekDouble> coord_x(npts);
2050  Array<OneD, NekDouble> coord_y(npts);
2051 
2052  //-> Dermine the point which is on coordinate (x -> -x + Lx/2, y-> -y)
2053  m_waveVelocities[0]->GetPlane(0)->GetCoords(coord_x, coord_y);
2054  NekDouble xmax = Vmath::Vmax(npts, coord_x, 1);
2055  // NekDouble tol =
2056  // NekConstants::kGeomFactorsTol*NekConstants::kGeomFactorsTol;
2057  NekDouble tol = 1e-5;
2058  NekDouble xnew, ynew;
2059 
2060  int start = npts - 1;
2061  int e_npts;
2062 
2063  bool useOnlyQuads = false;
2064  if (m_sessionVWI->DefinesSolverInfo("SymmetriseOnlyQuads"))
2065  {
2066  useOnlyQuads = true;
2067  }
2068 
2069  int cnt;
2070  for (int e = 0; e < nel; ++e)
2071  {
2072  e_npts = m_waveVelocities[0]->GetExp(e)->GetTotPoints();
2073  cnt = m_waveVelocities[0]->GetPhys_Offset(e);
2074 
2075  if (useOnlyQuads)
2076  {
2077  if (m_waveVelocities[0]->GetExp(e)->DetShapeType() ==
2078  LibUtilities::eTriangle)
2079  {
2080  for (i = 0; i < e_npts; ++i)
2081  {
2082  index[cnt + i] = -1;
2083  }
2084  continue;
2085  }
2086  }
2087 
2088  for (i = cnt; i < cnt + e_npts; ++i)
2089  {
2090  xnew = -coord_x[i] + xmax;
2091  ynew = -coord_y[i];
2092 
2093  for (j = start; j >= 0; --j)
2094  {
2095  if ((coord_x[j] - xnew) * (coord_x[j] - xnew) +
2096  (coord_y[j] - ynew) * (coord_y[j] - ynew) <
2097  tol)
2098  {
2099  index[i] = j;
2100  start = j;
2101  break;
2102  }
2103  }
2104 
2105  if (j == -1)
2106  {
2107 
2108  for (j = npts - 1; j > start; --j)
2109  {
2110 
2111  if ((coord_x[j] - xnew) * (coord_x[j] - xnew) +
2112  (coord_y[j] - ynew) * (coord_y[j] - ynew) <
2113  tol)
2114  {
2115  index[i] = j;
2116  break;
2117  }
2118  }
2119  ASSERTL0(j != start, "Failed to find matching point");
2120  }
2121  }
2122  }
2123  return index;
2124 }
2125 
2126 void VortexWaveInteraction::FileRelaxation(int reg)
2127 {
2128  cout << "relaxation..." << endl;
2129  static int cnt = 0;
2130  Array<OneD, MultiRegions::ExpListSharedPtr> Iexp =
2131  m_rollField[0]->GetBndCondExpansions();
2132  // cast to 1D explist (otherwise appenddata doesn't work)
2133  MultiRegions::ExpListSharedPtr Ilayer;
2134  Ilayer = MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(
2135  *std::static_pointer_cast<MultiRegions::ExpList>(Iexp[reg]));
2136  int nq = Ilayer->GetTotPoints();
2137  if (cnt == 0)
2138  {
2139  m_bcsForcing = Array<OneD, Array<OneD, NekDouble>>(4);
2140  m_bcsForcing[0] = Array<OneD, NekDouble>(4 * nq);
2141  for (int i = 1; i < 4; ++i)
2142  {
2143  m_bcsForcing[i] = m_bcsForcing[i - 1] + nq;
2144  }
2145  }
2146 
2147  // Read in mesh from input file
2148 
2149  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef_u;
2150  std::vector<std::vector<NekDouble>> FieldData_u;
2151  string file = m_sessionName;
2152 
2153  file += "_u_5.bc";
2154  LibUtilities::FieldIOSharedPtr fld =
2155  LibUtilities::FieldIO::CreateDefault(m_sessionVWI);
2156  fld->Import(file, FieldDef_u, FieldData_u);
2157  Ilayer->ExtractDataToCoeffs(FieldDef_u[0], FieldData_u[0],
2158  FieldDef_u[0]->m_fields[0],
2159  Ilayer->UpdateCoeffs());
2160  Ilayer->BwdTrans(Ilayer->GetCoeffs(), Ilayer->UpdatePhys());
2161 
2162  if (cnt == 0)
2163  {
2164  Vmath::Vcopy(nq, Ilayer->UpdatePhys(), 1, m_bcsForcing[2], 1);
2165  }
2166  Vmath::Vcopy(nq, Ilayer->UpdatePhys(), 1, m_bcsForcing[0], 1);
2167 
2168  // case relaxation==0 an additional array is needed
2169  Array<OneD, NekDouble> tmp_forcing(nq, 0.0);
2170 
2171  if (cnt != 0)
2172  {
2173  if (m_vwiRelaxation == 1.0)
2174  {
2175  Vmath::Vcopy(nq, m_bcsForcing[0], 1, tmp_forcing, 1);
2176  }
2177  Vmath::Smul(nq, 1.0 - m_vwiRelaxation, m_bcsForcing[0], 1,
2178  m_bcsForcing[0], 1);
2179  Vmath::Svtvp(nq, m_vwiRelaxation, m_bcsForcing[2], 1, m_bcsForcing[0],
2180  1, Ilayer->UpdatePhys(), 1);
2181  // generate again the bcs files:
2182 
2183  Array<OneD, Array<OneD, NekDouble>> fieldcoeffs(1);
2184  Ilayer->FwdTransLocalElmt(Ilayer->GetPhys(), Ilayer->UpdateCoeffs());
2185  fieldcoeffs[0] = Ilayer->UpdateCoeffs();
2186  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef1 =
2187  Ilayer->GetFieldDefinitions();
2188  std::vector<std::vector<NekDouble>> FieldData_1(FieldDef1.size());
2189  ;
2190  FieldDef1[0]->m_fields.push_back("u");
2191  Ilayer->AppendFieldData(FieldDef1[0], FieldData_1[0]);
2192  fld->Write(file, FieldDef1, FieldData_1);
2193  // save the bcs for the next iteration
2194  if (m_vwiRelaxation != 1.0)
2195  {
2196  Vmath::Smul(nq, 1. / (1.0 - m_vwiRelaxation), m_bcsForcing[0], 1,
2197  m_bcsForcing[0], 1);
2198  Vmath::Vcopy(nq, m_bcsForcing[0], 1, m_bcsForcing[2], 1);
2199  }
2200  else
2201  {
2202  Vmath::Vcopy(nq, tmp_forcing, 1, m_bcsForcing[2], 1);
2203  }
2204  }
2205 
2206  file = m_sessionName + "_v_5.bc";
2207 
2208  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef_v;
2209  std::vector<std::vector<NekDouble>> FieldData_v;
2210  fld->Import(file, FieldDef_v, FieldData_v);
2211  Ilayer->ExtractDataToCoeffs(FieldDef_v[0], FieldData_v[0],
2212  FieldDef_v[0]->m_fields[0],
2213  Ilayer->UpdateCoeffs());
2214  Ilayer->BwdTrans(Ilayer->GetCoeffs(), Ilayer->UpdatePhys());
2215  if (cnt == 0)
2216  {
2217  Vmath::Vcopy(nq, Ilayer->UpdatePhys(), 1, m_bcsForcing[3], 1);
2218  }
2219  Vmath::Vcopy(nq, Ilayer->UpdatePhys(), 1, m_bcsForcing[1], 1);
2220  if (cnt != 0)
2221  {
2222  if (m_vwiRelaxation == 1.0)
2223  {
2224  Vmath::Vcopy(nq, m_bcsForcing[1], 1, tmp_forcing, 1);
2225  }
2226  Vmath::Smul(nq, 1.0 - m_vwiRelaxation, m_bcsForcing[1], 1,
2227  m_bcsForcing[1], 1);
2228  Vmath::Svtvp(nq, m_vwiRelaxation, m_bcsForcing[3], 1, m_bcsForcing[1],
2229  1, Ilayer->UpdatePhys(), 1);
2230  // generate again the bcs files:
2231  Array<OneD, Array<OneD, NekDouble>> fieldcoeffs(1);
2232  Ilayer->FwdTransLocalElmt(Ilayer->GetPhys(), Ilayer->UpdateCoeffs());
2233  fieldcoeffs[0] = Ilayer->UpdateCoeffs();
2234  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef2 =
2235  Ilayer->GetFieldDefinitions();
2236  std::vector<std::vector<NekDouble>> FieldData_2(FieldDef2.size());
2237  ;
2238  FieldDef2[0]->m_fields.push_back("v");
2239  Ilayer->AppendFieldData(FieldDef2[0], FieldData_2[0]);
2240  fld->Write(file, FieldDef2, FieldData_2);
2241  // save the bcs for the next iteration
2242  if (m_vwiRelaxation != 1.0)
2243  {
2244  Vmath::Smul(nq, 1. / (1.0 - m_vwiRelaxation), m_bcsForcing[1], 1,
2245  m_bcsForcing[1], 1);
2246  Vmath::Vcopy(nq, m_bcsForcing[1], 1, m_bcsForcing[3], 1);
2247  }
2248  else
2249  {
2250  Vmath::Vcopy(nq, tmp_forcing, 1, m_bcsForcing[3], 1);
2251  }
2252  }
2253 
2254  cnt++;
2255 }
2256 } // 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:1131
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:197
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:111
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:1701
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:1109
Nektar::VortexWaveInteraction::CopyFile
void CopyFile(std::string file1end, std::string file2end)
Definition: VortexWaveInteraction.cpp:1150
Nektar::VortexWaveInteraction::GetReflectionIndex
Array< OneD, int > GetReflectionIndex(void)
Definition: VortexWaveInteraction.cpp:2041
Nektar::VortexWaveInteraction::SaveLoopDetails
void SaveLoopDetails(std::string dir, int i)
Definition: VortexWaveInteraction.cpp:1180
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:1199
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:1130
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:1162
Nektar::VortexWaveInteraction::m_sessionWave
LibUtilities::SessionReaderSharedPtr m_sessionWave
Definition: VortexWaveInteraction.h:267
Nektar::VortexWaveInteraction::UpdateDAlphaDWaveForceMag
void UpdateDAlphaDWaveForceMag(NekDouble alphainit)
Definition: VortexWaveInteraction.cpp:1921
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:327
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:2
Nektar::IncNavierStokesSharedPtr
std::shared_ptr< IncNavierStokes > IncNavierStokesSharedPtr
Definition: IncNavierStokes.h:292
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 scalar 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:303
tinysimd::abs
scalarT< T > abs(scalarT< T > in)
Definition: scalar.hpp:298
tinysimd::sqrt
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:294
Generated on Sat Mar 25 2023 12:35:53 for Nektar++ by doxygen 1.9.1