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
 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: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 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: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 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: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 355 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, Nektar::MultiRegions::GlobalLinSys::m_verbose, Nektar::LibUtilities::ReduceSum, 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  // If input residual is less than tolerance skip solve.
423  {
424  if (m_verbose && m_root)
425  {
426  cout << "CG iterations made = " << m_totalIterations
427  << " using tolerance of " << m_tolerance
428  << " (error = " << sqrt(eps/m_rhs_magnitude)
429  << ", rhs_mag = " << sqrt(m_rhs_magnitude) << ")"
430  << endl;
431  }
432  return;
433  }
434 
435  m_totalIterations = 1;
436  m_precon->DoPreconditioner(r_A, tmp = w_A + nDir);
437 
438  v_DoMatrixMultiply(w_A, s_A);
439 
440  k = 0;
441 
442  vExchange[0] = Vmath::Dot2(nNonDir,
443  r_A,
444  w_A + nDir,
445  m_map + nDir);
446 
447  vExchange[1] = Vmath::Dot2(nNonDir,
448  s_A + nDir,
449  w_A + nDir,
450  m_map + nDir);
451 
452  vComm->AllReduce(vExchange, Nektar::LibUtilities::ReduceSum);
453 
454  rho = vExchange[0];
455  mu = vExchange[1];
456  min_resid = m_rhs_magnitude;
457  beta = 0.0;
458  alpha = rho/mu;
459 
460  // Continue until convergence
461  while (true)
462  {
463  if(k >= m_maxiter)
464  {
465  if (m_root)
466  {
467  cout << "CG iterations made = " << m_totalIterations
468  << " using tolerance of " << m_tolerance
469  << " (error = " << sqrt(eps/m_rhs_magnitude)
470  << ", rhs_mag = " << sqrt(m_rhs_magnitude) << ")"
471  << endl;
472  }
473  ASSERTL0(false,
474  "Exceeded maximum number of iterations");
475  }
476 
477  // Compute new search direction p_k, q_k
478  Vmath::Svtvp(nNonDir, beta, &p_A[0], 1, &w_A[nDir], 1, &p_A[0], 1);
479  Vmath::Svtvp(nNonDir, beta, &q_A[0], 1, &s_A[nDir], 1, &q_A[0], 1);
480 
481  // Update solution x_{k+1}
482  Vmath::Svtvp(nNonDir, alpha, &p_A[0], 1, &pOutput[nDir], 1, &pOutput[nDir], 1);
483 
484  // Update residual vector r_{k+1}
485  Vmath::Svtvp(nNonDir, -alpha, &q_A[0], 1, &r_A[0], 1, &r_A[0], 1);
486 
487  // Apply preconditioner
488  m_precon->DoPreconditioner(r_A, tmp = w_A + nDir);
489 
490  // Perform the method-specific matrix-vector multiply operation.
491  v_DoMatrixMultiply(w_A, s_A);
492 
493  // <r_{k+1}, w_{k+1}>
494  vExchange[0] = Vmath::Dot2(nNonDir,
495  r_A,
496  w_A + nDir,
497  m_map + nDir);
498  // <s_{k+1}, w_{k+1}>
499  vExchange[1] = Vmath::Dot2(nNonDir,
500  s_A + nDir,
501  w_A + nDir,
502  m_map + nDir);
503 
504  // <r_{k+1}, r_{k+1}>
505  vExchange[2] = Vmath::Dot2(nNonDir,
506  r_A,
507  r_A,
508  m_map + nDir);
509 
510  // Perform inner-product exchanges
511  vComm->AllReduce(vExchange, Nektar::LibUtilities::ReduceSum);
512 
513  rho_new = vExchange[0];
514  mu = vExchange[1];
515  eps = vExchange[2];
516 
518  // test if norm is within tolerance
519  if (eps < m_tolerance * m_tolerance * m_rhs_magnitude)
520  {
521  if (m_verbose && m_root)
522  {
523  cout << "CG iterations made = " << m_totalIterations
524  << " using tolerance of " << m_tolerance
525  << " (error = " << sqrt(eps/m_rhs_magnitude)
526  << ", rhs_mag = " << sqrt(m_rhs_magnitude) << ")"
527  << endl;
528  }
529  break;
530  }
531  min_resid = min(min_resid, eps);
532 
533  // Compute search direction and solution coefficients
534  beta = rho_new/rho;
535  alpha = rho_new/(mu - rho_new*beta/alpha);
536  rho = rho_new;
537  k++;
538  }
539  }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
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: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 541 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().

543  {
544  Array<OneD, NekDouble> vExchange(1);
545  vExchange[0] = Vmath::Dot(pIn.GetDimension(),&pIn[0],&pIn[0]);
546 
547  m_expList.lock()->GetComm()->GetRowComm()->AllReduce(
549 
550  // To ensure that very different rhs values are not being
551  // used in subsequent solvers such as the velocit solve in
552  // INC NS. If this works we then need to work out a better
553  // way to control this.
554  NekDouble new_rhs_mag = (vExchange[0] > 1e-6)? vExchange[0] : 1.0;
555 
557  {
558  m_rhs_magnitude = new_rhs_mag;
559  }
560  else
561  {
563  (1.0-m_rhs_mag_sm)*new_rhs_mag);
564  }
565  }
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
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 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: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 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().