Nektar++
VortexWaveInteraction.cpp
Go to the documentation of this file.
1//////////////////////////////////////////////////////////////////////////////
2//
3// File: VortexWaveInteraction.cpp
4//
5// For more information, please see: http://www.nektar.info
6//
7// The MIT License
8//
9// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10// Department of Aeronautics, Imperial College London (UK), and Scientific
11// Computing and Imaging Institute, University of Utah (USA).
12//
13// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description: Vortex Wave Interaction class
32//
33///////////////////////////////////////////////////////////////////////////////
34
35#include <IncNavierStokesSolver/EquationSystems/IncNavierStokes.h>
36#include <MultiRegions/ExpList.h>
37#include <MultiRegions/GlobalLinSysKey.h>
38#include <SolverUtils/DriverModifiedArnoldi.h>
39#include <SpatialDomains/MeshGraphIO.h>
40#include <VortexWaveInteraction/VortexWaveInteraction.h>
41
42using namespace std;
43
44namespace Nektar
45{
46
47VortexWaveInteraction::VortexWaveInteraction(int argc, char *argv[])
48 : m_nOuterIterations(0)
49{
50
51 int storesize; // number of previous iterations to store
52
53 m_sessionName = argv[argc - 1];
54 string meshfile = m_sessionName + ".xml";
55
56 LibUtilities::SessionReaderSharedPtr session;
57
58 std::string VWICondFile(argv[argc - 1]);
59 VWICondFile += "_VWI.xml";
60 std::vector<std::string> VWIFilenames;
61 VWIFilenames.push_back(meshfile);
62 VWIFilenames.push_back(VWICondFile);
63
64 // Create Incompressible NavierStokesSolver session reader.
65 m_sessionVWI =
66 LibUtilities::SessionReader::CreateInstance(argc, argv, VWIFilenames);
67 m_sessionVWI->InitSession();
68
69 m_sessionVWI->LoadParameter("AlphaStep", m_alphaStep, 0.05);
70 m_sessionVWI->LoadParameter("OuterIterationStoreSize", storesize, 10);
71
72 // set a low tol for interfaceVWI
73 m_sessionVWI->LoadParameter("EigenvalueRelativeTol", m_eigRelTol, 1e-3);
74
75 m_sessionVWI->LoadParameter("NeutralPointTolerance", m_neutralPointTol,
76 1e-4);
77 m_sessionVWI->LoadParameter("MaxOuterIterations", m_maxOuterIterations,
78 100);
79
80 m_sessionVWI->LoadParameter("StartIteration", m_iterStart, 0);
81 m_sessionVWI->LoadParameter("EndIteration", m_iterEnd, 0);
82
83 m_sessionVWI->LoadParameter("WaveForceMagStep", m_waveForceMagStep, 0.01);
84 m_sessionVWI->LoadParameter("DAlphaDWaveForceMag", m_dAlphaDWaveForceMag,
85 0.0);
86 m_sessionVWI->LoadParameter("MaxWaveForceMagIter", m_maxWaveForceMagIter,
87 1);
88 m_sessionVWI->LoadParameter("RollForceScale", m_rollForceScale, 1.0);
89
90 if (m_sessionVWI->DefinesSolverInfo("DeltaFcnApprox"))
91 {
92 m_deltaFcnApprox = true;
93 m_sessionVWI->LoadParameter("DeltaFcnDecay", m_deltaFcnDecay,
94 1.0 / 500);
95 }
96 else
97 {
98 m_deltaFcnApprox = false;
99 m_deltaFcnDecay = 0.0;
100 }
101
102 if (m_sessionVWI->DefinesSolverInfo("LinfPressureNorm"))
103 {
104 m_useLinfPressureNorm = true;
105 }
106 else
107 {
108 m_useLinfPressureNorm = false;
109 }
110
111 if (m_sessionVWI->DefinesSolverInfo("INTERFACE"))
112 {
113 m_iterinterface = true;
114 }
115 else
116 {
117 m_iterinterface = false;
118 }
119
120 if (m_sessionVWI->DefinesSolverInfo("MoveMeshToCriticalLayer"))
121 {
122 m_moveMeshToCriticalLayer = true;
123 }
124 else
125 {
126 m_moveMeshToCriticalLayer = false;
127 }
128
129 m_alpha = Array<OneD, NekDouble>(storesize);
130 m_alpha[0] = m_sessionVWI->GetParameter("Alpha");
131 m_waveForceMag = Array<OneD, NekDouble>(storesize);
132 m_waveForceMag[0] = m_sessionVWI->GetParameter("WaveForceMag");
133
134 m_leading_real_evl = Array<OneD, NekDouble>(storesize);
135 m_leading_imag_evl = Array<OneD, NekDouble>(storesize);
136
137 if (m_sessionVWI->DefinesParameter("Relaxation"))
138 {
139 m_vwiRelaxation = m_sessionVWI->GetParameter("Relaxation");
140 // fix minimum number of iterations to be number of
141 // iterations required to make contribution of innitial
142 // forcing to 0.1
143 m_minInnerIterations = (int)(log(0.1) / log(m_vwiRelaxation));
144 }
145 else
146 {
147 m_vwiRelaxation = 0.0;
148 m_minInnerIterations = 1;
149 }
150
151 // Initialise NS Roll solver
152 std::string IncCondFile(argv[argc - 1]);
153 IncCondFile += "_IncNSCond.xml";
154 std::vector<std::string> IncNSFilenames;
155 IncNSFilenames.push_back(meshfile);
156 IncNSFilenames.push_back(IncCondFile);
157
158 // Create Incompressible NavierStokesSolver session reader.
159 m_sessionRoll = LibUtilities::SessionReader::CreateInstance(
160 argc, argv, IncNSFilenames, m_sessionVWI->GetComm());
161 m_graphRoll = SpatialDomains::MeshGraphIO::Read(m_sessionRoll);
162 std::string vEquation = m_sessionRoll->GetSolverInfo("SolverType");
163 m_solverRoll = GetEquationSystemFactory().CreateInstance(
164 vEquation, m_sessionRoll, m_graphRoll);
165 m_solverRoll->PrintSummary(cout);
166
167 if (m_sessionRoll->DefinesSolverInfo("INTERFACE"))
168 {
169 m_sessionVWI->LoadParameter("EigenvalueRelativeTol", m_eigRelTol, 1e-2);
170 }
171
172 int ncoeffs = m_solverRoll->UpdateFields()[0]->GetNcoeffs();
173 m_vwiForcing = Array<OneD, Array<OneD, NekDouble>>(4);
174 m_vwiForcing[0] = Array<OneD, NekDouble>(4 * ncoeffs);
175 for (int i = 1; i < 4; ++i)
176 {
177 m_vwiForcing[i] = m_vwiForcing[i - 1] + ncoeffs;
178 }
179
180 // Fill forcing into m_vwiForcing
181 // Has forcing even been initialised yet?
182 // Vmath::Vcopy(ncoeffs,m_solverRoll->UpdateForces()[0]->GetCoeffs(),1,m_vwiForcing[0],1);
183 // Vmath::Vcopy(ncoeffs,m_solverRoll->UpdateForces()[1]->GetCoeffs(),1,m_vwiForcing[1],1);
184
185 // Create AdvDiff Streak solver
186 std::string AdvDiffCondFile(argv[argc - 1]);
187 AdvDiffCondFile += "_AdvDiffCond.xml";
188 std::vector<std::string> AdvDiffFilenames;
189 AdvDiffFilenames.push_back(meshfile);
190 AdvDiffFilenames.push_back(AdvDiffCondFile);
191
192 // Create AdvDiffusion session reader.
193 m_sessionStreak = LibUtilities::SessionReader::CreateInstance(
194 argc, argv, AdvDiffFilenames, m_sessionVWI->GetComm());
195 m_graphStreak = SpatialDomains::MeshGraphIO::Read(m_sessionStreak);
196
197 // Initialise LinNS solver
198 std::string LinNSCondFile(argv[argc - 1]);
199 LinNSCondFile += "_LinNSCond.xml";
200 std::vector<std::string> LinNSFilenames;
201 LinNSFilenames.push_back(meshfile);
202 LinNSFilenames.push_back(LinNSCondFile);
203
204 // Create Linearised NS stability session reader.
205 m_sessionWave = LibUtilities::SessionReader::CreateInstance(
206 argc, argv, LinNSFilenames, m_sessionVWI->GetComm());
207 m_graphWave = SpatialDomains::MeshGraphIO::Read(m_sessionWave);
208
209 // Set the initial beta value in stability to be equal to VWI file
210 std::string LZstr("LZ");
211 NekDouble LZ = 2 * M_PI / m_alpha[0];
212 cout << "Setting LZ in Linearised solver to " << LZ << endl;
213 m_sessionWave->SetParameter(LZstr, LZ);
214
215 // Check for iteration type
216 if (m_sessionVWI->DefinesSolverInfo("VWIIterationType"))
217 {
218 std::string IterationTypeStr =
219 m_sessionVWI->GetSolverInfo("VWIIterationType");
220 for (int i = 0; i < (int)eVWIIterationTypeSize; ++i)
221 {
222 if (boost::iequals(VWIIterationTypeMap[i], IterationTypeStr))
223 {
224 m_VWIIterationType = (VWIIterationType)i;
225 break;
226 }
227 }
228 }
229 else
230 {
231 m_VWIIterationType = eFixedAlphaWaveForcing;
232 }
233
234 // Check for restart
235 bool restart;
236 m_sessionVWI->MatchSolverInfo("RestartIteration", "True", restart, false);
237 if (restart)
238 {
239 switch (m_VWIIterationType)
240 {
241 case eFixedAlpha:
242 break;
243 case eFixedWaveForcing:
244 {
245 FILE *fp;
246 // Check for OuterIter.his file to read
247 if ((fp = fopen("OuterIter.his", "r")))
248 {
249 char buf[BUFSIZ];
250 std::vector<NekDouble> Alpha, Growth, Phase;
251 NekDouble alpha, growth, phase;
252 while (fgets(buf, BUFSIZ, fp))
253 {
254 sscanf(buf, "%*d:%lf%lf%lf", &alpha, &growth, &phase);
255 Alpha.push_back(alpha);
256 Growth.push_back(growth);
257 Phase.push_back(phase);
258 }
259
260 m_nOuterIterations = Alpha.size();
261
262 int nvals =
263 std::min(m_nOuterIterations, (int)m_alpha.size());
264
265 for (int i = 0; i < nvals; ++i)
266 {
267 m_alpha[nvals - 1 - i] =
268 Alpha[m_nOuterIterations - nvals + i];
269 m_leading_real_evl[nvals - 1 - i] =
270 Growth[m_nOuterIterations - nvals + i];
271 m_leading_imag_evl[nvals - 1 - i] =
272 Phase[m_nOuterIterations - nvals + i];
273 }
274
275 UpdateAlpha(m_nOuterIterations++);
276 }
277 else
278 {
279 cout << " No File OuterIter.his to restart from" << endl;
280 }
281 }
282 break;
283 case eFixedAlphaWaveForcing:
284 {
285 string nstr = std::to_string(m_iterStart);
286 cout << "Restarting from iteration " << m_iterStart << endl;
287 std::string rstfile = "cp -f Save/" + m_sessionName + ".rst." +
288 nstr + " " + m_sessionName + ".rst";
289 cout << " " << rstfile << endl;
290 if (system(rstfile.c_str()))
291 {
292 ASSERTL0(false, rstfile.c_str());
293 }
294 std::string vwifile = "cp -f Save/" + m_sessionName + ".vwi." +
295 nstr + " " + m_sessionName + ".vwi";
296 cout << " " << vwifile << endl;
297 if (system(vwifile.c_str()))
298 {
299 ASSERTL0(false, vwifile.c_str());
300 }
301 }
302 break;
303 default:
304 ASSERTL0(false, "Unknown VWIITerationType in restart");
305 }
306 }
307
308 // Check for ConveredSoln to update DAlphaDWaveForce
309 {
310 FILE *fp;
311 if ((fp = fopen("ConvergedSolns", "r")))
312 {
313 char buf[BUFSIZ];
314 std::vector<NekDouble> WaveForce, Alpha;
315 NekDouble waveforce, alpha;
316 while (fgets(buf, BUFSIZ, fp))
317 {
318 sscanf(buf, "%*d:%lf%lf", &waveforce, &alpha);
319 WaveForce.push_back(waveforce);
320 Alpha.push_back(alpha);
321 }
322
323 if (Alpha.size() > 1)
324 {
325 int min_i = 0;
326 NekDouble min_alph = fabs(m_alpha[0] - Alpha[min_i]);
327 // find nearest point
328 for (int i = 1; i < Alpha.size(); ++i)
329 {
330 if (fabs(m_alpha[0] - Alpha[min_i]) < min_alph)
331 {
332 min_i = i;
333 min_alph = fabs(m_alpha[0] - Alpha[min_i]);
334 }
335 }
336
337 // find next nearest point
338 int min_j = (min_i == 0) ? 1 : 0;
339 min_alph = fabs(m_alpha[0] - Alpha[min_j]);
340 for (int i = 0; i < Alpha.size(); ++i)
341 {
342 if (i != min_i)
343 {
344 if (fabs(m_alpha[0] - Alpha[min_j]) < min_alph)
345 {
346 min_j = i;
347 min_alph = fabs(m_alpha[0] - Alpha[min_j]);
348 }
349 }
350 }
351
352 if (fabs(Alpha[min_i] - Alpha[min_j]) > 1e-4)
353 {
354 m_dAlphaDWaveForceMag =
355 (Alpha[min_i] - Alpha[min_j]) /
356 (WaveForce[min_i] - WaveForce[min_j]);
357 }
358 }
359 }
360 }
361}
362
363VortexWaveInteraction::~VortexWaveInteraction()
364{
365 m_sessionVWI->Finalise();
366}
367
368void VortexWaveInteraction::ExecuteRoll(void)
369{
370 // set up the equation system to update the mesh
371 if (m_sessionRoll->DefinesSolverInfo("INTERFACE"))
372 {
373 string vEquation = m_sessionRoll->GetSolverInfo("solvertype");
374 EquationSystemSharedPtr solverRoll =
375 GetEquationSystemFactory().CreateInstance(vEquation, m_sessionRoll,
376 m_graphRoll);
377 // The forcing terms are inserted as N bcs
378 // Execute Roll
379 cout << "Executing Roll solver" << endl;
380 solverRoll->DoInitialise();
381 solverRoll->DoSolve();
382 solverRoll->Output();
383 m_rollField = solverRoll->UpdateFields();
384 for (int g = 0; g < solverRoll->GetNvariables(); ++g)
385 {
386 NekDouble vL2Error = solverRoll->L2Error(g, false);
387 NekDouble vLinfError = solverRoll->LinfError(g);
388 cout << "L 2 error (variable " << solverRoll->GetVariable(g)
389 << ") : " << vL2Error << endl;
390 cout << "L inf error (variable " << solverRoll->GetVariable(g)
391 << ") : " << vLinfError << endl;
392 }
393 }
394 else
395 {
396 if (m_moveMeshToCriticalLayer)
397 {
398 string vEquation = m_sessionRoll->GetSolverInfo("solvertype");
399 EquationSystemSharedPtr solverRoll =
400 GetEquationSystemFactory().CreateInstance(
401 vEquation, m_sessionRoll, m_graphRoll);
402 }
403 else
404 {
405 const int npoints = m_solverRoll->GetNpoints();
406 const int ncoeffs = m_solverRoll->GetNcoeffs();
407
408 static int init = 1;
409 if (init)
410 {
411 m_solverRoll->DoInitialise();
412 m_vwiForcingObj =
413 std::dynamic_pointer_cast<SolverUtils::ForcingProgrammatic>(
414 GetForcingFactory().CreateInstance(
415 "Programmatic", m_sessionRoll, m_solverRoll,
416 m_solverRoll->UpdateFields(),
417 m_solverRoll->UpdateFields().size() - 1, 0));
418
419 std::vector<std::string> vFieldNames =
420 m_sessionRoll->GetVariables();
421 vFieldNames.erase(vFieldNames.end() - 1);
422
423 m_solverRoll->GetFunction("BodyForce")
424 ->Evaluate(vFieldNames, m_vwiForcingObj->UpdateForces());
425
426 // Scale forcing
427 for (int i = 0; i < m_vwiForcingObj->UpdateForces().size(); ++i)
428 {
429 m_solverRoll->UpdateFields()[0]->FwdTrans(
430 m_vwiForcingObj->UpdateForces()[i],
431 m_vwiForcing[2 + i]);
432 Vmath::Smul(npoints, m_rollForceScale,
433 m_vwiForcingObj->UpdateForces()[i], 1,
434 m_vwiForcingObj->UpdateForces()[i], 1);
435 }
436
437 IncNavierStokesSharedPtr ins =
438 m_solverRoll->as<IncNavierStokes>();
439 ins->AddForcing(m_vwiForcingObj);
440
441 init = 0;
442 }
443 else // use internal definition of forcing in m_vwiForcing
444 {
445 // Scale forcing
446 for (int i = 0; i < m_vwiForcingObj->UpdateForces().size(); ++i)
447 {
448 m_solverRoll->UpdateFields()[i]->BwdTrans(
449 m_vwiForcing[i], m_vwiForcingObj->UpdateForces()[i]);
450 Vmath::Smul(npoints, m_rollForceScale,
451 m_vwiForcingObj->UpdateForces()[i], 1,
452 m_vwiForcingObj->UpdateForces()[i], 1);
453 }
454
455 // Shift m_vwiForcing for new restart in case of relaxation
456 Vmath::Vcopy(ncoeffs, m_vwiForcing[0], 1, m_vwiForcing[2], 1);
457 Vmath::Vcopy(ncoeffs, m_vwiForcing[1], 1, m_vwiForcing[3], 1);
458 }
459 }
460
461 // Execute Roll
462 cout << "Executing Roll solver" << endl;
463 m_solverRoll->DoSolve();
464 m_solverRoll->Output();
465 m_rollField = m_solverRoll->UpdateFields();
466 for (int g = 0; g < m_solverRoll->GetNvariables(); ++g)
467 {
468 NekDouble vL2Error = m_solverRoll->L2Error(g, false);
469 NekDouble vLinfError = m_solverRoll->LinfError(g);
470 cout << "L 2 error (variable " << m_solverRoll->GetVariable(g)
471 << ") : " << vL2Error << endl;
472 cout << "L inf error (variable " << m_solverRoll->GetVariable(g)
473 << ") : " << vLinfError << endl;
474 }
475 }
476
477 // Copy .fld file to .rst and base.fld
478 cout << "Executing cp -f session.fld session.rst" << endl;
479 CopyFile(".fld", ".rst");
480
481 // Write out data into base flow with variable Vx,Vy
482 cout << "Writing data to session-Base.fld" << endl;
483
484 std::vector<std::string> variables(2);
485 variables[0] = "Vx";
486 variables[1] = "Vy";
487 std::vector<Array<OneD, NekDouble>> outfield(2);
488 outfield[0] = m_solverRoll->UpdateFields()[0]->UpdateCoeffs();
489 outfield[1] = m_solverRoll->UpdateFields()[1]->UpdateCoeffs();
490 std::string outname = m_sessionName + "-Base.fld";
491 m_solverRoll->WriteFld(outname, m_solverRoll->UpdateFields()[0], outfield,
492 variables);
493}
494
495void VortexWaveInteraction::ExecuteStreak(void)
496{
497 // Create driver
498#if 1
499 std::string vDriverModule;
500 m_sessionStreak->LoadSolverInfo("Driver", vDriverModule, "Standard");
501
502 DriverSharedPtr solverStreak = GetDriverFactory().CreateInstance(
503 vDriverModule, m_sessionStreak, m_graphStreak);
504 solverStreak->Execute();
505
506 m_streakField = solverStreak->GetEqu()[0]->UpdateFields();
507#else
508 // Setup and execute Advection Diffusion solver
509 string vEquation = m_sessionStreak->GetSolverInfo("EqType");
510 EquationSystemSharedPtr solverStreak =
511 GetEquationSystemFactory().CreateInstance(vEquation, m_sessionStreak,
512 m_graphStreak);
513
514 cout << "Executing Streak Solver" << endl;
515 solverStreak->DoInitialise();
516 solverStreak->DoSolve();
517 solverStreak->Output();
518
519 m_streakField = solverStreak->UpdateFields();
520
521 if (m_sessionVWI->DefinesSolverInfo("INTERFACE"))
522 {
523 for (int g = 0; g < solverStreak->GetNvariables(); ++g)
524 {
525 NekDouble vL2Error = solverStreak->L2Error(g, false);
526 NekDouble vLinfError = solverStreak->LinfError(g);
527 cout << "L 2 error (variable " << solverStreak->GetVariable(g)
528 << ") : " << vL2Error << endl;
529 cout << "L inf error (variable " << solverStreak->GetVariable(g)
530 << ") : " << vLinfError << endl;
531 }
532 }
533#endif
534
535 cout << "Executing cp -f session.fld session_streak.fld" << endl;
536 CopyFile(".fld", "_streak.fld");
537}
538
539void VortexWaveInteraction::ExecuteWave(void)
540{
541
542 // Set the initial beta value in stability to be equal to VWI file
543 std::string LZstr("LZ");
544 NekDouble LZ = 2 * M_PI / m_alpha[0];
545 cout << "Setting LZ in Linearised solver to " << LZ << endl;
546 m_sessionWave->SetParameter(LZstr, LZ);
547
548 // Create driver
549 std::string vDriverModule;
550 m_sessionWave->LoadSolverInfo("Driver", vDriverModule, "ModifiedArnoldi");
551 cout << "Setting up linearised NS Solver" << endl;
552 std::shared_ptr<DriverArnoldi> solverWave =
553 std::dynamic_pointer_cast<SolverUtils::DriverArnoldi>(
554 GetDriverFactory().CreateInstance(vDriverModule, m_sessionWave,
555 m_graphWave));
556
557 /// Do linearised NavierStokes Session with Modified Arnoldi
558 cout << "Executing wave solution " << endl;
559 solverWave->Execute();
560
561 // Copy file to a rst location for next restart
562 cout << "Executing cp -f session_eig_0 session_eig_0.rst" << endl;
563 CopyFile("_eig_0", "_eig_0.rst");
564
565 // Store data relevant to other operations
566 m_leading_real_evl[0] = solverWave->GetRealEvl()[0];
567 m_leading_imag_evl[0] = solverWave->GetImagEvl()[0];
568
569 // note this will only be true for modified Arnoldi
570 NekDouble realShift = 0.0;
571 if (m_sessionWave->DefinesParameter("RealShift"))
572 {
573 bool defineshift;
574 // only use shift in modifiedArnoldi solver since
575 // implicitly handled in Arpack.
576 m_sessionWave->MatchSolverInfo("Driver", "ModifiedArnoldi", defineshift,
577 true);
578 if (defineshift)
579 {
580 realShift = m_sessionWave->GetParameter("RealShift");
581 }
582 }
583
584 // Set up leading eigenvalue from inverse
585 NekDouble invmag = 1.0 / (m_leading_real_evl[0] * m_leading_real_evl[0] +
586 m_leading_imag_evl[0] * m_leading_imag_evl[0]);
587 m_leading_real_evl[0] *= -invmag;
588 m_leading_real_evl[0] += realShift;
589 m_leading_imag_evl[0] *= invmag;
590
591 m_waveVelocities = solverWave->GetEqu()[0]->UpdateFields();
592 m_wavePressure = solverWave->GetEqu()[0]->GetPressure();
593
594 if (m_sessionVWI->DefinesSolverInfo("INTERFACE"))
595 {
596 cout << "Growth =" << m_leading_real_evl[0] << endl;
597 cout << "Phase =" << m_leading_imag_evl[0] << endl;
598 }
599}
600
601void VortexWaveInteraction::CalcNonLinearWaveForce(void)
602{
603 if (m_sessionRoll->DefinesSolverInfo("INTERFACE"))
604 {
605 static int cnt = 0;
606 string wavefile = m_sessionName + ".fld";
607 string movedmesh = m_sessionName + "_advPost_moved.xml";
608 string filestreak = m_sessionName + "_streak.fld";
609 string syscall = "";
610 string c = std::to_string(cnt);
611 string c_alpha = std::to_string(m_alpha[0]);
612 if (m_sessionVWI->GetSolverInfo("INTERFACE") == "phase")
613 {
614 string filePost = m_sessionName + "_advPost.xml";
615 syscall = "../../utilities/PostProcessing/Extras/FldCalcBCs " +
616 filePost + " " +
617 "meshhalf_pos_Spen_stability_moved.fld "
618 "meshhalf_pos_Spen_advPost_moved.fld " +
619 c_alpha + " > data_alpha0";
620 cout << syscall.c_str() << endl;
621 if (system(syscall.c_str()))
622 {
623 ASSERTL0(false, syscall.c_str());
624 }
625
626 syscall = "cp -f meshhalf_pos_Spen_stability_moved_u_5.bc " +
627 m_sessionName + "_u_5.bc";
628 cout << syscall.c_str() << endl;
629 if (system(syscall.c_str()))
630 {
631 ASSERTL0(false, syscall.c_str());
632 }
633 syscall = "cp -f meshhalf_pos_Spen_stability_moved_v_5.bc " +
634 m_sessionName + "_v_5.bc";
635 cout << syscall.c_str() << endl;
636 if (system(syscall.c_str()))
637 {
638 ASSERTL0(false, syscall.c_str());
639 }
640 }
641 else
642 {
643 syscall = "../../utilities/PostProcessing/Extras/FldCalcBCs " +
644 movedmesh + " " + wavefile + " " + filestreak + " " +
645 c_alpha + " > datasub_" + c;
646 cout << syscall.c_str() << endl;
647 if (system(syscall.c_str()))
648 {
649 ASSERTL0(false, syscall.c_str());
650 }
651 }
652
653 // relaxation for different alpha values? does it make sense?
654
655 // save the wave
656 string wave_subalp = m_sessionName + "_wave_subalp_" + c + ".fld";
657 syscall = "cp -f " + wavefile + " " + wave_subalp;
658 cout << syscall.c_str() << endl;
659 if (system(syscall.c_str()))
660 {
661 ASSERTL0(false, syscall.c_str());
662 }
663 // FileRelaxation(3);
664 cnt++;
665 }
666 else
667 {
668 int npts = m_waveVelocities[0]->GetPlane(0)->GetNpoints();
669 int ncoeffs = m_waveVelocities[0]->GetPlane(0)->GetNcoeffs();
670 Array<OneD, NekDouble> val(npts), der1(2 * npts);
671 Array<OneD, NekDouble> der2 = der1 + npts;
672 Array<OneD, NekDouble> streak;
673 static int projectfield = -1;
674
675 if (m_deltaFcnApprox)
676 {
677 streak = Array<OneD, NekDouble>(npts);
678 m_streakField[0]->BwdTrans(m_streakField[0]->GetCoeffs(), streak);
679 }
680
681 // Set project field to be first field that has a Neumann
682 // boundary since this not impose any condition on the vertical
683 // boundaries Othersise set to zero.
684 if (projectfield == -1)
685 {
686 Array<OneD, const SpatialDomains::BoundaryConditionShPtr> BndConds;
687
688 for (int i = 0; i < m_waveVelocities.size(); ++i)
689 {
690 BndConds = m_waveVelocities[i]->GetBndConditions();
691 for (int j = 0; j < BndConds.size(); ++j)
692 {
693 if (BndConds[j]->GetBoundaryConditionType() ==
694 SpatialDomains::eNeumann)
695 {
696 projectfield = i;
697 break;
698 }
699 }
700 if (projectfield != -1)
701 {
702 break;
703 }
704 }
705 if (projectfield == -1)
706 {
707 cout << "using first field to project non-linear forcing which "
708 "imposes a Dirichlet condition"
709 << endl;
710 projectfield = 0;
711 }
712 }
713
714 // Shift m_vwiForcing in case of relaxation
715 Vmath::Vcopy(ncoeffs, m_vwiForcing[0], 1, m_vwiForcing[2], 1);
716 Vmath::Vcopy(ncoeffs, m_vwiForcing[1], 1, m_vwiForcing[3], 1);
717
718 // determine inverse of area normalised field.
719 m_wavePressure->GetPlane(0)->BwdTrans(
720 m_wavePressure->GetPlane(0)->GetCoeffs(),
721 m_wavePressure->GetPlane(0)->UpdatePhys());
722 m_wavePressure->GetPlane(1)->BwdTrans(
723 m_wavePressure->GetPlane(1)->GetCoeffs(),
724 m_wavePressure->GetPlane(1)->UpdatePhys());
725 NekDouble invnorm;
726
727 if (m_useLinfPressureNorm)
728 {
729 Vmath::Vmul(npts, m_wavePressure->GetPlane(0)->UpdatePhys(), 1,
730 m_wavePressure->GetPlane(0)->UpdatePhys(), 1, der1, 1);
731 Vmath::Vvtvp(npts, m_wavePressure->GetPlane(1)->UpdatePhys(), 1,
732 m_wavePressure->GetPlane(1)->UpdatePhys(), 1, der1, 1,
733 der1, 1);
734 Vmath::Vsqrt(npts, der1, 1, der1, 1);
735
736 NekDouble Linf = Vmath::Vmax(npts, der1, 1);
737
738 invnorm = 1.0 / Linf;
739 }
740 else
741 {
742 // Determine normalisation of pressure so that |P|/A = 1
743 NekDouble l2;
744 l2 = m_wavePressure->GetPlane(0)->L2(
745 m_wavePressure->GetPlane(0)->GetPhys());
746 invnorm = l2 * l2;
747 l2 = m_wavePressure->GetPlane(1)->L2(
748 m_wavePressure->GetPlane(1)->GetPhys());
749 invnorm += l2 * l2;
750 Vmath::Fill(2 * npts, 1.0, der1, 1);
751 NekDouble area = m_waveVelocities[0]->GetPlane(0)->Integral(der1);
752 cout << "Area: " << area << endl;
753 invnorm = sqrt(area / invnorm);
754 }
755
756 // Get hold of arrays.
757 m_waveVelocities[0]->GetPlane(0)->BwdTrans(
758 m_waveVelocities[0]->GetPlane(0)->GetCoeffs(),
759 m_waveVelocities[0]->GetPlane(0)->UpdatePhys());
760 Array<OneD, NekDouble> u_real =
761 m_waveVelocities[0]->GetPlane(0)->UpdatePhys();
762 Vmath::Smul(npts, invnorm, u_real, 1, u_real, 1);
763 m_waveVelocities[0]->GetPlane(1)->BwdTrans(
764 m_waveVelocities[0]->GetPlane(1)->GetCoeffs(),
765 m_waveVelocities[0]->GetPlane(1)->UpdatePhys());
766 Array<OneD, NekDouble> u_imag =
767 m_waveVelocities[0]->GetPlane(1)->UpdatePhys();
768 Vmath::Smul(npts, invnorm, u_imag, 1, u_imag, 1);
769 m_waveVelocities[1]->GetPlane(0)->BwdTrans(
770 m_waveVelocities[1]->GetPlane(0)->GetCoeffs(),
771 m_waveVelocities[1]->GetPlane(0)->UpdatePhys());
772 Array<OneD, NekDouble> v_real =
773 m_waveVelocities[1]->GetPlane(0)->UpdatePhys();
774 Vmath::Smul(npts, invnorm, v_real, 1, v_real, 1);
775 m_waveVelocities[1]->GetPlane(1)->BwdTrans(
776 m_waveVelocities[1]->GetPlane(1)->GetCoeffs(),
777 m_waveVelocities[1]->GetPlane(1)->UpdatePhys());
778 Array<OneD, NekDouble> v_imag =
779 m_waveVelocities[1]->GetPlane(1)->UpdatePhys();
780 Vmath::Smul(npts, invnorm, v_imag, 1, v_imag, 1);
781
782 // normalise wave for output
783 Vmath::Smul(2 * ncoeffs, invnorm, m_waveVelocities[0]->UpdateCoeffs(),
784 1, m_waveVelocities[0]->UpdateCoeffs(), 1);
785 Vmath::Smul(2 * ncoeffs, invnorm, m_waveVelocities[1]->UpdateCoeffs(),
786 1, m_waveVelocities[1]->UpdateCoeffs(), 1);
787 Vmath::Smul(2 * ncoeffs, invnorm, m_waveVelocities[2]->UpdateCoeffs(),
788 1, m_waveVelocities[2]->UpdateCoeffs(), 1);
789
790 // dump field
791 {
792 std::vector<std::string> variables(3);
793 std::vector<Array<OneD, NekDouble>> outfield(3);
794 variables[0] = "u_w";
795 variables[1] = "v_w";
796 variables[2] = "w_w";
797 outfield[0] = m_waveVelocities[0]->UpdateCoeffs();
798 outfield[1] = m_waveVelocities[1]->UpdateCoeffs();
799 outfield[2] = m_waveVelocities[2]->UpdateCoeffs();
800 std::string outname = m_sessionName + "_wave.fld";
801 m_solverRoll->WriteFld(outname, m_waveVelocities[0], outfield,
802 variables);
803 }
804
805#if 1
806 int ncoeffs_p = m_wavePressure->GetPlane(0)->GetNcoeffs();
807 Vmath::Smul(ncoeffs_p, invnorm,
808 m_wavePressure->GetPlane(0)->UpdateCoeffs(), 1,
809 m_wavePressure->GetPlane(0)->UpdateCoeffs(), 1);
810 Vmath::Smul(ncoeffs_p, invnorm,
811 m_wavePressure->GetPlane(1)->UpdateCoeffs(), 1,
812 m_wavePressure->GetPlane(1)->UpdateCoeffs(), 1);
813#else
814 m_wavePressure->GetPlane(0)->BwdTrans(
815 m_wavePressure->GetPlane(0)->GetCoeffs(),
816 m_wavePressure->GetPlane(0)->UpdatePhys());
817 Vmath::Smul(npts, invnorm, m_wavePressure->GetPlane(0)->UpdatePhys(), 1,
818 m_wavePressure->GetPlane(0)->UpdatePhys(), 1);
819 m_wavePressure->GetPlane(0)->FwdTrans(
820 m_wavePressure->GetPlane(0)->UpdatePhys(),
821 m_wavePressure->GetPlane(0)->UpdateCoeffs());
822
823 m_wavePressure->GetPlane(1)->BwdTrans(
824 m_wavePressure->GetPlane(1)->GetCoeffs(),
825 m_wavePressure->GetPlane(1)->UpdatePhys());
826 Vmath::Smul(npts, invnorm, m_wavePressure->GetPlane(1)->UpdatePhys(), 1,
827 m_wavePressure->GetPlane(1)->UpdatePhys(), 1);
828 m_wavePressure->GetPlane(1)->FwdTrans(
829 m_wavePressure->GetPlane(1)->UpdatePhys(),
830 m_wavePressure->GetPlane(1)->UpdateCoeffs());
831#endif
832
833 // Calculate non-linear terms for x and y directions
834 // d/dx(u u* + u* u)
835 Vmath::Vmul(npts, u_real, 1, u_real, 1, val, 1);
836 Vmath::Vvtvp(npts, u_imag, 1, u_imag, 1, val, 1, val, 1);
837 Vmath::Smul(npts, 2.0, val, 1, val, 1);
838 m_waveVelocities[0]->GetPlane(0)->PhysDeriv(0, val, der1);
839
840 // d/dy(v u* + v* u)
841 Vmath::Vmul(npts, u_real, 1, v_real, 1, val, 1);
842 Vmath::Vvtvp(npts, u_imag, 1, v_imag, 1, val, 1, val, 1);
843 Vmath::Smul(npts, 2.0, val, 1, val, 1);
844 m_waveVelocities[0]->GetPlane(0)->PhysDeriv(1, val, der2);
845
846 Vmath::Vadd(npts, der1, 1, der2, 1, der1, 1);
847#if 1
848 m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(der1,
849 m_vwiForcing[0]);
850#else
851 m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(
852 der1, m_vwiForcing[0]);
853#endif
854 Vmath::Smul(ncoeffs, -m_waveForceMag[0], m_vwiForcing[0], 1,
855 m_vwiForcing[0], 1);
856
857 // d/dx(u v* + u* v)
858 m_waveVelocities[0]->GetPlane(0)->PhysDeriv(0, val, der1);
859
860 // d/dy(v v* + v* v)
861 Vmath::Vmul(npts, v_real, 1, v_real, 1, val, 1);
862 Vmath::Vvtvp(npts, v_imag, 1, v_imag, 1, val, 1, val, 1);
863 Vmath::Smul(npts, 2.0, val, 1, val, 1);
864 m_waveVelocities[0]->GetPlane(0)->PhysDeriv(1, val, der2);
865
866 Vmath::Vadd(npts, der1, 1, der2, 1, der1, 1);
867
868#if 1
869 m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(der1,
870 m_vwiForcing[1]);
871#else
872 m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(
873 der1, m_vwiForcing[1]);
874#endif
875
876 Vmath::Smul(ncoeffs, -m_waveForceMag[0], m_vwiForcing[1], 1,
877 m_vwiForcing[1], 1);
878
879 // by default the symmetrization is on
880 bool symm = true;
881 m_sessionVWI->MatchSolverInfo("Symmetrization", "True", symm, true);
882#if 0
883 if(symm== true )
884 {
885
886 // Symmetrise forcing
887 //-> Get coordinates
888 Array<OneD, NekDouble> coord(2);
889 Array<OneD, NekDouble> coord_x(npts);
890 Array<OneD, NekDouble> coord_y(npts);
891
892 //-> Impose symmetry (x -> -x + Lx/2, y-> -y) on coordinates
893 m_waveVelocities[0]->GetPlane(0)->GetCoords(coord_x,coord_y);
894 NekDouble xmax = Vmath::Vmax(npts,coord_x,1);
895 Vmath::Neg(npts,coord_x,1);
896 Vmath::Sadd(npts,xmax,coord_x,1,coord_x,1);
897 Vmath::Neg(npts,coord_y,1);
898
899 int i, physoffset;
900
901 //-> Obtain list of expansion element ids for each point.
902 Array<OneD, int> Eid(npts);
903 // This search may not be necessary every iteration
904 for(i = 0; i < npts; ++i)
905 {
906 coord[0] = coord_x[i];
907 coord[1] = coord_y[i];
908
909 // Note this will not quite be symmetric.
910 Eid[i] = m_waveVelocities[0]->GetPlane(0)->GetExpIndex(coord,1e-6);
911 }
912
913 // Interpolate field 0
914 m_waveVelocities[0]->GetPlane(0)->BwdTrans(m_vwiForcing[0],der1);
915 for(i = 0; i < npts; ++i)
916 {
917 physoffset = m_waveVelocities[0]->GetPlane(0)->GetPhys_Offset(Eid[i]);
918 coord[0] = coord_x[i];
919 coord[1] = coord_y[i];
920 der2 [i] = m_waveVelocities[0]->GetPlane(0)->GetExp(Eid[i])->PhysEvaluate(coord, der1 + physoffset);
921 }
922 //-> Average field 0
923 Vmath::Vsub(npts,der1,1,der2,1,der2,1);
924 Vmath::Smul(npts,0.5,der2,1,der2,1);
925#if 1
926 m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(der2,m_vwiForcing[0]);
927#else
928 m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(der2, m_vwiForcing[0]);
929#endif
930
931 //-> Interpoloate field 1
932 m_waveVelocities[0]->GetPlane(0)->BwdTrans(m_vwiForcing[1],der1);
933 for(i = 0; i < npts; ++i)
934 {
935 physoffset = m_waveVelocities[0]->GetPlane(0)->GetPhys_Offset(Eid[i]);
936 coord[0] = coord_x[i];
937 coord[1] = coord_y[i];
938 der2[i] = m_waveVelocities[0]->GetPlane(0)->GetExp(Eid[i])->PhysEvaluate(coord, der1 + physoffset);
939 }
940
941 //-> Average field 1
942 Vmath::Vsub(npts,der1,1,der2,1,der2,1);
943 Vmath::Smul(npts,0.5,der2,1,der2,1);
944#if 1
945 m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(der2,m_vwiForcing[1]);
946#else
947 m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(der2, m_vwiForceing[1]);
948#endif
949 }
950#else
951 int i;
952 if (symm == true)
953 {
954 cout << "symmetrization is active" << endl;
955 static Array<OneD, int> index = GetReflectionIndex();
956
957 m_waveVelocities[0]->GetPlane(0)->BwdTrans(m_vwiForcing[0], der1);
958
959 for (i = 0; i < npts; ++i)
960 {
961 if (index[i] != -1)
962 {
963 val[i] = 0.5 * (der1[i] - der1[index[i]]);
964 }
965 }
966#if 1
967 m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(
968 val, m_vwiForcing[0]);
969#else
970 m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(
971 val, m_vwiForcing[0]);
972#endif
973
974 m_waveVelocities[0]->GetPlane(0)->BwdTrans(m_vwiForcing[1], der2);
975 for (i = 0; i < npts; ++i)
976 {
977 if (index[i] != -1)
978 {
979 val[i] = 0.5 * (der2[i] - der2[index[i]]);
980 }
981 }
982#if 1
983 m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(
984 val, m_vwiForcing[1]);
985#else
986 m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(
987 val, m_vwiForcing[1]);
988#endif
989 }
990
991 Vmath::Vmul(npts, der1, 1, der1, 1, val, 1);
992 Vmath::Vvtvp(npts, der2, 1, der2, 1, val, 1, val, 1);
993 Vmath::Vsqrt(npts, val, 1, val, 1);
994 cout << "F_Linf: " << Vmath::Vmax(npts, val, 1) << endl;
995
996#endif
997
998 if (m_vwiRelaxation)
999 {
1000 Vmath::Smul(ncoeffs, 1.0 - m_vwiRelaxation, m_vwiForcing[0], 1,
1001 m_vwiForcing[0], 1);
1002 Vmath::Svtvp(ncoeffs, m_vwiRelaxation, m_vwiForcing[2], 1,
1003 m_vwiForcing[0], 1, m_vwiForcing[0], 1);
1004
1005 Vmath::Smul(ncoeffs, 1.0 - m_vwiRelaxation, m_vwiForcing[1], 1,
1006 m_vwiForcing[1], 1);
1007 Vmath::Svtvp(ncoeffs, m_vwiRelaxation, m_vwiForcing[3], 1,
1008 m_vwiForcing[1], 1, m_vwiForcing[1], 1);
1009 }
1010
1011 if (m_deltaFcnApprox)
1012 {
1013 for (int j = 0; j < 2; ++j)
1014 {
1015
1016 m_waveVelocities[projectfield]->GetPlane(0)->BwdTrans(
1017 m_vwiForcing[j], der1);
1018 for (int i = 0; i < npts; ++i)
1019 {
1020 der1[i] *= exp(-streak[i] * streak[i] / m_deltaFcnDecay);
1021 }
1022 m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(
1023 der1, m_vwiForcing[j]);
1024 }
1025 }
1026
1027 // dump output
1028 std::vector<std::string> variables(4);
1029 std::vector<Array<OneD, NekDouble>> outfield(4);
1030 variables[0] = "u";
1031 variables[1] = "v";
1032 variables[2] = "pr";
1033 variables[3] = "pi";
1034 outfield[0] = m_vwiForcing[0];
1035 outfield[1] = m_vwiForcing[1];
1036 Array<OneD, NekDouble> soln(npts, 0.0);
1037 m_wavePressure->GetPlane(0)->BwdTrans(
1038 m_wavePressure->GetPlane(0)->GetCoeffs(),
1039 m_wavePressure->GetPlane(0)->UpdatePhys());
1040 outfield[2] = Array<OneD, NekDouble>(ncoeffs);
1041 m_waveVelocities[0]->GetPlane(0)->FwdTransLocalElmt(
1042 m_wavePressure->GetPlane(0)->GetPhys(), outfield[2]);
1043 m_wavePressure->GetPlane(1)->BwdTrans(
1044 m_wavePressure->GetPlane(1)->GetCoeffs(),
1045 m_wavePressure->GetPlane(1)->UpdatePhys());
1046
1047 Vmath::Vmul(npts, m_wavePressure->GetPlane(0)->UpdatePhys(), 1,
1048 m_wavePressure->GetPlane(0)->UpdatePhys(), 1, val, 1);
1049 Vmath::Vvtvp(npts, m_wavePressure->GetPlane(1)->UpdatePhys(), 1,
1050 m_wavePressure->GetPlane(1)->UpdatePhys(), 1, val, 1, val,
1051 1);
1052 cout << "int P^2: " << m_wavePressure->GetPlane(0)->Integral(val)
1053 << endl;
1054 Vmath::Vsqrt(npts, val, 1, val, 1);
1055 cout << "PLinf: " << Vmath::Vmax(npts, val, 1) << endl;
1056
1057 outfield[3] = Array<OneD, NekDouble>(ncoeffs);
1058 m_waveVelocities[1]->GetPlane(0)->FwdTransLocalElmt(
1059 m_wavePressure->GetPlane(1)->GetPhys(), outfield[3]);
1060
1061 std::string outname = m_sessionName + ".vwi";
1062
1063 m_solverRoll->WriteFld(outname, m_waveVelocities[0]->GetPlane(0),
1064 outfield, variables);
1065 }
1066}
1067
1068void VortexWaveInteraction::CalcL2ToLinfPressure(void)
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}
1110
1111void VortexWaveInteraction::SaveFile(string file, string dir, int n)
1112{
1113 static map<string, int> opendir;
1114
1115 if (opendir.find(dir) == opendir.end())
1116 {
1117 // make directory and presume will fail if it already exists
1118 string mkdir = "mkdir " + dir;
1119 ASSERTL0(system(mkdir.c_str()) == 0,
1120 "Failed to make directory '" + dir + "'");
1121
1122 opendir[dir] = 1;
1123 }
1124
1125 string savefile = dir + "/" + file + "." + std::to_string(n);
1126 string syscall = "cp -f " + file + " " + savefile;
1127
1128 ASSERTL0(system(syscall.c_str()) == 0, syscall.c_str());
1129}
1130
1131void VortexWaveInteraction::MoveFile(string file, string dir, int n)
1132{
1133 static map<string, int> opendir;
1134
1135 if (opendir.find(dir) == opendir.end())
1136 {
1137 // make directory and presume will fail if it already exists
1138 string mkdir = "mkdir " + dir;
1139 ASSERTL0(system(mkdir.c_str()) == 0,
1140 "Failed to make directory '" + dir + "'");
1141 opendir[dir] = 1;
1142 }
1143
1144 string savefile = dir + "/" + file + "." + std::to_string(n);
1145 string syscall = "mv -f " + file + " " + savefile;
1146
1147 ASSERTL0(system(syscall.c_str()) == 0, syscall.c_str());
1148}
1149
1150void VortexWaveInteraction::CopyFile(string file1end, string file2end)
1151{
1152 string cpfile1 = m_sessionName + file1end;
1153 string cpfile2 = m_sessionName + file2end;
1154 string syscall = "cp -f " + cpfile1 + " " + cpfile2;
1155
1156 if (system(syscall.c_str()))
1157 {
1158 ASSERTL0(false, syscall.c_str());
1159 }
1160}
1161
1162void VortexWaveInteraction::AppendEvlToFile(string file, int n)
1163{
1164 FILE *fp;
1165 fp = fopen(file.c_str(), "a");
1166 fprintf(fp, "%d: %lf %16.12le %16.12le\n", n, m_alpha[0],
1167 m_leading_real_evl[0], m_leading_imag_evl[0]);
1168 fclose(fp);
1169}
1170
1171void VortexWaveInteraction::AppendEvlToFile(string file, NekDouble WaveForceMag)
1172{
1173 FILE *fp;
1174 fp = fopen(file.c_str(), "a");
1175 fprintf(fp, "%lf %lf %16.12le %16.12le\n", WaveForceMag, m_alpha[0],
1176 m_leading_real_evl[0], m_leading_imag_evl[0]);
1177 fclose(fp);
1178}
1179
1180void VortexWaveInteraction::SaveLoopDetails(std::string SaveDir, int i)
1181
1182{
1183 // Save NS restart file
1184 SaveFile(m_sessionName + ".rst", SaveDir, i + 1);
1185 // Save Streak Solution
1186 SaveFile(m_sessionName + "_streak.fld", SaveDir, i);
1187 // Save Wave solution output
1188 SaveFile(m_sessionName + ".evl", SaveDir, i);
1189 SaveFile(m_sessionName + "_eig_0", SaveDir, i);
1190 // Save new field file of eigenvalue
1191 SaveFile(m_sessionName + ".fld", SaveDir, i);
1192 if (!(m_sessionVWI->DefinesSolverInfo("INTERFACE")))
1193 {
1194 // Save new forcing file
1195 SaveFile(m_sessionName + ".vwi", SaveDir, i + 1);
1196 }
1197}
1198
1199void VortexWaveInteraction::ExecuteLoop(bool CalcWaveForce)
1200{
1201
1202 // the global info has to be written in the
1203 // roll session file to use the interface loop
1204 if (m_sessionRoll->DefinesSolverInfo("INTERFACE"))
1205 {
1206 static int cnt = 0;
1207 bool skiprollstreak = false;
1208 if (cnt == 0 && m_sessionVWI->GetParameter("rollstreakfromit") == 1)
1209 {
1210 skiprollstreak = true;
1211 cout << "skip roll-streak at the first iteration" << endl;
1212 }
1213
1214 if (skiprollstreak != true)
1215 {
1216
1217 LibUtilities::NekManager<
1218 MultiRegions::GlobalLinSysKey,
1219 MultiRegions::GlobalLinSys>::EnableManagement("GlobalLinSys");
1220 ExecuteRoll();
1221 LibUtilities::NekManager<
1222 MultiRegions::GlobalLinSysKey,
1223 MultiRegions::GlobalLinSys>::DisableManagement("GlobalLinSys");
1224#ifndef _WIN32
1225 sleep(3);
1226#endif
1227 ExecuteStreak();
1228#ifndef _WIN32
1229 sleep(3);
1230#endif
1231 }
1232
1233 string syscall;
1234 string movedmesh = m_sessionName + "_advPost_moved.xml";
1235 string movedinterpmesh = m_sessionName + "_interp_moved.xml";
1236 // rewrite the Rollsessionfile (we start from the waleffe forcing)
1237 // string meshbndjumps = m_sessionName +"_bndjumps.xml";
1238 // if(cnt==0)
1239 //{
1240 // take the conditions tag from meshbndjumps and copy into
1241 // the rolls session file
1242 //}
1243
1244 string c = std::to_string(cnt);
1245
1246 // save old roll solution
1247 string oldroll = m_sessionName + "_roll_" + c + ".fld";
1248 syscall = "cp -f " + m_sessionName + "-Base.fld" + " " + oldroll;
1249 cout << syscall.c_str() << endl;
1250 if (system(syscall.c_str()))
1251 {
1252 ASSERTL0(false, syscall.c_str());
1253 }
1254 // define file names
1255 string filePost = m_sessionName + "_advPost.xml";
1256 string filestreak = m_sessionName + "_streak.fld";
1257 string filewave = m_sessionName + "_wave.fld";
1258 string filewavepressure = m_sessionName + "_wave_p_split.fld";
1259 string fileinterp = m_sessionName + "_interp.xml";
1260 string interpstreak = m_sessionName + "_interpstreak_" + c + ".fld";
1261 string interwavepressure =
1262 m_sessionName + "_wave_p_split_interp_" + c + ".fld";
1263 string c_alpha = std::to_string(m_alpha[0]);
1264 cout << "alpha = " << m_alpha[0] << endl;
1265
1266 if (m_sessionVWI->GetSolverInfo("INTERFACE") != "phase")
1267 {
1268 cout << "zerophase" << endl;
1269
1270 syscall = "../../utilities/PostProcessing/Extras/MoveMesh " +
1271 filePost + " " + filestreak + " " + fileinterp + " " +
1272 c_alpha;
1273
1274 cout << syscall.c_str() << endl;
1275 if (system(syscall.c_str()))
1276 {
1277 ASSERTL0(false, syscall.c_str());
1278 }
1279
1280 // move the advPost mesh (remark update alpha!!!)
1281 syscall = "../../utilities/PostProcessing/Extras/MoveMesh " +
1282 filePost + " " + filestreak + " " + filePost + " " +
1283 c_alpha;
1284 cout << syscall.c_str() << endl;
1285 if (system(syscall.c_str()))
1286 {
1287 ASSERTL0(false, syscall.c_str());
1288 }
1289
1290 // save oldstreak
1291 string oldstreak = m_sessionName + "_streak_" + c + ".fld";
1292 syscall = "cp -f " + filestreak + " " + oldstreak;
1293 cout << syscall.c_str() << endl;
1294 if (system(syscall.c_str()))
1295 {
1296 ASSERTL0(false, syscall.c_str());
1297 }
1298
1299 // interpolate the streak field into the new mesh
1300 string movedmesh = m_sessionName + "_advPost_moved.xml";
1301 string movedinterpmesh = m_sessionName + "_interp_moved.xml";
1302
1303 // create the interp streak
1304
1305 syscall = "../../utilities/PostProcessing/Extras/FieldToField " +
1306 fileinterp + " " + filestreak + " " + movedinterpmesh +
1307 " " + interpstreak;
1308
1309 cout << syscall.c_str() << endl;
1310 if (system(syscall.c_str()))
1311 {
1312 ASSERTL0(false, syscall.c_str());
1313 }
1314
1315 // save the old mesh
1316 string meshfile = m_sessionName + ".xml";
1317 string meshold = m_sessionName + "_" + c + ".xml";
1318 syscall = "cp -f " + meshfile + " " + meshold;
1319 cout << syscall.c_str() << endl;
1320 if (system(syscall.c_str()))
1321 {
1322 ASSERTL0(false, syscall.c_str());
1323 }
1324
1325 // overwriting the meshfile with the new mesh
1326 syscall = "cp -f " + movedmesh + " " + meshfile;
1327 cout << syscall.c_str() << endl;
1328 if (system(syscall.c_str()))
1329 {
1330 ASSERTL0(false, syscall.c_str());
1331 }
1332
1333 // overwriting the streak file!!
1334 syscall = "cp -f " + interpstreak + " " + filestreak;
1335 cout << syscall.c_str() << endl;
1336 if (system(syscall.c_str()))
1337 {
1338 ASSERTL0(false, syscall.c_str());
1339 }
1340
1341 // calculate the wave
1342 ExecuteWave();
1343
1344 // save the wave field:
1345 string oldwave = m_sessionName + "_wave_" + c + ".fld";
1346 syscall = "cp -f " + m_sessionName + ".fld" + " " + oldwave;
1347 cout << syscall.c_str() << endl;
1348 if (system(syscall.c_str()))
1349 {
1350 ASSERTL0(false, syscall.c_str());
1351 }
1352
1353 // save old jump conditions:
1354 string ujump = m_sessionName + "_u_5.bc";
1355 syscall = "cp -f " + ujump + " " + m_sessionName + "_u_5.bc_" + c;
1356 cout << syscall.c_str() << endl;
1357 if (system(syscall.c_str()))
1358 {
1359 ASSERTL0(false, syscall.c_str());
1360 }
1361
1362 string vjump = m_sessionName + "_v_5.bc";
1363 syscall = "cp -f " + vjump + " " + m_sessionName + "_v_5.bc_" + c;
1364 cout << syscall.c_str() << endl;
1365 if (system(syscall.c_str()))
1366 {
1367 ASSERTL0(false, syscall.c_str());
1368 }
1369 cnt++;
1370 c = std::to_string(cnt);
1371
1372 // use relaxation
1373 if (GetVWIIterationType() !=
1374 eFixedWaveForcingWithSubIterationOnAlpha)
1375 {
1376 // the critical layer should be the bnd region 3
1377 // int reg =3;
1378 // FileRelaxation(reg);
1379 }
1380 // calculate the jump conditions
1381 string wavefile = m_sessionName + ".fld";
1382 syscall = "../../utilities/PostProcessing/Extras/FldCalcBCs " +
1383 movedmesh + " " + wavefile + " " + interpstreak +
1384 "> data" + c;
1385 cout << syscall.c_str() << endl;
1386 if (system(syscall.c_str()))
1387 {
1388 ASSERTL0(false, syscall.c_str());
1389 }
1390
1391 // move the new name_interp_moved.xml into name_interp.xml
1392 syscall = "cp -f " + movedinterpmesh + " " + fileinterp;
1393 cout << syscall.c_str() << endl;
1394 if (system(syscall.c_str()))
1395 {
1396 ASSERTL0(false, syscall.c_str());
1397 }
1398 // move the new name_advPost_moved.xml into name_advPost.xml
1399 syscall = "cp -f " + movedmesh + " " + filePost;
1400 cout << syscall.c_str() << endl;
1401 if (system(syscall.c_str()))
1402 {
1403 ASSERTL0(false, syscall.c_str());
1404 }
1405 }
1406 else if (m_sessionVWI->GetSolverInfo("INTERFACE") == "phase")
1407 {
1408 cout << "phase" << endl;
1409 // determine cr:
1410 NekDouble cr;
1411 string cr_str;
1412 stringstream st;
1413
1414 // calculate the wave
1415 ExecuteWave();
1416
1417 // save oldstreak
1418 string oldstreak = m_sessionName + "_streak_" + c + ".fld";
1419 syscall = "cp -f " + filestreak + " " + oldstreak;
1420 cout << syscall.c_str() << endl;
1421 if (system(syscall.c_str()))
1422 {
1423 ASSERTL0(false, syscall.c_str());
1424 }
1425
1426 // save wave
1427 syscall = "cp -f " + m_sessionName + ".fld" + " " + filewave;
1428 cout << syscall.c_str() << endl;
1429 if (system(syscall.c_str()))
1430 {
1431 ASSERTL0(false, syscall.c_str());
1432 }
1433 // save the old mesh
1434 string meshfile = m_sessionName + ".xml";
1435 string meshold = m_sessionName + "_" + c + ".xml";
1436 syscall = "cp -f " + meshfile + " " + meshold;
1437 cout << syscall.c_str() << endl;
1438 if (system(syscall.c_str()))
1439 {
1440 ASSERTL0(false, syscall.c_str());
1441 }
1442
1443 // save the oldwave field:
1444 string oldwave = m_sessionName + "_wave_" + c + ".fld";
1445 syscall = "cp -f " + m_sessionName + ".fld" + " " + oldwave;
1446 cout << syscall.c_str() << endl;
1447 if (system(syscall.c_str()))
1448 {
1449 ASSERTL0(false, syscall.c_str());
1450 }
1451
1452 // save old jump conditions:
1453 string ujump = m_sessionName + "_u_5.bc";
1454 syscall = "cp -f " + ujump + " " + m_sessionName + "_u_5.bc_" + c;
1455 cout << syscall.c_str() << endl;
1456 if (system(syscall.c_str()))
1457 {
1458 ASSERTL0(false, syscall.c_str());
1459 }
1460
1461 string vjump = m_sessionName + "_v_5.bc";
1462 syscall = "cp -f " + vjump + " " + m_sessionName + "_v_5.bc_" + c;
1463 cout << syscall.c_str() << endl;
1464 if (system(syscall.c_str()))
1465 {
1466 ASSERTL0(false, syscall.c_str());
1467 }
1468 cnt++;
1469
1470 cr = m_leading_imag_evl[0] / m_alpha[0];
1471 st << cr;
1472 cr_str = st.str();
1473 cout << "cr=" << cr_str << endl;
1474 // NB -g or NOT!!!
1475 // move the mesh around the critical layer
1476 syscall = "../../utilities/PostProcessing/Extras/MoveMesh " +
1477 filePost + " " + filestreak + " " + fileinterp + " " +
1478 c_alpha + " " + cr_str;
1479
1480 cout << syscall.c_str() << endl;
1481 if (system(syscall.c_str()))
1482 {
1483 ASSERTL0(false, syscall.c_str());
1484 }
1485 // NB -g or NOT!!!
1486 // move the advPost mesh (remark update alpha!!!)
1487 syscall = "../../utilities/PostProcessing/Extras/MoveMesh " +
1488 filePost + " " + filestreak + " " + filePost + " " +
1489 c_alpha + " " + cr_str;
1490 cout << syscall.c_str() << endl;
1491 if (system(syscall.c_str()))
1492 {
1493 ASSERTL0(false, syscall.c_str());
1494 }
1495
1496 // interp streak into the new mesh
1497 syscall = "../../utilities/PostProcessing/Extras/FieldToField " +
1498 fileinterp + " " + filestreak + " " + movedinterpmesh +
1499 " " + interpstreak;
1500
1501 cout << syscall.c_str() << endl;
1502 if (system(syscall.c_str()))
1503 {
1504 ASSERTL0(false, syscall.c_str());
1505 }
1506
1507 // split wave sol
1508 syscall = "../../utilities/PostProcessing/Extras/SplitFld " +
1509 filePost + " " + filewave;
1510
1511 cout << syscall.c_str() << endl;
1512 if (system(syscall.c_str()))
1513 {
1514 ASSERTL0(false, syscall.c_str());
1515 }
1516 // interp wave
1517 syscall = "../../utilities/PostProcessing/Extras/FieldToField " +
1518 filePost + " " + filewavepressure + " " + movedmesh +
1519 " " + interwavepressure;
1520
1521 cout << syscall.c_str() << endl;
1522 if (system(syscall.c_str()))
1523 {
1524 ASSERTL0(false, syscall.c_str());
1525 }
1526
1527 // use relaxation
1528 if (GetVWIIterationType() !=
1529 eFixedWaveForcingWithSubIterationOnAlpha)
1530 {
1531 // the critical layer should be the bnd region 3
1532 // int reg =3;
1533 // FileRelaxation(reg);
1534 }
1535
1536 // cp wavepressure to m_sessionName.fld(to get
1537 // the right bcs names using FldCalcBCs)
1538 syscall =
1539 "cp -f " + interwavepressure + " " + m_sessionName + ".fld";
1540 cout << syscall.c_str() << endl;
1541 if (system(syscall.c_str()))
1542 {
1543 ASSERTL0(false, syscall.c_str());
1544 }
1545
1546 // calculate the jump conditions
1547 // NB -g or NOT!!!
1548 syscall = "../../utilities/PostProcessing/Extras/FldCalcBCs " +
1549 movedmesh + " " + m_sessionName + ".fld" + " " +
1550 interpstreak + "> data" + c;
1551 cout << syscall.c_str() << endl;
1552 if (system(syscall.c_str()))
1553 {
1554 ASSERTL0(false, syscall.c_str());
1555 }
1556
1557 // overwriting the meshfile with the new mesh
1558 syscall = "cp -f " + movedmesh + " " + meshfile;
1559 cout << syscall.c_str() << endl;
1560 if (system(syscall.c_str()))
1561 {
1562 ASSERTL0(false, syscall.c_str());
1563 }
1564
1565 // overwriting the streak file!! (maybe is useless)
1566 syscall = "cp -f " + interpstreak + " " + filestreak;
1567 cout << syscall.c_str() << endl;
1568 if (system(syscall.c_str()))
1569 {
1570 ASSERTL0(false, syscall.c_str());
1571 }
1572 // move the new name_interp_moved.xml into name_interp.xml
1573 syscall = "cp -f " + movedinterpmesh + " " + fileinterp;
1574 cout << syscall.c_str() << endl;
1575 if (system(syscall.c_str()))
1576 {
1577 ASSERTL0(false, syscall.c_str());
1578 }
1579 // move the new name_advPost_moved.xml into name_advPost.xml
1580 syscall = "cp -f " + movedmesh + " " + filePost;
1581 cout << syscall.c_str() << endl;
1582 if (system(syscall.c_str()))
1583 {
1584 ASSERTL0(false, syscall.c_str());
1585 }
1586 }
1587 }
1588 else
1589 {
1590 LibUtilities::NekManager<
1591 MultiRegions::GlobalLinSysKey,
1592 MultiRegions::GlobalLinSys>::EnableManagement("GlobalLinSys");
1593 ExecuteRoll();
1594 LibUtilities::NekManager<
1595 MultiRegions::GlobalLinSysKey,
1596 MultiRegions::GlobalLinSys>::DisableManagement("GlobalLinSys");
1597
1598#ifndef _WIN32
1599 sleep(3);
1600#endif
1601 ExecuteStreak();
1602#ifndef _WIN32
1603 sleep(3);
1604#endif
1605
1606 if (m_moveMeshToCriticalLayer)
1607 {
1608 string syscall;
1609 char alpchar[16] = "";
1610 snprintf(alpchar, 16, "%f", m_alpha[0]);
1611
1612 string filePost = m_sessionName + "_advPost.xml";
1613 string filestreak = m_sessionName + "_streak.fld";
1614 string filewave = m_sessionName + "_wave.fld";
1615 string filewavepressure = m_sessionName + "_wave_p_split.fld";
1616 string fileinterp = m_sessionName + "_interp.xml";
1617 string interpstreak = m_sessionName + "_interpstreak.fld";
1618 string interwavepressure =
1619 m_sessionName + "_wave_p_split_interp.fld";
1620 syscall = "../../utilities/PostProcessing/Extras/MoveMesh " +
1621 filePost + " " + filestreak + " " + fileinterp + " " +
1622 alpchar;
1623
1624 cout << syscall.c_str() << endl;
1625 if (system(syscall.c_str()))
1626 {
1627 ASSERTL0(false, syscall.c_str());
1628 }
1629
1630 // move the advPost mesh (remark update alpha!!!)
1631 syscall = "../../utilities/PostProcessing/Extras/MoveMesh " +
1632 filePost + " " + filestreak + " " + filePost + " " +
1633 alpchar;
1634 cout << syscall.c_str() << endl;
1635 if (system(syscall.c_str()))
1636 {
1637 ASSERTL0(false, syscall.c_str());
1638 }
1639
1640 // save oldstreak
1641 string oldstreak = m_sessionName + "_streak.fld";
1642 syscall = "cp -f " + filestreak + " " + oldstreak;
1643 cout << syscall.c_str() << endl;
1644 if (system(syscall.c_str()))
1645 {
1646 ASSERTL0(false, syscall.c_str());
1647 }
1648
1649 // interpolate the streak field into the new mesh
1650 string movedmesh = m_sessionName + "_advPost_moved.xml";
1651 string movedinterpmesh = m_sessionName + "_interp_moved.xml";
1652
1653 // create the interp streak
1654
1655 syscall = "../../utilities/PostProcessing/Extras/FieldToField " +
1656 fileinterp + " " + filestreak + " " + movedinterpmesh +
1657 " " + interpstreak;
1658
1659 cout << syscall.c_str() << endl;
1660 if (system(syscall.c_str()))
1661 {
1662 ASSERTL0(false, syscall.c_str());
1663 }
1664
1665 // save the old mesh
1666 string meshfile = m_sessionName + ".xml";
1667 string meshold = m_sessionName + ".xml";
1668 syscall = "cp -f " + meshfile + " " + meshold;
1669 cout << syscall.c_str() << endl;
1670 if (system(syscall.c_str()))
1671 {
1672 ASSERTL0(false, syscall.c_str());
1673 }
1674
1675 // overwriting the meshfile with the new mesh
1676 syscall = "cp -f " + movedmesh + " " + meshfile;
1677 cout << syscall.c_str() << endl;
1678 if (system(syscall.c_str()))
1679 {
1680 ASSERTL0(false, syscall.c_str());
1681 }
1682
1683 // overwriting the streak file!!
1684 syscall = "cp -f " + interpstreak + " " + filestreak;
1685 cout << syscall.c_str() << endl;
1686 if (system(syscall.c_str()))
1687 {
1688 ASSERTL0(false, syscall.c_str());
1689 }
1690 }
1691
1692 ExecuteWave();
1693
1694 if (CalcWaveForce)
1695 {
1696 CalcNonLinearWaveForce();
1697 }
1698 }
1699}
1700
1701bool VortexWaveInteraction::CheckEigIsStationary(bool reset)
1702{
1703 static NekDouble previous_real_evl = -1.0;
1704 static NekDouble previous_imag_evl = -1.0;
1705 static int min_iter = 0;
1706
1707 if (reset)
1708 {
1709 previous_real_evl = -1.0;
1710 min_iter = 0;
1711 }
1712
1713 if (previous_real_evl == -1.0 || min_iter < m_minInnerIterations)
1714 {
1715 previous_real_evl = m_leading_real_evl[0];
1716 previous_imag_evl = m_leading_imag_evl[0];
1717 min_iter++;
1718 return false;
1719 }
1720
1721 cout << "Growth tolerance: "
1722 << fabs((m_leading_real_evl[0] - previous_real_evl) /
1723 m_leading_real_evl[0])
1724 << endl;
1725 cout << "Phase tolerance: "
1726 << fabs((m_leading_imag_evl[0] - previous_imag_evl) /
1727 m_leading_imag_evl[0])
1728 << endl;
1729
1730 // See if real and imaginary growth have converged to with m_eigRelTol
1731 if ((fabs((m_leading_real_evl[0] - previous_real_evl) /
1732 m_leading_real_evl[0]) < m_eigRelTol) ||
1733 (fabs(m_leading_real_evl[0]) < 1e-6))
1734 {
1735 previous_real_evl = m_leading_real_evl[0];
1736 previous_imag_evl = m_leading_imag_evl[0];
1737 if ((fabs((m_leading_imag_evl[0] - previous_imag_evl) /
1738 m_leading_imag_evl[0]) < m_eigRelTol) ||
1739 (fabs(m_leading_imag_evl[0]) < 1e-6))
1740 {
1741 return true;
1742 }
1743 }
1744 else
1745 {
1746 if (fabs(m_leading_imag_evl[0]) > 1e-2)
1747 {
1748 cout << "Warning: imaginary eigenvalue is greater than 1e-2"
1749 << endl;
1750 }
1751 previous_real_evl = m_leading_real_evl[0];
1752 previous_imag_evl = m_leading_imag_evl[0];
1753 return false;
1754 }
1755 return false;
1756}
1757
1758// Check to see if leading eigenvalue is within tolerance defined
1759// in m_neutralPointTol
1760bool VortexWaveInteraction::CheckIfAtNeutralPoint(void)
1761{
1762 bool returnval = false;
1763 if (m_sessionRoll->DefinesSolverInfo("INTERFACE"))
1764 {
1765 if (m_sessionVWI->GetSolverInfo("INTERFACE") == "phase")
1766 {
1767 if (abs(m_leading_real_evl[0]) < 1e-4)
1768 {
1769 returnval = true;
1770 }
1771 }
1772 else
1773 {
1774 if (abs(m_leading_real_evl[0]) < 1e-4 &&
1775 abs(m_leading_imag_evl[0]) < 2e-6)
1776 {
1777 returnval = true;
1778 }
1779 }
1780 }
1781 else
1782 {
1783 if ((m_leading_real_evl[0] * m_leading_real_evl[0] +
1784 m_leading_imag_evl[0] * m_leading_imag_evl[0]) <
1785 m_neutralPointTol * m_neutralPointTol)
1786 {
1787 returnval = true;
1788 }
1789 }
1790 if (fabs(m_leading_imag_evl[0]) > 1e-2)
1791 {
1792 cout << "Warning: imaginary eigenvalue is greater than 1e-2" << endl;
1793 }
1794
1795 return returnval;
1796}
1797
1798// Similar routine to UpdateAlpha
1799
1800void VortexWaveInteraction::UpdateWaveForceMag(int outeriter)
1801{
1802 NekDouble wavef_new = 0.0;
1803
1804 if (outeriter == 1)
1805 {
1806 m_waveForceMag[1] = m_waveForceMag[0];
1807 if (m_leading_real_evl[0] > 0.0)
1808 {
1809 wavef_new = m_waveForceMag[0] - m_waveForceMagStep;
1810 }
1811 else
1812 {
1813 wavef_new = m_waveForceMag[0] + m_waveForceMagStep;
1814 }
1815 }
1816 else
1817 {
1818 int i;
1819 int nstore = (m_waveForceMag.size() < outeriter) ? m_waveForceMag.size()
1820 : outeriter;
1821 Array<OneD, NekDouble> WaveForce(nstore);
1822 Array<OneD, NekDouble> Growth(nstore);
1823
1824 Vmath::Vcopy(nstore, m_waveForceMag, 1, WaveForce, 1);
1825 Vmath::Vcopy(nstore, m_leading_real_evl, 1, Growth, 1);
1826
1827 // Sort WaveForce Growth values;
1828 NekDouble store;
1829 int k;
1830 for (i = 0; i < nstore; ++i)
1831 {
1832 k = Vmath::Imin(nstore - i, &WaveForce[i], 1);
1833
1834 store = WaveForce[i];
1835 WaveForce[i] = WaveForce[i + k];
1836 WaveForce[i + k] = store;
1837
1838 store = Growth[i];
1839 Growth[i] = Growth[i + k];
1840 Growth[i + k] = store;
1841 }
1842
1843 // See if we have any values that cross zero
1844 for (i = 0; i < nstore - 1; ++i)
1845 {
1846 if (Growth[i] * Growth[i + 1] < 0.0)
1847 {
1848 break;
1849 }
1850 }
1851
1852 if (i != nstore - 1)
1853 {
1854 if (nstore == 2)
1855 {
1856 wavef_new =
1857 (WaveForce[0] * Growth[1] - WaveForce[1] * Growth[0]) /
1858 (Growth[1] - Growth[0]);
1859 }
1860 else
1861 {
1862 // use a quadratic fit and step through 10000 points
1863 // to find zero.
1864 int j;
1865 int nsteps = 10000;
1866 int idx = (i == 0) ? 1 : i;
1867 NekDouble da = WaveForce[idx + 1] - WaveForce[idx - 1];
1868 NekDouble gval_m1 = Growth[idx - 1], a, gval;
1869 NekDouble c1 = Growth[idx - 1] /
1870 (WaveForce[idx - 1] - WaveForce[idx]) /
1871 (WaveForce[idx - 1] - WaveForce[idx + 1]);
1872 NekDouble c2 = Growth[idx] /
1873 (WaveForce[idx] - WaveForce[idx - 1]) /
1874 (WaveForce[idx] - WaveForce[idx + 1]);
1875 NekDouble c3 = Growth[idx + 1] /
1876 (WaveForce[idx + 1] - WaveForce[idx - 1]) /
1877 (WaveForce[idx + 1] - WaveForce[idx]);
1878
1879 for (j = 1; j < nsteps + 1; ++j)
1880 {
1881 a = WaveForce[i] + j * da / (NekDouble)nsteps;
1882 gval =
1883 c1 * (a - WaveForce[idx]) * (a - WaveForce[idx + 1]) +
1884 c2 * (a - WaveForce[idx - 1]) *
1885 (a - WaveForce[idx + 1]) +
1886 c3 * (a - WaveForce[idx - 1]) * (a - WaveForce[idx]);
1887
1888 if (gval * gval_m1 < 0.0)
1889 {
1890 wavef_new = ((a + da / (NekDouble)nsteps) * gval -
1891 a * gval_m1) /
1892 (gval - gval_m1);
1893 break;
1894 }
1895 }
1896 }
1897 }
1898 else // step backward/forward
1899 {
1900 if (Growth[i] > 0.0)
1901 {
1902 wavef_new = m_waveForceMag[0] - m_waveForceMagStep;
1903 }
1904 else
1905 {
1906 wavef_new = m_waveForceMag[0] + m_waveForceMagStep;
1907 }
1908 }
1909 }
1910
1911 for (int i = m_waveForceMag.size() - 1; i > 0; --i)
1912 {
1913 m_waveForceMag[i] = m_waveForceMag[i - 1];
1914 m_leading_real_evl[i] = m_leading_real_evl[i - 1];
1915 m_leading_imag_evl[i] = m_leading_imag_evl[i - 1];
1916 }
1917
1918 m_waveForceMag[0] = wavef_new;
1919}
1920
1921void VortexWaveInteraction::UpdateDAlphaDWaveForceMag(NekDouble alpha_init)
1922{
1923 m_dAlphaDWaveForceMag = (m_alpha[0] - alpha_init) / m_waveForceMagStep;
1924}
1925
1926void VortexWaveInteraction::UpdateAlpha(int outeriter)
1927{
1928 NekDouble alp_new = 0.0;
1929
1930 if (outeriter == 1)
1931 {
1932 m_alpha[1] = m_alpha[0];
1933 if (m_leading_real_evl[0] > 0.0)
1934 {
1935 alp_new = m_alpha[0] + m_alphaStep;
1936 }
1937 else
1938 {
1939 alp_new = m_alpha[0] - m_alphaStep;
1940 }
1941 }
1942 else
1943 {
1944 int i;
1945 int nstore = (m_alpha.size() < outeriter) ? m_alpha.size() : outeriter;
1946 Array<OneD, NekDouble> Alpha(nstore);
1947 Array<OneD, NekDouble> Growth(nstore);
1948
1949 Vmath::Vcopy(nstore, m_alpha, 1, Alpha, 1);
1950 Vmath::Vcopy(nstore, m_leading_real_evl, 1, Growth, 1);
1951
1952 // Sort Alpha Growth values;
1953 NekDouble store;
1954 int k;
1955 for (i = 0; i < nstore; ++i)
1956 {
1957 k = Vmath::Imin(nstore - i, &Alpha[i], 1);
1958
1959 store = Alpha[i];
1960 Alpha[i] = Alpha[i + k];
1961 Alpha[i + k] = store;
1962
1963 store = Growth[i];
1964 Growth[i] = Growth[i + k];
1965 Growth[i + k] = store;
1966 }
1967
1968 // See if we have any values that cross zero
1969 for (i = 0; i < nstore - 1; ++i)
1970 {
1971 if (Growth[i] * Growth[i + 1] < 0.0)
1972 {
1973 break;
1974 }
1975 }
1976
1977 if (i != nstore - 1)
1978 {
1979 if (nstore == 2)
1980 {
1981 alp_new = (Alpha[0] * Growth[1] - Alpha[1] * Growth[0]) /
1982 (Growth[1] - Growth[0]);
1983 }
1984 else
1985 {
1986 // use a quadratic fit and step through 10000 points
1987 // to find zero.
1988 int j;
1989 int nsteps = 10000;
1990 int idx = (i == 0) ? 1 : i;
1991 NekDouble da = Alpha[idx + 1] - Alpha[idx - 1];
1992 NekDouble gval_m1 = Growth[idx - 1], a, gval;
1993 NekDouble c1 = Growth[idx - 1] / (Alpha[idx - 1] - Alpha[idx]) /
1994 (Alpha[idx - 1] - Alpha[idx + 1]);
1995 NekDouble c2 = Growth[idx] / (Alpha[idx] - Alpha[idx - 1]) /
1996 (Alpha[idx] - Alpha[idx + 1]);
1997 NekDouble c3 = Growth[idx + 1] /
1998 (Alpha[idx + 1] - Alpha[idx - 1]) /
1999 (Alpha[idx + 1] - Alpha[idx]);
2000
2001 for (j = 1; j < nsteps + 1; ++j)
2002 {
2003 a = Alpha[i] + j * da / (NekDouble)nsteps;
2004 gval = c1 * (a - Alpha[idx]) * (a - Alpha[idx + 1]) +
2005 c2 * (a - Alpha[idx - 1]) * (a - Alpha[idx + 1]) +
2006 c3 * (a - Alpha[idx - 1]) * (a - Alpha[idx]);
2007
2008 if (gval * gval_m1 < 0.0)
2009 {
2010 alp_new = ((a + da / (NekDouble)nsteps) * gval -
2011 a * gval_m1) /
2012 (gval - gval_m1);
2013 break;
2014 }
2015 }
2016 }
2017 }
2018 else // step backward/forward
2019 {
2020 if (Growth[i] > 0.0)
2021 {
2022 alp_new = m_alpha[0] + m_alphaStep;
2023 }
2024 else
2025 {
2026 alp_new = m_alpha[0] - m_alphaStep;
2027 }
2028 }
2029 }
2030
2031 for (int i = m_alpha.size() - 1; i > 0; --i)
2032 {
2033 m_alpha[i] = m_alpha[i - 1];
2034 m_leading_real_evl[i] = m_leading_real_evl[i - 1];
2035 m_leading_imag_evl[i] = m_leading_imag_evl[i - 1];
2036 }
2037
2038 m_alpha[0] = alp_new;
2039}
2040
2041Array<OneD, int> VortexWaveInteraction::GetReflectionIndex(void)
2042{
2043 int i, j;
2044 int npts = m_waveVelocities[0]->GetPlane(0)->GetNpoints();
2045 int nel = m_waveVelocities[0]->GetNumElmts();
2046 Array<OneD, int> index(npts);
2047
2048 Array<OneD, NekDouble> coord(2);
2049 Array<OneD, NekDouble> coord_x(npts);
2050 Array<OneD, NekDouble> coord_y(npts);
2051
2052 //-> Dermine the point which is on coordinate (x -> -x + Lx/2, y-> -y)
2053 m_waveVelocities[0]->GetPlane(0)->GetCoords(coord_x, coord_y);
2054 NekDouble xmax = Vmath::Vmax(npts, coord_x, 1);
2055 NekDouble tol = 1e-5;
2056 NekDouble xnew, ynew;
2057
2058 int start = npts - 1;
2059 int e_npts;
2060
2061 bool useOnlyQuads = false;
2062 if (m_sessionVWI->DefinesSolverInfo("SymmetriseOnlyQuads"))
2063 {
2064 useOnlyQuads = true;
2065 }
2066
2067 int cnt;
2068 for (int e = 0; e < nel; ++e)
2069 {
2070 e_npts = m_waveVelocities[0]->GetExp(e)->GetTotPoints();
2071 cnt = m_waveVelocities[0]->GetPhys_Offset(e);
2072
2073 if (useOnlyQuads)
2074 {
2075 if (m_waveVelocities[0]->GetExp(e)->DetShapeType() ==
2076 LibUtilities::eTriangle)
2077 {
2078 for (i = 0; i < e_npts; ++i)
2079 {
2080 index[cnt + i] = -1;
2081 }
2082 continue;
2083 }
2084 }
2085
2086 for (i = cnt; i < cnt + e_npts; ++i)
2087 {
2088 xnew = -coord_x[i] + xmax;
2089 ynew = -coord_y[i];
2090
2091 for (j = start; j >= 0; --j)
2092 {
2093 if ((coord_x[j] - xnew) * (coord_x[j] - xnew) +
2094 (coord_y[j] - ynew) * (coord_y[j] - ynew) <
2095 tol)
2096 {
2097 index[i] = j;
2098 start = j;
2099 break;
2100 }
2101 }
2102
2103 if (j == -1)
2104 {
2105
2106 for (j = npts - 1; j > start; --j)
2107 {
2108
2109 if ((coord_x[j] - xnew) * (coord_x[j] - xnew) +
2110 (coord_y[j] - ynew) * (coord_y[j] - ynew) <
2111 tol)
2112 {
2113 index[i] = j;
2114 break;
2115 }
2116 }
2117 ASSERTL0(j != start, "Failed to find matching point");
2118 }
2119 }
2120 }
2121 return index;
2122}
2123
2124void VortexWaveInteraction::FileRelaxation(int reg)
2125{
2126 cout << "relaxation..." << endl;
2127 static int cnt = 0;
2128 Array<OneD, MultiRegions::ExpListSharedPtr> Iexp =
2129 m_rollField[0]->GetBndCondExpansions();
2130 // cast to 1D explist (otherwise appenddata doesn't work)
2131 MultiRegions::ExpListSharedPtr Ilayer;
2132 Ilayer = MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(
2133 *std::static_pointer_cast<MultiRegions::ExpList>(Iexp[reg]));
2134 int nq = Ilayer->GetTotPoints();
2135 if (cnt == 0)
2136 {
2137 m_bcsForcing = Array<OneD, Array<OneD, NekDouble>>(4);
2138 m_bcsForcing[0] = Array<OneD, NekDouble>(4 * nq);
2139 for (int i = 1; i < 4; ++i)
2140 {
2141 m_bcsForcing[i] = m_bcsForcing[i - 1] + nq;
2142 }
2143 }
2144
2145 // Read in mesh from input file
2146
2147 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef_u;
2148 std::vector<std::vector<NekDouble>> FieldData_u;
2149 string file = m_sessionName;
2150
2151 file += "_u_5.bc";
2152 LibUtilities::FieldIOSharedPtr fld =
2153 LibUtilities::FieldIO::CreateDefault(m_sessionVWI);
2154 fld->Import(file, FieldDef_u, FieldData_u);
2155 Ilayer->ExtractDataToCoeffs(FieldDef_u[0], FieldData_u[0],
2156 FieldDef_u[0]->m_fields[0],
2157 Ilayer->UpdateCoeffs());
2158 Ilayer->BwdTrans(Ilayer->GetCoeffs(), Ilayer->UpdatePhys());
2159
2160 if (cnt == 0)
2161 {
2162 Vmath::Vcopy(nq, Ilayer->UpdatePhys(), 1, m_bcsForcing[2], 1);
2163 }
2164 Vmath::Vcopy(nq, Ilayer->UpdatePhys(), 1, m_bcsForcing[0], 1);
2165
2166 // case relaxation==0 an additional array is needed
2167 Array<OneD, NekDouble> tmp_forcing(nq, 0.0);
2168
2169 if (cnt != 0)
2170 {
2171 if (m_vwiRelaxation == 1.0)
2172 {
2173 Vmath::Vcopy(nq, m_bcsForcing[0], 1, tmp_forcing, 1);
2174 }
2175 Vmath::Smul(nq, 1.0 - m_vwiRelaxation, m_bcsForcing[0], 1,
2176 m_bcsForcing[0], 1);
2177 Vmath::Svtvp(nq, m_vwiRelaxation, m_bcsForcing[2], 1, m_bcsForcing[0],
2178 1, Ilayer->UpdatePhys(), 1);
2179 // generate again the bcs files:
2180
2181 Array<OneD, Array<OneD, NekDouble>> fieldcoeffs(1);
2182 Ilayer->FwdTransLocalElmt(Ilayer->GetPhys(), Ilayer->UpdateCoeffs());
2183 fieldcoeffs[0] = Ilayer->UpdateCoeffs();
2184 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef1 =
2185 Ilayer->GetFieldDefinitions();
2186 std::vector<std::vector<NekDouble>> FieldData_1(FieldDef1.size());
2187 ;
2188 FieldDef1[0]->m_fields.push_back("u");
2189 Ilayer->AppendFieldData(FieldDef1[0], FieldData_1[0]);
2190 fld->Write(file, FieldDef1, FieldData_1);
2191 // save the bcs for the next iteration
2192 if (m_vwiRelaxation != 1.0)
2193 {
2194 Vmath::Smul(nq, 1. / (1.0 - m_vwiRelaxation), m_bcsForcing[0], 1,
2195 m_bcsForcing[0], 1);
2196 Vmath::Vcopy(nq, m_bcsForcing[0], 1, m_bcsForcing[2], 1);
2197 }
2198 else
2199 {
2200 Vmath::Vcopy(nq, tmp_forcing, 1, m_bcsForcing[2], 1);
2201 }
2202 }
2203
2204 file = m_sessionName + "_v_5.bc";
2205
2206 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef_v;
2207 std::vector<std::vector<NekDouble>> FieldData_v;
2208 fld->Import(file, FieldDef_v, FieldData_v);
2209 Ilayer->ExtractDataToCoeffs(FieldDef_v[0], FieldData_v[0],
2210 FieldDef_v[0]->m_fields[0],
2211 Ilayer->UpdateCoeffs());
2212 Ilayer->BwdTrans(Ilayer->GetCoeffs(), Ilayer->UpdatePhys());
2213 if (cnt == 0)
2214 {
2215 Vmath::Vcopy(nq, Ilayer->UpdatePhys(), 1, m_bcsForcing[3], 1);
2216 }
2217 Vmath::Vcopy(nq, Ilayer->UpdatePhys(), 1, m_bcsForcing[1], 1);
2218 if (cnt != 0)
2219 {
2220 if (m_vwiRelaxation == 1.0)
2221 {
2222 Vmath::Vcopy(nq, m_bcsForcing[1], 1, tmp_forcing, 1);
2223 }
2224 Vmath::Smul(nq, 1.0 - m_vwiRelaxation, m_bcsForcing[1], 1,
2225 m_bcsForcing[1], 1);
2226 Vmath::Svtvp(nq, m_vwiRelaxation, m_bcsForcing[3], 1, m_bcsForcing[1],
2227 1, Ilayer->UpdatePhys(), 1);
2228 // generate again the bcs files:
2229 Array<OneD, Array<OneD, NekDouble>> fieldcoeffs(1);
2230 Ilayer->FwdTransLocalElmt(Ilayer->GetPhys(), Ilayer->UpdateCoeffs());
2231 fieldcoeffs[0] = Ilayer->UpdateCoeffs();
2232 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef2 =
2233 Ilayer->GetFieldDefinitions();
2234 std::vector<std::vector<NekDouble>> FieldData_2(FieldDef2.size());
2235 ;
2236 FieldDef2[0]->m_fields.push_back("v");
2237 Ilayer->AppendFieldData(FieldDef2[0], FieldData_2[0]);
2238 fld->Write(file, FieldDef2, FieldData_2);
2239 // save the bcs for the next iteration
2240 if (m_vwiRelaxation != 1.0)
2241 {
2242 Vmath::Smul(nq, 1. / (1.0 - m_vwiRelaxation), m_bcsForcing[1], 1,
2243 m_bcsForcing[1], 1);
2244 Vmath::Vcopy(nq, m_bcsForcing[1], 1, m_bcsForcing[3], 1);
2245 }
2246 else
2247 {
2248 Vmath::Vcopy(nq, tmp_forcing, 1, m_bcsForcing[3], 1);
2249 }
2250 }
2251
2252 cnt++;
2253}
2254} // namespace Nektar
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
Nektar::IncNavierStokes
This class is the base class for Navier Stokes problems.
Definition: IncNavierStokes.h:148
Nektar::IncNavierStokes::AddForcing
void AddForcing(const SolverUtils::ForcingSharedPtr &pForce)
Definition: IncNavierStokes.cpp:801
Nektar::LibUtilities::FieldIO::CreateDefault
static std::shared_ptr< FieldIO > CreateDefault(const LibUtilities::SessionReaderSharedPtr session)
Returns an object for the default FieldIO method.
Definition: FieldIO.cpp:194
Nektar::LibUtilities::NekFactory::CreateInstance
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
Definition: BasicUtils/NekFactory.hpp:143
Nektar::LibUtilities::SessionReader::CreateInstance
static SessionReaderSharedPtr CreateInstance(int argc, char *argv[])
Creates an instance of the SessionReader class.
Definition: SessionReader.h:138
Nektar::MemoryManager::AllocateSharedPtr
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:166
Nektar::MultiRegions::GlobalLinSys
A global linear system.
Definition: GlobalLinSys.h:70
Nektar::SpatialDomains::MeshGraphIO::Read
static MeshGraphSharedPtr Read(const LibUtilities::SessionReaderSharedPtr pSession, LibUtilities::DomainRangeShPtr rng=LibUtilities::NullDomainRangeShPtr, bool fillGraph=true, SpatialDomains::MeshGraphSharedPtr partitionedGraph=nullptr)
Definition: MeshGraphIO.cpp:51
Nektar::VortexWaveInteraction::m_streakField
Array< OneD, MultiRegions::ExpListSharedPtr > m_streakField
Definition: VortexWaveInteraction.h:254
Nektar::VortexWaveInteraction::m_vwiForcing
Array< OneD, Array< OneD, NekDouble > > m_vwiForcing
Definition: VortexWaveInteraction.h:249
Nektar::VortexWaveInteraction::CheckEigIsStationary
bool CheckEigIsStationary(bool reset=false)
Definition: VortexWaveInteraction.cpp:1701
Nektar::VortexWaveInteraction::m_sessionVWI
LibUtilities::SessionReaderSharedPtr m_sessionVWI
Definition: VortexWaveInteraction.h:259
Nektar::VortexWaveInteraction::m_sessionRoll
LibUtilities::SessionReaderSharedPtr m_sessionRoll
Definition: VortexWaveInteraction.h:261
Nektar::VortexWaveInteraction::m_alpha
Array< OneD, NekDouble > m_alpha
< Leading imaginary eigenvalue
Definition: VortexWaveInteraction.h:233
Nektar::VortexWaveInteraction::m_sessionStreak
LibUtilities::SessionReaderSharedPtr m_sessionStreak
Definition: VortexWaveInteraction.h:265
Nektar::VortexWaveInteraction::SaveFile
void SaveFile(std::string fileend, std::string dir, int n)
Definition: VortexWaveInteraction.cpp:1111
Nektar::VortexWaveInteraction::CopyFile
void CopyFile(std::string file1end, std::string file2end)
Definition: VortexWaveInteraction.cpp:1150
Nektar::VortexWaveInteraction::GetReflectionIndex
Array< OneD, int > GetReflectionIndex(void)
Definition: VortexWaveInteraction.cpp:2041
Nektar::VortexWaveInteraction::VortexWaveInteraction
VortexWaveInteraction(int argc, char *argv[])
Definition: VortexWaveInteraction.cpp:47
Nektar::VortexWaveInteraction::SaveLoopDetails
void SaveLoopDetails(std::string dir, int i)
Definition: VortexWaveInteraction.cpp:1180
Nektar::VortexWaveInteraction::m_vwiForcingObj
SolverUtils::ForcingProgrammaticSharedPtr m_vwiForcingObj
Definition: VortexWaveInteraction.h:250
Nektar::VortexWaveInteraction::m_waveVelocities
Array< OneD, MultiRegions::ExpListSharedPtr > m_waveVelocities
Definition: VortexWaveInteraction.h:246
Nektar::VortexWaveInteraction::ExecuteLoop
void ExecuteLoop(bool CalcWaveForce=true)
Definition: VortexWaveInteraction.cpp:1199
Nektar::VortexWaveInteraction::m_graphStreak
SpatialDomains::MeshGraphSharedPtr m_graphStreak
Definition: VortexWaveInteraction.h:266
Nektar::VortexWaveInteraction::m_bcsForcing
Array< OneD, Array< OneD, NekDouble > > m_bcsForcing
Definition: VortexWaveInteraction.h:252
Nektar::VortexWaveInteraction::m_graphRoll
SpatialDomains::MeshGraphSharedPtr m_graphRoll
Definition: VortexWaveInteraction.h:262
Nektar::VortexWaveInteraction::m_rollField
Array< OneD, MultiRegions::ExpListSharedPtr > m_rollField
Definition: VortexWaveInteraction.h:256
Nektar::VortexWaveInteraction::MoveFile
void MoveFile(std::string fileend, std::string dir, int n)
Definition: VortexWaveInteraction.cpp:1131
Nektar::VortexWaveInteraction::m_graphWave
SpatialDomains::MeshGraphSharedPtr m_graphWave
Definition: VortexWaveInteraction.h:268
Nektar::VortexWaveInteraction::m_waveForceMag
Array< OneD, NekDouble > m_waveForceMag
Definition: VortexWaveInteraction.h:224
Nektar::VortexWaveInteraction::AppendEvlToFile
void AppendEvlToFile(std::string file, int n)
Definition: VortexWaveInteraction.cpp:1162
Nektar::VortexWaveInteraction::m_sessionWave
LibUtilities::SessionReaderSharedPtr m_sessionWave
Definition: VortexWaveInteraction.h:267
Nektar::VortexWaveInteraction::UpdateDAlphaDWaveForceMag
void UpdateDAlphaDWaveForceMag(NekDouble alphainit)
Definition: VortexWaveInteraction.cpp:1921
Nektar::VortexWaveInteraction::m_wavePressure
MultiRegions::ExpListSharedPtr m_wavePressure
Definition: VortexWaveInteraction.h:247
Nektar::VortexWaveInteraction::m_leading_real_evl
Array< OneD, NekDouble > m_leading_real_evl
Definition: VortexWaveInteraction.h:229
Nektar::VortexWaveInteraction::m_leading_imag_evl
Array< OneD, NekDouble > m_leading_imag_evl
< Leading real eigenvalue
Definition: VortexWaveInteraction.h:231
Nektar::VortexWaveInteraction::GetVWIIterationType
VWIIterationType GetVWIIterationType(void)
Definition: VortexWaveInteraction.h:106
Nektar::VortexWaveInteraction::m_solverRoll
EquationSystemSharedPtr m_solverRoll
Definition: VortexWaveInteraction.h:263
CG_Iterations.savefile
string savefile
Definition: CG_Iterations.py:312
Nektar::LibUtilities::FieldIOSharedPtr
std::shared_ptr< FieldIO > FieldIOSharedPtr
Definition: FieldIO.h:322
Nektar::LibUtilities::SessionReaderSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: SessionReader.h:119
Nektar::MultiRegions::ExpListSharedPtr
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Definition: LocTraceToTraceMap.h:56
Nektar::SolverUtils::DriverSharedPtr
std::shared_ptr< Driver > DriverSharedPtr
A shared pointer to a Driver object.
Definition: Driver.h:52
Nektar::SolverUtils::GetDriverFactory
DriverFactory & GetDriverFactory()
Definition: Driver.cpp:64
Nektar::SolverUtils::EquationSystemSharedPtr
std::shared_ptr< EquationSystem > EquationSystemSharedPtr
A shared pointer to an EquationSystem object.
Definition: EquationSystem.h:65
Nektar::SolverUtils::GetEquationSystemFactory
EquationSystemFactory & GetEquationSystemFactory()
Definition: EquationSystem.cpp:99
Nektar::SolverUtils::GetForcingFactory
ForcingFactory & GetForcingFactory()
Declaration of the forcing factory singleton.
Definition: Forcing.cpp:41
Nektar::IncNavierStokesSharedPtr
std::shared_ptr< IncNavierStokes > IncNavierStokesSharedPtr
Definition: IncNavierStokes.h:294
Nektar::VWIIterationTypeMap
const std::string VWIIterationTypeMap[]
Definition: VortexWaveInteraction.h:62
Nektar::eFixedWaveForcingWithSubIterationOnAlpha
@ eFixedWaveForcingWithSubIterationOnAlpha
Definition: VortexWaveInteraction.h:58
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Vmath::Vsqrt
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.hpp:340
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.hpp:72
Vmath::Svtvp
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Svtvp (scalar times vector plus vector): z = alpha*x + y.
Definition: Vmath.hpp:396
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.hpp:292
Vmath::Vvtvp
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.hpp:366
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.hpp:180
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.hpp:100
Vmath::Imin
int Imin(int n, const T *x, const int incx)
Return the index of the minimum element in x.
Definition: Vmath.hpp:704
Vmath::Fill
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.hpp:54
Vmath::Sadd
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
Definition: Vmath.hpp:194
Vmath::Vmax
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
Definition: Vmath.hpp:644
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825
Vmath::Vsub
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.hpp:220
std
STL namespace.
tinysimd::abs
scalarT< T > abs(scalarT< T > in)
Definition: scalar.hpp:289
tinysimd::log
scalarT< T > log(scalarT< T > in)
Definition: scalar.hpp:294
tinysimd::sqrt
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:285
Generated on Sun Jun 15 2025 07:52:35 for Nektar++ by doxygen 1.9.4