Nektar++
Public Member Functions | Private Attributes | List of all members
Nektar::VortexWaveInteraction Class Reference

#include <VortexWaveInteraction.h>

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 (std::string dir, int i)
 
void CalcNonLinearWaveForce (void)
 
void CalcL2ToLinfPressure (void)
 
void SaveFile (std::string fileend, std::string dir, int n)
 
void MoveFile (std::string fileend, std::string dir, int n)
 
void CopyFile (std::string file1end, std::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::ExpListSharedPtrm_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::ExpListSharedPtrm_streakField
 
Array< OneD, MultiRegions::ExpListSharedPtrm_rollField
 
std::string m_sessionName
 
LibUtilities::SessionReaderSharedPtr m_sessionVWI
 
LibUtilities::SessionReaderSharedPtr m_sessionRoll
 
SpatialDomains::MeshGraphSharedPtr m_graphRoll
 
EquationSystemSharedPtr m_solverRoll
 
LibUtilities::SessionReaderSharedPtr m_sessionStreak
 
SpatialDomains::MeshGraphSharedPtr m_graphStreak
 
LibUtilities::SessionReaderSharedPtr m_sessionWave
 
SpatialDomains::MeshGraphSharedPtr m_graphWave
 

Detailed Description

Definition at line 66 of file VortexWaveInteraction.h.

Constructor & Destructor Documentation

◆ VortexWaveInteraction()

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

Definition at line 46 of file VortexWaveInteraction.cpp.

48 {
49 
50  int storesize; // number of previous iterations to store
51 
52  m_sessionName = argv[argc - 1];
53  string meshfile = m_sessionName + ".xml";
54 
56 
57  std::string VWICondFile(argv[argc - 1]);
58  VWICondFile += "_VWI.xml";
59  std::vector<std::string> VWIFilenames;
60  VWIFilenames.push_back(meshfile);
61  VWIFilenames.push_back(VWICondFile);
62 
63  // Create Incompressible NavierStokesSolver session reader.
64  m_sessionVWI =
65  LibUtilities::SessionReader::CreateInstance(argc, argv, VWIFilenames);
66  m_sessionVWI->InitSession();
67 
68  m_sessionVWI->LoadParameter("AlphaStep", m_alphaStep, 0.05);
69  m_sessionVWI->LoadParameter("OuterIterationStoreSize", storesize, 10);
70 
71  // set a low tol for interfaceVWI
72  m_sessionVWI->LoadParameter("EigenvalueRelativeTol", m_eigRelTol, 1e-3);
73 
74  m_sessionVWI->LoadParameter("NeutralPointTolerance", m_neutralPointTol,
75  1e-4);
76  m_sessionVWI->LoadParameter("MaxOuterIterations", m_maxOuterIterations,
77  100);
78 
79  m_sessionVWI->LoadParameter("StartIteration", m_iterStart, 0);
80  m_sessionVWI->LoadParameter("EndIteration", m_iterEnd, 0);
81 
82  m_sessionVWI->LoadParameter("WaveForceMagStep", m_waveForceMagStep, 0.01);
83  m_sessionVWI->LoadParameter("DAlphaDWaveForceMag", m_dAlphaDWaveForceMag,
84  0.0);
85  m_sessionVWI->LoadParameter("MaxWaveForceMagIter", m_maxWaveForceMagIter,
86  1);
87  m_sessionVWI->LoadParameter("RollForceScale", m_rollForceScale, 1.0);
88 
89  if (m_sessionVWI->DefinesSolverInfo("DeltaFcnApprox"))
90  {
91  m_deltaFcnApprox = true;
92  m_sessionVWI->LoadParameter("DeltaFcnDecay", m_deltaFcnDecay,
93  1.0 / 500);
94  }
95  else
96  {
97  m_deltaFcnApprox = false;
98  m_deltaFcnDecay = 0.0;
99  }
100 
101  if (m_sessionVWI->DefinesSolverInfo("LinfPressureNorm"))
102  {
103  m_useLinfPressureNorm = true;
104  }
105  else
106  {
107  m_useLinfPressureNorm = false;
108  }
109 
110  if (m_sessionVWI->DefinesSolverInfo("INTERFACE"))
111  {
112  m_iterinterface = true;
113  }
114  else
115  {
116  m_iterinterface = false;
117  }
118 
119  if (m_sessionVWI->DefinesSolverInfo("MoveMeshToCriticalLayer"))
120  {
122  }
123  else
124  {
126  }
127 
128  m_alpha = Array<OneD, NekDouble>(storesize);
129  m_alpha[0] = m_sessionVWI->GetParameter("Alpha");
130  m_waveForceMag = Array<OneD, NekDouble>(storesize);
131  m_waveForceMag[0] = m_sessionVWI->GetParameter("WaveForceMag");
132 
133  m_leading_real_evl = Array<OneD, NekDouble>(storesize);
134  m_leading_imag_evl = Array<OneD, NekDouble>(storesize);
135 
136  if (m_sessionVWI->DefinesParameter("Relaxation"))
137  {
138  m_vwiRelaxation = m_sessionVWI->GetParameter("Relaxation");
139  // fix minimum number of iterations to be number of
140  // iterations required to make contribution of innitial
141  // forcing to 0.1
142  m_minInnerIterations = (int)(log(0.1) / log(m_vwiRelaxation));
143  }
144  else
145  {
146  m_vwiRelaxation = 0.0;
148  }
149 
150  // Initialise NS Roll solver
151  std::string IncCondFile(argv[argc - 1]);
152  IncCondFile += "_IncNSCond.xml";
153  std::vector<std::string> IncNSFilenames;
154  IncNSFilenames.push_back(meshfile);
155  IncNSFilenames.push_back(IncCondFile);
156 
157  // Create Incompressible NavierStokesSolver session reader.
159  argc, argv, IncNSFilenames, m_sessionVWI->GetComm());
161  std::string vEquation = m_sessionRoll->GetSolverInfo("SolverType");
163  vEquation, m_sessionRoll, m_graphRoll);
164  m_solverRoll->PrintSummary(cout);
165 
166  if (m_sessionRoll->DefinesSolverInfo("INTERFACE"))
167  {
168  m_sessionVWI->LoadParameter("EigenvalueRelativeTol", m_eigRelTol, 1e-2);
169  }
170 
171  int ncoeffs = m_solverRoll->UpdateFields()[0]->GetNcoeffs();
172  m_vwiForcing = Array<OneD, Array<OneD, NekDouble>>(4);
173  m_vwiForcing[0] = Array<OneD, NekDouble>(4 * ncoeffs);
174  for (int i = 1; i < 4; ++i)
175  {
176  m_vwiForcing[i] = m_vwiForcing[i - 1] + ncoeffs;
177  }
178 
179  // Fill forcing into m_vwiForcing
180  // Has forcing even been initialised yet?
181  // Vmath::Vcopy(ncoeffs,m_solverRoll->UpdateForces()[0]->GetCoeffs(),1,m_vwiForcing[0],1);
182  // Vmath::Vcopy(ncoeffs,m_solverRoll->UpdateForces()[1]->GetCoeffs(),1,m_vwiForcing[1],1);
183 
184  // Create AdvDiff Streak solver
185  std::string AdvDiffCondFile(argv[argc - 1]);
186  AdvDiffCondFile += "_AdvDiffCond.xml";
187  std::vector<std::string> AdvDiffFilenames;
188  AdvDiffFilenames.push_back(meshfile);
189  AdvDiffFilenames.push_back(AdvDiffCondFile);
190 
191  // Create AdvDiffusion session reader.
193  argc, argv, AdvDiffFilenames, m_sessionVWI->GetComm());
195 
196  // Initialise LinNS solver
197  std::string LinNSCondFile(argv[argc - 1]);
198  LinNSCondFile += "_LinNSCond.xml";
199  std::vector<std::string> LinNSFilenames;
200  LinNSFilenames.push_back(meshfile);
201  LinNSFilenames.push_back(LinNSCondFile);
202 
203  // Create Linearised NS stability session reader.
205  argc, argv, LinNSFilenames, m_sessionVWI->GetComm());
207 
208  // Set the initial beta value in stability to be equal to VWI file
209  std::string LZstr("LZ");
210  NekDouble LZ = 2 * M_PI / m_alpha[0];
211  cout << "Setting LZ in Linearised solver to " << LZ << endl;
212  m_sessionWave->SetParameter(LZstr, LZ);
213 
214  // Check for iteration type
215  if (m_sessionVWI->DefinesSolverInfo("VWIIterationType"))
216  {
217  std::string IterationTypeStr =
218  m_sessionVWI->GetSolverInfo("VWIIterationType");
219  for (int i = 0; i < (int)eVWIIterationTypeSize; ++i)
220  {
221  if (boost::iequals(VWIIterationTypeMap[i], IterationTypeStr))
222  {
224  break;
225  }
226  }
227  }
228  else
229  {
231  }
232 
233  // Check for restart
234  bool restart;
235  m_sessionVWI->MatchSolverInfo("RestartIteration", "True", restart, false);
236  if (restart)
237  {
238  switch (m_VWIIterationType)
239  {
240  case eFixedAlpha:
241  break;
242  case eFixedWaveForcing:
243  {
244  FILE *fp;
245  // Check for OuterIter.his file to read
246  if ((fp = fopen("OuterIter.his", "r")))
247  {
248  char buf[BUFSIZ];
249  std::vector<NekDouble> Alpha, Growth, Phase;
250  NekDouble alpha, growth, phase;
251  while (fgets(buf, BUFSIZ, fp))
252  {
253  sscanf(buf, "%*d:%lf%lf%lf", &alpha, &growth, &phase);
254  Alpha.push_back(alpha);
255  Growth.push_back(growth);
256  Phase.push_back(phase);
257  }
258 
259  m_nOuterIterations = Alpha.size();
260 
261  int nvals =
262  std::min(m_nOuterIterations, (int)m_alpha.size());
263 
264  for (int i = 0; i < nvals; ++i)
265  {
266  m_alpha[nvals - 1 - i] =
267  Alpha[m_nOuterIterations - nvals + i];
268  m_leading_real_evl[nvals - 1 - i] =
269  Growth[m_nOuterIterations - nvals + i];
270  m_leading_imag_evl[nvals - 1 - i] =
271  Phase[m_nOuterIterations - nvals + i];
272  }
273 
275  }
276  else
277  {
278  cout << " No File OuterIter.his to restart from" << endl;
279  }
280  }
281  break;
283  {
284  string nstr = boost::lexical_cast<std::string>(m_iterStart);
285  cout << "Restarting from iteration " << m_iterStart << endl;
286  std::string rstfile = "cp -f Save/" + m_sessionName + ".rst." +
287  nstr + " " + m_sessionName + ".rst";
288  cout << " " << rstfile << endl;
289  if (system(rstfile.c_str()))
290  {
291  ASSERTL0(false, rstfile.c_str());
292  }
293  std::string vwifile = "cp -f Save/" + m_sessionName + ".vwi." +
294  nstr + " " + m_sessionName + ".vwi";
295  cout << " " << vwifile << endl;
296  if (system(vwifile.c_str()))
297  {
298  ASSERTL0(false, vwifile.c_str());
299  }
300  }
301  break;
302  default:
303  ASSERTL0(false, "Unknown VWIITerationType in restart");
304  }
305  }
306 
307  // Check for ConveredSoln to update DAlphaDWaveForce
308  {
309  FILE *fp;
310  if ((fp = fopen("ConvergedSolns", "r")))
311  {
312  char buf[BUFSIZ];
313  std::vector<NekDouble> WaveForce, Alpha;
314  NekDouble waveforce, alpha;
315  while (fgets(buf, BUFSIZ, fp))
316  {
317  sscanf(buf, "%*d:%lf%lf", &waveforce, &alpha);
318  WaveForce.push_back(waveforce);
319  Alpha.push_back(alpha);
320  }
321 
322  if (Alpha.size() > 1)
323  {
324  int min_i = 0;
325  NekDouble min_alph = fabs(m_alpha[0] - Alpha[min_i]);
326  // find nearest point
327  for (int i = 1; i < Alpha.size(); ++i)
328  {
329  if (fabs(m_alpha[0] - Alpha[min_i]) < min_alph)
330  {
331  min_i = i;
332  min_alph = fabs(m_alpha[0] - Alpha[min_i]);
333  }
334  }
335 
336  // find next nearest point
337  int min_j = (min_i == 0) ? 1 : 0;
338  min_alph = fabs(m_alpha[0] - Alpha[min_j]);
339  for (int i = 0; i < Alpha.size(); ++i)
340  {
341  if (i != min_i)
342  {
343  if (fabs(m_alpha[0] - Alpha[min_j]) < min_alph)
344  {
345  min_j = i;
346  min_alph = fabs(m_alpha[0] - Alpha[min_j]);
347  }
348  }
349  }
350 
351  if (fabs(Alpha[min_i] - Alpha[min_j]) > 1e-4)
352  {
354  (Alpha[min_i] - Alpha[min_j]) /
355  (WaveForce[min_i] - WaveForce[min_j]);
356  }
357  }
358  }
359  }
360 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:144
static SessionReaderSharedPtr CreateInstance(int argc, char *argv[])
Creates an instance of the SessionReader class.
static MeshGraphSharedPtr Read(const LibUtilities::SessionReaderSharedPtr pSession, LibUtilities::DomainRangeShPtr rng=LibUtilities::NullDomainRangeShPtr, bool fillGraph=true)
Definition: MeshGraph.cpp:110
Array< OneD, Array< OneD, NekDouble > > m_vwiForcing
LibUtilities::SessionReaderSharedPtr m_sessionVWI
LibUtilities::SessionReaderSharedPtr m_sessionRoll
Array< OneD, NekDouble > m_alpha
< Leading imaginary eigenvalue
LibUtilities::SessionReaderSharedPtr m_sessionStreak
SpatialDomains::MeshGraphSharedPtr m_graphStreak
SpatialDomains::MeshGraphSharedPtr m_graphRoll
SpatialDomains::MeshGraphSharedPtr m_graphWave
Array< OneD, NekDouble > m_waveForceMag
LibUtilities::SessionReaderSharedPtr m_sessionWave
Array< OneD, NekDouble > m_leading_real_evl
Array< OneD, NekDouble > m_leading_imag_evl
< Leading real eigenvalue
EquationSystemSharedPtr m_solverRoll
std::shared_ptr< SessionReader > SessionReaderSharedPtr
EquationSystemFactory & GetEquationSystemFactory()
const std::string VWIIterationTypeMap[]
double NekDouble
scalarT< T > log(scalarT< T > in)
Definition: scalar.hpp:300

