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 exp = expList->GetExp(n);
169 m_linsys.lock()->GetStaticCondBlock(n)->GetBlock(0,0);
173 for (i = 0; i < exp->GetNverts(); ++i)
175 meshVertId = exp->GetGeom()->GetVid(i);
176 int locId = exp->GetVertexMap(i);
180 cnt + locId) - nDirBnd;
191 NekDouble vertVal = (*schurMat)(locId,locId);
195 gIt = idMats[0].find(gId);
197 if (gIt == idMats[0].end())
200 idMats[0][gId] = vector<NekDouble>(1, vertVal);
206 gIt->second[0] += vertVal;
213 pIt = periodicVerts.find(meshVertId);
214 if (pIt != periodicVerts.end())
216 for (j = 0; j < pIt->second.size(); ++j)
218 meshVertId = min(meshVertId, pIt->second[j].id);
224 gidMeshIds[0][gId] = meshVertId;
225 maxVertIds[0] = max(maxVertIds[0], meshVertId);
231 for (i = 0; i < exp->GetNedges(); ++i)
233 meshEdgeId = exp->GetGeom()->GetEid(i);
238 exp->GetGeom()->GetEorient(i);
242 pIt = periodicEdges.find(meshEdgeId);
243 if (pIt != periodicEdges.end())
245 pair<int, StdRegions::Orientation> idOrient =
247 meshEdgeId, edgeOrient, pIt->second);
248 meshEdgeId = idOrient.first;
249 edgeOrient = idOrient.second;
255 exp->GetEdgeInteriorMap(i, edgeOrient, bmap, sign);
256 bmap2 = exp->GetEdgeInverseBoundaryMap(i);
259 const int nEdgeCoeffs = bmap.num_elements();
260 vector<NekDouble> tmpStore(nEdgeCoeffs*nEdgeCoeffs);
264 for (j = 0; j < nEdgeCoeffs; ++j)
286 for (k = 0; k < nEdgeCoeffs; ++k)
288 tmpStore[k+j*nEdgeCoeffs] =
289 sign1*sign[k]*(*schurMat)(bmap2[j], bmap2[k]);
298 gidDofs[1][gId] = nEdgeCoeffs;
301 gIt = idMats[1].find(gId);
302 if (gIt == idMats[1].end())
304 idMats[1][gId] = tmpStore;
308 ASSERTL1(tmpStore.size() == gIt->second.size(),
309 "Number of modes mismatch");
310 Vmath::Vadd(nEdgeCoeffs*nEdgeCoeffs, &gIt->second[0], 1,
311 &tmpStore[0], 1, &gIt->second[0], 1);
314 gidMeshIds[1][gId] = meshEdgeId;
315 maxVertIds[2] = max(maxVertIds[2], meshEdgeId);
316 maxVertIds[3] = max(maxVertIds[3], nEdgeCoeffs);
321 for (i = 0; i < exp->GetNfaces(); ++i)
323 meshFaceId = exp->GetGeom()->GetFid(i);
328 exp->GetGeom()->GetForient(i);
332 pIt = periodicFaces.find(meshFaceId);
333 if (pIt != periodicFaces.end())
335 meshFaceId = min(meshFaceId, pIt->second[0].id);
337 faceOrient, pIt->second[0].orient);
340 exp->GetFaceInteriorMap(i, faceOrient, bmap, sign);
341 bmap2 = exp->GetFaceInverseBoundaryMap(i);
344 const int nFaceCoeffs = bmap.num_elements();
345 vector<NekDouble> tmpStore(nFaceCoeffs*nFaceCoeffs);
349 for (j = 0; j < nFaceCoeffs; ++j)
366 for (k = 0; k < nFaceCoeffs; ++k)
368 tmpStore[k+j*nFaceCoeffs] =
369 sign1*sign[k]*(*schurMat)(bmap2[j], bmap2[k]);
378 gidDofs[2][gId] = nFaceCoeffs;
381 gIt = idMats[2].find(gId);
382 if (gIt == idMats[2].end())
384 idMats[2][gId] = tmpStore;
388 ASSERTL1(tmpStore.size() == gIt->second.size(),
389 "Number of modes mismatch");
390 Vmath::Vadd(nFaceCoeffs*nFaceCoeffs, &gIt->second[0], 1,
391 &tmpStore[0], 1, &gIt->second[0], 1);
394 gidMeshIds[2][gId] = meshFaceId;
395 maxVertIds[4] = max(maxVertIds[4], meshFaceId);
396 maxVertIds[5] = max(maxVertIds[5], nFaceCoeffs);
399 cnt += exp->GetNcoeffs();
404 m_comm = expList->GetSession()->GetComm()->GetRowComm();
409 vector<NekDouble> storageBuf;
410 vector<long> globalToUniversal;
412 for (i = 0, cnt = 1; i < 3; ++i)
414 const int maxDofs = maxVertIds[2*i+1];
419 for (gIt = idMats[i].begin(); gIt != idMats[i].end(); ++gIt)
422 storageBuf.insert(storageBuf.end(),
423 gIt->second.begin(), gIt->second.end());
426 ASSERTL1(gidMeshIds[i].count(gIt->first) > 0,
427 "Unable to find global ID " +
428 boost::lexical_cast<string>(gIt->first) +
430 meshVertId = gidMeshIds[i][gIt->first];
432 for (j = 0; j < gIt->second.size(); ++j)
434 globalToUniversal.push_back(
435 cnt + meshVertId*maxDofs*maxDofs + j);
442 cnt += (maxVertIds[2*i]+1)*maxDofs*maxDofs;
445 ASSERTL1(storageBuf.size() == globalToUniversal.size(),
446 "Storage buffer and global to universal map size does "
450 storageBuf.size(), &storageBuf[0]);
452 globalToUniversal.size(), &globalToUniversal[0]);
456 globalToUniversalMap,
m_comm,
457 expList->GetSession()->DefinesCmdLineArgument(
"verbose"));
462 1 + idMats[1].size() + idMats[2].size());
465 n_blks[0] = idMats[0].size();
470 for (i = 1; i < 3; ++i)
472 for (gIt = idMats[i].begin(); gIt != idMats[i].end(); ++gIt)
474 ASSERTL1(gidDofs[i].count(gIt->first) > 0,
475 "Unable to find number of degrees of freedom for "
476 "global ID " + boost::lexical_cast<string>(
479 n_blks[cnt++] = gidDofs[i][gIt->first];
494 for (gIt = idMats[0].begin(); gIt != idMats[0].end(); ++gIt, ++cnt)
496 (*vertMat)(cnt, cnt) = 1.0/storageData[cnt];
505 for (i = 1; i < 3; ++i)
507 for (gIt = idMats[i].begin(); gIt != idMats[i].end();
510 int nDofs = gidDofs[i][gIt->first];
515 for (j = 0; j < nDofs; ++j)
517 for (k = 0; k < nDofs; ++k)
519 (*tmp)(j,k) = storageData[k+j*nDofs + cnt];
526 m_blkMat->SetBlock(cnt2, cnt2, tmp);
546 boost::shared_ptr<MultiRegions::ExpList>
547 expList = ((
m_linsys.lock())->GetLocMat()).lock();
548 boost::shared_ptr<MultiRegions::ExpList> trace = expList->GetTrace();
556 int i, j, k, n, cnt, cnt2;
559 int nTrace = expList->GetTrace()->GetExpSize();
562 for (cnt = n = 0; n < nTrace; ++n)
569 cnt += trace->GetExp(n)->GetNcoeffs();
576 int maxTraceSize = -1;
577 map<int, int> arrayOffsets;
579 for (cnt = 0, n = nDir; n < nTrace; ++n)
581 int nCoeffs = trace->GetExp(n)->GetNcoeffs();
582 int nCoeffs2 = nCoeffs * nCoeffs;
584 arrayOffsets[n] = cnt;
587 if (maxTraceSize < nCoeffs2)
589 maxTraceSize = nCoeffs2;
594 m_comm = expList->GetSession()->GetComm()->GetRowComm();
601 for (cnt = n = 0; n < expList->GetExpSize(); ++n)
604 locExpansion = expList->GetExp(elmt);
607 asmMap->GetElmtToTrace()[elmt];
610 loc_mat = (
m_linsys.lock())->GetStaticCondBlock(n);
611 bnd_mat = loc_mat->GetBlock(0,0);
613 int nFacets = locExpansion->GetNumBases() == 2 ?
614 locExpansion->GetNedges() : locExpansion->GetNfaces();
616 for (cnt2 = i = 0; i < nFacets; ++i)
618 int nCoeffs = elmtToTraceMap[i]->GetNcoeffs();
619 int elmtId = elmtToTraceMap[i]->GetElmtId ();
629 NekDouble *off = &tmpStore[arrayOffsets[elmtId]];
631 for (j = 0; j < nCoeffs; ++j)
633 NekDouble sign1 = asmMap->GetLocalToGlobalBndSign(
635 for (k = 0; k < nCoeffs; ++k)
637 NekDouble sign2 = asmMap->GetLocalToGlobalBndSign(
639 off[j*nCoeffs + k] +=
640 (*bnd_mat)(cnt2+j, cnt2+k) * sign1 * sign2;
651 for (cnt = 0, n = nDir; n < nTrace; ++n)
654 int nCoeffs = traceExp->GetNcoeffs();
655 int geomId = traceExp->GetGeom()->GetGlobalID();
657 for (i = 0; i < nCoeffs*nCoeffs; ++i)
659 uniIds[cnt++] = geomId * maxTraceSize + i + 1;
666 expList->GetSession()->DefinesCmdLineArgument(
"verbose"));
671 for (n = 0; n < nTrace - nDir; ++n)
673 n_blks[n] = trace->GetExp(n + nDir)->GetNcoeffs();
679 for (n = 0; n < nTrace - nDir; ++n)
681 int nCoeffs = trace->GetExp(n + nDir)->GetNcoeffs();
684 NekDouble *off = &tmpStore[arrayOffsets[n+nDir]];
686 for (i = 0; i < nCoeffs; ++i)
688 for (j = 0; j < nCoeffs; ++j)
690 (*tmp)(i,j) = off[i*nCoeffs + j];
709 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
DNekBlkMatSharedPtr m_blkMat
PreconFactory & GetPreconFactory()
boost::shared_ptr< DNekMat > DNekMatSharedPtr
boost::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
static gs_data * Init(const Nektar::Array< OneD, long > pId, const LibUtilities::CommSharedPtr &pComm, bool verbose=true)
Initialise Gather-Scatter map.
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.