#include <FilterHilbertFFTPhase.h>

Inheritance diagram for Nektar::FilterHilbertFFTPhase:

Public Member Functions
	FilterHilbertFFTPhase (const LibUtilities::SessionReaderSharedPtr &pSession, const std::weak_ptr< SolverUtils::EquationSystem > &pEquation, const ParamMap &pParams)

	~FilterHilbertFFTPhase () override

Public Member Functions inherited from Nektar::SolverUtils::Filter
SOLVER_UTILS_EXPORT	Filter (const LibUtilities::SessionReaderSharedPtr &pSession, const std::weak_ptr< EquationSystem > &pEquation)

virtual SOLVER_UTILS_EXPORT	~Filter ()

SOLVER_UTILS_EXPORT void	Initialise (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)

SOLVER_UTILS_EXPORT void	Update (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)

SOLVER_UTILS_EXPORT void	Finalise (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)

SOLVER_UTILS_EXPORT bool	IsTimeDependent ()

Static Public Member Functions
static SolverUtils::FilterSharedPtr	create (const LibUtilities::SessionReaderSharedPtr &pSession, const std::weak_ptr< SolverUtils::EquationSystem > &pEquation, const ParamMap &pParams)
	Creates an instance of this class. More...

Static Public Attributes
static std::string	className
	Name of the class. More...

Protected Member Functions
void	v_Initialise (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time) override

void	v_Update (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time) override

void	v_Finalise (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time) override

bool	v_IsTimeDependent () override

virtual void	v_Initialise (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)=0

virtual void	v_Update (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)=0

virtual void	v_Finalise (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)=0

virtual bool	v_IsTimeDependent ()=0

Private Attributes
unsigned int	m_index

unsigned int	m_outputIndex

unsigned int	vCounter

unsigned int	m_outputFrequency

int	m_window

unsigned int	m_overlap

bool	linear_trans = true

bool	provided_mean = false

std::string	m_outputFile

LibUtilities::FieldIOSharedPtr	m_fld

unsigned int	npoints

NekDouble	m_TimeStep

double	m_mean

Array< OneD, NekDouble >	oldV

Array< OneD, NekDouble >	oldV_zero_mean

Array< OneD, NekDouble >	fcoef

Array< OneD, NekDouble >	h

Array< OneD, NekDouble >	overlap_phase

fftw_plan	plan_forward

fftw_plan	plan_backward

Friends
class	MemoryManager< FilterHilbertFFTPhase >

Additional Inherited Members
Public Types inherited from Nektar::SolverUtils::Filter
typedef std::map< std::string, std::string >	ParamMap

Protected Attributes inherited from Nektar::SolverUtils::Filter
LibUtilities::SessionReaderSharedPtr	m_session

const std::weak_ptr< EquationSystem >	m_equ

Detailed Description

Definition at line 45 of file FilterHilbertFFTPhase.h.

Constructor & Destructor Documentation

◆ FilterHilbertFFTPhase()

Nektar::FilterHilbertFFTPhase::FilterHilbertFFTPhase	(	const LibUtilities::SessionReaderSharedPtr &	pSession,
		const std::weak_ptr< SolverUtils::EquationSystem > &	pEquation,
		const ParamMap &	pParams
	)

