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