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

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

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::MeshGraphIO::Read(), UpdateAlpha(), and Nektar::VWIIterationTypeMap.

◆ ~VortexWaveInteraction()

Nektar::VortexWaveInteraction::~VortexWaveInteraction ( void  )

Definition at line 363 of file VortexWaveInteraction.cpp.

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

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

1069{
1070 ExecuteWave();
1071
1072 m_wavePressure->GetPlane(0)->BwdTrans(
1073 m_wavePressure->GetPlane(0)->GetCoeffs(),
1074 m_wavePressure->GetPlane(0)->UpdatePhys());
1075 m_wavePressure->GetPlane(1)->BwdTrans(
1076 m_wavePressure->GetPlane(1)->GetCoeffs(),
1077 m_wavePressure->GetPlane(1)->UpdatePhys());
1078
1079 int npts = m_waveVelocities[0]->GetPlane(0)->GetNpoints();
1080 NekDouble Linf;
1081 Array<OneD, NekDouble> val(2 * npts, 0.0);
1082
1083 Vmath::Vmul(npts, m_wavePressure->GetPlane(0)->UpdatePhys(), 1,
1084 m_wavePressure->GetPlane(0)->UpdatePhys(), 1, val, 1);
1085 Vmath::Vvtvp(npts, m_wavePressure->GetPlane(1)->UpdatePhys(), 1,
1086 m_wavePressure->GetPlane(1)->UpdatePhys(), 1, val, 1, val, 1);
1087 cout << "int P^2 " << m_wavePressure->GetPlane(0)->Integral(val) << endl;
1088 Vmath::Vsqrt(npts, val, 1, val, 1);
1089
1090 Linf = Vmath::Vmax(npts, val, 1);
1091 cout << "Linf: " << Linf << endl;
1092
1093 NekDouble l2, norm;
1094 l2 =
1095 m_wavePressure->GetPlane(0)->L2(m_wavePressure->GetPlane(0)->GetPhys());
1096 norm = l2 * l2;
1097 l2 =
1098 m_wavePressure->GetPlane(1)->L2(m_wavePressure->GetPlane(1)->GetPhys());
1099 norm += l2 * l2;
1100
1101 Vmath::Fill(npts, 1.0, val, 1);
1102 NekDouble area = m_waveVelocities[0]->GetPlane(0)->Integral(val);
1103
1104 l2 = sqrt(norm / area);
1105
1106 cout << "L2: " << l2 << endl;
1107
1108 cout << "Ratio Linf/L2: " << Linf / l2 << endl;
1109}
Array< OneD, MultiRegions::ExpListSharedPtr > m_waveVelocities
MultiRegions::ExpListSharedPtr m_wavePressure
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.hpp:340
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.hpp:72
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.hpp:366
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.hpp:54
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.hpp:644
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:285

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

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

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) /
1724 << endl;
1725 cout << "Phase tolerance: "
1726 << fabs((m_leading_imag_evl[0] - previous_imag_evl) /
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 {
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:289

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

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

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

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

References CopyFile(), 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 2124 of file VortexWaveInteraction.cpp.

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

References m_alpha.

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 }

References m_alphaStep.

Referenced by main().

◆ GetDAlphaDWaveForceMag()

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

Definition at line 141 of file VortexWaveInteraction.h.

142 {
144 }

References m_dAlphaDWaveForceMag.

◆ GetEigRelTol()

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

Definition at line 151 of file VortexWaveInteraction.h.

153 {
154 return m_eigRelTol;
155 }

References m_eigRelTol.

Referenced by DoFixedForcingIteration().

◆ GetIterEnd()

int Nektar::VortexWaveInteraction::GetIterEnd ( )
inline

Definition at line 101 of file VortexWaveInteraction.h.

102 {
103 return m_iterEnd;
104 }

References m_iterEnd.

Referenced by DoFixedForcingIteration().

◆ GetIterStart()

int Nektar::VortexWaveInteraction::GetIterStart ( )
inline

Definition at line 96 of file VortexWaveInteraction.h.

97 {
98 return m_iterStart;
99 }

References m_iterStart.

Referenced by DoFixedForcingIteration().

◆ GetMaxOuterIterations()

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

Definition at line 116 of file VortexWaveInteraction.h.

117 {
119 }

References m_maxOuterIterations.

Referenced by DoFixedForcingIteration().

◆ GetMaxWaveForceMagIter()

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

Definition at line 146 of file VortexWaveInteraction.h.

147 {
149 }

References m_maxWaveForceMagIter.

Referenced by main().

◆ GetMinInnerIterations()

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

Definition at line 157 of file VortexWaveInteraction.h.

158 {
160 }

References m_minInnerIterations.

Referenced by DoFixedForcingIteration().

◆ GetNOuterIterations()

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

Definition at line 111 of file VortexWaveInteraction.h.

112 {
113 return m_nOuterIterations;
114 }

References m_nOuterIterations.

Referenced by DoFixedForcingIteration().

◆ GetPrevAlpha()

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

Definition at line 162 of file VortexWaveInteraction.h.

163 {
164 return m_prevAlpha;
165 }

References m_prevAlpha.

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 = 1e-5;
2056 NekDouble xnew, ynew;
2057
2058 int start = npts - 1;
2059 int e_npts;
2060
2061 bool useOnlyQuads = false;
2062 if (m_sessionVWI->DefinesSolverInfo("SymmetriseOnlyQuads"))
2063 {
2064 useOnlyQuads = true;
2065 }
2066
2067 int cnt;
2068 for (int e = 0; e < nel; ++e)
2069 {
2070 e_npts = m_waveVelocities[0]->GetExp(e)->GetTotPoints();
2071 cnt = m_waveVelocities[0]->GetPhys_Offset(e);
2072
2073 if (useOnlyQuads)
2074 {
2075 if (m_waveVelocities[0]->GetExp(e)->DetShapeType() ==
2077 {
2078 for (i = 0; i < e_npts; ++i)
2079 {
2080 index[cnt + i] = -1;
2081 }
2082 continue;
2083 }
2084 }
2085
2086 for (i = cnt; i < cnt + e_npts; ++i)
2087 {
2088 xnew = -coord_x[i] + xmax;
2089 ynew = -coord_y[i];
2090
2091 for (j = start; j >= 0; --j)
2092 {
2093 if ((coord_x[j] - xnew) * (coord_x[j] - xnew) +
2094 (coord_y[j] - ynew) * (coord_y[j] - ynew) <
2095 tol)
2096 {
2097 index[i] = j;
2098 start = j;
2099 break;
2100 }
2101 }
2102
2103 if (j == -1)
2104 {
2105
2106 for (j = npts - 1; j > start; --j)
2107 {
2108
2109 if ((coord_x[j] - xnew) * (coord_x[j] - xnew) +
2110 (coord_y[j] - ynew) * (coord_y[j] - ynew) <
2111 tol)
2112 {
2113 index[i] = j;
2114 break;
2115 }
2116 }
2117 ASSERTL0(j != start, "Failed to find matching point");
2118 }
2119 }
2120 }
2121 return index;
2122}

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 }

