Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | Static Public Member Functions | Static Public Attributes | Protected Attributes | Private Member Functions | List of all members
Nektar::MultiRegions::PreconditionerBlock Class Reference

#include <PreconditionerBlock.h>

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

Public Member Functions

 PreconditionerBlock (const boost::shared_ptr< GlobalLinSys > &plinsys, const AssemblyMapSharedPtr &pLocToGloMap)
 
virtual ~PreconditionerBlock ()
 
- Public Member Functions inherited from Nektar::MultiRegions::Preconditioner
 Preconditioner (const boost::shared_ptr< GlobalLinSys > &plinsys, const AssemblyMapSharedPtr &pLocToGloMap)
 
virtual ~Preconditioner ()
 
void DoPreconditioner (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
 
void DoPreconditionerWithNonVertOutput (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const Array< OneD, NekDouble > &pNonVertOutput, Array< OneD, NekDouble > &pVertForce=NullNekDouble1DArray)
 
void DoTransformToLowEnergy (Array< OneD, NekDouble > &pInOut, int offset)
 
void DoTransformToLowEnergy (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
 
void DoTransformFromLowEnergy (Array< OneD, NekDouble > &pInOut)
 
void DoMultiplybyInverseTransformationMatrix (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
 
void DoMultiplybyInverseTransposedTransformationMatrix (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, const boost::shared_ptr< DNekScalMat > &loc_mat)
 

Static Public Member Functions

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

Static Public Attributes

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

Protected Attributes

const boost::weak_ptr
< GlobalLinSys
m_linsys
 
PreconditionerType m_preconType
 
DNekBlkMatSharedPtr m_blkMat
 
boost::shared_ptr< AssemblyMapm_locToGloMap
 
- Protected Attributes inherited from Nektar::MultiRegions::Preconditioner
const boost::weak_ptr
< GlobalLinSys
m_linsys
 
PreconditionerType m_preconType
 
DNekMatSharedPtr m_preconditioner
 
boost::shared_ptr< AssemblyMapm_locToGloMap
 
LibUtilities::CommSharedPtr m_comm
 

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...
 
virtual void v_InitObject ()
 
virtual void v_DoPreconditioner (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
 Apply preconditioner to pInput and store the result in pOutput. More...
 
virtual void v_BuildPreconditioner ()
 

Additional Inherited Members

- Protected Member Functions inherited from Nektar::MultiRegions::Preconditioner
virtual DNekScalMatSharedPtr v_TransformedSchurCompl (int offset, const boost::shared_ptr< DNekScalMat > &loc_mat)
 Get block elemental transposed transformation matrix $\mathbf{R}^{T}$. More...
 

Detailed Description

Definition at line 52 of file PreconditionerBlock.h.

Constructor & Destructor Documentation

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

Definition at line 67 of file PreconditionerBlock.cpp.

70  : Preconditioner(plinsys, pLocToGloMap),
71  m_linsys(plinsys),
72  m_preconType(pLocToGloMap->GetPreconType()),
73  m_locToGloMap(pLocToGloMap)
74  {
75  }
Preconditioner(const boost::shared_ptr< GlobalLinSys > &plinsys, const AssemblyMapSharedPtr &pLocToGloMap)
boost::shared_ptr< AssemblyMap > m_locToGloMap
const boost::weak_ptr< GlobalLinSys > m_linsys
virtual Nektar::MultiRegions::PreconditionerBlock::~PreconditionerBlock ( )
inlinevirtual

Definition at line 75 of file PreconditionerBlock.h.

75 {}

Member Function Documentation

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

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

The preconditioner 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, edge and face blocks respectively.

Definition at line 116 of file PreconditionerBlock.cpp.

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

Referenced by v_BuildPreconditioner().

117  {
118  ExpListSharedPtr expList = m_linsys.lock()->GetLocMat().lock();
120  DNekScalBlkMatSharedPtr loc_mat;
121  DNekScalMatSharedPtr bnd_mat;
122 
123  int nel, i, j, k, n, cnt, gId;
124  int meshVertId, meshEdgeId, meshFaceId;
125 
126  const int nExp = expList->GetExpSize();
127  const int nDirBnd = m_locToGloMap->GetNumGlobalDirBndCoeffs();
128 
129  // Grab periodic geometry information.
131  PeriodicMap periodicVerts, periodicEdges, periodicFaces;
132  expList->GetPeriodicEntities(
133  periodicVerts, periodicEdges, periodicFaces);
134 
135  // The vectors below are of size 3 to have separate storage for
136  // vertices, edges and faces.
137 
138  // Maps from geometry ID to the matrix representing the extracted
139  // portion of S_1. For example idMats[2] folds the S_1 face blocks.
140  vector<map<int, vector<NekDouble> > > idMats(3);
141 
142  // Maps from the global ID, as obtained from AssemblyMapCG's
143  // localToGlobalMap, to the geometry ID.
144  vector<map<int, int> > gidMeshIds(3);
145 
146  // Maps from the global ID to the number of degrees of freedom for
147  // this geometry object.
148  vector<map<int, int> > gidDofs(3);
149 
150  // Array containing maximum information needed for the universal
151  // numbering later. For i = 0,1,2 for each geometry dimension:
152  //
153  // maxVertIds[2*i] = maximum geometry ID at dimension i
154  // maxVertIds[2*i+1] = maximum number of degrees of freedom for all
155  // elements of dimension i.
156  Array<OneD, int> maxVertIds(6, -1);
157 
158  // Iterator for idMats members.
159  map<int, vector<NekDouble> >::iterator gIt;
160 
161  // Figure out mapping from each elemental contribution to offset in
162  // (vert,edge,face) triples.
163  for (cnt = n = 0; n < nExp; ++n)
164  {
165  nel = expList->GetOffset_Elmt_Id(n);
166  exp = expList->GetExp(nel);
167 
168  // Grab reference to local Schur complement matrix.
169  DNekScalMatSharedPtr schurMat =
170  m_linsys.lock()->GetStaticCondBlock(n)->GetBlock(0,0);
171 
172  // Process vertices to extract relevant portion of the Schur
173  // complement matrix.
174  for (i = 0; i < exp->GetNverts(); ++i)
175  {
176  meshVertId = exp->GetGeom()->GetVid(i);
177  int locId = exp->GetVertexMap(i);
178 
179  // Get the global ID of this vertex.
180  gId = m_locToGloMap->GetLocalToGlobalMap(
181  cnt + locId) - nDirBnd;
182 
183  // Ignore all Dirichlet vertices.
184  if (gId < 0)
185  {
186  continue;
187  }
188 
189  gidDofs[0][gId] = 1;
190 
191  // Extract vertex value from Schur complement matrix.
192  NekDouble vertVal = (*schurMat)(locId,locId);
193 
194  // See if we have processed this vertex from another
195  // element.
196  gIt = idMats[0].find(gId);
197 
198  if (gIt == idMats[0].end())
199  {
200  // If not then put our 'matrix' inside idMats.
201  idMats[0][gId] = vector<NekDouble>(1, vertVal);
202  }
203  else
204  {
205  // Otherwise combine with the value that is already
206  // there (i.e. do assembly on this degree of freedom).
207  gIt->second[0] += vertVal;
208  }
209 
210  // Now check to see if the vertex is periodic. If it is,
211  // then we change meshVertId to be the minimum of all the
212  // other periodic vertices, so that we don't end up
213  // duplicating the matrix in our final block matrix.
214  pIt = periodicVerts.find(meshVertId);
215  if (pIt != periodicVerts.end())
216  {
217  for (j = 0; j < pIt->second.size(); ++j)
218  {
219  meshVertId = min(meshVertId, pIt->second[j].id);
220  }
221  }
222 
223  // Finally record the other information we need into the
224  // other matrices.
225  gidMeshIds[0][gId] = meshVertId;
226  maxVertIds[0] = max(maxVertIds[0], meshVertId);
227  maxVertIds[1] = 1;
228  }
229 
230  // Process edges. This logic is mostly the same as the previous
231  // block.
232  for (i = 0; i < exp->GetNedges(); ++i)
233  {
234  meshEdgeId = exp->GetGeom()->GetEid(i);
235 
236  Array<OneD, unsigned int> bmap, bmap2;
237  Array<OneD, int> sign;
238  StdRegions::Orientation edgeOrient =
239  exp->GetGeom()->GetEorient(i);
240 
241  // Check if this edge is periodic. We may need to flip
242  // orientation if it is.
243  pIt = periodicEdges.find(meshEdgeId);
244  if (pIt != periodicEdges.end())
245  {
246  pair<int, StdRegions::Orientation> idOrient =
248  meshEdgeId, edgeOrient, pIt->second);
249  meshEdgeId = idOrient.first;
250  edgeOrient = idOrient.second;
251  }
252 
253  // Grab edge interior map, and the edge inverse boundary
254  // map, so that we can extract this edge from the Schur
255  // complement matrix.
256  exp->GetEdgeInteriorMap(i, edgeOrient, bmap, sign);
257  bmap2 = exp->GetEdgeInverseBoundaryMap(i);
258 
259  // Allocate temporary storage for the extracted edge matrix.
260  const int nEdgeCoeffs = bmap.num_elements();
261  vector<NekDouble> tmpStore(nEdgeCoeffs*nEdgeCoeffs);
262 
263  gId = m_locToGloMap->GetLocalToGlobalMap(cnt + bmap[0]);
264 
265  for (j = 0; j < nEdgeCoeffs; ++j)
266  {
267  // We record the minimum ID from the edge for our
268  // maps. This follows the logic that the assembly map
269  // ordering will always give us a contiguous ordering of
270  // global degrees of freedom for edge interior
271  // coefficients.
272  gId = min(gId,
273  m_locToGloMap->GetLocalToGlobalMap(
274  cnt + bmap[j])
275  - nDirBnd);
276 
277  // Ignore Dirichlet edges.
278  if (gId < 0)
279  {
280  continue;
281  }
282 
283  const NekDouble sign1 = sign[j];
284 
285  // Extract this edge, along with sign array for assembly
286  // later.
287  for (k = 0; k < nEdgeCoeffs; ++k)
288  {
289  tmpStore[k+j*nEdgeCoeffs] =
290  sign1*sign[k]*(*schurMat)(bmap2[j], bmap2[k]);
291  }
292  }
293 
294  if (gId < 0)
295  {
296  continue;
297  }
298 
299  gidDofs[1][gId] = nEdgeCoeffs;
300 
301  // Assemble this edge matrix with another one, if it exists.
302  gIt = idMats[1].find(gId);
303  if (gIt == idMats[1].end())
304  {
305  idMats[1][gId] = tmpStore;
306  }
307  else
308  {
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);
313  }
314 
315  gidMeshIds[1][gId] = meshEdgeId;
316  maxVertIds[2] = max(maxVertIds[2], meshEdgeId);
317  maxVertIds[3] = max(maxVertIds[3], nEdgeCoeffs);
318  }
319 
320  // Process faces. This logic is mostly the same as the previous
321  // block.
322  for (i = 0; i < exp->GetNfaces(); ++i)
323  {
324  meshFaceId = exp->GetGeom()->GetFid(i);
325 
326  Array<OneD, unsigned int> bmap, bmap2;
327  Array<OneD, int> sign;
328  StdRegions::Orientation faceOrient =
329  exp->GetGeom()->GetForient(i);
330 
331  // Check if this face is periodic. We may need to flip
332  // orientation if it is.
333  pIt = periodicFaces.find(meshFaceId);
334  if (pIt != periodicFaces.end())
335  {
336  meshFaceId = min(meshFaceId, pIt->second[0].id);
337  faceOrient = DeterminePeriodicFaceOrient(
338  faceOrient, pIt->second[0].orient);
339  }
340 
341  exp->GetFaceInteriorMap(i, faceOrient, bmap, sign);
342  bmap2 = exp->GetFaceInverseBoundaryMap(i);
343 
344  // Allocate temporary storage for the extracted face matrix.
345  const int nFaceCoeffs = bmap.num_elements();
346  vector<NekDouble> tmpStore(nFaceCoeffs*nFaceCoeffs);
347 
348  gId = m_locToGloMap->GetLocalToGlobalMap(cnt + bmap[0]);
349 
350  for (j = 0; j < nFaceCoeffs; ++j)
351  {
352  gId = min(gId,
353  m_locToGloMap->GetLocalToGlobalMap(
354  cnt + bmap[j])
355  - nDirBnd);
356 
357  // Ignore Dirichlet faces.
358  if (gId < 0)
359  {
360  continue;
361  }
362 
363  const NekDouble sign1 = sign[j];
364 
365  // Extract this face, along with sign array for assembly
366  // later.
367  for (k = 0; k < nFaceCoeffs; ++k)
368  {
369  tmpStore[k+j*nFaceCoeffs] =
370  sign1*sign[k]*(*schurMat)(bmap2[j], bmap2[k]);
371  }
372  }
373 
374  if (gId < 0)
375  {
376  continue;
377  }
378 
379  gidDofs[2][gId] = nFaceCoeffs;
380 
381  // Assemble this face matrix with another one, if it exists.
382  gIt = idMats[2].find(gId);
383  if (gIt == idMats[2].end())
384  {
385  idMats[2][gId] = tmpStore;
386  }
387  else
388  {
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);
393  }
394 
395  gidMeshIds[2][gId] = meshFaceId;
396  maxVertIds[4] = max(maxVertIds[4], meshFaceId);
397  maxVertIds[5] = max(maxVertIds[5], nFaceCoeffs);
398  }
399 
400  cnt += exp->GetNcoeffs();
401  }
402 
403  // Perform a reduction to find maximum vertex, edge and face
404  // geometry IDs.
405  m_comm = expList->GetSession()->GetComm()->GetRowComm();
406  m_comm->AllReduce(maxVertIds, LibUtilities::ReduceMax);
407 
408  // Concatenate all matrices into contiguous storage and figure out
409  // universal ID numbering.
410  vector<NekDouble> storageBuf;
411  vector<long> globalToUniversal;
412 
413  for (i = 0, cnt = 1; i < 3; ++i)
414  {
415  const int maxDofs = maxVertIds[2*i+1];
416 
417  // Note that iterating over the map uses the property that keys
418  // are accessed in order of ascending order, putting everything
419  // in the right order for the global system.
420  for (gIt = idMats[i].begin(); gIt != idMats[i].end(); ++gIt)
421  {
422  // Copy matrix into storage.
423  storageBuf.insert(storageBuf.end(),
424  gIt->second.begin(), gIt->second.end());
425 
426  // Get mesh ID from global ID number.
427  ASSERTL1(gidMeshIds[i].count(gIt->first) > 0,
428  "Unable to find global ID " +
429  boost::lexical_cast<string>(gIt->first) +
430  " inside map");
431  meshVertId = gidMeshIds[i][gIt->first];
432 
433  for (j = 0; j < gIt->second.size(); ++j)
434  {
435  globalToUniversal.push_back(
436  cnt + meshVertId*maxDofs*maxDofs + j);
437  }
438 
439  // Free up the temporary storage.
440  gIt->second.clear();
441  }
442 
443  cnt += (maxVertIds[2*i]+1)*maxDofs*maxDofs;
444  }
445 
446  ASSERTL1(storageBuf.size() == globalToUniversal.size(),
447  "Storage buffer and global to universal map size does "
448  "not match");
449 
450  Array<OneD, NekDouble> storageData(
451  storageBuf.size(), &storageBuf[0]);
452  Array<OneD, long> globalToUniversalMap(
453  globalToUniversal.size(), &globalToUniversal[0]);
454 
455  // Use GS to assemble data between processors.
456  Gs::gs_data *tmpGs = Gs::Init(globalToUniversalMap, m_comm);
457  Gs::Gather(storageData, Gs::gs_add, tmpGs);
458 
459  // Figure out what storage we need in the block matrix.
460  Array<OneD, unsigned int> n_blks(
461  1 + idMats[1].size() + idMats[2].size());
462 
463  // Vertex block is a diagonal matrix.
464  n_blks[0] = idMats[0].size();
465 
466  // Now extract number of rows in each edge and face block from the
467  // gidDofs map.
468  cnt = 1;
469  for (i = 1; i < 3; ++i)
470  {
471  for (gIt = idMats[i].begin(); gIt != idMats[i].end(); ++gIt)
472  {
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>(
476  gIt->first));
477 
478  n_blks[cnt++] = gidDofs[i][gIt->first];
479  }
480  }
481 
482  // Allocate storage for the block matrix.
484  ::AllocateSharedPtr(n_blks, n_blks, eDIAGONAL);
485 
486  // We deal with the vertex matrix separately since all vertices can
487  // be condensed into a single, block-diagonal matrix.
489  ::AllocateSharedPtr(n_blks[0], n_blks[0], 0.0, eDIAGONAL);
490 
491  // Fill the vertex matrix with the inverse of each vertex value.
492  cnt = 0;
493  for (gIt = idMats[0].begin(); gIt != idMats[0].end(); ++gIt, ++cnt)
494  {
495  (*vertMat)(cnt, cnt) = 1.0/storageData[cnt];
496  }
497 
498  // Put the vertex matrix in the block matrix.
499  m_blkMat->SetBlock(0,0,vertMat);
500 
501  // Finally, grab the matrices from the block storage, invert them
502  // and place them in the correct position inside the block matrix.
503  int cnt2 = 1;
504  for (i = 1; i < 3; ++i)
505  {
506  for (gIt = idMats[i].begin(); gIt != idMats[i].end();
507  ++gIt, ++cnt2)
508  {
509  int nDofs = gidDofs[i][gIt->first];
510 
512  ::AllocateSharedPtr(nDofs, nDofs);
513 
514  for (j = 0; j < nDofs; ++j)
515  {
516  for (k = 0; k < nDofs; ++k)
517  {
518  (*tmp)(j,k) = storageData[k+j*nDofs + cnt];
519  }
520  }
521 
522  cnt += nDofs*nDofs;
523 
524  tmp->Invert();
525  m_blkMat->SetBlock(cnt2, cnt2, tmp);
526  }
527  }
528  }
LibUtilities::CommSharedPtr m_comm
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:224
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:22
static gs_data * Init(const Nektar::Array< OneD, long > pId, const LibUtilities::CommSharedPtr &pComm)
Initialise Gather-Scatter map.
Definition: GsLib.hpp:151
boost::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:70
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
Definition: NekTypeDefs.hpp:74
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
double NekDouble
boost::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
std::map< int, std::vector< PeriodicEntity > > PeriodicMap
boost::shared_ptr< AssemblyMap > m_locToGloMap
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
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...
Definition: ErrorUtil.hpp:191
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.cpp:285
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 543 of file PreconditionerBlock.cpp.

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

