Nektar::VortexWaveInteraction Class Reference

#include <VortexWaveInteraction.h>

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 (std::string dir, int i)
void	CalcNonLinearWaveForce (void)
void	CalcL2ToLinfPressure (void)
void	SaveFile (std::string fileend, std::string dir, int n)
void	MoveFile (std::string fileend, std::string dir, int n)
void	CopyFile (std::string file1end, std::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
std::string	m_sessionName
LibUtilities::SessionReaderSharedPtr	m_sessionVWI
LibUtilities::SessionReaderSharedPtr	m_sessionRoll
SpatialDomains::MeshGraphSharedPtr	m_graphRoll
EquationSystemSharedPtr	m_solverRoll
LibUtilities::SessionReaderSharedPtr	m_sessionStreak
SpatialDomains::MeshGraphSharedPtr	m_graphStreak
LibUtilities::SessionReaderSharedPtr	m_sessionWave
SpatialDomains::MeshGraphSharedPtr	m_graphWave

Detailed Description

Definition at line 66 of file VortexWaveInteraction.h.

Constructor & Destructor Documentation

VortexWaveInteraction()

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

Definition at line 47 of file VortexWaveInteraction.cpp.

    : 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->InitSession();
 
    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());
    m_graphRoll           = SpatialDomains::MeshGraphIO::Read(m_sessionRoll);
    std::string vEquation = m_sessionRoll->GetSolverInfo("SolverType");
    m_solverRoll          = GetEquationSystemFactory().CreateInstance(
        vEquation, m_sessionRoll, m_graphRoll);
    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());
    m_graphStreak = SpatialDomains::MeshGraphIO::Read(m_sessionStreak);
 
    // 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());
    m_graphWave = SpatialDomains::MeshGraphIO::Read(m_sessionWave);
 
    // 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 (boost::iequals(VWIIterationTypeMap[i], IterationTypeStr))
            {
                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.size());
 
                    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 = std::to_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]);
                }
            }
        }
    }
}

References ASSERTL0, Nektar::LibUtilities::SessionReader::CreateInstance(), Nektar::LibUtilities::NekFactory< tKey, tBase, tParam >::CreateInstance(), Nektar::eFixedAlpha, Nektar::eFixedAlphaWaveForcing, Nektar::eFixedWaveForcing, Nektar::eVWIIterationTypeSize, Nektar::SolverUtils::GetEquationSystemFactory(), tinysimd::log(), m_alpha, m_alphaStep, m_dAlphaDWaveForceMag, m_deltaFcnApprox, m_deltaFcnDecay, m_eigRelTol, m_graphRoll, m_graphStreak, m_graphWave, 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, Nektar::SpatialDomains::MeshGraphIO::Read(), UpdateAlpha(), and Nektar::VWIIterationTypeMap.

~VortexWaveInteraction()

Nektar::VortexWaveInteraction::~VortexWaveInteraction

(

void

)

Definition at line 363 of file VortexWaveInteraction.cpp.

{
    m_sessionVWI->Finalise();
}

References m_sessionVWI.

Member Function Documentation

AppendEvlToFile() [1/2]

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

Definition at line 1162 of file VortexWaveInteraction.cpp.

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

References m_alpha, m_leading_imag_evl, and m_leading_real_evl.

Referenced by DoFixedForcingIteration(), and main().

AppendEvlToFile() [2/2]

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

Definition at line 1171 of file VortexWaveInteraction.cpp.

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

References m_alpha, m_leading_imag_evl, and m_leading_real_evl.

CalcL2ToLinfPressure()

void Nektar::VortexWaveInteraction::CalcL2ToLinfPressure

(

void

)

Definition at line 1068 of file VortexWaveInteraction.cpp.

{
    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)->Integral(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)->Integral(val);
 
    l2 = sqrt(norm / area);
 
    cout << "L2:   " << l2 << endl;
 
    cout << "Ratio Linf/L2: " << Linf / l2 << endl;
}

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

Referenced by main().

CalcNonLinearWaveForce()

