#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 More...

Array< OneD, NekDouble >	m_alpha
	< Leading imaginary eigenvalue More...

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]);
                     }
                 }
             }
         }
     }

Nektar::VortexWaveInteraction::~VortexWaveInteraction ( void )

Definition at line 346 of file VortexWaveInteraction.cpp.

References m_sessionVWI.

     {
         m_sessionVWI->Finalise();
     }

Member Function Documentation

void Nektar::VortexWaveInteraction::AppendEvlToFile	(	std::string	file,
		int	n
	)

Definition at line 1060 of file VortexWaveInteraction.cpp.

References m_alpha, m_leading_imag_evl, and m_leading_real_evl.

Referenced by DoFixedForcingIteration(), and main().

      {
          FILE *fp;
          fp = fopen(file.c_str(),"a");
          fprintf(fp, "%d: %lf %16.12le  %16.12le\n",n, m_alpha[0], m_leading_real_evl[0],m_leading_imag_evl[0]);
          fclose(fp);
      }

void Nektar::VortexWaveInteraction::AppendEvlToFile	(	std::string	file,
		NekDouble	WaveForceMag
	)

Definition at line 1068 of file VortexWaveInteraction.cpp.

References m_alpha, m_leading_imag_evl, and m_leading_real_evl.

      {
          FILE *fp;
          fp = fopen(file.c_str(),"a");
          fprintf(fp, "%lf %lf %16.12le  %16.12le\n",WaveForceMag, m_alpha[0], m_leading_real_evl[0],m_leading_imag_evl[0]);
          fclose(fp);
      }

void Nektar::VortexWaveInteraction::CalcL2ToLinfPressure ( void )

Definition at line 963 of file VortexWaveInteraction.cpp.

References ExecuteWave(), Vmath::Fill(), m_wavePressure, m_waveVelocities, npts, Vmath::Vmax(), Vmath::Vmul(), Vmath::Vsqrt(), and Vmath::Vvtvp().

Referenced by main().

     {
 
         ExecuteWave();
 
         m_wavePressure->GetPlane(0)->BwdTrans(m_wavePressure->GetPlane(0)->GetCoeffs(),
                                               m_wavePressure->GetPlane(0)->UpdatePhys());
         m_wavePressure->GetPlane(1)->BwdTrans(m_wavePressure->GetPlane(1)->GetCoeffs(),
                                               m_wavePressure->GetPlane(1)->UpdatePhys());
 
         int npts    = m_waveVelocities[0]->GetPlane(0)->GetNpoints();
         NekDouble Linf;
         Array<OneD, NekDouble> val(2*npts,0.0);
         
         Vmath::Vmul(npts,m_wavePressure->GetPlane(0)->UpdatePhys(),1,m_wavePressure->GetPlane(0)->UpdatePhys(),1,val,1);
         Vmath::Vvtvp(npts,m_wavePressure->GetPlane(1)->UpdatePhys(),1,m_wavePressure->GetPlane(1)->UpdatePhys(),1,val,1,val,1);
         cout << "int P^2 " << m_wavePressure->GetPlane(0)->PhysIntegral(val) << endl;
         Vmath::Vsqrt(npts,val,1,val,1);
         
                
         Linf = Vmath::Vmax(npts,val,1);
         cout << "Linf: " << Linf << endl;
 
         NekDouble l2,norm;
         l2    = m_wavePressure->GetPlane(0)->L2(m_wavePressure->GetPlane(0)->GetPhys());
         norm  = l2*l2;
         l2    = m_wavePressure->GetPlane(1)->L2(m_wavePressure->GetPlane(1)->GetPhys());
         norm += l2*l2;
 
 
         Vmath::Fill(npts,1.0,val,1);
         NekDouble area = m_waveVelocities[0]->GetPlane(0)->PhysIntegral(val);
 
         l2 = sqrt(norm/area);
 
         cout << "L2:   " << l2 << endl;
         
         cout << "Ratio Linf/L2: "<< Linf/l2 << endl;
     }

void Nektar::VortexWaveInteraction::CalcNonLinearWaveForce ( void )

Definition at line 558 of file VortexWaveInteraction.cpp.

References ASSERTL0, Nektar::SpatialDomains::eNeumann, Vmath::Fill(), GetReflectionIndex(), m_alpha, m_deltaFcnApprox, m_deltaFcnDecay, m_sessionName, m_sessionRoll, m_sessionVWI, m_solverRoll, m_streakField, m_useLinfPressureNorm, m_vwiForcing, m_vwiRelaxation, m_waveForceMag, m_wavePressure, m_waveVelocities, Vmath::Neg(), npts, Vmath::Sadd(), Vmath::Smul(), Vmath::Svtvp(), Vmath::Vadd(), Vmath::Vcopy(), Vmath::Vmax(), Vmath::Vmul(), Vmath::Vsqrt(), Vmath::Vsub(), and Vmath::Vvtvp().

