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