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]
73 string fname = std::string(argv[2]);
74 int fdot = fname.find_last_of(
'.');
75 if (fdot != std::string::npos)
77 string ending = fname.substr(fdot);
82 if (ending ==
".chk" || ending ==
".fld")
84 fname = fname.substr(0,fdot);
88 fname = fname +
".txt";
96 int nBndEdgePts, nBndEdges, nBndRegions;
101 "Usage: ExtractSurface2DCFS meshfile fieldFile\n");
103 "Extracts a surface from a 2D fld file"
104 "(only for CompressibleFlowSolver and purely 2D .fld files)\n");
109 = LibUtilities::SessionReader::CreateInstance(3, argv);
111 std::string m_ViscosityType;
125 int nDimensions = m_spacedim;
129 ASSERTL0(vSession->DefinesParameter(
"Gamma"),
130 "Compressible flow sessions must define a Gamma parameter.");
131 vSession->LoadParameter(
"Gamma", m_gamma, 1.4);
134 ASSERTL0(vSession->DefinesParameter(
"pInf"),
135 "Compressible flow sessions must define a pInf parameter.");
136 vSession->LoadParameter(
"pInf", m_pInf, 101325);
139 ASSERTL0(vSession->DefinesParameter(
"rhoInf"),
140 "Compressible flow sessions must define a rhoInf parameter.");
141 vSession->LoadParameter(
"rhoInf", m_rhoInf, 1.225);
144 ASSERTL0(vSession->DefinesParameter(
"uInf"),
145 "Compressible flow sessions must define a uInf parameter.");
146 vSession->LoadParameter(
"uInf", m_uInf, 0.1);
149 if (m_spacedim == 2 || m_spacedim == 3)
151 ASSERTL0(vSession->DefinesParameter(
"vInf"),
152 "Compressible flow sessions must define a vInf parameter"
153 "for 2D/3D problems.");
154 vSession->LoadParameter(
"vInf", m_vInf, 0.0);
160 ASSERTL0(vSession->DefinesParameter(
"wInf"),
161 "Compressible flow sessions must define a wInf parameter"
163 vSession->LoadParameter(
"wInf", m_wInf, 0.0);
166 vSession->LoadParameter (
"GasConstant", m_gasConstant, 287.058);
167 vSession->LoadParameter (
"Twall", m_Twall, 300.15);
168 vSession->LoadSolverInfo(
"ViscosityType", m_ViscosityType,
"Constant");
169 vSession->LoadParameter (
"mu", m_mu, 1.78e-05);
170 vSession->LoadParameter (
"thermalConductivity",
171 m_thermalConductivity, 0.0257);
175 string meshfile(argv[1]);
177 SpatialDomains::MeshGraph::Read(vSession);
182 string fieldFile(argv[2]);
183 vector<LibUtilities::FieldDefinitionsSharedPtr> fieldDef;
184 vector<vector<NekDouble> > fieldData;
191 vector< vector<LibUtilities::PointsType> > pointsType;
192 for (i = 0; i < fieldDef.size(); ++i)
194 vector<LibUtilities::PointsType> ptype;
195 for (j = 0; j < 2; ++j)
199 pointsType.push_back(ptype);
201 graphShPt->SetExpansions(fieldDef, pointsType);
208 int nfields = fieldDef[0]->m_fields.size();
212 for(i = 0; i < pFields.num_elements(); i++)
216 vSession, graphShPt, vSession->GetVariable(i));
225 for (i = 1; i < nfields; ++i)
231 int nSolutionPts = pFields[0]->GetNpoints();
232 int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
233 int nElements = pFields[0]->GetExpSize();
249 pFields[0]->GetCoords(x, y, z);
251 pFields[0]->ExtractTracePhys(x, traceX);
252 pFields[0]->ExtractTracePhys(y, traceY);
253 pFields[0]->ExtractTracePhys(z, traceZ);
263 for (j = 0; j < nfields; ++j)
270 for (i = 0; i < fieldData.size(); ++i)
272 Exp[j]->ExtractDataToCoeffs(fieldDef[i], fieldData[i],
273 fieldDef[i]->m_fields[j],
274 Exp[j]->UpdateCoeffs());
276 Exp[j]->BwdTrans(Exp[j]->GetCoeffs(), Exp[j]->UpdatePhys());
277 Vmath::Vcopy(nSolutionPts, Exp[j]->GetPhys(), 1, uFields[j], 1);
278 pFields[0]->ExtractTracePhys(uFields[j], traceFields[j]);
283 int nfieldsAdded = 20;
287 for (j = 0; j < nfieldsAdded; ++j)
301 for(i = 0; i < nDimensions; ++i)
305 pFields[0]->GetTrace()->GetNormals(m_traceNormals);
308 for(i = 0; i < nDimensions; ++i)
316 &m_traceNormals[0][0], 1,
317 &traceFieldsAdded[0][0], 1);
321 &m_traceNormals[1][0], 1,
322 &traceFieldsAdded[1][0], 1);
326 &m_traceNormals[1][0], 1,
327 &m_traceTangents[0][0], 1);
328 Vmath::Neg(nTracePts, &m_traceTangents[0][0], 1);
331 &m_traceTangents[0][0], 1,
332 &traceFieldsAdded[2][0], 1);
336 &m_traceNormals[0][0], 1,
337 &m_traceTangents[1][0], 1);
340 &m_traceTangents[1][0], 1,
341 &traceFieldsAdded[3][0], 1);
351 for (i = 0; i < m_spacedim; i++)
354 &uFields[i + 1][0], 1,
355 &uFields[i + 1][0], 1,
375 &uFields[nfields - 1][0], 1,
384 pFields[0]->ExtractTracePhys(pressure, traceFieldsAdded[4]);
398 NekDouble GasConstantInv = 1.0/m_gasConstant;
404 pFields[0]->ExtractTracePhys(temperature, traceFieldsAdded[5]);
413 for (i = 0; i < nDimensions; ++ i)
419 for (i = 0; i < nDimensions; ++ i)
421 for (n = 0; n < nElements; n++)
423 phys_offset = pFields[0]->GetPhys_Offset(n);
425 pFields[i]->GetExp(n)->PhysDeriv(
426 i, temperature + phys_offset,
427 auxArray = Dtemperature[i] + phys_offset);
430 pFields[0]->ExtractTracePhys(Dtemperature[i], traceDtemperature[i]);
433 for(i = 0; i < nDimensions; ++i)
436 &m_traceNormals[i][0], 1,
437 &traceDtemperature[i][0], 1,
441 &traceFieldsAdded[6][0], 1,
443 &traceFieldsAdded[6][0], 1);
455 for (i = 0; i < nDimensions; ++ i)
461 for (i = 0; i < nDimensions; ++ i)
463 for (n = 0; n < nElements; n++)
465 phys_offset = pFields[0]->GetPhys_Offset(n);
467 pFields[i]->GetExp(n)->PhysDeriv(
468 i, pressure + phys_offset,
469 auxArray = Dpressure[i] + phys_offset);
472 pFields[0]->ExtractTracePhys(Dpressure[i], traceDpressure[i]);
476 for(i = 0; i < nDimensions; ++i)
479 &m_traceTangents[i][0], 1,
480 &traceDpressure[i][0], 1,
484 &traceFieldsAdded[7][0], 1,
486 &traceFieldsAdded[7][0], 1);
491 &traceDpressure[0][0], 1,
492 &traceFieldsAdded[8][0], 1);
496 &traceDpressure[1][0], 1,
497 &traceFieldsAdded[9][0], 1);
509 for (i = 0; i < nDimensions; ++ i)
515 Vmath::Vdiv(nSolutionPts, uFields[i+1], 1, uFields[0], 1,
518 for (j = 0; j < nDimensions; ++j)
525 for (i = 0; i < nDimensions; ++i)
527 for (j = 0; j < nDimensions; ++j)
529 for (n = 0; n < nElements; n++)
531 phys_offset = pFields[0]->GetPhys_Offset(n);
533 pFields[i]->GetExp(n)->PhysDeriv(
534 j, velocity[i] + phys_offset,
535 auxArray = Dvelocity[i][j] + phys_offset);
539 pFields[0]->ExtractTracePhys(Dvelocity[i][j], traceDvelocity[i][j]);
544 &traceDvelocity[0][0][0], 1,
545 &traceFieldsAdded[10][0], 1);
547 &traceDvelocity[0][1][0], 1,
548 &traceFieldsAdded[11][0], 1);
550 &traceDvelocity[1][0][0], 1,
551 &traceFieldsAdded[12][0], 1);
553 &traceDvelocity[1][1][0], 1,
554 &traceFieldsAdded[13][0], 1);
572 if (m_ViscosityType ==
"Variable")
575 NekDouble T_star = m_pInf / (m_rhoInf * m_gasConstant);
578 for (
int i = 0; i < nSolutionPts; ++i)
580 ratio = temperature[i] / T_star;
581 mu[i] = mu_star * ratio * sqrt(ratio) *
582 (T_star + 110.0) / (temperature[i] + 110.0);
595 Vmath::Smul(nSolutionPts, 2.0, &mu[0], 1, &mu2[0], 1);
599 &Dvelocity[0][0][0], 1, &divVel[0], 1);
601 &Dvelocity[1][1][0], 1, &divVel[0], 1);
604 Vmath::Smul(nSolutionPts, lambda, &divVel[0], 1, &divVel[0], 1);
605 Vmath::Vmul(nSolutionPts, &mu[0], 1, &divVel[0], 1, &divVel[0], 1);
609 for (j = 0; j < m_spacedim; ++j)
614 Vmath::Vmul(nSolutionPts, &mu2[0], 1, &Dvelocity[j][j][0], 1,
617 Vmath::Vadd(nSolutionPts, &temp[j][0], 1, &divVel[0], 1, &Sgg[j][0], 1);
625 &Dvelocity[1][0][0], 1, &Sxy[0], 1);
628 Vmath::Vmul(nSolutionPts, &mu[0], 1, &Sxy[0], 1, &Sxy[0], 1);
630 pFields[0]->ExtractTracePhys(Sgg[0], traceFieldsAdded[14]);
631 pFields[0]->ExtractTracePhys(Sgg[1], traceFieldsAdded[15]);
632 pFields[0]->ExtractTracePhys(Sxy, traceFieldsAdded[16]);
647 Vmath::Vcopy(nTracePts, &m_traceNormals[0][0], 1, &cosTeta[0], 1);
650 Vmath::Vcopy(nTracePts, &m_traceNormals[1][0], 1, &sinTeta[0], 1);
653 Vmath::Vsub(nTracePts, &traceFieldsAdded[14][0], 1,
654 &traceFieldsAdded[15][0], 1, &sigma_diff[0], 1);
657 Vmath::Vmul(nTracePts, &cosTeta[0], 1, &sinTeta[0], 1, &tmpTeta[0], 1);
658 Vmath::Smul(nTracePts, 2.0, &tmpTeta[0], 1, &sin2Teta[0], 1);
661 Vmath::Vmul(nTracePts, &cosTeta[0], 1, &cosTeta[0], 1, &cos2Teta[0], 1);
662 Vmath::Vmul(nTracePts, &sinTeta[0], 1, &sinTeta[0], 1, &tmpTeta[0], 1);
663 Vmath::Vsub(nTracePts, &cos2Teta[0], 1, &tmpTeta[0], 1, &cos2Teta[0], 1);
666 Vmath::Smul(nTracePts, -0.5, &sigma_diff[0], 1, &sigma_diff[0], 1);
667 Vmath::Vmul(nTracePts, &sigma_diff[0], 1, &sin2Teta[0], 1, &tau_t[0], 1);
668 Vmath::Vmul(nTracePts, &traceFieldsAdded[16][0], 1, &cos2Teta[0], 1,
670 Vmath::Vadd(nTracePts, &tau_t[0], 1, &tmpTeta[0], 1, &tau_t[0], 1);
672 Vmath::Vcopy(nTracePts, &tau_t[0], 1, &traceFieldsAdded[17][0], 1);
678 pFields[0]->ExtractTracePhys(mu, traceFieldsAdded[18]);
688 Vmath::Vdiv (nSolutionPts, pressure, 1, uFields[0], 1, soundspeed, 1);
689 Vmath::Smul (nSolutionPts, gamma, soundspeed, 1, soundspeed, 1);
690 Vmath::Vsqrt(nSolutionPts, soundspeed, 1, soundspeed, 1);
695 for (
int i = 0; i < m_spacedim; ++i)
697 Vmath::Vvtvp(nSolutionPts, uFields[i + 1], 1, uFields[i + 1], 1,
701 Vmath::Vdiv(nSolutionPts, mach, 1, uFields[0], 1, mach, 1);
702 Vmath::Vdiv(nSolutionPts, mach, 1, uFields[0], 1, mach, 1);
704 Vmath::Vdiv(nSolutionPts, mach, 1, soundspeed, 1, mach, 1);
706 pFields[0]->ExtractTracePhys(mach, traceFieldsAdded[19]);
711 if (pFields[0]->GetBndCondExpansions().num_elements())
715 nBndRegions = pFields[0]->GetBndCondExpansions().num_elements();
716 for (b = 0; b < nBndRegions; ++b)
718 nBndEdges = pFields[0]->GetBndCondExpansions()[b]->GetExpSize();
719 for (e = 0; e < nBndEdges; ++e)
721 nBndEdgePts = pFields[0]->
722 GetBndCondExpansions()[b]->GetExp(e)->GetNumPoints(0);
724 id2 = pFields[0]->GetTrace()->
725 GetPhys_Offset(pFields[0]->GetTraceMap()->
726 GetBndCondTraceToGlobalTraceMap(cnt++));
728 if (pFields[0]->GetBndConditions()[b]->
729 GetUserDefined() ==
"WallViscous" ||
730 pFields[0]->GetBndConditions()[b]->
731 GetUserDefined() ==
"WallAdiabatic" ||
732 pFields[0]->GetBndConditions()[b]->
733 GetUserDefined() ==
"Wall")
751 if (pFields[0]->GetBndCondExpansions().num_elements())
753 for (j = 0; j < nfields; ++j)
757 nBndRegions = pFields[j]->GetBndCondExpansions().num_elements();
758 for (b = 0; b < nBndRegions; ++b)
760 nBndEdges = pFields[j]->GetBndCondExpansions()[b]->GetExpSize();
761 for (e = 0; e < nBndEdges; ++e)
763 nBndEdgePts = pFields[j]->
764 GetBndCondExpansions()[b]->GetExp(e)->GetNumPoints(0);
766 id2 = pFields[j]->GetTrace()->
767 GetPhys_Offset(pFields[j]->GetTraceMap()->
768 GetBndCondTraceToGlobalTraceMap(cnt++));
770 if (pFields[j]->GetBndConditions()[b]->
771 GetUserDefined() ==
"WallViscous" ||
772 pFields[j]->GetBndConditions()[b]->
773 GetUserDefined() ==
"WallAdiabatic" ||
774 pFields[j]->GetBndConditions()[b]->
775 GetUserDefined() ==
"Wall")
778 &surfaceFields[j][id1], 1);
788 if (pFields[0]->GetBndCondExpansions().num_elements())
790 for (j = 0; j < nfieldsAdded; ++j)
794 nBndRegions = pFields[0]->GetBndCondExpansions().num_elements();
795 for (b = 0; b < nBndRegions; ++b)
797 nBndEdges = pFields[0]->GetBndCondExpansions()[b]->GetExpSize();
798 for (e = 0; e < nBndEdges; ++e)
800 nBndEdgePts = pFields[0]->
801 GetBndCondExpansions()[b]->GetExp(e)->GetNumPoints(0);
803 id2 = pFields[0]->GetTrace()->
804 GetPhys_Offset(pFields[0]->GetTraceMap()->
805 GetBndCondTraceToGlobalTraceMap(cnt++));
807 if (pFields[0]->GetBndConditions()[b]->
808 GetUserDefined() ==
"WallViscous" ||
809 pFields[0]->GetBndConditions()[b]->
810 GetUserDefined() ==
"WallAdiabatic" ||
811 pFields[0]->GetBndConditions()[b]->
812 GetUserDefined() ==
"Wall")
815 &surfaceFieldsAdded[j][id1], 1);
826 std::string vEquation = vSession->GetSolverInfo(
"EQType");
829 BndExp = pFields[0]->GetBndCondExpansions();
840 for(
int i = 0; i < BndExp[0]->GetExpSize(); ++i)
863 for(
int j = 0; j < nbc; ++j)
866 nxOnBnd[j] = surfaceFieldsAdded[0][GlobalIndex];
867 nyOnBnd[j] = surfaceFieldsAdded[1][GlobalIndex];
868 txOnBnd[j] = surfaceFieldsAdded[2][GlobalIndex];
869 tyOnBnd[j] = surfaceFieldsAdded[3][GlobalIndex];
871 PressurOnBnd[j] = surfaceFieldsAdded[4][GlobalIndex];
873 if (vEquation ==
"NavierStokesCFE")
875 ShearStressOnBnd[j] = surfaceFieldsAdded[17][GlobalIndex];
892 Vmath::Vmul(nbc,PressurOnBnd,1,nxOnBnd,1, drag_p,1);
893 Vmath::Vmul(nbc,PressurOnBnd,1,nyOnBnd,1, lift_p,1);
898 Fxp += bc->Integral(drag_p);
899 Fyp += bc->Integral(lift_p);
901 if (vEquation ==
"NavierStokesCFE")
903 Vmath::Vmul(nbc,ShearStressOnBnd,1,txOnBnd,1, drag_v,1);
904 Vmath::Vmul(nbc,ShearStressOnBnd,1,tyOnBnd,1, lift_v,1);
909 Fxv += bc->Integral(drag_v);
910 Fyv += bc->Integral(lift_v);
913 Sref += bc->Integral(Unity);
917 cout <<
"\n Sref = " << Sref << endl;
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;
936 outfile.open(fname.c_str());
937 outfile <<
"% x[m] " <<
" \t"
944 <<
"rho[kg/m^3] " <<
" \t"
945 <<
"rhou[kg/(m^2 s)] " <<
" \t"
946 <<
"rhov[kg/(m^2 s)] " <<
" \t"
950 <<
"dT/dn[k/m] " <<
" \t"
951 <<
"dp/dT[Pa/m] " <<
" \t"
952 <<
"dp/dx[Pa/m] " <<
" \t"
953 <<
"dp/dy[Pa/m] " <<
" \t"
954 <<
"du/dx[s^-1] " <<
" \t"
955 <<
"du/dy[s^-1] " <<
" \t"
956 <<
"dv/dx[s^-1] " <<
" \t"
957 <<
"dv/dy[s^-1] " <<
" \t"
958 <<
"tau_xx[Pa] " <<
" \t"
959 <<
"tau_yy[Pa] " <<
" \t"
960 <<
"tau_xy[Pa] " <<
" \t"
961 <<
"tau_t[Pa] " <<
" \t"
962 <<
"mu[Pa s] " <<
" \t"
966 for (i = 0; i < id1; ++i)
968 outfile << scientific
971 << surfaceX[i] <<
" \t "
972 << surfaceY[i] <<
" \t "
973 << surfaceZ[i] <<
" \t "
974 << surfaceFieldsAdded[0][i] <<
" \t "
975 << surfaceFieldsAdded[1][i] <<
" \t "
976 << surfaceFieldsAdded[2][i] <<
" \t "
977 << surfaceFieldsAdded[3][i] <<
" \t "
978 << surfaceFields[0][i] <<
" \t "
979 << surfaceFields[1][i] <<
" \t "
980 << surfaceFields[2][i] <<
" \t "
981 << surfaceFields[3][i] <<
" \t "
982 << surfaceFieldsAdded[4][i] <<
" \t "
983 << surfaceFieldsAdded[5][i] <<
" \t "
984 << surfaceFieldsAdded[6][i] <<
" \t "
985 << surfaceFieldsAdded[7][i] <<
" \t "
986 << surfaceFieldsAdded[8][i] <<
" \t "
987 << surfaceFieldsAdded[9][i] <<
" \t "
988 << surfaceFieldsAdded[10][i] <<
" \t "
989 << surfaceFieldsAdded[11][i] <<
" \t "
990 << surfaceFieldsAdded[12][i] <<
" \t "
991 << surfaceFieldsAdded[13][i] <<
" \t "
992 << surfaceFieldsAdded[14][i] <<
" \t "
993 << surfaceFieldsAdded[15][i] <<
" \t "
994 << surfaceFieldsAdded[16][i] <<
" \t "
995 << surfaceFieldsAdded[17][i] <<
" \t "
996 << surfaceFieldsAdded[18][i] <<
" \t "
997 << surfaceFieldsAdded[19][i] <<
" \t "
1001 outfile << endl << endl;
#define ASSERTL0(condition, msg)
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
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
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.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
void Import(const std::string &infilename, std::vector< FieldDefinitionsSharedPtr > &fielddefs, std::vector< std::vector< NekDouble > > &fielddata, FieldMetaDataMap &fieldinfomap, const Array< OneD, int > ElementiDs)
Imports an FLD file.
1D Evenly-spaced points using Lagrange polynomial
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
boost::shared_ptr< StdExpansion1D > StdExpansion1DSharedPtr
void Neg(int n, T *x, const int incx)
Negate x = -x.
boost::shared_ptr< ExpList2D > ExpList2DSharedPtr
Shared pointer to an ExpList2D object.
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.
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
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.
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.