Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | Private Attributes | List of all members
Nektar::VortexWaveInteraction Class Reference

#include <VortexWaveInteraction.h>

Collaboration diagram for Nektar::VortexWaveInteraction:
Collaboration graph
[legend]

Public Member Functions

 VortexWaveInteraction (int argc, char *argv[])
 
 ~VortexWaveInteraction (void)
 
void ExecuteRoll (void)
 
void ExecuteStreak (void)
 
void ExecuteWave (void)
 
void ExecuteLoop (bool CalcWaveForce=true)
 
void SaveLoopDetails (string dir, int i)
 
void CalcNonLinearWaveForce (void)
 
void CalcL2ToLinfPressure (void)
 
void SaveFile (string fileend, string dir, int n)
 
void MoveFile (string fileend, string dir, int n)
 
void CopyFile (string file1end, string file2end)
 
bool CheckEigIsStationary (bool reset=false)
 
bool CheckIfAtNeutralPoint (void)
 
void UpdateAlpha (int n)
 
void UpdateWaveForceMag (int n)
 
void UpdateDAlphaDWaveForceMag (NekDouble alphainit)
 
void AppendEvlToFile (std::string file, int n)
 
void AppendEvlToFile (std::string file, NekDouble WaveForceMag)
 
int GetIterStart ()
 
int GetIterEnd ()
 
VWIIterationType GetVWIIterationType (void)
 
int GetNOuterIterations (void)
 
int GetMaxOuterIterations (void)
 
NekDouble GetAlpha (void)
 
NekDouble GetAlphaStep (void)
 
NekDouble GetWaveForceMag (void)
 
NekDouble GetWaveForceMagStep (void)
 
NekDouble GetDAlphaDWaveForceMag (void)
 
int GetMaxWaveForceMagIter (void)
 
NekDouble GetEigRelTol (void)
 
int GetMinInnerIterations (void)
 
NekDouble GetPrevAlpha (void)
 
void SetAlpha (NekDouble alpha)
 
void SetWaveForceMag (NekDouble mag)
 
void SetEigRelTol (NekDouble tol)
 
void SetAlphaStep (NekDouble step)
 
void SetMinInnerIterations (int niter)
 
void SetPrevAlpha (NekDouble alpha)
 
bool IfIterInterface (void)
 
Array< OneD, int > GetReflectionIndex (void)
 
void FileRelaxation (int reg)
 

Private Attributes

int m_iterStart
 
int m_iterEnd
 
int m_nOuterIterations
 
int m_maxOuterIterations
 
int m_minInnerIterations
 
int m_maxWaveForceMagIter
 
bool m_deltaFcnApprox
 
bool m_useLinfPressureNorm
 
bool m_moveMeshToCriticalLayer
 
NekDouble m_deltaFcnDecay
 
Array< OneD, NekDoublem_waveForceMag
 
NekDouble m_waveForceMagStep
 
NekDouble m_rollForceScale
 
Array< OneD, NekDoublem_leading_real_evl
 
Array< OneD, NekDoublem_leading_imag_evl
 < Leading real eigenvalue More...
 
Array< OneD, NekDoublem_alpha
 < Leading imaginary eigenvalue More...
 
NekDouble m_alphaStep
 
NekDouble m_neutralPointTol
 
NekDouble m_eigRelTol
 
NekDouble m_vwiRelaxation
 
NekDouble m_dAlphaDWaveForceMag
 
NekDouble m_prevAlpha
 
bool m_iterinterface
 
VWIIterationType m_VWIIterationType
 
Array< OneD,
MultiRegions::ExpListSharedPtr
m_waveVelocities
 
MultiRegions::ExpListSharedPtr m_wavePressure
 
Array< OneD, Array< OneD,
NekDouble > > 
m_vwiForcing
 
SolverUtils::ForcingProgrammaticSharedPtr m_vwiForcingObj
 
Array< OneD, Array< OneD,
NekDouble > > 
m_bcsForcing
 
Array< OneD,
MultiRegions::ExpListSharedPtr
m_streakField
 
Array< OneD,
MultiRegions::ExpListSharedPtr
m_rollField
 
string m_sessionName
 
LibUtilities::SessionReaderSharedPtr m_sessionVWI
 
LibUtilities::SessionReaderSharedPtr m_sessionRoll
 
EquationSystemSharedPtr m_solverRoll
 
LibUtilities::SessionReaderSharedPtr m_sessionStreak
 
LibUtilities::SessionReaderSharedPtr m_sessionWave
 

Detailed Description

Definition at line 72 of file VortexWaveInteraction.h.

Constructor & Destructor Documentation

Nektar::VortexWaveInteraction::VortexWaveInteraction ( int  argc,
char *  argv[] 
)

Definition at line 45 of file VortexWaveInteraction.cpp.

