Nektar::VortexWaveInteraction Class Reference

#include <VortexWaveInteraction.h>

Collaboration diagram for Nektar::VortexWaveInteraction:

Public Member Functions
	VortexWaveInteraction (int argc, char *argv[])
	~VortexWaveInteraction (void)
void	ExecuteRoll (void)
void	ExecuteStreak (void)
void	ExecuteWave (void)
void	ExecuteLoop (bool CalcWaveForce=true)
void	SaveLoopDetails (string dir, int i)
void	CalcNonLinearWaveForce (void)
void	CalcL2ToLinfPressure (void)
void	SaveFile (string fileend, string dir, int n)
void	MoveFile (string fileend, string dir, int n)
void	CopyFile (string file1end, string file2end)
bool	CheckEigIsStationary (bool reset=false)
bool	CheckIfAtNeutralPoint (void)
void	UpdateAlpha (int n)
void	UpdateWaveForceMag (int n)
void	UpdateDAlphaDWaveForceMag (NekDouble alphainit)
void	AppendEvlToFile (std::string file, int n)
void	AppendEvlToFile (std::string file, NekDouble WaveForceMag)
int	GetIterStart ()
int	GetIterEnd ()
VWIIterationType	GetVWIIterationType (void)
int	GetNOuterIterations (void)
int	GetMaxOuterIterations (void)
NekDouble	GetAlpha (void)
NekDouble	GetAlphaStep (void)
NekDouble	GetWaveForceMag (void)
NekDouble	GetWaveForceMagStep (void)
NekDouble	GetDAlphaDWaveForceMag (void)
int	GetMaxWaveForceMagIter (void)
NekDouble	GetEigRelTol (void)
int	GetMinInnerIterations (void)
NekDouble	GetPrevAlpha (void)
void	SetAlpha (NekDouble alpha)
void	SetWaveForceMag (NekDouble mag)
void	SetEigRelTol (NekDouble tol)
void	SetAlphaStep (NekDouble step)
void	SetMinInnerIterations (int niter)
void	SetPrevAlpha (NekDouble alpha)
bool	IfIterInterface (void)
Array< OneD, int >	GetReflectionIndex (void)
void	FileRelaxation (int reg)

Private Attributes
int	m_iterStart
int	m_iterEnd
int	m_nOuterIterations
int	m_maxOuterIterations
int	m_minInnerIterations
int	m_maxWaveForceMagIter
bool	m_deltaFcnApprox
bool	m_useLinfPressureNorm
bool	m_moveMeshToCriticalLayer
NekDouble	m_deltaFcnDecay
Array< OneD, NekDouble >	m_waveForceMag
NekDouble	m_waveForceMagStep
NekDouble	m_rollForceScale
Array< OneD, NekDouble >	m_leading_real_evl
Array< OneD, NekDouble >	m_leading_imag_evl
	< Leading real eigenvalue
Array< OneD, NekDouble >	m_alpha
	< Leading imaginary eigenvalue
NekDouble	m_alphaStep
NekDouble	m_neutralPointTol
NekDouble	m_eigRelTol
NekDouble	m_vwiRelaxation
NekDouble	m_dAlphaDWaveForceMag
NekDouble	m_prevAlpha
bool	m_iterinterface
VWIIterationType	m_VWIIterationType
Array< OneD, MultiRegions::ExpListSharedPtr >	m_waveVelocities
MultiRegions::ExpListSharedPtr	m_wavePressure
Array< OneD, Array< OneD, NekDouble > >	m_vwiForcing
SolverUtils::ForcingProgrammaticSharedPtr	m_vwiForcingObj
Array< OneD, Array< OneD, NekDouble > >	m_bcsForcing
Array< OneD, MultiRegions::ExpListSharedPtr >	m_streakField
Array< OneD, MultiRegions::ExpListSharedPtr >	m_rollField
string	m_sessionName
LibUtilities::SessionReaderSharedPtr	m_sessionVWI
LibUtilities::SessionReaderSharedPtr	m_sessionRoll
EquationSystemSharedPtr	m_solverRoll
LibUtilities::SessionReaderSharedPtr	m_sessionStreak
LibUtilities::SessionReaderSharedPtr	m_sessionWave

Detailed Description

Definition at line 72 of file VortexWaveInteraction.h.

Constructor & Destructor Documentation

Nektar::VortexWaveInteraction::VortexWaveInteraction	(	int	argc,
		char *	argv[]
	)

