Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
CalcDivergence.cpp
Go to the documentation of this file.
1 /**
2  * This function calculate the divergence starting from an .fld file.
3  * It is meant to be used with solutions produced by the
4  * incompressible Navier-Stokes solver. To use it with solutions
5  * coming form another solver further generalisations are required.
6  */
7 #include <cstdio>
8 #include <cstdlib>
9 
10 #include <MultiRegions/ExpList.h>
11 #include <MultiRegions/ExpList1D.h>
12 #include <MultiRegions/ExpList2D.h>
13 #include <MultiRegions/ExpList3D.h>
17 
18 using namespace std;
19 using namespace Nektar;
20 
21 int main(int argc, char *argv[])
22 {
23  int i,j;
24 
25  if(argc != 3)
26  {
27  fprintf(stderr,"Usage: ./CalcDivergence file.xml file.fld\n");
28  exit(1);
29  }
30 
32  = LibUtilities::SessionReader::CreateInstance(argc, argv);
33 
34 
35  //----------------------------------------------
36  // Read in mesh from input file
37  string meshfile(argv[argc-2]);
38  SpatialDomains::MeshGraphSharedPtr graphShPt = SpatialDomains::MeshGraph::Read(vSession);//meshfile);
39  //----------------------------------------------
40 
41  //----------------------------------------------
42  // Import field file.
43  string fieldfile(argv[argc-1]);
44  vector<LibUtilities::FieldDefinitionsSharedPtr> fielddef;
45  vector<vector<NekDouble> > fielddata;
46  LibUtilities::Import(fieldfile,fielddef,fielddata);
47  bool useFFT = false;
48  bool dealiasing = false;
49  //----------------------------------------------
50 
51  //----------------------------------------------
52  // Define Expansion
53  int expdim = graphShPt->GetMeshDimension();
54  int nfields = fielddef[0]->m_fields.size();
55 
57 
58  switch(expdim)
59  {
60  case 1:
61  {
62  ASSERTL0(fielddef[0]->m_numHomogeneousDir <= 2,"Quasi-3D approach is only set up for 1 or 2 homogeneous directions");
63 
64  if(fielddef[0]->m_numHomogeneousDir == 1)
65  {
67 
68  // Define Homogeneous expansion
69  //int nplanes = fielddef[0]->m_numModes[1];
70 
71  int nplanes;
72  vSession->LoadParameter("HomModesZ",nplanes,fielddef[0]->m_numModes[1]);
73 
74  // choose points to be at evenly spaced points at
76  const LibUtilities::BasisKey Bkey(fielddef[0]->m_basis[1],nplanes,Pkey);
77  NekDouble ly = fielddef[0]->m_homogeneousLengths[0];
78 
79  Exp2DH1 = MemoryManager<MultiRegions::ExpList2DHomogeneous1D>::AllocateSharedPtr(vSession,Bkey,ly,useFFT,dealiasing,graphShPt);
80  Exp[0] = Exp2DH1;
81 
82  for(i = 1; i < nfields + 1; ++i)
83  {
85  }
86  }
87  else if(fielddef[0]->m_numHomogeneousDir == 2)
88  {
90 
91  // Define Homogeneous expansion
92  //int nylines = fielddef[0]->m_numModes[1];
93  //int nzlines = fielddef[0]->m_numModes[2];
94 
95  int nylines;
96  int nzlines;
97  vSession->LoadParameter("HomModesY",nylines,fielddef[0]->m_numModes[1]);
98  vSession->LoadParameter("HomModesZ",nzlines,fielddef[0]->m_numModes[2]);
99 
100  // choose points to be at evenly spaced points at
102  const LibUtilities::BasisKey BkeyY(fielddef[0]->m_basis[1],nylines,PkeyY);
103 
105  const LibUtilities::BasisKey BkeyZ(fielddef[0]->m_basis[2],nzlines,PkeyZ);
106 
107  NekDouble ly = fielddef[0]->m_homogeneousLengths[0];
108  NekDouble lz = fielddef[0]->m_homogeneousLengths[1];
109 
110  Exp3DH2 = MemoryManager<MultiRegions::ExpList3DHomogeneous2D>::AllocateSharedPtr(vSession,BkeyY,BkeyZ,ly,lz,useFFT,dealiasing,graphShPt);
111  Exp[0] = Exp3DH2;
112 
113  for(i = 1; i < nfields + 1; ++i)
114  {
116  }
117  }
118  else
119  {
122  ::AllocateSharedPtr(vSession,graphShPt);
123  Exp[0] = Exp1D;
124  for(i = 1; i < nfields + 1; ++i)
125  {
127  ::AllocateSharedPtr(*Exp1D);
128  }
129  }
130  }
131  break;
132  case 2:
133  {
134  ASSERTL0(fielddef[0]->m_numHomogeneousDir <= 1,"NumHomogeneousDir is only set up for 1");
135 
136  if(fielddef[0]->m_numHomogeneousDir == 1)
137  {
139 
140  // Define Homogeneous expansion
141  //int nplanes = fielddef[0]->m_numModes[2];
142 
143  int nplanes;
144  vSession->LoadParameter("HomModesZ",nplanes,fielddef[0]->m_numModes[2]);
145 
146  // choose points to be at evenly spaced points at
147  // nplanes + 1 points
149  const LibUtilities::BasisKey Bkey(fielddef[0]->m_basis[2],nplanes,Pkey);
150  NekDouble lz = fielddef[0]->m_homogeneousLengths[0];
151 
152  Exp3DH1 = MemoryManager<MultiRegions::ExpList3DHomogeneous1D>::AllocateSharedPtr(vSession,Bkey,lz,useFFT,dealiasing,graphShPt);
153  Exp[0] = Exp3DH1;
154 
155  for(i = 1; i < nfields + 1; ++i)
156  {
158  }
159  }
160  else
161  {
164  ::AllocateSharedPtr(vSession,graphShPt);
165  Exp[0] = Exp2D;
166 
167  for(i = 1; i < nfields + 1; ++i)
168  {
170  ::AllocateSharedPtr(*Exp2D);
171  }
172  }
173  }
174  break;
175  case 3:
176  {
179  ::AllocateSharedPtr(vSession,graphShPt);
180  Exp[0] = Exp3D;
181 
182  for(i = 1; i < nfields + 1; ++i)
183  {
185  ::AllocateSharedPtr(*Exp3D);
186  }
187  }
188  break;
189  default:
190  ASSERTL0(false,"Expansion dimension not recognised");
191  break;
192  }
193 
194  //----------------------------------------------
195  // Copy data from field file
196  for(j = 0; j < nfields; ++j)
197  {
198  for(unsigned int i = 0; i < fielddata.size(); ++i)
199  {
200  Exp[j]->ExtractDataToCoeffs(fielddef [i],
201  fielddata[i],
202  fielddef [i]->m_fields[j],
203  Exp[j]->UpdateCoeffs());
204  }
205  Exp[j]->BwdTrans(Exp[j]->GetCoeffs(),Exp[j]->UpdatePhys());
206  }
207  //----------------------------------------------
208 
209  int nq = Exp[0]->GetNpoints();
210 
211  Array<OneD, NekDouble> Ux(nq);
212  Array<OneD, NekDouble> Vy(nq);
213  Array<OneD, NekDouble> Div(nq);
214 
215  switch(expdim)
216  {
217  case 1:
218  {
219  if(fielddef[0]->m_numHomogeneousDir == 1)
220  {
221  ASSERTL0(false,"Not implemented yet");
222  }
223  else if(fielddef[0]->m_numHomogeneousDir == 2)
224  {
225  ASSERTL0(false,"Not implemented yet");
226  }
227  else
228  {
229  ASSERTL0(false,"Not implemented yet");
230  }
231  }
232  break;
233  case 2:
234  {
235  if(fielddef[0]->m_numHomogeneousDir == 1)
236  {
237  ASSERTL0(false,"Not implemented yet");
238  }
239  else
240  {
241  Exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[0],Exp[0]->GetPhys(),Ux);
242  Exp[0]->PhysDeriv(MultiRegions::DirCartesianMap[1],Exp[1]->GetPhys(),Vy);
243 
244  Vmath::Vsub(nq,Ux,1,Vy,1,Div,1);
245 
246  Exp[3]->FwdTrans(Div,Exp[3]->UpdateCoeffs());
247  }
248  }
249  break;
250  case 3:
251  {
252  ASSERTL0(false,"Not implemented yet");
253  }
254  break;
255  default:
256  {
257  ASSERTL0(false,"Expansion dimension not recognised");
258  }
259  break;
260  }
261 
262  //-----------------------------------------------
263  // Write solution to file with additional computed fields
264  string fldfilename(argv[2]);
265  string out = fldfilename.substr(0, fldfilename.find_last_of("."));
266  string endfile("_with_divergence.fld");
267  out += endfile;
268  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
269  = Exp[0]->GetFieldDefinitions();
270  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
271 
272  for(j = 0; j < nfields + 1; ++j)
273  {
274  for(i = 0; i < FieldDef.size(); ++i)
275  {
276  if (j >= nfields)
277  {
278  FieldDef[i]->m_fields.push_back("Div");
279  }
280  else
281  {
282  FieldDef[i]->m_fields.push_back(fielddef[i]->m_fields[j]);
283  }
284  Exp[j]->AppendFieldData(FieldDef[i], FieldData[i]);
285  }
286  }
287  LibUtilities::Write(out, FieldDef, FieldData);
288  //-----------------------------------------------
289 
290  return 0;
291 }
292 
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
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
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
STL namespace.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
boost::shared_ptr< ExpList3DHomogeneous2D > ExpList3DHomogeneous2DSharedPtr
Shared pointer to an ExpList3DHomogeneous2D object.
1D Evenly-spaced points using Lagrange polynomial
Definition: PointsType.h:65
boost::shared_ptr< ExpList1D > ExpList1DSharedPtr
Shared pointer to an ExpList1D object.
Definition: ExpList1D.h:50
Defines a specification for a set of points.
Definition: Points.h:58
double NekDouble
boost::shared_ptr< ExpList2D > ExpList2DSharedPtr
Shared pointer to an ExpList2D object.
Definition: ExpList2D.h:49
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
boost::shared_ptr< ExpList2DHomogeneous1D > ExpList2DHomogeneous1DSharedPtr
Shared pointer to an ExpList2DHomogeneous1D object.
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
boost::shared_ptr< ExpList3D > ExpList3DSharedPtr
Shared pointer to an ExpList3D object.
Definition: ExpList3D.h:115
MultiRegions::Direction const DirCartesianMap[]
Definition: ExpList.h:86
boost::shared_ptr< ExpList3DHomogeneous1D > ExpList3DHomogeneous1DSharedPtr
Shared pointer to an ExpList3DHomogeneous1D object.
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:442
int main(int argc, char *argv[])
Describes the specification for a Basis.
Definition: Basis.h:50