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:111
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:303

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 1162 of file VortexWaveInteraction.cpp.

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

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 1171 of file VortexWaveInteraction.cpp.

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

References m_alpha, m_leading_imag_evl, and m_leading_real_evl.

◆ CalcL2ToLinfPressure()

void Nektar::VortexWaveInteraction::CalcL2ToLinfPressure ( void  )

Definition at line 1065 of file VortexWaveInteraction.cpp.

1066 {
1067 
1068  ExecuteWave();
1069 
1070  m_wavePressure->GetPlane(0)->BwdTrans(
1071  m_wavePressure->GetPlane(0)->GetCoeffs(),
1072  m_wavePressure->GetPlane(0)->UpdatePhys());
1073  m_wavePressure->GetPlane(1)->BwdTrans(
1074  m_wavePressure->GetPlane(1)->GetCoeffs(),
1075  m_wavePressure->GetPlane(1)->UpdatePhys());
1076 
1077  int npts = m_waveVelocities[0]->GetPlane(0)->GetNpoints();
1078  NekDouble Linf;
1079  Array<OneD, NekDouble> val(2 * npts, 0.0);
1080 
1081  Vmath::Vmul(npts, m_wavePressure->GetPlane(0)->UpdatePhys(), 1,
1082  m_wavePressure->GetPlane(0)->UpdatePhys(), 1, val, 1);
1083  Vmath::Vvtvp(npts, m_wavePressure->GetPlane(1)->UpdatePhys(), 1,
1084  m_wavePressure->GetPlane(1)->UpdatePhys(), 1, val, 1, val, 1);
1085  cout << "int P^2 " << m_wavePressure->GetPlane(0)->Integral(val) << endl;
1086  Vmath::Vsqrt(npts, val, 1, val, 1);
1087 
1088  Linf = Vmath::Vmax(npts, val, 1);
1089  cout << "Linf: " << Linf << endl;
1090 
1091  NekDouble l2, norm;
1092  l2 =
1093  m_wavePressure->GetPlane(0)->L2(m_wavePressure->GetPlane(0)->GetPhys());
1094  norm = l2 * l2;
1095  l2 =
1096  m_wavePressure->GetPlane(1)->L2(m_wavePressure->GetPlane(1)->GetPhys());
1097  norm += l2 * l2;
1098 
1099  Vmath::Fill(npts, 1.0, val, 1);
1100  NekDouble area = m_waveVelocities[0]->GetPlane(0)->Integral(val);
1101 
1102  l2 = sqrt(norm / area);
1103 
1104  cout << "L2: " << l2 << endl;
1105 
1106  cout << "Ratio Linf/L2: " << Linf / l2 << endl;
1107 }
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:294

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

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

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 1760 of file VortexWaveInteraction.cpp.

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

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 1150 of file VortexWaveInteraction.cpp.

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

References ASSERTL0, and m_sessionName.

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

◆ ExecuteLoop()

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

Definition at line 1199 of file VortexWaveInteraction.cpp.

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

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

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

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 1130 of file VortexWaveInteraction.cpp.

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

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 1109 of file VortexWaveInteraction.cpp.

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

References ASSERTL0, and CG_Iterations::savefile.

Referenced by SaveLoopDetails().

◆ SaveLoopDetails()

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

Definition at line 1180 of file VortexWaveInteraction.cpp.

1182 {
1183  // Save NS restart file
1184  SaveFile(m_sessionName + ".rst", SaveDir, i + 1);
1185  // Save Streak Solution
1186  SaveFile(m_sessionName + "_streak.fld", SaveDir, i);
1187  // Save Wave solution output
1188  SaveFile(m_sessionName + ".evl", SaveDir, i);
1189  SaveFile(m_sessionName + "_eig_0", SaveDir, i);
1190  // Save new field file of eigenvalue
1191  SaveFile(m_sessionName + ".fld", SaveDir, i);
1192  if (!(m_sessionVWI->DefinesSolverInfo("INTERFACE")))
1193  {
1194  // Save new forcing file
1195  SaveFile(m_sessionName + ".vwi", SaveDir, i + 1);
1196  }
1197 }
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 1926 of file VortexWaveInteraction.cpp.

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

1922 {
1923  m_dAlphaDWaveForceMag = (m_alpha[0] - alpha_init) / m_waveForceMagStep;
1924 }

References m_alpha, m_dAlphaDWaveForceMag, and m_waveForceMagStep.

Referenced by DoFixedForcingIteration().

◆ UpdateWaveForceMag()

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

Definition at line 1800 of file VortexWaveInteraction.cpp.

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

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