35 #include <boost/core/ignore_unused.hpp>
64 m_session->LoadParameter(
"AdvectionJacFlag", ntmp, 1);
67 m_session->LoadParameter(
"ViscousJacFlag", ntmp, 1);
80 std::string SolverType =
"Newton";
81 if (
m_session->DefinesSolverInfo(
"NonlinSysIterSolver"))
83 SolverType =
m_session->GetSolverInfo(
"NonlinSysIterSolver");
86 "NekNonlinSys '" + SolverType +
"' is not defined.\n");
108 m_fields[0]->GetLocTraceToTraceMap();
110 locTraceToTraceMap->CalcLocTracePhysToTraceIDMap(
m_fields[0]->GetTrace(),
112 for (
int i = 1; i < nvariables; i++)
114 m_fields[i]->GetLocTraceToTraceMap()->SetLocTracePhysToTraceIDMap(
115 locTraceToTraceMap->GetLocTracephysToTraceIDMap());
154 boost::ignore_unused(flag);
155 unsigned int nvariables =
m_fields.size();
156 unsigned int npoints =
m_fields[0]->GetNcoeffs();
160 for (
int i = 0; i < nvariables; ++i)
162 int offset = i * npoints;
163 in2D[i] = inarray + offset;
164 out2D[i] = out + offset;
165 source2D[i] = source + offset;
180 unsigned int nvariable = inarray.size();
181 unsigned int ncoeffs = inarray[nvariable - 1].size();
182 unsigned int npoints =
m_fields[0]->GetNpoints();
186 for (
int i = 0; i < nvariable; ++i)
189 m_fields[i]->BwdTrans(inarray[i], inpnts[i]);
202 for (
int i = 0; i < nvariable; ++i)
210 for (
int i = 0; i < nvariable; ++i)
212 Vmath::Vsub(ncoeffs, out[i], 1, source[i], 1, out[i], 1);
227 int nvariables = inarray.size();
241 for (
int i = 0; i < nvariables; ++i)
245 m_fields[i]->GetFwdBwdTracePhys(inarray[i], Fwd[i], Bwd[i]);
256 for (
int i = 0; i < nvariables; ++i)
269 x->ApplyCoeff(
m_fields, inarray, outarray, time);
274 int nElements =
m_fields[0]->GetExpSize();
283 for (
int n = 0; n < nElements; ++n)
285 nq =
m_fields[0]->GetExp(n)->GetNcoeffs();
286 offset =
m_fields[0]->GetCoeff_Offset(n);
288 for (
int i = 0; i < nvariables; ++i)
291 tmp = outarray[i] + offset, 1);
306 int nvariables = inarray.size();
318 unsigned int nvariables = inpnts.size();
319 unsigned int ncoeffs =
m_fields[0]->GetNcoeffs();
320 unsigned int ntotal = nvariables * ncoeffs;
326 for (
int i = 0; i < nvariables; ++i)
328 int noffset = i * ncoeffs;
329 tmpArray = inarray + noffset;
330 m_fields[i]->FwdTrans(inpnts[i], tmpArray);
335 for (
int i = 0; i < nvariables; ++i)
337 int noffset = i * ncoeffs;
338 tmpArray = out + noffset;
339 m_fields[i]->BwdTrans(tmpArray, outpnt[i]);
348 boost::ignore_unused(inpnts);
353 unsigned int ntotal = inarray.size();
363 m_nonlinsol->v_SetupNekNonlinSystem(ntotal, inarray, inarray, 0);
375 unsigned int nvariables =
m_fields.size();
376 unsigned int ntotal = inarray.size();
377 unsigned int npoints = ntotal / nvariables;
379 unsigned int nTotalGlobal = ntotal;
381 unsigned int nTotalDOF = nTotalGlobal / nvariables;
387 for (
int i = 0; i < nvariables; ++i)
389 int offset = i * npoints;
391 Vmath::Dot(npoints, inarray + offset, inarray + offset);
395 for (
int i = 0; i < nvariables; ++i)
400 for (
int i = 2; i < nvariables - 1; ++i)
404 for (
int i = 2; i < nvariables - 1; ++i)
409 for (
int i = 0; i < nvariables; ++i)
415 for (
int i = 0; i < nvariables; ++i)
439 for (
int i = 0; i < nvariables; i++)
467 template <
typename DataType,
typename TypeNekBlkMatSharedPtr>
475 if (StdMatDataDBB.size() == 0)
480 int nSpaceDim =
m_graph->GetSpaceDimension();
481 int nvariable = inarray.size();
482 int npoints =
m_fields[0]->GetTotPoints();
483 int nVar2 = nvariable * nvariable;
484 std::shared_ptr<LocalRegions::ExpansionVector> expvect =
486 int nTotElmt = (*expvect).size();
497 for (
int i = 0; i < 3; i++)
501 normal3D[0][0] = 1.0;
502 normal3D[1][1] = 1.0;
503 normal3D[2][2] = 1.0;
509 nvariable - 1, nvariable, 0.0);
546 int nElmtCoefOld = -1;
547 for (
int ne = 0; ne < nTotElmt; ne++)
549 int nElmtCoef = (*expvect)[ne]->GetNcoeffs();
550 int nElmtCoef2 = nElmtCoef * nElmtCoef;
551 int nElmtPnt = (*expvect)[ne]->GetTotPoints();
555 int nQuotPlus = nQuot;
562 if (nElmtPnt > PntJacCons[0].size() || nElmtCoef > nElmtCoefOld)
564 nElmtCoefOld = nElmtCoef;
565 for (
int ndir = 0; ndir < 3; ndir++)
577 PntJacConsStd[ndir] =
579 for (
int i = 0; i < nElmtPnt; i++)
585 for (
int ndir1 = 0; ndir1 <
m_spacedim; ndir1++)
587 PntJacDerv[ndir][ndir1] =
589 PntJacDervStd[ndir][ndir1] =
593 for (
int i = 0; i < nElmtPnt; i++)
595 PntJacDerv[ndir][ndir1][i] =
597 PntJacDervStd[ndir][ndir1][i] =
605 for (
int j = 0; j < nvariable; j++)
607 locVars[j] = inarray[j] + noffset;
612 for (
int nFluxDir = 0; nFluxDir < nSpaceDim; nFluxDir++)
614 normals = normal3D[nFluxDir];
616 wspMat, PntJacCons[nFluxDir]);
622 for (
int j = 0; j < nSpaceDim; j++)
624 for (
int k = 0; k < nvariable; k++)
626 locDerv[j][k] = qfield[j][k] + noffset;
629 locmu = mu + noffset;
630 locDmuDT = DmuDT + noffset;
631 for (
int nFluxDir = 0; nFluxDir < nSpaceDim; nFluxDir++)
633 normals = normal3D[nFluxDir];
635 locDerv, locmu, locDmuDT, normals,
636 wspMatDrv, PntJacCons[nFluxDir]);
642 locmu = mu + noffset;
643 for (
int nFluxDir = 0; nFluxDir < nSpaceDim; nFluxDir++)
646 for (
int nDervDir = 0; nDervDir < nSpaceDim; nDervDir++)
649 nvariable, nElmtPnt, nDervDir, locVars, locmu,
650 normalPnt, wspMatDrv, PntJacDerv[nFluxDir][nDervDir]);
656 for (
int n = 0; n < nvariable; n++)
658 for (
int m = 0; m < nvariable; m++)
660 int nVarOffset = m + n * nvariable;
661 GmatxData = gmtxarray[m][n]->GetBlock(ne, ne)->GetPtr();
663 for (
int ndStd0 = 0; ndStd0 <
m_spacedim; ndStd0++)
669 for (
int i = 0; i < nElmtPnt; i++)
671 tmppnts[i] = PntJacCons[ndir][i][nVarOffset];
673 (*expvect)[ne]->AlignVectorToCollapsedDir(ndir, tmppnts,
677 Vmath::Vadd(nElmtPnt, ConsCurv[nd], 1, ConsStdd[nd], 1,
684 (*expvect)[ne]->MultiplyByQuadratureMetric(
685 ConsStdd[ndir], ConsStdd[ndir]);
686 for (
int i = 0; i < nElmtPnt; i++)
688 PntJacConsStd[ndir][i][nVarOffset] =
689 DataType(ConsStdd[ndir][i]);
697 for (
int m = 0; m < nvariable; m++)
699 for (
int n = 0; n < nvariable; n++)
701 int nVarOffset = m + n * nvariable;
702 for (
int ndStd0 = 0; ndStd0 <
m_spacedim; ndStd0++)
704 for (
int ndStd1 = 0; ndStd1 <
m_spacedim; ndStd1++)
706 Vmath::Zero(nElmtPnt, DervStdd[ndStd0][ndStd1], 1);
713 for (
int i = 0; i < nElmtPnt; i++)
716 PntJacDerv[nd0][nd1][i][nVarOffset];
719 (*expvect)[ne]->AlignVectorToCollapsedDir(
720 nd0, tmppnts, ConsCurv);
723 (*expvect)[ne]->AlignVectorToCollapsedDir(
724 nd1, ConsCurv[nd], DervCurv[nd]);
727 for (
int ndStd0 = 0; ndStd0 <
m_spacedim; ndStd0++)
733 DervCurv[ndStd0][ndStd1], 1,
734 DervStdd[ndStd0][ndStd1], 1,
735 DervStdd[ndStd0][ndStd1], 1);
744 (*expvect)[ne]->MultiplyByQuadratureMetric(
747 for (
int i = 0; i < nElmtPnt; i++)
749 PntJacDervStd[nd0][nd1][i][nVarOffset] =
750 -DataType(DervStdd[nd0][nd1][i]);
762 for (
int i = 0; i < nElmtPnt; i++)
765 &StdMatDataDBB[ne][ndir][i][0], 1,
766 &PntJacConsStd[ndir][i][0], 1, &MatData[0],
777 for (
int i = 0; i < nElmtPnt; i++)
780 &StdMatDataDBDB[ne][nd0][nd1][i][0], 1,
781 &PntJacDervStd[nd0][nd1][i][0], 1,
782 &MatData[0], nElmtCoef2Paded);
790 for (
int n = 0; n < nvariable; n++)
792 for (
int m = 0; m < nvariable; m++)
794 int nVarOffset = m + n * nvariable;
795 GmatxData = gmtxarray[m][n]->GetBlock(ne, ne)->GetPtr();
797 tmpA = MatData + nVarOffset * nElmtCoef2Paded, 1,
804 template <
typename DataType>
808 std::shared_ptr<LocalRegions::ExpansionVector> expvect =
810 int nTotElmt = (*expvect).size();
815 vector<DNekMatSharedPtr> VectStdDerivBase0;
816 vector<TensorOfArray3D<DataType>> VectStdDerivBase_Base;
817 vector<TensorOfArray4D<DataType>> VectStdDervBase_DervBase;
821 for (
int ne = 0; ne < nTotElmt; ne++)
824 stdExp = (*expvect)[ne]->GetStdExp();
826 stdExp->DetShapeType(), *stdExp);
827 MatStdDerivBase0 = stdExp->GetStdMatrix(matkey);
829 int nTotStdExp = VectStdDerivBase0.size();
830 int nFoundStdExp = -1;
831 for (
int i = 0; i < nTotStdExp; i++)
833 if ((*VectStdDerivBase0[i]) == (*MatStdDerivBase0))
838 if (nFoundStdExp >= 0)
840 StdMatDataDBB[ne] = VectStdDerivBase_Base[nFoundStdExp];
841 StdMatDataDBDB[ne] = VectStdDervBase_DervBase[nFoundStdExp];
845 int nElmtCoef = (*expvect)[ne]->GetNcoeffs();
846 int nElmtCoef2 = nElmtCoef * nElmtCoef;
847 int nElmtPnt = (*expvect)[ne]->GetTotPoints();
851 int nQuotPlus = nQuot;
858 ArrayStdMat[0] = MatStdDerivBase0;
863 ArrayStdMat[1] = stdExp->GetStdMatrix(matkey);
868 stdExp->DetShapeType(),
870 ArrayStdMat[2] = stdExp->GetStdMatrix(matkey);
875 ArrayStdMatData[nd0] = ArrayStdMat[nd0]->GetPtr();
879 stdExp->DetShapeType(), *stdExp);
889 for (
int i = 0; i < nElmtPnt; i++)
892 for (
int nc1 = 0; nc1 < nElmtCoef; nc1++)
894 int noffset = nc1 * nElmtCoef;
895 for (
int nc0 = 0; nc0 < nElmtCoef; nc0++)
897 tmpStdDBB[nd0][i][nc0 + noffset] = DataType(
898 ArrayStdMatData[nd0][i * nElmtCoef + nc0] *
899 BwdMatData[i * nElmtCoef + nc1]);
907 tmpStdDBDB[nd0][nd1] =
909 for (
int i = 0; i < nElmtPnt; i++)
911 tmpStdDBDB[nd0][nd1][i] =
913 for (
int nc1 = 0; nc1 < nElmtCoef; nc1++)
915 int noffset = nc1 * nElmtCoef;
916 for (
int nc0 = 0; nc0 < nElmtCoef; nc0++)
918 tmpStdDBDB[nd0][nd1][i][nc0 + noffset] =
921 [i * nElmtCoef + nc0] *
923 [i * nElmtCoef + nc1]);
929 VectStdDerivBase0.push_back(MatStdDerivBase0);
930 VectStdDerivBase_Base.push_back(tmpStdDBB);
931 VectStdDervBase_DervBase.push_back(tmpStdDBDB);
933 StdMatDataDBB[ne] = tmpStdDBB;
934 StdMatDataDBDB[ne] = tmpStdDBDB;
939 template <
typename DataType,
typename TypeNekBlkMatSharedPtr>
949 int nvariables = inarray.size();
953 GetTraceJac(inarray, qfield, TraceJac, TraceJacDeriv, TraceJacDerivSign,
963 gmtxarray, tmpJac, tmpSign);
968 template <
typename DataType,
typename TypeNekBlkMatSharedPtr>
971 TypeNekBlkMatSharedPtr &gmtVar,
const DataType &tmpDataType)
973 boost::ignore_unused(tmpDataType);
975 int n1d = gmtxarray.size();
976 int n2d = gmtxarray[0].size();
977 int nConvectiveFields = n1d;
979 ASSERTL0(n1d == n2d,
"ElmtVarInvMtrx requires n1d==n2d");
984 gmtxarray[0][0]->GetBlockSizes(rowSizes, colSizes);
985 int nTotElmt = rowSizes.size();
986 int nElmtCoef = rowSizes[0] - 1;
991 TypeNekBlkMatSharedPtr tmpGmtx;
996 for (
int nelmt = 0; nelmt < nTotElmt; nelmt++)
998 int nrows = gmtxarray[0][0]->GetBlock(nelmt, nelmt)->GetRows();
999 int ncols = gmtxarray[0][0]->GetBlock(nelmt, nelmt)->GetColumns();
1000 ASSERTL0(nrows == ncols,
"ElmtVarInvMtrx requires nrows==ncols");
1004 if (nElmtCoef0 != nElmtCoef)
1006 nElmtCoef0 = nElmtCoef;
1007 int nElmtCoefVar = nElmtCoef0 * nConvectiveFields;
1008 blocksize = nElmtCoefVar * nElmtCoefVar;
1009 tmprow[0] = nElmtCoefVar;
1011 GMatData = tmpGmtx->GetBlock(0, 0)->GetPtr();
1014 for (
int n = 0; n < nConvectiveFields; n++)
1016 for (
int m = 0; m < nConvectiveFields; m++)
1018 ElmtMatData = gmtxarray[m][n]->GetBlock(nelmt, nelmt)->GetPtr();
1020 for (
int ncl = 0; ncl < nElmtCoef; ncl++)
1023 (n * nElmtCoef + ncl) * nConvectiveFields * nElmtCoef +
1025 int Eoffset = ncl * nElmtCoef;
1027 Vmath::Vcopy(nElmtCoef, tmpArray1 = ElmtMatData + Eoffset,
1028 1, tmpArray2 = GMatData + Goffset, 1);
1033 tmpGmtx->GetBlock(0, 0)->Invert();
1035 for (
int m = 0; m < nConvectiveFields; m++)
1037 for (
int n = 0; n < nConvectiveFields; n++)
1039 ElmtMatData = gmtxarray[m][n]->GetBlock(nelmt, nelmt)->GetPtr();
1041 for (
int ncl = 0; ncl < nElmtCoef; ncl++)
1044 (n * nElmtCoef + ncl) * nConvectiveFields * nElmtCoef +
1046 int Eoffset = ncl * nElmtCoef;
1048 Vmath::Vcopy(nElmtCoef, tmpArray1 = GMatData + Goffset, 1,
1049 tmpArray2 = ElmtMatData + Eoffset, 1);
1053 ElmtMatData = gmtVar->GetBlock(nelmt, nelmt)->GetPtr();
1059 template <
typename DataType,
typename TypeNekBlkMatSharedPtr>
1068 boost::ignore_unused(TraceJacDeriv, TraceJacDerivSign);
1070 int nvariables = inarray.size();
1077 TypeNekBlkMatSharedPtr FJac, BJac;
1080 if (TraceJac.size() > 0)
1100 for (
int i = 0; i < nvariables; ++i)
1104 m_fields[i]->GetFwdBwdTracePhys(inarray[i], Fwd[i], Bwd[i]);
1112 TraceIPSymJacArray);
1118 template <
typename DataType,
typename TypeNekBlkMatSharedPtr>
1120 const int nConvectiveFields,
1127 TypeNekBlkMatSharedPtr &FJac, TypeNekBlkMatSharedPtr &BJac,
1130 boost::ignore_unused(TraceIPSymJacArray);
1133 int nvariables = nConvectiveFields;
1142 for (
int i = 0; i < nvariables; i++)
1154 for (
int i = 0; i < nvariables; ++i)
1160 fields[0]->GetTraceMap();
1165 for (
int nd = 0; nd < nDim; ++nd)
1169 for (
int i = 0; i < nConvectiveFields; ++i)
1174 fields[i]->GetFwdBwdTracePhys(qfield[nd][i], qFwd[nd][i],
1175 qBwd[nd][i],
true,
true,
false);
1176 TraceMap->GetAssemblyCommDG()->PerformExchange(qFwd[nd][i],
1183 PenaltyFactor2, fields, AdvVel, inarray, time,
1184 qfield, Fwd, Bwd, qFwd, qBwd, MuVarTrace,
1185 nonZeroIndex, numflux);
1187 int nFields = nvariables;
1192 for (
int i = 0; i < nFields; i++)
1201 for (
int i = 0; i < nFields; i++)
1211 for (
int i = 0; i < nFields; i++)
1215 Vmath::Sadd(nTracePts, epsvar, Fwd[i], 1, plusFwd[i], 1);
1219 for (
int i = 0; i < nFields; i++)
1226 x->Apply(FwdBnd, tmpinarry, time);
1230 for (
int j = 0; j < nFields; j++)
1232 m_fields[j]->FillBwdWithBoundCond(plusFwd[j], plusBwd[j]);
1236 PenaltyFactor2, fields, AdvVel, inarray, time,
1237 qfield, plusFwd, plusBwd, qFwd, qBwd, MuVarTrace,
1238 nonZeroIndex, plusflux);
1240 for (
int n = 0; n < nFields; n++)
1242 Vmath::Vsub(nTracePts, plusflux[n], 1, numflux[n], 1, Jacvect[n],
1244 Vmath::Smul(nTracePts, oepsvar, Jacvect[n], 1, Jacvect[n], 1);
1246 for (
int j = 0; j < nTracePts; j++)
1248 tmpMatData = FJac->GetBlock(j, j)->GetPtr();
1249 for (
int n = 0; n < nFields; n++)
1251 tmpMatData[n + i * nFields] = DataType(Jacvect[n][j]);
1261 for (
int i = 0; i < nFields; i++)
1268 x->Apply(FwdBnd, tmpinarry, time);
1272 for (
int i = 0; i < nFields; i++)
1277 for (
int i = 0; i < nFields; i++)
1282 Vmath::Sadd(nTracePts, epsvar, Bwd[i], 1, plusBwd[i], 1);
1284 for (
int j = 0; j < nFields; j++)
1286 m_fields[j]->FillBwdWithBoundCond(Fwd[j], plusBwd[j]);
1290 PenaltyFactor2, fields, AdvVel, inarray, time,
1291 qfield, Fwd, plusBwd, qFwd, qBwd, MuVarTrace,
1292 nonZeroIndex, plusflux);
1294 for (
int n = 0; n < nFields; n++)
1296 Vmath::Vsub(nTracePts, plusflux[n], 1, numflux[n], 1, Jacvect[n],
1298 Vmath::Smul(nTracePts, oepsvar, Jacvect[n], 1, Jacvect[n], 1);
1300 for (
int j = 0; j < nTracePts; j++)
1302 tmpMatData = BJac->GetBlock(j, j)->GetPtr();
1303 for (
int n = 0; n < nFields; n++)
1305 tmpMatData[n + i * nFields] = DataType(Jacvect[n][j]);
1314 const int nConvectiveFields,
const int nDim,
const int nPts,
1315 const int nTracePts,
const NekDouble PenaltyFactor2,
1327 boost::ignore_unused(nDim, nPts, PenaltyFactor2, time, qFwd, qBwd,
1338 for (
int i = 0; i < nConvectiveFields; i++)
1347 for (
int i = 0; i < nConvectiveFields; i++)
1352 m_diffusion->DiffuseTraceFlux(fields, inarray, qfield,
1354 vFwd, vBwd, nonZeroIndex);
1355 for (
int i = 0; i < nConvectiveFields; i++)
1357 Vmath::Vsub(nTracePts, traceflux[i], 1, visflux[i], 1, traceflux[i],
1392 TraceJacDerivSign, TraceIPSymJacArray);
1404 TraceJacDerivArray);
1408 const int nConvectiveFields,
const int nElmtPnt,
1415 boost::ignore_unused(nConvectiveFields, nElmtPnt, locVars, locDerv, locmu,
1416 locDmuDT, normals, wspMat, PntJacArray);
1420 template <
typename DataType,
typename TypeNekBlkMatSharedPtr>
1423 const NekDouble dtlamda,
const DataType tmpDataType)
1425 boost::ignore_unused(tmpDataType);
1428 std::shared_ptr<LocalRegions::ExpansionVector> pexp = explist->GetExp();
1429 int nTotElmt = (*pexp).size();
1430 int nConvectiveFields =
m_fields.size();
1435 Vmath::Fill(nTotElmt, Negdtlamda, pseudotimefactor, 1);
1438 for (
int m = 0; m < nConvectiveFields; m++)
1440 for (
int n = 0; n < nConvectiveFields; n++)
1442 for (
int nelmt = 0; nelmt < nTotElmt; nelmt++)
1444 GMatData = gmtxarray[m][n]->GetBlock(nelmt, nelmt)->GetPtr();
1445 DataType factor = DataType(pseudotimefactor[nelmt]);
1458 for (
int nelmt = 0; nelmt < nTotElmt; nelmt++)
1463 explist->GetExp(nelmt)->DetShapeType();
1465 if (tmp.size() != nelmtcoef || (ElmtTypeNow != ElmtTypePrevious))
1468 stdExp = explist->GetExp(nelmt)->GetStdExp();
1470 stdExp->DetShapeType(), *stdExp);
1473 BwdMatData = BwdMat->GetPtr();
1475 if (nelmtcoef != tmp.size())
1479 nelmtcoef, nelmtcoef, 0.0);
1480 MassMatData = MassMat->GetPtr();
1481 MassMatDataDataType =
1485 ElmtTypePrevious = ElmtTypeNow;
1488 for (
int np = 0; np < nelmtcoef; np++)
1490 explist->GetExp(nelmt)->IProductWRTBase(BwdMatData + np * nelmtpnts,
1492 Vmath::Vcopy(nelmtcoef, tmp, 1, tmp2 = MassMatData + np * nelmtcoef,
1495 for (
int i = 0; i < MassMatData.size(); i++)
1497 MassMatDataDataType[i] = DataType(MassMatData[i]);
1500 for (
int m = 0; m < nConvectiveFields; m++)
1502 GMatData = gmtxarray[m][m]->GetBlock(nelmt, nelmt)->GetPtr();
1503 Vmath::Vadd(MassMatData.size(), MassMatDataDataType, 1, GMatData, 1,
1510 template <
typename DataType,
typename TypeNekBlkMatSharedPtr>
1517 boost::ignore_unused(TraceJacArray, TraceJacDeriv, TraceJacDerivArray);
1519 int nFwdBwd, nDiagBlks, nvar0Jac, nvar1Jac;
1523 nFwdBwd = TraceJac.size();
1524 TraceJac[0]->GetBlockSizes(rowSizes, colSizes);
1525 nDiagBlks = rowSizes.size();
1526 nvar0Jac = rowSizes[1] - rowSizes[0];
1527 nvar1Jac = colSizes[1] - colSizes[0];
1529 if (0 == TraceJacArray.size())
1532 for (
int nlr = 0; nlr < nFwdBwd; nlr++)
1535 for (
int m = 0; m < nvar0Jac; m++)
1537 TraceJacArray[nlr][m] =
1539 for (
int n = 0; n < nvar1Jac; n++)
1547 for (
int nlr = 0; nlr < nFwdBwd; nlr++)
1549 const TypeNekBlkMatSharedPtr tmpMat = TraceJac[nlr];
1564 boost::ignore_unused(explist, normals, nDervDir, inarray, ElmtJacArray,
1570 const int nConvectiveFields,
const int nElmtPnt,
const int nDervDir,
1576 boost::ignore_unused(nConvectiveFields, nElmtPnt, nDervDir, locVars, locmu,
1577 locnormal, wspMat, PntJacArray);
1588 boost::ignore_unused(explist, normals, nDervDir, inarray, ElmtJac);
1592 const int nConvectiveFields,
const int nElmtPnt,
1599 int matsize = nConvectiveFields * nConvectiveFields;
1603 for (
int npnt = 0; npnt < nElmtPnt; npnt++)
1605 for (
int j = 0; j < nConvectiveFields; j++)
1607 pointVar[j] = locVars[j][npnt];
1612 Vmath::Vcopy(matsize, wspMatData, 1, PntJacArray[npnt], 1);
1621 int nvariables = conservVar.size();
1622 const int nvariables3D = 5;
1626 efix_StegerWarming = 0.0;
1628 if (nvariables > expDim + 2)
1635 fluxJacData = fluxJac->GetPtr();
1637 if (nConvectiveFields == nvariables3D)
1649 PointFJac3DData = PointFJac3D->GetPtr();
1655 index[nvariables - 1] = 4;
1656 for (
int i = 0; i < nvariables - 1; i++)
1663 for (
int j = 0; j < nvariables; j++)
1666 PointFwd[nj] = conservVar[j];
1670 efix_StegerWarming, fsw);
1672 for (
int j = 0; j < nvariables; j++)
1675 for (
int k = 0; k < nvariables; k++)
1678 fluxJacData[j + k * nConvectiveFields] =
1679 PointFJac3DData[nj + nk * nvariables3D];
1702 boost::ignore_unused(flag);
1709 unsigned int nTotalGlobal = inarray.size();
1715 unsigned int ntotal = inarray.size();
1719 Vmath::Svtvp(ntotal, eps, inarray, 1, refsol, 1, solplus, 1);
1723 Vmath::Vsub(ntotal, resplus, 1, refres, 1, out, 1);
1729 template <
typename DataType,
typename TypeNekBlkMatSharedPtr>
1735 BlkMat->GetBlockSizes(rowSizes, colSizes);
1736 int nDiagBlks = rowSizes.size();
1737 int nvar0 = rowSizes[1] - rowSizes[0];
1738 int nvar1 = colSizes[1] - colSizes[0];
1742 for (
int i = 0; i < nDiagBlks; i++)
1744 ElmtMatData = BlkMat->GetBlock(i, i)->GetPtr();
1745 for (
int n = 0; n < nvar1; n++)
1747 int noffset = n * nvar0;
1748 for (
int m = 0; m < nvar0; m++)
1750 MatArray[m][n][i] = ElmtMatData[m + noffset];
1756 template <
typename DataType,
typename TypeNekBlkMatSharedPtr>
1759 const DataType valu)
1762 int n1d = gmtxarray.size();
1764 for (
int n1 = 0; n1 < n1d; ++n1)
1770 template <
typename DataType,
typename TypeNekBlkMatSharedPtr>
1774 int n1d = gmtxarray.size();
1781 for (
int n1 = 0; n1 < n1d; ++n1)
1783 gmtxarray[n1]->GetBlockSizes(rowSizes, colSizes);
1784 int nelmts = rowSizes.size();
1786 for (
int i = 0; i < nelmts; ++i)
1788 loc_mat_arr = gmtxarray[n1]->GetBlock(i, i)->GetPtr();
1790 int nrows = gmtxarray[n1]->GetBlock(i, i)->GetRows();
1791 int ncols = gmtxarray[n1]->GetBlock(i, i)->GetColumns();
1809 const int nvariables3D = 5;
1812 NekDouble a, a2, h, h0, v2, vn, eps, eps2;
1815 NekDouble l1, l4, l5, al1, al4, al5, x1, x2, x3, y1;
1816 NekDouble c1, d1, c2, d2, c3, d3, c4, d4, c5, d5;
1833 NekDouble eL = (EL - 0.5 * (rhouL * vx + rhovL * vy + rhowL * vz)) / rhoL;
1835 ps =
m_varConv->Geteos()->GetPressure(rhoL, eL);
1839 v2 = vx * vx + vy * vy + vz * vz;
1840 a2 = gama * ps / ro;
1849 vn = nx * vx + ny * vy + nz * vz;
1850 sn = std::max(
sqrt(nx * nx + ny * ny + nz * nz), sml_ssf);
1864 al1 =
sqrt(l1 * l1 + eps2);
1865 al4 =
sqrt(l4 * l4 + eps2);
1866 al5 =
sqrt(l5 * l5 + eps2);
1868 l1 = 0.5 * (fExactorSplt * l1 + fsw * al1);
1869 l4 = 0.5 * (fExactorSplt * l4 + fsw * al4);
1870 l5 = 0.5 * (fExactorSplt * l5 + fsw * al5);
1872 x1 = 0.5 * (l4 + l5);
1873 x2 = 0.5 * (l4 - l5);
1880 int nVar1 = nvariables3D;
1881 int nVar2 = 2 * nvariables3D;
1882 int nVar3 = 3 * nvariables3D;
1883 int nVar4 = 4 * nvariables3D;
1884 FJacData[nVar0] = c1 * y1 - d1 * vna + l1;
1885 FJacData[nVar1] = -c1 * vx + d1 * nxa;
1886 FJacData[nVar2] = -c1 * vy + d1 * nya;
1887 FJacData[nVar3] = -c1 * vz + d1 * nza;
1888 FJacData[nVar4] = c1;
1889 c2 = c1 * vx + d1 * nxa * ae;
1890 d2 = x3 * nxa + d1 * vx;
1891 FJacData[1 + nVar0] = c2 * y1 - d2 * vna;
1892 FJacData[1 + nVar1] = -c2 * vx + d2 * nxa + l1;
1893 FJacData[1 + nVar2] = -c2 * vy + d2 * nya;
1894 FJacData[1 + nVar3] = -c2 * vz + d2 * nza;
1895 FJacData[1 + nVar4] = c2;
1896 c3 = c1 * vy + d1 * nya * ae;
1897 d3 = x3 * nya + d1 * vy;
1898 FJacData[2 + nVar0] = c3 * y1 - d3 * vna;
1899 FJacData[2 + nVar1] = -c3 * vx + d3 * nxa;
1900 FJacData[2 + nVar2] = -c3 * vy + d3 * nya + l1;
1901 FJacData[2 + nVar3] = -c3 * vz + d3 * nza;
1902 FJacData[2 + nVar4] = c3;
1903 c4 = c1 * vz + d1 * nza * ae;
1904 d4 = x3 * nza + d1 * vz;
1905 FJacData[3 + nVar0] = c4 * y1 - d4 * vna;
1906 FJacData[3 + nVar1] = -c4 * vx + d4 * nxa;
1907 FJacData[3 + nVar2] = -c4 * vy + d4 * nya;
1908 FJacData[3 + nVar3] = -c4 * vz + d4 * nza + l1;
1909 FJacData[3 + nVar4] = c4;
1910 c5 = c1 * h0 + d1 * vna * ae;
1911 d5 = x3 * vna + d1 * h0;
1912 FJacData[4 + nVar0] = c5 * y1 - d5 * vna;
1913 FJacData[4 + nVar1] = -c5 * vx + d5 * nxa;
1914 FJacData[4 + nVar2] = -c5 * vy + d5 * nya;
1915 FJacData[4 + nVar3] = -c5 * vz + d5 * nza;
1916 FJacData[4 + nVar4] = c5 + l1;
#define ASSERTL0(condition, msg)
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
size_type size() const
Returns the array's size.
void DoDiffusionCoeff(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const Array< OneD, const Array< OneD, NekDouble >> &pFwd, const Array< OneD, const Array< OneD, NekDouble >> &pBwd)
Add the diffusions terms to the right-hand side Similar to DoDiffusion() but with outarray in coeffic...
void CalcVolJacStdMat(TensorOfArray4D< DataType > &StdMatDataDBB, TensorOfArray5D< DataType > &StdMatDataDBDB)
void ElmtVarInvMtrx(Array< OneD, Array< OneD, TypeNekBlkMatSharedPtr >> &gmtxarray, TypeNekBlkMatSharedPtr &gmtVar, const DataType &tmpDatatype)
void DoAdvectionCoeff(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble time, const Array< OneD, const Array< OneD, NekDouble >> &pFwd, const Array< OneD, const Array< OneD, NekDouble >> &pBwd)
Compute the advection terms for the right-hand side.
virtual void v_DoDiffusionCoeff(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const Array< OneD, const Array< OneD, NekDouble >> &pFwd, const Array< OneD, const Array< OneD, NekDouble >> &pBwd)
LibUtilities::NekNonlinSysSharedPtr m_nonlinsol
void AddMatNSBlkDiagVol(const Array< OneD, const Array< OneD, NekDouble >> &inarray, const Array< OneD, const TensorOfArray2D< NekDouble >> &qfield, Array< OneD, Array< OneD, TypeNekBlkMatSharedPtr >> &gmtxarray, TensorOfArray4D< DataType > &StdMatDataDBB, TensorOfArray5D< DataType > &StdMatDataDBDB)
void PointFluxJacobianPoint(const Array< OneD, NekDouble > &Fwd, const Array< OneD, NekDouble > &normals, DNekMatSharedPtr &FJac, const NekDouble efix, const NekDouble fsw)
NekDouble m_jacobiFreeEps
void PreconCoeff(const Array< OneD, NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const bool &flag)
void MultiplyElmtInvMassPlusSource(Array< OneD, Array< OneD, TypeNekBlkMatSharedPtr >> &gmtxarray, const NekDouble dtlamda, const DataType tmpDataType)
void NonlinSysEvaluatorCoeff1D(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out, const bool &flag)
void CalcRefValues(const Array< OneD, const NekDouble > &inarray)
virtual bool UpdateTimeStepCheck()
virtual void v_GetFluxDerivJacDirctn(const MultiRegions::ExpListSharedPtr &explist, const Array< OneD, const Array< OneD, NekDouble >> &normals, const int nDervDir, const Array< OneD, const Array< OneD, NekDouble >> &inarray, TensorOfArray5D< NekDouble > &ElmtJacArray, const int nFluxDir)
void GetFluxVectorJacDirElmt(const int nConvectiveFields, const int nElmtPnt, const Array< OneD, const Array< OneD, NekDouble >> &locVars, const Array< OneD, NekDouble > &normals, DNekMatSharedPtr &wspMat, Array< OneD, Array< OneD, NekDouble >> &PntJacArray)
void Fill2DArrayOfBlkDiagonalMat(Array< OneD, Array< OneD, TypeNekBlkMatSharedPtr >> &gmtxarray, const DataType valu)
void DoImplicitSolveCoeff(const Array< OneD, const Array< OneD, NekDouble >> &inpnts, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out, const NekDouble time, const NekDouble lambda)
void CalcTraceNumericalFlux(const int nConvectiveFields, const int nDim, const int nPts, const int nTracePts, const NekDouble PenaltyFactor2, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const Array< OneD, NekDouble >> &AdvVel, const Array< OneD, const Array< OneD, NekDouble >> &inarray, const NekDouble time, TensorOfArray3D< NekDouble > &qfield, const Array< OneD, const Array< OneD, NekDouble >> &vFwd, const Array< OneD, const Array< OneD, NekDouble >> &vBwd, const Array< OneD, const TensorOfArray2D< NekDouble >> &qFwd, const Array< OneD, const TensorOfArray2D< NekDouble >> &qBwd, const Array< OneD, NekDouble > &MuVarTrace, Array< OneD, int > &nonZeroIndex, Array< OneD, Array< OneD, NekDouble >> &traceflux)
Array< OneD, Array< OneD, NekDouble > > m_solutionPhys
void AllocateNekBlkMatDig(SNekBlkMatSharedPtr &mat, const Array< OneD, unsigned int > nrow, const Array< OneD, unsigned int > ncol)
TensorOfArray4D< NekSingle > m_stdSMatDataDBB
void NonlinSysEvaluatorCoeff(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out, const bool &flag=true, const Array< OneD, const NekDouble > &source=NullNekDouble1DArray)
void TranSamesizeBlkDiagMatIntoArray(const TypeNekBlkMatSharedPtr &BlkMat, TensorOfArray3D< DataType > &MatArray)
void NumCalcRiemFluxJac(const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const Array< OneD, NekDouble >> &AdvVel, const Array< OneD, const Array< OneD, NekDouble >> &inarray, TensorOfArray3D< NekDouble > &qfield, const NekDouble &time, const Array< OneD, const Array< OneD, NekDouble >> &Fwd, const Array< OneD, const Array< OneD, NekDouble >> &Bwd, TypeNekBlkMatSharedPtr &FJac, TypeNekBlkMatSharedPtr &BJac, TensorOfArray5D< DataType > &TraceIPSymJacArray)
PreconCfsOpSharedPtr m_preconCfs
NekDouble m_newtonAbsoluteIteTol
virtual ~CFSImplicit()
Destructor for CFSImplicit class.
void DoOdeRhsCoeff(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble time)
Compute the right-hand side.
void GetTraceJac(const Array< OneD, const Array< OneD, NekDouble >> &inarray, TensorOfArray3D< NekDouble > &qfield, Array< OneD, TypeNekBlkMatSharedPtr > &TraceJac, Array< OneD, TypeNekBlkMatSharedPtr > &TraceJacDeriv, Array< OneD, Array< OneD, DataType >> &TraceJacDerivSign, TensorOfArray5D< DataType > &TraceIPSymJacArray)
void GetFluxVectorJacPoint(const int nConvectiveFields, const Array< OneD, NekDouble > &conservVar, const Array< OneD, NekDouble > &normals, DNekMatSharedPtr &fluxJac)
TensorOfArray5D< NekSingle > m_stdSMatDataDBDB
void GetFluxDerivJacDirctnElmt(const int nConvectiveFields, const int nElmtPnt, const int nDervDir, const Array< OneD, const Array< OneD, NekDouble >> &locVars, const Array< OneD, NekDouble > &locmu, const Array< OneD, const Array< OneD, NekDouble >> &locnormal, DNekMatSharedPtr &wspMat, Array< OneD, Array< OneD, NekDouble >> &PntJacArray)
void MatrixMultiplyMatrixFreeCoeff(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out, const bool &flag=false)
virtual void v_MinusDiffusionFluxJacPoint(const int nConvectiveFields, const int nElmtPnt, const Array< OneD, const Array< OneD, NekDouble >> &locVars, const TensorOfArray3D< NekDouble > &locDerv, const Array< OneD, NekDouble > &locmu, const Array< OneD, NekDouble > &locDmuDT, const Array< OneD, NekDouble > &normals, DNekMatSharedPtr &wspMat, Array< OneD, Array< OneD, NekDouble >> &PntJacArray)
void CalcPreconMatBRJCoeff(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, SNekBlkMatSharedPtr >> &gmtxarray, SNekBlkMatSharedPtr &gmtVar, Array< OneD, SNekBlkMatSharedPtr > &TraceJac, Array< OneD, SNekBlkMatSharedPtr > &TraceJacDeriv, Array< OneD, Array< OneD, NekSingle >> &TraceJacDerivSign, TensorOfArray4D< NekSingle > &TraceJacArray, TensorOfArray4D< NekSingle > &TraceJacDerivArray, TensorOfArray5D< NekSingle > &TraceIPSymJacArray)
void TransTraceJacMatToArray(const Array< OneD, TypeNekBlkMatSharedPtr > &TraceJac, const Array< OneD, TypeNekBlkMatSharedPtr > &TraceJacDeriv, TensorOfArray4D< DataType > &TraceJacArray, TensorOfArray4D< DataType > &TraceJacDerivArray)
void AddMatNSBlkDiagBnd(const Array< OneD, const Array< OneD, NekDouble >> &inarray, TensorOfArray3D< NekDouble > &qfield, TensorOfArray2D< TypeNekBlkMatSharedPtr > &gmtxarray, Array< OneD, TypeNekBlkMatSharedPtr > &TraceJac, Array< OneD, TypeNekBlkMatSharedPtr > &TraceJacDeriv, Array< OneD, Array< OneD, DataType >> &TraceJacDerivSign, TensorOfArray5D< DataType > &TraceIPSymJacArray)
void CalcPhysDeriv(const Array< OneD, const Array< OneD, NekDouble >> &inarray, TensorOfArray3D< NekDouble > &qfield)
void MinusDiffusionFluxJacPoint(const int nConvectiveFields, const int nElmtPnt, const Array< OneD, const Array< OneD, NekDouble >> &locVars, const TensorOfArray3D< NekDouble > &locDerv, const Array< OneD, NekDouble > &locmu, const Array< OneD, NekDouble > &locDmuDT, const Array< OneD, NekDouble > &normals, DNekMatSharedPtr &wspMat, Array< OneD, Array< OneD, NekDouble >> &PntJacArray)
void DoImplicitSolve(const Array< OneD, const Array< OneD, NekDouble >> &inpnts, Array< OneD, Array< OneD, NekDouble >> &outpnt, const NekDouble time, const NekDouble lambda)
virtual void v_InitObject(bool DeclareFields=true)
Initialization object for CFSImplicit class.
void Fill1DArrayOfBlkDiagonalMat(Array< OneD, TypeNekBlkMatSharedPtr > &gmtxarray, const DataType valu)
void InitialiseNonlinSysSolver()
void CalcMuDmuDT(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, NekDouble > &mu, Array< OneD, NekDouble > &DmuDT)
virtual void v_GetFluxDerivJacDirctnElmt(const int nConvectiveFields, const int nElmtPnt, const int nDervDir, const Array< OneD, const Array< OneD, NekDouble >> &locVars, const Array< OneD, NekDouble > &locmu, const Array< OneD, const Array< OneD, NekDouble >> &locnormal, DNekMatSharedPtr &wspMat, Array< OneD, Array< OneD, NekDouble >> &PntJacArray)
virtual void v_InitObject(bool DeclareFields=true)
Initialization object for CompressibleFlowSystem class.
NekDouble m_bndEvaluateTime
SolverUtils::DiffusionSharedPtr m_diffusion
std::vector< CFSBndCondSharedPtr > m_bndConds
void GetElmtTimeStep(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, NekDouble > &tstep)
Calculate the maximum timestep on each element subject to CFL restrictions.
void DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble time)
Compute the projection and call the method for imposing the boundary conditions in case of discontinu...
VariableConverterSharedPtr m_varConv
std::vector< SolverUtils::ForcingSharedPtr > m_forcing
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
NekDouble m_NonlinIterTolRelativeL2
int m_NekNonlinSysMaxIterations
NekDouble m_LinSysRelativeTolInNonlin
int m_NekLinSysMaxIterations
void DefineNekSysResEval(FuncPointerT func, ObjectPointerT obj)
void DefineNekSysLhsEval(FuncPointerT func, ObjectPointerT obj)
void DefineNekSysPrecon(FuncPointerT func, ObjectPointerT obj)
void DefineImplicitSolve(FuncPointerT func, ObjectPointerT obj)
void AccumulateRegion(std::string, int iolevel=0)
Accumulate elapsed time for a region.
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
void DefineCalcPreconMatBRJCoeff(FuncPointerT func, ObjectPointerT obj)
SolverUtils::AdvectionSharedPtr m_advObject
Advection term.
int m_spacedim
Spatial dimension (>= expansion dim).
SpatialDomains::MeshGraphSharedPtr m_graph
Pointer to graph defining mesh.
LibUtilities::CommSharedPtr m_comm
Communicator.
NekDouble m_timestep
Time step size.
NekDouble m_time
Current time of simulation.
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
NekDouble m_fintime
Finish time of the simulation.
NekDouble m_lastCheckTime
NekDouble m_checktime
Time between checkpoints.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
enum HomogeneousType m_HomogeneousType
SOLVER_UTILS_EXPORT int GetNpoints()
SOLVER_UTILS_EXPORT int GetNcoeffs()
SOLVER_UTILS_EXPORT int GetPhys_Offset(int n)
SOLVER_UTILS_EXPORT int GetTotPoints()
Base class for unsteady solvers.
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
NekDouble m_TimeIntegLambda
coefff of spacial derivatives(rhs or m_F in GLM) in calculating the residual of the whole equation(us...
bool m_explicitAdvection
Indicates if explicit or implicit treatment of advection is used.
bool m_explicitDiffusion
Indicates if explicit or implicit treatment of diffusion is used.
bool m_flagImplicitSolver
Array< OneD, NekDouble > m_magnitdEstimat
estimate the magnitude of each conserved varibles
bool m_flagUpdatePreconMat
static void Ger(const int &m, const int &n, const double &alpha, const double *x, const int &incx, const double *y, const int &incy, double *a, const int &lda)
BLAS level 2: Matrix vector multiply A = alpha*x*y**T + A where A[m x n].
NekNonlinSysFactory & GetNekNonlinSysFactory()
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< AssemblyMapDG > AssemblyMapDGSharedPtr
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
std::shared_ptr< LocTraceToTraceMap > LocTraceToTraceMapSharedPtr
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
The above copyright notice and this permission notice shall be included.
std::shared_ptr< SNekBlkMat > SNekBlkMatSharedPtr
static Array< OneD, Array< OneD, NekDouble > > NullNekDoubleArrayOfArray
static Array< OneD, Array< OneD, Array< OneD, NekDouble > > > NullNekDoubleTensorOfArray3D
static Array< OneD, NekDouble > NullNekDouble1DArray
std::shared_ptr< DNekMat > DNekMatSharedPtr
PreconCfsOpFactory & GetPreconCfsOpFactory()
Declaration of the boundary condition factory singleton.
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
void Neg(int n, T *x, const int incx)
Negate x = -x.
T Dot(int n, const T *w, const T *x)
vvtvp (vector times vector times vector): z = w*x*y
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 Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
void Zero(int n, T *x, const int incx)
Zero vector.
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha - x.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
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.
scalarT< T > abs(scalarT< T > in)
scalarT< T > sqrt(scalarT< T > in)