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