Nektar++
Public Member Functions | Static Public Member Functions | Static Public Attributes | Protected Member Functions | Protected Attributes | Private Member Functions | List of all members
Nektar::MultiRegions::PreconditionerBlock Class Reference

#include <PreconditionerBlock.h>

Inheritance diagram for Nektar::MultiRegions::PreconditionerBlock:
[legend]

Public Member Functions

 PreconditionerBlock (const std::shared_ptr< GlobalLinSys > &plinsys, const AssemblyMapSharedPtr &pLocToGloMap)
 
 ~PreconditionerBlock () override
 
- Public Member Functions inherited from Nektar::MultiRegions::Preconditioner
 Preconditioner (const std::shared_ptr< GlobalLinSys > &plinsys, const AssemblyMapSharedPtr &pLocToGloMap)
 
virtual ~Preconditioner ()
 
void DoPreconditioner (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const bool &IsLocal=false)
 
void DoAssembleLoc (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const bool &ZeroDir)
 Apply an assembly and scatter back to lcoal array. More...
 
void DoPreconditionerWithNonVertOutput (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const Array< OneD, NekDouble > &pNonVertOutput, Array< OneD, NekDouble > &pVertForce=NullNekDouble1DArray)
 
void DoTransformBasisToLowEnergy (Array< OneD, NekDouble > &pInOut)
 
void DoTransformCoeffsFromLowEnergy (Array< OneD, NekDouble > &pInOut)
 