Referenced by DoFixedForcingIteration(), ExecuteLoop(), and main().

     {
         if(m_sessionRoll->DefinesSolverInfo("INTERFACE") )
         {
             static int cnt=0;
             string wavefile  = m_sessionName +".fld"; 
             string movedmesh = m_sessionName + "_advPost_moved.xml";
             string filestreak   = m_sessionName + "_streak.fld";
             char c[16]="";
             sprintf(c,"%d",cnt);               
             char c_alpha[16]="";
             sprintf(c_alpha,"%f",m_alpha[0]);    
             string syscall;
             if( m_sessionVWI->GetSolverInfo("INTERFACE")=="phase" )
             {
                 string filePost = m_sessionName + "_advPost.xml";
                 syscall = "../../utilities/PostProcessing/Extras/FldCalcBCs  "
                     + filePost +"     "+
                     "meshhalf_pos_Spen_stability_moved.fld  meshhalf_pos_Spen_advPost_moved.fld "
                     +c_alpha +"  > data_alpha0";
                 cout<<syscall.c_str()<<endl;
                 if(system(syscall.c_str()))
                 {
                     ASSERTL0(false,syscall.c_str());
                 }
              
                 syscall = "cp -f meshhalf_pos_Spen_stability_moved_u_5.bc  "+m_sessionName+"_u_5.bc";  
                 cout<<syscall.c_str()<<endl;
                 if(system(syscall.c_str()))
                 {
                     ASSERTL0(false,syscall.c_str());
                 }
                 syscall = "cp -f meshhalf_pos_Spen_stability_moved_v_5.bc  "+m_sessionName+"_v_5.bc";  
                 cout<<syscall.c_str()<<endl;
                 if(system(syscall.c_str()))
                 {
                     ASSERTL0(false,syscall.c_str());
                 }
             }
             else
             {
                 syscall =  "../../utilities/PostProcessing/Extras/FldCalcBCs  "
                     + movedmesh + "  " + wavefile + "  " + filestreak + "   "+c_alpha +"  >  datasub_"+c;
                 cout<<syscall.c_str()<<endl;
                 if(system(syscall.c_str()))
                 {
                     ASSERTL0(false,syscall.c_str());
                 }
             }
 
              
              
             //relaxation for different alpha values? does it make sense?
 
             //save the wave
             string wave_subalp = m_sessionName + "_wave_subalp_"+c+".fld";
             syscall = "cp -f " + wavefile + "  " + wave_subalp;
             cout<<syscall.c_str()<<endl;
             if(system(syscall.c_str()))
             {
                 ASSERTL0(false,syscall.c_str());
             }
             //FileRelaxation(3);
             cnt++;
      
         }
         else
         {
             int npts    = m_waveVelocities[0]->GetPlane(0)->GetNpoints();
             int ncoeffs = m_waveVelocities[0]->GetPlane(0)->GetNcoeffs();
             Array<OneD, NekDouble> val(npts), der1(2*npts);
             Array<OneD, NekDouble> der2 = der1 + npts; 
             Array<OneD, NekDouble> streak;
             static int projectfield = -1;
 
             if(m_deltaFcnApprox)
             {
                 streak = Array<OneD, NekDouble> (npts);
                 m_streakField[0]->BwdTrans(m_streakField[0]->GetCoeffs(), streak);
             }
  
             // Set project field to be first field that has a Neumann
             // boundary since this not impose any condition on the vertical boundaries
             // Othersise set to zero. 
             if(projectfield == -1)
             {
                 Array<OneD, const SpatialDomains::BoundaryConditionShPtr > BndConds;
                  
                 for(int i = 0; i < m_waveVelocities.num_elements(); ++i)
                 {
                     BndConds = m_waveVelocities[i]->GetBndConditions();
                     for(int j = 0; j < BndConds.num_elements(); ++j)
                     {
                         if(BndConds[j]->GetBoundaryConditionType() == SpatialDomains::eNeumann)
                         {
                             projectfield = i;
                             break;
                         }
                     }
                     if(projectfield != -1)
                     {
                         break;
                     }
                 }
                 if(projectfield == -1)
                 {
                     cout << "using first field to project non-linear forcing which imposes a Dirichlet condition" << endl;
                     projectfield = 0;
                 }
             }
         
             // Shift m_vwiForcing in case of relaxation 
             Vmath::Vcopy(ncoeffs,m_vwiForcing[0],1,m_vwiForcing[2],1);
             Vmath::Vcopy(ncoeffs,m_vwiForcing[1],1,m_vwiForcing[3],1);
         
             // determine inverse of area normalised field. 
             m_wavePressure->GetPlane(0)->BwdTrans(m_wavePressure->GetPlane(0)->GetCoeffs(), m_wavePressure->GetPlane(0)->UpdatePhys());
             m_wavePressure->GetPlane(1)->BwdTrans(m_wavePressure->GetPlane(1)->GetCoeffs(), m_wavePressure->GetPlane(1)->UpdatePhys());
             NekDouble invnorm;
 
             if(m_useLinfPressureNorm)
             {
                 Vmath::Vmul(npts,m_wavePressure->GetPlane(0)->UpdatePhys(),1,m_wavePressure->GetPlane(0)->UpdatePhys(),1,der1,1);
                 Vmath::Vvtvp(npts,m_wavePressure->GetPlane(1)->UpdatePhys(),1,m_wavePressure->GetPlane(1)->UpdatePhys(),1,der1,1,der1,1);
                 Vmath::Vsqrt(npts,der1,1,der1,1);
                 
                 NekDouble Linf = Vmath::Vmax(npts,der1,1);
                 
                 invnorm = 1.0/Linf;
             }
             else
             {
                 // Determine normalisation of pressure so that |P|/A = 1
                 NekDouble l2;
                 l2    = m_wavePressure->GetPlane(0)->L2(m_wavePressure->GetPlane(0)->GetPhys());
                 invnorm  = l2*l2;
                 l2    = m_wavePressure->GetPlane(1)->L2(m_wavePressure->GetPlane(1)->GetPhys());
                 invnorm += l2*l2;
                 Vmath::Fill(2*npts,1.0,der1,1);
                 NekDouble area = m_waveVelocities[0]->GetPlane(0)->PhysIntegral(der1);
                 cout << "Area: " << area << endl;
                 invnorm = sqrt(area/invnorm);
             }
         
             // Get hold of arrays. 
             m_waveVelocities[0]->GetPlane(0)->BwdTrans(m_waveVelocities[0]->GetPlane(0)->GetCoeffs(),m_waveVelocities[0]->GetPlane(0)->UpdatePhys());
             Array<OneD, NekDouble> u_real = m_waveVelocities[0]->GetPlane(0)->UpdatePhys();
             Vmath::Smul(npts,invnorm,u_real,1,u_real,1);
             m_waveVelocities[0]->GetPlane(1)->BwdTrans(m_waveVelocities[0]->GetPlane(1)->GetCoeffs(),m_waveVelocities[0]->GetPlane(1)->UpdatePhys());
             Array<OneD, NekDouble> u_imag = m_waveVelocities[0]->GetPlane(1)->UpdatePhys();
             Vmath::Smul(npts,invnorm,u_imag,1,u_imag,1);
             m_waveVelocities[1]->GetPlane(0)->BwdTrans(m_waveVelocities[1]->GetPlane(0)->GetCoeffs(),m_waveVelocities[1]->GetPlane(0)->UpdatePhys());
             Array<OneD, NekDouble> v_real = m_waveVelocities[1]->GetPlane(0)->UpdatePhys(); 
             Vmath::Smul(npts,invnorm,v_real,1,v_real,1);
             m_waveVelocities[1]->GetPlane(1)->BwdTrans(m_waveVelocities[1]->GetPlane(1)->GetCoeffs(),m_waveVelocities[1]->GetPlane(1)->UpdatePhys());
             Array<OneD, NekDouble> v_imag = m_waveVelocities[1]->GetPlane(1)->UpdatePhys();
             Vmath::Smul(npts,invnorm,v_imag,1,v_imag,1);
 
             // normalise wave for output
             Vmath::Smul(2*ncoeffs,invnorm,m_waveVelocities[0]->UpdateCoeffs(),1,m_waveVelocities[0]->UpdateCoeffs(),1);
             Vmath::Smul(2*ncoeffs,invnorm,m_waveVelocities[1]->UpdateCoeffs(),1,m_waveVelocities[1]->UpdateCoeffs(),1);
             Vmath::Smul(2*ncoeffs,invnorm,m_waveVelocities[2]->UpdateCoeffs(),1,m_waveVelocities[2]->UpdateCoeffs(),1);
 
             // dump field
             {
                 std::vector<std::string> variables(3);
                 std::vector<Array<OneD, NekDouble> > outfield(3);
                 variables[0] = "u_w"; 
                 variables[1] = "v_w"; 
                 variables[2] = "w_w"; 
                 outfield[0] = m_waveVelocities[0]->UpdateCoeffs();
                 outfield[1] = m_waveVelocities[1]->UpdateCoeffs();
                 outfield[2] = m_waveVelocities[2]->UpdateCoeffs();
                 std::string outname = m_sessionName  + "_wave.fld";
                 m_solverRoll->WriteFld(outname, m_waveVelocities[0], outfield, variables);
             }
 
 #if 1
             int ncoeffs_p = m_wavePressure->GetPlane(0)->GetNcoeffs();
             Vmath::Smul(ncoeffs_p,invnorm,m_wavePressure->GetPlane(0)->UpdateCoeffs(),1,m_wavePressure->GetPlane(0)->UpdateCoeffs(),1);
             Vmath::Smul(ncoeffs_p,invnorm,m_wavePressure->GetPlane(1)->UpdateCoeffs(),1,m_wavePressure->GetPlane(1)->UpdateCoeffs(),1);
 #else
             m_wavePressure->GetPlane(0)->BwdTrans(m_wavePressure->GetPlane(0)->GetCoeffs(),m_wavePressure->GetPlane(0)->UpdatePhys());
             Vmath::Smul(npts,invnorm,m_wavePressure->GetPlane(0)->UpdatePhys(),1,m_wavePressure->GetPlane(0)->UpdatePhys(),1);
             m_wavePressure->GetPlane(0)->FwdTrans(m_wavePressure->GetPlane(0)->UpdatePhys(),m_wavePressure->GetPlane(0)->UpdateCoeffs());
             
             m_wavePressure->GetPlane(1)->BwdTrans(m_wavePressure->GetPlane(1)->GetCoeffs(),m_wavePressure->GetPlane(1)->UpdatePhys());
             Vmath::Smul(npts,invnorm,m_wavePressure->GetPlane(1)->UpdatePhys(),1,m_wavePressure->GetPlane(1)->UpdatePhys(),1);
             m_wavePressure->GetPlane(1)->FwdTrans(m_wavePressure->GetPlane(1)->UpdatePhys(),m_wavePressure->GetPlane(1)->UpdateCoeffs());
 #endif
        
             // Calculate non-linear terms for x and y directions
             // d/dx(u u* + u* u)
             Vmath::Vmul (npts,u_real,1,u_real,1,val,1);
             Vmath::Vvtvp(npts,u_imag,1,u_imag,1,val,1,val,1);
             Vmath::Smul (npts,2.0,val,1,val,1);
             m_waveVelocities[0]->GetPlane(0)->PhysDeriv(0,val,der1);
         
             // d/dy(v u* + v* u)
             Vmath::Vmul (npts,u_real,1,v_real,1,val,1);
             Vmath::Vvtvp(npts,u_imag,1,v_imag,1,val,1,val,1);
             Vmath::Smul (npts,2.0,val,1,val,1);
             m_waveVelocities[0]->GetPlane(0)->PhysDeriv(1,val,der2);
             
             Vmath::Vadd(npts,der1,1,der2,1,der1,1);
 #if 1
             m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(der1,m_vwiForcing[0]);
 #else
             m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(der1,m_vwiForcing[0]);
 #endif
             Vmath::Smul(ncoeffs,-m_waveForceMag[0],m_vwiForcing[0],1,m_vwiForcing[0],1);
 
             // d/dx(u v* + u* v)
             m_waveVelocities[0]->GetPlane(0)->PhysDeriv(0,val,der1);
         
             // d/dy(v v* + v* v)
             Vmath::Vmul(npts,v_real,1,v_real,1,val,1);
             Vmath::Vvtvp(npts,v_imag,1,v_imag,1,val,1,val,1);
             Vmath::Smul (npts,2.0,val,1,val,1);
             m_waveVelocities[0]->GetPlane(0)->PhysDeriv(1,val,der2);
         
             Vmath::Vadd(npts,der1,1,der2,1,der1,1);
 
 #if 1
             m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(der1,m_vwiForcing[1]);
 #else
             m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(der1,m_vwiForcing[1]);
 #endif
 
             Vmath::Smul(ncoeffs,-m_waveForceMag[0],m_vwiForcing[1],1,m_vwiForcing[1],1);
         
             //by default the symmetrization is on
             bool symm=true;
             m_sessionVWI->MatchSolverInfo("Symmetrization","True",symm,true);
 #if 0
             if(symm== true )
             {
                 
                 // Symmetrise forcing
                 //-> Get coordinates 
                 Array<OneD, NekDouble> coord(2);
                 Array<OneD, NekDouble> coord_x(npts);
                 Array<OneD, NekDouble> coord_y(npts);
                 
                 //-> Impose symmetry (x -> -x + Lx/2, y-> -y) on coordinates
                 m_waveVelocities[0]->GetPlane(0)->GetCoords(coord_x,coord_y);
                 NekDouble xmax = Vmath::Vmax(npts,coord_x,1);
                 Vmath::Neg(npts,coord_x,1);
                 Vmath::Sadd(npts,xmax,coord_x,1,coord_x,1);
                 Vmath::Neg(npts,coord_y,1);
                 
                 int i, physoffset;
                 
                 //-> Obtain list of expansion element ids for each point. 
                 Array<OneD, int> Eid(npts);
                 // This search may not be necessary every iteration
                 for(i = 0; i < npts; ++i)
                 {
                     coord[0] = coord_x[i];
                     coord[1] = coord_y[i];
                     
                     // Note this will not quite be symmetric. 
                     Eid[i] = m_waveVelocities[0]->GetPlane(0)->GetExpIndex(coord,1e-6);
                 }
                 
                 // Interpolate field 0 
                 m_waveVelocities[0]->GetPlane(0)->BwdTrans_IterPerExp(m_vwiForcing[0],der1);
                 for(i = 0; i < npts; ++i)
                 {
                     physoffset = m_waveVelocities[0]->GetPlane(0)->GetPhys_Offset(Eid[i]);
                     coord[0] = coord_x[i];
                     coord[1] = coord_y[i];
                     der2 [i] = m_waveVelocities[0]->GetPlane(0)->GetExp(Eid[i])->PhysEvaluate(coord, der1 + physoffset);
                 }
                 //-> Average field 0 
                 Vmath::Vsub(npts,der1,1,der2,1,der2,1);
                 Vmath::Smul(npts,0.5,der2,1,der2,1);
 #if 1
                 m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(der2,m_vwiForcing[0]);
 #else
                 m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(der2, m_vwiForcing[0]);
 #endif
                 
                 //-> Interpoloate field 1
                 m_waveVelocities[0]->GetPlane(0)->BwdTrans_IterPerExp(m_vwiForcing[1],der1);
                 for(i = 0; i < npts; ++i)
                 {
                     physoffset = m_waveVelocities[0]->GetPlane(0)->GetPhys_Offset(Eid[i]);
                     coord[0] = coord_x[i];
                     coord[1] = coord_y[i];
                     der2[i]  = m_waveVelocities[0]->GetPlane(0)->GetExp(Eid[i])->PhysEvaluate(coord,                                                                         der1 + physoffset);
                 }
                 
                 //-> Average field 1
                 Vmath::Vsub(npts,der1,1,der2,1,der2,1);
                 Vmath::Smul(npts,0.5,der2,1,der2,1);
 #if 1
                 m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(der2,m_vwiForcing[1]);
 #else
                 m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(der2, m_vwiForceing[1]);
 #endif
             }
 #else
             int i;
             if(symm== true )
             {
                 cout<<"symmetrization is active"<<endl;              
                 static Array<OneD, int> index = GetReflectionIndex();
                 
                 m_waveVelocities[0]->GetPlane(0)->BwdTrans_IterPerExp(m_vwiForcing[0],der1);
                 for(i = 0; i < npts; ++i)
                 {
                     if(index[i] != -1)
                     {
                         val[i] = 0.5*(der1[i] - der1[index[i]]);
                     }
                 }
 #if 1 
                 m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(val,m_vwiForcing[0]);
 #else
                 m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(val, m_vwiForcing[0]);
 #endif
 
                 m_waveVelocities[0]->GetPlane(0)->BwdTrans_IterPerExp(m_vwiForcing[1],der2);
                 for(i = 0; i < npts; ++i)
                 {
                     if(index[i] != -1)
                     {
                         val[i] = 0.5*(der2[i] - der2[index[i]]);
                     }
                 }        
 #if 1 
                 m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(val,m_vwiForcing[1]);
 #else
                 m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(val, m_vwiForcing[1]);
 #endif
             }
 
 
             Vmath::Vmul(npts,der1,1,der1,1,val,1);
             Vmath::Vvtvp(npts,der2,1,der2,1,val,1,val,1);
             Vmath::Vsqrt(npts,val,1,val,1);
             cout << "F_Linf: " <<  Vmath::Vmax(npts,val,1) << endl;
 
 #endif
 
             if(m_vwiRelaxation)
             {
                  Vmath::Smul(ncoeffs,1.0-m_vwiRelaxation,
                         m_vwiForcing[0],1,m_vwiForcing[0],1);
                  Vmath::Svtvp(ncoeffs,m_vwiRelaxation,m_vwiForcing[2],1,
                          m_vwiForcing[0],1,m_vwiForcing[0],1);
             
                  Vmath::Smul(ncoeffs,1.0-m_vwiRelaxation,
                         m_vwiForcing[1],1,m_vwiForcing[1],1);
                  Vmath::Svtvp(ncoeffs,m_vwiRelaxation,m_vwiForcing[3],1,
                          m_vwiForcing[1],1,m_vwiForcing[1],1);
             }
 
             if(m_deltaFcnApprox)
             {
                  for(int j = 0; j < 2; ++j)
                  {
                 
                       m_waveVelocities[projectfield]->GetPlane(0)->BwdTrans(m_vwiForcing[j],der1);
                       for(int i = 0; i < npts; ++i)
                       {
                            der1[i] *= exp(-streak[i]*streak[i]/m_deltaFcnDecay);
                       }
                       m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(der1,m_vwiForcing[j]);
                  }
             }
             
             
             // dump output
             std::vector<std::string> variables(4);
             std::vector<Array<OneD, NekDouble> > outfield(4);
             variables[0] = "u";  variables[1] = "v"; 
             variables[2] = "pr"; variables[3] = "pi";
             outfield[0]  = m_vwiForcing[0];
             outfield[1]  = m_vwiForcing[1];
             Array<OneD,NekDouble> soln(npts,0.0);
             m_wavePressure->GetPlane(0)->BwdTrans(m_wavePressure->GetPlane(0)->GetCoeffs(),m_wavePressure->GetPlane(0)->UpdatePhys());
             outfield[2] = Array<OneD,NekDouble>(ncoeffs);
             m_waveVelocities[0]->GetPlane(0)->FwdTrans_IterPerExp(m_wavePressure->GetPlane(0)->GetPhys(),outfield[2]);
             m_wavePressure->GetPlane(1)->BwdTrans(m_wavePressure->GetPlane(1)->GetCoeffs(),m_wavePressure->GetPlane(1)->UpdatePhys());
             
             Vmath::Vmul(npts,m_wavePressure->GetPlane(0)->UpdatePhys(),1,
                         m_wavePressure->GetPlane(0)->UpdatePhys(),1,val,1);
             Vmath::Vvtvp(npts,m_wavePressure->GetPlane(1)->UpdatePhys(),1,
                          m_wavePressure->GetPlane(1)->UpdatePhys(),1,val,1,val,1);
             cout << "int P^2: " << m_wavePressure->GetPlane(0)->PhysIntegral(val) << endl;
             Vmath::Vsqrt(npts,val,1,val,1);
             cout << "PLinf: " <<  Vmath::Vmax(npts,val,1) << endl;
 
             outfield[3] = Array<OneD,NekDouble>(ncoeffs);
             m_waveVelocities[1]->GetPlane(0)->FwdTrans_IterPerExp(m_wavePressure->GetPlane(1)->GetPhys(),outfield[3]);
             
             std::string outname = m_sessionName  + ".vwi";
             
             m_solverRoll->WriteFld(outname, m_waveVelocities[0]->GetPlane(0), outfield, variables);
 
         }
     }

