Nektar++
ForcingMovingBody.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: ForcingMovingBody.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: Moving Body m_forcing (movement of a body in a domain is
32// achieved via a m_forcing term which is the results of a coordinate system
33// motion)
34//
35///////////////////////////////////////////////////////////////////////////////
36
39
40namespace Nektar
41{
42
45 "MovingBody", ForcingMovingBody::create, "Moving Body Forcing");
46
49 const std::weak_ptr<SolverUtils::EquationSystem> &pEquation)
50 : Forcing(pSession, pEquation)
51{
52}
53
56 [[maybe_unused]] const unsigned int &pNumForcingFields,
57 const TiXmlElement *pForce)
58{
59 // Just 3D homogenous 1D problems can use this techinque
60 ASSERTL0(pFields[0]->GetExpType() == MultiRegions::e3DH1D,
61 "Moving body implemented just for 3D Homogenous 1D expansions.");
62
63 // At this point we know in the xml file where those quantities
64 // are declared (equation or file) - via a function name which is now
65 // stored in funcNameD etc. We need now to fill in with this info the
66 // m_zta and m_eta vectors (actuallythey are matrices) Array to control if
67 // the motion is determined by an equation or is from a file.(not Nektar++)
68 // check if we need to load a file or we have an equation
69 CheckIsFromFile(pForce);
70
71 // Initialise movingbody filter
72 InitialiseFilter(m_session, pFields, pForce);
73
74 // Initialise the cable model
76
77 // Load mapping
79 m_mapping->SetTimeDependent(true);
80
81 if (m_vdim > 0)
82 {
83 m_mapping->SetFromFunction(false);
84 }
85
88 // What are this bi-dimensional vectors
89 // ------------------------------------------ m_zta[0] = m_zta | m_eta[0] =
90 // m_eta | m_zta[1] = d(m_zta)/dt |
91 // m_eta[1] = d(m_eta)/dt | m_zta[2] = dd(m_zta)/ddtt |
92 // m_eta[2] = dd(m_eta)/ddtt |
93 //--------------------------------------------------------------------------------
94 size_t phystot = pFields[0]->GetTotPoints();
95 for (size_t i = 0; i < m_zta.size(); i++)
96 {
97 m_zta[i] = Array<OneD, NekDouble>(phystot, 0.0);
98 m_eta[i] = Array<OneD, NekDouble>(phystot, 0.0);
99 }
100}
101
104 [[maybe_unused]] const Array<OneD, Array<OneD, NekDouble>> &inarray,
105 [[maybe_unused]] Array<OneD, Array<OneD, NekDouble>> &outarray,
106 const NekDouble &time)
107{
108 // Update the forces from the calculation of fluid field, which is
109 // implemented in the movingbody filter
110 Array<OneD, NekDouble> Hydroforces(2 * m_np, 0.0);
111 m_MovBodyfilter->UpdateForce(m_session, pFields, Hydroforces, time);
112
113 // for "free" type (m_vdim = 2), the cable vibrates both in streamwise and
114 // crossflow dimections, for "constrained" type (m_vdim = 1), the cable only
115 // vibrates in crossflow dimection, and for "forced" type (m_vdim = 0), the
116 // calbe vibrates specifically along a given function or file
117 if (m_vdim == 1 || m_vdim == 2)
118 {
119 // For free vibration case, displacements, velocities and acceleartions
120 // are obtained through solving structure dynamic model
121 EvaluateStructDynModel(pFields, Hydroforces, time);
122
123 // Convert result to format required by mapping
124 size_t physTot = pFields[0]->GetTotPoints();
127 for (size_t i = 0; i < 3; i++)
128 {
129 coords[i] = Array<OneD, NekDouble>(physTot, 0.0);
130 coordsVel[i] = Array<OneD, NekDouble>(physTot, 0.0);
131 }
132 // Get original coordinates
133 pFields[0]->GetCoords(coords[0], coords[1], coords[2]);
134
135 // Add displacement to coordinates
136 Vmath::Vadd(physTot, coords[0], 1, m_zta[0], 1, coords[0], 1);
137 Vmath::Vadd(physTot, coords[1], 1, m_eta[0], 1, coords[1], 1);
138 // Copy velocities
139 Vmath::Vcopy(physTot, m_zta[1], 1, coordsVel[0], 1);
140 Vmath::Vcopy(physTot, m_eta[1], 1, coordsVel[1], 1);
141
142 // Update mapping
143 m_mapping->UpdateMapping(time, coords, coordsVel);
144 }
145 else if (m_vdim == 0)
146 {
147 // For forced vibration case, load from specified file or function
148 size_t cnt = 0;
149 for (size_t j = 0; j < m_funcName.size(); j++)
150 {
151 if (m_IsFromFile[cnt] && m_IsFromFile[cnt + 1])
152 {
153 ASSERTL0(false, "Motion loading from file needs specific "
154 "implementation: Work in Progress!");
155 }
156 else
157 {
158 GetFunction(pFields, m_session, m_funcName[j], true)
159 ->Evaluate(m_motion[0], m_zta[j], time);
160 GetFunction(pFields, m_session, m_funcName[j], true)
161 ->Evaluate(m_motion[1], m_eta[j], time);
162 cnt = cnt + 2;
163 }
164 }
165
166 // Update mapping
167 m_mapping->UpdateMapping(time);
168
169 // Convert result from mapping
170 size_t physTot = pFields[0]->GetTotPoints();
173 for (size_t i = 0; i < 3; i++)
174 {
175 coords[i] = Array<OneD, NekDouble>(physTot, 0.0);
176 coordsVel[i] = Array<OneD, NekDouble>(physTot, 0.0);
177 }
178 // Get original coordinates
179 pFields[0]->GetCoords(coords[0], coords[1], coords[2]);
180
181 // Get Coordinates and coord velocity from mapping
182 m_mapping->GetCartesianCoordinates(m_zta[0], m_eta[0], coords[2]);
183 m_mapping->GetCoordVelocity(coordsVel);
184
185 // Calculate displacement
186 Vmath::Vsub(physTot, m_zta[0], 1, coords[0], 1, m_zta[0], 1);
187 Vmath::Vsub(physTot, m_eta[0], 1, coords[1], 1, m_eta[0], 1);
188
189 Vmath::Vcopy(physTot, coordsVel[0], 1, m_zta[1], 1);
190 Vmath::Vcopy(physTot, coordsVel[1], 1, m_eta[1], 1);
191
192 for (size_t var = 0; var < 3; var++)
193 {
194 for (size_t plane = 0; plane < m_np; plane++)
195 {
196 size_t n = pFields[0]->GetPlane(plane)->GetTotPoints();
197 size_t offset = plane * n;
198 size_t Offset = var * m_np + plane;
199
200 m_MotionVars[0][Offset] = m_zta[var][offset];
201 m_MotionVars[1][Offset] = m_eta[var][offset];
202 }
203 }
204 }
205 else
206 {
207 ASSERTL0(false, "Unrecogized vibration type for cable's dynamic model");
208 }
209
210 LibUtilities::CommSharedPtr vcomm = pFields[0]->GetComm();
211 size_t colrank = vcomm->GetColumnComm()->GetRank();
212 // Pass the variables of the cable's motion to the movingbody filter
213 if (colrank == 0)
214 {
215 size_t n = m_MotionVars[0].size();
216 Array<OneD, NekDouble> tmpArray(2 * n), tmp(n);
217 Vmath::Vcopy(n, m_MotionVars[0], 1, tmpArray, 1);
218 Vmath::Vcopy(n, m_MotionVars[1], 1, tmp = tmpArray + n, 1);
219 m_MovBodyfilter->UpdateMotion(m_session, pFields, tmpArray, time);
220 }
221}
222
223/**
224 *
225 */
228 Array<OneD, NekDouble> &Hydroforces, [[maybe_unused]] NekDouble time)
229{
230 LibUtilities::CommSharedPtr vcomm = pFields[0]->GetComm();
231 size_t colrank = vcomm->GetColumnComm()->GetRank();
232 size_t nproc = vcomm->GetColumnComm()->GetSize();
233
234 bool homostrip;
235 m_session->MatchSolverInfo("HomoStrip", "True", homostrip, false);
236
237 // number of structral modes and number of strips
238 size_t npts, nstrips;
239 if (!homostrip) // full resolutions
240 {
241 npts = m_session->GetParameter("HomModesZ");
242 }
243 else
244 {
245 m_session->LoadParameter("HomStructModesZ", npts);
246 m_session->LoadParameter("Strip_Z", nstrips);
247
248 // ASSERTL0(nstrips < = npts,
249 // "the number of struct. modes must be larger than that of the
250 // strips.");
251 }
252
253 // the hydrodynamic forces
255 // forces in x-direction
256 fces[0] = Array<OneD, NekDouble>(npts, 0.0);
257 // forces in y-direction
258 fces[1] = Array<OneD, NekDouble>(npts, 0.0);
259
260 // fill the force vectors
261 if (colrank == 0)
262 {
263 if (!homostrip) // full resolutions
264 {
265 Vmath::Vcopy(m_np, Hydroforces, 1, fces[0], 1);
266 Vmath::Vcopy(m_np, Hydroforces + m_np, 1, fces[1], 1);
267 }
268 else // strip modelling
269 {
270 fces[0][0] = Hydroforces[0];
271 fces[1][0] = Hydroforces[m_np];
272 }
273
274 if (!homostrip) // full resolutions
275 {
277 for (size_t i = 1; i < nproc; ++i)
278 {
279 vcomm->GetColumnComm()->Recv(i, tmp);
280 for (size_t n = 0; n < m_np; n++)
281 {
282 for (size_t j = 0; j < 2; ++j)
283 {
284 fces[j][i * m_np + n] = tmp[j * m_np + n];
285 }
286 }
287 }
288 }
289 else // strip modelling
290 // if the body is submerged partly, the fces are filled partly
291 // by the flow induced forces
292 {
294 for (size_t i = 1; i < nstrips; ++i)
295 {
296 vcomm->GetColumnComm()->Recv(i, tmp);
297
298 for (size_t j = 0; j < 2; ++j)
299 {
300 fces[j][i] = tmp[j];
301 }
302 }
303 }
304 }
305 else
306 {
307 if (!homostrip) // full resolutions
308 {
309 vcomm->GetColumnComm()->Send(0, Hydroforces);
310 }
311 else // strip modelling
312 {
313 for (size_t i = 1; i < nstrips; ++i)
314 {
315 if (colrank == i)
316 {
318 tmp[0] = Hydroforces[0];
319 tmp[1] = Hydroforces[m_np];
320 vcomm->GetColumnComm()->Send(0, tmp);
321 }
322 }
323 }
324 }
325
326 if (colrank == 0)
327 {
328 // Fictitious mass method used to stablize the explicit coupling at
329 // relatively lower mass ratio
330 bool fictmass;
331 m_session->MatchSolverInfo("FictitiousMassMethod", "True", fictmass,
332 false);
333 if (fictmass)
334 {
335 NekDouble fictrho, fictdamp;
336 m_session->LoadParameter("FictMass", fictrho);
337 m_session->LoadParameter("FictDamp", fictdamp);
338
339 static NekDouble Betaq_Coeffs[2][2] = {{1.0, 0.0}, {2.0, -1.0}};
340
341 // only consider second order approximation for fictitious variables
342 size_t intOrder = 2;
343 size_t nint = std::min(m_movingBodyCalls + 1, intOrder);
344 size_t nlevels = m_fV[0].size();
345
346 for (size_t i = 0; i < m_motion.size(); ++i)
347 {
348 RollOver(m_fV[i]);
349 RollOver(m_fA[i]);
350
351 size_t Voffset = npts;
352 size_t Aoffset = 2 * npts;
353
354 Vmath::Vcopy(npts, m_MotionVars[i] + Voffset, 1, m_fV[i][0], 1);
355 Vmath::Vcopy(npts, m_MotionVars[i] + Aoffset, 1, m_fA[i][0], 1);
356
357 // Extrapolate to n+1
358 Vmath::Smul(npts, Betaq_Coeffs[nint - 1][nint - 1],
359 m_fV[i][nint - 1], 1, m_fV[i][nlevels - 1], 1);
360 Vmath::Smul(npts, Betaq_Coeffs[nint - 1][nint - 1],
361 m_fA[i][nint - 1], 1, m_fA[i][nlevels - 1], 1);
362
363 for (size_t n = 0; n < nint - 1; ++n)
364 {
365 Vmath::Svtvp(npts, Betaq_Coeffs[nint - 1][n], m_fV[i][n], 1,
366 m_fV[i][nlevels - 1], 1, m_fV[i][nlevels - 1],
367 1);
368 Vmath::Svtvp(npts, Betaq_Coeffs[nint - 1][n], m_fA[i][n], 1,
369 m_fA[i][nlevels - 1], 1, m_fA[i][nlevels - 1],
370 1);
371 }
372
373 // Add the fictitious forces on the RHS of the equation
374 Vmath::Svtvp(npts, fictdamp, m_fV[i][nlevels - 1], 1, fces[i],
375 1, fces[i], 1);
376 Vmath::Svtvp(npts, fictrho, m_fA[i][nlevels - 1], 1, fces[i], 1,
377 fces[i], 1);
378 }
379 }
380 }
381 // structural solver is implemented on the root process
382 if (colrank == 0)
383 {
384 // Tensioned cable model is evaluated in wave space
385 for (size_t n = 0, cn = 1; n < m_vdim; n++, cn--)
386 {
387 Newmark_betaSolver(pFields, fces[cn], m_MotionVars[cn]);
388 }
389 }
390
391 Array<OneD, NekDouble> Motvars(2 * 2 * m_np);
392 // send physical coeffients to all planes of each processor
393 if (!homostrip) // full resolutions
394 {
395 Array<OneD, NekDouble> tmp(2 * 2 * m_np);
396
397 if (colrank != 0)
398 {
399 vcomm->GetColumnComm()->Recv(0, tmp);
400 Vmath::Vcopy(2 * 2 * m_np, tmp, 1, Motvars, 1);
401 }
402 else
403 {
404 for (size_t i = 1; i < nproc; ++i)
405 {
406 for (size_t j = 0; j < 2; j++) // moving dimensions
407 {
408 for (size_t k = 0; k < 2; k++) // disp. and vel.
409 {
410 for (size_t n = 0; n < m_np; n++)
411 {
412 tmp[j * 2 * m_np + k * m_np + n] =
413 m_MotionVars[j][k * npts + i * m_np + n];
414 }
415 }
416 }
417 vcomm->GetColumnComm()->Send(i, tmp);
418 }
419
420 for (size_t j = 0; j < 2; j++)
421 {
422 for (size_t k = 0; k < 2; k++)
423 {
424 for (size_t n = 0; n < m_np; n++)
425 {
426 tmp[j * 2 * m_np + k * m_np + n] =
427 m_MotionVars[j][k * npts + n];
428 }
429 Vmath::Vcopy(2 * 2 * m_np, tmp, 1, Motvars, 1);
430 }
431 }
432 }
433 }
434 else // strip modelling
435 {
436 Array<OneD, NekDouble> tmp(2 * 2);
437
438 if (colrank != 0)
439 {
440 for (size_t j = 1; j < nproc / nstrips; j++)
441 {
442 if (colrank == j * nstrips)
443 {
444 vcomm->GetColumnComm()->Recv(0, tmp);
445
446 for (size_t plane = 0; plane < m_np; plane++)
447 {
448 for (size_t var = 0; var < 2; var++)
449 {
450 for (size_t k = 0; k < 2; k++)
451 {
452 Motvars[var * 2 * m_np + k * m_np + plane] =
453 tmp[var * 2 + k];
454 }
455 }
456 }
457 }
458 }
459
460 for (size_t i = 1; i < nstrips; i++)
461 {
462 for (size_t j = 0; j < nproc / nstrips; j++)
463 {
464 if (colrank == i + j * nstrips)
465 {
466 vcomm->GetColumnComm()->Recv(0, tmp);
467
468 for (size_t plane = 0; plane < m_np; plane++)
469 {
470 for (size_t var = 0; var < 2; var++)
471 {
472 for (size_t k = 0; k < 2; k++)
473 {
474 Motvars[var * 2 * m_np + k * m_np + plane] =
475 tmp[var * 2 + k];
476 }
477 }
478 }
479 }
480 }
481 }
482 }
483 else
484 {
485 for (size_t j = 0; j < 2; ++j)
486 {
487 for (size_t k = 0; k < 2; ++k)
488 {
489 tmp[j * 2 + k] = m_MotionVars[j][k * npts];
490 }
491 }
492
493 for (size_t j = 1; j < nproc / nstrips; j++)
494 {
495 vcomm->GetColumnComm()->Send(j * nstrips, tmp);
496 }
497
498 for (size_t plane = 0; plane < m_np; plane++)
499 {
500 for (size_t j = 0; j < 2; j++)
501 {
502 for (size_t k = 0; k < 2; ++k)
503 {
504 Motvars[j * 2 * m_np + k * m_np + plane] =
505 m_MotionVars[j][k * npts];
506 }
507 }
508 }
509
510 for (size_t i = 1; i < nstrips; ++i)
511 {
512 for (size_t j = 0; j < 2; ++j)
513 {
514 for (size_t k = 0; k < 2; ++k)
515 {
516 tmp[j * 2 + k] = m_MotionVars[j][i + k * npts];
517 }
518 }
519
520 for (size_t j = 0; j < nproc / nstrips; j++)
521 {
522 vcomm->GetColumnComm()->Send(i + j * nstrips, tmp);
523 }
524 }
525 }
526 }
527
528 // Set the m_forcing term based on the motion of the cable
529 for (size_t var = 0; var < 2; var++)
530 {
531 for (size_t plane = 0; plane < m_np; plane++)
532 {
533 size_t n = pFields[0]->GetPlane(plane)->GetTotPoints();
534
536
537 size_t offset = plane * n;
538 size_t xoffset = var * m_np + plane;
539 size_t yoffset = 2 * m_np + xoffset;
540
541 Vmath::Fill(n, Motvars[xoffset], tmp = m_zta[var] + offset, 1);
542 Vmath::Fill(n, Motvars[yoffset], tmp = m_eta[var] + offset, 1);
543 }
544 }
545}
546
547/**
548 *
549 */
551 [[maybe_unused]] const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
552 Array<OneD, NekDouble> &HydroForces, Array<OneD, NekDouble> &BodyMotions)
553{
554 std::string supptype = m_session->GetSolverInfo("SupportType");
555
556 size_t npts = HydroForces.size();
557
560
561 for (size_t i = 0; i < 4; i++)
562 {
563 fft_i[i] = Array<OneD, NekDouble>(npts, 0.0);
564 fft_o[i] = Array<OneD, NekDouble>(npts, 0.0);
565 }
566
567 Vmath::Vcopy(npts, HydroForces, 1, fft_i[0], 1);
568 for (size_t i = 0; i < 3; i++)
569 {
570 Vmath::Vcopy(npts, BodyMotions + i * npts, 1, fft_i[i + 1], 1);
571 }
572
573 // Implement Fourier transformation of the motion variables
574 if (boost::iequals(supptype, "Free-Free"))
575 {
576 for (size_t j = 0; j < 4; ++j)
577 {
578 m_FFT->FFTFwdTrans(fft_i[j], fft_o[j]);
579 }
580 }
581 else if (boost::iequals(supptype, "Pinned-Pinned"))
582 {
583 // TODO:
584 size_t N = fft_i[0].size();
585
586 for (size_t var = 0; var < 4; var++)
587 {
588 for (size_t i = 0; i < N; i++)
589 {
590 fft_o[var][i] = 0;
591
592 for (size_t k = 0; k < N; k++)
593 {
594 fft_o[var][i] +=
595 fft_i[var][k] * sin(M_PI / (N) * (k + 1 / 2) * (i + 1));
596 }
597 }
598 }
599 }
600 else
601 {
602 ASSERTL0(false, "Unrecognized support type for cable's motion");
603 }
604
605 // solve the ODE in the wave space
606 for (size_t i = 0; i < npts; ++i)
607 {
608 size_t nrows = 3;
609
610 Array<OneD, NekDouble> tmp0, tmp1, tmp2;
611 tmp0 = Array<OneD, NekDouble>(3, 0.0);
612 tmp1 = Array<OneD, NekDouble>(3, 0.0);
613 tmp2 = Array<OneD, NekDouble>(3, 0.0);
614
615 for (size_t var = 0; var < 3; var++)
616 {
617 tmp0[var] = fft_o[var + 1][i];
618 }
619
620 tmp2[0] = fft_o[0][i];
621
622 Blas::Dgemv('N', nrows, nrows, 1.0, &(m_CoeffMat_B[i]->GetPtr())[0],
623 nrows, &tmp0[0], 1, 0.0, &tmp1[0], 1);
624 Blas::Dgemv('N', nrows, nrows, 1.0 / m_structrho,
625 &(m_CoeffMat_A[i]->GetPtr())[0], nrows, &tmp2[0], 1, 1.0,
626 &tmp1[0], 1);
627
628 for (size_t var = 0; var < 3; var++)
629 {
630 fft_i[var][i] = tmp1[var];
631 }
632 }
633
634 // get physical coeffients via Backward fourier transformation of wave
635 // coefficients
636 if (boost::iequals(supptype, "Free-Free"))
637 {
638 for (size_t var = 0; var < 3; var++)
639 {
640 m_FFT->FFTBwdTrans(fft_i[var], fft_o[var]);
641 }
642 }
643 else if (boost::iequals(supptype, "Pinned-Pinned"))
644 {
645 // TODO:
646 size_t N = fft_i[0].size();
647
648 for (size_t var = 0; var < 3; var++)
649 {
650 for (size_t i = 0; i < N; i++)
651 {
652 fft_o[var][i] = 0;
653
654 for (size_t k = 0; k < N; k++)
655 {
656 fft_o[var][i] += fft_i[var][k] *
657 sin(M_PI / (N) * (k + 1) * (i + 1 / 2)) *
658 2 / N;
659 }
660 }
661 }
662 }
663 else
664 {
665 ASSERTL0(false, "Unrecognized support type for cable's motion");
666 }
667
668 for (size_t i = 0; i < 3; i++)
669 {
670 Array<OneD, NekDouble> tmp(npts, 0.0);
671 Vmath::Vcopy(npts, fft_o[i], 1, tmp = BodyMotions + i * npts, 1);
672 }
673}
674
675/**
676 *
677 */
679 [[maybe_unused]] const LibUtilities::SessionReaderSharedPtr &pSession,
680 [[maybe_unused]] const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields)
681{
683 m_session->LoadParameter("Kinvis", m_kinvis);
684 m_session->LoadParameter("TimeStep", m_timestep, 0.01);
685
686 LibUtilities::CommSharedPtr vcomm = pFields[0]->GetComm();
687 size_t colrank = vcomm->GetColumnComm()->GetRank();
688 size_t nproc = vcomm->GetColumnComm()->GetSize();
689
690 // number of structral modes
691 size_t npts;
692 bool homostrip;
693 m_session->MatchSolverInfo("HomoStrip", "True", homostrip, false);
694
695 if (!homostrip) // full resolutions
696 {
697 npts = m_session->GetParameter("HomModesZ");
698 }
699 else
700 {
701 m_session->LoadParameter("HomStructModesZ", npts);
702 }
703
705 m_MotionVars[0] = Array<OneD, NekDouble>(3 * npts, 0.0);
706 m_MotionVars[1] = Array<OneD, NekDouble>(3 * npts, 0.0);
707
709 ZIDs = pFields[0]->GetZIDs();
710 m_np = ZIDs.size();
711
712 std::string vibtype = m_session->GetSolverInfo("VibrationType");
713
714 if (boost::iequals(vibtype, "Constrained"))
715 {
716 m_vdim = 1;
717 }
718 else if (boost::iequals(vibtype, "Free"))
719 {
720 m_vdim = 2;
721 }
722 else if (boost::iequals(vibtype, "Forced"))
723 {
724 m_vdim = 0;
725 return;
726 }
727
728 if (!homostrip)
729 {
730 m_session->LoadParameter("LZ", m_lhom);
731 int nplanes = m_session->GetParameter("HomModesZ");
733 nplanes);
734 }
735 else
736 {
737 int nstrips;
738 NekDouble DistStrip;
739
740 m_session->LoadParameter("DistStrip", DistStrip);
741 m_session->LoadParameter("Strip_Z", nstrips);
742 m_lhom = nstrips * DistStrip;
744 nstrips);
745 }
746
747 // load the structural dynamic parameters from xml file
748 m_session->LoadParameter("StructRho", m_structrho);
749 m_session->LoadParameter("StructDamp", m_structdamp, 0.0);
750
751 // Identify whether the fictitious mass method is active for explicit
752 // coupling of fluid solver and structural dynamics solver
753 bool fictmass;
754 m_session->MatchSolverInfo("FictitiousMassMethod", "True", fictmass, false);
755 if (fictmass)
756 {
757 NekDouble fictrho, fictdamp;
758 m_session->LoadParameter("FictMass", fictrho);
759 m_session->LoadParameter("FictDamp", fictdamp);
760 m_structrho += fictrho;
761 m_structdamp += fictdamp;
762
763 // Storage array of Struct Velocity and Acceleration used for
764 // extrapolation of fictitious force
767 for (size_t i = 0; i < m_motion.size(); ++i)
768 {
771
772 for (size_t n = 0; n < 2; ++n)
773 {
774 m_fV[i][n] = Array<OneD, NekDouble>(npts, 0.0);
775 m_fA[i][n] = Array<OneD, NekDouble>(npts, 0.0);
776 }
777 }
778 }
779
780 // Setting the coefficient matrices for solving structural dynamic ODEs
781 SetDynEqCoeffMatrix(pFields);
782
783 // Set initial condition for cable's motion
784 size_t cnt = 0;
785
786 for (size_t j = 0; j < m_funcName.size(); j++)
787 {
788 // loading from the specified files through inputstream
789 if (m_IsFromFile[cnt] && m_IsFromFile[cnt + 1])
790 {
791 std::ifstream inputStream;
792 size_t nzpoints;
793
794 if (homostrip)
795 {
796 m_session->LoadParameter("HomStructModesZ", nzpoints);
797 }
798 else
799 {
800 nzpoints = pFields[0]->GetHomogeneousBasis()->GetNumModes();
801 }
802
803 if (vcomm->GetRank() == 0)
804 {
805 std::string filename =
806 m_session->GetFunctionFilename(m_funcName[j], m_motion[0]);
807
808 // Open intputstream for cable motions
809 inputStream.open(filename.c_str());
810
811 // Import the head string from the file
813 for (size_t n = 0; n < tmp.size(); n++)
814 {
815 inputStream >> tmp[n];
816 }
817
818 NekDouble time, z_cds;
819 // Import the motion variables from the file
820 for (size_t n = 0; n < nzpoints; n++)
821 {
822 inputStream >> std::setprecision(6) >> time;
823 inputStream >> std::setprecision(6) >> z_cds;
824 inputStream >> std::setprecision(8) >> m_MotionVars[0][n];
825 inputStream >> std::setprecision(8) >>
826 m_MotionVars[0][n + nzpoints];
827 inputStream >> std::setprecision(8) >>
828 m_MotionVars[0][n + 2 * nzpoints];
829 inputStream >> std::setprecision(8) >> m_MotionVars[1][n];
830 inputStream >> std::setprecision(8) >>
831 m_MotionVars[1][n + nzpoints];
832 inputStream >> std::setprecision(8) >>
833 m_MotionVars[1][n + 2 * nzpoints];
834 }
835 // Close inputstream for cable motions
836 inputStream.close();
837 }
838 cnt = cnt + 2;
839 }
840 else // Evaluate from the functions specified in xml file
841 {
842 if (!homostrip)
843 {
844 size_t nzpoints =
845 pFields[0]->GetHomogeneousBasis()->GetNumModes();
846 Array<OneD, NekDouble> z_coords(nzpoints, 0.0);
848 pFields[0]->GetHomogeneousBasis()->GetZ();
849
850 Vmath::Smul(nzpoints, m_lhom / 2.0, pts, 1, z_coords, 1);
851 Vmath::Sadd(nzpoints, m_lhom / 2.0, z_coords, 1, z_coords, 1);
852
856 for (size_t plane = 0; plane < m_np; plane++)
857 {
858 x2[plane] = z_coords[ZIDs[plane]];
859 }
860
862 Array<OneD, NekDouble> tmp0(m_np, 0.0);
863 Array<OneD, NekDouble> tmp1(m_np, 0.0);
864 LibUtilities::EquationSharedPtr ffunc0, ffunc1;
865
866 ffunc0 = m_session->GetFunction(m_funcName[j], m_motion[0]);
867 ffunc1 = m_session->GetFunction(m_funcName[j], m_motion[1]);
868 ffunc0->Evaluate(x0, x1, x2, 0.0, tmp0);
869 ffunc1->Evaluate(x0, x1, x2, 0.0, tmp1);
870
871 size_t offset = j * npts;
872 Vmath::Vcopy(m_np, tmp0, 1, tmp = m_MotionVars[0] + offset, 1);
873 Vmath::Vcopy(m_np, tmp1, 1, tmp = m_MotionVars[1] + offset, 1);
874
875 if (colrank == 0)
876 {
877 for (size_t i = 1; i < nproc; ++i)
878 {
879 vcomm->GetColumnComm()->Recv(i, tmp0);
880 vcomm->GetColumnComm()->Recv(i, tmp1);
881 Vmath::Vcopy(m_np, tmp0, 1,
882 tmp = m_MotionVars[0] + offset + i * m_np,
883 1);
884 Vmath::Vcopy(m_np, tmp1, 1,
885 tmp = m_MotionVars[1] + offset + i * m_np,
886 1);
887 }
888 }
889 else
890 {
891 vcomm->GetColumnComm()->Send(0, tmp0);
892 vcomm->GetColumnComm()->Send(0, tmp1);
893 }
894 }
895 else
896 {
897 if (colrank == 0)
898 {
899 size_t nstrips;
900 m_session->LoadParameter("Strip_Z", nstrips);
901
902 ASSERTL0(
903 m_session->DefinesSolverInfo("USEFFT"),
904 "Fourier transformation of cable motion is currently "
905 "implemented only for FFTW module.");
906
907 NekDouble DistStrip;
908 m_session->LoadParameter("DistStrip", DistStrip);
909
910 Array<OneD, NekDouble> x0(npts, 0.0);
911 Array<OneD, NekDouble> x1(npts, 0.0);
912 Array<OneD, NekDouble> x2(npts, 0.0);
913 Array<OneD, NekDouble> tmp(npts, 0.0);
914 Array<OneD, NekDouble> tmp0(npts, 0.0);
915 Array<OneD, NekDouble> tmp1(npts, 0.0);
916 for (size_t plane = 0; plane < npts; plane++)
917 {
918 x2[plane] = plane * DistStrip;
919 }
920 LibUtilities::EquationSharedPtr ffunc0, ffunc1;
921 ffunc0 = m_session->GetFunction(m_funcName[j], m_motion[0]);
922 ffunc1 = m_session->GetFunction(m_funcName[j], m_motion[1]);
923 ffunc0->Evaluate(x0, x1, x2, 0.0, tmp0);
924 ffunc1->Evaluate(x0, x1, x2, 0.0, tmp1);
925
926 size_t offset = j * npts;
927 Vmath::Vcopy(npts, tmp0, 1, tmp = m_MotionVars[0] + offset,
928 1);
929 Vmath::Vcopy(npts, tmp1, 1, tmp = m_MotionVars[1] + offset,
930 1);
931 }
932 }
933
934 cnt = cnt + 2;
935 }
936 }
937}
938
939/**
940 *
941 */
943 [[maybe_unused]] const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields)
944{
945 size_t nplanes;
946
947 bool homostrip;
948 m_session->MatchSolverInfo("HomoStrip", "True", homostrip, false);
949
950 if (!homostrip)
951 {
952 nplanes = m_session->GetParameter("HomModesZ");
953 }
954 else
955 {
956 m_session->LoadParameter("Strip_Z", nplanes);
957 }
958
961
962 NekDouble tmp1, tmp2, tmp3;
963 NekDouble tmp4, tmp5, tmp6, tmp7;
964
965 // load the structural dynamic parameters from xml file
966 NekDouble cabletension;
967 NekDouble bendingstiff;
968 NekDouble structstiff;
969 m_session->LoadParameter("StructStiff", structstiff, 0.0);
970 m_session->LoadParameter("CableTension", cabletension, 0.0);
971 m_session->LoadParameter("BendingStiff", bendingstiff, 0.0);
972
973 tmp1 = m_timestep * m_timestep;
974 tmp2 = structstiff / m_structrho;
975 tmp3 = m_structdamp / m_structrho;
976 tmp4 = cabletension / m_structrho;
977 tmp5 = bendingstiff / m_structrho;
978
979 // solve the ODE in the wave space for cable motion to obtain disp, vel and
980 // accel
981
982 std::string supptype = m_session->GetSolverInfo("SupportType");
983
984 for (size_t plane = 0; plane < nplanes; plane++)
985 {
986 int nel = 3;
987 m_CoeffMat_A[plane] =
989 m_CoeffMat_B[plane] =
991
992 // Initialised to avoid compiler warnings.
993 unsigned int K = 0;
994 NekDouble beta = 0.0;
995
996 if (boost::iequals(supptype, "Free-Free"))
997 {
998 K = plane / 2;
999 beta = 2.0 * M_PI / m_lhom;
1000 }
1001 else if (boost::iequals(supptype, "Pinned-Pinned"))
1002 {
1003 K = plane + 1;
1004 beta = M_PI / m_lhom;
1005 }
1006 else
1007 {
1008 ASSERTL0(false, "Unrecognized support type for cable's motion");
1009 }
1010
1011 tmp6 = beta * K;
1012 tmp6 = tmp6 * tmp6;
1013 tmp7 = tmp6 * tmp6;
1014 tmp7 = tmp2 + tmp4 * tmp6 + tmp5 * tmp7;
1015
1016 (*m_CoeffMat_A[plane])(0, 0) = tmp7;
1017 (*m_CoeffMat_A[plane])(0, 1) = tmp3;
1018 (*m_CoeffMat_A[plane])(0, 2) = 1.0;
1019 (*m_CoeffMat_A[plane])(1, 0) = 0.0;
1020 (*m_CoeffMat_A[plane])(1, 1) = 1.0;
1021 (*m_CoeffMat_A[plane])(1, 2) = -m_timestep / 2.0;
1022 (*m_CoeffMat_A[plane])(2, 0) = 1.0;
1023 (*m_CoeffMat_A[plane])(2, 1) = 0.0;
1024 (*m_CoeffMat_A[plane])(2, 2) = -tmp1 / 4.0;
1025 (*m_CoeffMat_B[plane])(0, 0) = 0.0;
1026 (*m_CoeffMat_B[plane])(0, 1) = 0.0;
1027 (*m_CoeffMat_B[plane])(0, 2) = 0.0;
1028 (*m_CoeffMat_B[plane])(1, 0) = 0.0;
1029 (*m_CoeffMat_B[plane])(1, 1) = 1.0;
1030 (*m_CoeffMat_B[plane])(1, 2) = m_timestep / 2.0;
1031 (*m_CoeffMat_B[plane])(2, 0) = 1.0;
1032 (*m_CoeffMat_B[plane])(2, 1) = m_timestep;
1033 (*m_CoeffMat_B[plane])(2, 2) = tmp1 / 4.0;
1034
1035 m_CoeffMat_A[plane]->Invert();
1036 (*m_CoeffMat_B[plane]) =
1037 (*m_CoeffMat_A[plane]) * (*m_CoeffMat_B[plane]);
1038 }
1039}
1040
1041/**
1042 * Function to roll time-level storages to the next step layout.
1043 * The stored data associated with the oldest time-level
1044 * (not required anymore) are moved to the top, where they will
1045 * be overwritten as the solution process progresses.
1046 */
1048{
1049 size_t nlevels = input.size();
1050
1052 tmp = input[nlevels - 1];
1053
1054 for (size_t n = nlevels - 1; n > 0; --n)
1055 {
1056 input[n] = input[n - 1];
1057 }
1058
1059 input[0] = tmp;
1060}
1061
1062/**
1063 *
1064 */
1065void ForcingMovingBody::CheckIsFromFile(const TiXmlElement *pForce)
1066{
1067
1070 m_motion[0] = "x";
1071 m_motion[1] = "y";
1072
1074 // Loading the x-dispalcement (m_zta) and the y-displacement (m_eta)
1075 // Those two variables are bith functions of z and t and the may come
1076 // from an equation (forced vibration) or from another solver which, given
1077 // the aerodynamic forces at the previous step, calculates the
1078 // displacements.
1079
1080 // Get the body displacement: m_zta and m_eta
1081 const TiXmlElement *funcNameElmt_D =
1082 pForce->FirstChildElement("DISPLACEMENTS");
1083 ASSERTL0(funcNameElmt_D,
1084 "MOVINGBODYFORCE tag has to specify a function name which "
1085 "prescribes the body displacement as d(z,t).");
1086
1087 m_funcName[0] = funcNameElmt_D->GetText();
1088 ASSERTL0(m_session->DefinesFunction(m_funcName[0]),
1089 "Function '" + m_funcName[0] + "' not defined.");
1090
1091 // Get the body velocity of movement: d(m_zta)/dt and d(m_eta)/dt
1092 const TiXmlElement *funcNameElmt_V =
1093 pForce->FirstChildElement("VELOCITIES");
1094 ASSERTL0(funcNameElmt_D,
1095 "MOVINGBODYFORCE tag has to specify a function name which "
1096 "prescribes the body velocity of movement as v(z,t).");
1097
1098 m_funcName[1] = funcNameElmt_V->GetText();
1099 ASSERTL0(m_session->DefinesFunction(m_funcName[1]),
1100 "Function '" + m_funcName[1] + "' not defined.");
1101
1102 // Get the body acceleration: dd(m_zta)/ddt and dd(m_eta)/ddt
1103 const TiXmlElement *funcNameElmt_A =
1104 pForce->FirstChildElement("ACCELERATIONS");
1105 ASSERTL0(funcNameElmt_A,
1106 "MOVINGBODYFORCE tag has to specify a function name which "
1107 "prescribes the body acceleration as a(z,t).");
1108
1109 m_funcName[2] = funcNameElmt_A->GetText();
1110 ASSERTL0(m_session->DefinesFunction(m_funcName[2]),
1111 "Function '" + m_funcName[2] + "' not defined.");
1112
1114
1115 // Check Displacement x
1116 vType = m_session->GetFunctionType(m_funcName[0], m_motion[0]);
1118 {
1119 m_IsFromFile[0] = false;
1120 }
1121 else if (vType == LibUtilities::eFunctionTypeFile)
1122 {
1123 m_IsFromFile[0] = true;
1124 }
1125 else
1126 {
1127 ASSERTL0(false, "The displacements in x must be specified via an "
1128 "equation <E> or a file <F>");
1129 }
1130
1131 // Check Displacement y
1132 vType = m_session->GetFunctionType(m_funcName[0], m_motion[1]);
1134 {
1135 m_IsFromFile[1] = false;
1136 }
1137 else if (vType == LibUtilities::eFunctionTypeFile)
1138 {
1139 m_IsFromFile[1] = true;
1140 }
1141 else
1142 {
1143 ASSERTL0(false, "The displacements in y must be specified via an "
1144 "equation <E> or a file <F>");
1145 }
1146
1147 // Check Velocity x
1148 vType = m_session->GetFunctionType(m_funcName[1], m_motion[0]);
1150 {
1151 m_IsFromFile[2] = false;
1152 }
1153 else if (vType == LibUtilities::eFunctionTypeFile)
1154 {
1155 m_IsFromFile[2] = true;
1156 }
1157 else
1158 {
1159 ASSERTL0(false, "The velocities in x must be specified via an equation "
1160 "<E> or a file <F>");
1161 }
1162
1163 // Check Velocity y
1164 vType = m_session->GetFunctionType(m_funcName[1], m_motion[1]);
1166 {
1167 m_IsFromFile[3] = false;
1168 }
1169 else if (vType == LibUtilities::eFunctionTypeFile)
1170 {
1171 m_IsFromFile[3] = true;
1172 }
1173 else
1174 {
1175 ASSERTL0(false, "The velocities in y must be specified via an equation "
1176 "<E> or a file <F>");
1177 }
1178
1179 // Check Acceleration x
1180 vType = m_session->GetFunctionType(m_funcName[2], m_motion[0]);
1182 {
1183 m_IsFromFile[4] = false;
1184 }
1185 else if (vType == LibUtilities::eFunctionTypeFile)
1186 {
1187 m_IsFromFile[4] = true;
1188 }
1189 else
1190 {
1191 ASSERTL0(false, "The accelerations in x must be specified via an "
1192 "equation <E> or a file <F>");
1193 }
1194
1195 // Check Acceleration y
1196 vType = m_session->GetFunctionType(m_funcName[2], m_motion[1]);
1198 {
1199 m_IsFromFile[5] = false;
1200 }
1201 else if (vType == LibUtilities::eFunctionTypeFile)
1202 {
1203 m_IsFromFile[5] = true;
1204 }
1205 else
1206 {
1207 ASSERTL0(false, "The accelerations in y must be specified via an "
1208 "equation <E> or a file <F>");
1209 }
1210}
1211
1212/**
1213 *
1214 */
1218 const TiXmlElement *pForce)
1219{
1220 // Get the outputfile name, output frequency and
1221 // the boundary's ID for the cable's wall
1222 std::string typeStr = pForce->Attribute("TYPE");
1223 std::map<std::string, std::string> vParams;
1224
1225 const TiXmlElement *param = pForce->FirstChildElement("PARAM");
1226 while (param)
1227 {
1228 ASSERTL0(param->Attribute("NAME"),
1229 "Missing attribute 'NAME' for parameter in filter " + typeStr +
1230 "'.");
1231 std::string nameStr = param->Attribute("NAME");
1232
1233 ASSERTL0(param->GetText(), "Empty value string for param.");
1234 std::string valueStr = param->GetText();
1235
1236 vParams[nameStr] = valueStr;
1237
1238 param = param->NextSiblingElement("PARAM");
1239 }
1240
1241 // Creat the filter for MovingBody, where we performed the calculation of
1242 // fluid forces and write both the aerodynamic forces and motion variables
1243 // into the output files
1245 pSession, m_equ.lock(), vParams);
1246
1247 // Initialise the object of MovingBody filter
1248 m_MovBodyfilter->Initialise(pFields, 0.0);
1249}
1250
1251} // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
GlobalMapping::MappingSharedPtr m_mapping
Array< OneD, std::string > m_funcName
[0] is displacements, [1] is velocities, [2] is accelerations
ForcingMovingBody(const LibUtilities::SessionReaderSharedPtr &pSession, const std::weak_ptr< SolverUtils::EquationSystem > &pEquation)
Array< OneD, Array< OneD, NekDouble > > m_zta
Store the derivatives of motion variables in x-direction.
NekDouble m_lhom
length ratio of the cable
NekDouble m_kinvis
fluid viscous
size_t m_np
number of planes per processors
Array< OneD, DNekMatSharedPtr > m_CoeffMat_A
matrices in Newmart-beta method
Array< OneD, Array< OneD, NekDouble > > m_eta
Store the derivatives of motion variables in y-direction.
void CheckIsFromFile(const TiXmlElement *pForce)
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_fV
fictitious velocity storage
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_fA
fictitious acceleration storage
void SetDynEqCoeffMatrix(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
void RollOver(Array< OneD, Array< OneD, NekDouble > > &input)
void InitialiseCableModel(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
NekDouble m_structrho
mass of the cable per unit length
size_t m_vdim
vibration dimension
NekDouble m_timestep
time step
void v_Apply(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble &time) override
Array< OneD, bool > m_IsFromFile
do determine if the the body motion come from an extern file
FilterMovingBodySharedPtr m_MovBodyfilter
Array< OneD, Array< OneD, NekDouble > > m_MotionVars
storage for the cable's motion(x,y) variables
static SolverUtils::ForcingSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const std::weak_ptr< SolverUtils::EquationSystem > &pEquation, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields, const TiXmlElement *pForce)
Creates an instance of this class.
Array< OneD, DNekMatSharedPtr > m_CoeffMat_B
matrices in Newmart-beta method
void EvaluateStructDynModel(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, Array< OneD, NekDouble > &Hydroforces, NekDouble time)
size_t m_movingBodyCalls
number of times the movbody have been called
void Newmark_betaSolver(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, Array< OneD, NekDouble > &FcePhysinArray, Array< OneD, NekDouble > &MotPhysinArray)
void InitialiseFilter(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const TiXmlElement *pForce)
Array< OneD, std::string > m_motion
motion direction: [0] is 'x' and [1] is 'y'
LibUtilities::NektarFFTSharedPtr m_FFT
static std::string className
Name of the class.
void v_InitObject(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const unsigned int &pNumForcingFields, const TiXmlElement *pForce) override
NekDouble m_structdamp
damping ratio of the cable
static GLOBAL_MAPPING_EXPORT MappingSharedPtr Load(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
Return a pointer to the mapping, creating it on first call.
Definition: Mapping.cpp:264
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Defines a forcing term to be explicitly applied.
Definition: Forcing.h:71
const std::weak_ptr< EquationSystem > m_equ
Weak pointer to equation system using this forcing.
Definition: Forcing.h:125
SOLVER_UTILS_EXPORT SessionFunctionSharedPtr GetFunction(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const LibUtilities::SessionReaderSharedPtr &pSession, std::string pName, bool pCache=false)
Get a SessionFunction by name.
Definition: Forcing.cpp:150
LibUtilities::SessionReaderSharedPtr m_session
Session reader.
Definition: Forcing.h:123
static void Dgemv(const char &trans, const int &m, const int &n, const double &alpha, const double *a, const int &lda, const double *x, const int &incx, const double &beta, double *y, const int &incy)
BLAS level 2: Matrix vector multiply y = alpha A x plus beta y where A[m x n].
Definition: Blas.hpp:211
NektarFFTFactory & GetNektarFFTFactory()
Definition: NektarFFT.cpp:65
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< Equation > EquationSharedPtr
Definition: Equation.h:125
@ beta
Gauss Radau pinned at x=-1,.
Definition: PointsType.h:59
std::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:55
ForcingFactory & GetForcingFactory()
Declaration of the forcing factory singleton.
Definition: Forcing.cpp:41
double NekDouble
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
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
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
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.hpp:54
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
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825
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