48 namespace MultiRegions
57 "Block Diagonal Preconditioning");
66 const boost::shared_ptr<GlobalLinSys> &plinsys,
70 m_preconType(pLocToGloMap->GetPreconType()),
71 m_locToGloMap(pLocToGloMap)
79 "Solver type not valid");
119 int nel, i, j, k, n, cnt, gId;
120 int meshVertId, meshEdgeId, meshFaceId;
122 const int nExp = expList->GetExpSize();
123 const int nDirBnd =
m_locToGloMap->GetNumGlobalDirBndCoeffs();
127 PeriodicMap periodicVerts, periodicEdges, periodicFaces;
128 expList->GetPeriodicEntities(
129 periodicVerts, periodicEdges, periodicFaces);
136 vector<map<int, vector<NekDouble> > > idMats(3);
140 vector<map<int, int> > gidMeshIds(3);
144 vector<map<int, int> > gidDofs(3);
159 for (cnt = n = 0; n < nExp; ++n)
161 nel = expList->GetOffset_Elmt_Id(n);
162 exp = expList->GetExp(nel);
166 m_linsys.lock()->GetStaticCondBlock(n)->GetBlock(0,0);
170 for (i = 0; i < exp->GetNverts(); ++i)
172 meshVertId = exp->GetGeom()->GetVid(i);
173 int locId = exp->GetVertexMap(i);
177 cnt + locId) - nDirBnd;
188 NekDouble vertVal = (*schurMat)(locId,locId);
192 gIt = idMats[0].find(gId);
194 if (gIt == idMats[0].end())
197 idMats[0][gId] = vector<NekDouble>(1, vertVal);
203 gIt->second[0] += vertVal;
210 pIt = periodicVerts.find(meshVertId);
211 if (pIt != periodicVerts.end())
213 for (j = 0; j < pIt->second.size(); ++j)
215 meshVertId = min(meshVertId, pIt->second[j].id);
221 gidMeshIds[0][gId] = meshVertId;
222 maxVertIds[0] = max(maxVertIds[0], meshVertId);
228 for (i = 0; i < exp->GetNedges(); ++i)
230 meshEdgeId = exp->GetGeom()->GetEid(i);
235 exp->GetGeom()->GetEorient(i);
239 pIt = periodicEdges.find(meshEdgeId);
240 if (pIt != periodicEdges.end())
242 pair<int, StdRegions::Orientation> idOrient =
244 meshEdgeId, edgeOrient, pIt->second);
245 meshEdgeId = idOrient.first;
246 edgeOrient = idOrient.second;
252 exp->GetEdgeInteriorMap(i, edgeOrient, bmap, sign);
253 bmap2 = exp->GetEdgeInverseBoundaryMap(i);
256 const int nEdgeCoeffs = bmap.num_elements();
257 vector<NekDouble> tmpStore(nEdgeCoeffs*nEdgeCoeffs);
261 for (j = 0; j < nEdgeCoeffs; ++j)
283 for (k = 0; k < nEdgeCoeffs; ++k)
285 tmpStore[k+j*nEdgeCoeffs] =
286 sign1*sign[k]*(*schurMat)(bmap2[j], bmap2[k]);
295 gidDofs[1][gId] = nEdgeCoeffs;
298 gIt = idMats[1].find(gId);
299 if (gIt == idMats[1].end())
301 idMats[1][gId] = tmpStore;
305 ASSERTL1(tmpStore.size() == gIt->second.size(),
306 "Number of modes mismatch");
307 Vmath::Vadd(nEdgeCoeffs*nEdgeCoeffs, &gIt->second[0], 1,
308 &tmpStore[0], 1, &gIt->second[0], 1);
311 gidMeshIds[1][gId] = meshEdgeId;
312 maxVertIds[2] = max(maxVertIds[2], meshEdgeId);
313 maxVertIds[3] = max(maxVertIds[3], nEdgeCoeffs);
318 for (i = 0; i < exp->GetNfaces(); ++i)
320 meshFaceId = exp->GetGeom()->GetFid(i);
325 exp->GetGeom()->GetForient(i);
329 pIt = periodicFaces.find(meshFaceId);
330 if (pIt != periodicFaces.end())
332 meshFaceId = min(meshFaceId, pIt->second[0].id);
334 faceOrient, pIt->second[0].orient);
337 exp->GetFaceInteriorMap(i, faceOrient, bmap, sign);
338 bmap2 = exp->GetFaceInverseBoundaryMap(i);
341 const int nFaceCoeffs = bmap.num_elements();
342 vector<NekDouble> tmpStore(nFaceCoeffs*nFaceCoeffs);
346 for (j = 0; j < nFaceCoeffs; ++j)
363 for (k = 0; k < nFaceCoeffs; ++k)
365 tmpStore[k+j*nFaceCoeffs] =
366 sign1*sign[k]*(*schurMat)(bmap2[j], bmap2[k]);
375 gidDofs[2][gId] = nFaceCoeffs;
378 gIt = idMats[2].find(gId);
379 if (gIt == idMats[2].end())
381 idMats[2][gId] = tmpStore;
385 ASSERTL1(tmpStore.size() == gIt->second.size(),
386 "Number of modes mismatch");
387 Vmath::Vadd(nFaceCoeffs*nFaceCoeffs, &gIt->second[0], 1,
388 &tmpStore[0], 1, &gIt->second[0], 1);
391 gidMeshIds[2][gId] = meshFaceId;
392 maxVertIds[4] = max(maxVertIds[4], meshFaceId);
393 maxVertIds[5] = max(maxVertIds[5], nFaceCoeffs);
396 cnt += exp->GetNcoeffs();
401 m_comm = expList->GetSession()->GetComm()->GetRowComm();
406 vector<NekDouble> storageBuf;
407 vector<long> globalToUniversal;
409 for (i = 0, cnt = 1; i < 3; ++i)
411 const int maxDofs = maxVertIds[2*i+1];
416 for (gIt = idMats[i].begin(); gIt != idMats[i].end(); ++gIt)
419 storageBuf.insert(storageBuf.end(),
420 gIt->second.begin(), gIt->second.end());
423 ASSERTL1(gidMeshIds[i].count(gIt->first) > 0,
424 "Unable to find global ID " +
425 boost::lexical_cast<string>(gIt->first) +
427 meshVertId = gidMeshIds[i][gIt->first];
429 for (j = 0; j < gIt->second.size(); ++j)
431 globalToUniversal.push_back(
432 cnt + meshVertId*maxDofs*maxDofs + j);
439 cnt += (maxVertIds[2*i]+1)*maxDofs*maxDofs;
442 ASSERTL1(storageBuf.size() == globalToUniversal.size(),
443 "Storage buffer and global to universal map size does "
447 storageBuf.size(), &storageBuf[0]);
449 globalToUniversal.size(), &globalToUniversal[0]);
457 1 + idMats[1].size() + idMats[2].size());
460 n_blks[0] = idMats[0].size();
465 for (i = 1; i < 3; ++i)
467 for (gIt = idMats[i].begin(); gIt != idMats[i].end(); ++gIt)
469 ASSERTL1(gidDofs[i].count(gIt->first) > 0,
470 "Unable to find number of degrees of freedom for "
471 "global ID " + boost::lexical_cast<string>(
474 n_blks[cnt++] = gidDofs[i][gIt->first];
489 for (gIt = idMats[0].begin(); gIt != idMats[0].end(); ++gIt, ++cnt)
491 (*vertMat)(cnt, cnt) = 1.0/storageData[cnt];
500 for (i = 1; i < 3; ++i)
502 for (gIt = idMats[i].begin(); gIt != idMats[i].end();
505 int nDofs = gidDofs[i][gIt->first];
510 for (j = 0; j < nDofs; ++j)
512 for (k = 0; k < nDofs; ++k)
514 (*tmp)(j,k) = storageData[k+j*nDofs + cnt];
521 m_blkMat->SetBlock(cnt2, cnt2, tmp);
541 boost::shared_ptr<MultiRegions::ExpList>
542 expList = ((
m_linsys.lock())->GetLocMat()).lock();
543 boost::shared_ptr<MultiRegions::ExpList> trace = expList->GetTrace();
551 int i, j, k, n, cnt, cnt2;
554 int nTrace = expList->GetTrace()->GetExpSize();
557 for (cnt = n = 0; n < nTrace; ++n)
564 cnt += trace->GetExp(n)->GetNcoeffs();
571 int maxTraceSize = -1;
572 map<int, int> arrayOffsets;
574 for (cnt = 0, n = nDir; n < nTrace; ++n)
576 int nCoeffs = trace->GetExp(n)->GetNcoeffs();
577 int nCoeffs2 = nCoeffs * nCoeffs;
579 arrayOffsets[n] = cnt;
582 if (maxTraceSize < nCoeffs2)
584 maxTraceSize = nCoeffs2;
589 m_comm = expList->GetSession()->GetComm()->GetRowComm();
596 for (cnt = n = 0; n < expList->GetExpSize(); ++n)
598 int elmt = expList->GetOffset_Elmt_Id(n);
599 locExpansion = expList->GetExp(elmt);
602 asmMap->GetElmtToTrace()[elmt];
605 loc_mat = (
m_linsys.lock())->GetStaticCondBlock(n);
606 bnd_mat = loc_mat->GetBlock(0,0);
608 int nFacets = locExpansion->GetNumBases() == 2 ?
609 locExpansion->GetNedges() : locExpansion->GetNfaces();
611 for (cnt2 = i = 0; i < nFacets; ++i)
613 int nCoeffs = elmtToTraceMap[i]->GetNcoeffs();
614 int elmtId = elmtToTraceMap[i]->GetElmtId ();
624 NekDouble *off = &tmpStore[arrayOffsets[elmtId]];
626 for (j = 0; j < nCoeffs; ++j)
628 NekDouble sign1 = asmMap->GetLocalToGlobalBndSign(
630 for (k = 0; k < nCoeffs; ++k)
632 NekDouble sign2 = asmMap->GetLocalToGlobalBndSign(
634 off[j*nCoeffs + k] +=
635 (*bnd_mat)(cnt2+j, cnt2+k) * sign1 * sign2;
646 for (cnt = 0, n = nDir; n < nTrace; ++n)
649 int nCoeffs = traceExp->GetNcoeffs();
650 int geomId = traceExp->GetGeom()->GetGlobalID();
652 for (i = 0; i < nCoeffs*nCoeffs; ++i)
654 uniIds[cnt++] = geomId * maxTraceSize + i + 1;
664 for (n = 0; n < nTrace - nDir; ++n)
666 n_blks[n] = trace->GetExp(n + nDir)->GetNcoeffs();
672 for (n = 0; n < nTrace - nDir; ++n)
674 int nCoeffs = trace->GetExp(n + nDir)->GetNcoeffs();
677 NekDouble *off = &tmpStore[arrayOffsets[n+nDir]];
679 for (i = 0; i < nCoeffs; ++i)
681 for (j = 0; j < nCoeffs; ++j)
683 (*tmp)(i,j) = off[i*nCoeffs + j];
702 int nNonDir = nGlobal-nDir;
virtual void v_BuildPreconditioner()
boost::shared_ptr< AssemblyMapDG > AssemblyMapDGSharedPtr
LibUtilities::CommSharedPtr m_comm
#define ASSERTL0(condition, msg)
virtual void v_DoPreconditioner(const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
Apply preconditioner to pInput and store the result in pOutput.
static void Gather(Nektar::Array< OneD, NekDouble > pU, gs_op pOp, gs_data *pGsh, Nektar::Array< OneD, NekDouble > pBuffer=NullNekDouble1DArray)
Performs a gather-scatter operation of the provided values.
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
std::map< int, vector< PeriodicEntity > > PeriodicMap
#define sign(a, b)
return the sign(b)*a
boost::shared_ptr< AssemblyMap > AssemblyMapSharedPtr
static gs_data * Init(const Nektar::Array< OneD, long > pId, const LibUtilities::CommSharedPtr &pComm)
Initialise Gather-Scatter map.
static std::string className
Name of class.
DNekBlkMatSharedPtr m_blkMat
PreconFactory & GetPreconFactory()
boost::shared_ptr< DNekMat > DNekMatSharedPtr
boost::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
StdRegions::Orientation DeterminePeriodicFaceOrient(StdRegions::Orientation faceOrient, StdRegions::Orientation perFaceOrient)
Determine relative orientation between two faces.
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
boost::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
pair< int, StdRegions::Orientation > DeterminePeriodicEdgeOrientId(int meshEdgeId, StdRegions::Orientation edgeOrient, const vector< PeriodicEntity > &periodicEdges)
Determine orientation of an edge to its periodic equivalents, as well as the ID of the representative...
void BlockPreconditionerHDG(void)
Construct a block preconditioner for the hybridized discontinuous Galerkin system.
static PreconditionerSharedPtr create(const boost::shared_ptr< GlobalLinSys > &plinsys, const boost::shared_ptr< AssemblyMap > &pLocToGloMap)
Creates an instance of this class.
Describe a linear system.
boost::shared_ptr< AssemblyMap > m_locToGloMap
StdRegions::MatrixType GetMatrixType() const
Return the matrix type.
virtual void v_InitObject()
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
boost::shared_ptr< Expansion > ExpansionSharedPtr
void BlockPreconditionerCG(void)
Construct a block preconditioner from for the continuous Galerkin system.
const boost::weak_ptr< GlobalLinSys > m_linsys
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
PreconditionerBlock(const boost::shared_ptr< GlobalLinSys > &plinsys, const AssemblyMapSharedPtr &pLocToGloMap)
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.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.