References ASSERTL0, Nektar::LibUtilities::SessionReader::CreateInstance(), Nektar::LibUtilities::NekFactory< tKey, tBase, >::CreateInstance(), Nektar::eFixedAlpha, Nektar::eFixedAlphaWaveForcing, Nektar::eFixedWaveForcing, Nektar::eVWIIterationTypeSize, Nektar::SolverUtils::GetEquationSystemFactory(), m_alpha, m_alphaStep, m_dAlphaDWaveForceMag, m_deltaFcnApprox, m_deltaFcnDecay, m_eigRelTol, m_iterEnd, m_iterinterface, m_iterStart, m_leading_imag_evl, m_leading_real_evl, m_maxOuterIterations, m_maxWaveForceMagIter, m_minInnerIterations, m_moveMeshToCriticalLayer, m_neutralPointTol, m_nOuterIterations, m_rollForceScale, m_sessionName, m_sessionRoll, m_sessionStreak, m_sessionVWI, m_sessionWave, m_solverRoll, m_useLinfPressureNorm, m_vwiForcing, m_VWIIterationType, m_vwiRelaxation, m_waveForceMag, m_waveForceMagStep, UpdateAlpha(), and Nektar::VWIIterationTypeMap.

45  :
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 
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.
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  {
116  }
117  else
118  {
120  }
121 
122  m_alpha = Array<OneD, NekDouble> (storesize);
123  m_alpha[0] = m_sessionVWI->GetParameter("Alpha");
125  m_waveForceMag[0] = m_sessionVWI->GetParameter("WaveForceMag");
126 
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;
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");
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();
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  {
213  break;
214  }
215  }
216  }
217  else
218  {
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 
262  }
263  else
264  {
265  cout << " No File OuterIter.his to restart from" << endl;
266  }
267  }
268  break;
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  }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
LibUtilities::SessionReaderSharedPtr m_sessionRoll
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
Array< OneD, NekDouble > m_leading_imag_evl
< Leading real eigenvalue
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
EquationSystemSharedPtr m_solverRoll
Array< OneD, NekDouble > m_leading_real_evl
static SessionReaderSharedPtr CreateInstance(int argc, char *argv[])
Creates an instance of the SessionReader class.
LibUtilities::SessionReaderSharedPtr m_sessionWave
const std::string VWIIterationTypeMap[]
Array< OneD, NekDouble > m_waveForceMag
double NekDouble
EquationSystemFactory & GetEquationSystemFactory()
Array< OneD, Array< OneD, NekDouble > > m_vwiForcing
LibUtilities::SessionReaderSharedPtr m_sessionVWI
Array< OneD, NekDouble > m_alpha
< Leading imaginary eigenvalue
LibUtilities::SessionReaderSharedPtr m_sessionStreak
Nektar::VortexWaveInteraction::~VortexWaveInteraction ( void  )

Definition at line 346 of file VortexWaveInteraction.cpp.

References m_sessionVWI.

347  {
348  m_sessionVWI->Finalise();
349  }
LibUtilities::SessionReaderSharedPtr m_sessionVWI

Member Function Documentation

void Nektar::VortexWaveInteraction::AppendEvlToFile ( std::string  file,
int  n 
)

Definition at line 1060 of file VortexWaveInteraction.cpp.

References m_alpha, m_leading_imag_evl, and m_leading_real_evl.

Referenced by DoFixedForcingIteration(), and main().

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  }
Array< OneD, NekDouble > m_leading_imag_evl
< Leading real eigenvalue
Array< OneD, NekDouble > m_leading_real_evl
Array< OneD, NekDouble > m_alpha
< Leading imaginary eigenvalue
void Nektar::VortexWaveInteraction::AppendEvlToFile ( std::string  file,
NekDouble  WaveForceMag 
)

Definition at line 1068 of file VortexWaveInteraction.cpp.

References m_alpha, m_leading_imag_evl, and m_leading_real_evl.

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  }
Array< OneD, NekDouble > m_leading_imag_evl
< Leading real eigenvalue
Array< OneD, NekDouble > m_leading_real_evl
Array< OneD, NekDouble > m_alpha
< Leading imaginary eigenvalue
void Nektar::VortexWaveInteraction::CalcL2ToLinfPressure ( void  )

Definition at line 963 of file VortexWaveInteraction.cpp.

References ExecuteWave(), Vmath::Fill(), m_wavePressure, m_waveVelocities, npts, Vmath::Vmax(), Vmath::Vmul(), Vmath::Vsqrt(), and Vmath::Vvtvp().

Referenced by main().

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  }
MultiRegions::ExpListSharedPtr m_wavePressure
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 Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:46
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
Array< OneD, MultiRegions::ExpListSharedPtr > m_waveVelocities
static std::string npts
Definition: InputFld.cpp:43
double NekDouble
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
void Nektar::VortexWaveInteraction::CalcNonLinearWaveForce ( void  )

Definition at line 558 of file VortexWaveInteraction.cpp.

