Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator 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
 Provide verbose output and root if parallel. More...
 
bool m_verbose
 
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 map< int,
RobinBCInfoSharedPtr
m_robinBCInfo
 Robin boundary info. More...
 

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 49 of file GlobalLinSysIterative.cpp.

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

54  : GlobalLinSys(pKey, pExpList, pLocToGloMap),
56  m_rhs_mag_sm(0.9),
59  m_useProjection(false),
60  m_numPrevSols(0)
61  {
63  = pExpList.lock()->GetSession();
64 
65  m_tolerance = pLocToGloMap->GetIterativeTolerance();
66  m_maxiter = pLocToGloMap->GetMaxIterations();
67 
68  LibUtilities::CommSharedPtr vComm = m_expList.lock()->GetComm()->GetRowComm();
69  m_root = (vComm->GetRank())? false : true;
70  m_verbose = (vSession->DefinesCmdLineArgument("verbose"))? true :false;
71 
72  int successiveRHS;
73 
74  if((successiveRHS = pLocToGloMap->GetSuccessiveRHS()))
75  {
76  m_prevLinSol.set_capacity(successiveRHS);
77  m_useProjection = true;
78  }
79  else
80  {
81  m_useProjection = false;
82  }
83  }
NekDouble m_rhs_mag_sm
cnt to how many times rhs_magnitude is called
bool m_root
Provide verbose output and root if parallel.
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< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:53
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 85 of file GlobalLinSysIterative.cpp.

86  {
87  }

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 227 of file GlobalLinSysIterative.cpp.

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

Referenced by UpdateKnownSolutions().

231  {
232  // Get the communicator for performing data exchanges
234  = m_expList.lock()->GetComm()->GetRowComm();
235 
236  // Get vector sizes
237  int nNonDir = nGlobal - nDir;
238 
239  // Allocate array storage
240  Array<OneD, NekDouble> tmpAx_s (nGlobal, 0.0);
241 
242  v_DoMatrixMultiply(in, tmpAx_s);
243 
244  NekDouble anorm_sq = Vmath::Dot2(nNonDir,
245  in + nDir,
246  tmpAx_s + nDir,
247  m_map + nDir);
248  vComm->AllReduce(anorm_sq, Nektar::LibUtilities::ReduceSum);
249  return std::sqrt(anorm_sq);
250  }
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:53
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:940
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 118 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().

124  {
125  // Get the communicator for performing data exchanges
127  = m_expList.lock()->GetComm()->GetRowComm();
128 
129  // Get vector sizes
130  int nNonDir = nGlobal - nDir;
132 
133  if (0 == m_numPrevSols)
134  {
135  // no previous solutions found, call CG
136 
137  DoConjugateGradient(nGlobal, pInput, pOutput, plocToGloMap, nDir);
138 
139  UpdateKnownSolutions(nGlobal, pOutput, nDir);
140  }
141  else
142  {
143  // Create NekVector wrappers for linear algebra operations
144  NekVector<NekDouble> b (nNonDir, pInput + nDir, eWrapper);
145  NekVector<NekDouble> x (nNonDir, tmp = pOutput + nDir, eWrapper);
146 
147  // check the input vector (rhs) is not zero
148 
149  NekDouble rhsNorm = Vmath::Dot2(nNonDir,
150  pInput + nDir,
151  pInput + nDir,
152  m_map + nDir);
153 
154  vComm->AllReduce(rhsNorm, Nektar::LibUtilities::ReduceSum);
155 
156  if (rhsNorm < NekConstants::kNekZeroTol)
157  {
158  Array<OneD, NekDouble> tmp = pOutput+nDir;
159  Vmath::Zero(nNonDir, tmp, 1);
160  return;
161  }
162 
163  // Allocate array storage
164  Array<OneD, NekDouble> px_s (nGlobal, 0.0);
165  Array<OneD, NekDouble> pb_s (nGlobal, 0.0);
166  Array<OneD, NekDouble> tmpAx_s (nGlobal, 0.0);
167  Array<OneD, NekDouble> tmpx_s (nGlobal, 0.0);
168 
169  NekVector<NekDouble> pb (nNonDir, tmp = pb_s + nDir, eWrapper);
170  NekVector<NekDouble> px (nNonDir, tmp = px_s + nDir, eWrapper);
171  NekVector<NekDouble> tmpAx (nNonDir, tmp = tmpAx_s + nDir, eWrapper);
172  NekVector<NekDouble> tmpx (nNonDir, tmp = tmpx_s + nDir, eWrapper);
173 
174 
175  // notation follows the paper cited:
176  // \alpha_i = \tilda{x_i}^T b^n
177  // projected x, px = \sum \alpha_i \tilda{x_i}
178 
179  Array<OneD, NekDouble> alpha (m_prevLinSol.size(), 0.0);
180  for (int i = 0; i < m_prevLinSol.size(); i++)
181  {
182  alpha[i] = Vmath::Dot2(nNonDir,
183  m_prevLinSol[i],
184  pInput + nDir,
185  m_map + nDir);
186  }
187  vComm->AllReduce(alpha, Nektar::LibUtilities::ReduceSum);
188 
189  for (int i = 0; i < m_prevLinSol.size(); i++)
190  {
191  if (alpha[i] < NekConstants::kNekZeroTol)
192  {
193  continue;
194  }
195 
196  NekVector<NekDouble> xi (nNonDir, m_prevLinSol[i], eWrapper);
197  px += alpha[i] * xi;
198  }
199 
200  // pb = b^n - A px
201  Vmath::Vcopy(nNonDir,
202  pInput.get() + nDir, 1,
203  pb_s.get() + nDir, 1);
204 
205  v_DoMatrixMultiply(px_s, tmpAx_s);
206 
207  pb -= tmpAx;
208 
209 
210  // solve the system with projected rhs
211  DoConjugateGradient(nGlobal, pb_s, tmpx_s, plocToGloMap, nDir);
212 
213 
214  // remainder solution + projection of previous solutions
215  x = tmpx + px;
216 
217  // save the auxiliary solution to prev. known solutions
218  UpdateKnownSolutions(nGlobal, tmpx_s, nDir);
219  }
220  }
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:53
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:940
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
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 357 of file GlobalLinSysIterative.cpp.

