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 
36 #include <IncNavierStokesSolver/AdvectionTerms/AdjointAdvection.h>
37 #include <StdRegions/StdSegExp.h>
38 
39 #include <MultiRegions/ContField1D.h>
40 #include <MultiRegions/ContField2D.h>
41 #include <MultiRegions/ContField3D.h>
42 #include <MultiRegions/DisContField1D.h>
43 #include <MultiRegions/DisContField2D.h>
44 #include <MultiRegions/ContField3DHomogeneous1D.h>
45 #include <MultiRegions/DisContField3DHomogeneous1D.h>
46 #include <MultiRegions/ContField3DHomogeneous2D.h>
47 #include <MultiRegions/DisContField3DHomogeneous2D.h>
48 
49 
50 namespace Nektar
51 {
52  string AdjointAdvection::className = GetAdvectionTermFactory().RegisterCreatorFunction("Adjoint", AdjointAdvection::create);
53 
54  /**
55  * Constructor. Creates ...
56  *
57  * \param
58  * \param
59  */
60 
61  AdjointAdvection::AdjointAdvection(
62  const LibUtilities::SessionReaderSharedPtr& pSession,
63  const SpatialDomains::MeshGraphSharedPtr& pGraph):
64  AdvectionTerm(pSession, pGraph)
65  {
66  }
67 
68 
69  void AdjointAdvection::v_InitObject()
70  {
71  AdvectionTerm::v_InitObject();
72  m_boundaryConditions = MemoryManager<SpatialDomains::BoundaryConditions>
73  ::AllocateSharedPtr(m_session, m_graph);
74 
75  //Setting parameters for homogeneous problems
76  m_HomoDirec = 0;
77  m_useFFT = false;
78  m_HomogeneousType = eNotHomogeneous;
79  m_SingleMode =false;
80  m_HalfMode =false;
81  m_MultipleModes =false;
82 
83  if(m_session->DefinesSolverInfo("HOMOGENEOUS"))
84  {
85  std::string HomoStr = m_session->GetSolverInfo("HOMOGENEOUS");
86  m_spacedim = 3;
87 
88  if((HomoStr == "HOMOGENEOUS1D")||(HomoStr == "Homogeneous1D")||
89  (HomoStr == "1D")||(HomoStr == "Homo1D"))
90  {
91  m_HomogeneousType = eHomogeneous1D;
92  m_LhomZ = m_session->GetParameter("LZ");
93  m_HomoDirec = 1;
94 
95  if(m_session->DefinesSolverInfo("ModeType"))
96  {
97  m_session->MatchSolverInfo("ModeType","SingleMode",m_SingleMode,false);
98  m_session->MatchSolverInfo("ModeType","HalfMode",m_HalfMode,false);
99  m_session->MatchSolverInfo("ModeType","MultipleModes",m_MultipleModes,false);
100  }
101  if(m_session->DefinesSolverInfo("ModeType"))
102  {
103  if(m_SingleMode)
104  {
105  m_npointsZ=2;
106 
107  }
108  else if(m_HalfMode)
109  {
110  m_npointsZ=1;
111 
112  }
113  else if(m_MultipleModes)
114  {
115  m_npointsZ = m_session->GetParameter("HomModesZ");
116  }
117  else
118  {
119  ASSERTL0(false, "SolverInfo ModeType not valid");
120 
121 
122  }
123  }
124  else
125  {
126  m_session->LoadParameter("HomModesZ",m_npointsZ);
127 
128  }
129 
130  }
131 
132  if((HomoStr == "HOMOGENEOUS2D")||(HomoStr == "Homogeneous2D")||
133  (HomoStr == "2D")||(HomoStr == "Homo2D"))
134  {
135  m_HomogeneousType = eHomogeneous2D;
136  m_session->LoadParameter("HomModesY", m_npointsY);
137  m_session->LoadParameter("LY", m_LhomY);
138  m_session->LoadParameter("HomModesZ", m_npointsZ);
139  m_session->LoadParameter("LZ", m_LhomZ);
140  m_HomoDirec = 2;
141  }
142 
143  if((HomoStr == "HOMOGENEOUS3D")||(HomoStr == "Homogeneous3D")||
144  (HomoStr == "3D")||(HomoStr == "Homo3D"))
145  {
146  m_HomogeneousType = eHomogeneous3D;
147  m_session->LoadParameter("HomModesX",m_npointsX);
148  m_session->LoadParameter("LX", m_LhomX );
149  m_session->LoadParameter("HomModesY",m_npointsY);
150  m_session->LoadParameter("LY", m_LhomY );
151  m_session->LoadParameter("HomModesZ",m_npointsZ);
152  m_session->LoadParameter("LZ", m_LhomZ );
153  m_HomoDirec = 3;
154  }
155 
156  if(m_session->DefinesSolverInfo("USEFFT"))
157  {
158  m_useFFT = true;
159  }
160  }
161  else
162  {
163  m_npointsZ = 1; // set to default value so can use to identify 2d or 3D (homogeneous) expansions
164  }
165 
166  if(m_session->DefinesSolverInfo("PROJECTION"))
167  {
168  std::string ProjectStr
169  = m_session->GetSolverInfo("PROJECTION");
170 
171  if((ProjectStr == "Continuous")||(ProjectStr == "Galerkin")||
172  (ProjectStr == "CONTINUOUS")||(ProjectStr == "GALERKIN"))
173  {
174  m_projectionType = MultiRegions::eGalerkin;
175  }
176  else if(ProjectStr == "DisContinuous")
177  {
178  m_projectionType = MultiRegions::eDiscontinuous;
179  }
180  else
181  {
182  ASSERTL0(false,"PROJECTION value not recognised");
183  }
184  }
185  else
186  {
187  cerr << "Projection type not specified in SOLVERINFO,"
188  "defaulting to continuous Galerkin" << endl;
189  m_projectionType = MultiRegions::eGalerkin;
190  }
191 
192  SetUpBaseFields(m_graph);
193  ASSERTL0(m_session->DefinesFunction("BaseFlow"),
194  "Base flow must be defined for linearised forms.");
195  string file = m_session->GetFunctionFilename("BaseFlow", 0);
196 
197 
198  //Periodic base flows
199  if(m_session->DefinesParameter("N_slices"))
200  {
201  m_session->LoadParameter("N_slices",m_slices);
202  if(m_slices>1)
203  {
204  DFT(file,m_slices);
205  }
206  else
207  {
208  ASSERTL0(false,"Number of slices must be a positive number greater than 1");
209  }
210  }
211  //Steady base-flow
212  else
213  {
214  m_slices=1;
215 
216  //BaseFlow from file
217  if (m_session->GetFunctionType("BaseFlow", m_session->GetVariable(0))
218  == LibUtilities::eFunctionTypeFile)
219  {
220  ImportFldBase(file,m_graph,1);
221 
222  }
223  //analytic base flow
224  else
225  {
226  int nq = m_base[0]->GetNpoints();
227  Array<OneD,NekDouble> x0(nq);
228  Array<OneD,NekDouble> x1(nq);
229  Array<OneD,NekDouble> x2(nq);
230 
231  // get the coordinates (assuming all fields have the same
232  // discretisation)
233  m_base[0]->GetCoords(x0,x1,x2);
234  for(unsigned int i = 0 ; i < m_base.num_elements(); i++)
235  {
236  LibUtilities::EquationSharedPtr ifunc
237  = m_session->GetFunction("BaseFlow", i);
238 
239  ifunc->Evaluate(x0,x1,x2,m_base[i]->UpdatePhys());
240 
241  m_base[i]->SetPhysState(true);
242  m_base[i]->FwdTrans_IterPerExp(m_base[i]->GetPhys(),
243  m_base[i]->UpdateCoeffs());
244  }
245 
246  }
247 
248  }
249 }
250 
251 
252  AdjointAdvection::~AdjointAdvection()
253  {
254  }
255 
256 
257  void AdjointAdvection::SetUpBaseFields(SpatialDomains::MeshGraphSharedPtr &mesh)
258  {
259  int nvariables = m_session->GetVariables().size();
260  int i;
261  m_base = Array<OneD, MultiRegions::ExpListSharedPtr>(nvariables);
262  if (m_projectionType == MultiRegions::eGalerkin)
263  {
264  switch (m_expdim)
265  {
266  case 1:
267  {
268  if(m_HomogeneousType == eHomogeneous2D)
269  {
270  const LibUtilities::PointsKey PkeyY(m_npointsY,LibUtilities::eFourierEvenlySpaced);
271  const LibUtilities::BasisKey BkeyY(LibUtilities::eFourier,m_npointsY,PkeyY);
272  const LibUtilities::PointsKey PkeyZ(m_npointsZ,LibUtilities::eFourierEvenlySpaced);
273  const LibUtilities::BasisKey BkeyZ(LibUtilities::eFourier,m_npointsZ,PkeyZ);
274 
275  for(i = 0 ; i < m_base.num_elements(); i++)
276  {
277  m_base[i] = MemoryManager<MultiRegions::ContField3DHomogeneous2D>
278  ::AllocateSharedPtr(m_session,BkeyY,BkeyZ,m_LhomY,m_LhomZ,m_useFFT,m_homogen_dealiasing,m_graph,m_session->GetVariable(i));
279  }
280  }
281 
282  else {
283 
284  for(i = 0 ; i < m_base.num_elements(); i++)
285  {
286  m_base[i] = MemoryManager<MultiRegions::ContField1D>
287  ::AllocateSharedPtr(m_session,mesh,
288  m_session->GetVariable(i));
289  }
290  }
291 
292  }
293  break;
294  case 2:
295  {
296  if(m_HomogeneousType == eHomogeneous1D)
297  {
298  if(m_SingleMode)
299  {
300  const LibUtilities::PointsKey PkeyZ(m_npointsZ,LibUtilities::eFourierSingleModeSpaced);
301  const LibUtilities::BasisKey BkeyZ(LibUtilities::eFourier,m_npointsZ,PkeyZ);
302 
303  for(i = 0 ; i < m_base.num_elements(); i++)
304  {
305  m_base[i] = MemoryManager<MultiRegions::ContField3DHomogeneous1D>
306  ::AllocateSharedPtr(m_session,BkeyZ,m_LhomZ,m_useFFT,m_homogen_dealiasing,m_graph,m_session->GetVariable(i));
307 
308  }
309  }
310  else if(m_HalfMode)
311  {
312  //1 plane field (half mode expansion)
313  const LibUtilities::PointsKey PkeyZ(m_npointsZ,LibUtilities::eFourierSingleModeSpaced);
314  const LibUtilities::BasisKey BkeyZ(LibUtilities::eFourierHalfModeRe,m_npointsZ,PkeyZ);
315 
316  for(i = 0 ; i < m_base.num_elements(); i++)
317  {
318  m_base[i] = MemoryManager<MultiRegions::ContField3DHomogeneous1D>
319  ::AllocateSharedPtr(m_session,BkeyZ,m_LhomZ,m_useFFT,m_homogen_dealiasing,m_graph,m_session->GetVariable(i));
320  }
321 
322  }
323  else
324  {
325  const LibUtilities::PointsKey PkeyZ(m_npointsZ,LibUtilities::eFourierEvenlySpaced);
326  const LibUtilities::BasisKey BkeyZ(LibUtilities::eFourier,m_npointsZ,PkeyZ);
327 
328 
329  for(i = 0 ; i < m_base.num_elements(); i++)
330  {
331  m_base[i] = MemoryManager<MultiRegions::ContField3DHomogeneous1D>
332  ::AllocateSharedPtr(m_session,BkeyZ,m_LhomZ,m_useFFT,m_homogen_dealiasing,m_graph,m_session->GetVariable(i));
333  m_base[i]->SetWaveSpace(false);
334  }
335 
336  }
337  }
338  else
339  {
340  i = 0;
341  MultiRegions::ContField2DSharedPtr firstbase =
342  MemoryManager<MultiRegions::ContField2D>
343  ::AllocateSharedPtr(m_session,mesh,
344  m_session->GetVariable(i));
345  m_base[0]=firstbase;
346 
347  for(i = 1 ; i < m_base.num_elements(); i++)
348  {
349  m_base[i] = MemoryManager<MultiRegions::ContField2D>
350  ::AllocateSharedPtr(*firstbase,mesh,
351  m_session->GetVariable(i));
352  }
353  }
354  }
355  break;
356  case 3:
357  {
358  if(m_HomogeneousType == eHomogeneous3D)
359  {
360  ASSERTL0(false,"3D fully periodic problems not implemented yet");
361  }
362  else
363  {
364  MultiRegions::ContField3DSharedPtr firstbase =
365  MemoryManager<MultiRegions::ContField3D>
366  ::AllocateSharedPtr(m_session,mesh,
367  m_session->GetVariable(0));
368  m_base[0] = firstbase;
369 
370  for(i = 1 ; i < m_base.num_elements(); i++)
371  {
372  m_base[i] = MemoryManager<MultiRegions::ContField3D>
373  ::AllocateSharedPtr(*firstbase,mesh,
374  m_session->GetVariable(i));
375  }
376  }
377  }
378  break;
379  default:
380  ASSERTL0(false,"Expansion dimension not recognised");
381  break;
382  }
383  }
384  else
385  {
386  switch(m_expdim)
387  {
388  case 1:
389  {
390  if(m_HomogeneousType == eHomogeneous2D)
391  {
392  const LibUtilities::PointsKey PkeyY(m_npointsY,LibUtilities::eFourierEvenlySpaced);
393  const LibUtilities::BasisKey BkeyY(LibUtilities::eFourier,m_npointsY,PkeyY);
394  const LibUtilities::PointsKey PkeyZ(m_npointsZ,LibUtilities::eFourierEvenlySpaced);
395  const LibUtilities::BasisKey BkeyZ(LibUtilities::eFourier,m_npointsZ,PkeyZ);
396 
397  for(i = 0 ; i < m_base.num_elements(); i++)
398  {
399  m_base[i] = MemoryManager<MultiRegions::DisContField3DHomogeneous2D>
400  ::AllocateSharedPtr(m_session,BkeyY,BkeyZ,m_LhomY,m_LhomZ,m_useFFT,m_homogen_dealiasing,m_graph,m_session->GetVariable(i));
401  }
402  }
403  else
404  {
405  for(i = 0 ; i < m_base.num_elements(); i++)
406  {
407  m_base[i] = MemoryManager<MultiRegions
408  ::DisContField1D>::AllocateSharedPtr(m_session,mesh,
409  m_session->GetVariable(i));
410  }
411  }
412  break;
413  }
414  case 2:
415  {
416  if(m_HomogeneousType == eHomogeneous1D)
417  {
418 
419  const LibUtilities::PointsKey PkeyZ(m_npointsZ,LibUtilities::eFourierEvenlySpaced);
420  const LibUtilities::BasisKey BkeyZ(LibUtilities::eFourier,m_npointsZ,PkeyZ);
421 
422  for(i = 0 ; i < m_base.num_elements(); i++)
423  {
424  m_base[i] = MemoryManager<MultiRegions::DisContField3DHomogeneous1D>
425  ::AllocateSharedPtr(m_session,BkeyZ,m_LhomZ,m_useFFT,m_homogen_dealiasing,m_graph,m_session->GetVariable(i));
426  }
427 
428 
429 
430  }
431  else
432  {
433  for(i = 0 ; i < m_base.num_elements(); i++)
434  {
435  m_base[i] = MemoryManager<MultiRegions
436  ::DisContField2D>::AllocateSharedPtr(m_session, mesh,
437  m_session->GetVariable(i));
438  }
439  }
440  break;
441 
442  }
443  case 3:
444  ASSERTL0(false,"3 D not set up");
445  default:
446  ASSERTL0(false,"Expansion dimension not recognised");
447  break;
448  }
449  }
450 
451  }
452 
453  /**
454  * Import field from infile and load into \a m_fields. This routine will
455  * also perform a \a BwdTrans to ensure data is in both the physical and
456  * coefficient storage.
457  * @param infile Filename to read.
458  */
459  void AdjointAdvection::ImportFldBase(std::string pInfile,
460  SpatialDomains::MeshGraphSharedPtr pGraph, int cnt)
461  {
462  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
463  std::vector<std::vector<NekDouble> > FieldData;
464  int nqtot = m_base[0]->GetTotPoints();
465 
466  //Get Homogeneous
467  LibUtilities::FieldIOSharedPtr fld =
468  MemoryManager<LibUtilities::FieldIO>::AllocateSharedPtr(
469  m_session->GetComm());
470  fld->Import(pInfile, FieldDef, FieldData);
471 
472  int nvar = m_session->GetVariables().size();
473  int s;
474 
475  if(m_session->DefinesSolverInfo("HOMOGENEOUS"))
476  {
477  std::string HomoStr = m_session->GetSolverInfo("HOMOGENEOUS");
478  }
479 
480  // copy FieldData into m_fields
481  for(int j = 0; j < nvar; ++j)
482  {
483  for(int i = 0; i < FieldDef.size(); ++i)
484  {
485  if ((m_session->DefinesSolverInfo("HOMOGENEOUS") &&
486  (m_session->GetSolverInfo("HOMOGENEOUS")=="HOMOGENEOUS1D" ||
487  m_session->GetSolverInfo("HOMOGENEOUS")=="1D" ||
488  m_session->GetSolverInfo("HOMOGENEOUS")=="Homo1D")) &&
489  nvar==4)
490  {
491  // w-component must be ignored and set to zero.
492  if (j != nvar - 2)
493  {
494  // p component (it is 4th variable of the 3D and corresponds 3nd variable of 2D)
495  s = (j == nvar - 1) ? 2 : j;
496 
497  //extraction of the 2D
498  m_base[j]->ExtractDataToCoeffs(
499  FieldDef[i],
500  FieldData[i],
501  FieldDef[i]->m_fields[s],
502  m_base[j]->UpdateCoeffs());
503  }
504 
505  // Put zero on higher modes
506  int ncplane = (m_base[0]->GetNcoeffs()) / m_npointsZ;
507 
508  if (m_npointsZ > 2)
509  {
510  Vmath::Zero(ncplane*(m_npointsZ-2),
511  &m_base[j]->UpdateCoeffs()[2*ncplane],1);
512  }
513  }
514  //2D cases and Homogeneous1D Base Flows
515  else
516  {
517  bool flag = FieldDef[i]->m_fields[j] ==
518  m_session->GetVariable(j);
519 
520  ASSERTL1(flag, (std::string("Order of ") + pInfile
521  + std::string(" data and that defined in "
522  "m_boundaryconditions differs")).c_str());
523 
524  m_base[j]->ExtractDataToCoeffs(FieldDef[i], FieldData[i],
525  FieldDef[i]->m_fields[j],
526  m_base[j]->UpdateCoeffs());
527  }
528  }
529 
530  if(m_SingleMode || m_HalfMode)
531  {
532  m_base[j]->SetWaveSpace(true);
533 
534  m_base[j]->BwdTrans(m_base[j]->GetCoeffs(),
535  m_base[j]->UpdatePhys());
536 
537  if(m_SingleMode)
538  {
539  //copy the bwd into the second plane for single Mode Analysis
540  int ncplane=(m_base[0]->GetNpoints())/m_npointsZ;
541  Vmath::Vcopy(ncplane,&m_base[j]->GetPhys()[0],1,&m_base[j]->UpdatePhys()[ncplane],1);
542  }
543  }
544  else
545  {
546  m_base[j]->BwdTrans(m_base[j]->GetCoeffs(),
547  m_base[j]->UpdatePhys());
548 
549  }
550 
551  }
552 
553  //std::string outname ="BaseFlow.bse";
554  //WriteFldBase(outname);
555 
556  if(m_session->DefinesParameter("N_slices"))
557  {
558  m_nConvectiveFields = m_base.num_elements()-1;
559 
560  for(int i=0; i<m_nConvectiveFields;++i)
561  {
562 
563  Vmath::Vcopy(nqtot, &m_base[i]->GetPhys()[0], 1, &m_interp[i][cnt*nqtot], 1);
564  }
565 
566  }
567  }
568 
569 
570  //Evaluation of the advective terms
571  void AdjointAdvection::v_ComputeAdvectionTerm(
572  Array<OneD, MultiRegions::ExpListSharedPtr > &pFields,
573  const Array<OneD, Array<OneD, NekDouble> > &pVelocity,
574  const Array<OneD, const NekDouble> &pU,
575  Array<OneD, NekDouble> &pOutarray,
576  int pVelocityComponent,
577  NekDouble m_time,
578  Array<OneD, NekDouble> &pWk)
579  {
580  int ndim = m_nConvectiveFields;
581  int nPointsTot = pFields[0]->GetNpoints();
582  Array<OneD, NekDouble> grad0,grad1,grad2;
583 
584  //Evaluation of the gradiend of each component of the base flow
585  //\nabla U
586  Array<OneD, NekDouble> grad_base_u0,grad_base_u1,grad_base_u2;
587  // \nabla V
588  Array<OneD, NekDouble> grad_base_v0,grad_base_v1,grad_base_v2;
589  // \nabla W
590  Array<OneD, NekDouble> grad_base_w0,grad_base_w1,grad_base_w2;
591 
592 
593  grad0 = Array<OneD, NekDouble> (nPointsTot);
594  grad_base_u0 = Array<OneD, NekDouble> (nPointsTot);
595  grad_base_v0 = Array<OneD, NekDouble> (nPointsTot);
596  grad_base_w0 = Array<OneD, NekDouble> (nPointsTot);
597 
598 
599  //Evaluation of the base flow for periodic cases
600  //(it requires fld files)
601 
602  if(m_slices>1)
603  {
604  if (m_session->GetFunctionType("BaseFlow", 0)
605  == LibUtilities::eFunctionTypeFile)
606  {
607  for(int i=0; i<m_nConvectiveFields;++i)
608  {
609  UpdateBase(m_slices,m_interp[i],m_base[i]->UpdatePhys(),m_time,m_period);
610  }
611  }
612  else
613  {
614  ASSERTL0(false, "Periodic Base flow requires .fld files");
615  }
616  }
617 
618  //Evaluate the Adjoint advection term
619  switch(ndim)
620  {
621  // 1D
622  case 1:
623  pFields[0]->PhysDeriv(pU,grad0);
624  pFields[0]->PhysDeriv(m_base[0]->GetPhys(),grad_base_u0);
625  //Evaluate U du'/dx
626  Vmath::Vmul(nPointsTot,grad0,1,m_base[0]->GetPhys(),1,pOutarray,1);
627  //Evaluate U du'/dx+ u' dU/dx
628  Vmath::Vvtvp(nPointsTot,grad_base_u0,1,pVelocity[0],1,pOutarray,1,pOutarray,1);
629  break;
630 
631  //2D
632  case 2:
633 
634  grad1 = Array<OneD, NekDouble> (nPointsTot);
635  grad_base_u1 = Array<OneD, NekDouble> (nPointsTot);
636  grad_base_v1 = Array<OneD, NekDouble> (nPointsTot);
637 
638  pFields[0]->PhysDeriv(pU,grad0,grad1);
639 
640  //Derivates of the base flow
641  pFields[0]-> PhysDeriv(m_base[0]->GetPhys(), grad_base_u0, grad_base_u1);
642  pFields[0]-> PhysDeriv(m_base[1]->GetPhys(), grad_base_v0, grad_base_v1);
643 
644  //Since the components of the velocity are passed one by one, it is necessary to distinguish which
645  //term is consider
646  switch (pVelocityComponent)
647  {
648  //x-equation
649  case 0:
650  // Evaluate U du'/dx
651  Vmath::Vmul (nPointsTot,grad0,1,m_base[0]->GetPhys(),1,pOutarray,1);
652  //Evaluate U du'/dx+ V du'/dy
653  Vmath::Vvtvp(nPointsTot,grad1,1,m_base[1]->GetPhys(),1,pOutarray,1,pOutarray,1);
654  //Evaluate - (U du'/dx+ V du'/dy)
655  Vmath::Neg(nPointsTot,pOutarray,1);
656  //Evaluate -(U du'/dx+ V du'/dy)+u' dU/dx
657  Vmath::Vvtvp(nPointsTot,grad_base_u0,1,pVelocity[0],1,pOutarray,1,pOutarray,1);
658  //Evaluate -(U du'/dx+ V du'/dy) +u' dU/dx +v' dV/dx
659  Vmath::Vvtvp(nPointsTot,grad_base_v0,1,pVelocity[1],1,pOutarray,1,pOutarray,1);
660  break;
661 
662  //y-equation
663  case 1:
664  // Evaluate U dv'/dx
665  Vmath::Vmul (nPointsTot,grad0,1,m_base[0]->GetPhys(),1,pOutarray,1);
666  //Evaluate U dv'/dx+ V dv'/dy
667  Vmath::Vvtvp(nPointsTot,grad1,1,m_base[1]->GetPhys(),1,pOutarray,1,pOutarray,1);
668  //Evaluate -(U dv'/dx+ V dv'/dy)
669  Vmath::Neg(nPointsTot,pOutarray,1);
670  //Evaluate (U dv'/dx+ V dv'/dy)+u' dU/dy
671  Vmath::Vvtvp(nPointsTot,grad_base_u1,1,pVelocity[0],1,pOutarray,1,pOutarray,1);
672  //Evaluate (U dv'/dx+ V dv'/dy +u' dv/dx)+v' dV/dy
673  Vmath::Vvtvp(nPointsTot,grad_base_v1,1,pVelocity[1],1,pOutarray,1,pOutarray,1);
674  break;
675  }
676  break;
677 
678 
679  //3D
680  case 3:
681 
682  grad1 = Array<OneD, NekDouble> (nPointsTot);
683  grad2 = Array<OneD, NekDouble> (nPointsTot);
684  grad_base_u1 = Array<OneD, NekDouble> (nPointsTot);
685  grad_base_v1 = Array<OneD, NekDouble> (nPointsTot);
686  grad_base_w1 = Array<OneD, NekDouble> (nPointsTot);
687 
688  grad_base_u2 = Array<OneD, NekDouble> (nPointsTot);
689  grad_base_v2 = Array<OneD, NekDouble> (nPointsTot);
690  grad_base_w2 = Array<OneD, NekDouble> (nPointsTot);
691 
692  m_base[0]->PhysDeriv(m_base[0]->GetPhys(), grad_base_u0, grad_base_u1,grad_base_u2);
693  m_base[0]->PhysDeriv(m_base[1]->GetPhys(), grad_base_v0, grad_base_v1,grad_base_v2);
694  m_base[0]->PhysDeriv(m_base[2]->GetPhys(), grad_base_w0, grad_base_w1, grad_base_w2);
695 
696  //HalfMode has W(x,y,t)=0
697  if(m_HalfMode)
698  {
699  for(int i=0; i<grad_base_u2.num_elements();++i)
700  {
701  grad_base_u2[i]=0;
702  grad_base_v2[i]=0;
703  grad_base_w2[i]=0;
704 
705  }
706  }
707 
708  pFields[0]->PhysDeriv(pU, grad0, grad1, grad2);
709 
710  switch (pVelocityComponent)
711  {
712  //x-equation
713  case 0:
714  //Evaluate U du'/dx
715  Vmath::Vmul (nPointsTot,grad0,1,m_base[0]->GetPhys(),1,pOutarray,1);
716  //Evaluate U du'/dx+ V du'/dy
717  Vmath::Vvtvp(nPointsTot,grad1,1,m_base[1]->GetPhys(),1,pOutarray,1,pOutarray,1);
718  //Evaluate U du'/dx+ V du'/dy+W du'/dz
719  Vmath::Vvtvp(nPointsTot,grad2,1,m_base[2]->GetPhys(),1,pOutarray,1,pOutarray,1);
720  //Evaluate -(U du'/dx+ V du'/dy+W du'/dz)
721  Vmath::Neg(nPointsTot,pOutarray,1);
722  //Evaluate -(U du'/dx+ V du'/dy+W du'/dz)+u' dU/dx
723  Vmath::Vvtvp(nPointsTot,grad_base_u0,1,pVelocity[0],1,pOutarray,1,pOutarray,1);
724  //Evaluate -(U du'/dx+ V du'/dy+W du'/dz)+u'dU/dx+ v' dV/dx
725  Vmath::Vvtvp(nPointsTot,grad_base_v0,1,pVelocity[1],1,pOutarray,1,pOutarray,1);
726  //Evaluate -(U du'/dx+ V du'/dy+W du'/dz)+u'dU/dx+ v' dV/dx+ w' dW/dz
727  Vmath::Vvtvp(nPointsTot,grad_base_w0,1,pVelocity[2],1,pOutarray,1,pOutarray,1);
728  break;
729  //y-equation
730  case 1:
731  //Evaluate U dv'/dx
732  Vmath::Vmul (nPointsTot,grad0,1,m_base[0]->GetPhys(),1,pOutarray,1);
733  //Evaluate U dv'/dx+ V dv'/dy
734  Vmath::Vvtvp(nPointsTot,grad1,1,m_base[1]->GetPhys(),1,pOutarray,1,pOutarray,1);
735  //Evaluate U dv'/dx+ V dv'/dy+W dv'/dz
736  Vmath::Vvtvp(nPointsTot,grad2,1,m_base[2]->GetPhys(),1,pOutarray,1,pOutarray,1);
737  //Evaluate -(U dv'/dx+ V dv'/dy+W dv'/dz)
738  Vmath::Neg(nPointsTot,pOutarray,1);
739  //Evaluate -(U dv'/dx+ V dv'/dy+W dv'/dz)+u' dU/dy
740  Vmath::Vvtvp(nPointsTot,grad_base_u1,1,pVelocity[0],1,pOutarray,1,pOutarray,1);
741  //Evaluate -(U dv'/dx+ V dv'/dy+W dv'/dz)+u' dU/dy +v' dV/dy
742  Vmath::Vvtvp(nPointsTot,grad_base_v1,1,pVelocity[1],1,pOutarray,1,pOutarray,1);
743  //Evaluate -(U dv'/dx+ V dv'/dy+W dv'/dz)+u' dU/dy +v' dV/dy+ w' dW/dy
744  Vmath::Vvtvp(nPointsTot,grad_base_w1,1,pVelocity[2],1,pOutarray,1,pOutarray,1);
745  break;
746 
747  //z-equation
748  case 2:
749  //Evaluate U dw'/dx
750  Vmath::Vmul (nPointsTot,grad0,1,m_base[0]->GetPhys(),1,pOutarray,1);
751  //Evaluate U dw'/dx+ V dw'/dx
752  Vmath::Vvtvp(nPointsTot,grad1,1,m_base[1]->GetPhys(),1,pOutarray,1,pOutarray,1);
753  //Evaluate U dw'/dx+ V dw'/dx+ W dw'/dz
754  Vmath::Vvtvp(nPointsTot,grad2,1,m_base[2]->GetPhys(),1,pOutarray,1,pOutarray,1);
755  //Evaluate -(U dw'/dx+ V dw'/dx+ W dw'/dz)
756  Vmath::Neg(nPointsTot,pOutarray,1);
757  //Evaluate -(U dw'/dx+ V dw'/dx+ W dw'/dz)+u' dU/dz
758  Vmath::Vvtvp(nPointsTot,grad_base_u2,1,pVelocity[0],1,pOutarray,1,pOutarray,1);
759  //Evaluate -(U dw'/dx+ V dw'/dx+ W dw'/dz)+u' dU/dz+v'dV/dz
760  Vmath::Vvtvp(nPointsTot,grad_base_v2,1,pVelocity[1],1,pOutarray,1,pOutarray,1);
761  //Evaluate -(U dw'/dx+ V dw'/dx+ W dw'/dz)+u' dU/dz+v'dV/dz + w' dW/dz
762  Vmath::Vvtvp(nPointsTot,grad_base_w2,1,pVelocity[2],1,pOutarray,1,pOutarray,1);
763  break;
764  }
765  break;
766 
767 
768  default:
769  ASSERTL0(false,"dimension unknown");
770  }
771  }
772 
773  void AdjointAdvection::UpdateBase( const NekDouble m_slices,
774  Array<OneD, const NekDouble> &inarray,
775  Array<OneD, NekDouble> &outarray,
776  const NekDouble m_time,
777  const NekDouble m_period)
778  {
779 
780  int npoints=m_base[0]->GetTotPoints();
781 
782  NekDouble BetaT=2*M_PI*fmod (m_time, m_period) / m_period;
783  NekDouble phase;
784  Array<OneD, NekDouble> auxiliary(npoints);
785 
786  Vmath::Vcopy(npoints,&inarray[0],1,&outarray[0],1);
787  Vmath::Svtvp(npoints, cos(0.5*m_slices*BetaT),&inarray[npoints],1,&outarray[0],1,&outarray[0],1);
788 
789  for (int i = 2; i < m_slices; i += 2)
790  {
791  phase = (i>>1) * BetaT;
792 
793  Vmath::Svtvp(npoints, cos(phase),&inarray[i*npoints],1,&outarray[0],1,&outarray[0],1);
794  Vmath::Svtvp(npoints, sin(phase), &inarray[(i+1)*npoints], 1, &outarray[0], 1,&outarray[0],1);
795  }
796 
797  }
798 
799  DNekBlkMatSharedPtr AdjointAdvection::GetFloquetBlockMatrix(FloquetMatType mattype, bool UseContCoeffs) const
800  {
801  DNekMatSharedPtr loc_mat;
802  DNekBlkMatSharedPtr BlkMatrix;
803  int n_exp = 0;
804 
805  n_exp = m_base[0]->GetTotPoints(); // will operatore on m_phys
806 
807  Array<OneD,unsigned int> nrows(n_exp);
808  Array<OneD,unsigned int> ncols(n_exp);
809 
810  nrows = Array<OneD, unsigned int>(n_exp,m_slices);
811  ncols = Array<OneD, unsigned int>(n_exp,m_slices);
812 
813  MatrixStorage blkmatStorage = eDIAGONAL;
814  BlkMatrix = MemoryManager<DNekBlkMat>
815  ::AllocateSharedPtr(nrows,ncols,blkmatStorage);
816 
817 
818  const LibUtilities::PointsKey Pkey(m_slices,LibUtilities::eFourierEvenlySpaced);
819  const LibUtilities::BasisKey BK(LibUtilities::eFourier,m_slices,Pkey);
820  StdRegions::StdSegExp StdSeg(BK);
821 
822  StdRegions::StdMatrixKey matkey(StdRegions::eFwdTrans,
823  StdSeg.DetShapeType(),
824  StdSeg);
825 
826  loc_mat = StdSeg.GetStdMatrix(matkey);
827 
828  // set up array of block matrices.
829  for(int i = 0; i < n_exp; ++i)
830  {
831  BlkMatrix->SetBlock(i,i,loc_mat);
832  }
833 
834  return BlkMatrix;
835  }
836 
837  //Discrete Fourier Transform for Floquet analysis
838  void AdjointAdvection::DFT(const string file, const NekDouble m_slices)
839  {
840  int npoints=m_base[0]->GetTotPoints();
841 
842  //Convected fields
843  int ConvectedFields=m_base.num_elements()-1;
844 
845  m_interp= Array<OneD, Array<OneD, NekDouble> > (ConvectedFields);
846  for(int i=0; i<ConvectedFields;++i)
847  {
848  m_interp[i]=Array<OneD,NekDouble>(npoints*m_slices);
849  }
850 
851  //Import the slides into auxiliary vector
852  //The base flow should be stored in the form filename_i.bse
853  for (int i=0; i< m_slices; ++i)
854  {
855  char chkout[16] = "";
856  sprintf(chkout, "%d", i);
857  ImportFldBase(file+"_"+chkout+".bse",m_graph,i);
858  }
859 
860 
861  // Discrete Fourier Transform of the fields
862  for(int k=0; k< ConvectedFields;++k)
863  {
864 #ifdef NEKTAR_USING_FFTW
865 
866  //Discrete Fourier Transform using FFTW
867 
868 
869  Array<OneD, NekDouble> fft_in(npoints*m_slices);
870  Array<OneD, NekDouble> fft_out(npoints*m_slices);
871 
872  Array<OneD, NekDouble> m_tmpIN(m_slices);
873  Array<OneD, NekDouble> m_tmpOUT(m_slices);
874 
875  //Shuffle the data
876  for(int j= 0; j < m_slices; ++j)
877  {
878  Vmath::Vcopy(npoints,&m_interp[k][j*npoints],1,&(fft_in[j]),m_slices);
879  }
880 
881  m_FFT = LibUtilities::GetNektarFFTFactory().CreateInstance("NekFFTW", m_slices);
882 
883  //FFT Transform
884  for(int i=0; i<npoints; i++)
885  {
886  m_FFT->FFTFwdTrans(m_tmpIN =fft_in + i*m_slices, m_tmpOUT =fft_out + i*m_slices);
887 
888  }
889 
890  //Reshuffle data
891  for(int s = 0; s < m_slices; ++s)
892  {
893  Vmath::Vcopy(npoints,&fft_out[s],m_slices,&m_interp[k][s*npoints],1);
894 
895  }
896 
897  Vmath::Zero(fft_in.num_elements(),&fft_in[0],1);
898  Vmath::Zero(fft_out.num_elements(),&fft_out[0],1);
899 #else
900  //Discrete Fourier Transform using MVM
901 
902 
903  DNekBlkMatSharedPtr blkmat;
904  blkmat = GetFloquetBlockMatrix(eForwardsPhys);
905 
906  int nrows = blkmat->GetRows();
907  int ncols = blkmat->GetColumns();
908 
909  Array<OneD, NekDouble> sortedinarray(ncols);
910  Array<OneD, NekDouble> sortedoutarray(nrows);
911 
912  //Shuffle the data
913  for(int j= 0; j < m_slices; ++j)
914  {
915  Vmath::Vcopy(npoints,&m_interp[k][j*npoints],1,&(sortedinarray[j]),m_slices);
916  }
917 
918  // Create NekVectors from the given data arrays
919  NekVector<NekDouble> in (ncols,sortedinarray,eWrapper);
920  NekVector<NekDouble> out(nrows,sortedoutarray,eWrapper);
921 
922  // Perform matrix-vector multiply.
923  out = (*blkmat)*in;
924 
925  //Reshuffle data
926  for(int s = 0; s < m_slices; ++s)
927  {
928  Vmath::Vcopy(npoints,&sortedoutarray[s],m_slices,&m_interp[k][s*npoints],1);
929  }
930 
931  Vmath::Zero(sortedinarray.num_elements(),&sortedinarray[0],1);
932  Vmath::Zero(sortedoutarray.num_elements(),&sortedoutarray[0],1);
933 
934 #endif
935 
936  //scaling of the Fourier coefficients
937  NekDouble j=-1;
938  for (int i = 2; i < m_slices; i += 2)
939  {
940  Vmath::Smul(2*npoints,j,&m_interp[k][i*npoints],1,&m_interp[k][i*npoints],1);
941  j=-j;
942 
943  }
944 
945  }
946 
947  if(m_session->DefinesParameter("period"))
948  {
949  m_period=m_session->GetParameter("period");
950  }
951  else
952  {
953  m_period=(m_session->GetParameter("TimeStep")*m_slices)/(m_slices-1.);
954  }
955 
956 
957  }
958 
959 
960 } //end of namespace
961 
Generated on Sun Mar 15 2015 00:11:29 for Nektar++ by   doxygen 1.8.1.2