Nektar++
Functions
ExtractSurface3DCFS.cpp File Reference
#include <cstdio>
#include <cstdlib>
#include <iomanip>
#include <iostream>
#include <string>
#include <MultiRegions/AssemblyMap/AssemblyMapDG.h>
#include <MultiRegions/ExpList.h>
#include <LocalRegions/Expansion.h>
#include <LocalRegions/Expansion3D.h>
#include <LocalRegions/MatrixKey.h>
#include <MultiRegions/DisContField.h>
#include <LibUtilities/BasicUtils/FieldIO.h>
#include <LibUtilities/BasicUtils/NekFactory.hpp>
#include <LibUtilities/BasicUtils/SessionReader.h>
#include <LibUtilities/BasicUtils/SharedArray.hpp>
#include <LibUtilities/Communication/Comm.h>
#include <LibUtilities/Memory/NekMemoryManager.hpp>
#include <MultiRegions/ContField.h>
#include <SpatialDomains/MeshGraph.h>
#include <SolverUtils/SolverUtilsDeclspec.h>

Go to the source code of this file.

Functions

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

Function Documentation

◆ main()

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

Evaluation of the velocity gradient in the cartesian directions Du_x: traceFieldsAdded[17] Du_y: traceFieldsAdded[18] Du_z: traceFieldsAdded[19] Dv_x: traceFieldsAdded[20] Dv_y: traceFieldsAdded[21] Dv_z: traceFieldsAdded[22] Dw_x: traceFieldsAdded[23] Dw_y: traceFieldsAdded[24] Dw_z: traceFieldsAdded[25]

Definition at line 65 of file ExtractSurface3DCFS.cpp.

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

References ASSERTL0, Nektar::LibUtilities::ePolyEvenlySpaced, Vmath::Fill(), Nektar::LibUtilities::Import(), m_mu, m_rhoInf, m_Twall, m_uInf, m_vInf, CG_Iterations::pressure, Vmath::Smul(), tinysimd::sqrt(), Vmath::Vadd(), Vmath::Vcopy(), Vmath::Vdiv(), Vmath::Vmul(), Vmath::Vsqrt(), Vmath::Vsub(), Vmath::Vvtvp(), and Nektar::UnitTests::z().

Generated on Wed Jan 31 2024 08:26:51 for Nektar++ by doxygen 1.9.4