References ASSERTL0, Nektar::MultiRegions::GlobalLinSys::CreatePrecon(), Vmath::Dot2(), 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, m_verbose, Nektar::LibUtilities::ReduceSum, Set_Rhs_Magnitude(), Vmath::Svtvp(), v_DoMatrixMultiply(), v_UniqueMap(), and Vmath::Zero().

Referenced by DoAconjugateProjection(), and v_SolveLinearSystem().

363  {
364  if (!m_precon)
365  {
366  v_UniqueMap();
367  m_precon = CreatePrecon(plocToGloMap);
368  m_precon->BuildPreconditioner();
369  }
370 
371  // Get the communicator for performing data exchanges
373  = m_expList.lock()->GetComm()->GetRowComm();
374 
375  // Get vector sizes
376  int nNonDir = nGlobal - nDir;
377 
378  // Allocate array storage
379  Array<OneD, NekDouble> w_A (nGlobal, 0.0);
380  Array<OneD, NekDouble> s_A (nGlobal, 0.0);
381  Array<OneD, NekDouble> p_A (nNonDir, 0.0);
382  Array<OneD, NekDouble> r_A (nNonDir, 0.0);
383  Array<OneD, NekDouble> q_A (nNonDir, 0.0);
385 
386  // Create NekVector wrappers for linear algebra operations
387  NekVector<NekDouble> in (nNonDir,pInput + nDir, eWrapper);
388  NekVector<NekDouble> out(nNonDir,tmp = pOutput + nDir,eWrapper);
389  NekVector<NekDouble> w (nNonDir,tmp = w_A + nDir, eWrapper);
390  NekVector<NekDouble> s (nNonDir,tmp = s_A + nDir, eWrapper);
391  NekVector<NekDouble> p (nNonDir,p_A, eWrapper);
392  NekVector<NekDouble> r (nNonDir,r_A, eWrapper);
393  NekVector<NekDouble> q (nNonDir,q_A, eWrapper);
394 
395  int k;
396  NekDouble alpha, beta, rho, rho_new, mu, eps, min_resid;
397  Array<OneD, NekDouble> vExchange(3,0.0);
398 
399  // Copy initial residual from input
400  r = in;
401  // zero homogeneous out array ready for solution updates
402  // Should not be earlier in case input vector is same as
403  // output and above copy has been peformed
404  Vmath::Zero(nNonDir,tmp = pOutput + nDir,1);
405 
406 
407  // evaluate initial residual error for exit check
408  vExchange[2] = Vmath::Dot2(nNonDir,
409  r_A,
410  r_A,
411  m_map + nDir);
412 
413  vComm->AllReduce(vExchange, Nektar::LibUtilities::ReduceSum);
414 
415  eps = vExchange[2];
416 
418  {
419  NekVector<NekDouble> inGlob (nGlobal, pInput, eWrapper);
420  Set_Rhs_Magnitude(inGlob);
421  }
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_totalIterations = 1;
438  m_precon->DoPreconditioner(r_A, tmp = w_A + nDir);
439 
440  v_DoMatrixMultiply(w_A, s_A);
441 
442  k = 0;
443 
444  vExchange[0] = Vmath::Dot2(nNonDir,
445  r_A,
446  w_A + nDir,
447  m_map + nDir);
448 
449  vExchange[1] = Vmath::Dot2(nNonDir,
450  s_A + nDir,
451  w_A + nDir,
452  m_map + nDir);
453 
454  vComm->AllReduce(vExchange, Nektar::LibUtilities::ReduceSum);
455 
456  rho = vExchange[0];
457  mu = vExchange[1];
458  min_resid = m_rhs_magnitude;
459  beta = 0.0;
460  alpha = rho/mu;
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  ASSERTL0(false,
475  "Exceeded maximum number of iterations");
476  }
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  }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
void Set_Rhs_Magnitude(const NekVector< NekDouble > &pIn)
bool m_root
Provide verbose output and root if parallel.
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:471
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:53
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:940
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
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::Dot(), Nektar::NekVector< DataType >::GetDimension(), Nektar::NekConstants::kNekUnsetDouble, Nektar::MultiRegions::GlobalLinSys::m_expList, 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);
547  vExchange[0] = Vmath::Dot(pIn.GetDimension(),&pIn[0],&pIn[0]);
548 
549  m_expList.lock()->GetComm()->GetRowComm()->AllReduce(
551 
552  // To ensure that very different rhs values are not being
553  // used in subsequent solvers such as the velocit solve in
554  // INC NS. If this works we then need to work out a better
555  // way to control this.
556  NekDouble new_rhs_mag = (vExchange[0] > 1e-6)? vExchange[0] : 1.0;
557 
559  {
560  m_rhs_magnitude = new_rhs_mag;
561  }
562  else
563  {
565  (1.0-m_rhs_mag_sm)*new_rhs_mag);
566  }
567  }
NekDouble m_rhs_mag_sm
cnt to how many times rhs_magnitude is called
double NekDouble
static const NekDouble kNekUnsetDouble
T Dot(int n, const T *w, const T *x)
vvtvp (vector times vector times vector): z = w*x*y
Definition: Vmath.cpp:900
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
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 256 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().

