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