50 "Block Relaxed Jacobi Preconditioner for CFS.");
61 size_t nvariables = pFields.size();
65 for (
size_t i = 0; i < nvariables; i++)
87 [[maybe_unused]]
const bool &flag)
89 ASSERTL0(inarray.size() == outarray.size(),
90 "In and Out not the same size in DoPreconCfs");
100 const NekDouble OmBRJParam = 1.0 - BRJParam;
102 size_t nvariables = pFields.size();
103 size_t npoints = pFields[0]->GetNcoeffs();
104 size_t ntotpnt = inarray.size();
106 ASSERTL0(nvariables * npoints == ntotpnt,
107 "nvariables*npoints!=ntotpnt in PreconCoeff");
116 for (
size_t m = 0; m < nvariables; m++)
118 size_t moffset = m * npoints;
119 rhs2d[m] = rhs + moffset;
120 out_2d[m] = outarray + moffset;
121 outTmp_2d[m] = outTmp + moffset;
122 pFields[m]->MultiplyByMassMatrix(inarray + moffset, rhs2d[m]);
125 size_t nphysic = pFields[0]->GetNpoints();
126 size_t nTracePts = pFields[0]->GetTrace()->GetNpoints();
131 for (
size_t j = 0; j < nvariables; j++)
138 for (
size_t i = 0; i < ntmpTrace; i++)
141 for (
size_t j = 0; j < nvariables; j++)
148 for (
size_t j = 0; j < nvariables; j++)
154 const size_t nwspTraceDataType = nvariables + 1;
156 for (
size_t m = 0; m < nwspTraceDataType; m++)
167 for (
size_t nrelax = 0; nrelax < nBRJIterTot - 1; nrelax++)
169 Vmath::Smul(ntotpnt, OmBRJParam, outarray, 1, outN, 1);
182 Vmath::Svtvp(ntotpnt, BRJParam, outTmp, 1, outN, 1, outarray, 1);
201 int nvariables = pFields.size();
202 int nelmts = pFields[0]->GetNumElmts();
204 for (
int i = 0; i < nelmts; i++)
206 matdim[i] = pFields[0]->GetExp(i)->GetNcoeffs() * nvariables;
218 std::cout <<
" ## CalcuPreconMat " << std::endl;
224 const auto vecwidth = vec_t::width;
226 alignas(vec_t::alignment) std::array<NekSingle, vec_t::width> tmp;
228 for (
int ne = 0; ne < nelmts; ne++)
230 const auto nElmtDof = matdim[ne];
231 const auto nblocks = nElmtDof / vecwidth;
234 PreconMatSingle->GetBlockPtr(ne, ne)->GetRawPtr();
236 for (
int i1 = 0; i1 < nblocks; ++i1)
238 for (
int j = 0; j < nElmtDof; ++j)
240 for (
int i = 0; i < vecwidth; ++i)
242 tmp[i] = mmat[j + (i1 * vecwidth + i) * nElmtDof];
249 const auto endwidth = nElmtDof - nblocks * vecwidth;
254 for (
int j = 0; j < nElmtDof; ++j)
256 for (
int i = 0; i < endwidth; ++i)
258 tmp[i] = mmat[j + (nblocks * vecwidth + i) * nElmtDof];
261 for (
int i = endwidth; i < vecwidth; ++i)
294 [[maybe_unused]]
const bool &flag)
306 unsigned int nvariables = pFields.size();
308 int nTotElmt = pFields[0]->GetNumElmts();
311 const auto vecwidth = vec_t::width;
314 std::vector<vec_t, tinysimd::allocator<vec_t>> Sinarray(
m_max_nblocks);
315 std::vector<vec_t, tinysimd::allocator<vec_t>> Soutarray(
m_max_nElmtDof);
317 alignas(vec_t::alignment) std::array<NekSingle, vec_t::width> tmp;
319 for (
int ne = 0, cnt = 0, icnt = 0, icnt1 = 0; ne < nTotElmt; ne++)
321 const auto nElmtCoef = pFields[0]->GetNcoeffs(ne);
322 const auto nElmtDof = nElmtCoef * nvariables;
323 const auto nblocks = (nElmtDof % vecwidth) ? nElmtDof / vecwidth + 1
324 : nElmtDof / vecwidth;
328 for (
int j = 0; j < nblocks; ++j, icnt += vecwidth)
330 for (
int i = 0; i < vecwidth; ++i)
335 Sinarray[j].load(tmp.data());
340 vec_t in = Sinarray[0];
341 for (
int i = 0; i < nElmtDof; ++i)
347 for (
int n = 1; n < nblocks; ++n)
350 for (
int i = 0; i < nElmtDof; ++i)
358 for (
int i = 0; i < nElmtDof; ++i)
361 Soutarray[i].store(tmp.data());
365 for (
int j = 1; j < vecwidth; ++j)
373 icnt1 += nblocks * vecwidth;
380template <
typename DataType>
383 const size_t nvariables,
const size_t nCoeffs,
390 size_t nTracePts = pFields[0]->GetTrace()->GetNpoints();
391 size_t npoints = pFields[0]->GetNpoints();
395 for (
size_t i = 0; i < nvariables; i++)
398 pFields[i]->BwdTrans(outarray[i], outpnts[i]);
408 size_t indexwspTrace = 0;
409 Fwd = wspTrace[indexwspTrace], indexwspTrace++;
410 Bwd = wspTrace[indexwspTrace], indexwspTrace++;
411 FwdFlux = wspTrace[indexwspTrace], indexwspTrace++;
412 BwdFlux = wspTrace[indexwspTrace], indexwspTrace++;
415 for (
size_t i = 0; i < nvariables; ++i)
418 pFields[i]->GetFwdBwdTracePhys(outpnts[i], Fwd[i], Bwd[i]);
423 size_t indexwspTraceDataType = 0;
425 for (
size_t m = 0; m < nvariables; ++m)
427 Fwdarray[m] = wspTraceDataType[indexwspTraceDataType],
428 indexwspTraceDataType++;
431 Fwdreslt = wspTraceDataType[indexwspTraceDataType], indexwspTraceDataType++;
433 for (
size_t m = 0; m < nvariables; ++m)
435 for (
size_t i = 0; i < nTracePts; ++i)
437 Fwdarray[m][i] = DataType(Fwd[m][i]);
440 for (
size_t m = 0; m < nvariables; ++m)
443 for (
size_t n = 0; n < nvariables; ++n)
445 for (
size_t p = 0;
p < nTracePts; ++
p)
447 Fwdreslt[
p] += TraceJacArray[0][m][n][
p] * Fwdarray[n][
p];
450 for (
size_t i = 0; i < nTracePts; ++i)
456 for (
size_t m = 0; m < nvariables; ++m)
458 for (
size_t i = 0; i < nTracePts; ++i)
460 Fwdarray[m][i] = DataType(Bwd[m][i]);
463 for (
size_t m = 0; m < nvariables; ++m)
466 for (
size_t n = 0; n < nvariables; ++n)
468 for (
size_t p = 0;
p < nTracePts; ++
p)
470 Fwdreslt[
p] += TraceJacArray[1][m][n][
p] * Fwdarray[n][
p];
473 for (
size_t i = 0; i < nTracePts; ++i)
479 for (
size_t i = 0; i < nvariables; ++i)
483 pFields[i]->AddTraceIntegralToOffDiag(FwdFlux[i], BwdFlux[i],
489 for (
size_t i = 0; i < nvariables; ++i)
491 for (
size_t p = 0;
p < nCoeffs; ++
p)
502template <
typename TypeNekBlkMatSharedPtr>
509 size_t nvars = pFields.size();
510 size_t nelmts = pFields[0]->GetNumElmts();
513 for (
size_t i = 0; i < nelmts; i++)
515 nelmtcoef = pFields[0]->GetExp(i)->GetNcoeffs();
516 nelmtmatdim[i] = nelmtcoef * nscale;
519 for (
size_t i = 0; i < nvars; i++)
521 for (
size_t j = 0; j < nvars; j++)
#define ASSERTL0(condition, msg)
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
void AccumulateRegion(std::string, int iolevel=0)
Accumulate elapsed time for a region.
void DoCalcPreconMatBRJCoeff(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)
Array< OneD, SNekBlkMatSharedPtr > m_TraceJacSingle
void v_BuildPreconCfs(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const Array< OneD, const Array< OneD, NekDouble > > &intmp, const NekDouble time, const NekDouble lambda) override
void v_DoPreconCfs(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const bool &flag) override
TensorOfArray4D< NekSingle > m_TraceJacArraySingle
void MinusOffDiag2Rhs(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const size_t nvariables, const size_t nCoeffs, const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, TensorOfArray3D< NekDouble > &wspTrace, Array< OneD, Array< OneD, DataType > > &wspTraceDataType, const TensorOfArray4D< DataType > &TraceJacArray)
Array< OneD, Array< OneD, SNekBlkMatSharedPtr > > m_PreconMatVarsSingle
static std::string className
Name of the class.
Array< OneD, SNekBlkMatSharedPtr > m_TraceJacDerivSingle
unsigned int m_max_nblocks
Array< OneD, Array< OneD, NekSingle > > m_TraceJacDerivSignSingle
static PreconCfsSharedPtr create(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const LibUtilities::SessionReaderSharedPtr &pSession, const LibUtilities::CommSharedPtr &vComm)
Creates an instance of this class.
void v_InitObject() override
TensorOfArray4D< NekSingle > m_TraceJacDerivArraySingle
PreconCfsBRJ(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const LibUtilities::SessionReaderSharedPtr &pSession, const LibUtilities::CommSharedPtr &vComm)
unsigned int m_max_nElmtDof
void AllocatePreconBlkDiagCoeff(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, Array< OneD, Array< OneD, TypeNekBlkMatSharedPtr > > &gmtxarray, const int &nscale=1)
TensorOfArray5D< NekSingle > m_TraceIPSymJacArraySingle
std::vector< int > m_inputIdx
void AllocateNekBlkMatDig(SNekBlkMatSharedPtr &mat, const Array< OneD, unsigned int > nrow, const Array< OneD, unsigned int > ncol)
void AllocateSIMDPreconBlkMatDiag(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
This function creates the matrix structure for the block diagonal operator. It organizes the way that...
void PreconBlkDiag(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const Array< OneD, NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
bool v_UpdatePreconMatCheck(const Array< OneD, const NekDouble > &res, const NekDouble dtLambda) override
std::vector< simd< NekSingle >, tinysimd::allocator< simd< NekSingle > > > m_sBlkDiagMat
void DoNullPrecon(const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const bool &flag)
NekPreconCfsOperators m_operator
NekDouble m_DtLambdaPreconMat
LibUtilities::CommSharedPtr m_Comm
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
PreconCfsFactory & GetPreconCfsFactory()
Declaration of the boundary condition factory singleton.
std::shared_ptr< SNekBlkMat > SNekBlkMatSharedPtr
tinysimd::simd< NekDouble > vec_t
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 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 Vcopy(int n, const T *x, const int incx, T *y, const int incy)
typename abi< ScalarType, width >::type simd