References ASSERTL0, Nektar::SpatialDomains::eNeumann, Vmath::Fill(), GetReflectionIndex(), m_alpha, m_deltaFcnApprox, m_deltaFcnDecay, m_sessionName, m_sessionRoll, m_sessionVWI, m_solverRoll, m_streakField, m_useLinfPressureNorm, m_vwiForcing, m_vwiRelaxation, m_waveForceMag, m_wavePressure, m_waveVelocities, Vmath::Neg(), npts, Vmath::Sadd(), Vmath::Smul(), Vmath::Svtvp(), Vmath::Vadd(), Vmath::Vcopy(), Vmath::Vmax(), Vmath::Vmul(), Vmath::Vsqrt(), Vmath::Vsub(), and Vmath::Vvtvp().

Referenced by DoFixedForcingIteration(), ExecuteLoop(), and main().

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  {
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 
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);
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);
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  }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
LibUtilities::SessionReaderSharedPtr m_sessionRoll
MultiRegions::ExpListSharedPtr m_wavePressure
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 Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:46
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
Array< OneD, MultiRegions::ExpListSharedPtr > m_waveVelocities
EquationSystemSharedPtr m_solverRoll
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
static std::string npts
Definition: InputFld.cpp:43
Array< OneD, NekDouble > m_waveForceMag
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:382
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
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
Array< OneD, MultiRegions::ExpListSharedPtr > m_streakField
Array< OneD, NekDouble > m_alpha
< Leading imaginary eigenvalue
Array< OneD, int > GetReflectionIndex(void)
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
bool Nektar::VortexWaveInteraction::CheckEigIsStationary ( bool  reset = false)

Definition at line 1603 of file VortexWaveInteraction.cpp.

References m_eigRelTol, m_leading_imag_evl, m_leading_real_evl, and m_minInnerIterations.

Referenced by DoFixedForcingIteration().

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  }
Array< OneD, NekDouble > m_leading_imag_evl
< Leading real eigenvalue
Array< OneD, NekDouble > m_leading_real_evl
double NekDouble
bool Nektar::VortexWaveInteraction::CheckIfAtNeutralPoint ( void  )

Definition at line 1651 of file VortexWaveInteraction.cpp.

References m_leading_imag_evl, m_leading_real_evl, m_neutralPointTol, m_sessionRoll, and m_sessionVWI.

Referenced by DoFixedForcingIteration().

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  {
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  }
LibUtilities::SessionReaderSharedPtr m_sessionRoll
Array< OneD, NekDouble > m_leading_imag_evl
< Leading real eigenvalue
Array< OneD, NekDouble > m_leading_real_evl
LibUtilities::SessionReaderSharedPtr m_sessionVWI
void Nektar::VortexWaveInteraction::CopyFile ( string  file1end,
string  file2end 
)

Definition at line 1048 of file VortexWaveInteraction.cpp.

References ASSERTL0, and m_sessionName.

Referenced by ExecuteRoll(), ExecuteStreak(), and ExecuteWave().

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  }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
void Nektar::VortexWaveInteraction::ExecuteLoop ( bool  CalcWaveForce = true)

Definition at line 1096 of file VortexWaveInteraction.cpp.

References ASSERTL0, CalcNonLinearWaveForce(), Nektar::LibUtilities::NekManager< KeyType, ValueT, opLessCreator >::DisableManagement(), Nektar::eFixedWaveForcingWithSubIterationOnAlpha, Nektar::LibUtilities::NekManager< KeyType, ValueT, opLessCreator >::EnableManagement(), ExecuteRoll(), ExecuteStreak(), ExecuteWave(), GetVWIIterationType(), m_alpha, m_leading_imag_evl, m_moveMeshToCriticalLayer, m_sessionName, m_sessionRoll, and m_sessionVWI.

Referenced by DoFixedForcingIteration(), and main().

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 
1117  ExecuteRoll();
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
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
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  {
1500  ExecuteRoll();
1502 
1503 #ifndef _WIN32
1504  sleep(3);
1505 #endif
1506  ExecuteStreak();
1507 #ifndef _WIN32
1508  sleep(3);
1509 #endif
1510 
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  {
1599  }
1600  }
1601  }
static void DisableManagement(std::string whichPool="")
Definition: NekManager.hpp:280
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
LibUtilities::SessionReaderSharedPtr m_sessionRoll
Array< OneD, NekDouble > m_leading_imag_evl
< Leading real eigenvalue
VWIIterationType GetVWIIterationType(void)
static void EnableManagement(std::string whichPool="")
Definition: NekManager.hpp:261
double NekDouble
LibUtilities::SessionReaderSharedPtr m_sessionVWI
Array< OneD, NekDouble > m_alpha
< Leading imaginary eigenvalue
void Nektar::VortexWaveInteraction::ExecuteRoll ( void  )

Definition at line 351 of file VortexWaveInteraction.cpp.