void Nektar::VortexWaveInteraction::CalcNonLinearWaveForce

(

void

)

Definition at line 601 of file VortexWaveInteraction.cpp.

{
    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";
        string syscall    = "";
        string c          = std::to_string(cnt);
        string c_alpha    = std::to_string(m_alpha[0]);
        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.size(); ++i)
            {
                BndConds = m_waveVelocities[i]->GetBndConditions();
                for (int j = 0; j < BndConds.size(); ++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)->Integral(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(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(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(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(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)->FwdTransLocalElmt(
            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)->Integral(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)->FwdTransLocalElmt(
            m_wavePressure->GetPlane(1)->GetPhys(), outfield[3]);
 
        std::string outname = m_sessionName + ".vwi";
 
        m_solverRoll->WriteFld(outname, m_waveVelocities[0]->GetPlane(0),
                               outfield, variables);
    }
}

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(), Vmath::Sadd(), Vmath::Smul(), tinysimd::sqrt(), Vmath::Svtvp(), Vmath::Vadd(), Vmath::Vcopy(), Vmath::Vmax(), Vmath::Vmul(), Vmath::Vsqrt(), Vmath::Vsub(), and Vmath::Vvtvp().

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

CheckEigIsStationary()

bool Nektar::VortexWaveInteraction::CheckEigIsStationary

(

bool

reset = false

)

Definition at line 1701 of file VortexWaveInteraction.cpp.

{
    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;
}

References m_eigRelTol, m_leading_imag_evl, m_leading_real_evl, and m_minInnerIterations.

Referenced by DoFixedForcingIteration().

CheckIfAtNeutralPoint()

bool Nektar::VortexWaveInteraction::CheckIfAtNeutralPoint

(

void

)

Definition at line 1760 of file VortexWaveInteraction.cpp.

{
    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;
}

References tinysimd::abs(), m_leading_imag_evl, m_leading_real_evl, m_neutralPointTol, m_sessionRoll, and m_sessionVWI.

Referenced by DoFixedForcingIteration().

CopyFile()

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

Definition at line 1150 of file VortexWaveInteraction.cpp.

{
    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());
    }
}

References ASSERTL0, and m_sessionName.

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

ExecuteLoop()

void Nektar::VortexWaveInteraction::ExecuteLoop

(

bool

CalcWaveForce = true

)

Definition at line 1199 of file VortexWaveInteraction.cpp.

{
 
    // 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;
        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
        //}
 
        string c = std::to_string(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";
        string c_alpha = std::to_string(m_alpha[0]);
        cout << "alpha = " << m_alpha[0] << endl;
 
        if (m_sessionVWI->GetSolverInfo("INTERFACE") != "phase")
        {
            cout << "zerophase" << endl;
 
            syscall = "../../utilities/PostProcessing/Extras/MoveMesh  " +
                      filePost + "  " + filestreak + "  " + fileinterp + "   " +
                      c_alpha;
 
            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 + "    " +
                      c_alpha;
            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++;
            c = std::to_string(cnt);
 
            // use relaxation
            if (GetVWIIterationType() !=
                eFixedWaveForcingWithSubIterationOnAlpha)
            {
                // the critical layer should be the bnd region 3
                // int reg =3;
                // FileRelaxation(reg);
            }
            // calculate the jump conditions
            string wavefile = m_sessionName + ".fld";
            syscall = "../../utilities/PostProcessing/Extras/FldCalcBCs  " +
                      movedmesh + "  " + wavefile + "  " + interpstreak +
                      ">  data" + c;
            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 + "   " +
                      c_alpha + "      " + 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 + "    " +
                      c_alpha + "      " + 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);
            }
 
            // 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" + c;
            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] = "";
            snprintf(alpchar, 16, "%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();
        }
    }
}

References ASSERTL0, CalcNonLinearWaveForce(), Nektar::eFixedWaveForcingWithSubIterationOnAlpha, ExecuteRoll(), ExecuteStreak(), ExecuteWave(), GetVWIIterationType(), m_alpha, m_leading_imag_evl, m_moveMeshToCriticalLayer, m_sessionName, m_sessionRoll, and m_sessionVWI.

