Forked from
Nektar / Nektar
9675 commits behind the upstream repository.
-
Chris Cantwell authoredChris Cantwell authored
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
DriverSteadyState.cpp 22.11 KiB
///////////////////////////////////////////////////////////////////////////////
//
// File DriverSteadyState.cpp
//
// For more information, please see: http://www.nektar.info
//
// The MIT License
//
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// License for the specific language governing rights and limitations under
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: Incompressible Navier Stokes solver
//
///////////////////////////////////////////////////////////////////////////////
#include <SolverUtils/DriverSteadyState.h>
#include <SolverUtils/AdvectionSystem.h>
namespace Nektar
{
namespace SolverUtils
{
string DriverSteadyState::className
= GetDriverFactory().RegisterCreatorFunction(
"SteadyState", DriverSteadyState::create);
string DriverSteadyState::driverLookupId
= LibUtilities::SessionReader::RegisterEnumValue(
"Driver","SteadyState",0);
/**
*
*/
DriverSteadyState::DriverSteadyState(
const LibUtilities::SessionReaderSharedPtr pSession)
: DriverModifiedArnoldi(pSession)
{
}
/**
*
*/
DriverSteadyState:: ~DriverSteadyState()
{
}
/**
* */
void DriverSteadyState::v_InitObject(ostream &out)
{
DriverModifiedArnoldi::v_InitObject(out);
}
void DriverSteadyState::v_Execute(ostream &out)
{
// With a loop over "DoSolve", this Driver implements the
// "encaplulated" Selective Frequency Damping method
// to find the steady state of a flow above the critical Reynolds
// number.
m_equ[m_nequ - 1]->PrintSummary(out);
m_equ[m_nequ - 1]->DoInitialise();
m_session->LoadParameter("IO_InfoSteps", m_infosteps, 1000);
m_session->LoadParameter("IO_CheckSteps", m_checksteps, 100000);
m_session->LoadParameter("ControlCoeff",m_X, 1);
m_session->LoadParameter("FilterWidth", m_Delta, 2);
// To evaluate optimum SFD parameters if growth rate provided in the
// xml file
m_session->LoadParameter("GrowthRateEV", GrowthRateEV, 0.0);
// To evaluate optimum SFD parameters if frequency provided in the xml
// file
m_session->LoadParameter("FrequencyEV", FrequencyEV, 0.0);
// Determine when SFD method is converged
m_session->LoadParameter("TOL", TOL, 1.0e-08);
// Used only for the Adaptive SFD method
m_session->LoadParameter("AdaptiveTOL", AdaptiveTOL, 1.0e-02);
// Used only for the Adaptive SFD method
m_session->LoadParameter("AdaptiveTime", AdaptiveTime, 50.0*m_Delta);
if (m_comm->GetRank() == 0)
{
PrintSummarySFD();
}
timer.Start();
// Definition of shared pointer (used only for the Adaptive SFD method)
AdvectionSystemSharedPtr A = m_equ[0]->as<AdvectionSystem>();
// Condition necessary to run SFD for the compressible case
NumVar_SFD = m_equ[m_nequ - 1]->UpdateFields()[0]->GetCoordim(0);
if (m_session->GetSolverInfo("EqType") == "EulerCFE" ||
m_session->GetSolverInfo("EqType") == "NavierStokesCFE")
{
// Number of variables for the compressible equations
NumVar_SFD += 2;
}
if(m_session->DefinesSolverInfo("HOMOGENEOUS"))
{
if (m_session->GetSolverInfo("HOMOGENEOUS") == "1D")
{
NumVar_SFD += 1;
}
}
// We store the time step
m_dt = m_equ[m_nequ - 1]->GetTimeStep();
// Evaluate optimum SFD parameters if dominent EV given by xml file
if (GrowthRateEV != 0.0 && FrequencyEV != 0.0) {
cout << "Besed on the dominant EV given in the xml file,"
<< "a 1D model is used to evaluate the optumum parameters"
<< "of the SFD method:" << endl;
complex<NekDouble> EV = polar(exp(GrowthRateEV), FrequencyEV);
GradientDescentMethod(EV, m_X, m_Delta);
}
// We set up the elements of the operator of the encapsulated
// formulation of the selective frequencive damping method
SetSFDOperator(m_X, m_Delta);
// m_steps is set to 1. Then "m_equ[m_nequ - 1]->DoSolve()" will run
// for only one time step
m_equ[m_nequ - 1]->SetStepsToOne();
ofstream m_file("ConvergenceHistory.txt", ios::out | ios::trunc);
Array<OneD, Array<OneD, NekDouble> > q0(NumVar_SFD);
Array<OneD, Array<OneD, NekDouble> > q1(NumVar_SFD);
Array<OneD, Array<OneD, NekDouble> > qBar0(NumVar_SFD);
Array<OneD, Array<OneD, NekDouble> > qBar1(NumVar_SFD);
for(int i = 0; i < NumVar_SFD; ++i)
{
q0[i] = Array<OneD, NekDouble> (m_equ[m_nequ-1]->GetTotPoints(),
0.0); //q0 is initialised
qBar0[i] = Array<OneD, NekDouble> (m_equ[m_nequ-1]->GetTotPoints(),
0.0);
m_equ[m_nequ - 1]->CopyFromPhysField(i, qBar0[i]);
}
///Definition of variables used in this algorithm
m_stepCounter = 0;
m_Check = 0;
m_Check_BaseFlow = 1;
m_NonConvergingStepsCounter = 0;
Diff_q_qBar = 1.0;
Diff_q1_q0 = 1.0;
cpuTime = 0.0;
elapsed = 0.0;
totalTime = 0.0;
FlowPartiallyConverged = false;
while (max(Diff_q_qBar, Diff_q1_q0) > TOL)
{
///Call the Navier-Stokes solver for one time step
m_equ[m_nequ - 1]->DoSolve();
for(int i = 0; i < NumVar_SFD; ++i)
{
///Copy the current flow field into q0
m_equ[m_nequ - 1]->CopyFromPhysField(i, q0[i]);
///Apply the linear operator to the outcome of the solver
ComputeSFD(i, q0, qBar0, q1, qBar1);
///Update qBar
qBar0[i] = qBar1[i];
///Copy the output of the SFD method into the current flow field
m_equ[m_nequ - 1]->CopyToPhysField(i, q1[i]);
}
if(m_infosteps && !((m_stepCounter+1)%m_infosteps))
{
ConvergenceHistory(qBar1, q0, Diff_q_qBar, Diff_q1_q0);
///Loop for the adaptive SFD method
if (m_EvolutionOperator == eAdaptiveSFD && FlowPartiallyConverged == false)
{
m_NonConvergingStepsCounter++;
if (Diff_q_qBar < AdaptiveTOL)
{
if (m_comm->GetRank() == 0)
{
cout << "\n\t The SFD method is converging: we compute "
<< "stability analysis using the 'partially "
<< "converged' steady state as base flow:\n" << endl;
}
m_equ[m_nequ - 1]->Checkpoint_BaseFlow(m_Check_BaseFlow);
m_Check_BaseFlow++;
A->GetAdvObject()->SetBaseFlow(q0);
DriverModifiedArnoldi::v_Execute(out);
if (m_comm->GetRank() == 0)
{
// Compute the update of the SFD parameters only on
// one processor
ComputeOptimization();
}
else
{
// On all the other processors, the parameters are set
// to 0
m_X = 0;
m_Delta = 0;
}
// The we give to all the processors the value of X and
// Delta of the first processor
m_comm->AllReduce(m_X, Nektar::LibUtilities::ReduceSum);
m_comm->AllReduce(m_Delta, Nektar::LibUtilities::ReduceSum);
FlowPartiallyConverged = true;
}
else if (m_NonConvergingStepsCounter * m_dt * m_infosteps
>= AdaptiveTime)
{
if (m_comm->GetRank() == 0)
{
cout << "\n\t We compute stability analysis using"
<< " the current flow field as base flow:\n"
<< endl;
}
m_equ[m_nequ - 1]->Checkpoint_BaseFlow(m_Check_BaseFlow);
m_Check_BaseFlow++;
A->GetAdvObject()->SetBaseFlow(q0);
DriverModifiedArnoldi::v_Execute(out);
if (m_comm->GetRank() == 0)
{
// Compute the update of the SFD parameters only on
// one processor
ComputeOptimization();
}
else
{
// On all the other processors, the parameters are set
// to 0
m_X = 0;
m_Delta = 0;
}
// The we give to all the processors the value of X and // Delta of the first processor
m_comm->AllReduce(m_X, Nektar::LibUtilities::ReduceSum);
m_comm->AllReduce(m_Delta, Nektar::LibUtilities::ReduceSum);
m_NonConvergingStepsCounter = 0;
}
}
}
if(m_checksteps && m_stepCounter&&(!((m_stepCounter+1)%m_checksteps)))
{
m_Check++;
m_equ[m_nequ - 1]->Checkpoint_Output(m_Check);
}
m_stepCounter++;
}
m_file.close();
///We save the final solution into a .fld file
m_equ[m_nequ - 1]->Output();
for(int j = 0; j < m_equ[m_nequ - 1]->GetNvariables(); ++j)
{
NekDouble vL2Error = m_equ[m_nequ - 1]->L2Error(j,false);
NekDouble vLinfError = m_equ[m_nequ - 1]->LinfError(j);
if (m_comm->GetRank() == 0)
{
out << "L 2 error (variable " << m_equ[m_nequ - 1]->GetVariable(j)
<< ") : " << vL2Error << endl;
out << "L inf error (variable " << m_equ[m_nequ - 1]->GetVariable(j)
<< ") : " << vLinfError << endl;
}
}
}
/**
* This routine defines the encapsulated SFD operator with first-order
* splitting and exact resolution of the second subproblem.
*(See http://scitation.aip.org/content/aip/journal/pof2/26/3/10.1063/1.4867482 for details)
*/
void DriverSteadyState::SetSFDOperator(const NekDouble X_input,
const NekDouble Delta_input)
{
NekDouble X = X_input*m_dt;
NekDouble Delta = Delta_input/m_dt;
NekDouble coeff = 1.0/(1.0 + X*Delta);
M11 = coeff*(1.0 + X*Delta*exp(-(X + 1.0/Delta)));
M12 = coeff*(X*Delta*(1.0 - exp(-(X + 1.0/Delta))));
M21 = coeff*(1.0 - exp(-(X + 1.0/Delta)));
M22 = coeff*(X*Delta + exp(-(X + 1.0/Delta)));
}
void DriverSteadyState::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)
{
q1[i] = Array<OneD, NekDouble> (m_equ[m_nequ - 1]->GetTotPoints(),0.0);
qBar1[i] = Array<OneD, NekDouble> (m_equ[m_nequ - 1]->GetTotPoints(),0.0);
///Encapsulated SFD method
Vmath::Svtvp(q1[i].num_elements(), M11, q0[i], 1, q1[i], 1, q1[i], 1 );
Vmath::Svtvp(q1[i].num_elements(), M12, qBar0[i], 1, q1[i], 1, q1[i], 1 );
Vmath::Svtvp(qBar1[i].num_elements(), M21, q0[i], 1, qBar1[i], 1,
qBar1[i], 1 ); Vmath::Svtvp(qBar1[i].num_elements(), M22, qBar0[i], 1, qBar1[i], 1,
qBar1[i], 1 );
}
void DriverSteadyState::ComputeOptimization()
{
NekDouble growthEV(0.0);
NekDouble frequencyEV(0.0);
// m_kdim is the dimension of Krylov subspace (defined in the xml file and
// used in DriverArnoldi.cpp)
ReadEVfile(m_kdim, growthEV, frequencyEV);
cout << "\n\tgrowthEV = " << growthEV << endl;
cout << "\tfrequencyEV = " << frequencyEV << endl;
complex<NekDouble> ApproxEV = polar(exp(growthEV), frequencyEV);
NekDouble X_new = m_X;
NekDouble Delta_new = m_Delta;
GradientDescentMethod(ApproxEV, X_new, Delta_new);
m_X = X_new;
m_Delta = Delta_new;
SetSFDOperator(m_X, m_Delta);
}
/**
* This routine implements a gradient descent method to find the parameters X
* end Delta which give the minimum eigenlavue of the SFD problem applied to
* the scalar case u(n+1) = \alpha*u(n).
*/
void DriverSteadyState::GradientDescentMethod(
const complex<NekDouble> &alpha,
NekDouble &X_output,
NekDouble &Delta_output)
{
cout << "\n\tWe enter the Gradient Descent Method [...]" << endl;
bool OptParmFound = false;
bool Descending = true;
NekDouble X_input = X_output;
NekDouble Delta_input = Delta_output;
NekDouble X_init = X_output;
NekDouble Delta_init = Delta_output;
int stepCounter(0);
NekDouble F0(0.0);
NekDouble F0x(0.0);
NekDouble F0y(0.0);
NekDouble F1(0.0);
NekDouble dx = 0.00000001;
NekDouble dirx(0.0);
NekDouble diry(0.0);
NekDouble s(0.0);
NekDouble CurrentMin = 1.0;
while (OptParmFound == false)
{
Descending = true;
EvalEV_ScalarSFD(X_input, Delta_input, alpha, F0);
EvalEV_ScalarSFD(X_input + dx, Delta_input, alpha, F0x);
EvalEV_ScalarSFD(X_input, Delta_input + dx, alpha, F0y);
dirx = (F0 - F0x);
diry = (F0 - F0y);
s = abs(0.000001/dirx);
X_output = X_input + s*dirx; Delta_output = Delta_input + s*diry;
F1 = F0;
while (Descending == true)
{
CurrentMin = F1;
X_input = X_output;
Delta_input = Delta_output;
EvalEV_ScalarSFD(X_output, Delta_output, alpha, F1);
if (F1 > CurrentMin)
{
Descending = false;
}
else
{
s = s + s*0.01;
X_output = X_input + s*dirx;
Delta_output = Delta_input + s*diry;
}
if(stepCounter > 9999999)
{
//We are stuck in this loop..
//Then we restart it with different initail conditions
Descending = false;
X_input = X_init;
Delta_init = Delta_init + Delta_init*0.1;
Delta_input = Delta_init;
stepCounter = 0;
}
stepCounter++;
}
if (abs(F0-F1) < dx)
{
cout << "\tThe Gradient Descent Method has converged!" << endl;
EvalEV_ScalarSFD(X_output, Delta_output, alpha, F1);
cout << "\n\tThe updated parameters are: X_tilde = " << X_output
<< " and Delta_tilde = " << Delta_output << endl;
OptParmFound = true;
}
}
}
/**
* This routine evaluates the maximum eigenvalue of the SFD system when applied
* to the 1D model u(n+1) = alpha*u(n)
*/
void DriverSteadyState::EvalEV_ScalarSFD(
const NekDouble &X_input,
const NekDouble &Delta_input,
const complex<NekDouble> &alpha,
NekDouble &MaxEV)
{
NekDouble A11 = ( 1.0 + X_input * Delta_input *
exp(-(X_input + 1.0/Delta_input)) )/(1.0 + X_input*Delta_input);
NekDouble A12 = ( X_input*Delta_input - X_input * Delta_input *
exp(-(X_input + 1.0/Delta_input)) )/(1.0 + X_input*Delta_input);
NekDouble A21 = ( 1.0 - 1.0 *
exp(-(X_input + 1.0/Delta_input)) )/(1 + X_input*Delta_input);
NekDouble A22 = ( X_input*Delta_input + 1.0 *
exp(-(X_input + 1.0/Delta_input)) )/(1.0 + X_input*Delta_input);
complex<NekDouble> B11 = alpha;
NekDouble B12 = 0.0;
NekDouble B21 = 0.0;
NekDouble B22 = 1.0;
complex<NekDouble> a = A11*B11 + A12*B21;
NekDouble b = A11*B12 + A12*B22;
complex<NekDouble> c = A21*B11 + A22*B21;
NekDouble d = A21*B12 + A22*B22;
complex<NekDouble> delt = sqrt((a-d)*(a-d) + 4.0*b*c);
complex<NekDouble> lambda_1 = 0.5*(a+d + delt);
complex<NekDouble> lambda_2 = 0.5*(a+d - delt);
NekDouble NORM_1 = abs(lambda_1);
NekDouble NORM_2 = abs(lambda_2);
MaxEV = max(NORM_1, NORM_2);
}
void DriverSteadyState::ReadEVfile(
int &KrylovSubspaceDim,
NekDouble &growthEV,
NekDouble &frequencyEV)
{
// This routine reads the .evl file written by the Arnoldi algorithm
// (written in September 2014)
std::string line;
int NumLinesInFile = 0;
std::string EVfileName = m_session->GetSessionName() + ".evl";
std::ifstream EVfile(EVfileName.c_str());
ASSERTL0(EVfile.good(), "Cannot open .evl file.");
if(EVfile)
{
// This block counts the total number of lines of the .evl file
// We keep going util we reach the end of the file
while(getline(EVfile, line))
{
++NumLinesInFile;
}
// It may happen that the Stability method that have produced the .elv
// file converges in less than m_kdim iterations. In this case,
// KrylovSubspaceDim has to be changed here
if(NumLinesInFile < KrylovSubspaceDim*2.0
+ KrylovSubspaceDim*(KrylovSubspaceDim+1.0)/2.0)
{
for(int i = 1; i <= KrylovSubspaceDim; ++i)
{
if(NumLinesInFile == i*2.0 + i*(i+1.0)/2.0)
{
KrylovSubspaceDim = i;
}
}
}
// go back to the beginning of the file
EVfile.clear();
EVfile.seekg(0, ios::beg);
// We now want to go to the line where the most unstable eigenlavue was
// written
for(int i = 0; i < (NumLinesInFile - KrylovSubspaceDim); ++i)
{
std::getline(EVfile, line);
cout << "Discard line: " << line << endl;
}
std::vector<std::string> tokens;
std::getline(EVfile, line);
boost::algorithm::split(tokens, line,
boost::is_any_of("\t "),boost::token_compress_on);
ASSERTL0(tokens.size() >= 5,
"Unexpected formatting of .evl file while reading line:\n"
+ line);
growthEV = boost::lexical_cast<NekDouble>(tokens[4]);
frequencyEV = boost::lexical_cast<NekDouble>(tokens[5]);
}
else
{
cout << "An error occurred when opening the .evl file" << endl;
}
EVfile.close();
}
/**
* This routine evaluates |q-qBar|_inf (and |q1-q0|_inf) and writes the values
* in "ConvergenceHistory.txt"
*/
void DriverSteadyState::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)
{
Array<OneD, NekDouble > NormDiff_q_qBar(NumVar_SFD, 1.0);
Array<OneD, NekDouble > NormDiff_q1_q0(NumVar_SFD, 1.0);
MaxNormDiff_q_qBar=0.0;
MaxNormDiff_q1_q0=0.0;
for(int i = 0; i < NumVar_SFD; ++i)
{
NormDiff_q_qBar[i] = m_equ[m_nequ - 1]->LinfError(i, qBar1[i]);
NormDiff_q1_q0[i] = m_equ[m_nequ - 1]->LinfError(i, q0[i]);
if (MaxNormDiff_q_qBar < NormDiff_q_qBar[i])
{
MaxNormDiff_q_qBar = NormDiff_q_qBar[i];
}
if (MaxNormDiff_q1_q0 < NormDiff_q1_q0[i])
{
MaxNormDiff_q1_q0 = NormDiff_q1_q0[i];
}
}
timer.Stop();
elapsed = timer.TimePerTest(1);
cpuTime += elapsed;
totalTime += elapsed;
if (m_comm->GetRank() == 0)
{
cout << "SFD - Step: " << left << m_stepCounter+1
<< ";\tTime: " << left << m_equ[m_nequ - 1]->GetFinalTime()
<< ";\tCPU time = " << left << cpuTime << " s"
<< ";\tTot time = " << left << totalTime << " s"
<< ";\tX = " << left << m_X
<< ";\tDelta = " << left << m_Delta
<< ";\t|q-qBar|inf = " << left << MaxNormDiff_q_qBar << endl;
std::ofstream m_file( "ConvergenceHistory.txt", std::ios_base::app);
m_file << m_stepCounter+1 << "\t"
<< m_equ[m_nequ - 1]->GetFinalTime() << "\t"
<< totalTime << "\t"
<< MaxNormDiff_q_qBar << "\t"
<< MaxNormDiff_q1_q0 << endl;
m_file.close();
}
cpuTime = 0.0;
timer.Start();}
void DriverSteadyState::PrintSummarySFD()
{
cout << "\n====================================="
"=================================="
<< endl;
cout << "Parameters for the SFD method:" << endl;
cout << "\tControl Coefficient: X = " << m_X << endl;
cout << "\tFilter Width: Delta = " << m_Delta << endl;
cout << "The simulation is stopped when |q-qBar|inf < " << TOL << endl;
if (m_EvolutionOperator == eAdaptiveSFD)
{
cout << "\nWe use the adaptive SFD method:" << endl;
cout << " The parameters are updated every " << AdaptiveTime
<< " time units;" << endl;
cout << " until |q-qBar|inf becomes smaller than " << AdaptiveTOL
<< endl;
}
cout << "====================================="
"==================================\n" << endl;
}
}
}