References ASSERTL0, Nektar::LibUtilities::SessionReader::CreateInstance(), Nektar::LibUtilities::NekFactory< tKey, tBase, tParam >::CreateInstance(), Nektar::eFixedAlpha, Nektar::eFixedAlphaWaveForcing, Nektar::eFixedWaveForcing, Nektar::eVWIIterationTypeSize, Nektar::SolverUtils::GetEquationSystemFactory(), tinysimd::log(), m_alpha, m_alphaStep, m_dAlphaDWaveForceMag, m_deltaFcnApprox, m_deltaFcnDecay, m_eigRelTol, m_graphRoll, m_graphStreak, m_graphWave, 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, Nektar::SpatialDomains::MeshGraph::Read(), UpdateAlpha(), and Nektar::VWIIterationTypeMap.

◆ ~VortexWaveInteraction()

Nektar::VortexWaveInteraction::~VortexWaveInteraction ( void  )

Definition at line 362 of file VortexWaveInteraction.cpp.

363 {
364  m_sessionVWI->Finalise();
365 }

References m_sessionVWI.

Member Function Documentation

◆ AppendEvlToFile() [1/2]

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

Definition at line 1164 of file VortexWaveInteraction.cpp.

1165 {
1166  FILE *fp;
1167  fp = fopen(file.c_str(), "a");
1168  fprintf(fp, "%d: %lf %16.12le %16.12le\n", n, m_alpha[0],
1170  fclose(fp);
1171 }

References m_alpha, m_leading_imag_evl, and m_leading_real_evl.

Referenced by DoFixedForcingIteration(), and main().

◆ AppendEvlToFile() [2/2]

void Nektar::VortexWaveInteraction::AppendEvlToFile ( std::string  file,
NekDouble  WaveForceMag 
)

Definition at line 1173 of file VortexWaveInteraction.cpp.

1174 {
1175  FILE *fp;
1176  fp = fopen(file.c_str(), "a");
1177  fprintf(fp, "%lf %lf %16.12le %16.12le\n", WaveForceMag, m_alpha[0],
1179  fclose(fp);
1180 }

References m_alpha, m_leading_imag_evl, and m_leading_real_evl.

◆ CalcL2ToLinfPressure()

void Nektar::VortexWaveInteraction::CalcL2ToLinfPressure ( void  )

Definition at line 1067 of file VortexWaveInteraction.cpp.

1068 {
1069 
1070  ExecuteWave();
1071 
1072  m_wavePressure->GetPlane(0)->BwdTrans(
1073  m_wavePressure->GetPlane(0)->GetCoeffs(),
1074  m_wavePressure->GetPlane(0)->UpdatePhys());
1075  m_wavePressure->GetPlane(1)->BwdTrans(
1076  m_wavePressure->GetPlane(1)->GetCoeffs(),
1077  m_wavePressure->GetPlane(1)->UpdatePhys());
1078 
1079  int npts = m_waveVelocities[0]->GetPlane(0)->GetNpoints();
1080  NekDouble Linf;
1081  Array<OneD, NekDouble> val(2 * npts, 0.0);
1082 
1083  Vmath::Vmul(npts, m_wavePressure->GetPlane(0)->UpdatePhys(), 1,
1084  m_wavePressure->GetPlane(0)->UpdatePhys(), 1, val, 1);
1085  Vmath::Vvtvp(npts, m_wavePressure->GetPlane(1)->UpdatePhys(), 1,
1086  m_wavePressure->GetPlane(1)->UpdatePhys(), 1, val, 1, val, 1);
1087  cout << "int P^2 " << m_wavePressure->GetPlane(0)->Integral(val) << endl;
1088  Vmath::Vsqrt(npts, val, 1, val, 1);
1089 
1090  Linf = Vmath::Vmax(npts, val, 1);
1091  cout << "Linf: " << Linf << endl;
1092 
1093  NekDouble l2, norm;
1094  l2 =
1095  m_wavePressure->GetPlane(0)->L2(m_wavePressure->GetPlane(0)->GetPhys());
1096  norm = l2 * l2;
1097  l2 =
1098  m_wavePressure->GetPlane(1)->L2(m_wavePressure->GetPlane(1)->GetPhys());
1099  norm += l2 * l2;
1100 
1101  Vmath::Fill(npts, 1.0, val, 1);
1102  NekDouble area = m_waveVelocities[0]->GetPlane(0)->Integral(val);
1103 
1104  l2 = sqrt(norm / area);
1105 
1106  cout << "L2: " << l2 << endl;
1107 
1108  cout << "Ratio Linf/L2: " << Linf / l2 << endl;
1109 }
Array< OneD, MultiRegions::ExpListSharedPtr > m_waveVelocities
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:534
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:209
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:574
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:45
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
Definition: Vmath.cpp:945
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:291

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

Referenced by main().

◆ CalcNonLinearWaveForce()

void Nektar::VortexWaveInteraction::CalcNonLinearWaveForce ( void  )

Definition at line 598 of file VortexWaveInteraction.cpp.

599 {
600  if (m_sessionRoll->DefinesSolverInfo("INTERFACE"))
601  {
602  static int cnt = 0;
603  string wavefile = m_sessionName + ".fld";
604  string movedmesh = m_sessionName + "_advPost_moved.xml";
605  string filestreak = m_sessionName + "_streak.fld";
606  char c[16] = "";
607  sprintf(c, "%d", cnt);
608  char c_alpha[16] = "";
609  sprintf(c_alpha, "%f", m_alpha[0]);
610  string syscall;
611  if (m_sessionVWI->GetSolverInfo("INTERFACE") == "phase")
612  {
613  string filePost = m_sessionName + "_advPost.xml";
614  syscall = "../../utilities/PostProcessing/Extras/FldCalcBCs " +
615  filePost + " " +
616  "meshhalf_pos_Spen_stability_moved.fld "
617  "meshhalf_pos_Spen_advPost_moved.fld " +
618  c_alpha + " > data_alpha0";
619  cout << syscall.c_str() << endl;
620  if (system(syscall.c_str()))
621  {
622  ASSERTL0(false, syscall.c_str());
623  }
624 
625  syscall = "cp -f meshhalf_pos_Spen_stability_moved_u_5.bc " +
626  m_sessionName + "_u_5.bc";
627  cout << syscall.c_str() << endl;
628  if (system(syscall.c_str()))
629  {
630  ASSERTL0(false, syscall.c_str());
631  }
632  syscall = "cp -f meshhalf_pos_Spen_stability_moved_v_5.bc " +
633  m_sessionName + "_v_5.bc";
634  cout << syscall.c_str() << endl;
635  if (system(syscall.c_str()))
636  {
637  ASSERTL0(false, syscall.c_str());
638  }
639  }
640  else
641  {
642  syscall = "../../utilities/PostProcessing/Extras/FldCalcBCs " +
643  movedmesh + " " + wavefile + " " + filestreak + " " +
644  c_alpha + " > datasub_" + c;
645  cout << syscall.c_str() << endl;
646  if (system(syscall.c_str()))
647  {
648  ASSERTL0(false, syscall.c_str());
649  }
650  }
651 
652  // relaxation for different alpha values? does it make sense?
653 
654  // save the wave
655  string wave_subalp = m_sessionName + "_wave_subalp_" + c + ".fld";
656  syscall = "cp -f " + wavefile + " " + wave_subalp;
657  cout << syscall.c_str() << endl;
658  if (system(syscall.c_str()))
659  {
660  ASSERTL0(false, syscall.c_str());
661  }
662  // FileRelaxation(3);
663  cnt++;
664  }
665  else
666  {
667  int npts = m_waveVelocities[0]->GetPlane(0)->GetNpoints();
668  int ncoeffs = m_waveVelocities[0]->GetPlane(0)->GetNcoeffs();
669  Array<OneD, NekDouble> val(npts), der1(2 * npts);
670  Array<OneD, NekDouble> der2 = der1 + npts;
671  Array<OneD, NekDouble> streak;
672  static int projectfield = -1;
673 
674  if (m_deltaFcnApprox)
675  {
676  streak = Array<OneD, NekDouble>(npts);
677  m_streakField[0]->BwdTrans(m_streakField[0]->GetCoeffs(), streak);
678  }
679 
680  // Set project field to be first field that has a Neumann
681  // boundary since this not impose any condition on the vertical
682  // boundaries Othersise set to zero.
683  if (projectfield == -1)
684  {
685  Array<OneD, const SpatialDomains::BoundaryConditionShPtr> BndConds;
686 
687  for (int i = 0; i < m_waveVelocities.size(); ++i)
688  {
689  BndConds = m_waveVelocities[i]->GetBndConditions();
690  for (int j = 0; j < BndConds.size(); ++j)
691  {
692  if (BndConds[j]->GetBoundaryConditionType() ==
694  {
695  projectfield = i;
696  break;
697  }
698  }
699  if (projectfield != -1)
700  {
701  break;
702  }
703  }
704  if (projectfield == -1)
705  {
706  cout << "using first field to project non-linear forcing which "
707  "imposes a Dirichlet condition"
708  << endl;
709  projectfield = 0;
710  }
711  }
712 
713  // Shift m_vwiForcing in case of relaxation
714  Vmath::Vcopy(ncoeffs, m_vwiForcing[0], 1, m_vwiForcing[2], 1);
715  Vmath::Vcopy(ncoeffs, m_vwiForcing[1], 1, m_vwiForcing[3], 1);
716 
717  // determine inverse of area normalised field.
718  m_wavePressure->GetPlane(0)->BwdTrans(
719  m_wavePressure->GetPlane(0)->GetCoeffs(),
720  m_wavePressure->GetPlane(0)->UpdatePhys());
721  m_wavePressure->GetPlane(1)->BwdTrans(
722  m_wavePressure->GetPlane(1)->GetCoeffs(),
723  m_wavePressure->GetPlane(1)->UpdatePhys());
724  NekDouble invnorm;
725 
727  {
728  Vmath::Vmul(npts, m_wavePressure->GetPlane(0)->UpdatePhys(), 1,
729  m_wavePressure->GetPlane(0)->UpdatePhys(), 1, der1, 1);
730  Vmath::Vvtvp(npts, m_wavePressure->GetPlane(1)->UpdatePhys(), 1,
731  m_wavePressure->GetPlane(1)->UpdatePhys(), 1, der1, 1,
732  der1, 1);
733  Vmath::Vsqrt(npts, der1, 1, der1, 1);
734 
735  NekDouble Linf = Vmath::Vmax(npts, der1, 1);
736 
737  invnorm = 1.0 / Linf;
738  }
739  else
740  {
741  // Determine normalisation of pressure so that |P|/A = 1
742  NekDouble l2;
743  l2 = m_wavePressure->GetPlane(0)->L2(
744  m_wavePressure->GetPlane(0)->GetPhys());
745  invnorm = l2 * l2;
746  l2 = m_wavePressure->GetPlane(1)->L2(
747  m_wavePressure->GetPlane(1)->GetPhys());
748  invnorm += l2 * l2;
749  Vmath::Fill(2 * npts, 1.0, der1, 1);
750  NekDouble area = m_waveVelocities[0]->GetPlane(0)->Integral(der1);
751  cout << "Area: " << area << endl;
752  invnorm = sqrt(area / invnorm);
753  }
754 
755  // Get hold of arrays.
756  m_waveVelocities[0]->GetPlane(0)->BwdTrans(
757  m_waveVelocities[0]->GetPlane(0)->GetCoeffs(),
758  m_waveVelocities[0]->GetPlane(0)->UpdatePhys());
759  Array<OneD, NekDouble> u_real =
760  m_waveVelocities[0]->GetPlane(0)->UpdatePhys();
761  Vmath::Smul(npts, invnorm, u_real, 1, u_real, 1);
762  m_waveVelocities[0]->GetPlane(1)->BwdTrans(
763  m_waveVelocities[0]->GetPlane(1)->GetCoeffs(),
764  m_waveVelocities[0]->GetPlane(1)->UpdatePhys());
765  Array<OneD, NekDouble> u_imag =
766  m_waveVelocities[0]->GetPlane(1)->UpdatePhys();
767  Vmath::Smul(npts, invnorm, u_imag, 1, u_imag, 1);
768  m_waveVelocities[1]->GetPlane(0)->BwdTrans(
769  m_waveVelocities[1]->GetPlane(0)->GetCoeffs(),
770  m_waveVelocities[1]->GetPlane(0)->UpdatePhys());
771  Array<OneD, NekDouble> v_real =
772  m_waveVelocities[1]->GetPlane(0)->UpdatePhys();
773  Vmath::Smul(npts, invnorm, v_real, 1, v_real, 1);
774  m_waveVelocities[1]->GetPlane(1)->BwdTrans(
775  m_waveVelocities[1]->GetPlane(1)->GetCoeffs(),
776  m_waveVelocities[1]->GetPlane(1)->UpdatePhys());
777  Array<OneD, NekDouble> v_imag =
778  m_waveVelocities[1]->GetPlane(1)->UpdatePhys();
779  Vmath::Smul(npts, invnorm, v_imag, 1, v_imag, 1);
780 
781  // normalise wave for output
782  Vmath::Smul(2 * ncoeffs, invnorm, m_waveVelocities[0]->UpdateCoeffs(),
783  1, m_waveVelocities[0]->UpdateCoeffs(), 1);
784  Vmath::Smul(2 * ncoeffs, invnorm, m_waveVelocities[1]->UpdateCoeffs(),
785  1, m_waveVelocities[1]->UpdateCoeffs(), 1);
786  Vmath::Smul(2 * ncoeffs, invnorm, m_waveVelocities[2]->UpdateCoeffs(),
787  1, m_waveVelocities[2]->UpdateCoeffs(), 1);
788 
789  // dump field
790  {
791  std::vector<std::string> variables(3);
792  std::vector<Array<OneD, NekDouble>> outfield(3);
793  variables[0] = "u_w";
794  variables[1] = "v_w";
795  variables[2] = "w_w";
796  outfield[0] = m_waveVelocities[0]->UpdateCoeffs();
797  outfield[1] = m_waveVelocities[1]->UpdateCoeffs();
798  outfield[2] = m_waveVelocities[2]->UpdateCoeffs();
799  std::string outname = m_sessionName + "_wave.fld";
800  m_solverRoll->WriteFld(outname, m_waveVelocities[0], outfield,
801  variables);
802  }
803 
804 #if 1
805  int ncoeffs_p = m_wavePressure->GetPlane(0)->GetNcoeffs();
806  Vmath::Smul(ncoeffs_p, invnorm,
807  m_wavePressure->GetPlane(0)->UpdateCoeffs(), 1,
808  m_wavePressure->GetPlane(0)->UpdateCoeffs(), 1);
809  Vmath::Smul(ncoeffs_p, invnorm,
810  m_wavePressure->GetPlane(1)->UpdateCoeffs(), 1,
811  m_wavePressure->GetPlane(1)->UpdateCoeffs(), 1);
812 #else
813  m_wavePressure->GetPlane(0)->BwdTrans(
814  m_wavePressure->GetPlane(0)->GetCoeffs(),
815  m_wavePressure->GetPlane(0)->UpdatePhys());
816  Vmath::Smul(npts, invnorm, m_wavePressure->GetPlane(0)->UpdatePhys(), 1,
817  m_wavePressure->GetPlane(0)->UpdatePhys(), 1);
818  m_wavePressure->GetPlane(0)->FwdTrans(
819  m_wavePressure->GetPlane(0)->UpdatePhys(),
820  m_wavePressure->GetPlane(0)->UpdateCoeffs());
821 
822  m_wavePressure->GetPlane(1)->BwdTrans(
823  m_wavePressure->GetPlane(1)->GetCoeffs(),
824  m_wavePressure->GetPlane(1)->UpdatePhys());
825  Vmath::Smul(npts, invnorm, m_wavePressure->GetPlane(1)->UpdatePhys(), 1,
826  m_wavePressure->GetPlane(1)->UpdatePhys(), 1);
827  m_wavePressure->GetPlane(1)->FwdTrans(
828  m_wavePressure->GetPlane(1)->UpdatePhys(),
829  m_wavePressure->GetPlane(1)->UpdateCoeffs());
830 #endif
831 
832  // Calculate non-linear terms for x and y directions
833  // d/dx(u u* + u* u)
834  Vmath::Vmul(npts, u_real, 1, u_real, 1, val, 1);
835  Vmath::Vvtvp(npts, u_imag, 1, u_imag, 1, val, 1, val, 1);
836  Vmath::Smul(npts, 2.0, val, 1, val, 1);
837  m_waveVelocities[0]->GetPlane(0)->PhysDeriv(0, val, der1);
838 
839  // d/dy(v u* + v* u)
840  Vmath::Vmul(npts, u_real, 1, v_real, 1, val, 1);
841  Vmath::Vvtvp(npts, u_imag, 1, v_imag, 1, val, 1, val, 1);
842  Vmath::Smul(npts, 2.0, val, 1, val, 1);
843  m_waveVelocities[0]->GetPlane(0)->PhysDeriv(1, val, der2);
844 
845  Vmath::Vadd(npts, der1, 1, der2, 1, der1, 1);
846 #if 1
847  m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(der1,
848  m_vwiForcing[0]);
849 #else
850  m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(
851  der1, m_vwiForcing[0]);
852 #endif
853  Vmath::Smul(ncoeffs, -m_waveForceMag[0], m_vwiForcing[0], 1,
854  m_vwiForcing[0], 1);
855 
856  // d/dx(u v* + u* v)
857  m_waveVelocities[0]->GetPlane(0)->PhysDeriv(0, val, der1);
858 
859  // d/dy(v v* + v* v)
860  Vmath::Vmul(npts, v_real, 1, v_real, 1, val, 1);
861  Vmath::Vvtvp(npts, v_imag, 1, v_imag, 1, val, 1, val, 1);
862  Vmath::Smul(npts, 2.0, val, 1, val, 1);
863  m_waveVelocities[0]->GetPlane(0)->PhysDeriv(1, val, der2);
864 
865  Vmath::Vadd(npts, der1, 1, der2, 1, der1, 1);
866 
867 #if 1
868  m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(der1,
869  m_vwiForcing[1]);
870 #else
871  m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(
872  der1, m_vwiForcing[1]);
873 #endif
874 
875  Vmath::Smul(ncoeffs, -m_waveForceMag[0], m_vwiForcing[1], 1,
876  m_vwiForcing[1], 1);
877 
878  // by default the symmetrization is on
879  bool symm = true;
880  m_sessionVWI->MatchSolverInfo("Symmetrization", "True", symm, true);
881 #if 0
882  if(symm== true )
883  {
884 
885  // Symmetrise forcing
886  //-> Get coordinates
887  Array<OneD, NekDouble> coord(2);
888  Array<OneD, NekDouble> coord_x(npts);
889  Array<OneD, NekDouble> coord_y(npts);
890 
891  //-> Impose symmetry (x -> -x + Lx/2, y-> -y) on coordinates
892  m_waveVelocities[0]->GetPlane(0)->GetCoords(coord_x,coord_y);
893  NekDouble xmax = Vmath::Vmax(npts,coord_x,1);
894  Vmath::Neg(npts,coord_x,1);
895  Vmath::Sadd(npts,xmax,coord_x,1,coord_x,1);
896  Vmath::Neg(npts,coord_y,1);
897 
898  int i, physoffset;
899 
900  //-> Obtain list of expansion element ids for each point.
901  Array<OneD, int> Eid(npts);
902  // This search may not be necessary every iteration
903  for(i = 0; i < npts; ++i)
904  {
905  coord[0] = coord_x[i];
906  coord[1] = coord_y[i];
907 
908  // Note this will not quite be symmetric.
909  Eid[i] = m_waveVelocities[0]->GetPlane(0)->GetExpIndex(coord,1e-6);
910  }
911 
912  // Interpolate field 0
913  m_waveVelocities[0]->GetPlane(0)->BwdTrans(m_vwiForcing[0],der1);
914  for(i = 0; i < npts; ++i)
915  {
916  physoffset = m_waveVelocities[0]->GetPlane(0)->GetPhys_Offset(Eid[i]);
917  coord[0] = coord_x[i];
918  coord[1] = coord_y[i];
919  der2 [i] = m_waveVelocities[0]->GetPlane(0)->GetExp(Eid[i])->PhysEvaluate(coord, der1 + physoffset);
920  }
921  //-> Average field 0
922  Vmath::Vsub(npts,der1,1,der2,1,der2,1);
923  Vmath::Smul(npts,0.5,der2,1,der2,1);
924 #if 1
925  m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(der2,m_vwiForcing[0]);
926 #else
927  m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(der2, m_vwiForcing[0]);
928 #endif
929 
930  //-> Interpoloate field 1
931  m_waveVelocities[0]->GetPlane(0)->BwdTrans(m_vwiForcing[1],der1);
932  for(i = 0; i < npts; ++i)
933  {
934  physoffset = m_waveVelocities[0]->GetPlane(0)->GetPhys_Offset(Eid[i]);
935  coord[0] = coord_x[i];
936  coord[1] = coord_y[i];
937  der2[i] = m_waveVelocities[0]->GetPlane(0)->GetExp(Eid[i])->PhysEvaluate(coord, der1 + physoffset);
938  }
939 
940  //-> Average field 1
941  Vmath::Vsub(npts,der1,1,der2,1,der2,1);
942  Vmath::Smul(npts,0.5,der2,1,der2,1);
943 #if 1
944  m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(der2,m_vwiForcing[1]);
945 #else
946  m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(der2, m_vwiForceing[1]);
947 #endif
948  }
949 #else
950  int i;
951  if (symm == true)
952  {
953  cout << "symmetrization is active" << endl;
954  static Array<OneD, int> index = GetReflectionIndex();
955 
956  m_waveVelocities[0]->GetPlane(0)->BwdTrans(m_vwiForcing[0], der1);
957 
958  for (i = 0; i < npts; ++i)
959  {
960  if (index[i] != -1)
961  {
962  val[i] = 0.5 * (der1[i] - der1[index[i]]);
963  }
964  }
965 #if 1
966  m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(
967  val, m_vwiForcing[0]);
968 #else
969  m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(
970  val, m_vwiForcing[0]);
971 #endif
972 
973  m_waveVelocities[0]->GetPlane(0)->BwdTrans(m_vwiForcing[1], der2);
974  for (i = 0; i < npts; ++i)
975  {
976  if (index[i] != -1)
977  {
978  val[i] = 0.5 * (der2[i] - der2[index[i]]);
979  }
980  }
981 #if 1
982  m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(
983  val, m_vwiForcing[1]);
984 #else
985  m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(
986  val, m_vwiForcing[1]);
987 #endif
988  }
989 
990  Vmath::Vmul(npts, der1, 1, der1, 1, val, 1);
991  Vmath::Vvtvp(npts, der2, 1, der2, 1, val, 1, val, 1);
992  Vmath::Vsqrt(npts, val, 1, val, 1);
993  cout << "F_Linf: " << Vmath::Vmax(npts, val, 1) << endl;
994 
995 #endif
996 
997  if (m_vwiRelaxation)
998  {
999  Vmath::Smul(ncoeffs, 1.0 - m_vwiRelaxation, m_vwiForcing[0], 1,
1000  m_vwiForcing[0], 1);
1001  Vmath::Svtvp(ncoeffs, m_vwiRelaxation, m_vwiForcing[2], 1,
1002  m_vwiForcing[0], 1, m_vwiForcing[0], 1);
1003 
1004  Vmath::Smul(ncoeffs, 1.0 - m_vwiRelaxation, m_vwiForcing[1], 1,
1005  m_vwiForcing[1], 1);
1006  Vmath::Svtvp(ncoeffs, m_vwiRelaxation, m_vwiForcing[3], 1,
1007  m_vwiForcing[1], 1, m_vwiForcing[1], 1);
1008  }
1009 
1010  if (m_deltaFcnApprox)
1011  {
1012  for (int j = 0; j < 2; ++j)
1013  {
1014 
1015  m_waveVelocities[projectfield]->GetPlane(0)->BwdTrans(
1016  m_vwiForcing[j], der1);
1017  for (int i = 0; i < npts; ++i)
1018  {
1019  der1[i] *= exp(-streak[i] * streak[i] / m_deltaFcnDecay);
1020  }
1021  m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(
1022  der1, m_vwiForcing[j]);
1023  }
1024  }
1025 
1026  // dump output
1027  std::vector<std::string> variables(4);
1028  std::vector<Array<OneD, NekDouble>> outfield(4);
1029  variables[0] = "u";
1030  variables[1] = "v";
1031  variables[2] = "pr";
1032  variables[3] = "pi";
1033  outfield[0] = m_vwiForcing[0];
1034  outfield[1] = m_vwiForcing[1];
1035  Array<OneD, NekDouble> soln(npts, 0.0);
1036  m_wavePressure->GetPlane(0)->BwdTrans(
1037  m_wavePressure->GetPlane(0)->GetCoeffs(),
1038  m_wavePressure->GetPlane(0)->UpdatePhys());
1039  outfield[2] = Array<OneD, NekDouble>(ncoeffs);
1040  m_waveVelocities[0]->GetPlane(0)->FwdTransLocalElmt(
1041  m_wavePressure->GetPlane(0)->GetPhys(), outfield[2]);
1042  m_wavePressure->GetPlane(1)->BwdTrans(
1043  m_wavePressure->GetPlane(1)->GetCoeffs(),
1044  m_wavePressure->GetPlane(1)->UpdatePhys());
1045 
1046  Vmath::Vmul(npts, m_wavePressure->GetPlane(0)->UpdatePhys(), 1,
1047  m_wavePressure->GetPlane(0)->UpdatePhys(), 1, val, 1);
1048  Vmath::Vvtvp(npts, m_wavePressure->GetPlane(1)->UpdatePhys(), 1,
1049  m_wavePressure->GetPlane(1)->UpdatePhys(), 1, val, 1, val,
1050  1);
1051  cout << "int P^2: " << m_wavePressure->GetPlane(0)->Integral(val)
1052  << endl;
1053  Vmath::Vsqrt(npts, val, 1, val, 1);
1054  cout << "PLinf: " << Vmath::Vmax(npts, val, 1) << endl;
1055 
1056  outfield[3] = Array<OneD, NekDouble>(ncoeffs);
1057  m_waveVelocities[1]->GetPlane(0)->FwdTransLocalElmt(
1058  m_wavePressure->GetPlane(1)->GetPhys(), outfield[3]);
1059 
1060  std::string outname = m_sessionName + ".vwi";
1061 
1062  m_solverRoll->WriteFld(outname, m_waveVelocities[0]->GetPlane(0),
1063  outfield, variables);
1064  }
1065 }
Array< OneD, MultiRegions::ExpListSharedPtr > m_streakField
Array< OneD, int > GetReflectionIndex(void)
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
Definition: Vmath.cpp:622
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:518
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:359
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.cpp:248
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha - x.
Definition: Vmath.cpp:384
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1255
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:419

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(), Vmath::Sadd(), Vmath::Smul(), tinysimd::sqrt(), Vmath::Svtvp(), Vmath::Vadd(), Vmath::Vcopy(), Vmath::Vmax(), Vmath::Vmul(), Vmath::Vsqrt(), Vmath::Vsub(), and Vmath::Vvtvp().

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

◆ CheckEigIsStationary()

bool Nektar::VortexWaveInteraction::CheckEigIsStationary ( bool  reset = false)

Definition at line 1707 of file VortexWaveInteraction.cpp.

1708 {
1709  static NekDouble previous_real_evl = -1.0;
1710  static NekDouble previous_imag_evl = -1.0;
1711  static int min_iter = 0;
1712 
1713  if (reset)
1714  {
1715  previous_real_evl = -1.0;
1716  min_iter = 0;
1717  }
1718 
1719  if (previous_real_evl == -1.0 || min_iter < m_minInnerIterations)
1720  {
1721  previous_real_evl = m_leading_real_evl[0];
1722  previous_imag_evl = m_leading_imag_evl[0];
1723  min_iter++;
1724  return false;
1725  }
1726 
1727  cout << "Growth tolerance: "
1728  << fabs((m_leading_real_evl[0] - previous_real_evl) /
1729  m_leading_real_evl[0])
1730  << endl;
1731  cout << "Phase tolerance: "
1732  << fabs((m_leading_imag_evl[0] - previous_imag_evl) /
1733  m_leading_imag_evl[0])
1734  << endl;
1735 
1736  // See if real and imaginary growth have converged to with m_eigRelTol
1737  if ((fabs((m_leading_real_evl[0] - previous_real_evl) /
1739  (fabs(m_leading_real_evl[0]) < 1e-6))
1740  {
1741  previous_real_evl = m_leading_real_evl[0];
1742  previous_imag_evl = m_leading_imag_evl[0];
1743  if ((fabs((m_leading_imag_evl[0] - previous_imag_evl) /
1745  (fabs(m_leading_imag_evl[0]) < 1e-6))
1746  {
1747  return true;
1748  }
1749  }
1750  else
1751  {
1752  if (fabs(m_leading_imag_evl[0]) > 1e-2)
1753  {
1754  cout << "Warning: imaginary eigenvalue is greater than 1e-2"
1755  << endl;
1756  }
1757  previous_real_evl = m_leading_real_evl[0];
1758  previous_imag_evl = m_leading_imag_evl[0];
1759  return false;
1760  }
1761  return false;
1762 }

References m_eigRelTol, m_leading_imag_evl, m_leading_real_evl, and m_minInnerIterations.

Referenced by DoFixedForcingIteration().

◆ CheckIfAtNeutralPoint()

bool Nektar::VortexWaveInteraction::CheckIfAtNeutralPoint ( void  )

Definition at line 1766 of file VortexWaveInteraction.cpp.

1767 {
1768  bool returnval = false;
1769  if (m_sessionRoll->DefinesSolverInfo("INTERFACE"))
1770  {
1771  if (m_sessionVWI->GetSolverInfo("INTERFACE") == "phase")
1772  {
1773  if (abs(m_leading_real_evl[0]) < 1e-4)
1774  {
1775  returnval = true;
1776  }
1777  }
1778  else
1779  {
1780  if (abs(m_leading_real_evl[0]) < 1e-4 &&
1781  abs(m_leading_imag_evl[0]) < 2e-6)
1782  {
1783  returnval = true;
1784  }
1785  }
1786  }
1787  else
1788  {
1789  if ((m_leading_real_evl[0] * m_leading_real_evl[0] +
1792  {
1793  returnval = true;
1794  }
1795  }
1796  if (fabs(m_leading_imag_evl[0]) > 1e-2)
1797  {
1798  cout << "Warning: imaginary eigenvalue is greater than 1e-2" << endl;
1799  }
1800 
1801  return returnval;
1802 }
scalarT< T > abs(scalarT< T > in)
Definition: scalar.hpp:295

References tinysimd::abs(), m_leading_imag_evl, m_leading_real_evl, m_neutralPointTol, m_sessionRoll, and m_sessionVWI.

Referenced by DoFixedForcingIteration().

◆ CopyFile()

void Nektar::VortexWaveInteraction::CopyFile ( std::string  file1end,
std::string  file2end 
)

Definition at line 1152 of file VortexWaveInteraction.cpp.

1153 {
1154  string cpfile1 = m_sessionName + file1end;
1155  string cpfile2 = m_sessionName + file2end;
1156  string syscall = "cp -f " + cpfile1 + " " + cpfile2;
1157 
1158  if (system(syscall.c_str()))
1159  {
1160  ASSERTL0(false, syscall.c_str());
1161  }
1162 }

References ASSERTL0, and m_sessionName.

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

◆ ExecuteLoop()

void Nektar::VortexWaveInteraction::ExecuteLoop ( bool  CalcWaveForce = true)

Definition at line 1201 of file VortexWaveInteraction.cpp.

1202 {
1203 
1204  // the global info has to be written in the
1205  // roll session file to use the interface loop
1206  if (m_sessionRoll->DefinesSolverInfo("INTERFACE"))
1207  {
1208  static int cnt = 0;
1209  bool skiprollstreak = false;
1210  if (cnt == 0 && m_sessionVWI->GetParameter("rollstreakfromit") == 1)
1211  {
1212  skiprollstreak = true;
1213  cout << "skip roll-streak at the first iteration" << endl;
1214  }
1215 
1216  if (skiprollstreak != true)
1217  {
1218 
1219  LibUtilities::NekManager<
1220  MultiRegions::GlobalLinSysKey,
1221  MultiRegions::GlobalLinSys>::EnableManagement("GlobalLinSys");
1222  ExecuteRoll();
1223  LibUtilities::NekManager<
1224  MultiRegions::GlobalLinSysKey,
1225  MultiRegions::GlobalLinSys>::DisableManagement("GlobalLinSys");
1226 #ifndef _WIN32
1227  sleep(3);
1228 #endif
1229  ExecuteStreak();
1230 #ifndef _WIN32
1231  sleep(3);
1232 #endif
1233  }
1234 
1235  string syscall;
1236  char c[16] = "";
1237  string movedmesh = m_sessionName + "_advPost_moved.xml";
1238  string movedinterpmesh = m_sessionName + "_interp_moved.xml";
1239  // rewrite the Rollsessionfile (we start from the waleffe forcing)
1240  // string meshbndjumps = m_sessionName +"_bndjumps.xml";
1241  // if(cnt==0)
1242  //{
1243  // take the conditions tag from meshbndjumps and copy into
1244  // the rolls session file
1245  //}
1246 
1247  sprintf(c, "%d", cnt);
1248  // save old roll solution
1249  string oldroll = m_sessionName + "_roll_" + c + ".fld";
1250  syscall = "cp -f " + m_sessionName + "-Base.fld" + " " + oldroll;
1251  cout << syscall.c_str() << endl;
1252  if (system(syscall.c_str()))
1253  {
1254  ASSERTL0(false, syscall.c_str());
1255  }
1256  // define file names
1257  string filePost = m_sessionName + "_advPost.xml";
1258  string filestreak = m_sessionName + "_streak.fld";
1259  string filewave = m_sessionName + "_wave.fld";
1260  string filewavepressure = m_sessionName + "_wave_p_split.fld";
1261  string fileinterp = m_sessionName + "_interp.xml";
1262  string interpstreak = m_sessionName + "_interpstreak_" + c + ".fld";
1263  string interwavepressure =
1264  m_sessionName + "_wave_p_split_interp_" + c + ".fld";
1265  char alpchar[16] = "";
1266  cout << "alpha = " << m_alpha[0] << endl;
1267  sprintf(alpchar, "%f", m_alpha[0]);
1268 
1269  if (m_sessionVWI->GetSolverInfo("INTERFACE") != "phase")
1270  {
1271  cout << "zerophase" << endl;
1272 
1273  syscall = "../../utilities/PostProcessing/Extras/MoveMesh " +
1274  filePost + " " + filestreak + " " + fileinterp + " " +
1275  alpchar;
1276 
1277  cout << syscall.c_str() << endl;
1278  if (system(syscall.c_str()))
1279  {
1280  ASSERTL0(false, syscall.c_str());
1281  }
1282 
1283  // move the advPost mesh (remark update alpha!!!)
1284  syscall = "../../utilities/PostProcessing/Extras/MoveMesh " +
1285  filePost + " " + filestreak + " " + filePost + " " +
1286  alpchar;
1287  cout << syscall.c_str() << endl;
1288  if (system(syscall.c_str()))
1289  {
1290  ASSERTL0(false, syscall.c_str());
1291  }
1292 
1293  // save oldstreak
1294  string oldstreak = m_sessionName + "_streak_" + c + ".fld";
1295  syscall = "cp -f " + filestreak + " " + oldstreak;
1296  cout << syscall.c_str() << endl;
1297  if (system(syscall.c_str()))
1298  {
1299  ASSERTL0(false, syscall.c_str());
1300  }
1301 
1302  // interpolate the streak field into the new mesh
1303  string movedmesh = m_sessionName + "_advPost_moved.xml";
1304  string movedinterpmesh = m_sessionName + "_interp_moved.xml";
1305 
1306  // create the interp streak
1307 
1308  syscall = "../../utilities/PostProcessing/Extras/FieldToField " +
1309  fileinterp + " " + filestreak + " " + movedinterpmesh +
1310  " " + interpstreak;
1311 
1312  cout << syscall.c_str() << endl;
1313  if (system(syscall.c_str()))
1314  {
1315  ASSERTL0(false, syscall.c_str());
1316  }
1317 
1318  // save the old mesh
1319  string meshfile = m_sessionName + ".xml";
1320  string meshold = m_sessionName + "_" + c + ".xml";
1321  syscall = "cp -f " + meshfile + " " + meshold;
1322  cout << syscall.c_str() << endl;
1323  if (system(syscall.c_str()))
1324  {
1325  ASSERTL0(false, syscall.c_str());
1326  }
1327 
1328  // overwriting the meshfile with the new mesh
1329  syscall = "cp -f " + movedmesh + " " + meshfile;
1330  cout << syscall.c_str() << endl;
1331  if (system(syscall.c_str()))
1332  {
1333  ASSERTL0(false, syscall.c_str());
1334  }
1335 
1336  // overwriting the streak file!!
1337  syscall = "cp -f " + interpstreak + " " + filestreak;
1338  cout << syscall.c_str() << endl;
1339  if (system(syscall.c_str()))
1340  {
1341  ASSERTL0(false, syscall.c_str());
1342  }
1343 
1344  // calculate the wave
1345  ExecuteWave();
1346 
1347  // save the wave field:
1348  string oldwave = m_sessionName + "_wave_" + c + ".fld";
1349  syscall = "cp -f " + m_sessionName + ".fld" + " " + oldwave;
1350  cout << syscall.c_str() << endl;
1351  if (system(syscall.c_str()))
1352  {
1353  ASSERTL0(false, syscall.c_str());
1354  }
1355 
1356  // save old jump conditions:
1357  string ujump = m_sessionName + "_u_5.bc";
1358  syscall = "cp -f " + ujump + " " + m_sessionName + "_u_5.bc_" + c;
1359  cout << syscall.c_str() << endl;
1360  if (system(syscall.c_str()))
1361  {
1362  ASSERTL0(false, syscall.c_str());
1363  }
1364 
1365  string vjump = m_sessionName + "_v_5.bc";
1366  syscall = "cp -f " + vjump + " " + m_sessionName + "_v_5.bc_" + c;
1367  cout << syscall.c_str() << endl;
1368  if (system(syscall.c_str()))
1369  {
1370  ASSERTL0(false, syscall.c_str());
1371  }
1372  cnt++;
1373 
1374  // use relaxation
1375  if (GetVWIIterationType() !=
1377  {
1378  // the critical layer should be the bnd region 3
1379  // int reg =3;
1380  // FileRelaxation(reg);
1381  }
1382  char c1[16] = "";
1383  sprintf(c1, "%d", cnt);
1384  // calculate the jump conditions
1385  string wavefile = m_sessionName + ".fld";
1386  syscall = "../../utilities/PostProcessing/Extras/FldCalcBCs " +
1387  movedmesh + " " + wavefile + " " + interpstreak +
1388  "> data" + c1;
1389  cout << syscall.c_str() << endl;
1390  if (system(syscall.c_str()))
1391  {
1392  ASSERTL0(false, syscall.c_str());
1393  }
1394 
1395  // move the new name_interp_moved.xml into name_interp.xml
1396  syscall = "cp -f " + movedinterpmesh + " " + fileinterp;
1397  cout << syscall.c_str() << endl;
1398  if (system(syscall.c_str()))
1399  {
1400  ASSERTL0(false, syscall.c_str());
1401  }
1402  // move the new name_advPost_moved.xml into name_advPost.xml
1403  syscall = "cp -f " + movedmesh + " " + filePost;
1404  cout << syscall.c_str() << endl;
1405  if (system(syscall.c_str()))
1406  {
1407  ASSERTL0(false, syscall.c_str());
1408  }
1409  }
1410  else if (m_sessionVWI->GetSolverInfo("INTERFACE") == "phase")
1411  {
1412  cout << "phase" << endl;
1413  // determine cr:
1414  NekDouble cr;
1415  string cr_str;
1416  stringstream st;
1417 
1418  // calculate the wave
1419  ExecuteWave();
1420 
1421  // save oldstreak
1422  string oldstreak = m_sessionName + "_streak_" + c + ".fld";
1423  syscall = "cp -f " + filestreak + " " + oldstreak;
1424  cout << syscall.c_str() << endl;
1425  if (system(syscall.c_str()))
1426  {
1427  ASSERTL0(false, syscall.c_str());
1428  }
1429 
1430  // save wave
1431  syscall = "cp -f " + m_sessionName + ".fld" + " " + filewave;
1432  cout << syscall.c_str() << endl;
1433  if (system(syscall.c_str()))
1434  {
1435  ASSERTL0(false, syscall.c_str());
1436  }
1437  // save the old mesh
1438  string meshfile = m_sessionName + ".xml";
1439  string meshold = m_sessionName + "_" + c + ".xml";
1440  syscall = "cp -f " + meshfile + " " + meshold;
1441  cout << syscall.c_str() << endl;
1442  if (system(syscall.c_str()))
1443  {
1444  ASSERTL0(false, syscall.c_str());
1445  }
1446 
1447  // save the oldwave field:
1448  string oldwave = m_sessionName + "_wave_" + c + ".fld";
1449  syscall = "cp -f " + m_sessionName + ".fld" + " " + oldwave;
1450  cout << syscall.c_str() << endl;
1451  if (system(syscall.c_str()))
1452  {
1453  ASSERTL0(false, syscall.c_str());
1454  }
1455 
1456  // save old jump conditions:
1457  string ujump = m_sessionName + "_u_5.bc";
1458  syscall = "cp -f " + ujump + " " + m_sessionName + "_u_5.bc_" + c;
1459  cout << syscall.c_str() << endl;
1460  if (system(syscall.c_str()))
1461  {
1462  ASSERTL0(false, syscall.c_str());
1463  }
1464 
1465  string vjump = m_sessionName + "_v_5.bc";
1466  syscall = "cp -f " + vjump + " " + m_sessionName + "_v_5.bc_" + c;
1467  cout << syscall.c_str() << endl;
1468  if (system(syscall.c_str()))
1469  {
1470  ASSERTL0(false, syscall.c_str());
1471  }
1472  cnt++;
1473 
1474  cr = m_leading_imag_evl[0] / m_alpha[0];
1475  st << cr;
1476  cr_str = st.str();
1477  cout << "cr=" << cr_str << endl;
1478  // NB -g or NOT!!!
1479  // move the mesh around the critical layer
1480  syscall = "../../utilities/PostProcessing/Extras/MoveMesh " +
1481  filePost + " " + filestreak + " " + fileinterp + " " +
1482  alpchar + " " + cr_str;
1483 
1484  cout << syscall.c_str() << endl;
1485  if (system(syscall.c_str()))
1486  {
1487  ASSERTL0(false, syscall.c_str());
1488  }
1489  // NB -g or NOT!!!
1490  // move the advPost mesh (remark update alpha!!!)
1491  syscall = "../../utilities/PostProcessing/Extras/MoveMesh " +
1492  filePost + " " + filestreak + " " + filePost + " " +
1493  alpchar + " " + cr_str;
1494  cout << syscall.c_str() << endl;
1495  if (system(syscall.c_str()))
1496  {
1497  ASSERTL0(false, syscall.c_str());
1498  }
1499 
1500  // interp streak into the new mesh
1501  syscall = "../../utilities/PostProcessing/Extras/FieldToField " +
1502  fileinterp + " " + filestreak + " " + movedinterpmesh +
1503  " " + interpstreak;
1504 
1505  cout << syscall.c_str() << endl;
1506  if (system(syscall.c_str()))
1507  {
1508  ASSERTL0(false, syscall.c_str());
1509  }
1510 
1511  // split wave sol
1512  syscall = "../../utilities/PostProcessing/Extras/SplitFld " +
1513  filePost + " " + filewave;
1514 
1515  cout << syscall.c_str() << endl;
1516  if (system(syscall.c_str()))
1517  {
1518  ASSERTL0(false, syscall.c_str());
1519  }
1520  // interp wave
1521  syscall = "../../utilities/PostProcessing/Extras/FieldToField " +
1522  filePost + " " + filewavepressure + " " + movedmesh +
1523  " " + interwavepressure;
1524 
1525  cout << syscall.c_str() << endl;
1526  if (system(syscall.c_str()))
1527  {
1528  ASSERTL0(false, syscall.c_str());
1529  }
1530 
1531  // use relaxation
1532  if (GetVWIIterationType() !=
1534  {
1535  // the critical layer should be the bnd region 3
1536  // int reg =3;
1537  // FileRelaxation(reg);
1538  }
1539  char c1[16] = "";
1540  sprintf(c1, "%d", cnt);
1541 
1542  // cp wavepressure to m_sessionName.fld(to get
1543  // the right bcs names using FldCalcBCs)
1544  syscall =
1545  "cp -f " + interwavepressure + " " + m_sessionName + ".fld";
1546  cout << syscall.c_str() << endl;
1547  if (system(syscall.c_str()))
1548  {
1549  ASSERTL0(false, syscall.c_str());
1550  }
1551 
1552  // calculate the jump conditions
1553  // NB -g or NOT!!!
1554  syscall = "../../utilities/PostProcessing/Extras/FldCalcBCs " +
1555  movedmesh + " " + m_sessionName + ".fld" + " " +
1556  interpstreak + "> data" + c1;
1557  cout << syscall.c_str() << endl;
1558  if (system(syscall.c_str()))
1559  {
1560  ASSERTL0(false, syscall.c_str());
1561  }
1562 
1563  // overwriting the meshfile with the new mesh
1564  syscall = "cp -f " + movedmesh + " " + meshfile;
1565  cout << syscall.c_str() << endl;
1566  if (system(syscall.c_str()))
1567  {
1568  ASSERTL0(false, syscall.c_str());
1569  }
1570 
1571  // overwriting the streak file!! (maybe is useless)
1572  syscall = "cp -f " + interpstreak + " " + filestreak;
1573  cout << syscall.c_str() << endl;
1574  if (system(syscall.c_str()))
1575  {
1576  ASSERTL0(false, syscall.c_str());
1577  }
1578  // move the new name_interp_moved.xml into name_interp.xml
1579  syscall = "cp -f " + movedinterpmesh + " " + fileinterp;
1580  cout << syscall.c_str() << endl;
1581  if (system(syscall.c_str()))
1582  {
1583  ASSERTL0(false, syscall.c_str());
1584  }
1585  // move the new name_advPost_moved.xml into name_advPost.xml
1586  syscall = "cp -f " + movedmesh + " " + filePost;
1587  cout << syscall.c_str() << endl;
1588  if (system(syscall.c_str()))
1589  {
1590  ASSERTL0(false, syscall.c_str());
1591  }
1592  }
1593  }
1594  else
1595  {
1596  LibUtilities::NekManager<
1597  MultiRegions::GlobalLinSysKey,
1598  MultiRegions::GlobalLinSys>::EnableManagement("GlobalLinSys");
1599  ExecuteRoll();
1600  LibUtilities::NekManager<
1601  MultiRegions::GlobalLinSysKey,
1602  MultiRegions::GlobalLinSys>::DisableManagement("GlobalLinSys");
1603 
1604 #ifndef _WIN32
1605  sleep(3);
1606 #endif
1607  ExecuteStreak();
1608 #ifndef _WIN32
1609  sleep(3);
1610 #endif
1611 
1613  {
1614  string syscall;
1615  char alpchar[16] = "";
1616  sprintf(alpchar, "%f", m_alpha[0]);
1617 
1618  string filePost = m_sessionName + "_advPost.xml";
1619  string filestreak = m_sessionName + "_streak.fld";
1620  string filewave = m_sessionName + "_wave.fld";
1621  string filewavepressure = m_sessionName + "_wave_p_split.fld";
1622  string fileinterp = m_sessionName + "_interp.xml";
1623  string interpstreak = m_sessionName + "_interpstreak.fld";
1624  string interwavepressure =
1625  m_sessionName + "_wave_p_split_interp.fld";
1626  syscall = "../../utilities/PostProcessing/Extras/MoveMesh " +
1627  filePost + " " + filestreak + " " + fileinterp + " " +
1628  alpchar;
1629 
1630  cout << syscall.c_str() << endl;
1631  if (system(syscall.c_str()))
1632  {
1633  ASSERTL0(false, syscall.c_str());
1634  }
1635 
1636  // move the advPost mesh (remark update alpha!!!)
1637  syscall = "../../utilities/PostProcessing/Extras/MoveMesh " +
1638  filePost + " " + filestreak + " " + filePost + " " +
1639  alpchar;
1640  cout << syscall.c_str() << endl;
1641  if (system(syscall.c_str()))
1642  {
1643  ASSERTL0(false, syscall.c_str());
1644  }
1645 
1646  // save oldstreak
1647  string oldstreak = m_sessionName + "_streak.fld";
1648  syscall = "cp -f " + filestreak + " " + oldstreak;
1649  cout << syscall.c_str() << endl;
1650  if (system(syscall.c_str()))
1651  {
1652  ASSERTL0(false, syscall.c_str());
1653  }
1654 
1655  // interpolate the streak field into the new mesh
1656  string movedmesh = m_sessionName + "_advPost_moved.xml";
1657  string movedinterpmesh = m_sessionName + "_interp_moved.xml";
1658 
1659  // create the interp streak
1660 
1661  syscall = "../../utilities/PostProcessing/Extras/FieldToField " +
1662  fileinterp + " " + filestreak + " " + movedinterpmesh +
1663  " " + interpstreak;
1664 
1665  cout << syscall.c_str() << endl;
1666  if (system(syscall.c_str()))
1667  {
1668  ASSERTL0(false, syscall.c_str());
1669  }
1670 
1671  // save the old mesh
1672  string meshfile = m_sessionName + ".xml";
1673  string meshold = m_sessionName + ".xml";
1674  syscall = "cp -f " + meshfile + " " + meshold;
1675  cout << syscall.c_str() << endl;
1676  if (system(syscall.c_str()))
1677  {
1678  ASSERTL0(false, syscall.c_str());
1679  }
1680 
1681  // overwriting the meshfile with the new mesh
1682  syscall = "cp -f " + movedmesh + " " + meshfile;
1683  cout << syscall.c_str() << endl;
1684  if (system(syscall.c_str()))
1685  {
1686  ASSERTL0(false, syscall.c_str());
1687  }
1688 
1689  // overwriting the streak file!!
1690  syscall = "cp -f " + interpstreak + " " + filestreak;
1691  cout << syscall.c_str() << endl;
1692  if (system(syscall.c_str()))
1693  {
1694  ASSERTL0(false, syscall.c_str());
1695  }
1696  }
1697 
1698  ExecuteWave();
1699 
1700  if (CalcWaveForce)
1701  {
1703  }
1704  }
1705 }
VWIIterationType GetVWIIterationType(void)
@ eFixedWaveForcingWithSubIterationOnAlpha

References ASSERTL0, CalcNonLinearWaveForce(), Nektar::eFixedWaveForcingWithSubIterationOnAlpha, ExecuteRoll(), ExecuteStreak(), ExecuteWave(), GetVWIIterationType(), m_alpha, m_leading_imag_evl, m_moveMeshToCriticalLayer, m_sessionName, m_sessionRoll, and m_sessionVWI.

Referenced by DoFixedForcingIteration(), and main().

◆ ExecuteRoll()

void Nektar::VortexWaveInteraction::ExecuteRoll ( void  )

Definition at line 367 of file VortexWaveInteraction.cpp.

368 {
369  // set up the equation system to update the mesh
370  if (m_sessionRoll->DefinesSolverInfo("INTERFACE"))
371  {
372  string vEquation = m_sessionRoll->GetSolverInfo("solvertype");
373  EquationSystemSharedPtr solverRoll =
375  m_graphRoll);
376  // The forcing terms are inserted as N bcs
377  // Execute Roll
378  cout << "Executing Roll solver" << endl;
379  solverRoll->DoInitialise();
380  solverRoll->DoSolve();
381  solverRoll->Output();
382  m_rollField = solverRoll->UpdateFields();
383  for (int g = 0; g < solverRoll->GetNvariables(); ++g)
384  {
385  NekDouble vL2Error = solverRoll->L2Error(g, false);
386  NekDouble vLinfError = solverRoll->LinfError(g);
387  cout << "L 2 error (variable " << solverRoll->GetVariable(g)
388  << ") : " << vL2Error << endl;
389  cout << "L inf error (variable " << solverRoll->GetVariable(g)
390  << ") : " << vLinfError << endl;
391  }
392  }
393  else
394  {
396  {
397  string vEquation = m_sessionRoll->GetSolverInfo("solvertype");
398  EquationSystemSharedPtr solverRoll =
400  vEquation, m_sessionRoll, m_graphRoll);
401  }
402  else
403  {
404  const int npoints = m_solverRoll->GetNpoints();
405  const int ncoeffs = m_solverRoll->GetNcoeffs();
406 
407  static int init = 1;
408  if (init)
409  {
410  m_solverRoll->DoInitialise();
412  std::dynamic_pointer_cast<SolverUtils::ForcingProgrammatic>(
413  GetForcingFactory().CreateInstance(
414  "Programmatic", m_sessionRoll, m_solverRoll,
415  m_solverRoll->UpdateFields(),
416  m_solverRoll->UpdateFields().size() - 1, 0));
417 
418  std::vector<std::string> vFieldNames =
419  m_sessionRoll->GetVariables();
420  vFieldNames.erase(vFieldNames.end() - 1);
421 
422  m_solverRoll->GetFunction("BodyForce")
423  ->Evaluate(vFieldNames, m_vwiForcingObj->UpdateForces());
424 
425  // Scale forcing
426  for (int i = 0; i < m_vwiForcingObj->UpdateForces().size(); ++i)
427  {
428  m_solverRoll->UpdateFields()[0]->FwdTrans(
429  m_vwiForcingObj->UpdateForces()[i],
430  m_vwiForcing[2 + i]);
431  Vmath::Smul(npoints, m_rollForceScale,
432  m_vwiForcingObj->UpdateForces()[i], 1,
433  m_vwiForcingObj->UpdateForces()[i], 1);
434  }
435 
437  m_solverRoll->as<IncNavierStokes>();
438  ins->AddForcing(m_vwiForcingObj);
439 
440  init = 0;
441  }
442  else // use internal definition of forcing in m_vwiForcing
443  {
444  // Scale forcing
445  for (int i = 0; i < m_vwiForcingObj->UpdateForces().size(); ++i)
446  {
447  m_solverRoll->UpdateFields()[i]->BwdTrans(
448  m_vwiForcing[i], m_vwiForcingObj->UpdateForces()[i]);
449  Vmath::Smul(npoints, m_rollForceScale,
450  m_vwiForcingObj->UpdateForces()[i], 1,
451  m_vwiForcingObj->UpdateForces()[i], 1);
452  }
453 
454  // Shift m_vwiForcing for new restart in case of relaxation
455  Vmath::Vcopy(ncoeffs, m_vwiForcing[0], 1, m_vwiForcing[2], 1);
456  Vmath::Vcopy(ncoeffs, m_vwiForcing[1], 1, m_vwiForcing[3], 1);
457  }
458  }
459 
460  // Execute Roll
461  cout << "Executing Roll solver" << endl;
462  m_solverRoll->DoSolve();
463  m_solverRoll->Output();
464  m_rollField = m_solverRoll->UpdateFields();
465  for (int g = 0; g < m_solverRoll->GetNvariables(); ++g)
466  {
467  NekDouble vL2Error = m_solverRoll->L2Error(g, false);
468  NekDouble vLinfError = m_solverRoll->LinfError(g);
469  cout << "L 2 error (variable " << m_solverRoll->GetVariable(g)
470  << ") : " << vL2Error << endl;
471  cout << "L inf error (variable " << m_solverRoll->GetVariable(g)
472  << ") : " << vLinfError << endl;
473  }
474  }
475 
476  // Copy .fld file to .rst and base.fld
477  cout << "Executing cp -f session.fld session.rst" << endl;
478  CopyFile(".fld", ".rst");
479 
480  // Write out data into base flow with variable Vx,Vy
481  cout << "Writing data to session-Base.fld" << endl;
482 
483  std::vector<std::string> variables(2);
484  variables[0] = "Vx";
485  variables[1] = "Vy";
486  std::vector<Array<OneD, NekDouble>> outfield(2);
487  outfield[0] = m_solverRoll->UpdateFields()[0]->UpdateCoeffs();
488  outfield[1] = m_solverRoll->UpdateFields()[1]->UpdateCoeffs();
489  std::string outname = m_sessionName + "-Base.fld";
490  m_solverRoll->WriteFld(outname, m_solverRoll->UpdateFields()[0], outfield,
491  variables);
492 }
void CopyFile(std::string file1end, std::string file2end)
SolverUtils::ForcingProgrammaticSharedPtr m_vwiForcingObj
Array< OneD, MultiRegions::ExpListSharedPtr > m_rollField
std::shared_ptr< EquationSystem > EquationSystemSharedPtr
A shared pointer to an EquationSystem object.
ForcingFactory & GetForcingFactory()
Declaration of the forcing factory singleton.
Definition: Forcing.cpp:44
std::shared_ptr< IncNavierStokes > IncNavierStokesSharedPtr

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

Referenced by ExecuteLoop().

◆ ExecuteStreak()

void Nektar::VortexWaveInteraction::ExecuteStreak ( void  )

Definition at line 494 of file VortexWaveInteraction.cpp.

495 {
496  // Create driver
497 #if 1
498  std::string vDriverModule;
499  m_sessionStreak->LoadSolverInfo("Driver", vDriverModule, "Standard");
500 
502  vDriverModule, m_sessionStreak, m_graphStreak);
503  solverStreak->Execute();
504 
505  m_streakField = solverStreak->GetEqu()[0]->UpdateFields();
506 #else
507  // Setup and execute Advection Diffusion solver
508  string vEquation = m_sessionStreak->GetSolverInfo("EqType");
509  EquationSystemSharedPtr solverStreak =
511  m_graphStreak);
512 
513  cout << "Executing Streak Solver" << endl;
514  solverStreak->DoInitialise();
515  solverStreak->DoSolve();
516  solverStreak->Output();
517 
518  m_streakField = solverStreak->UpdateFields();
519 
520  if (m_sessionVWI->DefinesSolverInfo("INTERFACE"))
521  {
522  for (int g = 0; g < solverStreak->GetNvariables(); ++g)
523  {
524  NekDouble vL2Error = solverStreak->L2Error(g, false);
525  NekDouble vLinfError = solverStreak->LinfError(g);
526  cout << "L 2 error (variable " << solverStreak->GetVariable(g)
527  << ") : " << vL2Error << endl;
528  cout << "L inf error (variable " << solverStreak->GetVariable(g)
529  << ") : " << vLinfError << endl;
530  }
531  }
532 #endif
533 
534  cout << "Executing cp -f session.fld session_streak.fld" << endl;
535  CopyFile(".fld", "_streak.fld");
536 }
std::shared_ptr< Driver > DriverSharedPtr
A shared pointer to a Driver object.
Definition: Driver.h:51
DriverFactory & GetDriverFactory()
Definition: Driver.cpp:64

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

Referenced by ExecuteLoop().

◆ ExecuteWave()

void Nektar::VortexWaveInteraction::ExecuteWave ( void  )

Do linearised NavierStokes Session with Modified Arnoldi

Definition at line 538 of file VortexWaveInteraction.cpp.

539 {
540 
541  // Set the initial beta value in stability to be equal to VWI file
542  std::string LZstr("LZ");
543  NekDouble LZ = 2 * M_PI / m_alpha[0];
544  cout << "Setting LZ in Linearised solver to " << LZ << endl;
545  m_sessionWave->SetParameter(LZstr, LZ);
546 
547  // Create driver
548  std::string vDriverModule;
549  m_sessionWave->LoadSolverInfo("Driver", vDriverModule, "ModifiedArnoldi");
550  cout << "Setting up linearised NS Solver" << endl;
552  vDriverModule, m_sessionWave, m_graphWave);
553 
554  /// Do linearised NavierStokes Session with Modified Arnoldi
555  cout << "Executing wave solution " << endl;
556  solverWave->Execute();
557 
558  // Copy file to a rst location for next restart
559  cout << "Executing cp -f session_eig_0 session_eig_0.rst" << endl;
560  CopyFile("_eig_0", "_eig_0.rst");
561 
562  // Store data relevant to other operations
563  m_leading_real_evl[0] = solverWave->GetRealEvl()[0];
564  m_leading_imag_evl[0] = solverWave->GetImagEvl()[0];
565 
566  // note this will only be true for modified Arnoldi
567  NekDouble realShift = 0.0;
568  if (m_sessionWave->DefinesParameter("RealShift"))
569  {
570  bool defineshift;
571  // only use shift in modifiedArnoldi solver since
572  // implicitly handled in Arpack.
573  m_sessionWave->MatchSolverInfo("Driver", "ModifiedArnoldi", defineshift,
574  true);
575  if (defineshift)
576  {
577  realShift = m_sessionWave->GetParameter("RealShift");
578  }
579  }
580 
581  // Set up leading eigenvalue from inverse
582  NekDouble invmag = 1.0 / (m_leading_real_evl[0] * m_leading_real_evl[0] +
584  m_leading_real_evl[0] *= -invmag;
585  m_leading_real_evl[0] += realShift;
586  m_leading_imag_evl[0] *= invmag;
587 
588  m_waveVelocities = solverWave->GetEqu()[0]->UpdateFields();
589  m_wavePressure = solverWave->GetEqu()[0]->GetPressure();
590 
591  if (m_sessionVWI->DefinesSolverInfo("INTERFACE"))
592  {
593  cout << "Growth =" << m_leading_real_evl[0] << endl;
594  cout << "Phase =" << m_leading_imag_evl[0] << endl;
595  }
596 }

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

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

◆ FileRelaxation()

void Nektar::VortexWaveInteraction::FileRelaxation ( int  reg)

Definition at line 2132 of file VortexWaveInteraction.cpp.

2133 {
2134  cout << "relaxation..." << endl;
2135  static int cnt = 0;
2136  Array<OneD, MultiRegions::ExpListSharedPtr> Iexp =
2137  m_rollField[0]->GetBndCondExpansions();
2138  // cast to 1D explist (otherwise appenddata doesn't work)
2141  *std::static_pointer_cast<MultiRegions::ExpList>(Iexp[reg]));
2142  int nq = Ilayer->GetTotPoints();
2143  if (cnt == 0)
2144  {
2145  m_bcsForcing = Array<OneD, Array<OneD, NekDouble>>(4);
2146  m_bcsForcing[0] = Array<OneD, NekDouble>(4 * nq);
2147  for (int i = 1; i < 4; ++i)
2148  {
2149  m_bcsForcing[i] = m_bcsForcing[i - 1] + nq;
2150  }
2151  }
2152 
2153  // Read in mesh from input file
2154 
2155  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef_u;
2156  std::vector<std::vector<NekDouble>> FieldData_u;
2157  string file = m_sessionName;
2158 
2159  file += "_u_5.bc";
2162  fld->Import(file, FieldDef_u, FieldData_u);
2163  Ilayer->ExtractDataToCoeffs(FieldDef_u[0], FieldData_u[0],
2164  FieldDef_u[0]->m_fields[0],
2165  Ilayer->UpdateCoeffs());
2166  Ilayer->BwdTrans(Ilayer->GetCoeffs(), Ilayer->UpdatePhys());
2167 
2168  if (cnt == 0)
2169  {
2170  Vmath::Vcopy(nq, Ilayer->UpdatePhys(), 1, m_bcsForcing[2], 1);
2171  }
2172  Vmath::Vcopy(nq, Ilayer->UpdatePhys(), 1, m_bcsForcing[0], 1);
2173 
2174  // case relaxation==0 an additional array is needed
2175  Array<OneD, NekDouble> tmp_forcing(nq, 0.0);
2176 
2177  if (cnt != 0)
2178  {
2179  if (m_vwiRelaxation == 1.0)
2180  {
2181  Vmath::Vcopy(nq, m_bcsForcing[0], 1, tmp_forcing, 1);
2182  }
2183  Vmath::Smul(nq, 1.0 - m_vwiRelaxation, m_bcsForcing[0], 1,
2184  m_bcsForcing[0], 1);
2186  1, Ilayer->UpdatePhys(), 1);
2187  // generate again the bcs files:
2188 
2189  Array<OneD, Array<OneD, NekDouble>> fieldcoeffs(1);
2190  Ilayer->FwdTransLocalElmt(Ilayer->GetPhys(), Ilayer->UpdateCoeffs());
2191  fieldcoeffs[0] = Ilayer->UpdateCoeffs();
2192  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef1 =
2193  Ilayer->GetFieldDefinitions();
2194  std::vector<std::vector<NekDouble>> FieldData_1(FieldDef1.size());
2195  ;
2196  FieldDef1[0]->m_fields.push_back("u");
2197  Ilayer->AppendFieldData(FieldDef1[0], FieldData_1[0]);
2198  fld->Write(file, FieldDef1, FieldData_1);
2199  // save the bcs for the next iteration
2200  if (m_vwiRelaxation != 1.0)
2201  {
2202  Vmath::Smul(nq, 1. / (1.0 - m_vwiRelaxation), m_bcsForcing[0], 1,
2203  m_bcsForcing[0], 1);
2204  Vmath::Vcopy(nq, m_bcsForcing[0], 1, m_bcsForcing[2], 1);
2205  }
2206  else
2207  {
2208  Vmath::Vcopy(nq, tmp_forcing, 1, m_bcsForcing[2], 1);
2209  }
2210  }
2211 
2212  file = m_sessionName + "_v_5.bc";
2213 
2214  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef_v;
2215  std::vector<std::vector<NekDouble>> FieldData_v;
2216  fld->Import(file, FieldDef_v, FieldData_v);
2217  Ilayer->ExtractDataToCoeffs(FieldDef_v[0], FieldData_v[0],
2218  FieldDef_v[0]->m_fields[0],
2219  Ilayer->UpdateCoeffs());
2220  Ilayer->BwdTrans(Ilayer->GetCoeffs(), Ilayer->UpdatePhys());
2221  if (cnt == 0)
2222  {
2223  Vmath::Vcopy(nq, Ilayer->UpdatePhys(), 1, m_bcsForcing[3], 1);
2224  }
2225  Vmath::Vcopy(nq, Ilayer->UpdatePhys(), 1, m_bcsForcing[1], 1);
2226  if (cnt != 0)
2227  {
2228  if (m_vwiRelaxation == 1.0)
2229  {
2230  Vmath::Vcopy(nq, m_bcsForcing[1], 1, tmp_forcing, 1);
2231  }
2232  Vmath::Smul(nq, 1.0 - m_vwiRelaxation, m_bcsForcing[1], 1,
2233  m_bcsForcing[1], 1);
2235  1, Ilayer->UpdatePhys(), 1);
2236  // generate again the bcs files:
2237  Array<OneD, Array<OneD, NekDouble>> fieldcoeffs(1);
2238  Ilayer->FwdTransLocalElmt(Ilayer->GetPhys(), Ilayer->UpdateCoeffs());
2239  fieldcoeffs[0] = Ilayer->UpdateCoeffs();
2240  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef2 =
2241  Ilayer->GetFieldDefinitions();
2242  std::vector<std::vector<NekDouble>> FieldData_2(FieldDef2.size());
2243  ;
2244  FieldDef2[0]->m_fields.push_back("v");
2245  Ilayer->AppendFieldData(FieldDef2[0], FieldData_2[0]);
2246  fld->Write(file, FieldDef2, FieldData_2);
2247  // save the bcs for the next iteration
2248  if (m_vwiRelaxation != 1.0)
2249  {
2250  Vmath::Smul(nq, 1. / (1.0 - m_vwiRelaxation), m_bcsForcing[1], 1,
2251  m_bcsForcing[1], 1);
2252  Vmath::Vcopy(nq, m_bcsForcing[1], 1, m_bcsForcing[3], 1);
2253  }
2254  else
2255  {
2256  Vmath::Vcopy(nq, tmp_forcing, 1, m_bcsForcing[3], 1);
2257  }
2258  }
2259 
2260  cnt++;
2261 }
static std::shared_ptr< FieldIO > CreateDefault(const LibUtilities::SessionReaderSharedPtr session)
Returns an object for the default FieldIO method.
Definition: FieldIO.cpp:198
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Array< OneD, Array< OneD, NekDouble > > m_bcsForcing
std::shared_ptr< FieldIO > FieldIOSharedPtr
Definition: FieldIO.h:301
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::LibUtilities::FieldIO::CreateDefault(), m_bcsForcing, m_rollField, m_sessionName, m_sessionVWI, m_vwiRelaxation, Vmath::Smul(), Vmath::Svtvp(), and Vmath::Vcopy().

◆ GetAlpha()

NekDouble Nektar::VortexWaveInteraction::GetAlpha ( void  )
inline

Definition at line 121 of file VortexWaveInteraction.h.

122  {
123  return m_alpha[0];
124  }

Referenced by DoFixedForcingIteration(), and main().

◆ GetAlphaStep()

NekDouble Nektar::VortexWaveInteraction::GetAlphaStep ( void  )
inline

Definition at line 126 of file VortexWaveInteraction.h.

127  {
128  return m_alphaStep;
129  }

Referenced by main().

◆ GetDAlphaDWaveForceMag()

NekDouble Nektar::VortexWaveInteraction::GetDAlphaDWaveForceMag ( void  )
inline

Definition at line 141 of file VortexWaveInteraction.h.

142  {
143  return m_dAlphaDWaveForceMag;
144  }

◆ GetEigRelTol()

NekDouble Nektar::VortexWaveInteraction::GetEigRelTol ( void  )
inline

Definition at line 151 of file VortexWaveInteraction.h.

153  {
154  return m_eigRelTol;
155  }

Referenced by DoFixedForcingIteration().

◆ GetIterEnd()

int Nektar::VortexWaveInteraction::GetIterEnd ( )
inline

Definition at line 101 of file VortexWaveInteraction.h.

102  {
103  return m_iterEnd;
104  }

Referenced by DoFixedForcingIteration().

◆ GetIterStart()

int Nektar::VortexWaveInteraction::GetIterStart ( )
inline

Definition at line 96 of file VortexWaveInteraction.h.

97  {
98  return m_iterStart;
99  }

Referenced by DoFixedForcingIteration().

◆ GetMaxOuterIterations()

int Nektar::VortexWaveInteraction::GetMaxOuterIterations ( void  )
inline

Definition at line 116 of file VortexWaveInteraction.h.

117  {
118  return m_maxOuterIterations;
119  }

Referenced by DoFixedForcingIteration().

◆ GetMaxWaveForceMagIter()

int Nektar::VortexWaveInteraction::GetMaxWaveForceMagIter ( void  )
inline

Definition at line 146 of file VortexWaveInteraction.h.

147  {
148  return m_maxWaveForceMagIter;
149  }

Referenced by main().

◆ GetMinInnerIterations()

int Nektar::VortexWaveInteraction::GetMinInnerIterations ( void  )
inline

Definition at line 157 of file VortexWaveInteraction.h.

158  {
159  return m_minInnerIterations;
160  }

Referenced by DoFixedForcingIteration().

◆ GetNOuterIterations()

int Nektar::VortexWaveInteraction::GetNOuterIterations ( void  )
inline

Definition at line 111 of file VortexWaveInteraction.h.

112  {
113  return m_nOuterIterations;
114  }

Referenced by DoFixedForcingIteration().

◆ GetPrevAlpha()

NekDouble Nektar::VortexWaveInteraction::GetPrevAlpha ( void  )
inline

Definition at line 162 of file VortexWaveInteraction.h.

163  {
164  return m_prevAlpha;
165  }

Referenced by DoFixedForcingIteration().

◆ GetReflectionIndex()

Array< OneD, int > Nektar::VortexWaveInteraction::GetReflectionIndex ( void  )

Definition at line 2047 of file VortexWaveInteraction.cpp.

2048 {
2049  int i, j;
2050  int npts = m_waveVelocities[0]->GetPlane(0)->GetNpoints();
2051  int nel = m_waveVelocities[0]->GetNumElmts();
2052  Array<OneD, int> index(npts);
2053 
2054  Array<OneD, NekDouble> coord(2);
2055  Array<OneD, NekDouble> coord_x(npts);
2056  Array<OneD, NekDouble> coord_y(npts);
2057 
2058  //-> Dermine the point which is on coordinate (x -> -x + Lx/2, y-> -y)
2059  m_waveVelocities[0]->GetPlane(0)->GetCoords(coord_x, coord_y);
2060  NekDouble xmax = Vmath::Vmax(npts, coord_x, 1);
2061  // NekDouble tol =
2062  // NekConstants::kGeomFactorsTol*NekConstants::kGeomFactorsTol;
2063  NekDouble tol = 1e-5;
2064  NekDouble xnew, ynew;
2065 
2066  int start = npts - 1;
2067  int e_npts;
2068 
2069  bool useOnlyQuads = false;
2070  if (m_sessionVWI->DefinesSolverInfo("SymmetriseOnlyQuads"))
2071  {
2072  useOnlyQuads = true;
2073  }
2074 
2075  int cnt;
2076  for (int e = 0; e < nel; ++e)
2077  {
2078  e_npts = m_waveVelocities[0]->GetExp(e)->GetTotPoints();
2079  cnt = m_waveVelocities[0]->GetPhys_Offset(e);
2080 
2081  if (useOnlyQuads)
2082  {
2083  if (m_waveVelocities[0]->GetExp(e)->DetShapeType() ==
2085  {
2086  for (i = 0; i < e_npts; ++i)
2087  {
2088  index[cnt + i] = -1;
2089  }
2090  continue;
2091  }
2092  }
2093 
2094  for (i = cnt; i < cnt + e_npts; ++i)
2095  {
2096  xnew = -coord_x[i] + xmax;
2097  ynew = -coord_y[i];
2098 
2099  for (j = start; j >= 0; --j)
2100  {
2101  if ((coord_x[j] - xnew) * (coord_x[j] - xnew) +
2102  (coord_y[j] - ynew) * (coord_y[j] - ynew) <
2103  tol)
2104  {
2105  index[i] = j;
2106  start = j;
2107  break;
2108  }
2109  }
2110 
2111  if (j == -1)
2112  {
2113 
2114  for (j = npts - 1; j > start; --j)
2115  {
2116 
2117  if ((coord_x[j] - xnew) * (coord_x[j] - xnew) +
2118  (coord_y[j] - ynew) * (coord_y[j] - ynew) <
2119  tol)
2120  {
2121  index[i] = j;
2122  break;
2123  }
2124  }
2125  ASSERTL0(j != start, "Failed to find matching point");
2126  }
2127  }
2128  }
2129  return index;
2130 }

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

Referenced by CalcNonLinearWaveForce().

◆ GetVWIIterationType()

VWIIterationType Nektar::VortexWaveInteraction::GetVWIIterationType ( void  )
inline

Definition at line 106 of file VortexWaveInteraction.h.

107  {
108  return m_VWIIterationType;
109  }

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

◆ GetWaveForceMag()

NekDouble Nektar::VortexWaveInteraction::GetWaveForceMag ( void  )
inline

Definition at line 131 of file VortexWaveInteraction.h.

132  {
133  return m_waveForceMag[0];
134  }

Referenced by main().

◆ GetWaveForceMagStep()

NekDouble Nektar::VortexWaveInteraction::GetWaveForceMagStep ( void  )
inline

Definition at line 136 of file VortexWaveInteraction.h.

137  {
138  return m_waveForceMagStep;
139  }

Referenced by main().

◆ IfIterInterface()

bool Nektar::VortexWaveInteraction::IfIterInterface ( void  )
inline

Definition at line 197 of file VortexWaveInteraction.h.

198  {
199  return m_iterinterface;
200  }

Referenced by DoFixedForcingIteration().

◆ MoveFile()

void Nektar::VortexWaveInteraction::MoveFile ( std::string  fileend,
std::string  dir,
int  n 
)

Definition at line 1132 of file VortexWaveInteraction.cpp.

1133 {
1134  static map<string, int> opendir;
1135 
1136  if (opendir.find(dir) == opendir.end())
1137  {
1138  // make directory and presume will fail if it already exists
1139  string mkdir = "mkdir " + dir;
1140  ASSERTL0(system(mkdir.c_str()) == 0,
1141  "Failed to make directory '" + dir + "'");
1142  opendir[dir] = 1;
1143  }
1144 
1145  string savefile =
1146  dir + "/" + file + "." + boost::lexical_cast<std::string>(n);
1147  string syscall = "mv -f " + file + " " + savefile;
1148 
1149  ASSERTL0(system(syscall.c_str()) == 0, syscall.c_str());
1150 }

References ASSERTL0, and CG_Iterations::savefile.

Referenced by DoFixedForcingIteration().

◆ SaveFile()

void Nektar::VortexWaveInteraction::SaveFile ( std::string  fileend,
std::string  dir,
int  n 
)

Definition at line 1111 of file VortexWaveInteraction.cpp.

1112 {
1113  static map<string, int> opendir;
1114 
1115  if (opendir.find(dir) == opendir.end())
1116  {
1117  // make directory and presume will fail if it already exists
1118  string mkdir = "mkdir " + dir;
1119  ASSERTL0(system(mkdir.c_str()) == 0,
1120  "Failed to make directory '" + dir + "'");
1121 
1122  opendir[dir] = 1;
1123  }
1124 
1125  string savefile =
1126  dir + "/" + file + "." + boost::lexical_cast<std::string>(n);
1127  string syscall = "cp -f " + file + " " + savefile;
1128 
1129  ASSERTL0(system(syscall.c_str()) == 0, syscall.c_str());
1130 }

References ASSERTL0, and CG_Iterations::savefile.

Referenced by SaveLoopDetails().

◆ SaveLoopDetails()

void Nektar::VortexWaveInteraction::SaveLoopDetails ( std::string  dir,
int  i 
)

Definition at line 1182 of file VortexWaveInteraction.cpp.

1184 {
1185  // Save NS restart file
1186  SaveFile(m_sessionName + ".rst", SaveDir, i + 1);
1187  // Save Streak Solution
1188  SaveFile(m_sessionName + "_streak.fld", SaveDir, i);
1189  // Save Wave solution output
1190  SaveFile(m_sessionName + ".evl", SaveDir, i);
1191  SaveFile(m_sessionName + "_eig_0", SaveDir, i);
1192  // Save new field file of eigenvalue
1193  SaveFile(m_sessionName + ".fld", SaveDir, i);
1194  if (!(m_sessionVWI->DefinesSolverInfo("INTERFACE")))
1195  {
1196  // Save new forcing file
1197  SaveFile(m_sessionName + ".vwi", SaveDir, i + 1);
1198  }
1199 }
void SaveFile(std::string fileend, std::string dir, int n)

References m_sessionName, m_sessionVWI, and SaveFile().

Referenced by DoFixedForcingIteration().

◆ SetAlpha()

void Nektar::VortexWaveInteraction::SetAlpha ( NekDouble  alpha)
inline

Definition at line 167 of file VortexWaveInteraction.h.

168  {
169  m_alpha[0] = alpha;
170  }

Referenced by main().

◆ SetAlphaStep()

void Nektar::VortexWaveInteraction::SetAlphaStep ( NekDouble  step)
inline

Definition at line 182 of file VortexWaveInteraction.h.

183  {
184  m_alphaStep = step;
185  }

◆ SetEigRelTol()

void Nektar::VortexWaveInteraction::SetEigRelTol ( NekDouble  tol)
inline

Definition at line 177 of file VortexWaveInteraction.h.

178  {
179  m_eigRelTol = tol;
180  }

Referenced by DoFixedForcingIteration().

◆ SetMinInnerIterations()

void Nektar::VortexWaveInteraction::SetMinInnerIterations ( int  niter)
inline

Definition at line 187 of file VortexWaveInteraction.h.

188  {
189  m_minInnerIterations = niter;
190  }

Referenced by DoFixedForcingIteration().

◆ SetPrevAlpha()

void Nektar::VortexWaveInteraction::SetPrevAlpha ( NekDouble  alpha)
inline

Definition at line 192 of file VortexWaveInteraction.h.

193  {
194  m_prevAlpha = alpha;
195  }

Referenced by main().

◆ SetWaveForceMag()

void Nektar::VortexWaveInteraction::SetWaveForceMag ( NekDouble  mag)
inline

Definition at line 172 of file VortexWaveInteraction.h.

173  {
174  m_waveForceMag[0] = mag;
175  }

Referenced by main().

◆ UpdateAlpha()

void Nektar::VortexWaveInteraction::UpdateAlpha ( int  n)

Definition at line 1932 of file VortexWaveInteraction.cpp.

1933 {
1934  NekDouble alp_new = 0.0;
1935 
1936  if (outeriter == 1)
1937  {
1938  m_alpha[1] = m_alpha[0];
1939  if (m_leading_real_evl[0] > 0.0)
1940  {
1941  alp_new = m_alpha[0] + m_alphaStep;
1942  }
1943  else
1944  {
1945  alp_new = m_alpha[0] - m_alphaStep;
1946  }
1947  }
1948  else
1949  {
1950  int i;
1951  int nstore = (m_alpha.size() < outeriter) ? m_alpha.size() : outeriter;
1952  Array<OneD, NekDouble> Alpha(nstore);
1953  Array<OneD, NekDouble> Growth(nstore);
1954 
1955  Vmath::Vcopy(nstore, m_alpha, 1, Alpha, 1);
1956  Vmath::Vcopy(nstore, m_leading_real_evl, 1, Growth, 1);
1957 
1958  // Sort Alpha Growth values;
1959  NekDouble store;
1960  int k;
1961  for (i = 0; i < nstore; ++i)
1962  {
1963  k = Vmath::Imin(nstore - i, &Alpha[i], 1);
1964 
1965  store = Alpha[i];
1966  Alpha[i] = Alpha[i + k];
1967  Alpha[i + k] = store;
1968 
1969  store = Growth[i];
1970  Growth[i] = Growth[i + k];
1971  Growth[i + k] = store;
1972  }
1973 
1974  // See if we have any values that cross zero
1975  for (i = 0; i < nstore - 1; ++i)
1976  {
1977  if (Growth[i] * Growth[i + 1] < 0.0)
1978  {
1979  break;
1980  }
1981  }
1982 
1983  if (i != nstore - 1)
1984  {
1985  if (nstore == 2)
1986  {
1987  alp_new = (Alpha[0] * Growth[1] - Alpha[1] * Growth[0]) /
1988  (Growth[1] - Growth[0]);
1989  }
1990  else
1991  {
1992  // use a quadratic fit and step through 10000 points
1993  // to find zero.
1994  int j;
1995  int nsteps = 10000;
1996  int idx = (i == 0) ? 1 : i;
1997  NekDouble da = Alpha[idx + 1] - Alpha[idx - 1];
1998  NekDouble gval_m1 = Growth[idx - 1], a, gval;
1999  NekDouble c1 = Growth[idx - 1] / (Alpha[idx - 1] - Alpha[idx]) /
2000  (Alpha[idx - 1] - Alpha[idx + 1]);
2001  NekDouble c2 = Growth[idx] / (Alpha[idx] - Alpha[idx - 1]) /
2002  (Alpha[idx] - Alpha[idx + 1]);
2003  NekDouble c3 = Growth[idx + 1] /
2004  (Alpha[idx + 1] - Alpha[idx - 1]) /
2005  (Alpha[idx + 1] - Alpha[idx]);
2006 
2007  for (j = 1; j < nsteps + 1; ++j)
2008  {
2009  a = Alpha[i] + j * da / (NekDouble)nsteps;
2010  gval = c1 * (a - Alpha[idx]) * (a - Alpha[idx + 1]) +
2011  c2 * (a - Alpha[idx - 1]) * (a - Alpha[idx + 1]) +
2012  c3 * (a - Alpha[idx - 1]) * (a - Alpha[idx]);
2013 
2014  if (gval * gval_m1 < 0.0)
2015  {
2016  alp_new = ((a + da / (NekDouble)nsteps) * gval -
2017  a * gval_m1) /
2018  (gval - gval_m1);
2019  break;
2020  }
2021  }
2022  }
2023  }
2024  else // step backward/forward
2025  {
2026  if (Growth[i] > 0.0)
2027  {
2028  alp_new = m_alpha[0] + m_alphaStep;
2029  }
2030  else
2031  {
2032  alp_new = m_alpha[0] - m_alphaStep;
2033  }
2034  }
2035  }
2036 
2037  for (int i = m_alpha.size() - 1; i > 0; --i)
2038  {
2039  m_alpha[i] = m_alpha[i - 1];
2042  }
2043 
2044  m_alpha[0] = alp_new;
2045 }
int Imin(int n, const T *x, const int incx)
Return the index of the minimum element in x.
Definition: Vmath.cpp:1023

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

Referenced by DoFixedForcingIteration(), and VortexWaveInteraction().

◆ UpdateDAlphaDWaveForceMag()

void Nektar::VortexWaveInteraction::UpdateDAlphaDWaveForceMag ( NekDouble  alphainit)

Definition at line 1927 of file VortexWaveInteraction.cpp.

1928 {
1929  m_dAlphaDWaveForceMag = (m_alpha[0] - alpha_init) / m_waveForceMagStep;
1930 }

References m_alpha, m_dAlphaDWaveForceMag, and m_waveForceMagStep.

Referenced by DoFixedForcingIteration().

◆ UpdateWaveForceMag()

void Nektar::VortexWaveInteraction::UpdateWaveForceMag ( int  n)

Definition at line 1806 of file VortexWaveInteraction.cpp.

1807 {
1808  NekDouble wavef_new = 0.0;
1809 
1810  if (outeriter == 1)
1811  {
1813  if (m_leading_real_evl[0] > 0.0)
1814  {
1815  wavef_new = m_waveForceMag[0] - m_waveForceMagStep;
1816  }
1817  else
1818  {
1819  wavef_new = m_waveForceMag[0] + m_waveForceMagStep;
1820  }
1821  }
1822  else
1823  {
1824  int i;
1825  int nstore = (m_waveForceMag.size() < outeriter) ? m_waveForceMag.size()
1826  : outeriter;
1827  Array<OneD, NekDouble> WaveForce(nstore);
1828  Array<OneD, NekDouble> Growth(nstore);
1829 
1830  Vmath::Vcopy(nstore, m_waveForceMag, 1, WaveForce, 1);
1831  Vmath::Vcopy(nstore, m_leading_real_evl, 1, Growth, 1);
1832 
1833  // Sort WaveForce Growth values;
1834  NekDouble store;
1835  int k;
1836  for (i = 0; i < nstore; ++i)
1837  {
1838  k = Vmath::Imin(nstore - i, &WaveForce[i], 1);
1839 
1840  store = WaveForce[i];
1841  WaveForce[i] = WaveForce[i + k];
1842  WaveForce[i + k] = store;
1843 
1844  store = Growth[i];
1845  Growth[i] = Growth[i + k];
1846  Growth[i + k] = store;
1847  }
1848 
1849  // See if we have any values that cross zero
1850  for (i = 0; i < nstore - 1; ++i)
1851  {
1852  if (Growth[i] * Growth[i + 1] < 0.0)
1853  {
1854  break;
1855  }
1856  }
1857 
1858  if (i != nstore - 1)
1859  {
1860  if (nstore == 2)
1861  {
1862  wavef_new =
1863  (WaveForce[0] * Growth[1] - WaveForce[1] * Growth[0]) /
1864  (Growth[1] - Growth[0]);
1865  }
1866  else
1867  {
1868  // use a quadratic fit and step through 10000 points
1869  // to find zero.
1870  int j;
1871  int nsteps = 10000;
1872  int idx = (i == 0) ? 1 : i;
1873  NekDouble da = WaveForce[idx + 1] - WaveForce[idx - 1];
1874  NekDouble gval_m1 = Growth[idx - 1], a, gval;
1875  NekDouble c1 = Growth[idx - 1] /
1876  (WaveForce[idx - 1] - WaveForce[idx]) /
1877  (WaveForce[idx - 1] - WaveForce[idx + 1]);
1878  NekDouble c2 = Growth[idx] /
1879  (WaveForce[idx] - WaveForce[idx - 1]) /
1880  (WaveForce[idx] - WaveForce[idx + 1]);
1881  NekDouble c3 = Growth[idx + 1] /
1882  (WaveForce[idx + 1] - WaveForce[idx - 1]) /
1883  (WaveForce[idx + 1] - WaveForce[idx]);
1884 
1885  for (j = 1; j < nsteps + 1; ++j)
1886  {
1887  a = WaveForce[i] + j * da / (NekDouble)nsteps;
1888  gval =
1889  c1 * (a - WaveForce[idx]) * (a - WaveForce[idx + 1]) +
1890  c2 * (a - WaveForce[idx - 1]) *
1891  (a - WaveForce[idx + 1]) +
1892  c3 * (a - WaveForce[idx - 1]) * (a - WaveForce[idx]);
1893 
1894  if (gval * gval_m1 < 0.0)
1895  {
1896  wavef_new = ((a + da / (NekDouble)nsteps) * gval -
1897  a * gval_m1) /
1898  (gval - gval_m1);
1899  break;
1900  }
1901  }
1902  }
1903  }
1904  else // step backward/forward
1905  {
1906  if (Growth[i] > 0.0)
1907  {
1908  wavef_new = m_waveForceMag[0] - m_waveForceMagStep;
1909  }
1910  else
1911  {
1912  wavef_new = m_waveForceMag[0] + m_waveForceMagStep;
1913  }
1914  }
1915  }
1916 
1917  for (int i = m_waveForceMag.size() - 1; i > 0; --i)
1918  {
1919  m_waveForceMag[i] = m_waveForceMag[i - 1];
1922  }
1923 
1924  m_waveForceMag[0] = wavef_new;
1925 }

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

Referenced by DoFixedForcingIteration().

Member Data Documentation

◆ m_alpha

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

◆ m_alphaStep

NekDouble Nektar::VortexWaveInteraction::m_alphaStep
private

Definition at line 235 of file VortexWaveInteraction.h.

Referenced by UpdateAlpha(), and VortexWaveInteraction().

◆ m_bcsForcing

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

Definition at line 252 of file VortexWaveInteraction.h.

Referenced by FileRelaxation().

◆ m_dAlphaDWaveForceMag

NekDouble Nektar::VortexWaveInteraction::m_dAlphaDWaveForceMag
private

Definition at line 239 of file VortexWaveInteraction.h.

Referenced by UpdateDAlphaDWaveForceMag(), and VortexWaveInteraction().

◆ m_deltaFcnApprox

bool Nektar::VortexWaveInteraction::m_deltaFcnApprox
private

Definition at line 217 of file VortexWaveInteraction.h.

Referenced by CalcNonLinearWaveForce(), and VortexWaveInteraction().

◆ m_deltaFcnDecay

NekDouble Nektar::VortexWaveInteraction::m_deltaFcnDecay
private

Definition at line 222 of file VortexWaveInteraction.h.

Referenced by CalcNonLinearWaveForce(), and VortexWaveInteraction().

◆ m_eigRelTol

NekDouble Nektar::VortexWaveInteraction::m_eigRelTol
private

Definition at line 237 of file VortexWaveInteraction.h.

Referenced by CheckEigIsStationary(), and VortexWaveInteraction().

◆ m_graphRoll

SpatialDomains::MeshGraphSharedPtr Nektar::VortexWaveInteraction::m_graphRoll
private

Definition at line 262 of file VortexWaveInteraction.h.

Referenced by ExecuteRoll(), and VortexWaveInteraction().

◆ m_graphStreak

SpatialDomains::MeshGraphSharedPtr Nektar::VortexWaveInteraction::m_graphStreak
private

Definition at line 266 of file VortexWaveInteraction.h.

Referenced by ExecuteStreak(), and VortexWaveInteraction().

◆ m_graphWave

SpatialDomains::MeshGraphSharedPtr Nektar::VortexWaveInteraction::m_graphWave
private

Definition at line 268 of file VortexWaveInteraction.h.

Referenced by ExecuteWave(), and VortexWaveInteraction().

◆ m_iterEnd

int Nektar::VortexWaveInteraction::m_iterEnd
private

Definition at line 209 of file VortexWaveInteraction.h.

Referenced by VortexWaveInteraction().

◆ m_iterinterface

bool Nektar::VortexWaveInteraction::m_iterinterface
private

Definition at line 242 of file VortexWaveInteraction.h.

Referenced by VortexWaveInteraction().

◆ m_iterStart

int Nektar::VortexWaveInteraction::m_iterStart
private

Definition at line 208 of file VortexWaveInteraction.h.

Referenced by VortexWaveInteraction().

◆ m_leading_imag_evl

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

◆ m_leading_real_evl

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

◆ m_maxOuterIterations

int Nektar::VortexWaveInteraction::m_maxOuterIterations
private

Definition at line 212 of file VortexWaveInteraction.h.

Referenced by VortexWaveInteraction().

◆ m_maxWaveForceMagIter

int Nektar::VortexWaveInteraction::m_maxWaveForceMagIter
private

Definition at line 215 of file VortexWaveInteraction.h.

Referenced by VortexWaveInteraction().

◆ m_minInnerIterations

int Nektar::VortexWaveInteraction::m_minInnerIterations
private

Definition at line 213 of file VortexWaveInteraction.h.

Referenced by CheckEigIsStationary(), and VortexWaveInteraction().

◆ m_moveMeshToCriticalLayer

bool Nektar::VortexWaveInteraction::m_moveMeshToCriticalLayer
private

Definition at line 220 of file VortexWaveInteraction.h.

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

◆ m_neutralPointTol

NekDouble Nektar::VortexWaveInteraction::m_neutralPointTol
private

Definition at line 236 of file VortexWaveInteraction.h.

Referenced by CheckIfAtNeutralPoint(), and VortexWaveInteraction().

◆ m_nOuterIterations

int Nektar::VortexWaveInteraction::m_nOuterIterations
private

Definition at line 211 of file VortexWaveInteraction.h.

Referenced by VortexWaveInteraction().

◆ m_prevAlpha

NekDouble Nektar::VortexWaveInteraction::m_prevAlpha
private

Definition at line 240 of file VortexWaveInteraction.h.

◆ m_rollField

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

Definition at line 256 of file VortexWaveInteraction.h.

Referenced by ExecuteRoll(), and FileRelaxation().

◆ m_rollForceScale

NekDouble Nektar::VortexWaveInteraction::m_rollForceScale
private

Definition at line 227 of file VortexWaveInteraction.h.

Referenced by ExecuteRoll(), and VortexWaveInteraction().

◆ m_sessionName

std::string Nektar::VortexWaveInteraction::m_sessionName
private

◆ m_sessionRoll

LibUtilities::SessionReaderSharedPtr Nektar::VortexWaveInteraction::m_sessionRoll
private

◆ m_sessionStreak

LibUtilities::SessionReaderSharedPtr Nektar::VortexWaveInteraction::m_sessionStreak
private

Definition at line 265 of file VortexWaveInteraction.h.

Referenced by ExecuteStreak(), and VortexWaveInteraction().

◆ m_sessionVWI

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

◆ m_sessionWave

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

Definition at line 267 of file VortexWaveInteraction.h.

Referenced by ExecuteWave(), and VortexWaveInteraction().

◆ m_solverRoll

EquationSystemSharedPtr Nektar::VortexWaveInteraction::m_solverRoll
private

◆ m_streakField

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

Definition at line 254 of file VortexWaveInteraction.h.

Referenced by CalcNonLinearWaveForce(), and ExecuteStreak().

◆ m_useLinfPressureNorm

bool Nektar::VortexWaveInteraction::m_useLinfPressureNorm
private

Definition at line 218 of file VortexWaveInteraction.h.

Referenced by CalcNonLinearWaveForce(), and VortexWaveInteraction().

◆ m_vwiForcing

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

◆ m_vwiForcingObj

SolverUtils::ForcingProgrammaticSharedPtr Nektar::VortexWaveInteraction::m_vwiForcingObj
private

Definition at line 250 of file VortexWaveInteraction.h.

Referenced by ExecuteRoll().

◆ m_VWIIterationType

VWIIterationType Nektar::VortexWaveInteraction::m_VWIIterationType
private

Definition at line 244 of file VortexWaveInteraction.h.

Referenced by VortexWaveInteraction().

◆ m_vwiRelaxation

NekDouble Nektar::VortexWaveInteraction::m_vwiRelaxation
private

◆ m_waveForceMag

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

◆ m_waveForceMagStep

NekDouble Nektar::VortexWaveInteraction::m_waveForceMagStep
private

◆ m_wavePressure

MultiRegions::ExpListSharedPtr Nektar::VortexWaveInteraction::m_wavePressure
private

◆ m_waveVelocities

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