Nektar++
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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 = boost::lexical_cast<std::string>(m_iterStart);
285 cout << "Restarting from iteration " << m_iterStart << endl;
286 std::string rstfile = "cp -f Save/" + m_sessionName + ".rst." +
287 nstr + " " + m_sessionName + ".rst";
288 cout << " " << rstfile << endl;
289 if (system(rstfile.c_str()))
290 {
291 ASSERTL0(false, rstfile.c_str());
292 }
293 std::string vwifile = "cp -f Save/" + m_sessionName + ".vwi." +
294 nstr + " " + m_sessionName + ".vwi";
295 cout << " " << vwifile << endl;
296 if (system(vwifile.c_str()))
297 {
298 ASSERTL0(false, vwifile.c_str());
299 }
300 }
301 break;
302 default:
303 ASSERTL0(false, "Unknown VWIITerationType in restart");
304 }
305 }
306
307 // Check for ConveredSoln to update DAlphaDWaveForce
308 {
309 FILE *fp;
310 if ((fp = fopen("ConvergedSolns", "r")))
311 {
312 char buf[BUFSIZ];
313 std::vector<NekDouble> WaveForce, Alpha;
314 NekDouble waveforce, alpha;
315 while (fgets(buf, BUFSIZ, fp))
316 {
317 sscanf(buf, "%*d:%lf%lf", &waveforce, &alpha);
318 WaveForce.push_back(waveforce);
319 Alpha.push_back(alpha);
320 }
321
322 if (Alpha.size() > 1)
323 {
324 int min_i = 0;
325 NekDouble min_alph = fabs(m_alpha[0] - Alpha[min_i]);
326 // find nearest point
327 for (int i = 1; i < Alpha.size(); ++i)
328 {
329 if (fabs(m_alpha[0] - Alpha[min_i]) < min_alph)
330 {
331 min_i = i;
332 min_alph = fabs(m_alpha[0] - Alpha[min_i]);
333 }
334 }
335
336 // find next nearest point
337 int min_j = (min_i == 0) ? 1 : 0;
338 min_alph = fabs(m_alpha[0] - Alpha[min_j]);
339 for (int i = 0; i < Alpha.size(); ++i)
340 {
341 if (i != min_i)
342 {
343 if (fabs(m_alpha[0] - Alpha[min_j]) < min_alph)
344 {
345 min_j = i;
346 min_alph = fabs(m_alpha[0] - Alpha[min_j]);
347 }
348 }
349 }
350
351 if (fabs(Alpha[min_i] - Alpha[min_j]) > 1e-4)
352 {
354 (Alpha[min_i] - Alpha[min_j]) /
355 (WaveForce[min_i] - WaveForce[min_j]);
356 }
357 }
358 }
359 }
360}
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:144
static SessionReaderSharedPtr CreateInstance(int argc, char *argv[])
Creates an instance of the SessionReader class.
static MeshGraphSharedPtr Read(const LibUtilities::SessionReaderSharedPtr pSession, LibUtilities::DomainRangeShPtr rng=LibUtilities::NullDomainRangeShPtr, bool fillGraph=true, SpatialDomains::MeshGraphSharedPtr partitionedGraph=nullptr)
Definition: MeshGraph.cpp:116
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 1164 of file VortexWaveInteraction.cpp.

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

References m_alpha, m_leading_imag_evl, and m_leading_real_evl.

Referenced by DoFixedForcingIteration(), and main().

◆ AppendEvlToFile() [2/2]

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

Definition at line 1173 of file VortexWaveInteraction.cpp.

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

References m_alpha, m_leading_imag_evl, and m_leading_real_evl.

◆ CalcL2ToLinfPressure()

void Nektar::VortexWaveInteraction::CalcL2ToLinfPressure ( void  )

Definition at line 1067 of file VortexWaveInteraction.cpp.

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

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

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

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

1763{
1764 bool returnval = false;
1765 if (m_sessionRoll->DefinesSolverInfo("INTERFACE"))
1766 {
1767 if (m_sessionVWI->GetSolverInfo("INTERFACE") == "phase")
1768 {
1769 if (abs(m_leading_real_evl[0]) < 1e-4)
1770 {
1771 returnval = true;
1772 }
1773 }
1774 else
1775 {
1776 if (abs(m_leading_real_evl[0]) < 1e-4 &&
1777 abs(m_leading_imag_evl[0]) < 2e-6)
1778 {
1779 returnval = true;
1780 }
1781 }
1782 }
1783 else
1784 {
1788 {
1789 returnval = true;
1790 }
1791 }
1792 if (fabs(m_leading_imag_evl[0]) > 1e-2)
1793 {
1794 cout << "Warning: imaginary eigenvalue is greater than 1e-2" << endl;
1795 }
1796
1797 return returnval;
1798}
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 1152 of file VortexWaveInteraction.cpp.

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

References ASSERTL0, and m_sessionName.

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

◆ ExecuteLoop()

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

Definition at line 1201 of file VortexWaveInteraction.cpp.

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

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

Referenced by ExecuteLoop().

◆ ExecuteStreak()

void Nektar::VortexWaveInteraction::ExecuteStreak ( void  )

Definition at line 494 of file VortexWaveInteraction.cpp.

495{
496 // Create driver
497#if 1
498 std::string vDriverModule;
499 m_sessionStreak->LoadSolverInfo("Driver", vDriverModule, "Standard");
500
502 vDriverModule, m_sessionStreak, m_graphStreak);
503 solverStreak->Execute();
504
505 m_streakField = solverStreak->GetEqu()[0]->UpdateFields();
506#else
507 // Setup and execute Advection Diffusion solver
508 string vEquation = m_sessionStreak->GetSolverInfo("EqType");
509 EquationSystemSharedPtr solverStreak =
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:54
DriverFactory & GetDriverFactory()
Definition: Driver.cpp:67

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

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

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

◆ GetAlpha()

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

Definition at line 121 of file VortexWaveInteraction.h.

122 {
123 return m_alpha[0];
124 }

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

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

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

Referenced by CalcNonLinearWaveForce().

◆ GetVWIIterationType()

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

Definition at line 106 of file VortexWaveInteraction.h.

107 {
108 return m_VWIIterationType;
109 }

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

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

References ASSERTL0, and CG_Iterations::savefile.

Referenced by DoFixedForcingIteration().

◆ SaveFile()

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

Definition at line 1111 of file VortexWaveInteraction.cpp.

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

References ASSERTL0, and CG_Iterations::savefile.

Referenced by SaveLoopDetails().

◆ SaveLoopDetails()

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

Definition at line 1182 of file VortexWaveInteraction.cpp.

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

References m_sessionName, m_sessionVWI, and SaveFile().

Referenced by DoFixedForcingIteration().

◆ SetAlpha()

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

Definition at line 167 of file VortexWaveInteraction.h.

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

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

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

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

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

References m_alpha, m_dAlphaDWaveForceMag, and m_waveForceMagStep.

Referenced by DoFixedForcingIteration().

◆ UpdateWaveForceMag()

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

Definition at line 1802 of file VortexWaveInteraction.cpp.

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

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