bool Nektar::VortexWaveInteraction::CheckEigIsStationary ( bool reset = false )

Definition at line 1603 of file VortexWaveInteraction.cpp.

References m_eigRelTol, m_leading_imag_evl, m_leading_real_evl, and m_minInnerIterations.

Referenced by DoFixedForcingIteration().

     {
         static NekDouble previous_real_evl = -1.0; 
         static NekDouble previous_imag_evl = -1.0; 
         static int min_iter = 0;
         
         if(reset)
         {
             previous_real_evl = -1.0;
             min_iter = 0;
         }
         
         if(previous_real_evl == -1.0 || min_iter < m_minInnerIterations)
         {
             previous_real_evl = m_leading_real_evl[0];
             previous_imag_evl = m_leading_imag_evl[0];
             min_iter++;
             return false;
         }
 
         cout << "Growth tolerance: " << fabs((m_leading_real_evl[0] - previous_real_evl)/m_leading_real_evl[0]) << endl; 
         cout << "Phase tolerance: " << fabs((m_leading_imag_evl[0] - previous_imag_evl)/m_leading_imag_evl[0]) << endl; 
 
         // See if real and imaginary growth have converged to with m_eigRelTol
         if((fabs((m_leading_real_evl[0] - previous_real_evl)/m_leading_real_evl[0]) < m_eigRelTol)||(fabs(m_leading_real_evl[0]) < 1e-6))
         {
             previous_real_evl = m_leading_real_evl[0];
             previous_imag_evl = m_leading_imag_evl[0];
             if((fabs((m_leading_imag_evl[0] - previous_imag_evl)/m_leading_imag_evl[0]) < m_eigRelTol)||(fabs(m_leading_imag_evl[0]) < 1e-6))
             {
                 return true;
             }
         }
         else
         {
             if(fabs(m_leading_imag_evl[0]) > 1e-2)
             {
                 cout << "Warning: imaginary eigenvalue is greater than 1e-2" << endl;
             }
             previous_real_evl = m_leading_real_evl[0];
             previous_imag_evl = m_leading_imag_evl[0];
             return false;
         }
         return false;
     }

bool Nektar::VortexWaveInteraction::CheckIfAtNeutralPoint ( void )

Definition at line 1651 of file VortexWaveInteraction.cpp.

References m_leading_imag_evl, m_leading_real_evl, m_neutralPointTol, m_sessionRoll, and m_sessionVWI.

Referenced by DoFixedForcingIteration().

     {
         bool returnval = false;
         if(m_sessionRoll->DefinesSolverInfo("INTERFACE"))
         {
               if(m_sessionVWI->GetSolverInfo("INTERFACE")=="phase")
               {
                   if( abs(m_leading_real_evl[0]) < 1e-4 )
                   {
                         returnval = true;
                   }
               }
               else
               {
                   if( abs(m_leading_real_evl[0]) < 1e-4  &&  abs(m_leading_imag_evl[0]) <2e-6 )
                   {
                         returnval = true;
                   }
               }
               
         }
         else
         {
                if((m_leading_real_evl[0]*m_leading_real_evl[0] + 
                   m_leading_imag_evl[0]*m_leading_imag_evl[0]) < 
                   m_neutralPointTol*m_neutralPointTol)
                {
                     returnval = true;
                }
         }
         if(fabs(m_leading_imag_evl[0]) > 1e-2)
         {
             cout << "Warning: imaginary eigenvalue is greater than 1e-2" << endl;
         }
 
         return returnval;
     }

void Nektar::VortexWaveInteraction::CopyFile	(	string	file1end,
		string	file2end
	)

Definition at line 1048 of file VortexWaveInteraction.cpp.

References ASSERTL0, and m_sessionName.

Referenced by ExecuteRoll(), ExecuteStreak(), and ExecuteWave().

      {
          string cpfile1   = m_sessionName + file1end;
          string cpfile2   = m_sessionName + file2end;
          string syscall  = "cp -f "  + cpfile1 + " " + cpfile2; 
 
          if(system(syscall.c_str()))
          {
              ASSERTL0(false,syscall.c_str());
          }
      }

void Nektar::VortexWaveInteraction::ExecuteLoop ( bool CalcWaveForce = true )

Definition at line 1096 of file VortexWaveInteraction.cpp.

References ASSERTL0, CalcNonLinearWaveForce(), Nektar::LibUtilities::NekManager< KeyType, ValueT, opLessCreator >::DisableManagement(), Nektar::eFixedWaveForcingWithSubIterationOnAlpha, Nektar::LibUtilities::NekManager< KeyType, ValueT, opLessCreator >::EnableManagement(), ExecuteRoll(), ExecuteStreak(), ExecuteWave(), GetVWIIterationType(), m_alpha, m_leading_imag_evl, m_moveMeshToCriticalLayer, m_sessionName, m_sessionRoll, and m_sessionVWI.

