Virtual function for solve implementation.
121{
122 Array<OneD, NekDouble> tmpworkd;
123
125 ->UpdateFields()[0]
126 ->GetNcoeffs();
129 int ido;
130
131 int info;
132
133
134 int iparam[11];
135 int ipntr[14];
136
137 Array<OneD, int> ritzSelect;
138 Array<OneD, NekDouble> dr;
139 Array<OneD, NekDouble> di;
140 Array<OneD, NekDouble> workev;
141 Array<OneD, NekDouble>
z;
143
144 Array<OneD, NekDouble> resid(n);
145 Array<OneD, NekDouble> v(n *
m_kdim);
146 Array<OneD, NekDouble> workl(lworkl, 0.0);
147 Array<OneD, NekDouble> workd(3 * n, 0.0);
148
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;
166 iparam[1] = 0;
168 iparam[3] = 1;
169 iparam[4] = 0;
170 iparam[5] = 0;
171
172
174 {
175 iparam[6] = 1;
176 B = 'I';
177 }
178 else
179 {
180 iparam[6] = 3;
181 B = 'G';
182 }
183#if 0
186 {
187 iparam[6] = 3;
188 B = 'G';
189 }
190 else
191 {
192 iparam[6] = 1;
193 B = 'I';
194 }
195#endif
196 iparam[7] = 0;
197 iparam[8] = 0;
198 iparam[9] = 0;
199 iparam[10] = 0;
200
201 int cycle = 0;
202 const char *problem =
204 "ArpackProblemType")]
205 .c_str();
206
208 std::ofstream pFile(
name.c_str());
209
210 ido = 0;
211
212 while (ido != 99)
213 {
214
217 iparam, ipntr, &workd[0], &workl[0], lworkl, info);
218
219
220
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
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:
259
260
261
263
264 m_equ[0]->TransCoeffToPhys();
265
268 {
269
271
272 m_equ[1]->TransCoeffToPhys();
274 }
275
277 {
278 out << std::endl;
280 }
281
282
284
285 cycle++;
286 break;
287 case 2:
288 {
289
291
292 m_equ[0]->TransCoeffToPhys();
293
294 Array<OneD, MultiRegions::ExpListSharedPtr> fields =
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
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
323 {
325 }
326 else
327 {
329 }
330
331
332
333
334 sigmai = 0;
335
336
337 Arpack::Dneupd(1,
"A", ritzSelect.data(), dr.data(), di.data(),
z.data(), n,
338 sigmar, sigmai, workev.data(), &B, n, problem,
m_nvec,
340 workd.data(), workl.data(), lworkl, info);
341
342 ASSERTL0(info == 0,
" Error with Dneupd");
343
344 int nconv = iparam[4];
345
346
348 {
350 }
351 else
352 {
354 }
355
357 "Need to implement Ritz re-evaluation of"
358 "eigenvalue. Only one half of complex "
359 "value will be correct");
360
361 Array<OneD, MultiRegions::ExpListSharedPtr> fields =
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 {
377
378 std::string file =
m_session->GetSessionName() +
"_eig_" +
379 std::to_string(i) + ".fld";
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 {
390 false);
392 false);
393
394 std::string file =
m_session->GetSessionName() +
"_eig_" +
395 std::to_string(i) + ".fld";
397 }
398 }
399
402
403 pFile.close();
404
405 for (
int j = 0; j <
m_equ[0]->GetNvariables(); ++j)
406 {
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}
#define ASSERTL0(condition, msg)
#define WARNINGL0(condition, msg)
void CopyFwdToAdj()
Copy the forward field to the adjoint system in transient growth calculations.
void WriteFld(std::string file, std::vector< Array< OneD, NekDouble > > coeffs)
Write coefficients to file.
int m_infosteps
underlying operator is time stepping
void CopyFieldToArnoldiArray(Array< OneD, NekDouble > &array)
Copy fields to Arnoldi storage.
int m_nvec
Dimension of Krylov subspace.
bool m_timeSteppingAlgorithm
Period of time stepping algorithm.
int m_nits
Number of vectors to test.
Array< OneD, NekDouble > m_imag_evl
void CopyArnoldiArrayToField(Array< OneD, NekDouble > &array)
Copy Arnoldi storage to fields.
NekDouble m_evtol
Maxmum number of iterations.
int m_nfields
interval to dump information if required.
Array< OneD, NekDouble > m_real_evl
Operator in solve call is negated.
void WriteEvs(std::ostream &evlout, const int k, const NekDouble real, const NekDouble imag, NekDouble resid=NekConstants::kNekUnsetDouble, bool DumpInverse=true)
static std::string ArpackProblemTypeTrans[]
LibUtilities::SessionReaderSharedPtr m_session
Session reader object.
LibUtilities::CommSharedPtr m_comm
Communication object.
enum EvolutionOperatorType m_EvolutionOperator
Evolution Operator.
Array< OneD, EquationSystemSharedPtr > m_equ
Equation system to solve.
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.
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.
static const NekDouble kNekUnsetDouble
static const NekDouble kNekZeroTol
std::vector< double > z(NPUPPER)
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.