50 namespace MultiRegions
55 string PreconditionerBlock::className
58 PreconditionerBlock::create,
59 "Block Diagonal Preconditioning");
67 PreconditionerBlock::PreconditionerBlock(
68 const boost::shared_ptr<GlobalLinSys> &plinsys,
72 m_preconType(pLocToGloMap->GetPreconType()),
73 m_locToGloMap(pLocToGloMap)
82 "Solver type not valid");
123 int nel, i, j, k, n, cnt, gId;
124 int meshVertId, meshEdgeId, meshFaceId;
126 const int nExp = expList->GetExpSize();
127 const int nDirBnd =
m_locToGloMap->GetNumGlobalDirBndCoeffs();
131 PeriodicMap periodicVerts, periodicEdges, periodicFaces;
132 expList->GetPeriodicEntities(
133 periodicVerts, periodicEdges, periodicFaces);
140 vector<map<int, vector<NekDouble> > > idMats(3);
144 vector<map<int, int> > gidMeshIds(3);
148 vector<map<int, int> > gidDofs(3);
163 for (cnt = n = 0; n < nExp; ++n)
165 nel = expList->GetOffset_Elmt_Id(n);
166 exp = expList->GetExp(nel);
170 m_linsys.lock()->GetStaticCondBlock(n)->GetBlock(0,0);
174 for (i = 0; i < exp->GetNverts(); ++i)
176 meshVertId = exp->GetGeom()->GetVid(i);
177 int locId = exp->GetVertexMap(i);
181 cnt + locId) - nDirBnd;
192 NekDouble vertVal = (*schurMat)(locId,locId);
196 gIt = idMats[0].find(gId);
198 if (gIt == idMats[0].end())
201 idMats[0][gId] = vector<NekDouble>(1, vertVal);
207 gIt->second[0] += vertVal;
214 pIt = periodicVerts.find(meshVertId);
215 if (pIt != periodicVerts.end())
217 for (j = 0; j < pIt->second.size(); ++j)
219 meshVertId = min(meshVertId, pIt->second[j].id);
225 gidMeshIds[0][gId] = meshVertId;
226 maxVertIds[0] = max(maxVertIds[0], meshVertId);
232 for (i = 0; i < exp->GetNedges(); ++i)
234 meshEdgeId = exp->GetGeom()->GetEid(i);
239 exp->GetGeom()->GetEorient(i);
243 pIt = periodicEdges.find(meshEdgeId);
244 if (pIt != periodicEdges.end())
246 pair<int, StdRegions::Orientation> idOrient =
248 meshEdgeId, edgeOrient, pIt->second);
249 meshEdgeId = idOrient.first;
250 edgeOrient = idOrient.second;
256 exp->GetEdgeInteriorMap(i, edgeOrient, bmap, sign);
257 bmap2 = exp->GetEdgeInverseBoundaryMap(i);
260 const int nEdgeCoeffs = bmap.num_elements();
261 vector<NekDouble> tmpStore(nEdgeCoeffs*nEdgeCoeffs);
265 for (j = 0; j < nEdgeCoeffs; ++j)
287 for (k = 0; k < nEdgeCoeffs; ++k)
289 tmpStore[k+j*nEdgeCoeffs] =
290 sign1*sign[k]*(*schurMat)(bmap2[j], bmap2[k]);
299 gidDofs[1][gId] = nEdgeCoeffs;
302 gIt = idMats[1].find(gId);
303 if (gIt == idMats[1].end())
305 idMats[1][gId] = tmpStore;
309 ASSERTL1(tmpStore.size() == gIt->second.size(),
310 "Number of modes mismatch");
311 Vmath::Vadd(nEdgeCoeffs*nEdgeCoeffs, &gIt->second[0], 1,
312 &tmpStore[0], 1, &gIt->second[0], 1);
315 gidMeshIds[1][gId] = meshEdgeId;
316 maxVertIds[2] = max(maxVertIds[2], meshEdgeId);
317 maxVertIds[3] = max(maxVertIds[3], nEdgeCoeffs);
322 for (i = 0; i < exp->GetNfaces(); ++i)
324 meshFaceId = exp->GetGeom()->GetFid(i);
329 exp->GetGeom()->GetForient(i);
333 pIt = periodicFaces.find(meshFaceId);
334 if (pIt != periodicFaces.end())
336 meshFaceId = min(meshFaceId, pIt->second[0].id);
338 faceOrient, pIt->second[0].orient);
341 exp->GetFaceInteriorMap(i, faceOrient, bmap, sign);
342 bmap2 = exp->GetFaceInverseBoundaryMap(i);
345 const int nFaceCoeffs = bmap.num_elements();
346 vector<NekDouble> tmpStore(nFaceCoeffs*nFaceCoeffs);
350 for (j = 0; j < nFaceCoeffs; ++j)
367 for (k = 0; k < nFaceCoeffs; ++k)
369 tmpStore[k+j*nFaceCoeffs] =
370 sign1*sign[k]*(*schurMat)(bmap2[j], bmap2[k]);
379 gidDofs[2][gId] = nFaceCoeffs;
382 gIt = idMats[2].find(gId);
383 if (gIt == idMats[2].end())
385 idMats[2][gId] = tmpStore;
389 ASSERTL1(tmpStore.size() == gIt->second.size(),
390 "Number of modes mismatch");
391 Vmath::Vadd(nFaceCoeffs*nFaceCoeffs, &gIt->second[0], 1,
392 &tmpStore[0], 1, &gIt->second[0], 1);
395 gidMeshIds[2][gId] = meshFaceId;
396 maxVertIds[4] = max(maxVertIds[4], meshFaceId);
397 maxVertIds[5] = max(maxVertIds[5], nFaceCoeffs);
400 cnt += exp->GetNcoeffs();
405 m_comm = expList->GetSession()->GetComm()->GetRowComm();
410 vector<NekDouble> storageBuf;
411 vector<long> globalToUniversal;
413 for (i = 0, cnt = 1; i < 3; ++i)
415 const int maxDofs = maxVertIds[2*i+1];
420 for (gIt = idMats[i].begin(); gIt != idMats[i].end(); ++gIt)
423 storageBuf.insert(storageBuf.end(),
424 gIt->second.begin(), gIt->second.end());
427 ASSERTL1(gidMeshIds[i].count(gIt->first) > 0,
428 "Unable to find global ID " +
429 boost::lexical_cast<string>(gIt->first) +
431 meshVertId = gidMeshIds[i][gIt->first];
433 for (j = 0; j < gIt->second.size(); ++j)
435 globalToUniversal.push_back(
436 cnt + meshVertId*maxDofs*maxDofs + j);
443 cnt += (maxVertIds[2*i]+1)*maxDofs*maxDofs;
446 ASSERTL1(storageBuf.size() == globalToUniversal.size(),
447 "Storage buffer and global to universal map size does "
451 storageBuf.size(), &storageBuf[0]);
453 globalToUniversal.size(), &globalToUniversal[0]);
461 1 + idMats[1].size() + idMats[2].size());
464 n_blks[0] = idMats[0].size();
469 for (i = 1; i < 3; ++i)
471 for (gIt = idMats[i].begin(); gIt != idMats[i].end(); ++gIt)
473 ASSERTL1(gidDofs[i].count(gIt->first) > 0,
474 "Unable to find number of degrees of freedom for "
475 "global ID " + boost::lexical_cast<string>(
478 n_blks[cnt++] = gidDofs[i][gIt->first];
493 for (gIt = idMats[0].begin(); gIt != idMats[0].end(); ++gIt, ++cnt)
495 (*vertMat)(cnt, cnt) = 1.0/storageData[cnt];
504 for (i = 1; i < 3; ++i)
506 for (gIt = idMats[i].begin(); gIt != idMats[i].end();
509 int nDofs = gidDofs[i][gIt->first];
514 for (j = 0; j < nDofs; ++j)
516 for (k = 0; k < nDofs; ++k)
518 (*tmp)(j,k) = storageData[k+j*nDofs + cnt];
525 m_blkMat->SetBlock(cnt2, cnt2, tmp);
545 boost::shared_ptr<MultiRegions::ExpList>
546 expList = ((
m_linsys.lock())->GetLocMat()).lock();
547 boost::shared_ptr<MultiRegions::ExpList> trace = expList->GetTrace();
555 int i, j, k, n, cnt, cnt2;
558 int nTrace = expList->GetTrace()->GetExpSize();
561 for (cnt = n = 0; n < nTrace; ++n)
568 cnt += trace->GetExp(n)->GetNcoeffs();
575 int maxTraceSize = -1;
576 map<int, int> arrayOffsets;
578 for (cnt = 0, n = nDir; n < nTrace; ++n)
580 int nCoeffs = trace->GetExp(n)->GetNcoeffs();
581 int nCoeffs2 = nCoeffs * nCoeffs;
583 arrayOffsets[n] = cnt;
586 if (maxTraceSize < nCoeffs2)
588 maxTraceSize = nCoeffs2;
593 m_comm = expList->GetSession()->GetComm()->GetRowComm();
600 for (cnt = n = 0; n < expList->GetExpSize(); ++n)
602 int elmt = expList->GetOffset_Elmt_Id(n);
603 locExpansion = expList->GetExp(elmt);
606 asmMap->GetElmtToTrace()[elmt];
609 loc_mat = (
m_linsys.lock())->GetStaticCondBlock(n);
610 bnd_mat = loc_mat->GetBlock(0,0);
612 int nFacets = locExpansion->GetNumBases() == 2 ?
613 locExpansion->GetNedges() : locExpansion->GetNfaces();
615 for (cnt2 = i = 0; i < nFacets; ++i)
617 int nCoeffs = elmtToTraceMap[i]->GetNcoeffs();
618 int elmtId = elmtToTraceMap[i]->GetElmtId ();
628 NekDouble *off = &tmpStore[arrayOffsets[elmtId]];
630 for (j = 0; j < nCoeffs; ++j)
632 NekDouble sign1 = asmMap->GetLocalToGlobalBndSign(
634 for (k = 0; k < nCoeffs; ++k)
636 NekDouble sign2 = asmMap->GetLocalToGlobalBndSign(
638 off[j*nCoeffs + k] +=
639 (*bnd_mat)(cnt2+j, cnt2+k) * sign1 * sign2;
650 for (cnt = 0, n = nDir; n < nTrace; ++n)
653 int nCoeffs = traceExp->GetNcoeffs();
654 int geomId = traceExp->GetGeom()->GetGlobalID();
656 for (i = 0; i < nCoeffs*nCoeffs; ++i)
658 uniIds[cnt++] = geomId * maxTraceSize + i + 1;
668 for (n = 0; n < nTrace - nDir; ++n)
670 n_blks[n] = trace->GetExp(n + nDir)->GetNcoeffs();
676 for (n = 0; n < nTrace - nDir; ++n)
678 int nCoeffs = trace->GetExp(n + nDir)->GetNcoeffs();
681 NekDouble *off = &tmpStore[arrayOffsets[n+nDir]];
683 for (i = 0; i < nCoeffs; ++i)
685 for (j = 0; j < nCoeffs; ++j)
687 (*tmp)(i,j) = off[i*nCoeffs + j];
706 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.
#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.
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< DNekScalBlkMat > DNekScalBlkMatSharedPtr
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
boost::shared_ptr< Expansion > ExpansionSharedPtr
void BlockPreconditionerHDG(void)
Construct a block preconditioner for the hybridized discontinuous Galerkin system.
std::map< int, std::vector< PeriodicEntity > > PeriodicMap
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
void BlockPreconditionerCG(void)
Construct a block preconditioner from for the continuous Galerkin system.
const boost::weak_ptr< GlobalLinSys > m_linsys
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...
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
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.