Nektar++
Functions
ExtractSurface2DCFS.cpp File Reference
#include <cstdio>
#include <cstdlib>
#include <iomanip>
#include <iostream>
#include <string>
#include <MultiRegions/AssemblyMap/AssemblyMapDG.h>
#include <MultiRegions/ExpList2DHomogeneous1D.h>
#include <MultiRegions/ExpList3DHomogeneous1D.h>
#include <MultiRegions/ExpList3DHomogeneous2D.h>
#include <LocalRegions/Expansion2D.h>
#include <LocalRegions/MatrixKey.h>
#include <MultiRegions/DisContField.h>
#include <LibUtilities/BasicUtils/ErrorUtil.hpp>
#include <LibUtilities/BasicUtils/FieldIO.h>
#include <LibUtilities/BasicUtils/NekFactory.hpp>
#include <LibUtilities/BasicUtils/SessionReader.h>
#include <LibUtilities/BasicUtils/SharedArray.hpp>
#include <LibUtilities/Memory/NekMemoryManager.hpp>
#include <MultiRegions/ContField.h>
#include <SpatialDomains/MeshGraphIO.h>
#include <SolverUtils/SolverUtilsDeclspec.h>
#include <boost/program_options.hpp>

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[10] Du_y: traceFieldsAdded[11] Dv_x: traceFieldsAdded[12] Dv_y: traceFieldsAdded[13]

Definition at line 70 of file ExtractSurface2DCFS.cpp.

