Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GlobalLinSysXxtStaticCond.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File GlobalLinSysIterativeStaticCond.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: GlobalLinSysIterativeStaticCond definition
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
40 
41 using namespace std;
42 
43 namespace Nektar
44 {
45  namespace MultiRegions
46  {
47  /**
48  * @class GlobalLinSysIterativeStaticCond
49  *
50  * Solves a linear system iteratively using single- or multi-level
51  * static condensation.
52  */
53 
54  /**
55  * Registers the class with the Factory.
56  */
57  string GlobalLinSysXxtStaticCond::className
59  "XxtStaticCond",
60  GlobalLinSysXxtStaticCond::create,
61  "Iterative static condensation.");
62 
63  string GlobalLinSysXxtStaticCond::className2
65  "XxtMultiLevelStaticCond",
66  GlobalLinSysXxtStaticCond::create,
67  "Xxt multi-level static condensation.");
68 
69  /**
70  * For a matrix system of the form @f[
71  * \left[ \begin{array}{cc}
72  * \boldsymbol{A} & \boldsymbol{B}\\
73  * \boldsymbol{C} & \boldsymbol{D}
74  * \end{array} \right]
75  * \left[ \begin{array}{c} \boldsymbol{x_1}\\ \boldsymbol{x_2}
76  * \end{array}\right]
77  * = \left[ \begin{array}{c} \boldsymbol{y_1}\\ \boldsymbol{y_2}
78  * \end{array}\right],
79  * @f]
80  * where @f$\boldsymbol{D}@f$ and
81  * @f$(\boldsymbol{A-BD^{-1}C})@f$ are invertible, store and assemble
82  * a static condensation system, according to a given local to global
83  * mapping. #m_linSys is constructed by AssembleSchurComplement().
84  * @param mKey Associated matrix key.
85  * @param pLocMatSys LocalMatrixSystem
86  * @param locToGloMap Local to global mapping.
87  */
88  GlobalLinSysXxtStaticCond::GlobalLinSysXxtStaticCond(
89  const GlobalLinSysKey &pKey,
90  const boost::weak_ptr<ExpList> &pExpList,
91  const boost::shared_ptr<AssemblyMap>
92  &pLocToGloMap)
93  : GlobalLinSys (pKey, pExpList, pLocToGloMap),
94  GlobalLinSysXxt (pKey, pExpList, pLocToGloMap),
95  GlobalLinSysStaticCond(pKey, pExpList, pLocToGloMap)
96  {
99  "This constructor is only valid when using static "
100  "condensation");
102  == pLocToGloMap->GetGlobalSysSolnType(),
103  "The local to global map is not set up for the requested "
104  "solution type");
105  }
106 
107 
108  /**
109  *
110  */
112  const GlobalLinSysKey &pKey,
113  const boost::weak_ptr<ExpList> &pExpList,
114  const DNekScalBlkMatSharedPtr pSchurCompl,
115  const DNekScalBlkMatSharedPtr pBinvD,
116  const DNekScalBlkMatSharedPtr pC,
117  const DNekScalBlkMatSharedPtr pInvD,
118  const boost::shared_ptr<AssemblyMap>
119  &pLocToGloMap)
120  : GlobalLinSys (pKey, pExpList, pLocToGloMap),
121  GlobalLinSysXxt (pKey, pExpList, pLocToGloMap),
122  GlobalLinSysStaticCond(pKey, pExpList, pLocToGloMap)
123  {
124  m_schurCompl = pSchurCompl;
125  m_BinvD = pBinvD;
126  m_C = pC;
127  m_invD = pInvD;
128  m_locToGloMap = pLocToGloMap;
129  }
130 
131 
132  /**
133  *
134  */
136  {
137 
138  }
139 
140  /**
141  * Create the inverse multiplicity map.
142  * @param locToGloMap Local to global mapping information.
143  */
145  const boost::shared_ptr<AssemblyMap> &pLocToGloMap)
146  {
147  const Array<OneD, const int> &vMap
148  = pLocToGloMap->GetLocalToGlobalBndMap();
149  unsigned int nGlo = pLocToGloMap->GetNumGlobalBndCoeffs();
150  unsigned int nEntries = pLocToGloMap->GetNumLocalBndCoeffs();
151  unsigned int i;
152 
153  // Count the multiplicity of each global DOF on this process
154  Array<OneD, NekDouble> vCounts(nGlo, 0.0);
155  for (i = 0; i < nEntries; ++i)
156  {
157  vCounts[vMap[i]] += 1.0;
158  }
159 
160  // Get universal multiplicity by globally assembling counts
161  pLocToGloMap->UniversalAssembleBnd(vCounts);
162 
163  // Construct a map of 1/multiplicity for use in XXT solve
165  for (i = 0; i < nEntries; ++i)
166  {
167  m_locToGloSignMult[i] = 1.0/vCounts[vMap[i]];
168  }
169 
170  m_map = pLocToGloMap->GetLocalToGlobalBndMap();
171  }
172 
173  /**
174  * Construct the local matrix row index, column index and value index
175  * arrays and initialize the XXT data structure with this information.
176  * @param locToGloMap Local to global mapping information.
177  */
179  boost::shared_ptr<AssemblyMap> pLocToGloMap)
180  {
181  CreateMap(pLocToGloMap);
182 
183  ExpListSharedPtr vExp = m_expList.lock();
184  unsigned int nElmt = m_schurCompl->GetNumberOfBlockRows();
185  DNekScalMatSharedPtr loc_mat;
186  unsigned int iCount = 0;
187  unsigned int rCount = 0;
188  unsigned int nRows = 0;
189  unsigned int nEntries = 0;
190  unsigned int numDirBnd = pLocToGloMap->GetNumGlobalDirBndCoeffs();
191  unsigned int nLocal = pLocToGloMap->GetNumLocalBndCoeffs();
192  const Array<OneD, NekDouble> &vMapSign
193  = pLocToGloMap->GetLocalToGlobalBndSign();
194  bool doSign = pLocToGloMap->GetSignChange();
195  unsigned int i = 0, j = 0, k = 0, n = 0;
196  int gid1;
197  Array<OneD, unsigned int> vSizes(nElmt);
198 
199  // First construct a map of the number of local DOFs in each block
200  // and the number of matrix entries for each block
201  for (n = 0; n < nElmt; ++n)
202  {
203  loc_mat = m_schurCompl->GetBlock(n,n);
204  vSizes[n] = loc_mat->GetRows();
205  nEntries += vSizes[n]*vSizes[n];
206  }
207 
208  // Set up i-index, j-index and value arrays
209  m_Ai = Array<OneD, unsigned int>(nEntries);
210  m_Aj = Array<OneD, unsigned int>(nEntries);
211  m_Ar = Array<OneD, double>(nEntries, 0.0);
212 
213  // Set up the universal ID array for XXT
214  Array<OneD, unsigned long> vId(nLocal);
215 
216  // Loop over each elemental block, extract matrix indices and value
217  // and set the universal ID array
218  for(n = iCount = 0; n < nElmt; ++n)
219  {
220  loc_mat = m_schurCompl->GetBlock(n,n);
221  nRows = loc_mat->GetRows();
222 
223  for(i = 0; i < nRows; ++i)
224  {
225  gid1 = pLocToGloMap->GetLocalToGlobalBndMap(iCount + i);
226  for(j = 0; j < nRows; ++j)
227  {
228  k = rCount + i*vSizes[n] + j;
229  m_Ai[k] = iCount + i;
230  m_Aj[k] = iCount + j;
231  m_Ar[k] = (*loc_mat)(i,j);
232  if (doSign)
233  {
234  m_Ar[k] *= vMapSign[iCount+i]*vMapSign[iCount+j];
235  }
236  }
237 
238  // Dirichlet DOFs are not included in the solve, so we set
239  // these to the special XXT id=0.
240  if (gid1 < numDirBnd)
241  {
242  vId[iCount + i] = 0;
243  }
244  else
245  {
246  vId[iCount + i]
247  = pLocToGloMap->GetGlobalToUniversalBndMap()[gid1];
248  }
249  }
250  iCount += vSizes[n];
251  rCount += vSizes[n]*vSizes[n];
252  }
253 
254  // Set up XXT and output some stats
255  LibUtilities::CommSharedPtr vComm = pLocToGloMap->GetComm()->GetRowComm();
256  m_crsData = Xxt::Init(nLocal, vId, m_Ai, m_Aj, m_Ar, vComm);
258  }
259 
261  const GlobalLinSysKey &mkey,
262  const boost::weak_ptr<ExpList> &pExpList,
263  const DNekScalBlkMatSharedPtr pSchurCompl,
264  const DNekScalBlkMatSharedPtr pBinvD,
265  const DNekScalBlkMatSharedPtr pC,
266  const DNekScalBlkMatSharedPtr pInvD,
267  const boost::shared_ptr<AssemblyMap> &l2gMap)
268  {
270  GlobalLinSysXxtStaticCond>::AllocateSharedPtr(
271  mkey, pExpList, pSchurCompl, pBinvD, pC, pInvD, l2gMap);
272  sys->Initialise(l2gMap);
273  return sys;
274  }
275  }
276 }
GlobalSysSolnType GetGlobalSysSolnType() const
Return the associated solution type.
DNekScalBlkMatSharedPtr m_schurCompl
Block Schur complement matrix.
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
DNekScalBlkMatSharedPtr m_C
Block matrix.
boost::shared_ptr< GlobalLinSysXxtStaticCond > GlobalLinSysXxtStaticCondSharedPtr
STL namespace.
DNekScalBlkMatSharedPtr m_invD
Block matrix.
virtual void v_AssembleSchurComplement(boost::shared_ptr< AssemblyMap > locToGloMap)
Assemble the Schur complement matrix.
boost::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:53
boost::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:74
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Describe a linear system.
static struct crs_data * Init(unsigned int pRank, const Nektar::Array< OneD, unsigned long > pId, const Nektar::Array< OneD, unsigned int > pAi, const Nektar::Array< OneD, unsigned int > pAj, const Nektar::Array< OneD, NekDouble > pAr, const LibUtilities::CommSharedPtr &pComm)
Initialise the matrix-solve.
Definition: Xxt.hpp:158
GlobalLinSysXxtStaticCond(const GlobalLinSysKey &mkey, const boost::weak_ptr< ExpList > &pExpList, const boost::shared_ptr< AssemblyMap > &locToGloMap)
Constructor for full direct matrix solve.
Array< OneD, unsigned int > m_Aj
A global linear system.
Definition: GlobalLinSys.h:74
virtual GlobalLinSysStaticCondSharedPtr v_Recurse(const GlobalLinSysKey &mkey, const boost::weak_ptr< ExpList > &pExpList, const DNekScalBlkMatSharedPtr pSchurCompl, const DNekScalBlkMatSharedPtr pBinvD, const DNekScalBlkMatSharedPtr pC, const DNekScalBlkMatSharedPtr pInvD, const boost::shared_ptr< AssemblyMap > &locToGloMap)
Array< OneD, unsigned int > m_Ai
boost::shared_ptr< AssemblyMap > m_locToGloMap
Local to global map.
GlobalLinSysFactory & GetGlobalLinSysFactory()
void nektar_crs_stats(struct crs_data *data)
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:218
void CreateMap(const boost::shared_ptr< AssemblyMap > &pLocToGloMap)
DNekScalBlkMatSharedPtr m_BinvD
Block matrix.
Array< OneD, NekDouble > m_locToGloSignMult
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215
boost::shared_ptr< GlobalLinSysStaticCond > GlobalLinSysStaticCondSharedPtr
const boost::weak_ptr< ExpList > m_expList
Local Matrix System.
Definition: GlobalLinSys.h:129