Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Public Member Functions | Protected Member Functions | Protected Attributes | Private Member Functions | List of all members
Nektar::MultiRegions::GlobalLinSysIterative Class Referenceabstract

A global linear system. More...

#include <GlobalLinSysIterative.h>

Inheritance diagram for Nektar::MultiRegions::GlobalLinSysIterative:
Inheritance graph
[legend]
Collaboration diagram for Nektar::MultiRegions::GlobalLinSysIterative:
Collaboration graph
[legend]

Public Member Functions

 GlobalLinSysIterative (const GlobalLinSysKey &pKey, const boost::weak_ptr< ExpList > &pExpList, const boost::shared_ptr< AssemblyMap > &pLocToGloMap)
 Constructor for full direct matrix solve. More...
 
virtual ~GlobalLinSysIterative ()
 
- Public Member Functions inherited from Nektar::MultiRegions::GlobalLinSys
 GlobalLinSys (const GlobalLinSysKey &pKey, const boost::weak_ptr< ExpList > &pExpList, const boost::shared_ptr< AssemblyMap > &pLocToGloMap)
 Constructor for full direct matrix solve. More...
 
virtual ~GlobalLinSys ()
 
const GlobalLinSysKeyGetKey (void) const
 Returns the key associated with the system. More...
 
const boost::weak_ptr< ExpList > & GetLocMat (void) const
 
void InitObject ()
 
void Initialise (const boost::shared_ptr< AssemblyMap > &pLocToGloMap)
 
void 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. More...
 
boost::shared_ptr< GlobalLinSysGetSharedThisPtr ()
 Returns a shared pointer to the current object. More...
 
int GetNumBlocks ()
 
DNekScalMatSharedPtr GetBlock (unsigned int n)
 
DNekScalBlkMatSharedPtr GetStaticCondBlock (unsigned int n)
 
void DropStaticCondBlock (unsigned int n)
 
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. More...
 

Protected Member Functions

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. More...
 
void DoConjugateGradient (const int pNumRows, const Array< OneD, const NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const AssemblyMapSharedPtr &locToGloMap, const int pNumDir)
 Actual iterative solve. More...
 
void Set_Rhs_Magnitude (const NekVector< NekDouble > &pIn)
 
virtual void v_UniqueMap ()=0
 
- Protected Member Functions inherited from Nektar::MultiRegions::GlobalLinSys
virtual int v_GetNumBlocks ()
 Get the number of blocks in this system. More...
 
virtual DNekScalMatSharedPtr v_GetBlock (unsigned int n)
 Retrieves the block matrix from n-th expansion using the matrix key provided by the m_linSysKey. More...
 
virtual DNekScalBlkMatSharedPtr v_GetStaticCondBlock (unsigned int n)
 Retrieves a the static condensation block matrices from n-th expansion using the matrix key provided by the m_linSysKey. More...
 
virtual void v_DropStaticCondBlock (unsigned int n)
 Releases the static condensation block matrices from NekManager of n-th expansion using the matrix key provided by the m_linSysKey. More...
 
PreconditionerSharedPtr CreatePrecon (AssemblyMapSharedPtr asmMap)
 Create a preconditioner object from the parameters defined in the supplied assembly map. More...
 

Protected Attributes

Array< OneD, int > m_map
 Global to universal unique map. More...
 
int m_maxiter
 maximum iterations More...
 
NekDouble m_tolerance
 Tolerance of iterative solver. More...
 
NekDouble m_rhs_magnitude
 dot product of rhs to normalise stopping criterion More...
 
NekDouble m_rhs_mag_sm
 cnt to how many times rhs_magnitude is called More...
 
PreconditionerSharedPtr m_precon
 
MultiRegions::PreconditionerType m_precontype
 
int m_totalIterations
 
bool m_useProjection
 Whether to apply projection technique. More...
 
bool m_root
 Root if parallel. More...
 
boost::circular_buffer< Array
< OneD, NekDouble > > 
m_prevLinSol
 Storage for solutions to previous linear problems. More...
 
int m_numPrevSols
 Total counter of previous solutions. More...
 
- Protected Attributes inherited from Nektar::MultiRegions::GlobalLinSys
const GlobalLinSysKey m_linSysKey
 Key associated with this linear system. More...
 
const boost::weak_ptr< ExpListm_expList
 Local Matrix System. More...
 
const std::map< int,
RobinBCInfoSharedPtr
m_robinBCInfo
 Robin boundary info. More...
 
bool m_verbose
 

