Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 using namespace std;
50 
51 namespace Nektar
52 {
53 
54 string AdjointAdvection::className = SolverUtils
56  AdjointAdvection::create);
57 
58 /**
59  *
60  */
61 AdjointAdvection::AdjointAdvection()
62  : Advection()
63 {
64 }
65 
66 
70 {
71 
72  Advection::v_InitObject(pSession, pFields);
73 
74  m_session = pSession;
76  ::AllocateSharedPtr(m_session, pFields[0]->GetGraph());
77  m_spacedim = pFields[0]->GetGraph()->GetSpaceDimension();
78  m_expdim = pFields[0]->GetGraph()->GetMeshDimension();
80 
81  //Setting parameters for homogeneous problems
82  m_HomoDirec = 0;
83  m_useFFT = false;
85  m_SingleMode = false;
86  m_HalfMode = false;
87  m_MultipleModes = false;
88  m_homogen_dealiasing = false;
89 
90  if(m_session->DefinesSolverInfo("HOMOGENEOUS"))
91  {
92  std::string HomoStr = m_session->GetSolverInfo("HOMOGENEOUS");
93  m_spacedim = 3;
94 
95  if((HomoStr == "HOMOGENEOUS1D")||(HomoStr == "Homogeneous1D")||
96  (HomoStr == "1D")||(HomoStr == "Homo1D"))
97  {
99  m_LhomZ = m_session->GetParameter("LZ");
100  m_HomoDirec = 1;
101 
102  m_homogen_dealiasing = pSession->DefinesSolverInfo("DEALIASING");
103 
104  if(m_session->DefinesSolverInfo("ModeType"))
105  {
106  m_session->MatchSolverInfo("ModeType","SingleMode",m_SingleMode,false);
107  m_session->MatchSolverInfo("ModeType","HalfMode",m_HalfMode,false);
108  m_session->MatchSolverInfo("ModeType","MultipleModes",m_MultipleModes,false);
109  }
110  if(m_session->DefinesSolverInfo("ModeType"))
111  {
112  if(m_SingleMode)
113  {
114  m_npointsZ=2;
115  }
116  else if(m_HalfMode)
117  {
118  m_npointsZ=1;
119  }
120  else if(m_MultipleModes)
121  {
122  m_npointsZ = m_session->GetParameter("HomModesZ");
123  }
124  else
125  {
126  ASSERTL0(false, "SolverInfo ModeType not valid");
127  }
128  }
129  else
130  {
131  m_session->LoadParameter("HomModesZ",m_npointsZ);
132 
133  }
134 
135  }
136 
137  if((HomoStr == "HOMOGENEOUS2D")||(HomoStr == "Homogeneous2D")||
138  (HomoStr == "2D")||(HomoStr == "Homo2D"))
139  {
141  m_session->LoadParameter("HomModesY", m_npointsY);
142  m_session->LoadParameter("LY", m_LhomY);
143  m_session->LoadParameter("HomModesZ", m_npointsZ);
144  m_session->LoadParameter("LZ", m_LhomZ);
145  m_HomoDirec = 2;
146  }
147 
148  if((HomoStr == "HOMOGENEOUS3D")||(HomoStr == "Homogeneous3D")||
149  (HomoStr == "3D")||(HomoStr == "Homo3D"))
150  {
152  m_session->LoadParameter("HomModesX",m_npointsX);
153  m_session->LoadParameter("LX", m_LhomX );
154  m_session->LoadParameter("HomModesY",m_npointsY);
155  m_session->LoadParameter("LY", m_LhomY );
156  m_session->LoadParameter("HomModesZ",m_npointsZ);
157  m_session->LoadParameter("LZ", m_LhomZ );
158  m_HomoDirec = 3;
159  }
160 
161  if(m_session->DefinesSolverInfo("USEFFT"))
162  {
163  m_useFFT = true;
164  }
165  }
166  else
167  {
168  m_npointsZ = 1; // set to default value so can use to identify 2d or 3D (homogeneous) expansions
169  }
170 
171  if(m_session->DefinesSolverInfo("PROJECTION"))
172  {
173  std::string ProjectStr
174  = m_session->GetSolverInfo("PROJECTION");
175 
176  if((ProjectStr == "Continuous")||(ProjectStr == "Galerkin")||
177  (ProjectStr == "CONTINUOUS")||(ProjectStr == "GALERKIN"))
178  {
180  }
181  else if(ProjectStr == "DisContinuous")
182  {
184  }
185  else
186  {
187  ASSERTL0(false,"PROJECTION value not recognised");
188  }
189  }
190  else
191  {
192  cerr << "Projection type not specified in SOLVERINFO,"
193  "defaulting to continuous Galerkin" << endl;
195  }
196 
197  int nvar = m_session->GetVariables().size();
199  for (int i = 0; i < nvar; ++i)
200  {
201  m_baseflow[i] = Array<OneD, NekDouble>(pFields[i]->GetTotPoints(), 0.0);
202  }
203  //SetUpBaseFields(pFields[0]->GetGraph());
204  ASSERTL0(m_session->DefinesFunction("BaseFlow"),
205  "Base flow must be defined for linearised forms.");
206  string file = m_session->GetFunctionFilename("BaseFlow", 0);
207 
208 
209  //Periodic base flows
210  if(m_session->DefinesParameter("N_slices"))
211  {
212  m_session->LoadParameter("N_slices",m_slices);
213  if(m_slices>1)
214  {
215  DFT(file,pFields,m_slices);
216  }
217  else
218  {
219  ASSERTL0(false,"Number of slices must be a positive number greater than 1");
220  }
221  }
222  //Steady base-flow
223  else
224  {
225  m_slices=1;
226 
227  //BaseFlow from file
228  if (m_session->GetFunctionType("BaseFlow", m_session->GetVariable(0))
230  {
231  ImportFldBase(file,pFields,1);
232 
233  }
234  //analytic base flow
235  else
236  {
237  int nq = pFields[0]->GetNpoints();
238  Array<OneD,NekDouble> x0(nq);
239  Array<OneD,NekDouble> x1(nq);
240  Array<OneD,NekDouble> x2(nq);
241 
242  // get the coordinates (assuming all fields have the same
243  // discretisation)
244  pFields[0]->GetCoords(x0,x1,x2);
245  for(unsigned int i = 0 ; i < m_baseflow.num_elements(); i++)
246  {
248  = m_session->GetFunction("BaseFlow", i);
249 
250  ifunc->Evaluate(x0,x1,x2,m_baseflow[i]);
251  }
252 
253  }
254 
255  }
256 }
257 
258 
260 {
261 }
262 
263 
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,
270  const NekDouble &time)
271 {
272  int nqtot = fields[0]->GetTotPoints();
273  ASSERTL1(nConvectiveFields == inarray.num_elements(),"Number of convective fields and Inarray are not compatible");
274 
275  Array<OneD, NekDouble > Deriv = Array<OneD, NekDouble> (nqtot*nConvectiveFields);
276 
277  for(int n = 0; n < nConvectiveFields; ++n)
278  {
279  //v_ComputeAdvectionTerm(fields,advVel,inarray[i],outarray[i],i,time,Deriv);
280  int ndim = advVel.num_elements();
281  int nPointsTot = fields[0]->GetNpoints();
282  Array<OneD, NekDouble> grad0,grad1,grad2;
283 
284  //Evaluation of the gradiend of each component of the base flow
285  //\nabla U
286  Array<OneD, NekDouble> grad_base_u0,grad_base_u1,grad_base_u2;
287  // \nabla V
288  Array<OneD, NekDouble> grad_base_v0,grad_base_v1,grad_base_v2;
289  // \nabla W
290  Array<OneD, NekDouble> grad_base_w0,grad_base_w1,grad_base_w2;
291 
292 
293  grad0 = Array<OneD, NekDouble> (nPointsTot);
294  grad_base_u0 = Array<OneD, NekDouble> (nPointsTot);
295  grad_base_v0 = Array<OneD, NekDouble> (nPointsTot);
296  grad_base_w0 = Array<OneD, NekDouble> (nPointsTot);
297 
298 
299  //Evaluation of the base flow for periodic cases
300  //(it requires fld files)
301  if(m_slices>1)
302  {
303  if (m_session->GetFunctionType("BaseFlow", 0)
305  {
306  for(int i=0; i<ndim;++i)
307  {
309  }
310  }
311  else
312  {
313  ASSERTL0(false, "Periodic Base flow requires .fld files");
314  }
315  }
316 
317  //Evaluate the Adjoint advection term
318  switch(ndim)
319  {
320  // 1D
321  case 1:
322  fields[0]->PhysDeriv(inarray[n],grad0);
323  fields[0]->PhysDeriv(m_baseflow[0],grad_base_u0);
324  //Evaluate U du'/dx
325  Vmath::Vmul(nPointsTot,grad0,1,m_baseflow[0],1,outarray[n],1);
326  //Evaluate U du'/dx+ u' dU/dx
327  Vmath::Vvtvp(nPointsTot,grad_base_u0,1,advVel[0],1,outarray[n],1,outarray[n],1);
328  break;
329 
330  //2D
331  case 2:
332 
333  grad1 = Array<OneD, NekDouble> (nPointsTot);
334  grad_base_u1 = Array<OneD, NekDouble> (nPointsTot);
335  grad_base_v1 = Array<OneD, NekDouble> (nPointsTot);
336 
337  fields[0]->PhysDeriv(inarray[n],grad0,grad1);
338 
339  //Derivates of the base flow
340  fields[0]-> PhysDeriv(m_baseflow[0], grad_base_u0, grad_base_u1);
341  fields[0]-> PhysDeriv(m_baseflow[1], grad_base_v0, grad_base_v1);
342 
343  //Since the components of the velocity are passed one by one, it is necessary to distinguish which
344  //term is consider
345  switch (n)
346  {
347  //x-equation
348  case 0:
349  // Evaluate U du'/dx
350  Vmath::Vmul (nPointsTot,grad0,1,m_baseflow[0],1,outarray[n],1);
351  //Evaluate U du'/dx+ V du'/dy
352  Vmath::Vvtvp(nPointsTot,grad1,1,m_baseflow[1],1,outarray[n],1,outarray[n],1);
353  //Evaluate - (U du'/dx+ V du'/dy)
354  Vmath::Neg(nPointsTot,outarray[n],1);
355  //Evaluate -(U du'/dx+ V du'/dy)+u' dU/dx
356  Vmath::Vvtvp(nPointsTot,grad_base_u0,1,advVel[0],1,outarray[n],1,outarray[n],1);
357  //Evaluate -(U du'/dx+ V du'/dy) +u' dU/dx +v' dV/dx
358  Vmath::Vvtvp(nPointsTot,grad_base_v0,1,advVel[1],1,outarray[n],1,outarray[n],1);
359  break;
360 
361  //y-equation
362  case 1:
363  // Evaluate U dv'/dx
364  Vmath::Vmul (nPointsTot,grad0,1,m_baseflow[0],1,outarray[n],1);
365  //Evaluate U dv'/dx+ V dv'/dy
366  Vmath::Vvtvp(nPointsTot,grad1,1,m_baseflow[1],1,outarray[n],1,outarray[n],1);
367  //Evaluate -(U dv'/dx+ V dv'/dy)
368  Vmath::Neg(nPointsTot,outarray[n],1);
369  //Evaluate (U dv'/dx+ V dv'/dy)+u' dU/dy
370  Vmath::Vvtvp(nPointsTot,grad_base_u1,1,advVel[0],1,outarray[n],1,outarray[n],1);
371  //Evaluate (U dv'/dx+ V dv'/dy +u' dv/dx)+v' dV/dy
372  Vmath::Vvtvp(nPointsTot,grad_base_v1,1,advVel[1],1,outarray[n],1,outarray[n],1);
373  break;
374  }
375  break;
376 
377 
378  //3D
379  case 3:
380 
381  grad1 = Array<OneD, NekDouble> (nPointsTot);
382  grad2 = Array<OneD, NekDouble> (nPointsTot);
383  grad_base_u1 = Array<OneD, NekDouble> (nPointsTot);
384  grad_base_v1 = Array<OneD, NekDouble> (nPointsTot);
385  grad_base_w1 = Array<OneD, NekDouble> (nPointsTot);
386 
387  grad_base_u2 = Array<OneD, NekDouble> (nPointsTot);
388  grad_base_v2 = Array<OneD, NekDouble> (nPointsTot);
389  grad_base_w2 = Array<OneD, NekDouble> (nPointsTot);
390 
391  fields[0]->PhysDeriv(m_baseflow[0], grad_base_u0, grad_base_u1,grad_base_u2);
392  fields[0]->PhysDeriv(m_baseflow[1], grad_base_v0, grad_base_v1,grad_base_v2);
393  fields[0]->PhysDeriv(m_baseflow[2], grad_base_w0, grad_base_w1, grad_base_w2);
394 
395  //HalfMode has W(x,y,t)=0
396  if(m_HalfMode || m_SingleMode)
397  {
398  for(int i=0; i<grad_base_u2.num_elements();++i)
399  {
400  grad_base_u2[i]=0;
401  grad_base_v2[i]=0;
402  grad_base_w2[i]=0;
403 
404  }
405  }
406 
407  fields[0]->PhysDeriv(inarray[n], grad0, grad1, grad2);
408 
409  switch (n)
410  {
411  //x-equation
412  case 0:
413  //Evaluate U du'/dx
414  Vmath::Vmul (nPointsTot,grad0,1,m_baseflow[0],1,outarray[n],1);
415  //Evaluate U du'/dx+ V du'/dy
416  Vmath::Vvtvp(nPointsTot,grad1,1,m_baseflow[1],1,outarray[n],1,outarray[n],1);
417  //Evaluate U du'/dx+ V du'/dy+W du'/dz
418  Vmath::Vvtvp(nPointsTot,grad2,1,m_baseflow[2],1,outarray[n],1,outarray[n],1);
419  //Evaluate -(U du'/dx+ V du'/dy+W du'/dz)
420  Vmath::Neg(nPointsTot,outarray[n],1);
421  //Evaluate -(U du'/dx+ V du'/dy+W du'/dz)+u' dU/dx
422  Vmath::Vvtvp(nPointsTot,grad_base_u0,1,advVel[0],1,outarray[n],1,outarray[n],1);
423  //Evaluate -(U du'/dx+ V du'/dy+W du'/dz)+u'dU/dx+ v' dV/dx
424  Vmath::Vvtvp(nPointsTot,grad_base_v0,1,advVel[1],1,outarray[n],1,outarray[n],1);
425  //Evaluate -(U du'/dx+ V du'/dy+W du'/dz)+u'dU/dx+ v' dV/dx+ w' dW/dz
426  Vmath::Vvtvp(nPointsTot,grad_base_w0,1,advVel[2],1,outarray[n],1,outarray[n],1);
427  break;
428  //y-equation
429  case 1:
430  //Evaluate U dv'/dx
431  Vmath::Vmul (nPointsTot,grad0,1,m_baseflow[0],1,outarray[n],1);
432  //Evaluate U dv'/dx+ V dv'/dy
433  Vmath::Vvtvp(nPointsTot,grad1,1,m_baseflow[1],1,outarray[n],1,outarray[n],1);
434  //Evaluate U dv'/dx+ V dv'/dy+W dv'/dz
435  Vmath::Vvtvp(nPointsTot,grad2,1,m_baseflow[2],1,outarray[n],1,outarray[n],1);
436  //Evaluate -(U dv'/dx+ V dv'/dy+W dv'/dz)
437  Vmath::Neg(nPointsTot,outarray[n],1);
438  //Evaluate -(U dv'/dx+ V dv'/dy+W dv'/dz)+u' dU/dy
439  Vmath::Vvtvp(nPointsTot,grad_base_u1,1,advVel[0],1,outarray[n],1,outarray[n],1);
440  //Evaluate -(U dv'/dx+ V dv'/dy+W dv'/dz)+u' dU/dy +v' dV/dy
441  Vmath::Vvtvp(nPointsTot,grad_base_v1,1,advVel[1],1,outarray[n],1,outarray[n],1);
442  //Evaluate -(U dv'/dx+ V dv'/dy+W dv'/dz)+u' dU/dy +v' dV/dy+ w' dW/dy
443  Vmath::Vvtvp(nPointsTot,grad_base_w1,1,advVel[2],1,outarray[n],1,outarray[n],1);
444  break;
445 
446  //z-equation
447  case 2:
448  //Evaluate U dw'/dx
449  Vmath::Vmul (nPointsTot,grad0,1,m_baseflow[0],1,outarray[n],1);
450  //Evaluate U dw'/dx+ V dw'/dx
451  Vmath::Vvtvp(nPointsTot,grad1,1,m_baseflow[1],1,outarray[n],1,outarray[n],1);
452  //Evaluate U dw'/dx+ V dw'/dx+ W dw'/dz
453  Vmath::Vvtvp(nPointsTot,grad2,1,m_baseflow[2],1,outarray[n],1,outarray[n],1);
454  //Evaluate -(U dw'/dx+ V dw'/dx+ W dw'/dz)
455  Vmath::Neg(nPointsTot,outarray[n],1);
456  //Evaluate -(U dw'/dx+ V dw'/dx+ W dw'/dz)+u' dU/dz
457  Vmath::Vvtvp(nPointsTot,grad_base_u2,1,advVel[0],1,outarray[n],1,outarray[n],1);
458  //Evaluate -(U dw'/dx+ V dw'/dx+ W dw'/dz)+u' dU/dz+v'dV/dz
459  Vmath::Vvtvp(nPointsTot,grad_base_v2,1,advVel[1],1,outarray[n],1,outarray[n],1);
460  //Evaluate -(U dw'/dx+ V dw'/dx+ W dw'/dz)+u' dU/dz+v'dV/dz + w' dW/dz
461  Vmath::Vvtvp(nPointsTot,grad_base_w2,1,advVel[2],1,outarray[n],1,outarray[n],1);
462  break;
463  }
464  break;
465  default:
466  ASSERTL0(false,"dimension unknown");
467  }
468 
469  Vmath::Neg(nqtot,outarray[n],1);
470  }
471 
472 }
473 
474 
476  const Array<OneD, Array<OneD, NekDouble> > &inarray)
477 {
478  ASSERTL1(inarray.num_elements() == m_baseflow.num_elements(),
479  "Number of base flow variables does not match what is"
480  "expected.");
481 
482  int npts = inarray[0].num_elements();
483  for (int i = 0; i < inarray.num_elements(); ++i)
484  {
485  ASSERTL1(npts == m_baseflow.num_elements(),
486  "Size of base flow array does not match expected.");
487  Vmath::Vcopy(npts, inarray[i], 1, m_baseflow[i], 1);
488  }
489 }
490 
491 
492 /**
493  * Import field from infile and load into \a m_fields. This routine will
494  * also perform a \a BwdTrans to ensure data is in both the physical and
495  * coefficient storage.
496  * @param infile Filename to read.
497  */
499  std::string pInfile,
501  int pSlice)
502 {
503  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
504  std::vector<std::vector<NekDouble> > FieldData;
505 
506  int nqtot = m_baseflow[0].num_elements();
507  Array<OneD, NekDouble> tmp_coeff(pFields[0]->GetNcoeffs(), 0.0);
508 
509  int numexp = pFields[0]->GetExpSize();
510  Array<OneD,int> ElementGIDs(numexp);
511 
512  // Define list of global element ids
513  for(int i = 0; i < numexp; ++i)
514  {
515  ElementGIDs[i] = pFields[0]->GetExp(i)->GetGeom()->GetGlobalID();
516  }
517 
520  m_session->GetComm());
521  fld->Import(pInfile, FieldDef, FieldData,
523  ElementGIDs);
524 
525  int nSessionVar = m_session->GetVariables().size();
526  int nSessionConvVar = nSessionVar - 1;
527  int nFileVar = FieldDef[0]->m_fields.size();
528  int nFileConvVar = nFileVar - 1; // Ignore pressure
529  if (m_HalfMode)
530  {
531  ASSERTL0(nFileVar == 3, "For half mode, expect 2D2C base flow.");
532  nFileConvVar = 2;
533  }
534 
535  for(int j = 0; j < nFileConvVar; ++j)
536  {
537  for(int i = 0; i < FieldDef.size(); ++i)
538  {
539  bool flag = FieldDef[i]->m_fields[j] ==
540  m_session->GetVariable(j);
541 
542  ASSERTL0(flag, (std::string("Order of ") + pInfile
543  + std::string(" data and that defined in "
544  "the session file differs")).c_str());
545 
546  pFields[j]->ExtractDataToCoeffs(
547  FieldDef[i],
548  FieldData[i],
549  FieldDef[i]->m_fields[j],
550  tmp_coeff);
551  }
552 
553  if(m_SingleMode || m_HalfMode)
554  {
555  pFields[j]->GetPlane(0)->BwdTrans(tmp_coeff, m_baseflow[j]);
556 
557  if(m_SingleMode)
558  {
559  //copy the bwd trans into the second plane for single
560  //Mode Analysis
561  int ncplane=(pFields[0]->GetNpoints())/m_npointsZ;
562  Vmath::Vcopy(ncplane,&m_baseflow[j][0],1,&m_baseflow[j][ncplane],1);
563  }
564  }
565  else // fully 3D base flow - put in physical space.
566  {
567  bool oldwavespace = pFields[j]->GetWaveSpace();
568  pFields[j]->SetWaveSpace(false);
569  pFields[j]->BwdTrans(tmp_coeff, m_baseflow[j]);
570  pFields[j]->SetWaveSpace(oldwavespace);
571  }
572  }
573 
574  // Zero unused fields (e.g. w in a 2D2C base flow).
575  for (int j = nFileConvVar; j < nSessionConvVar; ++j) {
576  Vmath::Fill(nqtot, 0.0, m_baseflow[j], 1);
577  }
578 
579  // If time-periodic, put loaded data into the slice storage.
580  if(m_session->DefinesParameter("N_slices"))
581  {
582  for(int i = 0; i < nSessionConvVar; ++i)
583  {
584  Vmath::Vcopy(nqtot, &m_baseflow[i][0], 1, &m_interp[i][pSlice*nqtot], 1);
585  }
586  }
587 }
588 
589 
591  const NekDouble m_slices,
592  const Array<OneD, const NekDouble> &inarray,
593  Array<OneD, NekDouble> &outarray,
594  const NekDouble m_time,
595  const NekDouble m_period)
596 {
597 
598  int npoints=m_baseflow[0].num_elements();
599 
600  NekDouble BetaT=2*M_PI*fmod (m_time, m_period) / m_period;
601  NekDouble phase;
602  Array<OneD, NekDouble> auxiliary(npoints);
603 
604  Vmath::Vcopy(npoints,&inarray[0],1,&outarray[0],1);
605  Vmath::Svtvp(npoints, cos(0.5*m_slices*BetaT),&inarray[npoints],1,&outarray[0],1,&outarray[0],1);
606 
607  for (int i = 2; i < m_slices; i += 2)
608  {
609  phase = (i>>1) * BetaT;
610 
611  Vmath::Svtvp(npoints, cos(phase),&inarray[i*npoints],1,&outarray[0],1,&outarray[0],1);
612  Vmath::Svtvp(npoints, sin(phase), &inarray[(i+1)*npoints], 1, &outarray[0], 1,&outarray[0],1);
613  }
614 
615 }
616 
618 {
619  DNekMatSharedPtr loc_mat;
620  DNekBlkMatSharedPtr BlkMatrix;
621  int n_exp = 0;
622 
623  n_exp = m_baseflow[0].num_elements(); // will operatore on m_phys
624 
625  Array<OneD,unsigned int> nrows(n_exp);
626  Array<OneD,unsigned int> ncols(n_exp);
627 
628  nrows = Array<OneD, unsigned int>(n_exp,m_slices);
629  ncols = Array<OneD, unsigned int>(n_exp,m_slices);
630 
631  MatrixStorage blkmatStorage = eDIAGONAL;
632  BlkMatrix = MemoryManager<DNekBlkMat>
633  ::AllocateSharedPtr(nrows,ncols,blkmatStorage);
634 
635 
638  StdRegions::StdSegExp StdSeg(BK);
639 
641  StdSeg.DetShapeType(),
642  StdSeg);
643 
644  loc_mat = StdSeg.GetStdMatrix(matkey);
645 
646  // set up array of block matrices.
647  for(int i = 0; i < n_exp; ++i)
648  {
649  BlkMatrix->SetBlock(i,i,loc_mat);
650  }
651 
652  return BlkMatrix;
653 }
654 
655 //Discrete Fourier Transform for Floquet analysis
656 void AdjointAdvection::DFT(const string file,
658  const NekDouble m_slices)
659 {
660  int npoints=m_baseflow[0].num_elements();
661 
662  //Convected fields
663  int ConvectedFields=m_baseflow.num_elements()-1;
664 
665  m_interp= Array<OneD, Array<OneD, NekDouble> > (ConvectedFields);
666  for(int i=0; i<ConvectedFields;++i)
667  {
669  }
670 
671  // Import the slides into auxiliary vector
672  // The base flow should be stored in the form "filename_%d.ext"
673  // A subdirectory can also be included, such as "dir/filename_%d.ext"
674  size_t found = file.find("%d");
675  ASSERTL0(found != string::npos && file.find("%d", found+1) == string::npos,
676  "Since N_slices is specified, the filename provided for function "
677  "'BaseFlow' must include exactly one instance of the format "
678  "specifier '%d', to index the time-slices.");
679  char* buffer = new char[file.length() + 8];
680  for (int i = 0; i < m_slices; ++i)
681  {
682  sprintf(buffer, file.c_str(), i);
683  ImportFldBase(buffer,pFields,i);
684  }
685  delete[] buffer;
686 
687 
688  // Discrete Fourier Transform of the fields
689  for(int k=0; k< ConvectedFields;++k)
690  {
691 #ifdef NEKTAR_USING_FFTW
692 
693  //Discrete Fourier Transform using FFTW
694 
695 
696  Array<OneD, NekDouble> fft_in(npoints*m_slices);
697  Array<OneD, NekDouble> fft_out(npoints*m_slices);
698 
701 
702  //Shuffle the data
703  for(int j= 0; j < m_slices; ++j)
704  {
705  Vmath::Vcopy(npoints,&m_interp[k][j*npoints],1,&(fft_in[j]),m_slices);
706  }
707 
708  m_FFT = LibUtilities::GetNektarFFTFactory().CreateInstance("NekFFTW", m_slices);
709 
710  //FFT Transform
711  for(int i=0; i<npoints; i++)
712  {
713  m_FFT->FFTFwdTrans(m_tmpIN =fft_in + i*m_slices, m_tmpOUT =fft_out + i*m_slices);
714 
715  }
716 
717  //Reshuffle data
718  for(int s = 0; s < m_slices; ++s)
719  {
720  Vmath::Vcopy(npoints,&fft_out[s],m_slices,&m_interp[k][s*npoints],1);
721 
722  }
723 
724  Vmath::Zero(fft_in.num_elements(),&fft_in[0],1);
725  Vmath::Zero(fft_out.num_elements(),&fft_out[0],1);
726 #else
727  //Discrete Fourier Transform using MVM
728 
729 
730  DNekBlkMatSharedPtr blkmat;
732 
733  int nrows = blkmat->GetRows();
734  int ncols = blkmat->GetColumns();
735 
736  Array<OneD, NekDouble> sortedinarray(ncols);
737  Array<OneD, NekDouble> sortedoutarray(nrows);
738 
739  //Shuffle the data
740  for(int j= 0; j < m_slices; ++j)
741  {
742  Vmath::Vcopy(npoints,&m_interp[k][j*npoints],1,&(sortedinarray[j]),m_slices);
743  }
744 
745  // Create NekVectors from the given data arrays
746  NekVector<NekDouble> in (ncols,sortedinarray,eWrapper);
747  NekVector<NekDouble> out(nrows,sortedoutarray,eWrapper);
748 
749  // Perform matrix-vector multiply.
750  out = (*blkmat)*in;
751 
752  //Reshuffle data
753  for(int s = 0; s < m_slices; ++s)
754  {
755  Vmath::Vcopy(npoints,&sortedoutarray[s],m_slices,&m_interp[k][s*npoints],1);
756  }
757 
758  Vmath::Zero(sortedinarray.num_elements(),&sortedinarray[0],1);
759  Vmath::Zero(sortedoutarray.num_elements(),&sortedoutarray[0],1);
760 
761 #endif
762 
763  //scaling of the Fourier coefficients
764  NekDouble j=-1;
765  for (int i = 2; i < m_slices; i += 2)
766  {
767  Vmath::Smul(2*npoints,j,&m_interp[k][i*npoints],1,&m_interp[k][i*npoints],1);
768  j=-j;
769 
770  }
771 
772  }
773 
774  if(m_session->DefinesParameter("period"))
775  {
776  m_period=m_session->GetParameter("period");
777  }
778  else
779  {
780  m_period=(m_session->GetParameter("TimeStep")*m_slices)/(m_slices-1.);
781  }
782 }
783 
784 
785 } //end of namespace
786 
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:470
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
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 Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:46
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:57
STL namespace.
NekDouble m_LhomY
physical length in Y direction (if homogeneous)
Fourier Expansion .
Definition: BasisType.h:52
LibUtilities::NektarFFTSharedPtr m_FFT
void DFT(const std::string file, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const NekDouble m_slices)
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:51
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:700
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:225
double NekDouble
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
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:1047
static FieldMetaDataMap NullFieldMetaDataMap
Definition: FieldIO.h:54
Describes the specification for a Basis.
Definition: Basis.h:50
virtual SOLVER_UTILS_EXPORT void v_InitObject(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Initialises the advection object.
Definition: Advection.cpp:97
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
Defines a callback function which evaluates the flux vector.
Definition: Advection.h:69
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