Definition at line 44 of file FilterHilbertFFTPhase.cpp.

    : Filter(pSession, pEquation)
{
    if (pParams.find("OutputFile") == pParams.end())
    {
        m_outputFile = m_session->GetSessionName();
    }
    else
    {
        ASSERTL0(!(pParams.find("OutputFile")->second.empty()),
                 "Missing parameter 'OutputFile'.");
        m_outputFile = pParams.find("OutputFile")->second;
    }
 
    ASSERTL0(pParams.find("OutputFrequency") != pParams.end(),
             "Missing parameter 'OutputFrequency'."); // Sampling frequency =
                                                      // per number of timesteps
    LibUtilities::Equation equ(m_session->GetInterpreter(),
                               pParams.find("OutputFrequency")->second);
    m_outputFrequency = floor(equ.Evaluate());
 
    ASSERTL0(pParams.find("WindowSize") != pParams.end(),
             "Missing parameter 'WindowSize'.");
    LibUtilities::Equation equ1(m_session->GetInterpreter(),
                                pParams.find("WindowSize")->second);
    m_window = floor(equ1.Evaluate()); // No. of output steps taken in each FFT
    ASSERTL0(m_window % 2 == 0, "'WindowSize' must be even.");
 
    ASSERTL0(pParams.find("OverlapSize") != pParams.end(),
             "Missing parameter 'OverlapSize'.");
    LibUtilities::Equation equ2(m_session->GetInterpreter(),
                                pParams.find("OverlapSize")->second);
    m_overlap = floor(equ2.Evaluate());
    ASSERTL0(m_window > m_overlap,
             "'OverlapSize' must be smaller than 'WindowSize'.");
 
    if (pParams.find("MeanV") != pParams.end())
    {
        provided_mean = true;
        LibUtilities::Equation equ3(m_session->GetInterpreter(),
                                    pParams.find("MeanV")->second);
        m_mean = equ3.Evaluate();
    }
 
    if (pParams.find("LinearTransitionOverlap") != pParams.end())
    {
        linear_trans =
            (pParams.find("LinearTransitionOverlap")->second == "true");
    }
 
    m_outputIndex = 0;
    m_index       = 0;
 
    m_fld = LibUtilities::FieldIO::CreateDefault(pSession);
}

References ASSERTL0, Nektar::LibUtilities::FieldIO::CreateDefault(), Nektar::LibUtilities::Equation::Evaluate(), linear_trans, m_fld, m_index, m_mean, m_outputFile, m_outputFrequency, m_outputIndex, m_overlap, Nektar::SolverUtils::Filter::m_session, m_window, and provided_mean.

◆ ~FilterHilbertFFTPhase()

Nektar::FilterHilbertFFTPhase::~FilterHilbertFFTPhase ( )

override

Definition at line 103 of file FilterHilbertFFTPhase.cpp.

104{

105}

Member Function Documentation

◆ create()

static SolverUtils::FilterSharedPtr Nektar::FilterHilbertFFTPhase::create	(	const LibUtilities::SessionReaderSharedPtr &	pSession,
		const std::weak_ptr< SolverUtils::EquationSystem > &	pEquation,
		const ParamMap &	pParams
	)

inlinestatic

Creates an instance of this class.

Definition at line 51 of file FilterHilbertFFTPhase.h.

    {
        SolverUtils::FilterSharedPtr p =
            MemoryManager<FilterHilbertFFTPhase>::AllocateSharedPtr(
                pSession, pEquation, pParams);
        return p;
    }

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), and CellMLToNektar.cellml_metadata::p.

◆ v_Finalise()

void Nektar::FilterHilbertFFTPhase::v_Finalise	(	const Array< OneD, const MultiRegions::ExpListSharedPtr > &	pFields,
		const NekDouble &	time
	)

overrideprotectedvirtual

Implements Nektar::SolverUtils::Filter.

Definition at line 293 of file FilterHilbertFFTPhase.cpp.

297{

298}

◆ v_Initialise()

void Nektar::FilterHilbertFFTPhase::v_Initialise	(	const Array< OneD, const MultiRegions::ExpListSharedPtr > &	pFields,
		const NekDouble &	time
	)

overrideprotectedvirtual

Implements Nektar::SolverUtils::Filter.

Definition at line 107 of file FilterHilbertFFTPhase.cpp.

{
    m_index       = 0;
    m_outputIndex = 0;
    vCounter      = 0;
 
    m_TimeStep = pFields[0]->GetSession()->GetParameter("TimeStep");
 
    npoints        = pFields[0]->GetPhys().size();
    oldV           = Array<OneD, NekDouble>(npoints * m_window);
    oldV_zero_mean = Array<OneD, NekDouble>(npoints * m_window);
    fcoef          = Array<OneD, NekDouble>(npoints * m_window);
    h              = Array<OneD, NekDouble>(npoints * m_window);
    overlap_phase  = Array<OneD, NekDouble>(npoints * m_overlap);
    std::fill_n(&overlap_phase[0], npoints * m_overlap, 0.0);
 
    // Initialise FFTW
    fftw_r2r_kind *fwd_kind = new fftw_r2r_kind[npoints];
    std::fill_n(fwd_kind, npoints, FFTW_R2HC);
    plan_forward = fftw_plan_many_r2r(
        1, &m_window, npoints, &oldV_zero_mean[0], &m_window, npoints, 1,
        &fcoef[0], &m_window, npoints, 1, &fwd_kind[0], FFTW_ESTIMATE);
 
    fftw_r2r_kind *bwd_kind = new fftw_r2r_kind[npoints];
    std::fill_n(bwd_kind, npoints, FFTW_HC2R);
    plan_backward = fftw_plan_many_r2r(1, &m_window, npoints, &fcoef[0],
                                       &m_window, npoints, 1, &h[0], &m_window,
                                       npoints, 1, &bwd_kind[0], FFTW_ESTIMATE);
 
    delete[] fwd_kind;
    delete[] bwd_kind;
    v_Update(pFields, 0.0);
}