References Nektar::IncNavierStokes::AddForcing(), CopyFile(), Nektar::LibUtilities::NekFactory< tKey, tBase, >::CreateInstance(), Nektar::SolverUtils::GetEquationSystemFactory(), Nektar::SolverUtils::GetForcingFactory(), m_moveMeshToCriticalLayer, m_rollField, m_rollForceScale, m_sessionName, m_sessionRoll, m_solverRoll, m_vwiForcing, m_vwiForcingObj, Vmath::Smul(), and Vmath::Vcopy().

Referenced by ExecuteLoop().

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");
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  {
376  {
377  string vEquation = m_sessionRoll->GetSolverInfo("solvertype");
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 
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  }
LibUtilities::SessionReaderSharedPtr m_sessionRoll
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
boost::shared_ptr< IncNavierStokes > IncNavierStokesSharedPtr
Array< OneD, MultiRegions::ExpListSharedPtr > m_rollField
ForcingFactory & GetForcingFactory()
Declaration of the forcing factory singleton.
Definition: Forcing.cpp:42
EquationSystemSharedPtr m_solverRoll
SolverUtils::ForcingProgrammaticSharedPtr m_vwiForcingObj
boost::shared_ptr< EquationSystem > EquationSystemSharedPtr
A shared pointer to an EquationSystem object.
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
double NekDouble
EquationSystemFactory & GetEquationSystemFactory()
Array< OneD, Array< OneD, NekDouble > > m_vwiForcing
void AddForcing(const SolverUtils::ForcingSharedPtr &pForce)
void CopyFile(string file1end, string file2end)
This class is the base class for Navier Stokes problems.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
void Nektar::VortexWaveInteraction::ExecuteStreak ( void  )

Definition at line 457 of file VortexWaveInteraction.cpp.

References CopyFile(), Nektar::LibUtilities::NekFactory< tKey, tBase, >::CreateInstance(), Nektar::SolverUtils::GetDriverFactory(), Nektar::SolverUtils::GetEquationSystemFactory(), m_sessionStreak, m_sessionVWI, and m_streakField.

Referenced by ExecuteLoop().

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");
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  }
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
boost::shared_ptr< EquationSystem > EquationSystemSharedPtr
A shared pointer to an EquationSystem object.
double NekDouble
EquationSystemFactory & GetEquationSystemFactory()
LibUtilities::SessionReaderSharedPtr m_sessionVWI
void CopyFile(string file1end, string file2end)
DriverFactory & GetDriverFactory()
Definition: Driver.cpp:64
Array< OneD, MultiRegions::ExpListSharedPtr > m_streakField
LibUtilities::SessionReaderSharedPtr m_sessionStreak
void Nektar::VortexWaveInteraction::ExecuteWave ( void  )

Do linearised NavierStokes Session with Modified Arnoldi

Definition at line 496 of file VortexWaveInteraction.cpp.

References CopyFile(), Nektar::LibUtilities::NekFactory< tKey, tBase, >::CreateInstance(), Nektar::SolverUtils::GetDriverFactory(), m_alpha, m_leading_imag_evl, m_leading_real_evl, m_sessionVWI, m_sessionWave, m_wavePressure, and m_waveVelocities.

Referenced by CalcL2ToLinfPressure(), DoFixedForcingIteration(), ExecuteLoop(), and main().

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] +
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  }
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
Array< OneD, NekDouble > m_leading_imag_evl
< Leading real eigenvalue
Array< OneD, MultiRegions::ExpListSharedPtr > m_waveVelocities
Array< OneD, NekDouble > m_leading_real_evl
LibUtilities::SessionReaderSharedPtr m_sessionWave
double NekDouble
LibUtilities::SessionReaderSharedPtr m_sessionVWI
void CopyFile(string file1end, string file2end)
DriverFactory & GetDriverFactory()
Definition: Driver.cpp:64
Array< OneD, NekDouble > m_alpha
< Leading imaginary eigenvalue
void Nektar::VortexWaveInteraction::FileRelaxation ( int  reg)

====================================================

Todo:

Please update to use MeshGraph::Read(vSession)

Definition at line 2003 of file VortexWaveInteraction.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), m_bcsForcing, m_rollField, m_sessionName, m_sessionVWI, m_vwiRelaxation, Nektar::SpatialDomains::MeshGraph::Read(), Vmath::Smul(), Vmath::Svtvp(), and Vmath::Vcopy().

2004  {
2005  cout<<"relaxation..."<<endl;
2006  static int cnt=0;
2008  =m_rollField[0]->GetBndCondExpansions();
2009  //cast to 1D explist (otherwise appenddata doesn't work)
2013  *boost::static_pointer_cast<MultiRegions::ExpList1D>(Iexp[reg]));
2014  int nq = Ilayer->GetTotPoints();
2015  if( cnt==0)
2016  {
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  /// ====================================================
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";
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  }
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  }
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  }
static boost::shared_ptr< MeshGraph > Read(const LibUtilities::SessionReaderSharedPtr &pSession, DomainRangeShPtr &rng=NullDomainRangeShPtr)
Definition: MeshGraph.cpp:119
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
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
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
boost::shared_ptr< ExpList1D > ExpList1DSharedPtr
Shared pointer to an ExpList1D object.
Definition: ExpList1D.h:50
boost::shared_ptr< FieldIO > FieldIOSharedPtr
Definition: FieldIO.h:225
LibUtilities::SessionReaderSharedPtr m_sessionVWI
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
NekDouble Nektar::VortexWaveInteraction::GetAlpha ( void  )
inline