Referenced by DoFixedForcingIteration(), and main().

ExecuteRoll()

void Nektar::VortexWaveInteraction::ExecuteRoll

(

void

)

Definition at line 368 of file VortexWaveInteraction.cpp.

{
    // 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,
                                                      m_graphRoll);
        // 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, m_graphRoll);
        }
        else
        {
            const int npoints = m_solverRoll->GetNpoints();
            const int ncoeffs = m_solverRoll->GetNcoeffs();
 
            static int init = 1;
            if (init)
            {
                m_solverRoll->DoInitialise();
                m_vwiForcingObj =
                    std::dynamic_pointer_cast<SolverUtils::ForcingProgrammatic>(
                        GetForcingFactory().CreateInstance(
                            "Programmatic", m_sessionRoll, m_solverRoll,
                            m_solverRoll->UpdateFields(),
                            m_solverRoll->UpdateFields().size() - 1, 0));
 
                std::vector<std::string> vFieldNames =
                    m_sessionRoll->GetVariables();
                vFieldNames.erase(vFieldNames.end() - 1);
 
                m_solverRoll->GetFunction("BodyForce")
                    ->Evaluate(vFieldNames, m_vwiForcingObj->UpdateForces());
 
                // Scale forcing
                for (int i = 0; i < m_vwiForcingObj->UpdateForces().size(); ++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().size(); ++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);
}

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

Referenced by ExecuteLoop().

ExecuteStreak()

void Nektar::VortexWaveInteraction::ExecuteStreak

(

void

)

Definition at line 495 of file VortexWaveInteraction.cpp.

{
    // Create driver
#if 1
    std::string vDriverModule;
    m_sessionStreak->LoadSolverInfo("Driver", vDriverModule, "Standard");
 
    DriverSharedPtr solverStreak = GetDriverFactory().CreateInstance(
        vDriverModule, m_sessionStreak, m_graphStreak);
    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,
                                                  m_graphStreak);
 
    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");
}

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

Referenced by ExecuteLoop().

ExecuteWave()

void Nektar::VortexWaveInteraction::ExecuteWave

(

void

)

Do linearised NavierStokes Session with Modified Arnoldi

Definition at line 539 of file VortexWaveInteraction.cpp.

{
 
    // 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 Solver" << endl;
    std::shared_ptr<DriverArnoldi> solverWave =
        std::dynamic_pointer_cast<SolverUtils::DriverArnoldi>(
            GetDriverFactory().CreateInstance(vDriverModule, m_sessionWave,
                                              m_graphWave));
 
    /// 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;
    }
}

References CopyFile(), Nektar::SolverUtils::GetDriverFactory(), m_alpha, m_graphWave, m_leading_imag_evl, m_leading_real_evl, m_sessionVWI, m_sessionWave, m_wavePressure, and m_waveVelocities.

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

FileRelaxation()

void Nektar::VortexWaveInteraction::FileRelaxation

(

int

reg

)

Definition at line 2124 of file VortexWaveInteraction.cpp.

