Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Block Preconditioner definition
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <LibUtilities/BasicUtils/VDmathArray.hpp>
37 #include <MultiRegions/PreconditionerBlock.h>
38 #include <MultiRegions/GlobalMatrixKey.h>
39 #include <MultiRegions/GlobalLinSysIterativeStaticCond.h>
40 #include <MultiRegions/GlobalLinSys.h>
41 #include <MultiRegions/AssemblyMap/AssemblyMapDG.h>
42 #include <LocalRegions/MatrixKey.h>
43 #include <LocalRegions/SegExp.h>
44 #include <math.h>
45 
46 namespace Nektar
47 {
48  namespace MultiRegions
49  {
50  /**
51  * Registers the class with the Factory.
52  */
53  string PreconditionerBlock::className
54  = GetPreconFactory().RegisterCreatorFunction(
55  "Block",
56  PreconditionerBlock::create,
57  "Block Diagonal Preconditioning");
58  /**
59  * @class Block Preconditioner
60  *
61  * This class implements Block preconditioning for the conjugate
62  * gradient matrix solver.
63  */
64 
65  PreconditionerBlock::PreconditionerBlock(
66  const boost::shared_ptr<GlobalLinSys> &plinsys,
67  const AssemblyMapSharedPtr &pLocToGloMap)
68  : Preconditioner(plinsys, pLocToGloMap),
69  m_linsys(plinsys),
70  m_preconType(pLocToGloMap->GetPreconType()),
71  m_locToGloMap(pLocToGloMap)
72  {
73  }
74 
75  void PreconditionerBlock::v_InitObject()
76  {
77  GlobalSysSolnType solvertype=m_locToGloMap->GetGlobalSysSolnType();
78  ASSERTL0(solvertype == MultiRegions::eIterativeStaticCond,"Solver type not valid");
79  }
80 
81 
82  void PreconditionerBlock::v_BuildPreconditioner()
83  {
84  // Different setup for HDG
85  GlobalLinSysKey key = m_linsys.lock()->GetKey();
86  if (key.GetMatrixType() == StdRegions::eHybridDGHelmBndLam)
87  {
88  BlockPreconditionerHDG();
89  return;
90  }
91 
92  boost::shared_ptr<MultiRegions::ExpList>
93  expList=((m_linsys.lock())->GetLocMat()).lock();
94  StdRegions::StdExpansionSharedPtr locExpansion;
95  locExpansion = expList->GetExp(0);
96  int nDim = locExpansion->GetShapeDimension();
97 
98  if (nDim == 1)
99  {
100  ASSERTL0(0,"Unknown preconditioner");
101  }
102  else if (nDim == 2)
103  {
104  BlockPreconditioner2D();
105  }
106  else if (nDim == 3)
107  {
108  BlockPreconditioner3D();
109  }
110  }
111 
112  /**
113  * \brief Construct a block preconditioner from
114  * \f$\mathbf{S}_{1}\f$
115  *
116  *\f[\mathbf{M}^{-1}=\left[\begin{array}{ccc}
117  * Diag[(\mathbf{S_{1}})_{vv}] & & \\ & (\mathbf{S}_{1})_{eb} & \\ & &
118  * (\mathbf{S}_{1})_{fb} \end{array}\right] \f]
119  *
120  * where \f$\mathbf{S}_{1}\f$ is the local schur complement matrix for
121  * each element.
122  */
123  void PreconditionerBlock::BlockPreconditioner2D()
124  {
125  boost::shared_ptr<MultiRegions::ExpList>
126  expList=((m_linsys.lock())->GetLocMat()).lock();
127  LocalRegions::ExpansionSharedPtr locExpansion;
128  GlobalLinSysKey m_linSysKey=(m_linsys.lock())->GetKey();
129  StdRegions::VarCoeffMap vVarCoeffMap;
130  int i, j, k, nel;
131  int nVerts, nEdges;
132  int eid, n, cnt, nedgemodes;
133  NekDouble zero = 0.0;
134 
135  int vMap1, vMap2, sign1, sign2;
136  int m, v, eMap1, eMap2;
137  int offset, globalrow, globalcol;
138 
139  //matrix storage
140  MatrixStorage storage = eFULL;
141  MatrixStorage vertstorage = eDIAGONAL;
142  MatrixStorage blkmatStorage = eDIAGONAL;
143 
144  //local element static condensed matrices
145  DNekScalBlkMatSharedPtr loc_mat;
146  DNekScalMatSharedPtr bnd_mat;
147 
148  DNekMatSharedPtr VertBlk;
149 
150  int nDirBnd = m_locToGloMap->GetNumGlobalDirBndCoeffs();
151  int nNonDirVerts = m_locToGloMap->GetNumNonDirVertexModes();
152 
153  //Vertex and edge preconditioner matrices
154  VertBlk = MemoryManager<DNekMat>::
155  AllocateSharedPtr(nNonDirVerts,nNonDirVerts,zero,vertstorage);
156 
157  Array<OneD, NekDouble> vertArray(nNonDirVerts,0.0);
158  Array<OneD, long> VertBlockToUniversalMap(nNonDirVerts,-1);
159 
160  int n_exp = expList->GetNumElmts();
161  int nNonDirEdgeIDs=m_locToGloMap->GetNumNonDirEdges();
162 
163  //set the number of blocks in the matrix
164  Array<OneD,unsigned int> n_blks(1+nNonDirEdgeIDs);
165  n_blks[0]=nNonDirVerts;
166 
167  map<int,int> edgeDirMap;
168  map<int,int> uniqueEdgeMap;
169 
170  //this should be of size total number of local edges
171  Array<OneD, int> edgemodeoffset(nNonDirEdgeIDs,0);
172  Array<OneD, int> edgeglobaloffset(nNonDirEdgeIDs,0);
173 
174  const Array<OneD, const ExpListSharedPtr>& bndCondExp = expList->GetBndCondExpansions();
175  StdRegions::StdExpansion1DSharedPtr bndCondFaceExp;
176  LocalRegions::SegExpSharedPtr bndSegExp;
177  const Array<OneD, const SpatialDomains::BoundaryConditionShPtr>&
178  bndConditions = expList->GetBndConditions();
179 
180  int meshVertId;
181  int meshEdgeId;
182 
183  //Determine which boundary edges and faces have dirichlet values
184  for(i = 0; i < bndCondExp.num_elements(); i++)
185  {
186  cnt = 0;
187  for(j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
188  {
189  bndSegExp = bndCondExp[i]->GetExp(j)
190  ->as<LocalRegions::SegExp>();
191  if (bndConditions[i]->GetBoundaryConditionType() ==
192  SpatialDomains::eDirichlet)
193  {
194  meshEdgeId = (bndSegExp->GetGeom1D())->GetEid();
195  edgeDirMap[meshEdgeId] = 1;
196  }
197  }
198  }
199 
200  int dof=0;
201  int maxEdgeDof=0;
202  int nlocalNonDirEdges=0;
203 
204  int edgematrixlocation=0;
205  int ntotaledgeentries=0;
206 
207  // Loop over all the elements in the domain and compute max edge
208  // DOF. Reduce across all processes to get universal maximum.
209  for(cnt=n=0; n < n_exp; ++n)
210  {
211  nel = expList->GetOffset_Elmt_Id(n);
212 
213  locExpansion = expList->GetExp(nel);
214 
215  for (j = 0; j < locExpansion->GetNedges(); ++j)
216  {
217  dof = locExpansion->GetEdgeNcoeffs(j)-2;
218  maxEdgeDof = (dof > maxEdgeDof ? dof : maxEdgeDof);
219  meshEdgeId = locExpansion->as<LocalRegions::Expansion2D>()
220  ->GetGeom2D()->GetEid(j);
221 
222  if(edgeDirMap.count(meshEdgeId)==0)
223  {
224  if(uniqueEdgeMap.count(meshEdgeId)==0)
225  {
226  uniqueEdgeMap[meshEdgeId]=edgematrixlocation;
227 
228  edgeglobaloffset[edgematrixlocation]+=ntotaledgeentries;
229 
230  edgemodeoffset[edgematrixlocation]=dof*dof;
231 
232  ntotaledgeentries+=dof*dof;
233 
234  n_blks[1+edgematrixlocation++]=dof;
235 
236  }
237 
238  nlocalNonDirEdges+=dof*dof;
239  }
240  }
241  }
242 
243  m_comm = expList->GetComm()->GetRowComm();
244  m_comm->AllReduce(maxEdgeDof, LibUtilities::ReduceMax);
245 
246  //Allocate arrays for block to universal map (number of expansions * p^2)
247  Array<OneD, long> EdgeBlockToUniversalMap(ntotaledgeentries,-1);
248 
249  Array<OneD, int> localEdgeToGlobalMatrixMap(nlocalNonDirEdges,-1);
250 
251  //Allocate arrays to store matrices (number of expansions * p^2)
252  Array<OneD, NekDouble> EdgeBlockArray(nlocalNonDirEdges,-1);
253 
254  int edgematrixoffset=0;
255  int vGlobal;
256 
257  for(n=0; n < n_exp; ++n)
258  {
259  nel = expList->GetOffset_Elmt_Id(n);
260 
261  locExpansion = expList->GetExp(nel);
262 
263  //loop over the edges of the expansion
264  for(j = 0; j < locExpansion->GetNedges(); ++j)
265  {
266  //get mesh edge id
267  meshEdgeId = locExpansion->as<LocalRegions::Expansion2D>()
268  ->GetGeom2D()->GetEid(j);
269 
270  nedgemodes=locExpansion->GetEdgeNcoeffs(j)-2;
271 
272  if(edgeDirMap.count(meshEdgeId)==0)
273  {
274  for(k=0; k<nedgemodes*nedgemodes; ++k)
275  {
276  vGlobal=edgeglobaloffset[uniqueEdgeMap[meshEdgeId]]+k;
277 
278 
279  localEdgeToGlobalMatrixMap[edgematrixoffset+k]=vGlobal;
280 
281  EdgeBlockToUniversalMap[vGlobal]
282  = meshEdgeId * maxEdgeDof * maxEdgeDof + k + 1;
283  }
284  edgematrixoffset+=nedgemodes*nedgemodes;
285  }
286  }
287  }
288 
289  edgematrixoffset=0;
290 
291  m_blkMat = MemoryManager<DNekBlkMat>
292  ::AllocateSharedPtr(n_blks, n_blks, blkmatStorage);
293 
294  //Here we loop over the expansion and build the block low energy
295  //preconditioner as well as the block versions of the transformation
296  //matrices.
297  for(cnt=n=0; n < n_exp; ++n)
298  {
299  nel = expList->GetOffset_Elmt_Id(n);
300 
301  locExpansion = expList->GetExp(nel);
302 
303  nVerts=locExpansion->GetGeom()->GetNumVerts();
304  nEdges=locExpansion->GetGeom()->GetNumEdges();
305 
306  //Get statically condensed matrix
307  loc_mat = (m_linsys.lock())->GetStaticCondBlock(n);
308 
309  //Extract boundary block (elemental S1)
310  bnd_mat=loc_mat->GetBlock(0,0);
311 
312  //offset by number of rows
313  offset = bnd_mat->GetRows();
314 
315  DNekScalMat &S=(*bnd_mat);
316 
317  //loop over vertices of the element and return the vertex map
318  //for each vertex
319  for (v=0; v<nVerts; ++v)
320  {
321  vMap1=locExpansion->GetVertexMap(v);
322 
323  //Get vertex map
324  globalrow = m_locToGloMap->
325  GetLocalToGlobalBndMap(cnt+vMap1)-nDirBnd;
326 
327  if(globalrow >= 0)
328  {
329  for (m=0; m<nVerts; ++m)
330  {
331  vMap2=locExpansion->GetVertexMap(m);
332 
333  //global matrix location (without offset due to
334  //dirichlet values)
335  globalcol = m_locToGloMap->
336  GetLocalToGlobalBndMap(cnt+vMap2)-nDirBnd;
337 
338  //offset for dirichlet conditions
339  if (globalcol == globalrow)
340  {
341  meshVertId = locExpansion->as<LocalRegions::Expansion2D>()->GetGeom2D()->GetVid(v);
342 
343  //modal connectivity between elements
344  sign1 = m_locToGloMap->
345  GetLocalToGlobalBndSign(cnt + vMap1);
346  sign2 = m_locToGloMap->
347  GetLocalToGlobalBndSign(cnt + vMap2);
348 
349  vertArray[globalrow]
350  += sign1*sign2*S(vMap1,vMap2);
351 
352  VertBlockToUniversalMap[globalrow]
353  = meshVertId * maxEdgeDof * maxEdgeDof + 1;
354  }
355  }
356  }
357  }
358 
359  //loop over edges of the element and return the edge map
360  for (eid=0; eid<nEdges; ++eid)
361  {
362  nedgemodes=locExpansion->GetEdgeNcoeffs(eid)-2;
363 
364  DNekMatSharedPtr locMat =
365  MemoryManager<DNekMat>::AllocateSharedPtr
366  (nedgemodes,nedgemodes,zero,storage);
367 
368  meshEdgeId = locExpansion->as<LocalRegions::Expansion2D>()->GetGeom2D()->GetEid(eid);
369  Array<OneD, unsigned int> edgemodearray =
370  locExpansion->GetEdgeInverseBoundaryMap(eid);
371 
372  if(edgeDirMap.count(meshEdgeId)==0)
373  {
374  for (v=0; v<nedgemodes; ++v)
375  {
376  eMap1=edgemodearray[v];
377 
378  for (m=0; m<nedgemodes; ++m)
379  {
380  eMap2=edgemodearray[m];
381 
382  //modal connectivity between elements
383  sign1 = m_locToGloMap->
384  GetLocalToGlobalBndSign(cnt + eMap1);
385  sign2 = m_locToGloMap->
386  GetLocalToGlobalBndSign(cnt + eMap2);
387 
388  NekDouble globalEdgeValue =
389  sign1*sign2*S(eMap1,eMap2);
390 
391  EdgeBlockArray[edgematrixoffset+v*nedgemodes+m]=
392  globalEdgeValue;
393  }
394  }
395  edgematrixoffset+=nedgemodes*nedgemodes;
396  }
397  }
398 
399  //offset for the expansion
400  cnt+=offset;
401  }
402 
403  //Assemble edge matrices of each process
404  Array<OneD, NekDouble> GlobalEdgeBlock(ntotaledgeentries);
405  Vmath::Zero(ntotaledgeentries, GlobalEdgeBlock.get(), 1);
406  Vmath::Assmb(EdgeBlockArray.num_elements(),
407  EdgeBlockArray.get(),
408  localEdgeToGlobalMatrixMap.get(),
409  GlobalEdgeBlock.get());
410 
411  //Exchange vertex data over different processes
412  if(nNonDirVerts!=0)
413  {
414  Gs::gs_data *tmp = Gs::Init(VertBlockToUniversalMap, m_comm);
415  Gs::Gather(vertArray, Gs::gs_add, tmp);
416  }
417 
418  //Exchange edge data over different processes
419  Gs::gs_data *tmp1 = Gs::Init(EdgeBlockToUniversalMap, m_comm);
420  Gs::Gather(GlobalEdgeBlock, Gs::gs_add, tmp1);
421 
422  // Populate vertex block
423  for (int i = 0; i < nNonDirVerts; ++i)
424  {
425  VertBlk->SetValue(i,i,1.0/vertArray[i]);
426  }
427 
428  //Set the first block to be the diagonal of the vertex space
429  m_blkMat->SetBlock(0,0, VertBlk);
430 
431  offset=0;
432  //Build the edge matrices from the vector
433  for(int loc=0; loc<nNonDirEdgeIDs; ++loc)
434  {
435  DNekMatSharedPtr gmat =
436  MemoryManager<DNekMat>::AllocateSharedPtr
437  (nedgemodes,nedgemodes,zero,storage);
438 
439  for (v=0; v<nedgemodes; ++v)
440  {
441  for (m=0; m<nedgemodes; ++m)
442  {
443  NekDouble EdgeValue = GlobalEdgeBlock[offset+v*nedgemodes+m];
444  gmat->SetValue(v,m,EdgeValue);
445  }
446  }
447 
448  m_blkMat->SetBlock(1+loc,1+loc, gmat);
449 
450  offset+=edgemodeoffset[loc];
451  }
452 
453  int totblks=m_blkMat->GetNumberOfBlockRows();
454  for (i=1; i< totblks; ++i)
455  {
456  unsigned int nmodes=m_blkMat->GetNumberOfRowsInBlockRow(i);
457  DNekMatSharedPtr tmp_mat =
458  MemoryManager<DNekMat>::AllocateSharedPtr
459  (nmodes,nmodes,zero,storage);
460 
461  tmp_mat=m_blkMat->GetBlock(i,i);
462  tmp_mat->Invert();
463  m_blkMat->SetBlock(i,i,tmp_mat);
464  }
465  }
466 
467  /**
468  *
469  */
470  void PreconditionerBlock::BlockPreconditioner3D()
471  {
472  boost::shared_ptr<MultiRegions::ExpList>
473  expList=((m_linsys.lock())->GetLocMat()).lock();
474  LocalRegions::ExpansionSharedPtr locExpansion;
475  GlobalLinSysKey m_linSysKey=(m_linsys.lock())->GetKey();
476  StdRegions::VarCoeffMap vVarCoeffMap;
477  int i, j, k, nel;
478  int nVerts, nEdges,nFaces;
479  int eid, fid, n, cnt, nedgemodes, nfacemodes;
480  NekDouble zero = 0.0;
481 
482  int vMap1, vMap2, sign1, sign2;
483  int m, v, eMap1, eMap2, fMap1, fMap2;
484  int offset, globalrow, globalcol;
485 
486  //matrix storage
487  MatrixStorage storage = eFULL;
488  MatrixStorage vertstorage = eDIAGONAL;
489  MatrixStorage blkmatStorage = eDIAGONAL;
490 
491  //local element static condensed matrices
492  DNekScalBlkMatSharedPtr loc_mat;
493  DNekScalMatSharedPtr bnd_mat;
494 
495  DNekMatSharedPtr VertBlk;
496 
497  int nDirBnd = m_locToGloMap->GetNumGlobalDirBndCoeffs();
498  int nNonDirVerts = m_locToGloMap->GetNumNonDirVertexModes();
499 
500  //Vertex, edge and face preconditioner matrices
501  VertBlk = MemoryManager<DNekMat>::
502  AllocateSharedPtr(nNonDirVerts,nNonDirVerts,zero,vertstorage);
503 
504  Array<OneD, NekDouble> vertArray(nNonDirVerts,0.0);
505  Array<OneD, long> VertBlockToUniversalMap(nNonDirVerts,-1);
506 
507  int n_exp = expList->GetNumElmts();
508  int nNonDirEdgeIDs=m_locToGloMap->GetNumNonDirEdges();
509  int nNonDirFaceIDs=m_locToGloMap->GetNumNonDirFaces();
510 
511  //set the number of blocks in the matrix
512  Array<OneD,unsigned int> n_blks(1+nNonDirEdgeIDs+nNonDirFaceIDs);
513  n_blks[0]=nNonDirVerts;
514 
515  map<int,int> edgeDirMap;
516  map<int,int> faceDirMap;
517  map<int,int> uniqueEdgeMap;
518  map<int,int> uniqueFaceMap;
519 
520  //this should be of size total number of local edges
521  Array<OneD, int> edgemodeoffset(nNonDirEdgeIDs,0);
522  Array<OneD, int> facemodeoffset(nNonDirFaceIDs,0);
523 
524  Array<OneD, int> edgeglobaloffset(nNonDirEdgeIDs,0);
525  Array<OneD, int> faceglobaloffset(nNonDirFaceIDs,0);
526 
527  const Array<OneD, const ExpListSharedPtr>& bndCondExp = expList->GetBndCondExpansions();
528  StdRegions::StdExpansion2DSharedPtr bndCondFaceExp;
529  const Array<OneD, const SpatialDomains::BoundaryConditionShPtr>& bndConditions = expList->GetBndConditions();
530 
531  int meshVertId;
532  int meshEdgeId;
533  int meshFaceId;
534 
535  const Array<OneD, const int> &extradiredges
536  = m_locToGloMap->GetExtraDirEdges();
537  for(i=0; i<extradiredges.num_elements(); ++i)
538  {
539  meshEdgeId=extradiredges[i];
540  edgeDirMap[meshEdgeId] = 1;
541  }
542 
543  //Determine which boundary edges and faces have dirichlet values
544  for(i = 0; i < bndCondExp.num_elements(); i++)
545  {
546  cnt = 0;
547  for(j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
548  {
549  bndCondFaceExp = bndCondExp[i]->GetExp(j)->as<StdRegions::StdExpansion2D>();
550  if (bndConditions[i]->GetBoundaryConditionType() ==
551  SpatialDomains::eDirichlet)
552  {
553  for(k = 0; k < bndCondFaceExp->GetNedges(); k++)
554  {
555  meshEdgeId = bndCondFaceExp->as<LocalRegions::Expansion2D>()->GetGeom2D()->GetEid(k);
556  if(edgeDirMap.count(meshEdgeId) == 0)
557  {
558  edgeDirMap[meshEdgeId] = 1;
559  }
560  }
561  meshFaceId = bndCondFaceExp->as<LocalRegions::Expansion2D>()->GetGeom2D()->GetFid();
562  faceDirMap[meshFaceId] = 1;
563  }
564  }
565  }
566 
567  int dof=0;
568  int maxFaceDof=0;
569  int maxEdgeDof=0;
570  int nlocalNonDirEdges=0;
571  int nlocalNonDirFaces=0;
572 
573  int edgematrixlocation=0;
574  int ntotaledgeentries=0;
575 
576  // Loop over all the elements in the domain and compute max edge
577  // DOF. Reduce across all processes to get universal maximum.
578  for(cnt=n=0; n < n_exp; ++n)
579  {
580  nel = expList->GetOffset_Elmt_Id(n);
581 
582  locExpansion = expList->GetExp(nel);
583 
584  for (j = 0; j < locExpansion->GetNedges(); ++j)
585  {
586  dof = locExpansion->GetEdgeNcoeffs(j)-2;
587  maxEdgeDof = (dof > maxEdgeDof ? dof : maxEdgeDof);
588  meshEdgeId = locExpansion->as<LocalRegions::Expansion3D>()->GetGeom3D()->GetEid(j);
589 
590  if(edgeDirMap.count(meshEdgeId)==0)
591  {
592  if(uniqueEdgeMap.count(meshEdgeId)==0 && dof > 0)
593  {
594  uniqueEdgeMap[meshEdgeId]=edgematrixlocation;
595 
596  edgeglobaloffset[edgematrixlocation]+=ntotaledgeentries;
597 
598  edgemodeoffset[edgematrixlocation]=dof*dof;
599 
600  ntotaledgeentries+=dof*dof;
601 
602  n_blks[1+edgematrixlocation++]=dof;
603  }
604 
605  nlocalNonDirEdges+=dof*dof;
606  }
607  }
608  }
609 
610  int facematrixlocation=0;
611  int ntotalfaceentries=0;
612 
613  // Loop over all the elements in the domain and compute max face
614  // DOF. Reduce across all processes to get universal maximum.
615  for(cnt=n=0; n < n_exp; ++n)
616  {
617  nel = expList->GetOffset_Elmt_Id(n);
618 
619  locExpansion = expList->GetExp(nel);
620 
621  for (j = 0; j < locExpansion->GetNfaces(); ++j)
622  {
623  dof = locExpansion->GetFaceIntNcoeffs(j);
624  maxFaceDof = (dof > maxFaceDof ? dof : maxFaceDof);
625 
626  meshFaceId = locExpansion->as<LocalRegions::Expansion3D>()->GetGeom3D()->GetFid(j);
627 
628  if(faceDirMap.count(meshFaceId)==0)
629  {
630  if(uniqueFaceMap.count(meshFaceId)==0 && dof > 0)
631  {
632  uniqueFaceMap[meshFaceId]=facematrixlocation;
633 
634  facemodeoffset[facematrixlocation]=dof*dof;
635 
636  faceglobaloffset[facematrixlocation]+=ntotalfaceentries;
637 
638  ntotalfaceentries+=dof*dof;
639 
640  n_blks[1+nNonDirEdgeIDs+facematrixlocation++]=dof;
641  }
642  nlocalNonDirFaces+=dof*dof;
643  }
644 
645  }
646  }
647 
648  m_comm = expList->GetComm();
649  m_comm->AllReduce(maxEdgeDof, LibUtilities::ReduceMax);
650  m_comm->AllReduce(maxFaceDof, LibUtilities::ReduceMax);
651 
652  //Allocate arrays for block to universal map (number of expansions * p^2)
653  Array<OneD, long> EdgeBlockToUniversalMap(ntotaledgeentries,-1);
654  Array<OneD, long> FaceBlockToUniversalMap(ntotalfaceentries,-1);
655 
656  Array<OneD, int> localEdgeToGlobalMatrixMap(nlocalNonDirEdges,-1);
657  Array<OneD, int> localFaceToGlobalMatrixMap(nlocalNonDirFaces,-1);
658 
659  //Allocate arrays to store matrices (number of expansions * p^2)
660  Array<OneD, NekDouble> EdgeBlockArray(nlocalNonDirEdges,-1);
661  Array<OneD, NekDouble> FaceBlockArray(nlocalNonDirFaces,-1);
662 
663  int edgematrixoffset=0;
664  int facematrixoffset=0;
665  int vGlobal;
666 
667  for(n=0; n < n_exp; ++n)
668  {
669  nel = expList->GetOffset_Elmt_Id(n);
670 
671  locExpansion = expList->GetExp(nel);
672 
673  //loop over the edges of the expansion
674  for(j = 0; j < locExpansion->GetNedges(); ++j)
675  {
676  //get mesh edge id
677  meshEdgeId = locExpansion->as<LocalRegions::Expansion3D>()->GetGeom3D()->GetEid(j);
678 
679  nedgemodes=locExpansion->GetEdgeNcoeffs(j)-2;
680 
681  if(edgeDirMap.count(meshEdgeId)==0)
682  {
683  for(k=0; k<nedgemodes*nedgemodes; ++k)
684  {
685  vGlobal=edgeglobaloffset[uniqueEdgeMap[meshEdgeId]]+k;
686 
687 
688  localEdgeToGlobalMatrixMap[edgematrixoffset+k]=vGlobal;
689 
690  EdgeBlockToUniversalMap[vGlobal]
691  = meshEdgeId * maxEdgeDof * maxEdgeDof + k + 1;
692  }
693  edgematrixoffset+=nedgemodes*nedgemodes;
694  }
695  }
696 
697  //loop over the faces of the expansion
698  for(j = 0; j < locExpansion->GetNfaces(); ++j)
699  {
700  //get mesh face id
701  meshFaceId = locExpansion->as<LocalRegions::Expansion3D>()->GetGeom3D()->GetFid(j);
702 
703  nfacemodes = locExpansion->GetFaceIntNcoeffs(j);
704 
705  //Check if face is has dirichlet values
706  if(faceDirMap.count(meshFaceId)==0)
707  {
708  for(k=0; k<nfacemodes*nfacemodes; ++k)
709  {
710  vGlobal=faceglobaloffset[uniqueFaceMap[meshFaceId]]+k;
711 
712  localFaceToGlobalMatrixMap[facematrixoffset+k]
713  = vGlobal;
714 
715  FaceBlockToUniversalMap[vGlobal]
716  = meshFaceId * maxFaceDof * maxFaceDof + k + 1;
717  }
718  facematrixoffset+=nfacemodes*nfacemodes;
719  }
720  }
721  }
722 
723  edgematrixoffset=0;
724  facematrixoffset=0;
725 
726  m_blkMat = MemoryManager<DNekBlkMat>
727  ::AllocateSharedPtr(n_blks, n_blks, blkmatStorage);
728 
729  //Here we loop over the expansion and build the block low energy
730  //preconditioner as well as the block versions of the transformation
731  //matrices.
732  for(cnt=n=0; n < n_exp; ++n)
733  {
734  nel = expList->GetOffset_Elmt_Id(n);
735 
736  locExpansion = expList->GetExp(nel);
737 
738  nVerts=locExpansion->GetGeom()->GetNumVerts();
739  nEdges=locExpansion->GetGeom()->GetNumEdges();
740  nFaces=locExpansion->GetGeom()->GetNumFaces();
741 
742  //Get statically condensed matrix
743  loc_mat = (m_linsys.lock())->GetStaticCondBlock(n);
744 
745  //Extract boundary block (elemental S1)
746  bnd_mat=loc_mat->GetBlock(0,0);
747 
748  //offset by number of rows
749  offset = bnd_mat->GetRows();
750 
751  DNekScalMat &S=(*bnd_mat);
752 
753  //loop over vertices of the element and return the vertex map
754  //for each vertex
755  for (v=0; v<nVerts; ++v)
756  {
757  vMap1=locExpansion->GetVertexMap(v);
758 
759  //Get vertex map
760  globalrow = m_locToGloMap->
761  GetLocalToGlobalBndMap(cnt+vMap1)-nDirBnd;
762 
763  if(globalrow >= 0)
764  {
765  for (m=0; m<nVerts; ++m)
766  {
767  vMap2=locExpansion->GetVertexMap(m);
768 
769  //global matrix location (without offset due to
770  //dirichlet values)
771  globalcol = m_locToGloMap->
772  GetLocalToGlobalBndMap(cnt+vMap2)-nDirBnd;
773 
774  //offset for dirichlet conditions
775  if (globalcol == globalrow)
776  {
777  meshVertId = locExpansion->as<LocalRegions::Expansion3D>()->GetGeom3D()->GetVid(v);
778 
779  //modal connectivity between elements
780  sign1 = m_locToGloMap->
781  GetLocalToGlobalBndSign(cnt + vMap1);
782  sign2 = m_locToGloMap->
783  GetLocalToGlobalBndSign(cnt + vMap2);
784 
785  vertArray[globalrow]
786  += sign1*sign2*S(vMap1,vMap2);
787 
788  VertBlockToUniversalMap[globalrow]
789  = meshVertId * maxEdgeDof * maxEdgeDof + 1;
790  }
791  }
792  }
793  }
794 
795  //loop over edges of the element and return the edge map
796  for (eid=0; eid<nEdges; ++eid)
797  {
798  nedgemodes=locExpansion->GetEdgeNcoeffs(eid)-2;
799 
800  DNekMatSharedPtr locMat =
801  MemoryManager<DNekMat>::AllocateSharedPtr
802  (nedgemodes,nedgemodes,zero,storage);
803 
804  meshEdgeId = locExpansion->as<LocalRegions::Expansion3D>()->GetGeom3D()->GetEid(eid);
805  Array<OneD, unsigned int> edgemodearray = locExpansion->GetEdgeInverseBoundaryMap(eid);
806 
807  if(edgeDirMap.count(meshEdgeId)==0)
808  {
809  for (v=0; v<nedgemodes; ++v)
810  {
811  eMap1=edgemodearray[v];
812 
813  for (m=0; m<nedgemodes; ++m)
814  {
815  eMap2=edgemodearray[m];
816 
817  //modal connectivity between elements
818  sign1 = m_locToGloMap->
819  GetLocalToGlobalBndSign(cnt + eMap1);
820  sign2 = m_locToGloMap->
821  GetLocalToGlobalBndSign(cnt + eMap2);
822 
823  NekDouble globalEdgeValue = sign1*sign2*S(eMap1,eMap2);
824 
825  EdgeBlockArray[edgematrixoffset+v*nedgemodes+m]=globalEdgeValue;
826  }
827  }
828  edgematrixoffset+=nedgemodes*nedgemodes;
829  }
830  }
831 
832  //loop over faces of the element and return the face map
833  for (fid=0; fid<nFaces; ++fid)
834  {
835  nfacemodes = locExpansion->GetFaceIntNcoeffs(fid);
836 
837  DNekMatSharedPtr locMat =
838  MemoryManager<DNekMat>::AllocateSharedPtr
839  (nfacemodes,nfacemodes,zero,storage);
840 
841  meshFaceId = locExpansion->as<LocalRegions::Expansion3D>()->GetGeom3D()->GetFid(fid);
842 
843  Array<OneD, unsigned int> facemodearray = locExpansion->GetFaceInverseBoundaryMap(fid);
844 
845  if(faceDirMap.count(meshFaceId)==0)
846  {
847  for (v=0; v<nfacemodes; ++v)
848  {
849  fMap1=facemodearray[v];
850 
851  for (m=0; m<nfacemodes; ++m)
852  {
853  fMap2=facemodearray[m];
854 
855  //modal connectivity between elements
856  sign1 = m_locToGloMap->
857  GetLocalToGlobalBndSign(cnt + fMap1);
858  sign2 = m_locToGloMap->
859  GetLocalToGlobalBndSign(cnt + fMap2);
860 
861  // Get the face-face value from the low energy matrix (S2)
862  NekDouble globalFaceValue = sign1*sign2*S(fMap1,fMap2);
863 
864  //local face value to global face value
865  FaceBlockArray[facematrixoffset+v*nfacemodes+m]=globalFaceValue;
866  }
867  }
868  facematrixoffset+=nfacemodes*nfacemodes;
869  }
870  }
871 
872  //offset for the expansion
873  cnt+=offset;
874  }
875 
876  //Assemble edge matrices of each process
877  Array<OneD, NekDouble> GlobalEdgeBlock(ntotaledgeentries);
878  Vmath::Zero(ntotaledgeentries, GlobalEdgeBlock.get(), 1);
879  Vmath::Assmb(EdgeBlockArray.num_elements(),
880  EdgeBlockArray.get(),
881  localEdgeToGlobalMatrixMap.get(),
882  GlobalEdgeBlock.get());
883 
884  //Assemble face matrices of each process
885  Array<OneD, NekDouble> GlobalFaceBlock(ntotalfaceentries);
886  Vmath::Zero(ntotalfaceentries, GlobalFaceBlock.get(), 1);
887  Vmath::Assmb(FaceBlockArray.num_elements(),
888  FaceBlockArray.get(),
889  localFaceToGlobalMatrixMap.get(),
890  GlobalFaceBlock.get());
891 
892  //Exchange vertex data over different processes
893  if(nNonDirVerts!=0)
894  {
895  Gs::gs_data *tmp = Gs::Init(VertBlockToUniversalMap, m_comm);
896  Gs::Gather(vertArray, Gs::gs_add, tmp);
897  }
898 
899  //Exchange edge data over different processes
900  Gs::gs_data *tmp1 = Gs::Init(EdgeBlockToUniversalMap, m_comm);
901  Gs::Gather(GlobalEdgeBlock, Gs::gs_add, tmp1);
902 
903  //Exchange face data over different processes
904  Gs::gs_data *tmp2 = Gs::Init(FaceBlockToUniversalMap, m_comm);
905  Gs::Gather(GlobalFaceBlock, Gs::gs_add, tmp2);
906 
907  // Populate vertex block
908  for (int i = 0; i < nNonDirVerts; ++i)
909  {
910  VertBlk->SetValue(i,i,1.0/vertArray[i]);
911  }
912 
913  //Set the first block to be the diagonal of the vertex space
914  m_blkMat->SetBlock(0,0, VertBlk);
915 
916  offset=0;
917  //Build the edge matrices from the vector
918  for(int loc=0; loc<nNonDirEdgeIDs; ++loc)
919  {
920  DNekMatSharedPtr gmat =
921  MemoryManager<DNekMat>::AllocateSharedPtr
922  (nedgemodes,nedgemodes,zero,storage);
923 
924  for (v=0; v<nedgemodes; ++v)
925  {
926  for (m=0; m<nedgemodes; ++m)
927  {
928  NekDouble EdgeValue = GlobalEdgeBlock[offset+v*nedgemodes+m];
929  gmat->SetValue(v,m,EdgeValue);
930  }
931  }
932 
933  m_blkMat->SetBlock(1+loc,1+loc, gmat);
934 
935  offset+=edgemodeoffset[loc];
936  }
937 
938  offset=0;
939  //Build the face matrices from the vector
940  for(int loc=0; loc<nNonDirFaceIDs; ++loc)
941  {
942  nfacemodes=n_blks[1+nNonDirEdgeIDs+loc];
943 
944  DNekMatSharedPtr gmat =
945  MemoryManager<DNekMat>::AllocateSharedPtr
946  (nfacemodes,nfacemodes,zero,storage);
947 
948  for (v=0; v<nfacemodes; ++v)
949  {
950  for (m=0; m<nfacemodes; ++m)
951  {
952  NekDouble FaceValue = GlobalFaceBlock[offset+v*nfacemodes+m];
953  gmat->SetValue(v,m,FaceValue);
954  }
955  }
956 
957  m_blkMat->SetBlock(1+nNonDirEdgeIDs+loc,1+nNonDirEdgeIDs+loc, gmat);
958 
959  offset+=facemodeoffset[loc];
960  }
961 
962 
963  int totblks=m_blkMat->GetNumberOfBlockRows();
964  for (i=1; i< totblks; ++i)
965  {
966  unsigned int nmodes=m_blkMat->GetNumberOfRowsInBlockRow(i);
967  DNekMatSharedPtr tmp_mat =
968  MemoryManager<DNekMat>::AllocateSharedPtr
969  (nmodes,nmodes,zero,storage);
970 
971  tmp_mat=m_blkMat->GetBlock(i,i);
972  tmp_mat->Invert();
973 
974  m_blkMat->SetBlock(i,i,tmp_mat);
975  }
976  }
977 
978  void PreconditionerBlock::BlockPreconditionerHDG()
979  {
980  boost::shared_ptr<MultiRegions::ExpList>
981  expList=((m_linsys.lock())->GetLocMat()).lock();
982  boost::shared_ptr<MultiRegions::ExpList> trace = expList->GetTrace();
983  LocalRegions::ExpansionSharedPtr locExpansion;
984  DNekScalBlkMatSharedPtr loc_mat;
985  DNekScalMatSharedPtr bnd_mat;
986 
987  AssemblyMapDGSharedPtr asmMap = boost::dynamic_pointer_cast<
988  AssemblyMapDG>(m_locToGloMap);
989 
990  int i, j, k, n, cnt, cnt2;
991 
992  // Figure out number of Dirichlet trace elements
993  int nTrace = expList->GetTrace()->GetExpSize();
994  int nDir = m_locToGloMap->GetNumGlobalDirBndCoeffs();
995 
996  for (cnt = n = 0; n < nTrace; ++n)
997  {
998  if (cnt >= nDir)
999  {
1000  break;
1001  }
1002 
1003  cnt += trace->GetExp(n)->GetNcoeffs();
1004  }
1005 
1006  nDir = n;
1007 
1008  // Allocate storage for block matrix. Need number of unique faces in
1009  // trace space.
1010  int maxTraceSize = -1;
1011  map<int, int> arrayOffsets;
1012 
1013  for (cnt = 0, n = nDir; n < nTrace; ++n)
1014  {
1015  int nCoeffs = trace->GetExp(n)->GetNcoeffs();
1016  int nCoeffs2 = nCoeffs * nCoeffs;
1017 
1018  arrayOffsets[n] = cnt;
1019  cnt += nCoeffs2;
1020 
1021  if (maxTraceSize < nCoeffs2)
1022  {
1023  maxTraceSize = nCoeffs2;
1024  }
1025  }
1026 
1027  // Find maximum trace size.
1028  LibUtilities::CommSharedPtr comm =
1029  expList->GetSession()->GetComm()->GetRowComm();
1030  comm->AllReduce(maxTraceSize, LibUtilities::ReduceMax);
1031 
1032  // Zero matrix storage.
1033  Array<OneD, NekDouble> tmpStore(cnt, 0.0);
1034 
1035  // Assemble block matrices for each trace element.
1036  for (cnt = n = 0; n < expList->GetExpSize(); ++n)
1037  {
1038  int elmt = expList->GetOffset_Elmt_Id(n);
1039  locExpansion = expList->GetExp(elmt);
1040 
1041  Array<OneD, LocalRegions::ExpansionSharedPtr> &elmtToTraceMap =
1042  asmMap->GetElmtToTrace()[elmt];
1043 
1044  // Block matrix (lambda)
1045  loc_mat = (m_linsys.lock())->GetStaticCondBlock(n);
1046  bnd_mat = loc_mat->GetBlock(0,0);
1047 
1048  int nFacets = locExpansion->GetNumBases() == 2 ?
1049  locExpansion->GetNedges() : locExpansion->GetNfaces();
1050 
1051  for (cnt2 = i = 0; i < nFacets; ++i)
1052  {
1053  int nCoeffs = elmtToTraceMap[i]->GetNcoeffs();
1054  int elmtId = elmtToTraceMap[i]->GetElmtId ();
1055 
1056  // Ignore Dirichlet edges.
1057  if (elmtId < nDir)
1058  {
1059  cnt += nCoeffs;
1060  cnt2 += nCoeffs;
1061  continue;
1062  }
1063 
1064  NekDouble *off = &tmpStore[arrayOffsets[elmtId]];
1065 
1066  for (j = 0; j < nCoeffs; ++j)
1067  {
1068  NekDouble sign1 = asmMap->GetLocalToGlobalBndSign(
1069  cnt + j);
1070  for (k = 0; k < nCoeffs; ++k)
1071  {
1072  NekDouble sign2 = asmMap->GetLocalToGlobalBndSign(
1073  cnt + k);
1074  off[j*nCoeffs + k] +=
1075  (*bnd_mat)(cnt2+j, cnt2+k) * sign1 * sign2;
1076  }
1077  }
1078 
1079  cnt += nCoeffs;
1080  cnt2 += nCoeffs;
1081  }
1082  }
1083 
1084  // Set up IDs for universal numbering.
1085  Array<OneD, long> uniIds(tmpStore.num_elements());
1086  for (cnt = 0, n = nDir; n < nTrace; ++n)
1087  {
1088  LocalRegions::ExpansionSharedPtr traceExp = trace->GetExp(n);
1089  int nCoeffs = traceExp->GetNcoeffs();
1090  int geomId = traceExp->GetGeom()->GetGlobalID();
1091 
1092  for (i = 0; i < nCoeffs*nCoeffs; ++i)
1093  {
1094  uniIds[cnt++] = geomId * maxTraceSize + i + 1;
1095  }
1096  }
1097 
1098  // Assemble matrices across partitions.
1099  Gs::gs_data *gsh = Gs::Init(uniIds, comm);
1100  Gs::Gather(tmpStore, Gs::gs_add, gsh);
1101 
1102  // Set up diagonal block matrix
1103  Array<OneD, unsigned int> n_blks(nTrace - nDir);
1104  for (n = 0; n < nTrace - nDir; ++n)
1105  {
1106  n_blks[n] = trace->GetExp(n + nDir)->GetNcoeffs();
1107  }
1108 
1109  m_blkMat = MemoryManager<DNekBlkMat>
1110  ::AllocateSharedPtr(n_blks, n_blks, eDIAGONAL);
1111 
1112  for (n = 0; n < nTrace - nDir; ++n)
1113  {
1114  int nCoeffs = trace->GetExp(n + nDir)->GetNcoeffs();
1115  DNekMatSharedPtr tmp = MemoryManager<DNekMat>
1116  ::AllocateSharedPtr(nCoeffs, nCoeffs);
1117  NekDouble *off = &tmpStore[arrayOffsets[n+nDir]];
1118 
1119  for (i = 0; i < nCoeffs; ++i)
1120  {
1121  for (j = 0; j < nCoeffs; ++j)
1122  {
1123  (*tmp)(i,j) = off[i*nCoeffs + j];
1124  }
1125  }
1126 
1127  tmp->Invert();
1128  m_blkMat->SetBlock(n, n, tmp);
1129  }
1130  }
1131 
1132  /**
1133  *
1134  */
1135  void PreconditionerBlock::v_DoPreconditioner(
1136  const Array<OneD, NekDouble>& pInput,
1137  Array<OneD, NekDouble>& pOutput)
1138  {
1139  int nDir = m_locToGloMap->GetNumGlobalDirBndCoeffs();
1140  int nGlobal = m_locToGloMap->GetNumGlobalBndCoeffs();
1141  int nNonDir = nGlobal-nDir;
1142  DNekBlkMat &M = (*m_blkMat);
1143  NekVector<NekDouble> r(nNonDir,pInput,eWrapper);
1144  NekVector<NekDouble> z(nNonDir,pOutput,eWrapper);
1145  z = M * r;
1146  }
1147  }
1148 }
Generated on Wed Mar 4 2015 16:33:16 for Nektar++ by   doxygen 1.8.1.2