Definition at line 130 of file VortexWaveInteraction.h.

Referenced by DoFixedForcingIteration(), and main().

131  {
132  return m_alpha[0];
133  }
Array< OneD, NekDouble > m_alpha
< Leading imaginary eigenvalue
NekDouble Nektar::VortexWaveInteraction::GetAlphaStep ( void  )
inline

Definition at line 136 of file VortexWaveInteraction.h.

Referenced by main().

137  {
138  return m_alphaStep;
139  }
NekDouble Nektar::VortexWaveInteraction::GetDAlphaDWaveForceMag ( void  )
inline

Definition at line 151 of file VortexWaveInteraction.h.

152  {
153  return m_dAlphaDWaveForceMag;
154  }
NekDouble Nektar::VortexWaveInteraction::GetEigRelTol ( void  )
inline

Definition at line 162 of file VortexWaveInteraction.h.

Referenced by DoFixedForcingIteration().

164  {
165  return m_eigRelTol;
166  }
int Nektar::VortexWaveInteraction::GetIterEnd ( )
inline

Definition at line 110 of file VortexWaveInteraction.h.

Referenced by DoFixedForcingIteration().

111  {
112  return m_iterEnd;
113  }
int Nektar::VortexWaveInteraction::GetIterStart ( )
inline

Definition at line 105 of file VortexWaveInteraction.h.

Referenced by DoFixedForcingIteration().

106  {
107  return m_iterStart;
108  }
int Nektar::VortexWaveInteraction::GetMaxOuterIterations ( void  )
inline

Definition at line 125 of file VortexWaveInteraction.h.

Referenced by DoFixedForcingIteration().

126  {
127  return m_maxOuterIterations;
128  }
int Nektar::VortexWaveInteraction::GetMaxWaveForceMagIter ( void  )
inline

Definition at line 156 of file VortexWaveInteraction.h.

Referenced by main().

157  {
158  return m_maxWaveForceMagIter;
159  }
int Nektar::VortexWaveInteraction::GetMinInnerIterations ( void  )
inline

Definition at line 168 of file VortexWaveInteraction.h.

Referenced by DoFixedForcingIteration().

169  {
170  return m_minInnerIterations;
171  }
int Nektar::VortexWaveInteraction::GetNOuterIterations ( void  )
inline

Definition at line 120 of file VortexWaveInteraction.h.

Referenced by DoFixedForcingIteration().

121  {
122  return m_nOuterIterations;
123  }
NekDouble Nektar::VortexWaveInteraction::GetPrevAlpha ( void  )
inline

Definition at line 173 of file VortexWaveInteraction.h.

Referenced by DoFixedForcingIteration().

174  {
175  return m_prevAlpha;
176  }
Array< OneD, int > Nektar::VortexWaveInteraction::GetReflectionIndex ( void  )

Definition at line 1923 of file VortexWaveInteraction.cpp.

References ASSERTL0, Nektar::LibUtilities::eTriangle, m_sessionVWI, m_waveVelocities, npts, and Vmath::Vmax().

Referenced by CalcNonLinearWaveForce().

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  }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
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
Array< OneD, MultiRegions::ExpListSharedPtr > m_waveVelocities
static std::string npts
Definition: InputFld.cpp:43
double NekDouble
LibUtilities::SessionReaderSharedPtr m_sessionVWI
VWIIterationType Nektar::VortexWaveInteraction::GetVWIIterationType ( void  )
inline

Definition at line 115 of file VortexWaveInteraction.h.

Referenced by DoFixedForcingIteration(), ExecuteLoop(), and main().

116  {
117  return m_VWIIterationType;
118  }
NekDouble Nektar::VortexWaveInteraction::GetWaveForceMag ( void  )
inline

Definition at line 141 of file VortexWaveInteraction.h.

Referenced by main().

142  {
143  return m_waveForceMag[0];
144  }
Array< OneD, NekDouble > m_waveForceMag
NekDouble Nektar::VortexWaveInteraction::GetWaveForceMagStep ( void  )
inline

Definition at line 146 of file VortexWaveInteraction.h.

Referenced by main().

147  {
148  return m_waveForceMagStep;
149  }
bool Nektar::VortexWaveInteraction::IfIterInterface ( void  )
inline

Definition at line 209 of file VortexWaveInteraction.h.

Referenced by DoFixedForcingIteration().

