Nektar++
Preconditioner.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File Preconditioner.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: Preconditioner definition
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
38 #include <LocalRegions/MatrixKey.h>
39 #include <math.h>
40 
41 namespace Nektar
42 {
43  namespace MultiRegions
44  {
45  std::string Preconditioner::lookupIds[8] = {
47  "Preconditioner", "Null", eNull),
49  "Preconditioner", "Diagonal", eDiagonal),
51  "Preconditioner", "FullLinearSpaceWithDiagonal",
54  "Preconditioner", "FullLinearSpace",eLinear),
56  "Preconditioner", "LowEnergyBlock",eLowEnergy),
58  "Preconditioner", "FullLinearSpaceWithLowEnergyBlock",
61  "Preconditioner", "Block",eBlock),
63  "Preconditioner", "FullLinearSpaceWithBlock",eLinearWithBlock),
64  };
65  std::string Preconditioner::def =
67  "Preconditioner", "Diagonal");
68 
69  /**
70  * @class Preconditioner
71  *
72  * This class implements preconditioning for the conjugate
73  * gradient matrix solver.
74  */
75 
77  const boost::shared_ptr<GlobalLinSys> &plinsys,
78  const AssemblyMapSharedPtr &pLocToGloMap)
79  : m_linsys(plinsys),
80  m_preconType(pLocToGloMap->GetPreconType()),
81  m_locToGloMap(pLocToGloMap)
82  {
83  }
84 
85  /**
86  *
87  */
89  {
90  typedef Loki::SingletonHolder<PreconFactory,
91  Loki::CreateUsingNew,
92  Loki::NoDestroy,
93  Loki::SingleThreaded> Type;
94  return Type::Instance();
95  }
96 
98  {
99  NEKERROR(ErrorUtil::efatal,"Method does not exist" );
100  }
101 
102  /**
103  * \brief Apply a preconditioner to the conjugate gradient method
104  */
106  const Array<OneD, NekDouble>& pInput,
107  Array<OneD, NekDouble>& pOutput)
108  {
109  NEKERROR(ErrorUtil::efatal,"Method does not exist" );
110  }
111 
112  /**
113  * \brief Apply a preconditioner to the conjugate gradient method with
114  * an output for non-vertex degrees of freedom.
115  */
117  const Array<OneD, NekDouble> &pInput,
118  Array<OneD, NekDouble> &pOutput,
119  const Array<OneD, NekDouble> &pNonVertOutput,
120  Array<OneD, NekDouble>& pVertForce)
121  {
122  NEKERROR(ErrorUtil::efatal, "Method does not exist");
123  }
124 
125  /**
126  * \brief Transform from original basis to low energy basis
127  */
129  Array<OneD, NekDouble>& pInOut,
130  int offset)
131  {
132  }
133 
134  /**
135  * \brief Transform from original basis to low energy basis
136  */
138  const Array<OneD, NekDouble> &pInOut,
139  Array<OneD, NekDouble> &pOutput)
140  {
141  }
142 
143  /**
144  * \brief Transform from low energy basis to orignal basis
145  */
147  Array<OneD, NekDouble>& pInput)
148  {
149  Vmath::Smul(pInput.num_elements(), 1.0, pInput, 1, pInput, 1);
150  }
151 
152  /**
153  * \brief Multiply by the block inverse transformation matrix
154  */
156  const Array<OneD, NekDouble> &pInput,
157  Array<OneD, NekDouble> &pOutput)
158  {
159  NEKERROR(ErrorUtil::efatal,"Method does not exist" );
160  }
161 
162  /**
163  * \brief Multiply by the block transposed inverse transformation matrix
164  */
166  const Array<OneD, NekDouble> &pInput,
167  Array<OneD, NekDouble> &pOutput)
168  {
169  NEKERROR(ErrorUtil::efatal,"Method does not exist" );
170  }
171 
173  {
174  }
175 
176  /**
177  * \brief Get block elemental transposed transformation matrix
178  * \f$\mathbf{R}^{T}\f$
179  */
181  int offset, const boost::shared_ptr<DNekScalMat> &loc_mat)
182  {
183  return loc_mat;
184  }
185 
186  /**
187  * @brief Performs global assembly of diagonal entries to global Schur
188  * complement matrix.
189  */
192  {
193  int nGlobalBnd = m_locToGloMap->GetNumGlobalBndCoeffs();
194  int nDirBnd = m_locToGloMap->GetNumGlobalDirBndCoeffs();
195  int rows = nGlobalBnd - nDirBnd;
196 
197  DNekScalBlkMatSharedPtr loc_mat;
198  DNekScalMatSharedPtr bnd_mat;
199  int sign1, sign2, gid1, gid2, i, j, n, cnt;
200  Array<OneD, NekDouble> diagonals(rows,0.0);
201 
202  // Extract diagonal contributions of globally assembled
203  // schur complement matrix
204  for (cnt = n = 0; n < m_linsys.lock()->GetNumBlocks(); ++n)
205  {
206  // Get statically condensed local matrix.
207  loc_mat = (m_linsys.lock())->GetStaticCondBlock(n);
208 
209  // Extract boundary block.
210  bnd_mat = loc_mat->GetBlock(0, 0);
211 
212  // Offset by number of rows.
213  int bnd_row = bnd_mat->GetRows();
214 
215  for (i = 0; i < bnd_row; ++i)
216  {
217  gid1 = m_locToGloMap->GetLocalToGlobalBndMap (cnt + i)
218  - nDirBnd;
219  sign1 = m_locToGloMap->GetLocalToGlobalBndSign(cnt + i);
220 
221  if (gid1 < 0)
222  {
223  continue;
224  }
225 
226  for (j = 0; j < bnd_row; ++j)
227  {
228  gid2 = m_locToGloMap->GetLocalToGlobalBndMap (cnt + j)
229  - nDirBnd;
230  sign2 = m_locToGloMap->GetLocalToGlobalBndSign(cnt + j);
231 
232  if (gid2 == gid1)
233  {
234  diagonals[gid1] += sign1 * sign2 * (*bnd_mat)(i, j);
235  }
236  }
237  }
238  cnt += bnd_row;
239  }
240 
241  return diagonals;
242  }
243  }
244 }
245 
246 
247 
248 
249 
250 
boost::shared_ptr< AssemblyMap > m_locToGloMap
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
Definition: ErrorUtil.hpp:158
virtual void v_DoMultiplybyInverseTransposedTransformationMatrix(const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
Multiply by the block transposed inverse transformation matrix.
static std::string RegisterEnumValue(std::string pEnum, std::string pString, int pEnumValue)
Registers an enumeration value.
boost::shared_ptr< AssemblyMap > AssemblyMapSharedPtr
Definition: AssemblyMap.h:53
Preconditioner(const boost::shared_ptr< GlobalLinSys > &plinsys, const AssemblyMapSharedPtr &pLocToGloMap)
virtual void v_DoPreconditioner(const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
Apply a preconditioner to the conjugate gradient method.
LibUtilities::NekFactory< std::string, Preconditioner, const boost::shared_ptr< GlobalLinSys > &, const boost::shared_ptr< AssemblyMap > & > PreconFactory
PreconFactory & GetPreconFactory()
virtual void v_DoPreconditionerWithNonVertOutput(const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const Array< OneD, NekDouble > &pNonVertOutput, Array< OneD, NekDouble > &pVertForce)
Apply a preconditioner to the conjugate gradient method with an output for non-vertex degrees of free...
virtual void v_DoMultiplybyInverseTransformationMatrix(const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
Multiply by the block inverse transformation matrix.
boost::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
const boost::weak_ptr< GlobalLinSys > m_linsys
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:199
virtual void v_DoTransformFromLowEnergy(Array< OneD, NekDouble > &pInput)
Transform from low energy basis to orignal basis.
Array< OneD, NekDouble > AssembleStaticCondGlobalDiagonals()
Performs global assembly of diagonal entries to global Schur complement matrix.
boost::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:74
virtual DNekScalMatSharedPtr v_TransformedSchurCompl(int offset, const boost::shared_ptr< DNekScalMat > &loc_mat)
Get block elemental transposed transformation matrix .
No Solution type specified.
virtual void v_DoTransformToLowEnergy(Array< OneD, NekDouble > &pInOut, int offset)
Transform from original basis to low energy basis.
static std::string RegisterDefaultSolverInfo(const std::string &pName, const std::string &pValue)
Registers the default string value of a solver info property.
Provides a generic Factory class.
Definition: NekFactory.hpp:116