Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties 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  m_verbose(m_expList.lock()->GetSession()->
197  DefinesCmdLineArgument("verbose"))
198  {
199  }
200 
201  /**
202  *
203  */
205  {
206  typedef Loki::SingletonHolder<GlobalLinSysFactory,
207  Loki::CreateUsingNew,
208  Loki::NoDestroy,
209  Loki::SingleThreaded> Type;
210  return Type::Instance();
211  }
212 
213  /**
214  * @brief Create a preconditioner object from the parameters defined in
215  * the supplied assembly map.
216  *
217  * @param asmMap Assembly map used to construct the global system.
218  */
220  {
221  PreconditionerType pType = asmMap->GetPreconType();
222  std::string PreconType = MultiRegions::PreconditionerTypeMap[pType];
224  PreconType, GetSharedThisPtr(), asmMap);
225  }
226 
227 
228  /**
229  * @brief Get the number of blocks in this system.
230  *
231  * At the top level this corresponds to the number of elements in the
232  * expansion list.
233  */
235  {
236  return m_expList.lock()->GetExpSize();
237  }
238 
239  /**
240  * @brief Retrieves the block matrix from n-th expansion using the
241  * matrix key provided by the #m_linSysKey.
242  *
243  * @param n Number of the expansion.
244  * @return Block matrix for the specified expansion.
245  */
247  {
248  boost::shared_ptr<MultiRegions::ExpList> expList = m_expList.lock();
249  int cnt = 0;
250  DNekScalMatSharedPtr loc_mat;
251 
253  boost::dynamic_pointer_cast<LocalRegions::Expansion>(
254  expList->GetExp(n));
255 
256  // need to be initialised with zero size for non variable
257  // coefficient case
258  StdRegions::VarCoeffMap vVarCoeffMap;
259 
260  // retrieve variable coefficients
261  if(m_linSysKey.GetNVarCoeffs() > 0)
262  {
263  StdRegions::VarCoeffMap::const_iterator x;
264  cnt = expList->GetPhys_Offset(n);
265 
266  for (x = m_linSysKey.GetVarCoeffs().begin();
267  x != m_linSysKey.GetVarCoeffs().end(); ++x)
268  {
269  vVarCoeffMap[x->first] = x->second + cnt;
270  }
271  }
272 
274  vExp->DetShapeType(),
275  *vExp, m_linSysKey.GetConstFactors(),
276  vVarCoeffMap);
277  loc_mat = vExp->GetLocMatrix(matkey);
278 
279  // apply robin boundary conditions to the matrix.
280  if(m_robinBCInfo.count(n) != 0) // add robin mass matrix
281  {
283 
284  // declare local matrix from scaled matrix.
285  int rows = loc_mat->GetRows();
286  int cols = loc_mat->GetColumns();
287  const NekDouble *dat = loc_mat->GetRawPtr();
289  AllocateSharedPtr(rows,cols,dat);
290  Blas::Dscal(rows*cols,loc_mat->Scale(),new_mat->GetRawPtr(),1);
291 
292  // add local matrix contribution
293  for(rBC = m_robinBCInfo.find(n)->second;rBC; rBC = rBC->next)
294  {
295  vExp->AddRobinMassMatrix(
296  rBC->m_robinID, rBC->m_robinPrimitiveCoeffs, new_mat);
297  }
298 
299  // redeclare loc_mat to point to new_mat plus the scalar.
301  1.0, new_mat);
302  }
303 
304  // finally return the matrix.
305  return loc_mat;
306  }
307 
308  /**
309  * @brief Retrieves a the static condensation block matrices from n-th
310  * expansion using the matrix key provided by the #m_linSysKey.
311  *
312  * @param n Number of the expansion
313  * @return 2x2 Block matrix holding the static condensation
314  * matrices for the n-th expansion.
315  */
317  unsigned int n)
318  {
319  boost::shared_ptr<MultiRegions::ExpList> expList = m_expList.lock();
320  int cnt = 0;
321  DNekScalBlkMatSharedPtr loc_mat;
322  DNekScalMatSharedPtr tmp_mat;
323 
324  StdRegions::StdExpansionSharedPtr vExp = expList->GetExp(n);
325 
326  // need to be initialised with zero size for non variable
327  // coefficient case
328  StdRegions::VarCoeffMap vVarCoeffMap;
329 
330  // retrieve variable coefficients
331  if(m_linSysKey.GetNVarCoeffs() > 0)
332  {
333  StdRegions::VarCoeffMap::const_iterator x;
334  cnt = expList->GetPhys_Offset(n);
335  for (x = m_linSysKey.GetVarCoeffs().begin();
336  x != m_linSysKey.GetVarCoeffs().end (); ++x)
337  {
338  vVarCoeffMap[x->first] = x->second + cnt;
339  }
340  }
341 
343  vExp->DetShapeType(),
344  *vExp,
346  vVarCoeffMap);
347 
348  loc_mat = vExp->GetLocStaticCondMatrix(matkey);
349 
350  if(m_robinBCInfo.count(n) != 0) // add robin mass matrix
351  {
353 
354  tmp_mat = loc_mat->GetBlock(0,0);
355 
356  // declare local matrix from scaled matrix.
357  int rows = tmp_mat->GetRows();
358  int cols = tmp_mat->GetColumns();
359  const NekDouble *dat = tmp_mat->GetRawPtr();
361  AllocateSharedPtr(rows, cols, dat);
362  Blas::Dscal(rows*cols,tmp_mat->Scale(),new_mat->GetRawPtr(),1);
363 
364  // add local matrix contribution
365  for(rBC = m_robinBCInfo.find(n)->second;rBC; rBC = rBC->next)
366  {
367  vExp->AddRobinMassMatrix(
368  rBC->m_robinID, rBC->m_robinPrimitiveCoeffs, new_mat);
369  }
370 
371  // redeclare loc_mat to point to new_mat plus the scalar.
373  1.0, new_mat);
374  DNekScalBlkMatSharedPtr new_loc_mat;
375  unsigned int exp_size[] = {tmp_mat->GetRows(), loc_mat->GetBlock(1,1)->GetRows()};
376  unsigned int nblks = 2;
377  new_loc_mat = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(nblks, nblks, exp_size, exp_size);
378 
379 
380  new_loc_mat->SetBlock(0,0,tmp_mat);
381  new_loc_mat->SetBlock(0,1,loc_mat->GetBlock(0,1));
382  new_loc_mat->SetBlock(1,0,loc_mat->GetBlock(1,0));
383  new_loc_mat->SetBlock(1,1,loc_mat->GetBlock(1,1));
384  loc_mat = new_loc_mat;
385  }
386 
387  return loc_mat;
388  }
389 
390  /**
391  * @brief Releases the static condensation block matrices from NekManager
392  * of n-th expansion using the matrix key provided by the #m_linSysKey.
393  *
394  * @param n Number of the expansion
395  */
397  {
398  boost::shared_ptr<MultiRegions::ExpList> expList = m_expList.lock();
399 
400  StdRegions::StdExpansionSharedPtr vExp = expList->GetExp(n);
401 
402  // need to be initialised with zero size for non variable
403  // coefficient case
404  StdRegions::VarCoeffMap vVarCoeffMap;
405 
406  // retrieve variable coefficients
407  if(m_linSysKey.GetNVarCoeffs() > 0)
408  {
409  StdRegions::VarCoeffMap::const_iterator x;
410  int cnt = expList->GetPhys_Offset(n);
411  for (x = m_linSysKey.GetVarCoeffs().begin();
412  x != m_linSysKey.GetVarCoeffs().end (); ++x)
413  {
414  vVarCoeffMap[x->first] = x->second + cnt;
415  }
416  }
417 
419  vExp->DetShapeType(),
420  *vExp,
422  vVarCoeffMap);
423 
424  vExp->DropLocStaticCondMatrix(matkey);
425  }
426 
428  {
429  NEKERROR(ErrorUtil::efatal, "Method does not exist" );
430  }
431 
433  const boost::shared_ptr<AssemblyMap>& pLocToGloMap)
434  {
435  NEKERROR(ErrorUtil::efatal, "Method does not exist" );
436  }
437  } //end of namespace
438 } //end of namespace
439 
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:191
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:227
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:163
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