Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GlobalLinSys.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File GlobalLinSys.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: GlobalLinSys definition
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <MultiRegions/GlobalLinSys.h>
37 #include <LocalRegions/MatrixKey.h>
38 #include <LocalRegions/Expansion.h>
39 #include <LibUtilities/BasicUtils/SessionReader.h>
40 
41 #include <LibUtilities/LinearAlgebra/NekTypeDefs.hpp>
42 #include <LibUtilities/LinearAlgebra/NekMatrix.hpp>
43 
44 namespace Nektar
45 {
46  namespace MultiRegions
47  {
48  std::string GlobalLinSys::lookupIds[12] = {
49  LibUtilities::SessionReader::RegisterEnumValue(
50  "GlobalSysSoln", "DirectFull",
51  MultiRegions::eDirectFullMatrix),
52  LibUtilities::SessionReader::RegisterEnumValue(
53  "GlobalSysSoln", "DirectStaticCond",
54  MultiRegions::eDirectStaticCond),
55  LibUtilities::SessionReader::RegisterEnumValue(
56  "GlobalSysSoln", "DirectMultiLevelStaticCond",
57  MultiRegions::eDirectMultiLevelStaticCond),
58  LibUtilities::SessionReader::RegisterEnumValue(
59  "GlobalSysSoln", "IterativeFull",
60  MultiRegions::eIterativeFull),
61  LibUtilities::SessionReader::RegisterEnumValue(
62  "GlobalSysSoln", "IterativeStaticCond",
63  MultiRegions::eIterativeStaticCond),
64  LibUtilities::SessionReader::RegisterEnumValue(
65  "GlobalSysSoln", "IterativeMultiLevelStaticCond",
66  MultiRegions::eIterativeMultiLevelStaticCond),
67  LibUtilities::SessionReader::RegisterEnumValue(
68  "GlobalSysSoln", "XxtFull",
69  MultiRegions::eXxtFullMatrix),
70  LibUtilities::SessionReader::RegisterEnumValue(
71  "GlobalSysSoln", "XxtStaticCond",
72  MultiRegions::eXxtStaticCond),
73  LibUtilities::SessionReader::RegisterEnumValue(
74  "GlobalSysSoln", "XxtMultiLevelStaticCond",
75  MultiRegions::eXxtMultiLevelStaticCond),
76  LibUtilities::SessionReader::RegisterEnumValue(
77  "GlobalSysSoln", "PETScFull",
78  MultiRegions::ePETScFullMatrix),
79  LibUtilities::SessionReader::RegisterEnumValue(
80  "GlobalSysSoln", "PETScStaticCond",
81  MultiRegions::ePETScStaticCond),
82  LibUtilities::SessionReader::RegisterEnumValue(
83  "GlobalSysSoln", "PETScMultiLevelStaticCond",
84  MultiRegions::ePETScMultiLevelStaticCond)
85  };
86 
87  std::string GlobalLinSys::def = LibUtilities::SessionReader::
88  RegisterDefaultSolverInfo("GlobalSysSoln",
89  "DirectMultiLevelStaticCond");
90 
91  /**
92  * @class GlobalLinSys
93  *
94  * Consider the linear system
95  * \f$\boldsymbol{M\hat{u}}_g=\boldsymbol{\hat{f}}\f$.
96  * Distinguishing between the boundary and interior components of
97  * \f$\boldsymbol{\hat{u}}_g\f$ and \f$\boldsymbol{\hat{f}}\f$ using
98  * \f$\boldsymbol{\hat{u}}_b\f$,\f$\boldsymbol{\hat{u}}_i\f$ and
99  * \f$\boldsymbol{\hat{f}}_b\f$,\f$\boldsymbol{\hat{f}}_i\f$
100  * respectively, this system can be split into its constituent parts as
101  * \f[\left[\begin{array}{cc}
102  * \boldsymbol{M}_b&\boldsymbol{M}_{c1}\\
103  * \boldsymbol{M}_{c2}&\boldsymbol{M}_i\\
104  * \end{array}\right]
105  * \left[\begin{array}{c}
106  * \boldsymbol{\hat{u}_b}\\
107  * \boldsymbol{\hat{u}_i}\\
108  * \end{array}\right]=
109  * \left[\begin{array}{c}
110  * \boldsymbol{\hat{f}_b}\\
111  * \boldsymbol{\hat{f}_i}\\
112  * \end{array}\right]\f]
113  * where \f$\boldsymbol{M}_b\f$ represents the components of
114  * \f$\boldsymbol{M}\f$ resulting from boundary-boundary mode
115  * interactions,
116  * \f$\boldsymbol{M}_{c1}\f$ and \f$\boldsymbol{M}_{c2}\f$ represent the
117  * components resulting from coupling between the boundary-interior
118  * modes, and \f$\boldsymbol{M}_i\f$ represents the components of
119  * \f$\boldsymbol{M}\f$ resulting from interior-interior mode
120  * interactions.
121  *
122  * The solution of the linear system can now be determined in two steps:
123  * \f{eqnarray*}
124  * \mathrm{step 1:}&\quad&(\boldsymbol{M}_b-\boldsymbol{M}_{c1}
125  * \boldsymbol{M}_i^{-1}\boldsymbol{M}_{c2}) \boldsymbol{\hat{u}_b} =
126  * \boldsymbol{\hat{f}}_b - \boldsymbol{M}_{c1}\boldsymbol{M}_i^{-1}
127  * \boldsymbol{\hat{f}}_i,\nonumber \\
128  * \mathrm{step 2:}&\quad&\boldsymbol{\hat{u}_i}=\boldsymbol{M}_i^{-1}
129  * \left( \boldsymbol{\hat{f}}_i
130  * - \boldsymbol{M}_{c2}\boldsymbol{\hat{u}_b}
131  * \right). \nonumber \\ \f}
132  * As the inverse of \f$\boldsymbol{M}_i^{-1}\f$ is
133  * \f[ \boldsymbol{M}_i^{-1} = \left [\underline{\boldsymbol{M}^e_i}
134  * \right ]^{-1} = \underline{[\boldsymbol{M}^e_i]}^{-1} \f]
135  * and the following operations can be evaluated as,
136  * \f{eqnarray*}
137  * \boldsymbol{M}_{c1}\boldsymbol{M}_i^{-1}\boldsymbol{\hat{f}}_i &
138  * =& \boldsymbol{\mathcal{A}}_b^T \underline{\boldsymbol{M}^e_{c1}}
139  * \underline{[\boldsymbol{M}^e_i]}^{-1} \boldsymbol{\hat{f}}_i \\
140  * \boldsymbol{M}_{c2} \boldsymbol{\hat{u}_b} &=&
141  * \underline{\boldsymbol{M}^e_{c2}} \boldsymbol{\mathcal{A}}_b
142  * \boldsymbol{\hat{u}_b}.\f}
143  * where \f$\boldsymbol{\mathcal{A}}_b \f$ is the permutation matrix
144  * which scatters from global to local degrees of freedom, only the
145  * following four matrices should be constructed:
146  * - \f$\underline{[\boldsymbol{M}^e_i]}^{-1}\f$
147  * - \f$\underline{\boldsymbol{M}^e_{c1}}
148  * \underline{[\boldsymbol{M}^e_i]}^{-1}\f$
149  * - \f$\underline{\boldsymbol{M}^e_{c2}}\f$
150  * - The Schur complement: \f$\boldsymbol{M}_{\mathrm{Schur}}=
151  * \quad\boldsymbol{M}_b-\boldsymbol{M}_{c1}\boldsymbol{M}_i^{-1}
152  * \boldsymbol{M}_{c2}\f$
153  *
154  * The first three matrices are just a concatenation of the
155  * corresponding local matrices and they can be created as such. They
156  * also allow for an elemental evaluation of the operations concerned.
157  *
158  * The global Schur complement however should be assembled from the
159  * concatenation of the local elemental Schur complements, that is,
160  * \f[ \boldsymbol{M}_{\mathrm{Schur}}=\boldsymbol{M}_b
161  * - \boldsymbol{M}_{c1}
162  * \boldsymbol{M}_i^{-1} \boldsymbol{M}_{c2} =
163  * \boldsymbol{\mathcal{A}}_b^T \left [\underline{\boldsymbol{M}^e_b -
164  * \boldsymbol{M}^e_{c1} [\boldsymbol{M}^e_i]^{-1}
165  * (\boldsymbol{M}^e_{c2})} \right ] \boldsymbol{\mathcal{A}}_b \f]
166  * and it is the only matrix operation that need to be evaluated on a
167  * global level when using static condensation.
168  * However, due to the size and sparsity of the matrix
169  * \f$\boldsymbol{\mathcal{A}}_b\f$, it is more efficient to assemble
170  * the global Schur matrix using the mapping array bmap\f$[e][i]\f$
171  * contained in the input argument \a locToGloMap. The global Schur
172  * complement is then constructed as:
173  * \f[\boldsymbol{M}_{\mathrm{Schur}}\left[\mathrm{\a bmap}[e][i]\right]
174  * \left[\mathrm{\a bmap}[e][j]\right]=\mathrm{\a bsign}[e][i]\cdot
175  * \mathrm{\a bsign}[e][j]
176  * \cdot\boldsymbol{M}^e_{\mathrm{Schur}}[i][j]\f]
177  * All four matrices are stored in the \a GlobalLinSys returned by this
178  * function.
179  */
180 
181  /**
182  * Given a block matrix, construct a global matrix system according to
183  * a local to global mapping. #m_linSys is constructed by
184  * AssembleFullMatrix().
185  * @param pkey Associated linear system key.
186  * @param locToGloMap Local to global mapping.
187  */
188  GlobalLinSys::GlobalLinSys(const GlobalLinSysKey &pKey,
189  const boost::weak_ptr<ExpList> &pExpList,
190  const boost::shared_ptr<AssemblyMap>
191  &pLocToGloMap):
192  m_linSysKey(pKey),
193  m_expList(pExpList),
194  m_robinBCInfo(m_expList.lock()->GetRobinBCInfo())
195  {
196  }
197 
198  /**
199  *
200  */
201  GlobalLinSysFactory& GetGlobalLinSysFactory()
202  {
203  typedef Loki::SingletonHolder<GlobalLinSysFactory,
204  Loki::CreateUsingNew,
205  Loki::NoDestroy > Type;
206  return Type::Instance();
207  }
208 
209  /**
210  * @brief Get the number of blocks in this system.
211  *
212  * At the top level this corresponds to the number of elements in the
213  * expansion list.
214  */
215  int GlobalLinSys::v_GetNumBlocks()
216  {
217  return m_expList.lock()->GetExpSize();
218  }
219 
220  /**
221  * @brief Retrieves the block matrix from n-th expansion using the
222  * matrix key provided by the #m_linSysKey.
223  *
224  * @param n Number of the expansion.
225  * @return Block matrix for the specified expansion.
226  */
227  DNekScalMatSharedPtr GlobalLinSys::v_GetBlock(unsigned int n)
228  {
229  boost::shared_ptr<MultiRegions::ExpList> expList = m_expList.lock();
230  int cnt = 0;
231  DNekScalMatSharedPtr loc_mat;
232 
233  LocalRegions::ExpansionSharedPtr vExp =
234  boost::dynamic_pointer_cast<LocalRegions::Expansion>(
235  expList->GetExp(n));
236 
237  // need to be initialised with zero size for non variable
238  // coefficient case
239  StdRegions::VarCoeffMap vVarCoeffMap;
240 
241  // retrieve variable coefficients
242  if(m_linSysKey.GetNVarCoeffs() > 0)
243  {
244  StdRegions::VarCoeffMap::const_iterator x;
245  cnt = expList->GetPhys_Offset(n);
246 
247  for (x = m_linSysKey.GetVarCoeffs().begin();
248  x != m_linSysKey.GetVarCoeffs().end(); ++x)
249  {
250  vVarCoeffMap[x->first] = x->second + cnt;
251  }
252  }
253 
254  LocalRegions::MatrixKey matkey(m_linSysKey.GetMatrixType(),
255  vExp->DetShapeType(),
256  *vExp, m_linSysKey.GetConstFactors(),
257  vVarCoeffMap);
258  loc_mat = vExp->GetLocMatrix(matkey);
259 
260  // apply robin boundary conditions to the matrix.
261  if(m_robinBCInfo.count(n) != 0) // add robin mass matrix
262  {
263  RobinBCInfoSharedPtr rBC;
264 
265  // declare local matrix from scaled matrix.
266  int rows = loc_mat->GetRows();
267  int cols = loc_mat->GetColumns();
268  const NekDouble *dat = loc_mat->GetRawPtr();
269  DNekMatSharedPtr new_mat = MemoryManager<DNekMat>::
270  AllocateSharedPtr(rows,cols,dat);
271  Blas::Dscal(rows*cols,loc_mat->Scale(),new_mat->GetRawPtr(),1);
272 
273  // add local matrix contribution
274  for(rBC = m_robinBCInfo.find(n)->second;rBC; rBC = rBC->next)
275  {
276  vExp->AddRobinMassMatrix(
277  rBC->m_robinID, rBC->m_robinPrimitiveCoeffs, new_mat);
278  }
279 
280  // redeclare loc_mat to point to new_mat plus the scalar.
281  loc_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(
282  1.0, new_mat);
283  }
284 
285  // finally return the matrix.
286  return loc_mat;
287  }
288 
289  /**
290  * @brief Retrieves a the static condensation block matrices from n-th
291  * expansion using the matrix key provided by the #m_linSysKey.
292  *
293  * @param n Number of the expansion
294  * @return 2x2 Block matrix holding the static condensation
295  * matrices for the n-th expansion.
296  */
297  DNekScalBlkMatSharedPtr GlobalLinSys::v_GetStaticCondBlock(
298  unsigned int n)
299  {
300  boost::shared_ptr<MultiRegions::ExpList> expList = m_expList.lock();
301  int cnt = 0;
302  DNekScalBlkMatSharedPtr loc_mat;
303  DNekScalMatSharedPtr tmp_mat;
304 
305  StdRegions::StdExpansionSharedPtr vExp = expList->GetExp(n);
306 
307  // need to be initialised with zero size for non variable
308  // coefficient case
309  StdRegions::VarCoeffMap vVarCoeffMap;
310 
311  // retrieve variable coefficients
312  if(m_linSysKey.GetNVarCoeffs() > 0)
313  {
314  StdRegions::VarCoeffMap::const_iterator x;
315  cnt = expList->GetPhys_Offset(n);
316  for (x = m_linSysKey.GetVarCoeffs().begin();
317  x != m_linSysKey.GetVarCoeffs().end (); ++x)
318  {
319  vVarCoeffMap[x->first] = x->second + cnt;
320  }
321  }
322 
323  LocalRegions::MatrixKey matkey(m_linSysKey.GetMatrixType(),
324  vExp->DetShapeType(),
325  *vExp,
326  m_linSysKey.GetConstFactors(),
327  vVarCoeffMap);
328 
329  loc_mat = vExp->GetLocStaticCondMatrix(matkey);
330 
331  if(m_robinBCInfo.count(n) != 0) // add robin mass matrix
332  {
333  RobinBCInfoSharedPtr rBC;
334 
335  tmp_mat = loc_mat->GetBlock(0,0);
336 
337  // declare local matrix from scaled matrix.
338  int rows = tmp_mat->GetRows();
339  int cols = tmp_mat->GetColumns();
340  const NekDouble *dat = tmp_mat->GetRawPtr();
341  DNekMatSharedPtr new_mat = MemoryManager<DNekMat>::
342  AllocateSharedPtr(rows, cols, dat);
343  Blas::Dscal(rows*cols,tmp_mat->Scale(),new_mat->GetRawPtr(),1);
344 
345  // add local matrix contribution
346  for(rBC = m_robinBCInfo.find(n)->second;rBC; rBC = rBC->next)
347  {
348  vExp->AddRobinMassMatrix(
349  rBC->m_robinID, rBC->m_robinPrimitiveCoeffs, new_mat);
350  }
351 
352  // redeclare loc_mat to point to new_mat plus the scalar.
353  tmp_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(
354  1.0, new_mat);
355  loc_mat->SetBlock(0,0,tmp_mat);
356  }
357 
358  return loc_mat;
359  }
360 
361  /**
362  * @brief Releases the static condensation block matrices from NekManager
363  * of n-th expansion using the matrix key provided by the #m_linSysKey.
364  *
365  * @param n Number of the expansion
366  */
367  void GlobalLinSys::v_DropStaticCondBlock(unsigned int n)
368  {
369  boost::shared_ptr<MultiRegions::ExpList> expList = m_expList.lock();
370 
371  StdRegions::StdExpansionSharedPtr vExp = expList->GetExp(n);
372 
373  // need to be initialised with zero size for non variable
374  // coefficient case
375  StdRegions::VarCoeffMap vVarCoeffMap;
376 
377  // retrieve variable coefficients
378  if(m_linSysKey.GetNVarCoeffs() > 0)
379  {
380  StdRegions::VarCoeffMap::const_iterator x;
381  int cnt = expList->GetPhys_Offset(n);
382  for (x = m_linSysKey.GetVarCoeffs().begin();
383  x != m_linSysKey.GetVarCoeffs().end (); ++x)
384  {
385  vVarCoeffMap[x->first] = x->second + cnt;
386  }
387  }
388 
389  LocalRegions::MatrixKey matkey(m_linSysKey.GetMatrixType(),
390  vExp->DetShapeType(),
391  *vExp,
392  m_linSysKey.GetConstFactors(),
393  vVarCoeffMap);
394 
395  vExp->DropLocStaticCondMatrix(matkey);
396  }
397 
398  void GlobalLinSys::v_InitObject()
399  {
400  NEKERROR(ErrorUtil::efatal, "Method does not exist" );
401  }
402 
403  void GlobalLinSys::v_Initialise(
404  const boost::shared_ptr<AssemblyMap>& pLocToGloMap)
405  {
406  NEKERROR(ErrorUtil::efatal, "Method does not exist" );
407  }
408  } //end of namespace
409 } //end of namespace
410 
Generated on Sun Mar 15 2015 00:11:32 for Nektar++ by   doxygen 1.8.1.2