{
    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::ExpListSharedPtr Ilayer;
    Ilayer = MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(
        *std::static_pointer_cast<MultiRegions::ExpList>(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
 
    std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef_u;
    std::vector<std::vector<NekDouble>> FieldData_u;
    string file = m_sessionName;
 
    file += "_u_5.bc";
    LibUtilities::FieldIOSharedPtr fld =
        LibUtilities::FieldIO::CreateDefault(m_sessionVWI);
    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(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->FwdTransLocalElmt(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(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->FwdTransLocalElmt(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++;
}

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::LibUtilities::FieldIO::CreateDefault(), m_bcsForcing, m_rollField, m_sessionName, m_sessionVWI, m_vwiRelaxation, Vmath::Smul(), Vmath::Svtvp(), and Vmath::Vcopy().

GetAlpha()

NekDouble Nektar::VortexWaveInteraction::GetAlpha ( void  )

inline

Definition at line 121 of file VortexWaveInteraction.h.

    {
        return m_alpha[0];
    }

References m_alpha.

Referenced by DoFixedForcingIteration(), and main().

GetAlphaStep()

NekDouble Nektar::VortexWaveInteraction::GetAlphaStep ( void  )

inline

Definition at line 126 of file VortexWaveInteraction.h.

    {
        return m_alphaStep;
    }

References m_alphaStep.

Referenced by main().

GetDAlphaDWaveForceMag()

NekDouble Nektar::VortexWaveInteraction::GetDAlphaDWaveForceMag ( void  )

inline

Definition at line 141 of file VortexWaveInteraction.h.

    {
        return m_dAlphaDWaveForceMag;
    }

References m_dAlphaDWaveForceMag.

GetEigRelTol()

NekDouble Nektar::VortexWaveInteraction::GetEigRelTol ( void  )

inline

Definition at line 151 of file VortexWaveInteraction.h.

    {
        return m_eigRelTol;
    }

References m_eigRelTol.

Referenced by DoFixedForcingIteration().

GetIterEnd()

int Nektar::VortexWaveInteraction::GetIterEnd ( )

inline

Definition at line 101 of file VortexWaveInteraction.h.

    {
        return m_iterEnd;
    }

References m_iterEnd.

Referenced by DoFixedForcingIteration().

GetIterStart()

int Nektar::VortexWaveInteraction::GetIterStart ( )

inline

Definition at line 96 of file VortexWaveInteraction.h.

    {
        return m_iterStart;
    }

References m_iterStart.

Referenced by DoFixedForcingIteration().

GetMaxOuterIterations()

int Nektar::VortexWaveInteraction::GetMaxOuterIterations ( void  )

inline

Definition at line 116 of file VortexWaveInteraction.h.

    {
        return m_maxOuterIterations;
    }

References m_maxOuterIterations.

Referenced by DoFixedForcingIteration().

GetMaxWaveForceMagIter()

int Nektar::VortexWaveInteraction::GetMaxWaveForceMagIter ( void  )

inline

Definition at line 146 of file VortexWaveInteraction.h.

    {
        return m_maxWaveForceMagIter;
    }

References m_maxWaveForceMagIter.

Referenced by main().

GetMinInnerIterations()

int Nektar::VortexWaveInteraction::GetMinInnerIterations ( void  )

inline

Definition at line 157 of file VortexWaveInteraction.h.

    {
        return m_minInnerIterations;
    }

References m_minInnerIterations.

Referenced by DoFixedForcingIteration().

GetNOuterIterations()

int Nektar::VortexWaveInteraction::GetNOuterIterations ( void  )

inline

Definition at line 111 of file VortexWaveInteraction.h.

    {
        return m_nOuterIterations;
    }

References m_nOuterIterations.

Referenced by DoFixedForcingIteration().

GetPrevAlpha()

NekDouble Nektar::VortexWaveInteraction::GetPrevAlpha ( void  )

inline

Definition at line 162 of file VortexWaveInteraction.h.

    {
        return m_prevAlpha;
    }

References m_prevAlpha.

Referenced by DoFixedForcingIteration().

GetReflectionIndex()

Array< OneD, int > Nektar::VortexWaveInteraction::GetReflectionIndex

(

void

)

Definition at line 2041 of file VortexWaveInteraction.cpp.

{
    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  = 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;
}

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

Referenced by CalcNonLinearWaveForce().

GetVWIIterationType()

VWIIterationType Nektar::VortexWaveInteraction::GetVWIIterationType ( void  )

inline

Definition at line 106 of file VortexWaveInteraction.h.

    {
        return m_VWIIterationType;
    }

References m_VWIIterationType.

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

GetWaveForceMag()

NekDouble Nektar::VortexWaveInteraction::GetWaveForceMag ( void  )

inline

Definition at line 131 of file VortexWaveInteraction.h.

    {
        return m_waveForceMag[0];
    }

References m_waveForceMag.

Referenced by main().

GetWaveForceMagStep()

NekDouble Nektar::VortexWaveInteraction::GetWaveForceMagStep ( void  )

inline

Definition at line 136 of file VortexWaveInteraction.h.

    {
        return m_waveForceMagStep;
    }

References m_waveForceMagStep.

Referenced by main().

IfIterInterface()

bool Nektar::VortexWaveInteraction::IfIterInterface ( void  )

inline

Definition at line 197 of file VortexWaveInteraction.h.

    {
        return m_iterinterface;
    }

References m_iterinterface.

Referenced by DoFixedForcingIteration().

MoveFile()

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

Definition at line 1131 of file VortexWaveInteraction.cpp.

{
    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;
        ASSERTL0(system(mkdir.c_str()) == 0,
                 "Failed to make directory '" + dir + "'");
        opendir[dir] = 1;
    }
 
    string savefile = dir + "/" + file + "." + std::to_string(n);
    string syscall  = "mv -f " + file + " " + savefile;
 
    ASSERTL0(system(syscall.c_str()) == 0, syscall.c_str());
}

References ASSERTL0, and CG_Iterations::savefile.

Referenced by DoFixedForcingIteration().

SaveFile()

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

Definition at line 1111 of file VortexWaveInteraction.cpp.

{
    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;
        ASSERTL0(system(mkdir.c_str()) == 0,
                 "Failed to make directory '" + dir + "'");
 
        opendir[dir] = 1;
    }
 
    string savefile = dir + "/" + file + "." + std::to_string(n);
    string syscall  = "cp -f " + file + " " + savefile;
 
    ASSERTL0(system(syscall.c_str()) == 0, syscall.c_str());
}

References ASSERTL0, and CG_Iterations::savefile.

Referenced by SaveLoopDetails().

SaveLoopDetails()

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

Definition at line 1180 of file VortexWaveInteraction.cpp.

{
    // 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);
    }
}

References m_sessionName, m_sessionVWI, and SaveFile().

Referenced by DoFixedForcingIteration().

SetAlpha()

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

inline

Definition at line 167 of file VortexWaveInteraction.h.

    {
        m_alpha[0] = alpha;
    }

References m_alpha.

Referenced by main().

SetAlphaStep()

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

inline

Definition at line 182 of file VortexWaveInteraction.h.

    {
        m_alphaStep = step;
    }

References m_alphaStep.

SetEigRelTol()

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

inline

Definition at line 177 of file VortexWaveInteraction.h.

    {
        m_eigRelTol = tol;
    }

References m_eigRelTol.

Referenced by DoFixedForcingIteration().

SetMinInnerIterations()

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

inline

Definition at line 187 of file VortexWaveInteraction.h.

    {
        m_minInnerIterations = niter;
    }

References m_minInnerIterations.

Referenced by DoFixedForcingIteration().

SetPrevAlpha()

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

inline

Definition at line 192 of file VortexWaveInteraction.h.

    {
        m_prevAlpha = alpha;
    }

References m_prevAlpha.

Referenced by main().

SetWaveForceMag()

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

inline

Definition at line 172 of file VortexWaveInteraction.h.

    {
        m_waveForceMag[0] = mag;
    }

References m_waveForceMag.

Referenced by main().

UpdateAlpha()

void Nektar::VortexWaveInteraction::UpdateAlpha

(

int

n

)

Definition at line 1926 of file VortexWaveInteraction.cpp.

{
    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.size() < outeriter) ? m_alpha.size() : 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.size() - 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;
}

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

Referenced by DoFixedForcingIteration(), and VortexWaveInteraction().

UpdateDAlphaDWaveForceMag()

void Nektar::VortexWaveInteraction::UpdateDAlphaDWaveForceMag

(

NekDouble

alphainit

)

Definition at line 1921 of file VortexWaveInteraction.cpp.

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

References m_alpha, m_dAlphaDWaveForceMag, and m_waveForceMagStep.

Referenced by DoFixedForcingIteration().

UpdateWaveForceMag()

void Nektar::VortexWaveInteraction::UpdateWaveForceMag

(

int

n

)

Definition at line 1800 of file VortexWaveInteraction.cpp.

{
    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.size() < outeriter) ? m_waveForceMag.size()
                                                         : 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.size() - 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;
}

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

Referenced by DoFixedForcingIteration().

Member Data Documentation

m_alpha

Array<OneD, NekDouble> Nektar::VortexWaveInteraction::m_alpha

private

< Leading imaginary eigenvalue

Definition at line 233 of file VortexWaveInteraction.h.

Referenced by AppendEvlToFile(), CalcNonLinearWaveForce(), ExecuteLoop(), ExecuteWave(), GetAlpha(), SetAlpha(), UpdateAlpha(), UpdateDAlphaDWaveForceMag(), and VortexWaveInteraction().

m_alphaStep

NekDouble Nektar::VortexWaveInteraction::m_alphaStep

private

Definition at line 235 of file VortexWaveInteraction.h.

Referenced by GetAlphaStep(), SetAlphaStep(), UpdateAlpha(), and VortexWaveInteraction().

m_bcsForcing

Array<OneD, Array<OneD, NekDouble> > Nektar::VortexWaveInteraction::m_bcsForcing

private

Definition at line 252 of file VortexWaveInteraction.h.

Referenced by FileRelaxation().

m_dAlphaDWaveForceMag

NekDouble Nektar::VortexWaveInteraction::m_dAlphaDWaveForceMag

private

Definition at line 239 of file VortexWaveInteraction.h.

Referenced by GetDAlphaDWaveForceMag(), UpdateDAlphaDWaveForceMag(), and VortexWaveInteraction().

m_deltaFcnApprox

bool Nektar::VortexWaveInteraction::m_deltaFcnApprox

private

Definition at line 217 of file VortexWaveInteraction.h.

Referenced by CalcNonLinearWaveForce(), and VortexWaveInteraction().

m_deltaFcnDecay

NekDouble Nektar::VortexWaveInteraction::m_deltaFcnDecay

private

Definition at line 222 of file VortexWaveInteraction.h.

Referenced by CalcNonLinearWaveForce(), and VortexWaveInteraction().

m_eigRelTol

NekDouble Nektar::VortexWaveInteraction::m_eigRelTol

private

Definition at line 237 of file VortexWaveInteraction.h.

Referenced by CheckEigIsStationary(), GetEigRelTol(), SetEigRelTol(), and VortexWaveInteraction().

m_graphRoll

SpatialDomains::MeshGraphSharedPtr Nektar::VortexWaveInteraction::m_graphRoll

private

Definition at line 262 of file VortexWaveInteraction.h.

Referenced by ExecuteRoll(), and VortexWaveInteraction().

m_graphStreak

SpatialDomains::MeshGraphSharedPtr Nektar::VortexWaveInteraction::m_graphStreak

private

Definition at line 266 of file VortexWaveInteraction.h.

Referenced by ExecuteStreak(), and VortexWaveInteraction().

m_graphWave

SpatialDomains::MeshGraphSharedPtr Nektar::VortexWaveInteraction::m_graphWave

private

Definition at line 268 of file VortexWaveInteraction.h.

Referenced by ExecuteWave(), and VortexWaveInteraction().

m_iterEnd

int Nektar::VortexWaveInteraction::m_iterEnd

private

Definition at line 209 of file VortexWaveInteraction.h.

Referenced by GetIterEnd(), and VortexWaveInteraction().

m_iterinterface

bool Nektar::VortexWaveInteraction::m_iterinterface

private

Definition at line 242 of file VortexWaveInteraction.h.

Referenced by IfIterInterface(), and VortexWaveInteraction().

m_iterStart

int Nektar::VortexWaveInteraction::m_iterStart

private

Definition at line 208 of file VortexWaveInteraction.h.

Referenced by GetIterStart(), and VortexWaveInteraction().

m_leading_imag_evl

Array<OneD, NekDouble> Nektar::VortexWaveInteraction::m_leading_imag_evl

private

< Leading real eigenvalue

Definition at line 231 of file VortexWaveInteraction.h.

Referenced by AppendEvlToFile(), CheckEigIsStationary(), CheckIfAtNeutralPoint(), ExecuteLoop(), ExecuteWave(), UpdateAlpha(), UpdateWaveForceMag(), and VortexWaveInteraction().

m_leading_real_evl

Array<OneD, NekDouble> Nektar::VortexWaveInteraction::m_leading_real_evl

private

Definition at line 229 of file VortexWaveInteraction.h.

Referenced by AppendEvlToFile(), CheckEigIsStationary(), CheckIfAtNeutralPoint(), ExecuteWave(), UpdateAlpha(), UpdateWaveForceMag(), and VortexWaveInteraction().

m_maxOuterIterations

int Nektar::VortexWaveInteraction::m_maxOuterIterations

private

Definition at line 212 of file VortexWaveInteraction.h.

Referenced by GetMaxOuterIterations(), and VortexWaveInteraction().

m_maxWaveForceMagIter

int Nektar::VortexWaveInteraction::m_maxWaveForceMagIter

private

Definition at line 215 of file VortexWaveInteraction.h.

Referenced by GetMaxWaveForceMagIter(), and VortexWaveInteraction().

m_minInnerIterations

int Nektar::VortexWaveInteraction::m_minInnerIterations

private

Definition at line 213 of file VortexWaveInteraction.h.

Referenced by CheckEigIsStationary(), GetMinInnerIterations(), SetMinInnerIterations(), and VortexWaveInteraction().

m_moveMeshToCriticalLayer

bool Nektar::VortexWaveInteraction::m_moveMeshToCriticalLayer

private

Definition at line 220 of file VortexWaveInteraction.h.

Referenced by ExecuteLoop(), ExecuteRoll(), and VortexWaveInteraction().

m_neutralPointTol

NekDouble Nektar::VortexWaveInteraction::m_neutralPointTol

private

Definition at line 236 of file VortexWaveInteraction.h.

Referenced by CheckIfAtNeutralPoint(), and VortexWaveInteraction().

m_nOuterIterations

int Nektar::VortexWaveInteraction::m_nOuterIterations

private

Definition at line 211 of file VortexWaveInteraction.h.

Referenced by GetNOuterIterations(), and VortexWaveInteraction().

m_prevAlpha

NekDouble Nektar::VortexWaveInteraction::m_prevAlpha

private

Definition at line 240 of file VortexWaveInteraction.h.

Referenced by GetPrevAlpha(), and SetPrevAlpha().

m_rollField

Array<OneD, MultiRegions::ExpListSharedPtr> Nektar::VortexWaveInteraction::m_rollField

private

Definition at line 256 of file VortexWaveInteraction.h.

Referenced by ExecuteRoll(), and FileRelaxation().

m_rollForceScale

NekDouble Nektar::VortexWaveInteraction::m_rollForceScale

private

Definition at line 227 of file VortexWaveInteraction.h.

Referenced by ExecuteRoll(), and VortexWaveInteraction().

m_sessionName

std::string Nektar::VortexWaveInteraction::m_sessionName

private

Definition at line 258 of file VortexWaveInteraction.h.

Referenced by CalcNonLinearWaveForce(), CopyFile(), ExecuteLoop(), ExecuteRoll(), FileRelaxation(), SaveLoopDetails(), and VortexWaveInteraction().

m_sessionRoll

LibUtilities::SessionReaderSharedPtr Nektar::VortexWaveInteraction::m_sessionRoll

private

Definition at line 261 of file VortexWaveInteraction.h.

Referenced by CalcNonLinearWaveForce(), CheckIfAtNeutralPoint(), ExecuteLoop(), ExecuteRoll(), and VortexWaveInteraction().

m_sessionStreak

LibUtilities::SessionReaderSharedPtr Nektar::VortexWaveInteraction::m_sessionStreak

private

Definition at line 265 of file VortexWaveInteraction.h.

Referenced by ExecuteStreak(), and VortexWaveInteraction().

m_sessionVWI

LibUtilities::SessionReaderSharedPtr Nektar::VortexWaveInteraction::m_sessionVWI

private

Definition at line 259 of file VortexWaveInteraction.h.

Referenced by CalcNonLinearWaveForce(), CheckIfAtNeutralPoint(), ExecuteLoop(), ExecuteStreak(), ExecuteWave(), FileRelaxation(), GetReflectionIndex(), SaveLoopDetails(), VortexWaveInteraction(), and ~VortexWaveInteraction().

m_sessionWave

LibUtilities::SessionReaderSharedPtr Nektar::VortexWaveInteraction::m_sessionWave

private

Definition at line 267 of file VortexWaveInteraction.h.

Referenced by ExecuteWave(), and VortexWaveInteraction().

m_solverRoll

EquationSystemSharedPtr Nektar::VortexWaveInteraction::m_solverRoll

private

Definition at line 263 of file VortexWaveInteraction.h.

Referenced by CalcNonLinearWaveForce(), ExecuteRoll(), and VortexWaveInteraction().

m_streakField

Array<OneD, MultiRegions::ExpListSharedPtr> Nektar::VortexWaveInteraction::m_streakField

private

Definition at line 254 of file VortexWaveInteraction.h.

Referenced by CalcNonLinearWaveForce(), and ExecuteStreak().

m_useLinfPressureNorm

bool Nektar::VortexWaveInteraction::m_useLinfPressureNorm

private

Definition at line 218 of file VortexWaveInteraction.h.

Referenced by CalcNonLinearWaveForce(), and VortexWaveInteraction().

m_vwiForcing

Array<OneD, Array<OneD, NekDouble> > Nektar::VortexWaveInteraction::m_vwiForcing

private

Definition at line 249 of file VortexWaveInteraction.h.

Referenced by CalcNonLinearWaveForce(), ExecuteRoll(), and VortexWaveInteraction().

m_vwiForcingObj

SolverUtils::ForcingProgrammaticSharedPtr Nektar::VortexWaveInteraction::m_vwiForcingObj

private

Definition at line 250 of file VortexWaveInteraction.h.

Referenced by ExecuteRoll().

m_VWIIterationType

VWIIterationType Nektar::VortexWaveInteraction::m_VWIIterationType

private

Definition at line 244 of file VortexWaveInteraction.h.

Referenced by GetVWIIterationType(), and VortexWaveInteraction().

m_vwiRelaxation

NekDouble Nektar::VortexWaveInteraction::m_vwiRelaxation

private

Definition at line 238 of file VortexWaveInteraction.h.

Referenced by CalcNonLinearWaveForce(), FileRelaxation(), and VortexWaveInteraction().

m_waveForceMag

Array<OneD, NekDouble> Nektar::VortexWaveInteraction::m_waveForceMag

private

Definition at line 224 of file VortexWaveInteraction.h.

Referenced by CalcNonLinearWaveForce(), GetWaveForceMag(), SetWaveForceMag(), UpdateWaveForceMag(), and VortexWaveInteraction().

m_waveForceMagStep

NekDouble Nektar::VortexWaveInteraction::m_waveForceMagStep

private

Definition at line 225 of file VortexWaveInteraction.h.

Referenced by GetWaveForceMagStep(), UpdateDAlphaDWaveForceMag(), UpdateWaveForceMag(), and VortexWaveInteraction().

m_wavePressure

MultiRegions::ExpListSharedPtr Nektar::VortexWaveInteraction::m_wavePressure

private

Definition at line 247 of file VortexWaveInteraction.h.

Referenced by CalcL2ToLinfPressure(), CalcNonLinearWaveForce(), and ExecuteWave().

m_waveVelocities

Array<OneD, MultiRegions::ExpListSharedPtr> Nektar::VortexWaveInteraction::m_waveVelocities

private

Definition at line 246 of file VortexWaveInteraction.h.

Referenced by CalcL2ToLinfPressure(), CalcNonLinearWaveForce(), ExecuteWave(), and GetReflectionIndex().