Nektar++
GlobalLinSysIterativeFull.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 
36 #include <map>
39 
40 namespace Nektar
41 {
42  namespace MultiRegions
43  {
44  /**
45  * @class GlobalLinSysIterativeCG
46  *
47  * This class implements a conjugate gradient matrix solver.
48  * Preconditioning is implemented using a Jacobi (diagonal)
49  * preconditioner.
50  */
51 
52  /**
53  * Registers the class with the Factory.
54  */
57  "IterativeFull",
59  "Iterative solver for full matrix system.");
60 
61 
62  /**
63  * Constructor for full direct matrix solve.
64  * @param pKey Key specifying matrix to solve.
65  * @param pExp Shared pointer to expansion list for applying
66  * matrix evaluations.
67  * @param pLocToGloMap Local to global mapping.
68  */
70  const GlobalLinSysKey &pKey,
71  const boost::weak_ptr<ExpList> &pExp,
72  const boost::shared_ptr<AssemblyMap> &pLocToGloMap)
73  : GlobalLinSys (pKey, pExp, pLocToGloMap),
74  GlobalLinSysIterative(pKey, pExp, pLocToGloMap)
75  {
77  "This routine should only be used when using an Iterative "
78  "conjugate gradient matrix solve.");
79  }
80 
81 
82  /**
83  *
84  */
86  {
87 
88  }
89 
90 
91  /**
92  * Solve a global linear system with Dirichlet forcing using a
93  * conjugate gradient method. This routine performs handling of the
94  * Dirichlet forcing terms and wraps the underlying iterative solver
95  * used for the remaining degrees of freedom.
96  *
97  * Consider solving for \f$x\f$, the matrix system \f$Ax=b\f$, where
98  * \f$b\f$ is known. To enforce the Dirichlet terms we instead solve
99  * \f[A(x-x_0) = b - Ax_0 \f]
100  * where \f$x_0\f$ is the Dirichlet forcing.
101  *
102  * @param pInput RHS of linear system, \f$b\f$.
103  * @param pOutput On input, values of dirichlet degrees
104  * of freedom with initial guess on other values.
105  * On output, the solution \f$x\f$.
106  * @param pLocToGloMap Local to global mapping.
107  * @param pDirForcing Precalculated Dirichlet forcing.
108  */
110  const Array<OneD, const NekDouble> &pInput,
111  Array<OneD, NekDouble> &pOutput,
112  const AssemblyMapSharedPtr &pLocToGloMap,
113  const Array<OneD, const NekDouble> &pDirForcing)
114  {
115  boost::shared_ptr<MultiRegions::ExpList> expList = m_expList.lock();
116  bool vCG;
117  if ((m_locToGloMap = boost::dynamic_pointer_cast<AssemblyMapCG>(
118  pLocToGloMap)))
119  {
120  vCG = true;
121  }
122  else if ((m_locToGloMap = boost::dynamic_pointer_cast<
123  AssemblyMapDG>(pLocToGloMap)))
124  {
125  vCG = false;
126  }
127  else
128  {
129  ASSERTL0(false, "Unknown map type");
130  }
131 
132  bool dirForcCalculated = (bool) pDirForcing.num_elements();
133  int nDirDofs = pLocToGloMap->GetNumGlobalDirBndCoeffs();
134  int nGlobDofs = pLocToGloMap->GetNumGlobalCoeffs();
135  int nDirTotal = nDirDofs;
136 
137  expList->GetComm()->AllReduce(nDirTotal, LibUtilities::ReduceSum);
138 
139  Array<OneD, NekDouble> tmp(nGlobDofs), tmp2;
140 
141  if(nDirTotal)
142  {
143  // calculate the Dirichlet forcing
144  if(dirForcCalculated)
145  {
146  Vmath::Vsub(nGlobDofs, pInput.get(), 1,
147  pDirForcing.get(), 1,
148  tmp.get(), 1);
149  }
150  else
151  {
152  // Calculate the dirichlet forcing B_b (== X_b) and
153  // substract it from the rhs
154  expList->GeneralMatrixOp(
155  m_linSysKey, pOutput, tmp, eGlobal);
156 
157  Vmath::Vsub(nGlobDofs, pInput.get(), 1,
158  tmp.get(), 1,
159  tmp.get(), 1);
160  }
161  if (vCG)
162  {
163  Array<OneD, NekDouble> out(nGlobDofs,0.0);
164 
165  // solve for perturbation from intiial guess in pOutput
167  nGlobDofs, tmp, out, pLocToGloMap, nDirDofs);
168  Vmath::Vadd(nGlobDofs-nDirDofs, &out [nDirDofs], 1,
169  &pOutput[nDirDofs], 1, &pOutput[nDirDofs], 1);
170  }
171  else
172  {
173  ASSERTL0(false, "Need DG solve if using Dir BCs");
174  }
175  }
176  else
177  {
178  Vmath::Vcopy(nGlobDofs, pInput, 1, tmp, 1);
179  SolveLinearSystem(nGlobDofs, tmp, pOutput, pLocToGloMap);
180  }
181  }
182 
183 
184  /**
185  *
186  */
188  const Array<OneD, NekDouble>& pInput,
189  Array<OneD, NekDouble>& pOutput)
190  {
191  boost::shared_ptr<MultiRegions::ExpList> expList = m_expList.lock();
192  // Perform matrix-vector operation A*d_i
193  expList->GeneralMatrixOp(m_linSysKey,
194  pInput, pOutput, eGlobal);
195 
196  // retrieve robin boundary condition information and apply robin
197  // boundary conditions to the solution.
198  const std::map<int, RobinBCInfoSharedPtr> vRobinBCInfo
199  = expList->GetRobinBCInfo();
200  if(vRobinBCInfo.size() > 0)
201  {
202  ASSERTL0(false,
203  "Robin boundaries not set up in IterativeFull solver.");
204  int nGlobal = m_locToGloMap->GetNumGlobalCoeffs();
205  int nLocal = m_locToGloMap->GetNumLocalCoeffs();
206  int nDir = m_locToGloMap->GetNumGlobalDirBndCoeffs();
207  int nNonDir = nGlobal - nDir;
208  Array<OneD, NekDouble> robin_A(nGlobal, 0.0);
209  Array<OneD, NekDouble> robin_l(nLocal, 0.0);
211  NekVector<NekDouble> robin(nNonDir,
212  tmp = robin_A + nDir, eWrapper);
213 
214  // Operation: p_A = A * d_A
215  // First map d_A to local solution
216  m_locToGloMap->GlobalToLocal(pInput, robin_l);
217 
218  // Iterate over all the elements computing Robin BCs where
219  // necessary
220  for (int n = 0; n < expList->GetNumElmts(); ++n)
221  {
222  int nel = expList->GetOffset_Elmt_Id(n);
223  int offset = expList->GetCoeff_Offset(n);
224  int ncoeffs = expList->GetExp(nel)->GetNcoeffs();
225 
226  if(vRobinBCInfo.count(nel) != 0) // add robin mass matrix
227  {
230  StdRegions::StdExpansionSharedPtr vExp = expList->GetExp(nel);
231 
232  // add local matrix contribution
233  for(rBC = vRobinBCInfo.find(nel)->second;rBC; rBC = rBC->next)
234  {
235  vExp->AddRobinEdgeContribution(rBC->m_robinID,rBC->m_robinPrimitiveCoeffs, tmp = robin_l + offset);
236  }
237  }
238  else
239  {
240  Vmath::Zero(ncoeffs, &robin_l[offset], 1);
241  }
242  }
243 
244  // Map local Robin contribution back to global coefficients
245  m_locToGloMap->LocalToGlobal(robin_l, robin_A);
246  // Add them to the output of the GeneralMatrixOp
247  Vmath::Vadd(nGlobal, pOutput, 1, robin_A, 1, pOutput, 1);
248  }
249 
250  }
251 
252  /**
253  *
254  */
256  {
257  m_map = m_locToGloMap->GetGlobalToUniversalMapUnique();
258  }
259 
260  }
261 }
GlobalSysSolnType GetGlobalSysSolnType() const
Return the associated solution type.
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
boost::shared_ptr< AssemblyMap > AssemblyMapSharedPtr
Definition: AssemblyMap.h:53
boost::shared_ptr< RobinBCInfo > RobinBCInfoSharedPtr
void SolveLinearSystem(const int pNumRows, const Array< OneD, const NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const AssemblyMapSharedPtr &locToGloMap, const int pNumDir=0)
Solve the linear system for given input and output vectors.
Definition: GlobalLinSys.h:193
Global coefficients.
Array< OneD, int > m_map
Global to universal unique map.
Describe a linear system.
const GlobalLinSysKey m_linSysKey
Key associated with this linear system.
Definition: GlobalLinSys.h:123
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:329
A global linear system.
Definition: GlobalLinSys.h:70
GlobalLinSysIterativeFull(const GlobalLinSysKey &pLinSysKey, const boost::weak_ptr< ExpList > &pExpList, const boost::shared_ptr< AssemblyMap > &pLocToGloMap)
Constructor for full direct matrix solve.
static GlobalLinSysSharedPtr create(const GlobalLinSysKey &pLinSysKey, const boost::weak_ptr< ExpList > &pExpList, const boost::shared_ptr< AssemblyMap > &pLocToGloMap)
Creates an instance of this class.
virtual void v_Solve(const Array< OneD, const NekDouble > &in, Array< OneD, NekDouble > &out, const AssemblyMapSharedPtr &locToGloMap, const Array< OneD, const NekDouble > &dirForcing=NullNekDouble1DArray)
Solve the linear system for given input and output vectors using a specified local to global map...
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
GlobalLinSysFactory & GetGlobalLinSysFactory()
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:191
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1038
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:285
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215
virtual void v_DoMatrixMultiply(const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)
const boost::weak_ptr< ExpList > m_expList
Local Matrix System.
Definition: GlobalLinSys.h:125