Referenced by DoFixedForcingIteration(), and main().

     {
 
 
     //the global info has to be written in the 
     //roll session file to use the interface loop
         if(m_sessionRoll->DefinesSolverInfo("INTERFACE"))
     {
              static int cnt=0;     
              bool skiprollstreak=false;
              if(cnt==0 && m_sessionVWI->GetParameter("rollstreakfromit")==1)
              {
                   skiprollstreak =true;
                   cout<<"skip roll-streak at the first iteration"<<endl;
              }
              
              
              if(skiprollstreak != true)
              {
 
                  LibUtilities::NekManager<MultiRegions::GlobalLinSysKey,MultiRegions::GlobalLinSys>::EnableManagement("GlobalLinSys");
                  ExecuteRoll();
                  LibUtilities::NekManager<MultiRegions::GlobalLinSysKey,MultiRegions::GlobalLinSys>::DisableManagement("GlobalLinSys");
 #ifndef _WIN32
              sleep(3);
 #endif
                  ExecuteStreak();
 #ifndef _WIN32
              sleep(3);
 #endif
              }
               
              string syscall;
              char c[16]="";
              string movedmesh = m_sessionName + "_advPost_moved.xml";
              string movedinterpmesh = m_sessionName + "_interp_moved.xml";
              //rewrite the Rollsessionfile (we start from the waleffe forcing)
              //string meshbndjumps = m_sessionName +"_bndjumps.xml";             
              //if(cnt==0)
              //{
                  //take the conditions tag from meshbndjumps and copy into 
                  // the rolls session file
              //}
 
 
              sprintf(c,"%d",cnt);  
              //save old roll solution
              string oldroll = m_sessionName +"_roll_"+c +".fld";    
              syscall = "cp -f " + m_sessionName+"-Base.fld" + "  " + oldroll;
              cout<<syscall.c_str()<<endl;
              if(system(syscall.c_str()))
              {
                   ASSERTL0(false,syscall.c_str());
              } 
              //define file names
              string filePost   = m_sessionName + "_advPost.xml";
              string filestreak   = m_sessionName + "_streak.fld";
              string filewave    = m_sessionName + "_wave.fld";
              string filewavepressure = m_sessionName + "_wave_p_split.fld";
              string fileinterp = m_sessionName + "_interp.xml";
              string interpstreak = m_sessionName +"_interpstreak_"+ c +".fld";  
              string interwavepressure  = m_sessionName +"_wave_p_split_interp_"+ c +".fld";
              char alpchar[16]="";
 cout<<"alpha = "<<m_alpha[0]<<endl;
              sprintf(alpchar, "%f", m_alpha[0]);
 
 
              if( m_sessionVWI->GetSolverInfo("INTERFACE")!="phase" )
              {
                  cout<<"zerophase"<<endl;
 
                  syscall  = "../../utilities/PostProcessing/Extras/MoveMesh  "
                              + filePost +"  "+ filestreak +"  "+ fileinterp + "   "+ alpchar; 
 
                  cout<<syscall.c_str()<<endl;
                  if(system(syscall.c_str()))
                  {
                       ASSERTL0(false,syscall.c_str());
                  }
 
                  //move the advPost mesh (remark update alpha!!!)
                  syscall  =  "../../utilities/PostProcessing/Extras/MoveMesh  "
                        + filePost + "  " + filestreak + "  " + filePost + "    "+ alpchar;
                  cout<<syscall.c_str()<<endl;
                  if(system(syscall.c_str()))
                  {
                      ASSERTL0(false,syscall.c_str());
                  }
 
 
 
                  //save oldstreak
                  string oldstreak = m_sessionName +"_streak_"+ c +".fld";            
              syscall = "cp -f " + filestreak + "  " + oldstreak;
                  cout<<syscall.c_str()<<endl;
                  if(system(syscall.c_str()))
                  {
                       ASSERTL0(false,syscall.c_str());
                  } 
 
                  //interpolate the streak field into the new mesh
                  string movedmesh = m_sessionName + "_advPost_moved.xml";
                  string movedinterpmesh = m_sessionName + "_interp_moved.xml";
 
                  //create the interp streak             
 
                  syscall  =  "../../utilities/PostProcessing/Extras/FieldToField  "
                       + fileinterp + "  " + filestreak + "  " + movedinterpmesh + "  " 
                   + interpstreak;
 
                  cout<<syscall.c_str()<<endl;
                  if(system(syscall.c_str()))
                  {
                       ASSERTL0(false,syscall.c_str());
                  } 
 
 
                  //save the old mesh     
                  string meshfile = m_sessionName + ".xml";                  
                  string meshold = m_sessionName +"_"+ c +".xml";
                  syscall = "cp -f " + meshfile + "  " + meshold;
                  cout<<syscall.c_str()<<endl;
                  if(system(syscall.c_str()))
                  {
                      ASSERTL0(false,syscall.c_str());
                  }  
 
                  //overwriting the meshfile with the new mesh
              syscall = "cp -f " + movedmesh + "  " + meshfile;
                  cout<<syscall.c_str()<<endl;
                  if(system(syscall.c_str()))
                  {
                      ASSERTL0(false,syscall.c_str());
                  }
 
                  //overwriting the streak file!!          
                  syscall = "cp -f " + interpstreak + "  " + filestreak;
                  cout<<syscall.c_str()<<endl;
                  if(system(syscall.c_str()))
                  {
                      ASSERTL0(false,syscall.c_str());
                  } 
 
                  //calculate the wave
                  ExecuteWave();
 
                  //save the wave field:
                  string oldwave = m_sessionName +"_wave_"+c +".fld";    
              syscall = "cp -f " + m_sessionName+".fld" + "  " + oldwave;
                  cout<<syscall.c_str()<<endl;
                  if(system(syscall.c_str()))
                  {
                       ASSERTL0(false,syscall.c_str());
                  } 
 
                  //save old jump conditions:
                  string ujump = m_sessionName+"_u_5.bc";
              syscall = "cp -f " + ujump + "  " + m_sessionName+"_u_5.bc_"+c;
                  cout<<syscall.c_str()<<endl;
                  if(system(syscall.c_str()))
                  {
                       ASSERTL0(false,syscall.c_str());
                  }              
 
                  string vjump = m_sessionName+"_v_5.bc";
                  syscall = "cp -f " + vjump + "  " + m_sessionName+"_v_5.bc_"+c;
                  cout<<syscall.c_str()<<endl;
                  if(system(syscall.c_str()))
                  {
                      ASSERTL0(false,syscall.c_str());
                  }    
                  cnt++;
 
 
 
 
                   //use relaxation
                   if(GetVWIIterationType()!=eFixedWaveForcingWithSubIterationOnAlpha 
                     )
                   {
                        // the critical layer should be the bnd region 3
                        //int reg =3;
                        //FileRelaxation(reg);
                   }
                   char c1[16]="";
                   sprintf(c1,"%d",cnt);   
                   //calculate the jump conditions
                   string wavefile  = m_sessionName +".fld"; 
                   syscall =  "../../utilities/PostProcessing/Extras/FldCalcBCs  "
                         + movedmesh + "  " + wavefile + "  " + interpstreak + ">  data"+c1;
                   cout<<syscall.c_str()<<endl;
                   if(system(syscall.c_str()))
                   {
                        ASSERTL0(false,syscall.c_str());
                   }
 
                   //move the new name_interp_moved.xml into name_interp.xml
               syscall = "cp -f " + movedinterpmesh + "  " + fileinterp;
                   cout<<syscall.c_str()<<endl;
                   if(system(syscall.c_str()))
                   {
                        ASSERTL0(false,syscall.c_str());
                   } 
                   //move the new name_advPost_moved.xml into name_advPost.xml
               syscall = "cp -f " + movedmesh + "  " + filePost;
                   cout<<syscall.c_str()<<endl;
                   if(system(syscall.c_str()))
                   {
                        ASSERTL0(false,syscall.c_str());
                   } 
  
 
              }
              else if(  m_sessionVWI->GetSolverInfo("INTERFACE")=="phase" )
              {    
 cout<<"phase"<<endl;
                   //determine cr:
                   NekDouble cr;
                   string cr_str;
                   stringstream st;
 
                   //calculate the wave
                   ExecuteWave();
    
                   //save oldstreak
                   string oldstreak = m_sessionName +"_streak_"+ c +".fld";            
               syscall = "cp -f " + filestreak + "  " + oldstreak;
                   cout<<syscall.c_str()<<endl;
                   if(system(syscall.c_str()))
                   {
                        ASSERTL0(false,syscall.c_str());
                   } 
 
                   //save wave
               syscall = "cp -f " + m_sessionName+".fld" + "  " + filewave;
                   cout<<syscall.c_str()<<endl;
                   if(system(syscall.c_str()))
                   {
                        ASSERTL0(false,syscall.c_str());
                   } 
                   //save the old mesh     
                   string meshfile = m_sessionName + ".xml";                  
                   string meshold = m_sessionName +"_"+ c +".xml";
               syscall = "cp -f " + meshfile + "  " + meshold;
                   cout<<syscall.c_str()<<endl;
                   if(system(syscall.c_str()))
                   {
                        ASSERTL0(false,syscall.c_str());
                   }  
 
                   //save the oldwave field:
                   string oldwave = m_sessionName +"_wave_"+c +".fld";    
                   syscall = "cp -f " + m_sessionName+".fld" + "  " + oldwave;
                   cout<<syscall.c_str()<<endl;
                   if(system(syscall.c_str()))
                   {
                        ASSERTL0(false,syscall.c_str());
                   }       
 
                   //save old jump conditions:
                   string ujump = m_sessionName+"_u_5.bc";
               syscall = "cp -f " + ujump + "  " + m_sessionName+"_u_5.bc_"+c;
                   cout<<syscall.c_str()<<endl;
                   if(system(syscall.c_str()))
                   {
                       ASSERTL0(false,syscall.c_str());
                   }              
 
                   string vjump = m_sessionName+"_v_5.bc";
                   syscall = "cp -f " + vjump + "  " + m_sessionName+"_v_5.bc_"+c;
                   cout<<syscall.c_str()<<endl;
                   if(system(syscall.c_str()))
                   {
                      ASSERTL0(false,syscall.c_str());
                   }
                   cnt++;
 
 
                   cr = m_leading_imag_evl[0]/m_alpha[0];
                   st << cr; 
                   cr_str = st.str();
 cout<<"cr="<<cr_str<<endl;
                   //NB -g or NOT!!!
                   //move the mesh around the critical layer    
                   syscall  = "../../utilities/PostProcessing/Extras/MoveMesh  "
                                + filePost +"  "+ filestreak +"  "+ fileinterp + "   "+ alpchar
                                +"      "+cr_str; 
   
                   cout<<syscall.c_str()<<endl;
                   if(system(syscall.c_str()))
                   {
                        ASSERTL0(false,syscall.c_str());
                   }
                   //NB -g or NOT!!!
                   //move the advPost mesh (remark update alpha!!!)
                   syscall  =  "../../utilities/PostProcessing/Extras/MoveMesh  "
                         + filePost + "  " + filestreak + "  " + filePost + "    "+ alpchar
                                +"      "+cr_str; 
                   cout<<syscall.c_str()<<endl;
                   if(system(syscall.c_str()))
                   {
                       ASSERTL0(false,syscall.c_str());
                   }
 
                   //interp streak into the new mesh
                   syscall  =  "../../utilities/PostProcessing/Extras/FieldToField  "
                         + fileinterp + "  " + filestreak + "  " + movedinterpmesh + "  " 
                     + interpstreak;
 
                   cout<<syscall.c_str()<<endl;
                   if(system(syscall.c_str()))
                   {
                       ASSERTL0(false,syscall.c_str());
                   }
 
                   //split wave sol
                   syscall  =  "../../utilities/PostProcessing/Extras/SplitFld  "
                         + filePost + "  " + filewave;
 
                   cout<<syscall.c_str()<<endl;
                   if(system(syscall.c_str()))
                   {
                       ASSERTL0(false,syscall.c_str());
                   }                
                   //interp wave
                   syscall  =  "../../utilities/PostProcessing/Extras/FieldToField  "
                         + filePost + "  " + filewavepressure + "  " + movedmesh + "  " 
                     + interwavepressure;
 
                   cout<<syscall.c_str()<<endl;
                   if(system(syscall.c_str()))
                   {
                       ASSERTL0(false,syscall.c_str());
                   }
 
 
 
 
                   //use relaxation
                   if(GetVWIIterationType()!=eFixedWaveForcingWithSubIterationOnAlpha 
                     )
                   {
                       // the critical layer should be the bnd region 3
                       //int reg =3;
                       //FileRelaxation(reg);
                   }
                   char c1[16]="";
                   sprintf(c1,"%d",cnt);   
 
                   //cp wavepressure to m_sessionName.fld(to get
                   // the right bcs names using FldCalcBCs)
                   syscall = "cp -f "+ interwavepressure +"  "+m_sessionName+".fld";
                   cout<<syscall.c_str()<<endl;
                   if(system(syscall.c_str()))
                   {
                       ASSERTL0(false,syscall.c_str());
                   }
 
                   //calculate the jump conditions
                   //NB -g or NOT!!!
                   syscall =  "../../utilities/PostProcessing/Extras/FldCalcBCs  "
                         + movedmesh + "  " +m_sessionName+".fld"  + "  " 
                         + interpstreak + ">  data"+c1;
                   cout<<syscall.c_str()<<endl;
                   if(system(syscall.c_str()))
                   {
                       ASSERTL0(false,syscall.c_str());
                   }
 
 
                   //overwriting the meshfile with the new mesh
               syscall = "cp -f " + movedmesh + "  " + meshfile;
                   cout<<syscall.c_str()<<endl;
                   if(system(syscall.c_str()))
                   {
                       ASSERTL0(false,syscall.c_str());
                   }
 
                   //overwriting the streak file!!    (maybe is useless)      
                   syscall = "cp -f " + interpstreak + "  " + filestreak;
                   cout<<syscall.c_str()<<endl;
                   if(system(syscall.c_str()))
                   {
                     ASSERTL0(false,syscall.c_str());
                   } 
                   //move the new name_interp_moved.xml into name_interp.xml
               syscall = "cp -f " + movedinterpmesh + "  " + fileinterp;
                   cout<<syscall.c_str()<<endl;
                   if(system(syscall.c_str()))
                   {
                        ASSERTL0(false,syscall.c_str());
                   } 
                   //move the new name_advPost_moved.xml into name_advPost.xml
               syscall = "cp -f " + movedmesh + "  " + filePost;
                   cout<<syscall.c_str()<<endl;
                   if(system(syscall.c_str()))
                   {
                        ASSERTL0(false,syscall.c_str());
                   }
              }
     }
     else
     {
             LibUtilities::NekManager<MultiRegions::GlobalLinSysKey,MultiRegions::GlobalLinSys>::EnableManagement("GlobalLinSys");
             ExecuteRoll();
             LibUtilities::NekManager<MultiRegions::GlobalLinSysKey,MultiRegions::GlobalLinSys>::DisableManagement("GlobalLinSys");
 
 #ifndef _WIN32
         sleep(3);
 #endif
             ExecuteStreak();
 #ifndef _WIN32
             sleep(3);
 #endif
             
             if(m_moveMeshToCriticalLayer)
             {
                 string syscall;
                 char alpchar[16]="";
                 sprintf(alpchar, "%f", m_alpha[0]);
 
                 string filePost          = m_sessionName + "_advPost.xml";
                 string filestreak        = m_sessionName + "_streak.fld";
                 string filewave          = m_sessionName + "_wave.fld";
                 string filewavepressure  = m_sessionName + "_wave_p_split.fld";
                 string fileinterp        = m_sessionName + "_interp.xml";
                 string interpstreak      = m_sessionName +"_interpstreak.fld";  
                 string interwavepressure = m_sessionName +"_wave_p_split_interp.fld";
                 syscall  = "../../utilities/PostProcessing/Extras/MoveMesh  "
                     + filePost +"  "+ filestreak +"  "+ fileinterp + "   "+ alpchar; 
                 
                 cout<<syscall.c_str()<<endl;
                 if(system(syscall.c_str()))
                 {
                     ASSERTL0(false,syscall.c_str());
                 }
                 
                 //move the advPost mesh (remark update alpha!!!)
                 syscall  =  "../../utilities/PostProcessing/Extras/MoveMesh  "
                     + filePost + "  " + filestreak + "  " + filePost + "    "+ alpchar;
                 cout<<syscall.c_str()<<endl;
                 if(system(syscall.c_str()))
                 {
                     ASSERTL0(false,syscall.c_str());
                 }
                 
                 //save oldstreak
                 string oldstreak = m_sessionName +"_streak.fld";
                 syscall = "cp -f " + filestreak + "  " + oldstreak;
                 cout<<syscall.c_str()<<endl;
                 if(system(syscall.c_str()))
                 {
                     ASSERTL0(false,syscall.c_str());
                 } 
                 
                 //interpolate the streak field into the new mesh
                 string movedmesh = m_sessionName + "_advPost_moved.xml";
                 string movedinterpmesh = m_sessionName + "_interp_moved.xml";
                 
                 //create the interp streak             
                 
                 syscall  =  "../../utilities/PostProcessing/Extras/FieldToField  "
                     + fileinterp + "  " + filestreak + "  " + movedinterpmesh 
                     + "  "  + interpstreak;
                 
                 cout<<syscall.c_str()<<endl;
                 if(system(syscall.c_str()))
                 {
                     ASSERTL0(false,syscall.c_str());
                 } 
                 
                 //save the old mesh     
                 string meshfile = m_sessionName + ".xml";                  
                 string meshold = m_sessionName + ".xml";
                 syscall = "cp -f " + meshfile + "  " + meshold;
                 cout<<syscall.c_str()<<endl;
                 if(system(syscall.c_str()))
                 {
                     ASSERTL0(false,syscall.c_str());
                 }  
                 
                 //overwriting the meshfile with the new mesh
                 syscall = "cp -f " + movedmesh + "  " + meshfile;
                 cout<<syscall.c_str()<<endl;
                 if(system(syscall.c_str()))
                 {
                     ASSERTL0(false,syscall.c_str());
                 }
                 
                 //overwriting the streak file!!          
                 syscall = "cp -f " + interpstreak + "  " + filestreak;
                 cout<<syscall.c_str()<<endl;
                 if(system(syscall.c_str()))
                 {
                     ASSERTL0(false,syscall.c_str());
                 } 
             }
             
             ExecuteWave();
 
             if(CalcWaveForce)
             {
                 CalcNonLinearWaveForce();
             }
     }
     }

