Nektar++
ExtractSurface3DCFS.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: ExtractSurface3DCFS.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// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description: Extract 2D surfaces from 3D file and output relevant quantities
32// for compressible flow solver.
33//
34///////////////////////////////////////////////////////////////////////////////
35
36#include <cstdio>
37#include <cstdlib>
38#include <iomanip>
39#include <iostream>
40#include <string>
41
42#include <MultiRegions/AssemblyMap/AssemblyMapDG.h>
43#include <MultiRegions/ExpList.h>
44
45#include <LocalRegions/Expansion.h>
46#include <LocalRegions/Expansion3D.h>
47#include <LocalRegions/MatrixKey.h>
48#include <MultiRegions/DisContField.h>
49
50#include <LibUtilities/BasicUtils/FieldIO.h>
51#include <LibUtilities/BasicUtils/NekFactory.hpp>
52#include <LibUtilities/BasicUtils/SessionReader.h>
53#include <LibUtilities/BasicUtils/SharedArray.hpp>
54
55#include <LibUtilities/Memory/NekMemoryManager.hpp>
56#include <MultiRegions/ContField.h>
57#include <SpatialDomains/MeshGraphIO.h>
58
59#include <SolverUtils/SolverUtilsDeclspec.h>
60
61using namespace std;
62using namespace Nektar;
63
64int main(int argc, char *argv[])
65{
66 string fname = std::string(argv[2]);
67 int fdot = fname.find_last_of('.');
68 if (fdot != std::string::npos)
69 {
70 string ending = fname.substr(fdot);
71
72 // If .chk or .fld we exchange the extension in the output file.
73 // For all other files (e.g. .bse) we append the extension to avoid
74 // conflicts.
75 if (ending == ".chk" || ending == ".fld")
76 {
77 fname = fname.substr(0, fdot);
78 }
79 }
80
81 fname = fname + ".txt";
82
83 int cnt;
84 int id1, id2;
85 int i, j, n, e, b;
86 Array<OneD, NekDouble> auxArray;
87
88 int nBndEdgePts, nBndEdges, nBndRegions;
89
90 if (argc < 3)
91 {
92 fprintf(stderr, "Usage: ExtractSurface3DCFS meshfile fieldFile\n");
93 fprintf(stderr,
94 "Extracts a surface from a 3D fld file"
95 "(only for CompressibleFlowSolver and purely 3D .fld files)\n");
96 exit(1);
97 }
98
99 LibUtilities::SessionReaderSharedPtr vSession =
100 LibUtilities::SessionReader::CreateInstance(3, argv);
101 SpatialDomains::MeshGraphSharedPtr graphShPt =
102 SpatialDomains::MeshGraphIO::Read(vSession);
103
104 std::string m_ViscosityType;
105
106 NekDouble m_gamma;
107 NekDouble m_pInf;
108 NekDouble m_rhoInf;
109 NekDouble m_uInf;
110 NekDouble m_vInf;
111 NekDouble m_wInf;
112 NekDouble m_gasConstant;
113 NekDouble m_Twall;
114 NekDouble m_mu;
115
116 int m_spacedim = 3;
117 int nDimensions = m_spacedim;
118 int phys_offset;
119
120 // Get gamma parameter from session file.
121 ASSERTL0(vSession->DefinesParameter("Gamma"),
122 "Compressible flow sessions must define a Gamma parameter.");
123 vSession->LoadParameter("Gamma", m_gamma, 1.4);
124
125 // Get E0 parameter from session file.
126 ASSERTL0(vSession->DefinesParameter("pInf"),
127 "Compressible flow sessions must define a pInf parameter.");
128 vSession->LoadParameter("pInf", m_pInf, 101325);
129
130 // Get rhoInf parameter from session file.
131 ASSERTL0(vSession->DefinesParameter("rhoInf"),
132 "Compressible flow sessions must define a rhoInf parameter.");
133 vSession->LoadParameter("rhoInf", m_rhoInf, 1.225);
134
135 // Get uInf parameter from session file.
136 ASSERTL0(vSession->DefinesParameter("uInf"),
137 "Compressible flow sessions must define a uInf parameter.");
138 vSession->LoadParameter("uInf", m_uInf, 0.1);
139
140 // Get vInf parameter from session file.
141 if (m_spacedim == 2 || m_spacedim == 3)
142 {
143 ASSERTL0(vSession->DefinesParameter("vInf"),
144 "Compressible flow sessions must define a vInf parameter"
145 "for 2D/3D problems.");
146 vSession->LoadParameter("vInf", m_vInf, 0.0);
147 }
148
149 // Get wInf parameter from session file.
150 if (m_spacedim == 3)
151 {
152 ASSERTL0(vSession->DefinesParameter("wInf"),
153 "Compressible flow sessions must define a wInf parameter"
154 "for 3D problems.");
155 vSession->LoadParameter("wInf", m_wInf, 0.0);
156 }
157
158 vSession->LoadParameter("GasConstant", m_gasConstant, 287.058);
159 vSession->LoadParameter("Twall", m_Twall, 300.15);
160 vSession->LoadSolverInfo("ViscosityType", m_ViscosityType, "Constant");
161 vSession->LoadParameter("mu", m_mu, 1.78e-05);
162
163 //--------------------------------------------------------------------------
164 // Import field file
165 string fieldFile(argv[2]);
166 vector<LibUtilities::FieldDefinitionsSharedPtr> fieldDef;
167 vector<vector<NekDouble>> fieldData;
168
169 LibUtilities::Import(fieldFile, fieldDef, fieldData);
170 //--------------------------------------------------------------------------
171
172 //--------------------------------------------------------------------------
173 // Set up Expansion information
174 vector<vector<LibUtilities::PointsType>> pointsType;
175 for (i = 0; i < fieldDef.size(); ++i)
176 {
177 vector<LibUtilities::PointsType> ptype;
178 for (j = 0; j < 3; ++j)
179 {
180 ptype.push_back(LibUtilities::ePolyEvenlySpaced);
181 }
182 pointsType.push_back(ptype);
183 }
184 graphShPt->SetExpansionInfo(fieldDef, pointsType);
185
186 //--------------------------------------------------------------------------
187
188 //--------------------------------------------------------------------------
189 // Define Expansion
190 int nfields = fieldDef[0]->m_fields.size();
191 Array<OneD, MultiRegions::ExpListSharedPtr> Exp(nfields);
192 Array<OneD, MultiRegions::ExpListSharedPtr> pFields(nfields);
193
194 for (i = 0; i < pFields.size(); i++)
195 {
196 pFields[i] =
197 MemoryManager<MultiRegions ::DisContField>::AllocateSharedPtr(
198 vSession, graphShPt, vSession->GetVariable(i));
199 }
200
201 MultiRegions::ExpListSharedPtr Exp3D;
202 Exp3D = MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(vSession,
203 graphShPt);
204
205 Exp[0] = Exp3D;
206
207 for (i = 1; i < nfields; ++i)
208 {
209 Exp[i] =
210 MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(*Exp3D);
211 }
212
213 // Count of the point on the surface
214 int nSurfacePts = 0;
215 if (pFields[0]->GetBndCondExpansions().size())
216 {
217 nSurfacePts = 0;
218 cnt = 0;
219 nBndRegions = pFields[0]->GetBndCondExpansions().size();
220 for (b = 0; b < nBndRegions; ++b)
221 {
222 nBndEdges = pFields[0]->GetBndCondExpansions()[b]->GetExpSize();
223 for (e = 0; e < nBndEdges; ++e)
224 {
225 nBndEdgePts = pFields[0]
226 ->GetBndCondExpansions()[b]
227 ->GetExp(e)
228 ->GetTotPoints();
229
230 if (pFields[0]->GetBndConditions()[b]->GetUserDefined() ==
231 "WallViscous" ||
232 pFields[0]->GetBndConditions()[b]->GetUserDefined() ==
233 "WallAdiabatic" ||
234 pFields[0]->GetBndConditions()[b]->GetUserDefined() ==
235 "Wall")
236 {
237 nSurfacePts += nBndEdgePts;
238 }
239 }
240 }
241 }
242
243 int nSolutionPts = pFields[0]->GetNpoints();
244 int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
245 int nElements = pFields[0]->GetExpSize();
246
247 Array<OneD, NekDouble> tmp(nSolutionPts, 0.0);
248
249 Array<OneD, NekDouble> x(nSolutionPts);
250 Array<OneD, NekDouble> y(nSolutionPts);
251 Array<OneD, NekDouble> z(nSolutionPts);
252
253 Array<OneD, NekDouble> traceX(nTracePts);
254 Array<OneD, NekDouble> traceY(nTracePts);
255 Array<OneD, NekDouble> traceZ(nTracePts);
256
257 Array<OneD, NekDouble> surfaceX(nSurfacePts);
258 Array<OneD, NekDouble> surfaceY(nSurfacePts);
259 Array<OneD, NekDouble> surfaceZ(nSurfacePts);
260
261 pFields[0]->GetCoords(x, y, z);
262
263 pFields[0]->ExtractTracePhys(x, traceX);
264 pFields[0]->ExtractTracePhys(y, traceY);
265 pFields[0]->ExtractTracePhys(z, traceZ);
266 //--------------------------------------------------------------------------
267
268 //--------------------------------------------------------------------------
269 // Copy data from field file
270 Array<OneD, Array<OneD, NekDouble>> uFields(nfields);
271 Array<OneD, Array<OneD, NekDouble>> traceFields(nfields);
272 Array<OneD, Array<OneD, NekDouble>> surfaceFields(nfields);
273
274 // Extract the physical values of the solution at the boundaries
275 for (j = 0; j < nfields; ++j)
276 {
277 uFields[j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
278 traceFields[j] = Array<OneD, NekDouble>(nTracePts, 0.0);
279 surfaceFields[j] = Array<OneD, NekDouble>(nSurfacePts, 0.0);
280
281 for (i = 0; i < fieldData.size(); ++i)
282 {
283 Exp[j]->ExtractDataToCoeffs(fieldDef[i], fieldData[i],
284 fieldDef[i]->m_fields[j],
285 Exp[j]->UpdateCoeffs());
286 }
287 Exp[j]->BwdTrans(Exp[j]->GetCoeffs(), Exp[j]->UpdatePhys());
288 Vmath::Vcopy(nSolutionPts, Exp[j]->GetPhys(), 1, uFields[j], 1);
289 pFields[0]->ExtractTracePhys(uFields[j], traceFields[j]);
290 }
291
292 // Fields to add in the output file
293 int nfieldsAdded = 34;
294 Array<OneD, Array<OneD, NekDouble>> traceFieldsAdded(nfieldsAdded);
295 Array<OneD, Array<OneD, NekDouble>> surfaceFieldsAdded(nfieldsAdded);
296
297 for (j = 0; j < nfieldsAdded; ++j)
298 {
299 traceFieldsAdded[j] = Array<OneD, NekDouble>(nTracePts, 0.0);
300 surfaceFieldsAdded[j] = Array<OneD, NekDouble>(nSurfacePts, 0.0);
301 }
302
303 /******** Evaluation of normals and tangents on the trace *****************
304 * nx -> traceFieldsAdded[0];
305 * ny -> traceFieldsAdded[1];
306 * nz -> traceFieldsAdded[2];
307 * bx -> traceFieldsAdded[3];
308 * by -> traceFieldsAdded[4];
309 * bz -> traceFieldsAdded[5];
310 * tx -> traceFieldsAdded[6];
311 * ty -> traceFieldsAdded[7];
312 * tz -> traceFieldsAdded[8];
313
314 ***************************************************************************/
315
316 Array<OneD, Array<OneD, NekDouble>> m_traceNormals(nDimensions);
317 for (i = 0; i < nDimensions; ++i)
318 {
319 m_traceNormals[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
320 }
321 pFields[0]->GetTrace()->GetNormals(m_traceNormals);
322
323 Array<OneD, Array<OneD, NekDouble>> m_traceTangents(nDimensions);
324 Array<OneD, Array<OneD, NekDouble>> m_traceBinormals(nDimensions);
325 Array<OneD, Array<OneD, NekDouble>> h(nDimensions);
326 Array<OneD, NekDouble> tmpNorm(nTracePts, 1.0);
327 Array<OneD, NekDouble> NormH(nTracePts, 0.0);
328 Array<OneD, NekDouble> tmpTrace(nTracePts, 0.0);
329
330 for (i = 0; i < nDimensions; ++i)
331 {
332 m_traceTangents[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
333 m_traceBinormals[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
334 h[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
335 }
336
337 // Normals
338
339 // nx
340 Vmath::Vcopy(nTracePts, &m_traceNormals[0][0], 1, &traceFieldsAdded[0][0],
341 1);
342
343 // ny
344 Vmath::Vcopy(nTracePts, &m_traceNormals[1][0], 1, &traceFieldsAdded[1][0],
345 1);
346
347 // nz
348 Vmath::Vcopy(nTracePts, &m_traceNormals[2][0], 1, &traceFieldsAdded[2][0],
349 1);
350
351 // Tangents and Binormals
352 // h1
353 Vmath::Vadd(nTracePts, &m_traceNormals[0][0], 1, &tmpNorm[0], 1, &h[0][0],
354 1);
355 // h2
356 Vmath::Vcopy(nTracePts, &m_traceNormals[1][0], 1, &h[1][0], 1);
357 // h3
358 Vmath::Vcopy(nTracePts, &m_traceNormals[2][0], 1, &h[2][0], 1);
359
360 // Norm of h
361 for (i = 0; i < m_spacedim; i++)
362 {
363 Vmath::Vvtvp(nTracePts, &h[i][0], 1, &h[i][0], 1, &NormH[0], 1,
364 &NormH[0], 1);
365 }
366
367 // b1
368 Vmath::Vmul(nTracePts, &h[0][0], 1, &h[1][0], 1, &tmpTrace[0], 1);
369
370 Vmath::Vdiv(nTracePts, &tmpTrace[0], 1, &NormH[0], 1, &tmpTrace[0], 1);
371
372 Vmath::Smul(nTracePts, -2.0, &tmpTrace[0], 1, &m_traceBinormals[0][0], 1);
373
374 Vmath::Vcopy(nTracePts, &m_traceBinormals[0][0], 1, &traceFieldsAdded[3][0],
375 1);
376
377 // b2
378 Vmath::Vmul(nTracePts, &h[1][0], 1, &h[1][0], 1, &tmpTrace[0], 1);
379
380 Vmath::Vdiv(nTracePts, &tmpTrace[0], 1, &NormH[0], 1, &tmpTrace[0], 1);
381
382 Vmath::Smul(nTracePts, -2.0, &tmpTrace[0], 1, &tmpTrace[0], 1);
383
384 Vmath::Vadd(nTracePts, &tmpTrace[0], 1, &tmpNorm[0], 1,
385 &m_traceBinormals[1][0], 1);
386
387 Vmath::Vcopy(nTracePts, &m_traceBinormals[1][0], 1, &traceFieldsAdded[4][0],
388 1);
389
390 // b3
391 Vmath::Vmul(nTracePts, &h[1][0], 1, &h[2][0], 1, &tmpTrace[0], 1);
392
393 Vmath::Vdiv(nTracePts, &tmpTrace[0], 1, &NormH[0], 1, &tmpTrace[0], 1);
394
395 Vmath::Smul(nTracePts, -2.0, &tmpTrace[0], 1, &m_traceBinormals[2][0], 1);
396
397 Vmath::Vcopy(nTracePts, &m_traceBinormals[2][0], 1, &traceFieldsAdded[5][0],
398 1);
399
400 // t1
401 Vmath::Vmul(nTracePts, &h[0][0], 1, &h[2][0], 1, &tmpTrace[0], 1);
402
403 Vmath::Vdiv(nTracePts, &tmpTrace[0], 1, &NormH[0], 1, &tmpTrace[0], 1);
404
405 Vmath::Smul(nTracePts, -2.0, &tmpTrace[0], 1, &m_traceTangents[0][0], 1);
406
407 Vmath::Vcopy(nTracePts, &m_traceTangents[0][0], 1, &traceFieldsAdded[6][0],
408 1);
409
410 // t2
411 Vmath::Vcopy(nTracePts, &m_traceBinormals[2][0], 1, &m_traceTangents[1][0],
412 1);
413
414 Vmath::Vcopy(nTracePts, &m_traceTangents[1][0], 1, &traceFieldsAdded[7][0],
415 1);
416
417 // t3
418 Vmath::Vmul(nTracePts, &h[2][0], 1, &h[2][0], 1, &tmpTrace[0], 1);
419
420 Vmath::Vdiv(nTracePts, &tmpTrace[0], 1, &NormH[0], 1, &tmpTrace[0], 1);
421
422 Vmath::Smul(nTracePts, -2.0, &tmpTrace[0], 1, &tmpTrace[0], 1);
423
424 Vmath::Vadd(nTracePts, &tmpTrace[0], 1, &tmpNorm[0], 1,
425 &m_traceTangents[2][0], 1);
426
427 Vmath::Vcopy(nTracePts, &m_traceTangents[2][0], 1, &traceFieldsAdded[8][0],
428 1);
429
430 /******** Evaluation of the pressure ***************************************
431 * P = (E-1/2.*rho.*((rhou./rho).^2+(rhov./rho).^2))*(gamma - 1);
432 * P -> traceFieldsAdded[9];
433 ***************************************************************************/
434
435 Array<OneD, NekDouble> pressure(nSolutionPts, 0.0);
436 NekDouble gammaMinusOne = m_gamma - 1.0;
437
438 for (i = 0; i < m_spacedim; i++)
439 {
440 Vmath::Vmul(nSolutionPts, &uFields[i + 1][0], 1, &uFields[i + 1][0], 1,
441 &tmp[0], 1);
442
443 Vmath::Smul(nSolutionPts, 0.5, &tmp[0], 1, &tmp[0], 1);
444
445 Vmath::Vadd(nSolutionPts, &pressure[0], 1, &tmp[0], 1, &pressure[0], 1);
446 }
447
448 Vmath::Vdiv(nSolutionPts, &pressure[0], 1, &uFields[0][0], 1, &pressure[0],
449 1);
450
451 Vmath::Vsub(nSolutionPts, &uFields[nfields - 1][0], 1, &pressure[0], 1,
452 &pressure[0], 1);
453
454 Vmath::Smul(nSolutionPts, gammaMinusOne, &pressure[0], 1, &pressure[0], 1);
455
456 // Extract trace
457 pFields[0]->ExtractTracePhys(pressure, traceFieldsAdded[9]);
458
459 /******** Evaluation of the temperature ************************************
460 * T = P/(R*rho);
461 * T -> traceFieldsAdded[10];
462 ***************************************************************************/
463
464 Array<OneD, NekDouble> temperature(nSolutionPts, 0.0);
465
466 Vmath::Vdiv(nSolutionPts, &pressure[0], 1, &uFields[0][0], 1,
467 &temperature[0], 1);
468
469 NekDouble GasConstantInv = 1.0 / m_gasConstant;
470 Vmath::Smul(nSolutionPts, GasConstantInv, &temperature[0], 1,
471 &temperature[0], 1);
472
473 // Extract trace
474 pFields[0]->ExtractTracePhys(temperature, traceFieldsAdded[10]);
475
476 /*** Evaluation of the temperature gradient in the normal direction ********
477 * DT_n -> traceFieldsAdded[11]
478 ***************************************************************************/
479
480 Array<OneD, Array<OneD, NekDouble>> Dtemperature(nDimensions);
481 Array<OneD, Array<OneD, NekDouble>> traceDtemperature(nDimensions);
482
483 for (i = 0; i < nDimensions; ++i)
484 {
485 Dtemperature[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
486 traceDtemperature[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
487 }
488
489 for (i = 0; i < nDimensions; ++i)
490 {
491 for (n = 0; n < nElements; n++)
492 {
493 phys_offset = pFields[0]->GetPhys_Offset(n);
494
495 pFields[i]->GetExp(n)->PhysDeriv(i, temperature + phys_offset,
496 auxArray =
497 Dtemperature[i] + phys_offset);
498 }
499 // Extract trace
500 pFields[0]->ExtractTracePhys(Dtemperature[i], traceDtemperature[i]);
501 }
502
503 for (i = 0; i < nDimensions; ++i)
504 {
505 Vmath::Vmul(nTracePts, &m_traceNormals[i][0], 1,
506 &traceDtemperature[i][0], 1, &tmp[0], 1);
507
508 Vmath::Vadd(nTracePts, &traceFieldsAdded[11][0], 1, &tmp[0], 1,
509 &traceFieldsAdded[11][0], 1);
510 }
511
512 /*** Evaluation of the pressure gradient ***********************************
513 * DP_t -> traceFieldsAdded[12] tangent direction
514 * DP_b -> traceFieldsAdded[13] binormal direction
515 * DP_x -> traceFieldsAdded[14]
516 * DP_y -> traceFieldsAdded[15]
517 * DP_z -> traceFieldsAdded[16]
518 ***************************************************************************/
519
520 Array<OneD, Array<OneD, NekDouble>> Dpressure(nDimensions);
521 Array<OneD, Array<OneD, NekDouble>> traceDpressure(nDimensions);
522
523 for (i = 0; i < nDimensions; ++i)
524 {
525 Dpressure[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
526 traceDpressure[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
527 }
528
529 for (i = 0; i < nDimensions; ++i)
530 {
531 for (n = 0; n < nElements; n++)
532 {
533 phys_offset = pFields[0]->GetPhys_Offset(n);
534
535 pFields[i]->GetExp(n)->PhysDeriv(i, pressure + phys_offset,
536 auxArray =
537 Dpressure[i] + phys_offset);
538 }
539 // Extract trace
540 pFields[0]->ExtractTracePhys(Dpressure[i], traceDpressure[i]);
541 }
542
543 // Dp_t
544 for (i = 0; i < nDimensions; ++i)
545 {
546 Vmath::Vmul(nTracePts, &m_traceTangents[i][0], 1, &traceDpressure[i][0],
547 1, &tmp[0], 1);
548
549 Vmath::Vadd(nTracePts, &traceFieldsAdded[12][0], 1, &tmp[0], 1,
550 &traceFieldsAdded[12][0], 1);
551 }
552
553 // Dp_b
554 for (i = 0; i < nDimensions; ++i)
555 {
556 Vmath::Vmul(nTracePts, &m_traceBinormals[i][0], 1,
557 &traceDpressure[i][0], 1, &tmp[0], 1);
558
559 Vmath::Vadd(nTracePts, &traceFieldsAdded[13][0], 1, &tmp[0], 1,
560 &traceFieldsAdded[13][0], 1);
561 }
562
563 // Dp_x
564 Vmath::Vcopy(nTracePts, &traceDpressure[0][0], 1, &traceFieldsAdded[14][0],
565 1);
566
567 // Dp_y
568 Vmath::Vcopy(nTracePts, &traceDpressure[1][0], 1, &traceFieldsAdded[15][0],
569 1);
570
571 // Dp_z
572 Vmath::Vcopy(nTracePts, &traceDpressure[2][0], 1, &traceFieldsAdded[16][0],
573 1);
574
575 /** Evaluation of the velocity gradient in the cartesian directions
576 * Du_x: traceFieldsAdded[17]
577 * Du_y: traceFieldsAdded[18]
578 * Du_z: traceFieldsAdded[19]
579 * Dv_x: traceFieldsAdded[20]
580 * Dv_y: traceFieldsAdded[21]
581 * Dv_z: traceFieldsAdded[22]
582 * Dw_x: traceFieldsAdded[23]
583 * Dw_y: traceFieldsAdded[24]
584 * Dw_z: traceFieldsAdded[25]
585 **/
586
587 Array<OneD, Array<OneD, Array<OneD, NekDouble>>> Dvelocity(nDimensions);
588 Array<OneD, Array<OneD, Array<OneD, NekDouble>>> traceDvelocity(
589 nDimensions);
590 Array<OneD, Array<OneD, NekDouble>> velocity(nDimensions);
591
592 for (i = 0; i < nDimensions; ++i)
593 {
594 Dvelocity[i] = Array<OneD, Array<OneD, NekDouble>>(nDimensions);
595 traceDvelocity[i] = Array<OneD, Array<OneD, NekDouble>>(nDimensions);
596 velocity[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
597
598 Vmath::Vdiv(nSolutionPts, uFields[i + 1], 1, uFields[0], 1, velocity[i],
599 1);
600
601 for (j = 0; j < nDimensions; ++j)
602 {
603 Dvelocity[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
604 traceDvelocity[i][j] = Array<OneD, NekDouble>(nTracePts, 0.0);
605 }
606 }
607
608 for (i = 0; i < nDimensions; ++i)
609 {
610 for (j = 0; j < nDimensions; ++j)
611 {
612 for (n = 0; n < nElements; n++)
613 {
614 phys_offset = pFields[0]->GetPhys_Offset(n);
615
616 pFields[i]->GetExp(n)->PhysDeriv(j, velocity[i] + phys_offset,
617 auxArray = Dvelocity[i][j] +
618 phys_offset);
619 }
620
621 // Extract trace
622 pFields[0]->ExtractTracePhys(Dvelocity[i][j], traceDvelocity[i][j]);
623 }
624 }
625
626 Vmath::Vcopy(nTracePts, &traceDvelocity[0][0][0], 1,
627 &traceFieldsAdded[17][0], 1);
628 Vmath::Vcopy(nTracePts, &traceDvelocity[0][1][0], 1,
629 &traceFieldsAdded[18][0], 1);
630 Vmath::Vcopy(nTracePts, &traceDvelocity[0][2][0], 1,
631 &traceFieldsAdded[19][0], 1);
632 Vmath::Vcopy(nTracePts, &traceDvelocity[1][0][0], 1,
633 &traceFieldsAdded[20][0], 1);
634 Vmath::Vcopy(nTracePts, &traceDvelocity[1][1][0], 1,
635 &traceFieldsAdded[21][0], 1);
636 Vmath::Vcopy(nTracePts, &traceDvelocity[1][2][0], 1,
637 &traceFieldsAdded[22][0], 1);
638 Vmath::Vcopy(nTracePts, &traceDvelocity[2][0][0], 1,
639 &traceFieldsAdded[23][0], 1);
640 Vmath::Vcopy(nTracePts, &traceDvelocity[2][1][0], 1,
641 &traceFieldsAdded[24][0], 1);
642 Vmath::Vcopy(nTracePts, &traceDvelocity[2][2][0], 1,
643 &traceFieldsAdded[25][0], 1);
644
645 /*** Evaluation of shear stresses ******************************************
646 * tau_xx -> traceFieldsAdded[26]
647 * tau_yy -> traceFieldsAdded[27]
648 * tau_zz -> traceFieldsAdded[28]
649 * tau_xy -> traceFieldsAdded[29]
650 * tau_xz -> traceFieldsAdded[30]
651 * tau_yz -> traceFieldsAdded[31]
652 ***************************************************************************/
653
654 // Stokes hypotesis
655 const NekDouble lambda = -2.0 / 3.0;
656
657 // Auxiliary variables
658 Array<OneD, NekDouble> mu(nSolutionPts, 0.0);
659 Array<OneD, NekDouble> mu2(nSolutionPts, 0.0);
660 Array<OneD, NekDouble> divVel(nSolutionPts, 0.0);
661
662 // Variable viscosity through the Sutherland's law
663 if (m_ViscosityType == "Variable")
664 {
665 NekDouble mu_star = m_mu;
666 NekDouble T_star = m_pInf / (m_rhoInf * m_gasConstant);
667 NekDouble ratio;
668
669 for (int i = 0; i < nSolutionPts; ++i)
670 {
671 ratio = temperature[i] / T_star;
672 mu[i] = mu_star * ratio * sqrt(ratio) * (T_star + 110.0) /
673 (temperature[i] + 110.0);
674 }
675 }
676 else
677 {
678 Vmath::Fill(nSolutionPts, m_mu, &mu[0], 1);
679 }
680
681 // Computing diagonal terms of viscous stress tensor
682 Array<OneD, Array<OneD, NekDouble>> temp(m_spacedim);
683 Array<OneD, Array<OneD, NekDouble>> Sgg(m_spacedim);
684
685 // mu2 = 2 * mu
686 Vmath::Smul(nSolutionPts, 2.0, &mu[0], 1, &mu2[0], 1);
687
688 // Velocity divergence
689 Vmath::Vadd(nSolutionPts, &divVel[0], 1, &Dvelocity[0][0][0], 1, &divVel[0],
690 1);
691 Vmath::Vadd(nSolutionPts, &divVel[0], 1, &Dvelocity[1][1][0], 1, &divVel[0],
692 1);
693
694 // Velocity divergence scaled by lambda * mu
695 Vmath::Smul(nSolutionPts, lambda, &divVel[0], 1, &divVel[0], 1);
696 Vmath::Vmul(nSolutionPts, &mu[0], 1, &divVel[0], 1, &divVel[0], 1);
697
698 // Diagonal terms of viscous stress tensor (Sxx, Syy)
699 // Sjj = 2 * mu * du_j/dx_j - (2 / 3) * mu * sum_j(du_j/dx_j)
700 for (j = 0; j < m_spacedim; ++j)
701 {
702 temp[j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
703 Sgg[j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
704
705 Vmath::Vmul(nSolutionPts, &mu2[0], 1, &Dvelocity[j][j][0], 1,
706 &temp[j][0], 1);
707
708 Vmath::Vadd(nSolutionPts, &temp[j][0], 1, &divVel[0], 1, &Sgg[j][0], 1);
709 }
710
711 // Extra diagonal terms of viscous stress tensor
712 Array<OneD, NekDouble> Sxy(nSolutionPts, 0.0);
713 Array<OneD, NekDouble> Sxz(nSolutionPts, 0.0);
714 Array<OneD, NekDouble> Syz(nSolutionPts, 0.0);
715
716 // Sxy = (du/dy + dv/dx)
717 Vmath::Vadd(nSolutionPts, &Dvelocity[0][1][0], 1, &Dvelocity[1][0][0], 1,
718 &Sxy[0], 1);
719
720 // Sxz = (du/dz + dw/dx)
721 Vmath::Vadd(nSolutionPts, &Dvelocity[0][2][0], 1, &Dvelocity[2][0][0], 1,
722 &Sxz[0], 1);
723
724 // Syz = (dv/dz + dw/dy)
725 Vmath::Vadd(nSolutionPts, &Dvelocity[1][2][0], 1, &Dvelocity[2][1][0], 1,
726 &Syz[0], 1);
727
728 // Sxy = mu * (du/dy + dv/dx)
729 Vmath::Vmul(nSolutionPts, &mu[0], 1, &Sxy[0], 1, &Sxy[0], 1);
730
731 // Sxz = mu * (du/dy + dv/dx)
732 Vmath::Vmul(nSolutionPts, &mu[0], 1, &Sxz[0], 1, &Sxz[0], 1);
733
734 // Syz = mu * (du/dy + dv/dx)
735 Vmath::Vmul(nSolutionPts, &mu[0], 1, &Syz[0], 1, &Syz[0], 1);
736
737 pFields[0]->ExtractTracePhys(Sgg[0], traceFieldsAdded[26]);
738 pFields[0]->ExtractTracePhys(Sgg[1], traceFieldsAdded[27]);
739 pFields[0]->ExtractTracePhys(Sgg[2], traceFieldsAdded[28]);
740 pFields[0]->ExtractTracePhys(Sxy, traceFieldsAdded[29]);
741 pFields[0]->ExtractTracePhys(Sxz, traceFieldsAdded[30]);
742 pFields[0]->ExtractTracePhys(Syz, traceFieldsAdded[31]);
743
744 /*** Evaluation of dinamic viscosity ***************************************
745 * mu -> traceFieldsAdded[32]
746 ***************************************************************************/
747
748 pFields[0]->ExtractTracePhys(mu, traceFieldsAdded[32]);
749
750 /*** Evaluation of Mach number *********************************************
751 * M -> traceFieldsAdded[33]
752 ***************************************************************************/
753 NekDouble gamma = m_gamma;
754
755 // Speed of sound
756 Array<OneD, NekDouble> soundspeed(nSolutionPts, 0.0);
757
758 Vmath::Vdiv(nSolutionPts, pressure, 1, uFields[0], 1, soundspeed, 1);
759 Vmath::Smul(nSolutionPts, gamma, soundspeed, 1, soundspeed, 1);
760 Vmath::Vsqrt(nSolutionPts, soundspeed, 1, soundspeed, 1);
761
762 // Mach
763 Array<OneD, NekDouble> mach(nSolutionPts, 0.0);
764
765 for (int i = 0; i < m_spacedim; ++i)
766 {
767 Vmath::Vvtvp(nSolutionPts, uFields[i + 1], 1, uFields[i + 1], 1, mach,
768 1, mach, 1);
769 }
770
771 Vmath::Vdiv(nSolutionPts, mach, 1, uFields[0], 1, mach, 1);
772 Vmath::Vdiv(nSolutionPts, mach, 1, uFields[0], 1, mach, 1);
773 Vmath::Vsqrt(nSolutionPts, mach, 1, mach, 1);
774 Vmath::Vdiv(nSolutionPts, mach, 1, soundspeed, 1, mach, 1);
775
776 pFields[0]->ExtractTracePhys(mach, traceFieldsAdded[33]);
777
778 /**************************************************************************/
779 // Extract coordinates
780
781 if (pFields[0]->GetBndCondExpansions().size())
782 {
783 id1 = 0;
784 cnt = 0;
785 nBndRegions = pFields[0]->GetBndCondExpansions().size();
786 for (b = 0; b < nBndRegions; ++b)
787 {
788 nBndEdges = pFields[0]->GetBndCondExpansions()[b]->GetExpSize();
789 for (e = 0; e < nBndEdges; ++e)
790 {
791 nBndEdgePts = pFields[0]
792 ->GetBndCondExpansions()[b]
793 ->GetExp(e)
794 ->GetTotPoints();
795
796 id2 = pFields[0]->GetTrace()->GetPhys_Offset(
797 pFields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(
798 cnt++));
799
800 if (pFields[0]->GetBndConditions()[b]->GetUserDefined() ==
801 "WallViscous" ||
802 pFields[0]->GetBndConditions()[b]->GetUserDefined() ==
803 "WallAdiabatic" ||
804 pFields[0]->GetBndConditions()[b]->GetUserDefined() ==
805 "Wall")
806 {
807
808 Vmath::Vcopy(nBndEdgePts, &traceX[id2], 1, &surfaceX[id1],
809 1);
810
811 Vmath::Vcopy(nBndEdgePts, &traceY[id2], 1, &surfaceY[id1],
812 1);
813
814 Vmath::Vcopy(nBndEdgePts, &traceZ[id2], 1, &surfaceZ[id1],
815 1);
816
817 id1 += nBndEdgePts;
818 }
819 }
820 }
821 }
822
823 // Extract fields
824 if (pFields[0]->GetBndCondExpansions().size())
825 {
826
827 for (j = 0; j < nfields; ++j)
828 {
829 cout << "field " << j << endl;
830
831 id1 = 0;
832 cnt = 0;
833 nBndRegions = pFields[j]->GetBndCondExpansions().size();
834 for (b = 0; b < nBndRegions; ++b)
835 {
836 nBndEdges = pFields[j]->GetBndCondExpansions()[b]->GetExpSize();
837 for (e = 0; e < nBndEdges; ++e)
838 {
839 nBndEdgePts = pFields[j]
840 ->GetBndCondExpansions()[b]
841 ->GetExp(e)
842 ->GetTotPoints();
843
844 id2 = pFields[j]->GetTrace()->GetPhys_Offset(
845 pFields[j]->GetTraceMap()->GetBndCondIDToGlobalTraceID(
846 cnt++));
847
848 if (pFields[j]->GetBndConditions()[b]->GetUserDefined() ==
849 "WallViscous" ||
850 pFields[j]->GetBndConditions()[b]->GetUserDefined() ==
851 "WallAdiabatic" ||
852 pFields[j]->GetBndConditions()[b]->GetUserDefined() ==
853 "Wall")
854 {
855 Vmath::Vcopy(nBndEdgePts, &traceFields[j][id2], 1,
856 &surfaceFields[j][id1], 1);
857
858 id1 += nBndEdgePts;
859 }
860 }
861 }
862 }
863 }
864
865 // Extract fields added
866 if (pFields[0]->GetBndCondExpansions().size())
867 {
868 for (j = 0; j < nfieldsAdded; ++j)
869 {
870 cout << "field added " << j << endl;
871
872 id1 = 0;
873 cnt = 0;
874 nBndRegions = pFields[0]->GetBndCondExpansions().size();
875 for (b = 0; b < nBndRegions; ++b)
876 {
877 nBndEdges = pFields[0]->GetBndCondExpansions()[b]->GetExpSize();
878 for (e = 0; e < nBndEdges; ++e)
879 {
880 nBndEdgePts = pFields[0]
881 ->GetBndCondExpansions()[b]
882 ->GetExp(e)
883 ->GetTotPoints();
884
885 id2 = pFields[0]->GetTrace()->GetPhys_Offset(
886 pFields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(
887 cnt++));
888
889 if (pFields[0]->GetBndConditions()[b]->GetUserDefined() ==
890 "WallViscous" ||
891 pFields[0]->GetBndConditions()[b]->GetUserDefined() ==
892 "WallAdiabatic" ||
893 pFields[0]->GetBndConditions()[b]->GetUserDefined() ==
894 "Wall")
895 {
896 Vmath::Vcopy(nBndEdgePts, &traceFieldsAdded[j][id2], 1,
897 &surfaceFieldsAdded[j][id1], 1);
898
899 id1 += nBndEdgePts;
900 }
901 }
902 }
903 }
904 }
905
906 //==========================================================================
907 //==========================================================================
908 //==========================================================================
909
910 // Print the surface coordinates and the surface solution in a .txt file
911 ofstream outfile;
912 outfile.open(fname.c_str());
913 outfile << "% x[m] "
914 << " \t"
915 << "y[m] "
916 << " \t"
917 << "z[m] "
918 << " \t"
919 << "nx[] "
920 << " \t"
921 << "ny[] "
922 << " \t"
923 << "nz[] "
924 << " \t"
925 << "bx[] "
926 << " \t"
927 << "by[] "
928 << " \t"
929 << "bz[] "
930 << " \t"
931 << "tx[] "
932 << " \t"
933 << "ty[] "
934 << " \t"
935 << "tz[] "
936 << " \t"
937 << "rho[kg/m^3] "
938 << " \t"
939 << "rhou[kg/(m^2 s)] "
940 << " \t"
941 << "rhov[kg/(m^2 s)] "
942 << " \t"
943 << "rhow[kg/(m^2 s)] "
944 << " \t"
945 << "E[Pa] "
946 << " \t"
947 << "p[Pa] "
948 << " \t"
949 << "T[k] "
950 << " \t"
951 << "dT/dn[k/m] "
952 << " \t"
953 << "dp/dT[Pa/m] "
954 << " \t"
955 << "dp/dB[Pa/m] "
956 << " \t"
957 << "dp/dx[Pa/m] "
958 << " \t"
959 << "dp/dy[Pa/m] "
960 << " \t"
961 << "dp/dz[Pa/m] "
962 << " \t"
963 << "du/dx[s^-1] "
964 << " \t"
965 << "du/dy[s^-1] "
966 << " \t"
967 << "du/dz[s^-1] "
968 << " \t"
969 << "dv/dx[s^-1] "
970 << " \t"
971 << "dv/dy[s^-1] "
972 << " \t"
973 << "dv/dz[s^-1] "
974 << " \t"
975 << "dw/dx[s^-1] "
976 << " \t"
977 << "dw/dy[s^-1] "
978 << " \t"
979 << "dw/dz[s^-1] "
980 << " \t"
981 << "tau_xx[Pa] "
982 << " \t"
983 << "tau_yy[Pa] "
984 << " \t"
985 << "tau_zz[Pa] "
986 << " \t"
987 << "tau_xy[Pa] "
988 << " \t"
989 << "tau_xz[Pa] "
990 << " \t"
991 << "tau_yz[Pa] "
992 << " \t"
993 << "mu[Pa s] "
994 << " \t"
995 << "M[] "
996 << " \t" << endl;
997 for (i = 0; i < nSurfacePts; ++i)
998 {
999 outfile
1000 << scientific << setw(17) << setprecision(16) << surfaceX[i]
1001 << " \t " << surfaceY[i] << " \t " << surfaceZ[i] << " \t "
1002 << surfaceFieldsAdded[0][i] << " \t " << surfaceFieldsAdded[1][i]
1003 << " \t " << surfaceFieldsAdded[2][i] << " \t "
1004 << surfaceFieldsAdded[3][i] << " \t " << surfaceFieldsAdded[4][i]
1005 << " \t " << surfaceFieldsAdded[5][i] << " \t "
1006 << surfaceFieldsAdded[6][i] << " \t " << surfaceFieldsAdded[7][i]
1007 << " \t " << surfaceFieldsAdded[8][i] << " \t "
1008 << surfaceFields[0][i] << " \t " << surfaceFields[1][i] << " \t "
1009 << surfaceFields[2][i] << " \t " << surfaceFields[3][i] << " \t "
1010 << surfaceFields[4][i] << " \t " << surfaceFieldsAdded[9][i]
1011 << " \t " << surfaceFieldsAdded[10][i] << " \t "
1012 << surfaceFieldsAdded[11][i] << " \t " << surfaceFieldsAdded[12][i]
1013 << " \t " << surfaceFieldsAdded[13][i] << " \t "
1014 << surfaceFieldsAdded[14][i] << " \t " << surfaceFieldsAdded[15][i]
1015 << " \t " << surfaceFieldsAdded[16][i] << " \t "
1016 << surfaceFieldsAdded[17][i] << " \t " << surfaceFieldsAdded[18][i]
1017 << " \t " << surfaceFieldsAdded[19][i] << " \t "
1018 << surfaceFieldsAdded[20][i] << " \t " << surfaceFieldsAdded[21][i]
1019 << " \t " << surfaceFieldsAdded[22][i] << " \t "
1020 << surfaceFieldsAdded[23][i] << " \t " << surfaceFieldsAdded[24][i]
1021 << " \t " << surfaceFieldsAdded[25][i] << " \t "
1022 << surfaceFieldsAdded[26][i] << " \t " << surfaceFieldsAdded[27][i]
1023 << " \t " << surfaceFieldsAdded[28][i] << " \t "
1024 << surfaceFieldsAdded[29][i] << " \t " << surfaceFieldsAdded[30][i]
1025 << " \t " << surfaceFieldsAdded[31][i] << " \t "
1026 << surfaceFieldsAdded[32][i] << " \t " << surfaceFieldsAdded[33][i]
1027 << " \t " << endl;
1028 }
1029 outfile << endl << endl;
1030 outfile.close();
1031
1032 return 0;
1033}
m_mu
NekDouble m_mu
Definition: CompressibleBL.cpp:81
m_Twall
NekDouble m_Twall
Definition: CompressibleBL.cpp:73
m_rhoInf
NekDouble m_rhoInf
Definition: CompressibleBL.cpp:78
m_vInf
NekDouble m_vInf
Definition: CompressibleBL.cpp:80
m_uInf
NekDouble m_uInf
Definition: CompressibleBL.cpp:77
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
main
int main(int argc, char *argv[])
Definition: ExtractSurface3DCFS.cpp:64
Nektar::MemoryManager
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
Definition: NekMemoryManager.hpp:81
CG_Iterations.pressure
pressure
Definition: CG_Iterations.py:229
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:288
Nektar::LibUtilities::SessionReaderSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: SessionReader.h:113
Nektar::LibUtilities::ePolyEvenlySpaced
@ ePolyEvenlySpaced
1D Evenly-spaced points using Lagrange polynomial
Definition: PointsType.h:73
Nektar::MultiRegions::ExpListSharedPtr
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Definition: LocTraceToTraceMap.h:56
Nektar::SpatialDomains::MeshGraphSharedPtr
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:174
Nektar::UnitTests::z
std::vector< double > z(NPUPPER)
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Vmath::Vsqrt
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.hpp:340
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.hpp:72
Vmath::Vvtvp
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.hpp:366
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.hpp:180
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.hpp:100
Vmath::Vdiv
void Vdiv(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x/y.
Definition: Vmath.hpp:126
Vmath::Fill
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.hpp:54
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825
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.hpp:220
std
STL namespace.
tinysimd::sqrt
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:294
Generated on Wed Oct 30 2024 16:07:33 for Nektar++ by doxygen 1.9.4