Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
DriverModifiedArnoldi.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File DriverModifiedArnoldi.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: Driver for eigenvalue analysis using the modified Arnoldi
33 // method.
34 //
35 ///////////////////////////////////////////////////////////////////////////////
36 
37 #include <iomanip>
38 
39 #include <SolverUtils/DriverModifiedArnoldi.h>
40 
41 namespace Nektar
42 {
43  namespace SolverUtils
44  {
45  string DriverModifiedArnoldi::className = GetDriverFactory().RegisterCreatorFunction("ModifiedArnoldi", DriverModifiedArnoldi::create);
46  string DriverModifiedArnoldi::driverLookupId = LibUtilities::SessionReader::RegisterEnumValue("Driver","ModifiedArnoldi",0);
47 
48  /**
49  *
50  */
51  DriverModifiedArnoldi::DriverModifiedArnoldi(
52  const LibUtilities::SessionReaderSharedPtr pSession)
53  : DriverArnoldi(pSession)
54  {
55  }
56 
57 
58  /**
59  *
60  */
61  DriverModifiedArnoldi::~DriverModifiedArnoldi()
62  {
63  }
64 
65 
66  /**
67  *
68  */
69  void DriverModifiedArnoldi::v_InitObject(ostream &out)
70  {
71  DriverArnoldi::v_InitObject(out);
72 
73  m_equ[0]->PrintSummary(out);
74 
75  // Print session parameters
76  out << "\tArnoldi solver type : Modified Arnoldi" << endl;
77 
78  DriverArnoldi::ArnoldiSummary(out);
79 
80  m_equ[m_nequ - 1]->DoInitialise();
81 
82  //FwdTrans Initial conditions to be in Coefficient Space
83  m_equ[m_nequ-1] ->TransPhysToCoeff();
84 
85  }
86 
87 
88  /**
89  *
90  */
91  void DriverModifiedArnoldi::v_Execute(ostream &out)
92  {
93  int i = 0;
94  int j = 0;
95  int nq = m_equ[0]->UpdateFields()[0]->GetNcoeffs();
96  int ntot = m_nfields*nq;
97  int converged = 0;
98  NekDouble resnorm = 0.0;
99  std::string evlFile = m_session->GetFilename().substr(0,m_session->GetFilename().find_last_of('.')) + ".evl";
100  ofstream evlout(evlFile.c_str());
101 
102  // Allocate memory
103  Array<OneD, NekDouble> alpha = Array<OneD, NekDouble> (m_kdim+1, 0.0);
104  Array<OneD, NekDouble> wr = Array<OneD, NekDouble> (m_kdim, 0.0);
105  Array<OneD, NekDouble> wi = Array<OneD, NekDouble> (m_kdim, 0.0);
106  Array<OneD, NekDouble> zvec = Array<OneD, NekDouble> (m_kdim*m_kdim, 0.0);
107 
108  Array<OneD, Array<OneD, NekDouble> > Kseq
109  = Array<OneD, Array<OneD, NekDouble> > (m_kdim + 1);
110  Array<OneD, Array<OneD, NekDouble> > Tseq
111  = Array<OneD, Array<OneD, NekDouble> > (m_kdim + 1);
112  for (i = 0; i < m_kdim + 1; ++i)
113  {
114  Kseq[i] = Array<OneD, NekDouble>(ntot, 0.0);
115  Tseq[i] = Array<OneD, NekDouble>(ntot, 0.0);
116  }
117 
118 
119  // Copy starting vector into second sequence element (temporary).
120  if(m_session->DefinesFunction("InitialConditions"))
121  {
122  out << "\tInital vector : specified in input file " << endl;
123  m_equ[0]->SetInitialConditions(0.0,false);
124 
125  CopyFieldToArnoldiArray(Kseq[1]);
126  }
127  else
128  {
129  out << "\tInital vector : random " << endl;
130 
131  NekDouble eps=1;
132 
133  Vmath::FillWhiteNoise(ntot, eps , &Kseq[1][0], 1);
134 
135  }
136 
137  // Perform one iteration to enforce boundary conditions.
138  // Set this as the initial value in the sequence.
139  EV_update(Kseq[1], Kseq[0]);
140  out << "Iteration: " << 0 << endl;
141 
142  // Normalise first vector in sequence
143  alpha[0] = std::sqrt(Blas::Ddot(ntot, &Kseq[0][0], 1,
144  &Kseq[0][0], 1));
145  m_comm->AllReduce(alpha[0], Nektar::LibUtilities::ReduceSum);
146 
147  //alpha[0] = std::sqrt(alpha[0]);
148  Vmath::Smul(ntot, 1.0/alpha[0], Kseq[0], 1, Kseq[0], 1);
149 
150  // Fill initial krylov sequence
151  NekDouble resid0;
152  for (i = 1; !converged && i <= m_kdim; ++i)
153  {
154  // Compute next vector
155  EV_update(Kseq[i-1], Kseq[i]);
156 
157  // Normalise
158  alpha[i] = std::sqrt(Blas::Ddot(ntot, &Kseq[i][0], 1,
159  &Kseq[i][0], 1));
160  m_comm->AllReduce(alpha[i], Nektar::LibUtilities::ReduceSum);
161 
162 
163  //alpha[i] = std::sqrt(alpha[i]);
164  Vmath::Smul(ntot, 1.0/alpha[i], Kseq[i], 1, Kseq[i], 1);
165 
166  // Copy Krylov sequence into temporary storage
167  for (int k = 0; k < i + 1; ++k)
168  {
169  Vmath::Vcopy(ntot, Kseq[k], 1, Tseq[k], 1);
170  }
171 
172  // Generate Hessenberg matrix and compute eigenvalues of it.
173  EV_small(Tseq, ntot, alpha, i, zvec, wr, wi, resnorm);
174 
175  // Test for convergence.
176  converged = EV_test(i,i,zvec,wr,wi,resnorm,std::min(i,m_nvec),evlout,resid0);
177  converged = max (converged, 0);
178  out << "Iteration: " << i << " (residual : " << resid0 << ")" <<endl;
179  }
180 
181  // Continue with full sequence
182  if (!converged)
183  {
184  for (i = m_kdim + 1; !converged && i <= m_nits; ++i)
185  {
186  // Shift all the vectors in the sequence.
187  // First vector is removed.
188  //NekDouble invnorm = 1.0/sqrt(Blas::Ddot(ntot,Kseq[1],1,Kseq[1],1));
189  for (int j = 1; j <= m_kdim; ++j)
190  {
191  alpha[j-1] = alpha[j];
192  //Vmath::Smul(ntot,invnorm,Kseq[j],1,Kseq[j],1);
193  Vmath::Vcopy(ntot, Kseq[j], 1, Kseq[j-1], 1);
194  }
195 
196  // Compute next vector
197  EV_update(Kseq[m_kdim - 1], Kseq[m_kdim]);
198 
199  // Compute new scale factor
200  alpha[m_kdim] = std::sqrt(Vmath::Dot(ntot, &Kseq[m_kdim][0], 1, &Kseq[m_kdim][0], 1));
201  //alpha[m_kdim] = std::sqrt(alpha[m_kdim]);
202  Vmath::Smul(ntot, 1.0/alpha[m_kdim], Kseq[m_kdim], 1, Kseq[m_kdim], 1);
203 
204  // Copy Krylov sequence into temporary storage
205  for (int k = 0; k < m_kdim + 1; ++k)
206  {
207  Vmath::Vcopy(ntot, Kseq[k], 1, Tseq[k], 1);
208  }
209 
210  // Generate Hessenberg matrix and compute eigenvalues of it
211  EV_small(Tseq, ntot, alpha, m_kdim, zvec, wr, wi, resnorm);
212 
213  // Test for convergence.
214  converged = EV_test(i,m_kdim,zvec,wr,wi,resnorm,m_nvec,evlout,resid0);
215  out << "Iteration: " << i << " (residual : " << resid0 << ")" <<endl;
216  }
217  }
218 
219  m_equ[0]->Output();
220 
221  // Evaluate and output computation time and solution accuracy.
222  // The specific format of the error output is essential for the
223  // regression tests to work.
224  // Evaluate L2 Error
225  for(j = 0; j < m_equ[0]->GetNvariables(); ++j)
226  {
227  NekDouble vL2Error = m_equ[0]->L2Error(j,false);
228  NekDouble vLinfError = m_equ[0]->LinfError(j);
229  if (m_comm->GetRank() == 0)
230  {
231  out << "L 2 error (variable " << m_equ[0]->GetVariable(j) << ") : " << vL2Error << endl;
232  out << "L inf error (variable " << m_equ[0]->GetVariable(j) << ") : " << vLinfError << endl;
233  }
234  }
235 
236  // Process eigenvectors and write out.
237  EV_post(Tseq, Kseq, ntot, min(--i, m_kdim), m_nvec, zvec, wr, wi, converged);
238 
239  // store eigenvalues so they can be access from driver class
240  m_real_evl = wr;
241  m_imag_evl = wi;
242 
243  // Close the runtime info file.
244  evlout.close();
245  }
246 
247 
248  /**
249  *
250  */
251  void DriverModifiedArnoldi::EV_update(
252  Array<OneD, NekDouble> &src,
253  Array<OneD, NekDouble> &tgt)
254  {
255  // Copy starting vector into first sequence element.
256  CopyArnoldiArrayToField(src);
257 
258  m_equ[0]->TransCoeffToPhys();
259 
260 
261  m_equ[0]->DoSolve();
262 
263  if(m_EvolutionOperator == eTransientGrowth)
264  {
265  Array<OneD, MultiRegions::ExpListSharedPtr> fields;
266  fields = m_equ[0]->UpdateFields();
267 
268  //start Adjoint with latest fields of direct
269  CopyFwdToAdj();
270 
271  m_equ[1]->TransCoeffToPhys();
272 
273  m_equ[1]->DoSolve();
274  }
275 
276  // Copy starting vector into first sequence element.
277  CopyFieldToArnoldiArray(tgt);
278  }
279 
280 
281  /**
282  *
283  */
284  void DriverModifiedArnoldi::EV_small(
285  Array<OneD, Array<OneD, NekDouble> > &Kseq,
286  const int ntot,
287  const Array<OneD, NekDouble> &alpha,
288  const int kdim,
289  Array<OneD, NekDouble> &zvec,
290  Array<OneD, NekDouble> &wr,
291  Array<OneD, NekDouble> &wi,
292  NekDouble &resnorm)
293  {
294  int kdimp = kdim + 1;
295  int lwork = 10*kdim;
296  int ier;
297  Array<OneD, NekDouble> R(kdimp * kdimp, 0.0);
298  Array<OneD, NekDouble> H(kdimp * kdim, 0.0);
299  Array<OneD, NekDouble> rwork(lwork, 0.0);
300 
301  // Modified G-S orthonormalisation
302  for (int i = 0; i < kdimp; ++i)
303  {
304  NekDouble gsc = std::sqrt(Vmath::Dot(ntot, &Kseq[i][0], 1, &Kseq[i][0], 1));
305  ASSERTL0(gsc != 0.0, "Vectors are linearly independent.");
306  R[i*kdimp+i] = gsc;
307  Vmath::Smul(ntot, 1.0/gsc, Kseq[i], 1, Kseq[i], 1);
308  for (int j = i + 1; j < kdimp; ++j)
309  {
310  gsc = Vmath::Dot(ntot, Kseq[i], 1, Kseq[j], 1);
311  Vmath::Svtvp(ntot, -gsc, Kseq[i], 1, Kseq[j], 1, Kseq[j], 1);
312  R[j*kdimp+i] = gsc;
313  }
314  }
315 
316  // Compute H matrix
317  for (int i = 0; i < kdim; ++i)
318  {
319  for (int j = 0; j < kdim; ++j)
320  {
321  H[j*kdim+i] = alpha[j+1] * R[(j+1)*kdimp+i]
322  - Vmath::Dot(j, &H[0] + i, kdim, &R[0] + j*kdimp, 1);
323  H[j*kdim+i] /= R[j*kdimp+j];
324  }
325  }
326 
327  H[(kdim-1)*kdim+kdim] = alpha[kdim]
328  * std::fabs(R[kdim*kdimp+kdim] / R[(kdim-1)*kdimp + kdim-1]);
329 
330  Lapack::dgeev_('N','V',kdim,&H[0],kdim,&wr[0],&wi[0],0,1,&zvec[0],kdim,&rwork[0],lwork,ier);
331 
332  ASSERTL0(!ier, "Error with dgeev");
333 
334  resnorm = H[(kdim-1)*kdim + kdim];
335  }
336 
337 
338  /**
339  *
340  */
341  int DriverModifiedArnoldi::EV_test(
342  const int itrn,
343  const int kdim,
344  Array<OneD, NekDouble> &zvec,
345  Array<OneD, NekDouble> &wr,
346  Array<OneD, NekDouble> &wi,
347  const NekDouble resnorm,
348  const int nvec,
349  ofstream &evlout,
350  NekDouble &resid0)
351  {
352  int idone = 0;
353  // NekDouble period = 0.1;
354 
355  Array<OneD, NekDouble> resid(kdim);
356  for (int i = 0; i < kdim; ++i)
357  {
358  resid[i] = resnorm * std::fabs(zvec[kdim - 1 + i*kdim]) /
359  std::sqrt(Vmath::Dot(kdim, &zvec[0] + i*kdim, 1, &zvec[0] + i*kdim, 1));
360  if (wi[i] < 0.0) resid[i-1] = resid[i] = hypot(resid[i-1], resid[i]);
361  }
362  EV_sort(zvec, wr, wi, resid, kdim);
363 
364  if (resid[nvec-1] < m_evtol) idone = nvec;
365 
366  evlout << "-- Iteration = " << itrn << ", H(k+1, k) = " << resnorm << endl;
367 
368  evlout.precision(4);
369  evlout.setf(ios::scientific, ios::floatfield);
370  if(m_timeSteppingAlgorithm)
371  {
372  evlout << "EV Magnitude Angle Growth Frequency Residual"
373  << endl;
374  }
375  else
376  {
377  evlout << "EV Real Imaginary inverse real inverse imag Residual"
378  << endl;
379  }
380 
381  for (int i = 0; i < kdim; i++)
382  {
383  WriteEvs(evlout,i,wr[i],wi[i],resid[i]);
384  }
385 
386 
387  resid0 = resid[nvec-1];
388  return idone;
389  }
390 
391 
392  /**
393  *
394  */
395  void DriverModifiedArnoldi::EV_sort(
396  Array<OneD, NekDouble> &evec,
397  Array<OneD, NekDouble> &wr,
398  Array<OneD, NekDouble> &wi,
399  Array<OneD, NekDouble> &test,
400  const int dim)
401  {
402  Array<OneD, NekDouble> z_tmp(dim,0.0);
403  NekDouble wr_tmp, wi_tmp, te_tmp;
404  for (int j = 1; j < dim; ++j)
405  {
406  wr_tmp = wr[j];
407  wi_tmp = wi[j];
408  te_tmp = test[j];
409  Vmath::Vcopy(dim, &evec[0] + j*dim, 1, &z_tmp[0], 1);
410  int i = j - 1;
411  while (i >= 0 && test[i] > te_tmp)
412  {
413  wr[i+1] = wr[i];
414  wi[i+1] = wi[i];
415  test[i+1] = test[i];
416  Vmath::Vcopy(dim, &evec[0] + i*dim, 1, &evec[0] + (i+1)*dim, 1);
417  i--;
418  }
419  wr[i+1] = wr_tmp;
420  wi[i+1] = wi_tmp;
421  test[i+1] = te_tmp;
422  Vmath::Vcopy(dim, &z_tmp[0], 1, &evec[0] + (i+1)*dim, 1);
423  }
424  };
425 
426 
427  /**
428  *
429  */
430  void DriverModifiedArnoldi::EV_post(
431  Array<OneD, Array<OneD, NekDouble> > &Tseq,
432  Array<OneD, Array<OneD, NekDouble> > &Kseq,
433  const int ntot,
434  const int kdim,
435  const int nvec,
436  Array<OneD, NekDouble> &zvec,
437  Array<OneD, NekDouble> &wr,
438  Array<OneD, NekDouble> &wi,
439  const int icon)
440  {
441  if (icon == 0)
442  {
443  // Not converged, write final Krylov vector
444  ASSERTL0(false, "Convergence was not achieved within the prescribed number of iterations.");
445  }
446  else if (icon < 0)
447  {
448  // Minimum residual reached
449  ASSERTL0(false, "Minimum residual reached.");
450  }
451  else if (icon == nvec)
452  {
453  // Converged, write out eigenvectors
454  EV_big(Tseq, Kseq, ntot, kdim, icon, zvec, wr, wi);
455  Array<OneD, MultiRegions::ExpListSharedPtr> fields
456  = m_equ[0]->UpdateFields();
457 
458  for (int j = 0; j < icon; ++j)
459  {
460  std::string file = m_session->GetFilename().substr(0,m_session->GetFilename().find_last_of('.')) + "_eig_" + boost::lexical_cast<std::string>(j);
461 
462  WriteFld(file,Kseq[j]);
463  }
464  }
465  else
466  {
467  // Not recognised
468  ASSERTL0(false, "Unrecognised value.");
469  }
470  }
471 
472 
473  /**
474  *
475  */
476  void DriverModifiedArnoldi::EV_big(
477  Array<OneD, Array<OneD, NekDouble> > &bvecs,
478  Array<OneD, Array<OneD, NekDouble> > &evecs,
479  const int ntot,
480  const int kdim,
481  const int nvec,
482  Array<OneD, NekDouble> &zvec,
483  Array<OneD, NekDouble> &wr,
484  Array<OneD, NekDouble> &wi)
485  {
486  NekDouble wgt, norm;
487 
488  // Generate big eigenvectors
489  for (int j = 0; j < nvec; ++j)
490  {
491  Vmath::Zero(ntot, evecs[j], 1);
492  for (int i = 0; i < kdim; ++i)
493  {
494  wgt = zvec[i + j*kdim];
495  Vmath::Svtvp(ntot, wgt, bvecs[i], 1, evecs[j], 1, evecs[j], 1);
496  }
497  }
498 
499  // Normalise the big eigenvectors
500  for (int i = 0; i < nvec; ++i)
501  {
502  if (wi[i] == 0.0) // Real mode
503  {
504  norm = std::sqrt(Vmath::Dot(ntot, evecs[i], evecs[i]));
505  Vmath::Smul(ntot, 1.0/norm, evecs[i], 1, evecs[i], 1);
506  }
507  else
508  {
509  norm = Vmath::Dot(ntot, evecs[i], 1, evecs[i], 1);
510  norm += Vmath::Dot(ntot, evecs[i+1], 1, evecs[i+1], 1);
511  norm = std::sqrt(norm);
512  Vmath::Smul(ntot, 1.0/norm, evecs[i], 1, evecs[i], 1);
513  Vmath::Smul(ntot, 1.0/norm, evecs[i+1], 1, evecs[i+1], 1);
514  i++;
515  }
516  }
517  }
518  }
519 }
Generated on Sun Mar 15 2015 00:11:30 for Nektar++ by   doxygen 1.8.1.2