void Nektar::VortexWaveInteraction::ExecuteRoll ( void )

Definition at line 351 of file VortexWaveInteraction.cpp.

References Nektar::IncNavierStokes::AddForcing(), CopyFile(), Nektar::LibUtilities::NekFactory< tKey, tBase, >::CreateInstance(), Nektar::SolverUtils::GetEquationSystemFactory(), Nektar::SolverUtils::GetForcingFactory(), m_moveMeshToCriticalLayer, m_rollField, m_rollForceScale, m_sessionName, m_sessionRoll, m_solverRoll, m_vwiForcing, m_vwiForcingObj, Vmath::Smul(), and Vmath::Vcopy().

Referenced by ExecuteLoop().

     {
         //set up the equation system to update the mesh
         if(m_sessionRoll->DefinesSolverInfo("INTERFACE"))
         {
             string vEquation = m_sessionRoll->GetSolverInfo("solvertype");
             EquationSystemSharedPtr solverRoll = GetEquationSystemFactory().CreateInstance(vEquation,m_sessionRoll);
             // The forcing terms are inserted as N bcs
             // Execute Roll 
             cout << "Executing Roll solver" << endl;
             solverRoll->DoInitialise();
             solverRoll->DoSolve();
             solverRoll->Output();
             m_rollField = solverRoll->UpdateFields();
             for(int g=0; g< solverRoll->GetNvariables(); ++g)
             {
                 NekDouble vL2Error = solverRoll->L2Error(g,false);
                 NekDouble vLinfError = solverRoll->LinfError(g);
                 cout << "L 2 error (variable " << solverRoll->GetVariable(g) << ") : " << vL2Error << endl;
                 cout << "L inf error (variable " << solverRoll->GetVariable(g) << ") : " << vLinfError << endl;
             }  
         }
         else
         {
             if(m_moveMeshToCriticalLayer)
             {
                 string vEquation = m_sessionRoll->GetSolverInfo("solvertype");
                 EquationSystemSharedPtr solverRoll = GetEquationSystemFactory().CreateInstance(vEquation,m_sessionRoll);
             }
             else
             {
                 const int npoints = m_solverRoll->GetNpoints();
                 const int ncoeffs = m_solverRoll->GetNcoeffs();
 
                 static int init = 1;
                 if(init)
                 {
                     m_solverRoll->DoInitialise();
                     m_vwiForcingObj = boost::dynamic_pointer_cast<SolverUtils::ForcingProgrammatic>(GetForcingFactory().CreateInstance("Programmatic", m_sessionRoll, m_solverRoll->UpdateFields(), m_solverRoll->UpdateFields().num_elements() - 1, 0));
 
                     std::vector<std::string> vFieldNames = m_sessionRoll->GetVariables();
                     vFieldNames.erase(vFieldNames.end()-1);
 
                     m_solverRoll->EvaluateFunction(vFieldNames, m_vwiForcingObj->UpdateForces(), "BodyForce");
 
                     // Scale forcing
                     for(int i = 0; i < m_vwiForcingObj->UpdateForces().num_elements(); ++i)
                     {
                         m_solverRoll->UpdateFields()[0]->FwdTrans(m_vwiForcingObj->UpdateForces()[i], m_vwiForcing[2+i]);
                         Vmath::Smul(npoints,m_rollForceScale,m_vwiForcingObj->UpdateForces()[i],1,m_vwiForcingObj->UpdateForces()[i],1);
                     }
 
                     IncNavierStokesSharedPtr ins = m_solverRoll->as<IncNavierStokes>();
                     ins->AddForcing(m_vwiForcingObj);
 
                     init = 0;
                 }
                 else // use internal definition of forcing in m_vwiForcing
                 {
                     // Scale forcing
                     for(int i = 0; i < m_vwiForcingObj->UpdateForces().num_elements(); ++i)
                     {
                         m_solverRoll->UpdateFields()[i]->BwdTrans(m_vwiForcing[i],m_vwiForcingObj->UpdateForces()[i]);
                         Vmath::Smul(npoints,m_rollForceScale,m_vwiForcingObj->UpdateForces()[i],1,m_vwiForcingObj->UpdateForces()[i],1);
                     }
                     
                     // Shift m_vwiForcing for new restart in case of relaxation 
                     Vmath::Vcopy(ncoeffs,m_vwiForcing[0],1,m_vwiForcing[2],1);
                     Vmath::Vcopy(ncoeffs,m_vwiForcing[1],1,m_vwiForcing[3],1);
                 }
             }
 
             // Execute Roll 
             cout << "Executing Roll solver" << endl;
             m_solverRoll->DoSolve();
             m_solverRoll->Output();
             m_rollField = m_solverRoll->UpdateFields();
             for(int g=0; g< m_solverRoll->GetNvariables(); ++g)
             {
                 NekDouble vL2Error = m_solverRoll->L2Error(g,false);
                 NekDouble vLinfError = m_solverRoll->LinfError(g);
                 cout << "L 2 error (variable " << m_solverRoll->GetVariable(g) << ") : " << vL2Error << endl;
                 cout << "L inf error (variable " << m_solverRoll->GetVariable(g) << ") : " << vLinfError << endl;
             }  
 
 
         }
         
         // Copy .fld file to .rst and base.fld
         cout << "Executing cp -f session.fld session.rst" << endl;
         CopyFile(".fld",".rst");
         
         // Write out data into base flow with variable Vx,Vy
         cout << "Writing data to session-Base.fld" << endl;
         
         std::vector<std::string> variables(2);
         variables[0] = "Vx";   variables[1] = "Vy";
         std::vector<Array<OneD, NekDouble> > outfield(2);
         outfield[0]  = m_solverRoll->UpdateFields()[0]->UpdateCoeffs(); 
         outfield[1]  = m_solverRoll->UpdateFields()[1]->UpdateCoeffs(); 
         std::string outname = m_sessionName  + "-Base.fld";
         m_solverRoll->WriteFld(outname, m_solverRoll->UpdateFields()[0], 
                                outfield, variables);
     }

