Nektar++
LinearisedAdvection.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: LinearisedAdvection.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: Evaluation of the linearised advective term
32//
33///////////////////////////////////////////////////////////////////////////////
34
37#include <boost/format.hpp>
38
39namespace Nektar
40{
41
44 "Linearised", LinearisedAdvection::create,
45 "Linearised Non-Conservative");
46
47/**
48 * Constructor. Creates ...
49 *
50 * \param
51 * \param
52 */
53
55{
56}
57
61{
62 Advection::v_InitObject(pSession, pFields);
63
64 m_session = pSession;
67 m_session, pFields[0]->GetGraph());
68 m_spacedim = pFields[0]->GetGraph()->GetSpaceDimension();
69 m_expdim = pFields[0]->GetGraph()->GetMeshDimension();
70
71 // Setting parameters for homogeneous problems
72 m_HomoDirec = 0;
73 m_useFFT = false;
75 m_singleMode = false;
76 m_halfMode = false;
77 m_multipleModes = false;
78
79 if (m_session->DefinesSolverInfo("HOMOGENEOUS"))
80 {
81 std::string HomoStr = m_session->GetSolverInfo("HOMOGENEOUS");
82 m_spacedim = 3;
83
84 if ((HomoStr == "HOMOGENEOUS1D") || (HomoStr == "Homogeneous1D") ||
85 (HomoStr == "1D") || (HomoStr == "Homo1D"))
86 {
88 m_LhomZ = m_session->GetParameter("LZ");
89 m_HomoDirec = 1;
90
91 ASSERTL0(m_session->DefinesSolverInfo("ModeType"),
92 "Need to specify ModeType as HalfMode,SingleMode or "
93 "MultipleModes");
94
95 m_session->MatchSolverInfo("ModeType", "SingleMode", m_singleMode,
96 false);
97 m_session->MatchSolverInfo("ModeType", "HalfMode", m_halfMode,
98 false);
99 m_session->MatchSolverInfo("ModeType", "MultipleModes",
100 m_multipleModes, false);
101
102 if (m_singleMode)
103 {
104 m_npointsZ = 2;
105 }
106 else if (m_halfMode)
107 {
108 m_npointsZ = 1;
109 }
110 else if (m_multipleModes)
111 {
112 m_npointsZ = m_session->GetParameter("HomModesZ");
113 }
114 }
115
116 if ((HomoStr == "HOMOGENEOUS2D") || (HomoStr == "Homogeneous2D") ||
117 (HomoStr == "2D") || (HomoStr == "Homo2D"))
118 {
120 m_session->LoadParameter("HomModesY", m_npointsY);
121 m_session->LoadParameter("LY", m_LhomY);
122 m_session->LoadParameter("HomModesZ", m_npointsZ);
123 m_session->LoadParameter("LZ", m_LhomZ);
124 m_HomoDirec = 2;
125 }
126
127 if ((HomoStr == "HOMOGENEOUS3D") || (HomoStr == "Homogeneous3D") ||
128 (HomoStr == "3D") || (HomoStr == "Homo3D"))
129 {
131 m_session->LoadParameter("HomModesX", m_npointsX);
132 m_session->LoadParameter("LX", m_LhomX);
133 m_session->LoadParameter("HomModesY", m_npointsY);
134 m_session->LoadParameter("LY", m_LhomY);
135 m_session->LoadParameter("HomModesZ", m_npointsZ);
136 m_session->LoadParameter("LZ", m_LhomZ);
137 m_HomoDirec = 3;
138 }
139
140 if (m_session->DefinesSolverInfo("USEFFT"))
141 {
142 m_useFFT = true;
143 }
144 }
145 else
146 {
147 m_npointsZ = 1; // set to default value so can use to identify 2d or 3D
148 // (homogeneous) expansions
149 }
150
151 size_t nvar = m_session->GetVariables().size();
153 for (size_t i = 0; i < nvar; ++i)
154 {
155 m_baseflow[i] = Array<OneD, NekDouble>(pFields[i]->GetTotPoints(), 0.0);
156 }
157
158 size_t nBaseDerivs = (m_halfMode || m_singleMode) ? 2 : m_spacedim;
159 m_gradBase = Array<OneD, Array<OneD, NekDouble>>(nvar * nBaseDerivs);
160 for (size_t i = 0; i < nvar; ++i)
161 {
162 for (size_t j = 0; j < nBaseDerivs; ++j)
163 {
164 m_gradBase[i * nBaseDerivs + j] =
165 Array<OneD, NekDouble>(pFields[i]->GetTotPoints(), 0.0);
166 }
167 }
168
169 ASSERTL0(m_session->DefinesFunction("BaseFlow"),
170 "Base flow must be defined for linearised forms.");
171 std::string file = m_session->GetFunctionFilename("BaseFlow", 0);
172
173 // Periodic base flows
174 if (m_session->DefinesParameter("N_slices"))
175 {
176 m_session->LoadParameter("N_slices", m_slices);
177 if (m_slices > 1)
178 {
179 ASSERTL0(m_session->GetFunctionType("BaseFlow", 0) ==
181 "Base flow should be a sequence of files.");
182 m_session->LoadParameter("BaseFlow_interporder", m_interporder, 0);
184 "BaseFlow_interporder should be smaller than or equal to "
185 "N_slices.");
187 m_session->LoadParameter("N_start", m_start, 0);
188 m_session->LoadParameter("N_skip", m_skip, 1);
189 DFT(file, pFields);
190 }
191 else
192 {
193 ASSERTL0(false, "Number of slices must be a positive number "
194 "greater than 1");
195 }
196 }
197 // Steady base-flow
198 else
199 {
200 m_slices = 1;
201
202 // BaseFlow from file
203 if (m_session->GetFunctionType("BaseFlow", m_session->GetVariable(0)) ==
205 {
206 ImportFldBase(file, pFields, 1);
207 }
208 // analytic base flow
209 else
210 {
211 size_t nq = pFields[0]->GetNpoints();
215
216 // get the coordinates (assuming all fields have the same
217 // discretisation)
218 pFields[0]->GetCoords(x0, x1, x2);
219
220 for (size_t i = 0; i < pFields.size(); i++)
221 {
223 m_session->GetFunction("BaseFlow", i);
224
225 ifunc->Evaluate(x0, x1, x2, m_baseflow[i]);
226 }
227 }
228 }
229
230 for (size_t i = 0; i < nvar; ++i)
231 {
232 UpdateGradBase(i, pFields[i]);
233 }
234
235 if (m_session->DefinesParameter("period"))
236 {
237 m_period = m_session->GetParameter("period");
238 }
239 else
240 {
241 m_period =
242 (m_session->GetParameter("TimeStep") * m_slices) / (m_slices - 1.);
243 }
244 if (m_session->GetComm()->GetRank() == 0)
245 {
246 std::cout << "baseflow info : interpolation order " << m_interporder
247 << ", period " << m_period << ", periodicity ";
248 if (m_isperiodic)
249 {
250 std::cout << "yes\n";
251 }
252 else
253 {
254 std::cout << "no\n";
255 }
256 std::cout << "baseflow info : files from " << m_start << " to "
257 << (m_start + (m_slices - 1) * m_skip) << " (skip " << m_skip
258 << ") with " << (m_slices - (m_interporder > 1))
259 << " time intervals" << std::endl;
260 }
261}
262
263// Advection function
265 const int nConvectiveFields,
267 const Array<OneD, Array<OneD, NekDouble>> &advVel,
268 const Array<OneD, Array<OneD, NekDouble>> &inarray,
269 Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble &time,
270 [[maybe_unused]] const Array<OneD, Array<OneD, NekDouble>> &pFwd,
271 [[maybe_unused]] const Array<OneD, Array<OneD, NekDouble>> &pBwd)
272{
273 ASSERTL1(nConvectiveFields == inarray.size(),
274 "Number of convective fields and Inarray are not compatible");
275
276 size_t nPointsTot = fields[0]->GetNpoints();
277 size_t ndim = advVel.size();
278 size_t nBaseDerivs = (m_halfMode || m_singleMode) ? 2 : m_spacedim;
279 size_t nDerivs = (m_halfMode) ? 2 : m_spacedim;
280
282 for (size_t i = 0; i < ndim; ++i)
283 {
284 if (fields[i]->GetWaveSpace() && !m_singleMode && !m_halfMode)
285 {
286 velocity[i] = Array<OneD, NekDouble>(nPointsTot, 0.0);
287 fields[i]->HomogeneousBwdTrans(nPointsTot, advVel[i], velocity[i]);
288 }
289 else
290 {
291 velocity[i] = advVel[i];
292 }
293 }
294
296 for (size_t i = 0; i < nDerivs; ++i)
297 {
298 grad[i] = Array<OneD, NekDouble>(nPointsTot);
299 }
300
301 // Evaluation of the base flow for periodic cases
302 if (m_slices > 1)
303 {
304 for (size_t i = 0; i < ndim; ++i)
305 {
306 UpdateBase(m_interp[i], m_baseflow[i], time);
307 UpdateGradBase(i, fields[i]);
308 }
309 }
310
311 // Evaluate the linearised advection term
312 for (size_t i = 0; i < (size_t)nConvectiveFields; ++i)
313 {
314 // Calculate gradient
315 switch (nDerivs)
316 {
317 case 1:
318 {
319 fields[i]->PhysDeriv(inarray[i], grad[0]);
320 }
321 break;
322 case 2:
323 {
324 fields[i]->PhysDeriv(inarray[i], grad[0], grad[1]);
325 }
326 break;
327 case 3:
328 {
329 fields[i]->PhysDeriv(inarray[i], grad[0], grad[1], grad[2]);
330 if (m_multipleModes)
331 {
332 // transform gradients into physical Fourier space
333 fields[i]->HomogeneousBwdTrans(nPointsTot, grad[0],
334 grad[0]);
335 fields[i]->HomogeneousBwdTrans(nPointsTot, grad[1],
336 grad[1]);
337 fields[i]->HomogeneousBwdTrans(nPointsTot, grad[2],
338 grad[2]);
339 }
340 }
341 break;
342 }
343
344 // Calculate U_j du'_i/dx_j
345 Vmath::Vmul(nPointsTot, grad[0], 1, m_baseflow[0], 1, outarray[i], 1);
346 for (size_t j = 1; j < nDerivs; ++j)
347 {
348 Vmath::Vvtvp(nPointsTot, grad[j], 1, m_baseflow[j], 1, outarray[i],
349 1, outarray[i], 1);
350 }
351
352 // Add u'_j dU_i/dx_j
353 size_t lim = (m_halfMode || m_singleMode) ? 2 : ndim;
354 if (m_halfMode && i == 2)
355 {
356 lim = 0;
357 }
358 for (size_t j = 0; j < lim; ++j)
359 {
360 Vmath::Vvtvp(nPointsTot, m_gradBase[i * nBaseDerivs + j], 1,
361 velocity[j], 1, outarray[i], 1, outarray[i], 1);
362 }
363
364 if (m_multipleModes)
365 {
366 fields[i]->HomogeneousFwdTrans(nPointsTot, outarray[i],
367 outarray[i]);
368 }
369 Vmath::Neg(nPointsTot, outarray[i], 1);
370 }
371}
372
374 const Array<OneD, Array<OneD, NekDouble>> &inarray,
376{
377 if (m_session->GetSolverInfo("EqType") == "UnsteadyNavierStokes")
378 {
379 // The SFD method is only applied to the velocity variables in
380 // incompressible
381 ASSERTL1(inarray.size() == (m_baseflow.size() - 1),
382 "Number of base flow variables does not match what is "
383 "expected.");
384 }
385 else
386 {
387 ASSERTL1(
388 inarray.size() == (m_baseflow.size()),
389 "Number of base flow variables does not match what is expected.");
390 }
391
392 size_t npts = inarray[0].size();
393
394 for (size_t i = 0; i < inarray.size(); ++i)
395 {
396 ASSERTL1(npts == m_baseflow[i].size(),
397 "Size of base flow array does not match expected.");
398 Vmath::Vcopy(npts, inarray[i], 1, m_baseflow[i], 1);
399 UpdateGradBase(i, fields[i]);
400 }
401}
402
403/**
404 * Import field from infile and load into \a m_fields. This routine will
405 * also perform a \a BwdTrans to ensure data is in both the physical and
406 * coefficient storage.
407 * @param pInFile Filename to read.
408 * @param pFields Array of expansion lists
409 */
411 std::string pInfile, Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
412 int pSlice)
413{
414 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
415 std::vector<std::vector<NekDouble>> FieldData;
416
417 size_t nqtot = m_baseflow[0].size();
418 Array<OneD, NekDouble> tmp_coeff(pFields[0]->GetNcoeffs(), 0.0);
419
420 size_t numexp = pFields[0]->GetExpSize();
421 Array<OneD, int> ElementGIDs(numexp);
422
423 // Define list of global element ids
424 for (size_t i = 0; i < numexp; ++i)
425 {
426 ElementGIDs[i] = pFields[0]->GetExp(i)->GetGeom()->GetGlobalID();
427 }
428
429 // Get Homogeneous
432 fld->Import(pInfile, FieldDef, FieldData,
434
435 size_t nSessionVar = m_session->GetVariables().size();
436 size_t nSessionConvVar = nSessionVar - 1;
437 size_t nFileVar = FieldDef[0]->m_fields.size();
438
439 std::unordered_map<int, int> zIdToPlane;
441 {
442 zIdToPlane[0] = 0;
443 }
444
445 for (size_t j = 0; j < nFileVar; ++j)
446 {
447 size_t k = 0;
448 for (; k < nSessionConvVar; ++k)
449 {
450 if (m_session->GetVariable(k) == FieldDef[0]->m_fields[j])
451 {
452 break;
453 }
454 }
455 if (k == nSessionConvVar)
456 {
457 continue;
458 }
459 for (size_t i = 0; i < FieldDef.size(); ++i)
460 {
461 pFields[j]->ExtractDataToCoeffs(FieldDef[i], FieldData[i],
462 FieldDef[i]->m_fields[j], tmp_coeff,
463 zIdToPlane);
464 }
465
467 {
468 pFields[j]->GetPlane(0)->BwdTrans(tmp_coeff, m_baseflow[k]);
469
470 if (m_singleMode)
471 {
472 // copy the bwd trans into the second plane for single
473 // Mode Analysis
474 int ncplane = (pFields[0]->GetNpoints()) / m_npointsZ;
475 Vmath::Vcopy(ncplane, &m_baseflow[k][0], 1,
476 &m_baseflow[k][ncplane], 1);
477 }
478 }
479 else // fully 3D base flow - put in physical space.
480 {
481 bool oldwavespace = pFields[j]->GetWaveSpace();
482 pFields[j]->SetWaveSpace(false);
483 pFields[j]->BwdTrans(tmp_coeff, m_baseflow[k]);
484 pFields[j]->SetWaveSpace(oldwavespace);
485 }
486 }
487
488 // If time-periodic, put loaded data into the slice storage.
489 if (m_slices > 1)
490 {
491 for (size_t i = 0; i < nSessionConvVar; ++i)
492 {
493 Vmath::Vcopy(nqtot, &m_baseflow[i][0], 1,
494 &m_interp[i][pSlice * nqtot], 1);
495 }
496 }
497}
498
500 const Array<OneD, const NekDouble> &inarray,
501 Array<OneD, NekDouble> &outarray, const NekDouble time)
502{
503 int npoints = m_baseflow[0].size();
504 if (m_isperiodic)
505 {
506 NekDouble BetaT = 2 * M_PI * fmod(time, m_period) / m_period;
507 NekDouble phase;
508 Array<OneD, NekDouble> auxiliary(npoints);
509
510 Vmath::Vcopy(npoints, &inarray[0], 1, &outarray[0], 1);
511 Vmath::Svtvp(npoints, cos(0.5 * m_slices * BetaT), &inarray[npoints], 1,
512 &outarray[0], 1, &outarray[0], 1);
513
514 for (int i = 2; i < m_slices; i += 2)
515 {
516 phase = (i >> 1) * BetaT;
517
518 Vmath::Svtvp(npoints, cos(phase), &inarray[i * npoints], 1,
519 &outarray[0], 1, &outarray[0], 1);
520 Vmath::Svtvp(npoints, -sin(phase), &inarray[(i + 1) * npoints], 1,
521 &outarray[0], 1, &outarray[0], 1);
522 }
523 }
524 else
525 {
526 NekDouble x = time;
527 x = x / m_period * (m_slices - 1);
528 int ix = x;
529 if (ix < 0)
530 {
531 ix = 0;
532 }
533 if (ix > m_slices - 2)
534 {
535 ix = m_slices - 2;
536 }
537 int padleft = m_interporder / 2 - 1;
538 if (padleft > ix)
539 {
540 padleft = ix;
541 }
542 int padright = m_interporder - 1 - padleft;
543 if (padright > m_slices - 1 - ix)
544 {
545 padright = m_slices - 1 - ix;
546 }
547 padleft = m_interporder - 1 - padright;
549 for (int i = 0; i < m_interporder; ++i)
550 {
551 for (int j = 0; j < m_interporder; ++j)
552 {
553 if (i != j)
554 {
555 coeff[i] *= (x - ix + padleft - (NekDouble)j) /
556 ((NekDouble)i - (NekDouble)j);
557 }
558 }
559 }
560 Vmath::Zero(npoints, &outarray[0], 1);
561 for (int i = ix - padleft; i < ix + padright + 1; ++i)
562 {
563 Vmath::Svtvp(npoints, coeff[i - ix + padleft],
564 &inarray[i * npoints], 1, &outarray[0], 1,
565 &outarray[0], 1);
566 }
567 }
568}
569
571 const int var, const MultiRegions::ExpListSharedPtr &field)
572{
573 int nBaseDerivs = (m_halfMode || m_singleMode) ? 2 : m_spacedim;
574 switch (m_spacedim)
575 {
576 case 1: // 1D
577 {
578 field->PhysDeriv(m_baseflow[var],
579 m_gradBase[var * nBaseDerivs + 0]);
580 }
581 break;
582 case 2: // 2D
583 {
584 field->PhysDeriv(m_baseflow[var], m_gradBase[var * nBaseDerivs + 0],
585 m_gradBase[var * nBaseDerivs + 1]);
586 }
587 break;
588 case 3:
589 {
590 if (m_halfMode) // can assume W = 0 and d/dz == 0
591 {
592 if (var < 2)
593 {
594 field->PhysDeriv(m_baseflow[var],
595 m_gradBase[var * nBaseDerivs + 0],
596 m_gradBase[var * nBaseDerivs + 1]);
597 }
598 }
599 else if (m_singleMode) // single mode where d/dz = 0
600 {
601 field->PhysDeriv(m_baseflow[var],
602 m_gradBase[var * nBaseDerivs + 0],
603 m_gradBase[var * nBaseDerivs + 1]);
604 }
605 else
606 {
607 // Differentiate base flow in physical space
608 bool oldwavespace = field->GetWaveSpace();
609 field->SetWaveSpace(false);
610 field->PhysDeriv(m_baseflow[var],
611 m_gradBase[var * nBaseDerivs + 0],
612 m_gradBase[var * nBaseDerivs + 1],
613 m_gradBase[var * nBaseDerivs + 2]);
614 field->SetWaveSpace(oldwavespace);
615 }
616 }
617 break;
618 }
619}
620
622 [[maybe_unused]] FloquetMatType mattype,
623 [[maybe_unused]] bool UseContCoeffs) const
624{
625 DNekMatSharedPtr loc_mat;
626 DNekBlkMatSharedPtr BlkMatrix;
627 size_t n_exp = 0;
628
629 n_exp = m_baseflow[0].size(); // will operatore on m_phys
630
631 Array<OneD, unsigned int> nrows(n_exp);
632 Array<OneD, unsigned int> ncols(n_exp);
633
634 nrows = Array<OneD, unsigned int>(n_exp, m_slices);
635 ncols = Array<OneD, unsigned int>(n_exp, m_slices);
636
637 MatrixStorage blkmatStorage = eDIAGONAL;
638 BlkMatrix = MemoryManager<DNekBlkMat>::AllocateSharedPtr(nrows, ncols,
639 blkmatStorage);
640
644 StdRegions::StdSegExp StdSeg(BK);
645
647 StdSeg.DetShapeType(), StdSeg);
648
649 loc_mat = StdSeg.GetStdMatrix(matkey);
650
651 // set up array of block matrices.
652 for (size_t i = 0; i < n_exp; ++i)
653 {
654 BlkMatrix->SetBlock(i, i, loc_mat);
655 }
656
657 return BlkMatrix;
658}
659
660// Discrete Fourier Transform for Floquet analysis
662 const std::string file,
664{
665 size_t ConvectedFields = m_baseflow.size() - 1;
666 size_t npoints = m_baseflow[0].size();
668
669 for (size_t i = 0; i < ConvectedFields; ++i)
670 {
671 m_interp[i] = Array<OneD, NekDouble>(npoints * m_slices, 0.0);
672 }
673
674 // Import the slides into auxiliary vector
675 // The base flow should be stored in the form "filename_%d.ext"
676 // A subdirectory can also be included, such as "dir/filename_%d.ext"
677 size_t found = file.find("%d");
678 ASSERTL0(found != std::string::npos &&
679 file.find("%d", found + 1) == std::string::npos,
680 "Since N_slices is specified, the filename provided for function "
681 "'BaseFlow' must include exactly one instance of the format "
682 "specifier '%d', to index the time-slices.");
683 size_t nstart = m_start;
684 for (size_t i = nstart; i < nstart + m_slices * m_skip; i += m_skip)
685 {
686 boost::format filename(file);
687 filename % i;
688 ImportFldBase(filename.str(), pFields, (i - nstart) / m_skip);
689 if (m_session->GetComm()->GetRank() == 0)
690 {
691 std::cout << "read base flow file " << filename.str() << std::endl;
692 }
693 }
694 if (!m_isperiodic)
695 {
696 return;
697 }
698
699 // Discrete Fourier Transform of the fields
700 for (size_t k = 0; k < ConvectedFields; ++k)
701 {
702#ifdef NEKTAR_USING_FFTW
703
704 // Discrete Fourier Transform using FFTW
705 Array<OneD, NekDouble> fft_in(npoints * m_slices);
706 Array<OneD, NekDouble> fft_out(npoints * m_slices);
707
710
711 // Shuffle the data
712 for (int j = 0; j < m_slices; ++j)
713 {
714 Vmath::Vcopy(npoints, &m_interp[k][j * npoints], 1, &(fft_in[j]),
715 m_slices);
716 }
717
719 m_slices);
720
721 // FFT Transform
722 for (size_t i = 0; i < npoints; i++)
723 {
724 m_FFT->FFTFwdTrans(m_tmpIN = fft_in + i * m_slices,
725 m_tmpOUT = fft_out + i * m_slices);
726 }
727
728 // Reshuffle data
729 for (int s = 0; s < m_slices; ++s)
730 {
731 Vmath::Vcopy(npoints, &fft_out[s], m_slices,
732 &m_interp[k][s * npoints], 1);
733 }
734
735 Vmath::Zero(fft_in.size(), &fft_in[0], 1);
736 Vmath::Zero(fft_out.size(), &fft_out[0], 1);
737#else
738 // Discrete Fourier Transform using MVM
739 DNekBlkMatSharedPtr blkmat;
741
742 int nrows = blkmat->GetRows();
743 int ncols = blkmat->GetColumns();
744
745 Array<OneD, NekDouble> sortedinarray(ncols);
746 Array<OneD, NekDouble> sortedoutarray(nrows);
747
748 // Shuffle the data
749 for (int j = 0; j < m_slices; ++j)
750 {
751 Vmath::Vcopy(npoints, &m_interp[k][j * npoints], 1,
752 &(sortedinarray[j]), m_slices);
753 }
754
755 // Create NekVectors from the given data arrays
756 NekVector<NekDouble> in(ncols, sortedinarray, eWrapper);
757 NekVector<NekDouble> out(nrows, sortedoutarray, eWrapper);
758
759 // Perform matrix-vector multiply.
760 out = (*blkmat) * in;
761
762 // Reshuffle data
763 for (int s = 0; s < m_slices; ++s)
764 {
765 Vmath::Vcopy(npoints, &sortedoutarray[s], m_slices,
766 &m_interp[k][s * npoints], 1);
767 }
768
769 for (size_t r = 0; r < sortedinarray.size(); ++r)
770 {
771 sortedinarray[0] = 0;
772 sortedoutarray[0] = 0;
773 }
774
775#endif
776 }
777}
778
779} // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:242
Describes the specification for a Basis.
Definition: Basis.h:45
static std::shared_ptr< FieldIO > CreateForFile(const LibUtilities::SessionReaderSharedPtr session, const std::string &filename)
Construct a FieldIO object for a given input filename.
Definition: FieldIO.cpp:223
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.
Defines a specification for a set of points.
Definition: Points.h:50
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
int m_npointsZ
number of points in Z direction (if homogeneous)
NekDouble m_LhomZ
physical length in Z direction (if homogeneous)
void v_InitObject(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields) override
Initialises the advection object.
enum HomogeneousType m_HomogeneousType
void UpdateBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const NekDouble time)
NekDouble m_period
period length
void v_Advect(const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &advVel, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble &time, const Array< OneD, Array< OneD, NekDouble > > &pFwd=NullNekDoubleArrayOfArray, const Array< OneD, Array< OneD, NekDouble > > &pBwd=NullNekDoubleArrayOfArray) override
Advects a vector field.
int m_npointsX
number of points in X direction (if homogeneous)
int m_npointsY
number of points in Y direction (if homogeneous)
int m_HomoDirec
number of homogenous directions
static SolverUtils::AdvectionSharedPtr create(std::string)
Creates an instance of this class.
bool m_useFFT
flag to determine if use or not the FFT for transformations
void v_SetBaseFlow(const Array< OneD, Array< OneD, NekDouble > > &inarray, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields) override
Overrides the base flow used during linearised advection.
LibUtilities::SessionReaderSharedPtr m_session
Array< OneD, NekDouble > m_tmpIN
void UpdateGradBase(const int var, const MultiRegions::ExpListSharedPtr &field)
DNekBlkMatSharedPtr GetFloquetBlockMatrix(FloquetMatType mattype, bool UseContCoeffs=false) const
bool m_singleMode
flag to determine if use single mode or not
Array< OneD, Array< OneD, NekDouble > > m_baseflow
Storage for base flow.
Array< OneD, NekDouble > m_tmpOUT
Array< OneD, Array< OneD, NekDouble > > m_gradBase
bool m_multipleModes
flag to determine if use multiple mode or not
void DFT(const std::string file, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
LibUtilities::NektarFFTSharedPtr m_FFT
auxiliary variables
bool m_halfMode
flag to determine if use half mode or not
NekDouble m_LhomY
physical length in Y direction (if homogeneous)
void ImportFldBase(std::string pInfile, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, int slice)
NekDouble m_LhomX
physical length in X direction (if homogeneous)
Array< OneD, Array< OneD, NekDouble > > m_interp
interpolation vector
static std::string className
Name of class.
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
An abstract base class encapsulating the concept of advection of a vector field.
Definition: Advection.h:81
virtual SOLVER_UTILS_EXPORT void v_InitObject(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Initialises the advection object.
Definition: Advection.cpp:264
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:612
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition: StdExpansion.h:370
Class representing a segment element in reference space All interface of this class sits in StdExpans...
Definition: StdSegExp.h:45
NektarFFTFactory & GetNektarFFTFactory()
Definition: NektarFFT.cpp:65
std::shared_ptr< FieldIO > FieldIOSharedPtr
Definition: FieldIO.h:322
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< Equation > EquationSharedPtr
Definition: Equation.h:125
static FieldMetaDataMap NullFieldMetaDataMap
Definition: FieldIO.h:51
@ eFourierEvenlySpaced
1D Evenly-spaced points using Fourier Fit
Definition: PointsType.h:74
@ eFourier
Fourier Expansion .
Definition: BasisType.h:55
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
Definition: Advection.cpp:43
std::shared_ptr< DNekBlkMat > DNekBlkMatSharedPtr
Definition: NekTypeDefs.hpp:77
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:75
double NekDouble
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.hpp:72
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 Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.hpp:292
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.hpp:366
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.hpp:273
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825