void DoTransformCoeffsToLowEnergy (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
 
void DoTransformBasisFromLowEnergy (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
 
void BuildPreconditioner ()
 
void InitObject ()
 
Array< OneD, NekDoubleAssembleStaticCondGlobalDiagonals ()
 Performs global assembly of diagonal entries to global Schur complement matrix. More...
 
const DNekScalBlkMatSharedPtrGetBlockTransformedSchurCompl () const
 
const DNekScalBlkMatSharedPtrGetBlockCMatrix () const
 
const DNekScalBlkMatSharedPtrGetBlockInvDMatrix () const
 
const DNekScalBlkMatSharedPtrGetBlockSchurCompl () const
 
const DNekScalBlkMatSharedPtrGetBlockTransformationMatrix () const
 
const DNekScalBlkMatSharedPtrGetBlockTransposedTransformationMatrix () const
 
DNekScalMatSharedPtr TransformedSchurCompl (int offset, int bndoffset, const std::shared_ptr< DNekScalMat > &loc_mat)
 

Static Public Member Functions

static PreconditionerSharedPtr create (const std::shared_ptr< GlobalLinSys > &plinsys, const std::shared_ptr< AssemblyMap > &pLocToGloMap)
 Creates an instance of this class. More...
 

Static Public Attributes

static std::string className
 Name of class. More...
 

Protected Member Functions

void v_InitObject () override
 
void v_DoPreconditioner (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const bool &isLocal=false) override
 Apply preconditioner to pInput and store the result in pOutput. More...
 
void v_BuildPreconditioner () override
 
- Protected Member Functions inherited from Nektar::MultiRegions::Preconditioner
virtual DNekScalMatSharedPtr v_TransformedSchurCompl (int offset, int bndoffset, const std::shared_ptr< DNekScalMat > &loc_mat)
 Get block elemental transposed transformation matrix \(\mathbf{R}^{T}\). More...
 
virtual void v_InitObject ()
 
virtual void v_DoPreconditioner (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const bool &isLocal=false)
 Apply a preconditioner to the conjugate gradient method. More...
 
virtual void v_DoPreconditionerWithNonVertOutput (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const Array< OneD, NekDouble > &pNonVertOutput, Array< OneD, NekDouble > &pVertForce)
 Apply a preconditioner to the conjugate gradient method with an output for non-vertex degrees of freedom. More...
 
virtual void v_DoTransformBasisToLowEnergy (Array< OneD, NekDouble > &pInOut)
 Transform from original basis to low energy basis. More...
 
virtual void v_DoTransformCoeffsFromLowEnergy (Array< OneD, NekDouble > &pInOut)
 Transform from low energy coeffs to orignal basis. More...
 
virtual void v_DoTransformCoeffsToLowEnergy (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
 Multiply by the block inverse transformation matrix. More...
 
virtual void v_DoTransformBasisFromLowEnergy (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
 Multiply by the block transposed inverse transformation matrix. More...
 
virtual void v_BuildPreconditioner ()
 

Protected Attributes

bool m_isFull
 
DNekBlkMatSharedPtr m_blkMat
 
Gs::gs_datam_gs_mapping
 
- Protected Attributes inherited from Nektar::MultiRegions::Preconditioner
const std::weak_ptr< GlobalLinSysm_linsys
 
std::string m_preconType
 
DNekMatSharedPtr m_preconditioner
 
std::weak_ptr< AssemblyMapm_locToGloMap
 

Private Member Functions

void BlockPreconditionerCG (void)
 Construct a block preconditioner from \(\mathbf{S}_{1}\) for the continuous Galerkin system. More...
 
void BlockPreconditionerHDG (void)
 Construct a block preconditioner for the hybridized discontinuous Galerkin system. More...
 

Detailed Description

Definition at line 48 of file PreconditionerBlock.h.

Constructor & Destructor Documentation

◆ PreconditionerBlock()

Nektar::MultiRegions::PreconditionerBlock::PreconditionerBlock ( const std::shared_ptr< GlobalLinSys > &  plinsys,
const AssemblyMapSharedPtr pLocToGloMap 
)

Definition at line 61 of file PreconditionerBlock.cpp.

64 : Preconditioner(plinsys, pLocToGloMap)
65{
66}
Preconditioner(const std::shared_ptr< GlobalLinSys > &plinsys, const AssemblyMapSharedPtr &pLocToGloMap)

◆ ~PreconditionerBlock()

Nektar::MultiRegions::PreconditionerBlock::~PreconditionerBlock ( )
inlineoverride

Definition at line 71 of file PreconditionerBlock.h.

72 {
74 }
static void Free(gs_data *pGsh)
Deallocates GSLib mapping data without finalising MPI.
Definition: GsLib.hpp:263

References Gs::Free(), and m_gs_mapping.

Member Function Documentation

◆ BlockPreconditionerCG()

void Nektar::MultiRegions::PreconditionerBlock::BlockPreconditionerCG ( void  )
private

Construct a block preconditioner from \(\mathbf{S}_{1}\) for the continuous Galerkin system.

The preconditioner for the statically-condensed system is defined as:

\[\mathbf{M}^{-1}=\left[\begin{array}{ccc} \mathrm{Diag}[(\mathbf{S_{1}})_{vv}] & & \\ & (\mathbf{S}_{1})_{eb} & \\ & & (\mathbf{S}_{1})_{fb} \end{array}\right] \]

where \(\mathbf{S}_{1}\) is the local Schur complement matrix for each element and the subscript denotes the portion of the Schur complement associated with the vertex (vv), edge (eb) and face blocks (fb) respectively.

In case of the full system \(\mathbf{A}\), the preconditioner includes the interior blocks as well. The full block preconditioner is defined as:

\[\mathbf{M}^{-1}=\left[\begin{array}{cccc} \mathrm{Diag}[(\mathbf{A})_{vv}] & & & \\ & (\mathbf{A})_{eb} & & \\ & & (\mathbf{A})_{fb} & \\ & & & (\mathbf{A})_{int} \end{array}\right] \]

where the interior portion is added at the end following the ordering inside data structures.

Process boundary matrices

Definition at line 123 of file PreconditionerBlock.cpp.

124{
125 ExpListSharedPtr expList = m_linsys.lock()->GetLocMat().lock();
128 DNekScalMatSharedPtr bnd_mat;
129
130 int i, j, k, n, cnt, gId;
131 int meshVertId, meshEdgeId, meshFaceId;
132
133 auto asmMap = m_locToGloMap.lock();
134
135 const int nExp = expList->GetExpSize();
136 const int nDirBnd = asmMap->GetNumGlobalDirBndCoeffs();
137
138 // Grab periodic geometry information.
139 PeriodicMap periodicVerts, periodicEdges, periodicFaces;
140 expList->GetPeriodicEntities(periodicVerts, periodicEdges, periodicFaces);
141
142 // The vectors below are of size 3(+1) to have separate storage for
143 // vertices, edges and faces.
144 // For the full matrix approach additional storage for the
145 // interior matrices is required
146 int nStorage = m_isFull ? 4 : 3;
147
148 // Maps from geometry ID to the matrix representing the extracted
149 // portion of S_1. For example idMats[2] folds the S_1 face blocks.
150 vector<map<int, vector<NekDouble>>> idMats(nStorage);
151
152 // Maps from the global ID, as obtained from AssemblyMapCG's
153 // localToGlobalMap, to the geometry ID.
154 vector<map<int, int>> gidMeshIds(3);
155
156 // Maps from the global ID to the number of degrees of freedom for
157 // this geometry object.
158 vector<map<int, int>> gidDofs(nStorage);
159
160 // Array containing maximum information needed for the universal
161 // numbering later. For i = 0,1,2 for each geometry dimension:
162 //
163 // maxVertIds[2*i] = maximum geometry ID at dimension i
164 // maxVertIds[2*i+1] = maximum number of degrees of freedom for all
165 // elements of dimension i.
166 Array<OneD, int> maxVertIds(6, -1);
167
168 // Figure out mapping from each elemental contribution to offset in
169 // (vert,edge,face) triples.
170 for (cnt = n = 0; n < nExp; ++n)
171 {
172 exp = expList->GetExp(n);
173
174 // Grab reference to
175 // StaticCond: Local Schur complement matrix
176 // Full: Boundary and interior matrix
177 DNekScalMatSharedPtr schurMat;
178
179 if (m_isFull)
180 {
181 // Get local matrix
182 auto tmpMat = m_linsys.lock()->GetBlock(n);
183
184 // Get number of boundary dofs
185 int nBndDofs = exp->NumBndryCoeffs();
186 int nIntDofs = exp->GetNcoeffs() - nBndDofs;
187
188 /// Process boundary matrices
189 // Get boundary map
190 Array<OneD, unsigned int> bndMap(nBndDofs);
191 exp->GetBoundaryMap(bndMap);
192
193 // Create temporary matrices
194 DNekMatSharedPtr bndryMat =
196
197 // Extract boundary and send to StaticCond (Schur)
198 // framework
199 for (int i = 0; i < nBndDofs; ++i)
200 {
201 for (int j = 0; j < nBndDofs; ++j)
202 {
203 (*bndryMat)(i, j) = (*tmpMat)(bndMap[i], bndMap[j]);
204 }
205 }
206
207 schurMat =
209
210 // Extract interior matrix and save for block matrix setting
211 if (nIntDofs)
212 {
213 Array<OneD, unsigned int> intMap(nIntDofs);
214 exp->GetInteriorMap(intMap);
215
216 vector<NekDouble> tmpStore(nIntDofs * nIntDofs);
217 for (int i = 0; i < nIntDofs; ++i)
218 {
219 for (int j = 0; j < nIntDofs; ++j)
220 {
221 tmpStore[j + i * nIntDofs] =
222 (*tmpMat)(intMap[i], intMap[j]);
223 }
224 }
225
226 // Save interior matrices and nDofs
227 idMats[3][n] = tmpStore;
228 gidDofs[3][n] = nIntDofs;
229 // then intMat gets added to a block diagonal matrix that
230 // handles the extra interior dofs no need to do any
231 // communication for that because all interior dofs are
232 // local to each element
233 }
234 }
235 else
236 {
237 schurMat = m_linsys.lock()->GetStaticCondBlock(n)->GetBlock(0, 0);
238 }
239
240 // Process vertices to extract relevant portion of the Schur
241 // complement matrix.
242 for (i = 0; i < exp->GetNverts(); ++i)
243 {
244 meshVertId = exp->GetGeom()->GetVid(i);
245 int locId = exp->GetVertexMap(i);
246
247 // Get the global ID of this vertex.
248 gId = asmMap->GetLocalToGlobalMap(cnt + locId) - nDirBnd;
249
250 // Ignore all Dirichlet vertices.
251 if (gId < 0)
252 {
253 continue;
254 }
255
256 gidDofs[0][gId] = 1;
257
258 // Extract vertex value from Schur complement matrix.
259 NekDouble vertVal = (*schurMat)(locId, locId);
260
261 // See if we have processed this vertex from another
262 // element.
263 auto gIt = idMats[0].find(gId);
264
265 if (gIt == idMats[0].end())
266 {
267 // If not then put our 'matrix' inside idMats.
268 idMats[0][gId] = vector<NekDouble>(1, vertVal);
269 }
270 else
271 {
272 // Otherwise combine with the value that is already
273 // there (i.e. do assembly on this degree of freedom).
274 gIt->second[0] += vertVal;
275 }
276
277 // Now check to see if the vertex is periodic. If it is,
278 // then we change meshVertId to be the minimum of all the
279 // other periodic vertices, so that we don't end up
280 // duplicating the matrix in our final block matrix.
281 auto pIt = periodicVerts.find(meshVertId);
282 if (pIt != periodicVerts.end())
283 {
284 for (j = 0; j < pIt->second.size(); ++j)
285 {
286 meshVertId = min(meshVertId, pIt->second[j].id);
287 }
288 }
289
290 // Finally record the other information we need into the
291 // other matrices.
292 gidMeshIds[0][gId] = meshVertId;
293 maxVertIds[0] = max(maxVertIds[0], meshVertId);
294 maxVertIds[1] = 1;
295 }
296
297 // Process edges. This logic is mostly the same as the previous
298 // block.
299 for (i = 0; i < exp->GetGeom()->GetNumEdges(); ++i)
300 {
301 meshEdgeId = exp->GetGeom()->GetEid(i);
302
303 Array<OneD, unsigned int> bmap, bmap2;
304 Array<OneD, int> sign;
305 StdRegions::Orientation edgeOrient = exp->GetGeom()->GetEorient(i);
306
307 // Check if this edge is periodic. We may need to flip
308 // orientation if it is.
309 auto pIt = periodicEdges.find(meshEdgeId);
310 if (pIt != periodicEdges.end())
311 {
312 pair<int, StdRegions::Orientation> idOrient =
313 DeterminePeriodicEdgeOrientId(meshEdgeId, edgeOrient,
314 pIt->second);
315 meshEdgeId = idOrient.first;
316 edgeOrient = idOrient.second;
317 }
318
319 // Grab edge interior map, and the edge inverse boundary
320 // map, so that we can extract this edge from the Schur
321 // complement matrix.
322 if (exp->GetGeom()->GetNumFaces()) // 3D Element calls
323 {
324 exp->as<LocalRegions::Expansion3D>()
325 ->GetEdgeInteriorToElementMap(i, bmap, sign, edgeOrient);
326 bmap2 = exp->as<LocalRegions::Expansion3D>()
327 ->GetEdgeInverseBoundaryMap(i);
328 }
329 else
330 {
331 exp->GetTraceInteriorToElementMap(i, bmap, sign, edgeOrient);
332 bmap2 = exp->as<LocalRegions::Expansion2D>()
333 ->GetTraceInverseBoundaryMap(i);
334 }
335
336 // Check that edge has edge-interior dofs
337 const int nEdgeCoeffs = bmap.size();
338 if (nEdgeCoeffs == 0)
339 {
340 continue;
341 }
342
343 // Allocate temporary storage for the extracted edge matrix.
344 vector<NekDouble> tmpStore(nEdgeCoeffs * nEdgeCoeffs);
345
346 gId = asmMap->GetLocalToGlobalMap(cnt + bmap[0]);
347
348 for (j = 0; j < nEdgeCoeffs; ++j)
349 {
350 // We record the minimum ID from the edge for our
351 // maps. This follows the logic that the assembly map
352 // ordering will always give us a contiguous ordering of
353 // global degrees of freedom for edge interior
354 // coefficients.
355 gId = min(gId,
356 asmMap->GetLocalToGlobalMap(cnt + bmap[j]) - nDirBnd);
357
358 // Ignore Dirichlet edges.
359 if (gId < 0)
360 {
361 continue;
362 }
363
364 const NekDouble sign1 = sign[j];
365
366 // Extract this edge, along with sign array for assembly
367 // later.
368 for (k = 0; k < nEdgeCoeffs; ++k)
369 {
370 tmpStore[k + j * nEdgeCoeffs] =
371 sign1 * sign[k] * (*schurMat)(bmap2[j], bmap2[k]);
372 }
373 }
374
375 if (gId < 0)
376 {
377 continue;
378 }
379
380 gidDofs[1][gId] = nEdgeCoeffs;
381
382 // Assemble this edge matrix with another one, if it exists.
383 auto gIt = idMats[1].find(gId);
384 if (gIt == idMats[1].end())
385 {
386 idMats[1][gId] = tmpStore;
387 }
388 else
389 {
390 ASSERTL1(tmpStore.size() == gIt->second.size(),
391 "Number of modes mismatch");
392 Vmath::Vadd(nEdgeCoeffs * nEdgeCoeffs, &gIt->second[0], 1,
393 &tmpStore[0], 1, &gIt->second[0], 1);
394 }
395
396 gidMeshIds[1][gId] = meshEdgeId;
397 maxVertIds[2] = max(maxVertIds[2], meshEdgeId);
398 maxVertIds[3] = max(maxVertIds[3], nEdgeCoeffs);
399 }
400
401 // Process faces. This logic is mostly the same as the previous
402 // block.
403 for (i = 0; i < exp->GetGeom()->GetNumFaces(); ++i)
404 {
405 meshFaceId = exp->GetGeom()->GetFid(i);
406
407 Array<OneD, unsigned int> bmap, bmap2;
408 Array<OneD, int> sign;
409 StdRegions::Orientation faceOrient = exp->GetGeom()->GetForient(i);
410
411 // Check if this face is periodic. We may need to flip
412 // orientation if it is.
413 auto pIt = periodicFaces.find(meshFaceId);
414 if (pIt != periodicFaces.end())
415 {
416 meshFaceId = min(meshFaceId, pIt->second[0].id);
417 faceOrient = DeterminePeriodicFaceOrient(faceOrient,
418 pIt->second[0].orient);
419 }
420
421 // Trace Interior to Element Map
422 exp->GetTraceInteriorToElementMap(i, bmap, sign, faceOrient);
423 bmap2 = exp->as<LocalRegions::Expansion3D>()
424 ->GetTraceInverseBoundaryMap(i);
425
426 // Check that face has face-interior dofs
427 const int nFaceCoeffs = bmap.size();
428 if (nFaceCoeffs == 0)
429 {
430 continue;
431 }
432
433 // Allocate temporary storage for the extracted face matrix.
434 vector<NekDouble> tmpStore(nFaceCoeffs * nFaceCoeffs);
435
436 gId = asmMap->GetLocalToGlobalMap(cnt + bmap[0]);
437
438 for (j = 0; j < nFaceCoeffs; ++j)
439 {
440 gId = min(gId,
441 asmMap->GetLocalToGlobalMap(cnt + bmap[j]) - nDirBnd);
442
443 // Ignore Dirichlet faces.
444 if (gId < 0)
445 {
446 continue;
447 }
448
449 const NekDouble sign1 = sign[j];
450
451 // Extract this face, along with sign array for assembly
452 // later.
453 for (k = 0; k < nFaceCoeffs; ++k)
454 {
455 tmpStore[k + j * nFaceCoeffs] =
456 sign1 * sign[k] * (*schurMat)(bmap2[j], bmap2[k]);
457 }
458 }
459
460 if (gId < 0)
461 {
462 continue;
463 }
464
465 gidDofs[2][gId] = nFaceCoeffs;
466
467 // Assemble this face matrix with another one, if it exists.
468 auto gIt = idMats[2].find(gId);
469 if (gIt == idMats[2].end())
470 {
471 idMats[2][gId] = tmpStore;
472 }
473 else
474 {
475 ASSERTL1(tmpStore.size() == gIt->second.size(),
476 "Number of modes mismatch");
477 Vmath::Vadd(nFaceCoeffs * nFaceCoeffs, &gIt->second[0], 1,
478 &tmpStore[0], 1, &gIt->second[0], 1);
479 }
480
481 gidMeshIds[2][gId] = meshFaceId;
482 maxVertIds[4] = max(maxVertIds[4], meshFaceId);
483 maxVertIds[5] = max(maxVertIds[5], nFaceCoeffs);
484 }
485
486 cnt += exp->GetNcoeffs();
487 }
488
489 // Perform a reduction to find maximum vertex, edge and face
490 // geometry IDs.
492 expList->GetSession()->GetComm()->GetRowComm();
493 vComm->AllReduce(maxVertIds, LibUtilities::ReduceMax);
494
495 // Concatenate all matrices into contiguous storage and figure out
496 // universal ID numbering.
497 vector<NekDouble> storageBuf;
498 vector<long> globalToUniversal;
499
500 for (i = 0, cnt = 1; i < 3; ++i)
501 {
502 const int maxDofs = maxVertIds[2 * i + 1];
503
504 // Note that iterating over the map uses the property that keys
505 // are accessed in order of ascending order, putting everything
506 // in the right order for the global system.
507 for (auto &gIt : idMats[i])
508 {
509 // Copy matrix into storage.
510 storageBuf.insert(storageBuf.end(), gIt.second.begin(),
511 gIt.second.end());
512
513 // Get mesh ID from global ID number.
514 ASSERTL1(gidMeshIds[i].count(gIt.first) > 0,
515 "Unable to find global ID " + std::to_string(gIt.first) +
516 " inside map");
517 meshVertId = gidMeshIds[i][gIt.first];
518
519 for (j = 0; j < gIt.second.size(); ++j)
520 {
521 globalToUniversal.push_back(cnt +
522 meshVertId * maxDofs * maxDofs + j);
523 }
524
525 // Free up the temporary storage.
526 gIt.second.clear();
527 }
528
529 cnt += (maxVertIds[2 * i] + 1) * maxDofs * maxDofs;
530 }
531
532 ASSERTL1(storageBuf.size() == globalToUniversal.size(),
533 "Storage buffer and global to universal map size does "
534 "not match");
535
536 Array<OneD, NekDouble> storageData(storageBuf.size(), &storageBuf[0]);
537 Array<OneD, long> globalToUniversalMap(globalToUniversal.size(),
538 &globalToUniversal[0]);
539
540 // Initialise Gs mapping
542 Gs::Init(globalToUniversalMap, vComm,
543 expList->GetSession()->DefinesCmdLineArgument("verbose"));
544
545 // Use GSlib to assemble data between processors.
546 Gs::Gather(storageData, Gs::gs_add, m_gs_mapping);
547
548 // Figure out what storage we need in the block matrix.
549 int nblksSize = 1 + idMats[1].size() + idMats[2].size();
550 if (m_isFull)
551 {
552 nblksSize += idMats[3].size();
553 }
554 Array<OneD, unsigned int> n_blks(nblksSize);
555
556 // Vertex block is a diagonal matrix.
557 n_blks[0] = idMats[0].size();
558
559 // Now extract number of rows in each edge and face block from the
560 // gidDofs map.
561 cnt = 1;
562 for (i = 1; i < nStorage; ++i)
563 {
564 for (auto &gIt : idMats[i])
565 {
566 ASSERTL1(gidDofs[i].count(gIt.first) > 0,
567 "Unable to find number of degrees of freedom for "
568 "global ID " +
569 std::to_string(gIt.first));
570
571 // Check for non-zero DoF count
572 if (gidDofs[i][gIt.first])
573 {
574 n_blks[cnt++] = gidDofs[i][gIt.first];
575 }
576 }
577 }
578
579 // Allocate storage for the block matrix.
580 m_blkMat =
582
583 // We deal with the vertex matrix separately since all vertices can
584 // be condensed into a single, block-diagonal matrix.
586 n_blks[0], n_blks[0], 0.0, eDIAGONAL);
587
588 // Fill the vertex matrix with the inverse of each vertex value.
589 cnt = 0;
590 for (auto gIt = idMats[0].begin(); gIt != idMats[0].end(); ++gIt, ++cnt)
591 {
592 (*vertMat)(cnt, cnt) = 1.0 / storageData[cnt];
593 }
594
595 // Put the vertex matrix in the block matrix.
596 m_blkMat->SetBlock(0, 0, vertMat);
597
598 // Finally, grab the matrices from the block storage, invert them
599 // and place them in the correct position inside the block matrix.
600 int cnt2 = 1;
601 for (i = 1; i < nStorage; ++i)
602 {
603 for (auto &gIt : idMats[i])
604 {
605 int nDofs = gidDofs[i][gIt.first];
606
607 // Skip, if zero DoFs
608 if (nDofs == 0)
609 {
610 continue;
611 }
612
613 DNekMatSharedPtr tmp =
615
616 for (j = 0; j < nDofs; ++j)
617 {
618 for (k = 0; k < nDofs; ++k)
619 {
620 // Interior matrices do not need mapping
621 // and are stored in idMats[3]
622 if (m_isFull && i == 3)
623 {
624 (*tmp)(j, k) = gIt.second[k + j * nDofs];
625 }
626 // Boundary matrices
627 else
628 {
629 (*tmp)(j, k) = storageData[k + j * nDofs + cnt];
630 }
631 }
632 }
633
634 cnt += nDofs * nDofs;
635
636 tmp->Invert();
637 m_blkMat->SetBlock(cnt2, cnt2, tmp);
638 ++cnt2;
639 }
640 }
641}
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:242
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:47
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
const std::weak_ptr< GlobalLinSys > m_linsys
std::weak_ptr< AssemblyMap > m_locToGloMap
static gs_data * Init(const Nektar::Array< OneD, long > &pId, const LibUtilities::CommSharedPtr &pComm, bool verbose=true)
Initialise Gather-Scatter map.
Definition: GsLib.hpp:190
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.
Definition: GsLib.hpp:278
@ gs_add
Definition: GsLib.hpp:60
std::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:55
std::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:66
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
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...
std::map< int, std::vector< PeriodicEntity > > PeriodicMap
StdRegions::Orientation DeterminePeriodicFaceOrient(StdRegions::Orientation faceOrient, StdRegions::Orientation perFaceOrient)
Determine relative orientation between two faces.
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
std::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:79
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:75
double NekDouble
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.
Definition: Vmath.hpp:180

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL1, Nektar::MultiRegions::DeterminePeriodicEdgeOrientId(), Nektar::MultiRegions::DeterminePeriodicFaceOrient(), Nektar::eDIAGONAL, Gs::Gather(), Gs::gs_add, Gs::Init(), m_blkMat, m_gs_mapping, m_isFull, Nektar::MultiRegions::Preconditioner::m_linsys, Nektar::MultiRegions::Preconditioner::m_locToGloMap, Nektar::LibUtilities::ReduceMax, sign, and Vmath::Vadd().

Referenced by v_BuildPreconditioner().

◆ BlockPreconditionerHDG()

void Nektar::MultiRegions::PreconditionerBlock::BlockPreconditionerHDG ( void  )
private

Construct a block preconditioner for the hybridized discontinuous Galerkin system.

This system is built in a similar fashion to its continuous variant found in PreconditionerBlock::BlockPreconditionerCG. In this setting however, the matrix is constructed as:

\[ M^{-1} = \mathrm{Diag}[ (\mathbf{S_{1}})_{f}^{-1} ] \]

where each matrix is the Schur complement system restricted to a single face of the trace system.

Definition at line 656 of file PreconditionerBlock.cpp.

657{
658 std::shared_ptr<MultiRegions::ExpList> expList =
659 ((m_linsys.lock())->GetLocMat()).lock();
660 std::shared_ptr<MultiRegions::ExpList> trace = expList->GetTrace();
663 DNekScalMatSharedPtr bnd_mat;
664
666 std::dynamic_pointer_cast<AssemblyMapDG>(m_locToGloMap.lock());
667
668 int i, j, k, n, cnt, cnt2;
669
670 // Figure out number of Dirichlet trace elements
671 int nTrace = expList->GetTrace()->GetExpSize();
672 int nDir = asmMap->GetNumGlobalDirBndCoeffs();
673
674 for (cnt = n = 0; n < nTrace; ++n)
675 {
676 if (cnt >= nDir)
677 {
678 break;
679 }
680
681 cnt += trace->GetExp(n)->GetNcoeffs();
682 }
683
684 nDir = n;
685
686 // Allocate storage for block matrix. Need number of unique faces in
687 // trace space.
688 int maxTraceSize = -1;
689 map<int, int> arrayOffsets;
690
691 for (cnt = 0, n = nDir; n < nTrace; ++n)
692 {
693 int nCoeffs = trace->GetExp(n)->GetNcoeffs();
694 int nCoeffs2 = nCoeffs * nCoeffs;
695
696 arrayOffsets[n] = cnt;
697 cnt += nCoeffs2;
698
699 if (maxTraceSize < nCoeffs2)
700 {
701 maxTraceSize = nCoeffs2;
702 }
703 }
704
705 // Find maximum trace size.
707 expList->GetSession()->GetComm()->GetRowComm();
708 vComm->AllReduce(maxTraceSize, LibUtilities::ReduceMax);
709
710 // Zero matrix storage.
711 Array<OneD, NekDouble> tmpStore(cnt, 0.0);
712
713 // Assemble block matrices for each trace element.
714 for (cnt = n = 0; n < expList->GetExpSize(); ++n)
715 {
716 int elmt = n;
717 locExpansion = expList->GetExp(elmt);
718
719 Array<OneD, LocalRegions::ExpansionSharedPtr> &elmtToTraceMap =
720 asmMap->GetElmtToTrace()[elmt];
721
722 // Block matrix (lambda)
723 loc_mat = (m_linsys.lock())->GetStaticCondBlock(n);
724 bnd_mat = loc_mat->GetBlock(0, 0);
725
726 int nFacets = locExpansion->GetNtraces();
727
728 for (cnt2 = i = 0; i < nFacets; ++i)
729 {
730 int nCoeffs = elmtToTraceMap[i]->GetNcoeffs();
731 int elmtId = elmtToTraceMap[i]->GetElmtId();
732
733 // Ignore Dirichlet edges.
734 if (elmtId < nDir)
735 {
736 cnt += nCoeffs;
737 cnt2 += nCoeffs;
738 continue;
739 }
740
741 NekDouble *off = &tmpStore[arrayOffsets[elmtId]];
742
743 for (j = 0; j < nCoeffs; ++j)
744 {
745 NekDouble sign1 = asmMap->GetLocalToGlobalBndSign(cnt + j);
746 for (k = 0; k < nCoeffs; ++k)
747 {
748 NekDouble sign2 = asmMap->GetLocalToGlobalBndSign(cnt + k);
749 off[j * nCoeffs + k] +=
750 (*bnd_mat)(cnt2 + j, cnt2 + k) * sign1 * sign2;
751 }
752 }
753
754 cnt += nCoeffs;
755 cnt2 += nCoeffs;
756 }
757 }
758
759 // Set up IDs for universal numbering.
760 Array<OneD, long> uniIds(tmpStore.size());
761 for (cnt = 0, n = nDir; n < nTrace; ++n)
762 {
763 LocalRegions::ExpansionSharedPtr traceExp = trace->GetExp(n);
764 int nCoeffs = traceExp->GetNcoeffs();
765 int geomId = traceExp->GetGeom()->GetGlobalID();
766
767 for (i = 0; i < nCoeffs * nCoeffs; ++i)
768 {
769 uniIds[cnt++] = geomId * maxTraceSize + i + 1;
770 }
771 }
772
773 // Assemble matrices across partitions.
774 Gs::gs_data *gsh =
775 Gs::Init(uniIds, vComm,
776 expList->GetSession()->DefinesCmdLineArgument("verbose"));
777 Gs::Gather(tmpStore, Gs::gs_add, gsh);
778
779 // Set up diagonal block matrix
780 Array<OneD, unsigned int> n_blks(nTrace - nDir);
781 for (n = 0; n < nTrace - nDir; ++n)
782 {
783 n_blks[n] = trace->GetExp(n + nDir)->GetNcoeffs();
784 }
785
786 m_blkMat =
788
789 for (n = 0; n < nTrace - nDir; ++n)
790 {
791 int nCoeffs = trace->GetExp(n + nDir)->GetNcoeffs();
792 DNekMatSharedPtr tmp =
794 NekDouble *off = &tmpStore[arrayOffsets[n + nDir]];
795
796 for (i = 0; i < nCoeffs; ++i)
797 {
798 for (j = 0; j < nCoeffs; ++j)
799 {
800 (*tmp)(i, j) = off[i * nCoeffs + j];
801 }
802 }
803
804 tmp->Invert();
805 m_blkMat->SetBlock(n, n, tmp);
806 }
807}
std::shared_ptr< AssemblyMapDG > AssemblyMapDGSharedPtr
Definition: AssemblyMapDG.h:46

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::eDIAGONAL, Gs::Gather(), Gs::gs_add, Gs::Init(), m_blkMat, Nektar::MultiRegions::Preconditioner::m_linsys, Nektar::MultiRegions::Preconditioner::m_locToGloMap, and Nektar::LibUtilities::ReduceMax.

Referenced by v_BuildPreconditioner().

◆ create()

static PreconditionerSharedPtr Nektar::MultiRegions::PreconditionerBlock::create ( const std::shared_ptr< GlobalLinSys > &  plinsys,
const std::shared_ptr< AssemblyMap > &  pLocToGloMap 
)
inlinestatic

Creates an instance of this class.

Definition at line 52 of file PreconditionerBlock.h.

55 {
58 pLocToGloMap);
59 p->InitObject();
60 return p;
61 }
std::shared_ptr< Preconditioner > PreconditionerSharedPtr
Definition: GlobalLinSys.h:58

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), and CellMLToNektar.cellml_metadata::p.

◆ v_BuildPreconditioner()

void Nektar::MultiRegions::PreconditionerBlock::v_BuildPreconditioner ( )
overrideprotectedvirtual

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 79 of file PreconditionerBlock.cpp.

80{
81 GlobalLinSysKey key = m_linsys.lock()->GetKey();
82
83 // Different setup for HDG and CG.
84 if (key.GetMatrixType() == StdRegions::eHybridDGHelmBndLam)
85 {
87 }
88 else
89 {
91 }
92}
void BlockPreconditionerCG(void)
Construct a block preconditioner from for the continuous Galerkin system.
void BlockPreconditionerHDG(void)
Construct a block preconditioner for the hybridized discontinuous Galerkin system.

References BlockPreconditionerCG(), BlockPreconditionerHDG(), Nektar::StdRegions::eHybridDGHelmBndLam, Nektar::MultiRegions::GlobalMatrixKey::GetMatrixType(), and Nektar::MultiRegions::Preconditioner::m_linsys.

◆ v_DoPreconditioner()

void Nektar::MultiRegions::PreconditionerBlock::v_DoPreconditioner ( const Array< OneD, NekDouble > &  pInput,
Array< OneD, NekDouble > &  pOutput,
const bool &  isLocal = false 
)
overrideprotectedvirtual

Apply preconditioner to pInput and store the result in pOutput.

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 813 of file PreconditionerBlock.cpp.

816{
817 // Get assembly map and solver type
818 auto asmMap = m_locToGloMap.lock();
819
820 // Determine matrix size
821 int nGlobal = m_isFull ? asmMap->GetNumGlobalCoeffs()
822 : asmMap->GetNumGlobalBndCoeffs();
823 int nDir = asmMap->GetNumGlobalDirBndCoeffs();
824 int nNonDir = nGlobal - nDir;
825
826 // Apply preconditioner
827 DNekBlkMat &M = (*m_blkMat);
828
829 if (isLocal)
830 {
831 Array<OneD, NekDouble> wk(nGlobal);
832 Array<OneD, NekDouble> wk1(nGlobal);
833 asmMap->Assemble(pInput, wk);
834
835 NekVector<NekDouble> r(nNonDir, wk + nDir, eWrapper);
836 NekVector<NekDouble> z(nNonDir, wk1 + nDir, eWrapper);
837
838 z = M * r;
839 Vmath::Zero(nDir, wk1, 1);
840
841 asmMap->GlobalToLocal(wk1, pOutput);
842 }
843 else
844 {
845 NekVector<NekDouble> r(nNonDir, pInput, eWrapper);
846 NekVector<NekDouble> z(nNonDir, pOutput, eWrapper);
847 z = M * r;
848 }
849}
std::vector< double > z(NPUPPER)
NekMatrix< NekMatrix< NekDouble, StandardMatrixTag >, BlockMatrixTag > DNekBlkMat
Definition: NekTypeDefs.hpp:58
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.hpp:273

References Nektar::eWrapper, m_isFull, Nektar::MultiRegions::Preconditioner::m_locToGloMap, Nektar::UnitTests::z(), and Vmath::Zero().

◆ v_InitObject()

void Nektar::MultiRegions::PreconditionerBlock::v_InitObject ( )
overrideprotectedvirtual

Member Data Documentation

◆ className

string Nektar::MultiRegions::PreconditionerBlock::className
static
Initial value:
=
"Block", PreconditionerBlock::create, "Block Diagonal Preconditioning")
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
static PreconditionerSharedPtr create(const std::shared_ptr< GlobalLinSys > &plinsys, const std::shared_ptr< AssemblyMap > &pLocToGloMap)
Creates an instance of this class.
PreconFactory & GetPreconFactory()

Name of class.

Registers the class with the Factory.

Definition at line 64 of file PreconditionerBlock.h.

◆ m_blkMat

DNekBlkMatSharedPtr Nektar::MultiRegions::PreconditionerBlock::m_blkMat
protected

Definition at line 78 of file PreconditionerBlock.h.

Referenced by BlockPreconditionerCG(), and BlockPreconditionerHDG().

◆ m_gs_mapping

Gs::gs_data* Nektar::MultiRegions::PreconditionerBlock::m_gs_mapping
protected

Definition at line 81 of file PreconditionerBlock.h.

Referenced by BlockPreconditionerCG(), and ~PreconditionerBlock().

◆ m_isFull

bool Nektar::MultiRegions::PreconditionerBlock::m_isFull
protected

Definition at line 77 of file PreconditionerBlock.h.

Referenced by BlockPreconditionerCG(), v_DoPreconditioner(), and v_InitObject().