void Nektar::VortexWaveInteraction::ExecuteStreak ( void )

Definition at line 457 of file VortexWaveInteraction.cpp.

References CopyFile(), Nektar::LibUtilities::NekFactory< tKey, tBase, >::CreateInstance(), Nektar::SolverUtils::GetDriverFactory(), Nektar::SolverUtils::GetEquationSystemFactory(), m_sessionStreak, m_sessionVWI, and m_streakField.

Referenced by ExecuteLoop().

     {
         // Create driver
 #if 1
         std::string vDriverModule;
         m_sessionStreak->LoadSolverInfo("Driver", vDriverModule, "Standard");
         
         DriverSharedPtr solverStreak = GetDriverFactory().CreateInstance(vDriverModule, m_sessionStreak); 
         solverStreak->Execute();
 
         m_streakField = solverStreak->GetEqu()[0]->UpdateFields();
 #else        
         // Setup and execute Advection Diffusion solver 
         string vEquation = m_sessionStreak->GetSolverInfo("EqType");
         EquationSystemSharedPtr solverStreak = GetEquationSystemFactory().CreateInstance(vEquation,m_sessionStreak);
 
         cout << "Executing Streak Solver" << endl;
         solverStreak->DoInitialise();
         solverStreak->DoSolve();
         solverStreak->Output();
 
         m_streakField = solverStreak->UpdateFields();
 
         if(m_sessionVWI->DefinesSolverInfo("INTERFACE"))
         {
             for(int g=0; g< solverStreak->GetNvariables(); ++g)
             {
                 NekDouble vL2Error = solverStreak->L2Error(g,false);
                 NekDouble vLinfError = solverStreak->LinfError(g);
                 cout << "L 2 error (variable " << solverStreak->GetVariable(g) << ") : " << vL2Error << endl;
                 cout << "L inf error (variable " << solverStreak->GetVariable(g) << ") : " << vLinfError << endl;
             }
         }
 #endif
 
         cout << "Executing cp -f session.fld session_streak.fld" << endl;
         CopyFile(".fld","_streak.fld");
     }

void Nektar::VortexWaveInteraction::ExecuteWave ( void )

Do linearised NavierStokes Session with Modified Arnoldi

Definition at line 496 of file VortexWaveInteraction.cpp.

References CopyFile(), Nektar::LibUtilities::NekFactory< tKey, tBase, >::CreateInstance(), Nektar::SolverUtils::GetDriverFactory(), m_alpha, m_leading_imag_evl, m_leading_real_evl, m_sessionVWI, m_sessionWave, m_wavePressure, and m_waveVelocities.

Referenced by CalcL2ToLinfPressure(), DoFixedForcingIteration(), ExecuteLoop(), and main().

     {
 
         // 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);
 
         // Create driver
         std::string vDriverModule;
         m_sessionWave->LoadSolverInfo("Driver", vDriverModule, "ModifiedArnoldi");
         cout << "Setting up linearised NS sovler" << endl;
         DriverSharedPtr solverWave = GetDriverFactory().CreateInstance(vDriverModule, m_sessionWave);  
 
         /// Do linearised NavierStokes Session  with Modified Arnoldi
         cout << "Executing wave solution " << endl;
         solverWave->Execute();
 
         // Copy file to a rst location for next restart
         cout << "Executing cp -f session_eig_0 session_eig_0.rst" << endl;
         CopyFile("_eig_0","_eig_0.rst");
 
         // Store data relevant to other operations 
         m_leading_real_evl[0] = solverWave->GetRealEvl()[0];
         m_leading_imag_evl[0] = solverWave->GetImagEvl()[0];
 
         // note this will only be true for modified Arnoldi
         NekDouble realShift = 0.0;
         if(m_sessionWave->DefinesParameter("RealShift"))
         {
             bool defineshift;
             // only use shift in modifiedArnoldi solver since
             // implicitly handled in Arpack.
             m_sessionWave->MatchSolverInfo("Driver","ModifiedArnoldi",defineshift,true);
             if(defineshift)
             {
                 realShift = m_sessionWave->GetParameter("RealShift");
             }
 
         }
 
         // Set up leading eigenvalue from inverse 
         NekDouble invmag = 1.0/(m_leading_real_evl[0]*m_leading_real_evl[0] + 
                                 m_leading_imag_evl[0]*m_leading_imag_evl[0]);
         m_leading_real_evl[0] *= -invmag;
         m_leading_real_evl[0] += realShift;
         m_leading_imag_evl[0] *= invmag;
 
 
         m_waveVelocities = solverWave->GetEqu()[0]->UpdateFields();
         m_wavePressure   = solverWave->GetEqu()[0]->GetPressure();
 
 
         if( m_sessionVWI->DefinesSolverInfo("INTERFACE")  )
         {
             cout << "Growth =" <<m_leading_real_evl[0]<<endl; 
             cout << "Phase =" <<m_leading_imag_evl[0]<<endl; 
         }
         
     }

void Nektar::VortexWaveInteraction::FileRelaxation ( int reg )

====================================================

Todo:: Please update to use MeshGraph::Read(vSession)

Definition at line 2003 of file VortexWaveInteraction.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), m_bcsForcing, m_rollField, m_sessionName, m_sessionVWI, m_vwiRelaxation, Nektar::SpatialDomains::MeshGraph::Read(), Vmath::Smul(), Vmath::Svtvp(), and Vmath::Vcopy().

     {
           cout<<"relaxation..."<<endl;
           static int cnt=0;
           Array<OneD, MultiRegions::ExpListSharedPtr> Iexp 
                                            =m_rollField[0]->GetBndCondExpansions();
           //cast to 1D explist (otherwise appenddata doesn't work)
           MultiRegions::ExpList1DSharedPtr Ilayer;  
           Ilayer = MemoryManager<MultiRegions::ExpList1D>::
                           AllocateSharedPtr(  
                           *boost::static_pointer_cast<MultiRegions::ExpList1D>(Iexp[reg]));
           int nq = Ilayer->GetTotPoints();
           if( cnt==0)
           {
                 m_bcsForcing = Array<OneD, Array<OneD, NekDouble> > (4);
                 m_bcsForcing[0] = Array<OneD, NekDouble> (4*nq);
                 for(int i = 1; i < 4; ++i)
                 {
                       m_bcsForcing[i] = m_bcsForcing[i-1] + nq;
                 }           
           }
 
           // Read in mesh from input file
 
           /// ====================================================
           /// \todo Please update to use MeshGraph::Read(vSession)
           /// ====================================================
           SpatialDomains::MeshGraphSharedPtr graphShPt = 
                                      SpatialDomains::MeshGraph::Read(m_sessionName+".xml");
 
 
           std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef_u;
           std::vector<std::vector<NekDouble> > FieldData_u;
           string file = m_sessionName;
 
 
           file += "_u_5.bc";
           LibUtilities::FieldIOSharedPtr fld =
               MemoryManager<LibUtilities::FieldIO>::AllocateSharedPtr(m_sessionVWI->GetComm());
           fld->Import(file,FieldDef_u, FieldData_u);
           Ilayer->ExtractDataToCoeffs(FieldDef_u[0], FieldData_u[0], FieldDef_u[0]->m_fields[0],Ilayer->UpdateCoeffs());
           Ilayer->BwdTrans_IterPerExp(Ilayer->GetCoeffs(), Ilayer->UpdatePhys());
           
           if(cnt==0)
           {
                Vmath::Vcopy(nq,Ilayer->UpdatePhys(),1,m_bcsForcing[2],1);
           }
           Vmath::Vcopy(nq,Ilayer->UpdatePhys(),1,m_bcsForcing[0],1);
 
           //case relaxation==0 an additional array is needed
           Array<OneD, NekDouble> tmp_forcing(nq, 0.0);
 
           if(cnt!=0)
           {
               if(m_vwiRelaxation==1.0)
               {
                  Vmath::Vcopy(nq, m_bcsForcing[0],1, tmp_forcing,1);
               }
               Vmath::Smul(nq,1.0-m_vwiRelaxation,
                         m_bcsForcing[0],1,m_bcsForcing[0],1);
               Vmath::Svtvp(nq,m_vwiRelaxation,m_bcsForcing[2],1,
                          m_bcsForcing[0],1,Ilayer->UpdatePhys(),1);
               //generate again the bcs files:
               
               Array<OneD, Array<OneD, NekDouble> > fieldcoeffs(1);   
               Ilayer->FwdTrans_IterPerExp(Ilayer->GetPhys(),Ilayer->UpdateCoeffs()); 
               fieldcoeffs[0] = Ilayer->UpdateCoeffs();      
           std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef1  = Ilayer->GetFieldDefinitions();               
               std::vector<std::vector<NekDouble> > FieldData_1(FieldDef1.size());;
               FieldDef1[0]->m_fields.push_back("u");                    
               Ilayer->AppendFieldData(FieldDef1[0], FieldData_1[0]);                    
               fld->Write(file,FieldDef1,FieldData_1);
               //save the bcs for the next iteration
               if(m_vwiRelaxation!=1.0)
               {
                    Vmath::Smul(nq,1./(1.0-m_vwiRelaxation),
                         m_bcsForcing[0],1,m_bcsForcing[0],1);              
                    Vmath::Vcopy(nq,m_bcsForcing[0],1,m_bcsForcing[2],1);
               }
               else
               {
                    Vmath::Vcopy(nq, tmp_forcing,1, m_bcsForcing[2],1);                   
               }
           }
                    
 
 
           file = m_sessionName+ "_v_5.bc"; 
 
           std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef_v;
           std::vector<std::vector<NekDouble> > FieldData_v;
           fld->Import(file,FieldDef_v, FieldData_v);
           Ilayer->ExtractDataToCoeffs(FieldDef_v[0], FieldData_v[0], FieldDef_v[0]->m_fields[0],Ilayer->UpdateCoeffs());
           Ilayer->BwdTrans_IterPerExp(Ilayer->GetCoeffs(), Ilayer->UpdatePhys());
           if(cnt==0)
           {
                Vmath::Vcopy(nq,Ilayer->UpdatePhys(),1,m_bcsForcing[3],1);
           }
           Vmath::Vcopy(nq,Ilayer->UpdatePhys(),1,m_bcsForcing[1],1);
           if(cnt!=0)
           {
               if(m_vwiRelaxation==1.0)
               {
                  Vmath::Vcopy(nq, m_bcsForcing[1],1, tmp_forcing,1);
               }
               Vmath::Smul(nq,1.0-m_vwiRelaxation,
                         m_bcsForcing[1],1,m_bcsForcing[1],1);
               Vmath::Svtvp(nq,m_vwiRelaxation,m_bcsForcing[3],1,
                          m_bcsForcing[1],1,Ilayer->UpdatePhys(),1);
               //generate again the bcs files:
               Array<OneD, Array<OneD, NekDouble> > fieldcoeffs(1);   
               Ilayer->FwdTrans_IterPerExp(Ilayer->GetPhys(),Ilayer->UpdateCoeffs()); 
               fieldcoeffs[0] = Ilayer->UpdateCoeffs();      
           std::vector<LibUtilities::FieldDefinitionsSharedPtr>  FieldDef2  = Ilayer->GetFieldDefinitions();         
               std::vector<std::vector<NekDouble> > FieldData_2(FieldDef2.size());;      
               FieldDef2[0]->m_fields.push_back("v");                    
               Ilayer->AppendFieldData(FieldDef2[0], FieldData_2[0]);                                
               fld->Write(file,FieldDef2,FieldData_2);
               //save the bcs for the next iteration
               if(m_vwiRelaxation!=1.0)
               {
                   Vmath::Smul(nq,1./(1.0-m_vwiRelaxation),
                               m_bcsForcing[1],1,m_bcsForcing[1],1);              
                   Vmath::Vcopy(nq,m_bcsForcing[1],1,m_bcsForcing[3],1);
               }
               else
               {
                   Vmath::Vcopy(nq, tmp_forcing,1, m_bcsForcing[3],1);                   
               }
 
 
           }
 
 
            cnt++;
     }