References m_VWIIterationType.

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 }

References m_waveForceMag.

Referenced by main().

◆ GetWaveForceMagStep()

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

Definition at line 136 of file VortexWaveInteraction.h.

137 {
138 return m_waveForceMagStep;
139 }

References m_waveForceMagStep.

Referenced by main().

◆ IfIterInterface()

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

Definition at line 197 of file VortexWaveInteraction.h.

198 {
199 return m_iterinterface;
200 }

References m_iterinterface.

Referenced by DoFixedForcingIteration().

◆ MoveFile()

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

Definition at line 1131 of file VortexWaveInteraction.cpp.

1132{
1133 static map<string, int> opendir;
1134
1135 if (opendir.find(dir) == opendir.end())
1136 {
1137 // make directory and presume will fail if it already exists
1138 string mkdir = "mkdir " + dir;
1139 ASSERTL0(system(mkdir.c_str()) == 0,
1140 "Failed to make directory '" + dir + "'");
1141 opendir[dir] = 1;
1142 }
1143
1144 string savefile = dir + "/" + file + "." + std::to_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 1111 of file VortexWaveInteraction.cpp.

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

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 }

References m_alpha.

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 }

References m_alphaStep.

◆ SetEigRelTol()

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

Definition at line 177 of file VortexWaveInteraction.h.

178 {
179 m_eigRelTol = tol;
180 }

References m_eigRelTol.

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 }

References m_minInnerIterations.

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 }

References m_prevAlpha.

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 }

References m_waveForceMag.

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.hpp:704

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

◆ 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

◆ 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

◆ 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 GetIterEnd(), and VortexWaveInteraction().

◆ m_iterinterface

bool Nektar::VortexWaveInteraction::m_iterinterface
private

Definition at line 242 of file VortexWaveInteraction.h.

Referenced by IfIterInterface(), and VortexWaveInteraction().

◆ m_iterStart

int Nektar::VortexWaveInteraction::m_iterStart
private

Definition at line 208 of file VortexWaveInteraction.h.

Referenced by GetIterStart(), and 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 GetMaxOuterIterations(), and VortexWaveInteraction().

◆ m_maxWaveForceMagIter

int Nektar::VortexWaveInteraction::m_maxWaveForceMagIter
private

Definition at line 215 of file VortexWaveInteraction.h.

Referenced by GetMaxWaveForceMagIter(), and VortexWaveInteraction().

◆ m_minInnerIterations

int Nektar::VortexWaveInteraction::m_minInnerIterations
private

◆ 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 GetNOuterIterations(), and VortexWaveInteraction().

◆ m_prevAlpha

NekDouble Nektar::VortexWaveInteraction::m_prevAlpha
private

Definition at line 240 of file VortexWaveInteraction.h.

Referenced by GetPrevAlpha(), and SetPrevAlpha().

◆ 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 GetVWIIterationType(), and 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