260  {
261  // Get vector sizes
262  int nNonDir = nGlobal - nDir;
263 
264  // Get the communicator for performing data exchanges
266  = m_expList.lock()->GetComm()->GetRowComm();
267 
268  // Check the solution is non-zero
269  NekDouble solNorm = Vmath::Dot2(nNonDir,
270  newX + nDir,
271  newX + nDir,
272  m_map + nDir);
273  vComm->AllReduce(solNorm, Nektar::LibUtilities::ReduceSum);
274 
275  if (solNorm < NekConstants::kNekZeroTol)
276  {
277  return;
278  }
279 
280 
281  // Allocate array storage
282  Array<OneD, NekDouble> tmpAx_s (nGlobal, 0.0);
283  Array<OneD, NekDouble> px_s (nGlobal, 0.0);
284  Array<OneD, NekDouble> tmp1, tmp2;
285 
286  // Create NekVector wrappers for linear algebra operations
287  NekVector<NekDouble> px (nNonDir, tmp1 = px_s + nDir, eWrapper);
288  NekVector<NekDouble> tmpAx (nNonDir, tmp2 = tmpAx_s + nDir, eWrapper);
289 
290 
291  // calculating \tilda{x} - sum \alpha_i\tilda{x}_i
292 
293  Vmath::Vcopy(nNonDir,
294  tmp1 = newX + nDir, 1,
295  tmp2 = px_s + nDir, 1);
296 
297  if (m_prevLinSol.size() > 0)
298  {
299  v_DoMatrixMultiply(newX, tmpAx_s);
300  }
301 
302  Array<OneD, NekDouble> alpha (m_prevLinSol.size(), 0.0);
303  for (int i = 0; i < m_prevLinSol.size(); i++)
304  {
305  alpha[i] = Vmath::Dot2(nNonDir,
306  m_prevLinSol[i],
307  tmpAx_s + nDir,
308  m_map + nDir);
309  }
310  vComm->AllReduce(alpha, Nektar::LibUtilities::ReduceSum);
311 
312  for (int i = 0; i < m_prevLinSol.size(); i++)
313  {
314  if (alpha[i] < NekConstants::kNekZeroTol)
315  {
316  continue;
317  }
318 
319  NekVector<NekDouble> xi (nNonDir, m_prevLinSol[i], eWrapper);
320  px -= alpha[i] * xi;
321  }
322 
323 
324  // Some solutions generated by CG are identical zeros, see
325  // solutions generated for Test_Tet_equitri.xml (IncNavierStokesSolver).
326  // Not going to store identically zero solutions.
327 
328  NekDouble anorm = CalculateAnorm(nGlobal, px_s, nDir);
329  if (anorm < NekConstants::kNekZeroTol)
330  {
331  return;
332  }
333 
334  // normalisation of new solution
335  Vmath::Smul(nNonDir, 1.0/anorm, px_s.get() + nDir, 1, px_s.get() + nDir, 1);
336 
337  // updating storage with non-Dirichlet-dof part of new solution vector
338  m_prevLinSol.push_back(px_s + nDir);
339  m_numPrevSols++;
340  }
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:53
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:199
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:940
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
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 93 of file GlobalLinSysIterative.cpp.

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

99  {
100  if (m_useProjection)
101  {
102  DoAconjugateProjection(nGlobal, pInput, pOutput, plocToGloMap, nDir);
103  }
104  else
105  {
106  // applying plain Conjugate Gradient
107  DoConjugateGradient(nGlobal, pInput, pOutput, plocToGloMap, nDir);
108  }
109  }
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 96 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 93 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

Provide verbose output and 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().

bool Nektar::MultiRegions::GlobalLinSysIterative::m_verbose
protected

Definition at line 90 of file GlobalLinSysIterative.h.

Referenced by DoConjugateGradient(), and GlobalLinSysIterative().