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 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 //
20 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description: GlobalLinSys definition
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 #include <boost/core/ignore_unused.hpp>
36 
38 #include <LocalRegions/Expansion.h>
39 #include <LocalRegions/MatrixKey.h>
42 
45 
46 namespace Nektar
47 {
48 namespace MultiRegions
49 {
50 std::string GlobalLinSys::lookupIds[12] = {
52  "GlobalSysSoln", "DirectFull", MultiRegions::eDirectFullMatrix),
54  "GlobalSysSoln", "DirectStaticCond", MultiRegions::eDirectStaticCond),
56  "GlobalSysSoln", "DirectMultiLevelStaticCond",
59  "GlobalSysSoln", "IterativeFull", MultiRegions::eIterativeFull),
61  "GlobalSysSoln", "IterativeStaticCond",
64  "GlobalSysSoln", "IterativeMultiLevelStaticCond",
67  "GlobalSysSoln", "XxtFull", MultiRegions::eXxtFullMatrix),
69  "GlobalSysSoln", "XxtStaticCond", MultiRegions::eXxtStaticCond),
71  "GlobalSysSoln", "XxtMultiLevelStaticCond",
74  "GlobalSysSoln", "PETScFull", MultiRegions::ePETScFullMatrix),
76  "GlobalSysSoln", "PETScStaticCond", MultiRegions::ePETScStaticCond),
78  "GlobalSysSoln", "PETScMultiLevelStaticCond",
80 
81 #ifdef NEKTAR_USE_SCOTCH
82 std::string GlobalLinSys::def =
84  "GlobalSysSoln", "DirectMultiLevelStaticCond");
85 #else
86 std::string GlobalLinSys::def =
88  "DirectStaticCond");
89 #endif
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 std::weak_ptr<ExpList> &pExpList,
190  const std::shared_ptr<AssemblyMap> &pLocToGloMap)
191  : m_linSysKey(pKey), m_expList(pExpList),
192  m_robinBCInfo(m_expList.lock()->GetRobinBCInfo()),
193  m_verbose(
194  m_expList.lock()->GetSession()->DefinesCmdLineArgument("verbose"))
195 {
196  boost::ignore_unused(pLocToGloMap);
197 }
198 
199 /**
200  *
201  */
203 {
204  static GlobalLinSysFactory instance;
205  return instance;
206 }
207 
208 /**
209  * @brief Create a preconditioner object from the parameters defined in
210  * the supplied assembly map.
211  *
212  * @param asmMap Assembly map used to construct the global system.
213  */
215 {
216  PreconditionerType pType = asmMap->GetPreconType();
217  std::string PreconType = MultiRegions::PreconditionerTypeMap[pType];
218  return GetPreconFactory().CreateInstance(PreconType, GetSharedThisPtr(),
219  asmMap);
220 }
221 
222 /**
223  * @brief Get the number of blocks in this system.
224  *
225  * At the top level this corresponds to the number of elements in the
226  * expansion list.
227  */
229 {
230  return m_expList.lock()->GetExpSize();
231 }
232 
233 /**
234  Assemble the matrix key for each block n
235 **/
236 
238 {
239 
240  std::shared_ptr<MultiRegions::ExpList> expList = m_expList.lock();
241  int cnt = 0;
242 
243  LocalRegions::ExpansionSharedPtr vExp = expList->GetExp(n);
244 
245  // need to be initialised with zero size for non variable
246  // coefficient case
247  StdRegions::VarCoeffMap vVarCoeffMap;
248 
250 
251  // setup variable factors
252  if (m_linSysKey.GetNVarFactors() > 0)
253  {
254  if (m_linSysKey.GetVarFactors().count(
256  {
257  vConstFactorMap[StdRegions::eFactorSVVDiffCoeff] =
259 
262  "VarCoeffSVVCuroffRatio is set but "
263  " not FactorSVVCutoffRatio");
264 
265  vConstFactorMap[StdRegions::eFactorSVVCutoffRatio] =
267  }
268 
269  if (m_linSysKey.GetVarFactors().count(
271  {
272  vConstFactorMap[StdRegions::eFactorSVVPowerKerDiffCoeff] =
275  }
276 
277  if (m_linSysKey.GetVarFactors().count(
279  {
280  vConstFactorMap[StdRegions::eFactorSVVDGKerDiffCoeff] =
283  }
284  }
285 
286  // retrieve variable coefficients
287  if (m_linSysKey.GetNVarCoeffs() > 0)
288  {
289  cnt = expList->GetPhys_Offset(n);
290 
291  for (auto &x : m_linSysKey.GetVarCoeffs())
292  {
293  vVarCoeffMap[x.first] = x.second + cnt;
294  }
295  }
296 
298  vExp->DetShapeType(), *vExp, vConstFactorMap,
299  vVarCoeffMap);
300  return matkey;
301 }
302 
303 /**
304  * @brief Retrieves the block matrix from n-th expansion using the
305  * matrix key provided by the #m_linSysKey.
306  *
307  * @param n Number of the expansion.
308  * @return Block matrix for the specified expansion.
309  */
311 {
312  LocalRegions::ExpansionSharedPtr vExp = m_expList.lock()->GetExp(n);
313  DNekScalMatSharedPtr loc_mat;
314  loc_mat = vExp->GetLocMatrix(GetBlockMatrixKey(n));
315 
316  // apply robin boundary conditions to the matrix.
317  if (m_robinBCInfo.count(n) != 0) // add robin mass matrix
318  {
320 
321  // declare local matrix from scaled matrix.
322  int rows = loc_mat->GetRows();
323  int cols = loc_mat->GetColumns();
324  const NekDouble *dat = loc_mat->GetRawPtr();
325  DNekMatSharedPtr new_mat =
327  Blas::Dscal(rows * cols, loc_mat->Scale(), new_mat->GetRawPtr(), 1);
328 
329  // add local matrix contribution
330  for (rBC = m_robinBCInfo.find(n)->second; rBC; rBC = rBC->next)
331  {
332  vExp->AddRobinMassMatrix(rBC->m_robinID,
333  rBC->m_robinPrimitiveCoeffs, new_mat);
334  }
335 
336  // redeclare loc_mat to point to new_mat plus the scalar.
337  loc_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0, new_mat);
338  }
339 
340  // finally return the matrix.
341  return loc_mat;
342 }
343 
344 /**
345  * @brief Retrieves a the static condensation block matrices from n-th
346  * expansion using the matrix key provided by the #m_linSysKey.
347  *
348  * @param n Number of the expansion
349  * @return 2x2 Block matrix holding the static condensation
350  * matrices for the n-th expansion.
351  */
353 {
354 
355  LocalRegions::ExpansionSharedPtr vExp = m_expList.lock()->GetExp(n);
356  DNekScalBlkMatSharedPtr loc_mat;
357  loc_mat = vExp->GetLocStaticCondMatrix(GetBlockMatrixKey(n));
358 
359  if (m_robinBCInfo.count(n) != 0) // add robin mass matrix
360  {
361  DNekScalMatSharedPtr tmp_mat;
363 
364  tmp_mat = loc_mat->GetBlock(0, 0);
365 
366  // declare local matrix from scaled matrix.
367  int rows = tmp_mat->GetRows();
368  int cols = tmp_mat->GetColumns();
369  const NekDouble *dat = tmp_mat->GetRawPtr();
370  DNekMatSharedPtr new_mat =
372  Blas::Dscal(rows * cols, tmp_mat->Scale(), new_mat->GetRawPtr(), 1);
373 
374  // add local matrix contribution
375  for (rBC = m_robinBCInfo.find(n)->second; rBC; rBC = rBC->next)
376  {
377  vExp->AddRobinMassMatrix(rBC->m_robinID,
378  rBC->m_robinPrimitiveCoeffs, new_mat);
379  }
380 
381  // redeclare loc_mat to point to new_mat plus the scalar.
382  tmp_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0, new_mat);
383  DNekScalBlkMatSharedPtr new_loc_mat;
384  unsigned int exp_size[] = {tmp_mat->GetRows(),
385  loc_mat->GetBlock(1, 1)->GetRows()};
386  unsigned int nblks = 2;
388  nblks, nblks, exp_size, exp_size);
389 
390  new_loc_mat->SetBlock(0, 0, tmp_mat);
391  new_loc_mat->SetBlock(0, 1, loc_mat->GetBlock(0, 1));
392  new_loc_mat->SetBlock(1, 0, loc_mat->GetBlock(1, 0));
393  new_loc_mat->SetBlock(1, 1, loc_mat->GetBlock(1, 1));
394  loc_mat = new_loc_mat;
395  }
396 
397  return loc_mat;
398 }
399 
400 /**
401  * @brief Releases the static condensation block matrices from NekManager
402  * of n-th expansion using the matrix key provided by the #m_linSysKey.
403  *
404  * @param n Number of the expansion
405  */
407 {
408  LocalRegions::ExpansionSharedPtr vExp = m_expList.lock()->GetExp(n);
409  vExp->DropLocStaticCondMatrix(GetBlockMatrixKey(n));
410 }
411 
413 {
414  NEKERROR(ErrorUtil::efatal, "Method does not exist");
415 }
416 
418  const std::shared_ptr<AssemblyMap> &pLocToGloMap)
419 {
420  boost::ignore_unused(pLocToGloMap);
421  NEKERROR(ErrorUtil::efatal, "Method does not exist");
422 }
423 } // namespace MultiRegions
424 } // namespace Nektar
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition: ErrorUtil.hpp:209
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:249
Provides a generic Factory class.
Definition: NekFactory.hpp:105
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:144
static std::string RegisterEnumValue(std::string pEnum, std::string pString, int pEnumValue)
Registers an enumeration value.
static std::string RegisterDefaultSolverInfo(const std::string &pName, const std::string &pValue)
Registers the default string value of a solver info property.
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
virtual int v_GetNumBlocks()
Get the number of blocks in this system.
const std::weak_ptr< ExpList > m_expList
Local Matrix System.
Definition: GlobalLinSys.h:123
static std::string lookupIds[]
Definition: GlobalLinSys.h:156
std::shared_ptr< GlobalLinSys > GetSharedThisPtr()
Returns a shared pointer to the current object.
Definition: GlobalLinSys.h:102
LocalRegions::MatrixKey GetBlockMatrixKey(unsigned int n)
const std::map< int, RobinBCInfoSharedPtr > m_robinBCInfo
Robin boundary info.
Definition: GlobalLinSys.h:125
GlobalLinSys(const GlobalLinSysKey &pKey, const std::weak_ptr< ExpList > &pExpList, const std::shared_ptr< AssemblyMap > &pLocToGloMap)
Constructor for full direct matrix solve.
const GlobalLinSysKey m_linSysKey
Key associated with this linear system.
Definition: GlobalLinSys.h:121
virtual DNekScalMatSharedPtr v_GetBlock(unsigned int n)
Retrieves the block matrix from n-th expansion using the matrix key provided by the m_linSysKey.
PreconditionerSharedPtr CreatePrecon(AssemblyMapSharedPtr asmMap)
Create a preconditioner object from the parameters defined in the supplied assembly map.
virtual DNekScalBlkMatSharedPtr v_GetStaticCondBlock(unsigned int n)
Retrieves a the static condensation block matrices from n-th expansion using the matrix key provided ...
virtual void v_Initialise(const std::shared_ptr< AssemblyMap > &pLocToGloMap)
virtual void v_DropStaticCondBlock(unsigned int n)
Releases the static condensation block matrices from NekManager of n-th expansion using the matrix ke...
const Array< OneD, const NekDouble > & GetVarFactors(const StdRegions::ConstFactorType &coeff) const
const StdRegions::ConstFactorMap & GetConstFactors() const
Returns all the constants.
const StdRegions::VarCoeffMap & GetVarCoeffs() const
StdRegions::MatrixType GetMatrixType() const
Return the matrix type.
static void Dscal(const int &n, const double &alpha, double *x, const int &incx)
BLAS level 1: x = alpha x.
Definition: Blas.hpp:168
std::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
const char *const PreconditionerTypeMap[]
std::shared_ptr< RobinBCInfo > RobinBCInfoSharedPtr
GlobalLinSysFactory & GetGlobalLinSysFactory()
std::shared_ptr< Preconditioner > PreconditionerSharedPtr
Definition: GlobalLinSys.h:60
PreconFactory & GetPreconFactory()
std::shared_ptr< AssemblyMap > AssemblyMapSharedPtr
Definition: AssemblyMap.h:51
std::map< StdRegions::VarCoeffType, Array< OneD, NekDouble > > VarCoeffMap
Definition: StdRegions.hpp:240
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:282
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
std::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:79
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:75
double NekDouble