References fcoef, h, m_index, m_outputIndex, m_overlap, m_TimeStep, m_window, npoints, oldV, oldV_zero_mean, overlap_phase, plan_backward, plan_forward, v_Update(), and vCounter.

◆ v_IsTimeDependent()

bool Nektar::FilterHilbertFFTPhase::v_IsTimeDependent ( )

overrideprotectedvirtual

Implements Nektar::SolverUtils::Filter.

Definition at line 300 of file FilterHilbertFFTPhase.cpp.

{
    return true;
}

◆ v_Update()

void Nektar::FilterHilbertFFTPhase::v_Update	(	const Array< OneD, const MultiRegions::ExpListSharedPtr > &	pFields,
		const NekDouble &	time
	)

overrideprotectedvirtual

Implements Nektar::SolverUtils::Filter.

Definition at line 143 of file FilterHilbertFFTPhase.cpp.

{
    if (m_index++ % m_outputFrequency > 0)
    {
        return;
    }
 
    // Store old data
    Array<OneD, NekDouble> currentV = pFields[0]->GetPhys();
    memcpy(&oldV[vCounter * npoints], &currentV[0],
           npoints * sizeof(NekDouble));
    vCounter++;
 
    // Do Hilbert Transform and output phase
    if (vCounter == m_window)
    {
        // Zero mean
        if (!provided_mean)
        {
            for (int j = 0; j < npoints; ++j)
            {
                m_mean = 0.0;
                for (int i = 0; i < m_window; ++i)
                {
                    m_mean += oldV[i * npoints + j];
                }
                m_mean = m_mean / (double)m_window; // time-average
 
                for (int i = 0; i < m_window; ++i)
                {
                    oldV_zero_mean[i * npoints + j] =
                        oldV[i * npoints + j] - m_mean;
                }
            }
        }
        else // minus mean directly
        {
            for (int i = 0; i < npoints * m_window; ++i)
            {
                oldV_zero_mean[i] = oldV[i] - m_mean;
            }
        }
 
        // Do Hilbert Transform:
        fftw_execute(plan_forward);
        // Rearrange out: imag->real, -real->imag
        double tmp_fcoef;
        for (int i = 1; i < m_window / 2; ++i)
        {
            for (int j = 0; j < npoints; ++j)
            {
                tmp_fcoef = -fcoef[i * npoints + j] / m_window;
                fcoef[i * npoints + j] =
                    fcoef[(m_window - i) * npoints + j] / m_window;
                fcoef[(m_window - i) * npoints + j] = tmp_fcoef;
            }
        }
        std::fill_n(&fcoef[0], npoints, 0.0);
        std::fill_n(&fcoef[(m_window / 2) * npoints], npoints, 0.0);
        fftw_execute(plan_backward);
 
        // Output batch
        Array<OneD, NekDouble> phase_coeff(pFields[0]->GetNcoeffs());
        Array<OneD, NekDouble> phase(npoints);
 
        for (int t = 0; t < m_window - m_overlap; ++t)
        {
            std::stringstream vOutputFilename;
            vOutputFilename << m_outputFile << "_" << m_outputIndex << ".chk";
 
            // Get FieldDefinitions and allocate Fielddata
            std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef =
                pFields[0]->GetFieldDefinitions();
            std::vector<std::vector<NekDouble>> FieldData(FieldDef.size());
 
            // Calculate phase
            if (t < m_overlap)
            {
                if (linear_trans)
                {
                    for (int i = 0; i < npoints; ++i)
                    {
                        phase[i] = ((double)(m_overlap - t) / m_overlap) *
                                       overlap_phase[t * npoints + i] +
                                   ((double)t / m_overlap) *
                                       atan2(h[t * npoints + i],
                                             oldV_zero_mean[t * npoints + i]);
                    }
                }
                else // simple average
                {
                    for (int i = 0; i < npoints; ++i)
                    {
                        phase[i] = 0.5 * overlap_phase[t * npoints + i] +
                                   0.5 * atan2(h[t * npoints + i],
                                               oldV_zero_mean[t * npoints + i]);
                    }
                }
            }
            else
            {
                for (int i = 0; i < npoints; ++i)
                {
                    phase[i] = atan2(h[t * npoints + i],
                                     oldV_zero_mean[t * npoints + i]);
                }
            }
 
            // Convert to coefficient space
            pFields[0]->FwdTransLocalElmt(phase, phase_coeff);
 
            // for each item in FieldDef
            //      push "phase" to FieldDef[i]->m_fields
            //      call AppendFieldData(FieldDef[i], FieldData[i], data)
            for (int i = 0; i < FieldDef.size(); ++i)
            {
                FieldDef[i]->m_fields.push_back("phase");
                pFields[0]->AppendFieldData(FieldDef[i], FieldData[i],
                                            phase_coeff);
            }
 
            // Update time in field info if required
            LibUtilities::FieldMetaDataMap fieldMetaDataMap;
            fieldMetaDataMap["Time"] = boost::lexical_cast<std::string>(
                m_outputIndex * m_outputFrequency * m_TimeStep);
 
            m_fld->Write(vOutputFilename.str(), FieldDef, FieldData,
                         fieldMetaDataMap);
            m_outputIndex++;
        }
 
        // Save phase for next iter
        for (int t = m_window - m_overlap; t < m_window; ++t)
        {
            for (int i = 0; i < npoints; ++i)
            {
                overlap_phase[(t - m_window + m_overlap) * npoints + i] =
                    atan2(h[t * npoints + i], oldV_zero_mean[t * npoints + i]);
            }
        }
 
        // Move overlap part to front
        memcpy(&oldV[0], &oldV[(m_window - m_overlap) * npoints],
               npoints * m_overlap * sizeof(NekDouble));
        vCounter = m_overlap;
    }
}

