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

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

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 703 of file PreconditionerBlock.cpp.

References Nektar::eWrapper, and m_locToGloMap.

706  {
707  int nDir = m_locToGloMap->GetNumGlobalDirBndCoeffs();
708  int nGlobal = m_locToGloMap->GetNumGlobalBndCoeffs();
709  int nNonDir = nGlobal-nDir;
710  DNekBlkMat &M = (*m_blkMat);
711  NekVector<NekDouble> r(nNonDir,pInput,eWrapper);
712  NekVector<NekDouble> z(nNonDir,pOutput,eWrapper);
713  z = M * r;
714  }
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:198
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.