Private Member Functions

void UpdateKnownSolutions (const int pGlobalBndDofs, const Array< OneD, const NekDouble > &pSolution, const int pNumDirBndDofs)
 
NekDouble CalculateAnorm (const int nGlobal, const Array< OneD, const NekDouble > &in, const int nDir)
 
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. More...
 
virtual void v_DoMatrixMultiply (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)=0
 

Detailed Description

A global linear system.

Solves a linear system using iterative methods.

Definition at line 52 of file GlobalLinSysIterative.h.

Constructor & Destructor Documentation

Nektar::MultiRegions::GlobalLinSysIterative::GlobalLinSysIterative ( const GlobalLinSysKey pKey,
const boost::weak_ptr< ExpList > &  pExpList,
const boost::shared_ptr< AssemblyMap > &  pLocToGloMap 
)

Constructor for full direct matrix solve.

Definition at line 51 of file GlobalLinSysIterative.cpp.

References Nektar::MultiRegions::GlobalLinSys::m_expList, m_maxiter, m_prevLinSol, m_root, m_tolerance, and m_useProjection.

56  : GlobalLinSys(pKey, pExpList, pLocToGloMap),
58  m_rhs_mag_sm(0.9),
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  }
NekDouble m_rhs_mag_sm
cnt to how many times rhs_magnitude is called
bool m_useProjection
Whether to apply projection technique.
boost::circular_buffer< Array< OneD, NekDouble > > m_prevLinSol
Storage for solutions to previous linear problems.
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:55
static PreconditionerSharedPtr NullPreconditionerSharedPtr
NekDouble m_tolerance
Tolerance of iterative solver.
int m_numPrevSols
Total counter of previous solutions.
GlobalLinSys(const GlobalLinSysKey &pKey, const boost::weak_ptr< ExpList > &pExpList, const boost::shared_ptr< AssemblyMap > &pLocToGloMap)
Constructor for full direct matrix solve.
static const NekDouble kNekUnsetDouble
NekDouble m_rhs_magnitude
dot product of rhs to normalise stopping criterion
const boost::weak_ptr< ExpList > m_expList
Local Matrix System.
Definition: GlobalLinSys.h:129
Nektar::MultiRegions::GlobalLinSysIterative::~GlobalLinSysIterative ( )
virtual

Definition at line 83 of file GlobalLinSysIterative.cpp.

84  {
85  }

Member Function Documentation

NekDouble Nektar::MultiRegions::GlobalLinSysIterative::CalculateAnorm ( const int  nGlobal,
const Array< OneD, const NekDouble > &  in,
const int  nDir 
)
private

Calculating A-norm of an input vector, A-norm(x) := sqrt( < x, Ax > )

Definition at line 225 of file GlobalLinSysIterative.cpp.

References Vmath::Dot2(), Nektar::MultiRegions::GlobalLinSys::m_expList, m_map, Nektar::LibUtilities::ReduceSum, and v_DoMatrixMultiply().

Referenced by UpdateKnownSolutions().

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  }
virtual void v_DoMatrixMultiply(const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)=0
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:55
Array< OneD, int > m_map
Global to universal unique map.
double NekDouble
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
const boost::weak_ptr< ExpList > m_expList
Local Matrix System.
Definition: GlobalLinSys.h:129
void Nektar::MultiRegions::GlobalLinSysIterative::DoAconjugateProjection ( const int  nGlobal,
const Array< OneD, const NekDouble > &  pInput,
Array< OneD, NekDouble > &  pOutput,
const AssemblyMapSharedPtr plocToGloMap,
const int  nDir 
)
protected

A-conjugate projection technique.

This method implements A-conjugate projection technique in order to speed up successive linear solves with right-hand sides arising from time-dependent discretisations. (P.F.Fischer, Comput. Methods Appl. Mech. Engrg. 163, 1998)

Definition at line 116 of file GlobalLinSysIterative.cpp.

References DoConjugateGradient(), Vmath::Dot2(), Nektar::eWrapper, Nektar::NekConstants::kNekZeroTol, Nektar::MultiRegions::GlobalLinSys::m_expList, m_map, m_numPrevSols, m_prevLinSol, Nektar::LibUtilities::ReduceSum, UpdateKnownSolutions(), v_DoMatrixMultiply(), Vmath::Vcopy(), and Vmath::Zero().