71{
72 string fname;
73 po::options_description desc("Available options");
74 desc.add_options()("help,h", "Produce this help message.");
75
76 po::options_description hidden("Hidden options");
77 hidden.add_options()("input-file", po::value<vector<string>>(),
78 "input filename");
79
80 po::options_description cmdline_options;
81 cmdline_options.add(desc).add(hidden);
82
83 po::options_description visible("Allowed options");
84 visible.add(desc);
85
86 po::positional_options_description p;
87 p.add("input-file", -1);
88
89 po::variables_map vm;
90
91 try
92 {
93 po::store(po::command_line_parser(argc, argv)
94 .options(cmdline_options)
95 .positional(p)
96 .run(),
97 vm);
98 po::notify(vm);
99 }
100 catch (const std::exception &e)
101 {
102 cerr << e.what() << endl;
103 cerr << desc;
104 return 1;
105 }
106
107 std::vector<std::string> filenames;
108 if (vm.count("input-file"))
109 {
110 filenames = vm["input-file"].as<std::vector<std::string>>();
111 }
112
113 if (vm.count("help") || vm.count("input-file") != 1)
114 {
115 cerr << "Description: Extracts a surface from a 2D .fld file created "
116 "by the CompressibleFlowSolver and places it into a .cfs file"
117 << endl;
118 cerr << "Usage: ExtractSurface2DCFS [options] meshFile fieldFile"
119 << endl;
120 cout << desc;
121 return 1;
122 }
123
124 LibUtilities::SessionReaderSharedPtr vSession =
125 LibUtilities::SessionReader::CreateInstance(argc, argv);
126 SpatialDomains::MeshGraphSharedPtr graphShPt =
127 SpatialDomains::MeshGraphIO::Read(vSession);
128
129 fname = vSession->GetSessionName() + ".cfs";
130
131 // Find .fld or .chk file name
132 std::string fieldFile;
133 for (auto &file : filenames)
134 {
135 if (file.size() > 4 && (file.substr(file.size() - 4, 4) == ".fld" ||
136 file.substr(file.size() - 4, 4) == ".chk"))
137 {
138 fieldFile = file;
139 break;
140 }
141 }
142
143 int cnt;
144 int id1 = 0;
145 int id2 = 0;
146 int i, j, n, e, b;
147 Array<OneD, NekDouble> auxArray;
148
149 int nBndEdgePts, nBndEdges, nBndRegions;
150
151 std::string m_ViscosityType;
152 NekDouble m_gamma;
153 NekDouble m_pInf;
154 NekDouble m_rhoInf;
155 NekDouble m_uInf;
156 NekDouble m_vInf;
157 NekDouble m_wInf;
158 NekDouble m_gasConstant;
159 NekDouble m_Twall;
160 NekDouble m_mu;
161
162 int m_spacedim = 2;
163 int nDimensions = m_spacedim;
164 int phys_offset;
165
166 // Get gamma parameter from session file.
167 ASSERTL0(vSession->DefinesParameter("Gamma"),
168 "Compressible flow sessions must define a Gamma parameter.");
169 vSession->LoadParameter("Gamma", m_gamma, 1.4);
170
171 // Get E0 parameter from session file.
172 ASSERTL0(vSession->DefinesParameter("pInf"),
173 "Compressible flow sessions must define a pInf parameter.");
174 vSession->LoadParameter("pInf", m_pInf, 101325);
175
176 // Get rhoInf parameter from session file.
177 ASSERTL0(vSession->DefinesParameter("rhoInf"),
178 "Compressible flow sessions must define a rhoInf parameter.");
179 vSession->LoadParameter("rhoInf", m_rhoInf, 1.225);
180
181 // Get uInf parameter from session file.
182 ASSERTL0(vSession->DefinesParameter("uInf"),
183 "Compressible flow sessions must define a uInf parameter.");
184 vSession->LoadParameter("uInf", m_uInf, 0.1);
185
186 // Get vInf parameter from session file.
187 if (m_spacedim == 2 || m_spacedim == 3)
188 {
189 ASSERTL0(vSession->DefinesParameter("vInf"),
190 "Compressible flow sessions must define a vInf parameter"
191 "for 2D/3D problems.");
192 vSession->LoadParameter("vInf", m_vInf, 0.0);
193 }
194
195 // Get wInf parameter from session file.
196 if (m_spacedim == 3)
197 {
198 ASSERTL0(vSession->DefinesParameter("wInf"),
199 "Compressible flow sessions must define a wInf parameter"
200 "for 3D problems.");
201 vSession->LoadParameter("wInf", m_wInf, 0.0);
202 }
203
204 vSession->LoadParameter("GasConstant", m_gasConstant, 287.058);
205 vSession->LoadParameter("Twall", m_Twall, 300.15);
206 vSession->LoadSolverInfo("ViscosityType", m_ViscosityType, "Constant");
207 vSession->LoadParameter("mu", m_mu, 1.78e-05);
208
209 //--------------------------------------------------------------------------
210 // Import field file
211 vector<LibUtilities::FieldDefinitionsSharedPtr> fieldDef;
212 vector<vector<NekDouble>> fieldData;
213
214 LibUtilities::FieldMetaDataMap fieldMetaDataMap;
215 LibUtilities::Import(fieldFile, fieldDef, fieldData, fieldMetaDataMap);
216 //--------------------------------------------------------------------------
217
218 //--------------------------------------------------------------------------
219 // Set up Expansion information
220 vector<vector<LibUtilities::PointsType>> pointsType;
221 for (i = 0; i < fieldDef.size(); ++i)
222 {
223 vector<LibUtilities::PointsType> ptype;
224 for (j = 0; j < 2; ++j)
225 {
226 ptype.push_back(LibUtilities::ePolyEvenlySpaced);
227 }
228 pointsType.push_back(ptype);
229 }
230 graphShPt->SetExpansionInfo(fieldDef, pointsType);
231
232 //--------------------------------------------------------------------------
233
234 //--------------------------------------------------------------------------
235 // Define Expansion
236 int nfields = vSession->GetVariables().size();
237 Array<OneD, MultiRegions::ExpListSharedPtr> Exp(nfields);
238 Array<OneD, MultiRegions::ExpListSharedPtr> pFields(nfields);
239
240 for (i = 0; i < pFields.size(); i++)
241 {
242 pFields[i] =
243 MemoryManager<MultiRegions::DisContField>::AllocateSharedPtr(
244 vSession, graphShPt, vSession->GetVariable(i));
245 }
246
247 //@TODO: Might need this to rotate mesh based on time
248 for (auto &fld : pFields)
249 {
250 if (fld->GetGraph()->GetMovement() != nullptr)
251 {
252 fld->GetGraph()->GetMovement()->PerformMovement(
253 std::stod(fieldMetaDataMap["Time"]));
254 fld->Reset();
255 fld->SetUpPhysNormals();
256 }
257 }
258
259 MultiRegions::ExpListSharedPtr Exp2D;
260 Exp2D = MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(vSession,
261 graphShPt);
262
263 Exp[0] = Exp2D;
264
265 for (i = 1; i < nfields; ++i)
266 {
267 Exp[i] =
268 MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(*Exp2D);
269 }
270
271 int nSolutionPts = pFields[0]->GetNpoints();
272 int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
273 int nElements = pFields[0]->GetExpSize();
274
275 Array<OneD, NekDouble> tmp(nSolutionPts, 0.0);
276
277 Array<OneD, NekDouble> x(nSolutionPts);
278 Array<OneD, NekDouble> y(nSolutionPts);
279 Array<OneD, NekDouble> z(nSolutionPts);
280
281 Array<OneD, NekDouble> traceX(nTracePts);
282 Array<OneD, NekDouble> traceY(nTracePts);
283 Array<OneD, NekDouble> traceZ(nTracePts);
284
285 Array<OneD, NekDouble> surfaceX(nTracePts);
286 Array<OneD, NekDouble> surfaceY(nTracePts);
287 Array<OneD, NekDouble> surfaceZ(nTracePts);
288
289 pFields[0]->GetCoords(x, y, z);
290
291 pFields[0]->ExtractTracePhys(x, traceX);
292 pFields[0]->ExtractTracePhys(y, traceY);
293 pFields[0]->ExtractTracePhys(z, traceZ);
294 //--------------------------------------------------------------------------
295
296 //--------------------------------------------------------------------------
297 // Copy data from field file
298 Array<OneD, Array<OneD, NekDouble>> uFields(nfields);
299 Array<OneD, Array<OneD, NekDouble>> traceFields(nfields);
300 Array<OneD, Array<OneD, NekDouble>> surfaceFields(nfields);
301
302 // Extract the physical values of the solution at the boundaries
303 for (j = 0; j < nfields; ++j)
304 {
305 uFields[j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
306 traceFields[j] = Array<OneD, NekDouble>(nTracePts, 0.0);
307 surfaceFields[j] = Array<OneD, NekDouble>(nTracePts, 0.0);
308
309 for (i = 0; i < fieldData.size(); ++i)
310 {
311 Exp[j]->ExtractDataToCoeffs(fieldDef[i], fieldData[i],
312 fieldDef[i]->m_fields[j],
313 Exp[j]->UpdateCoeffs());
314 }
315 Exp[j]->BwdTrans(Exp[j]->GetCoeffs(), Exp[j]->UpdatePhys());
316 Vmath::Vcopy(nSolutionPts, Exp[j]->GetPhys(), 1, uFields[j], 1);
317 pFields[0]->ExtractTracePhys(uFields[j], traceFields[j]);
318 }
319
320 // Fields to add in the output file
321
322 int nfieldsAdded = 20;
323 Array<OneD, Array<OneD, NekDouble>> traceFieldsAdded(nfieldsAdded);
324 Array<OneD, Array<OneD, NekDouble>> surfaceFieldsAdded(nfieldsAdded);
325
326 for (j = 0; j < nfieldsAdded; ++j)
327 {
328 traceFieldsAdded[j] = Array<OneD, NekDouble>(nTracePts, 0.0);
329 surfaceFieldsAdded[j] = Array<OneD, NekDouble>(nTracePts, 0.0);
330 }
331
332 /******** Evaluation of normals and tangents on the trace *****************
333 * nx -> traceFieldsAdded[0];
334 * ny -> traceFieldsAdded[1];
335 * tx -> traceFieldsAdded[2];
336 * ty -> traceFieldsAdded[3];
337 ***************************************************************************/
338
339 Array<OneD, Array<OneD, NekDouble>> m_traceNormals(nDimensions);
340 for (i = 0; i < nDimensions; ++i)
341 {
342 m_traceNormals[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
343 }
344 pFields[0]->GetTrace()->GetNormals(m_traceNormals);
345
346 Array<OneD, Array<OneD, NekDouble>> m_traceTangents(nDimensions);
347 for (i = 0; i < nDimensions; ++i)
348 {
349 m_traceTangents[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
350 }
351
352 // nx
353 Vmath::Vcopy(nTracePts, &m_traceNormals[0][0], 1, &traceFieldsAdded[0][0],
354 1);
355
356 // ny
357 Vmath::Vcopy(nTracePts, &m_traceNormals[1][0], 1, &traceFieldsAdded[1][0],
358 1);
359
360 // t_x = - n_y
361 Vmath::Vcopy(nTracePts, &m_traceNormals[1][0], 1, &m_traceTangents[0][0],
362 1);
363 Vmath::Neg(nTracePts, &m_traceTangents[0][0], 1);
364
365 Vmath::Vcopy(nTracePts, &m_traceTangents[0][0], 1, &traceFieldsAdded[2][0],
366 1);
367
368 // t_y = n_x
369 Vmath::Vcopy(nTracePts, &m_traceNormals[0][0], 1, &m_traceTangents[1][0],
370 1);
371
372 Vmath::Vcopy(nTracePts, &m_traceTangents[1][0], 1, &traceFieldsAdded[3][0],
373 1);
374
375 /******** Evaluation of the pressure ***************************************
376 * P = (E-1/2.*rho.*((rhou./rho).^2+(rhov./rho).^2))*(gamma - 1);
377 * P -> traceFieldsAdded[4];
378 ***************************************************************************/
379
380 Array<OneD, NekDouble> pressure(nSolutionPts, 0.0);
381 NekDouble gammaMinusOne = m_gamma - 1.0;
382
383 for (i = 0; i < m_spacedim; i++)
384 {
385 Vmath::Vmul(nSolutionPts, &uFields[i + 1][0], 1, &uFields[i + 1][0], 1,
386 &tmp[0], 1);
387
388 Vmath::Smul(nSolutionPts, 0.5, &tmp[0], 1, &tmp[0], 1);
389
390 Vmath::Vadd(nSolutionPts, &pressure[0], 1, &tmp[0], 1, &pressure[0], 1);
391 }
392
393 Vmath::Vdiv(nSolutionPts, &pressure[0], 1, &uFields[0][0], 1, &pressure[0],
394 1);
395
396 Vmath::Vsub(nSolutionPts, &uFields[nfields - 1][0], 1, &pressure[0], 1,
397 &pressure[0], 1);
398
399 Vmath::Smul(nSolutionPts, gammaMinusOne, &pressure[0], 1, &pressure[0], 1);
400
401 // Extract trace
402 pFields[0]->ExtractTracePhys(pressure, traceFieldsAdded[4]);
403
404 /******** Evaluation of the temperature ************************************
405 * T = P/(R*rho);
406 * T -> traceFieldsAdded[5];
407 ***************************************************************************/
408
409 Array<OneD, NekDouble> temperature(nSolutionPts, 0.0);
410
411 Vmath::Vdiv(nSolutionPts, &pressure[0], 1, &uFields[0][0], 1,
412 &temperature[0], 1);
413
414 NekDouble GasConstantInv = 1.0 / m_gasConstant;
415 Vmath::Smul(nSolutionPts, GasConstantInv, &temperature[0], 1,
416 &temperature[0], 1);
417
418 // Extract trace
419 pFields[0]->ExtractTracePhys(temperature, traceFieldsAdded[5]);
420
421 /*** Evaluation of the temperature gradient in the normal direction ********
422 * DT_n -> traceFieldsAdded[6]
423 ***************************************************************************/
424
425 Array<OneD, Array<OneD, NekDouble>> Dtemperature(nDimensions);
426 Array<OneD, Array<OneD, NekDouble>> traceDtemperature(nDimensions);
427
428 for (i = 0; i < nDimensions; ++i)
429 {
430 Dtemperature[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
431 traceDtemperature[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
432 }
433
434 for (i = 0; i < nDimensions; ++i)
435 {
436 for (n = 0; n < nElements; n++)
437 {
438 phys_offset = pFields[0]->GetPhys_Offset(n);
439
440 pFields[i]->GetExp(n)->PhysDeriv(i, temperature + phys_offset,
441 auxArray =
442 Dtemperature[i] + phys_offset);
443 }
444 // Extract trace
445 pFields[0]->ExtractTracePhys(Dtemperature[i], traceDtemperature[i]);
446 }
447
448 for (i = 0; i < nDimensions; ++i)
449 {
450 Vmath::Vmul(nTracePts, &m_traceNormals[i][0], 1,
451 &traceDtemperature[i][0], 1, &tmp[0], 1);
452
453 Vmath::Vadd(nTracePts, &traceFieldsAdded[6][0], 1, &tmp[0], 1,
454 &traceFieldsAdded[6][0], 1);
455 }
456
457 /*** Evaluation of the pressure gradient ***********************************
458 * DP_t -> traceFieldsAdded[7] tangent direction
459 * DP_x -> traceFieldsAdded[8]
460 * DP_y -> traceFieldsAdded[9]
461 ***************************************************************************/
462
463 Array<OneD, Array<OneD, NekDouble>> Dpressure(nDimensions);
464 Array<OneD, Array<OneD, NekDouble>> traceDpressure(nDimensions);
465
466 for (i = 0; i < nDimensions; ++i)
467 {
468 Dpressure[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
469 traceDpressure[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
470 }
471
472 for (i = 0; i < nDimensions; ++i)
473 {
474 for (n = 0; n < nElements; n++)
475 {
476 phys_offset = pFields[0]->GetPhys_Offset(n);
477
478 pFields[i]->GetExp(n)->PhysDeriv(i, pressure + phys_offset,
479 auxArray =
480 Dpressure[i] + phys_offset);
481 }
482 // Extract trace
483 pFields[0]->ExtractTracePhys(Dpressure[i], traceDpressure[i]);
484 }
485
486 // Dp_t
487 for (i = 0; i < nDimensions; ++i)
488 {
489 Vmath::Vmul(nTracePts, &m_traceTangents[i][0], 1, &traceDpressure[i][0],
490 1, &tmp[0], 1);
491
492 Vmath::Vadd(nTracePts, &traceFieldsAdded[7][0], 1, &tmp[0], 1,
493 &traceFieldsAdded[7][0], 1);
494 }
495
496 // Dp_x
497 Vmath::Vcopy(nTracePts, &traceDpressure[0][0], 1, &traceFieldsAdded[8][0],
498 1);
499
500 // Dp_y
501 Vmath::Vcopy(nTracePts, &traceDpressure[1][0], 1, &traceFieldsAdded[9][0],
502 1);
503
504 /** Evaluation of the velocity gradient in the cartesian directions
505 * Du_x: traceFieldsAdded[10]
506 * Du_y: traceFieldsAdded[11]
507 * Dv_x: traceFieldsAdded[12]
508 * Dv_y: traceFieldsAdded[13]
509 **/
510 Array<OneD, Array<OneD, Array<OneD, NekDouble>>> Dvelocity(nDimensions);
511 Array<OneD, Array<OneD, Array<OneD, NekDouble>>> traceDvelocity(
512 nDimensions);
513 Array<OneD, Array<OneD, NekDouble>> velocity(nDimensions);
514
515 for (i = 0; i < nDimensions; ++i)
516 {
517 Dvelocity[i] = Array<OneD, Array<OneD, NekDouble>>(nDimensions);
518 traceDvelocity[i] = Array<OneD, Array<OneD, NekDouble>>(nDimensions);
519 velocity[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
520
521 Vmath::Vdiv(nSolutionPts, uFields[i + 1], 1, uFields[0], 1, velocity[i],
522 1);
523
524 for (j = 0; j < nDimensions; ++j)
525 {
526 Dvelocity[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
527 traceDvelocity[i][j] = Array<OneD, NekDouble>(nTracePts, 0.0);
528 }
529 }
530
531 for (i = 0; i < nDimensions; ++i)
532 {
533 for (j = 0; j < nDimensions; ++j)
534 {
535 for (n = 0; n < nElements; n++)
536 {
537 phys_offset = pFields[0]->GetPhys_Offset(n);
538
539 pFields[i]->GetExp(n)->PhysDeriv(j, velocity[i] + phys_offset,
540 auxArray = Dvelocity[i][j] +
541 phys_offset);
542 }
543
544 // Extract trace
545 pFields[0]->ExtractTracePhys(Dvelocity[i][j], traceDvelocity[i][j]);
546 }
547 }
548
549 Vmath::Vcopy(nTracePts, &traceDvelocity[0][0][0], 1,
550 &traceFieldsAdded[10][0], 1);
551 Vmath::Vcopy(nTracePts, &traceDvelocity[0][1][0], 1,
552 &traceFieldsAdded[11][0], 1);
553 Vmath::Vcopy(nTracePts, &traceDvelocity[1][0][0], 1,
554 &traceFieldsAdded[12][0], 1);
555 Vmath::Vcopy(nTracePts, &traceDvelocity[1][1][0], 1,
556 &traceFieldsAdded[13][0], 1);
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) * (T_star + 110.0) /
583 (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, &Dvelocity[0][0][0], 1, &divVel[0],
600 1);
601 Vmath::Vadd(nSolutionPts, &divVel[0], 1, &Dvelocity[1][1][0], 1, &divVel[0],
602 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, &Dvelocity[1][0][0], 1,
626 &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, mach,
699 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().size())
713 {
714 id1 = 0;
715 cnt = 0;
716 nBndRegions = pFields[0]->GetBndCondExpansions().size();
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]
724 ->GetExp(e)
725 ->GetNumPoints(0);
726
727 id2 = pFields[0]->GetTrace()->GetPhys_Offset(
728 pFields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(
729 cnt++));
730
731 if (pFields[0]->GetBndConditions()[b]->GetUserDefined() ==
732 "WallViscous" ||
733 pFields[0]->GetBndConditions()[b]->GetUserDefined() ==
734 "WallAdiabatic" ||
735 pFields[0]->GetBndConditions()[b]->GetUserDefined() ==
736 "Wall")
737 {
738 Vmath::Vcopy(nBndEdgePts, &traceX[id2], 1, &surfaceX[id1],
739 1);
740
741 Vmath::Vcopy(nBndEdgePts, &traceY[id2], 1, &surfaceY[id1],
742 1);
743
744 Vmath::Vcopy(nBndEdgePts, &traceZ[id2], 1, &surfaceZ[id1],
745 1);
746
747 id1 += nBndEdgePts;
748 }
749 }
750 }
751 }
752
753 // Extract fields
754 if (pFields[0]->GetBndCondExpansions().size())
755 {
756 for (j = 0; j < nfields; ++j)
757 {
758 id1 = 0;
759 cnt = 0;
760 nBndRegions = pFields[j]->GetBndCondExpansions().size();
761 for (b = 0; b < nBndRegions; ++b)
762 {
763 nBndEdges = pFields[j]->GetBndCondExpansions()[b]->GetExpSize();
764 for (e = 0; e < nBndEdges; ++e)
765 {
766 nBndEdgePts = pFields[j]
767 ->GetBndCondExpansions()[b]
768 ->GetExp(e)
769 ->GetNumPoints(0);
770
771 id2 = pFields[j]->GetTrace()->GetPhys_Offset(
772 pFields[j]->GetTraceMap()->GetBndCondIDToGlobalTraceID(
773 cnt++));
774
775 if (pFields[j]->GetBndConditions()[b]->GetUserDefined() ==
776 "WallViscous" ||
777 pFields[j]->GetBndConditions()[b]->GetUserDefined() ==
778 "WallAdiabatic" ||
779 pFields[j]->GetBndConditions()[b]->GetUserDefined() ==
780 "Wall")
781 {
782 Vmath::Vcopy(nBndEdgePts, &traceFields[j][id2], 1,
783 &surfaceFields[j][id1], 1);
784
785 id1 += nBndEdgePts;
786 }
787 }
788 }
789 }
790 }
791
792 // Extract fields added
793 if (pFields[0]->GetBndCondExpansions().size())
794 {
795 for (j = 0; j < nfieldsAdded; ++j)
796 {
797 id1 = 0;
798 cnt = 0;
799 nBndRegions = pFields[0]->GetBndCondExpansions().size();
800 for (b = 0; b < nBndRegions; ++b)
801 {
802 nBndEdges = pFields[0]->GetBndCondExpansions()[b]->GetExpSize();
803 for (e = 0; e < nBndEdges; ++e)
804 {
805 nBndEdgePts = pFields[0]
806 ->GetBndCondExpansions()[b]
807 ->GetExp(e)
808 ->GetNumPoints(0);
809
810 id2 = pFields[0]->GetTrace()->GetPhys_Offset(
811 pFields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(
812 cnt++));
813
814 if (pFields[0]->GetBndConditions()[b]->GetUserDefined() ==
815 "WallViscous" ||
816 pFields[0]->GetBndConditions()[b]->GetUserDefined() ==
817 "WallAdiabatic" ||
818 pFields[0]->GetBndConditions()[b]->GetUserDefined() ==
819 "Wall")
820 {
821 Vmath::Vcopy(nBndEdgePts, &traceFieldsAdded[j][id2], 1,
822 &surfaceFieldsAdded[j][id1], 1);
823
824 id1 += nBndEdgePts;
825 }
826 }
827 }
828 }
829 }
830 //===================================================================================================
831 //===================================================================================================
832 //===================================================================================================
833 std::string vEquation = vSession->GetSolverInfo("EQType");
834
835 Array<OneD, MultiRegions::ExpListSharedPtr> BndExp;
836 BndExp = pFields[0]->GetBndCondExpansions();
837 StdRegions::StdExpansion1DSharedPtr bc;
838
839 NekDouble Fxp(0.0);
840 NekDouble Fyp(0.0);
841 NekDouble Fxv(0.0);
842 NekDouble Fyv(0.0);
843 NekDouble Sref(0.0);
844
845 int GlobalIndex(0);
846
847 for (int i = 0; i < BndExp[0]->GetExpSize(); ++i)
848 {
849 bc = BndExp[0]->GetExp(i)->as<LocalRegions::Expansion1D>();
850
851 int nbc = bc->GetTotPoints();
852
853 Array<OneD, NekDouble> nxOnBnd(nbc, 0.0);
854 Array<OneD, NekDouble> nyOnBnd(nbc, 0.0);
855 Array<OneD, NekDouble> txOnBnd(nbc, 0.0);
856 Array<OneD, NekDouble> tyOnBnd(nbc, 0.0);
857 // Array<OneD, NekDouble> CoeffAero(nbc,0.0);
858 // Array<OneD, NekDouble> tmp(nbc,0.0);
859
860 Array<OneD, NekDouble> drag_p(nbc, 0.0);
861 Array<OneD, NekDouble> lift_p(nbc, 0.0);
862 Array<OneD, NekDouble> PressurOnBnd(nbc, 0.0);
863
864 Array<OneD, NekDouble> drag_v(nbc, 0.0);
865 Array<OneD, NekDouble> lift_v(nbc, 0.0);
866 Array<OneD, NekDouble> ShearStressOnBnd(nbc, 0.0);
867
868 Array<OneD, NekDouble> Unity(nbc, 1.0);
869
870 for (int j = 0; j < nbc; ++j)
871 {
872
873 nxOnBnd[j] = surfaceFieldsAdded[0][GlobalIndex];
874 nyOnBnd[j] = surfaceFieldsAdded[1][GlobalIndex];
875 txOnBnd[j] = surfaceFieldsAdded[2][GlobalIndex];
876 tyOnBnd[j] = surfaceFieldsAdded[3][GlobalIndex];
877
878 PressurOnBnd[j] = surfaceFieldsAdded[4][GlobalIndex];
879
880 if (vEquation == "NavierStokesCFE")
881 {
882 ShearStressOnBnd[j] = surfaceFieldsAdded[17][GlobalIndex];
883 }
884
885 // CoeffAero[j] = surfaceFields[0][GlobalIndex];
886 // tmp[j] =
887 // surfaceFields[1][GlobalIndex]*surfaceFields[1][GlobalIndex];
888 // tmp[j] = tmp[j] +
889 // surfaceFields[2][GlobalIndex]*surfaceFields[2][GlobalIndex];
890 // tmp[j] = sqrt(tmp[j]);
891 // CoeffAero[j] = CoeffAero[j]*tmp[j];
892 // CoeffAero[j] = 1.0/CoeffAero[j];
893 //
894 // PressurOnBnd[j] =
895 // CoeffAero[j]*surfaceFieldsAdded[4][GlobalIndex];
896 //
897 // cout << "CoeffAero = " << CoeffAero[j] << endl;
898
899 GlobalIndex++;
900 }
901
902 Vmath::Vmul(nbc, PressurOnBnd, 1, nxOnBnd, 1, drag_p, 1);
903 Vmath::Vmul(nbc, PressurOnBnd, 1, nyOnBnd, 1, lift_p, 1);
904
905 // Vmath::Vmul(nbc,drag_p,1,CoeffAero,1, drag_p,1);
906 // Vmath::Vmul(nbc,lift_p,1,CoeffAero,1, lift_p,1);
907
908 Fxp += bc->Integral(drag_p);
909 Fyp += bc->Integral(lift_p);
910
911 if (vEquation == "NavierStokesCFE")
912 {
913 Vmath::Vmul(nbc, ShearStressOnBnd, 1, txOnBnd, 1, drag_v, 1);
914 Vmath::Vmul(nbc, ShearStressOnBnd, 1, tyOnBnd, 1, lift_v, 1);
915
916 // Vmath::Vdiv(nbc,drag_v,1,CoeffAero,1, drag_v,1);
917 // Vmath::Vdiv(nbc,lift_v,1,CoeffAero,1, lift_v,1);
918
919 Fxv += bc->Integral(drag_v);
920 Fyv += bc->Integral(lift_v);
921 }
922
923 Sref += bc->Integral(Unity);
924 }
925
926 cout << "\n Sref = " << Sref << endl;
927 Fxp = Fxp / Sref;
928 Fyp = Fyp / Sref;
929 Fxv = Fxv / Sref;
930 Fyv = Fyv / Sref;
931 cout << " Pressure drag (Fxp) = " << Fxp << endl;
932 cout << " Pressure lift (Fyp) = " << Fyp << endl;
933 cout << " Viscous drag (Fxv) = " << Fxv << endl;
934 cout << " Viscous lift (Fyv) = " << Fyv << endl;
935 cout << "\n ==> Total drag = " << Fxp + Fxv << endl;
936 cout << " ==> Total lift = " << Fyp + Fyv << "\n" << endl;
937
938 //===================================================================================================
939 //===================================================================================================
940 //===================================================================================================
941
942 // Print the surface coordinates and the surface solution in a .txt file
943 ofstream outfile;
944 outfile.open(fname.c_str());
945 outfile << "% x[m] "
946 << " \t"
947 << "y[m] "
948 << " \t"
949 << "z[m] "
950 << " \t"
951 << "nx[] "
952 << " \t"
953 << "ny[] "
954 << " \t"
955 << "tx[] "
956 << " \t"
957 << "ty[] "
958 << " \t"
959 << "rho[kg/m^3] "
960 << " \t"
961 << "rhou[kg/(m^2 s)] "
962 << " \t"
963 << "rhov[kg/(m^2 s)] "
964 << " \t"
965 << "E[Pa] "
966 << " \t"
967 << "p[Pa] "
968 << " \t"
969 << "T[k] "
970 << " \t"
971 << "dT/dn[k/m] "
972 << " \t"
973 << "dp/dT[Pa/m] "
974 << " \t"
975 << "dp/dx[Pa/m] "
976 << " \t"
977 << "dp/dy[Pa/m] "
978 << " \t"
979 << "du/dx[s^-1] "
980 << " \t"
981 << "du/dy[s^-1] "
982 << " \t"
983 << "dv/dx[s^-1] "
984 << " \t"
985 << "dv/dy[s^-1] "
986 << " \t"
987 << "tau_xx[Pa] "
988 << " \t"
989 << "tau_yy[Pa] "
990 << " \t"
991 << "tau_xy[Pa] "
992 << " \t"
993 << "tau_t[Pa] "
994 << " \t"
995 << "mu[Pa s] "
996 << " \t"
997 << "M[] "
998 << " \t"
999 // << "Fxp " << " \t"
1000 << endl;
1001 for (i = 0; i < id1; ++i)
1002 {
1003 outfile << scientific << setw(17) << setprecision(16) << surfaceX[i]
1004 << " \t " << surfaceY[i] << " \t " << surfaceZ[i] << " \t "
1005 << surfaceFieldsAdded[0][i] << " \t "
1006 << surfaceFieldsAdded[1][i] << " \t "
1007 << surfaceFieldsAdded[2][i] << " \t "
1008 << surfaceFieldsAdded[3][i] << " \t " << surfaceFields[0][i]
1009 << " \t " << surfaceFields[1][i] << " \t "
1010 << surfaceFields[2][i] << " \t " << surfaceFields[3][i]
1011 << " \t " << surfaceFieldsAdded[4][i] << " \t "
1012 << surfaceFieldsAdded[5][i] << " \t "
1013 << surfaceFieldsAdded[6][i] << " \t "
1014 << surfaceFieldsAdded[7][i] << " \t "
1015 << surfaceFieldsAdded[8][i] << " \t "
1016 << surfaceFieldsAdded[9][i] << " \t "
1017 << surfaceFieldsAdded[10][i] << " \t "
1018 << surfaceFieldsAdded[11][i] << " \t "
1019 << surfaceFieldsAdded[12][i] << " \t "
1020 << surfaceFieldsAdded[13][i] << " \t "
1021 << surfaceFieldsAdded[14][i] << " \t "
1022 << surfaceFieldsAdded[15][i] << " \t "
1023 << surfaceFieldsAdded[16][i] << " \t "
1024 << surfaceFieldsAdded[17][i] << " \t "
1025 << surfaceFieldsAdded[18][i] << " \t "
1026 << surfaceFieldsAdded[19][i]
1027 << " \t "
1028 // << Fxp << " \t "
1029 << endl;
1030 }
1031 outfile << endl << endl;
1032 outfile.close();
1033
1034 return 0;
1035}
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
Nektar::MemoryManager
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
Definition: NekMemoryManager.hpp:81
Nektar::StdRegions::StdExpansion::GetTotPoints
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:134
CG_Iterations.pressure
pressure
Definition: CG_Iterations.py:229
CellMLToNektar.cellml_metadata.p
p
Definition: cellml_metadata.py:350
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::FieldMetaDataMap
std::map< std::string, std::string > FieldMetaDataMap
Definition: FieldIO.h:50
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::StdRegions::StdExpansion1DSharedPtr
std::shared_ptr< StdExpansion1D > StdExpansion1DSharedPtr
Definition: StdExpansion1D.h:111
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::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.hpp:292
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
tinysimd::sqrt
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:294

References ASSERTL0, Nektar::LibUtilities::ePolyEvenlySpaced, Vmath::Fill(), Nektar::StdRegions::StdExpansion::GetTotPoints(), Nektar::LibUtilities::Import(), m_mu, m_rhoInf, m_Twall, m_uInf, m_vInf, Vmath::Neg(), CellMLToNektar.cellml_metadata::p, CG_Iterations::pressure, CellMLToNektar.translators::run(), 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 Oct 30 2024 16:07:35 for Nektar++ by doxygen 1.9.4