NekDouble Nektar::VortexWaveInteraction::GetAlpha ( void )

inline

Definition at line 130 of file VortexWaveInteraction.h.

Referenced by DoFixedForcingIteration(), and main().

         {
             return m_alpha[0];
         }

NekDouble Nektar::VortexWaveInteraction::GetAlphaStep ( void )

inline

Definition at line 136 of file VortexWaveInteraction.h.

Referenced by main().

         {
             return m_alphaStep;
         }

NekDouble Nektar::VortexWaveInteraction::GetDAlphaDWaveForceMag ( void )

inline

Definition at line 151 of file VortexWaveInteraction.h.

         {
             return m_dAlphaDWaveForceMag; 
         }

NekDouble Nektar::VortexWaveInteraction::GetEigRelTol ( void )

inline

Definition at line 162 of file VortexWaveInteraction.h.

Referenced by DoFixedForcingIteration().

     {
             return m_eigRelTol;
     }

int Nektar::VortexWaveInteraction::GetIterEnd ( )

inline

Definition at line 110 of file VortexWaveInteraction.h.

Referenced by DoFixedForcingIteration().

         {
             return m_iterEnd;
         }

int Nektar::VortexWaveInteraction::GetIterStart ( )

inline

Definition at line 105 of file VortexWaveInteraction.h.

Referenced by DoFixedForcingIteration().

         {
             return m_iterStart;
         }

int Nektar::VortexWaveInteraction::GetMaxOuterIterations ( void )

inline

Definition at line 125 of file VortexWaveInteraction.h.

Referenced by DoFixedForcingIteration().

         {
             return m_maxOuterIterations;
         }

int Nektar::VortexWaveInteraction::GetMaxWaveForceMagIter ( void )

inline

Definition at line 156 of file VortexWaveInteraction.h.

Referenced by main().

         {
             return m_maxWaveForceMagIter;
         }

int Nektar::VortexWaveInteraction::GetMinInnerIterations ( void )

inline

Definition at line 168 of file VortexWaveInteraction.h.

Referenced by DoFixedForcingIteration().

     {
             return  m_minInnerIterations;
     }

int Nektar::VortexWaveInteraction::GetNOuterIterations ( void )

inline

Definition at line 120 of file VortexWaveInteraction.h.

Referenced by DoFixedForcingIteration().

         {
             return m_nOuterIterations;
         }

NekDouble Nektar::VortexWaveInteraction::GetPrevAlpha ( void )

inline

Definition at line 173 of file VortexWaveInteraction.h.

Referenced by DoFixedForcingIteration().

         {
             return m_prevAlpha;
         }

Array< OneD, int > Nektar::VortexWaveInteraction::GetReflectionIndex ( void )

Definition at line 1923 of file VortexWaveInteraction.cpp.

References ASSERTL0, Nektar::LibUtilities::eTriangle, m_sessionVWI, m_waveVelocities, npts, and Vmath::Vmax().

Referenced by CalcNonLinearWaveForce().

     {
         int i,j;
         int npts = m_waveVelocities[0]->GetPlane(0)->GetNpoints();
         int nel  = m_waveVelocities[0]->GetNumElmts();
         Array<OneD, int> index(npts);
 
         Array<OneD, NekDouble> coord(2);
         Array<OneD, NekDouble> coord_x(npts);
         Array<OneD, NekDouble> coord_y(npts);
         
         //-> Dermine the point which is on coordinate (x -> -x + Lx/2, y-> -y)
         m_waveVelocities[0]->GetPlane(0)->GetCoords(coord_x,coord_y);
         NekDouble xmax = Vmath::Vmax(npts,coord_x,1);
         //NekDouble tol = NekConstants::kGeomFactorsTol*NekConstants::kGeomFactorsTol;
         NekDouble tol = 1e-5;
         NekDouble xnew,ynew;
 
         int start  = npts-1; 
         int e_npts;
 
         bool useOnlyQuads = false;
         if(m_sessionVWI->DefinesSolverInfo("SymmetriseOnlyQuads"))
         {
             useOnlyQuads = true;
         }
         
         int cnt;
         for(int e = 0; e < nel; ++e)
         {
             e_npts = m_waveVelocities[0]->GetExp(e)->GetTotPoints();
             cnt = m_waveVelocities[0]->GetPhys_Offset(e);
             
             if(useOnlyQuads)
             {
                 if(m_waveVelocities[0]->GetExp(e)->DetShapeType() == LibUtilities::eTriangle)
                 {
                     for(i = 0; i < e_npts; ++i)
                     {
                         index[cnt+i] = -1;
                     }
                     continue;
                 }
             }
             
             for(i = cnt; i < cnt+e_npts; ++i)
             {
                 xnew = - coord_x[i]  + xmax;
                 ynew = - coord_y[i];
                 
                 for(j = start; j >=0 ; --j)
                 {
                     if((coord_x[j]-xnew)*(coord_x[j]-xnew) + (coord_y[j]-ynew)*(coord_y[j]-ynew) < tol)
                     {
                         index[i] = j;
                         start = j;
                         break;
                     }
                 }
                 
                 if(j == -1)
                 {
                     
                     for(j = npts-1; j > start; --j)
                     {
                         
                         if((coord_x[j]-xnew)*(coord_x[j]-xnew) + (coord_y[j]-ynew)*(coord_y[j]-ynew) < tol)
                         {
                             index[i] = j;
                             break;
                         }
                     }
                     ASSERTL0(j != start,"Failed to find matching point");
                 }
             }
         }
         return index;
     }

VWIIterationType Nektar::VortexWaveInteraction::GetVWIIterationType ( void )

inline

Definition at line 115 of file VortexWaveInteraction.h.

Referenced by DoFixedForcingIteration(), ExecuteLoop(), and main().

         {
             return m_VWIIterationType;
         }

NekDouble Nektar::VortexWaveInteraction::GetWaveForceMag ( void )

inline

Definition at line 141 of file VortexWaveInteraction.h.

Referenced by main().

         {
             return m_waveForceMag[0];
         }

NekDouble Nektar::VortexWaveInteraction::GetWaveForceMagStep ( void )

inline

Definition at line 146 of file VortexWaveInteraction.h.

Referenced by main().

         {
             return m_waveForceMagStep;
         }

bool Nektar::VortexWaveInteraction::IfIterInterface ( void )

inline

Definition at line 209 of file VortexWaveInteraction.h.

Referenced by DoFixedForcingIteration().

         {
             return m_iterinterface;
         }

void Nektar::VortexWaveInteraction::MoveFile	(	string	fileend,
		string	dir,
		int	n
	)

Definition at line 1027 of file VortexWaveInteraction.cpp.

References ASSERTL0.

Referenced by DoFixedForcingIteration().

     {
         static map<string,int> opendir;
 
         if(opendir.find(dir) == opendir.end())
         {
             // make directory and presume will fail if it already exists
             string mkdir = "mkdir " + dir;
             system(mkdir.c_str());
             opendir[dir] = 1;
         }
         
         string savefile = dir + "/" + file + "." + boost::lexical_cast<std::string>(n);
         string syscall  = "mv -f "  + file + " " + savefile; 
 
         if(system(syscall.c_str()))
         {
             ASSERTL0(false,syscall.c_str());
         }
      }

void Nektar::VortexWaveInteraction::SaveFile	(	string	fileend,
		string	dir,
		int	n
	)

Definition at line 1003 of file VortexWaveInteraction.cpp.

References ASSERTL0.

Referenced by SaveLoopDetails().

     {
         static map<string,int> opendir;
 
         if(opendir.find(dir) == opendir.end())
         {
             // make directory and presume will fail if it already exists
             string mkdir = "mkdir " + dir;
             system(mkdir.c_str());
 
             opendir[dir] = 1;
         }
         
         string savefile = dir + "/" + file + "." + boost::lexical_cast<std::string>(n);
         string syscall  = "cp -f "  + file + " " + savefile; 
 
         if(system(syscall.c_str()))
         {
             ASSERTL0(false,syscall.c_str());
         }
 
      }

void Nektar::VortexWaveInteraction::SaveLoopDetails	(	std::string	SaveDir,
		int	i
	)

Definition at line 1076 of file VortexWaveInteraction.cpp.

References m_sessionName, m_sessionVWI, and SaveFile().

