46 string DriverSteadyState::className
48 "SteadyState", DriverSteadyState::create);
49 string DriverSteadyState::driverLookupId
50 = LibUtilities::SessionReader::RegisterEnumValue(
51 "Driver",
"SteadyState",0);
56 DriverSteadyState::DriverSteadyState(
111 if (
m_comm->GetRank() == 0)
123 if (
m_session->GetSolverInfo(
"EqType") ==
"EulerCFE" ||
124 m_session->GetSolverInfo(
"EqType") ==
"NavierStokesCFE")
129 if(
m_session->DefinesSolverInfo(
"HOMOGENEOUS"))
131 if (
m_session->GetSolverInfo(
"HOMOGENEOUS") ==
"1D")
143 cout <<
"Besed on the dominant EV given in the xml file,"
144 <<
"a 1D model is used to evaluate the optumum parameters"
145 <<
"of the SFD method:" << endl;
158 ofstream
m_file(
"ConvergenceHistory.txt", ios::out | ios::trunc);
219 if (
m_comm->GetRank() == 0)
221 cout <<
"\n\t The SFD method is converging: we compute "
222 <<
"stability analysis using the 'partially "
223 <<
"converged' steady state as base flow:\n" << endl;
229 A->GetAdvObject()->SetBaseFlow(q0,
m_equ[0]->UpdateFields());
232 if (
m_comm->GetRank() == 0)
255 if (
m_comm->GetRank() == 0)
257 cout <<
"\n\t We compute stability analysis using"
258 <<
" the current flow field as base flow:\n"
265 A->GetAdvObject()->SetBaseFlow(q0,
m_equ[0]->UpdateFields());
268 if (
m_comm->GetRank() == 0)
304 for(
int j = 0; j <
m_equ[
m_nequ - 1]->GetNvariables(); ++j)
308 if (
m_comm->GetRank() == 0)
310 out <<
"L 2 error (variable " <<
m_equ[
m_nequ - 1]->GetVariable(j)
311 <<
") : " << vL2Error << endl;
312 out <<
"L inf error (variable " <<
m_equ[
m_nequ - 1]->GetVariable(j)
313 <<
") : " << vLinfError << endl;
330 M11 = coeff*(1.0 + X*Delta*exp(-(X + 1.0/Delta)));
331 M12 = coeff*(X*Delta*(1.0 - exp(-(X + 1.0/Delta))));
332 M21 = coeff*(1.0 - exp(-(X + 1.0/Delta)));
333 M22 = coeff*(X*Delta + exp(-(X + 1.0/Delta)));
347 Vmath::Svtvp(q1[i].num_elements(),
M11, q0[i], 1, q1[i], 1, q1[i], 1 );
348 Vmath::Svtvp(q1[i].num_elements(),
M12, qBar0[i], 1, q1[i], 1, q1[i], 1 );
366 cout <<
"\n\tgrowthEV = " << growthEV << endl;
367 cout <<
"\tfrequencyEV = " << frequencyEV << endl;
369 complex<NekDouble> ApproxEV = polar(exp(growthEV), frequencyEV);
388 const complex<NekDouble> &alpha,
392 cout <<
"\n\tWe enter the Gradient Descent Method [...]" << endl;
393 bool OptParmFound =
false;
394 bool Descending =
true;
412 while (OptParmFound ==
false)
420 s = abs(0.000001/dirx);
421 X_output = X_input + s*dirx;
422 Delta_output = Delta_input + s*diry;
425 while (Descending ==
true)
429 Delta_input = Delta_output;
439 X_output = X_input + s*dirx;
440 Delta_output = Delta_input + s*diry;
443 if(stepCounter > 9999999)
449 Delta_init = Delta_init + Delta_init*0.1;
450 Delta_input = Delta_init;
458 cout <<
"\tThe Gradient Descent Method has converged!" << endl;
460 cout <<
"\n\tThe updated parameters are: X_tilde = " << X_output
461 <<
" and Delta_tilde = " << Delta_output << endl;
476 const complex<NekDouble> &alpha,
479 NekDouble A11 = ( 1.0 + X_input * Delta_input *
480 exp(-(X_input + 1.0/Delta_input)) )/(1.0 + X_input*Delta_input);
481 NekDouble A12 = ( X_input*Delta_input - X_input * Delta_input *
482 exp(-(X_input + 1.0/Delta_input)) )/(1.0 + X_input*Delta_input);
484 exp(-(X_input + 1.0/Delta_input)) )/(1 + X_input*Delta_input);
485 NekDouble A22 = ( X_input*Delta_input + 1.0 *
486 exp(-(X_input + 1.0/Delta_input)) )/(1.0 + X_input*Delta_input);
488 complex<NekDouble> B11 = alpha;
493 complex<NekDouble> a = A11*B11 + A12*B21;
495 complex<NekDouble> c = A21*B11 + A22*B21;
498 complex<NekDouble> delt = sqrt((a-d)*(a-d) + 4.0*b*c);
499 complex<NekDouble> lambda_1 = 0.5*(a+d + delt);
500 complex<NekDouble> lambda_2 = 0.5*(a+d - delt);
505 MaxEV = max(NORM_1, NORM_2);
510 int &KrylovSubspaceDim,
517 int NumLinesInFile = 0;
518 std::string EVfileName =
m_session->GetSessionName() +
".evl";
519 std::ifstream EVfile(EVfileName.c_str());
520 ASSERTL0(EVfile.good(),
"Cannot open .evl file.");
526 while(getline(EVfile, line))
534 if(NumLinesInFile < KrylovSubspaceDim*2.0
535 + KrylovSubspaceDim*(KrylovSubspaceDim+1.0)/2.0)
537 for(
int i = 1; i <= KrylovSubspaceDim; ++i)
539 if(NumLinesInFile == i*2.0 + i*(i+1.0)/2.0)
541 KrylovSubspaceDim = i;
548 EVfile.seekg(0, ios::beg);
552 for(
int i = 0; i < (NumLinesInFile - KrylovSubspaceDim); ++i)
554 std::getline(EVfile, line);
555 cout <<
"Discard line: " << line << endl;
558 std::vector<std::string> tokens;
559 std::getline(EVfile, line);
560 boost::algorithm::split(tokens, line,
561 boost::is_any_of(
"\t "),boost::token_compress_on);
564 "Unexpected formatting of .evl file while reading line:\n"
566 growthEV = boost::lexical_cast<
NekDouble>(tokens[4]);
567 frequencyEV = boost::lexical_cast<
NekDouble>(tokens[5]);
571 cout <<
"An error occurred when opening the .evl file" << endl;
589 MaxNormDiff_q_qBar=0.0;
590 MaxNormDiff_q1_q0=0.0;
594 NormDiff_q_qBar[i] =
m_equ[
m_nequ - 1]->LinfError(i, qBar1[i]);
595 NormDiff_q1_q0[i] =
m_equ[
m_nequ - 1]->LinfError(i, q0[i]);
597 if (MaxNormDiff_q_qBar < NormDiff_q_qBar[i])
599 MaxNormDiff_q_qBar = NormDiff_q_qBar[i];
601 if (MaxNormDiff_q1_q0 < NormDiff_q1_q0[i])
603 MaxNormDiff_q1_q0 = NormDiff_q1_q0[i];
612 if (
m_comm->GetRank() == 0)
615 <<
";\tTime: " << left <<
m_equ[
m_nequ - 1]->GetFinalTime()
616 <<
";\tCPU time = " << left <<
cpuTime <<
" s"
617 <<
";\tTot time = " << left <<
totalTime <<
" s"
618 <<
";\tX = " << left <<
m_X
619 <<
";\tDelta = " << left <<
m_Delta
620 <<
";\t|q-qBar|inf = " << left << MaxNormDiff_q_qBar << endl;
621 std::ofstream
m_file(
"ConvergenceHistory.txt", std::ios_base::app);
625 << MaxNormDiff_q_qBar <<
"\t"
626 << MaxNormDiff_q1_q0 << endl;
637 cout <<
"\n====================================="
638 "=================================="
640 cout <<
"Parameters for the SFD method:" << endl;
641 cout <<
"\tControl Coefficient: X = " <<
m_X << endl;
642 cout <<
"\tFilter Width: Delta = " <<
m_Delta << endl;
643 cout <<
"The simulation is stopped when |q-qBar|inf < " <<
TOL << endl;
646 cout <<
"\nWe use the adaptive SFD method:" << endl;
647 cout <<
" The parameters are updated every " <<
AdaptiveTime
648 <<
" time units;" << endl;
649 cout <<
" until |q-qBar|inf becomes smaller than " <<
AdaptiveTOL
652 cout <<
"====================================="
653 "==================================\n" << endl;
#define ASSERTL0(condition, msg)
boost::shared_ptr< AdvectionSystem > AdvectionSystemSharedPtr
Shared pointer to an AdvectionSystem class.
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 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
virtual SOLVER_UTILS_EXPORT void v_InitObject(std::ostream &out=std::cout)
Second-stage initialisation.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
virtual void v_InitObject(std::ostream &out=std::cout)
Virtual function for initialisation 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)
int m_NonConvergingStepsCounter
void ReadEVfile(int &KrylovSubspaceDim, NekDouble &growthEV, NekDouble &frequencyEV)
enum EvolutionOperatorType m_EvolutionOperator
Evolution Operator.
void GradientDescentMethod(const std::complex< NekDouble > &alpha, NekDouble &X_output, NekDouble &Delta_output)
void ComputeOptimization()
void EvalEV_ScalarSFD(const NekDouble &X_input, const NekDouble &Delta_input, const std::complex< NekDouble > &alpha, NekDouble &MaxEV)
virtual SOLVER_UTILS_EXPORT ~DriverSteadyState()
Destructor.
void SetSFDOperator(const NekDouble X_input, const NekDouble Delta_input)
Array< OneD, EquationSystemSharedPtr > m_equ
Equation system to solve.
NekDouble M11
Definition of the SFD operator.
DriverFactory & GetDriverFactory()
LibUtilities::CommSharedPtr m_comm
Communication object.
NekDouble m_Delta
Definition of the SFD parameters:
virtual void v_Execute(std::ostream &out=std::cout)
Virtual function for solve implementation.
virtual SOLVER_UTILS_EXPORT void v_Execute(std::ostream &out=std::cout)
Virtual function for solve implementation.
A base class for PDEs which include an advection component.
NekDouble TimePerTest(unsigned int n)
Returns amount of seconds per iteration in a test with n iterations.
LibUtilities::SessionReaderSharedPtr m_session
Session reader object.
int m_nequ
number of equations
bool FlowPartiallyConverged
For adaptive SFD method.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.