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