Nektar++
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 
37 #include <LocalRegions/MatrixKey.h>
38 #include <LocalRegions/Expansion.h>
40 
43 
44 namespace Nektar
45 {
46  namespace MultiRegions
47  {
48  std::string GlobalLinSys::lookupIds[12] = {
50  "GlobalSysSoln", "DirectFull",
53  "GlobalSysSoln", "DirectStaticCond",
56  "GlobalSysSoln", "DirectMultiLevelStaticCond",
59  "GlobalSysSoln", "IterativeFull",
62  "GlobalSysSoln", "IterativeStaticCond",
65  "GlobalSysSoln", "IterativeMultiLevelStaticCond",
68  "GlobalSysSoln", "XxtFull",
71  "GlobalSysSoln", "XxtStaticCond",
74  "GlobalSysSoln", "XxtMultiLevelStaticCond",
77  "GlobalSysSoln", "PETScFull",
80  "GlobalSysSoln", "PETScStaticCond",
83  "GlobalSysSoln", "PETScMultiLevelStaticCond",
85  };
86 
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  */
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  */
202  {
203  typedef Loki::SingletonHolder<GlobalLinSysFactory,
204  Loki::CreateUsingNew,
205  Loki::NoDestroy,
206  Loki::SingleThreaded> Type;
207  return Type::Instance();
208  }
209 
210  /**
211  * @brief Get the number of blocks in this system.
212  *
213  * At the top level this corresponds to the number of elements in the
214  * expansion list.
215  */
217  {
218  return m_expList.lock()->GetExpSize();
219  }
220 
221  /**
222  * @brief Retrieves the block matrix from n-th expansion using the
223  * matrix key provided by the #m_linSysKey.
224  *
225  * @param n Number of the expansion.
226  * @return Block matrix for the specified expansion.
227  */
229  {
230  boost::shared_ptr<MultiRegions::ExpList> expList = m_expList.lock();
231  int cnt = 0;
232  DNekScalMatSharedPtr loc_mat;
233 
235  boost::dynamic_pointer_cast<LocalRegions::Expansion>(
236  expList->GetExp(n));
237 
238  // need to be initialised with zero size for non variable
239  // coefficient case
240  StdRegions::VarCoeffMap vVarCoeffMap;
241 
242  // retrieve variable coefficients
243  if(m_linSysKey.GetNVarCoeffs() > 0)
244  {
245  StdRegions::VarCoeffMap::const_iterator x;
246  cnt = expList->GetPhys_Offset(n);
247 
248  for (x = m_linSysKey.GetVarCoeffs().begin();
249  x != m_linSysKey.GetVarCoeffs().end(); ++x)
250  {
251  vVarCoeffMap[x->first] = x->second + cnt;
252  }
253  }
254 
256  vExp->DetShapeType(),
257  *vExp, m_linSysKey.GetConstFactors(),
258  vVarCoeffMap);
259  loc_mat = vExp->GetLocMatrix(matkey);
260 
261  // apply robin boundary conditions to the matrix.
262  if(m_robinBCInfo.count(n) != 0) // add robin mass matrix
263  {
265 
266  // declare local matrix from scaled matrix.
267  int rows = loc_mat->GetRows();
268  int cols = loc_mat->GetColumns();
269  const NekDouble *dat = loc_mat->GetRawPtr();
271  AllocateSharedPtr(rows,cols,dat);
272  Blas::Dscal(rows*cols,loc_mat->Scale(),new_mat->GetRawPtr(),1);
273 
274  // add local matrix contribution
275  for(rBC = m_robinBCInfo.find(n)->second;rBC; rBC = rBC->next)
276  {
277  vExp->AddRobinMassMatrix(
278  rBC->m_robinID, rBC->m_robinPrimitiveCoeffs, new_mat);
279  }
280 
281  // redeclare loc_mat to point to new_mat plus the scalar.
283  1.0, new_mat);
284  }
285 
286  // finally return the matrix.
287  return loc_mat;
288  }
289 
290  /**
291  * @brief Retrieves a the static condensation block matrices from n-th
292  * expansion using the matrix key provided by the #m_linSysKey.
293  *
294  * @param n Number of the expansion
295  * @return 2x2 Block matrix holding the static condensation
296  * matrices for the n-th expansion.
297  */
299  unsigned int n)
300  {
301  boost::shared_ptr<MultiRegions::ExpList> expList = m_expList.lock();
302  int cnt = 0;
303  DNekScalBlkMatSharedPtr loc_mat;
304  DNekScalMatSharedPtr tmp_mat;
305 
306  StdRegions::StdExpansionSharedPtr vExp = expList->GetExp(n);
307 
308  // need to be initialised with zero size for non variable
309  // coefficient case
310  StdRegions::VarCoeffMap vVarCoeffMap;
311 
312  // retrieve variable coefficients
313  if(m_linSysKey.GetNVarCoeffs() > 0)
314  {
315  StdRegions::VarCoeffMap::const_iterator x;
316  cnt = expList->GetPhys_Offset(n);
317  for (x = m_linSysKey.GetVarCoeffs().begin();
318  x != m_linSysKey.GetVarCoeffs().end (); ++x)
319  {
320  vVarCoeffMap[x->first] = x->second + cnt;
321  }
322  }
323 
325  vExp->DetShapeType(),
326  *vExp,
328  vVarCoeffMap);
329 
330  loc_mat = vExp->GetLocStaticCondMatrix(matkey);
331 
332  if(m_robinBCInfo.count(n) != 0) // add robin mass matrix
333  {
335 
336  tmp_mat = loc_mat->GetBlock(0,0);
337 
338  // declare local matrix from scaled matrix.
339  int rows = tmp_mat->GetRows();
340  int cols = tmp_mat->GetColumns();
341  const NekDouble *dat = tmp_mat->GetRawPtr();
343  AllocateSharedPtr(rows, cols, dat);
344  Blas::Dscal(rows*cols,tmp_mat->Scale(),new_mat->GetRawPtr(),1);
345 
346  // add local matrix contribution
347  for(rBC = m_robinBCInfo.find(n)->second;rBC; rBC = rBC->next)
348  {
349  vExp->AddRobinMassMatrix(
350  rBC->m_robinID, rBC->m_robinPrimitiveCoeffs, new_mat);
351  }
352 
353  // redeclare loc_mat to point to new_mat plus the scalar.
355  1.0, new_mat);
356  loc_mat->SetBlock(0,0,tmp_mat);
357  }
358 
359  return loc_mat;
360  }
361 
362  /**
363  * @brief Releases the static condensation block matrices from NekManager
364  * of n-th expansion using the matrix key provided by the #m_linSysKey.
365  *
366  * @param n Number of the expansion
367  */
369  {
370  boost::shared_ptr<MultiRegions::ExpList> expList = m_expList.lock();
371 
372  StdRegions::StdExpansionSharedPtr vExp = expList->GetExp(n);
373 
374  // need to be initialised with zero size for non variable
375  // coefficient case
376  StdRegions::VarCoeffMap vVarCoeffMap;
377 
378  // retrieve variable coefficients
379  if(m_linSysKey.GetNVarCoeffs() > 0)
380  {
381  StdRegions::VarCoeffMap::const_iterator x;
382  int cnt = expList->GetPhys_Offset(n);
383  for (x = m_linSysKey.GetVarCoeffs().begin();
384  x != m_linSysKey.GetVarCoeffs().end (); ++x)
385  {
386  vVarCoeffMap[x->first] = x->second + cnt;
387  }
388  }
389 
391  vExp->DetShapeType(),
392  *vExp,
394  vVarCoeffMap);
395 
396  vExp->DropLocStaticCondMatrix(matkey);
397  }
398 
400  {
401  NEKERROR(ErrorUtil::efatal, "Method does not exist" );
402  }
403 
405  const boost::shared_ptr<AssemblyMap>& pLocToGloMap)
406  {
407  NEKERROR(ErrorUtil::efatal, "Method does not exist" );
408  }
409  } //end of namespace
410 } //end of namespace
411 
LibUtilities::NekFactory< std::string, GlobalLinSys, const GlobalLinSysKey &, const boost::weak_ptr< ExpList > &, const boost::shared_ptr< AssemblyMap > & > GlobalLinSysFactory
Datatype of the NekFactory used to instantiate classes derived from the EquationSystem class...
Definition: GlobalLinSys.h:65
const StdRegions::VarCoeffMap & GetVarCoeffs() const
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
Definition: ErrorUtil.hpp:158
static std::string RegisterEnumValue(std::string pEnum, std::string pString, int pEnumValue)
Registers an enumeration value.
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
boost::shared_ptr< RobinBCInfo > RobinBCInfoSharedPtr
virtual DNekScalBlkMatSharedPtr v_GetStaticCondBlock(unsigned int n)
Retrieves a the static condensation block matrices from n-th expansion using the matrix key provided ...
const StdRegions::ConstFactorMap & GetConstFactors() const
Returns all the constants.
boost::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:70
boost::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
virtual void v_Initialise(const boost::shared_ptr< AssemblyMap > &pLocToGloMap)
std::map< StdRegions::VarCoeffType, Array< OneD, NekDouble > > VarCoeffMap
Definition: StdRegions.hpp:225
virtual DNekScalMatSharedPtr v_GetBlock(unsigned int n)
Retrieves the block matrix from n-th expansion using the matrix key provided by the m_linSysKey...
const map< int, RobinBCInfoSharedPtr > m_robinBCInfo
Robin boundary info.
Definition: GlobalLinSys.h:127
boost::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:74
virtual int v_GetNumBlocks()
Get the number of blocks in this system.
double NekDouble
GlobalLinSys(const GlobalLinSysKey &pKey, const boost::weak_ptr< ExpList > &pExpList, const boost::shared_ptr< AssemblyMap > &pLocToGloMap)
Constructor for full direct matrix solve.
Describe a linear system.
StdRegions::MatrixType GetMatrixType() const
Return the matrix type.
const GlobalLinSysKey m_linSysKey
Key associated with this linear system.
Definition: GlobalLinSys.h:123
boost::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
virtual void v_DropStaticCondBlock(unsigned int n)
Releases the static condensation block matrices from NekManager of n-th expansion using the matrix ke...
static std::string lookupIds[]
Definition: GlobalLinSys.h:155
static std::string RegisterDefaultSolverInfo(const std::string &pName, const std::string &pValue)
Registers the default string value of a solver info property.
GlobalLinSysFactory & GetGlobalLinSysFactory()
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
Provides a generic Factory class.
Definition: NekFactory.hpp:116
const boost::weak_ptr< ExpList > m_expList
Local Matrix System.
Definition: GlobalLinSys.h:125