46 string DriverSteadyState::className
48 "SteadyState", DriverSteadyState::create);
49 string DriverSteadyState::driverLookupId
50 = LibUtilities::SessionReader::RegisterEnumValue(
51 "Driver",
"SteadyState",0);
56 DriverSteadyState::DriverSteadyState(
112 if (
m_comm->GetRank() == 0)
124 if (
m_session->GetSolverInfo(
"EqType") ==
"EulerCFE" ||
125 m_session->GetSolverInfo(
"EqType") ==
"NavierStokesCFE")
130 if(
m_session->DefinesSolverInfo(
"HOMOGENEOUS"))
132 if (
m_session->GetSolverInfo(
"HOMOGENEOUS") ==
"1D")
144 cout <<
"Besed on the dominant EV given in the xml file,"
145 <<
"a 1D model is used to evaluate the optumum parameters"
146 <<
"of the SFD method:" << endl;
159 ofstream
m_file(
"ConvergenceHistory.txt", ios::out | ios::trunc);
220 if (
m_comm->GetRank() == 0)
222 cout <<
"\n\t The SFD method is converging: we compute "
223 <<
"stability analysis using the 'partially "
224 <<
"converged' steady state as base flow:\n" << endl;
230 A->GetAdvObject()->SetBaseFlow(q0);
233 if (
m_comm->GetRank() == 0)
256 if (
m_comm->GetRank() == 0)
258 cout <<
"\n\t We compute stability analysis using"
259 <<
" the current flow field as base flow:\n"
266 A->GetAdvObject()->SetBaseFlow(q0);
269 if (
m_comm->GetRank() == 0)
305 for(
int j = 0; j <
m_equ[
m_nequ - 1]->GetNvariables(); ++j)
309 if (
m_comm->GetRank() == 0)
311 out <<
"L 2 error (variable " <<
m_equ[
m_nequ - 1]->GetVariable(j)
312 <<
") : " << vL2Error << endl;
313 out <<
"L inf error (variable " <<
m_equ[
m_nequ - 1]->GetVariable(j)
314 <<
") : " << vLinfError << endl;
331 M11 = coeff*(1.0 + X*Delta*exp(-(X + 1.0/Delta)));
332 M12 = coeff*(X*Delta*(1.0 - exp(-(X + 1.0/Delta))));
333 M21 = coeff*(1.0 - exp(-(X + 1.0/Delta)));
334 M22 = coeff*(X*Delta + exp(-(X + 1.0/Delta)));
348 Vmath::Svtvp(q1[i].num_elements(),
M11, q0[i], 1, q1[i], 1, q1[i], 1 );
349 Vmath::Svtvp(q1[i].num_elements(),
M12, qBar0[i], 1, q1[i], 1, q1[i], 1 );
367 cout <<
"\n\tgrowthEV = " << growthEV << endl;
368 cout <<
"\tfrequencyEV = " << frequencyEV << endl;
370 complex<NekDouble> ApproxEV = polar(exp(growthEV), frequencyEV);
389 const complex<NekDouble> &alpha,
393 cout <<
"\n\tWe enter the Gradient Descent Method [...]" << endl;
394 bool OptParmFound =
false;
395 bool Descending =
true;
413 while (OptParmFound ==
false)
421 s = abs(0.000001/dirx);
422 X_output = X_input + s*dirx;
423 Delta_output = Delta_input + s*diry;
426 while (Descending ==
true)
430 Delta_input = Delta_output;
440 X_output = X_input + s*dirx;
441 Delta_output = Delta_input + s*diry;
444 if(stepCounter > 9999999)
450 Delta_init = Delta_init + Delta_init*0.1;
451 Delta_input = Delta_init;
459 cout <<
"\tThe Gradient Descent Method has converged!" << endl;
461 cout <<
"\n\tThe updated parameters are: X_tilde = " << X_output
462 <<
" and Delta_tilde = " << Delta_output << endl;
477 const complex<NekDouble> &alpha,
480 NekDouble A11 = ( 1.0 + X_input * Delta_input *
481 exp(-(X_input + 1.0/Delta_input)) )/(1.0 + X_input*Delta_input);
482 NekDouble A12 = ( X_input*Delta_input - X_input * Delta_input *
483 exp(-(X_input + 1.0/Delta_input)) )/(1.0 + X_input*Delta_input);
485 exp(-(X_input + 1.0/Delta_input)) )/(1 + X_input*Delta_input);
486 NekDouble A22 = ( X_input*Delta_input + 1.0 *
487 exp(-(X_input + 1.0/Delta_input)) )/(1.0 + X_input*Delta_input);
489 complex<NekDouble> B11 = alpha;
494 complex<NekDouble> a = A11*B11 + A12*B21;
496 complex<NekDouble> c = A21*B11 + A22*B21;
499 complex<NekDouble> delt = sqrt((a-d)*(a-d) + 4.0*b*c);
500 complex<NekDouble> lambda_1 = 0.5*(a+d + delt);
501 complex<NekDouble> lambda_2 = 0.5*(a+d - delt);
506 MaxEV = max(NORM_1, NORM_2);
511 int &KrylovSubspaceDim,
518 int NumLinesInFile = 0;
519 std::string EVfileName =
m_session->GetSessionName() +
".evl";
520 std::ifstream EVfile(EVfileName.c_str());
521 ASSERTL0(EVfile.good(),
"Cannot open .evl file.");
527 while(getline(EVfile, line))
535 if(NumLinesInFile < KrylovSubspaceDim*2.0
536 + KrylovSubspaceDim*(KrylovSubspaceDim+1.0)/2.0)
538 for(
int i = 1; i <= KrylovSubspaceDim; ++i)
540 if(NumLinesInFile == i*2.0 + i*(i+1.0)/2.0)
542 KrylovSubspaceDim = i;
549 EVfile.seekg(0, ios::beg);
553 for(
int i = 0; i < (NumLinesInFile - KrylovSubspaceDim); ++i)
555 std::getline(EVfile, line);
556 cout <<
"Discard line: " << line << endl;
559 std::vector<std::string> tokens;
560 std::getline(EVfile, line);
561 boost::algorithm::split(tokens, line,
562 boost::is_any_of(
"\t "),boost::token_compress_on);
565 "Unexpected formatting of .evl file while reading line:\n"
567 growthEV = boost::lexical_cast<
NekDouble>(tokens[4]);
568 frequencyEV = boost::lexical_cast<
NekDouble>(tokens[5]);
572 cout <<
"An error occurred when opening the .evl file" << endl;
590 MaxNormDiff_q_qBar=0.0;
591 MaxNormDiff_q1_q0=0.0;
595 NormDiff_q_qBar[i] =
m_equ[
m_nequ - 1]->LinfError(i, qBar1[i]);
596 NormDiff_q1_q0[i] =
m_equ[
m_nequ - 1]->LinfError(i, q0[i]);
598 if (MaxNormDiff_q_qBar < NormDiff_q_qBar[i])
600 MaxNormDiff_q_qBar = NormDiff_q_qBar[i];
602 if (MaxNormDiff_q1_q0 < NormDiff_q1_q0[i])
604 MaxNormDiff_q1_q0 = NormDiff_q1_q0[i];
613 if (
m_comm->GetRank() == 0)
616 <<
";\tTime: " << left <<
m_equ[
m_nequ - 1]->GetFinalTime()
617 <<
";\tCPU time = " << left <<
cpuTime <<
" s"
618 <<
";\tTot time = " << left <<
totalTime <<
" s"
619 <<
";\tX = " << left <<
m_X
620 <<
";\tDelta = " << left <<
m_Delta
621 <<
";\t|q-qBar|inf = " << left << MaxNormDiff_q_qBar << endl;
622 std::ofstream
m_file(
"ConvergenceHistory.txt", std::ios_base::app);
626 << MaxNormDiff_q_qBar <<
"\t"
627 << MaxNormDiff_q1_q0 << endl;
638 cout <<
"\n====================================="
639 "=================================="
641 cout <<
"Parameters for the SFD method:" << endl;
642 cout <<
"\tControl Coefficient: X = " <<
m_X << endl;
643 cout <<
"\tFilter Width: Delta = " <<
m_Delta << endl;
644 cout <<
"The simulation is stopped when |q-qBar|inf < " <<
TOL << endl;
647 cout <<
"\nWe use the adaptive SFD method:" << endl;
648 cout <<
" The parameters are updated every " <<
AdaptiveTime
649 <<
" time units;" << endl;
650 cout <<
" until |q-qBar|inf becomes smaller than " <<
AdaptiveTOL
653 cout <<
"====================================="
654 "==================================\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.