Definition at line 45 of file VortexWaveInteraction.cpp.

References ASSERTL0, Nektar::LibUtilities::SessionReader::CreateInstance(), Nektar::LibUtilities::NekFactory< tKey, tBase, >::CreateInstance(), Nektar::eFixedAlpha, Nektar::eFixedAlphaWaveForcing, Nektar::eFixedWaveForcing, Nektar::eVWIIterationTypeSize, Nektar::SolverUtils::GetEquationSystemFactory(), m_alpha, m_alphaStep, m_dAlphaDWaveForceMag, m_deltaFcnApprox, m_deltaFcnDecay, m_eigRelTol, m_iterEnd, m_iterinterface, m_iterStart, m_leading_imag_evl, m_leading_real_evl, m_maxOuterIterations, m_maxWaveForceMagIter, m_minInnerIterations, m_moveMeshToCriticalLayer, m_neutralPointTol, m_nOuterIterations, m_rollForceScale, m_sessionName, m_sessionRoll, m_sessionStreak, m_sessionVWI, m_sessionWave, m_solverRoll, m_useLinfPressureNorm, m_vwiForcing, m_VWIIterationType, m_vwiRelaxation, m_waveForceMag, m_waveForceMagStep, UpdateAlpha(), and Nektar::VWIIterationTypeMap.

                                                                       :
        m_nOuterIterations(0)
    {

        int storesize;// number of previous iterations to store

        m_sessionName = argv[argc-1];
        string meshfile = m_sessionName + ".xml";

        LibUtilities::SessionReaderSharedPtr session; 

        std::string VWICondFile(argv[argc-1]); 
        VWICondFile += "_VWI.xml"; 
        std::vector<std::string> VWIFilenames;
        VWIFilenames.push_back(meshfile);   
        VWIFilenames.push_back(VWICondFile);

        // Create Incompressible NavierStokesSolver session reader.
        m_sessionVWI = LibUtilities::SessionReader::CreateInstance(argc, argv, 
                                                                   VWIFilenames);
        
        m_sessionVWI->LoadParameter("AlphaStep",              m_alphaStep,0.05);
        m_sessionVWI->LoadParameter("OuterIterationStoreSize",storesize,10);

        //set a low tol for interfaceVWI
    m_sessionVWI->LoadParameter("EigenvalueRelativeTol",  m_eigRelTol,1e-3);
       
        
        m_sessionVWI->LoadParameter("NeutralPointTolerance",  m_neutralPointTol,1e-4);
        m_sessionVWI->LoadParameter("MaxOuterIterations",     m_maxOuterIterations,100);
        
        m_sessionVWI->LoadParameter("StartIteration",m_iterStart, 0);
        m_sessionVWI->LoadParameter("EndIteration",  m_iterEnd, 0);

        m_sessionVWI->LoadParameter("WaveForceMagStep",   m_waveForceMagStep,0.01);
        m_sessionVWI->LoadParameter("DAlphaDWaveForceMag", m_dAlphaDWaveForceMag,0.0);
        m_sessionVWI->LoadParameter("MaxWaveForceMagIter",m_maxWaveForceMagIter,1);
        m_sessionVWI->LoadParameter("RollForceScale",     m_rollForceScale,1.0);
        
        if(m_sessionVWI->DefinesSolverInfo("DeltaFcnApprox"))
        {
            m_deltaFcnApprox = true;
            m_sessionVWI->LoadParameter("DeltaFcnDecay", m_deltaFcnDecay,1.0/500);
        }
        else
        {
            m_deltaFcnApprox = false;
            m_deltaFcnDecay = 0.0;
        }

        if(m_sessionVWI->DefinesSolverInfo("LinfPressureNorm"))
        {
            m_useLinfPressureNorm  = true;
        }
        else
        {
            m_useLinfPressureNorm  = false;
        }

        if( m_sessionVWI->DefinesSolverInfo("INTERFACE")  )
        {
            m_iterinterface = true;
        }
        else
        {
            m_iterinterface = false;
        }

        if(m_sessionVWI->DefinesSolverInfo("MoveMeshToCriticalLayer"))
        {
            m_moveMeshToCriticalLayer = true;
        }
        else
        {
            m_moveMeshToCriticalLayer = false;
        }

        m_alpha = Array<OneD, NekDouble> (storesize);
        m_alpha[0]         = m_sessionVWI->GetParameter("Alpha");
        m_waveForceMag     = Array<OneD, NekDouble> (storesize);
        m_waveForceMag[0]  = m_sessionVWI->GetParameter("WaveForceMag");
        
        m_leading_real_evl = Array<OneD, NekDouble> (storesize);
        m_leading_imag_evl = Array<OneD, NekDouble> (storesize);
        
        if(m_sessionVWI->DefinesParameter("Relaxation"))
        {
            m_vwiRelaxation = m_sessionVWI->GetParameter("Relaxation");
            // fix minimum number of iterations to be number of
            // iterations required to make contribution of innitial
            // forcing to 0.1
            m_minInnerIterations = (int) (log(0.1)/log(m_vwiRelaxation)); 
        }
        else
        {
            m_vwiRelaxation = 0.0;
            m_minInnerIterations = 1;
        }
        
        // Initialise NS Roll solver 
        std::string IncCondFile(argv[argc-1]); 
        IncCondFile += "_IncNSCond.xml"; 
        std::vector<std::string> IncNSFilenames;
        IncNSFilenames.push_back(meshfile);   
        IncNSFilenames.push_back(IncCondFile);

        // Create Incompressible NavierStokesSolver session reader.
        m_sessionRoll = LibUtilities::SessionReader::CreateInstance(argc, argv, IncNSFilenames, 
                                                                    m_sessionVWI->GetComm());
        std::string vEquation = m_sessionRoll->GetSolverInfo("SolverType");
        m_solverRoll = GetEquationSystemFactory().CreateInstance(vEquation, m_sessionRoll);
        m_solverRoll->PrintSummary(cout);


        if(m_sessionRoll->DefinesSolverInfo("INTERFACE"))
    {
            m_sessionVWI->LoadParameter("EigenvalueRelativeTol",  m_eigRelTol,1e-2);
    }

        int ncoeffs = m_solverRoll->UpdateFields()[0]->GetNcoeffs();
        m_vwiForcing = Array<OneD, Array<OneD, NekDouble> > (4);
        m_vwiForcing[0] = Array<OneD, NekDouble> (4*ncoeffs);
        for(int i = 1; i < 4; ++i)
        {
            m_vwiForcing[i] = m_vwiForcing[i-1] + ncoeffs;
        }

        // Fill forcing into m_vwiForcing
        // Has forcing even been initialised yet?
//        Vmath::Vcopy(ncoeffs,m_solverRoll->UpdateForces()[0]->GetCoeffs(),1,m_vwiForcing[0],1);
//        Vmath::Vcopy(ncoeffs,m_solverRoll->UpdateForces()[1]->GetCoeffs(),1,m_vwiForcing[1],1);
        

        // Create AdvDiff Streak solver 
        std::string AdvDiffCondFile(argv[argc-1]); 
        AdvDiffCondFile += "_AdvDiffCond.xml"; 
        std::vector<std::string> AdvDiffFilenames;
        AdvDiffFilenames.push_back(meshfile);   
        AdvDiffFilenames.push_back(AdvDiffCondFile);
        
        // Create AdvDiffusion session reader.
        m_sessionStreak = LibUtilities::SessionReader::CreateInstance(argc, argv, AdvDiffFilenames, m_sessionVWI->GetComm());

        // Initialise LinNS solver 
        std::string LinNSCondFile(argv[argc-1]); 
        LinNSCondFile += "_LinNSCond.xml"; 
        std::vector<std::string> LinNSFilenames;
        LinNSFilenames.push_back(meshfile);   
        LinNSFilenames.push_back(LinNSCondFile);
        
        // Create Linearised NS stability session reader.
        m_sessionWave = LibUtilities::SessionReader::CreateInstance(argc, argv, LinNSFilenames, m_sessionVWI->GetComm());

        // Set the initial beta value in stability to be equal to VWI file
        std::string LZstr("LZ");
        NekDouble LZ = 2*M_PI/m_alpha[0];
        cout << "Setting LZ in Linearised solver to " << LZ << endl;
        m_sessionWave->SetParameter(LZstr,LZ);

        // Check for iteration type 
        if(m_sessionVWI->DefinesSolverInfo("VWIIterationType"))
        {
            std::string IterationTypeStr = m_sessionVWI->GetSolverInfo("VWIIterationType");
            for(int i = 0; i < (int) eVWIIterationTypeSize; ++i)
            {
                if(m_solverRoll->NoCaseStringCompare(VWIIterationTypeMap[i],IterationTypeStr) == 0 )
                {
                    m_VWIIterationType = (VWIIterationType)i; 
                    break;
                }
            }
        }
        else
        {
            m_VWIIterationType = eFixedAlphaWaveForcing;
        }



        // Check for restart
        bool restart;
        m_sessionVWI->MatchSolverInfo("RestartIteration","True",restart,false);
        if(restart)
        {
            switch(m_VWIIterationType)
            {
            case  eFixedAlpha:
                break;
            case  eFixedWaveForcing:
                {
                    FILE *fp;
                    // Check for OuterIter.his file to read
                    if((fp = fopen("OuterIter.his","r")))
                    {
                        char buf[BUFSIZ];
                        std::vector<NekDouble> Alpha, Growth, Phase;
                        NekDouble alpha,growth,phase;
                        while(fgets(buf,BUFSIZ,fp))
                        {
                            sscanf(buf,"%*d:%lf%lf%lf",&alpha,&growth,&phase);
                            Alpha.push_back(alpha);
                            Growth.push_back(growth);
                            Phase.push_back(phase);
                        }

                        m_nOuterIterations = Alpha.size();
                        
                        int nvals = std::min(m_nOuterIterations,(int)m_alpha.num_elements());
                        
                        for(int i = 0; i < nvals; ++i)
                        {
                            m_alpha[nvals-1-i]            = Alpha[m_nOuterIterations-nvals+i];
                            m_leading_real_evl[nvals-1-i] = Growth[m_nOuterIterations-nvals+i];
                            m_leading_imag_evl[nvals-1-i] = Phase [m_nOuterIterations-nvals+i];
                        }

                        UpdateAlpha(m_nOuterIterations++);
                    }
                    else
                    {
                        cout << " No File OuterIter.his to restart from" << endl;
                    }
                }
                break;
            case eFixedAlphaWaveForcing:
                {
                    string nstr =  boost::lexical_cast<std::string>(m_iterStart);
                    cout << "Restarting from iteration " << m_iterStart << endl;
                    std::string rstfile = "cp -f Save/" + m_sessionName + ".rst." + nstr + " " + m_sessionName + ".rst"; 
                    cout << "      " << rstfile << endl;
                    if(system(rstfile.c_str()))
                    {
                        ASSERTL0(false,rstfile.c_str());
                    }
                    std::string vwifile = "cp -f Save/" + m_sessionName + ".vwi." + nstr + " " + m_sessionName + ".vwi"; 
                    cout << "      " << vwifile << endl;
                    if(system(vwifile.c_str()))
                    {
                        ASSERTL0(false,vwifile.c_str());
                    }
                }
                break;
            default:
                ASSERTL0(false,"Unknown VWIITerationType in restart");
            }
        }

        // Check for ConveredSoln to update DAlphaDWaveForce
        {
            FILE *fp;
            if((fp = fopen("ConvergedSolns","r")))
            {
                char buf[BUFSIZ];
                std::vector<NekDouble> WaveForce, Alpha;
                NekDouble waveforce,alpha;
                while(fgets(buf,BUFSIZ,fp))
                {
                    sscanf(buf,"%*d:%lf%lf",&waveforce,&alpha);
                    WaveForce.push_back(waveforce);
                    Alpha.push_back(alpha);
                }

                if(Alpha.size() > 1)
                {
                    int min_i = 0;
                    NekDouble min_alph = fabs(m_alpha[0]-Alpha[min_i]);
                    // find nearest point 
                    for(int i = 1; i < Alpha.size(); ++i)
                    {
                        if(fabs(m_alpha[0]-Alpha[min_i]) < min_alph)
                        {
                            min_i = i;
                            min_alph = fabs(m_alpha[0]-Alpha[min_i]);
                        }
                    }

                    // find next nearest point 
                    int min_j = (min_i == 0)? 1:0;
                    min_alph = fabs(m_alpha[0]-Alpha[min_j]);
                    for(int i = 0; i < Alpha.size(); ++i)
                    {
                        if(i != min_i)
                        {
                            if(fabs(m_alpha[0]-Alpha[min_j]) < min_alph)
                            {
                                min_j = i;
                                min_alph = fabs(m_alpha[0]-Alpha[min_j]);
                            }
                        }
                    }
                    
                    if(fabs(Alpha[min_i] - Alpha[min_j]) > 1e-4)
                    {
                        m_dAlphaDWaveForceMag = (Alpha[min_i]-Alpha[min_j])/(WaveForce[min_i]-WaveForce[min_j]);
                    }
                }
            }
        }