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.
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 {
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
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;
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 = std::to_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:208
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:143
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: MeshGraph.cpp:115
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 1067 of file VortexWaveInteraction.cpp.

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

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

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;
551 std::shared_ptr<DriverArnoldi> solverWave =
552 std::dynamic_pointer_cast<SolverUtils::DriverArnoldi>(
553 GetDriverFactory().CreateInstance(vDriverModule, m_sessionWave,
554 m_graphWave));
555
556 /// Do linearised NavierStokes Session with Modified Arnoldi
557 cout << "Executing wave solution " << endl;
558 solverWave->Execute();
559
560 // Copy file to a rst location for next restart
561 cout << "Executing cp -f session_eig_0 session_eig_0.rst" << endl;
562 CopyFile("_eig_0", "_eig_0.rst");
563
564 // Store data relevant to other operations
565 m_leading_real_evl[0] = solverWave->GetRealEvl()[0];
566 m_leading_imag_evl[0] = solverWave->GetImagEvl()[0];
567
568 // note this will only be true for modified Arnoldi
569 NekDouble realShift = 0.0;
570 if (m_sessionWave->DefinesParameter("RealShift"))
571 {
572 bool defineshift;
573 // only use shift in modifiedArnoldi solver since
574 // implicitly handled in Arpack.
575 m_sessionWave->MatchSolverInfo("Driver", "ModifiedArnoldi", defineshift,
576 true);
577 if (defineshift)
578 {
579 realShift = m_sessionWave->GetParameter("RealShift");
580 }
581 }
582
583 // Set up leading eigenvalue from inverse
584 NekDouble invmag = 1.0 / (m_leading_real_evl[0] * m_leading_real_evl[0] +
586 m_leading_real_evl[0] *= -invmag;
587 m_leading_real_evl[0] += realShift;
588 m_leading_imag_evl[0] *= invmag;
589
590 m_waveVelocities = solverWave->GetEqu()[0]->UpdateFields();
591 m_wavePressure = solverWave->GetEqu()[0]->GetPressure();
592
593 if (m_sessionVWI->DefinesSolverInfo("INTERFACE"))
594 {
595 cout << "Growth =" << m_leading_real_evl[0] << endl;
596 cout << "Phase =" << m_leading_imag_evl[0] << endl;
597 }
598}

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:195
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