Referenced by DoFixedForcingIteration().

     {
         // Save NS restart file
         SaveFile(m_sessionName + ".rst",SaveDir,i+1);
         // Save Streak Solution
         SaveFile(m_sessionName + "_streak.fld",SaveDir,i);
         // Save Wave solution output
         SaveFile(m_sessionName + ".evl",SaveDir,i);
         SaveFile(m_sessionName + "_eig_0",SaveDir,i);
         // Save new field file of eigenvalue
         SaveFile(m_sessionName + ".fld",SaveDir,i);
         if(!(m_sessionVWI->DefinesSolverInfo("INTERFACE")))
         {
             // Save new forcing file
             SaveFile(m_sessionName + ".vwi",SaveDir,i+1);
         }
 
     }

void Nektar::VortexWaveInteraction::SetAlpha ( NekDouble alpha )

inline

Definition at line 178 of file VortexWaveInteraction.h.

Referenced by main().

         {
             m_alpha[0] = alpha; 
         }

void Nektar::VortexWaveInteraction::SetAlphaStep ( NekDouble step )

inline

Definition at line 194 of file VortexWaveInteraction.h.

         {
         m_alphaStep = step;
         }

void Nektar::VortexWaveInteraction::SetEigRelTol ( NekDouble tol )

inline

Definition at line 189 of file VortexWaveInteraction.h.

Referenced by DoFixedForcingIteration().

     {
         m_eigRelTol = tol;
     }

void Nektar::VortexWaveInteraction::SetMinInnerIterations ( int niter )

inline

Definition at line 199 of file VortexWaveInteraction.h.

Referenced by DoFixedForcingIteration().

     {
       m_minInnerIterations = niter;
     }

void Nektar::VortexWaveInteraction::SetPrevAlpha ( NekDouble alpha )

inline

Definition at line 204 of file VortexWaveInteraction.h.

Referenced by main().

         {
             m_prevAlpha = alpha; 
         }

void Nektar::VortexWaveInteraction::SetWaveForceMag ( NekDouble mag )

inline

Definition at line 184 of file VortexWaveInteraction.h.

Referenced by main().

         {
             m_waveForceMag[0] = mag; 
         }

void Nektar::VortexWaveInteraction::UpdateAlpha ( int n )

Definition at line 1810 of file VortexWaveInteraction.cpp.

References Vmath::Imin(), m_alpha, m_alphaStep, m_leading_imag_evl, m_leading_real_evl, and Vmath::Vcopy().

Referenced by DoFixedForcingIteration(), and VortexWaveInteraction().

     {
         NekDouble alp_new = 0.0;
 
 
         if(outeriter == 1)
         {
             m_alpha[1] = m_alpha[0];
             if(m_leading_real_evl[0] > 0.0)
             {
                 alp_new = m_alpha[0] + m_alphaStep;
             }
             else
             {
                 alp_new = m_alpha[0] - m_alphaStep;
             }
         }
         else
         {
             int i;
             int nstore = (m_alpha.num_elements() < outeriter)? m_alpha.num_elements(): outeriter;
             Array<OneD, NekDouble> Alpha(nstore);
             Array<OneD, NekDouble> Growth(nstore);
             
             Vmath::Vcopy(nstore,m_alpha,1,Alpha,1);
             Vmath::Vcopy(nstore,m_leading_real_evl,1,Growth,1);
             
             // Sort Alpha Growth values; 
             NekDouble store;
             int k;
             for(i = 0; i < nstore; ++i)
             {
                 k = Vmath::Imin(nstore-i,&Alpha[i],1);
                 
                 store    = Alpha[i]; 
                 Alpha[i] = Alpha[i+k];
                 Alpha[i+k] = store;
 
                 store     = Growth[i];
                 Growth[i] = Growth[i+k];
                 Growth[i+k] = store; 
             }
 
             // See if we have any values that cross zero
             for(i = 0; i < nstore-1; ++i)
             {
                 if(Growth[i]*Growth[i+1] < 0.0)
                 {
                     break;
                 }
             }
             
             if(i != nstore-1)
             {
                 if(nstore == 2)
                 {
                     alp_new = (Alpha[0]*Growth[1] - Alpha[1]*Growth[0])/(Growth[1]-Growth[0]);
                 }
                 else
                 {
                     // use a quadratic fit and step through 10000 points
                     // to find zero.
                     int     j; 
                     int     nsteps = 10000;
                     int     idx = (i == 0)?1:i;
                     NekDouble  da = Alpha[idx+1] - Alpha[idx-1];
                     NekDouble  gval_m1 = Growth[idx-1],a,gval;
                     NekDouble  c1 = Growth[idx-1]/(Alpha[idx-1]-Alpha[idx])/
                         (Alpha[idx-1]-Alpha[idx+1]);
                     NekDouble  c2 = Growth[idx]/(Alpha[idx]-Alpha[idx-1])/
                         (Alpha[idx]-Alpha[idx+1]);
                     NekDouble  c3 = Growth[idx+1]/(Alpha[idx+1]-Alpha[idx-1])/
                         (Alpha[idx+1]-Alpha[idx]);
                     
                     for(j = 1; j < nsteps+1; ++j)
                     {
                         a = Alpha[i] + j*da/(NekDouble) nsteps;
                         gval = c1*(a - Alpha[idx  ])*(a - Alpha[idx+1]) 
                             +  c2*(a - Alpha[idx-1])*(a - Alpha[idx+1])
                             +  c3*(a - Alpha[idx-1])*(a - Alpha[idx]);
                         
                         if(gval*gval_m1 < 0.0)
                         {
                             alp_new = ((a+da/(NekDouble)nsteps)*gval - a*gval_m1)/
                                 (gval - gval_m1);
                             break;
                         }
                     }
                 }
             }
             else // step backward/forward 
             {
                 if(Growth[i] > 0.0)
                 {
                     alp_new = m_alpha[0] + m_alphaStep;
                 }
                 else
                 {
                     alp_new = m_alpha[0] - m_alphaStep;
                 }
             }
         }
         
         for(int i = m_alpha.num_elements()-1; i > 0; --i)
         {
             m_alpha[i] = m_alpha[i-1];
             m_leading_real_evl[i] = m_leading_real_evl[i-1];
             m_leading_imag_evl[i] = m_leading_imag_evl[i-1];
         }
 
         m_alpha[0] = alp_new;    
     }

void Nektar::VortexWaveInteraction::UpdateDAlphaDWaveForceMag ( NekDouble alphainit )

Definition at line 1805 of file VortexWaveInteraction.cpp.

References m_alpha, m_dAlphaDWaveForceMag, and m_waveForceMagStep.

Referenced by DoFixedForcingIteration().

     {
         m_dAlphaDWaveForceMag = (m_alpha[0]-alpha_init)/m_waveForceMagStep;
     }

void Nektar::VortexWaveInteraction::UpdateWaveForceMag ( int n )

Definition at line 1691 of file VortexWaveInteraction.cpp.

References Vmath::Imin(), m_leading_imag_evl, m_leading_real_evl, m_waveForceMag, m_waveForceMagStep, and Vmath::Vcopy().

Referenced by DoFixedForcingIteration().

     {
         NekDouble wavef_new = 0.0;
 
 
         if(outeriter == 1)
         {
             m_waveForceMag[1] = m_waveForceMag[0];
             if(m_leading_real_evl[0] > 0.0)
             {
                 wavef_new = m_waveForceMag[0] - m_waveForceMagStep;
             }
             else
             {
                 wavef_new = m_waveForceMag[0] + m_waveForceMagStep;
             }
         }
         else
         {
             int i;
             int nstore = (m_waveForceMag.num_elements() < outeriter)? m_waveForceMag.num_elements(): outeriter;
             Array<OneD, NekDouble> WaveForce(nstore);
             Array<OneD, NekDouble> Growth(nstore);
             
             Vmath::Vcopy(nstore,m_waveForceMag,1,WaveForce,1);
             Vmath::Vcopy(nstore,m_leading_real_evl,1,Growth,1);
             
             // Sort WaveForce Growth values; 
             NekDouble store;
             int k;
             for(i = 0; i < nstore; ++i)
             {
                 k = Vmath::Imin(nstore-i,&WaveForce[i],1);
                 
                 store    = WaveForce[i]; 
                 WaveForce[i] = WaveForce[i+k];
                 WaveForce[i+k] = store;
 
                 store     = Growth[i];
                 Growth[i] = Growth[i+k];
                 Growth[i+k] = store; 
             }
 
             // See if we have any values that cross zero
             for(i = 0; i < nstore-1; ++i)
             {
                 if(Growth[i]*Growth[i+1] < 0.0)
                 {
                     break;
                 }
             }
             
             if(i != nstore-1)
             {
                 if(nstore == 2)
                 {
                     wavef_new = (WaveForce[0]*Growth[1] - WaveForce[1]*Growth[0])/(Growth[1]-Growth[0]);
                 }
                 else
                 {
                     // use a quadratic fit and step through 10000 points
                     // to find zero.
                     int     j; 
                     int     nsteps = 10000;
                     int     idx = (i == 0)?1:i;
                     NekDouble  da = WaveForce[idx+1] - WaveForce[idx-1];
                     NekDouble  gval_m1 = Growth[idx-1],a,gval;
                     NekDouble  c1 = Growth[idx-1]/(WaveForce[idx-1]-WaveForce[idx])/
                         (WaveForce[idx-1]-WaveForce[idx+1]);
                     NekDouble  c2 = Growth[idx]/(WaveForce[idx]-WaveForce[idx-1])/
                         (WaveForce[idx]-WaveForce[idx+1]);
                     NekDouble  c3 = Growth[idx+1]/(WaveForce[idx+1]-WaveForce[idx-1])/
                         (WaveForce[idx+1]-WaveForce[idx]);
                     
                     for(j = 1; j < nsteps+1; ++j)
                     {
                         a = WaveForce[i] + j*da/(NekDouble) nsteps;
                         gval = c1*(a - WaveForce[idx  ])*(a - WaveForce[idx+1]) 
                             +  c2*(a - WaveForce[idx-1])*(a - WaveForce[idx+1])
                             +  c3*(a - WaveForce[idx-1])*(a - WaveForce[idx]);
                         
                         if(gval*gval_m1 < 0.0)
                         {
                             wavef_new = ((a+da/(NekDouble)nsteps)*gval - a*gval_m1)/
                                 (gval - gval_m1);
                             break;
                         }
                     }
                 }
             }
             else // step backward/forward 
             {
                 if(Growth[i] > 0.0)
                 {
                     wavef_new = m_waveForceMag[0] - m_waveForceMagStep;
                 }
                 else
                 {
                     wavef_new = m_waveForceMag[0] + m_waveForceMagStep;
                 }
             }
         }
         
         for(int i = m_waveForceMag.num_elements()-1; i > 0; --i)
         {
             m_waveForceMag[i] = m_waveForceMag[i-1];
             m_leading_real_evl[i] = m_leading_real_evl[i-1];
             m_leading_imag_evl[i] = m_leading_imag_evl[i-1];
         }
 
         m_waveForceMag[0] = wavef_new;
         
     }