Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
PreconditionerLinear.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 
42 
43 #ifdef NEKTAR_USING_PETSC
45 #endif
46 
47 #include <LocalRegions/MatrixKey.h>
48 #include <math.h>
49 
50 namespace Nektar
51 {
52  namespace MultiRegions
53  {
54  /**
55  * Registers the class with the Factory.
56  */
57 
60  "FullLinearSpace",
62  "Full Linear space inverse Preconditioning");
63 
66  "LinearPreconSolver",
67  "Xxt");
68  std::string PreconditionerLinear::solveTypeIds[] = {
70  "LinearPreconSolver",
71  "PETSc",
74  "LinearPreconSolver",
75  "Xxt",
77  };
78 
79  /**
80  * @class PreconditionerLinear
81  *
82  * This class implements preconditioning for the conjugate
83  * gradient matrix solver.
84  */
85 
87  const boost::shared_ptr<GlobalLinSys> &plinsys,
88  const AssemblyMapSharedPtr &pLocToGloMap)
89  : Preconditioner(plinsys, pLocToGloMap)
90  {
91  }
92 
94  {
95  }
96 
98  {
99  GlobalSysSolnType sType = m_locToGloMap->GetGlobalSysSolnType();
101  "This type of preconditioning is not implemented "
102  "for this solver");
103 
104  boost::shared_ptr<MultiRegions::ExpList>
105  expList=((m_linsys.lock())->GetLocMat()).lock();
106 
108  expList->GetSession()->GetSolverInfoAsEnum<LinearPreconSolver>(
109  "LinearPreconSolver");
110 
111  GlobalSysSolnType linSolveType;
112 
113  switch(solveType)
114  {
115  case eLinearPreconPETSc:
116  {
117  linSolveType = ePETScFullMatrix;
118 #ifndef NEKTAR_USING_PETSC
119  ASSERTL0(false, "Nektar++ has not been compiled with "
120  "PETSc support.");
121 #endif
122  }
123  case eLinearPreconXxt:
124  default:
125  {
126  linSolveType = eXxtFullMatrix;
127  break;
128  }
129  }
130 
131  m_vertLocToGloMap = m_locToGloMap->LinearSpaceMap(*expList, linSolveType);
132 
133  // Generate linear solve system.
134  StdRegions::MatrixType mType =
135  m_linsys.lock()->GetKey().GetMatrixType() == StdRegions::eMass ?
138 
139  GlobalLinSysKey preconKey(mType, m_vertLocToGloMap, (m_linsys.lock())->GetKey().GetConstFactors());
140 
141  switch(solveType)
142  {
143  case eLinearPreconXxt:
144  {
146  AllocateSharedPtr(preconKey,expList,m_vertLocToGloMap);
147  break;
148  }
149  case eLinearPreconPETSc:
150  {
151 #ifdef NEKTAR_USING_PETSC
153  AllocateSharedPtr(preconKey,expList,m_vertLocToGloMap);
154 #else
155  ASSERTL0(false, "Nektar++ has not been compiled with "
156  "PETSc support.");
157 #endif
158  }
159  }
160  }
161 
162  /**
163  *
164  */
166  const Array<OneD, NekDouble>& pInput,
167  Array<OneD, NekDouble>& pOutput)
168  {
171  }
172 
173  /**
174  *
175  */
177  const Array<OneD, NekDouble>& pInput,
178  Array<OneD, NekDouble>& pOutput,
179  const Array<OneD, NekDouble>& pNonVertOutput,
180  Array<OneD, NekDouble>& pVertForce)
181  {
182  GlobalSysSolnType solvertype=m_locToGloMap->GetGlobalSysSolnType();
183  switch(solvertype)
184  {
186  {
187  int i,val;
188  int nloc = m_vertLocToGloMap->GetNumLocalCoeffs();
189  int nglo = m_vertLocToGloMap->GetNumGlobalCoeffs();
190  // mapping from full space to vertices
191  Array<OneD, int> LocToGloBnd = m_vertLocToGloMap->GetLocalToGlobalBndMap();
192 
193  // Global to local for linear solver (different from above)
194  Array<OneD, int> LocToGlo = m_vertLocToGloMap->GetLocalToGlobalMap();
195 
196  // number of Dir coeffs in from full problem
197  int nDirFull = m_locToGloMap->GetNumGlobalDirBndCoeffs();
198 
199  Array<OneD,NekDouble> In(nglo,0.0);
200  Array<OneD,NekDouble> Out(nglo,0.0);
201 
202  // Gather rhs
203  for(i = 0; i < nloc; ++i)
204  {
205  val = LocToGloBnd[i];
206  if(val >= nDirFull)
207  {
208  In[LocToGlo[i]] = pInput[val-nDirFull];
209  }
210  }
211 
212  // Do solve without enforcing any boundary conditions.
213  m_vertLinsys->SolveLinearSystem(
214  m_vertLocToGloMap->GetNumGlobalCoeffs(),
215  In,Out,m_vertLocToGloMap,
216  m_vertLocToGloMap->GetNumGlobalDirBndCoeffs());
217 
218 
219  if(pNonVertOutput != NullNekDouble1DArray)
220  {
221  ASSERTL1(pNonVertOutput.num_elements() >= pOutput.num_elements(),"Non Vert output is not of sufficient length");
222  Vmath::Vcopy(pOutput.num_elements(),pNonVertOutput,1,pOutput,1);
223  }
224  else
225  {
226  //Copy input to output as a unit preconditioner on
227  //any other value
228  Vmath::Vcopy(pInput.num_elements(),pInput,1,pOutput,1);
229  }
230 
231  if(pVertForce != NullNekDouble1DArray)
232  {
233  Vmath::Zero(pVertForce.num_elements(),pVertForce,1);
234  // Scatter back soln from linear solve
235  for(i = 0; i < nloc; ++i)
236  {
237  val = LocToGloBnd[i];
238  if(val >= nDirFull)
239  {
240  pOutput[val-nDirFull] = Out[LocToGlo[i]];
241  // copy vertex forcing into this vector
242  pVertForce[val-nDirFull] = In[LocToGlo[i]];
243  }
244  }
245  }
246  else
247  {
248  // Scatter back soln from linear solve
249  for(i = 0; i < nloc; ++i)
250  {
251  val = LocToGloBnd[i];
252  if(val >= nDirFull)
253  {
254  pOutput[val-nDirFull] = Out[LocToGlo[i]];
255  }
256  }
257  }
258  }
259  break;
260  default:
261  ASSERTL0(0,"Unsupported solver type");
262  break;
263  }
264  }
265  }
266 }
267 
268 
269 
270 
271 
272