Referenced by v_BuildPreconditioner().

544  {
545  boost::shared_ptr<MultiRegions::ExpList>
546  expList = ((m_linsys.lock())->GetLocMat()).lock();
547  boost::shared_ptr<MultiRegions::ExpList> trace = expList->GetTrace();
549  DNekScalBlkMatSharedPtr loc_mat;
550  DNekScalMatSharedPtr bnd_mat;
551 
552  AssemblyMapDGSharedPtr asmMap = boost::dynamic_pointer_cast<
553  AssemblyMapDG>(m_locToGloMap);
554 
555  int i, j, k, n, cnt, cnt2;
556 
557  // Figure out number of Dirichlet trace elements
558  int nTrace = expList->GetTrace()->GetExpSize();
559  int nDir = m_locToGloMap->GetNumGlobalDirBndCoeffs();
560 
561  for (cnt = n = 0; n < nTrace; ++n)
562  {
563  if (cnt >= nDir)
564  {
565  break;
566  }
567 
568  cnt += trace->GetExp(n)->GetNcoeffs();
569  }
570 
571  nDir = n;
572 
573  // Allocate storage for block matrix. Need number of unique faces in
574  // trace space.
575  int maxTraceSize = -1;
576  map<int, int> arrayOffsets;
577 
578  for (cnt = 0, n = nDir; n < nTrace; ++n)
579  {
580  int nCoeffs = trace->GetExp(n)->GetNcoeffs();
581  int nCoeffs2 = nCoeffs * nCoeffs;
582 
583  arrayOffsets[n] = cnt;
584  cnt += nCoeffs2;
585 
586  if (maxTraceSize < nCoeffs2)
587  {
588  maxTraceSize = nCoeffs2;
589  }
590  }
591 
592  // Find maximum trace size.
593  m_comm = expList->GetSession()->GetComm()->GetRowComm();
594  m_comm->AllReduce(maxTraceSize, LibUtilities::ReduceMax);
595 
596  // Zero matrix storage.
597  Array<OneD, NekDouble> tmpStore(cnt, 0.0);
598 
599  // Assemble block matrices for each trace element.
600  for (cnt = n = 0; n < expList->GetExpSize(); ++n)
601  {
602  int elmt = expList->GetOffset_Elmt_Id(n);
603  locExpansion = expList->GetExp(elmt);
604 
605  Array<OneD, LocalRegions::ExpansionSharedPtr> &elmtToTraceMap =
606  asmMap->GetElmtToTrace()[elmt];
607 
608  // Block matrix (lambda)
609  loc_mat = (m_linsys.lock())->GetStaticCondBlock(n);
610  bnd_mat = loc_mat->GetBlock(0,0);
611 
612  int nFacets = locExpansion->GetNumBases() == 2 ?
613  locExpansion->GetNedges() : locExpansion->GetNfaces();
614 
615  for (cnt2 = i = 0; i < nFacets; ++i)
616  {
617  int nCoeffs = elmtToTraceMap[i]->GetNcoeffs();
618  int elmtId = elmtToTraceMap[i]->GetElmtId ();
619 
620  // Ignore Dirichlet edges.
621  if (elmtId < nDir)
622  {
623  cnt += nCoeffs;
624  cnt2 += nCoeffs;
625  continue;
626  }
627 
628  NekDouble *off = &tmpStore[arrayOffsets[elmtId]];
629 
630  for (j = 0; j < nCoeffs; ++j)
631  {
632  NekDouble sign1 = asmMap->GetLocalToGlobalBndSign(
633  cnt + j);
634  for (k = 0; k < nCoeffs; ++k)
635  {
636  NekDouble sign2 = asmMap->GetLocalToGlobalBndSign(
637  cnt + k);
638  off[j*nCoeffs + k] +=
639  (*bnd_mat)(cnt2+j, cnt2+k) * sign1 * sign2;
640  }
641  }
642 
643  cnt += nCoeffs;
644  cnt2 += nCoeffs;
645  }
646  }
647 
648  // Set up IDs for universal numbering.
649  Array<OneD, long> uniIds(tmpStore.num_elements());
650  for (cnt = 0, n = nDir; n < nTrace; ++n)
651  {
652  LocalRegions::ExpansionSharedPtr traceExp = trace->GetExp(n);
653  int nCoeffs = traceExp->GetNcoeffs();
654  int geomId = traceExp->GetGeom()->GetGlobalID();
655 
656  for (i = 0; i < nCoeffs*nCoeffs; ++i)
657  {
658  uniIds[cnt++] = geomId * maxTraceSize + i + 1;
659  }
660  }
661 
662  // Assemble matrices across partitions.
663  Gs::gs_data *gsh = Gs::Init(uniIds, m_comm);
664  Gs::Gather(tmpStore, Gs::gs_add, gsh);
665 
666  // Set up diagonal block matrix
667  Array<OneD, unsigned int> n_blks(nTrace - nDir);
668  for (n = 0; n < nTrace - nDir; ++n)
669  {
670  n_blks[n] = trace->GetExp(n + nDir)->GetNcoeffs();
671  }
672 
674  ::AllocateSharedPtr(n_blks, n_blks, eDIAGONAL);
675 
676  for (n = 0; n < nTrace - nDir; ++n)
677  {
678  int nCoeffs = trace->GetExp(n + nDir)->GetNcoeffs();
680  ::AllocateSharedPtr(nCoeffs, nCoeffs);
681  NekDouble *off = &tmpStore[arrayOffsets[n+nDir]];
682 
683  for (i = 0; i < nCoeffs; ++i)
684  {
685  for (j = 0; j < nCoeffs; ++j)
686  {
687  (*tmp)(i,j) = off[i*nCoeffs + j];
688  }
689  }
690 
691  tmp->Invert();
692  m_blkMat->SetBlock(n, n, tmp);
693  }
694  }
boost::shared_ptr< AssemblyMapDG > AssemblyMapDGSharedPtr
Definition: AssemblyMapDG.h:49
LibUtilities::CommSharedPtr m_comm
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:224
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
static gs_data * Init(const Nektar::Array< OneD, long > pId, const LibUtilities::CommSharedPtr &pComm)
Initialise Gather-Scatter map.
Definition: GsLib.hpp:151
boost::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:70
boost::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
boost::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:74
double NekDouble
boost::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
boost::shared_ptr< AssemblyMap > m_locToGloMap
const boost::weak_ptr< GlobalLinSys > m_linsys
static PreconditionerSharedPtr Nektar::MultiRegions::PreconditionerBlock::create ( const boost::shared_ptr< GlobalLinSys > &  plinsys,
const boost::shared_ptr< AssemblyMap > &  pLocToGloMap 
)
inlinestatic