210  {
211  return m_iterinterface;
212  }
void Nektar::VortexWaveInteraction::MoveFile ( string  fileend,
string  dir,
int  n 
)

Definition at line 1027 of file VortexWaveInteraction.cpp.

References ASSERTL0.

Referenced by DoFixedForcingIteration().

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  }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
void Nektar::VortexWaveInteraction::SaveFile ( string  fileend,
string  dir,
int  n 
)

Definition at line 1003 of file VortexWaveInteraction.cpp.

References ASSERTL0.

Referenced by SaveLoopDetails().

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  }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
void Nektar::VortexWaveInteraction::SaveLoopDetails ( std::string  SaveDir,
int  i 
)

Definition at line 1076 of file VortexWaveInteraction.cpp.

References m_sessionName, m_sessionVWI, and SaveFile().

Referenced by DoFixedForcingIteration().

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  }
LibUtilities::SessionReaderSharedPtr m_sessionVWI
void SaveFile(string fileend, string dir, int n)
void Nektar::VortexWaveInteraction::SetAlpha ( NekDouble  alpha)
inline

Definition at line 178 of file VortexWaveInteraction.h.

Referenced by main().

179  {
180  m_alpha[0] = alpha;
181  }
Array< OneD, NekDouble > m_alpha
< Leading imaginary eigenvalue
void Nektar::VortexWaveInteraction::SetAlphaStep ( NekDouble  step)
inline

Definition at line 194 of file VortexWaveInteraction.h.

195  {
196  m_alphaStep = step;
197  }
void Nektar::VortexWaveInteraction::SetEigRelTol ( NekDouble  tol)
inline

Definition at line 189 of file VortexWaveInteraction.h.

Referenced by DoFixedForcingIteration().

190  {
191  m_eigRelTol = tol;
192  }
void Nektar::VortexWaveInteraction::SetMinInnerIterations ( int  niter)
inline

Definition at line 199 of file VortexWaveInteraction.h.

Referenced by DoFixedForcingIteration().

200  {
201  m_minInnerIterations = niter;
202  }
void Nektar::VortexWaveInteraction::SetPrevAlpha ( NekDouble  alpha)
inline

Definition at line 204 of file VortexWaveInteraction.h.

Referenced by main().

205  {
206  m_prevAlpha = alpha;
207  }
void Nektar::VortexWaveInteraction::SetWaveForceMag ( NekDouble  mag)
inline

Definition at line 184 of file VortexWaveInteraction.h.

Referenced by main().

185  {
186  m_waveForceMag[0] = mag;
187  }
Array< OneD, NekDouble > m_waveForceMag
void Nektar::VortexWaveInteraction::UpdateAlpha ( int  n)

Definition at line 1810 of file VortexWaveInteraction.cpp.

References Vmath::Imin(), m_alpha, m_alphaStep, m_leading_imag_evl, m_leading_real_evl, and Vmath::Vcopy().

Referenced by DoFixedForcingIteration(), and VortexWaveInteraction().

