Nektar++
DriverSteadyState.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: DriverSteadyState.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 class for the steady-state solver
32//
33///////////////////////////////////////////////////////////////////////////////
34
37
38#include <boost/algorithm/string/classification.hpp>
39#include <boost/algorithm/string/predicate.hpp>
40#include <boost/algorithm/string/split.hpp>
41
42namespace Nektar::SolverUtils
43{
44
49 LibUtilities::SessionReader::RegisterEnumValue("Driver", "SteadyState", 0);
50
51/**
52 *
53 */
57 : DriverModifiedArnoldi(pSession, pGraph)
58{
59}
60
61/**
62 *
63 */
64void DriverSteadyState::v_InitObject(std::ostream &out)
65{
67}
68
69/**
70 *
71 */
72void DriverSteadyState::v_Execute(std::ostream &out)
73{
74 // With a loop over "DoSolve", this Driver implements the
75 // "encapsulated" Selective Frequency Damping method(SFD)
76 // to find the steady state of a flow above the critical Reynolds
77 // number.
78 m_equ[m_nequ - 1]->PrintSummary(out);
79
80 m_session->LoadParameter("IO_InfoSteps", m_infosteps, 1000);
81 m_session->LoadParameter("IO_CheckSteps", m_checksteps, 100000);
82 m_session->LoadParameter("ControlCoeff", m_X, 1);
83 m_session->LoadParameter("FilterWidth", m_Delta, 2);
84
85 // To evaluate optimum SFD parameters if growth rate provided in the
86 // xml file
87 m_session->LoadParameter("GrowthRateEV", GrowthRateEV, 0.0);
88
89 // To evaluate optimum SFD parameters if frequency provided in the xml
90 // file
91 m_session->LoadParameter("FrequencyEV", FrequencyEV, 0.0);
92
93 // Determine when SFD method is converged
94 m_session->LoadParameter("TOL", TOL, 1.0e-08);
95
96 // Used only for the Adaptive SFD method
97 m_session->LoadParameter("AdaptiveTOL", AdaptiveTOL, 1.0e-02);
98
99 // Used only for the Adaptive SFD method
100 m_session->LoadParameter("AdaptiveTime", AdaptiveTime, 50.0 * m_Delta);
101
102 if (m_comm->GetRank() == 0)
103 {
105 }
106
107 timer.Start();
108
109 // Definition of shared pointer (used only for the Adaptive SFD method)
111
112 NumVar_SFD = m_equ[m_nequ - 1]->UpdateFields()[0]->GetCoordim(0);
113 // SFD to run for incompressible case with scalar field
114 if (m_session->GetSolverInfo("EqType") == "UnsteadyNavierStokes" ||
115 m_session->GetSolverInfo("EqType") == "SteadyNavierStokes")
116 {
117 int nConvectiveFields = m_session->GetVariables().size();
118 if (boost::iequals(m_session->GetVariable(nConvectiveFields - 1), "p"))
119 {
120 nConvectiveFields -= 1;
121 }
122 NumVar_SFD = nConvectiveFields;
123 }
124 // Condition necessary to run SFD for the compressible case
125 if (m_session->GetSolverInfo("EqType") == "EulerCFE" ||
126 m_session->GetSolverInfo("EqType") == "NavierStokesCFE" ||
127 m_session->GetSolverInfo("EqType") == "NavierStokesImplicitCFE")
128 {
129 // Number of variables for the compressible equations
130 NumVar_SFD += 2;
131 }
132 if (m_session->DefinesSolverInfo("HOMOGENEOUS"))
133 {
134 if (m_session->GetSolverInfo("HOMOGENEOUS") == "1D")
135 {
136 NumVar_SFD += 1;
137 }
138 }
139
140 // We store the time step
141 m_dt = m_equ[m_nequ - 1]->GetTimeStep();
142
143 // Evaluate optimum SFD parameters if dominent EV given by xml file
144 if (GrowthRateEV != 0.0 && FrequencyEV != 0.0)
145 {
146 std::cout << "Besed on the dominant EV given in the xml file,"
147 << "a 1D model is used to evaluate the optumum parameters"
148 << "of the SFD method:" << std::endl;
149 std::complex<NekDouble> EV = std::polar(exp(GrowthRateEV), FrequencyEV);
151 }
152
153 // We set up the elements of the operator of the encapsulated
154 // formulation of the selective frequencive damping method
156
157 // m_steps is set to 1. Then "m_equ[m_nequ - 1]->DoSolve()" will run
158 // for only one time step
159 m_equ[m_nequ - 1]->SetSteps(1);
160 std::ofstream m_file("ConvergenceHistory.txt",
161 std::ios::out | std::ios::trunc);
162
167
168 for (int i = 0; i < NumVar_SFD; ++i)
169 {
170 q0[i] = Array<OneD, NekDouble>(m_equ[m_nequ - 1]->GetTotPoints(),
171 0.0); // q0 is initialised
172 qBar0[i] =
173 Array<OneD, NekDouble>(m_equ[m_nequ - 1]->GetTotPoints(), 0.0);
174 m_equ[m_nequ - 1]->CopyFromPhysField(i, qBar0[i]);
175 }
176
177 /// Definition of variables used in this algorithm
178 m_stepCounter = 0;
179 m_Check = 0;
182 Diff_q_qBar = 1.0;
183 Diff_q1_q0 = 1.0;
184 cpuTime = 0.0;
185 elapsed = 0.0;
186 totalTime = 0.0;
188
189 while (std::max(Diff_q_qBar, Diff_q1_q0) > TOL)
190 {
191 /// Call the Navier-Stokes solver for one time step
192 m_equ[m_nequ - 1]->DoSolve();
193
194 for (int i = 0; i < NumVar_SFD; ++i)
195 {
196 /// Copy the current flow field into q0
197 m_equ[m_nequ - 1]->CopyFromPhysField(i, q0[i]);
198
199 /// Apply the linear operator to the outcome of the solver
200 ComputeSFD(i, q0, qBar0, q1, qBar1);
201
202 /// Update qBar
203 qBar0[i] = qBar1[i];
204
205 /// Copy the output of the SFD method into the current flow field
206 m_equ[m_nequ - 1]->CopyToPhysField(i, q1[i]);
207 }
208
209 if (m_infosteps && !((m_stepCounter + 1) % m_infosteps))
210 {
212
213 /// Loop for the adaptive SFD method
215 FlowPartiallyConverged == false)
216 {
217
219
221 {
222 if (m_comm->GetRank() == 0)
223 {
224 std::cout
225 << "\n\t The SFD method is converging: we compute "
226 << "stability analysis using the 'partially "
227 << "converged' steady state as base flow:\n"
228 << std::endl;
229 }
230
231 m_equ[m_nequ - 1]->Checkpoint_BaseFlow(m_Check_BaseFlow);
233
234 A->GetAdvObject()->SetBaseFlow(q0,
235 m_equ[0]->UpdateFields());
237
238 if (m_comm->GetRank() == 0)
239 {
240 // Compute the update of the SFD parameters only on
241 // one processor
243 }
244 else
245 {
246 // On all the other processors, the parameters are set
247 // to 0
248 m_X = 0;
249 m_Delta = 0;
250 }
251 // The we give to all the processors the value of X and
252 // Delta of the first processor
255
257 }
260 {
261 if (m_comm->GetRank() == 0)
262 {
263 std::cout << "\n\t We compute stability analysis using"
264 << " the current flow field as base flow:\n"
265 << std::endl;
266 }
267
268 m_equ[m_nequ - 1]->Checkpoint_BaseFlow(m_Check_BaseFlow);
270
271 A->GetAdvObject()->SetBaseFlow(q0,
272 m_equ[0]->UpdateFields());
274
275 if (m_comm->GetRank() == 0)
276 {
277 // Compute the update of the SFD parameters only on
278 // one processor
280 }
281 else
282 {
283 // On all the other processors, the parameters are set
284 // to 0
285 m_X = 0;
286 m_Delta = 0;
287 }
288 // The we give to all the processors the value of X and
289 // Delta of the first processor
292
294 }
295 }
296 }
297
299 (!((m_stepCounter + 1) % m_checksteps)))
300 {
301 m_Check++;
302 m_equ[m_nequ - 1]->Checkpoint_Output(m_Check);
303 }
305 }
306
307 m_file.close();
308
309 /// We save the final solution into a .fld file
310 m_equ[m_nequ - 1]->Output();
311
312 for (int j = 0; j < m_equ[m_nequ - 1]->GetNvariables(); ++j)
313 {
314 NekDouble vL2Error = m_equ[m_nequ - 1]->L2Error(j, false);
315 NekDouble vLinfError = m_equ[m_nequ - 1]->LinfError(j);
316 if (m_comm->GetRank() == 0)
317 {
318 out << "L 2 error (variable " << m_equ[m_nequ - 1]->GetVariable(j)
319 << ") : " << vL2Error << std::endl;
320 out << "L inf error (variable " << m_equ[m_nequ - 1]->GetVariable(j)
321 << ") : " << vLinfError << std::endl;
322 }
323 }
324}
325
326/**
327 * This routine defines the encapsulated SFD operator with first-order
328 * splitting and exact resolution of the second subproblem.
329 *(See http://scitation.aip.org/content/aip/journal/pof2/26/3/10.1063/1.4867482
330 *for details)
331 */
333 const NekDouble Delta_input)
334{
335 NekDouble X = X_input * m_dt;
336 NekDouble Delta = Delta_input / m_dt;
337 NekDouble coeff = 1.0 / (1.0 + X * Delta);
338 M11 = coeff * (1.0 + X * Delta * exp(-(X + 1.0 / Delta)));
339 M12 = coeff * (X * Delta * (1.0 - exp(-(X + 1.0 / Delta))));
340 M21 = coeff * (1.0 - exp(-(X + 1.0 / Delta)));
341 M22 = coeff * (X * Delta + exp(-(X + 1.0 / Delta)));
342}
343
344/**
345 *
346 */
348 const int i, const Array<OneD, const Array<OneD, NekDouble>> &q0,
349 const Array<OneD, const Array<OneD, NekDouble>> &qBar0,
352{
353 q1[i] = Array<OneD, NekDouble>(m_equ[m_nequ - 1]->GetTotPoints(), 0.0);
354 qBar1[i] = Array<OneD, NekDouble>(m_equ[m_nequ - 1]->GetTotPoints(), 0.0);
355
356 /// Encapsulated SFD method
357 Vmath::Svtvp(q1[i].size(), M11, q0[i], 1, q1[i], 1, q1[i], 1);
358 Vmath::Svtvp(q1[i].size(), M12, qBar0[i], 1, q1[i], 1, q1[i], 1);
359
360 Vmath::Svtvp(qBar1[i].size(), M21, q0[i], 1, qBar1[i], 1, qBar1[i], 1);
361 Vmath::Svtvp(qBar1[i].size(), M22, qBar0[i], 1, qBar1[i], 1, qBar1[i], 1);
362}
363
364/**
365 *
366 */
368{
369 NekDouble growthEV(0.0);
370 NekDouble frequencyEV(0.0);
371
372 // m_kdim is the dimension of Krylov subspace (defined in the xml file and
373 // used in DriverArnoldi.cpp)
374 ReadEVfile(m_kdim, growthEV, frequencyEV);
375
376 std::cout << "\n\tgrowthEV = " << growthEV << std::endl;
377 std::cout << "\tfrequencyEV = " << frequencyEV << std::endl;
378
379 std::complex<NekDouble> ApproxEV = std::polar(exp(growthEV), frequencyEV);
380
381 NekDouble X_new = m_X;
382 NekDouble Delta_new = m_Delta;
383
384 GradientDescentMethod(ApproxEV, X_new, Delta_new);
385
386 m_X = X_new;
387 m_Delta = Delta_new;
388
390}
391
392/**
393 * This routine implements a gradient descent method to find the parameters X
394 * end Delta which give the minimum eigenlavue of the SFD problem applied to
395 * the scalar case u(n+1) = \alpha*u(n).
396 */
398 const std::complex<NekDouble> &alpha, NekDouble &X_output,
399 NekDouble &Delta_output)
400{
401 std::cout << "\n\tWe enter the Gradient Descent Method [...]" << std::endl;
402 bool OptParmFound = false;
403 bool Descending = true;
404 NekDouble X_input = X_output;
405 NekDouble Delta_input = Delta_output;
406
407 NekDouble X_init = X_output;
408 NekDouble Delta_init = Delta_output;
409 int stepCounter(0);
410
411 NekDouble F0(0.0);
412 NekDouble F0x(0.0);
413 NekDouble F0y(0.0);
414 NekDouble F1(0.0);
415 NekDouble dx = 0.00000001;
416 NekDouble dirx(0.0);
417 NekDouble diry(0.0);
418 NekDouble s(0.0);
419 NekDouble CurrentMin = 1.0;
420
421 while (OptParmFound == false)
422 {
423 Descending = true;
424 EvalEV_ScalarSFD(X_input, Delta_input, alpha, F0);
425 EvalEV_ScalarSFD(X_input + dx, Delta_input, alpha, F0x);
426 EvalEV_ScalarSFD(X_input, Delta_input + dx, alpha, F0y);
427 dirx = (F0 - F0x);
428 diry = (F0 - F0y);
429 s = abs(0.000001 / dirx);
430 X_output = X_input + s * dirx;
431 Delta_output = Delta_input + s * diry;
432 F1 = F0;
433
434 while (Descending == true)
435 {
436 CurrentMin = F1;
437 X_input = X_output;
438 Delta_input = Delta_output;
439 EvalEV_ScalarSFD(X_output, Delta_output, alpha, F1);
440
441 if (F1 > CurrentMin)
442 {
443 Descending = false;
444 }
445 else
446 {
447 s = s + s * 0.01;
448 X_output = X_input + s * dirx;
449 Delta_output = Delta_input + s * diry;
450 }
451
452 if (stepCounter > 9999999)
453 {
454 // We are stuck in this loop..
455 // Then we restart it with different initail conditions
456 Descending = false;
457 X_input = X_init;
458 Delta_init = Delta_init + Delta_init * 0.1;
459 Delta_input = Delta_init;
460 stepCounter = 0;
461 }
462 stepCounter++;
463 }
464
465 if (abs(F0 - F1) < dx)
466 {
467 std::cout << "\tThe Gradient Descent Method has converged!"
468 << std::endl;
469 EvalEV_ScalarSFD(X_output, Delta_output, alpha, F1);
470 std::cout << "\n\tThe updated parameters are: X_tilde = "
471 << X_output << " and Delta_tilde = " << Delta_output
472 << std::endl;
473 OptParmFound = true;
474 }
475 }
476}
477
478/**
479 * This routine evaluates the maximum eigenvalue of the SFD system when applied
480 * to the 1D model u(n+1) = alpha*u(n)
481 */
483 const NekDouble &Delta_input,
484 const std::complex<NekDouble> &alpha,
485 NekDouble &MaxEV)
486{
487 NekDouble A11 = (1.0 + X_input * Delta_input *
488 std::exp(-(X_input + 1.0 / Delta_input))) /
489 (1.0 + X_input * Delta_input);
490 NekDouble A12 =
491 (X_input * Delta_input -
492 X_input * Delta_input * exp(-(X_input + 1.0 / Delta_input))) /
493 (1.0 + X_input * Delta_input);
494 NekDouble A21 = (1.0 - 1.0 * exp(-(X_input + 1.0 / Delta_input))) /
495 (1 + X_input * Delta_input);
496 NekDouble A22 =
497 (X_input * Delta_input + 1.0 * exp(-(X_input + 1.0 / Delta_input))) /
498 (1.0 + X_input * Delta_input);
499
500 std::complex<NekDouble> B11 = alpha;
501 NekDouble B12 = 0.0;
502 NekDouble B21 = 0.0;
503 NekDouble B22 = 1.0;
504
505 std::complex<NekDouble> a = A11 * B11 + A12 * B21;
506 NekDouble b = A11 * B12 + A12 * B22;
507 std::complex<NekDouble> c = A21 * B11 + A22 * B21;
508 NekDouble d = A21 * B12 + A22 * B22;
509
510 std::complex<NekDouble> delt = sqrt((a - d) * (a - d) + 4.0 * b * c);
511 std::complex<NekDouble> lambda_1 = 0.5 * (a + d + delt);
512 std::complex<NekDouble> lambda_2 = 0.5 * (a + d - delt);
513
514 NekDouble NORM_1 = abs(lambda_1);
515 NekDouble NORM_2 = abs(lambda_2);
516
517 MaxEV = std::max(NORM_1, NORM_2);
518}
519
520/**
521 *
522 */
523void DriverSteadyState::ReadEVfile(int &KrylovSubspaceDim, NekDouble &growthEV,
524 NekDouble &frequencyEV)
525{
526 // This routine reads the .evl file written by the Arnoldi algorithm
527 // (written in September 2014)
528 std::string line;
529 int NumLinesInFile = 0;
530 std::string EVfileName = m_session->GetSessionName() + ".evl";
531 std::ifstream EVfile(EVfileName.c_str());
532 ASSERTL0(EVfile.good(), "Cannot open .evl file.");
533
534 if (EVfile)
535 {
536 // This block counts the total number of lines of the .evl file
537 // We keep going util we reach the end of the file
538 while (getline(EVfile, line))
539 {
540 ++NumLinesInFile;
541 }
542
543 // It may happen that the Stability method that have produced the .elv
544 // file converges in less than m_kdim iterations. In this case,
545 // KrylovSubspaceDim has to be changed here
546 if (NumLinesInFile <
547 KrylovSubspaceDim * 2.0 +
548 KrylovSubspaceDim * (KrylovSubspaceDim + 1.0) / 2.0)
549 {
550 for (int i = 1; i <= KrylovSubspaceDim; ++i)
551 {
552 if (NumLinesInFile == i * 2.0 + i * (i + 1.0) / 2.0)
553 {
554 KrylovSubspaceDim = i;
555 }
556 }
557 }
558
559 // go back to the beginning of the file
560 EVfile.clear();
561 EVfile.seekg(0, std::ios::beg);
562
563 // We now want to go to the line where the most unstable eigenlavue was
564 // written
565 for (int i = 0; i < (NumLinesInFile - KrylovSubspaceDim); ++i)
566 {
567 std::getline(EVfile, line);
568 std::cout << "Discard line: " << line << std::endl;
569 }
570
571 std::vector<std::string> tokens;
572 std::getline(EVfile, line);
573 boost::algorithm::split(tokens, line, boost::is_any_of("\t "),
574 boost::token_compress_on);
575
576 ASSERTL0(tokens.size() >= 5,
577 "Unexpected formatting of .evl file while reading line:\n" +
578 line);
579 growthEV = std::stod(tokens[4]);
580 frequencyEV = std::stod(tokens[5]);
581 }
582 else
583 {
584 std::cout << "An error occurred when opening the .evl file"
585 << std::endl;
586 }
587 EVfile.close();
588}
589
590/**
591 * This routine evaluates |q-qBar|_inf (and |q1-q0|_inf) and writes the values
592 * in "ConvergenceHistory.txt"
593 */
595 const Array<OneD, const Array<OneD, NekDouble>> &qBar1,
596 const Array<OneD, const Array<OneD, NekDouble>> &q0,
597 NekDouble &MaxNormDiff_q_qBar, NekDouble &MaxNormDiff_q1_q0)
598{
599 Array<OneD, NekDouble> NormDiff_q_qBar(NumVar_SFD, 1.0);
600 Array<OneD, NekDouble> NormDiff_q1_q0(NumVar_SFD, 1.0);
601 MaxNormDiff_q_qBar = 0.0;
602 MaxNormDiff_q1_q0 = 0.0;
603
604 for (int i = 0; i < NumVar_SFD; ++i)
605 {
606 NormDiff_q_qBar[i] = m_equ[m_nequ - 1]->LinfError(i, qBar1[i]);
607 NormDiff_q1_q0[i] = m_equ[m_nequ - 1]->LinfError(i, q0[i]);
608
609 if (MaxNormDiff_q_qBar < NormDiff_q_qBar[i])
610 {
611 MaxNormDiff_q_qBar = NormDiff_q_qBar[i];
612 }
613 if (MaxNormDiff_q1_q0 < NormDiff_q1_q0[i])
614 {
615 MaxNormDiff_q1_q0 = NormDiff_q1_q0[i];
616 }
617 }
618
619 timer.Stop();
621 cpuTime += elapsed;
623
624 if (m_comm->GetRank() == 0)
625 {
626 std::cout << "SFD - Step: " << std::left << m_stepCounter + 1
627 << ";\tTime: " << std::left << m_equ[m_nequ - 1]->GetTime()
628 << ";\tCPU time = " << std::left << cpuTime << " s"
629 << ";\tTot time = " << std::left << totalTime << " s"
630 << ";\tX = " << std::left << m_X << ";\tDelta = " << std::left
631 << m_Delta << ";\t|q-qBar|inf = " << std::left
632 << MaxNormDiff_q_qBar << std::endl;
633 std::ofstream m_file("ConvergenceHistory.txt", std::ios_base::app);
634 m_file << m_stepCounter + 1 << "\t" << m_equ[m_nequ - 1]->GetTime()
635 << "\t" << totalTime << "\t" << MaxNormDiff_q_qBar << "\t"
636 << MaxNormDiff_q1_q0 << std::endl;
637 m_file.close();
638 }
639
640 cpuTime = 0.0;
641 timer.Start();
642}
643
644/**
645 *
646 */
648{
649 std::cout << "\n====================================="
650 "=================================="
651 << std::endl;
652 std::cout << "Parameters for the SFD method:" << std::endl;
653 std::cout << "\tControl Coefficient: X = " << m_X << std::endl;
654 std::cout << "\tFilter Width: Delta = " << m_Delta << std::endl;
655 std::cout << "The simulation is stopped when |q-qBar|inf < " << TOL
656 << std::endl;
658 {
659 std::cout << "\nWe use the adaptive SFD method:" << std::endl;
660 std::cout << " The parameters are updated every " << AdaptiveTime
661 << " time units;" << std::endl;
662 std::cout << " until |q-qBar|inf becomes smaller than " << AdaptiveTOL
663 << std::endl;
664 }
665 std::cout << "====================================="
666 "==================================\n"
667 << std::endl;
668}
669
670} // namespace Nektar::SolverUtils
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
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.
NekDouble TimePerTest(unsigned int n)
Returns amount of seconds per iteration in a test with n iterations.
Definition: Timer.cpp:65
A base class for PDEs which include an advection component.
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 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 ComputeSFD(const int i, const Array< OneD, const Array< OneD, NekDouble > > &q0, const Array< OneD, const Array< OneD, NekDouble > > &qBar0, Array< OneD, Array< OneD, NekDouble > > &q1, Array< OneD, Array< OneD, NekDouble > > &qBar1)
SOLVER_UTILS_EXPORT void v_InitObject(std::ostream &out=std::cout) override
Initialises EquationSystem class members.
void ConvergenceHistory(const Array< OneD, const Array< OneD, NekDouble > > &qBar1, const Array< OneD, const Array< OneD, NekDouble > > &q0, NekDouble &MaxNormDiff_q_qBar, NekDouble &MaxNormDiff_q1_q0)
static DriverSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
SOLVER_UTILS_EXPORT DriverSteadyState(const LibUtilities::SessionReaderSharedPtr pSession, const SpatialDomains::MeshGraphSharedPtr pGraph)
Constructor.
NekDouble m_Delta
Definition of the SFD parameters:
void SetSFDOperator(const NekDouble X_input, const NekDouble Delta_input)
NekDouble M11
Definition of the SFD operator.
SOLVER_UTILS_EXPORT void v_Execute(std::ostream &out=std::cout) override
Virtual function for solve implementation.
void GradientDescentMethod(const std::complex< NekDouble > &alpha, NekDouble &X_output, NekDouble &Delta_output)
static std::string className
Name of the class.
void EvalEV_ScalarSFD(const NekDouble &X_input, const NekDouble &Delta_input, const std::complex< NekDouble > &alpha, NekDouble &MaxEV)
bool FlowPartiallyConverged
For adaptive SFD method.
void ReadEVfile(int &KrylovSubspaceDim, NekDouble &growthEV, NekDouble &frequencyEV)
std::shared_ptr< SessionReader > SessionReaderSharedPtr
DriverFactory & GetDriverFactory()
Definition: Driver.cpp:64
std::shared_ptr< AdvectionSystem > AdvectionSystemSharedPtr
Shared pointer to an AdvectionSystem class.
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:174
std::vector< double > d(NPUPPER *NPUPPER)
double NekDouble
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
scalarT< T > abs(scalarT< T > in)
Definition: scalar.hpp:289
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:285