43 NekDouble Extrapolate::StifflyStable_Betaq_Coeffs[3][3] = {
44 { 1.0, 0.0, 0.0},{ 2.0, -1.0, 0.0},{ 3.0, -3.0, 1.0}};
45 NekDouble Extrapolate::StifflyStable_Alpha_Coeffs[3][3] = {
46 { 1.0, 0.0, 0.0},{ 2.0, -0.5, 0.0},{ 3.0, -1.5, 1.0/3.0}};
47 NekDouble Extrapolate::StifflyStable_Gamma0_Coeffs[3] = {
55 Loki::SingleThreaded > Type;
56 return Type::Instance();
59 Extrapolate::Extrapolate(
65 : m_session(pSession),
67 m_pressure(pPressure),
69 m_advObject(advObject)
81 "StandardExtrapolate",
"StandardExtrapolate");
102 accelerationTerm, 1);
104 for(
int i = 0; i < acc_order; i++)
110 accelerationTerm, 1);
131 for(
int n = 0; n < nint-1; ++n)
146 for(cnt = n = 0; n <
m_PBndConds.num_elements(); ++n)
149 if(boost::iequals(
m_PBndConds[n]->GetUserDefined(),
"H"))
172 int pindex=N.num_elements();
186 for(
int j = 0 ; j <
m_HBCdata.num_elements() ; j++)
198 Velocity[i] = fields[i] +
m_HBCdata[j].m_physOffset;
199 Advection[i] = N[i] +
m_HBCdata[j].m_physOffset;
205 Velocity[2] = fields[2] +
m_HBCdata[j].m_assPhysOffset;
222 switch(
m_fields[pindex]->GetExpType())
227 elmt->GetEdgePhysVals(
m_HBCdata[j].m_elmtTraceID,Pbc,Q[0],BndValues[0]);
228 elmt->GetEdgePhysVals(
m_HBCdata[j].m_elmtTraceID,Pbc,Q[1],BndValues[1]);
229 Pbc->NormVectorIProductWRTBase(BndValues[0],BndValues[1],Pvals);
236 Pvals[0] = -1.0*Q[0][0];
238 else if (
m_HBCdata[j].m_elmtTraceID == 1)
240 Pvals[0] = Q[0][
m_HBCdata[j].m_ptsInElmt-1];
245 "In the 3D homogeneous 2D approach BCs edge "
246 "ID can be just 0 or 1 ");
252 elmt->GetFacePhysVals(
m_HBCdata[j].m_elmtTraceID,Pbc,Q[0],BndValues[0]);
253 elmt->GetFacePhysVals(
m_HBCdata[j].m_elmtTraceID,Pbc,Q[1],BndValues[1]);
254 elmt->GetFacePhysVals(
m_HBCdata[j].m_elmtTraceID,Pbc,Q[2],BndValues[2]);
255 Pbc->NormVectorIProductWRTBase(BndValues[0],BndValues[1],BndValues[2],Pvals);
259 ASSERTL0(0,
"Dimension not supported");
274 static bool init =
true;
275 static bool noHOBC =
false;
286 for(
int n = 0; n <
m_PBndConds.num_elements(); ++n)
288 if(boost::iequals(
m_PBndConds[n]->GetUserDefined(),
"HOutflow"))
290 totbndpts +=
m_PBndExp[n]->GetTotPoints();
336 int num_planes = planes.num_elements();
338 for(
int n = 0; n <
m_PBndConds.num_elements(); ++n)
340 int exp_size =
m_PBndExp[n]->GetExpSize();
344 for(
int n = 0; n <
m_PBndConds.num_elements(); ++n)
384 int nbc,cnt,cnt_start;
394 for(
int n = 0; n <
m_PBndConds.num_elements(); ++n)
397 if(boost::iequals(
m_PBndConds[n]->GetUserDefined(),
"HOutflow"))
399 for(
int i = 0; i <
m_PBndExp[n]->GetExpSize(); ++i,cnt++)
401 cnt = max(cnt,
m_PBndExp[n]->GetTotPoints());
412 m_session->LoadParameter(
"U0_HighOrderBC",U0,1.0);
413 m_session->LoadParameter(
"Delta_HighOrderBC",delta,1/20.0);
416 for(
int n = 0; n <
m_PBndConds.num_elements(); ++n)
421 if(boost::iequals(
m_PBndConds[n]->GetUserDefined(),
"HOutflow"))
427 for(
int i = 0; i <
m_PBndExp[n]->GetExpSize(); ++i, cnt++)
435 int offset =
m_fields[0]->GetPhys_Offset(elmtid);
440 int nq = elmt->GetTotPoints();
441 int nbc =
m_PBndExp[n]->GetExp(i)->GetTotPoints();
450 elmt->GetTracePhysVals(boundary,Bc,BndElmt[j],
468 for(
int i = 0; i <
m_PBndExp[n]->GetExpSize(); ++i, cnt++)
473 int nbc =
m_PBndExp[n]->GetExp(i)->GetTotPoints();
476 normDotu(nbc,0.0), utot(nbc,0.0);
478 normals=elmt->GetSurfaceNormal(boundary);
513 normDotu, 1, normDotu, 1);
516 int Offset =
m_PBndExp[n]->GetPhys_Offset(i);
518 for(
int k = 0; k < nbc; ++k)
522 NekDouble fac = 0.5*(1.0-tanh(normDotu[k]/(U0*delta)));
531 UBndExp[0][n]->HomogeneousFwdTrans(
544 UBndExp[j][n]->HomogeneousFwdTrans(
545 UBndExp[j][n]->UpdatePhys(),
552 for(
int i = 0; i <
m_PBndExp[n]->GetExpSize(); ++i,cnt++)
560 int offset =
m_fields[0]->GetPhys_Offset(elmtid);
563 int nq = elmt->GetTotPoints();
573 int nbc =
m_PBndExp[n]->GetExp(i)->GetTotPoints();
579 normals=elmt->GetSurfaceNormal(boundary);
588 BndElmt[j], gradtmp);
589 elmt->GetTracePhysVals(boundary, Bc, gradtmp,
592 nGradu[j], 1, nGradu[j],1);
609 int p_offset =
m_PBndExp[n]->GetPhys_Offset(i);
611 for(
int k = 0; k < nbc; ++k)
617 ptmp[k] = kinvis * ptmp[k]
622 int u_offset = UBndExp[0][n]->GetPhys_Offset(i);
626 for(
int k = 0; k < nbc; ++k)
628 ubc[j][k + u_offset] = (1.0 / kinvis)
636 for(
int k = 0; k < nbc; ++k)
638 ubc[
m_bnd_dim][k + u_offset] = (1.0 / kinvis)
642 u_offset = UBndExp[
m_bnd_dim][n]->GetPhys_Offset(i);
643 UBCvals = UBndExp[
m_bnd_dim][n]->UpdateCoeffs()
644 + UBndExp[
m_bnd_dim][n]->GetCoeff_Offset(i);
645 Bc->IProductWRTBase(ubc[m_bnd_dim] + u_offset, UBCvals);
655 elmt->GetTracePhysVals(boundary,Bc,BndElmt[j],
693 normDotu, 1, normDotu, 1);
701 for(
int k = 0; k < nbc; ++k)
703 NekDouble fac = 0.5*(1.0-tanh(normDotu[k]/(U0*delta)));
708 ptmp[k] = kinvis * ptmp[k] - 0.5 * utot[k] * fac
712 int u_offset = UBndExp[0][n]->GetPhys_Offset(i);
716 UBCvals = UBndExp[j][n]->GetPhys()
717 + UBndExp[j][n]->GetPhys_Offset(i);
719 for(
int k = 0; k < nbc; ++k)
721 NekDouble fac = 0.5 * (1.0 - tanh(normDotu[k]
723 ubc[j][k + u_offset] = (1.0 / kinvis)
724 * (UBCvals[k] + 0.5 * utot[k] * fac
733 m_PBndExp[n]->GetExp(i)->FwdTrans(ptmp,PBCvals);
741 if(boost::iequals(UBndConds[j][n]->GetUserDefined(),
"HOutflow"))
744 for(
int i = 0; i < UBndExp[0][n]->GetExpSize();
750 (UBndExp[0][n]->GetExp(i));
752 nbc = UBndExp[0][n]->GetExp(i)->GetTotPoints();
760 normals = elmt->GetSurfaceNormal(boundary);
765 Pbc->BwdTrans(PBCvals,pb);
767 int u_offset = UBndExp[j][n]->GetPhys_Offset(i);
769 for(
int k = 0; k < nbc; ++k)
771 ub[k] = ubc[j][k + u_offset]
772 + pb[k] * normals[j][k] / kinvis;
775 UBCvals = UBndExp[j][n]->UpdateCoeffs()
776 + UBndExp[j][n]->GetCoeff_Offset(i);
777 Bc->IProductWRTBase(ub,UBCvals);
814 elmt->PhysDeriv(Dummy,Q[1],Q[0]);
837 Vel[0], 1, Dummy1, 1);
846 Vel[1], 1, Dummy2, 1);
897 elmt->PhysDeriv(Vel[0], Dummy, Uy, Uz);
898 elmt->PhysDeriv(Vel[1], Vx, Dummy, Vz);
899 elmt->PhysDeriv(Vel[2], Wx, Wy, Dummy);
905 elmt->PhysDeriv(Q[0], Dummy, Wy, Vx);
906 elmt->PhysDeriv(Q[1], Wx, Dummy, Uz);
907 elmt->PhysDeriv(Q[2], Vz, Uy, Dummy);
915 ASSERTL0(0,
"Dimension not supported");
928 int pindex=
m_fields.num_elements()-1;
940 for(n = 0; n <
m_HBCdata.num_elements(); ++n)
951 Velmt = Vel[j] +
m_HBCdata[n].m_physOffset;
952 elmt->GetTracePhysVals(
m_HBCdata[n].m_elmtTraceID,Pbc,Velmt,velbc[j]);
955 IProdVnTmp = IProdVn +
m_HBCdata[n].m_coeffOffset;
958 Pbc->NormVectorIProductWRTBase(velbc,IProdVnTmp);
986 for(n = 0; n <
m_HBCdata.num_elements(); ++n)
994 VelBndExp[j][bndid]->GetExp(elmtid)->BwdTrans(VelBndExp[j][bndid]->GetCoeffs() +
995 VelBndExp[j][bndid]->GetCoeff_Offset(elmtid),
999 IProdVnTmp = IProdVn +
m_HBCdata[n].m_coeffOffset;
1002 Pbc->NormVectorIProductWRTBase(velbc,IProdVnTmp);
1014 int nlevels = input.num_elements();
1018 tmp = input[nlevels-1];
1020 for(
int n = nlevels-1; n > 0; --n)
1022 input[n] = input[n-1];
1035 int pindex=
m_fields.num_elements()-1;
1049 for(cnt = n = 0; n <
m_PBndConds.num_elements(); ++n)
1051 for(
int i = 0; i <
m_PBndExp[n]->GetExpSize(); ++i)
1054 m_PBndExp[n]->GetExp(i)->GetTotPoints());
1066 for(cnt = n = 0; n <
m_PBndConds.num_elements(); ++n)
1069 if(boost::iequals(
m_PBndConds[n]->GetUserDefined(),
"H"))
1072 HBCnumber +=
m_PBndExp[n]->GetExpSize();
1076 int checkHBC = HBCnumber;
1118 ASSERTL0(0,
"Dimension not supported");
1131 int coeff_count = 0;
1135 for(
int n = 0 ; n <
m_PBndConds.num_elements(); ++n)
1139 if(boost::iequals(
m_PBndConds[n]->GetUserDefined(),
"H"))
1141 for(
int i = 0; i < exp_size; ++i,cnt++)
1145 m_HBCdata[j].m_ptsInElmt = elmt->GetTotPoints();
1150 m_HBCdata[j].m_coeffOffset = coeff_count;
1151 coeff_count += elmt->GetTraceNcoeffs(
m_HBCdata[j].m_elmtTraceID);
1167 int num_planes = planes.num_elements();
1168 int num_elm_per_plane = (
m_pressure->GetExpSize())/num_planes;
1173 int exp_size, exp_size_per_plane;
1179 m_session->MatchSolverInfo(
"ModeType",
"SingleMode",
1181 m_session->MatchSolverInfo(
"ModeType",
"HalfMode",
1183 m_session->MatchSolverInfo(
"ModeType",
"MultipleModes",
1188 if(
m_session->DefinesSolverInfo(
"ModeType"))
1204 ASSERTL0(
false,
"SolverInfo ModeType not valid");
1212 int coeff_count = 0;
1214 for(n = 0, j= 0, cnt = 0; n <
m_PBndConds.num_elements(); ++n)
1217 exp_size_per_plane = exp_size/num_planes;
1219 if(boost::iequals(
m_PBndConds[n]->GetUserDefined(),
"H"))
1221 for(k = 0; k < num_planes; k++)
1224 for(i = 0; i < exp_size_per_plane; ++i, ++j, ++cnt)
1228 m_HBCdata[j].m_ptsInElmt = elmt->GetTotPoints();
1230 m_HBCdata[j].m_bndElmtID = i+k*exp_size_per_plane;
1233 m_HBCdata[j].m_coeffOffset = coeff_count;
1234 coeff_count += elmt->GetEdgeNcoeffs(
m_HBCdata[j].m_elmtTraceID);
1258 assElmtID =
m_HBCdata[j].m_globalElmtID;
1263 assElmtID =
m_HBCdata[j].m_globalElmtID + num_elm_per_plane;
1268 assElmtID =
m_HBCdata[j].m_globalElmtID - num_elm_per_plane;
1287 int exp_size, exp_size_per_line;
1294 for(
int n = 0 ; n <
m_PBndConds.num_elements(); ++n)
1298 exp_size_per_line = exp_size/(m_npointsZ*
m_npointsY);
1300 if(boost::iequals(
m_PBndConds[n]->GetUserDefined(),
"H"))
1302 for(
int i = 0; i < exp_size_per_line; ++i,cnt++)
1307 m_HBCdata[j].m_ptsInElmt = elmt->GetTotPoints();
1309 m_HBCdata[j].m_bndElmtID = i+(k1*m_npointsY+k2)*exp_size_per_line;
1318 cnt += exp_size_per_line;
1326 ASSERTL0(0,
"Dimension not supported");
1340 int n_points_0 =
m_fields[0]->GetExp(0)->GetTotPoints();
1341 int n_element =
m_fields[0]->GetExpSize();
1342 int nvel = inarray.num_elements();
1352 for (
int i = 0; i < nvel; ++i)
1360 for (
int el = 0; el < n_element; ++el)
1362 int n_points =
m_fields[0]->GetExp(el)->GetTotPoints();
1363 ptsKeys =
m_fields[0]->GetExp(el)->GetPointsKeys();
1366 if(n_points != n_points_0)
1368 for (
int j = 0; j < nvel; ++j)
1372 n_points_0 = n_points;
1376 m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo()->GetDerivFactors(ptsKeys);
1378 if (
m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo()->GetGtype()
1381 for (
int i = 0; i < n_points; i++)
1383 stdVelocity[0][i] = gmat[0][i]*inarray[0][i+cnt]
1384 + gmat[2][i]*inarray[1][i+cnt];
1386 stdVelocity[1][i] = gmat[1][i]*inarray[0][i+cnt]
1387 + gmat[3][i]*inarray[1][i+cnt];
1392 for (
int i = 0; i < n_points; i++)
1394 stdVelocity[0][i] = gmat[0][0]*inarray[0][i+cnt]
1395 + gmat[2][0]*inarray[1][i+cnt];
1397 stdVelocity[1][i] = gmat[1][0]*inarray[0][i+cnt]
1398 + gmat[3][0]*inarray[1][i+cnt];
1405 for (
int i = 0; i < n_points; i++)
1407 pntVelocity = stdVelocity[0][i]*stdVelocity[0][i]
1408 + stdVelocity[1][i]*stdVelocity[1][i];
1410 if (pntVelocity>maxV[el])
1412 maxV[el] = pntVelocity;
1415 maxV[el] = sqrt(maxV[el]);
1421 for (
int el = 0; el < n_element; ++el)
1424 int n_points =
m_fields[0]->GetExp(el)->GetTotPoints();
1425 ptsKeys =
m_fields[0]->GetExp(el)->GetPointsKeys();
1428 if(n_points != n_points_0)
1430 for (
int j = 0; j < nvel; ++j)
1434 n_points_0 = n_points;
1438 m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo()->GetDerivFactors(ptsKeys);
1440 if (
m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo()->GetGtype()
1443 for (
int i = 0; i < n_points; i++)
1445 stdVelocity[0][i] = gmat[0][i]*inarray[0][i+cnt]
1446 + gmat[3][i]*inarray[1][i+cnt]
1447 + gmat[6][i]*inarray[2][i+cnt];
1449 stdVelocity[1][i] = gmat[1][i]*inarray[0][i+cnt]
1450 + gmat[4][i]*inarray[1][i+cnt]
1451 + gmat[7][i]*inarray[2][i+cnt];
1453 stdVelocity[2][i] = gmat[2][i]*inarray[0][i+cnt]
1454 + gmat[5][i]*inarray[1][i+cnt]
1455 + gmat[8][i]*inarray[2][i+cnt];
1460 for (
int i = 0; i < n_points; i++)
1462 stdVelocity[0][i] = gmat[0][0]*inarray[0][i+cnt]
1463 + gmat[3][0]*inarray[1][i+cnt]
1464 + gmat[6][0]*inarray[2][i+cnt];
1466 stdVelocity[1][i] = gmat[1][0]*inarray[0][i+cnt]
1467 + gmat[4][0]*inarray[1][i+cnt]
1468 + gmat[7][0]*inarray[2][i+cnt];
1470 stdVelocity[2][i] = gmat[2][0]*inarray[0][i+cnt]
1471 + gmat[5][0]*inarray[1][i+cnt]
1472 + gmat[8][0]*inarray[2][i+cnt];
1478 for (
int i = 0; i < n_points; i++)
1480 pntVelocity = stdVelocity[0][i]*stdVelocity[0][i]
1481 + stdVelocity[1][i]*stdVelocity[1][i]
1482 + stdVelocity[2][i]*stdVelocity[2][i];
1484 if (pntVelocity > maxV[el])
1486 maxV[el] = pntVelocity;
1490 maxV[el] = sqrt(maxV[el]);
1513 int nPts = newarray.num_elements();
1521 oldarrays[nint-1], 1,
1524 for(
int n = 0; n < nint-1; ++n)
1527 oldarrays[n],1,outarray,1,
#define ASSERTL0(condition, msg)
std::vector< PointsKey > PointsKeyVector
#define sign(a, b)
return the sign(b)*a
ExtrapolateFactory & GetExtrapolateFactory()
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
boost::shared_ptr< Advection > AdvectionSharedPtr
A shared pointer to an Advection object.
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
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
The base class for all shapes.
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
LibUtilities::NekFactory< std::string, Extrapolate, const LibUtilities::SessionReaderSharedPtr &, Array< OneD, MultiRegions::ExpListSharedPtr > &, MultiRegions::ExpListSharedPtr &, const Array< OneD, int > &, const SolverUtils::AdvectionSharedPtr & > ExtrapolateFactory
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList 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.
MultiRegions::Direction const DirCartesianMap[]
static std::string RegisterDefaultSolverInfo(const std::string &pName, const std::string &pValue)
Registers the default string value of a solver info property.
void Zero(int n, T *x, const int incx)
Zero vector.
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Geometry is curved or has non-constant factors.
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.
Defines a callback function which evaluates the flux vector.
Provides a generic Factory class.