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 <VortexWaveInteraction/VortexWaveInteraction.h>
40
41using namespace std;
42
43namespace Nektar
44{
45
46VortexWaveInteraction::VortexWaveInteraction(int argc, char *argv[])
47 : m_nOuterIterations(0)
48{
49
50 int storesize; // number of previous iterations to store
51
52 m_sessionName = argv[argc - 1];
53 string meshfile = m_sessionName + ".xml";
54
55 LibUtilities::SessionReaderSharedPtr session;
56
57 std::string VWICondFile(argv[argc - 1]);
58 VWICondFile += "_VWI.xml";
59 std::vector<std::string> VWIFilenames;
60 VWIFilenames.push_back(meshfile);
61 VWIFilenames.push_back(VWICondFile);
62
63 // Create Incompressible NavierStokesSolver session reader.
64 m_sessionVWI =
65 LibUtilities::SessionReader::CreateInstance(argc, argv, VWIFilenames);
66 m_sessionVWI->InitSession();
67
68 m_sessionVWI->LoadParameter("AlphaStep", m_alphaStep, 0.05);
69 m_sessionVWI->LoadParameter("OuterIterationStoreSize", storesize, 10);
70
71 // set a low tol for interfaceVWI
72 m_sessionVWI->LoadParameter("EigenvalueRelativeTol", m_eigRelTol, 1e-3);
73
74 m_sessionVWI->LoadParameter("NeutralPointTolerance", m_neutralPointTol,
75 1e-4);
76 m_sessionVWI->LoadParameter("MaxOuterIterations", m_maxOuterIterations,
77 100);
78
79 m_sessionVWI->LoadParameter("StartIteration", m_iterStart, 0);
80 m_sessionVWI->LoadParameter("EndIteration", m_iterEnd, 0);
81
82 m_sessionVWI->LoadParameter("WaveForceMagStep", m_waveForceMagStep, 0.01);
83 m_sessionVWI->LoadParameter("DAlphaDWaveForceMag", m_dAlphaDWaveForceMag,
84 0.0);
85 m_sessionVWI->LoadParameter("MaxWaveForceMagIter", m_maxWaveForceMagIter,
86 1);
87 m_sessionVWI->LoadParameter("RollForceScale", m_rollForceScale, 1.0);
88
89 if (m_sessionVWI->DefinesSolverInfo("DeltaFcnApprox"))
90 {
91 m_deltaFcnApprox = true;
92 m_sessionVWI->LoadParameter("DeltaFcnDecay", m_deltaFcnDecay,
93 1.0 / 500);
94 }
95 else
96 {
97 m_deltaFcnApprox = false;
98 m_deltaFcnDecay = 0.0;
99 }
100
101 if (m_sessionVWI->DefinesSolverInfo("LinfPressureNorm"))
102 {
103 m_useLinfPressureNorm = true;
104 }
105 else
106 {
107 m_useLinfPressureNorm = false;
108 }
109
110 if (m_sessionVWI->DefinesSolverInfo("INTERFACE"))
111 {
112 m_iterinterface = true;
113 }
114 else
115 {
116 m_iterinterface = false;
117 }
118
119 if (m_sessionVWI->DefinesSolverInfo("MoveMeshToCriticalLayer"))
120 {
121 m_moveMeshToCriticalLayer = true;
122 }
123 else
124 {
125 m_moveMeshToCriticalLayer = false;
126 }
127
128 m_alpha = Array<OneD, NekDouble>(storesize);
129 m_alpha[0] = m_sessionVWI->GetParameter("Alpha");
130 m_waveForceMag = Array<OneD, NekDouble>(storesize);
131 m_waveForceMag[0] = m_sessionVWI->GetParameter("WaveForceMag");
132
133 m_leading_real_evl = Array<OneD, NekDouble>(storesize);
134 m_leading_imag_evl = Array<OneD, NekDouble>(storesize);
135
136 if (m_sessionVWI->DefinesParameter("Relaxation"))
137 {
138 m_vwiRelaxation = m_sessionVWI->GetParameter("Relaxation");
139 // fix minimum number of iterations to be number of
140 // iterations required to make contribution of innitial
141 // forcing to 0.1
142 m_minInnerIterations = (int)(log(0.1) / log(m_vwiRelaxation));
143 }
144 else
145 {
146 m_vwiRelaxation = 0.0;
147 m_minInnerIterations = 1;
148 }
149
150 // Initialise NS Roll solver
151 std::string IncCondFile(argv[argc - 1]);
152 IncCondFile += "_IncNSCond.xml";
153 std::vector<std::string> IncNSFilenames;
154 IncNSFilenames.push_back(meshfile);
155 IncNSFilenames.push_back(IncCondFile);
156
157 // Create Incompressible NavierStokesSolver session reader.
158 m_sessionRoll = LibUtilities::SessionReader::CreateInstance(
159 argc, argv, IncNSFilenames, m_sessionVWI->GetComm());
160 m_graphRoll = SpatialDomains::MeshGraph::Read(m_sessionRoll);
161 std::string vEquation = m_sessionRoll->GetSolverInfo("SolverType");
162 m_solverRoll = GetEquationSystemFactory().CreateInstance(
163 vEquation, m_sessionRoll, m_graphRoll);
164 m_solverRoll->PrintSummary(cout);
165
166 if (m_sessionRoll->DefinesSolverInfo("INTERFACE"))
167 {
168 m_sessionVWI->LoadParameter("EigenvalueRelativeTol", m_eigRelTol, 1e-2);
169 }
170
171 int ncoeffs = m_solverRoll->UpdateFields()[0]->GetNcoeffs();
172 m_vwiForcing = Array<OneD, Array<OneD, NekDouble>>(4);
173 m_vwiForcing[0] = Array<OneD, NekDouble>(4 * ncoeffs);
174 for (int i = 1; i < 4; ++i)
175 {
176 m_vwiForcing[i] = m_vwiForcing[i - 1] + ncoeffs;
177 }
178
179 // Fill forcing into m_vwiForcing
180 // Has forcing even been initialised yet?
181 // Vmath::Vcopy(ncoeffs,m_solverRoll->UpdateForces()[0]->GetCoeffs(),1,m_vwiForcing[0],1);
182 // Vmath::Vcopy(ncoeffs,m_solverRoll->UpdateForces()[1]->GetCoeffs(),1,m_vwiForcing[1],1);
183
184 // Create AdvDiff Streak solver
185 std::string AdvDiffCondFile(argv[argc - 1]);
186 AdvDiffCondFile += "_AdvDiffCond.xml";
187 std::vector<std::string> AdvDiffFilenames;
188 AdvDiffFilenames.push_back(meshfile);
189 AdvDiffFilenames.push_back(AdvDiffCondFile);
190
191 // Create AdvDiffusion session reader.
192 m_sessionStreak = LibUtilities::SessionReader::CreateInstance(
193 argc, argv, AdvDiffFilenames, m_sessionVWI->GetComm());
194 m_graphStreak = SpatialDomains::MeshGraph::Read(m_sessionStreak);
195
196 // Initialise LinNS solver
197 std::string LinNSCondFile(argv[argc - 1]);
198 LinNSCondFile += "_LinNSCond.xml";
199 std::vector<std::string> LinNSFilenames;
200 LinNSFilenames.push_back(meshfile);
201 LinNSFilenames.push_back(LinNSCondFile);
202
203 // Create Linearised NS stability session reader.
204 m_sessionWave = LibUtilities::SessionReader::CreateInstance(
205 argc, argv, LinNSFilenames, m_sessionVWI->GetComm());
206 m_graphWave = SpatialDomains::MeshGraph::Read(m_sessionWave);
207
208 // Set the initial beta value in stability to be equal to VWI file
209 std::string LZstr("LZ");
210 NekDouble LZ = 2 * M_PI / m_alpha[0];
211 cout << "Setting LZ in Linearised solver to " << LZ << endl;
212 m_sessionWave->SetParameter(LZstr, LZ);
213
214 // Check for iteration type
215 if (m_sessionVWI->DefinesSolverInfo("VWIIterationType"))
216 {
217 std::string IterationTypeStr =
218 m_sessionVWI->GetSolverInfo("VWIIterationType");
219 for (int i = 0; i < (int)eVWIIterationTypeSize; ++i)
220 {
221 if (boost::iequals(VWIIterationTypeMap[i], IterationTypeStr))
222 {
223 m_VWIIterationType = (VWIIterationType)i;
224 break;
225 }
226 }
227 }
228 else
229 {
230 m_VWIIterationType = eFixedAlphaWaveForcing;
231 }
232
233 // Check for restart
234 bool restart;
235 m_sessionVWI->MatchSolverInfo("RestartIteration", "True", restart, false);
236 if (restart)
237 {
238 switch (m_VWIIterationType)
239 {
240 case eFixedAlpha:
241 break;
242 case eFixedWaveForcing:
243 {
244 FILE *fp;
245 // Check for OuterIter.his file to read
246 if ((fp = fopen("OuterIter.his", "r")))
247 {
248 char buf[BUFSIZ];
249 std::vector<NekDouble> Alpha, Growth, Phase;
250 NekDouble alpha, growth, phase;
251 while (fgets(buf, BUFSIZ, fp))
252 {
253 sscanf(buf, "%*d:%lf%lf%lf", &alpha, &growth, &phase);
254 Alpha.push_back(alpha);
255 Growth.push_back(growth);
256 Phase.push_back(phase);
257 }
258
259 m_nOuterIterations = Alpha.size();
260
261 int nvals =
262 std::min(m_nOuterIterations, (int)m_alpha.size());
263
264 for (int i = 0; i < nvals; ++i)
265 {
266 m_alpha[nvals - 1 - i] =
267 Alpha[m_nOuterIterations - nvals + i];
268 m_leading_real_evl[nvals - 1 - i] =
269 Growth[m_nOuterIterations - nvals + i];
270 m_leading_imag_evl[nvals - 1 - i] =
271 Phase[m_nOuterIterations - nvals + i];
272 }
273
274 UpdateAlpha(m_nOuterIterations++);
275 }
276 else
277 {
278 cout << " No File OuterIter.his to restart from" << endl;
279 }
280 }
281 break;
282 case eFixedAlphaWaveForcing:
283 {
284 string nstr = std::to_string(m_iterStart);
285 cout << "Restarting from iteration " << m_iterStart << endl;
286 std::string rstfile = "cp -f Save/" + m_sessionName + ".rst." +
287 nstr + " " + m_sessionName + ".rst";
288 cout << " " << rstfile << endl;
289 if (system(rstfile.c_str()))
290 {
291 ASSERTL0(false, rstfile.c_str());
292 }
293 std::string vwifile = "cp -f Save/" + m_sessionName + ".vwi." +
294 nstr + " " + m_sessionName + ".vwi";
295 cout << " " << vwifile << endl;
296 if (system(vwifile.c_str()))
297 {
298 ASSERTL0(false, vwifile.c_str());
299 }
300 }
301 break;
302 default:
303 ASSERTL0(false, "Unknown VWIITerationType in restart");
304 }
305 }
306
307 // Check for ConveredSoln to update DAlphaDWaveForce
308 {
309 FILE *fp;
310 if ((fp = fopen("ConvergedSolns", "r")))
311 {
312 char buf[BUFSIZ];
313 std::vector<NekDouble> WaveForce, Alpha;
314 NekDouble waveforce, alpha;
315 while (fgets(buf, BUFSIZ, fp))
316 {
317 sscanf(buf, "%*d:%lf%lf", &waveforce, &alpha);
318 WaveForce.push_back(waveforce);
319 Alpha.push_back(alpha);
320 }
321
322 if (Alpha.size() > 1)
323 {
324 int min_i = 0;
325 NekDouble min_alph = fabs(m_alpha[0] - Alpha[min_i]);
326 // find nearest point
327 for (int i = 1; i < Alpha.size(); ++i)
328 {
329 if (fabs(m_alpha[0] - Alpha[min_i]) < min_alph)
330 {
331 min_i = i;
332 min_alph = fabs(m_alpha[0] - Alpha[min_i]);
333 }
334 }
335
336 // find next nearest point
337 int min_j = (min_i == 0) ? 1 : 0;
338 min_alph = fabs(m_alpha[0] - Alpha[min_j]);
339 for (int i = 0; i < Alpha.size(); ++i)
340 {
341 if (i != min_i)
342 {
343 if (fabs(m_alpha[0] - Alpha[min_j]) < min_alph)
344 {
345 min_j = i;
346 min_alph = fabs(m_alpha[0] - Alpha[min_j]);
347 }
348 }
349 }
350
351 if (fabs(Alpha[min_i] - Alpha[min_j]) > 1e-4)
352 {
353 m_dAlphaDWaveForceMag =
354 (Alpha[min_i] - Alpha[min_j]) /
355 (WaveForce[min_i] - WaveForce[min_j]);
356 }
357 }
358 }
359 }
360}
361
362VortexWaveInteraction::~VortexWaveInteraction()
363{
364 m_sessionVWI->Finalise();
365}
366
367void VortexWaveInteraction::ExecuteRoll(void)
368{
369 // set up the equation system to update the mesh
370 if (m_sessionRoll->DefinesSolverInfo("INTERFACE"))
371 {
372 string vEquation = m_sessionRoll->GetSolverInfo("solvertype");
373 EquationSystemSharedPtr solverRoll =
374 GetEquationSystemFactory().CreateInstance(vEquation, m_sessionRoll,
375 m_graphRoll);
376 // The forcing terms are inserted as N bcs
377 // Execute Roll
378 cout << "Executing Roll solver" << endl;
379 solverRoll->DoInitialise();
380 solverRoll->DoSolve();
381 solverRoll->Output();
382 m_rollField = solverRoll->UpdateFields();
383 for (int g = 0; g < solverRoll->GetNvariables(); ++g)
384 {
385 NekDouble vL2Error = solverRoll->L2Error(g, false);
386 NekDouble vLinfError = solverRoll->LinfError(g);
387 cout << "L 2 error (variable " << solverRoll->GetVariable(g)
388 << ") : " << vL2Error << endl;
389 cout << "L inf error (variable " << solverRoll->GetVariable(g)
390 << ") : " << vLinfError << endl;
391 }
392 }
393 else
394 {
395 if (m_moveMeshToCriticalLayer)
396 {
397 string vEquation = m_sessionRoll->GetSolverInfo("solvertype");
398 EquationSystemSharedPtr solverRoll =
399 GetEquationSystemFactory().CreateInstance(
400 vEquation, m_sessionRoll, m_graphRoll);
401 }
402 else
403 {
404 const int npoints = m_solverRoll->GetNpoints();
405 const int ncoeffs = m_solverRoll->GetNcoeffs();
406
407 static int init = 1;
408 if (init)
409 {
410 m_solverRoll->DoInitialise();
411 m_vwiForcingObj =
412 std::dynamic_pointer_cast<SolverUtils::ForcingProgrammatic>(
413 GetForcingFactory().CreateInstance(
414 "Programmatic", m_sessionRoll, m_solverRoll,
415 m_solverRoll->UpdateFields(),
416 m_solverRoll->UpdateFields().size() - 1, 0));
417
418 std::vector<std::string> vFieldNames =
419 m_sessionRoll->GetVariables();
420 vFieldNames.erase(vFieldNames.end() - 1);
421
422 m_solverRoll->GetFunction("BodyForce")
423 ->Evaluate(vFieldNames, m_vwiForcingObj->UpdateForces());
424
425 // Scale forcing
426 for (int i = 0; i < m_vwiForcingObj->UpdateForces().size(); ++i)
427 {
428 m_solverRoll->UpdateFields()[0]->FwdTrans(
429 m_vwiForcingObj->UpdateForces()[i],
430 m_vwiForcing[2 + i]);
431 Vmath::Smul(npoints, m_rollForceScale,
432 m_vwiForcingObj->UpdateForces()[i], 1,
433 m_vwiForcingObj->UpdateForces()[i], 1);
434 }
435
436 IncNavierStokesSharedPtr ins =
437 m_solverRoll->as<IncNavierStokes>();
438 ins->AddForcing(m_vwiForcingObj);
439
440 init = 0;
441 }
442 else // use internal definition of forcing in m_vwiForcing
443 {
444 // Scale forcing
445 for (int i = 0; i < m_vwiForcingObj->UpdateForces().size(); ++i)
446 {
447 m_solverRoll->UpdateFields()[i]->BwdTrans(
448 m_vwiForcing[i], m_vwiForcingObj->UpdateForces()[i]);
449 Vmath::Smul(npoints, m_rollForceScale,
450 m_vwiForcingObj->UpdateForces()[i], 1,
451 m_vwiForcingObj->UpdateForces()[i], 1);
452 }
453
454 // Shift m_vwiForcing for new restart in case of relaxation
455 Vmath::Vcopy(ncoeffs, m_vwiForcing[0], 1, m_vwiForcing[2], 1);
456 Vmath::Vcopy(ncoeffs, m_vwiForcing[1], 1, m_vwiForcing[3], 1);
457 }
458 }
459
460 // Execute Roll
461 cout << "Executing Roll solver" << endl;
462 m_solverRoll->DoSolve();
463 m_solverRoll->Output();
464 m_rollField = m_solverRoll->UpdateFields();
465 for (int g = 0; g < m_solverRoll->GetNvariables(); ++g)
466 {
467 NekDouble vL2Error = m_solverRoll->L2Error(g, false);
468 NekDouble vLinfError = m_solverRoll->LinfError(g);
469 cout << "L 2 error (variable " << m_solverRoll->GetVariable(g)
470 << ") : " << vL2Error << endl;
471 cout << "L inf error (variable " << m_solverRoll->GetVariable(g)
472 << ") : " << vLinfError << endl;
473 }
474 }
475
476 // Copy .fld file to .rst and base.fld
477 cout << "Executing cp -f session.fld session.rst" << endl;
478 CopyFile(".fld", ".rst");
479
480 // Write out data into base flow with variable Vx,Vy
481 cout << "Writing data to session-Base.fld" << endl;
482
483 std::vector<std::string> variables(2);
484 variables[0] = "Vx";
485 variables[1] = "Vy";
486 std::vector<Array<OneD, NekDouble>> outfield(2);
487 outfield[0] = m_solverRoll->UpdateFields()[0]->UpdateCoeffs();
488 outfield[1] = m_solverRoll->UpdateFields()[1]->UpdateCoeffs();
489 std::string outname = m_sessionName + "-Base.fld";
490 m_solverRoll->WriteFld(outname, m_solverRoll->UpdateFields()[0], outfield,
491 variables);
492}
493
494void VortexWaveInteraction::ExecuteStreak(void)
495{
496 // Create driver
497#if 1
498 std::string vDriverModule;
499 m_sessionStreak->LoadSolverInfo("Driver", vDriverModule, "Standard");
500
501 DriverSharedPtr solverStreak = GetDriverFactory().CreateInstance(
502 vDriverModule, m_sessionStreak, m_graphStreak);
503 solverStreak->Execute();
504
505 m_streakField = solverStreak->GetEqu()[0]->UpdateFields();
506#else
507 // Setup and execute Advection Diffusion solver
508 string vEquation = m_sessionStreak->GetSolverInfo("EqType");
509 EquationSystemSharedPtr solverStreak =
510 GetEquationSystemFactory().CreateInstance(vEquation, m_sessionStreak,
511 m_graphStreak);
512
513 cout << "Executing Streak Solver" << endl;
514 solverStreak->DoInitialise();
515 solverStreak->DoSolve();
516 solverStreak->Output();
517
518 m_streakField = solverStreak->UpdateFields();
519
520 if (m_sessionVWI->DefinesSolverInfo("INTERFACE"))
521 {
522 for (int g = 0; g < solverStreak->GetNvariables(); ++g)
523 {
524 NekDouble vL2Error = solverStreak->L2Error(g, false);
525 NekDouble vLinfError = solverStreak->LinfError(g);
526 cout << "L 2 error (variable " << solverStreak->GetVariable(g)
527 << ") : " << vL2Error << endl;
528 cout << "L inf error (variable " << solverStreak->GetVariable(g)
529 << ") : " << vLinfError << endl;
530 }
531 }
532#endif
533
534 cout << "Executing cp -f session.fld session_streak.fld" << endl;
535 CopyFile(".fld", "_streak.fld");
536}
537
538void VortexWaveInteraction::ExecuteWave(void)
539{
540
541 // Set the initial beta value in stability to be equal to VWI file
542 std::string LZstr("LZ");
543 NekDouble LZ = 2 * M_PI / m_alpha[0];
544 cout << "Setting LZ in Linearised solver to " << LZ << endl;
545 m_sessionWave->SetParameter(LZstr, LZ);
546
547 // Create driver
548 std::string vDriverModule;
549 m_sessionWave->LoadSolverInfo("Driver", vDriverModule, "ModifiedArnoldi");
550 cout << "Setting up linearised NS Solver" << endl;
551 std::shared_ptr<DriverArnoldi> solverWave =
552 std::dynamic_pointer_cast<SolverUtils::DriverArnoldi>(
553 GetDriverFactory().CreateInstance(vDriverModule, m_sessionWave,
554 m_graphWave));
555
556 /// Do linearised NavierStokes Session with Modified Arnoldi
557 cout << "Executing wave solution " << endl;
558 solverWave->Execute();
559
560 // Copy file to a rst location for next restart
561 cout << "Executing cp -f session_eig_0 session_eig_0.rst" << endl;
562 CopyFile("_eig_0", "_eig_0.rst");
563
564 // Store data relevant to other operations
565 m_leading_real_evl[0] = solverWave->GetRealEvl()[0];
566 m_leading_imag_evl[0] = solverWave->GetImagEvl()[0];
567
568 // note this will only be true for modified Arnoldi
569 NekDouble realShift = 0.0;
570 if (m_sessionWave->DefinesParameter("RealShift"))
571 {
572 bool defineshift;
573 // only use shift in modifiedArnoldi solver since
574 // implicitly handled in Arpack.
575 m_sessionWave->MatchSolverInfo("Driver", "ModifiedArnoldi", defineshift,
576 true);
577 if (defineshift)
578 {
579 realShift = m_sessionWave->GetParameter("RealShift");
580 }
581 }
582
583 // Set up leading eigenvalue from inverse
584 NekDouble invmag = 1.0 / (m_leading_real_evl[0] * m_leading_real_evl[0] +
585 m_leading_imag_evl[0] * m_leading_imag_evl[0]);
586 m_leading_real_evl[0] *= -invmag;
587 m_leading_real_evl[0] += realShift;
588 m_leading_imag_evl[0] *= invmag;
589
590 m_waveVelocities = solverWave->GetEqu()[0]->UpdateFields();
591 m_wavePressure = solverWave->GetEqu()[0]->GetPressure();
592
593 if (m_sessionVWI->DefinesSolverInfo("INTERFACE"))
594 {
595 cout << "Growth =" << m_leading_real_evl[0] << endl;
596 cout << "Phase =" << m_leading_imag_evl[0] << endl;
597 }
598}
599
600void VortexWaveInteraction::CalcNonLinearWaveForce(void)
601{
602 if (m_sessionRoll->DefinesSolverInfo("INTERFACE"))
603 {
604 static int cnt = 0;
605 string wavefile = m_sessionName + ".fld";
606 string movedmesh = m_sessionName + "_advPost_moved.xml";
607 string filestreak = m_sessionName + "_streak.fld";
608 string syscall = "";
609 string c = std::to_string(cnt);
610 string c_alpha = std::to_string(m_alpha[0]);
611 if (m_sessionVWI->GetSolverInfo("INTERFACE") == "phase")
612 {
613 string filePost = m_sessionName + "_advPost.xml";
614 syscall = "../../utilities/PostProcessing/Extras/FldCalcBCs " +
615 filePost + " " +
616 "meshhalf_pos_Spen_stability_moved.fld "
617 "meshhalf_pos_Spen_advPost_moved.fld " +
618 c_alpha + " > data_alpha0";
619 cout << syscall.c_str() << endl;
620 if (system(syscall.c_str()))
621 {
622 ASSERTL0(false, syscall.c_str());
623 }
624
625 syscall = "cp -f meshhalf_pos_Spen_stability_moved_u_5.bc " +
626 m_sessionName + "_u_5.bc";
627 cout << syscall.c_str() << endl;
628 if (system(syscall.c_str()))
629 {
630 ASSERTL0(false, syscall.c_str());
631 }
632 syscall = "cp -f meshhalf_pos_Spen_stability_moved_v_5.bc " +
633 m_sessionName + "_v_5.bc";
634 cout << syscall.c_str() << endl;
635 if (system(syscall.c_str()))
636 {
637 ASSERTL0(false, syscall.c_str());
638 }
639 }
640 else
641 {
642 syscall = "../../utilities/PostProcessing/Extras/FldCalcBCs " +
643 movedmesh + " " + wavefile + " " + filestreak + " " +
644 c_alpha + " > datasub_" + c;
645 cout << syscall.c_str() << endl;
646 if (system(syscall.c_str()))
647 {
648 ASSERTL0(false, syscall.c_str());
649 }
650 }
651
652 // relaxation for different alpha values? does it make sense?
653
654 // save the wave
655 string wave_subalp = m_sessionName + "_wave_subalp_" + c + ".fld";
656 syscall = "cp -f " + wavefile + " " + wave_subalp;
657 cout << syscall.c_str() << endl;
658 if (system(syscall.c_str()))
659 {
660 ASSERTL0(false, syscall.c_str());
661 }
662 // FileRelaxation(3);
663 cnt++;
664 }
665 else
666 {
667 int npts = m_waveVelocities[0]->GetPlane(0)->GetNpoints();
668 int ncoeffs = m_waveVelocities[0]->GetPlane(0)->GetNcoeffs();
669 Array<OneD, NekDouble> val(npts), der1(2 * npts);
670 Array<OneD, NekDouble> der2 = der1 + npts;
671 Array<OneD, NekDouble> streak;
672 static int projectfield = -1;
673
674 if (m_deltaFcnApprox)
675 {
676 streak = Array<OneD, NekDouble>(npts);
677 m_streakField[0]->BwdTrans(m_streakField[0]->GetCoeffs(), streak);
678 }
679
680 // Set project field to be first field that has a Neumann
681 // boundary since this not impose any condition on the vertical
682 // boundaries Othersise set to zero.
683 if (projectfield == -1)
684 {
685 Array<OneD, const SpatialDomains::BoundaryConditionShPtr> BndConds;
686
687 for (int i = 0; i < m_waveVelocities.size(); ++i)
688 {
689 BndConds = m_waveVelocities[i]->GetBndConditions();
690 for (int j = 0; j < BndConds.size(); ++j)
691 {
692 if (BndConds[j]->GetBoundaryConditionType() ==
693 SpatialDomains::eNeumann)
694 {
695 projectfield = i;
696 break;
697 }
698 }
699 if (projectfield != -1)
700 {
701 break;
702 }
703 }
704 if (projectfield == -1)
705 {
706 cout << "using first field to project non-linear forcing which "
707 "imposes a Dirichlet condition"
708 << endl;
709 projectfield = 0;
710 }
711 }
712
713 // Shift m_vwiForcing in case of relaxation
714 Vmath::Vcopy(ncoeffs, m_vwiForcing[0], 1, m_vwiForcing[2], 1);
715 Vmath::Vcopy(ncoeffs, m_vwiForcing[1], 1, m_vwiForcing[3], 1);
716
717 // determine inverse of area normalised field.
718 m_wavePressure->GetPlane(0)->BwdTrans(
719 m_wavePressure->GetPlane(0)->GetCoeffs(),
720 m_wavePressure->GetPlane(0)->UpdatePhys());
721 m_wavePressure->GetPlane(1)->BwdTrans(
722 m_wavePressure->GetPlane(1)->GetCoeffs(),
723 m_wavePressure->GetPlane(1)->UpdatePhys());
724 NekDouble invnorm;
725
726 if (m_useLinfPressureNorm)
727 {
728 Vmath::Vmul(npts, m_wavePressure->GetPlane(0)->UpdatePhys(), 1,
729 m_wavePressure->GetPlane(0)->UpdatePhys(), 1, der1, 1);
730 Vmath::Vvtvp(npts, m_wavePressure->GetPlane(1)->UpdatePhys(), 1,
731 m_wavePressure->GetPlane(1)->UpdatePhys(), 1, der1, 1,
732 der1, 1);
733 Vmath::Vsqrt(npts, der1, 1, der1, 1);
734
735 NekDouble Linf = Vmath::Vmax(npts, der1, 1);
736
737 invnorm = 1.0 / Linf;
738 }
739 else
740 {
741 // Determine normalisation of pressure so that |P|/A = 1
742 NekDouble l2;
743 l2 = m_wavePressure->GetPlane(0)->L2(
744 m_wavePressure->GetPlane(0)->GetPhys());
745 invnorm = l2 * l2;
746 l2 = m_wavePressure->GetPlane(1)->L2(
747 m_wavePressure->GetPlane(1)->GetPhys());
748 invnorm += l2 * l2;
749 Vmath::Fill(2 * npts, 1.0, der1, 1);
750 NekDouble area = m_waveVelocities[0]->GetPlane(0)->Integral(der1);
751 cout << "Area: " << area << endl;
752 invnorm = sqrt(area / invnorm);
753 }
754
755 // Get hold of arrays.
756 m_waveVelocities[0]->GetPlane(0)->BwdTrans(
757 m_waveVelocities[0]->GetPlane(0)->GetCoeffs(),
758 m_waveVelocities[0]->GetPlane(0)->UpdatePhys());
759 Array<OneD, NekDouble> u_real =
760 m_waveVelocities[0]->GetPlane(0)->UpdatePhys();
761 Vmath::Smul(npts, invnorm, u_real, 1, u_real, 1);
762 m_waveVelocities[0]->GetPlane(1)->BwdTrans(
763 m_waveVelocities[0]->GetPlane(1)->GetCoeffs(),
764 m_waveVelocities[0]->GetPlane(1)->UpdatePhys());
765 Array<OneD, NekDouble> u_imag =
766 m_waveVelocities[0]->GetPlane(1)->UpdatePhys();
767 Vmath::Smul(npts, invnorm, u_imag, 1, u_imag, 1);
768 m_waveVelocities[1]->GetPlane(0)->BwdTrans(
769 m_waveVelocities[1]->GetPlane(0)->GetCoeffs(),
770 m_waveVelocities[1]->GetPlane(0)->UpdatePhys());
771 Array<OneD, NekDouble> v_real =
772 m_waveVelocities[1]->GetPlane(0)->UpdatePhys();
773 Vmath::Smul(npts, invnorm, v_real, 1, v_real, 1);
774 m_waveVelocities[1]->GetPlane(1)->BwdTrans(
775 m_waveVelocities[1]->GetPlane(1)->GetCoeffs(),
776 m_waveVelocities[1]->GetPlane(1)->UpdatePhys());
777 Array<OneD, NekDouble> v_imag =
778 m_waveVelocities[1]->GetPlane(1)->UpdatePhys();
779 Vmath::Smul(npts, invnorm, v_imag, 1, v_imag, 1);
780
781 // normalise wave for output
782 Vmath::Smul(2 * ncoeffs, invnorm, m_waveVelocities[0]->UpdateCoeffs(),
783 1, m_waveVelocities[0]->UpdateCoeffs(), 1);
784 Vmath::Smul(2 * ncoeffs, invnorm, m_waveVelocities[1]->UpdateCoeffs(),
785 1, m_waveVelocities[1]->UpdateCoeffs(), 1);
786 Vmath::Smul(2 * ncoeffs, invnorm, m_waveVelocities[2]->UpdateCoeffs(),
787 1, m_waveVelocities[2]->UpdateCoeffs(), 1);
788
789 // dump field
790 {
791 std::vector<std::string> variables(3);
792 std::vector<Array<OneD, NekDouble>> outfield(3);
793 variables[0] = "u_w";
794 variables[1] = "v_w";
795 variables[2] = "w_w";
796 outfield[0] = m_waveVelocities[0]->UpdateCoeffs();
797 outfield[1] = m_waveVelocities[1]->UpdateCoeffs();
798 outfield[2] = m_waveVelocities[2]->UpdateCoeffs();
799 std::string outname = m_sessionName + "_wave.fld";
800 m_solverRoll->WriteFld(outname, m_waveVelocities[0], outfield,
801 variables);
802 }
803
804#if 1
805 int ncoeffs_p = m_wavePressure->GetPlane(0)->GetNcoeffs();
806 Vmath::Smul(ncoeffs_p, invnorm,
807 m_wavePressure->GetPlane(0)->UpdateCoeffs(), 1,
808 m_wavePressure->GetPlane(0)->UpdateCoeffs(), 1);
809 Vmath::Smul(ncoeffs_p, invnorm,
810 m_wavePressure->GetPlane(1)->UpdateCoeffs(), 1,
811 m_wavePressure->GetPlane(1)->UpdateCoeffs(), 1);
812#else
813 m_wavePressure->GetPlane(0)->BwdTrans(
814 m_wavePressure->GetPlane(0)->GetCoeffs(),
815 m_wavePressure->GetPlane(0)->UpdatePhys());
816 Vmath::Smul(npts, invnorm, m_wavePressure->GetPlane(0)->UpdatePhys(), 1,
817 m_wavePressure->GetPlane(0)->UpdatePhys(), 1);
818 m_wavePressure->GetPlane(0)->FwdTrans(
819 m_wavePressure->GetPlane(0)->UpdatePhys(),
820 m_wavePressure->GetPlane(0)->UpdateCoeffs());
821
822 m_wavePressure->GetPlane(1)->BwdTrans(
823 m_wavePressure->GetPlane(1)->GetCoeffs(),
824 m_wavePressure->GetPlane(1)->UpdatePhys());
825 Vmath::Smul(npts, invnorm, m_wavePressure->GetPlane(1)->UpdatePhys(), 1,
826 m_wavePressure->GetPlane(1)->UpdatePhys(), 1);
827 m_wavePressure->GetPlane(1)->FwdTrans(
828 m_wavePressure->GetPlane(1)->UpdatePhys(),
829 m_wavePressure->GetPlane(1)->UpdateCoeffs());
830#endif
831
832 // Calculate non-linear terms for x and y directions
833 // d/dx(u u* + u* u)
834 Vmath::Vmul(npts, u_real, 1, u_real, 1, val, 1);
835 Vmath::Vvtvp(npts, u_imag, 1, u_imag, 1, val, 1, val, 1);
836 Vmath::Smul(npts, 2.0, val, 1, val, 1);
837 m_waveVelocities[0]->GetPlane(0)->PhysDeriv(0, val, der1);
838
839 // d/dy(v u* + v* u)
840 Vmath::Vmul(npts, u_real, 1, v_real, 1, val, 1);
841 Vmath::Vvtvp(npts, u_imag, 1, v_imag, 1, val, 1, val, 1);
842 Vmath::Smul(npts, 2.0, val, 1, val, 1);
843 m_waveVelocities[0]->GetPlane(0)->PhysDeriv(1, val, der2);
844
845 Vmath::Vadd(npts, der1, 1, der2, 1, der1, 1);
846#if 1
847 m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(der1,
848 m_vwiForcing[0]);
849#else
850 m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(
851 der1, m_vwiForcing[0]);
852#endif
853 Vmath::Smul(ncoeffs, -m_waveForceMag[0], m_vwiForcing[0], 1,
854 m_vwiForcing[0], 1);
855
856 // d/dx(u v* + u* v)
857 m_waveVelocities[0]->GetPlane(0)->PhysDeriv(0, val, der1);
858
859 // d/dy(v v* + v* v)
860 Vmath::Vmul(npts, v_real, 1, v_real, 1, val, 1);
861 Vmath::Vvtvp(npts, v_imag, 1, v_imag, 1, val, 1, val, 1);
862 Vmath::Smul(npts, 2.0, val, 1, val, 1);
863 m_waveVelocities[0]->GetPlane(0)->PhysDeriv(1, val, der2);
864
865 Vmath::Vadd(npts, der1, 1, der2, 1, der1, 1);
866
867#if 1
868 m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(der1,
869 m_vwiForcing[1]);
870#else
871 m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(
872 der1, m_vwiForcing[1]);
873#endif
874
875 Vmath::Smul(ncoeffs, -m_waveForceMag[0], m_vwiForcing[1], 1,
876 m_vwiForcing[1], 1);
877
878 // by default the symmetrization is on
879 bool symm = true;
880 m_sessionVWI->MatchSolverInfo("Symmetrization", "True", symm, true);
881#if 0
882 if(symm== true )
883 {
884
885 // Symmetrise forcing
886 //-> Get coordinates
887 Array<OneD, NekDouble> coord(2);
888 Array<OneD, NekDouble> coord_x(npts);
889 Array<OneD, NekDouble> coord_y(npts);
890
891 //-> Impose symmetry (x -> -x + Lx/2, y-> -y) on coordinates
892 m_waveVelocities[0]->GetPlane(0)->GetCoords(coord_x,coord_y);
893 NekDouble xmax = Vmath::Vmax(npts,coord_x,1);
894 Vmath::Neg(npts,coord_x,1);
895 Vmath::Sadd(npts,xmax,coord_x,1,coord_x,1);
896 Vmath::Neg(npts,coord_y,1);
897
898 int i, physoffset;
899
900 //-> Obtain list of expansion element ids for each point.
901 Array<OneD, int> Eid(npts);
902 // This search may not be necessary every iteration
903 for(i = 0; i < npts; ++i)
904 {
905 coord[0] = coord_x[i];
906 coord[1] = coord_y[i];
907
908 // Note this will not quite be symmetric.
909 Eid[i] = m_waveVelocities[0]->GetPlane(0)->GetExpIndex(coord,1e-6);
910 }
911
912 // Interpolate field 0
913 m_waveVelocities[0]->GetPlane(0)->BwdTrans(m_vwiForcing[0],der1);
914 for(i = 0; i < npts; ++i)
915 {
916 physoffset = m_waveVelocities[0]->GetPlane(0)->GetPhys_Offset(Eid[i]);
917 coord[0] = coord_x[i];
918 coord[1] = coord_y[i];
919 der2 [i] = m_waveVelocities[0]->GetPlane(0)->GetExp(Eid[i])->PhysEvaluate(coord, der1 + physoffset);
920 }
921 //-> Average field 0
922 Vmath::Vsub(npts,der1,1,der2,1,der2,1);
923 Vmath::Smul(npts,0.5,der2,1,der2,1);
924#if 1
925 m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(der2,m_vwiForcing[0]);
926#else
927 m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(der2, m_vwiForcing[0]);
928#endif
929
930 //-> Interpoloate field 1
931 m_waveVelocities[0]->GetPlane(0)->BwdTrans(m_vwiForcing[1],der1);
932 for(i = 0; i < npts; ++i)
933 {
934 physoffset = m_waveVelocities[0]->GetPlane(0)->GetPhys_Offset(Eid[i]);
935 coord[0] = coord_x[i];
936 coord[1] = coord_y[i];
937 der2[i] = m_waveVelocities[0]->GetPlane(0)->GetExp(Eid[i])->PhysEvaluate(coord, der1 + physoffset);
938 }
939
940 //-> Average field 1
941 Vmath::Vsub(npts,der1,1,der2,1,der2,1);
942 Vmath::Smul(npts,0.5,der2,1,der2,1);
943#if 1
944 m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(der2,m_vwiForcing[1]);
945#else
946 m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(der2, m_vwiForceing[1]);
947#endif
948 }
949#else
950 int i;
951 if (symm == true)
952 {
953 cout << "symmetrization is active" << endl;
954 static Array<OneD, int> index = GetReflectionIndex();
955
956 m_waveVelocities[0]->GetPlane(0)->BwdTrans(m_vwiForcing[0], der1);
957
958 for (i = 0; i < npts; ++i)
959 {
960 if (index[i] != -1)
961 {
962 val[i] = 0.5 * (der1[i] - der1[index[i]]);
963 }
964 }
965#if 1
966 m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(
967 val, m_vwiForcing[0]);
968#else
969 m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(
970 val, m_vwiForcing[0]);
971#endif
972
973 m_waveVelocities[0]->GetPlane(0)->BwdTrans(m_vwiForcing[1], der2);
974 for (i = 0; i < npts; ++i)
975 {
976 if (index[i] != -1)
977 {
978 val[i] = 0.5 * (der2[i] - der2[index[i]]);
979 }
980 }
981#if 1
982 m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(
983 val, m_vwiForcing[1]);
984#else
985 m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(
986 val, m_vwiForcing[1]);
987#endif
988 }
989
990 Vmath::Vmul(npts, der1, 1, der1, 1, val, 1);
991 Vmath::Vvtvp(npts, der2, 1, der2, 1, val, 1, val, 1);
992 Vmath::Vsqrt(npts, val, 1, val, 1);
993 cout << "F_Linf: " << Vmath::Vmax(npts, val, 1) << endl;
994
995#endif
996
997 if (m_vwiRelaxation)
998 {
999 Vmath::Smul(ncoeffs, 1.0 - m_vwiRelaxation, m_vwiForcing[0], 1,
1000 m_vwiForcing[0], 1);
1001 Vmath::Svtvp(ncoeffs, m_vwiRelaxation, m_vwiForcing[2], 1,
1002 m_vwiForcing[0], 1, m_vwiForcing[0], 1);
1003
1004 Vmath::Smul(ncoeffs, 1.0 - m_vwiRelaxation, m_vwiForcing[1], 1,
1005 m_vwiForcing[1], 1);
1006 Vmath::Svtvp(ncoeffs, m_vwiRelaxation, m_vwiForcing[3], 1,
1007 m_vwiForcing[1], 1, m_vwiForcing[1], 1);
1008 }
1009
1010 if (m_deltaFcnApprox)
1011 {
1012 for (int j = 0; j < 2; ++j)
1013 {
1014
1015 m_waveVelocities[projectfield]->GetPlane(0)->BwdTrans(
1016 m_vwiForcing[j], der1);
1017 for (int i = 0; i < npts; ++i)
1018 {
1019 der1[i] *= exp(-streak[i] * streak[i] / m_deltaFcnDecay);
1020 }
1021 m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(
1022 der1, m_vwiForcing[j]);
1023 }
1024 }
1025
1026 // dump output
1027 std::vector<std::string> variables(4);
1028 std::vector<Array<OneD, NekDouble>> outfield(4);
1029 variables[0] = "u";
1030 variables[1] = "v";
1031 variables[2] = "pr";
1032 variables[3] = "pi";
1033 outfield[0] = m_vwiForcing[0];
1034 outfield[1] = m_vwiForcing[1];
1035 Array<OneD, NekDouble> soln(npts, 0.0);
1036 m_wavePressure->GetPlane(0)->BwdTrans(
1037 m_wavePressure->GetPlane(0)->GetCoeffs(),
1038 m_wavePressure->GetPlane(0)->UpdatePhys());
1039 outfield[2] = Array<OneD, NekDouble>(ncoeffs);
1040 m_waveVelocities[0]->GetPlane(0)->FwdTransLocalElmt(
1041 m_wavePressure->GetPlane(0)->GetPhys(), outfield[2]);
1042 m_wavePressure->GetPlane(1)->BwdTrans(
1043 m_wavePressure->GetPlane(1)->GetCoeffs(),
1044 m_wavePressure->GetPlane(1)->UpdatePhys());
1045
1046 Vmath::Vmul(npts, m_wavePressure->GetPlane(0)->UpdatePhys(), 1,
1047 m_wavePressure->GetPlane(0)->UpdatePhys(), 1, val, 1);
1048 Vmath::Vvtvp(npts, m_wavePressure->GetPlane(1)->UpdatePhys(), 1,
1049 m_wavePressure->GetPlane(1)->UpdatePhys(), 1, val, 1, val,
1050 1);
1051 cout << "int P^2: " << m_wavePressure->GetPlane(0)->Integral(val)
1052 << endl;
1053 Vmath::Vsqrt(npts, val, 1, val, 1);
1054 cout << "PLinf: " << Vmath::Vmax(npts, val, 1) << endl;
1055
1056 outfield[3] = Array<OneD, NekDouble>(ncoeffs);
1057 m_waveVelocities[1]->GetPlane(0)->FwdTransLocalElmt(
1058 m_wavePressure->GetPlane(1)->GetPhys(), outfield[3]);
1059
1060 std::string outname = m_sessionName + ".vwi";
1061
1062 m_solverRoll->WriteFld(outname, m_waveVelocities[0]->GetPlane(0),
1063 outfield, variables);
1064 }
1065}
1066
1067void VortexWaveInteraction::CalcL2ToLinfPressure(void)
1068{
1069
1070 ExecuteWave();
1071
1072 m_wavePressure->GetPlane(0)->BwdTrans(
1073 m_wavePressure->GetPlane(0)->GetCoeffs(),
1074 m_wavePressure->GetPlane(0)->UpdatePhys());
1075 m_wavePressure->GetPlane(1)->BwdTrans(
1076 m_wavePressure->GetPlane(1)->GetCoeffs(),
1077 m_wavePressure->GetPlane(1)->UpdatePhys());
1078
1079 int npts = m_waveVelocities[0]->GetPlane(0)->GetNpoints();
1080 NekDouble Linf;
1081 Array<OneD, NekDouble> val(2 * npts, 0.0);
1082
1083 Vmath::Vmul(npts, m_wavePressure->GetPlane(0)->UpdatePhys(), 1,
1084 m_wavePressure->GetPlane(0)->UpdatePhys(), 1, val, 1);
1085 Vmath::Vvtvp(npts, m_wavePressure->GetPlane(1)->UpdatePhys(), 1,
1086 m_wavePressure->GetPlane(1)->UpdatePhys(), 1, val, 1, val, 1);
1087 cout << "int P^2 " << m_wavePressure->GetPlane(0)->Integral(val) << endl;
1088 Vmath::Vsqrt(npts, val, 1, val, 1);
1089
1090 Linf = Vmath::Vmax(npts, val, 1);
1091 cout << "Linf: " << Linf << endl;
1092
1093 NekDouble l2, norm;
1094 l2 =
1095 m_wavePressure->GetPlane(0)->L2(m_wavePressure->GetPlane(0)->GetPhys());
1096 norm = l2 * l2;
1097 l2 =
1098 m_wavePressure->GetPlane(1)->L2(m_wavePressure->GetPlane(1)->GetPhys());
1099 norm += l2 * l2;
1100
1101 Vmath::Fill(npts, 1.0, val, 1);
1102 NekDouble area = m_waveVelocities[0]->GetPlane(0)->Integral(val);
1103
1104 l2 = sqrt(norm / area);
1105
1106 cout << "L2: " << l2 << endl;
1107
1108 cout << "Ratio Linf/L2: " << Linf / l2 << endl;
1109}
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:1057
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: 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::MeshGraph::Read
static MeshGraphSharedPtr Read(const LibUtilities::SessionReaderSharedPtr pSession, LibUtilities::DomainRangeShPtr rng=LibUtilities::NullDomainRangeShPtr, bool fillGraph=true, SpatialDomains::MeshGraphSharedPtr partitionedGraph=nullptr)
Definition: MeshGraph.cpp:115
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:46
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: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:65
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:306
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
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 Jan 31 2024 09:35:58 for Nektar++ by doxygen 1.9.4