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
1071 ExecuteWave();
1072
1073 m_wavePressure->GetPlane(0)->BwdTrans(
1074 m_wavePressure->GetPlane(0)->GetCoeffs(),
1075 m_wavePressure->GetPlane(0)->UpdatePhys());
1076 m_wavePressure->GetPlane(1)->BwdTrans(
1077 m_wavePressure->GetPlane(1)->GetCoeffs(),
1078 m_wavePressure->GetPlane(1)->UpdatePhys());
1079
1080 int npts = m_waveVelocities[0]->GetPlane(0)->GetNpoints();
1081 NekDouble Linf;
1082 Array<OneD, NekDouble> val(2 * npts, 0.0);
1083
1084 Vmath::Vmul(npts, m_wavePressure->GetPlane(0)->UpdatePhys(), 1,
1085 m_wavePressure->GetPlane(0)->UpdatePhys(), 1, val, 1);
1086 Vmath::Vvtvp(npts, m_wavePressure->GetPlane(1)->UpdatePhys(), 1,
1087 m_wavePressure->GetPlane(1)->UpdatePhys(), 1, val, 1, val, 1);
1088 cout << "int P^2 " << m_wavePressure->GetPlane(0)->Integral(val) << endl;
1089 Vmath::Vsqrt(npts, val, 1, val, 1);
1090
1091 Linf = Vmath::Vmax(npts, val, 1);
1092 cout << "Linf: " << Linf << endl;
1093
1094 NekDouble l2, norm;
1095 l2 =
1096 m_wavePressure->GetPlane(0)->L2(m_wavePressure->GetPlane(0)->GetPhys());
1097 norm = l2 * l2;
1098 l2 =
1099 m_wavePressure->GetPlane(1)->L2(m_wavePressure->GetPlane(1)->GetPhys());
1100 norm += l2 * l2;
1101
1102 Vmath::Fill(npts, 1.0, val, 1);
1103 NekDouble area = m_waveVelocities[0]->GetPlane(0)->Integral(val);
1104
1105 l2 = sqrt(norm / area);
1106
1107 cout << "L2: " << l2 << endl;
1108
1109 cout << "Ratio Linf/L2: " << Linf / l2 << endl;
1110}
1111
1112void VortexWaveInteraction::SaveFile(string file, string dir, int n)
1113{
1114 static map<string, int> opendir;
1115
1116 if (opendir.find(dir) == opendir.end())
1117 {
1118 // make directory and presume will fail if it already exists
1119 string mkdir = "mkdir " + dir;
1120 ASSERTL0(system(mkdir.c_str()) == 0,
1121 "Failed to make directory '" + dir + "'");
1122
1123 opendir[dir] = 1;
1124 }
1125
1126 string savefile = dir + "/" + file + "." + std::to_string(n);
1127 string syscall = "cp -f " + file + " " + savefile;
1128
1129 ASSERTL0(system(syscall.c_str()) == 0, syscall.c_str());
1130}
1131
1132void VortexWaveInteraction::MoveFile(string file, string dir, int n)
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 = dir + "/" + file + "." + std::to_string(n);
1146 string syscall = "mv -f " + file + " " + savefile;
1147
1148 ASSERTL0(system(syscall.c_str()) == 0, syscall.c_str());
1149}
1150
1151void VortexWaveInteraction::CopyFile(string file1end, string file2end)
1152{
1153 string cpfile1 = m_sessionName + file1end;
1154 string cpfile2 = m_sessionName + file2end;
1155 string syscall = "cp -f " + cpfile1 + " " + cpfile2;
1156
1157 if (system(syscall.c_str()))
1158 {
1159 ASSERTL0(false, syscall.c_str());
1160 }
1161}
1162
1163void VortexWaveInteraction::AppendEvlToFile(string file, int n)
1164{
1165 FILE *fp;
1166 fp = fopen(file.c_str(), "a");
1167 fprintf(fp, "%d: %lf %16.12le %16.12le\n", n, m_alpha[0],
1168 m_leading_real_evl[0], m_leading_imag_evl[0]);
1169 fclose(fp);
1170}
1171
1172void VortexWaveInteraction::AppendEvlToFile(string file, NekDouble WaveForceMag)
1173{
1174 FILE *fp;
1175 fp = fopen(file.c_str(), "a");
1176 fprintf(fp, "%lf %lf %16.12le %16.12le\n", WaveForceMag, m_alpha[0],
1177 m_leading_real_evl[0], m_leading_imag_evl[0]);
1178 fclose(fp);
1179}
1180
1181void VortexWaveInteraction::SaveLoopDetails(std::string SaveDir, int i)
1182
1183{
1184 // Save NS restart file
1185 SaveFile(m_sessionName + ".rst", SaveDir, i + 1);
1186 // Save Streak Solution
1187 SaveFile(m_sessionName + "_streak.fld", SaveDir, i);
1188 // Save Wave solution output
1189 SaveFile(m_sessionName + ".evl", SaveDir, i);
1190 SaveFile(m_sessionName + "_eig_0", SaveDir, i);
1191 // Save new field file of eigenvalue
1192 SaveFile(m_sessionName + ".fld", SaveDir, i);
1193 if (!(m_sessionVWI->DefinesSolverInfo("INTERFACE")))
1194 {
1195 // Save new forcing file
1196 SaveFile(m_sessionName + ".vwi", SaveDir, i + 1);
1197 }
1198}
1199
1200void VortexWaveInteraction::ExecuteLoop(bool CalcWaveForce)
1201{
1202
1203 // the global info has to be written in the
1204 // roll session file to use the interface loop
1205 if (m_sessionRoll->DefinesSolverInfo("INTERFACE"))
1206 {
1207 static int cnt = 0;
1208 bool skiprollstreak = false;
1209 if (cnt == 0 && m_sessionVWI->GetParameter("rollstreakfromit") == 1)
1210 {
1211 skiprollstreak = true;
1212 cout << "skip roll-streak at the first iteration" << endl;
1213 }
1214
1215 if (skiprollstreak != true)
1216 {
1217
1218 LibUtilities::NekManager<
1219 MultiRegions::GlobalLinSysKey,
1220 MultiRegions::GlobalLinSys>::EnableManagement("GlobalLinSys");
1221 ExecuteRoll();
1222 LibUtilities::NekManager<
1223 MultiRegions::GlobalLinSysKey,
1224 MultiRegions::GlobalLinSys>::DisableManagement("GlobalLinSys");
1225#ifndef _WIN32
1226 sleep(3);
1227#endif
1228 ExecuteStreak();
1229#ifndef _WIN32
1230 sleep(3);
1231#endif
1232 }
1233
1234 string syscall;
1235 string movedmesh = m_sessionName + "_advPost_moved.xml";
1236 string movedinterpmesh = m_sessionName + "_interp_moved.xml";
1237 // rewrite the Rollsessionfile (we start from the waleffe forcing)
1238 // string meshbndjumps = m_sessionName +"_bndjumps.xml";
1239 // if(cnt==0)
1240 //{
1241 // take the conditions tag from meshbndjumps and copy into
1242 // the rolls session file
1243 //}
1244
1245 string c = std::to_string(cnt);
1246
1247 // save old roll solution
1248 string oldroll = m_sessionName + "_roll_" + c + ".fld";
1249 syscall = "cp -f " + m_sessionName + "-Base.fld" + " " + oldroll;
1250 cout << syscall.c_str() << endl;
1251 if (system(syscall.c_str()))
1252 {
1253 ASSERTL0(false, syscall.c_str());
1254 }
1255 // define file names
1256 string filePost = m_sessionName + "_advPost.xml";
1257 string filestreak = m_sessionName + "_streak.fld";
1258 string filewave = m_sessionName + "_wave.fld";
1259 string filewavepressure = m_sessionName + "_wave_p_split.fld";
1260 string fileinterp = m_sessionName + "_interp.xml";
1261 string interpstreak = m_sessionName + "_interpstreak_" + c + ".fld";
1262 string interwavepressure =
1263 m_sessionName + "_wave_p_split_interp_" + c + ".fld";
1264 string c_alpha = std::to_string(m_alpha[0]);
1265 cout << "alpha = " << m_alpha[0] << endl;
1266
1267 if (m_sessionVWI->GetSolverInfo("INTERFACE") != "phase")
1268 {
1269 cout << "zerophase" << endl;
1270
1271 syscall = "../../utilities/PostProcessing/Extras/MoveMesh " +
1272 filePost + " " + filestreak + " " + fileinterp + " " +
1273 c_alpha;
1274
1275 cout << syscall.c_str() << endl;
1276 if (system(syscall.c_str()))
1277 {
1278 ASSERTL0(false, syscall.c_str());
1279 }
1280
1281 // move the advPost mesh (remark update alpha!!!)
1282 syscall = "../../utilities/PostProcessing/Extras/MoveMesh " +
1283 filePost + " " + filestreak + " " + filePost + " " +
1284 c_alpha;
1285 cout << syscall.c_str() << endl;
1286 if (system(syscall.c_str()))
1287 {
1288 ASSERTL0(false, syscall.c_str());
1289 }
1290
1291 // save oldstreak
1292 string oldstreak = m_sessionName + "_streak_" + c + ".fld";
1293 syscall = "cp -f " + filestreak + " " + oldstreak;
1294 cout << syscall.c_str() << endl;
1295 if (system(syscall.c_str()))
1296 {
1297 ASSERTL0(false, syscall.c_str());
1298 }
1299
1300 // interpolate the streak field into the new mesh
1301 string movedmesh = m_sessionName + "_advPost_moved.xml";
1302 string movedinterpmesh = m_sessionName + "_interp_moved.xml";
1303
1304 // create the interp streak
1305
1306 syscall = "../../utilities/PostProcessing/Extras/FieldToField " +
1307 fileinterp + " " + filestreak + " " + movedinterpmesh +
1308 " " + interpstreak;
1309
1310 cout << syscall.c_str() << endl;
1311 if (system(syscall.c_str()))
1312 {
1313 ASSERTL0(false, syscall.c_str());
1314 }
1315
1316 // save the old mesh
1317 string meshfile = m_sessionName + ".xml";
1318 string meshold = m_sessionName + "_" + c + ".xml";
1319 syscall = "cp -f " + meshfile + " " + meshold;
1320 cout << syscall.c_str() << endl;
1321 if (system(syscall.c_str()))
1322 {
1323 ASSERTL0(false, syscall.c_str());
1324 }
1325
1326 // overwriting the meshfile with the new mesh
1327 syscall = "cp -f " + movedmesh + " " + meshfile;
1328 cout << syscall.c_str() << endl;
1329 if (system(syscall.c_str()))
1330 {
1331 ASSERTL0(false, syscall.c_str());
1332 }
1333
1334 // overwriting the streak file!!
1335 syscall = "cp -f " + interpstreak + " " + filestreak;
1336 cout << syscall.c_str() << endl;
1337 if (system(syscall.c_str()))
1338 {
1339 ASSERTL0(false, syscall.c_str());
1340 }
1341
1342 // calculate the wave
1343 ExecuteWave();
1344
1345 // save the wave field:
1346 string oldwave = m_sessionName + "_wave_" + c + ".fld";
1347 syscall = "cp -f " + m_sessionName + ".fld" + " " + oldwave;
1348 cout << syscall.c_str() << endl;
1349 if (system(syscall.c_str()))
1350 {
1351 ASSERTL0(false, syscall.c_str());
1352 }
1353
1354 // save old jump conditions:
1355 string ujump = m_sessionName + "_u_5.bc";
1356 syscall = "cp -f " + ujump + " " + m_sessionName + "_u_5.bc_" + c;
1357 cout << syscall.c_str() << endl;
1358 if (system(syscall.c_str()))
1359 {
1360 ASSERTL0(false, syscall.c_str());
1361 }
1362
1363 string vjump = m_sessionName + "_v_5.bc";
1364 syscall = "cp -f " + vjump + " " + m_sessionName + "_v_5.bc_" + c;
1365 cout << syscall.c_str() << endl;
1366 if (system(syscall.c_str()))
1367 {
1368 ASSERTL0(false, syscall.c_str());
1369 }
1370 cnt++;
1371 c = std::to_string(cnt);
1372
1373 // use relaxation
1374 if (GetVWIIterationType() !=
1375 eFixedWaveForcingWithSubIterationOnAlpha)
1376 {
1377 // the critical layer should be the bnd region 3
1378 // int reg =3;
1379 // FileRelaxation(reg);
1380 }
1381 // calculate the jump conditions
1382 string wavefile = m_sessionName + ".fld";
1383 syscall = "../../utilities/PostProcessing/Extras/FldCalcBCs " +
1384 movedmesh + " " + wavefile + " " + interpstreak +
1385 "> data" + c;
1386 cout << syscall.c_str() << endl;
1387 if (system(syscall.c_str()))
1388 {
1389 ASSERTL0(false, syscall.c_str());
1390 }
1391
1392 // move the new name_interp_moved.xml into name_interp.xml
1393 syscall = "cp -f " + movedinterpmesh + " " + fileinterp;
1394 cout << syscall.c_str() << endl;
1395 if (system(syscall.c_str()))
1396 {
1397 ASSERTL0(false, syscall.c_str());
1398 }
1399 // move the new name_advPost_moved.xml into name_advPost.xml
1400 syscall = "cp -f " + movedmesh + " " + filePost;
1401 cout << syscall.c_str() << endl;
1402 if (system(syscall.c_str()))
1403 {
1404 ASSERTL0(false, syscall.c_str());
1405 }
1406 }
1407 else if (m_sessionVWI->GetSolverInfo("INTERFACE") == "phase")
1408 {
1409 cout << "phase" << endl;
1410 // determine cr:
1411 NekDouble cr;
1412 string cr_str;
1413 stringstream st;
1414
1415 // calculate the wave
1416 ExecuteWave();
1417
1418 // save oldstreak
1419 string oldstreak = m_sessionName + "_streak_" + c + ".fld";
1420 syscall = "cp -f " + filestreak + " " + oldstreak;
1421 cout << syscall.c_str() << endl;
1422 if (system(syscall.c_str()))
1423 {
1424 ASSERTL0(false, syscall.c_str());
1425 }
1426
1427 // save wave
1428 syscall = "cp -f " + m_sessionName + ".fld" + " " + filewave;
1429 cout << syscall.c_str() << endl;
1430 if (system(syscall.c_str()))
1431 {
1432 ASSERTL0(false, syscall.c_str());
1433 }
1434 // save the old mesh
1435 string meshfile = m_sessionName + ".xml";
1436 string meshold = m_sessionName + "_" + c + ".xml";
1437 syscall = "cp -f " + meshfile + " " + meshold;
1438 cout << syscall.c_str() << endl;
1439 if (system(syscall.c_str()))
1440 {
1441 ASSERTL0(false, syscall.c_str());
1442 }
1443
1444 // save the oldwave field:
1445 string oldwave = m_sessionName + "_wave_" + c + ".fld";
1446 syscall = "cp -f " + m_sessionName + ".fld" + " " + oldwave;
1447 cout << syscall.c_str() << endl;
1448 if (system(syscall.c_str()))
1449 {
1450 ASSERTL0(false, syscall.c_str());
1451 }
1452
1453 // save old jump conditions:
1454 string ujump = m_sessionName + "_u_5.bc";
1455 syscall = "cp -f " + ujump + " " + m_sessionName + "_u_5.bc_" + c;
1456 cout << syscall.c_str() << endl;
1457 if (system(syscall.c_str()))
1458 {
1459 ASSERTL0(false, syscall.c_str());
1460 }
1461
1462 string vjump = m_sessionName + "_v_5.bc";
1463 syscall = "cp -f " + vjump + " " + m_sessionName + "_v_5.bc_" + c;
1464 cout << syscall.c_str() << endl;
1465 if (system(syscall.c_str()))
1466 {
1467 ASSERTL0(false, syscall.c_str());
1468 }
1469 cnt++;
1470
1471 cr = m_leading_imag_evl[0] / m_alpha[0];
1472 st << cr;
1473 cr_str = st.str();
1474 cout << "cr=" << cr_str << endl;
1475 // NB -g or NOT!!!
1476 // move the mesh around the critical layer
1477 syscall = "../../utilities/PostProcessing/Extras/MoveMesh " +
1478 filePost + " " + filestreak + " " + fileinterp + " " +
1479 c_alpha + " " + cr_str;
1480
1481 cout << syscall.c_str() << endl;
1482 if (system(syscall.c_str()))
1483 {
1484 ASSERTL0(false, syscall.c_str());
1485 }
1486 // NB -g or NOT!!!
1487 // move the advPost mesh (remark update alpha!!!)
1488 syscall = "../../utilities/PostProcessing/Extras/MoveMesh " +
1489 filePost + " " + filestreak + " " + filePost + " " +
1490 c_alpha + " " + cr_str;
1491 cout << syscall.c_str() << endl;
1492 if (system(syscall.c_str()))
1493 {
1494 ASSERTL0(false, syscall.c_str());
1495 }
1496
1497 // interp streak into the new mesh
1498 syscall = "../../utilities/PostProcessing/Extras/FieldToField " +
1499 fileinterp + " " + filestreak + " " + movedinterpmesh +
1500 " " + interpstreak;
1501
1502 cout << syscall.c_str() << endl;
1503 if (system(syscall.c_str()))
1504 {
1505 ASSERTL0(false, syscall.c_str());
1506 }
1507
1508 // split wave sol
1509 syscall = "../../utilities/PostProcessing/Extras/SplitFld " +
1510 filePost + " " + filewave;
1511
1512 cout << syscall.c_str() << endl;
1513 if (system(syscall.c_str()))
1514 {
1515 ASSERTL0(false, syscall.c_str());
1516 }
1517 // interp wave
1518 syscall = "../../utilities/PostProcessing/Extras/FieldToField " +
1519 filePost + " " + filewavepressure + " " + movedmesh +
1520 " " + interwavepressure;
1521
1522 cout << syscall.c_str() << endl;
1523 if (system(syscall.c_str()))
1524 {
1525 ASSERTL0(false, syscall.c_str());
1526 }
1527
1528 // use relaxation
1529 if (GetVWIIterationType() !=
1530 eFixedWaveForcingWithSubIterationOnAlpha)
1531 {
1532 // the critical layer should be the bnd region 3
1533 // int reg =3;
1534 // FileRelaxation(reg);
1535 }
1536
1537 // cp wavepressure to m_sessionName.fld(to get
1538 // the right bcs names using FldCalcBCs)
1539 syscall =
1540 "cp -f " + interwavepressure + " " + m_sessionName + ".fld";
1541 cout << syscall.c_str() << endl;
1542 if (system(syscall.c_str()))
1543 {
1544 ASSERTL0(false, syscall.c_str());
1545 }
1546
1547 // calculate the jump conditions
1548 // NB -g or NOT!!!
1549 syscall = "../../utilities/PostProcessing/Extras/FldCalcBCs " +
1550 movedmesh + " " + m_sessionName + ".fld" + " " +
1551 interpstreak + "> data" + c;
1552 cout << syscall.c_str() << endl;
1553 if (system(syscall.c_str()))
1554 {
1555 ASSERTL0(false, syscall.c_str());
1556 }
1557
1558 // overwriting the meshfile with the new mesh
1559 syscall = "cp -f " + movedmesh + " " + meshfile;
1560 cout << syscall.c_str() << endl;
1561 if (system(syscall.c_str()))
1562 {
1563 ASSERTL0(false, syscall.c_str());
1564 }
1565
1566 // overwriting the streak file!! (maybe is useless)
1567 syscall = "cp -f " + interpstreak + " " + filestreak;
1568 cout << syscall.c_str() << endl;
1569 if (system(syscall.c_str()))
1570 {
1571 ASSERTL0(false, syscall.c_str());
1572 }
1573 // move the new name_interp_moved.xml into name_interp.xml
1574 syscall = "cp -f " + movedinterpmesh + " " + fileinterp;
1575 cout << syscall.c_str() << endl;
1576 if (system(syscall.c_str()))
1577 {
1578 ASSERTL0(false, syscall.c_str());
1579 }
1580 // move the new name_advPost_moved.xml into name_advPost.xml
1581 syscall = "cp -f " + movedmesh + " " + filePost;
1582 cout << syscall.c_str() << endl;
1583 if (system(syscall.c_str()))
1584 {
1585 ASSERTL0(false, syscall.c_str());
1586 }
1587 }
1588 }
1589 else
1590 {
1591 LibUtilities::NekManager<
1592 MultiRegions::GlobalLinSysKey,
1593 MultiRegions::GlobalLinSys>::EnableManagement("GlobalLinSys");
1594 ExecuteRoll();
1595 LibUtilities::NekManager<
1596 MultiRegions::GlobalLinSysKey,
1597 MultiRegions::GlobalLinSys>::DisableManagement("GlobalLinSys");
1598
1599#ifndef _WIN32
1600 sleep(3);
1601#endif
1602 ExecuteStreak();
1603#ifndef _WIN32
1604 sleep(3);
1605#endif
1606
1607 if (m_moveMeshToCriticalLayer)
1608 {
1609 string syscall;
1610 char alpchar[16] = "";
1611 snprintf(alpchar, 16, "%f", m_alpha[0]);
1612
1613 string filePost = m_sessionName + "_advPost.xml";
1614 string filestreak = m_sessionName + "_streak.fld";
1615 string filewave = m_sessionName + "_wave.fld";
1616 string filewavepressure = m_sessionName + "_wave_p_split.fld";
1617 string fileinterp = m_sessionName + "_interp.xml";
1618 string interpstreak = m_sessionName + "_interpstreak.fld";
1619 string interwavepressure =
1620 m_sessionName + "_wave_p_split_interp.fld";
1621 syscall = "../../utilities/PostProcessing/Extras/MoveMesh " +
1622 filePost + " " + filestreak + " " + fileinterp + " " +
1623 alpchar;
1624
1625 cout << syscall.c_str() << endl;
1626 if (system(syscall.c_str()))
1627 {
1628 ASSERTL0(false, syscall.c_str());
1629 }
1630
1631 // move the advPost mesh (remark update alpha!!!)
1632 syscall = "../../utilities/PostProcessing/Extras/MoveMesh " +
1633 filePost + " " + filestreak + " " + filePost + " " +
1634 alpchar;
1635 cout << syscall.c_str() << endl;
1636 if (system(syscall.c_str()))
1637 {
1638 ASSERTL0(false, syscall.c_str());
1639 }
1640
1641 // save oldstreak
1642 string oldstreak = m_sessionName + "_streak.fld";
1643 syscall = "cp -f " + filestreak + " " + oldstreak;
1644 cout << syscall.c_str() << endl;
1645 if (system(syscall.c_str()))
1646 {
1647 ASSERTL0(false, syscall.c_str());
1648 }
1649
1650 // interpolate the streak field into the new mesh
1651 string movedmesh = m_sessionName + "_advPost_moved.xml";
1652 string movedinterpmesh = m_sessionName + "_interp_moved.xml";
1653
1654 // create the interp streak
1655
1656 syscall = "../../utilities/PostProcessing/Extras/FieldToField " +
1657 fileinterp + " " + filestreak + " " + movedinterpmesh +
1658 " " + interpstreak;
1659
1660 cout << syscall.c_str() << endl;
1661 if (system(syscall.c_str()))
1662 {
1663 ASSERTL0(false, syscall.c_str());
1664 }
1665
1666 // save the old mesh
1667 string meshfile = m_sessionName + ".xml";
1668 string meshold = m_sessionName + ".xml";
1669 syscall = "cp -f " + meshfile + " " + meshold;
1670 cout << syscall.c_str() << endl;
1671 if (system(syscall.c_str()))
1672 {
1673 ASSERTL0(false, syscall.c_str());
1674 }
1675
1676 // overwriting the meshfile with the new mesh
1677 syscall = "cp -f " + movedmesh + " " + meshfile;
1678 cout << syscall.c_str() << endl;
1679 if (system(syscall.c_str()))
1680 {
1681 ASSERTL0(false, syscall.c_str());
1682 }
1683
1684 // overwriting the streak file!!
1685 syscall = "cp -f " + interpstreak + " " + filestreak;
1686 cout << syscall.c_str() << endl;
1687 if (system(syscall.c_str()))
1688 {
1689 ASSERTL0(false, syscall.c_str());
1690 }
1691 }
1692
1693 ExecuteWave();
1694
1695 if (CalcWaveForce)
1696 {
1697 CalcNonLinearWaveForce();
1698 }
1699 }
1700}
1701
1702bool VortexWaveInteraction::CheckEigIsStationary(bool reset)
1703{
1704 static NekDouble previous_real_evl = -1.0;
1705 static NekDouble previous_imag_evl = -1.0;
1706 static int min_iter = 0;
1707
1708 if (reset)
1709 {
1710 previous_real_evl = -1.0;
1711 min_iter = 0;
1712 }
1713
1714 if (previous_real_evl == -1.0 || min_iter < m_minInnerIterations)
1715 {
1716 previous_real_evl = m_leading_real_evl[0];
1717 previous_imag_evl = m_leading_imag_evl[0];
1718 min_iter++;
1719 return false;
1720 }
1721
1722 cout << "Growth tolerance: "
1723 << fabs((m_leading_real_evl[0] - previous_real_evl) /
1724 m_leading_real_evl[0])
1725 << endl;
1726 cout << "Phase tolerance: "
1727 << fabs((m_leading_imag_evl[0] - previous_imag_evl) /
1728 m_leading_imag_evl[0])
1729 << endl;
1730
1731 // See if real and imaginary growth have converged to with m_eigRelTol
1732 if ((fabs((m_leading_real_evl[0] - previous_real_evl) /
1733 m_leading_real_evl[0]) < m_eigRelTol) ||
1734 (fabs(m_leading_real_evl[0]) < 1e-6))
1735 {
1736 previous_real_evl = m_leading_real_evl[0];
1737 previous_imag_evl = m_leading_imag_evl[0];
1738 if ((fabs((m_leading_imag_evl[0] - previous_imag_evl) /
1739 m_leading_imag_evl[0]) < m_eigRelTol) ||
1740 (fabs(m_leading_imag_evl[0]) < 1e-6))
1741 {
1742 return true;
1743 }
1744 }
1745 else
1746 {
1747 if (fabs(m_leading_imag_evl[0]) > 1e-2)
1748 {
1749 cout << "Warning: imaginary eigenvalue is greater than 1e-2"
1750 << endl;
1751 }
1752 previous_real_evl = m_leading_real_evl[0];
1753 previous_imag_evl = m_leading_imag_evl[0];
1754 return false;
1755 }
1756 return false;
1757}
1758
1759// Check to see if leading eigenvalue is within tolerance defined
1760// in m_neutralPointTol
1761bool VortexWaveInteraction::CheckIfAtNeutralPoint(void)
1762{
1763 bool returnval = false;
1764 if (m_sessionRoll->DefinesSolverInfo("INTERFACE"))
1765 {
1766 if (m_sessionVWI->GetSolverInfo("INTERFACE") == "phase")
1767 {
1768 if (abs(m_leading_real_evl[0]) < 1e-4)
1769 {
1770 returnval = true;
1771 }
1772 }
1773 else
1774 {
1775 if (abs(m_leading_real_evl[0]) < 1e-4 &&
1776 abs(m_leading_imag_evl[0]) < 2e-6)
1777 {
1778 returnval = true;
1779 }
1780 }
1781 }
1782 else
1783 {
1784 if ((m_leading_real_evl[0] * m_leading_real_evl[0] +
1785 m_leading_imag_evl[0] * m_leading_imag_evl[0]) <
1786 m_neutralPointTol * m_neutralPointTol)
1787 {
1788 returnval = true;
1789 }
1790 }
1791 if (fabs(m_leading_imag_evl[0]) > 1e-2)
1792 {
1793 cout << "Warning: imaginary eigenvalue is greater than 1e-2" << endl;
1794 }
1795
1796 return returnval;
1797}
1798
1799// Similar routine to UpdateAlpha
1800
1801void VortexWaveInteraction::UpdateWaveForceMag(int outeriter)
1802{
1803 NekDouble wavef_new = 0.0;
1804
1805 if (outeriter == 1)
1806 {
1807 m_waveForceMag[1] = m_waveForceMag[0];
1808 if (m_leading_real_evl[0] > 0.0)
1809 {
1810 wavef_new = m_waveForceMag[0] - m_waveForceMagStep;
1811 }
1812 else
1813 {
1814 wavef_new = m_waveForceMag[0] + m_waveForceMagStep;
1815 }
1816 }
1817 else
1818 {
1819 int i;
1820 int nstore = (m_waveForceMag.size() < outeriter) ? m_waveForceMag.size()
1821 : outeriter;
1822 Array<OneD, NekDouble> WaveForce(nstore);
1823 Array<OneD, NekDouble> Growth(nstore);
1824
1825 Vmath::Vcopy(nstore, m_waveForceMag, 1, WaveForce, 1);
1826 Vmath::Vcopy(nstore, m_leading_real_evl, 1, Growth, 1);
1827
1828 // Sort WaveForce Growth values;
1829 NekDouble store;
1830 int k;
1831 for (i = 0; i < nstore; ++i)
1832 {
1833 k = Vmath::Imin(nstore - i, &WaveForce[i], 1);
1834
1835 store = WaveForce[i];
1836 WaveForce[i] = WaveForce[i + k];
1837 WaveForce[i + k] = store;
1838
1839 store = Growth[i];
1840 Growth[i] = Growth[i + k];
1841 Growth[i + k] = store;
1842 }
1843
1844 // See if we have any values that cross zero
1845 for (i = 0; i < nstore - 1; ++i)
1846 {
1847 if (Growth[i] * Growth[i + 1] < 0.0)
1848 {
1849 break;
1850 }
1851 }
1852
1853 if (i != nstore - 1)
1854 {
1855 if (nstore == 2)
1856 {
1857 wavef_new =
1858 (WaveForce[0] * Growth[1] - WaveForce[1] * Growth[0]) /
1859 (Growth[1] - Growth[0]);
1860 }
1861 else
1862 {
1863 // use a quadratic fit and step through 10000 points
1864 // to find zero.
1865 int j;
1866 int nsteps = 10000;
1867 int idx = (i == 0) ? 1 : i;
1868 NekDouble da = WaveForce[idx + 1] - WaveForce[idx - 1];
1869 NekDouble gval_m1 = Growth[idx - 1], a, gval;
1870 NekDouble c1 = Growth[idx - 1] /
1871 (WaveForce[idx - 1] - WaveForce[idx]) /
1872 (WaveForce[idx - 1] - WaveForce[idx + 1]);
1873 NekDouble c2 = Growth[idx] /
1874 (WaveForce[idx] - WaveForce[idx - 1]) /
1875 (WaveForce[idx] - WaveForce[idx + 1]);
1876 NekDouble c3 = Growth[idx + 1] /
1877 (WaveForce[idx + 1] - WaveForce[idx - 1]) /
1878 (WaveForce[idx + 1] - WaveForce[idx]);
1879
1880 for (j = 1; j < nsteps + 1; ++j)
1881 {
1882 a = WaveForce[i] + j * da / (NekDouble)nsteps;
1883 gval =
1884 c1 * (a - WaveForce[idx]) * (a - WaveForce[idx + 1]) +
1885 c2 * (a - WaveForce[idx - 1]) *
1886 (a - WaveForce[idx + 1]) +
1887 c3 * (a - WaveForce[idx - 1]) * (a - WaveForce[idx]);
1888
1889 if (gval * gval_m1 < 0.0)
1890 {
1891 wavef_new = ((a + da / (NekDouble)nsteps) * gval -
1892 a * gval_m1) /
1893 (gval - gval_m1);
1894 break;
1895 }
1896 }
1897 }
1898 }
1899 else // step backward/forward
1900 {
1901 if (Growth[i] > 0.0)
1902 {
1903 wavef_new = m_waveForceMag[0] - m_waveForceMagStep;
1904 }
1905 else
1906 {
1907 wavef_new = m_waveForceMag[0] + m_waveForceMagStep;
1908 }
1909 }
1910 }
1911
1912 for (int i = m_waveForceMag.size() - 1; i > 0; --i)
1913 {
1914 m_waveForceMag[i] = m_waveForceMag[i - 1];
1915 m_leading_real_evl[i] = m_leading_real_evl[i - 1];
1916 m_leading_imag_evl[i] = m_leading_imag_evl[i - 1];
1917 }
1918
1919 m_waveForceMag[0] = wavef_new;
1920}
1921
1922void VortexWaveInteraction::UpdateDAlphaDWaveForceMag(NekDouble alpha_init)
1923{
1924 m_dAlphaDWaveForceMag = (m_alpha[0] - alpha_init) / m_waveForceMagStep;
1925}
1926
1927void VortexWaveInteraction::UpdateAlpha(int outeriter)
1928{
1929 NekDouble alp_new = 0.0;
1930
1931 if (outeriter == 1)
1932 {
1933 m_alpha[1] = m_alpha[0];
1934 if (m_leading_real_evl[0] > 0.0)
1935 {
1936 alp_new = m_alpha[0] + m_alphaStep;
1937 }
1938 else
1939 {
1940 alp_new = m_alpha[0] - m_alphaStep;
1941 }
1942 }
1943 else
1944 {
1945 int i;
1946 int nstore = (m_alpha.size() < outeriter) ? m_alpha.size() : outeriter;
1947 Array<OneD, NekDouble> Alpha(nstore);
1948 Array<OneD, NekDouble> Growth(nstore);
1949
1950 Vmath::Vcopy(nstore, m_alpha, 1, Alpha, 1);
1951 Vmath::Vcopy(nstore, m_leading_real_evl, 1, Growth, 1);
1952
1953 // Sort Alpha Growth values;
1954 NekDouble store;
1955 int k;
1956 for (i = 0; i < nstore; ++i)
1957 {
1958 k = Vmath::Imin(nstore - i, &Alpha[i], 1);
1959
1960 store = Alpha[i];
1961 Alpha[i] = Alpha[i + k];
1962 Alpha[i + k] = store;
1963
1964 store = Growth[i];
1965 Growth[i] = Growth[i + k];
1966 Growth[i + k] = store;
1967 }
1968
1969 // See if we have any values that cross zero
1970 for (i = 0; i < nstore - 1; ++i)
1971 {
1972 if (Growth[i] * Growth[i + 1] < 0.0)
1973 {
1974 break;
1975 }
1976 }
1977
1978 if (i != nstore - 1)
1979 {
1980 if (nstore == 2)
1981 {
1982 alp_new = (Alpha[0] * Growth[1] - Alpha[1] * Growth[0]) /
1983 (Growth[1] - Growth[0]);
1984 }
1985 else
1986 {
1987 // use a quadratic fit and step through 10000 points
1988 // to find zero.
1989 int j;
1990 int nsteps = 10000;
1991 int idx = (i == 0) ? 1 : i;
1992 NekDouble da = Alpha[idx + 1] - Alpha[idx - 1];
1993 NekDouble gval_m1 = Growth[idx - 1], a, gval;
1994 NekDouble c1 = Growth[idx - 1] / (Alpha[idx - 1] - Alpha[idx]) /
1995 (Alpha[idx - 1] - Alpha[idx + 1]);
1996 NekDouble c2 = Growth[idx] / (Alpha[idx] - Alpha[idx - 1]) /
1997 (Alpha[idx] - Alpha[idx + 1]);
1998 NekDouble c3 = Growth[idx + 1] /
1999 (Alpha[idx + 1] - Alpha[idx - 1]) /
2000 (Alpha[idx + 1] - Alpha[idx]);
2001
2002 for (j = 1; j < nsteps + 1; ++j)
2003 {
2004 a = Alpha[i] + j * da / (NekDouble)nsteps;
2005 gval = c1 * (a - Alpha[idx]) * (a - Alpha[idx + 1]) +
2006 c2 * (a - Alpha[idx - 1]) * (a - Alpha[idx + 1]) +
2007 c3 * (a - Alpha[idx - 1]) * (a - Alpha[idx]);
2008
2009 if (gval * gval_m1 < 0.0)
2010 {
2011 alp_new = ((a + da / (NekDouble)nsteps) * gval -
2012 a * gval_m1) /
2013 (gval - gval_m1);
2014 break;
2015 }
2016 }
2017 }
2018 }
2019 else // step backward/forward
2020 {
2021 if (Growth[i] > 0.0)
2022 {
2023 alp_new = m_alpha[0] + m_alphaStep;
2024 }
2025 else
2026 {
2027 alp_new = m_alpha[0] - m_alphaStep;
2028 }
2029 }
2030 }
2031
2032 for (int i = m_alpha.size() - 1; i > 0; --i)
2033 {
2034 m_alpha[i] = m_alpha[i - 1];
2035 m_leading_real_evl[i] = m_leading_real_evl[i - 1];
2036 m_leading_imag_evl[i] = m_leading_imag_evl[i - 1];
2037 }
2038
2039 m_alpha[0] = alp_new;
2040}
2041
2042Array<OneD, int> VortexWaveInteraction::GetReflectionIndex(void)
2043{
2044 int i, j;
2045 int npts = m_waveVelocities[0]->GetPlane(0)->GetNpoints();
2046 int nel = m_waveVelocities[0]->GetNumElmts();
2047 Array<OneD, int> index(npts);
2048
2049 Array<OneD, NekDouble> coord(2);
2050 Array<OneD, NekDouble> coord_x(npts);
2051 Array<OneD, NekDouble> coord_y(npts);
2052
2053 //-> Dermine the point which is on coordinate (x -> -x + Lx/2, y-> -y)
2054 m_waveVelocities[0]->GetPlane(0)->GetCoords(coord_x, coord_y);
2055 NekDouble xmax = Vmath::Vmax(npts, coord_x, 1);
2056 NekDouble tol = 1e-5;
2057 NekDouble xnew, ynew;
2058
2059 int start = npts - 1;
2060 int e_npts;
2061
2062 bool useOnlyQuads = false;
2063 if (m_sessionVWI->DefinesSolverInfo("SymmetriseOnlyQuads"))
2064 {
2065 useOnlyQuads = true;
2066 }
2067
2068 int cnt;
2069 for (int e = 0; e < nel; ++e)
2070 {
2071 e_npts = m_waveVelocities[0]->GetExp(e)->GetTotPoints();
2072 cnt = m_waveVelocities[0]->GetPhys_Offset(e);
2073
2074 if (useOnlyQuads)
2075 {
2076 if (m_waveVelocities[0]->GetExp(e)->DetShapeType() ==
2077 LibUtilities::eTriangle)
2078 {
2079 for (i = 0; i < e_npts; ++i)
2080 {
2081 index[cnt + i] = -1;
2082 }
2083 continue;
2084 }
2085 }
2086
2087 for (i = cnt; i < cnt + e_npts; ++i)
2088 {
2089 xnew = -coord_x[i] + xmax;
2090 ynew = -coord_y[i];
2091
2092 for (j = start; j >= 0; --j)
2093 {
2094 if ((coord_x[j] - xnew) * (coord_x[j] - xnew) +
2095 (coord_y[j] - ynew) * (coord_y[j] - ynew) <
2096 tol)
2097 {
2098 index[i] = j;
2099 start = j;
2100 break;
2101 }
2102 }
2103
2104 if (j == -1)
2105 {
2106
2107 for (j = npts - 1; j > start; --j)
2108 {
2109
2110 if ((coord_x[j] - xnew) * (coord_x[j] - xnew) +
2111 (coord_y[j] - ynew) * (coord_y[j] - ynew) <
2112 tol)
2113 {
2114 index[i] = j;
2115 break;
2116 }
2117 }
2118 ASSERTL0(j != start, "Failed to find matching point");
2119 }
2120 }
2121 }
2122 return index;
2123}
2124
2125void VortexWaveInteraction::FileRelaxation(int reg)
2126{
2127 cout << "relaxation..." << endl;
2128 static int cnt = 0;
2129 Array<OneD, MultiRegions::ExpListSharedPtr> Iexp =
2130 m_rollField[0]->GetBndCondExpansions();
2131 // cast to 1D explist (otherwise appenddata doesn't work)
2132 MultiRegions::ExpListSharedPtr Ilayer;
2133 Ilayer = MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(
2134 *std::static_pointer_cast<MultiRegions::ExpList>(Iexp[reg]));
2135 int nq = Ilayer->GetTotPoints();
2136 if (cnt == 0)
2137 {
2138 m_bcsForcing = Array<OneD, Array<OneD, NekDouble>>(4);
2139 m_bcsForcing[0] = Array<OneD, NekDouble>(4 * nq);
2140 for (int i = 1; i < 4; ++i)
2141 {
2142 m_bcsForcing[i] = m_bcsForcing[i - 1] + nq;
2143 }
2144 }
2145
2146 // Read in mesh from input file
2147
2148 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef_u;
2149 std::vector<std::vector<NekDouble>> FieldData_u;
2150 string file = m_sessionName;
2151
2152 file += "_u_5.bc";
2153 LibUtilities::FieldIOSharedPtr fld =
2154 LibUtilities::FieldIO::CreateDefault(m_sessionVWI);
2155 fld->Import(file, FieldDef_u, FieldData_u);
2156 Ilayer->ExtractDataToCoeffs(FieldDef_u[0], FieldData_u[0],
2157 FieldDef_u[0]->m_fields[0],
2158 Ilayer->UpdateCoeffs());
2159 Ilayer->BwdTrans(Ilayer->GetCoeffs(), Ilayer->UpdatePhys());
2160
2161 if (cnt == 0)
2162 {
2163 Vmath::Vcopy(nq, Ilayer->UpdatePhys(), 1, m_bcsForcing[2], 1);
2164 }
2165 Vmath::Vcopy(nq, Ilayer->UpdatePhys(), 1, m_bcsForcing[0], 1);
2166
2167 // case relaxation==0 an additional array is needed
2168 Array<OneD, NekDouble> tmp_forcing(nq, 0.0);
2169
2170 if (cnt != 0)
2171 {
2172 if (m_vwiRelaxation == 1.0)
2173 {
2174 Vmath::Vcopy(nq, m_bcsForcing[0], 1, tmp_forcing, 1);
2175 }
2176 Vmath::Smul(nq, 1.0 - m_vwiRelaxation, m_bcsForcing[0], 1,
2177 m_bcsForcing[0], 1);
2178 Vmath::Svtvp(nq, m_vwiRelaxation, m_bcsForcing[2], 1, m_bcsForcing[0],
2179 1, Ilayer->UpdatePhys(), 1);
2180 // generate again the bcs files:
2181
2182 Array<OneD, Array<OneD, NekDouble>> fieldcoeffs(1);
2183 Ilayer->FwdTransLocalElmt(Ilayer->GetPhys(), Ilayer->UpdateCoeffs());
2184 fieldcoeffs[0] = Ilayer->UpdateCoeffs();
2185 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef1 =
2186 Ilayer->GetFieldDefinitions();
2187 std::vector<std::vector<NekDouble>> FieldData_1(FieldDef1.size());
2188 ;
2189 FieldDef1[0]->m_fields.push_back("u");
2190 Ilayer->AppendFieldData(FieldDef1[0], FieldData_1[0]);
2191 fld->Write(file, FieldDef1, FieldData_1);
2192 // save the bcs for the next iteration
2193 if (m_vwiRelaxation != 1.0)
2194 {
2195 Vmath::Smul(nq, 1. / (1.0 - m_vwiRelaxation), m_bcsForcing[0], 1,
2196 m_bcsForcing[0], 1);
2197 Vmath::Vcopy(nq, m_bcsForcing[0], 1, m_bcsForcing[2], 1);
2198 }
2199 else
2200 {
2201 Vmath::Vcopy(nq, tmp_forcing, 1, m_bcsForcing[2], 1);
2202 }
2203 }
2204
2205 file = m_sessionName + "_v_5.bc";
2206
2207 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef_v;
2208 std::vector<std::vector<NekDouble>> FieldData_v;
2209 fld->Import(file, FieldDef_v, FieldData_v);
2210 Ilayer->ExtractDataToCoeffs(FieldDef_v[0], FieldData_v[0],
2211 FieldDef_v[0]->m_fields[0],
2212 Ilayer->UpdateCoeffs());
2213 Ilayer->BwdTrans(Ilayer->GetCoeffs(), Ilayer->UpdatePhys());
2214 if (cnt == 0)
2215 {
2216 Vmath::Vcopy(nq, Ilayer->UpdatePhys(), 1, m_bcsForcing[3], 1);
2217 }
2218 Vmath::Vcopy(nq, Ilayer->UpdatePhys(), 1, m_bcsForcing[1], 1);
2219 if (cnt != 0)
2220 {
2221 if (m_vwiRelaxation == 1.0)
2222 {
2223 Vmath::Vcopy(nq, m_bcsForcing[1], 1, tmp_forcing, 1);
2224 }
2225 Vmath::Smul(nq, 1.0 - m_vwiRelaxation, m_bcsForcing[1], 1,
2226 m_bcsForcing[1], 1);
2227 Vmath::Svtvp(nq, m_vwiRelaxation, m_bcsForcing[3], 1, m_bcsForcing[1],
2228 1, Ilayer->UpdatePhys(), 1);
2229 // generate again the bcs files:
2230 Array<OneD, Array<OneD, NekDouble>> fieldcoeffs(1);
2231 Ilayer->FwdTransLocalElmt(Ilayer->GetPhys(), Ilayer->UpdateCoeffs());
2232 fieldcoeffs[0] = Ilayer->UpdateCoeffs();
2233 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef2 =
2234 Ilayer->GetFieldDefinitions();
2235 std::vector<std::vector<NekDouble>> FieldData_2(FieldDef2.size());
2236 ;
2237 FieldDef2[0]->m_fields.push_back("v");
2238 Ilayer->AppendFieldData(FieldDef2[0], FieldData_2[0]);
2239 fld->Write(file, FieldDef2, FieldData_2);
2240 // save the bcs for the next iteration
2241 if (m_vwiRelaxation != 1.0)
2242 {
2243 Vmath::Smul(nq, 1. / (1.0 - m_vwiRelaxation), m_bcsForcing[1], 1,
2244 m_bcsForcing[1], 1);
2245 Vmath::Vcopy(nq, m_bcsForcing[1], 1, m_bcsForcing[3], 1);
2246 }
2247 else
2248 {
2249 Vmath::Vcopy(nq, tmp_forcing, 1, m_bcsForcing[3], 1);
2250 }
2251 }
2252
2253 cnt++;
2254}
2255} // 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:810
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:195
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:132
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:53
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:1702
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:1112
Nektar::VortexWaveInteraction::CopyFile
void CopyFile(std::string file1end, std::string file2end)
Definition: VortexWaveInteraction.cpp:1151
Nektar::VortexWaveInteraction::GetReflectionIndex
Array< OneD, int > GetReflectionIndex(void)
Definition: VortexWaveInteraction.cpp:2042
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:1181
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:1200
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:1132
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:1163
Nektar::VortexWaveInteraction::m_sessionWave
LibUtilities::SessionReaderSharedPtr m_sessionWave
Definition: VortexWaveInteraction.h:267
Nektar::VortexWaveInteraction::UpdateDAlphaDWaveForceMag
void UpdateDAlphaDWaveForceMag(NekDouble alphainit)
Definition: VortexWaveInteraction.cpp:1922
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:113
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:66
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:42
Nektar::IncNavierStokesSharedPtr
std::shared_ptr< IncNavierStokes > IncNavierStokesSharedPtr
Definition: IncNavierStokes.h:300
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:298
tinysimd::log
scalarT< T > log(scalarT< T > in)
Definition: scalar.hpp:303
tinysimd::sqrt
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:294
Generated on Wed Oct 30 2024 16:07:34 for Nektar++ by doxygen 1.9.4