References fcoef, h, linear_trans, m_fld, m_index, m_mean, m_outputFile, m_outputFrequency, m_outputIndex, m_overlap, m_TimeStep, m_window, npoints, oldV, oldV_zero_mean, overlap_phase, plan_backward, plan_forward, provided_mean, and vCounter.

Referenced by v_Initialise().

Friends And Related Function Documentation

◆ MemoryManager< FilterHilbertFFTPhase >

friend class MemoryManager< FilterHilbertFFTPhase >

friend

Definition at line 1 of file FilterHilbertFFTPhase.h.

Member Data Documentation

◆ className

std::string Nektar::FilterHilbertFFTPhase::className

static

Initial value:

=
    SolverUtils::GetFilterFactory().RegisterCreatorFunction(
        "HilbertPhase", FilterHilbertFFTPhase::create)

Name of the class.

Definition at line 63 of file FilterHilbertFFTPhase.h.

◆ fcoef

Array<OneD, NekDouble> Nektar::FilterHilbertFFTPhase::fcoef

private

Definition at line 106 of file FilterHilbertFFTPhase.h.

Referenced by v_Initialise(), and v_Update().

◆ h

Array<OneD, NekDouble> Nektar::FilterHilbertFFTPhase::h

private

Definition at line 107 of file FilterHilbertFFTPhase.h.

Referenced by v_Initialise(), and v_Update().

◆ linear_trans

bool Nektar::FilterHilbertFFTPhase::linear_trans = true

private

Definition at line 93 of file FilterHilbertFFTPhase.h.

Referenced by FilterHilbertFFTPhase(), and v_Update().

◆ m_fld

LibUtilities::FieldIOSharedPtr Nektar::FilterHilbertFFTPhase::m_fld

private

Definition at line 96 of file FilterHilbertFFTPhase.h.

Referenced by FilterHilbertFFTPhase(), and v_Update().

◆ m_index

unsigned int Nektar::FilterHilbertFFTPhase::m_index

private

Definition at line 85 of file FilterHilbertFFTPhase.h.

Referenced by FilterHilbertFFTPhase(), v_Initialise(), and v_Update().

◆ m_mean

double Nektar::FilterHilbertFFTPhase::m_mean

private

Definition at line 103 of file FilterHilbertFFTPhase.h.