Referenced by v_SolveLinearSystem().

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;
129  Array<OneD, NekDouble> tmp;
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  }
void UpdateKnownSolutions(const int pGlobalBndDofs, const Array< OneD, const NekDouble > &pSolution, const int pNumDirBndDofs)
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
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:55
static const NekDouble kNekZeroTol
int m_numPrevSols
Total counter of previous solutions.
Array< OneD, int > m_map
Global to universal unique map.
double NekDouble
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
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.
const boost::weak_ptr< ExpList > m_expList
Local Matrix System.
Definition: GlobalLinSys.h:129
void Nektar::MultiRegions::GlobalLinSysIterative::DoConjugateGradient ( const int  nGlobal,
const Array< OneD, const NekDouble > &  pInput,
Array< OneD, NekDouble > &  pOutput,
const AssemblyMapSharedPtr plocToGloMap,
const int  nDir 
)
protected

Actual iterative solve.

  Solve a global linear system using the conjugate gradient method.   We solve only for the non-Dirichlet modes. The operator is evaluated   using an auxiliary function v_DoMatrixMultiply defined by the   specific solver. Distributed math routines are used to support   parallel execution of the solver.     The implemented algorithm uses a reduced-communication reordering of   the standard PCG method (Demmel, Heath and Vorst, 1993)    

Parameters
pInputInput residual of all DOFs.  
pOutputSolution vector of all DOFs.  

Definition at line 355 of file GlobalLinSysIterative.cpp.

References Nektar::MultiRegions::GlobalLinSys::CreatePrecon(), Vmath::Dot2(), ErrorUtil::efatal, Nektar::eWrapper, Nektar::NekConstants::kNekUnsetDouble, Nektar::MultiRegions::GlobalLinSys::m_expList, m_map, m_maxiter, m_precon, m_rhs_magnitude, m_root, m_tolerance, m_totalIterations, Nektar::MultiRegions::GlobalLinSys::m_verbose, CellMLToNektar.cellml_metadata::p, Nektar::LibUtilities::ReduceSum, ROOTONLY_NEKERROR, Set_Rhs_Magnitude(), Vmath::Svtvp(), v_DoMatrixMultiply(), v_UniqueMap(), and Vmath::Zero().

Referenced by DoAconjugateProjection(), and v_SolveLinearSystem().

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);
382  Array<OneD, NekDouble> tmp;
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  }
void Set_Rhs_Magnitude(const NekVector< NekDouble > &pIn)
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
virtual void v_DoMatrixMultiply(const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput)=0
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:55
NekDouble m_tolerance
Tolerance of iterative solver.
Array< OneD, int > m_map
Global to universal unique map.
double NekDouble
static const NekDouble kNekUnsetDouble
NekDouble m_rhs_magnitude
dot product of rhs to normalise stopping criterion
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
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:373
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
void Nektar::MultiRegions::GlobalLinSysIterative::Set_Rhs_Magnitude ( const NekVector< NekDouble > &  pIn)
protected

Definition at line 543 of file GlobalLinSysIterative.cpp.

References Vmath::Dot2(), Nektar::NekVector< DataType >::GetDimension(), Nektar::NekConstants::kNekUnsetDouble, Nektar::MultiRegions::GlobalLinSys::m_expList, m_map, m_rhs_mag_sm, m_rhs_magnitude, and Nektar::LibUtilities::ReduceSum.

Referenced by DoConjugateGradient(), and Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::v_PreSolve().

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  }
NekDouble m_rhs_mag_sm
cnt to how many times rhs_magnitude is called
Array< OneD, int > m_map
Global to universal unique map.
double NekDouble
static const NekDouble kNekUnsetDouble
unsigned int GetDimension() const
Returns the number of dimensions for the point.
Definition: NekVector.cpp:212
NekDouble m_rhs_magnitude
dot product of rhs to normalise stopping criterion
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
const boost::weak_ptr< ExpList > m_expList
Local Matrix System.
Definition: GlobalLinSys.h:129
void Nektar::MultiRegions::GlobalLinSysIterative::UpdateKnownSolutions ( const int  nGlobal,
const Array< OneD, const NekDouble > &  newX,
const int  nDir 
)
private

Updates the storage of previously known solutions. Performs normalisation of input vector wrt A-norm.

Definition at line 254 of file GlobalLinSysIterative.cpp.

References CalculateAnorm(), Vmath::Dot2(), Nektar::eWrapper, Nektar::NekConstants::kNekZeroTol, Nektar::MultiRegions::GlobalLinSys::m_expList, m_map, m_numPrevSols, m_prevLinSol, Nektar::LibUtilities::ReduceSum, Vmath::Smul(), v_DoMatrixMultiply(), and Vmath::Vcopy().

