Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
DriverArpack.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File DriverArpack.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Arnoldi solver using Arpack
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
37 
38 using namespace std;
39 
40 namespace Nektar
41 {
42 namespace SolverUtils
43 {
44 
45 std::string DriverArpack::arpackProblemTypeLookupIds[6] = {
46  LibUtilities::SessionReader::RegisterEnumValue("ArpackProblemType","LargestReal" ,0),
47  LibUtilities::SessionReader::RegisterEnumValue("ArpackProblemType","SmallestReal" ,1),
48  LibUtilities::SessionReader::RegisterEnumValue("ArpackProblemType","LargestImag" ,2),
49  LibUtilities::SessionReader::RegisterEnumValue("ArpackProblemType","SmallestImag" ,3),
50  LibUtilities::SessionReader::RegisterEnumValue("ArpackProblemType","LargestMag" ,4),
51  LibUtilities::SessionReader::RegisterEnumValue("ArpackProblemType","SmallestMag" ,5),
52 };
53 std::string DriverArpack::arpackProblemTypeDefault = LibUtilities::SessionReader::RegisterDefaultSolverInfo("ArpackProblemType","LargestMag");
54 std::string DriverArpack::driverLookupId = LibUtilities::SessionReader::RegisterEnumValue("Driver","Arpack",0);
55 
56 std::string DriverArpack::className = GetDriverFactory().RegisterCreatorFunction("Arpack", DriverArpack::create);
57 
58 std::string DriverArpack::ArpackProblemTypeTrans[6] =
59 { "LR", "SR", "LI", "SI", "LM", "SM" };
60 
61 
62 /**
63  *
64  */
65 DriverArpack::DriverArpack(const LibUtilities::SessionReaderSharedPtr pSession)
66  : DriverArnoldi(pSession)
67 {
68 }
69 
70 
71 /**
72  *
73  */
75 {
76 }
77 
78 /**
79  *
80  */
81 void DriverArpack::v_InitObject(ostream &out)
82 {
84 
85  //Initialisation of Arnoldi parameters
86  m_maxn = 1000000; // Maximum size of the problem
87  m_maxnev = 200; // maximum number of eigenvalues requested
88  m_maxncv = 500; // Largest number of basis vector used in Implicitly Restarted Arnoldi
89 
90  // Error alerts
91  ASSERTL0(m_nvec <= m_maxnev,"NEV is greater than MAXNEV");
92  ASSERTL0(m_kdim <= m_maxncv,"NEV is greater than MAXNEV");
93  ASSERTL0(2 <= m_kdim-m_nvec,"NCV-NEV is less than 2");
94 
95  m_equ[0]->PrintSummary(out);
96 
97  // Print session parameters
98  out << "\tArnoldi solver type : Arpack" << endl;
99 
100  out << "\tArpack problem type : ";
101  out << ArpackProblemTypeTrans[m_session->GetSolverInfoAsEnum<int>("ArpackProblemType")] << endl;
103 
104  m_equ[m_nequ - 1]->DoInitialise();
105 
106  // FwdTrans Initial conditions to be in Coefficient Space
107  m_equ[m_nequ-1] ->TransPhysToCoeff();
108 
109 }
110 
111 
112 void DriverArpack::v_Execute(ostream &out)
113 
114 {
115  Array<OneD, NekDouble> tmpworkd;
116 
117  int nq = m_equ[0]->UpdateFields()[0]->GetNcoeffs(); // Number of points in the mesh
118  int n = m_nfields*nq; // Number of points in eigenvalue calculation
119  int lworkl = 3*m_kdim*(m_kdim+2); // Size of work array
120  int ido ; //REVERSE COMMUNICATION parameter. At the first call must be initialised at 0
121  int info; // do not set initial vector (info=0 random initial vector, info=1 read initial vector from session file)
122 
123  int iparam[11];
124  int ipntr[14];
125 
126  Array<OneD, int> ritzSelect;
129  Array<OneD, NekDouble> workev;
131  NekDouble sigmar, sigmai;
132 
133  Array<OneD, NekDouble> resid(n);
135  Array<OneD, NekDouble> workl(lworkl, 0.0);
136  Array<OneD, NekDouble> workd(3*n, 0.0);
137 
138  ASSERTL0(n <= m_maxn, "N is greater than MAXN");
139 
140  if(m_session->DefinesFunction("InitialConditions"))
141  {
142  out << "\tInital vector : input file " << endl;
143  info = 1;
145  }
146  else
147  {
148  out << "\tInital vector : random " << endl;
149  info = 0;
150  }
151 
152  char B;
153 
154  iparam[0] = 1; // strategy for shift-invert
155  iparam[1] = 0; // (deprecated)
156  iparam[2] = m_nits; // maximum number of iterations allowed/taken
157  iparam[3] = 1; // blocksize to be used for recurrence
158  iparam[4] = 0; // number of converged ritz eigenvalues
159  iparam[5] = 0; // (deprecated)
160 
161  // Use generalized B matrix for coupled solver.
163  {
164  iparam[6] = 1; // computation mode 1=> matrix-vector prod
165  B = 'I';
166  }
167  else {
168  iparam[6] = 3; // computation mode 1=> matrix-vector prod
169  B = 'G';
170  }
171 #if 0
172  if((fabs(m_realShift) > NekConstants::kNekZeroTol)|| // use shift if m_realShift > 1e-12
174  {
175  iparam[6] = 3; // This was 3 need to know what to set it to
176  B = 'G';
177  }
178  else
179  {
180  iparam[6] = 1; // computation mode 1=> matrix-vector prod
181  B = 'I';
182  }
183 #endif
184  iparam[7] = 0; // (for shift-invert)
185  iparam[8] = 0; // number of MV operations
186  iparam[9] = 0; // number of BV operations
187  iparam[10]= 0; // number of reorthogonalisation steps
188 
189  int cycle = 0;
190  const char* problem = ArpackProblemTypeTrans[m_session->GetSolverInfoAsEnum<int>("ArpackProblemType")].c_str();
191 
192  std::string name = m_session->GetSessionName() + ".evl";
193  ofstream pFile(name.c_str());
194 
195  ido = 0; //At the first call must be initialisedat 0
196 
197  while(ido != 99)//ido==-1 || ido==1 || ido==0)
198  {
199  //Routine for eigenvalue evaluation for non-symmetric operators
200  Arpack::Dnaupd( ido, &B, // B='I' for std eval problem
201  n, problem, m_nvec,
202  m_evtol, &resid[0], m_kdim,
203  &v[0], n, iparam, ipntr, &workd[0],
204  &workl[0], lworkl, info);
205 
206  //Plotting of real and imaginary part of the
207  //eigenvalues from workl
208  out << "\rIteration " << cycle << ", output: " << info << ", ido=" << ido << " " << std::flush;
209 
210  if(!((cycle-1)%m_kdim)&&(cycle> m_kdim)&&(ido!=2))
211  {
212  pFile << "Krylov spectrum at iteration: " << cycle << endl;
213 
215  {
216  pFile << "EV Magnitude Angle Growth Frequency Residual" << endl;
217  }
218  else
219  {
220  pFile << "EV Real Imaginary inverse real inverse imag Residual" << endl;
221  }
222 
223  out << endl;
224  for(int k = 0; k < m_kdim; ++k)
225  {
226  // write m_kdim eigs to screen
227  WriteEvs(pFile,k, workl[ipntr[5]-1+k],workl[ipntr[6]-1+k]);
228  }
229  }
230 
231  if (ido == 99) break;
232 
233  switch(ido)
234  {
235  case -1:
236  case 1: // Note that ido=1 we are using input x
237  // (workd[inptr[0]-1]) rather than Mx as
238  // recommended in manual since it is not
239  // possible to impose forcing directly.
240  CopyArnoldiArrayToField(tmpworkd = workd + (ipntr[0]-1));
241 
242  m_equ[0]->TransCoeffToPhys();
243 
244  m_equ[0]->DoSolve();
246  {
247  //start Adjoint with latest fields of direct
248  CopyFwdToAdj();
249 
250  m_equ[1]->TransCoeffToPhys();
251  m_equ[1]->DoSolve();
252  }
253 
254  if(!(cycle%m_infosteps))
255  {
256  out << endl;
257  m_equ[0]->Output();
258  }
259 
260  // operated fields are copied into workd[inptr[1]-1]
261  CopyFieldToArnoldiArray(tmpworkd = workd + (ipntr[1]-1));
262 
263  cycle++;
264  break;
265  case 2: // provide y = M x (bwd trans and iproduct);
266  {
267  //workd[inptr[0]-1] copied into operator fields
268  CopyArnoldiArrayToField(tmpworkd = workd + (ipntr[0]-1));
269 
270  m_equ[0]->TransCoeffToPhys();
271 
272  Array<OneD, MultiRegions::ExpListSharedPtr> fields = m_equ[0]->UpdateFields();
273  for (int i = 0; i < fields.num_elements(); ++i)
274  {
275  fields[i]->IProductWRTBase(fields[i]->GetPhys(),
276  fields[i]->UpdateCoeffs());
277  }
278 
279  // operated fields are copied into workd[inptr[1]-1]
280  CopyFieldToArnoldiArray(tmpworkd = workd + (ipntr[1]-1));
281  break;
282  }
283  default:
284  ASSERTL0(false, "Unexpected reverse communication request.");
285  }
286 
287  }
288 
289  out<< endl << "Converged in " << iparam[8] << " iterations" << endl;
290 
291  ASSERTL0(info >= 0," Error with Dnaupd");
292 
293  ritzSelect = Array<OneD, int> (m_kdim,0);
294  dr = Array<OneD, NekDouble> (m_nvec+1,0.0);
295  di = Array<OneD, NekDouble> (m_nvec+1,0.0);
296  workev = Array<OneD, NekDouble> (3*m_kdim);
297  z = Array<OneD, NekDouble> (n*(m_nvec+1));
298 
299  if(m_negatedOp)
300  {
301  sigmar = -m_realShift;
302  }
303  else
304  {
305  sigmar = m_realShift;
306  }
307 
308  // Do not pass imaginary shift to Arpack since we have not
309  // used a Fortran complex number format and so processing
310  // is mucked up. Need to do some processing afterwards.
311  sigmai = 0;
312 
313  //Setting 'A', Ritz vectors are computed. 'S' for Shur vectors
314  Arpack::Dneupd(1, "A", ritzSelect.get(), dr.get(), di.get(),
315  z.get(), n, sigmar, sigmai, workev.get(), &B,
316  n, problem, m_nvec, m_evtol, resid.get(), m_kdim,
317  v.get(), n, iparam, ipntr, workd.get(),
318  workl.get(),lworkl,info);
319 
320  ASSERTL0(info == 0, " Error with Dneupd");
321 
322  int nconv=iparam[4];
323 
324  // Subtract off complex shift if it exists
325  if(m_negatedOp)
326  {
327  Vmath::Sadd(nconv,m_imagShift,di,1,di,1);
328  }
329  else
330  {
331  Vmath::Sadd(nconv,-m_imagShift,di,1,di,1);
332  }
333 
334  WARNINGL0(m_imagShift == 0,"Complex Shift applied. "
335  "Need to implement Ritz re-evaluation of"
336  "eigenvalue. Only one half of complex "
337  "value will be correct");
338 
339 
340  Array<OneD, MultiRegions::ExpListSharedPtr> fields = m_equ[0]->UpdateFields();
341 
342  out << "Converged Eigenvalues: " << nconv << endl;
343  pFile << "Converged Eigenvalues: " << nconv << endl;
344 
346  {
347  out << " Magnitude Angle Growth Frequency" << endl;
348  pFile << " Magnitude Angle Growth Frequency" << endl;
349  for(int i= 0; i< nconv; ++i)
350  {
351  WriteEvs(out,i,dr[i],di[i]);
352  WriteEvs(pFile,i,dr[i],di[i]);
353 
354  std::string file = m_session->GetSessionName() + "_eig_"
355  + boost::lexical_cast<std::string>(i)
356  + ".fld";
357  WriteFld(file,z + i*n);
358  }
359  }
360  else
361  {
362  out << " Real Imaginary " << endl;
363  pFile << " Real Imaginary " << endl;
364  for(int i= 0; i< nconv; ++i)
365  {
366  WriteEvs(out,i,dr[i],di[i],
368  WriteEvs(pFile,i,dr[i],di[i],
370 
371  std::string file = m_session->GetSessionName() + "_eig_"
372  + boost::lexical_cast<std::string>(i)
373  + ".fld";
374  WriteFld(file,z + i*n);
375  }
376  }
377 
378  m_real_evl = dr;
379  m_imag_evl = di;
380 
381  pFile.close();
382 
383  for(int j = 0; j < m_equ[0]->GetNvariables(); ++j)
384  {
385  NekDouble vL2Error = m_equ[0]->L2Error(j,false);
386  NekDouble vLinfError = m_equ[0]->LinfError(j);
387  if (m_comm->GetRank() == 0)
388  {
389  out << "L 2 error (variable " << m_equ[0]->GetVariable(j) << ") : " << vL2Error << endl;
390  out << "L inf error (variable " << m_equ[0]->GetVariable(j) << ") : " << vLinfError << endl;
391  }
392  }
393 }
394 
395 }
396 }
Array< OneD, NekDouble > m_imag_evl
Definition: DriverArnoldi.h:70
static void Dneupd(const int &rvec, const char *howmny, const int *select, double *dr, double *di, double *z, const int &ldz, const double &sigmar, const double &sigmai, double *workev, const char *bmat, const int &n, const char *which, const int &nev, const double &tol, double *resid, const int &ncv, double *v, const int &ldv, int *iparam, int *ipntr, double *workd, double *workl, const int &lworkl, int &info)
Post-processing routine to computed eigenvector of computed eigenvalues in Dnaupd.
Definition: Arpack.hpp:149
virtual ~DriverArpack()
Destructor.
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
void CopyFwdToAdj()
Copy the forward field to the adjoint system in transient growth calculations.
NekDouble m_evtol
Maxmum number of iterations.
Definition: DriverArnoldi.h:58
void CopyArnoldiArrayToField(Array< OneD, NekDouble > &array)
Copy Arnoldi storage to fields.
virtual void v_Execute(std::ostream &out=std::cout)
Virtual function for solve implementation.
SOLVER_UTILS_EXPORT void ArnoldiSummary(std::ostream &out)
STL namespace.
void CopyFieldToArnoldiArray(Array< OneD, NekDouble > &array)
Copy fields to Arnoldi storage.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
void WriteFld(std::string file, std::vector< Array< OneD, NekDouble > > coeffs)
Write coefficients to file.
virtual void v_InitObject(std::ostream &out=std::cout)
Virtual function for initialisation implementation.
static const NekDouble kNekZeroTol
int m_nvec
Dimension of Krylov subspace.
Definition: DriverArnoldi.h:56
Array< OneD, NekDouble > m_real_evl
Operator in solve call is negated.
Definition: DriverArnoldi.h:69
void WriteEvs(std::ostream &evlout, const int k, const NekDouble real, const NekDouble imag, NekDouble resid=NekConstants::kNekUnsetDouble, bool DumpInverse=true)
#define WARNINGL0(condition, msg)
Definition: ErrorUtil.hpp:167
virtual void v_InitObject(std::ostream &out=std::cout)
static void Dnaupd(int &ido, const char *bmat, const int &n, const char *which, const int &nev, const double &tol, double *resid, const int &ncv, double *v, const int &ldv, int *iparam, int *ipntr, double *workd, double *workl, const int &lworkl, int &info)
Top level reverse communication interface to solve real double-precision non-symmetric problems...
Definition: Arpack.hpp:135
enum EvolutionOperatorType m_EvolutionOperator
Evolution Operator.
Definition: Driver.h:100
double NekDouble
static const NekDouble kNekUnsetDouble
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
Definition: Vmath.cpp:301
int m_nfields
interval to dump information if required.
Definition: DriverArnoldi.h:64
bool m_timeSteppingAlgorithm
Period of time stepping algorithm.
Definition: DriverArnoldi.h:60
static std::string ArpackProblemTypeTrans[]
Definition: DriverArpack.h:88
Array< OneD, EquationSystemSharedPtr > m_equ
Equation system to solve.
Definition: Driver.h:94
int m_infosteps
underlying operator is time stepping
Definition: DriverArnoldi.h:62
DriverFactory & GetDriverFactory()
Definition: Driver.cpp:66
LibUtilities::CommSharedPtr m_comm
Communication object.
Definition: Driver.h:85
Base class for the development of solvers.
Definition: DriverArnoldi.h:47
LibUtilities::SessionReaderSharedPtr m_session
Session reader object.
Definition: Driver.h:88
int m_nequ
number of equations
Definition: Driver.h:97
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215
int m_nits
Number of vectors to test.
Definition: DriverArnoldi.h:57