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 
38 #include <LocalRegions/MatrixKey.h>
39 #include <LocalRegions/Expansion.h>
41 
44 
45 namespace Nektar
46 {
47  namespace MultiRegions
48  {
49  std::string GlobalLinSys::lookupIds[12] = {
51  "GlobalSysSoln", "DirectFull",
54  "GlobalSysSoln", "DirectStaticCond",
57  "GlobalSysSoln", "DirectMultiLevelStaticCond",
60  "GlobalSysSoln", "IterativeFull",
63  "GlobalSysSoln", "IterativeStaticCond",
66  "GlobalSysSoln", "IterativeMultiLevelStaticCond",
69  "GlobalSysSoln", "XxtFull",
72  "GlobalSysSoln", "XxtStaticCond",
75  "GlobalSysSoln", "XxtMultiLevelStaticCond",
78  "GlobalSysSoln", "PETScFull",
81  "GlobalSysSoln", "PETScStaticCond",
84  "GlobalSysSoln", "PETScMultiLevelStaticCond",
86  };
87 
89  RegisterDefaultSolverInfo("GlobalSysSoln",
90  "DirectMultiLevelStaticCond");
91 
92  /**
93  * @class GlobalLinSys
94  *
95  * Consider the linear system
96  * \f$\boldsymbol{M\hat{u}}_g=\boldsymbol{\hat{f}}\f$.
97  * Distinguishing between the boundary and interior components of
98  * \f$\boldsymbol{\hat{u}}_g\f$ and \f$\boldsymbol{\hat{f}}\f$ using
99  * \f$\boldsymbol{\hat{u}}_b\f$,\f$\boldsymbol{\hat{u}}_i\f$ and
100  * \f$\boldsymbol{\hat{f}}_b\f$,\f$\boldsymbol{\hat{f}}_i\f$
101  * respectively, this system can be split into its constituent parts as
102  * \f[\left[\begin{array}{cc}
103  * \boldsymbol{M}_b&\boldsymbol{M}_{c1}\\
104  * \boldsymbol{M}_{c2}&\boldsymbol{M}_i\\
105  * \end{array}\right]
106  * \left[\begin{array}{c}
107  * \boldsymbol{\hat{u}_b}\\
108  * \boldsymbol{\hat{u}_i}\\
109  * \end{array}\right]=
110  * \left[\begin{array}{c}
111  * \boldsymbol{\hat{f}_b}\\
112  * \boldsymbol{\hat{f}_i}\\
113  * \end{array}\right]\f]
114  * where \f$\boldsymbol{M}_b\f$ represents the components of
115  * \f$\boldsymbol{M}\f$ resulting from boundary-boundary mode
116  * interactions,
117  * \f$\boldsymbol{M}_{c1}\f$ and \f$\boldsymbol{M}_{c2}\f$ represent the
118  * components resulting from coupling between the boundary-interior
119  * modes, and \f$\boldsymbol{M}_i\f$ represents the components of
120  * \f$\boldsymbol{M}\f$ resulting from interior-interior mode
121  * interactions.
122  *
123  * The solution of the linear system can now be determined in two steps:
124  * \f{eqnarray*}
125  * \mathrm{step 1:}&\quad&(\boldsymbol{M}_b-\boldsymbol{M}_{c1}
126  * \boldsymbol{M}_i^{-1}\boldsymbol{M}_{c2}) \boldsymbol{\hat{u}_b} =
127  * \boldsymbol{\hat{f}}_b - \boldsymbol{M}_{c1}\boldsymbol{M}_i^{-1}
128  * \boldsymbol{\hat{f}}_i,\nonumber \\
129  * \mathrm{step 2:}&\quad&\boldsymbol{\hat{u}_i}=\boldsymbol{M}_i^{-1}
130  * \left( \boldsymbol{\hat{f}}_i
131  * - \boldsymbol{M}_{c2}\boldsymbol{\hat{u}_b}
132  * \right). \nonumber \\ \f}
133  * As the inverse of \f$\boldsymbol{M}_i^{-1}\f$ is
134  * \f[ \boldsymbol{M}_i^{-1} = \left [\underline{\boldsymbol{M}^e_i}
135  * \right ]^{-1} = \underline{[\boldsymbol{M}^e_i]}^{-1} \f]
136  * and the following operations can be evaluated as,
137  * \f{eqnarray*}
138  * \boldsymbol{M}_{c1}\boldsymbol{M}_i^{-1}\boldsymbol{\hat{f}}_i &
139  * =& \boldsymbol{\mathcal{A}}_b^T \underline{\boldsymbol{M}^e_{c1}}
140  * \underline{[\boldsymbol{M}^e_i]}^{-1} \boldsymbol{\hat{f}}_i \\
141  * \boldsymbol{M}_{c2} \boldsymbol{\hat{u}_b} &=&
142  * \underline{\boldsymbol{M}^e_{c2}} \boldsymbol{\mathcal{A}}_b
143  * \boldsymbol{\hat{u}_b}.\f}
144  * where \f$\boldsymbol{\mathcal{A}}_b \f$ is the permutation matrix
145  * which scatters from global to local degrees of freedom, only the
146  * following four matrices should be constructed:
147  * - \f$\underline{[\boldsymbol{M}^e_i]}^{-1}\f$
148  * - \f$\underline{\boldsymbol{M}^e_{c1}}
149  * \underline{[\boldsymbol{M}^e_i]}^{-1}\f$
150  * - \f$\underline{\boldsymbol{M}^e_{c2}}\f$
151  * - The Schur complement: \f$\boldsymbol{M}_{\mathrm{Schur}}=
152  * \quad\boldsymbol{M}_b-\boldsymbol{M}_{c1}\boldsymbol{M}_i^{-1}
153  * \boldsymbol{M}_{c2}\f$
154  *
155  * The first three matrices are just a concatenation of the
156  * corresponding local matrices and they can be created as such. They
157  * also allow for an elemental evaluation of the operations concerned.
158  *
159  * The global Schur complement however should be assembled from the
160  * concatenation of the local elemental Schur complements, that is,
161  * \f[ \boldsymbol{M}_{\mathrm{Schur}}=\boldsymbol{M}_b
162  * - \boldsymbol{M}_{c1}
163  * \boldsymbol{M}_i^{-1} \boldsymbol{M}_{c2} =
164  * \boldsymbol{\mathcal{A}}_b^T \left [\underline{\boldsymbol{M}^e_b -
165  * \boldsymbol{M}^e_{c1} [\boldsymbol{M}^e_i]^{-1}
166  * (\boldsymbol{M}^e_{c2})} \right ] \boldsymbol{\mathcal{A}}_b \f]
167  * and it is the only matrix operation that need to be evaluated on a
168  * global level when using static condensation.
169  * However, due to the size and sparsity of the matrix
170  * \f$\boldsymbol{\mathcal{A}}_b\f$, it is more efficient to assemble
171  * the global Schur matrix using the mapping array bmap\f$[e][i]\f$
172  * contained in the input argument \a locToGloMap. The global Schur
173  * complement is then constructed as:
174  * \f[\boldsymbol{M}_{\mathrm{Schur}}\left[\mathrm{\a bmap}[e][i]\right]
175  * \left[\mathrm{\a bmap}[e][j]\right]=\mathrm{\a bsign}[e][i]\cdot
176  * \mathrm{\a bsign}[e][j]
177  * \cdot\boldsymbol{M}^e_{\mathrm{Schur}}[i][j]\f]
178  * All four matrices are stored in the \a GlobalLinSys returned by this
179  * function.
180  */
181 
182  /**
183  * Given a block matrix, construct a global matrix system according to
184  * a local to global mapping. #m_linSys is constructed by
185  * AssembleFullMatrix().
186  * @param pkey Associated linear system key.
187  * @param locToGloMap Local to global mapping.
188  */
190  const boost::weak_ptr<ExpList> &pExpList,
191  const boost::shared_ptr<AssemblyMap>
192  &pLocToGloMap):
193  m_linSysKey(pKey),
194  m_expList(pExpList),
195  m_robinBCInfo(m_expList.lock()->GetRobinBCInfo())
196  {
197  }
198 
199  /**
200  *
201  */
203  {
204  typedef Loki::SingletonHolder<GlobalLinSysFactory,
205  Loki::CreateUsingNew,
206  Loki::NoDestroy,
207  Loki::SingleThreaded> Type;
208  return Type::Instance();
209  }
210 
211  /**
212  * @brief Create a preconditioner object from the parameters defined in
213  * the supplied assembly map.
214  *
215  * @param asmMap Assembly map used to construct the global system.
216  */
218  {
219  PreconditionerType pType = asmMap->GetPreconType();
220  std::string PreconType = MultiRegions::PreconditionerTypeMap[pType];
222  PreconType, GetSharedThisPtr(), asmMap);
223  }
224 
225 
226  /**
227  * @brief Get the number of blocks in this system.
228  *
229  * At the top level this corresponds to the number of elements in the
230  * expansion list.
231  */
233  {
234  return m_expList.lock()->GetExpSize();
235  }
236 
237  /**
238  * @brief Retrieves the block matrix from n-th expansion using the
239  * matrix key provided by the #m_linSysKey.
240  *
241  * @param n Number of the expansion.
242  * @return Block matrix for the specified expansion.
243  */
245  {
246  boost::shared_ptr<MultiRegions::ExpList> expList = m_expList.lock();
247  int cnt = 0;
248  DNekScalMatSharedPtr loc_mat;
249 
251  boost::dynamic_pointer_cast<LocalRegions::Expansion>(
252  expList->GetExp(n));
253 
254  // need to be initialised with zero size for non variable
255  // coefficient case
256  StdRegions::VarCoeffMap vVarCoeffMap;
257 
258  // retrieve variable coefficients
259  if(m_linSysKey.GetNVarCoeffs() > 0)
260  {
261  StdRegions::VarCoeffMap::const_iterator x;
262  cnt = expList->GetPhys_Offset(n);
263 
264  for (x = m_linSysKey.GetVarCoeffs().begin();
265  x != m_linSysKey.GetVarCoeffs().end(); ++x)
266  {
267  vVarCoeffMap[x->first] = x->second + cnt;
268  }
269  }
270 
272  vExp->DetShapeType(),
273  *vExp, m_linSysKey.GetConstFactors(),
274  vVarCoeffMap);
275  loc_mat = vExp->GetLocMatrix(matkey);
276 
277  // apply robin boundary conditions to the matrix.
278  if(m_robinBCInfo.count(n) != 0) // add robin mass matrix
279  {
281 
282  // declare local matrix from scaled matrix.
283  int rows = loc_mat->GetRows();
284  int cols = loc_mat->GetColumns();
285  const NekDouble *dat = loc_mat->GetRawPtr();
287  AllocateSharedPtr(rows,cols,dat);
288  Blas::Dscal(rows*cols,loc_mat->Scale(),new_mat->GetRawPtr(),1);
289 
290  // add local matrix contribution
291  for(rBC = m_robinBCInfo.find(n)->second;rBC; rBC = rBC->next)
292  {
293  vExp->AddRobinMassMatrix(
294  rBC->m_robinID, rBC->m_robinPrimitiveCoeffs, new_mat);
295  }
296 
297  // redeclare loc_mat to point to new_mat plus the scalar.
299  1.0, new_mat);
300  }
301 
302  // finally return the matrix.
303  return loc_mat;
304  }
305 
306  /**
307  * @brief Retrieves a the static condensation block matrices from n-th
308  * expansion using the matrix key provided by the #m_linSysKey.
309  *
310  * @param n Number of the expansion
311  * @return 2x2 Block matrix holding the static condensation
312  * matrices for the n-th expansion.
313  */
315  unsigned int n)
316  {
317  boost::shared_ptr<MultiRegions::ExpList> expList = m_expList.lock();
318  int cnt = 0;
319  DNekScalBlkMatSharedPtr loc_mat;
320  DNekScalMatSharedPtr tmp_mat;
321 
322  StdRegions::StdExpansionSharedPtr vExp = expList->GetExp(n);
323 
324  // need to be initialised with zero size for non variable
325  // coefficient case
326  StdRegions::VarCoeffMap vVarCoeffMap;
327 
328  // retrieve variable coefficients
329  if(m_linSysKey.GetNVarCoeffs() > 0)
330  {
331  StdRegions::VarCoeffMap::const_iterator x;
332  cnt = expList->GetPhys_Offset(n);
333  for (x = m_linSysKey.GetVarCoeffs().begin();
334  x != m_linSysKey.GetVarCoeffs().end (); ++x)
335  {
336  vVarCoeffMap[x->first] = x->second + cnt;
337  }
338  }
339 
341  vExp->DetShapeType(),
342  *vExp,
344  vVarCoeffMap);
345 
346  loc_mat = vExp->GetLocStaticCondMatrix(matkey);
347 
348  if(m_robinBCInfo.count(n) != 0) // add robin mass matrix
349  {
351 
352  tmp_mat = loc_mat->GetBlock(0,0);
353 
354  // declare local matrix from scaled matrix.
355  int rows = tmp_mat->GetRows();
356  int cols = tmp_mat->GetColumns();
357  const NekDouble *dat = tmp_mat->GetRawPtr();
359  AllocateSharedPtr(rows, cols, dat);
360  Blas::Dscal(rows*cols,tmp_mat->Scale(),new_mat->GetRawPtr(),1);
361 
362  // add local matrix contribution
363  for(rBC = m_robinBCInfo.find(n)->second;rBC; rBC = rBC->next)
364  {
365  vExp->AddRobinMassMatrix(
366  rBC->m_robinID, rBC->m_robinPrimitiveCoeffs, new_mat);
367  }
368 
369  // redeclare loc_mat to point to new_mat plus the scalar.
371  1.0, new_mat);
372  DNekScalBlkMatSharedPtr new_loc_mat;
373  unsigned int exp_size[] = {tmp_mat->GetRows(), loc_mat->GetBlock(1,1)->GetRows()};
374  unsigned int nblks = 2;
375  new_loc_mat = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(nblks, nblks, exp_size, exp_size);
376 
377 
378  new_loc_mat->SetBlock(0,0,tmp_mat);
379  new_loc_mat->SetBlock(0,1,loc_mat->GetBlock(0,1));
380  new_loc_mat->SetBlock(1,0,loc_mat->GetBlock(1,0));
381  new_loc_mat->SetBlock(1,1,loc_mat->GetBlock(1,1));
382  loc_mat = new_loc_mat;
383  }
384 
385  return loc_mat;
386  }
387 
388  /**
389  * @brief Releases the static condensation block matrices from NekManager
390  * of n-th expansion using the matrix key provided by the #m_linSysKey.
391  *
392  * @param n Number of the expansion
393  */
395  {
396  boost::shared_ptr<MultiRegions::ExpList> expList = m_expList.lock();
397 
398  StdRegions::StdExpansionSharedPtr vExp = expList->GetExp(n);
399 
400  // need to be initialised with zero size for non variable
401  // coefficient case
402  StdRegions::VarCoeffMap vVarCoeffMap;
403 
404  // retrieve variable coefficients
405  if(m_linSysKey.GetNVarCoeffs() > 0)
406  {
407  StdRegions::VarCoeffMap::const_iterator x;
408  int cnt = expList->GetPhys_Offset(n);
409  for (x = m_linSysKey.GetVarCoeffs().begin();
410  x != m_linSysKey.GetVarCoeffs().end (); ++x)
411  {
412  vVarCoeffMap[x->first] = x->second + cnt;
413  }
414  }
415 
417  vExp->DetShapeType(),
418  *vExp,
420  vVarCoeffMap);
421 
422  vExp->DropLocStaticCondMatrix(matkey);
423  }
424 
426  {
427  NEKERROR(ErrorUtil::efatal, "Method does not exist" );
428  }
429 
431  const boost::shared_ptr<AssemblyMap>& pLocToGloMap)
432  {
433  NEKERROR(ErrorUtil::efatal, "Method does not exist" );
434  }
435  } //end of namespace
436 } //end of namespace
437 
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:69
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
tBaseSharedPtr CreateInstance(tKey idKey BOOST_PP_COMMA_IF(MAX_PARAM) BOOST_PP_ENUM_BINARY_PARAMS(MAX_PARAM, tParam, x))
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:162
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< AssemblyMap > AssemblyMapSharedPtr
Definition: AssemblyMap.h:53
boost::shared_ptr< RobinBCInfo > RobinBCInfoSharedPtr
PreconFactory & GetPreconFactory()
boost::shared_ptr< Preconditioner > PreconditionerSharedPtr
Definition: GlobalLinSys.h:62
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:226
virtual DNekScalMatSharedPtr v_GetBlock(unsigned int n)
Retrieves the block matrix from n-th expansion using the matrix key provided by the m_linSysKey...
boost::shared_ptr< GlobalLinSys > GetSharedThisPtr()
Returns a shared pointer to the current object.
Definition: GlobalLinSys.h:107
boost::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:74
virtual int v_GetNumBlocks()
Get the number of blocks in this system.
double NekDouble
boost::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
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:127
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:161
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
const char *const PreconditionerTypeMap[]
const std::map< int, RobinBCInfoSharedPtr > m_robinBCInfo
Robin boundary info.
Definition: GlobalLinSys.h:131
PreconditionerSharedPtr CreatePrecon(AssemblyMapSharedPtr asmMap)
Create a preconditioner object from the parameters defined in the supplied assembly map...
Provides a generic Factory class.
Definition: NekFactory.hpp:116
const boost::weak_ptr< ExpList > m_expList
Local Matrix System.
Definition: GlobalLinSys.h:129