Referenced by FilterHilbertFFTPhase(), and v_Update().

◆ m_outputFile

std::string Nektar::FilterHilbertFFTPhase::m_outputFile

private

Definition at line 95 of file FilterHilbertFFTPhase.h.

Referenced by FilterHilbertFFTPhase(), and v_Update().

◆ m_outputFrequency

unsigned int Nektar::FilterHilbertFFTPhase::m_outputFrequency

private

Definition at line 90 of file FilterHilbertFFTPhase.h.

Referenced by FilterHilbertFFTPhase(), and v_Update().

◆ m_outputIndex

unsigned int Nektar::FilterHilbertFFTPhase::m_outputIndex

private

Definition at line 86 of file FilterHilbertFFTPhase.h.

Referenced by FilterHilbertFFTPhase(), v_Initialise(), and v_Update().

◆ m_overlap

unsigned int Nektar::FilterHilbertFFTPhase::m_overlap

private

Definition at line 92 of file FilterHilbertFFTPhase.h.

Referenced by FilterHilbertFFTPhase(), v_Initialise(), and v_Update().

◆ m_TimeStep

NekDouble Nektar::FilterHilbertFFTPhase::m_TimeStep

private

Definition at line 100 of file FilterHilbertFFTPhase.h.

Referenced by v_Initialise(), and v_Update().

◆ m_window

int Nektar::FilterHilbertFFTPhase::m_window

private

Definition at line 91 of file FilterHilbertFFTPhase.h.

Referenced by FilterHilbertFFTPhase(), v_Initialise(), and v_Update().

◆ npoints

unsigned int Nektar::FilterHilbertFFTPhase::npoints

private

Definition at line 99 of file FilterHilbertFFTPhase.h.

Referenced by v_Initialise(), and v_Update().

◆ oldV

Array<OneD, NekDouble> Nektar::FilterHilbertFFTPhase::oldV

private

Definition at line 104 of file FilterHilbertFFTPhase.h.

Referenced by v_Initialise(), and v_Update().

◆ oldV_zero_mean

Array<OneD, NekDouble> Nektar::FilterHilbertFFTPhase::oldV_zero_mean

private

Definition at line 105 of file FilterHilbertFFTPhase.h.

Referenced by v_Initialise(), and v_Update().

◆ overlap_phase

Array<OneD, NekDouble> Nektar::FilterHilbertFFTPhase::overlap_phase

private

Definition at line 108 of file FilterHilbertFFTPhase.h.

Referenced by v_Initialise(), and v_Update().

◆ plan_backward

fftw_plan Nektar::FilterHilbertFFTPhase::plan_backward

private

Definition at line 111 of file FilterHilbertFFTPhase.h.

Referenced by v_Initialise(), and v_Update().

◆ plan_forward

fftw_plan Nektar::FilterHilbertFFTPhase::plan_forward

private

Definition at line 110 of file FilterHilbertFFTPhase.h.

Referenced by v_Initialise(), and v_Update().

◆ provided_mean

bool Nektar::FilterHilbertFFTPhase::provided_mean = false

private

Definition at line 94 of file FilterHilbertFFTPhase.h.

Referenced by FilterHilbertFFTPhase(), and v_Update().

◆ vCounter

unsigned int Nektar::FilterHilbertFFTPhase::vCounter

private

Definition at line 87 of file FilterHilbertFFTPhase.h.

Referenced by v_Initialise(), and v_Update().

Public Member Functions

Static Public Member Functions

Static Public Attributes

Protected Member Functions

Private Attributes

Friends

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ FilterHilbertFFTPhase()

◆ ~FilterHilbertFFTPhase()

Member Function Documentation

◆ create()

◆ v_Finalise()

◆ v_Initialise()

◆ v_IsTimeDependent()

◆ v_Update()

Friends And Related Function Documentation

◆ MemoryManager< FilterHilbertFFTPhase >

Member Data Documentation

◆ className

◆ fcoef

◆ h

◆ linear_trans

◆ m_fld

◆ m_index

◆ m_mean

◆ m_outputFile

◆ m_outputFrequency

◆ m_outputIndex

◆ m_overlap

◆ m_TimeStep

◆ m_window

◆ npoints

◆ oldV

◆ oldV_zero_mean

◆ overlap_phase

◆ plan_backward

◆ plan_forward

◆ provided_mean

◆ vCounter