Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
GlobalLinSysIterative.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: GlobalLinSysIterative.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: GlobalLinSysIterative definition
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
37 
38 using namespace std;
39 
40 namespace Nektar
41 {
42  namespace MultiRegions
43  {
44  /**
45  * @class GlobalLinSysIterative
46  *
47  * Solves a linear system using iterative methods.
48  */
49 
50  /// Constructor for full direct matrix solve.
51  GlobalLinSysIterative::GlobalLinSysIterative(
52  const GlobalLinSysKey &pKey,
53  const boost::weak_ptr<ExpList> &pExpList,
54  const boost::shared_ptr<AssemblyMap>
55  &pLocToGloMap)
56  : GlobalLinSys(pKey, pExpList, pLocToGloMap),
57  m_rhs_magnitude(NekConstants::kNekUnsetDouble),
58  m_rhs_mag_sm(0.9),
60  m_totalIterations(0),
61  m_useProjection(false),
62  m_numPrevSols(0)
63  {
64  m_tolerance = pLocToGloMap->GetIterativeTolerance();
65  m_maxiter = pLocToGloMap->GetMaxIterations();
66 
67  LibUtilities::CommSharedPtr vComm = m_expList.lock()->GetComm()->GetRowComm();
68  m_root = (vComm->GetRank())? false : true;
69 
70  int successiveRHS;
71 
72  if((successiveRHS = pLocToGloMap->GetSuccessiveRHS()))
73  {
74  m_prevLinSol.set_capacity(successiveRHS);
75  m_useProjection = true;
76  }
77  else
78  {
79  m_useProjection = false;
80  }
81  }
82 
84  {
85  }
86 
87 
88  /**
89  *
90  */
92  const int nGlobal,
93  const Array<OneD,const NekDouble> &pInput,
94  Array<OneD, NekDouble> &pOutput,
95  const AssemblyMapSharedPtr &plocToGloMap,
96  const int nDir)
97  {
98  if (m_useProjection)
99  {
100  DoAconjugateProjection(nGlobal, pInput, pOutput, plocToGloMap, nDir);
101  }
102  else
103  {
104  // applying plain Conjugate Gradient
105  DoConjugateGradient(nGlobal, pInput, pOutput, plocToGloMap, nDir);
106  }
107  }
108 
109 
110  /**
111  * This method implements A-conjugate projection technique
112  * in order to speed up successive linear solves with
113  * right-hand sides arising from time-dependent discretisations.
114  * (P.F.Fischer, Comput. Methods Appl. Mech. Engrg. 163, 1998)
115  */
117  const int nGlobal,
118  const Array<OneD,const NekDouble> &pInput,
119  Array<OneD, NekDouble> &pOutput,
120  const AssemblyMapSharedPtr &plocToGloMap,
121  const int nDir)
122  {
123  // Get the communicator for performing data exchanges
125  = m_expList.lock()->GetComm()->GetRowComm();
126 
127  // Get vector sizes
128  int nNonDir = nGlobal - nDir;
130 
131  if (0 == m_numPrevSols)
132  {
133  // no previous solutions found, call CG
134 
135  DoConjugateGradient(nGlobal, pInput, pOutput, plocToGloMap, nDir);
136 
137  UpdateKnownSolutions(nGlobal, pOutput, nDir);
138  }
139  else
140  {
141  // Create NekVector wrappers for linear algebra operations
142  NekVector<NekDouble> b (nNonDir, pInput + nDir, eWrapper);
143  NekVector<NekDouble> x (nNonDir, tmp = pOutput + nDir, eWrapper);
144 
145  // check the input vector (rhs) is not zero
146 
147  NekDouble rhsNorm = Vmath::Dot2(nNonDir,
148  pInput + nDir,
149  pInput + nDir,
150  m_map + nDir);
151 
152  vComm->AllReduce(rhsNorm, Nektar::LibUtilities::ReduceSum);
153 
154  if (rhsNorm < NekConstants::kNekZeroTol)
155  {
156  Array<OneD, NekDouble> tmp = pOutput+nDir;
157  Vmath::Zero(nNonDir, tmp, 1);
158  return;
159  }
160 
161  // Allocate array storage
162  Array<OneD, NekDouble> px_s (nGlobal, 0.0);
163  Array<OneD, NekDouble> pb_s (nGlobal, 0.0);
164  Array<OneD, NekDouble> tmpAx_s (nGlobal, 0.0);
165  Array<OneD, NekDouble> tmpx_s (nGlobal, 0.0);
166 
167  NekVector<NekDouble> pb (nNonDir, tmp = pb_s + nDir, eWrapper);
168  NekVector<NekDouble> px (nNonDir, tmp = px_s + nDir, eWrapper);
169  NekVector<NekDouble> tmpAx (nNonDir, tmp = tmpAx_s + nDir, eWrapper);
170  NekVector<NekDouble> tmpx (nNonDir, tmp = tmpx_s + nDir, eWrapper);
171 
172 
173  // notation follows the paper cited:
174  // \alpha_i = \tilda{x_i}^T b^n
175  // projected x, px = \sum \alpha_i \tilda{x_i}
176 
177  Array<OneD, NekDouble> alpha (m_prevLinSol.size(), 0.0);
178  for (int i = 0; i < m_prevLinSol.size(); i++)
179  {
180  alpha[i] = Vmath::Dot2(nNonDir,
181  m_prevLinSol[i],
182  pInput + nDir,
183  m_map + nDir);
184  }
185  vComm->AllReduce(alpha, Nektar::LibUtilities::ReduceSum);
186 
187  for (int i = 0; i < m_prevLinSol.size(); i++)
188  {
189  if (alpha[i] < NekConstants::kNekZeroTol)
190  {
191  continue;
192  }
193 
194  NekVector<NekDouble> xi (nNonDir, m_prevLinSol[i], eWrapper);
195  px += alpha[i] * xi;
196  }
197 
198  // pb = b^n - A px
199  Vmath::Vcopy(nNonDir,
200  pInput.get() + nDir, 1,
201  pb_s.get() + nDir, 1);
202 
203  v_DoMatrixMultiply(px_s, tmpAx_s);
204 
205  pb -= tmpAx;
206 
207 
208  // solve the system with projected rhs
209  DoConjugateGradient(nGlobal, pb_s, tmpx_s, plocToGloMap, nDir);
210 
211 
212  // remainder solution + projection of previous solutions
213  x = tmpx + px;
214 
215  // save the auxiliary solution to prev. known solutions
216  UpdateKnownSolutions(nGlobal, tmpx_s, nDir);
217  }
218  }
219 
220 
221  /**
222  * Calculating A-norm of an input vector,
223  * A-norm(x) := sqrt( < x, Ax > )
224  */
226  const int nGlobal,
227  const Array<OneD,const NekDouble> &in,
228  const int nDir)
229  {
230  // Get the communicator for performing data exchanges
232  = m_expList.lock()->GetComm()->GetRowComm();
233 
234  // Get vector sizes
235  int nNonDir = nGlobal - nDir;
236 
237  // Allocate array storage
238  Array<OneD, NekDouble> tmpAx_s (nGlobal, 0.0);
239 
240  v_DoMatrixMultiply(in, tmpAx_s);
241 
242  NekDouble anorm_sq = Vmath::Dot2(nNonDir,
243  in + nDir,
244  tmpAx_s + nDir,
245  m_map + nDir);
246  vComm->AllReduce(anorm_sq, Nektar::LibUtilities::ReduceSum);
247  return std::sqrt(anorm_sq);
248  }
249 
250  /**
251  * Updates the storage of previously known solutions.
252  * Performs normalisation of input vector wrt A-norm.
253  */
255  const int nGlobal,
256  const Array<OneD,const NekDouble> &newX,
257  const int nDir)
258  {
259  // Get vector sizes
260  int nNonDir = nGlobal - nDir;
261 
262  // Get the communicator for performing data exchanges
264  = m_expList.lock()->GetComm()->GetRowComm();
265 
266  // Check the solution is non-zero
267  NekDouble solNorm = Vmath::Dot2(nNonDir,
268  newX + nDir,
269  newX + nDir,
270  m_map + nDir);
271  vComm->AllReduce(solNorm, Nektar::LibUtilities::ReduceSum);
272 
273  if (solNorm < NekConstants::kNekZeroTol)
274  {
275  return;
276  }
277 
278 
279  // Allocate array storage
280  Array<OneD, NekDouble> tmpAx_s (nGlobal, 0.0);
281  Array<OneD, NekDouble> px_s (nGlobal, 0.0);
282  Array<OneD, NekDouble> tmp1, tmp2;
283 
284  // Create NekVector wrappers for linear algebra operations
285  NekVector<NekDouble> px (nNonDir, tmp1 = px_s + nDir, eWrapper);
286  NekVector<NekDouble> tmpAx (nNonDir, tmp2 = tmpAx_s + nDir, eWrapper);
287 
288 
289  // calculating \tilda{x} - sum \alpha_i\tilda{x}_i
290 
291  Vmath::Vcopy(nNonDir,
292  tmp1 = newX + nDir, 1,
293  tmp2 = px_s + nDir, 1);
294 
295  if (m_prevLinSol.size() > 0)
296  {
297  v_DoMatrixMultiply(newX, tmpAx_s);
298  }
299 
300  Array<OneD, NekDouble> alpha (m_prevLinSol.size(), 0.0);
301  for (int i = 0; i < m_prevLinSol.size(); i++)
302  {
303  alpha[i] = Vmath::Dot2(nNonDir,
304  m_prevLinSol[i],
305  tmpAx_s + nDir,
306  m_map + nDir);
307  }
308  vComm->AllReduce(alpha, Nektar::LibUtilities::ReduceSum);
309 
310  for (int i = 0; i < m_prevLinSol.size(); i++)
311  {
312  if (alpha[i] < NekConstants::kNekZeroTol)
313  {
314  continue;
315  }
316 
317  NekVector<NekDouble> xi (nNonDir, m_prevLinSol[i], eWrapper);
318  px -= alpha[i] * xi;
319  }
320 
321 
322  // Some solutions generated by CG are identical zeros, see
323  // solutions generated for Test_Tet_equitri.xml (IncNavierStokesSolver).
324  // Not going to store identically zero solutions.
325 
326  NekDouble anorm = CalculateAnorm(nGlobal, px_s, nDir);
327  if (anorm < NekConstants::kNekZeroTol)
328  {
329  return;
330  }
331 
332  // normalisation of new solution
333  Vmath::Smul(nNonDir, 1.0/anorm, px_s.get() + nDir, 1, px_s.get() + nDir, 1);
334 
335  // updating storage with non-Dirichlet-dof part of new solution vector
336  m_prevLinSol.push_back(px_s + nDir);
337  m_numPrevSols++;
338  }
339 
340 
341 
342  /**  
343  * Solve a global linear system using the conjugate gradient method.  
344  * We solve only for the non-Dirichlet modes. The operator is evaluated  
345  * using an auxiliary function v_DoMatrixMultiply defined by the  
346  * specific solver. Distributed math routines are used to support  
347  * parallel execution of the solver.  
348  *  
349  * The implemented algorithm uses a reduced-communication reordering of  
350  * the standard PCG method (Demmel, Heath and Vorst, 1993)  
351  *  
352  * @param pInput Input residual of all DOFs.  
353  * @param pOutput Solution vector of all DOFs.  
354  */
356  const int nGlobal,
357  const Array<OneD,const NekDouble> &pInput,
358  Array<OneD, NekDouble> &pOutput,
359  const AssemblyMapSharedPtr &plocToGloMap,
360  const int nDir)
361  {
362  if (!m_precon)
363  {
364  v_UniqueMap();
365  m_precon = CreatePrecon(plocToGloMap);
366  m_precon->BuildPreconditioner();
367  }
368 
369  // Get the communicator for performing data exchanges
371  = m_expList.lock()->GetComm()->GetRowComm();
372 
373  // Get vector sizes
374  int nNonDir = nGlobal - nDir;
375 
376  // Allocate array storage
377  Array<OneD, NekDouble> w_A (nGlobal, 0.0);
378  Array<OneD, NekDouble> s_A (nGlobal, 0.0);
379  Array<OneD, NekDouble> p_A (nNonDir, 0.0);
380  Array<OneD, NekDouble> r_A (nNonDir, 0.0);
381  Array<OneD, NekDouble> q_A (nNonDir, 0.0);
383 
384  // Create NekVector wrappers for linear algebra operations
385  NekVector<NekDouble> in (nNonDir,pInput + nDir, eWrapper);
386  NekVector<NekDouble> out(nNonDir,tmp = pOutput + nDir,eWrapper);
387  NekVector<NekDouble> w (nNonDir,tmp = w_A + nDir, eWrapper);
388  NekVector<NekDouble> s (nNonDir,tmp = s_A + nDir, eWrapper);
389  NekVector<NekDouble> p (nNonDir,p_A, eWrapper);
390  NekVector<NekDouble> r (nNonDir,r_A, eWrapper);
391  NekVector<NekDouble> q (nNonDir,q_A, eWrapper);
392 
393  int k;
394  NekDouble alpha, beta, rho, rho_new, mu, eps, min_resid;
395  Array<OneD, NekDouble> vExchange(3,0.0);
396 
397  // Copy initial residual from input
398  r = in;
399  // zero homogeneous out array ready for solution updates
400  // Should not be earlier in case input vector is same as
401  // output and above copy has been peformed
402  Vmath::Zero(nNonDir,tmp = pOutput + nDir,1);
403 
404 
405  // evaluate initial residual error for exit check
406  vExchange[2] = Vmath::Dot2(nNonDir,
407  r_A,
408  r_A,
409  m_map + nDir);
410 
411  vComm->AllReduce(vExchange, Nektar::LibUtilities::ReduceSum);
412 
413  eps = vExchange[2];
414 
416  {
417  NekVector<NekDouble> inGlob (nGlobal, pInput, eWrapper);
418  Set_Rhs_Magnitude(inGlob);
419  }
420 
421  m_totalIterations = 0;
422 
423  // If input residual is less than tolerance skip solve.
425  {
426  if (m_verbose && m_root)
427  {
428  cout << "CG iterations made = " << m_totalIterations
429  << " using tolerance of " << m_tolerance
430  << " (error = " << sqrt(eps/m_rhs_magnitude)
431  << ", rhs_mag = " << sqrt(m_rhs_magnitude) << ")"
432  << endl;
433  }
434  return;
435  }
436 
437  m_precon->DoPreconditioner(r_A, tmp = w_A + nDir);
438 
439  v_DoMatrixMultiply(w_A, s_A);
440 
441  k = 0;
442 
443  vExchange[0] = Vmath::Dot2(nNonDir,
444  r_A,
445  w_A + nDir,
446  m_map + nDir);
447 
448  vExchange[1] = Vmath::Dot2(nNonDir,
449  s_A + nDir,
450  w_A + nDir,
451  m_map + nDir);
452 
453  vComm->AllReduce(vExchange, Nektar::LibUtilities::ReduceSum);
454 
455  rho = vExchange[0];
456  mu = vExchange[1];
457  min_resid = m_rhs_magnitude;
458  beta = 0.0;
459  alpha = rho/mu;
460  m_totalIterations = 1;
461 
462  // Continue until convergence
463  while (true)
464  {
465  if(k >= m_maxiter)
466  {
467  if (m_root)
468  {
469  cout << "CG iterations made = " << m_totalIterations
470  << " using tolerance of " << m_tolerance
471  << " (error = " << sqrt(eps/m_rhs_magnitude)
472  << ", rhs_mag = " << sqrt(m_rhs_magnitude) << ")"
473  << endl;
474  }
476  "Exceeded maximum number of iterations");
477  }
478 
479  // Compute new search direction p_k, q_k
480  Vmath::Svtvp(nNonDir, beta, &p_A[0], 1, &w_A[nDir], 1, &p_A[0], 1);
481  Vmath::Svtvp(nNonDir, beta, &q_A[0], 1, &s_A[nDir], 1, &q_A[0], 1);
482 
483  // Update solution x_{k+1}
484  Vmath::Svtvp(nNonDir, alpha, &p_A[0], 1, &pOutput[nDir], 1, &pOutput[nDir], 1);
485 
486  // Update residual vector r_{k+1}
487  Vmath::Svtvp(nNonDir, -alpha, &q_A[0], 1, &r_A[0], 1, &r_A[0], 1);
488 
489  // Apply preconditioner
490  m_precon->DoPreconditioner(r_A, tmp = w_A + nDir);
491 
492  // Perform the method-specific matrix-vector multiply operation.
493  v_DoMatrixMultiply(w_A, s_A);
494 
495  // <r_{k+1}, w_{k+1}>
496  vExchange[0] = Vmath::Dot2(nNonDir,
497  r_A,
498  w_A + nDir,
499  m_map + nDir);
500  // <s_{k+1}, w_{k+1}>
501  vExchange[1] = Vmath::Dot2(nNonDir,
502  s_A + nDir,
503  w_A + nDir,
504  m_map + nDir);
505 
506  // <r_{k+1}, r_{k+1}>
507  vExchange[2] = Vmath::Dot2(nNonDir,
508  r_A,
509  r_A,
510  m_map + nDir);
511 
512  // Perform inner-product exchanges
513  vComm->AllReduce(vExchange, Nektar::LibUtilities::ReduceSum);
514 
515  rho_new = vExchange[0];
516  mu = vExchange[1];
517  eps = vExchange[2];
518 
520  // test if norm is within tolerance
521  if (eps < m_tolerance * m_tolerance * m_rhs_magnitude)
522  {
523  if (m_verbose && m_root)
524  {
525  cout << "CG iterations made = " << m_totalIterations
526  << " using tolerance of " << m_tolerance
527  << " (error = " << sqrt(eps/m_rhs_magnitude)
528  << ", rhs_mag = " << sqrt(m_rhs_magnitude) << ")"
529  << endl;
530  }
531  break;
532  }
533  min_resid = min(min_resid, eps);
534 
535  // Compute search direction and solution coefficients
536  beta = rho_new/rho;
537  alpha = rho_new/(mu - rho_new*beta/alpha);
538  rho = rho_new;
539  k++;
540  }
541  }
542 
544  const NekVector<NekDouble> &pIn)
545  {
546  Array<OneD, NekDouble> vExchange(1, 0.0);
547  if (m_map.num_elements() > 0)
548  {
549  vExchange[0] = Vmath::Dot2(pIn.GetDimension(),
550  &pIn[0],&pIn[0],&m_map[0]);
551  }
552 
553  m_expList.lock()->GetComm()->GetRowComm()->AllReduce(
555 
556  // To ensure that very different rhs values are not being
557  // used in subsequent solvers such as the velocit solve in
558  // INC NS. If this works we then need to work out a better
559  // way to control this.
560  NekDouble new_rhs_mag = (vExchange[0] > 1e-6)? vExchange[0] : 1.0;
561 
563  {
564  m_rhs_magnitude = new_rhs_mag;
565  }
566  else
567  {
569  (1.0-m_rhs_mag_sm)*new_rhs_mag);
570  }
571  }
572 
573  }
574 }
boost::shared_ptr< AssemblyMap > AssemblyMapSharedPtr
Definition: AssemblyMap.h:53
NekDouble m_rhs_mag_sm
cnt to how many times rhs_magnitude is called
void Set_Rhs_Magnitude(const NekVector< NekDouble > &pIn)
void UpdateKnownSolutions(const int pGlobalBndDofs, const Array< OneD, const NekDouble > &pSolution, const int pNumDirBndDofs)
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
Definition: Vmath.cpp:485
bool m_useProjection
Whether to apply projection technique.
STL namespace.
boost::circular_buffer< Array< OneD, NekDouble > > m_prevLinSol
Storage for solutions to previous linear problems.
virtual void v_DoMatrixMultiply(const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)=0
NekDouble CalculateAnorm(const int nGlobal, const Array< OneD, const NekDouble > &in, const int nDir)
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:55
static const NekDouble kNekZeroTol
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:213
static PreconditionerSharedPtr NullPreconditionerSharedPtr
NekDouble m_tolerance
Tolerance of iterative solver.
int m_numPrevSols
Total counter of previous solutions.
Array< OneD, int > m_map
Global to universal unique map.
double NekDouble
static const NekDouble kNekUnsetDouble
Describe a linear system.
unsigned int GetDimension() const
Returns the number of dimensions for the point.
Definition: NekVector.cpp:212
void DoAconjugateProjection(const int pNumRows, const Array< OneD, const NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const AssemblyMapSharedPtr &locToGloMap, const int pNumDir)
A-conjugate projection technique.
NekDouble m_rhs_magnitude
dot product of rhs to normalise stopping criterion
A global linear system.
Definition: GlobalLinSys.h:74
T Dot2(int n, const T *w, const T *x, const int *y)
vvtvp (vector times vector times vector): z = w*x*y
Definition: Vmath.cpp:954
#define ROOTONLY_NEKERROR(type, msg)
Definition: ErrorUtil.hpp:195
virtual void v_SolveLinearSystem(const int pNumRows, const Array< OneD, const NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const AssemblyMapSharedPtr &locToGloMap, const int pNumDir)
Solve the matrix system.
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:373
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1061
void DoConjugateGradient(const int pNumRows, const Array< OneD, const NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const AssemblyMapSharedPtr &locToGloMap, const int pNumDir)
Actual iterative solve.
PreconditionerSharedPtr CreatePrecon(AssemblyMapSharedPtr asmMap)
Create a preconditioner object from the parameters defined in the supplied assembly map...
const boost::weak_ptr< ExpList > m_expList
Local Matrix System.
Definition: GlobalLinSys.h:129