Nektar++
PreconditionerBlock.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: PreconditionerBlock.cpp
4//
5// For more information, please see: http://www.nektar.info
6//
7// The MIT License
8//
9// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10// Department of Aeronautics, Imperial College London (UK), and Scientific
11// Computing and Imaging Institute, University of Utah (USA).
12//
13// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description: Block Preconditioner definition
32//
33///////////////////////////////////////////////////////////////////////////////
34
36#include <LocalRegions/SegExp.h>
42#include <cmath>
43
44using namespace std;
45
47{
48/**
49 * Registers the class with the Factory.
50 */
53 "Block", PreconditionerBlock::create, "Block Diagonal Preconditioning");
54/**
55 * @class Block Preconditioner
56 *
57 * This class implements Block preconditioning for the conjugate
58 * gradient matrix solver.
59 */
60
62 const std::shared_ptr<GlobalLinSys> &plinsys,
63 const AssemblyMapSharedPtr &pLocToGloMap)
64 : Preconditioner(plinsys, pLocToGloMap)
65{
66}
67
69{
70 GlobalSysSolnType solverType = m_locToGloMap.lock()->GetGlobalSysSolnType();
71 m_isFull =
72 (solverType == eIterativeFull) || (solverType == ePETScFullMatrix);
73
76 "Solver type not valid");
77}
78
80{
81 GlobalLinSysKey key = m_linsys.lock()->GetKey();
82
83 // Different setup for HDG and CG.
85 {
87 }
88 else
89 {
91 }
92}
93
94/**
95 * \brief Construct a block preconditioner from \f$\mathbf{S}_{1}\f$ for
96 * the continuous Galerkin system.
97 *
98 * The preconditioner for the statically-condensed system is defined as:
99 *
100 * \f[\mathbf{M}^{-1}=\left[\begin{array}{ccc}
101 * \mathrm{Diag}[(\mathbf{S_{1}})_{vv}] & & \\
102 * & (\mathbf{S}_{1})_{eb} & \\
103 * & & (\mathbf{S}_{1})_{fb} \end{array}\right] \f]
104 *
105 * where \f$\mathbf{S}_{1}\f$ is the local Schur complement matrix for
106 * each element and the subscript denotes the portion of the Schur
107 * complement associated with the vertex (vv), edge (eb) and face
108 * blocks (fb) respectively.
109 *
110 * In case of the full system \f$\mathbf{A}\f$, the preconditioner
111 * includes the interior blocks as well. The full block preconditioner
112 * is defined as:
113 *
114 * \f[\mathbf{M}^{-1}=\left[\begin{array}{cccc}
115 * \mathrm{Diag}[(\mathbf{A})_{vv}] & & & \\
116 * & (\mathbf{A})_{eb} & & \\
117 * & & (\mathbf{A})_{fb} & \\
118 * & & & (\mathbf{A})_{int} \end{array}\right] \f]
119 *
120 * where the interior portion is added at the end following the
121 * ordering inside data structures.
122 */
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;
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 {
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;
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}
642
643/**
644 * @brief Construct a block preconditioner for the hybridized
645 * discontinuous Galerkin system.
646 *
647 * This system is built in a similar fashion to its continuous variant
648 * found in PreconditionerBlock::BlockPreconditionerCG. In this setting
649 * however, the matrix is constructed as:
650 *
651 * \f[ M^{-1} = \mathrm{Diag}[ (\mathbf{S_{1}})_{f}^{-1} ] \f]
652 *
653 * where each matrix is the Schur complement system restricted to a
654 * single face of the trace system.
655 */
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
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}
808
809/**
810 * @brief Apply preconditioner to \p pInput and store the result in \p
811 * pOutput.
812 */
814 const Array<OneD, NekDouble> &pInput, Array<OneD, NekDouble> &pOutput,
815 const bool &isLocal)
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}
850} // namespace Nektar::MultiRegions
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
#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
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
StdRegions::MatrixType GetMatrixType() const
Return the matrix type.
void BlockPreconditionerCG(void)
Construct a block preconditioner from for the continuous Galerkin system.
static PreconditionerSharedPtr create(const std::shared_ptr< GlobalLinSys > &plinsys, const std::shared_ptr< AssemblyMap > &pLocToGloMap)
Creates an instance of this class.
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.
void BlockPreconditionerHDG(void)
Construct a block preconditioner for the hybridized discontinuous Galerkin system.
PreconditionerBlock(const std::shared_ptr< GlobalLinSys > &plinsys, const AssemblyMapSharedPtr &pLocToGloMap)
static std::string className
Name of class.
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< AssemblyMapDG > AssemblyMapDGSharedPtr
Definition: AssemblyMapDG.h:46
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
PreconFactory & GetPreconFactory()
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::shared_ptr< AssemblyMap > AssemblyMapSharedPtr
Definition: AssemblyMap.h:50
std::map< int, std::vector< PeriodicEntity > > PeriodicMap
StdRegions::Orientation DeterminePeriodicFaceOrient(StdRegions::Orientation faceOrient, StdRegions::Orientation perFaceOrient)
Determine relative orientation between two faces.
std::vector< double > z(NPUPPER)
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
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.hpp:273
STL namespace.