1811  {
1812  NekDouble alp_new = 0.0;
1813 
1814 
1815  if(outeriter == 1)
1816  {
1817  m_alpha[1] = m_alpha[0];
1818  if(m_leading_real_evl[0] > 0.0)
1819  {
1820  alp_new = m_alpha[0] + m_alphaStep;
1821  }
1822  else
1823  {
1824  alp_new = m_alpha[0] - m_alphaStep;
1825  }
1826  }
1827  else
1828  {
1829  int i;
1830  int nstore = (m_alpha.num_elements() < outeriter)? m_alpha.num_elements(): outeriter;
1831  Array<OneD, NekDouble> Alpha(nstore);
1832  Array<OneD, NekDouble> Growth(nstore);
1833 
1834  Vmath::Vcopy(nstore,m_alpha,1,Alpha,1);
1835  Vmath::Vcopy(nstore,m_leading_real_evl,1,Growth,1);
1836 
1837  // Sort Alpha Growth values;
1838  NekDouble store;
1839  int k;
1840  for(i = 0; i < nstore; ++i)
1841  {
1842  k = Vmath::Imin(nstore-i,&Alpha[i],1);
1843 
1844  store = Alpha[i];
1845  Alpha[i] = Alpha[i+k];
1846  Alpha[i+k] = store;
1847 
1848  store = Growth[i];
1849  Growth[i] = Growth[i+k];
1850  Growth[i+k] = store;
1851  }
1852 
1853  // See if we have any values that cross zero
1854  for(i = 0; i < nstore-1; ++i)
1855  {
1856  if(Growth[i]*Growth[i+1] < 0.0)
1857  {
1858  break;
1859  }
1860  }
1861 
1862  if(i != nstore-1)
1863  {
1864  if(nstore == 2)
1865  {
1866  alp_new = (Alpha[0]*Growth[1] - Alpha[1]*Growth[0])/(Growth[1]-Growth[0]);
1867  }
1868  else
1869  {
1870  // use a quadratic fit and step through 10000 points
1871  // to find zero.
1872  int j;
1873  int nsteps = 10000;
1874  int idx = (i == 0)?1:i;
1875  NekDouble da = Alpha[idx+1] - Alpha[idx-1];
1876  NekDouble gval_m1 = Growth[idx-1],a,gval;
1877  NekDouble c1 = Growth[idx-1]/(Alpha[idx-1]-Alpha[idx])/
1878  (Alpha[idx-1]-Alpha[idx+1]);
1879  NekDouble c2 = Growth[idx]/(Alpha[idx]-Alpha[idx-1])/
1880  (Alpha[idx]-Alpha[idx+1]);
1881  NekDouble c3 = Growth[idx+1]/(Alpha[idx+1]-Alpha[idx-1])/
1882  (Alpha[idx+1]-Alpha[idx]);
1883 
1884  for(j = 1; j < nsteps+1; ++j)
1885  {
1886  a = Alpha[i] + j*da/(NekDouble) nsteps;
1887  gval = c1*(a - Alpha[idx ])*(a - Alpha[idx+1])
1888  + c2*(a - Alpha[idx-1])*(a - Alpha[idx+1])
1889  + c3*(a - Alpha[idx-1])*(a - Alpha[idx]);
1890 
1891  if(gval*gval_m1 < 0.0)
1892  {
1893  alp_new = ((a+da/(NekDouble)nsteps)*gval - a*gval_m1)/
1894  (gval - gval_m1);
1895  break;
1896  }
1897  }
1898  }
1899  }
1900  else // step backward/forward
1901  {
1902  if(Growth[i] > 0.0)
1903  {
1904  alp_new = m_alpha[0] + m_alphaStep;
1905  }
1906  else
1907  {
1908  alp_new = m_alpha[0] - m_alphaStep;
1909  }
1910  }
1911  }
1912 
1913  for(int i = m_alpha.num_elements()-1; i > 0; --i)
1914  {
1915  m_alpha[i] = m_alpha[i-1];
1918  }
1919 
1920  m_alpha[0] = alp_new;
1921  }
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
Array< OneD, NekDouble > m_leading_real_evl
double NekDouble
Array< OneD, NekDouble > m_alpha
< Leading imaginary eigenvalue
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
void Nektar::VortexWaveInteraction::UpdateDAlphaDWaveForceMag ( NekDouble  alphainit)

Definition at line 1805 of file VortexWaveInteraction.cpp.

References m_alpha, m_dAlphaDWaveForceMag, and m_waveForceMagStep.

Referenced by DoFixedForcingIteration().

1806  {
1808  }
Array< OneD, NekDouble > m_alpha
< Leading imaginary eigenvalue
void Nektar::VortexWaveInteraction::UpdateWaveForceMag ( int  n)

Definition at line 1691 of file VortexWaveInteraction.cpp.

References Vmath::Imin(), m_leading_imag_evl, m_leading_real_evl, m_waveForceMag, m_waveForceMagStep, and Vmath::Vcopy().

Referenced by DoFixedForcingIteration().

1692  {
1693  NekDouble wavef_new = 0.0;
1694 
1695 
1696  if(outeriter == 1)
1697  {
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];
1799  }
1800 
1801  m_waveForceMag[0] = wavef_new;
1802 
1803  }
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
Array< OneD, NekDouble > m_leading_real_evl
Array< OneD, NekDouble > m_waveForceMag
double NekDouble
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047

Member Data Documentation

Array<OneD, NekDouble> Nektar::VortexWaveInteraction::m_alpha
private
NekDouble Nektar::VortexWaveInteraction::m_alphaStep
private

Definition at line 246 of file VortexWaveInteraction.h.

Referenced by UpdateAlpha(), and VortexWaveInteraction().

Array<OneD, Array<OneD, NekDouble > > Nektar::VortexWaveInteraction::m_bcsForcing
private

Definition at line 263 of file VortexWaveInteraction.h.

Referenced by FileRelaxation().

NekDouble Nektar::VortexWaveInteraction::m_dAlphaDWaveForceMag
private

Definition at line 250 of file VortexWaveInteraction.h.

Referenced by UpdateDAlphaDWaveForceMag(), and VortexWaveInteraction().

bool Nektar::VortexWaveInteraction::m_deltaFcnApprox
private

Definition at line 229 of file VortexWaveInteraction.h.

Referenced by CalcNonLinearWaveForce(), and VortexWaveInteraction().

NekDouble Nektar::VortexWaveInteraction::m_deltaFcnDecay
private

Definition at line 234 of file VortexWaveInteraction.h.

Referenced by CalcNonLinearWaveForce(), and VortexWaveInteraction().

NekDouble Nektar::VortexWaveInteraction::m_eigRelTol
private

Definition at line 248 of file VortexWaveInteraction.h.

Referenced by CheckEigIsStationary(), and VortexWaveInteraction().

int Nektar::VortexWaveInteraction::m_iterEnd
private

