Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ExtractSurface2DCFS.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File ExtractSurface2DCFS.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: Extract 1D surfaces from 2D file and output relevant quantities
33 // for compressible flow solver.
34 //
35 ///////////////////////////////////////////////////////////////////////////////
36 
37 #include <cstdio>
38 #include <cstdlib>
39 #include <string>
40 #include <iostream>
41 #include <iomanip>
42 
43 #include <MultiRegions/ExpList.h>
44 #include <MultiRegions/ExpList1D.h>
45 #include <MultiRegions/ExpList2D.h>
46 #include <MultiRegions/ExpList2DHomogeneous1D.h>
47 #include <MultiRegions/ExpList3DHomogeneous1D.h>
48 #include <MultiRegions/ExpList3DHomogeneous2D.h>
49 #include <MultiRegions/AssemblyMap/AssemblyMapDG.h>
50 
51 #include <MultiRegions/DisContField2D.h>
52 #include <LocalRegions/MatrixKey.h>
53 #include <LocalRegions/Expansion2D.h>
54 #include <LocalRegions/Expansion.h>
55 
56 #include <LibUtilities/BasicUtils/FieldIO.h>
57 #include <LibUtilities/BasicUtils/NekFactory.hpp>
58 #include <LibUtilities/BasicUtils/SessionReader.h>
59 #include <LibUtilities/BasicUtils/SharedArray.hpp>
60 #include <LibUtilities/Communication/Comm.h>
61 
62 #include <LibUtilities/Memory/NekMemoryManager.hpp>
63 #include <MultiRegions/ContField2D.h>
64 #include <SpatialDomains/MeshGraph2D.h>
65 
66 #include <SolverUtils/SolverUtilsDeclspec.h>
67 
68 // #include <SolverUtils/Filters/Filter.h>
69 
70 using namespace Nektar;
71 
72 int main(int argc, char *argv[])
73 {
74  string fname = std::string(argv[2]);
75  int fdot = fname.find_last_of('.');
76  if (fdot != std::string::npos)
77  {
78  string ending = fname.substr(fdot);
79 
80  // If .chk or .fld we exchange the extension in the output file.
81  // For all other files (e.g. .bse) we append the extension to avoid
82  // conflicts.
83  if (ending == ".chk" || ending == ".fld")
84  {
85  fname = fname.substr(0,fdot);
86  }
87  }
88 
89  fname = fname + ".txt";
90 
91  int cnt;
92  int id1 = 0;
93  int id2 = 0;
94  int i, j, n, e, b;
95  Array<OneD, NekDouble> auxArray;
96 
97  int nBndEdgePts, nBndEdges, nBndRegions;
98 
99  if (argc < 3)
100  {
101  fprintf(stderr,
102  "Usage: ExtractSurface2DCFS meshfile fieldFile\n");
103  fprintf(stderr,
104  "Extracts a surface from a 2D fld file"
105  "(only for CompressibleFlowSolver and purely 2D .fld files)\n");
106  exit(1);
107  }
108 
109  LibUtilities::SessionReaderSharedPtr vSession
110  = LibUtilities::SessionReader::CreateInstance(3, argv);
111 
112  std::string m_ViscosityType;
113 
114  NekDouble m_gamma;
115  NekDouble m_pInf;
116  NekDouble m_rhoInf;
117  NekDouble m_uInf;
118  NekDouble m_vInf;
119  NekDouble m_wInf;
120  NekDouble m_gasConstant;
121  NekDouble m_Twall;
122  NekDouble m_mu;
123  NekDouble m_thermalConductivity;
124 
125  int m_spacedim = 2;
126  int nDimensions = m_spacedim;
127  int phys_offset;
128 
129  // Get gamma parameter from session file.
130  ASSERTL0(vSession->DefinesParameter("Gamma"),
131  "Compressible flow sessions must define a Gamma parameter.");
132  vSession->LoadParameter("Gamma", m_gamma, 1.4);
133 
134  // Get E0 parameter from session file.
135  ASSERTL0(vSession->DefinesParameter("pInf"),
136  "Compressible flow sessions must define a pInf parameter.");
137  vSession->LoadParameter("pInf", m_pInf, 101325);
138 
139  // Get rhoInf parameter from session file.
140  ASSERTL0(vSession->DefinesParameter("rhoInf"),
141  "Compressible flow sessions must define a rhoInf parameter.");
142  vSession->LoadParameter("rhoInf", m_rhoInf, 1.225);
143 
144  // Get uInf parameter from session file.
145  ASSERTL0(vSession->DefinesParameter("uInf"),
146  "Compressible flow sessions must define a uInf parameter.");
147  vSession->LoadParameter("uInf", m_uInf, 0.1);
148 
149  // Get vInf parameter from session file.
150  if (m_spacedim == 2 || m_spacedim == 3)
151  {
152  ASSERTL0(vSession->DefinesParameter("vInf"),
153  "Compressible flow sessions must define a vInf parameter"
154  "for 2D/3D problems.");
155  vSession->LoadParameter("vInf", m_vInf, 0.0);
156  }
157 
158  // Get wInf parameter from session file.
159  if (m_spacedim == 3)
160  {
161  ASSERTL0(vSession->DefinesParameter("wInf"),
162  "Compressible flow sessions must define a wInf parameter"
163  "for 3D problems.");
164  vSession->LoadParameter("wInf", m_wInf, 0.0);
165  }
166 
167  vSession->LoadParameter ("GasConstant", m_gasConstant, 287.058);
168  vSession->LoadParameter ("Twall", m_Twall, 300.15);
169  vSession->LoadSolverInfo("ViscosityType", m_ViscosityType, "Constant");
170  vSession->LoadParameter ("mu", m_mu, 1.78e-05);
171  vSession->LoadParameter ("thermalConductivity",
172  m_thermalConductivity, 0.0257);
173 
174  //--------------------------------------------------------------------------
175  // Read in mesh from input file
176  string meshfile(argv[1]);
177  SpatialDomains::MeshGraphSharedPtr graphShPt =
178  SpatialDomains::MeshGraph::Read(vSession);
179  //--------------------------------------------------------------------------
180 
181  //--------------------------------------------------------------------------
182  // Import field file
183  string fieldFile(argv[2]);
184  vector<LibUtilities::FieldDefinitionsSharedPtr> fieldDef;
185  vector<vector<NekDouble> > fieldData;
186 
187  LibUtilities::Import(fieldFile, fieldDef, fieldData);
188  //--------------------------------------------------------------------------
189 
190  //--------------------------------------------------------------------------
191  // Set up Expansion information
192  vector< vector<LibUtilities::PointsType> > pointsType;
193  for (i = 0; i < fieldDef.size(); ++i)
194  {
195  vector<LibUtilities::PointsType> ptype;
196  for (j = 0; j < 2; ++j)
197  {
198  ptype.push_back(LibUtilities::ePolyEvenlySpaced);
199  }
200  pointsType.push_back(ptype);
201  }
202  graphShPt->SetExpansions(fieldDef, pointsType);
203 
204  //--------------------------------------------------------------------------
205 
206 
207  //--------------------------------------------------------------------------
208  // Define Expansion
209  int nfields = fieldDef[0]->m_fields.size();
210  Array<OneD, MultiRegions::ExpListSharedPtr> Exp(nfields);
211  Array<OneD, MultiRegions::ExpListSharedPtr> pFields(nfields);
212 
213  for(i = 0; i < pFields.num_elements(); i++)
214  {
215  pFields[i] = MemoryManager<
216  MultiRegions::DisContField2D>::AllocateSharedPtr(
217  vSession, graphShPt, vSession->GetVariable(i));
218  }
219 
220  MultiRegions::ExpList2DSharedPtr Exp2D;
221  Exp2D = MemoryManager<MultiRegions::ExpList2D>
222  ::AllocateSharedPtr(vSession, graphShPt);
223 
224  Exp[0] = Exp2D;
225 
226  for (i = 1; i < nfields; ++i)
227  {
228  Exp[i] = MemoryManager<MultiRegions::ExpList2D>
229  ::AllocateSharedPtr(*Exp2D);
230  }
231 
232  int nSolutionPts = pFields[0]->GetNpoints();
233  int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
234  int nElements = pFields[0]->GetExpSize();
235 
236  Array<OneD, NekDouble> tmp(nSolutionPts, 0.0);
237 
238  Array<OneD, NekDouble> x(nSolutionPts);
239  Array<OneD, NekDouble> y(nSolutionPts);
240  Array<OneD, NekDouble> z(nSolutionPts);
241 
242  Array<OneD, NekDouble> traceX(nTracePts);
243  Array<OneD, NekDouble> traceY(nTracePts);
244  Array<OneD, NekDouble> traceZ(nTracePts);
245 
246  Array<OneD, NekDouble> surfaceX(nTracePts);
247  Array<OneD, NekDouble> surfaceY(nTracePts);
248  Array<OneD, NekDouble> surfaceZ(nTracePts);
249 
250  pFields[0]->GetCoords(x, y, z);
251 
252  pFields[0]->ExtractTracePhys(x, traceX);
253  pFields[0]->ExtractTracePhys(y, traceY);
254  pFields[0]->ExtractTracePhys(z, traceZ);
255  //--------------------------------------------------------------------------
256 
257  //--------------------------------------------------------------------------
258  // Copy data from field file
259  Array<OneD, Array<OneD, NekDouble> > uFields(nfields);
260  Array<OneD, Array<OneD, NekDouble> > traceFields(nfields);
261  Array<OneD, Array<OneD, NekDouble> > surfaceFields(nfields);
262 
263  // Extract the physical values of the solution at the boundaries
264  for (j = 0; j < nfields; ++j)
265  {
266  uFields[j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
267  traceFields[j] = Array<OneD, NekDouble>(nTracePts, 0.0);
268  surfaceFields[j] = Array<OneD, NekDouble>(nTracePts, 0.0);
269 
270 
271  for (i = 0; i < fieldData.size(); ++i)
272  {
273  Exp[j]->ExtractDataToCoeffs(fieldDef[i], fieldData[i],
274  fieldDef[i]->m_fields[j],
275  Exp[j]->UpdateCoeffs());
276  }
277  Exp[j]->BwdTrans(Exp[j]->GetCoeffs(), Exp[j]->UpdatePhys());
278  Vmath::Vcopy(nSolutionPts, Exp[j]->GetPhys(), 1, uFields[j], 1);
279  pFields[0]->ExtractTracePhys(uFields[j], traceFields[j]);
280  }
281 
282  //Fields to add in the output file
283 
284  int nfieldsAdded = 20;
285  Array<OneD, Array<OneD, NekDouble> > traceFieldsAdded(nfieldsAdded);
286  Array<OneD, Array<OneD, NekDouble> > surfaceFieldsAdded(nfieldsAdded);
287 
288  for (j = 0; j < nfieldsAdded; ++j)
289  {
290  traceFieldsAdded[j] = Array<OneD, NekDouble>(nTracePts, 0.0);
291  surfaceFieldsAdded[j] = Array<OneD, NekDouble>(nTracePts, 0.0);
292  }
293 
294  /******** Evaluation of normals and tangents on the trace *****************
295  * nx -> traceFieldsAdded[0];
296  * ny -> traceFieldsAdded[1];
297  * tx -> traceFieldsAdded[2];
298  * ty -> traceFieldsAdded[3];
299  ***************************************************************************/
300 
301  Array<OneD, Array<OneD, NekDouble> > m_traceNormals (nDimensions);
302  for(i = 0; i < nDimensions; ++i)
303  {
304  m_traceNormals[i] = Array<OneD, NekDouble> (nTracePts, 0.0);
305  }
306  pFields[0]->GetTrace()->GetNormals(m_traceNormals);
307 
308  Array<OneD, Array<OneD, NekDouble> > m_traceTangents (nDimensions);
309  for(i = 0; i < nDimensions; ++i)
310  {
311  m_traceTangents[i] = Array<OneD, NekDouble> (nTracePts, 0.0);
312  }
313 
314 
315  // nx
316  Vmath::Vcopy(nTracePts,
317  &m_traceNormals[0][0], 1,
318  &traceFieldsAdded[0][0], 1);
319 
320  // ny
321  Vmath::Vcopy(nTracePts,
322  &m_traceNormals[1][0], 1,
323  &traceFieldsAdded[1][0], 1);
324 
325  // t_x = - n_y
326  Vmath::Vcopy(nTracePts,
327  &m_traceNormals[1][0], 1,
328  &m_traceTangents[0][0], 1);
329  Vmath::Neg(nTracePts, &m_traceTangents[0][0], 1);
330 
331  Vmath::Vcopy(nTracePts,
332  &m_traceTangents[0][0], 1,
333  &traceFieldsAdded[2][0], 1);
334 
335  // t_y = n_x
336  Vmath::Vcopy(nTracePts,
337  &m_traceNormals[0][0], 1,
338  &m_traceTangents[1][0], 1);
339 
340  Vmath::Vcopy(nTracePts,
341  &m_traceTangents[1][0], 1,
342  &traceFieldsAdded[3][0], 1);
343 
344  /******** Evaluation of the pressure ***************************************
345  * P = (E-1/2.*rho.*((rhou./rho).^2+(rhov./rho).^2))*(gamma - 1);
346  * P -> traceFieldsAdded[4];
347  ***************************************************************************/
348 
349  Array<OneD, NekDouble> pressure(nSolutionPts, 0.0);
350  NekDouble gammaMinusOne = m_gamma - 1.0;
351 
352  for (i = 0; i < m_spacedim; i++)
353  {
354  Vmath::Vmul(nSolutionPts,
355  &uFields[i + 1][0], 1,
356  &uFields[i + 1][0], 1,
357  &tmp[0],1);
358 
359 
360  Vmath::Smul(nSolutionPts, 0.5,
361  &tmp[0], 1,
362  &tmp[0], 1);
363 
364  Vmath::Vadd(nSolutionPts,
365  &pressure[0], 1,
366  &tmp[0], 1,
367  &pressure[0], 1);
368  }
369 
370  Vmath::Vdiv(nSolutionPts,
371  &pressure[0], 1,
372  &uFields[0][0], 1,
373  &pressure[0],1);
374 
375  Vmath::Vsub(nSolutionPts,
376  &uFields[nfields - 1][0], 1,
377  &pressure[0], 1,
378  &pressure[0],1);
379 
380  Vmath::Smul(nSolutionPts, gammaMinusOne,
381  &pressure[0], 1,
382  &pressure[0], 1);
383 
384  // Extract trace
385  pFields[0]->ExtractTracePhys(pressure, traceFieldsAdded[4]);
386 
387  /******** Evaluation of the temperature ************************************
388  * T = P/(R*rho);
389  * T -> traceFieldsAdded[5];
390  ***************************************************************************/
391 
392  Array<OneD, NekDouble> temperature(nSolutionPts, 0.0);
393 
394  Vmath::Vdiv(nSolutionPts,
395  &pressure[0], 1,
396  &uFields[0][0], 1,
397  &temperature[0],1);
398 
399  NekDouble GasConstantInv = 1.0/m_gasConstant;
400  Vmath::Smul(nSolutionPts, GasConstantInv,
401  &temperature[0], 1,
402  &temperature[0], 1);
403 
404  // Extract trace
405  pFields[0]->ExtractTracePhys(temperature, traceFieldsAdded[5]);
406 
407  /*** Evaluation of the temperature gradient in the normal direction ********
408  * DT_n -> traceFieldsAdded[6]
409  ***************************************************************************/
410 
411  Array<OneD, Array<OneD, NekDouble> > Dtemperature(nDimensions);
412  Array<OneD, Array<OneD, NekDouble> > traceDtemperature(nDimensions);
413 
414  for (i = 0; i < nDimensions; ++ i)
415  {
416  Dtemperature[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
417  traceDtemperature[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
418  }
419 
420  for (i = 0; i < nDimensions; ++ i)
421  {
422  for (n = 0; n < nElements; n++)
423  {
424  phys_offset = pFields[0]->GetPhys_Offset(n);
425 
426  pFields[i]->GetExp(n)->PhysDeriv(
427  i, temperature + phys_offset,
428  auxArray = Dtemperature[i] + phys_offset);
429  }
430  // Extract trace
431  pFields[0]->ExtractTracePhys(Dtemperature[i], traceDtemperature[i]);
432  }
433 
434  for(i = 0; i < nDimensions; ++i)
435  {
436  Vmath::Vmul(nTracePts,
437  &m_traceNormals[i][0], 1,
438  &traceDtemperature[i][0], 1,
439  &tmp[0],1);
440 
441  Vmath::Vadd(nTracePts,
442  &traceFieldsAdded[6][0], 1,
443  &tmp[0], 1,
444  &traceFieldsAdded[6][0], 1);
445  }
446 
447  /*** Evaluation of the pressure gradient ***********************************
448  * DP_t -> traceFieldsAdded[7] tangent direction
449  * DP_x -> traceFieldsAdded[8]
450  * DP_y -> traceFieldsAdded[9]
451  ***************************************************************************/
452 
453  Array<OneD, Array<OneD, NekDouble> > Dpressure(nDimensions);
454  Array<OneD, Array<OneD, NekDouble> > traceDpressure(nDimensions);
455 
456  for (i = 0; i < nDimensions; ++ i)
457  {
458  Dpressure[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
459  traceDpressure[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
460  }
461 
462  for (i = 0; i < nDimensions; ++ i)
463  {
464  for (n = 0; n < nElements; n++)
465  {
466  phys_offset = pFields[0]->GetPhys_Offset(n);
467 
468  pFields[i]->GetExp(n)->PhysDeriv(
469  i, pressure + phys_offset,
470  auxArray = Dpressure[i] + phys_offset);
471  }
472  // Extract trace
473  pFields[0]->ExtractTracePhys(Dpressure[i], traceDpressure[i]);
474  }
475 
476  // Dp_t
477  for(i = 0; i < nDimensions; ++i)
478  {
479  Vmath::Vmul(nTracePts,
480  &m_traceTangents[i][0], 1,
481  &traceDpressure[i][0], 1,
482  &tmp[0],1);
483 
484  Vmath::Vadd(nTracePts,
485  &traceFieldsAdded[7][0], 1,
486  &tmp[0], 1,
487  &traceFieldsAdded[7][0], 1);
488  }
489 
490  // Dp_x
491  Vmath::Vcopy(nTracePts,
492  &traceDpressure[0][0], 1,
493  &traceFieldsAdded[8][0], 1);
494 
495  // Dp_y
496  Vmath::Vcopy(nTracePts,
497  &traceDpressure[1][0], 1,
498  &traceFieldsAdded[9][0], 1);
499 
500  /** Evaluation of the velocity gradient in the cartesian directions
501  * Du_x: traceFieldsAdded[10]
502  * Du_y: traceFieldsAdded[11]
503  * Dv_x: traceFieldsAdded[12]
504  * Dv_y: traceFieldsAdded[13]
505  **/
506  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > Dvelocity(nDimensions);
507  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > traceDvelocity(nDimensions);
508  Array<OneD, Array<OneD, NekDouble> > velocity(nDimensions);
509 
510  for (i = 0; i < nDimensions; ++ i)
511  {
512  Dvelocity[i] = Array<OneD, Array<OneD, NekDouble> >(nDimensions);
513  traceDvelocity[i] = Array<OneD, Array<OneD, NekDouble> >(nDimensions);
514  velocity[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
515 
516  Vmath::Vdiv(nSolutionPts, uFields[i+1], 1, uFields[0], 1,
517  velocity[i], 1);
518 
519  for (j = 0; j < nDimensions; ++j)
520  {
521  Dvelocity[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
522  traceDvelocity[i][j] = Array<OneD, NekDouble>(nTracePts, 0.0);
523  }
524  }
525 
526  for (i = 0; i < nDimensions; ++i)
527  {
528  for (j = 0; j < nDimensions; ++j)
529  {
530  for (n = 0; n < nElements; n++)
531  {
532  phys_offset = pFields[0]->GetPhys_Offset(n);
533 
534  pFields[i]->GetExp(n)->PhysDeriv(
535  j, velocity[i] + phys_offset,
536  auxArray = Dvelocity[i][j] + phys_offset);
537  }
538 
539  // Extract trace
540  pFields[0]->ExtractTracePhys(Dvelocity[i][j], traceDvelocity[i][j]);
541  }
542  }
543 
544  Vmath::Vcopy(nTracePts,
545  &traceDvelocity[0][0][0], 1,
546  &traceFieldsAdded[10][0], 1);
547  Vmath::Vcopy(nTracePts,
548  &traceDvelocity[0][1][0], 1,
549  &traceFieldsAdded[11][0], 1);
550  Vmath::Vcopy(nTracePts,
551  &traceDvelocity[1][0][0], 1,
552  &traceFieldsAdded[12][0], 1);
553  Vmath::Vcopy(nTracePts,
554  &traceDvelocity[1][1][0], 1,
555  &traceFieldsAdded[13][0], 1);
556 
557 
558  /*** Evaluation of shear stresses ******************************************
559  * tau_xx -> traceFieldsAdded[14]
560  * tau_yy -> traceFieldsAdded[15]
561  * tau_xy -> traceFieldsAdded[16]
562  ***************************************************************************/
563 
564  // Stokes hypotesis
565  const NekDouble lambda = -2.0/3.0;
566 
567  // Auxiliary variables
568  Array<OneD, NekDouble > mu (nSolutionPts, 0.0);
569  Array<OneD, NekDouble > mu2 (nSolutionPts, 0.0);
570  Array<OneD, NekDouble > divVel(nSolutionPts, 0.0);
571 
572  // Variable viscosity through the Sutherland's law
573  if (m_ViscosityType == "Variable")
574  {
575  NekDouble mu_star = m_mu;
576  NekDouble T_star = m_pInf / (m_rhoInf * m_gasConstant);
577  NekDouble ratio;
578 
579  for (int i = 0; i < nSolutionPts; ++i)
580  {
581  ratio = temperature[i] / T_star;
582  mu[i] = mu_star * ratio * sqrt(ratio) *
583  (T_star + 110.0) / (temperature[i] + 110.0);
584  }
585  }
586  else
587  {
588  Vmath::Fill(nSolutionPts, m_mu, &mu[0], 1);
589  }
590 
591  // Computing diagonal terms of viscous stress tensor
592  Array<OneD, Array<OneD, NekDouble> > temp(m_spacedim);
593  Array<OneD, Array<OneD, NekDouble> > Sgg(m_spacedim);
594 
595  // mu2 = 2 * mu
596  Vmath::Smul(nSolutionPts, 2.0, &mu[0], 1, &mu2[0], 1);
597 
598  // Velocity divergence
599  Vmath::Vadd(nSolutionPts, &divVel[0], 1,
600  &Dvelocity[0][0][0], 1, &divVel[0], 1);
601  Vmath::Vadd(nSolutionPts, &divVel[0], 1,
602  &Dvelocity[1][1][0], 1, &divVel[0], 1);
603 
604  // Velocity divergence scaled by lambda * mu
605  Vmath::Smul(nSolutionPts, lambda, &divVel[0], 1, &divVel[0], 1);
606  Vmath::Vmul(nSolutionPts, &mu[0], 1, &divVel[0], 1, &divVel[0], 1);
607 
608  // Diagonal terms of viscous stress tensor (Sxx, Syy)
609  // Sjj = 2 * mu * du_j/dx_j - (2 / 3) * mu * sum_j(du_j/dx_j)
610  for (j = 0; j < m_spacedim; ++j)
611  {
612  temp[j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
613  Sgg[j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
614 
615  Vmath::Vmul(nSolutionPts, &mu2[0], 1, &Dvelocity[j][j][0], 1,
616  &temp[j][0], 1);
617 
618  Vmath::Vadd(nSolutionPts, &temp[j][0], 1, &divVel[0], 1, &Sgg[j][0], 1);
619  }
620 
621  // Extra diagonal terms of viscous stress tensor (Sxy = Syx)
622  Array<OneD, NekDouble > Sxy(nSolutionPts, 0.0);
623 
624  // Sxy = (du/dy + dv/dx)
625  Vmath::Vadd(nSolutionPts, &Dvelocity[0][1][0], 1,
626  &Dvelocity[1][0][0], 1, &Sxy[0], 1);
627 
628  // Sxy = mu * (du/dy + dv/dx)
629  Vmath::Vmul(nSolutionPts, &mu[0], 1, &Sxy[0], 1, &Sxy[0], 1);
630 
631  pFields[0]->ExtractTracePhys(Sgg[0], traceFieldsAdded[14]);
632  pFields[0]->ExtractTracePhys(Sgg[1], traceFieldsAdded[15]);
633  pFields[0]->ExtractTracePhys(Sxy, traceFieldsAdded[16]);
634 
635  /*** Evaluation of the shear stress in tangent direction *******************
636  * tau_t -> traceFieldsAdded[17]
637  ***************************************************************************/
638  Array<OneD, NekDouble > sigma_diff (nTracePts, 0.0);
639  Array<OneD, NekDouble > cosTeta (nTracePts, 0.0);
640  Array<OneD, NekDouble > sinTeta (nTracePts, 0.0);
641  Array<OneD, NekDouble > cos2Teta (nTracePts, 0.0);
642  Array<OneD, NekDouble > sin2Teta (nTracePts, 0.0);
643  Array<OneD, NekDouble > tau_t (nTracePts, 0.0);
644 
645  Array<OneD, NekDouble > tmpTeta (nTracePts, 0.0);
646 
647  // cos(teta) = nx
648  Vmath::Vcopy(nTracePts, &m_traceNormals[0][0], 1, &cosTeta[0], 1);
649 
650  // sin(teta) = ny
651  Vmath::Vcopy(nTracePts, &m_traceNormals[1][0], 1, &sinTeta[0], 1);
652 
653  // sigma_diff = sigma_x - sigma_y
654  Vmath::Vsub(nTracePts, &traceFieldsAdded[14][0], 1,
655  &traceFieldsAdded[15][0], 1, &sigma_diff[0], 1);
656 
657  // sin(2*teta)
658  Vmath::Vmul(nTracePts, &cosTeta[0], 1, &sinTeta[0], 1, &tmpTeta[0], 1);
659  Vmath::Smul(nTracePts, 2.0, &tmpTeta[0], 1, &sin2Teta[0], 1);
660 
661  // cos(2*teta)
662  Vmath::Vmul(nTracePts, &cosTeta[0], 1, &cosTeta[0], 1, &cos2Teta[0], 1);
663  Vmath::Vmul(nTracePts, &sinTeta[0], 1, &sinTeta[0], 1, &tmpTeta[0], 1);
664  Vmath::Vsub(nTracePts, &cos2Teta[0], 1, &tmpTeta[0], 1, &cos2Teta[0], 1);
665 
666  // tau_t = -0.5*sigma_diff * sin(2*teta) + tau_xy * cos(2*teta)
667  Vmath::Smul(nTracePts, -0.5, &sigma_diff[0], 1, &sigma_diff[0], 1);
668  Vmath::Vmul(nTracePts, &sigma_diff[0], 1, &sin2Teta[0], 1, &tau_t[0], 1);
669  Vmath::Vmul(nTracePts, &traceFieldsAdded[16][0], 1, &cos2Teta[0], 1,
670  &tmpTeta[0], 1);
671  Vmath::Vadd(nTracePts, &tau_t[0], 1, &tmpTeta[0], 1, &tau_t[0], 1);
672 
673  Vmath::Vcopy(nTracePts, &tau_t[0], 1, &traceFieldsAdded[17][0], 1);
674 
675  /*** Evaluation of dinamic viscosity ***************************************
676  * mu -> traceFieldsAdded[18]
677  ***************************************************************************/
678 
679  pFields[0]->ExtractTracePhys(mu, traceFieldsAdded[18]);
680 
681  /*** Evaluation of Mach number *********************************************
682  * M -> traceFieldsAdded[18]
683  ***************************************************************************/
684  NekDouble gamma = m_gamma;
685 
686  // Speed of sound
687  Array<OneD, NekDouble> soundspeed(nSolutionPts, 0.0);
688 
689  Vmath::Vdiv (nSolutionPts, pressure, 1, uFields[0], 1, soundspeed, 1);
690  Vmath::Smul (nSolutionPts, gamma, soundspeed, 1, soundspeed, 1);
691  Vmath::Vsqrt(nSolutionPts, soundspeed, 1, soundspeed, 1);
692 
693  // Mach
694  Array<OneD, NekDouble> mach(nSolutionPts, 0.0);
695 
696  for (int i = 0; i < m_spacedim; ++i)
697  {
698  Vmath::Vvtvp(nSolutionPts, uFields[i + 1], 1, uFields[i + 1], 1,
699  mach, 1, mach, 1);
700  }
701 
702  Vmath::Vdiv(nSolutionPts, mach, 1, uFields[0], 1, mach, 1);
703  Vmath::Vdiv(nSolutionPts, mach, 1, uFields[0], 1, mach, 1);
704  Vmath::Vsqrt(nSolutionPts, mach, 1, mach, 1);
705  Vmath::Vdiv(nSolutionPts, mach, 1, soundspeed, 1, mach, 1);
706 
707  pFields[0]->ExtractTracePhys(mach, traceFieldsAdded[19]);
708 
709  /**************************************************************************/
710  // Extract coordinates
711 
712  if (pFields[0]->GetBndCondExpansions().num_elements())
713  {
714  id1 = 0;
715  cnt = 0;
716  nBndRegions = pFields[0]->GetBndCondExpansions().num_elements();
717  for (b = 0; b < nBndRegions; ++b)
718  {
719  nBndEdges = pFields[0]->GetBndCondExpansions()[b]->GetExpSize();
720  for (e = 0; e < nBndEdges; ++e)
721  {
722  nBndEdgePts = pFields[0]->
723  GetBndCondExpansions()[b]->GetExp(e)->GetNumPoints(0);
724 
725  id2 = pFields[0]->GetTrace()->
726  GetPhys_Offset(pFields[0]->GetTraceMap()->
727  GetBndCondTraceToGlobalTraceMap(cnt++));
728 
729  if (pFields[0]->GetBndConditions()[b]->
730  GetUserDefined() == SpatialDomains::eWallViscous ||
731  pFields[0]->GetBndConditions()[b]->
732  GetUserDefined() == SpatialDomains::eWall)
733  {
734  Vmath::Vcopy(nBndEdgePts, &traceX[id2], 1,
735  &surfaceX[id1], 1);
736 
737  Vmath::Vcopy(nBndEdgePts, &traceY[id2], 1,
738  &surfaceY[id1], 1);
739 
740  Vmath::Vcopy(nBndEdgePts, &traceZ[id2], 1,
741  &surfaceZ[id1], 1);
742 
743  id1 += nBndEdgePts;
744  }
745  }
746  }
747  }
748 
749  // Extract fields
750  if (pFields[0]->GetBndCondExpansions().num_elements())
751  {
752  for (j = 0; j < nfields; ++j)
753  {
754  id1 = 0;
755  cnt = 0;
756  nBndRegions = pFields[j]->GetBndCondExpansions().num_elements();
757  for (b = 0; b < nBndRegions; ++b)
758  {
759  nBndEdges = pFields[j]->GetBndCondExpansions()[b]->GetExpSize();
760  for (e = 0; e < nBndEdges; ++e)
761  {
762  nBndEdgePts = pFields[j]->
763  GetBndCondExpansions()[b]->GetExp(e)->GetNumPoints(0);
764 
765  id2 = pFields[j]->GetTrace()->
766  GetPhys_Offset(pFields[j]->GetTraceMap()->
767  GetBndCondTraceToGlobalTraceMap(cnt++));
768 
769  if (pFields[j]->GetBndConditions()[b]->
770  GetUserDefined() == SpatialDomains::eWallViscous ||
771  pFields[j]->GetBndConditions()[b]->
772  GetUserDefined() == SpatialDomains::eWall)
773  {
774  Vmath::Vcopy(nBndEdgePts, &traceFields[j][id2], 1,
775  &surfaceFields[j][id1], 1);
776 
777  id1 += nBndEdgePts;
778  }
779  }
780  }
781  }
782  }
783 
784  // Extract fields added
785  if (pFields[0]->GetBndCondExpansions().num_elements())
786  {
787  for (j = 0; j < nfieldsAdded; ++j)
788  {
789  id1 = 0;
790  cnt = 0;
791  nBndRegions = pFields[0]->GetBndCondExpansions().num_elements();
792  for (b = 0; b < nBndRegions; ++b)
793  {
794  nBndEdges = pFields[0]->GetBndCondExpansions()[b]->GetExpSize();
795  for (e = 0; e < nBndEdges; ++e)
796  {
797  nBndEdgePts = pFields[0]->
798  GetBndCondExpansions()[b]->GetExp(e)->GetNumPoints(0);
799 
800  id2 = pFields[0]->GetTrace()->
801  GetPhys_Offset(pFields[0]->GetTraceMap()->
802  GetBndCondTraceToGlobalTraceMap(cnt++));
803 
804  if (pFields[0]->GetBndConditions()[b]->
805  GetUserDefined() == SpatialDomains::eWallViscous ||
806  pFields[0]->GetBndConditions()[b]->
807  GetUserDefined() == SpatialDomains::eWall)
808  {
809  Vmath::Vcopy(nBndEdgePts, &traceFieldsAdded[j][id2], 1,
810  &surfaceFieldsAdded[j][id1], 1);
811 
812  id1 += nBndEdgePts;
813  }
814  }
815  }
816  }
817  }
818  //===================================================================================================
819  //===================================================================================================
820  //===================================================================================================
821  std::string vEquation = vSession->GetSolverInfo("EQType");
822 
823  Array<OneD, MultiRegions::ExpListSharedPtr> BndExp;
824  BndExp = pFields[0]->GetBndCondExpansions();
825  StdRegions::StdExpansion1DSharedPtr bc;
826 
827  NekDouble Fxp(0.0);
828  NekDouble Fyp(0.0);
829  NekDouble Fxv(0.0);
830  NekDouble Fyv(0.0);
831  NekDouble Sref(0.0);
832 
833  int GlobalIndex(0);
834 
835  for(int i = 0; i < BndExp[0]->GetExpSize(); ++i)
836  {
837  bc = BndExp[0]->GetExp(i)->as<LocalRegions::Expansion1D> ();
838 
839  int nbc = bc->GetTotPoints();
840 
841  Array<OneD, NekDouble> nxOnBnd(nbc,0.0);
842  Array<OneD, NekDouble> nyOnBnd(nbc,0.0);
843  Array<OneD, NekDouble> txOnBnd(nbc,0.0);
844  Array<OneD, NekDouble> tyOnBnd(nbc,0.0);
845  // Array<OneD, NekDouble> CoeffAero(nbc,0.0);
846  // Array<OneD, NekDouble> tmp(nbc,0.0);
847 
848  Array<OneD, NekDouble> drag_p(nbc,0.0);
849  Array<OneD, NekDouble> lift_p(nbc,0.0);
850  Array<OneD, NekDouble> PressurOnBnd(nbc,0.0);
851 
852  Array<OneD, NekDouble> drag_v(nbc,0.0);
853  Array<OneD, NekDouble> lift_v(nbc,0.0);
854  Array<OneD, NekDouble> ShearStressOnBnd(nbc,0.0);
855 
856  Array<OneD, NekDouble> Unity(nbc,1.0);
857 
858  for(int j = 0; j < nbc; ++j)
859  {
860 
861  nxOnBnd[j] = surfaceFieldsAdded[0][GlobalIndex];
862  nyOnBnd[j] = surfaceFieldsAdded[1][GlobalIndex];
863  txOnBnd[j] = surfaceFieldsAdded[2][GlobalIndex];
864  tyOnBnd[j] = surfaceFieldsAdded[3][GlobalIndex];
865 
866  PressurOnBnd[j] = surfaceFieldsAdded[4][GlobalIndex];
867 
868  if (vEquation == "NavierStokesCFE")
869  {
870  ShearStressOnBnd[j] = surfaceFieldsAdded[17][GlobalIndex];
871  }
872 
873  // CoeffAero[j] = surfaceFields[0][GlobalIndex];
874  // tmp[j] = surfaceFields[1][GlobalIndex]*surfaceFields[1][GlobalIndex];
875  // tmp[j] = tmp[j] + surfaceFields[2][GlobalIndex]*surfaceFields[2][GlobalIndex];
876  // tmp[j] = sqrt(tmp[j]);
877  // CoeffAero[j] = CoeffAero[j]*tmp[j];
878  // CoeffAero[j] = 1.0/CoeffAero[j];
879  //
880  // PressurOnBnd[j] = CoeffAero[j]*surfaceFieldsAdded[4][GlobalIndex];
881  //
882  // cout << "CoeffAero = " << CoeffAero[j] << endl;
883 
884  GlobalIndex++;
885  }
886 
887  Vmath::Vmul(nbc,PressurOnBnd,1,nxOnBnd,1, drag_p,1);
888  Vmath::Vmul(nbc,PressurOnBnd,1,nyOnBnd,1, lift_p,1);
889 
890  // Vmath::Vmul(nbc,drag_p,1,CoeffAero,1, drag_p,1);
891  // Vmath::Vmul(nbc,lift_p,1,CoeffAero,1, lift_p,1);
892 
893  Fxp += bc->Integral(drag_p);
894  Fyp += bc->Integral(lift_p);
895 
896  if (vEquation == "NavierStokesCFE")
897  {
898  Vmath::Vmul(nbc,ShearStressOnBnd,1,txOnBnd,1, drag_v,1);
899  Vmath::Vmul(nbc,ShearStressOnBnd,1,tyOnBnd,1, lift_v,1);
900 
901  // Vmath::Vdiv(nbc,drag_v,1,CoeffAero,1, drag_v,1);
902  // Vmath::Vdiv(nbc,lift_v,1,CoeffAero,1, lift_v,1);
903 
904  Fxv += bc->Integral(drag_v);
905  Fyv += bc->Integral(lift_v);
906  }
907 
908  Sref += bc->Integral(Unity);
909 
910  }
911 
912  cout << "\n Sref = " << Sref << endl;
913  Fxp = Fxp/Sref;
914  Fyp = Fyp/Sref;
915  Fxv = Fxv/Sref;
916  Fyv = Fyv/Sref;
917  cout << " Pressure drag (Fxp) = " << Fxp << endl;
918  cout << " Pressure lift (Fyp) = " << Fyp << endl;
919  cout << " Viscous drag (Fxv) = " << Fxv << endl;
920  cout << " Viscous lift (Fyv) = " << Fyv << endl;
921  cout << "\n ==> Total drag = " << Fxp+Fxv << endl;
922  cout << " ==> Total lift = " << Fyp+Fyv << "\n" << endl;
923 
924  //===================================================================================================
925  //===================================================================================================
926  //===================================================================================================
927 
928 
929  // Print the surface coordinates and the surface solution in a .txt file
930  ofstream outfile;
931  outfile.open(fname.c_str());
932  outfile << "% x[m] " << " \t"
933  << "y[m] " << " \t"
934  << "z[m] " << " \t"
935  << "nx[] " << " \t"
936  << "ny[] " << " \t"
937  << "tx[] " << " \t"
938  << "ty[] " << " \t"
939  << "rho[kg/m^3] " << " \t"
940  << "rhou[kg/(m^2 s)] " << " \t"
941  << "rhov[kg/(m^2 s)] " << " \t"
942  << "E[Pa] " << " \t"
943  << "p[Pa] " << " \t"
944  << "T[k] " << " \t"
945  << "dT/dn[k/m] " << " \t"
946  << "dp/dT[Pa/m] " << " \t"
947  << "dp/dx[Pa/m] " << " \t"
948  << "dp/dy[Pa/m] " << " \t"
949  << "du/dx[s^-1] " << " \t"
950  << "du/dy[s^-1] " << " \t"
951  << "dv/dx[s^-1] " << " \t"
952  << "dv/dy[s^-1] " << " \t"
953  << "tau_xx[Pa] " << " \t"
954  << "tau_yy[Pa] " << " \t"
955  << "tau_xy[Pa] " << " \t"
956  << "tau_t[Pa] " << " \t"
957  << "mu[Pa s] " << " \t"
958  << "M[] " << " \t"
959  // << "Fxp " << " \t"
960  << endl;
961  for (i = 0; i < id1; ++i)
962  {
963  outfile << scientific
964  << setw (17)
965  << setprecision(16)
966  << surfaceX[i] << " \t "
967  << surfaceY[i] << " \t "
968  << surfaceZ[i] << " \t "
969  << surfaceFieldsAdded[0][i] << " \t "
970  << surfaceFieldsAdded[1][i] << " \t "
971  << surfaceFieldsAdded[2][i] << " \t "
972  << surfaceFieldsAdded[3][i] << " \t "
973  << surfaceFields[0][i] << " \t "
974  << surfaceFields[1][i] << " \t "
975  << surfaceFields[2][i] << " \t "
976  << surfaceFields[3][i] << " \t "
977  << surfaceFieldsAdded[4][i] << " \t "
978  << surfaceFieldsAdded[5][i] << " \t "
979  << surfaceFieldsAdded[6][i] << " \t "
980  << surfaceFieldsAdded[7][i] << " \t "
981  << surfaceFieldsAdded[8][i] << " \t "
982  << surfaceFieldsAdded[9][i] << " \t "
983  << surfaceFieldsAdded[10][i] << " \t "
984  << surfaceFieldsAdded[11][i] << " \t "
985  << surfaceFieldsAdded[12][i] << " \t "
986  << surfaceFieldsAdded[13][i] << " \t "
987  << surfaceFieldsAdded[14][i] << " \t "
988  << surfaceFieldsAdded[15][i] << " \t "
989  << surfaceFieldsAdded[16][i] << " \t "
990  << surfaceFieldsAdded[17][i] << " \t "
991  << surfaceFieldsAdded[18][i] << " \t "
992  << surfaceFieldsAdded[19][i] << " \t "
993  // << Fxp << " \t "
994  << endl;
995  }
996  outfile << endl << endl;
997  outfile.close();
998 
999  return 0;
1000 }
Generated on Wed Mar 4 2015 16:33:14 for Nektar++ by   doxygen 1.8.1.2