Nektar++
GlobalLinSysPETScStaticCond.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 
38 #include <petscsys.h>
39 #include <petscksp.h>
40 #include <petscmat.h>
41 
42 namespace Nektar
43 {
44  namespace MultiRegions
45  {
46  /**
47  * @class GlobalLinSysPETSc
48  *
49  * Solves a linear system using single- or multi-level static
50  * condensation.
51  */
52 
53  /**
54  * Registers the class with the Factory.
55  */
58  "PETScStaticCond",
60  "PETSc static condensation.");
61 
64  "PETScMultiLevelStaticCond",
66  "PETSc multi-level static condensation.");
67 
68  /**
69  * For a matrix system of the form @f[
70  * \left[ \begin{array}{cc}
71  * \boldsymbol{A} & \boldsymbol{B}\\
72  * \boldsymbol{C} & \boldsymbol{D}
73  * \end{array} \right]
74  * \left[ \begin{array}{c} \boldsymbol{x_1}\\ \boldsymbol{x_2}
75  * \end{array}\right]
76  * = \left[ \begin{array}{c} \boldsymbol{y_1}\\ \boldsymbol{y_2}
77  * \end{array}\right],
78  * @f]
79  * where @f$\boldsymbol{D}@f$ and
80  * @f$(\boldsymbol{A-BD^{-1}C})@f$ are invertible, store and assemble
81  * a static condensation system, according to a given local to global
82  * mapping. #m_linSys is constructed by AssembleSchurComplement().
83  * @param mKey Associated matrix key.
84  * @param pLocMatSys LocalMatrixSystem
85  * @param locToGloMap Local to global mapping.
86  */
88  const GlobalLinSysKey &pKey,
89  const boost::weak_ptr<ExpList> &pExpList,
90  const boost::shared_ptr<AssemblyMap>
91  &pLocToGloMap)
92  : GlobalLinSys (pKey, pExpList, pLocToGloMap),
93  GlobalLinSysPETSc (pKey, pExpList, pLocToGloMap),
94  GlobalLinSysStaticCond(pKey, pExpList, pLocToGloMap)
95  {
98  "This constructor is only valid when using static "
99  "condensation");
101  == pLocToGloMap->GetGlobalSysSolnType(),
102  "The local to global map is not set up for the requested "
103  "solution type");
104  }
105 
106  /**
107  *
108  */
110  const GlobalLinSysKey &pKey,
111  const boost::weak_ptr<ExpList> &pExpList,
112  const DNekScalBlkMatSharedPtr pSchurCompl,
113  const DNekScalBlkMatSharedPtr pBinvD,
114  const DNekScalBlkMatSharedPtr pC,
115  const DNekScalBlkMatSharedPtr pInvD,
116  const boost::shared_ptr<AssemblyMap>
117  &pLocToGloMap)
118  : GlobalLinSys (pKey, pExpList, pLocToGloMap),
119  GlobalLinSysPETSc (pKey, pExpList, pLocToGloMap),
120  GlobalLinSysStaticCond(pKey, pExpList, pLocToGloMap)
121  {
122  m_schurCompl = pSchurCompl;
123  m_BinvD = pBinvD;
124  m_C = pC;
125  m_invD = pInvD;
126  }
127 
128  /**
129  *
130  */
132  {
133 
134  }
135 
136  /**
137  * Assemble the schur complement matrix from the block matrices stored
138  * in #m_blkMatrices and the given local to global mapping information.
139  * @param locToGloMap Local to global mapping information.
140  */
142  AssemblyMapSharedPtr pLocToGloMap)
143  {
144  int i, j, n, cnt, gid1, gid2, loc_lda;
145  NekDouble sign1, sign2, value;
146 
147  const int nDirDofs = pLocToGloMap->GetNumGlobalDirBndCoeffs();
148 
153  DNekScalMatSharedPtr loc_mat;
154 
155  // CALCULATE REORDERING MAPPING
156  CalculateReordering(pLocToGloMap->GetGlobalToUniversalBndMap(),
157  pLocToGloMap->GetGlobalToUniversalBndMapUnique(),
158  pLocToGloMap);
159 
160  // SET UP VECTORS AND MATRIX
161  SetUpMatVec();
162 
163  for(n = cnt = 0; n < m_schurCompl->GetNumberOfBlockRows(); ++n)
164  {
165  loc_mat = m_schurCompl->GetBlock(n,n);
166  loc_lda = loc_mat->GetRows();
167 
168  for(i = 0; i < loc_lda; ++i)
169  {
170  gid1 = pLocToGloMap->GetLocalToGlobalBndMap(cnt + i)-nDirDofs;
171  sign1 = pLocToGloMap->GetLocalToGlobalBndSign(cnt + i);
172  if(gid1 >= 0)
173  {
174  int gid1ro = m_reorderedMap[gid1];
175  for(j = 0; j < loc_lda; ++j)
176  {
177  gid2 = pLocToGloMap->GetLocalToGlobalBndMap(cnt + j)
178  - nDirDofs;
179  sign2 = pLocToGloMap->GetLocalToGlobalBndSign(cnt + j);
180  if(gid2 >= 0)
181  {
182  int gid2ro = m_reorderedMap[gid2];
183  value = sign1*sign2*(*loc_mat)(i,j);
184  MatSetValue(
185  m_matrix, gid1ro, gid2ro, value, ADD_VALUES);
186  }
187  }
188  }
189  }
190  cnt += loc_lda;
191  }
192 
193  // ASSEMBLE MATRIX
194  MatAssemblyBegin(m_matrix, MAT_FINAL_ASSEMBLY);
195  MatAssemblyEnd (m_matrix, MAT_FINAL_ASSEMBLY);
196 
197  // SET UP SCATTER OBJECTS
198  SetUpScatter();
199 
200  // CONSTRUCT KSP OBJECT
201  SetUpSolver(pLocToGloMap->GetIterativeTolerance());
202  }
203 
205  const GlobalLinSysKey &mkey,
206  const boost::weak_ptr<ExpList> &pExpList,
207  const DNekScalBlkMatSharedPtr pSchurCompl,
208  const DNekScalBlkMatSharedPtr pBinvD,
209  const DNekScalBlkMatSharedPtr pC,
210  const DNekScalBlkMatSharedPtr pInvD,
211  const boost::shared_ptr<AssemblyMap> &l2gMap)
212  {
214  GlobalLinSysPETScStaticCond>::AllocateSharedPtr(
215  mkey, pExpList, pSchurCompl, pBinvD, pC, pInvD, l2gMap);
216  sys->Initialise(l2gMap);
217  return sys;
218  }
219  }
220 }
GlobalSysSolnType GetGlobalSysSolnType() const
Return the associated solution type.
boost::shared_ptr< GlobalLinSysPETScStaticCond > GlobalLinSysPETScStaticCondSharedPtr
void SetUpMatVec()
Construct PETSc matrix and vector handles.
boost::shared_ptr< AssemblyMap > AssemblyMapSharedPtr
Definition: AssemblyMap.h:53
DNekScalBlkMatSharedPtr m_schurCompl
Block Schur complement matrix.
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
DNekScalBlkMatSharedPtr m_C
Block matrix.
void CalculateReordering(const Array< OneD, const int > &glo2uniMap, const Array< OneD, const int > &glo2unique, const AssemblyMapSharedPtr &pLocToGloMap)
Calculate a reordering of universal IDs for PETSc.
GlobalLinSysPETScStaticCond(const GlobalLinSysKey &mkey, const boost::weak_ptr< ExpList > &pExpList, const boost::shared_ptr< AssemblyMap > &locToGloMap)
Constructor for full direct matrix solve.
void SetUpSolver(NekDouble tolerance)
Set up KSP solver object.
DNekScalBlkMatSharedPtr m_invD
Block matrix.
boost::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
virtual GlobalLinSysStaticCondSharedPtr v_Recurse(const GlobalLinSysKey &mkey, const boost::weak_ptr< ExpList > &pExpList, const DNekScalBlkMatSharedPtr pSchurCompl, const DNekScalBlkMatSharedPtr pBinvD, const DNekScalBlkMatSharedPtr pC, const DNekScalBlkMatSharedPtr pInvD, const boost::shared_ptr< AssemblyMap > &locToGloMap)
boost::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:74
double NekDouble
Describe a linear system.
void SetUpScatter()
Set up PETSc local (equivalent to Nektar++ global) and global (equivalent to universal) scatter maps...
A global linear system.
Definition: GlobalLinSys.h:70
virtual void v_AssembleSchurComplement(boost::shared_ptr< AssemblyMap > locToGloMap)
Assemble the Schur complement matrix.
GlobalLinSysFactory & GetGlobalLinSysFactory()
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:191
DNekScalBlkMatSharedPtr m_BinvD
Block matrix.
static GlobalLinSysSharedPtr create(const GlobalLinSysKey &pLinSysKey, const boost::weak_ptr< ExpList > &pExpList, const boost::shared_ptr< AssemblyMap > &pLocToGloMap)
Creates an instance of this class.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215
boost::shared_ptr< GlobalLinSysStaticCond > GlobalLinSysStaticCondSharedPtr