Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Functions
CalcVorticity.cpp File Reference
#include <cstdio>
#include <cstdlib>
#include <MultiRegions/ExpList.h>
#include <MultiRegions/ExpList1D.h>
#include <MultiRegions/ExpList2D.h>
#include <MultiRegions/ExpList3D.h>
#include <MultiRegions/ExpList2DHomogeneous1D.h>
#include <MultiRegions/ExpList3DHomogeneous1D.h>
#include <MultiRegions/ExpList3DHomogeneous2D.h>
Include dependency graph for CalcVorticity.cpp:

Go to the source code of this file.

Functions

int main (int argc, char *argv[])
 

Function Documentation

int main ( int  argc,
char *  argv[] 
)

Definition at line 20 of file CalcVorticity.cpp.

References ASSERTL0, Nektar::MultiRegions::DirCartesianMap, Nektar::LibUtilities::ePolyEvenlySpaced, Nektar::LibUtilities::Import(), Vmath::Vsub(), and Nektar::LibUtilities::Write().

21 {
22  int i,j;
23 
24  if(argc != 3)
25  {
26  fprintf(stderr,"Usage: ./CalcVorticity file.xml file.fld\n");
27  exit(1);
28  }
29 
30  LibUtilities::SessionReaderSharedPtr vSession
31  = LibUtilities::SessionReader::CreateInstance(argc, argv);
32 
33 
34  //----------------------------------------------
35  // Read in mesh from input file
36  string meshfile(argv[argc-2]);
37  SpatialDomains::MeshGraphSharedPtr graphShPt = SpatialDomains::MeshGraph::Read(vSession);//meshfile);
38  //----------------------------------------------
39 
40  //----------------------------------------------
41  // Import field file.
42  string fieldfile(argv[argc-1]);
43  vector<LibUtilities::FieldDefinitionsSharedPtr> fielddef;
44  vector<vector<NekDouble> > fielddata;
45  LibUtilities::Import(fieldfile,fielddef,fielddata);
46  bool useFFT = false;
47  bool dealiasing = false;
48  //----------------------------------------------
49 
50  //----------------------------------------------
51  // Define Expansion
52  int expdim = graphShPt->GetMeshDimension();
53  int nfields = fielddef[0]->m_fields.size();
54  int vorticitydim;
55  if(expdim == 1)
56  {
57  if(fielddef[0]->m_numHomogeneousDir == 2)//3D Homogeneous 2D
58  {
59  vorticitydim = 3;
60  }
61  else // 1D
62  {
63  vorticitydim = 0;
64  }
65  }
66  else if(expdim ==2)
67  {
68  if(fielddef[0]->m_numHomogeneousDir == 1)// 3D Homogeneous 1D
69  {
70  vorticitydim = 3;
71  }
72  else //2D
73  {
74  vorticitydim = 1;
75  }
76 
77  }
78  else // Full 3D
79  {
80  vorticitydim = 3;
81  }
82 
83 
84  Array<OneD, MultiRegions::ExpListSharedPtr> Exp(nfields + vorticitydim);
85 
86  switch(expdim)
87  {
88  case 1:
89  {
90  ASSERTL0(fielddef[0]->m_numHomogeneousDir <= 2,"Quasi-3D approach is only set up for 1 or 2 homogeneous directions");
91 
92  if(fielddef[0]->m_numHomogeneousDir == 1)
93  {
94  MultiRegions::ExpList2DHomogeneous1DSharedPtr Exp2DH1;
95 
96  // Define Homogeneous expansion
97  //int nplanes = fielddef[0]->m_numModes[1];
98 
99  int nplanes;
100  vSession->LoadParameter("HomModesZ",nplanes,fielddef[0]->m_numModes[1]);
101 
102  // choose points to be at evenly spaced points at
103  const LibUtilities::PointsKey Pkey(nplanes+1,LibUtilities::ePolyEvenlySpaced);
104  const LibUtilities::BasisKey Bkey(fielddef[0]->m_basis[1],nplanes,Pkey);
105  NekDouble ly = fielddef[0]->m_homogeneousLengths[0];
106 
107  Exp2DH1 = MemoryManager<MultiRegions::ExpList2DHomogeneous1D>::AllocateSharedPtr(vSession,Bkey,ly,useFFT,dealiasing,graphShPt);
108  Exp[0] = Exp2DH1;
109 
110  for(i = 1; i < nfields + vorticitydim; ++i)
111  {
112  Exp[i] = MemoryManager<MultiRegions::ExpList2DHomogeneous1D>::AllocateSharedPtr(*Exp2DH1);
113  }
114  }
115  else if(fielddef[0]->m_numHomogeneousDir == 2)
116  {
117  MultiRegions::ExpList3DHomogeneous2DSharedPtr Exp3DH2;
118 
119  // Define Homogeneous expansion
120  //int nylines = fielddef[0]->m_numModes[1];
121  //int nzlines = fielddef[0]->m_numModes[2];
122 
123  int nylines;
124  int nzlines;
125  vSession->LoadParameter("HomModesY",nylines,fielddef[0]->m_numModes[1]);
126  vSession->LoadParameter("HomModesZ",nzlines,fielddef[0]->m_numModes[2]);
127 
128  // choose points to be at evenly spaced points at
129  const LibUtilities::PointsKey PkeyY(nylines+1,LibUtilities::ePolyEvenlySpaced);
130  const LibUtilities::BasisKey BkeyY(fielddef[0]->m_basis[1],nylines,PkeyY);
131 
132  const LibUtilities::PointsKey PkeyZ(nzlines+1,LibUtilities::ePolyEvenlySpaced);
133  const LibUtilities::BasisKey BkeyZ(fielddef[0]->m_basis[2],nzlines,PkeyZ);
134 
135  NekDouble ly = fielddef[0]->m_homogeneousLengths[0];
136  NekDouble lz = fielddef[0]->m_homogeneousLengths[1];
137 
138  Exp3DH2 = MemoryManager<MultiRegions::ExpList3DHomogeneous2D>::AllocateSharedPtr(vSession,BkeyY,BkeyZ,ly,lz,useFFT,dealiasing,graphShPt);
139  Exp[0] = Exp3DH2;
140 
141  for(i = 1; i < nfields + vorticitydim; ++i)
142  {
143  Exp[i] = MemoryManager<MultiRegions::ExpList3DHomogeneous2D>::AllocateSharedPtr(*Exp3DH2);
144  }
145  }
146  else
147  {
148  MultiRegions::ExpList1DSharedPtr Exp1D;
149  Exp1D = MemoryManager<MultiRegions::ExpList1D>
150  ::AllocateSharedPtr(vSession,graphShPt);
151  Exp[0] = Exp1D;
152  for(i = 1; i < nfields + vorticitydim; ++i)
153  {
154  Exp[i] = MemoryManager<MultiRegions::ExpList1D>
155  ::AllocateSharedPtr(*Exp1D);
156  }
157  }
158  }
159  break;
160  case 2:
161  {
162  ASSERTL0(fielddef[0]->m_numHomogeneousDir <= 1,"NumHomogeneousDir is only set up for 1");
163 
164  if(fielddef[0]->m_numHomogeneousDir == 1)
165  {
166  MultiRegions::ExpList3DHomogeneous1DSharedPtr Exp3DH1;
167 
168  // Define Homogeneous expansion
169  //int nplanes = fielddef[0]->m_numModes[2];
170 
171  int nplanes;
172  vSession->LoadParameter("HomModesZ",nplanes,fielddef[0]->m_numModes[2]);
173 
174  // choose points to be at evenly spaced points at
175  // nplanes + 1 points
176  const LibUtilities::PointsKey Pkey(nplanes+1,LibUtilities::ePolyEvenlySpaced);
177  const LibUtilities::BasisKey Bkey(fielddef[0]->m_basis[2],nplanes,Pkey);
178  NekDouble lz = fielddef[0]->m_homogeneousLengths[0];
179 
180  Exp3DH1 = MemoryManager<MultiRegions::ExpList3DHomogeneous1D>::AllocateSharedPtr(vSession,Bkey,lz,useFFT,dealiasing,graphShPt);
181  Exp[0] = Exp3DH1;
182 
183  for(i = 1; i < nfields + vorticitydim; ++i)
184  {
185  Exp[i] = MemoryManager<MultiRegions::ExpList3DHomogeneous1D>::AllocateSharedPtr(*Exp3DH1);
186  }
187  }
188  else
189  {
190  MultiRegions::ExpList2DSharedPtr Exp2D;
191  Exp2D = MemoryManager<MultiRegions::ExpList2D>
192  ::AllocateSharedPtr(vSession,graphShPt);
193  Exp[0] = Exp2D;
194 
195  for(i = 1; i < nfields + vorticitydim; ++i)
196  {
197  Exp[i] = MemoryManager<MultiRegions::ExpList2D>
198  ::AllocateSharedPtr(*Exp2D);
199  }
200  }
201  }
202  break;
203  case 3:
204  {
205  MultiRegions::ExpList3DSharedPtr Exp3D;
206  Exp3D = MemoryManager<MultiRegions::ExpList3D>
207  ::AllocateSharedPtr(vSession,graphShPt);
208  Exp[0] = Exp3D;
209 
210  for(i = 1; i < nfields + vorticitydim; ++i)
211  {
212  Exp[i] = MemoryManager<MultiRegions::ExpList3D>
213  ::AllocateSharedPtr(*Exp3D);
214  }
215  }
216  break;
217  default:
218  ASSERTL0(false,"Expansion dimension not recognised");
219  break;
220  }
221 
222  //----------------------------------------------
223  // Copy data from field file
224  for(j = 0; j < nfields; ++j)
225  {
226  for(unsigned int i = 0; i < fielddata.size(); ++i)
227  {
228  Exp[j]->ExtractDataToCoeffs(fielddef [i],
229  fielddata[i],
230  fielddef [i]->m_fields[j],
231  Exp[j]->UpdateCoeffs());
232  }
233  Exp[j]->BwdTrans(Exp[j]->GetCoeffs(),Exp[j]->UpdatePhys());
234  }
235  //----------------------------------------------
236 
237  int nq = Exp[0]->GetNpoints();
238 
239  Array<OneD, NekDouble> Uy(nq);
240  Array<OneD, NekDouble> Uz(nq);
241  Array<OneD, NekDouble> Vx(nq);
242  Array<OneD, NekDouble> Vz(nq);
243  Array<OneD, NekDouble> Wx(nq);
244  Array<OneD, NekDouble> Wy(nq);
245  Array<OneD, NekDouble> Qx(nq);
246  Array<OneD, NekDouble> Qy(nq);
247  Array<OneD, NekDouble> Qz(nq);
248 
249  switch(expdim)
250  {
251  case 1:
252  {
253  if(fielddef[0]->m_numHomogeneousDir == 1)
254  {
255  ASSERTL0(false,"Not implemented yet");
256  }
257  else if(fielddef[0]->m_numHomogeneousDir == 2)
258  {
259  Exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],Exp[2]->GetPhys(),Wy);
260  Exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],Exp[1]->GetPhys(),Vz);
261  Exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],Exp[0]->GetPhys(),Uz);
262  Exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],Exp[2]->GetPhys(),Wx);
263  Exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],Exp[1]->GetPhys(),Vx);
264  Exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],Exp[0]->GetPhys(),Uy);
265 
266  Vmath::Vsub(nq,Wy,1,Vz,1,Qx,1);
267  Vmath::Vsub(nq,Uz,1,Wx,1,Qy,1);
268  Vmath::Vsub(nq,Vx,1,Uy,1,Qz,1);
269 
270  Exp[4]->FwdTrans(Qx,Exp[4]->UpdateCoeffs());
271  Exp[5]->FwdTrans(Qy,Exp[5]->UpdateCoeffs());
272  Exp[6]->FwdTrans(Qz,Exp[6]->UpdateCoeffs());
273  }
274  else
275  {
276  ASSERTL0(false,"Not implemented yet");
277  }
278  }
279  break;
280  case 2:
281  {
282  if(fielddef[0]->m_numHomogeneousDir == 1)
283  {
284  Exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],Exp[2]->GetPhys(),Wy);
285  Exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],Exp[1]->GetPhys(),Vz);
286  Exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],Exp[0]->GetPhys(),Uz);
287  Exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],Exp[2]->GetPhys(),Wx);
288  Exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],Exp[1]->GetPhys(),Vx);
289  Exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],Exp[0]->GetPhys(),Uy);
290 
291  Vmath::Vsub(nq,Wy,1,Vz,1,Qx,1);
292  Vmath::Vsub(nq,Uz,1,Wx,1,Qy,1);
293  Vmath::Vsub(nq,Vx,1,Uy,1,Qz,1);
294 
295  Exp[4]->FwdTrans(Qx,Exp[4]->UpdateCoeffs());
296  Exp[5]->FwdTrans(Qy,Exp[5]->UpdateCoeffs());
297  Exp[6]->FwdTrans(Qz,Exp[6]->UpdateCoeffs());
298  }
299  else
300  {
301  Exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],Exp[1]->GetPhys(),Vx);
302  Exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],Exp[0]->GetPhys(),Uy);
303 
304  Vmath::Vsub(nq,Vx,1,Uy,1,Qz,1);
305 
306  Exp[3]->FwdTrans(Qz,Exp[3]->UpdateCoeffs());
307  }
308  }
309  break;
310  case 3:
311  {
312  Exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],Exp[2]->GetPhys(),Wy);
313  Exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],Exp[1]->GetPhys(),Vz);
314  Exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[2],Exp[0]->GetPhys(),Uz);
315  Exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],Exp[2]->GetPhys(),Wx);
316  Exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],Exp[1]->GetPhys(),Vx);
317  Exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],Exp[0]->GetPhys(),Uy);
318 
319  Vmath::Vsub(nq,Wy,1,Vz,1,Qx,1);
320  Vmath::Vsub(nq,Uz,1,Wx,1,Qy,1);
321  Vmath::Vsub(nq,Vx,1,Uy,1,Qz,1);
322 
323  Exp[4]->FwdTrans(Qx,Exp[4]->UpdateCoeffs());
324  Exp[5]->FwdTrans(Qy,Exp[5]->UpdateCoeffs());
325  Exp[6]->FwdTrans(Qz,Exp[6]->UpdateCoeffs());
326  }
327  break;
328  default:
329  {
330  ASSERTL0(false,"Expansion dimension not recognised");
331  }
332  break;
333  }
334 
335  //-----------------------------------------------
336  // Write solution to file with additional computed fields
337  string fldfilename(argv[2]);
338  string out = fldfilename.substr(0, fldfilename.find_last_of("."));
339  string endfile("_with_vorticity.fld");
340  out += endfile;
341  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
342  = Exp[0]->GetFieldDefinitions();
343  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
344 
345  for(j = 0; j < nfields + vorticitydim; ++j)
346  {
347  for(i = 0; i < FieldDef.size(); ++i)
348  {
349  if (j >= nfields)
350  {
351  if(j == 4)
352  {
353  FieldDef[i]->m_fields.push_back("Qx");
354  }
355  else if(j == 5)
356  {
357  FieldDef[i]->m_fields.push_back("Qy");
358  }
359  else
360  {
361  FieldDef[i]->m_fields.push_back("Qz");
362  }
363  }
364  else
365  {
366  FieldDef[i]->m_fields.push_back(fielddef[i]->m_fields[j]);
367  }
368  Exp[j]->AppendFieldData(FieldDef[i], FieldData[i]);
369  }
370  }
371  LibUtilities::Write(out, FieldDef, FieldData);
372  //-----------------------------------------------
373 
374  return 0;
375 }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
Nektar::LibUtilities::Import
void Import(const std::string &infilename, std::vector< FieldDefinitionsSharedPtr > &fielddefs, std::vector< std::vector< NekDouble > > &fielddata, FieldMetaDataMap &fieldinfomap, const Array< OneD, int > &ElementIDs)
This function allows for data to be imported from an FLD file when a session and/or communicator is n...
Definition: FieldIO.cpp:279
Nektar::MemoryManager
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
Definition: NekMemoryManager.hpp:102
Nektar::LibUtilities::SessionReaderSharedPtr
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
Nektar::MultiRegions::ExpList3DHomogeneous2DSharedPtr
boost::shared_ptr< ExpList3DHomogeneous2D > ExpList3DHomogeneous2DSharedPtr
Shared pointer to an ExpList3DHomogeneous2D object.
Definition: ExpList3DHomogeneous2D.h:50
Nektar::LibUtilities::ePolyEvenlySpaced
1D Evenly-spaced points using Lagrange polynomial
Definition: PointsType.h:65
Nektar::MultiRegions::ExpList1DSharedPtr
boost::shared_ptr< ExpList1D > ExpList1DSharedPtr
Shared pointer to an ExpList1D object.
Definition: ExpList1D.h:50
Nektar::LibUtilities::PointsKey
Defines a specification for a set of points.
Definition: Points.h:58
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Nektar::MultiRegions::ExpList2DSharedPtr
boost::shared_ptr< ExpList2D > ExpList2DSharedPtr
Shared pointer to an ExpList2D object.
Definition: ExpList2D.h:49
Nektar::LibUtilities::Write
void Write(const std::string &outFile, std::vector< FieldDefinitionsSharedPtr > &fielddefs, std::vector< std::vector< NekDouble > > &fielddata, const FieldMetaDataMap &fieldinfomap, const bool backup)
This function allows for data to be written to an FLD file when a session and/or communicator is not ...
Definition: FieldIO.cpp:235
Nektar::MultiRegions::ExpList2DHomogeneous1DSharedPtr
boost::shared_ptr< ExpList2DHomogeneous1D > ExpList2DHomogeneous1DSharedPtr
Shared pointer to an ExpList2DHomogeneous1D object.
Definition: ExpList2DHomogeneous1D.h:50
Vmath::Vsub
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:343
Nektar::MultiRegions::ExpList3DSharedPtr
boost::shared_ptr< ExpList3D > ExpList3DSharedPtr
Shared pointer to an ExpList3D object.
Definition: ExpList3D.h:112
Nektar::MultiRegions::DirCartesianMap
MultiRegions::Direction const DirCartesianMap[]
Definition: ExpList.h:86
Nektar::MultiRegions::ExpList3DHomogeneous1DSharedPtr
boost::shared_ptr< ExpList3DHomogeneous1D > ExpList3DHomogeneous1DSharedPtr
Shared pointer to an ExpList3DHomogeneous1D object.
Definition: ExpList3DHomogeneous1D.h:50
Nektar::SpatialDomains::MeshGraphSharedPtr
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:442
Nektar::LibUtilities::BasisKey
Describes the specification for a Basis.
Definition: Basis.h:50
Generated on Sat Nov 25 2017 16:18:21 for Nektar++ by   doxygen 1.8.8