Creates an instance of this class.

Definition at line 56 of file PreconditionerBlock.h.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr().

60  {
62  ::AllocateSharedPtr(plinsys,pLocToGloMap);
63  p->InitObject();
64  return p;
65  }
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
boost::shared_ptr< Preconditioner > PreconditionerSharedPtr
Definition: GlobalLinSys.h:62
void Nektar::MultiRegions::PreconditionerBlock::v_BuildPreconditioner ( )
privatevirtual

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 85 of file PreconditionerBlock.cpp.

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

86  {
87  GlobalLinSysKey key = m_linsys.lock()->GetKey();
88 
89  // Different setup for HDG and CG.
90  if (key.GetMatrixType() == StdRegions::eHybridDGHelmBndLam)
91  {
93  }
94  else
95  {
97  }
98  }
void BlockPreconditionerHDG(void)
Construct a block preconditioner for the hybridized discontinuous Galerkin system.
void BlockPreconditionerCG(void)
Construct a block preconditioner from for the continuous Galerkin system.
const boost::weak_ptr< GlobalLinSys > m_linsys
void Nektar::MultiRegions::PreconditionerBlock::v_DoPreconditioner ( const Array< OneD, NekDouble > &  pInput,
Array< OneD, NekDouble > &  pOutput 
)
privatevirtual

