Nektar++
AdjointAdvection.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File AdjointAdvection.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 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Evaluation of the adjoint advective term
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
37 #include <StdRegions/StdSegExp.h>
38 
48 
49 
50 namespace Nektar
51 {
52 
56 
57 /**
58  *
59  */
61  : Advection()
62 {
63 }
64 
65 
69 {
70 
71  Advection::v_InitObject(pSession, pFields);
72 
73  m_session = pSession;
75  ::AllocateSharedPtr(m_session, pFields[0]->GetGraph());
76  m_spacedim = pFields[0]->GetGraph()->GetSpaceDimension();
77  m_expdim = pFields[0]->GetGraph()->GetMeshDimension();
79 
80  //Setting parameters for homogeneous problems
81  m_HomoDirec = 0;
82  m_useFFT = false;
84  m_SingleMode = false;
85  m_HalfMode = false;
86  m_MultipleModes = false;
87  m_homogen_dealiasing = false;
88 
89  if(m_session->DefinesSolverInfo("HOMOGENEOUS"))
90  {
91  std::string HomoStr = m_session->GetSolverInfo("HOMOGENEOUS");
92  m_spacedim = 3;
93 
94  if((HomoStr == "HOMOGENEOUS1D")||(HomoStr == "Homogeneous1D")||
95  (HomoStr == "1D")||(HomoStr == "Homo1D"))
96  {
98  m_LhomZ = m_session->GetParameter("LZ");
99  m_HomoDirec = 1;
100 
101  m_homogen_dealiasing = pSession->DefinesSolverInfo("DEALIASING");
102 
103  if(m_session->DefinesSolverInfo("ModeType"))
104  {
105  m_session->MatchSolverInfo("ModeType","SingleMode",m_SingleMode,false);
106  m_session->MatchSolverInfo("ModeType","HalfMode",m_HalfMode,false);
107  m_session->MatchSolverInfo("ModeType","MultipleModes",m_MultipleModes,false);
108  }
109  if(m_session->DefinesSolverInfo("ModeType"))
110  {
111  if(m_SingleMode)
112  {
113  m_npointsZ=2;
114  }
115  else if(m_HalfMode)
116  {
117  m_npointsZ=1;
118  }
119  else if(m_MultipleModes)
120  {
121  m_npointsZ = m_session->GetParameter("HomModesZ");
122  }
123  else
124  {
125  ASSERTL0(false, "SolverInfo ModeType not valid");
126  }
127  }
128  else
129  {
130  m_session->LoadParameter("HomModesZ",m_npointsZ);
131 
132  }
133 
134  }
135 
136  if((HomoStr == "HOMOGENEOUS2D")||(HomoStr == "Homogeneous2D")||
137  (HomoStr == "2D")||(HomoStr == "Homo2D"))
138  {
140  m_session->LoadParameter("HomModesY", m_npointsY);
141  m_session->LoadParameter("LY", m_LhomY);
142  m_session->LoadParameter("HomModesZ", m_npointsZ);
143  m_session->LoadParameter("LZ", m_LhomZ);
144  m_HomoDirec = 2;
145  }
146 
147  if((HomoStr == "HOMOGENEOUS3D")||(HomoStr == "Homogeneous3D")||
148  (HomoStr == "3D")||(HomoStr == "Homo3D"))
149  {
151  m_session->LoadParameter("HomModesX",m_npointsX);
152  m_session->LoadParameter("LX", m_LhomX );
153  m_session->LoadParameter("HomModesY",m_npointsY);
154  m_session->LoadParameter("LY", m_LhomY );
155  m_session->LoadParameter("HomModesZ",m_npointsZ);
156  m_session->LoadParameter("LZ", m_LhomZ );
157  m_HomoDirec = 3;
158  }
159 
160  if(m_session->DefinesSolverInfo("USEFFT"))
161  {
162  m_useFFT = true;
163  }
164  }
165  else
166  {
167  m_npointsZ = 1; // set to default value so can use to identify 2d or 3D (homogeneous) expansions
168  }
169 
170  if(m_session->DefinesSolverInfo("PROJECTION"))
171  {
172  std::string ProjectStr
173  = m_session->GetSolverInfo("PROJECTION");
174 
175  if((ProjectStr == "Continuous")||(ProjectStr == "Galerkin")||
176  (ProjectStr == "CONTINUOUS")||(ProjectStr == "GALERKIN"))
177  {
179  }
180  else if(ProjectStr == "DisContinuous")
181  {
183  }
184  else
185  {
186  ASSERTL0(false,"PROJECTION value not recognised");
187  }
188  }
189  else
190  {
191  cerr << "Projection type not specified in SOLVERINFO,"
192  "defaulting to continuous Galerkin" << endl;
194  }
195 
196  int nvar = m_session->GetVariables().size();
198  for (int i = 0; i < nvar; ++i)
199  {
200  m_baseflow[i] = Array<OneD, NekDouble>(pFields[i]->GetTotPoints(), 0.0);
201  }
202  //SetUpBaseFields(pFields[0]->GetGraph());
203  ASSERTL0(m_session->DefinesFunction("BaseFlow"),
204  "Base flow must be defined for linearised forms.");
205  string file = m_session->GetFunctionFilename("BaseFlow", 0);
206 
207 
208  //Periodic base flows
209  if(m_session->DefinesParameter("N_slices"))
210  {
211  m_session->LoadParameter("N_slices",m_slices);
212  if(m_slices>1)
213  {
214  DFT(file,pFields,m_slices);
215  }
216  else
217  {
218  ASSERTL0(false,"Number of slices must be a positive number greater than 1");
219  }
220  }
221  //Steady base-flow
222  else
223  {
224  m_slices=1;
225 
226  //BaseFlow from file
227  if (m_session->GetFunctionType("BaseFlow", m_session->GetVariable(0))
229  {
230  ImportFldBase(file,pFields,1);
231 
232  }
233  //analytic base flow
234  else
235  {
236  int nq = pFields[0]->GetNpoints();
237  Array<OneD,NekDouble> x0(nq);
238  Array<OneD,NekDouble> x1(nq);
239  Array<OneD,NekDouble> x2(nq);
240 
241  // get the coordinates (assuming all fields have the same
242  // discretisation)
243  pFields[0]->GetCoords(x0,x1,x2);
244  for(unsigned int i = 0 ; i < m_baseflow.num_elements(); i++)
245  {
247  = m_session->GetFunction("BaseFlow", i);
248 
249  ifunc->Evaluate(x0,x1,x2,m_baseflow[i]);
250  }
251 
252  }
253 
254  }
255 }
256 
257 
259 {
260 }
261 
262 
264  const int nConvectiveFields,
266  const Array<OneD, Array<OneD, NekDouble> > &advVel,
267  const Array<OneD, Array<OneD, NekDouble> > &inarray,
268  Array<OneD, Array<OneD, NekDouble> > &outarray,
269  const NekDouble &time)
270 {
271  int nqtot = fields[0]->GetTotPoints();
272  ASSERTL1(nConvectiveFields == inarray.num_elements(),"Number of convective fields and Inarray are not compatible");
273 
274  Array<OneD, NekDouble > Deriv = Array<OneD, NekDouble> (nqtot*nConvectiveFields);
275 
276  for(int n = 0; n < nConvectiveFields; ++n)
277  {
278  //v_ComputeAdvectionTerm(fields,advVel,inarray[i],outarray[i],i,time,Deriv);
279  int ndim = advVel.num_elements();
280  int nPointsTot = fields[0]->GetNpoints();
281  Array<OneD, NekDouble> grad0,grad1,grad2;
282 
283  //Evaluation of the gradiend of each component of the base flow
284  //\nabla U
285  Array<OneD, NekDouble> grad_base_u0,grad_base_u1,grad_base_u2;
286  // \nabla V
287  Array<OneD, NekDouble> grad_base_v0,grad_base_v1,grad_base_v2;
288  // \nabla W
289  Array<OneD, NekDouble> grad_base_w0,grad_base_w1,grad_base_w2;
290 
291 
292  grad0 = Array<OneD, NekDouble> (nPointsTot);
293  grad_base_u0 = Array<OneD, NekDouble> (nPointsTot);
294  grad_base_v0 = Array<OneD, NekDouble> (nPointsTot);
295  grad_base_w0 = Array<OneD, NekDouble> (nPointsTot);
296 
297 
298  //Evaluation of the base flow for periodic cases
299  //(it requires fld files)
300  if(m_slices>1)
301  {
302  if (m_session->GetFunctionType("BaseFlow", 0)
304  {
305  for(int i=0; i<ndim;++i)
306  {
308  }
309  }
310  else
311  {
312  ASSERTL0(false, "Periodic Base flow requires .fld files");
313  }
314  }
315 
316  //Evaluate the Adjoint advection term
317  switch(ndim)
318  {
319  // 1D
320  case 1:
321  fields[0]->PhysDeriv(inarray[n],grad0);
322  fields[0]->PhysDeriv(m_baseflow[0],grad_base_u0);
323  //Evaluate U du'/dx
324  Vmath::Vmul(nPointsTot,grad0,1,m_baseflow[0],1,outarray[n],1);
325  //Evaluate U du'/dx+ u' dU/dx
326  Vmath::Vvtvp(nPointsTot,grad_base_u0,1,advVel[0],1,outarray[n],1,outarray[n],1);
327  break;
328 
329  //2D
330  case 2:
331 
332  grad1 = Array<OneD, NekDouble> (nPointsTot);
333  grad_base_u1 = Array<OneD, NekDouble> (nPointsTot);
334  grad_base_v1 = Array<OneD, NekDouble> (nPointsTot);
335 
336  fields[0]->PhysDeriv(inarray[n],grad0,grad1);
337 
338  //Derivates of the base flow
339  fields[0]-> PhysDeriv(m_baseflow[0], grad_base_u0, grad_base_u1);
340  fields[0]-> PhysDeriv(m_baseflow[1], grad_base_v0, grad_base_v1);
341 
342  //Since the components of the velocity are passed one by one, it is necessary to distinguish which
343  //term is consider
344  switch (n)
345  {
346  //x-equation
347  case 0:
348  // Evaluate U du'/dx
349  Vmath::Vmul (nPointsTot,grad0,1,m_baseflow[0],1,outarray[n],1);
350  //Evaluate U du'/dx+ V du'/dy
351  Vmath::Vvtvp(nPointsTot,grad1,1,m_baseflow[1],1,outarray[n],1,outarray[n],1);
352  //Evaluate - (U du'/dx+ V du'/dy)
353  Vmath::Neg(nPointsTot,outarray[n],1);
354  //Evaluate -(U du'/dx+ V du'/dy)+u' dU/dx
355  Vmath::Vvtvp(nPointsTot,grad_base_u0,1,advVel[0],1,outarray[n],1,outarray[n],1);
356  //Evaluate -(U du'/dx+ V du'/dy) +u' dU/dx +v' dV/dx
357  Vmath::Vvtvp(nPointsTot,grad_base_v0,1,advVel[1],1,outarray[n],1,outarray[n],1);
358  break;
359 
360  //y-equation
361  case 1:
362  // Evaluate U dv'/dx
363  Vmath::Vmul (nPointsTot,grad0,1,m_baseflow[0],1,outarray[n],1);
364  //Evaluate U dv'/dx+ V dv'/dy
365  Vmath::Vvtvp(nPointsTot,grad1,1,m_baseflow[1],1,outarray[n],1,outarray[n],1);
366  //Evaluate -(U dv'/dx+ V dv'/dy)
367  Vmath::Neg(nPointsTot,outarray[n],1);
368  //Evaluate (U dv'/dx+ V dv'/dy)+u' dU/dy
369  Vmath::Vvtvp(nPointsTot,grad_base_u1,1,advVel[0],1,outarray[n],1,outarray[n],1);
370  //Evaluate (U dv'/dx+ V dv'/dy +u' dv/dx)+v' dV/dy
371  Vmath::Vvtvp(nPointsTot,grad_base_v1,1,advVel[1],1,outarray[n],1,outarray[n],1);
372  break;
373  }
374  break;
375 
376 
377  //3D
378  case 3:
379 
380  grad1 = Array<OneD, NekDouble> (nPointsTot);
381  grad2 = Array<OneD, NekDouble> (nPointsTot);
382  grad_base_u1 = Array<OneD, NekDouble> (nPointsTot);
383  grad_base_v1 = Array<OneD, NekDouble> (nPointsTot);
384  grad_base_w1 = Array<OneD, NekDouble> (nPointsTot);
385 
386  grad_base_u2 = Array<OneD, NekDouble> (nPointsTot);
387  grad_base_v2 = Array<OneD, NekDouble> (nPointsTot);
388  grad_base_w2 = Array<OneD, NekDouble> (nPointsTot);
389 
390  fields[0]->PhysDeriv(m_baseflow[0], grad_base_u0, grad_base_u1,grad_base_u2);
391  fields[0]->PhysDeriv(m_baseflow[1], grad_base_v0, grad_base_v1,grad_base_v2);
392  fields[0]->PhysDeriv(m_baseflow[2], grad_base_w0, grad_base_w1, grad_base_w2);
393 
394  //HalfMode has W(x,y,t)=0
395  if(m_HalfMode || m_SingleMode)
396  {
397  for(int i=0; i<grad_base_u2.num_elements();++i)
398  {
399  grad_base_u2[i]=0;
400  grad_base_v2[i]=0;
401  grad_base_w2[i]=0;
402 
403  }
404  }
405 
406  fields[0]->PhysDeriv(inarray[n], grad0, grad1, grad2);
407 
408  switch (n)
409  {
410  //x-equation
411  case 0:
412  //Evaluate U du'/dx
413  Vmath::Vmul (nPointsTot,grad0,1,m_baseflow[0],1,outarray[n],1);
414  //Evaluate U du'/dx+ V du'/dy
415  Vmath::Vvtvp(nPointsTot,grad1,1,m_baseflow[1],1,outarray[n],1,outarray[n],1);
416  //Evaluate U du'/dx+ V du'/dy+W du'/dz
417  Vmath::Vvtvp(nPointsTot,grad2,1,m_baseflow[2],1,outarray[n],1,outarray[n],1);
418  //Evaluate -(U du'/dx+ V du'/dy+W du'/dz)
419  Vmath::Neg(nPointsTot,outarray[n],1);
420  //Evaluate -(U du'/dx+ V du'/dy+W du'/dz)+u' dU/dx
421  Vmath::Vvtvp(nPointsTot,grad_base_u0,1,advVel[0],1,outarray[n],1,outarray[n],1);
422  //Evaluate -(U du'/dx+ V du'/dy+W du'/dz)+u'dU/dx+ v' dV/dx
423  Vmath::Vvtvp(nPointsTot,grad_base_v0,1,advVel[1],1,outarray[n],1,outarray[n],1);
424  //Evaluate -(U du'/dx+ V du'/dy+W du'/dz)+u'dU/dx+ v' dV/dx+ w' dW/dz
425  Vmath::Vvtvp(nPointsTot,grad_base_w0,1,advVel[2],1,outarray[n],1,outarray[n],1);
426  break;
427  //y-equation
428  case 1:
429  //Evaluate U dv'/dx
430  Vmath::Vmul (nPointsTot,grad0,1,m_baseflow[0],1,outarray[n],1);
431  //Evaluate U dv'/dx+ V dv'/dy
432  Vmath::Vvtvp(nPointsTot,grad1,1,m_baseflow[1],1,outarray[n],1,outarray[n],1);
433  //Evaluate U dv'/dx+ V dv'/dy+W dv'/dz
434  Vmath::Vvtvp(nPointsTot,grad2,1,m_baseflow[2],1,outarray[n],1,outarray[n],1);
435  //Evaluate -(U dv'/dx+ V dv'/dy+W dv'/dz)
436  Vmath::Neg(nPointsTot,outarray[n],1);
437  //Evaluate -(U dv'/dx+ V dv'/dy+W dv'/dz)+u' dU/dy
438  Vmath::Vvtvp(nPointsTot,grad_base_u1,1,advVel[0],1,outarray[n],1,outarray[n],1);
439  //Evaluate -(U dv'/dx+ V dv'/dy+W dv'/dz)+u' dU/dy +v' dV/dy
440  Vmath::Vvtvp(nPointsTot,grad_base_v1,1,advVel[1],1,outarray[n],1,outarray[n],1);
441  //Evaluate -(U dv'/dx+ V dv'/dy+W dv'/dz)+u' dU/dy +v' dV/dy+ w' dW/dy
442  Vmath::Vvtvp(nPointsTot,grad_base_w1,1,advVel[2],1,outarray[n],1,outarray[n],1);
443  break;
444 
445  //z-equation
446  case 2:
447  //Evaluate U dw'/dx
448  Vmath::Vmul (nPointsTot,grad0,1,m_baseflow[0],1,outarray[n],1);
449  //Evaluate U dw'/dx+ V dw'/dx
450  Vmath::Vvtvp(nPointsTot,grad1,1,m_baseflow[1],1,outarray[n],1,outarray[n],1);
451  //Evaluate U dw'/dx+ V dw'/dx+ W dw'/dz
452  Vmath::Vvtvp(nPointsTot,grad2,1,m_baseflow[2],1,outarray[n],1,outarray[n],1);
453  //Evaluate -(U dw'/dx+ V dw'/dx+ W dw'/dz)
454  Vmath::Neg(nPointsTot,outarray[n],1);
455  //Evaluate -(U dw'/dx+ V dw'/dx+ W dw'/dz)+u' dU/dz
456  Vmath::Vvtvp(nPointsTot,grad_base_u2,1,advVel[0],1,outarray[n],1,outarray[n],1);
457  //Evaluate -(U dw'/dx+ V dw'/dx+ W dw'/dz)+u' dU/dz+v'dV/dz
458  Vmath::Vvtvp(nPointsTot,grad_base_v2,1,advVel[1],1,outarray[n],1,outarray[n],1);
459  //Evaluate -(U dw'/dx+ V dw'/dx+ W dw'/dz)+u' dU/dz+v'dV/dz + w' dW/dz
460  Vmath::Vvtvp(nPointsTot,grad_base_w2,1,advVel[2],1,outarray[n],1,outarray[n],1);
461  break;
462  }
463  break;
464  default:
465  ASSERTL0(false,"dimension unknown");
466  }
467 
468  Vmath::Neg(nqtot,outarray[n],1);
469  }
470 
471 }
472 
473 
475  const Array<OneD, Array<OneD, NekDouble> > &inarray)
476 {
477  ASSERTL1(inarray.num_elements() == m_baseflow.num_elements(),
478  "Number of base flow variables does not match what is"
479  "expected.");
480 
481  int npts = inarray[0].num_elements();
482  for (int i = 0; i < inarray.num_elements(); ++i)
483  {
484  ASSERTL1(npts == m_baseflow.num_elements(),
485  "Size of base flow array does not match expected.");
486  Vmath::Vcopy(npts, inarray[i], 1, m_baseflow[i], 1);
487  }
488 }
489 
490 
491 /**
492  * Import field from infile and load into \a m_fields. This routine will
493  * also perform a \a BwdTrans to ensure data is in both the physical and
494  * coefficient storage.
495  * @param infile Filename to read.
496  */
497 void AdjointAdvection::ImportFldBase(std::string pInfile,
499 {
500  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
501  std::vector<std::vector<NekDouble> > FieldData;
502 
503  int nqtot = pFields[0]->GetTotPoints();
504  Array<OneD, NekDouble> tmp_coeff(pFields[0]->GetNcoeffs(), 0.0);
505 
506  //Get Homogeneous
509  m_session->GetComm());
510  fld->Import(pInfile, FieldDef, FieldData);
511 
512  int nvar = m_session->GetVariables().size();
513  int s;
514 
515  if(m_session->DefinesSolverInfo("HOMOGENEOUS"))
516  {
517  std::string HomoStr = m_session->GetSolverInfo("HOMOGENEOUS");
518  }
519 
520  // copy FieldData into m_fields
521  for(int j = 0; j < nvar; ++j)
522  {
523  for(int i = 0; i < FieldDef.size(); ++i)
524  {
525  if ((m_session->DefinesSolverInfo("HOMOGENEOUS") &&
526  (m_session->GetSolverInfo("HOMOGENEOUS")=="HOMOGENEOUS1D" ||
527  m_session->GetSolverInfo("HOMOGENEOUS")=="1D" ||
528  m_session->GetSolverInfo("HOMOGENEOUS")=="Homo1D")) &&
529  nvar==4)
530  {
531  // w-component must be ignored and set to zero.
532  if (j != nvar - 2)
533  {
534  // p component (it is 4th variable of the 3D and corresponds 3nd variable of 2D)
535  s = (j == nvar - 1) ? 2 : j;
536 
537  //extraction of the 2D
538  pFields[j]->ExtractDataToCoeffs(
539  FieldDef[i],
540  FieldData[i],
541  FieldDef[i]->m_fields[s],
542  tmp_coeff);
543  }
544 
545  // Put zero on higher modes
546  int ncplane = (pFields[0]->GetNcoeffs()) / m_npointsZ;
547 
548  if (m_npointsZ > 2)
549  {
550  Vmath::Zero(ncplane*(m_npointsZ-2),
551  &tmp_coeff[2*ncplane],1);
552  }
553  }
554  //2D cases and Homogeneous1D Base Flows
555  else
556  {
557  bool flag = FieldDef[i]->m_fields[j] ==
558  m_session->GetVariable(j);
559 
560  ASSERTL0(flag, (std::string("Order of ") + pInfile
561  + std::string(" data and that defined in "
562  "m_boundaryconditions differs")).c_str());
563 
564  pFields[j]->ExtractDataToCoeffs(FieldDef[i], FieldData[i],
565  FieldDef[i]->m_fields[j],
566  tmp_coeff);
567  }
568  }
569 
570  if(m_SingleMode || m_HalfMode)
571  {
572  pFields[j]->GetPlane(0)->BwdTrans(tmp_coeff, m_baseflow[j]);
573 
574  if(m_SingleMode)
575  {
576  //copy the bwd into the second plane for single Mode Analysis
577  int ncplane=(pFields[0]->GetNpoints())/m_npointsZ;
578  Vmath::Vcopy(ncplane,&m_baseflow[j][0],1,&m_baseflow[j][ncplane],1);
579  }
580  }
581  else
582  {
583  pFields[j]->BwdTrans(tmp_coeff, m_baseflow[j]);
584  }
585  }
586 
587  if(m_session->DefinesParameter("N_slices"))
588  {
589  int n = pFields.num_elements()-1;
590 
591  for(int i=0; i<n;++i)
592  {
593  Vmath::Vcopy(nqtot, &m_baseflow[i][0], 1, &m_interp[i][slice*nqtot], 1);
594  }
595 
596  }
597 }
598 
599 
601  const NekDouble m_slices,
602  const Array<OneD, const NekDouble> &inarray,
603  Array<OneD, NekDouble> &outarray,
604  const NekDouble m_time,
605  const NekDouble m_period)
606 {
607 
608  int npoints=m_baseflow[0].num_elements();
609 
610  NekDouble BetaT=2*M_PI*fmod (m_time, m_period) / m_period;
611  NekDouble phase;
612  Array<OneD, NekDouble> auxiliary(npoints);
613 
614  Vmath::Vcopy(npoints,&inarray[0],1,&outarray[0],1);
615  Vmath::Svtvp(npoints, cos(0.5*m_slices*BetaT),&inarray[npoints],1,&outarray[0],1,&outarray[0],1);
616 
617  for (int i = 2; i < m_slices; i += 2)
618  {
619  phase = (i>>1) * BetaT;
620 
621  Vmath::Svtvp(npoints, cos(phase),&inarray[i*npoints],1,&outarray[0],1,&outarray[0],1);
622  Vmath::Svtvp(npoints, sin(phase), &inarray[(i+1)*npoints], 1, &outarray[0], 1,&outarray[0],1);
623  }
624 
625 }
626 
628 {
629  DNekMatSharedPtr loc_mat;
630  DNekBlkMatSharedPtr BlkMatrix;
631  int n_exp = 0;
632 
633  n_exp = m_baseflow[0].num_elements(); // will operatore on m_phys
634 
635  Array<OneD,unsigned int> nrows(n_exp);
636  Array<OneD,unsigned int> ncols(n_exp);
637 
638  nrows = Array<OneD, unsigned int>(n_exp,m_slices);
639  ncols = Array<OneD, unsigned int>(n_exp,m_slices);
640 
641  MatrixStorage blkmatStorage = eDIAGONAL;
642  BlkMatrix = MemoryManager<DNekBlkMat>
643  ::AllocateSharedPtr(nrows,ncols,blkmatStorage);
644 
645 
648  StdRegions::StdSegExp StdSeg(BK);
649 
651  StdSeg.DetShapeType(),
652  StdSeg);
653 
654  loc_mat = StdSeg.GetStdMatrix(matkey);
655 
656  // set up array of block matrices.
657  for(int i = 0; i < n_exp; ++i)
658  {
659  BlkMatrix->SetBlock(i,i,loc_mat);
660  }
661 
662  return BlkMatrix;
663 }
664 
665 //Discrete Fourier Transform for Floquet analysis
666 void AdjointAdvection::DFT(const string file,
668  const NekDouble m_slices)
669 {
670  int npoints=m_baseflow[0].num_elements();
671 
672  //Convected fields
673  int ConvectedFields=m_baseflow.num_elements()-1;
674 
675  m_interp= Array<OneD, Array<OneD, NekDouble> > (ConvectedFields);
676  for(int i=0; i<ConvectedFields;++i)
677  {
679  }
680 
681  // Import the slides into auxiliary vector
682  // The base flow should be stored in the form "filename_%d.ext"
683  // A subdirectory can also be included, such as "dir/filename_%d.ext"
684  size_t found = file.find("%d");
685  ASSERTL0(found != string::npos && file.find("%d", found+1) == string::npos,
686  "Since N_slices is specified, the filename provided for function "
687  "'BaseFlow' must include exactly one instance of the format "
688  "specifier '%d', to index the time-slices.");
689  char* buffer = new char[file.length() + 8];
690  for (int i = 0; i < m_slices; ++i)
691  {
692  sprintf(buffer, file.c_str(), i);
693  ImportFldBase(buffer,pFields,i);
694  }
695  delete[] buffer;
696 
697 
698  // Discrete Fourier Transform of the fields
699  for(int k=0; k< ConvectedFields;++k)
700  {
701 #ifdef NEKTAR_USING_FFTW
702 
703  //Discrete Fourier Transform using FFTW
704 
705 
706  Array<OneD, NekDouble> fft_in(npoints*m_slices);
707  Array<OneD, NekDouble> fft_out(npoints*m_slices);
708 
711 
712  //Shuffle the data
713  for(int j= 0; j < m_slices; ++j)
714  {
715  Vmath::Vcopy(npoints,&m_interp[k][j*npoints],1,&(fft_in[j]),m_slices);
716  }
717 
718  m_FFT = LibUtilities::GetNektarFFTFactory().CreateInstance("NekFFTW", m_slices);
719 
720  //FFT Transform
721  for(int i=0; i<npoints; i++)
722  {
723  m_FFT->FFTFwdTrans(m_tmpIN =fft_in + i*m_slices, m_tmpOUT =fft_out + i*m_slices);
724 
725  }
726 
727  //Reshuffle data
728  for(int s = 0; s < m_slices; ++s)
729  {
730  Vmath::Vcopy(npoints,&fft_out[s],m_slices,&m_interp[k][s*npoints],1);
731 
732  }
733 
734  Vmath::Zero(fft_in.num_elements(),&fft_in[0],1);
735  Vmath::Zero(fft_out.num_elements(),&fft_out[0],1);
736 #else
737  //Discrete Fourier Transform using MVM
738 
739 
740  DNekBlkMatSharedPtr blkmat;
742 
743  int nrows = blkmat->GetRows();
744  int ncols = blkmat->GetColumns();
745 
746  Array<OneD, NekDouble> sortedinarray(ncols);
747  Array<OneD, NekDouble> sortedoutarray(nrows);
748 
749  //Shuffle the data
750  for(int j= 0; j < m_slices; ++j)
751  {
752  Vmath::Vcopy(npoints,&m_interp[k][j*npoints],1,&(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,&m_interp[k][s*npoints],1);
766  }
767 
768  Vmath::Zero(sortedinarray.num_elements(),&sortedinarray[0],1);
769  Vmath::Zero(sortedoutarray.num_elements(),&sortedoutarray[0],1);
770 
771 #endif
772 
773  //scaling of the Fourier coefficients
774  NekDouble j=-1;
775  for (int i = 2; i < m_slices; i += 2)
776  {
777  Vmath::Smul(2*npoints,j,&m_interp[k][i*npoints],1,&m_interp[k][i*npoints],1);
778  j=-j;
779 
780  }
781 
782  }
783 
784  if(m_session->DefinesParameter("period"))
785  {
786  m_period=m_session->GetParameter("period");
787  }
788  else
789  {
790  m_period=(m_session->GetParameter("TimeStep")*m_slices)/(m_slices-1.);
791  }
792 }
793 
794 
795 } //end of namespace
796 
Array< OneD, NekDouble > m_tmpOUT
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition: StdExpansion.h:454
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
static SolverUtils::AdvectionSharedPtr create(std::string)
Creates an instance of this class.
Local coefficients.
tBaseSharedPtr CreateInstance(tKey idKey BOOST_PP_COMMA_IF(MAX_PARAM) BOOST_PP_ENUM_BINARY_PARAMS(MAX_PARAM, tParam, x))
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:162
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
void ImportFldBase(std::string pInfile, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, int slice)
Import Base flow.
LibUtilities::SessionReaderSharedPtr m_session
int m_npointsX
number of points in X direction (if homogeneous)
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:471
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:428
array buffer
Definition: GsLib.hpp:56
NekDouble m_LhomY
physical length in Y direction (if homogeneous)
Fourier Expansion .
Definition: BasisType.h:52
LibUtilities::NektarFFTSharedPtr m_FFT
boost::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:70
NektarFFTFactory & GetNektarFFTFactory()
Definition: NektarFFT.cpp:69
Array< OneD, Array< OneD, NekDouble > > m_baseflow
Storage for base flow.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:50
virtual 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)
Advects a vector field.
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:684
virtual void v_InitObject(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Initialises the advection object.
Class representing a segment element in reference space.
Definition: StdSegExp.h:54
Array< OneD, Array< OneD, NekDouble > > m_interp
DNekBlkMatSharedPtr GetFloquetBlockMatrix(FloquetMatType mattype, bool UseContCoeffs=false) const
1D Evenly-spaced points using Fourier Fit
Definition: PointsType.h:64
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:199
int m_npointsY
number of points in Y direction (if homogeneous)
static std::string npts
Definition: InputFld.cpp:43
bool m_HalfMode
flag to determine if use half mode or not
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
Definition: Advection.cpp:46
int m_HomoDirec
number of homogenous directions
Defines a specification for a set of points.
Definition: Points.h:58
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:382
boost::shared_ptr< FieldIO > FieldIOSharedPtr
Definition: FieldIO.h:236
double NekDouble
static std::string className
Name of class.
boost::shared_ptr< Equation > EquationSharedPtr
NekDouble m_LhomX
physical length in X direction (if homogeneous)
void UpdateBase(const NekDouble m_slices, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const NekDouble m_time, const NekDouble m_period)
MultiRegions::CoeffState m_CoeffState
boost::shared_ptr< DNekBlkMat > DNekBlkMatSharedPtr
Definition: NekTypeDefs.hpp:72
enum HomogeneousType m_HomogeneousType
void DFT(const string file, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const NekDouble m_slices)
NekDouble m_LhomZ
physical length in Z direction (if homogeneous)
virtual void v_SetBaseFlow(const Array< OneD, Array< OneD, NekDouble > > &inarray)
Overrides the base flow used during linearised advection.
bool m_SingleMode
flag to determine if use single mode or not
int m_npointsZ
number of points in Z direction (if homogeneous)
Array< OneD, NekDouble > m_tmpIN
MultiRegions::ProjectionType m_projectionType
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
bool m_useFFT
flag to determine if use or not the FFT for transformations
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:191
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1038
Describes the specification for a Basis.
Definition: Basis.h:50
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:169
bool m_MultipleModes
flag to determine if use multiple mode or not
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215