Nektar++
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// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description: Driver for eigenvalue analysis using the modified Arnoldi
32// method.
33//
34///////////////////////////////////////////////////////////////////////////////
35
36#include <iomanip>
37
40
41namespace Nektar::SolverUtils
42{
43
45 GetDriverFactory().RegisterCreatorFunction("ModifiedArnoldi",
48 LibUtilities::SessionReader::RegisterEnumValue("Driver", "ModifiedArnoldi",
49 0);
50
51/**
52 *
53 */
57 : DriverArnoldi(pSession, pGraph)
58{
59}
60
61/**
62 *
63 */
65{
67
68 m_equ[0]->PrintSummary(out);
69
70 // Print session parameters
71 if (m_comm->GetRank() == 0)
72 {
73 out << "\tArnoldi solver type : Modified Arnoldi" << std::endl;
74 }
75
77
78 for (int i = 0; i < m_nequ; ++i)
79 {
80 m_equ[i]->DoInitialise();
81 }
82
83 // FwdTrans Initial conditions to be in Coefficient Space
84 m_equ[m_nequ - 1]->TransPhysToCoeff();
85}
86
87/**
88 *
89 */
90void DriverModifiedArnoldi::v_Execute(std::ostream &out)
91{
92 int i = 0;
93 int j = 0;
94 int nm = m_equ[0]->UpdateFields()[0]->GetNcoeffs();
95 int ntot = m_nfields * nm;
96 int converged = 0;
97 NekDouble resnorm = 0.0;
98 std::ofstream evlout;
99 std::string evlFile = m_session->GetSessionName() + ".evl";
100
101 if (m_comm->GetRank() == 0)
102 {
103 evlout.open(evlFile.c_str());
104 }
105
106 // Allocate memory
111
118 for (i = 0; i < m_kdim + 1; ++i)
119 {
120 Kseq[i] = Array<OneD, NekDouble>(ntot, 0.0);
121 if (m_useMask)
122 {
123 Kseqcopy[i] = Array<OneD, NekDouble>(ntot, 0.0);
124 }
125 else
126 {
127 Kseqcopy[i] = Kseq[i];
128 }
129 Tseq[i] = Array<OneD, NekDouble>(ntot, 0.0);
130 }
131
132 // Copy starting vector into second sequence element (temporary).
133 if (m_session->DefinesFunction("InitialConditions"))
134 {
135 if (m_comm->GetRank() == 0)
136 {
137 out << "\tInital vector : specified in input file "
138 << std::endl;
139 }
140 m_equ[0]->SetInitialConditions(0.0, false);
141
143 }
144 else
145 {
146 if (m_comm->GetRank() == 0)
147 {
148 out << "\tInital vector : random " << std::endl;
149 }
150
151 NekDouble eps = 0.0001;
152 Vmath::FillWhiteNoise(ntot, eps, &Kseq[1][0], 1);
153
154 auto contfield = std::dynamic_pointer_cast<MultiRegions::ContField>(
155 m_equ[0]->UpdateFields()[0]);
156 if (contfield)
157 {
158 int nGlobal =
159 contfield->GetLocalToGlobalMap()->GetNumGlobalCoeffs();
160 Array<OneD, NekDouble> global(nGlobal), tmp;
161 for (int n = 0; n < m_nfields; ++n)
162 {
163 contfield->LocalToGlobal(Kseq[1] + n * nm, global);
164 contfield->GlobalToLocal(global, tmp = Kseq[1] + n * nm);
165 }
166 }
167 if (m_useMask)
168 {
169 Vmath::Vmul(ntot, Kseq[1], 1, m_maskCoeffs, 1, Kseq[1], 1);
170 }
171 }
172
173 // Perform one iteration to enforce boundary conditions.
174 // Set this as the initial value in the sequence.
175 EV_update(Kseq[1], Kseq[0]);
176 if (m_comm->GetRank() == 0)
177 {
178 out << "Iteration: " << 0 << std::endl;
179 }
180
181 // Normalise first vector in sequence
182 if (m_useMask)
183 {
184 Vmath::Vmul(ntot, Kseq[0], 1, m_maskCoeffs, 1, Kseqcopy[0], 1);
185 }
186 alpha[0] = Blas::Ddot(ntot, &Kseqcopy[0][0], 1, &Kseqcopy[0][0], 1);
187 m_comm->AllReduce(alpha[0], Nektar::LibUtilities::ReduceSum);
188 alpha[0] = std::sqrt(alpha[0]);
189 Vmath::Smul(ntot, 1.0 / alpha[0], Kseq[0], 1, Kseq[0], 1);
190
191 // Fill initial krylov sequence
192 NekDouble resid0;
193 for (i = 1; !converged && i <= m_kdim; ++i)
194 {
195 // Compute next vector
196 EV_update(Kseq[i - 1], Kseq[i]);
197
198 // Normalise
199 if (m_useMask)
200 {
201 Vmath::Vmul(ntot, Kseq[i], 1, m_maskCoeffs, 1, Kseqcopy[i], 1);
202 }
203 alpha[i] = Blas::Ddot(ntot, &Kseqcopy[i][0], 1, &Kseqcopy[i][0], 1);
204 m_comm->AllReduce(alpha[i], Nektar::LibUtilities::ReduceSum);
205 alpha[i] = std::sqrt(alpha[i]);
206
207 // alpha[i] = std::sqrt(alpha[i]);
208 Vmath::Smul(ntot, 1.0 / alpha[i], Kseq[i], 1, Kseq[i], 1);
209
210 // Copy Krylov sequence into temporary storage
211 for (int k = 0; k < i + 1; ++k)
212 {
213 if (m_useMask)
214 {
215 Vmath::Vmul(ntot, Kseq[k], 1, m_maskCoeffs, 1, Tseq[k], 1);
216 Vmath::Vcopy(ntot, Kseq[k], 1, Kseqcopy[k], 1);
217 }
218 else
219 {
220 Vmath::Vcopy(ntot, Kseq[k], 1, Tseq[k], 1);
221 ;
222 }
223 }
224
225 // Generate Hessenberg matrix and compute eigenvalues of it.
226 EV_small(Tseq, Kseqcopy, ntot, alpha, i, zvec, wr, wi, resnorm);
227
228 // Test for convergence.
229 converged = EV_test(i, i, zvec, wr, wi, resnorm, std::min(i, m_nvec),
230 evlout, resid0);
231 if (i >= m_nvec)
232 {
233 converged = std::max(converged, 0);
234 }
235 else
236 {
237 converged = 0;
238 }
239
240 if (m_comm->GetRank() == 0)
241 {
242 out << "Iteration: " << i << " (residual : " << resid0 << ")"
243 << std::endl;
244 }
245 }
246
247 // Continue with full sequence
248 if (!converged)
249 {
250 for (i = m_kdim + 1; !converged && i <= m_nits; ++i)
251 {
252 // Shift all the vectors in the sequence.
253 // First vector is removed.
254 for (int j = 1; j <= m_kdim; ++j)
255 {
256 alpha[j - 1] = alpha[j];
257 Vmath::Vcopy(ntot, Kseq[j], 1, Kseq[j - 1], 1);
258 }
259
260 // Compute next vector
261 EV_update(Kseq[m_kdim - 1], Kseq[m_kdim]);
262
263 // Compute new scale factor
264 if (m_useMask)
265 {
266 Vmath::Vmul(ntot, Kseq[m_kdim], 1, m_maskCoeffs, 1,
267 Kseqcopy[m_kdim], 1);
268 }
269 alpha[m_kdim] = Blas::Ddot(ntot, &Kseqcopy[m_kdim][0], 1,
270 &Kseqcopy[m_kdim][0], 1);
272 alpha[m_kdim] = std::sqrt(alpha[m_kdim]);
273 Vmath::Smul(ntot, 1.0 / alpha[m_kdim], Kseq[m_kdim], 1,
274 Kseq[m_kdim], 1);
275
276 // Copy Krylov sequence into temporary storage
277 for (int k = 0; k < m_kdim + 1; ++k)
278 {
279 if (m_useMask)
280 {
281 Vmath::Vmul(ntot, Kseq[k], 1, m_maskCoeffs, 1, Tseq[k], 1);
282 Vmath::Vcopy(ntot, Kseq[k], 1, Kseqcopy[k], 1);
283 }
284 else
285 {
286 Vmath::Vcopy(ntot, Kseq[k], 1, Tseq[k], 1);
287 ;
288 }
289 }
290
291 // Generate Hessenberg matrix and compute eigenvalues of it
292 EV_small(Tseq, Kseqcopy, ntot, alpha, m_kdim, zvec, wr, wi,
293 resnorm);
294
295 // Test for convergence.
296 converged = EV_test(i, m_kdim, zvec, wr, wi, resnorm, m_nvec,
297 evlout, resid0);
298
299 if (m_comm->GetRank() == 0)
300 {
301 out << "Iteration: " << i << " (residual : " << resid0 << ")"
302 << std::endl;
303 }
304 }
305 }
306
307 m_equ[0]->Output();
308
309 // Evaluate and output computation time and solution accuracy.
310 // The specific format of the error output is essential for the
311 // regression tests to work.
312 // Evaluate L2 Error
313 for (j = 0; j < m_equ[0]->GetNvariables(); ++j)
314 {
315 NekDouble vL2Error = m_equ[0]->L2Error(j, false);
316 NekDouble vLinfError = m_equ[0]->LinfError(j);
317 if (m_comm->GetRank() == 0)
318 {
319 out << "L 2 error (variable " << m_equ[0]->GetVariable(j)
320 << ") : " << vL2Error << std::endl;
321 out << "L inf error (variable " << m_equ[0]->GetVariable(j)
322 << ") : " << vLinfError << std::endl;
323 }
324 }
325
326 // Process eigenvectors and write out.
327 m_nvec = converged;
328 EV_post(Tseq, Kseqcopy, ntot, std::min(--i, m_kdim), m_nvec, zvec, wr, wi,
329 converged);
330
331 WARNINGL0(m_imagShift == 0, "Complex Shift applied. "
332 "Need to implement Ritz re-evaluation of"
333 "eigenvalue. Only one half of complex "
334 "value will be correct");
335
336 // store eigenvalues so they can be accessed from driver class
337 m_real_evl = wr;
338 m_imag_evl = wi;
339
340 // Close the runtime info file.
341 if (m_comm->GetRank() == 0)
342 {
343 evlout.close();
344 }
345}
346
347/**
348 *
349 */
352{
353 // Copy starting vector into first sequence element.
355 m_equ[0]->TransCoeffToPhys();
356
357 m_equ[0]->SetTime(0.);
358 m_equ[0]->DoSolve();
359
361 {
363 fields = m_equ[0]->UpdateFields();
364
365 // start Adjoint with latest fields of direct
366 CopyFwdToAdj();
367 m_equ[1]->TransCoeffToPhys();
368
369 m_equ[1]->SetTime(0.);
370 m_equ[1]->DoSolve();
371 }
372
373 // Copy starting vector into first sequence element.
375}
376
377/**
378 * Computes the Ritz eigenvalues and eigenvectors of the Hessenberg matrix
379 * constructed from the Krylov sequence.
380 */
383 Array<OneD, Array<OneD, NekDouble>> &Kseqcopy, const int ntot,
384 const Array<OneD, NekDouble> &alpha, const int kdim,
386 Array<OneD, NekDouble> &wi, NekDouble &resnorm)
387{
388 int kdimp = kdim + 1;
389 int lwork = 10 * kdim;
390 int ier;
391 Array<OneD, NekDouble> R(kdimp * kdimp, 0.0);
392 Array<OneD, NekDouble> H(kdimp * kdim, 0.0);
393 Array<OneD, NekDouble> rwork(lwork, 0.0);
394
395 // Modified G-S orthonormalisation
396 for (int i = 0; i < kdimp; ++i)
397 {
398 NekDouble gsc = Blas::Ddot(ntot, &Kseq[i][0], 1, &Kseq[i][0], 1);
400 gsc = std::sqrt(gsc);
401 ASSERTL0(gsc != 0.0, "Vectors are linearly independent.");
402
403 R[i * kdimp + i] = gsc;
404 Vmath::Smul(ntot, 1.0 / gsc, Kseq[i], 1, Kseq[i], 1);
405 if (m_useMask)
406 {
407 Vmath::Smul(ntot, 1.0 / gsc, Kseqcopy[i], 1, Kseqcopy[i], 1);
408 }
409
410 for (int j = i + 1; j < kdimp; ++j)
411 {
412 gsc = Blas::Ddot(ntot, &Kseq[i][0], 1, &Kseq[j][0], 1);
414 Vmath::Svtvp(ntot, -gsc, Kseq[i], 1, Kseq[j], 1, Kseq[j], 1);
415 if (m_useMask)
416 {
417 Vmath::Svtvp(ntot, -gsc, Kseqcopy[i], 1, Kseqcopy[j], 1,
418 Kseqcopy[j], 1);
419 }
420 R[j * kdimp + i] = gsc;
421 }
422 }
423
424 // Compute H matrix
425 for (int i = 0; i < kdim; ++i)
426 {
427 for (int j = 0; j < kdim; ++j)
428 {
429 H[j * kdim + i] =
430 alpha[j + 1] * R[(j + 1) * kdimp + i] -
431 Vmath::Dot(j, &H[0] + i, kdim, &R[0] + j * kdimp, 1);
432 H[j * kdim + i] /= R[j * kdimp + j];
433 }
434 }
435
436 H[(kdim - 1) * kdim + kdim] =
437 alpha[kdim] *
438 std::fabs(R[kdim * kdimp + kdim] / R[(kdim - 1) * kdimp + kdim - 1]);
439
440 Lapack::dgeev_('N', 'V', kdim, &H[0], kdim, &wr[0], &wi[0], nullptr, 1,
441 &zvec[0], kdim, &rwork[0], lwork, ier);
442
443 ASSERTL0(!ier, "Error with dgeev");
444
445 resnorm = H[(kdim - 1) * kdim + kdim];
446}
447
448/**
449 * Check for convergence of the residuals of the eigenvalues of H.
450 */
451int DriverModifiedArnoldi::EV_test(const int itrn, const int kdim,
455 const NekDouble resnorm, int nvec,
456 std::ofstream &evlout, NekDouble &resid0)
457{
458 int idone = 0;
459 // NekDouble period = 0.1;
460
461 Array<OneD, NekDouble> resid(kdim);
462 for (int i = 0; i < kdim; ++i)
463 {
464 NekDouble tmp = std::sqrt(
465 Vmath::Dot(kdim, &zvec[0] + i * kdim, 1, &zvec[0] + i * kdim, 1));
466 resid[i] = resnorm * std::fabs(zvec[kdim - 1 + i * kdim]) / tmp;
467 if (wi[i] < 0.0)
468 {
469 resid[i - 1] = resid[i] = hypot(resid[i - 1], resid[i]);
470 }
471 }
472
473 EV_sort(zvec, wr, wi, resid, kdim);
474
475 while (nvec <= kdim && resid[nvec - 1] < m_evtol)
476 {
477 idone = nvec;
478 ++nvec;
479 }
480 nvec -= (idone > 0);
481
482 if (m_comm->GetRank() == 0)
483 {
484 evlout << "-- Iteration = " << itrn << ", H(k+1, k) = " << resnorm
485 << std::endl;
486 evlout.precision(4);
487 evlout.setf(std::ios::scientific, std::ios::floatfield);
489 {
490 evlout << " Magnitude Angle Growth "
491 << "Frequency Residual" << std::endl;
492 }
493 else
494 {
495 evlout << " Real Imaginary inverse real "
496 << "inverse imag Residual" << std::endl;
497 }
498
499 for (int i = 0; i < kdim; i++)
500 {
501 WriteEvs(evlout, i, wr[i], wi[i], resid[i]);
502 }
503 }
504
505 resid0 = resid[nvec - 1];
506 return idone;
507}
508
509/**
510 * Sorts the computed eigenvalues by smallest residual first.
511 */
515 Array<OneD, NekDouble> &test, const int dim)
516{
517 Array<OneD, NekDouble> z_tmp(dim, 0.0);
518 NekDouble wr_tmp, wi_tmp, te_tmp;
519 for (int j = 1; j < dim; ++j)
520 {
521 wr_tmp = wr[j];
522 wi_tmp = wi[j];
523 te_tmp = test[j];
524 Vmath::Vcopy(dim, &evec[0] + j * dim, 1, &z_tmp[0], 1);
525 int i = j - 1;
526 while (i >= 0 && test[i] > te_tmp)
527 {
528 wr[i + 1] = wr[i];
529 wi[i + 1] = wi[i];
530 test[i + 1] = test[i];
531 Vmath::Vcopy(dim, &evec[0] + i * dim, 1, &evec[0] + (i + 1) * dim,
532 1);
533 i--;
534 }
535 wr[i + 1] = wr_tmp;
536 wi[i + 1] = wi_tmp;
537 test[i + 1] = te_tmp;
538 Vmath::Vcopy(dim, &z_tmp[0], 1, &evec[0] + (i + 1) * dim, 1);
539 }
540}
541
542/**
543 * Post-process the Ritz eigenvalues/eigenvectors of the matrix H, to compute
544 * estimations of the leading eigenvalues and eigenvectors of the original
545 * matrix.
546 */
549 const int ntot, const int kdim,
550 const int nvec,
553 Array<OneD, NekDouble> &wi, const int icon)
554{
555 if (icon == 0)
556 {
557 // Not converged, write final Krylov vector
558 ASSERTL0(false, "Convergence was not achieved within the "
559 "prescribed number of iterations.");
560 }
561 else if (icon < 0)
562 {
563 // Minimum residual reached
564 ASSERTL0(false, "Minimum residual reached.");
565 }
566 else if (icon == nvec)
567 {
568 // Converged, write out eigenvectors
569 EV_big(Tseq, Kseq, ntot, kdim, icon, zvec, wr, wi);
571 m_equ[0]->UpdateFields();
572
573 for (int j = 0; j < icon; ++j)
574 {
575 std::string file = m_session->GetSessionName() + "_eig_" +
576 std::to_string(j) + ".fld";
577
578 if (m_comm->GetRank() == 0)
579 {
580 WriteEvs(std::cout, j, wr[j], wi[j]);
581 }
582 WriteFld(file, Kseq[j]);
583 if (m_useMask)
584 {
585 std::string fileunmask = m_session->GetSessionName() +
586 "_eig_masked_" + std::to_string(j) +
587 ".fld";
588 WriteFld(fileunmask, Tseq[j]);
589 }
590 }
591 }
592 else
593 {
594 // Not recognised
595 ASSERTL0(false, "Unrecognised value.");
596 }
597}
598
599/**
600 * Compute estimates of the eigenvalues/eigenvectors of the original system.
601 */
604 const int ntot, const int kdim,
605 const int nvec, Array<OneD, NekDouble> &zvec,
606 [[maybe_unused]] Array<OneD, NekDouble> &wr,
608{
609 NekDouble wgt, norm;
612 for (int i = 0; i < nvec; ++i)
613 {
614 if (m_useMask)
615 {
616 btmp[i] = Array<OneD, NekDouble>(ntot, 0.);
617 etmp[i] = Array<OneD, NekDouble>(ntot, 0.);
618 }
619 else
620 {
621 btmp[i] = evecs[i];
622 }
623 }
624
625 // Generate big eigenvectors
626 for (int j = 0; j < nvec; ++j)
627 {
628 Vmath::Zero(ntot, btmp[j], 1);
629 if (m_useMask)
630 {
631 Vmath::Zero(ntot, etmp[j], 1);
632 }
633 for (int i = 0; i < kdim; ++i)
634 {
635 wgt = zvec[i + j * kdim];
636 Vmath::Svtvp(ntot, wgt, bvecs[i], 1, btmp[j], 1, btmp[j], 1);
637 if (m_useMask)
638 {
639 Vmath::Svtvp(ntot, wgt, evecs[i], 1, etmp[j], 1, etmp[j], 1);
640 }
641 }
642 }
643
644 // Normalise the big eigenvectors
645 for (int i = 0; i < nvec; ++i)
646 {
647 if (wi[i] == 0.0) // Real mode
648 {
649 if (m_session->GetComm()->GetRank() == 0)
650 {
651 std::cout << "eigenvalue " << i << ": real mode" << std::endl;
652 }
653 norm = Blas::Ddot(ntot, &btmp[i][0], 1, &btmp[i][0], 1);
654 m_comm->AllReduce(norm, Nektar::LibUtilities::ReduceSum);
655 norm = std::sqrt(norm);
656 if (m_useMask)
657 {
658 Vmath::Smul(ntot, 1.0 / norm, btmp[i], 1, bvecs[i], 1);
659 Vmath::Smul(ntot, 1.0 / norm, etmp[i], 1, evecs[i], 1);
660 }
661 else
662 {
663 Vmath::Smul(ntot, 1.0 / norm, btmp[i], 1, evecs[i], 1);
664 }
665 }
666 else
667 {
668 if (m_session->GetComm()->GetRank() == 0)
669 {
670 std::cout << "eigenvalues " << i << ", " << i + 1
671 << ": complex modes" << std::endl;
672 }
673 norm = Blas::Ddot(ntot, &btmp[i][0], 1, &btmp[i][0], 1);
674 norm += Blas::Ddot(ntot, &btmp[i + 1][0], 1, &btmp[i + 1][0], 1);
675 m_comm->AllReduce(norm, Nektar::LibUtilities::ReduceSum);
676 norm = std::sqrt(norm);
677
678 if (m_useMask)
679 {
680 Vmath::Smul(ntot, 1.0 / norm, btmp[i], 1, bvecs[i], 1);
681 Vmath::Smul(ntot, 1.0 / norm, btmp[i + 1], 1, bvecs[i + 1], 1);
682 Vmath::Smul(ntot, 1.0 / norm, etmp[i], 1, evecs[i], 1);
683 Vmath::Smul(ntot, 1.0 / norm, etmp[i + 1], 1, evecs[i + 1], 1);
684 }
685 else
686 {
687 Vmath::Smul(ntot, 1.0 / norm, btmp[i], 1, evecs[i], 1);
688 Vmath::Smul(ntot, 1.0 / norm, btmp[i + 1], 1, evecs[i + 1], 1);
689 }
690
691 i++;
692 }
693 }
694}
695
696} // namespace Nektar::SolverUtils
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
#define WARNINGL0(condition, msg)
Definition: ErrorUtil.hpp:215
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
static std::string RegisterEnumValue(std::string pEnum, std::string pString, int pEnumValue)
Registers an enumeration value.
Base class for the development of solvers.
Definition: DriverArnoldi.h:45
void CopyFwdToAdj()
Copy the forward field to the adjoint system in transient growth calculations.
void v_InitObject(std::ostream &out=std::cout) override
Virtual function for initialisation implementation.
void WriteFld(std::string file, std::vector< Array< OneD, NekDouble > > coeffs)
Write coefficients to file.
Array< OneD, NekDouble > m_maskCoeffs
Definition: DriverArnoldi.h:78
void CopyFieldToArnoldiArray(Array< OneD, NekDouble > &array)
Copy fields to Arnoldi storage.
int m_nvec
Dimension of Krylov subspace.
Definition: DriverArnoldi.h:61
bool m_timeSteppingAlgorithm
Period of time stepping algorithm.
Definition: DriverArnoldi.h:65
int m_nits
Number of vectors to test.
Definition: DriverArnoldi.h:62
Array< OneD, NekDouble > m_imag_evl
Definition: DriverArnoldi.h:75
void CopyArnoldiArrayToField(Array< OneD, NekDouble > &array)
Copy Arnoldi storage to fields.
NekDouble m_evtol
Maxmum number of iterations.
Definition: DriverArnoldi.h:63
int m_nfields
interval to dump information if required.
Definition: DriverArnoldi.h:69
SOLVER_UTILS_EXPORT void ArnoldiSummary(std::ostream &out)
Array< OneD, NekDouble > m_real_evl
Operator in solve call is negated.
Definition: DriverArnoldi.h:74
void WriteEvs(std::ostream &evlout, const int k, const NekDouble real, const NekDouble imag, NekDouble resid=NekConstants::kNekUnsetDouble, bool DumpInverse=true)
LibUtilities::SessionReaderSharedPtr m_session
Session reader object.
Definition: Driver.h:83
LibUtilities::CommSharedPtr m_comm
Communication object.
Definition: Driver.h:80
enum EvolutionOperatorType m_EvolutionOperator
Evolution Operator.
Definition: Driver.h:98
Array< OneD, EquationSystemSharedPtr > m_equ
Equation system to solve.
Definition: Driver.h:92
int m_nequ
number of equations
Definition: Driver.h:95
void EV_update(Array< OneD, NekDouble > &src, Array< OneD, NekDouble > &tgt)
Generates a new vector in the sequence by applying the linear operator.
void EV_post(Array< OneD, Array< OneD, NekDouble > > &Tseq, Array< OneD, Array< OneD, NekDouble > > &Kseq, const int ntot, const int kdim, const int nvec, Array< OneD, NekDouble > &zvec, Array< OneD, NekDouble > &wr, Array< OneD, NekDouble > &wi, const int icon)
int EV_test(const int itrn, const int kdim, Array< OneD, NekDouble > &zvec, Array< OneD, NekDouble > &wr, Array< OneD, NekDouble > &wi, const NekDouble resnorm, int nvec, std::ofstream &evlout, NekDouble &resid0)
Tests for convergence of eigenvalues of H.
void EV_small(Array< OneD, Array< OneD, NekDouble > > &Kseq, Array< OneD, Array< OneD, NekDouble > > &Kseqcopy, const int ntot, const Array< OneD, NekDouble > &alpha, const int kdim, Array< OneD, NekDouble > &zvec, Array< OneD, NekDouble > &wr, Array< OneD, NekDouble > &wi, NekDouble &resnorm)
Generates the upper Hessenberg matrix H and computes its eigenvalues.
void EV_sort(Array< OneD, NekDouble > &evec, Array< OneD, NekDouble > &wr, Array< OneD, NekDouble > &wi, Array< OneD, NekDouble > &test, const int dim)
Sorts a sequence of eigenvectors/eigenvalues by magnitude.
static DriverSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
DriverModifiedArnoldi(const LibUtilities::SessionReaderSharedPtr pSession, const SpatialDomains::MeshGraphSharedPtr pGraph)
Constructor.
static std::string className
Name of the class.
void v_InitObject(std::ostream &out=std::cout) override
Virtual function for initialisation implementation.
void v_Execute(std::ostream &out=std::cout) override
Virtual function for solve implementation.
void EV_big(Array< OneD, Array< OneD, NekDouble > > &bvecs, Array< OneD, Array< OneD, NekDouble > > &evecs, const int ntot, const int kdim, const int nvec, Array< OneD, NekDouble > &zvec, Array< OneD, NekDouble > &wr, Array< OneD, NekDouble > &wi)
static double Ddot(const int &n, const double *x, const int &incx, const double *y, const int &incy)
BLAS level 1: output = .
Definition: Blas.hpp:163
std::shared_ptr< SessionReader > SessionReaderSharedPtr
DriverFactory & GetDriverFactory()
Definition: Driver.cpp:64
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:174
MultiRegions::ContFieldSharedPtr contfield
double NekDouble
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.hpp:72
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.hpp:396
T Dot(int n, const T *w, const T *x)
dot product
Definition: Vmath.hpp:761
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.hpp:100
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.hpp:273
void FillWhiteNoise(int n, const T eps, T *x, const int incx, int outseed)
Fills a vector with white noise.
Definition: Vmath.cpp:154
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:285