38#include <boost/algorithm/string/classification.hpp>
39#include <boost/algorithm/string/predicate.hpp>
40#include <boost/algorithm/string/split.hpp>
102 if (
m_comm->GetRank() == 0)
114 if (
m_session->GetSolverInfo(
"EqType") ==
"UnsteadyNavierStokes" ||
115 m_session->GetSolverInfo(
"EqType") ==
"SteadyNavierStokes")
117 int nConvectiveFields =
m_session->GetVariables().size();
118 if (boost::iequals(
m_session->GetVariable(nConvectiveFields - 1),
"p"))
120 nConvectiveFields -= 1;
125 if (
m_session->GetSolverInfo(
"EqType") ==
"EulerCFE" ||
126 m_session->GetSolverInfo(
"EqType") ==
"NavierStokesCFE" ||
127 m_session->GetSolverInfo(
"EqType") ==
"NavierStokesImplicitCFE")
132 if (
m_session->DefinesSolverInfo(
"HOMOGENEOUS"))
134 if (
m_session->GetSolverInfo(
"HOMOGENEOUS") ==
"1D")
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;
160 std::ofstream
m_file(
"ConvergenceHistory.txt",
161 std::ios::out | std::ios::trunc);
222 if (
m_comm->GetRank() == 0)
225 <<
"\n\t The SFD method is converging: we compute "
226 <<
"stability analysis using the 'partially "
227 <<
"converged' steady state as base flow:\n"
234 A->GetAdvObject()->SetBaseFlow(q0,
235 m_equ[0]->UpdateFields());
238 if (
m_comm->GetRank() == 0)
261 if (
m_comm->GetRank() == 0)
263 std::cout <<
"\n\t We compute stability analysis using"
264 <<
" the current flow field as base flow:\n"
271 A->GetAdvObject()->SetBaseFlow(q0,
272 m_equ[0]->UpdateFields());
275 if (
m_comm->GetRank() == 0)
312 for (
int j = 0; j <
m_equ[
m_nequ - 1]->GetNvariables(); ++j)
316 if (
m_comm->GetRank() == 0)
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;
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)));
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);
376 std::cout <<
"\n\tgrowthEV = " << growthEV << std::endl;
377 std::cout <<
"\tfrequencyEV = " << frequencyEV << std::endl;
379 std::complex<NekDouble> ApproxEV = std::polar(exp(growthEV), frequencyEV);
398 const std::complex<NekDouble> &alpha,
NekDouble &X_output,
401 std::cout <<
"\n\tWe enter the Gradient Descent Method [...]" << std::endl;
402 bool OptParmFound =
false;
403 bool Descending =
true;
421 while (OptParmFound ==
false)
429 s =
abs(0.000001 / dirx);
430 X_output = X_input + s * dirx;
431 Delta_output = Delta_input + s * diry;
434 while (Descending ==
true)
438 Delta_input = Delta_output;
448 X_output = X_input + s * dirx;
449 Delta_output = Delta_input + s * diry;
452 if (stepCounter > 9999999)
458 Delta_init = Delta_init + Delta_init * 0.1;
459 Delta_input = Delta_init;
465 if (
abs(F0 - F1) < dx)
467 std::cout <<
"\tThe Gradient Descent Method has converged!"
470 std::cout <<
"\n\tThe updated parameters are: X_tilde = "
471 << X_output <<
" and Delta_tilde = " << Delta_output
484 const std::complex<NekDouble> &alpha,
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);
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);
497 (X_input * Delta_input + 1.0 * exp(-(X_input + 1.0 / Delta_input))) /
498 (1.0 + X_input * Delta_input);
500 std::complex<NekDouble> B11 = alpha;
505 std::complex<NekDouble> a = A11 * B11 + A12 * B21;
507 std::complex<NekDouble> c = A21 * B11 + A22 * B21;
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);
517 MaxEV = std::max(NORM_1, NORM_2);
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.");
538 while (getline(EVfile, line))
547 KrylovSubspaceDim * 2.0 +
548 KrylovSubspaceDim * (KrylovSubspaceDim + 1.0) / 2.0)
550 for (
int i = 1; i <= KrylovSubspaceDim; ++i)
552 if (NumLinesInFile == i * 2.0 + i * (i + 1.0) / 2.0)
554 KrylovSubspaceDim = i;
561 EVfile.seekg(0, std::ios::beg);
565 for (
int i = 0; i < (NumLinesInFile - KrylovSubspaceDim); ++i)
567 std::getline(EVfile, line);
568 std::cout <<
"Discard line: " << line << std::endl;
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);
577 "Unexpected formatting of .evl file while reading line:\n" +
579 growthEV = std::stod(tokens[4]);
580 frequencyEV = std::stod(tokens[5]);
584 std::cout <<
"An error occurred when opening the .evl file"
601 MaxNormDiff_q_qBar = 0.0;
602 MaxNormDiff_q1_q0 = 0.0;
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]);
609 if (MaxNormDiff_q_qBar < NormDiff_q_qBar[i])
611 MaxNormDiff_q_qBar = NormDiff_q_qBar[i];
613 if (MaxNormDiff_q1_q0 < NormDiff_q1_q0[i])
615 MaxNormDiff_q1_q0 = NormDiff_q1_q0[i];
624 if (
m_comm->GetRank() == 0)
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);
635 <<
"\t" <<
totalTime <<
"\t" << MaxNormDiff_q_qBar <<
"\t"
636 << MaxNormDiff_q1_q0 << std::endl;
649 std::cout <<
"\n====================================="
650 "=================================="
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
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
665 std::cout <<
"====================================="
666 "==================================\n"
#define ASSERTL0(condition, msg)
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.
A base class for PDEs which include an advection component.
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.
int m_nequ
number of equations
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.
LibUtilities::Timer timer
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)
int m_NonConvergingStepsCounter
static std::string driverLookupId
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)
void ComputeOptimization()
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()
std::shared_ptr< AdvectionSystem > AdvectionSystemSharedPtr
Shared pointer to an AdvectionSystem class.
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
std::vector< double > d(NPUPPER *NPUPPER)
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.
scalarT< T > abs(scalarT< T > in)
scalarT< T > sqrt(scalarT< T > in)