Nektar++
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// 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: Arnoldi solver using Arpack
32//
33///////////////////////////////////////////////////////////////////////////////
34
36
37namespace Nektar::SolverUtils
38{
39
42 "LargestReal", 0),
44 "SmallestReal", 1),
46 "LargestImag", 2),
48 "SmallestImag", 3),
50 "LargestMag", 4),
52 "SmallestMag", 5),
53};
56 "LargestMag");
59
60std::string DriverArpack::className =
62
63std::string DriverArpack::ArpackProblemTypeTrans[6] = {"LR", "SR", "LI",
64 "SI", "LM", "SM"};
65
66/**
67 *
68 */
71 : DriverArnoldi(pSession, pGraph)
72{
73}
74
75/**
76 *
77 */
78void DriverArpack::v_InitObject(std::ostream &out)
79{
81
82 // Initialisation of Arnoldi parameters
83 m_maxn = 1000000; // Maximum size of the problem
84 m_maxnev = 200; // maximum number of eigenvalues requested
85 m_maxncv = 500; // Largest number of basis vector used in Implicitly
86 // Restarted Arnoldi
87
88 // Error alerts
89 ASSERTL0(m_nvec <= m_maxnev, "NEV is greater than MAXNEV");
90 ASSERTL0(m_kdim <= m_maxncv, "NEV is greater than MAXNEV");
91 ASSERTL0(2 <= m_kdim - m_nvec, "NCV-NEV is less than 2");
92
93 ASSERTL0(m_comm->GetSize() == 1,
94 "..ARPACK is not currently set-up for parallel execution...\n");
95
96 m_equ[0]->PrintSummary(out);
97
98 // Print session parameters
99 out << "\tArnoldi solver type : Arpack" << std::endl;
100 out << "\tArpack problem type : ";
101 out << ArpackProblemTypeTrans[m_session->GetSolverInfoAsEnum<int>(
102 "ArpackProblemType")]
103 << std::endl;
105
106 // Initialization
107 for (int i = 0; i < m_nequ; ++i)
108 {
109 m_equ[i]->DoInitialise();
110 }
111
112 // FwdTrans Initial conditions to be in Coefficient Space
113 m_equ[m_nequ - 1]->TransPhysToCoeff();
114}
115
116/**
117 *
118 */
119void DriverArpack::v_Execute(std::ostream &out)
120
121{
122 Array<OneD, NekDouble> tmpworkd;
123
124 int nq = m_equ[0]
125 ->UpdateFields()[0]
126 ->GetNcoeffs(); // Number of points in the mesh
127 int n = m_nfields * nq; // Number of points in eigenvalue calculation
128 int lworkl = 3 * m_kdim * (m_kdim + 2); // Size of work array
129 int ido; // REVERSE COMMUNICATION parameter. At the first call must be
130 // initialised at 0
131 int info; // do not set initial vector (info=0 random initial vector, info=1
132 // read initial vector from session file)
133
134 int iparam[11];
135 int ipntr[14];
136
137 Array<OneD, int> ritzSelect;
142 NekDouble sigmar, sigmai;
143
144 Array<OneD, NekDouble> resid(n);
146 Array<OneD, NekDouble> workl(lworkl, 0.0);
147 Array<OneD, NekDouble> workd(3 * n, 0.0);
148
149 ASSERTL0(n <= m_maxn, "N is greater than MAXN");
150
151 if (m_session->DefinesFunction("InitialConditions"))
152 {
153 out << "\tInital vector : input file " << std::endl;
154 info = 1;
156 }
157 else
158 {
159 out << "\tInital vector : random " << std::endl;
160 info = 0;
161 }
162
163 char B;
164
165 iparam[0] = 1; // strategy for shift-invert
166 iparam[1] = 0; // (deprecated)
167 iparam[2] = m_nits; // maximum number of iterations allowed/taken
168 iparam[3] = 1; // blocksize to be used for recurrence
169 iparam[4] = 0; // number of converged ritz eigenvalues
170 iparam[5] = 0; // (deprecated)
171
172 // Use generalized B matrix for coupled solver.
174 {
175 iparam[6] = 1; // computation mode 1=> matrix-vector prod
176 B = 'I';
177 }
178 else
179 {
180 iparam[6] = 3; // computation mode 1=> matrix-vector prod
181 B = 'G';
182 }
183#if 0
184 if((fabs(m_realShift) > NekConstants::kNekZeroTol)|| // use shift if m_realShift > 1e-12
186 {
187 iparam[6] = 3; // This was 3 need to know what to set it to
188 B = 'G';
189 }
190 else
191 {
192 iparam[6] = 1; // computation mode 1=> matrix-vector prod
193 B = 'I';
194 }
195#endif
196 iparam[7] = 0; // (for shift-invert)
197 iparam[8] = 0; // number of MV operations
198 iparam[9] = 0; // number of BV operations
199 iparam[10] = 0; // number of reorthogonalisation steps
200
201 int cycle = 0;
202 const char *problem =
203 ArpackProblemTypeTrans[m_session->GetSolverInfoAsEnum<int>(
204 "ArpackProblemType")]
205 .c_str();
206
207 std::string name = m_session->GetSessionName() + ".evl";
208 std::ofstream pFile(name.c_str());
209
210 ido = 0; // At the first call must be initialisedat 0
211
212 while (ido != 99) // ido==-1 || ido==1 || ido==0)
213 {
214 // Routine for eigenvalue evaluation for non-symmetric operators
215 Arpack::Dnaupd(ido, &B, // B='I' for std eval problem
216 n, problem, m_nvec, m_evtol, &resid[0], m_kdim, &v[0], n,
217 iparam, ipntr, &workd[0], &workl[0], lworkl, info);
218
219 // Plotting of real and imaginary part of the
220 // eigenvalues from workl
221 out << "\rIteration " << cycle << ", output: " << info
222 << ", ido=" << ido << " " << std::flush;
223
224 if (!((cycle - 1) % m_kdim) && (cycle > m_kdim) && (ido != 2))
225 {
226 pFile << "Krylov spectrum at iteration: " << cycle << std::endl;
227
229 {
230 pFile << "EV Magnitude Angle Growth Frequency "
231 "Residual"
232 << std::endl;
233 }
234 else
235 {
236 pFile << "EV Real Imaginary inverse real inverse "
237 "imag Residual"
238 << std::endl;
239 }
240
241 out << std::endl;
242 for (int k = 0; k < m_kdim; ++k)
243 {
244 // write m_kdim eigs to screen
245 WriteEvs(pFile, k, workl[ipntr[5] - 1 + k],
246 workl[ipntr[6] - 1 + k]);
247 }
248 }
249
250 if (ido == 99)
251 {
252 break;
253 }
254
255 switch (ido)
256 {
257 case -1:
258 case 1: // Note that ido=1 we are using input x
259 // (workd[inptr[0]-1]) rather than Mx as
260 // recommended in manual since it is not
261 // possible to impose forcing directly.
262 CopyArnoldiArrayToField(tmpworkd = workd + (ipntr[0] - 1));
263
264 m_equ[0]->TransCoeffToPhys();
265
266 m_equ[0]->DoSolve();
268 {
269 // start Adjoint with latest fields of direct
270 CopyFwdToAdj();
271
272 m_equ[1]->TransCoeffToPhys();
273 m_equ[1]->DoSolve();
274 }
275
276 if (!(cycle % m_infosteps))
277 {
278 out << std::endl;
279 m_equ[0]->Output();
280 }
281
282 // operated fields are copied into workd[inptr[1]-1]
283 CopyFieldToArnoldiArray(tmpworkd = workd + (ipntr[1] - 1));
284
285 cycle++;
286 break;
287 case 2: // provide y = M x (bwd trans and iproduct);
288 {
289 // workd[inptr[0]-1] copied into operator fields
290 CopyArnoldiArrayToField(tmpworkd = workd + (ipntr[0] - 1));
291
292 m_equ[0]->TransCoeffToPhys();
293
295 m_equ[0]->UpdateFields();
296 for (int i = 0; i < fields.size(); ++i)
297 {
298 fields[i]->IProductWRTBase(fields[i]->GetPhys(),
299 fields[i]->UpdateCoeffs());
300 }
301
302 // operated fields are copied into workd[inptr[1]-1]
303 CopyFieldToArnoldiArray(tmpworkd = workd + (ipntr[1] - 1));
304 break;
305 }
306 default:
307 ASSERTL0(false, "Unexpected reverse communication request.");
308 }
309 }
310
311 out << std::endl
312 << "Converged in " << iparam[8] << " iterations" << std::endl;
313
314 ASSERTL0(info >= 0, " Error with Dnaupd");
315
316 ritzSelect = Array<OneD, int>(m_kdim, 0);
317 dr = Array<OneD, NekDouble>(m_nvec + 1, 0.0);
318 di = Array<OneD, NekDouble>(m_nvec + 1, 0.0);
319 workev = Array<OneD, NekDouble>(3 * m_kdim);
320 z = Array<OneD, NekDouble>(n * (m_nvec + 1));
321
322 if (m_negatedOp)
323 {
324 sigmar = -m_realShift;
325 }
326 else
327 {
328 sigmar = m_realShift;
329 }
330
331 // Do not pass imaginary shift to Arpack since we have not
332 // used a Fortran complex number format and so processing
333 // is mucked up. Need to do some processing afterwards.
334 sigmai = 0;
335
336 // Setting 'A', Ritz vectors are computed. 'S' for Shur vectors
337 Arpack::Dneupd(1, "A", ritzSelect.data(), dr.data(), di.data(), z.data(), n,
338 sigmar, sigmai, workev.data(), &B, n, problem, m_nvec,
339 m_evtol, resid.data(), m_kdim, v.data(), n, iparam, ipntr,
340 workd.data(), workl.data(), lworkl, info);
341
342 ASSERTL0(info == 0, " Error with Dneupd");
343
344 int nconv = iparam[4];
345
346 // Subtract off complex shift if it exists
347 if (m_negatedOp)
348 {
349 Vmath::Sadd(nconv, m_imagShift, di, 1, di, 1);
350 }
351 else
352 {
353 Vmath::Sadd(nconv, -m_imagShift, di, 1, di, 1);
354 }
355
356 WARNINGL0(m_imagShift == 0, "Complex Shift applied. "
357 "Need to implement Ritz re-evaluation of"
358 "eigenvalue. Only one half of complex "
359 "value will be correct");
360
362 m_equ[0]->UpdateFields();
363
364 out << "Converged Eigenvalues: " << nconv << std::endl;
365 pFile << "Converged Eigenvalues: " << nconv << std::endl;
366
368 {
369 out << " Magnitude Angle Growth Frequency"
370 << std::endl;
371 pFile << " Magnitude Angle Growth Frequency"
372 << std::endl;
373 for (int i = 0; i < nconv; ++i)
374 {
375 WriteEvs(out, i, dr[i], di[i]);
376 WriteEvs(pFile, i, dr[i], di[i]);
377
378 std::string file = m_session->GetSessionName() + "_eig_" +
379 std::to_string(i) + ".fld";
380 WriteFld(file, z + i * n);
381 }
382 }
383 else
384 {
385 out << " Real Imaginary " << std::endl;
386 pFile << " Real Imaginary " << std::endl;
387 for (int i = 0; i < nconv; ++i)
388 {
389 WriteEvs(out, i, dr[i], di[i], NekConstants::kNekUnsetDouble,
390 false);
391 WriteEvs(pFile, i, dr[i], di[i], NekConstants::kNekUnsetDouble,
392 false);
393
394 std::string file = m_session->GetSessionName() + "_eig_" +
395 std::to_string(i) + ".fld";
396 WriteFld(file, z + i * n);
397 }
398 }
399
400 m_real_evl = dr;
401 m_imag_evl = di;
402
403 pFile.close();
404
405 for (int j = 0; j < m_equ[0]->GetNvariables(); ++j)
406 {
407 NekDouble vL2Error = m_equ[0]->L2Error(j, false);
408 NekDouble vLinfError = m_equ[0]->LinfError(j);
409 if (m_comm->GetRank() == 0)
410 {
411 out << "L 2 error (variable " << m_equ[0]->GetVariable(j)
412 << ") : " << vL2Error << std::endl;
413 out << "L inf error (variable " << m_equ[0]->GetVariable(j)
414 << ") : " << vLinfError << std::endl;
415 }
416 }
417}
418
419} // 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.
static std::string RegisterDefaultSolverInfo(const std::string &pName, const std::string &pValue)
Registers the default string value of a solver info property.
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.
int m_infosteps
underlying operator is time stepping
Definition: DriverArnoldi.h:67
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)
static std::string arpackProblemTypeDefault
Definition: DriverArpack.h:88
static DriverSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
Definition: DriverArpack.h:52
DriverArpack(const LibUtilities::SessionReaderSharedPtr pSession, const SpatialDomains::MeshGraphSharedPtr pGraph)
Constructor.
void v_InitObject(std::ostream &out=std::cout) override
Virtual function for initialisation implementation.
static std::string driverLookupId
Definition: DriverArpack.h:84
static std::string arpackProblemTypeLookupIds[]
Definition: DriverArpack.h:87
void v_Execute(std::ostream &out=std::cout) override
Virtual function for solve implementation.
static std::string className
Name of the class.
Definition: DriverArpack.h:63
static std::string ArpackProblemTypeTrans[]
Definition: DriverArpack.h:89
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
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:114
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:128
std::shared_ptr< SessionReader > SessionReaderSharedPtr
static const NekDouble kNekUnsetDouble
static const NekDouble kNekZeroTol
DriverFactory & GetDriverFactory()
Definition: Driver.cpp:64
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:174
std::vector< double > z(NPUPPER)
double NekDouble
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.hpp:194