Apply preconditioner to pInput and store the result in pOutput.

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 700 of file PreconditionerBlock.cpp.

References Nektar::eWrapper, and m_locToGloMap.

703  {
704  int nDir = m_locToGloMap->GetNumGlobalDirBndCoeffs();
705  int nGlobal = m_locToGloMap->GetNumGlobalBndCoeffs();
706  int nNonDir = nGlobal-nDir;
707  DNekBlkMat &M = (*m_blkMat);
708  NekVector<NekDouble> r(nNonDir,pInput,eWrapper);
709  NekVector<NekDouble> z(nNonDir,pOutput,eWrapper);
710  z = M * r;
711  }
NekMatrix< NekMatrix< NekDouble, StandardMatrixTag >, BlockMatrixTag > DNekBlkMat
Definition: NekTypeDefs.hpp:60
boost::shared_ptr< AssemblyMap > m_locToGloMap
void Nektar::MultiRegions::PreconditionerBlock::v_InitObject ( )
privatevirtual

Reimplemented from Nektar::MultiRegions::Preconditioner.

Definition at line 77 of file PreconditionerBlock.cpp.

References ASSERTL0, Nektar::MultiRegions::eIterativeStaticCond, Nektar::MultiRegions::ePETScStaticCond, and m_locToGloMap.

78  {
79  GlobalSysSolnType solvertype=m_locToGloMap->GetGlobalSysSolnType();
81  solvertype == MultiRegions::ePETScStaticCond,
82  "Solver type not valid");
83  }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
boost::shared_ptr< AssemblyMap > m_locToGloMap

Member Data Documentation

string Nektar::MultiRegions::PreconditionerBlock::className
static
Initial value:
"Block",
"Block Diagonal Preconditioning")

Name of class.

Registers the class with the Factory.

Definition at line 68 of file PreconditionerBlock.h.

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

Definition at line 80 of file PreconditionerBlock.h.

Referenced by BlockPreconditionerCG(), and BlockPreconditionerHDG().

const boost::weak_ptr<GlobalLinSys> Nektar::MultiRegions::PreconditionerBlock::m_linsys
protected
boost::shared_ptr<AssemblyMap> Nektar::MultiRegions::PreconditionerBlock::m_locToGloMap
protected
PreconditionerType Nektar::MultiRegions::PreconditionerBlock::m_preconType
protected

Definition at line 79 of file PreconditionerBlock.h.