Definition at line 222 of file VortexWaveInteraction.h.

Referenced by VortexWaveInteraction().

bool Nektar::VortexWaveInteraction::m_iterinterface
private

Definition at line 253 of file VortexWaveInteraction.h.

Referenced by VortexWaveInteraction().

int Nektar::VortexWaveInteraction::m_iterStart
private

Definition at line 221 of file VortexWaveInteraction.h.

Referenced by VortexWaveInteraction().

Array<OneD, NekDouble> Nektar::VortexWaveInteraction::m_leading_imag_evl
private
Array<OneD, NekDouble> Nektar::VortexWaveInteraction::m_leading_real_evl
private
int Nektar::VortexWaveInteraction::m_maxOuterIterations
private

Definition at line 225 of file VortexWaveInteraction.h.

Referenced by VortexWaveInteraction().

int Nektar::VortexWaveInteraction::m_maxWaveForceMagIter
private

Definition at line 227 of file VortexWaveInteraction.h.

Referenced by VortexWaveInteraction().

int Nektar::VortexWaveInteraction::m_minInnerIterations
private

Definition at line 226 of file VortexWaveInteraction.h.

Referenced by CheckEigIsStationary(), and VortexWaveInteraction().

bool Nektar::VortexWaveInteraction::m_moveMeshToCriticalLayer
private

Definition at line 232 of file VortexWaveInteraction.h.

Referenced by ExecuteLoop(), ExecuteRoll(), and VortexWaveInteraction().

NekDouble Nektar::VortexWaveInteraction::m_neutralPointTol
private

Definition at line 247 of file VortexWaveInteraction.h.

Referenced by CheckIfAtNeutralPoint(), and VortexWaveInteraction().

int Nektar::VortexWaveInteraction::m_nOuterIterations
private

Definition at line 224 of file VortexWaveInteraction.h.

Referenced by VortexWaveInteraction().

NekDouble Nektar::VortexWaveInteraction::m_prevAlpha
private

Definition at line 251 of file VortexWaveInteraction.h.

Array<OneD, MultiRegions::ExpListSharedPtr> Nektar::VortexWaveInteraction::m_rollField
private

Definition at line 267 of file VortexWaveInteraction.h.

Referenced by ExecuteRoll(), and FileRelaxation().

NekDouble Nektar::VortexWaveInteraction::m_rollForceScale
private

Definition at line 239 of file VortexWaveInteraction.h.

Referenced by ExecuteRoll(), and VortexWaveInteraction().

string Nektar::VortexWaveInteraction::m_sessionName
private
LibUtilities::SessionReaderSharedPtr Nektar::VortexWaveInteraction::m_sessionRoll
private
LibUtilities::SessionReaderSharedPtr Nektar::VortexWaveInteraction::m_sessionStreak
private

Definition at line 275 of file VortexWaveInteraction.h.

Referenced by ExecuteStreak(), and VortexWaveInteraction().

LibUtilities::SessionReaderSharedPtr Nektar::VortexWaveInteraction::m_sessionVWI
private
LibUtilities::SessionReaderSharedPtr Nektar::VortexWaveInteraction::m_sessionWave
private

Definition at line 276 of file VortexWaveInteraction.h.

Referenced by ExecuteWave(), and VortexWaveInteraction().

EquationSystemSharedPtr Nektar::VortexWaveInteraction::m_solverRoll
private
Array<OneD, MultiRegions::ExpListSharedPtr> Nektar::VortexWaveInteraction::m_streakField
private

Definition at line 265 of file VortexWaveInteraction.h.

Referenced by CalcNonLinearWaveForce(), and ExecuteStreak().

bool Nektar::VortexWaveInteraction::m_useLinfPressureNorm
private

Definition at line 230 of file VortexWaveInteraction.h.

Referenced by CalcNonLinearWaveForce(), and VortexWaveInteraction().

Array<OneD, Array<OneD, NekDouble > > Nektar::VortexWaveInteraction::m_vwiForcing
private
SolverUtils::ForcingProgrammaticSharedPtr Nektar::VortexWaveInteraction::m_vwiForcingObj
private

Definition at line 261 of file VortexWaveInteraction.h.

Referenced by ExecuteRoll().

VWIIterationType Nektar::VortexWaveInteraction::m_VWIIterationType
private

Definition at line 255 of file VortexWaveInteraction.h.

Referenced by VortexWaveInteraction().

NekDouble Nektar::VortexWaveInteraction::m_vwiRelaxation
private
Array<OneD, NekDouble> Nektar::VortexWaveInteraction::m_waveForceMag
private
NekDouble Nektar::VortexWaveInteraction::m_waveForceMagStep
private
MultiRegions::ExpListSharedPtr Nektar::VortexWaveInteraction::m_wavePressure
private
Array<OneD, MultiRegions::ExpListSharedPtr> Nektar::VortexWaveInteraction::m_waveVelocities
private