Referenced by DoAconjugateProjection().

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  }
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
int m_numPrevSols
Total counter of previous solutions.
Array< OneD, int > m_map
Global to universal unique map.
double NekDouble
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
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1061
const boost::weak_ptr< ExpList > m_expList
Local Matrix System.
Definition: GlobalLinSys.h:129
virtual void Nektar::MultiRegions::GlobalLinSysIterative::v_DoMatrixMultiply ( const Array< OneD, NekDouble > &  pInput,
Array< OneD, NekDouble > &  pOutput 
)
privatepure virtual
void Nektar::MultiRegions::GlobalLinSysIterative::v_SolveLinearSystem ( const int  pNumRows,
const Array< OneD, const NekDouble > &  pInput,
Array< OneD, NekDouble > &  pOutput,
const AssemblyMapSharedPtr locToGloMap,
const int  pNumDir 
)
privatevirtual

Solve the matrix system.

Implements Nektar::MultiRegions::GlobalLinSys.

Definition at line 91 of file GlobalLinSysIterative.cpp.

References DoAconjugateProjection(), DoConjugateGradient(), and m_useProjection.

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  }
bool m_useProjection
Whether to apply projection technique.
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.
void DoConjugateGradient(const int pNumRows, const Array< OneD, const NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const AssemblyMapSharedPtr &locToGloMap, const int pNumDir)
Actual iterative solve.
virtual void Nektar::MultiRegions::GlobalLinSysIterative::v_UniqueMap ( )
protectedpure virtual

Member Data Documentation

Array<OneD, int> Nektar::MultiRegions::GlobalLinSysIterative::m_map
protected
int Nektar::MultiRegions::GlobalLinSysIterative::m_maxiter
protected

maximum iterations

Definition at line 68 of file GlobalLinSysIterative.h.

Referenced by DoConjugateGradient(), and GlobalLinSysIterative().

int Nektar::MultiRegions::GlobalLinSysIterative::m_numPrevSols
protected

Total counter of previous solutions.

Definition at line 95 of file GlobalLinSysIterative.h.

Referenced by DoAconjugateProjection(), and UpdateKnownSolutions().

PreconditionerSharedPtr Nektar::MultiRegions::GlobalLinSysIterative::m_precon
protected
MultiRegions::PreconditionerType Nektar::MultiRegions::GlobalLinSysIterative::m_precontype
protected

Definition at line 81 of file GlobalLinSysIterative.h.

boost::circular_buffer<Array<OneD, NekDouble> > Nektar::MultiRegions::GlobalLinSysIterative::m_prevLinSol
protected

Storage for solutions to previous linear problems.

Definition at line 92 of file GlobalLinSysIterative.h.

Referenced by DoAconjugateProjection(), GlobalLinSysIterative(), and UpdateKnownSolutions().

NekDouble Nektar::MultiRegions::GlobalLinSysIterative::m_rhs_mag_sm
protected

cnt to how many times rhs_magnitude is called

Definition at line 77 of file GlobalLinSysIterative.h.

Referenced by Set_Rhs_Magnitude().

NekDouble Nektar::MultiRegions::GlobalLinSysIterative::m_rhs_magnitude
protected

dot product of rhs to normalise stopping criterion

Definition at line 74 of file GlobalLinSysIterative.h.

Referenced by DoConjugateGradient(), Set_Rhs_Magnitude(), and Nektar::MultiRegions::GlobalLinSysIterativeStaticCond::v_PreSolve().

bool Nektar::MultiRegions::GlobalLinSysIterative::m_root
protected

Root if parallel.

Definition at line 89 of file GlobalLinSysIterative.h.

Referenced by DoConjugateGradient(), and GlobalLinSysIterative().

NekDouble Nektar::MultiRegions::GlobalLinSysIterative::m_tolerance
protected

Tolerance of iterative solver.

Definition at line 71 of file GlobalLinSysIterative.h.

Referenced by DoConjugateGradient(), and GlobalLinSysIterative().

int Nektar::MultiRegions::GlobalLinSysIterative::m_totalIterations
protected

Definition at line 83 of file GlobalLinSysIterative.h.

Referenced by DoConjugateGradient().

bool Nektar::MultiRegions::GlobalLinSysIterative::m_useProjection
protected

Whether to apply projection technique.

Definition at line 86 of file GlobalLinSysIterative.h.

Referenced by GlobalLinSysIterative(), and v_SolveLinearSystem().